Dear Samyak,
You second expression
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)
seems to me perfectly correct. Implicitely, the zero pressure
corresponds to the pressure outside the body (you have zero Neumann
conditions) so that what you obtain is a the difference of pressure
relatively to the ambient pressure. Do you find this unconsistent ?
Yves.
Le 10/02/2017 à 05:47, samyak jain a écrit :
> 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]
> <mailto:[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] <mailto:[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] <mailto:[email protected]>
>> https://mail.gna.org/listinfo/getfem-users
>> <https://mail.gna.org/listinfo/getfem-users>
>>
>>
>>
>>
>> _______________________________________________
>> Getfem-users mailing list
>> [email protected] <mailto:[email protected]>
>> https://mail.gna.org/listinfo/getfem-users
>> <https://mail.gna.org/listinfo/getfem-users>
>
> --
>
> Yves Renard ([email protected]
> <mailto:[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
> <http://math.univ-lyon1.fr/%7Erenard>
>
> ---------
>
--
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
---------
_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
