Newmark algorithm for solving transient problems

For the sake of brevity in the rest of this page, we denote by $\ \underline{\delta}$ the vector of nodal unknowns, mixing nodal displacements and temperature difference in a known order. The Finite Element matrices are also assumed to follow the same numbering.

We assume the transient study starts at time $\ t = 0$ and goes on with time step $\Delta t$ for $N$ time steps. Subscript $\ \underline{\delta}_n$ will denote the time step numbering, with $\ n \in {1, \ldots, N}$.

As, in this package, mechanical and thermal properties are assumed constant through time, only the vector of nodal unknowns $\ \underline{\delta}$ is indexed by time.

Also, the global force vector $\ \underline{F}$ is assumed to already incorporate Dirichlet boundary conditions with respect to displacement and temperature. This means that nodal unknowns vector $\ \underline{\delta}$ only contains unconstrained degrees of freedom.

Finally, as this algorithm is sometimes called the Newmark-beta method, parametrized by $\ \gamma$ and $\ \beta$, we will use both these symbols in the equation, meaning that $\ \beta$ does not refer to the thermal modulus in this page. Common values for these parameters are $\ \gamma = \frac{1}{2}$ and $\ \beta = \frac{1}{4}$, which yields an implicit scheme called the average constant acceleration.

Updating rules of the Newmark algorithm

Given acceleration $\ \underline{\delta}_{n+1}$ at time step $\ n+1$, we have the following updating rules for the velocity and 'displacement' (which is a mix of mechanical displacement and temperature differences):

$$\underline{\delta}_{n+1} = \underline{\dot{\delta}}_n + (1-\gamma) \; \Delta t \; \underline{\ddot{\delta}}_n + \gamma \; \Delta t \; \underline{\ddot{u}}_{n+1}$$ $$\underline{\delta}_{n+1} = \underline{\delta}_n + \Delta t \; \underline{\dot{u}}_n + \frac{\Delta t^2}{2} \; \left((1 - 2\beta) \; \underline{\ddot{u}}_n + 2\beta \; \underline{\ddot{u}}_{n+1}\right)$$

Governing equation

As per this page, the discretized equation for classical linear thermoelasticity can be formulated in terms of nodal unknowns at time step $\ n+1$ as:

$$[M] \; \underline{\ddot{\delta}}_{n+1} + [D] \; \underline{\dot{\delta}}_{n+1} + [K] \; \underline{\delta}_{n+1} = \underline{F}$$

Calculation of acceleration

Injecting the updating rules in the governing equation gives us a system of equations where only $\ \underline{\ddot{u}}_{n+1}$ is unknown, the rest of the terms depending on the previous time steps, where the acceleration, velocity and 'displacement' are known. By reorganizing with terms related to $\ \underline{\delta}_{n+1}$ on the left-hand side and the rest on the right-hand side, we obtain:

$$\left([M] + \gamma \; \Delta t \; [D] + \beta \; \Delta t^2 \; [K]\right) \; \underline{\delta}_{n+1} = \underline{F} - [D] \; \left(\underline{\dot{\delta}}_n + (1 - \gamma) \; \Delta t \; \underline{\ddot{\delta}}_n\right) - [K] \; \left(\underline{\delta}_n + \Delta t \; \underline{\dot{\delta}}_n + \Delta t^2 \; \frac{1 - 2\beta}{2} \; \underline{\ddot{\delta}}_n\right)$$

For brevity, we rewrite this as:

$$[K^{\textrm{dyn}}] \; \underline{\ddot{\delta}}_{n+1} = \underline{F}^{\textrm{dyn}}$$

which yields $\underline{\ddot{\delta}}_{n+1}$ by inverting this system where $[K^{\textrm{dyn}}]$ is called the dynamic stiffness matrix and $\underline{F}^{\textrm{dyn}}$ is called the dynamic force vector as such:

$$\underline{\ddot{\delta}}_{n+1} = [K^{\textrm{dyn}}]^{-1} . \underline{F}^{\textrm{dyn}}$$

Algorithm step-by-step:

1. Initialize $\ \underline{\delta}_{1}$ and $\ \underline{\dot{\delta}}_{1}$ with known initial vectors, as well as $\ \underline{\ddot{\delta}}_{1}$ if it is known. Set $\ n=1$
2. Calculate $\ \underline{\ddot{\delta}}_{n+1}$
3. Use the updating rules to obtain $\ \underline{\dot{\delta}}_{n+1}$ and $\ \underline{\delta}_{n+1}$
4. Set $\ n \leftarrow n+1$ and go to step 2 until step $\ N$ has been performed