# FEM and deformable objects questions

Hello!
I’m currently using create to simulate deformable objects and have some questions about this, someone can help me? thanks! ^__^

-What are the fundamental equilibrium equations implemented and solved by Create? Are they default Continuum Mechanics Cauchy-like equations? Do you have any references for those equations?

-What kind of solvers do you use? Explicit dynamic ones? Where can we find the formulation or even the references for the equations that have been implemented? Even just references to the names of the methods would help us.

-What kind of finite elements are implemented? First-order finite element shape functions? Is there any different FE implemented?

-How boundary conditions are applied? And how it relates to the different types of meshes?

-How does the constitutive/material behavior is implemented?

-What constitutive models are implemented? What nonlinear constitutive models are available? Can we select the ones that we want to use?

-What material laws are available (the one-dimensional relation between strains and stresses)? Can we tweak them? For instance, could we enter an experimental material strength curve?

-Where can we find the definition, with references to the actual equations, of each physics parameter of a deformable body?

-Is there any material or structural failure model implemented?

-On the UI there are basically two elastic/material parameters, the Young Modulus, and the Poisson. What does Create mean by elasticity damping and damping scale? If the constitutive models are nonlinear, there is a lack of additional phenomenological parameters to set up the equations. How does the simulation inherit the physical behavior?

-Why does the mesh refinement impact so much the physical behavior of the solid? The refinement does add additional degrees of freedom in the discrete models and enhances the kinematics of the FE model, but the impact is so big that it seems that the constitutive behavior is also influenced. In the end, it seems that the analyst needs to identify a trade-off between mesh refinement and material constants, which is not something usual. Not only that, even if one sets really high values for the elastic constants (i.e., an undeformable material), the deformation still would be expressive if the mesh is really fine. So, how does the implementation works underneath this behavior?

-Can we extract strains and stresses from the models? Would Create display a default FEA color map or even general state and dual variables values of the initial/boundary value problem

-Mass and damping seem to be major components in Create. There are some damping and mass parameters to be input by the user. We didn’t find technical information about those parameters. Do you have any reference for the general equilibrium equations and where those parameters are used in the formulation/implementation?

-What contact models are implemented? There is a dynamic friction parameter, so is the contact a frictional contact model? Is there any reference for the equations?

Hi

I finally managed to put together some references that should answer most of the questions. Since this is a lot of content and takes a considerable amount of time to summarize, I decided to directly write it in a format such that we can add it to the online PhysX documentation in the future (Welcome to PhysX — physx 5.1.2 documentation). Therefore some markdown syntax is still in the text.
Please note that we try to keep the Create documentation a bit more lightweight (as opposed to the PhysX SDK documentation which is mainly for software developers) since many people using it are not strongly interested in the underlying formulas.

The mathematical model is based on XPBD which is a semi implicit time stepping scheme. If the TGS solver is active, then a timestep is split into multiple smaller timesteps (dt_small=1/numIterations) which one iteration each. For the PGS solver, a full time step is done while the solution is refined until the maximal number of solver iterations is reached. The system of all motion constraints is solved in a Non-Linear Gauss-Seidel fashion where position updates of already solved constraints are fed directly into the solution of constraints that get solved later in the same (sub-)time step.

The following formulas and more detailed explanations can be found in

Both PGS and TGS solver update the soft body’s vertices according to constraints that ensure the correct physical behavior. [3] describes how TGS solves the constraint equations compared to previous methods and how this strategy improves the convergence rate. PGS follows the original approach for XPBD. The following section focuses on the constraints that define the material behavior.

Note that for soft bodies with a very high stiffness and many tetrahedral elements on the simulation mesh, the solver might have trouble to fully converge, even at a high iteration count. Therefore a soft body with many small tetrahedra might behave much softer than one with a low resolution tetmesh. Making the time step arbitrarily small is not ideal too since at some point 32 bit floating point numbers run into cancellation problems. Other floating point types with higher precision are not supported since they run quite a bit slower on GPUs.

Material Constraints:

There are two models implemented for strain calculation, both of them work on tetrahedral meshes: The first one is described in ([1], chapter 5.8.1) and is called Saint-Venant Kirchhoff model for isotropic materials. Anisotropy is not supported yet. The constraint formula is given by:

:math:C(F)=V \Psi_s(F)

:math:\Psi_s(F)=\frac{1}{2} tr(\epsilon^T S)

:math:S(F)=C \epsilon

The :math:C matrix in the formula above should not be confused with the constraint function as it represents the elasticity tensor of Hooke’s law in 3D.

