public class AreaCompliantContact extends ContactForceBehavior
K = Ka*A in parallel with a damper of coefficient
D = max (stiffnessDamping, stepDamping*h) * Kwhere
Ka is a stiffness per unit area (with units of
force/length^3), A is the area associated with the contact, h is the simulation step size, stepDamping is a dimensionless
damping factor, and stiffnessDamping is a damping time constant.
Specifying stiffness per unit area, instead of per contact, makes the
contact pressure for a given penetration independent of the mesh
resolution.
This behavior is intended for vertex-based contacts (e.g., VERTEX_PENETRATION). The area A is given by the contactArea argument of computeResponse(double[], double, artisynth.core.modelbase.ContactPoint, artisynth.core.modelbase.ContactPoint, maspack.matrix.Vector3d, double, int), which for a vertex
contact is one third of the area of each face surrounding the penetrating
vertex, projected onto the plane perpendicular to the contact normal
(and halved for two-way contact).
Choosing Ka. Adding contact compliance is roughly equivalent to
adding springs to each contact, so that the resulting local deformation
delta (i.e., the penetration) at each contact point is about
p/Ka, where p is the contact pressure. If the maximum
contact pressure pmax and a maximum allowable delta are
known, one can choose Ka = pmax/deltaMax.
Alternatively, for contact involving a deformable body, Ka can
be set to try to control the ratio deltaB/delta, where deltaB is the deformation in the contact direction that the body itself
experiences under pressure p. If the body is assumed to be linearly
elastic with Young's modulus E and a Poisson's ratio of 0, then
deltaB is about p*Le/E, where Le is the
effective thickness of the body along the contact direction; i.e.,
the depth over which the contact stress stays concentrated over the contact
patch. Le lies between the true body thickness L and the
minimal width dp of the contact patch. With Ap the area of
the patch and AB the cross-sectional area of the body in the
contact plane, the following rules of thumb can be useful:
| body type | condition | Le |
|---|---|---|
| thin layer | L << dp | ~ L |
| compact | L >> dp, AB >> Ap | ~ dp |
| slender | L >> dp, AB ~ Ap | ~ L |
A nonzero Poisson's ratio will make the body stiffer, particularly for
a thin, nearly incompressible layer whose lateral expansion is constrained
(by a factor of about 3.8 for a Poisson's ratio of 0.45), in which case a
smaller value of Ka can be used. Given an estimate for Le,
Ka can then be chosen from
Ka = rho * E / Lewhere
rho is a stiffness ratio that equals deltaB/delta
under the assumptions given above. rho should be large enough that
delta is small relative to deltaB, but not so large that
the noisy, unevenly distributed contact pressures that can arise with
non-compliant vertex contact are reintroduced. Tests, while not
exhaustive, suggest that values of rho between 10 and 30 are
usually a good compromise. The default value of Ka corresponds to
rho = 20 for a cartilage-like layer with E = 5e6 and a
thickness of 2 mm (in SI units). It should usually be adjusted for other
materials, geometries and units.
Choosing the damping. Compliant contact forces are integrated
using the penetration estimated at the middle of the time step. As a
result, the damping time constant D/K must be at least h/2. If it is smaller, stiff contacts oscillate from one step to the
next, with contacts possibly switching on and off at every step, and the
oscillation becomes more persistent as Ka increases. The stepDamping property provides this minimum damping automatically, by
setting D/K to stepDamping*h, where h is the
current step size. This does not depend on the masses or units of the
model, and adjusts automatically if the step size is changed. A value of
0.5 damps the oscillation within a single step, while larger values give a
smoother but more sluggish response. The default value of 1 is suitable
for most models. Values less than 0.5 should be avoided unless stiffnessDamping is set to at least h/2.
Step damping is a numerical device and does not model physical
damping. In particular, when a light rigid body strikes a stiff compliant
contact, it rebounds with a coefficient of restitution of roughly 0.5/(0.5+stepDamping) (or about 1/3 for the default value), regardless
of the step size; increasing stepDamping reduces this rebound.
The stiffnessDamping property provides damping that is
independent of the step size, and is needed only when such damping is
required; for example, when the oscillation of a rigid body of mass
M on a soft contact, with frequency w = sqrt(Ktot/M)
(where Ktot is the total stiffness of all its contacts), is
resolved by the step size (i.e., w*h is less than about 1). A
damping ratio of about zeta can then be obtained by setting stiffnessDamping to 2*zeta/w. FEM models rarely require stiffnessDamping.
| Modifier and Type | Field and Description |
|---|---|
static PropertyList |
myProps |
TWO_WAY_CONTACT| Constructor and Description |
|---|
AreaCompliantContact()
Need no-args constructor for scanning
|
AreaCompliantContact(double kA)
Creates a new AreaCompliantContact with a specified stiffness per unit
area and default damping properties.
|
AreaCompliantContact(double kA,
double stepDamping,
double stiffnessDamping)
Creates a new AreaCompliantContact with specified properties.
|
| Modifier and Type | Method and Description |
|---|---|
void |
computeResponse(double[] fres,
double dist,
ContactPoint cpnt0,
ContactPoint cpnt1,
Vector3d normal,
double contactArea,
int flags)
Computes the force response for a given contact.
|
PropertyList |
getAllPropertyInfo()
Returns a list giving static information about all properties exported by
this object.
|
double |
getKa()
Queries the contact stiffness per unit area.
|
double |
getStepDamping()
Queries the step damping factor.
|
double |
getStiffnessDamping()
Queries the stiffness damping time constant.
|
void |
setKa(double ka)
Sets the contact stiffness per unit area, such that the stiffness of
each contact is
Ka times its associated area. |
void |
setStepDamping(double damping)
Sets the step damping factor, such that the damping coefficient of each
contact is at least
stepDamping*h times its stiffness, where
h is the current step size. |
void |
setStiffnessDamping(double tau)
Sets the stiffness damping time constant, such that the damping
coefficient of each contact is at least
stiffnessDamping times
its stiffness, independently of the step size. |
computeContactArea, getSubClasses, registerSubclassadvanceState, clone, createChangeEvent, createStateObject, getProperty, getPropertyHost, getPropertyInfo, hasProperty, hasState, hasSymmetricTangent, isWritable, postscan, scaleDistance, scaleMass, scan, setPropertyHost, setPropertyInfo, symmetryOrStateChanged, updateMaterial, write, writeItemspublic static PropertyList myProps
public AreaCompliantContact()
public AreaCompliantContact(double kA)
kA - stiffness per unit area (force/length^3)public AreaCompliantContact(double kA,
double stepDamping,
double stiffnessDamping)
kA - stiffness per unit area (force/length^3)stepDamping - step damping factorstiffnessDamping - stiffness damping time constantpublic PropertyList getAllPropertyInfo()
HasPropertiesgetAllPropertyInfo in interface HasPropertiesgetAllPropertyInfo in class MaterialBasepublic double getKa()
setKa(double).public void setKa(double ka)
Ka times its associated area. See the class
documentation for advice on choosing this value.ka - stiffness per unit area (force/length^3)public double getStepDamping()
setStepDamping(double).public void setStepDamping(double damping)
stepDamping*h times its stiffness, where
h is the current step size. The default value is 1, and values
less than 0.5 should usually be avoided; see the class documentation
for details.damping - step damping factorpublic double getStiffnessDamping()
setStiffnessDamping(double).public void setStiffnessDamping(double tau)
stiffnessDamping times
its stiffness, independently of the step size. The default value is 0;
see the class documentation for when a nonzero value is needed.tau - stiffness damping time constantpublic void computeResponse(double[] fres,
double dist,
ContactPoint cpnt0,
ContactPoint cpnt1,
Vector3d normal,
double contactArea,
int flags)
fres. Entries
that are not set default to 0. The force is defined as the restoring
force
force = f(dist)computed in response to the interpenetration distance
dist.
Note that in general, dist will be negative, and a
negative force will be expected in response. The compliance is the
reciprocal of the derivative of f(dist) evaluated at
dist. The damping is a positive damping parameter used to
generate damping forces in response to motion in the contact direction.
The step damping is an optional dimensionless factor that specifies a
minimum damping proportional to the step size h, such that the
effective damping d is given by
d = max (damping, stepDamping*h/compliance)This allows the damping time constant
d*compliance to be kept at
or above h/2, which is the minimum required for stable compliant
contact, without the behavior needing to know h. Damping and step
damping are applied only if the compliance is greater than 0.
Other information the can be used includes the two contact points associated the contact, and the contact normal, which is facing outward from the contact surface.
computeResponse in class ContactForceBehaviorfres - returns the force, compliance, damping and step damping for
a given contact distance.dist - contact interpenetration distance. This value is usually
negative, with greater negativity indicating greater interpenetration.cpnt0 - first contact pointcpnt1 - second contact pointnormal - contact normal, directed from cpnt0 to cpnt1contactArea - area associated with the contact. For vertex
penetration contacts, this is one third of the area of each face
surrounding the penetrating vertex; for edge-edge contacts, it is one
third of the area of each face adjacent to the first edge; in both
cases, the area is projected onto the plane perpendicular to the
normal. For contour region contacts, it is the area of the contact
region divided by the number of contact points. For two-way contact
(when vertex penetrations are computed for both collidables), the area
of vertex contacts is halved, since the vertices of both collidables
then cover the same contact patch.flags - flags providing additional information to the calculation.
No flags are currently defined.