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
- http://mmacklin.com/2017-EG-CourseNotes.pdf
- https://mmacklin.com/neohookean.pdf
- http://mmacklin.com/smallsteps.pdf
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.
Hope that answers most of your questions