Homec4science
Diffusion tamaas (feature/mpi_fft_engine) Browse

tamaas/tests/test_hertz_adhesion.py
87f2d0f8cfb1feature/mpi_fft_engine

Last Change
Raw File

test_hertz_adhesion.py
View Options

# -*- coding: utf-8 -*-
# -----------------------------------------------------------------------------
# @author Lucas Frérot <lucas.frerot@epfl.ch>
# @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/>.
# -----------------------------------------------------------------------------
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
#Other quantities
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 :
return 1
return 0
if __name__ == "__main__":
sys.exit(main())
c4science · Help