Metal Forging Simulation

Hi all,

I’m interested in using SOFA for the simulation of metal forging using a power hammer. My aim is to use this simulation within the initial stages of a deep reinforcement learning model. Therefore accuracy is less of a concern, it just needs to be fast. My background is in computer science and so I have little to no experience with finite element analysis or using simulation software. I’ve managed to piece together something which seems close to what I need but it doesn’t quite work correctly.

In my simulation I have a stationary bottom die, a rigid top die and work piece in the middle which is modelled using a TetrahedronFEMForceField for plastic behaviour. The top die is supposed to accelerate into the work piece in order to plastically deform it. I’m using parameters to try and model a piece of copper at room temperature (cold forming). I want the work piece to undergo plastic deformation due to the force applied by the top die. At the moment this does not work correctly. The work piece is showing significant elastic deformation when the top die has a mass of only 10kg. This is not realistic behaviour. Another issue is that any deformation seems to be happening in the y direction – there is no expansion of the work piece in the x or z directions which would be expected under a compressive load due to the Poisson effect.

I think my issues might be coming from an incorrect collision model however I’m unable to find out how to correct my problems using the documentation. I’m assuming all standard SI units so my work piece mesh has dimensions of 0.1m x 0.1m x 0.1m. Does anybody have any experience using SOFA for FEM simulations of deforming metals due to large amounts of stress or does anyone know if this is possible using SOFA?

Here is my code:

<Node name="root" dt="0.0001" gravity="0 -9.81 0">
    <RequiredPlugin name="SofaOpenglVisual"/>
	<RequiredPlugin name="SofaMiscCollision"/>
	<RequiredPlugin  name="SofaPython" printLog="1"  pluginName="SofaPython" />
	<VisualStyle displayFlags="showCollisionModels showForceFields" />
   
    <LCPConstraintSolver tolerance="1e-6" maxIt="1000" initial_guess="false" build_lcp="false"  printLog="0" mu="0.2"/>
    <FreeMotionAnimationLoop name="FreeMotionAnimationLoop" />
	<DefaultPipeline depth="6" verbose="0" draw="0" />
	<BruteForceDetection name="N2" />
	<LocalMinDistance name="Proximity" alarmDistance="0.01" contactDistance="0.005" />
	<DefaultContactManager response="FrictionContact" />
    
	<MeshGmshLoader name="top_loader" filename="top_die.msh" createSubelements="true" translation="0.0125 .10 0.025" />
	<MeshGmshLoader name="work_piece_loader" filename="work_piece.msh" createSubelements="true" translation="0 0.05 0"/>
	<MeshGmshLoader name="bottom_loader" filename="bottom_die.msh" createSubelements="true" translation="0 -.1 0" />
	<PythonScriptController name="Hammer" filename="hammer.py" classname="Hammer" printLog="true"/>

	<Node name="Top">
		<EulerImplicitSolver name="cg_odesolver" printLog="false"  rayleighStiffness="0.1" rayleighMass="0.1" />
		<CGLinearSolver iterations="500" name="linear solver" tolerance="1.0e-15" threshold="1.0e-15" />
		
		<MeshTopology src="@../top_loader"/>
		<MechanicalObject template="Rigid3d" name="mo_top"/>

		<UniformMass totalMass="10" />		
		<ConstantForceField name="const" force="0 0 0 0 0 0 0"/>
		<PartialFixedConstraint fixedDirections="0 0 0 1 1 1"/>
		<UncoupledConstraintCorrection compliance=".5" />
		
		<Node name="VisuTop" tags="Visual">
			<OglModel name="Visual" src="@../../top_loader" color="blue" />
			<RigidMapping input="@../mo_top" output="@Visual" />
		</Node>
		<Node name="Collision">
			<MeshTopology src="@../../top_loader"/>
			<MechanicalObject name = "mo_top_collision" />
			<LineCollisionModel/>
			<RigidMapping input="@../mo_top" output="@mo_top_collision"/>
		</Node>
	</Node>	

	<Node name="Work Piece">
		<EulerImplicitSolver name="cg_odesolver" printLog="false"  rayleighStiffness="0.1" rayleighMass="0.1" />
		<CGLinearSolver iterations="500" name="linear solver" tolerance="1.0e-15" threshold="1.0e-15"/>

		<TetrahedronSetTopologyContainer name="TetWorkPiece" src="@../work_piece_loader" />
		<MechanicalObject name="mo_work_piece" />
		<TetrahedronSetGeometryAlgorithms name="GeomAlgo" />
		<TetrahedronFEMForceField name="FEM" method="large" poissonRatio="0.36" youngModulus="117e9" computeGlobalMatrix="false" plasticYieldThreshold="0.0006" plasticMaxThreshold="0.0019" plasticCreep=".9"  />
		
		<UniformMass totalMass="8.96" />
		<UncoupledConstraintCorrection compliance=".5"/>
		
		<Node name="VisuWorkPiece" tags="Visual">
			<OglModel name="Visual" src="@../../work_piece_loader" color="blue" />
			<IdentityMapping input="@../mo_work_piece" output="@Visual" />
		</Node>
		<Node name="Collision">
			<MeshTopology src="@../../work_piece_loader"/>
			<MechanicalObject name = "mo_work_piece_collision"/>
			<LineCollisionModel />
			<IdentityMapping input="@../mo_work_piece" output="@mo_work_piece_collision"/>
		</Node>
	</Node>	

	<Node name="Bottom">
		<MeshTopology src="@../bottom_loader"/>
		<MechanicalObject name="mo_bot"/>

		<Node name="VisuBottom" tags="Visual">
			<OglModel name="Visual" src="@../../bottom_loader" color="blue" />
			<IdentityMapping input="@../mo_bot" output="@Visual" />
		</Node>
		<Node name="Collision">
			<MeshTopology src="@../../bottom_loader"/>
			<MechanicalObject name = "mo_bottom_collision" />
			<LineCollisionModel  simulated="0" moving="0" contactFriction="0.01" />
			<IdentityMapping input="@../mo_bot" output="@mo_bottom_collision"/>
		</Node>
	</Node>
</Node>

I have a python script called hammer.py which lifts up the top die after a certain amount of time. This is to crudely mimic the operation of a power hammer:

import Sofa
import sys    

class Hammer(Sofa.PythonScriptController):
     
    def initGraph(self,node):
        top_die_node = node.getChild('Top')
        self.constant_force_field = top_die_node.getObject('const')
        self.time = 0

        return 0

    def onBeginAnimationStep(self,dt):
        self.time += dt
        if self.time > .17:
            self.constant_force_field.findData('force').value = "0 5000 0"
      
        return 0

The 3 meshes have been created using gmsh and can be found along with the rest of the code at : files