Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/tests/test_hertz_disp.py
2316c9ebcdd7features/netcdf-mpi

Last Change
Raw File

test_hertz_disp.py
View Options

# -*- coding: utf-8 -*-
# @file
# @section LICENSE
#
# Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
from __future__ import print_function, division
import tamaas as tm
import numpy as np
def constructHertzProfile(size, curvature):
radius = 1. / curvature
x = np.linspace(-0.5, 0.5, size, dtype=tm.dtype)
y = np.linspace(-0.5, 0.5, size, dtype=tm.dtype)
x, y = np.meshgrid(x, y)
surface = np.sqrt(radius**2 - x**2 - y**2)
surface -= surface.min()
return surface.copy()
def computeHertzDisplacement(e_star, contact_size, max_pressure, size):
x = np.linspace(-0.5, 0.5, size)
y = np.linspace(-0.5, 0.5, size)
x, y = np.meshgrid(x, y)
disp = np.pi * max_pressure / (4 * contact_size * e_star) \
* (2 * contact_size**2 - (x**2 + y**2))
return disp.copy()
def test_hertz_disp(tamaas_fixture):
grid_size = 512
curvature = 0.5
effective_modulus = 1.
load = 0.12
surface = constructHertzProfile(grid_size, curvature)
model = tm.ModelFactory.createModel(tm.model_type_basic_2d,
[1., 1.],
[grid_size, grid_size])
model.setElasticity(1, 0)
solver = tm.PolonskyKeerRey(model, surface, 1e-12,
tm.PolonskyKeerRey.gap,
tm.PolonskyKeerRey.gap)
solver.solve(load - surface.mean())
tractions = model.getTraction()
true_displacements = model.getDisplacement()
true_gap = true_displacements - surface
pressure = np.mean(tractions)
contact_points = np.where(1.0e-10 >= true_gap)
contact_area = (np.size(contact_points)/2)/float(grid_size*grid_size)
print('The contact area computed with the gap map is '+str(contact_area))
hertz_contact_size = (3 * pressure
/ (4 * curvature * effective_modulus))**(1. / 3.)
print('Exact Hertz contact radius is '+str(hertz_contact_size))
hertz_area = np.pi * hertz_contact_size**2
print('Exact Hertz contact area is '+str(hertz_area))
area_error = np.abs(hertz_area - contact_area) / hertz_area
print("Area error: {}".format(area_error))
# Testing maximum pressure
max_pressure = tractions.max()
hertz_max_pressure = (6 * pressure * effective_modulus**2
* curvature ** 2)**(1. / 3.) / np.pi
pressure_error = np.abs(hertz_max_pressure - max_pressure) \
/ hertz_max_pressure
print('Exact Hertz max pressure is '+str(hertz_max_pressure))
print('max pressure is '+str(np.max(tractions)))
print("Max pressure error: {}".format(pressure_error))
# Testing displacements
hertz_displacements = computeHertzDisplacement(effective_modulus,
hertz_contact_size,
hertz_max_pressure,
grid_size)
# Selecing only the points that are in contact
contact_indexes = [(i, j)
for i in range(grid_size) for j in range(grid_size)
if true_gap[i, j] < 1e-12]
# Displacements of bem are centered around the mean of the whole surface
# and Hertz displacements are not centered, so we need to compute mean
# on the contact zone for both arrays
bem_mean = 0.
hertz_mean = 0.
for index in contact_indexes:
bem_mean += true_displacements[index]
hertz_mean += hertz_displacements[index]
bem_mean /= len(contact_indexes)
hertz_mean /= len(contact_indexes)
# Correction applied when computing error
correction = hertz_mean - bem_mean
# Computation of error of displacement in contact zone
error = 0.
hertz_norm = 0.
for index in contact_indexes:
error += (hertz_displacements[index] -
true_displacements[index] - correction)**2
hertz_norm += (hertz_displacements[index] - hertz_mean)**2
displacement_error = np.sqrt(error / hertz_norm)
print("Displacement error (in contact zone): {}".format(displacement_error))
assert area_error < 1e-2 \
and pressure_error < 1e-2 \
and displacement_error < 2e-3
if __name__ == "__main__":
test_hertz_disp()
c4science · Help