====== Implicit dynamic integration schemes ======

===== Description =====

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


==== Alpha-Generalized Family ====
 
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 [[doc:user:integration:general:parameters]]). 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 [[doc:user:integration:general:parameters]]).

^     Algorithm     ^  $\beta$  ^  $\gamma$  ^  Order  ^
|conservative Newmark |  $\frac{1}{4}$   |   $\frac{1}{2}$  |    2    |
|dissipative Newmark  |  $\displaystyle\frac{1}{(1+\rho_{\infty})^2}$ |  $\displaystyle\frac{3-\rho_{\infty}}{2+2\rho_{\infty}}$ |  1  | 

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 [[doc:user:integration:general:parameters]]). 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:

^     Algorithme     ^  $\beta$  ^  $\beta_0$  ^  $\gamma$  ^  $\gamma_0$  ^  $\alpha_M$  ^  $\alpha_F$  ^ Ordre  ^
| Wood - Bossak - Zienkiewicz |  $\beta_0(1-\alpha_M)^2$  |  0.25  |  $\gamma_0 (1-2\alpha_M)$  |  0.5  |  $\frac{\rho_\infty-1}{\rho_\infty+1}$  |  0  |  2  |
| Hilber-Hughes-Taylor |  $\beta_0(1+\alpha_F)^2$  |  $\frac{1}{4}$  |  $\gamma_0 (1+2\alpha_F)$  |  $\frac{1}{2}$  |  0  |  $\displaystyle\frac{1-\rho_\infty}{1+\rho_\infty}$  |  2  | 
| Chung-Hulbert |  $\beta_0(1-\alpha_M+\alpha_F)^2$  |  $\frac{1}{4}$  |  $\gamma_0 (1-2\alpha_M+2\alpha_F)$  |  $\frac{1}{2}$  |  $\displaystyle\frac{2\rho_\infty-1}{\rho_\infty+1}$  |  $\displaystyle\frac{\rho_\infty}{\rho_\infty+1}$  |  2  |
 

==== Damped Alpha-Generalized family ==== 
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.

<code>
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
</code>
=== Parameters to the scheme === 
  * 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 [[doc:user:elements:volumes:volumeelement?&#parameters|element properties]] (''DAMPSTIFF'',''DAMPMASS''). Parameters can depend of time.
        
==== Generalized Midpoint Rule ==== 

The equilibrium is computed for a time defined by the parameter $\theta$ (''MDR_THEM'') (see [[doc:user:integration:general:parameters]]):

$$\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 $


==== Conservative scheme ====

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$)

^ Algorithm   ^                       $\xi$                          ^  order ^
| EMCA         |                         0                            |    2    |
| EDMC         |  $\displaystyle\frac{1-\rho_\infty}{1+\rho_\infty}$  |    1    |

=== Management of conservative algorithm ===

It is **mandatory** to instanciate the ConsistentAlgorithmFunctions and then the parameters can adjusted if required :
<code>
    caf = ConsistentAlgorithmFunctions(metafor)
    caf.setDissipationParameter(diss)
    caf.setConservingContactMethod(consContactMeth)
</code>
where

|''consContactMeth ''| = CONSERVINGCONTACTMETHOD_ARMEROPETOCZ : the default contact method is the one from Armero and Petöcz. |
|        | = CONSERVINGCONTACTMETHOD_LAURSEN : the contact method is the one from Laursen. (not validated) |
| 	 | = CONSERVINGCONTACTMETHOD_LOVE : the contact method is the one from Love (discontinuities in velocities) . (Not validated)|
|''diss''| numerical dissipation (default value 0.)|

Contact methods can be combined with the following methods:

^     Contact method        ^  Penalty ^ Augmented Lagrangian ^
| Classical                |   X       |       X             |
| Consistent Armero-Petöcz |   X       |       0             |
| Consistent Laursen       |   X       |       X             | 
| Consistent Love          |   X       |       X             | 

Contact methods can be combined with the following materials:

^     Contact method    ^  Frictionless ^ VariablePenaltyFrictionless ^ Sticking ^ VariablePenaltySticking ^ Coulomb ^ Tresca ^
| Classical                 |   X  |  X   |  X  |   X     |  X       |  X  |
| Consistent Armero-Petöcz |   X  |  0   |  X  |   0     |  X       |  X  |
| Consistent Laursen       |   X  |  0   |  X  |   X     |  X       |  X  |         
| Consistent Love          |   X  |  0   |  0  |   0     |  0       |  0  |  

===== Data set =====

==== Defining material density ====

For the domain ''domain'', material number ''no'' and mass ''mass''

  matset[no].put(MASS_DENSITY, mass)

==== Defining initial velocites ====

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

==== Choosing the algorithm ====

=== Old Metafor Version <= 2422 ===

^          Scheme         ^  ''MDE_NDYN''  ^  ''MDR_ALPM''  ^  ''MDR_ALPF''  ^  ''MDR_BET0''  ^  ''MDR_GAM0''  ^  ''MDR_THEM''  ^   ''CONSERVINGMETHOD''  ^  Dissipation  ^
|Chung Hulbert            |       2        |        X       |       X        |       X        |       X        |                |                         |               |
|HHT                      |       2        |        0       |       X        |       X        |       X        |                |                         |               |
|Newmark                  |       2        |        0       |       0        |       X        |       X        |                |                         |               |      
|Generalized midpoint     |       4        |                |                |                |                |       X        |                         |               |
|Conservative              |       5        |                |                |                |                |                |            X            |        X      | 

(see [[doc:user:integration:general:parameters]])

=== New Metafor Version > 2422 ===

== Consistent Time Integration (EMCA or EDMC) ==

<code>
ti = ConsistentTimeIntegration(metafor)
metafor.setTimeIntegration(ti)
</code>

== Generalized Midpoint Rule ==

<code>
ti = MpgTimeIntegration(metafor)
ti.setTheta(_theta)
metafor.setTimeIntegration(ti)
</code>

The parameter "_theta" is defined above. The default value is 1.1.

== Alpha Generalized Family == 

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

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 [[quasistatique|quasi-statique]].