Damping Options

NairnMPM has several forms of artificial damping done during particle updates. Their most common use is to damp dynamic effects such that the solution converges to a static result. Some materials also support artificial viscosity as a damping mechanism done on the material model level. In addition, NairnMPM supports FLIP, PIC, and combinations of FLIP and PIC by representing PIC damping of FLIP.

Introduction

Many times it is useful to apply damping to results, such as when the goal to achieve a simulation of quasi-static results. This section explains various artificial damping options that are applied to the particle update in the MPM time step. The particle updates for position and velocity with damping become:

      [math]\displaystyle{ \vec x_p^{(n+1)} = \vec x_p^{(n)} + \left(\vec v_{g\to p}^{(n+1)} - {\Delta t\over 2}\biggl(\vec a_{g\to p}^{(n)} + \alpha_g \vec v_{g\to p}^{(n)} + \alpha_p \vec v_{p}^{(n)} + \alpha(\beta)\bigl(\vec v_{p}^{(n)} -\vec v_{g\to p}^{(n)} \bigr) + \alpha(\beta) \vec v_{*,g\to p}^{(n)}\biggr)\right)\Delta t }[/math]

      [math]\displaystyle{ \vec v_p^{(n+1)} = \vec v_p^{(n)} + \biggl(\vec a_{g\to p}^{(n)} - \alpha_g \vec v_{g\to p}^{(n)} - \alpha_p \vec v_{p}^{(n)} - \alpha(\beta)\bigl(\vec v_{p}^{(n)} -\vec v_{g\to p}^{(n)} \bigr) - \alpha(\beta) \vec v_{*,g\to p}^{(n)} \biggr)\Delta t }[/math]

where x is position, v is velocity, a is acceleration, and Δt is the time step. Subscript p indicates a particle property and subscript g→p indicates a grid property that has been extrapolated to the particle. Subscripts (n) and (n+1) indicate at the start or end of time step n. The damping terms are:

• α_g - grid damping that scales with the velocity extrapolated to the particle
• α_p - particle damping that scales with the particle velocity
• α(β) - PIC damping where β is the fraction FLIP in the simulation (from 0 to 1) and it scales with the extrapolation error between the particle velocity and the velocity extrapolated from the grid back to the particle.

Although these three damping options could be reduced to two that scale with either particle velocity or the velocity extrapolated from the grid to the particle, separating them as above adds physically insight into the damping mechanisms, especially for the PIC damping option. The [math]\displaystyle{ \vec v_{*,g\to p}^{(n)} }[/math] is a special velocity used to implement an extended for of PIC known as XPIC Damping (and available only in OSParticulas).

Grid and Particle Damping

Two types of grid and particle damping are available. Simple damping just specifies the damping value (which can be a function of time). Alternatively, feedback damping provides a variant of Nose-Hoover feedback, which implements a scheme to evolve the damping coefficient depending on the current total kinetic energy. The net grid or particle damping will be a sum of the two damping options:

      [math]\displaystyle{ \alpha_g = \alpha_{g,simple}(t) + \alpha_{g,feedback}(t) \qquad {\rm and}\qquad \alpha_p = \alpha_{p,simple}(t) + \alpha_{p,feedback}(t) }[/math]

A single simulation can use any combination of these damping options. The α_g,simple and α_p,simple are specified with Damping and PDamping commands, respeectively. The α_g,feedback and α_p,feedback are initialized to zero, but can separately evolve during a simulation based on settings in FeedbackDamping and PFeedbackDamping commands. Feedback damping will evolve such that kinetic energy approaches zero (for a static result) or approaches any non-zero target kinetic energy. The α_g,feedback and α_p,feedback terms evolve by kinetic energy evaluated on the grid:

      [math]\displaystyle{ {d\alpha\over dt} = {2K\over m_{tot}} \left(\sum_{{\rm node}\ i} {1\over 2} m_i |v_i|^2 - T_k \right) }[/math]

where K is a gain parameter, m_tot is total mass, and T_k is the target kinetic energy. The grid and particle feedback damping terms can set separate values for K and T_k. The (2/m_tot) is an arbitrary scaling used to make K similar to a different feedback damping method used in earlier versions of this code and an attempt to make K less dependent on problem size.

PIC Damping

A paper on using MPM for animation suggested that simulations could be improved by modifying the velocity update to include a fraction β of standard FLIP style MPM velocity updates and a fraction (1-β) of PIC style updates. In brief, an undamped FLIP update increments the particle velocity using grid acceleration:

      [math]\displaystyle{ \vec v_{p,FLIP}^{(n+1)} = \vec v_p^{(n)} + \vec a_{g\to p}^{(n)} \Delta t }[/math]

while an undamped PIC update sets particle velocity to the grid velocity extrapolated to the particle:

      [math]\displaystyle{ \vec v_{p,PIC}^{(n+1)} = \vec v_{g\to p}^{(n+1)} }[/math]

The suggestion is that an improved update combines these two updates using:

      [math]\displaystyle{ \vec v_{p}^{(n+1)} = \beta \vec v_{p,FLIP}^{(n+1)} + (1-\beta)\vec v_{p,PIC}^{(n+1)} }[/math]

