Damping Options

NairnMPM has two forms of grid damping. Their most common use to damp dynamic effects such that the solution converges to a static result.

Introduction

Grid Damping

Grid damping is a damping term applied to each node on the grid that alters the grid forces (or accelerations). In brief, the force on node i for step n with damping is modified to be:

      [math]\displaystyle{ \vec f_{i}^{(n)*} =\vec f_{i}^{(n)} - \alpha \vec p_i^{(n)} }[/math]

where f_i is total force on node i from all sources, p_i is nodal momentum, and α is the grid damping constant (in units of 1/sec). The nodal velocity update becomes:

      [math]\displaystyle{ \vec v_{i}(t+\Delta t) = \vec v_{i}(t) + {\vec f_{i}^{(n)*}\over m_i^{(n)}}\Delta t = \vec v_{i}(t) + (\vec a_i - \alpha \vec v_i) \Delta t }[/math]

where a_i = f_i⁽ⁿ⁾/m_i⁽ⁿ⁾ is nodal acceleration without damping, v_i = p_i⁽ⁿ⁾/m_i⁽ⁿ⁾ is nodal velocity, and m_i⁽ⁿ⁾ is nodal mass. Similarly, the particle velocity update (when done using grid accelerations) is

      [math]\displaystyle{ \vec v_{p}(t+\Delta t) = \vec v_{p}(t) + \vec a_{g\to p}^{(n)*}\Delta t = \vec v_{p}(t) + (\vec a_{g\to p}^{(n)} - \alpha\vec v_{g\to p}^{(n)})\Delta t }[/math]

where g→p indicates acceleration or velocity extrapolated from the grid to the particle, superscript * means extrapolation using force modified for damping, and no superscript means extrapolation using unmodified force.

Two types of grid damping are available in NairnMPM. 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 α will be a sum of the two damping options:

      [math]\displaystyle{ \alpha = \alpha_{simple}(t) + \alpha_{feedback}(t) }[/math]

A single simulation can use one or the other type or combine them. The α_simple is specified with a Damping command. The α_feedback is initialized to zero, but can evolve during a simulation based on settings in a FeedbackDamping command. Feedback damping will evolve such that kinetic energy approaches zero (for a static result) or approaches any non-zero target kinetic energy. The α_feedback term evolves 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 (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.

Particle Damping

Particle damping is applied during the particle's velocity update and the damping term is related to the initial particle velocity instead of the grid velocity (which is used by grid damping). It is currently only available in OSPaticulas. With particle damping, the particle velocity update becomes

      [math]\displaystyle{ \vec v_{p}(t+\Delta t) = \vec v_{p}(t) + (\vec a_{g\to p}^{(n)*}-\alpha_p\vec v_{p}(t))\Delta t = \vec v_{p}(t) + (\vec a_{g\to p}^{(n)} - \alpha\vec v_{g\to p}^{(n)}-\alpha_p\vec v_{p}(t))\Delta t }[/math]

where α_p is a particle damping constant (in units of 1/sec). Notice that this update shows both grid and particle damping, which means they be combined or used separately. As with grid damping, there are two types of particle damping methods resulting in net particle damping constant of:

      [math]\displaystyle{ \alpha_p = \alpha_{p,simple}(t) + \alpha_{p,feedback}(t) }[/math]

The α_p,simple is specified with a PDamping command. The α_p,feedback is initialized to zero, but can evolve during a simulation based on settings in a PFeedbackDamping command. The particle α_p,feedback constant evolves exactly like the grid damping α_feedback, except it can specify different values for the gain parameter (K_p) and target kinetic energy (T_p,k)

Damping Commands

In scripted files, grid damping is activated with the two commands

Damping (alphaVsT)
FeedbackDamping (gain),<(Tk)>,<(maxAlpha)>

Although both can be used, usually only one damping method is used in a simulation. The parameters are:

• (alphaVsT) is a constant damping factor or a user defined function of time that evaluates to a grid damping constant. The units are 1/sec and the default value is 0.
• (gain) is gain for evolution of α in feedback damping. It has units of 1/mm². The default is 0 or no damping.
• (Tk) in an optional user-defined function of time that evaluates to a target kinetic energy in μJ. The default value is 0.
• (maxAlpha) can optionally specify a maximum feedback damping coefficient (α) entered in units of 1/sec. 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.

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

<Damping function='(alphaVsT)'>(alphaNum)</Damping>
<FeedbackDamping target='(Tk)' max='(maxAlpha)'>(gain)</FeedbackDamping>
<PDamping function='(alphaVsT)'>(alphaNum)</PDamping>
<PFeedbackDamping target='(Tk)' max='(maxAlpha)'>(gain)</PFeedbackDamping>

where the (alphaNum) parameters are constant grid and particle damping factors (only numbers are allowed). These constant value are only used if no function attribute is provided. The units are 1/sec and the default value is 0. The Damping and FeedbackDamping commands set grid damping and particle damping terms, respectively.

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 = {1\over 2\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.1 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 vibration 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, α 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 time constant 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 expected to be minor. The two options are available, however, in case one type (or a combination) works better in specific problems.

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.