Metafor

ULiege - Aerospace & Mechanical Engineering

User Tools

Site Tools


doc:user:integration:scheme:dynimpl

This is an old revision of the document!


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

Generalized Midpoint Rule

The equilibrium i 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)
<\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 ====

=== Central Difference ===

<code>
ti = CentralDifferenceTimeIntegration(metafor)
metafor.setTimeIntegration(ti)

Chung Hulbert

ti = ChExplicitTimeIntegration(metafor)
ti.setRhoB(_rhoB)
metafor.setTimeIntegration(ti)

The parameter “_rhoB” is the spectral radius at bifurcation point ([0, 1]). The default value is 0.8182.

Tchamwa

ti = TchamwaTimeIntegration(metafor)
ti.setRhoB(_rhoB)
metafor.setTimeIntegration(ti)

The parameter “_rhoB” is the spectral radius at bifurcation point ([0, 1]). The default value is 0.8182.

Useful parameters : see quasi-statique.

doc/user/integration/scheme/dynimpl.1444402373.txt.gz · Last modified: 2016/03/30 15:22 (external edit)

Donate Powered by PHP Valid HTML5 Valid CSS Driven by DokuWiki