Modelling flexible inelastic materials

[rob]

Hi, what combination of Poisson Ratio and Young’s Modulus would you recommend for simulating a flexible but inelastic fabric?

Are there other parameters for the TetrahedronFEMForceField that would also be recommended?

In my experiments it always seems as though the material is stretching and so not quite what I’m after.

[Hugo]

Hi @robclark

Inelastic fabric is a knitting cloth right?
For this I must say I don’t know reference values.

Is the deformation of the tissue going to be large (>5-10%)? If so, you should pay attention to the field “method” of the TetrahedronFEMForceField and change it accordingly.

If you are not looking for stretch what are you looking for exactly?

Interesting topic anyway!
Best

Hugo

[rob]

Hi Hugo, many thanks for your reply.

The material I had in mind would be similar to the outer material of the inflatable kayak I have from Decathlon! For example: Decathlon Inflatable Kayak.

The idea was to have a flat envelope-like rectangular pouch which when inflated would take on a cylindrical shape. For the most part this works well but I see less contraction of the long edges of the envelope than expected. My guess is that the material is stretching a bit.

I will take a look at the “method” as you recommend.

Thank you.

[rob]

In case it’s helpful all my current experimental files are here:

https://drive.google.com/drive/folders/16M5DWp_7wWsTHtN3fDrcyuwlKBUwbuKp?usp=sharing

I tried changing from “large” to “small” and that caused over-expansion of the structure with almost no deformation of the rectangular perimeter.

[jnbrunet]

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

[rob]

Hi J-N,

Many thanks for your suggestion.

I have tried replacing my force-field with this, based on your XML:

pouch.createObject(‘TetrahedronHyperelasticityFEMForceField’, name = “FEM”,
ParameterSet = “3448.2759 31034.483”,
AnisotropyDirections = “0 1 0”,
materialName = “StVenantKirchhoff”)

However on a “Step” in SOFA I manage to crash the whole program. Obviously I’m doing something a little wrong… am I transliterating the XML->Python too naively?

[jnbrunet]

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

[rob]

Here is the code with the commented-out previous force-field and the new one. I’ve put the calculations in as you typed them, so hopefully they are correct. The Young’s Modulus and Poisson Ratio I had been using I had arrived upon through trial and error.

    young_modulus = 2
    poisson_ratio = 0.35
    # pouch.createObject('TetrahedronFEMForceField', template='Vec3d',
    #     name='FEM', method='large',
    #     poissonRatio = poisson_ratio,
    #     youngModulus = young_modulus,
    #     drawAsEdges="true")

    mu = young_modulus / (2. * (1. + poisson_ratio))
    lm = young_modulus * poisson_ratio / ((1. + poisson_ratio) * (1. - 2.*poisson_ratio))
    print(young_modulus, poisson_ratio, "->", mu, lm)

    pouch.createObject('TetrahedronHyperelasticityFEMForceField', name = "FEM",
        ParameterSet = [ mu, lm ],
        AnisotropyDirections = "0 1 0",
        materialName = "StVenantKirchhoff"
    )

I just get a SEGV on a single step, both from the SOFA I compiled on my MacBook and also the pre-compiled version I downloaded for the soft robotics toolkit (SOFA_v19.06.99_custom_MacOS_v11).

Have tried setting the AnisotropyDirections to all-zeros but the same SEGV.

Also changing the timestep to 0.5, 0.05, or 0.0001 all resulted in similar crashes.

Any ideas welcome!

[jnbrunet]

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
    )

[rob]

No specific warnings are shown in the console. The MacOS crash dump report that appears shows the stack location as:

Thread 0 Crashed:: Dispatch queue: com.apple.main-thread
0   libSofaMiscFem.19.06.99.dylib 	0x000000010117fe2c sofa::component::forcefield::TetrahedronHyperelasticityFEMForceField<sofa::defaulttype::StdVectorTypes<sofa::defaulttype::Vec<3, double>, sofa::defaulttype::Vec<3, double>, double> >::updateTangentMatrix() + 1132
1   libSofaMiscFem.19.06.99.dylib 	0x000000010117f816 sofa::component::forcefield::TetrahedronHyperelasticityFEMForceField<sofa::defaulttype::StdVectorTypes<sofa::defaulttype::Vec<3, double>, sofa::defaulttype::Vec<3, double>, double> >::addDForce(sofa::core::MechanicalParams const*, sofa::core::objectmodel::Data<sofa::helper::vector<sofa::defaulttype::Vec<3, double>, sofa::helper::CPUMemoryManager<sofa::defaulttype::Vec<3, double> > > >&, sofa::core::objectmodel::Data<sofa::helper::vector<sofa::defaulttype::Vec<3, double>, sofa::helper::CPUMemoryManager<sofa::defaulttype::Vec<3, double> > > > const&) + 310

Your example Python scene works perfectly on my build of SOFA. I will try to understand the differences between the two.

Thank you.

[rob]

The following items from your example are clearly missing from my code:

    root.meca.createObject('TetrahedronSetTopologyContainer', name='Container')
    root.meca.createObject('TetrahedronSetTopologyModifier')
    root.meca.createObject('TetrahedronSetTopologyAlgorithms')

If I just add them in I now get an error on loading the scene:

[ERROR] [TetrahedronHyperelasticityFEMForceField(FEM)] ERROR(TetrahedronHyperelasticityFEMForceField): object must have a Tetrahedral Set Topology.

In your example you use a Hexa2TetraTopologicalMapping whereas I thought I was already loading my object as tetras so don’t know whether I need a mapping or if not then what I need to do…

[jnbrunet]

Hey @rob,

The Hexa2TetraTopologicalMapping is there simply to create a tetrahedral topology from an hexahedral one (the regular grid). Hence you should not need it in your scene if you already have tetrahedrons.

You should however have a tetrahedral topology container (TetrahedronSetTopologyContainer). Make sure it is filled with tetrahedrons by double clicking on it and looking at the data field tetrahedra in the Property tab.

Also, open up the properties of the TetrahedronHyperelasticityFEMForceField (again by double clicking on it) and look at the data field *topology* in the Links tab. Does it point to your tetrahedron container? If not, you can force it with, for example:

root.meca.createObject('TetrahedronSetTopologyContainer', name='Container')
(...)
root.meca.createObject('TetrahedronHyperelasticityFEMForceField', 
    (...),
    topology='@./Container'
)

J-N

[Author]

Posts