A new model-dependent time integration scheme with effective numerical damping for dynamic analysis

Abstract

This paper presents a new semi-explicit dissipative model-dependent time integration algorithm for solving structural dynamics problems. Motivated by the superior properties of the composite time-stepping scheme, the proposed method is designed, so that it fully inherits the numerical characteristics of its parent algorithm, namely the Bathe method. The algorithm design procedure is carried out by assuming unknown integration parameters for the proposed method. Afterwards, by time discretization of an SDOF model equation, the unknown parameters can be obtained explicitly by solving nonlinear system of equations. Some numerical examples are analyzed by the presented technique and comparisons are also made with two other dissipative model-dependent time integration algorithms as well as the Bathe method. Results demonstrate that the suggested technique can effectively damp out the spurious oscillations of the high-frequency modes, while the other schemes exhibit significant overshoot in the calculated responses. Furthermore, it is also observed that numerical results of the presented method totally coincide with the parent algorithm. While the Bathe method subdivides each time increment into two sub-steps, the proposed algorithm is single-step, non-iterative and does not involve any time-step subdividing.

Appendices

Appendix 1

$$\begin{aligned} g_{0} & = \left( {\gamma - 2} \right)^{2} , \, g_{1} = 2\left( {\gamma - 2} \right)\left( {\gamma^{2} - 2} \right), \, g_{2} = 2\gamma \left( {\gamma - 2} \right)\left( {\gamma - 1} \right) \\ g_{3} & = \gamma^{2} - 2\gamma + 2, \, g_{4} = \gamma \left( {\gamma - 1} \right)\left( {\gamma^{2} - 2} \right), \, g_{5} = \gamma \left( {\gamma - 1} \right),g_{6} = \left( {\gamma - 2} \right)\left( {\gamma - 1} \right) \\ g_{7} & = \left( {\gamma - 1} \right)^{2} , \, g_{8} = 4 - \gamma^{2} \left( {6 + \gamma \left( {\gamma - 4} \right)} \right),g_{9} = \gamma \left( {1 - \gamma } \right)\left( {\gamma^{2} - 2\gamma + 2} \right) \\ g_{10} & = 2\left( {\gamma - 2} \right)\left( {\gamma^{2} - 2\gamma + 2} \right), \, g_{11} = \gamma \left( {4 + \gamma } \right) - 6,g_{12} = \gamma^{2} - 2. \\ \end{aligned}$$

(45)

Appendix 2

Table 2 summarizes the necessary steps for conducting time integration through the MDED technique. In this table, N represents the numerator of the integration parameters in MDOF space. Moreover, $\left[ {\mathbf{I}} \right]$ is the identity matrix of size q, where q is the number of degrees of freedom.

Table 2 Pseudocode for the proposed MDED method

Appendix 3

The ‘N’ notation is an operator which returns the numerator of a rational function. Similarly, operator ‘D’ could be defined which give the denominator of a rational function. For instance, considering α_i (integration parameter) as a rational function:

