Modelling flexible inelastic materials

Hey @rob,

The TetrahedronFEMForceField implements a linear isotropic material, which means that it offers the same resistance to strain in every directions. It also means that it can only represent small and linear displacements. Since rotations are non-linear rigid motions, it will generate wrong elastic forces if your object is rotated. This is why you have an over-expansion.

When you set the method to “large”, it uses a corotational approach, which is still isotropic linear elasticity, but tries to remove the rotation before computing the elastic forces).

If you want anisotropy (more resistance to stretch in some given direction), or if you want a true non-linear material that is not sensible to rigid motions, you might have a look at the TetrahedronHyperelasticityFEMForceField component. For example,

<TetrahedronHyperelasticityFEMForceField ParameterSet="3448.2759 31034.483" AnisotropyDirections="0 1 0" materialName="StVenantKirchhoff"/>

Here the ParameterSet here uses, in order, the mu and lambda parameters instead of the Young modulus and Poisson ratio used in the TetrahedronFEMForceField. You can use this formulae to pass from one parameter set to the other:

mu     = young_modulus / (2. * (1. + poisson_ratio))
lambda = young_modulus * poisson_ratio / ((1. + poisson_ratio) * (1. - 2.*poisson_ratio))

Other examples can be found in the file “examples/Components/forcefield/TetrahedronHyperelasticityFEMForceField.scn”

J-N

Hey @rob,

The python code looks fine. You will probably have to set up correctly the material parameters (mu and lambda) so that it makes senses with the magnitude of your applied external forces . If you had found a Young modulus and Poisson ratio that worked well with the corotational model, start with these values using the equation I gave you to get the equivalent Lamé parameters.

Remove the anisotropy to make sure the crash doesn’t come from that.

Hyperelastic materials might also be harder to reach convergence, you can play your time step to hopefully find what is going on.

J-N

Hey @rob,

Does runSofa displays any error or warning before the crash? Can you paste its log output here?

The following python scene works well on my side:

import Sofa

poisson_ratio = 0.45
young_modulus = 1e4
mu = young_modulus / (2. * (1. + poisson_ratio))
lm = young_modulus * poisson_ratio / ((1. + poisson_ratio) * (1. - 2.*poisson_ratio))
    
def createScene(root):
    root.createObject('VisualStyle', displayFlags=['showForceFields', 'showBehaviorModels'])
    root.createChild('meca')
    root.meca.createObject('EulerImplicitSolver')
    root.meca.createObject('CGLinearSolver')
    
    root.meca.createObject('RegularGridTopology', name='grid', min=[0,0,0], max=[1,1,2.7], n=[3,3,8])
    root.meca.createObject('MechanicalObject', name='mo')
    root.meca.createObject('TetrahedronSetTopologyContainer', name='Container')
    root.meca.createObject('TetrahedronSetTopologyModifier')
    root.meca.createObject('TetrahedronSetTopologyAlgorithms')
    root.meca.createObject('Hexa2TetraTopologicalMapping', input="@./grid", output="@./Container" )
    root.meca.createObject('TetrahedronHyperelasticityFEMForceField', 
        ParameterSet=[mu, lm],  
        materialName="StVenantKirchhoff", 
        AnisotropyDirections=[0, 0, 1]
    )
    root.meca.createObject('BoxROI', drawBoxes=True, box=[-0.05, -0.05, -0.05, 1.05, 1.05, 0.05], name="fixed_box")
    root.meca.createObject('FixedConstraint', indices="@fixed_box.indices")
    root.meca.createObject('BoxROI', drawBoxes=True, box=[-0.05, -0.05, +2.60, 1.05, 1.05, 2.7], name="traction_box")
    root.meca.createObject('TriangleSetGeometryAlgorithms')
    root.meca.createObject('TrianglePressureForceField', 
        pressure=[0, -2e2, 0], 
        triangleList='@traction_box.triangleIndices', 
        showForces=True
    )