Implicit dynamic integration schemes

This method can manage several sub-methods :Newmark, HHT, Chung Hulbert, generalized midpoint rule, conservative scheme.

Concerning Chung and Hulbert scheme, the equilibrium equation is rewritten pondering inertia forces by the parameter $\alpha_M$ (MDR_ALPM), and internal and external forces by the parameter $\alpha_F$(MDR_ALPF) (see Global Parameters [REMOVED]). For time step $t^{n+1}$, this leads to :

$$(1-\alpha_M) \boldsymbol{F}^{\text{inert}}(t^{n+1}) + \alpha_M \boldsymbol{F}^{\text{inert}}(t^n) + (1-\alpha_F) \boldsymbol{F}^{\text{int}}(t^{n+1}) + \alpha_F \boldsymbol{F}^{\text{int}}(t^n)$$ $$ = (1-\alpha_F) \boldsymbol{F}^{\text{ext}}(t^{n+1}) + \alpha_F \boldsymbol{F}^{\text{ext}}(t^n)$$

In addition, relations between the displacements $\boldsymbol{x}$, velocities $\boldsymbol{v}$ and accelerations $\boldsymbol{a}$ are:

$$\boldsymbol{x}(t^{n+1}) = \boldsymbol{x}(t^n) + [t^{n+1}-t^n] \boldsymbol{v}(t^n) + (t^{n+1}-t^n)^2 \{[0.5-\beta]\boldsymbol{a}(t^n)+\beta \boldsymbol{a}(t^{n+1})\} $$

$$\boldsymbol{v}(t^{n+1}) = \boldsymbol{v}(t^n) + [t^{n+1}-t^n] \{[1-\gamma]\boldsymbol{a}(t^n)+\gamma \boldsymbol{a}(t^{n+1})\} $$

Using particular values of these weighting parameters can lead to :

$\alpha_M = \alpha_F = 0 $: Newmark scheme
$\alpha_M = 0 $: Hilber-Hughes-Taylor scheme (HHT)
$\alpha_F >0, \alpha_M \neq 0 $: Chung-Hulbert scheme (CH)
$\alpha_F = 0 $: Wood-Bossak-Zienkiewicz scheme (WBZ)

Unconditional stability (stability for any time step) is guaranteed for

$$\alpha_M \leq \alpha_F \leq \frac{1}{2} $$

$$\gamma \geq \frac{1}{2} - \alpha_M + \alpha_F $$

$$\beta \geq \frac{1}{4} \left(\frac{1}{2} + \gamma\right)^2 $$

Second order accuracy is obtained for:

$$\gamma = \frac{1}{2} - \alpha_M + \alpha_F $$

Optimal numeric dissipation of high frequencies is obtained for:

$$3 \alpha_F - \alpha_M = 1 $$

$$\beta = \frac{1}{4} \left(1 - \alpha_M + \alpha_F\right)^2 $$

The only combination leading to a second order accuracy, unconditional stability and no numerical dissipation is:

$$\alpha_M = \alpha_F = 0 $$

$$\gamma = 0.5 $$

$$\beta = 0.25 $$

Newmark dissipative schemes ($\alpha_M = \alpha_F = 0$, $\gamma >\frac{1}{2}$, $\beta = \frac{1}{4} (\frac{1}{2} + \gamma)^2$) is only accurate to the first order. Integration parameters can then be computed as a function of the spectral radius for an infinite frequency ($\rho_{\infty}$). The parameter $\beta$ is entered by MDR_BET0 and $\gamma$ by MDR_GAM0 (see Global Parameters [REMOVED]).

Second order dissipative schemes (HHT, CH, WBZ) are activated by the parameters $\alpha_M$ (MDR_ALPM), $\alpha_F$ (MDR_ALPF), $\beta_0$ (MDR_BET0) and $\gamma_0$ (MDR_GAM0) (see Global Parameters [REMOVED]). Parameters $\beta$ and $\gamma$ are computed internally :

• $\gamma = \gamma_0 (1 - 2 \alpha_M + 2 \alpha_F) $
• $\beta = \beta_0 (1 - \alpha_M + \alpha_F)^2 $