$$\alpha_{i} = \frac{1 + 2\zeta \varOmega }{{1 + 2\zeta \varOmega + \varOmega^{2} }} \to \left\{ \begin{aligned} N_{{\left( {\alpha_{i} } \right)}} = 1 + 2\zeta \varOmega \hfill \\ D_{{\left( {\alpha_{i} } \right)}} = 1 + 2\zeta \varOmega + \varOmega^{2} \hfill \\ \end{aligned} \right..$$

(46)

Thus

$$\alpha_{i} = \frac{{N_{{\left( {\alpha_{i} } \right)}} }}{{D_{{\left( {\alpha_{i} } \right)}} }} = \left( {D_{{\left( {\alpha_{i} } \right)}} } \right)^{ - 1} N_{{\left( {\alpha_{i} } \right)}} .$$

(47)

Regarding the generalization of this concept to MDOF systems (Matrix form), a simple example is introduced to clarify the matter. Let us suppose one wants to generalize α_i for a system with three degrees of freedom. Since every mode of vibration has a unique natural frequency, there are three integration parameters (with different $\varOmega$ and ζ). Thus, one could write

$$\alpha_{i1} = \frac{{1 + 2\zeta_{1} \varOmega_{1} }}{{1 + 2\zeta_{1} \varOmega_{1} + \varOmega_{1}^{2} }},\quad \alpha_{i2} = \frac{{1 + 2\zeta_{2} \varOmega_{2} }}{{1 + 2\zeta_{2} \varOmega_{2} + \varOmega_{2}^{2} }},\quad \alpha_{i3} = \frac{{1 + 2\zeta_{3} \varOmega_{3} }}{{1 + 2\zeta_{3} \varOmega_{3} + \varOmega_{3}^{2} }}.$$

(48)

Here, subscripts 1, 2 and 3 represent the mode number. Equation (48) could be written, alternatively, in matrix form:

$$\left[ {\bar{\rm A}_{i} } \right] = \left[ {\begin{array}{*{20}c} {\alpha_{i1} } & 0 & 0 \\ 0 & {\alpha_{i2} } & 0 \\ 0 & 0 & {\alpha_{i3} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {\frac{{1 + 2\zeta_{1} \varOmega_{1} }}{{1 + 2\zeta_{1} \varOmega_{1} + \varOmega_{1}^{2} }}} & 0 & 0 \\ 0 & {\frac{{1 + 2\zeta_{2} \varOmega_{2} }}{{1 + 2\zeta_{2} \varOmega_{2} + \varOmega_{2}^{2} }}} & 0 \\ 0 & 0 & {\frac{{1 + 2\zeta_{3} \varOmega_{3} }}{{1 + 2\zeta_{3} \varOmega_{3} + \varOmega_{3}^{2} }}} \\ \end{array} } \right].$$

(49)

Owing to diagonality of (49), simple matrix decomposition can be performed as follows:

$$\left[ {{\mathbf{\bar{\rm A}}}_{i} } \right] = \left[ {\begin{array}{*{20}c} {\frac{1}{{1 + 2\zeta_{1} \varOmega_{1} + \varOmega_{1}^{2} }}} & 0 & 0 \\ 0 & {\frac{1}{{1 + 2\zeta_{2} \varOmega_{2} + \varOmega_{2}^{2} }}} & 0 \\ 0 & 0 & {\frac{1}{{1 + 2\zeta_{3} \varOmega_{3} + \varOmega_{3}^{2} }}} \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {1 + 2\zeta_{1} \varOmega_{1} } & 0 & 0 \\ 0 & {1 + 2\zeta_{2} \varOmega_{2} } & 0 \\ 0 & 0 & {1 + 2\zeta_{3} \varOmega_{3} } \\ \end{array} } \right].$$

(50)

Or

$$\left[ {{\mathbf{\bar{\rm A}}}_{i} } \right] = \left[ {\begin{array}{*{20}c} {1 + 2\zeta_{1} \varOmega_{1} + \varOmega_{1}^{2} } & 0 & 0 \\ 0 & {1 + 2\zeta_{2} \varOmega_{2} + \varOmega_{2}^{2} } & 0 \\ 0 & 0 & {1 + 2\zeta_{3} \varOmega_{3} + \varOmega_{3}^{2} } \\ \end{array} } \right]^{ - 1} \left[ {\begin{array}{*{20}c} {1 + 2\zeta_{1} \varOmega_{1} } & 0 & 0 \\ 0 & {1 + 2\zeta_{2} \varOmega_{2} } & 0 \\ 0 & 0 & {1 + 2\zeta_{3} \varOmega_{3} } \\ \end{array} } \right]$$

(51)

$$\left[ {{\mathbf{\bar{\rm A}}}_{i} } \right] = \left[ {\begin{array}{*{20}c} {D_{{\left( {\alpha_{i1} } \right)}} } & 0 & 0 \\ 0 & {D_{{\left( {\alpha_{i2} } \right)}} } & 0 \\ 0 & 0 & {D_{{\left( {\alpha_{i3} } \right)}} } \\ \end{array} } \right]^{ - 1} \left[ {\begin{array}{*{20}c} {N_{{\left( {\alpha_{i1} } \right)}} } & 0 & 0 \\ 0 & {N_{{\left( {\alpha_{i2} } \right)}} } & 0 \\ 0 & 0 & {N_{{\left( {\alpha_{i3} } \right)}} } \\ \end{array} } \right] \to \left[ {{\mathbf{\bar{\rm A}}}_{i} } \right] = \left[ {D_{{\left( {{\mathbf{\bar{\rm A}}}_{i} } \right)}} } \right]^{ - 1} \left[ {N_{{\left( {{\mathbf{\bar{\rm A}}}_{i} } \right)}} } \right].$$

(52)

However, Eq. (52) involves modal parameters and in practice, it should be transformed to standard finite-element coordinates. This can be done by employing the transformations introduced in this paper [Eqs. (32) and (33)]. For example, $\left[ {N_{{\left( {{\mathbf{\bar{\rm A}}}_{i} } \right)}} } \right]$ could be written as

$$\begin{aligned} \left[ {N_{{\left( {{\mathbf{\bar{\rm A}}}_{i} } \right)}} } \right] & = \left[ {\begin{array}{*{20}c} {1 + 2\zeta_{1} \varOmega_{1} } & 0 & 0 \\ 0 & {1 + 2\zeta_{2} \varOmega_{2} } & 0 \\ 0 & 0 & {1 + 2\zeta_{3} \varOmega_{3} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array} } \right] + \left[ {\begin{array}{*{20}c} {2\zeta_{1} \varOmega_{1} } & 0 & 0 \\ 0 & {2\zeta_{2} \varOmega_{2} } & 0 \\ 0 & 0 & {2\zeta_{3} \varOmega_{3} } \\ \end{array} } \right] \\ & = \left[ {\mathbf{I}} \right] + 2\left[ {{\bar{\mathbf{\zeta }}}} \right]\left[ {{\bar{\mathbf{\varOmega }}}} \right] = \left[ {\varvec{\Phi}} \right]^{\text{T}} \left[ {\varvec{\Phi}} \right]^{{ - {\text{T}}}} + \Delta t\left[ {\varvec{\Phi}} \right]^{\text{T}} \left( {\left[ {\mathbf{C}} \right]\left[ {\mathbf{M}} \right]^{ - 1} } \right)\left[ {\varvec{\Phi}} \right]^{{ - {\text{T}}}} \to \left[ {N_{{\left( {{\mathbf{\bar{\rm A}}}_{i} } \right)}} } \right] \\ & = \left[ {\varvec{\Phi}} \right]^{\text{T}} \left( {\left[ {\mathbf{I}} \right] + \Delta t\left[ {\mathbf{C}} \right]\left[ {\mathbf{M}} \right]^{ - 1} } \right)\left[ {\varvec{\Phi}} \right]^{{ - {\text{T}}}} . \\ \end{aligned}$$

(53)

The above procedure could be performed for all of the integration parameters. Thus, they can be described only in terms of original system matrices and in the form of $\left[ {\varvec{\Phi}} \right]^{T} \left( \ldots \right)\left[ {\varvec{\Phi}} \right]^{{ - {\text{T}}}}$.