:math:\epsilon(F)=\frac{1}{2}(F^T F - I)

where :math:F represents the deformation gradient of a tetrahedron, :math:S is Hooke’s law in 3D and :math:V is the volume of the tetrahedron in its rest state.

Note that the model is applied in a co-rotated local coordinate system to minimize artifacts due to linearization. The rotation matrix R can be obtained using polar decomposition. Check slides https://matthias-research.github.io/pages/publications/realtimeSlidesMatthias.pdf starting on page 43.

The second (recently) implemented model is the Neo-Hookean deviatoric constraint described in ([2], chapter 3.2.2). The constraint equation can be written as follows:

:math:C_D(F)=\sqrt{tr(F^T F)}

It is more difficult to achieve a high stiffness with this model due to convergence reasons. But it produces great results for soft rubber like materials and is incredible stable even under extreme deformations.

The second part to realistic soft body behavior is volume conservation. Especially rubber like materials keep their initial volume even if they get heavily deformed. A volume constraint based on a Neo-Hookean model as described using the term Hydrostatic Constraint in ([2], chapter 3.2.1) is used by PhysX. The constraint equation can be written as follows:

:math:C_H(F)=det(F)-1

where :math:F is the deformation gradient of a tetrahedron. Note that this formulation is independent of the orientation of the tetrahedron, therefore it does not matter if it is computed in global space or some rotate space (e. g. for a co-rotational model).

It is important to note that always a strain and a volume constraint are active per tetrahedron. Only the strain constraint can be chosen out of the two offered variants.

Stress reporting:

Optionally a stress measure can be requested for every tetrahedron of the collision mesh. Be aware that the stress that drives the deformation is actually computed on the simulation mesh which in most cases only coarsely approximates the mesh’s surface and can be of a completely different resolution than the collision tetmesh.
The stress measure is currently only implemented for the Saint-Venant Kirchhoff model. The stress calculation is based on the non-linear Green strain tensor:

:math:\epsilon(F)=\frac{1}{2}(F^T F - I)

Based on the strain tensor, Hooke’s law for isotropic materials in 3D can be applied. See Hooke's law - Wikipedia for more details and how to use the Lamé parameters (which can be computed out of Young’s modulus and Poisson’s ratio) to simplify its notation.

:math:\sigma=S(F)=2 \mu \epsilon + \lambda tr(\epsilon)I

where :math:\mu and :math:\lambda are the Lamé parameters. The stress reporting also computes the scalar von Mises stress value out of :math:\epsilon which is defined as follows:

:math:\sigma_v=\sqrt{3 J_2}

:math:J_2=\frac{1}{6} ((\sigma_{11}-\sigma_{22})^2+(\sigma_{22}-\sigma_{33})^2+(\sigma_{33}-\sigma_{11})^2)+\sigma_{12}^2+\sigma_{23}^2+\sigma_{31}^2

Note that :math:\sigma is a symmetrix 3x3 matrix.

The von Mises stress :math:\sigma_v can be used to evaluate the von Mises yield criterion (Cauchy stress tensor - Wikipedia) which is a common model to determine if a material can withstand the applied load or if it is likely to break.

Plasticity:

An experimental model for plasticity is implemented. It is not well tested and not guaranteed to be stable under all conditions. Material failure like tearing is not available.

Damping:

Damping values > 0 lead to energy dissipation inside of tetrahedral elements when they deform. Very high damping (e. g. 0.75) allows to create behaviors similar to memory foam.

DampingScale is an artificial (non physical) parameter in the range of 0 to 1 (default value 0) that allows to create certain special effects on a softbody like a sofbody that moves as it was filled with water.

Attachments:

Soft bodies can get attached to other objects like rigids. An explicit attachment constraint must be created to connect two objects. If the objects overlap, collision filtering needs to be active for the overlapping region because otherwise the collision detection routing will try to separate the two objects while the attachment constraint tries to keep the attached points at the same location leading to fighting constraints and jittering of the objects.

An alternative (and computationally very cheap) way which allows to fixate soft body vertices to a never changing location in world space is setting the inverse mass of every simulation mesh vertex that should never move to 0. The inverse mass is stored in the forth component (w) of every soft body vertex.

Contact Model:

Soft Bodies can collide with all other physics objects. The contact model is currently limited to support only dynamic friction. Static friction support might come in a future version. Dynamic friction uses a very simple model:

:math:f=\mu N

where :math:N is the normal force and :math:f a force that acts in the opposite direction of the tangential sliding motion between the two objects in contact.