# Metafor

ULiege - Aerospace & Mechanical Engineering

### Site Tools

doc:user:integration:scheme:dynimpl

# 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 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]).

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 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:

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.

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

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

### 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 :

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

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

#### New Metafor Version > 2422

##### Consistent Time Integration (EMCA or EDMC)
ti = ConsistentTimeIntegration(metafor)
metafor.setTimeIntegration(ti)
##### Generalized Midpoint Rule
ti = MpgTimeIntegration(metafor)
ti.setTheta(_theta)
metafor.setTimeIntegration(ti)

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

##### Alpha Generalized Family
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.