This FLIP/PIC option appears to provide benefit in some simulations, but its use is difficult to justify. We noticed, however, that the FLIP/PIC update is actually a special case of the above damping where the PIC damping constant is set equal to:

      [math]\displaystyle{ \alpha(\beta) = {1-\beta\over \Delta t} }[/math]

This view leads to two significant improvements for the use of PIC in simulations:

• It provides insight and justification for adding PIC. The PIC damping term is related to the extrapolation error between particle velocity and extrapolated grid velocity. This error will be large in elements that contain particles with widely varying velocity (as sometimes happens in MPM simulations) and small in areas with rather uniform velocity. Thus PIC damping will selectively damp out regions with large velocity variations and have little affect on regions with smooth velocities.
• By implementing PIC as a damping mechanism, it became clear how to modify the position update to include a fraction PIC as well. The original reference on this option did not alter the position update to account for addition of PIC. The above position update addresses that issue and greatly improves particle displacement results in simulations with a significant amount of PIC. Also note that the damping terms in the position update are second order, but because PIC damping is inversely proportional to Δt, the PIC damping terms become first order terms in the position update.

The development of PIC as a damping scheme and how it modifies position update is described in a paper on simulating cutting in MPM

PIC damping is added to simulations using the Damping or PDamping commands. For reference, the undamped FLIP and PIC position updates reduce to

      [math]\displaystyle{ \vec x_{p,FLIP}^{(n+1)} = \vec x_p^{(n)} + \left( \vec v_{g\to p}^{(n+1)} - {\Delta t\over 2}\vec a_{g\to p}^{(n)} \right)\Delta t }[/math]

      [math]\displaystyle{ \vec x_{p,PIC}^{(n+1)} = \vec x_p^{(n)} + \left( \vec v_{g\to p}^{(n+1)} - {1\over 2}(\vec v_p^{(n)} - \vec v_{g\to p}^{(n+1)})\right)\Delta t }[/math]

XPIC Damping

XPIC is an improved form of PIC, which is only available in OSParticulas. The PIC method can be described as applying a projection operator that modifies (and filters) particle velocities before updating them with the grid acceleration. The problem with PIC is that its projection operator filters many problems too heavily resulting in significant dissipation of energy. In contrast, XPIC defines a series of new projection operators that can significantly reduce the over damping of PIC simulations. XPIC is defined by an order m, where m=1 is PIC, m>1 is XPIC, and large m approaches a FLIP method (although may actually approach an improved FLIP method). In other words, m from 1 to infinity provides an interpolation between PIC and FLIP, but is probably a better interpolation than the simple linear combination provided by standard PIC damping.

The drawback of XPIC is that each higher order of XPIC requires two extra extrapolations. XPIC is therefore less efficient than PIC or FLIP (or a PIC damping linear combination of FLIP and PIC). In many problems, m=2 can provide sufficient improvement over PIC and reduce undesirable energy dissipation with minimal extra calculations. Larger m can be tried but typically become very inefficient and have diminishing returns. Very high values of m (e.g., m>40) are typcially unstable (due to too many additional extrapolations).

XPIC is activated by using the XPIC command in combination with activating PIC by having β<1 in a Damping or a PDamping command. Note that XPIC can be used in combination with PIC damping where that method will now provide linear combination of fraction β FLIP and (1-β) of XPIC.

Damping Commands

In scripted files, grid, particle, and PIC damping are activated with the following commands

Damping (alphagVsT),<(fractionPIC)>
FeedbackDamping (gaing),<(Tkg)>,<(maxAlphag)>
PDamping (alphapVsT),<(fractionPIC)>
PFeedbackDamping (gainp),<(Tkp)>,<(maxAlphap)>
XMLData "MPMHeader"
    <XPIC order='(XPICOrder)'/>
EndXMLData

The Damping and FeedbackDamping commands or the PDamping and PFeedbackDamping commands set grid or particle damping terms, respectively. Any combination of these damping options can be used in a simulation. The parameters are:

• (alphagVsT) and (alphapVsT) are constant damping factors or user defined functions of time that evaluate to a damping constant. The units are 1/time units and the default value is 0.
• (gaing) and (gainp) are gains for evolution of α_g or α_p in feedback damping. It has 1/length² units. The default is 0 or no damping.
• (Tkg) and (Tkp) are optional user-defined functions of time that evaluate to a target kinetic energy in μJ. The default value is 0.
• (maxAlphag) and (maxAlphap) can optionally specify a maximum feedback damping coefficient (α_g or α_p) entered in 1/time units. The default is to have no limit to damping, but a maximum can be used to prevent out-of-control overdamping. When a limit is needed, however, it might be a symptom of a poor simulation with too much kinetic energy.
• (fractionPIC) is the fraction PIC to use in the simulation (i.e., 1-β). It can vary from 1 for pure PIC to 0 for pure FLIP. The default is 0. Although this parameter can be set in either the Damping or the PDamping command, only one setting is allowed; a simulation will use whichever setting is last.
• (XPICOrder) to set the XPIC order or m to use for the XPIC method. It must be an integer and values less than 1 are set to 1. Note that this command is ignored unless a Damping or PDamping command activates use of PIC by having (fractionPIC) > 1. The default order is 2 (when the command is present) or 1 when it is not used.

