Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/tests/test_hertz_adhesion.py
586c9bd80954features/netcdf-mpi

Last Change
Raw File

test_hertz_adhesion.py
View Options

# -*- coding: utf-8 -*-
# -----------------------------------------------------------------------------
# @author Valentine Rey <valentine.rey@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/>.
# -----------------------------------------------------------------------------
from __future__ import print_function, division
import sys
import tamaas as tm
import numpy as np
import scipy.optimize as scio
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 = -(x**2 + y**2)/(2*radius)
return surface.copy()
def main():
grid_size = 512
curvature = 0.5
effective_modulus = 1.
surface_energy = 0.0001
rho = 0.002
radius = 1./curvature
K = 4*effective_modulus/3
KI = np.sqrt(surface_energy*2*effective_modulus)
tabor = surface_energy**(2/3)/(rho*effective_modulus**(2/3))*(1/curvature)**(1/3)
sig0 = surface_energy/(rho)
lam = 2 * sig0/(np.pi*surface_energy*K**2/radius)**(1/3)
surface = constructHertzProfile(grid_size, curvature)
bem = tm.BemGigipol(surface)
bem.setDumpFreq(100)
functional = tm.MaugisAdhesionFunctional(bem)
functional.setParameter('rho', rho)
functional.setParameter('surface_energy', surface_energy)
bem.setEffectiveModulus(effective_modulus)
bem.addFunctional(functional)
bem.computeEquilibrium(1e-8, 0.00116480372339)
tractions = bem.getTractions()
true_gap = bem.getGap()
true_displacements = bem.getTrueDisplacements()
adhesive_pressure = np.zeros((grid_size, grid_size))
for i in np.arange(0, grid_size, 1):
for j in np.arange(0, grid_size, 1):
if (rho > true_gap[i, j]):
adhesive_pressure[i, j] = surface_energy/rho
# computed solution
p0 = np.mean(tractions)
p_max = np.max(tractions)
contact_points = np.where(true_gap == 0)
contact_area = np.size(contact_points)/2/(float(grid_size)**2)
ag = (contact_area/np.pi)**0.5
adh_points = np.where(adhesive_pressure != 0.)
adh_area = np.size(adh_points)/2/(float(grid_size)**2)
cg = (adh_area/np.pi)**0.5
delta = ag**2/radius-8*sig0/(3*K)*np.sqrt(cg**2-ag**2)
# Exact solution
beta = (1-np.exp(lam/-0.924))/1.02
a_chap = 1.54+((2.28*lam**1.3-1)/(2.28*lam**1.3+1))*0.279
L_chap = -7./4.+((4.04*lam**1.4-1)/(4.04*lam**1.4+1))/4.
a0 = a_chap*(np.pi*surface_energy*radius**2/K)**(1./3.)
pc = L_chap*np.pi*surface_energy*radius
P = -0.00078
ac = a0*((beta+np.sqrt(1-P/pc))/(1+beta))**(2./3.)
def fonction1(m):
return np.sqrt(m**2-1)+m**2*np.arctan(np.sqrt(m**2-1))-np.pi*KI/(np.sqrt(np.pi*ac)*sig0)
m = scio.brenth(fonction1, 1., 2.)
cc = m*ac
delta_cos = ac**2/radius-8*sig0/(3*K)*np.sqrt(cc**2-ac**2)
x = np.linspace(-0.5, 0.5, grid_size)
y = np.linspace(-0.5, 0.5, grid_size)
x, y = np.meshgrid(x, y)
a = ac
c = cc
# stress
sig_hertz = np.zeros((grid_size, grid_size))
disp = np.zeros((grid_size, grid_size))
for i in np.arange(0, grid_size, 1):
for j in np.arange(0, grid_size, 1):
if a**2>(x[i,j]**2 + y[i,j]**2):
disp[i,j]=3*K/(2*np.pi*radius)*np.sqrt(a**2-(x[i,j]**2 + y[i,j]**2))-2*sig0/np.pi*np.arctan(np.sqrt((c**2-a**2)/(a**2-(x[i,j]**2 + y[i,j]**2))))
if c**2>(x[i,j]**2 + y[i,j]**2)>a**2:
disp[i,j]=-sig0
sig_hertz = disp
a_error = abs((ac-ag)/ac)
c_error = abs((cc-cg)/cc)
delta_error = abs((delta-delta_cos)/delta)
pressure_error = abs((p0-np.mean(sig_hertz))/np.mean(sig_hertz))
if a_error > 1.e-2 or c_error > 1.2e-2 or delta_error > 5.e-2 or pressure_error > 1.e-1:
for err in [a_error, c_error, delta_error, pressure_error]:
print(err)
return 1
return 0
if __name__ == "__main__":
sys.exit(main())
c4science · Help