Homec4science
Diffusion tamaas (archive/legacy_mpi) Browse

tamaas/tests/test_hertz_kato.py
7c0aff50a114archive/legacy_mpi

Last Change
Raw File

test_hertz_kato.py
View Options

#!/usr/bin/env python
# coding: utf-8
# -----------------------------------------------------------------------------
# @author Lucas Frérot <lucas.frerot@epfl.ch>
#
# @section LICENSE
#
# Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de
# Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
# Solides)
#
# Tamaas is free software: you can redistribute it and/or modify it under the
# terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# Tamaas 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 Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
# -----------------------------------------------------------------------------
import sys
import numpy as np
import tamaas as tm
def constructHertzProfile(size, curvature):
radius = 1. / curvature
x = np.linspace(-0.5, 0.5, size)
y = np.linspace(-0.5, 0.5, size)
x, y = np.meshgrid(x, y)
surface = np.sqrt(radius**2 - x**2 - y**2)
surface -= surface.mean()
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 main():
grid_size = 512
curvature = 0.1
effective_modulus = 1.
load = 0.0001
surface = constructHertzProfile(grid_size, curvature)
bem = tm.BemKato(surface)
bem.setEffectiveModulus(effective_modulus)
bem.setMaxIterations(100000)
bem.computeEquilibrium(1e-7, load)
tractions = bem.getTractions()
displacements = bem.getDisplacements()
# Testing contact area against Hertz solution for solids of revolution
contact_area = tm.SurfaceStatistics.computeContactArea(tractions)
hertz_contact_size = (3 * load / (4 * curvature * effective_modulus))**(1. / 3.)
hertz_area = np.pi * hertz_contact_size**2
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 * load * effective_modulus**2 * curvature ** 2)**(1. / 3.) / np.pi
pressure_error = np.abs(hertz_max_pressure - max_pressure) / hertz_max_pressure
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, tractions[i, j] > 0) for i in range(grid_size) for j in range(grid_size)]
contact_indexes = map(lambda x: x[0:2], filter(lambda x: x[2], contact_indexes))
# 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 += 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] - 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)
if area_error > 6e-3 or pressure_error > 4e-3 or displacement_error > 3e-4:
return 1
return 0
if __name__ == "__main__":
sys.exit(main())
c4science · Help