[SOLVED] How to stabilize the beam model using BeamFEMForceField?

[shang]

Hi, all
I using using the BeamFEMForceField component to simulate a beam-like catheter. However, it turns out that the model is not stable at all, could anyone help me figure out what the problem is?
Here is the scene file:

<Node name="root" dt="0.01">
    <VisualStyle displayFlags="hideBehaviorModels hideCollisionModels hideMappings hideForceFields" />
    <RequiredPlugin name="Geomagic plugin" pluginName="Geomagic" />
    <RequiredPlugin pluginName="SofaPython"/>
    <BackgroundSetting color="0.8 0.4 0.6"/>
    
    <CollisionPipeline name="pipeline" depth="6" verbose="0"/>
    <BruteForceDetection name="detection" />
    <CollisionResponse name="response" response="FrictionContact" />
    <LocalMinDistance name="proximity" alarmDistance="0.15" contactDistance="0.05" angleCone="0.0" />
    <FreeMotionAnimationLoop/>
    <LCPConstraintSolver tolerance="0.001" maxIt="1000"/>
    <GeomagicDriver name="GeomagicDevice" deviceName="Default Device" scale="1" drawDeviceFrame="1" positionBase="0 0 0" orientationBase="0 0.707 0 -0.707"  />

    <Node name="Omni">
	 <MechanicalObject template="Rigid" name="DOFs" position="@GeomagicDevice.positionDevice"/>
         <MechanicalStateController template="Rigid" listening="true" mainDirection="0.0 0.0 1.0" handleEventTriggersUpdate="true"/>
	 <PythonScriptController name="device" filename="E:/sofa-engine/sofa/applications/plugins/Geomagic/scenes/ScriptEvents.py" classname="geomagic"/>
    </Node>

    <Node name="beam" gravity="0 9.8 0">
       <EulerImplicit rayleighStiffness="0.1" printLog="false" />
       <CGLinearSolver template="GraphScattered" iterations="20" threshold="1e-008" tolerance="1e-5" />
       <MechanicalObject template="Rigid" name="catheterDOFs" rotation="0 -90 0"
		position="
					0 -1 0 0 0 0 1
					0.5 -1 0 0 0 0 1
					1 -1 0 0 0 0 1
					1.5 -1 0 0 0 0 1
					2 -1 0 0 0 0 1
					2.5 -1 0 0 0 0 1
					3 -1 0 0 0 0 1
					3.5 -1 0 0 0 0 1
					4 -1 0 0 0 0 1
					4.5 -1 0 0 0 0 1
					5 -1 0 0 0 0 1
					5.5 -1 0 0 0 0 1
					6 -1 0 0 0 0 1
					6.5 -1 0 0 0 0 1
					7 -1 0 0 0 0 1
					7.5 -1 0 0 0 0 1
					8 -1 0 0 0 0 1
					8.5 -1 0 0 0 0 1
					9 -1 0 0 0 0 1
					9.5 -1 0 0 0 0 1
					10 -1 0 0 0 0 1
					10.5 -1 0 0 0 0 1
					11 -1 0 0 0 0 1
					11.5 -1 0 0 0 0 1
					12 -1 0 0 0 0 1
					12.5 -1 0 0 0 0 1
					13 -1 0 0 0 0 1
					13.5 -1 0 0 0 0 1
					14 -1 0 0 0 0 1
					14.5 -1 0 0 0 0 1
					15 -1 0 0 0 0 1
					15.5 -1 0 0 0 0 1
					16 -1 0 0 0 0 1
					16.5 -1 0 0 0 0 1
					17 -1 0 0 0 0 1
					17.5 -1 0 0 0 0 1
					18 -1 0 0 0 0 1
					18.5 -1 0 0 0 0 1
					19 -1 0 0 0 0 1
					19.5 -1 0 0 0 0 1
					20 -1 0 0 0 0 1
					20.5 -1 0 0 0 0 1
					21 -1 0 0 0 0 1
					21.5 -1 0 0 0 0 1
					22 -1 0 0 0 0 1
					22.5 -1 0 0 0 0 1
					23 -1 0 0 0 0 1
					23.5 -1 0 0 0 0 1
					24 -1 0 0 0 0 1
					24.5 -1 0 0 0 0 1
					25 -1 0 0 0 0 1
					25.5 -1 0 0 0 0 1
					26 -1 0 0 0 0 1
					26.5 -1 0 0 0 0 1
					27 -1 0 0 0 0 1
					27.5 -1 0 0 0 0 1
					28 -1 0 0 0 0 1
					28.5 -1 0 0 0 0 1
					29 -1 0 0 0 0 1
					29.5 -1 0 0 0 0 1
					" />
        <Mesh name="lines" lines="0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 21 21 22 22 23 23 24 24 25 25 26 26 27 27 28 28 29 29 30 30 31 31 32 32 33 33 34 34 35 35 36 36 37 37 38 38 39 39 40 40 41 41 42 42 43 43 44 44 45 45 46 46 47 47 48 48 49 49 50 50 51 51 52 52 53 53 54 54 55 55 56 56 57 57 58 58 59" />
	<Monitor name="catheterInfo" template="Rigid" listening="1" indices="0 1 2 3 4" showPositions="0" PositionsColor="1 1 0 1"  showMinThreshold="0.01" TrajectoriesPrecision="0.1" TrajectoriesColor="1 1 0 1" sizeFactor="1" />
		<FixedConstraint name="FixedConstraint" indices="0" />
		<UniformMass totalmass="6"  />
        <BeamFEMForceField name="FEM" radius="10" youngModulus="2000000" poissonRatio="0.9"/>
		<UncoupledConstraintCorrection />
		<ConstantForceField />
        <Node name="Collision">
	    <CylinderGridTopology name="coli" nx="6" ny="6" nz="59" length="59" radius="0.3" axis="1 0 0" center="0 -1 0"/>		
            <MechanicalObject />
            <BeamLinearMapping isMechanical="true" />
            <Triangle />
            <Line />
	    <Point />
        </Node>
	<Node name="visu">
	   <CylinderGridTopology name="coli" nx="6" ny="6" nz="59" length="59" radius="0.3" axis="1 0 0" center="0 -1 0"/>
	   <OglModel color="gray" />
	   <BeamLinearMapping isMechanical="true" />
        </Node>
	<PythonScriptController name="catheter" filename="ScriptEvents.py" classname="catheter"/>
    </Node>	
    <Node name="clamper" gravity="0 9.8 0">		
        <Node name="CollisionModel" >
            <MeshObjLoader filename="mesh/clamper.obj"  name="loader" scale="10" />
            <Mesh src="@loader" name="colision"  />
            <MechanicalObject src="@loader" name="instrumentCollisionState" />
            <Triangle contactStiffness="10" />
            <Line contactStiffness="10"  />    
            <Point contactStiffness="10" />      
        </Node>
        <Node name="VisualModel" >
            <OglModel name="visual" fileMesh="mesh/clamper.obj" color="blue" scale="10" />       
        </Node>
</Node>

[Hugo]

Hey @xiaojuan

Are you choosing all numerical settings and model parameters on purpose?
Could you especially explain the Poisson ration = 0.9 ?

For sure, more stability would be obtained by using a mass that is based on the Beam FEM. The theory for the Beam in SOFA is the Timoshenko beam theory. Here you are using the UniformMass, i.e. the same mass value at each point of your beams which should not be the case. A mass matrix with a correct integration is a tri-diagonal matrix taking into account the section of the beam and the length of it. Unfortunately such a Beam mass does not exist yet in SOFA. I know a branch was working on this here.

Moreover you can find papers on Timoshenko beam theory online. Like here or any other article. We would be really glad to welcome your contribution on this.

Best

Hugo

[shang]

Nice,Hugo,
Thank you very much! I will read these articles carefully.

[Hugo]

Ok I close this topic.
In case, you have further question on Beam stability, do not hesitate to ask it here.
Or create another thread for any other question.

Best

Hugo

[Author]

Posts