When these parameters are written as functions of the spectral radius for infinite frequency:

Adding damping forces to the Alpha-Generalized time integration scheme family allow user to dissipate energy in a regulated way. According to the reference below, damping forces can be balanced between previous and current time as others (internal or external) forces through $\alpha_F$ parameter. The global system is then written by :

$$(1-\alpha_M) \boldsymbol{F}^{\text{inert}}(t^{n+1}) + \alpha_M \boldsymbol{F}^{\text{inert}}(t^n) + (1-\alpha_F) \boldsymbol{F}^{\text{damp}}(t^{n+1}) + \alpha_F \boldsymbol{F}^{\text{damp}}(t^n)$$ $$+ (1-\alpha_F) \boldsymbol{F}^{\text{int}}(t^{n+1}) + \alpha_F \boldsymbol{F}^{\text{int}}(t^n) = (1-\alpha_F) \boldsymbol{F}^{\text{ext}}(t^{n+1}) + \alpha_F \boldsymbol{F}^{\text{ext}}(t^n)$$

where the damping forces are computed proportional to velocities : $$ \boldsymbol{F}^{\text{damp}} = \boldsymbol{C} * v $$

According to ref2 (chapter 3), an easy way to build a diagonal damping matrix is to compute a ponderated sum of the mass and stiffness matrices : $$\boldsymbol{C} = a_k \boldsymbol{K} + a_m \boldsymbol{M}$$

The modal damping factor corresponding to each eigen pulsation $\omega_{0r} = 2 \pi \phi$: $$\epsilon_r=\frac{1}{2}(a_k \omega_{0r} + \frac{a_m}{\omega_{0r}} )$$

telling us that the mass damping factor $a_m$ will induice damping on lowest eigen frequencies and that the stiff damping factor $a_k$ will induce damping on the higher eigen frequencies.

ref 1 : "The analysis of the Generalized-$\alpha$ method for non linear dynamic problems" -
        S.Erlicher, L.Bonaventura, O.S.Bursi -  Computanional Mechanics 28 (2002) 83-104
ref 2 : "Théorie des vibrations - Application à la dynamique des structures" -
         M.Géradin, D.Rixen - Editions Masson

• Name of the scheme : DampedAlphaGeneralizedTimeIntegration
• Parameters :
  • AlphaGeneralizedTimeIntegration parameters has to be defined as usual
  • update of the damping matrix managed through ti.setDampingMatrixUpdate(DMUxxx)
    • DMUINIT : Damping matrix computed on initial configuration only
    • DMUPERSTAGE : Damping matrix computed at each stage change
    • DMUPERSTEP : Damping matrix computed at the beginning of each time step (base on previous equilibrated solution to avoid need of stiffness computation)
  • Mass Damping Factor and Stiffness Damping Factor are defined though the element properties (DAMPSTIFF,DAMPMASS). Parameters can depend of time.

The equilibrium is computed for a time defined by the parameter $\theta$ (MDR_THEM) (see Global Parameters [REMOVED]):

$$\boldsymbol{F}^{\text{inert}}(t^{n+\theta}) + \boldsymbol{F}^{\text{int}}(t^{n+\theta})= \boldsymbol{F}^{\text{ext}}(t^{n+\theta})$$

In addition, relations between the displacements $\boldsymbol{x}$, velocities $\boldsymbol{v}$ and accelerations $\boldsymbol{a}$ are:

$$\begin{eqnarray} \boldsymbol{x}(t^{n+\theta}) &=& \boldsymbol{x}(t^n) + \theta(t^{n+1} - t^n) \boldsymbol{v}(t^n) + \frac{1}{2} \theta^2(t^{n+1} - t^n)^2 \boldsymbol{a}(t^{n+\theta}) \\ \boldsymbol{v}(t^{n+\theta}) &=& \boldsymbol{v}(t^n) + \theta(t^{n+1} - t^n) \boldsymbol{a}(t^{n+\theta}) \end{eqnarray}$$ once the equilibrium is reached, values for $t^{n+1}$ are determined:

