Dear Yves/Kostas, I think I understand what you are trying to tell me. Please find below an example on contact between a hyperelastic hemispherical rubber and a very rigid base. The displacement I am getting makes perfect sense but the hydrostatic pressure values are 10-50 times higher than what I am expecting.
Could you please have a look and tell me what am I doing wrong. I have also
attached the mesh files.
############################### Perfect Incompressible Case
####################################
import getfem as gf
import numpy as np
########### Physical parameters #############
# For Rubber/tool
young = 0.7967
poisson = 0.4999
toolOffset = -0.4 # Rubber Offset 0.4 mm
E1 = young # Young Modulus (N/mm^2)
nu1 = poisson # Poisson ratio
clambda1 = E1*nu1/((1+nu1)*(1-2*nu1)) # First Lame coefficient (N/mm^2)
cmu1 = E1/(2*(1+nu1)) # Second Lame coefficient (N/mm^2)
C10 = 0.23 # R1C10 - Mooney Rivlin 1st
Parameter - 0.23 MPa = 0.23 N/mm^2
C01 = 0.057 # R1C01 - Mooney Rivlin 2nd
Parameter - 0.057 MPa = 0.057 N/mm^2
# For Workpiece
E2 = 1e15 # Very high Young Modulus
(N/mm^2) --> Close to rigid body
nu2 = 0.3
clambda2 = E2*nu2/((1+nu2)*(1-2*nu2)) # First Lame coefficient (N/mm^2)
cmu2 = E2/(2*(1+nu2)) # Second Lame coefficient (N/mm^2)
########### Numerical parameters #############
elements_degree = 3 # Degree of the finite element methods
######### Mesh Import ###############
file_msh1 = "tool_hemisphere.msh"
file_msh2 = "workpiece1.msh"
mesh1 = gf.Mesh('import','gmsh',file_msh1)
mesh1.set('optimize_structure')
mesh2 = gf.Mesh('import','gmsh',file_msh2)
mesh2.set('optimize_structure')
########## Boundary selection ###########
P1 = mesh1.pts()
P2 = mesh2.pts()
top = (P1[2,:] + 1e-6) > 0.0
workTop = (P2[2,:] + 1e-6) > -10.0 # Only top surface of the workpiece
work = (P2[2,:] + 1e-6) > -1001.0 ### Get the entire mesh ####
pidtop = np.compress(top,range(0,mesh1.nbpts()))
pidworkTop = np.compress(workTop,range(0,mesh2.nbpts()))
pidwork = np.compress(work,range(0,mesh2.nbpts()))
ftop = mesh1.faces_from_pid(pidtop) # Tool-hemisphere top surface
fcontactWork = mesh2.faces_from_pid(pidworkTop) # Workpiece's top face
fAllWork = mesh2.faces_from_pid(pidwork) # Entire Workpiece's
fcontact = mesh1.outer_faces_with_direction([0., 0., -1.], np.pi/4.5)
# Contact boundary of the tool (rubber)
fbot = mesh2.outer_faces_with_direction([0., 0., -1.], 0.05) #
Workpiece's bottom face - Constrained by ground below
TOOL_TOP_BOUND=1; CONTACT_BOUND=2; CONTACT_BOUND2=3; BOTTOM_BOUND=4; ALL_WORK=5;
########## Body - 1 ###############
# Definition of finite elements methods and integration method
mfu1 = gf.MeshFem(mesh1, 3)
mfu1.set_classical_fem(2)
pre_mflambda1 = gf.MeshFem(mesh1, 3)
pre_mflambda1.set_classical_fem(1)
# FEM for Stress Tensor
mfPKST = gf.MeshFem(mesh1)
mfPKST.set_classical_discontinuous_fem(2)
# Set few regions for the mesh
mesh1.set_region(TOOL_TOP_BOUND, ftop)
mesh1.set_region(CONTACT_BOUND, fcontact)
# Integration Methods
mim1 = gf.MeshIm(mesh1, 4)
mim1_contact = gf.MeshIm(mesh1, 4)
######## Body - 2 ##################
mfu2 = gf.MeshFem(mesh2, 3)
mfu2.set_classical_fem(2)
# Set few regions for the mesh
mesh2.set_region(CONTACT_BOUND2, fcontactWork)
mesh2.set_region(BOTTOM_BOUND, fbot)
mesh2.set_region(ALL_WORK, fAllWork)
# Integration Methods
mim2 = gf.MeshIm(mesh2, 4)
mim2_contact = gf.MeshIm(mesh2, 4)
######## Model definition ###############
md=gf.Model('real');
md.add_fem_variable('u1', mfu1) # Displacement of the structure 1
md.add_filtered_fem_variable('lambda1', pre_mflambda1, CONTACT_BOUND)
# Add data of the model and add a elasticity brick for the rubber
md.add_initialized_data('cmu1', [cmu1])
md.add_initialized_data('clambda1', [clambda1])
md.add_initialized_data('mrParams', [0.23, 0.057])
# md.add_initialized_data('mrParams', [0.23])
md.add_finite_strain_elasticity_brick(mim1, 'Incompressible Mooney
Rivlin', 'u1', '[0.23, 0.057]')
md.add_fem_variable('u2', mfu2) # Displacement of the structure 2
# Add the Lame coefficients as data of the model and add a linearized
elasticity brick for the workpiece
md.add_initialized_data('cmu2', [cmu2])
md.add_initialized_data('clambda2', [clambda2])
md.add_isotropic_linearized_elasticity_brick(mim2, 'u2', 'clambda2', 'cmu2')
# Large Sliding Contact condition
md.add_initialized_data('r', gamma0)
indbrick = md.add_integral_large_sliding_contact_brick_raytracing('r', 1.0)
md.add_slave_contact_boundary_to_large_sliding_contact_brick(indbrick,
mim1_contact, CONTACT_BOUND, 'u1', 'lambda1')
md.add_master_contact_boundary_to_large_sliding_contact_brick(indbrick,
mim2_contact, CONTACT_BOUND2, 'u2')
# md.add_rigid_obstacle_to_large_sliding_contact_brick(indbrick, 'z+6', 3)
# Add Incompressibility brick as we are using incompressible mooney-rivlin law
mfp = gf.MeshFem(mesh1,1)
mfp.set_classical_discontinuous_fem(1)
md.add_fem_variable('p', mfp) # Hydrostatic pressure multiplier term
md.add_finite_strain_incompressibility_brick(mim1, 'u1', 'p')
# Add dirichlet displacement boundary conditions on the rubber surface
(TOOL_TOP_BOUND) as there is no force condition
md.add_initialized_data('toolOffset', [0.0,0.0,toolOffset])
md.add_Dirichlet_condition_with_multipliers(mim1, 'u1', 1,
TOOL_TOP_BOUND, 'toolOffset')
# Add dirichlet boundary conditions on the workpiece surface (BOTTOM_BOUND)
md.add_initialized_data('fixedBottom', [0.0,0.0,0.0])
md.add_Dirichlet_condition_with_multipliers(mim2, 'u2', 1, -1,
'fixedBottom') # Entire Mesh
# md.add_Dirichlet_condition_with_multipliers(mim2, 'u2', 1,
BOTTOM_BOUND, 'fixedBottom')
######### Model solve ##############
print 'Solve problem with ', md.nbdof(), ' dofs'
md.solve('max_res', 1E-7, 'max_iter', 20, 'noisy' , 'lsearch',
'simplest', 'alpha min', 0.8)
# Solution export
U1 = md.variable('u1') # Values are correct
U2 = md.variable('u2') # Values are correct
Pr = md.variable('p') # Hydrostatic pressure - Values are 10 to 50
times large and both positive/negative
####### Pressure and Stress Calculations ##############
# Hydrostatic Pressure Calculations from Cauchy Stress (The values are zero)
pressureMethod1 =
md.interpolation("-0.33*Trace(Cauchy_stress_from_PK2(Incompressible_Mooney_Rivlin_sigma(Grad_u1,
[0.23, 0.057]),Grad_u1))", mfPKST, -1)
# Hydrostatic Pressure Calculations from Cauchy Stress after adding
the incompressible pressure part (The values are 10 times higher and
are both positive and negative)
pressureMethod2 =
md.interpolation("-0.33*Trace(Cauchy_stress_from_PK2(Incompressible_Mooney_Rivlin_sigma(Grad_u1,
[0.23, 0.057]),Grad_u1) - p*Id(3) )", mfPKST, -1)
slice = gf.Slice(('planar',0,[[0],[0],[-10+abs(toolOffset)]],[[0],[0],[1]]),
mfu1, 3)
pressureOnSlice = gf.compute_interpolate_on(mfPKST, pressure, slice)
slice.export_to_pos('slice.pos', mfPKST, pressure,
'Pressure(N/mm^2)', mfu1, U1, 'Displacement(mm)')
Thanking You
Yours sincerely
Samyak Jain
On Thu, Feb 9, 2017 at 5:40 PM, Yves Renard <[email protected]>
wrote:
>
> Dear Samyak,
>
> May be to be more precise, I would add that if you use the finite strain
> incompressible law, then the multiplier (call it lambda) used to prescribe
> the incompressiblity constraint is to be added to the spheric part of the
> Hyperelastic Law used
>
> sigma = Hyperelastic_law_used - lambda Id(3)
>
> So that the hydrostatic pressure is -1/3*Tr(sigma) contains two terms.
>
> Yves.
>
> Le 09/02/2017 à 09:26, Konstantinos Poulios a écrit :
>
> Dear Samyak,
>
> It is definitely not true that the hydrostatic term -1/3*Tr(sigma) should
> be zero in an incompressible material. If this is the case you are simply
> calculating the deviatoric part of sigma instead of sigma. In order to get
> sigma you need to add the term p*Id(3). Actually the hydrostatic term you
> are looking for is simply equal to p.
>
> Best regards
> Kostas
>
>
>
> On Thu, Feb 9, 2017 at 4:20 AM, samyak jain <[email protected]>
> wrote:
>
>> Dear getfem-users,
>>
>> I am currently trying to solve a contact problem between a hyperelastic
>> rubber and a rigid bosy and I need to calculate the pressure values on
>> either on the rubber.
>>
>> I am using Incompressible Mooney-Rivlin Hyperelastic law and if my model
>> I am adding also adding finite strain incompressibility brick.
>>
>> Now when I calculate Cauchy Stress from second piola kirchhoff stress, I
>> am getting the Hydrostatic term (-1/3Tr(sigma)) of the cauchy stress tensor
>> as zero which is what it should be as the material is incompressible.
>>
>> So, is there a way is getfem to calculate the hydrostatic pressure term
>> for such incompressible materials.I believe treating the material as
>> nearly incompressible (Poisson's ratio 0.499) is one way to solve it but I
>> don't know how it works or if it is implemented in the model.
>>
>> Could you guys please provide any help or suggestion to calculate the
>> hydrostatic pressure for such a case.
>>
>> Thanks a lot.
>>
>> Yours sincerely
>> Samyak
>>
>> _______________________________________________
>> Getfem-users mailing list
>> [email protected]
>> https://mail.gna.org/listinfo/getfem-users
>>
>>
>
>
> _______________________________________________
> Getfem-users mailing
> [email protected]https://mail.gna.org/listinfo/getfem-users
>
>
>
> --
>
> Yves Renard ([email protected]) tel : (33) 04.72.43.87.08
> Pole de Mathematiques, INSA-Lyon fax : (33) 04.72.43.85.29
> 20, rue Albert Einstein
> 69621 Villeurbanne Cedex, FRANCE
> http://math.univ-lyon1.fr/~renard
>
> ---------
>
>
test.py
Description: Binary data
tool_hemisphere.msh
Description: Mesh model
workpiece1.msh
Description: Mesh model
_______________________________________________ Getfem-users mailing list [email protected] https://mail.gna.org/listinfo/getfem-users