In XML input files, grid, particle, and/or PIC damping are activated with the the following commands, which must be within the <MPMHeader> element:

<Damping PIC='(fractionPIC)' function='(alphagVsT)'>(alphagNum)</Damping>
<FeedbackDamping target='(Tkg)' max='(maxAlphag)'>(gaing)</FeedbackDamping>
<PDamping function='(alphapVsT)'>(alphapNum)</PDamping>
<PFeedbackDamping target='(Tkp)' max='(maxAlphap)'>(gainp)</PFeedbackDamping>
<XPIC order='(XPICOrder)'/>

where

• (alphagNum) and (alphapNum) parameters are constant grid and particle damping factors (only numbers are allowed). These constant values are only used if no function attribute is provided. The units are 1/time units and the default value is 0.

Custom Particle and PIC Damping

The above commands set global damping parameters. If desired, particles for different materials can use custom particle and/or PIC damping. The custom damping for specific materials are set within the material properties when defining any material type.

Using Damping

Adding damping allows a dynamic calculation to approach static or quasi static results after sufficiently long times. Too much damping, however, will alter the results. The challenge is to select the appropriate amount of damping. The amount of damping to approach static results without over damping depends on the problem being damped. Under static conditions (no longer applying forces or moving particles), damping will cause the total kinetic energy on the grid to decay with a exponential time constant on the order of

      [math]\displaystyle{ \tau = {2\over \alpha} }[/math]

One set of tension simulations for an isotropic material with E = 2300 MPa and ρ = 1 g/cm³ were damped weill with α = 5000 sec^-1 or τ ≈ 0.4 msec. For other materials, an appropriate amount of damping will likely scale from this choice approximately with wave speed √(E/ρ) of the materials. When the problem involves bending, however, the damping coefficient may need to based on the transverse vibrational frequency rather than on material wave speed.

For problems with monotonic loading, you would not want to damp to zero kinetic energy, but rather to a non-zero kinetic energy appropriate for the problem. For example, a simple tension test the x direction on a bar with one end (at x=0) fixed and the other end (at x=l) moving at velocity v in quasi-static conditions, would have position-dependent velocity of vx/l. The quasi-static kinetic energy for this simulation would be

      [math]\displaystyle{ T_k = \int_V {\rho\over2} \left({vx\over l}\right)^2 dV = {1\over 6} m v^2 }[/math]

where m is total mass of the object. This problem is best damped using feedback damping with the target kinetic energy equal to this evaluated T_k. The (gain) parameter would have to be selected by trial and error to avoid too little or too much damping.

You can monitor total grid and particle damping constants, α_p and α_p, by using global archiving options. Ideally the simulation α will remain as a reasonable number for the dynamics of the problem (i.e., provided damping gain is appropriate for the problem). When using a target kinetic energy, α should decrease as kinetic energy moves toward the target result. If the damping cannot keep up with dynamic effects, the α may continue to increase and may eventually cause overdamping. The (maxAlpha) parameter can be used to prevent α growing too large, but typically a need to use this parameter might mean the simulations will not work well and other methods should be considered.

The differences between grid and particle damping are typically minor. The two options are available, however, in case one type (or a combination) works better in specific problems.

PIC damping can be used alone or in addition to grid and particle. It seems to have the most affect in problems involving non-zero velocity boundary conditions, problems involving contact, and problems with crack propagation. It has less affect on problems with traction boundary conditions. The current recommendation is purely trial and error. Try adding some fraction PIC and see if it helps the results without altering them to a invalid solution.

Notes

1. The use of feedback damping in MPM calculations was first described in Ayton, et al., 2002. This reference evolved α using
         [math]\displaystyle{ {d\alpha\over dt} = K \sum_p |v_p|^2 }[/math]
   In other words the update was based on sum of particle velocities and there was no target kinetic energy. This approach is only suitable for problems with static boundary conditions. The modifications in NairnMPM allow the use of damping with variable boundary conditions (assuming you can estimate the expected kinetic energy) and uses kinetic energy rather than a velocity sum. The use of kinetic energy makes the method more general and therefore better able to handle composite materials having particles of different masses. Furthermore, the kinetic energy for NairnMPM damping is found from grid masses and velocities rather than particle masses and velocities becuase in some problems particle kinetric energy can develop artifacts while grid kinetic energy is more robust estimate of total kinetic energy.
2. A Nose-Hoover thermostat works in molecular calculations by tying kinetic energy to temperature. In MPM the thermostat ties kinetic energy to an expected kinetic energy after vibrations are damped out. When boundary conditions are constant, either simple damping or feedback damping with no target kinetic energy will help the problem converge to the static solution. For non-constant boundary conditions, you need to use feedback damping, estimate the "static" (i.e., very slow loading) kinetic energy, and set (Tk) parameter to this estimate to have this thermostat work correctly.