$$\begin{eqnarray} \boldsymbol{x}(t^{n+1}) &= &\boldsymbol{x}(t^n) + (t^{n+1} - t^n) \boldsymbol{v}(t^n) + \frac{1}{2} (t^{n+1}-t^n)^2 \boldsymbol{a}(t^{n+\theta}) \\ \boldsymbol{v}(t^{n+1}) &=& \boldsymbol{v}(t^n) + (t^{n+1} - t^n) \boldsymbol{a}(t^{n+\theta}) \\ \boldsymbol{a}(t^{n+1}) &=& \boldsymbol{a}(t^{n+\theta}) \end{eqnarray}$$ This scheme is unconditionally stable and dissipative for : $ \theta>1 $

This is a midpoint where internal forces $\boldsymbol{F}^{\text{int}}$ are caclualted to conserve the energy. For a hypoelastic material with plasticity, an approximation of the forces is given for CONSERVINGMETHOD=VES_WITHOUTCORRECTION, when the exact formulation requires CONSERVINGMETHOD=VES_WITHCORRECTION. Numerical dissipation forces $\boldsymbol{F}^{\text{diss}}$ and numerical dissipation velocities $G^{\text{diss}}$ can be induced for a given dissipation parameter. The equilibrium is computed :

$$ \boldsymbol{F}^{\text{inert}}(t^{n+1/2}) + \boldsymbol{F}^{\text{int}} + \boldsymbol{F}^{\text{diss}} = \boldsymbol{F}^{\text{ext}}(t^{n+1/2}) $$

Relations between displacements $\boldsymbol{x}$, velocities $\boldsymbol{v}$ and accelerations $\boldsymbol{a}$ are:

$$\begin{eqnarray} \boldsymbol{x}(t^{n+1}) &=& \boldsymbol{x}(t^n) + 0.5 (t^{n+1}-t^n) \boldsymbol{v}(t^n) + 0.5 (t^{n+1}-t^n) \boldsymbol{v}(t^{n+1}) + (t^{n+1}-t^n) G^{\text{diss}} \\ \boldsymbol{v}(t^{n+1}) &=& \boldsymbol{v}(t^n) + 0.5 (t^{n+1}-t^n) \boldsymbol{a}(t^n) + 0.5 (t^{n+1}-t^n) \boldsymbol{a}(t^{n+1}) \\ \boldsymbol{a}(t^{n+1/2}) &=& 0.5 \boldsymbol{a}(t^n) + 0.5 \boldsymbol{a}(t^{n+1}) \end{eqnarray}$$

This scheme is unconditionally stable but requires programming an adequate formulation of internal forces for each element and material. When the scheme is done without numerical dissipation (EMCA), the accuracy is of the second order, but when there is numerical dissipation (EDMC) ($\xi=$ dissipation), accuracy is only of the first order. The dissipation parameter $\xi$ can be written as function of the spectral radius for infinite frequency (ideally $\rho_\infty=0.8$)

It is mandatory to instanciate the ConsistentAlgorithmFunctions and then the parameters can adjusted if required :

    caf = ConsistentAlgorithmFunctions(metafor)
    caf.setDissipationParameter(diss)
    caf.setConservingContactMethod(consContactMeth)

where

Contact methods can be combined with the following methods:

Contact methods can be combined with the following materials:

For the domain domain, material number no and mass mass

matset[no].put(MASS_DENSITY, mass)

metafor.getInitialConditionSet().define(numero, entite, lock, ampl)

(see Global Parameters [REMOVED])

ti = ConsistentTimeIntegration(metafor)
metafor.setTimeIntegration(ti)

ti = MpgTimeIntegration(metafor)
ti.setTheta(_theta)
metafor.setTimeIntegration(ti)

The parameter “_theta” is defined above. The default value is 1.1.

ti = AlphaGeneralizedTimeIntegration(metafor)
ti.setAlphaM(_AlphaM)
ti.setAlphaF(_AlphaF)
ti.setBeta0(_Beta0)
ti.setGamma0(_Gamma0)
metafor.setTimeIntegration(ti)

The parameters “_AlphaM”, “_AlphaF”, “_Beta0” and “_Gamma0” are defined above. The default values are -0.97, 0.01, 0.25 and 0.5 respectively.

Useful parameters : see quasi-statique.