Table of Contents
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.