Insight Into Feasibility of Structure-Dependent Methods for Dynamic Analysis

Shuenn-Yih Chang

Department of Civil Engineering, National Taipei University of Technology, Taipei, Taiwan, Republic of China
E-mail: changsy@ntut.edu.tw

Received 20 August 2022; Accepted 17 January 2023; Publication 29 April 2023

Abstract

The first family of structure-dependent integration methods have been successfully developed for nonlinear dynamic analysis. Although its numerical properties were evaluated and its performance was numerically corroborated for both linear and nonlinear systems, its feasibility is still under debate due to the lack of a theoretical background. It seems that an eigen-based theory can provide a fundamental basis for the proof of the feasibility of structure-dependent integration methods. This can be manifested from each major stage of the development of structure-dependent integration methods. Therefore, the development of the first family of structure-dependent integration methods will be presented and the correlation between each major stage and an eigen-based theory will be explored and explained. Besides, this developing sequence can lay a typical procedure for developing a general structure-dependent integration method.

Keywords: An eigen-based theory, low frequency modes, high frequency modes, eigen-dependent method, structure-dependent method.

1 Introduction

An unconditionally stable, explicit Structure-Dependent Integration Method (SDIM) was first developed by Chang in 2002 [1] for the pseudo-dynamic testing. Subsequently, some improved integration methods of the same type were further developed [2–6]. SDIMs are different from conventional integration methods [7–19] in the coefficients of displacement and/or velocity difference equations. The coefficients are scalar constants for conventional integration methods while for SDIMs they can be functions of initial structural properties for defining the problem under analysis [1–6]. The most promising property of SDIMs is that it can combine unconditional stability and explicit formulation at the same time and thus it can be very computationally efficient for solving inertial problems although explicit formulation and unconditional stability are the repellent properties for a numerical method based on the Dahlquist theorem [20], which states that there does not exist any explicit method that is absolutely stable in the linear multistep methods.

Numerical properties of the first family of SDIMs that is explicit and unconditionally stable [3] have been explored for both linear elastic and nonlinear systems in the early stage of development and the feasibility of nonlinear performance is also numerically validated. However, the development details are undisclosed. In addition, there is still lack of a solid fundamental base for this development. This family of SDIM will be referred as the Chang Family Method (CFM). This method is a semi-explicit method since it involves an explicit displacement difference equation and an implicit velocity difference equation. An eigen-based theory is applicable to disclose why CFM can combine unconditional stability and explicitness of each step. Hence, it is of interest to apply it to gain an insight into CFM. In addition, a useful procedure for developing such a SDIM can be also constructed from the use of an eigen-based theory to each major step. Thus, a summary of developing procedure is described next. Conceptually, a coupled equation of motion can be first decomposed into a system of uncoupled modal equations of motion. Subsequently, an Eigen-Dependent Integration Method (EDIM) can be developed to solve each modal equation of motion, where the low frequency modes must be accurately integrated while no instability must be guaranteed for high frequency modes. Finally, all the EDIMs will be converted to a SDIM since this conversion can enable the feasibility of this type of integration methods for practical applications. Besides, this procedure can be used to develop a general SDIM, such as CFM. The development details of CFM will be presented in this work.

2 Concept of Eigen-based Theory

In structural dynamics or earthquake engineering, an equation of motion for a multiple degree of freedom system is a second-order ordinary differential equation and is used to govern the dynamic behaviors of the system. Although an analytical solution might be theoretically obtained based on the fundamental theory of structural dynamics for some specific problems, it is a difficult task or impossible to solve a dynamic problem with an arbitrary dynamic loading, such as an earthquake ground shaking, or a nonlinear structural system [21]. As a consequence, an integration method is generally adopted as an alternative for solving such problems. A structural dynamic problem is often classified as an inertial problem since its total response is often dominated by low frequency modes while high frequency responses contribute insignificantly. As a result, an implicit integration method is best suited to solving such a problem since it can accurately achieve low frequency results while no abnormal amplitude growth is found for high frequency modes although a significant period distortion is often found for a high frequency mode [22].

Clearly, the critical drawback of an implicit integration method for solving a general structural dynamics problem is the implicitness of each step since it requires an iteration procedure for each step. Thus, an integration method will be promising for solving inertial problems if it can combine explicit formulation and unconditional stability together at the same time. However, it is forbidden by Dahlquist barrier: there exists no linear multistep methods that can have explicit formulation and absolute stability simultaneously. There is a motivation to develop a new type of integration methods by using matrix coefficients to replace scalar coefficients for determining the displacement and/or velocity difference equations. This concept is intended to overcome the Dahlquist barrier that has been proved for linear multistep methods with scalar coefficients but not for matrix coefficients.

A coupled equation of motion can be decomposed into a set of uncoupled equations of motion by exploiting an eigen-decomposition technique. Hence, a total solution can be accurately obtained if the responses for the modes of interest are accurately integrated. This indicates that each modal equation of motion of interest must be reliably integrated by an integration method. Besides, it is found that the characteristic equation of a conventional integration method for solving a single degree of freedom system is, in general, a function of $\mathrm{\Omega }=\omega (\mathrm{\Delta }t)$, where $\omega$ is a natural frequency and $\mathrm{\Delta }t$ is a step size. It is natural to assume that the coefficients of the displacement and/or velocity difference equation are functions of ${\mathrm{\Omega }}^{(k)}={\omega }^{(k)}(\mathrm{\Delta }t)$, where ${\omega }^{(k)}$ is introduced to denote the natural frequency of the k-th the mode, since this assumption will also result in a characteristic equation that is also a function of ${\mathrm{\Omega }}^{(k)}$ for using the integration method to solve the k-th modal equation of motion. As a consequence, an EDIM is proposed to solve the modal equation of motion for the k-th mode since the eigendata of the k-th mode ${\omega }^{(k)}$ is involved in the formulation. Each modal equation of motion will be solved by the corresponding EDIM. In addition, modal solutions for low frequency modes must be accurately obtained while no instability experiences for high frequency modes.

Clearly, it is impractical to solve each modal equation of motion of interest by using a different EDIM. This is because an eigen-analysis must be conducted so that the eigendata of each mode can be obtained and it is time consuming for a matrix of large order. Besides, after solving each modal equation of motion of interest, it is also impractical to use a modal superposition method to sum up the modal solutions of interest. Notice that the eigendata may vary per step, which implies that an eigen-analysis must be performed as the structural properties vary. It will be too complex and too time consuming to conduct an eigen-analysis for each time step. As an alternative, it is very important to transform all the EDIMs into a SDIM by using a reverse procedure of an eigen-decomposition technique. This integration method is no longer eigen dependent but structure dependent after this transformation. Thus, a SDIM can be developed.

3 Development of CFM

After a literature review, it can be found that CFM is the first family of the SDIMs [3]. It has been shown that it can have the same characteristic equation as that of the Newmark family method (NFM) [7] for linear elastic systems. Hence, it exhibits the same numerical properties as those of NFM since these numerical properties are derived from the same characteristic equation. Besides, its feasibility for solving nonlinear structural systems were numerically confirmed. However, the detailed development of this family of SDIMs is still covered and the simultaneous combination of the unconditional stability and explicitness of each step is still unclear. Thus, both topics will be explored next.

3.1 Eigen-Decomposition of Equation of Motion

In structural dynamics, the governing equation for an n degree of freedom system can be, in general, expressed as:

$$\mathbf{\text{M}}{\ddot{\mathbf{\text{u}}}}_{i+1}+\mathbf{\text{C}}{\dot{\mathbf{\text{u}}}}_{i+1}+{\mathbf{\text{Ku}}}_{i+1}={\mathbf{\text{f}}}_{i+1}.$$

(1)

where M, C and K are the mass, viscous damping and stiffness matrices, respectively; and ${\mathbf{\text{u}}}_{i+1}$, ${\dot{\mathbf{\text{u}}}}_{i+1}$, ${\ddot{\mathbf{\text{u}}}}_{i+1}$ and ${\mathbf{\text{f}}}_{i+1}$ are in correspondence to the vectors of the displacement, velocity, acceleration and dynamic loading at the end of the (i+1)-th time step. After conducting an eigen-analysis of Equation (1), the natural frequencies ${\omega }^{(k)}$ and corresponding eigenvectors ${\varphi }^{(k)}$ can be found for $k=1,2,3,\mathrm{\dots },n$, where ${\omega }^{(1)}\le {\omega }^{(2)}\le \mathrm{\cdots }\le {\omega }^{(k)}$. In general, Equation (1) is a coupled equation of motion through the mass or stiffness matrix or both matrices and it can be decoupled into a series of uncoupled modal equations of motion:

$$m^{(k)}{\ddot{u}}^{(k)}+c^{(k)}{\dot{u}}^{(k)}+k^{(k)}{u}^{(k)}=f^{(k)},k=1,2,3,\mathrm{\dots },n.$$

(2)

where

$$\begin{array}{c}{m}^{(k)}={[{\varphi }^{(k)}]}^{\mathrm{T}}\mathbf{\text{M}}{\varphi }^{(k)},c^{(k)}={[{\varphi }^{(k)}]}^{\mathrm{T}}\mathbf{\text{C}}{\varphi }^{(k)}\hfill \\ k^{(k)}={[{\varphi }^{(k)}]}^{\mathrm{T}}\mathbf{\text{K}}{\varphi }^{(k)},f^{(k)}={[{\varphi }^{(k)}]}^{\mathrm{T}}\mathbf{\text{f}}(t)\hfill \end{array}.$$

(3)

for $k=1,2,\mathrm{\dots },n$. In addition, $m^{(k)}$, $c^{(k)}$, $k^{(k)}$ and $f^{(k)}$ are the generalized terms of mass, viscous damping coefficient, stiffness and dynamic loading for the k-th mode, respectively [21]. Notice that a classical damping matrix is assumed so that the coupled equations of motion can be decoupled in the derivation. However, the application of this method is not limited and will be shown later. In the subsequent study, the superscript (k) for representing the k-th mode will be deleted from its original expression for brevity although the k-th mode is still indicated.

Equation (2) is an uncoupled modal equation of motion and an integration method can be applied to solve it. To develop the novel type of SDIMs, an EDIM must be developed to solve each modal equation of motion of Equation (2) at first. Based on the concept of an eigen-based theory, either the displacement and/or velocity difference equation can be chosen to be eigen dependent or structure dependent [1]. In developing CFM, the displacement difference equation is assumed to be eigen dependent and the velocity difference equation used by NFM is adopted. As a result, an EDIM for solving Equation (2) can be expressed as:

$$\begin{array}{c}m{a}_{i+1}+c{v}_{i+1}+k{d}_{i+1}=f_{i+1}\hfill \\ d_{i+1}=d_i+{\beta }_1(\mathrm{\Delta }t)v_i+{\beta }_2{(\mathrm{\Delta }t)}^2{a}_i+\frac{1}{m}{(\mathrm{\Delta }t)}^2{p}_{i+1}\hfill \\ v_{i+1}=v_i+(\mathrm{\Delta }t)[(1-\gamma )a_i+\gamma a_{i+1}]\hfill \end{array}.$$

(4)

where ${\beta }_1$ and ${\beta }_2$ are eigen-dependent coefficients and $p_{i+1}$ is a load-dependent term; and $\gamma$ is a scalar constant. Clearly, the second line of this equation is an explicit displacement difference equation since it involves no current step data. Although the load-dependent term $p_{i+1}$ has never been seen in the difference equations of a conventional integration method, it is important to include it since it can be applied to eliminate an adverse overshoot in high frequency steady-state responses [23].

3.2 An Eigen-Dependent Integration Method

The core to successfully develop an EDIM is that it can accurately solve low frequency modes while no instability occurs for high frequency modes. For this purpose, the basic assumptions for ${\beta }_1$ and ${\beta }_2$ play a key role since both will be involved in the characteristic equation, which controls numerical properties, such as stability, accuracy and numerical damping. In addition, the behaviors in the limit of ${\mathrm{\Omega }}_0\to 0$ and ${\mathrm{\Omega }}_0\to \mathrm{\infty }$ indicate the behaviors of low and high frequency modes, respectively, where ${\mathrm{\Omega }}_0={\omega }_0(\mathrm{\Delta }t)$ is defined; ${\omega }_0=\sqrt{k_0/m}$ is an initial natural frequency and $k_0$ is an initial stiffness. Two simple guidelines can be used to assume the formulations of ${\beta }_1$ and ${\beta }_2$. One is to assume that ${\beta }_1$ and ${\beta }_2$ are fractions of ${\mathrm{\Omega }}_0$ and the other is to assume that the asymptotic values of ${\beta }_1$ and ${\beta }_2$ are scalar constants in the limit of ${\mathrm{\Omega }}_0\to 0$ and ${\mathrm{\Omega }}_0\to \mathrm{\infty }$. For simplicity, the numerator of ${\beta }_1$ and ${\beta }_2$ is assumed to be a linear polynomial function of ${\mathrm{\Omega }}_0$ while the denominator is a quadratic polynomial function of ${\mathrm{\Omega }}_0$. As a consequence, the formulation of ${\beta }_1$ and ${\beta }_2$ can be written as:

$${\beta }_i=\frac{1}{D}(g_i+h_{i}2\xi {\mathrm{\Omega }}_0),i=1,2.$$

(5)

where $g_i$ and $h_i$ are scalar constants and the denominator of $D=1+\gamma 2\xi {\mathrm{\Omega }}_0+\beta {\mathrm{\Omega }}_0^2$ is chosen. This choice of $D$ can be referred to NFM for simplicity. Notice that it can be also assumed in a different polynomial form or even other higher order forms. However, the highest power of ${\mathrm{\Omega }}_0$ in numerator must be no more than that of the denominator. This is because numerical explosions will occur if it is greater than that of the denominator.

Clearly, $g_i$ and $h_i$ must be properly determined so that desired numerical properties can be yielded. Since any numerical method must be a convergent method, the satisfaction of the convergence can be used to determine the scalar coefficients of $g_i$ and $h_i$. Besides, consistency together with stability can imply the convergence of a numerical method based on the Lax or Dahlquist equivalence theory [24, 25]. To prove consistency, one can show that the proposed EDIM has at least a first order accuracy. The detailed derivation of its local truncation error for a free vibration response can be found in [22]. As a result, it is:

$E$	$={\displaystyle \frac{1}{BD}}(1-g_1){\ddot{u}}_i+{\displaystyle \frac{1}{BD}}[(1-g_1)+(\gamma -h_1)2\xi \mathrm{\Omega }]2\xi {\omega }_0{\dot{u}}_i$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{1}{BD}}\left[\left(2\gamma -{\displaystyle \frac{1}{2}}-h_1\right)+\left(\gamma -{\displaystyle \frac{1}{2}}\right)\gamma 2\xi {\mathrm{\Omega }}_0\right]2\xi {\mathrm{\Omega }}_0{\ddot{u}}_i$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{1}{BD}}\{[g_2-g_1(1-\gamma )]+[\beta +h_2-(1-\gamma )h_1]2\xi {\mathrm{\Omega }}_0\}{\mathrm{\Omega }}_0{\omega }_0{\dot{u}}_i$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{1}{BD}}\{[\beta -{\displaystyle \frac{1}{2}}{g}_2+{\displaystyle \frac{1}{2}}{g}_1(1-\gamma )]$
	$\mathrm{\hspace{1em}}+[\beta \gamma -{\displaystyle \frac{1}{2}}{h}_2+{\displaystyle \frac{1}{2}}{h}_1(1-\gamma )]2\xi {\mathrm{\Omega }}_0\}\mathrm{\Omega }{}_0{}^2\ddot{u}{}_i$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{1}{6}}{\displaystyle \frac{1}{BD}}2\xi {\mathrm{\Omega }}_0(\mathrm{\Delta }t){\stackrel{\dot{}\dot{}\dot{}}{u}}_i+{\displaystyle \frac{1}{24}}{\displaystyle \frac{1}{BD}}2{(\mathrm{\Delta }t)}^2{\stackrel{\mathrm{\dots }.}{u}}_i+O[{(\mathrm{\Delta }t)}^3].$	(6)

where $B=1+\gamma 2\xi {\mathrm{\Omega }}_0$. In general, a competitive integration method must possess at least a second order accuracy and thus the coefficients of the terms that have an order of accuracy less than 2 must be equal to zero. Utilizing the first three terms on the right side of this equation, one can have $g_1=1$ and $h_1=\gamma$. Evidently, the first two terms disappear and the third term becomes $(\gamma -\frac{1}{2})2\xi {\mathrm{\Omega }}_0{\ddot{u}}_i/D$, where both numerator and denominator have the factor of B and thus it can be removed. Similarly, one can also assume that the numerator of the fourth term also possesses the factor of B. As a result, the following equation can be found:

$$[g_2-g_1(1-\gamma )]\gamma =\beta +h_2-(1-\gamma )h_1.$$

(7)

After substituting $g_1=1$ and $h_1=\gamma$ into this equation, $g_2\gamma =\beta +h_2$ is achieved and the fourth term will become $(g_2-1+\gamma ){\mathrm{\Omega }}_0{\omega }_0{\dot{u}}_i/D$. As a result, based on the assumption of $g_2-1+\gamma =\gamma -1/2$, the results of $g_2=1/2$ and $h_2=\gamma /2-\beta$ can be found. Thus, the eigen-dependent coefficients of ${\beta }_1$ and ${\beta }_2$ are summarized to be:

$${\beta }_1=\frac{1}{D}(1+\gamma 2\xi {\mathrm{\Omega }}_0),{\beta }_2=\frac{1}{D}\left[\frac{1}{2}+\left(\frac{1}{2}\gamma -\beta \right)2\xi {\mathrm{\Omega }}_0\right].$$

(8)

Clearly, both ${\beta }_1$ and ${\beta }_2$ are eigen-dependent coefficients. At first glance, it seems that only consistency is required to determine the scalar constant coefficients of ${\beta }_1$ and ${\beta }_2$ while the condition of stability is not involved up to this time. In fact, it has been implicitly adopted in the assumed general formulation of ${\beta }_1$ and ${\beta }_2$, where the maximum order of numerator must be no more than that of denominator so that instability can be automatically avoided.

After determining ${\beta }_1$ and ${\beta }_2$, the load-dependent term $p_{i+1}$ must be further derived and it can be also determined from a local truncation error. In this case, a local truncation error is no longer determined from a free vibration response but a forced vibration response since the dynamic loading must be considered for determining the order of accuracy of the proposed EDIM. As a result, it is found to be:

$E$	$=-{\displaystyle \frac{1}{D}}\left(\gamma -{\displaystyle \frac{1}{2}}\right)(\mathrm{\Delta }t){\stackrel{\dot{}\dot{}\dot{}}{u}}_i+{\displaystyle \frac{1}{D}}\beta {\mathrm{\Omega }}_0^2{\ddot{u}}_i$
	$\mathrm{\hspace{1em}}-{\displaystyle \frac{1}{D}}{\displaystyle \frac{1}{2}}(\gamma -{\displaystyle \frac{1}{2}})[2\xi {\mathrm{\Omega }}_0(\mathrm{\Delta }t){\stackrel{\dot{}\dot{}\dot{}}{u}}_i+{(\mathrm{\Delta }t)}^2{\stackrel{\mathrm{\dots }.}{u}}_i]$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{1}{D}}\left({\displaystyle \frac{1}{6}}\right)2\xi {\mathrm{\Omega }}_0(\mathrm{\Delta }t){\stackrel{\dot{}\dot{}\dot{}}{u}}_i+{\displaystyle \frac{1}{D}}{\displaystyle \frac{1}{12}}{(\mathrm{\Delta }t)}^2{\stackrel{\mathrm{\dots }.}{u}}_i+{\displaystyle \frac{1}{BD}}\beta 2\xi {\mathrm{\Omega }}_0{\displaystyle \frac{1}{m}}(\mathrm{\Delta }t){\dot{f}}_i$
	$\mathrm{\hspace{1em}}-{\displaystyle \frac{1}{BD}}{\displaystyle \frac{1}{2}}\beta 2\xi {\mathrm{\Omega }}_0{\displaystyle \frac{1}{m}}{(\mathrm{\Delta }t)}^2{\ddot{f}}_i+\left(1-{\displaystyle \frac{1}{B}}2\xi {\mathrm{\Omega }}_0\right)p_i-p_{i+1}+O[{(\mathrm{\Delta }t)}^3].$	(9)

This equation reveals that CFM can generally possess a first order accuracy for zero dynamic load while a second order accuracy can be achieved for $\gamma =1/2$. However, it can only have a first order accuracy for nonzero dynamic loading even for $\gamma =1/2$. The term of $\beta {\mathrm{\Omega }}_0^2{\ddot{u}}_i/D$ is the only quadratic error term of ${\mathrm{\Omega }}_0$ for $\gamma =1/2$ and it has been shown that it will lead to an adverse overshoot in high frequency steady-state responses [23]. Thus, it must be removed from the local truncation error to eliminate this overshoot. Since the displacement and velocity difference equations are involved in the derivation of a local truncation error, the addition of the load-dependent term into the displacement difference equation is intended to remove the dominant error term for CFM.

For this purpose, the following relation will be applied to take out the adverse dominant error term of $\beta {\mathrm{\Omega }}_0^2{\ddot{u}}_i/D$:

$${(\mathrm{\Delta }t)}^2{\stackrel{\mathrm{\dots }.}{u}}_i+2\xi {\mathrm{\Omega }}_0(\mathrm{\Delta }t){\stackrel{\dot{}\dot{}\dot{}}{u}}_i=\frac{1}{m}{(\mathrm{\Delta }t)}^2{\ddot{f}}_i-{\mathrm{\Omega }}_0^2{\ddot{u}}_i.$$

(10)

Clearly, this equation is directly derived from the second time derivative of the equation of motion. It can be found that the maximum order of ${\mathrm{\Omega }}_0$ is 1 on the left side while on the right side it is 2. As a result, the left two terms can replace the right two terms and then the quadratic term of $\beta {\mathrm{\Omega }}_0^2{\ddot{u}}_i/D$ can be removed and the order of ${\mathrm{\Omega }}_0$ is reduced. As a result, an adverse overshoot can be eliminated. It is worth noting that the equation of motion and its first time derivative has been used in determining the scalar constant coefficients of ${\beta }_1$ and ${\beta }_2$. Finally, the load-dependent term can be derived from Equation (3.2) and is found to be:

$$p_{i+1}=\frac{1}{D}\beta (f_{i+1}-f_i).$$

(11)

Finally, after substituting Equation (11) into Equation (3.2), the final form of the local truncation error for the proposed EDIM will become:

$E$	$=-{\displaystyle \frac{1}{D}}\left(\gamma -{\displaystyle \frac{1}{2}}\right)(\mathrm{\Delta }t){\stackrel{\dot{}\dot{}\dot{}}{u}}_i+{\displaystyle \frac{1}{D}}\left[{\displaystyle \frac{1}{6}}-\beta +{\displaystyle \frac{1}{2}}\left(\gamma -{\displaystyle \frac{1}{2}}\right)\right]2\xi {\mathrm{\Omega }}_0(\mathrm{\Delta }t){\stackrel{\dot{}\dot{}\dot{}}{u}}_i$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{1}{D}}\left[{\displaystyle \frac{1}{12}}-\beta +{\displaystyle \frac{1}{2}}\left(\gamma -{\displaystyle \frac{1}{2}}\right)\right]{(\mathrm{\Delta }t)}^2{\stackrel{\mathrm{\dots }.}{u}}_i+O[{(\mathrm{\Delta }t)}^3].$	(12)

It is disclosed by this equation that the proposed EDIM can generally have a first order accuracy for classical damping and dynamic loading. A second order accuracy can be further achieved if $\gamma =1/2$ is adopted. As a result, an EDIM for solving the modal equation of motion of the k-th mode is developed. This EDIM also implies that a total number of n EDIMs are developed for solving the n modal equations of motion for $k=1,2,\mathrm{\dots },n$, respectively.

It is important to examine whether the proposed EDIM can satisfy the major goal to successfully develop an EDIM. Therefore, the limiting cases of ${\mathrm{\Omega }}_0\to 0$ and ${\mathrm{\Omega }}_0\to \mathrm{\infty }$ are considered. In the limit ${\mathrm{\Omega }}_0\to 0$, the displacement difference equation will reduce to be:

$$d_{i+1}=d_i+(\mathrm{\Delta }t)v_i+\frac{1}{2}{(\mathrm{\Delta }t)}^2{a}_i.$$

(13)

Clearly, it can be manifested from Equation (4) that the proposed EDIM becomes the Newmark explicit method for any $\beta$ and $\gamma =1/2$. Thus, the performance of the proposed EDIM will be the roughly same as the Newmark explicit method in the limit ${\mathrm{\Omega }}_0\to 0$ or low frequency modes. On the other hand, in the limit ${\mathrm{\Omega }}_0\to \mathrm{\infty }$, the displacement difference equation will become:

$$d_{i+1}=d_i.$$

(14)

It can be found that the characteristic equation of Equation (4) with the displacement difference equation given in Equation (14) is:

$$\left[{\lambda }^2-\left(2-\frac{2\xi {\mathrm{\Omega }}_0}{1+2\gamma \xi {\mathrm{\Omega }}_0}\right)\lambda +\left(1-\frac{2\xi {\mathrm{\Omega }}_0}{1+2\gamma \xi {\mathrm{\Omega }}_0}\right)\right]\lambda =0.$$

(15)

It is very important to find that the characteristic equation in the limit ${\mathrm{\Omega }}_0\to \mathrm{\infty }$ has three different eigenvalues. In fact, they are found to be:

$${\lambda }_1=0,{\lambda }_2=1-\frac{\gamma 2\xi {\mathrm{\Omega }}_0}{1+\gamma 2\xi {\mathrm{\Omega }}_0},{\lambda }_3=1.$$

(16)

Since all the three eigenvalues are less than or equal to 1, an unconditional stability can be achieved in the limit ${\mathrm{\Omega }}_0\to \mathrm{\infty }$ or high frequency modes. Hence, the proposed EDIM can meet the major goal to successfully develop an EDIM. In addition, it will exhibit no weak instability since it can have three linearly independent eigenvectors in the limiting case of ${\mathrm{\Omega }}_0\to \mathrm{\infty }$ [26].

To decompose the equation of motion into a set of uncoupled equations of motion, it is required to assume that the damping matrix is a classical damping matrix. However, the extreme cases of ${\mathrm{\Omega }}_0\to 0$ and ${\mathrm{\Omega }}_0\to \mathrm{\infty }$ seem to reveal that the assumption is unnecessary. In the limit ${\mathrm{\Omega }}_0\to 0$ or low frequency modes, $D=1+\gamma 2\xi {\mathrm{\Omega }}_0+\beta {\mathrm{\Omega }}_0^2$ is controlled by the constant of 1 and is unaffected by damping. On the other hand, in the limit ${\mathrm{\Omega }}_0\to \mathrm{\infty }$ or high frequency modes, $D=1+\gamma 2\xi {\mathrm{\Omega }}_0+\beta {\mathrm{\Omega }}_0^2$ is dominated by the quadratic term of $\beta {\mathrm{\Omega }}_0^2$ and is nothing to do with damping. As a consequence, this assumption is to give a rigorous fundamental basis for derivation but it is unnecessary for practical applications.

Using the fundamental theory of structural dynamics, one can have $(k_0/m){(\mathrm{\Delta }t)}^2={\mathrm{\Omega }}_0^2$ and $(c_0/m)(\mathrm{\Delta }t)=2\xi {\mathrm{\Omega }}_0$. They can be applied to convert an EDIM to a generalized modal equation of motion [21]. As a result, ${\beta }_1$, ${\beta }_2$ and $p_{i+1}$ become:

${\beta }_1$	$={\displaystyle \frac{1}{D}}[m+\gamma (\mathrm{\Delta }t)c_0]$
${\beta }_2$	$={\displaystyle \frac{1}{D}}\left[{\displaystyle \frac{1}{2}}m+\left({\displaystyle \frac{1}{2}}\gamma -\beta \right)(\mathrm{\Delta }t)c_0\right]$	(17)
$p_{i+1}$	$={\displaystyle \frac{1}{D}}\beta (f_{i+1}-f_i).$

for the EDIM for solving the modal equation of motion of the k-th mode for $k=1,2,3,\mathrm{\dots },n$. As a summary, a total number of n EDIMs are developed to solve the n uncoupled modal equations of motion correspondingly. Besides, the low frequency modal equations of motion can be solved almost as accurately as those by the Newmark explicit method while no abnormal amplitude growth in high frequency modal equations of motion is guaranteed.

3.3 A Structure-Dependent Integration Method

Although an EDIM is successfully developed for solving each modal equation of motion, it is still inconvenient and time consuming if using a modal superposition scheme to sum up each modal response to yield a complete solution. An eigen-analysis must be frequently conducted as the structural properties of the coupled equation of motion might vary during time integration. This is because that eigendata will be involved both in developing an EDIM and in adding up all modal responses. This seems impractical for practical applications.

It is learned from an eigen-decomposition technique that it can decompose a coupled equation of motion into a series of uncoupled equations of motion. As a consequence, there is a great motivation to apply the reverse procedure of the eigen-decomposition technique to transform all the developed EDIMs into a SDIM so that the time consuming eigen-analysis can be avoided. As a consequence, the following equation of Equation (18) can be yielded after the transformation. The reverse of the eigen-decomposition technique implies that Equation (4) can be obtained after Equation (18) is decomposed [2, 5, 27].

$$\begin{array}{c}{\mathbf{\text{Ma}}}_{i+1}+{\mathbf{\text{C}}}_{i+1}{\mathbf{\text{v}}}_{i+1}+{\mathbf{\text{K}}}_{i+1}{\mathbf{\text{d}}}_{i+1}={\mathbf{\text{f}}}_{i+1}\hfill \\ {\mathbf{\text{d}}}_{i+1}={\mathbf{\text{d}}}_i+{\mathbf{\text{B}}}_1(\mathrm{\Delta }t){\mathbf{\text{v}}}_i+{\mathbf{\text{B}}}_2{(\mathrm{\Delta }t)}^2{\mathbf{\text{a}}}_i+{(\mathrm{\Delta }t)}^2{\mathbf{\text{p}}}_{i+1}\hfill \\ {\mathbf{\text{v}}}_{i+1}={\mathbf{\text{v}}}_i+(\mathrm{\Delta }t)[(1-\gamma ){\mathbf{\text{a}}}_i+\gamma {\mathbf{\text{a}}}_{i+1}]\hfill \end{array}.$$

(18)

where

$$\begin{array}{c}{\mathbf{\text{B}}}_1={\mathbf{\text{D}}}^{-1}[\mathbf{\text{M}}+\gamma (\mathrm{\Delta }t){\mathbf{\text{C}}}_0]\hfill \\ {\mathbf{\text{B}}}_2={\mathbf{\text{D}}}^{-1}\left[\frac{1}{2}\mathbf{\text{M}}+\left(\frac{1}{2}\gamma -\beta \right)(\mathrm{\Delta }t){\mathbf{\text{C}}}_0\right]\hfill \\ {\mathbf{\text{p}}}_{i+1}=\beta {\mathbf{\text{D}}}^{-1}({\mathbf{\text{f}}}_{i+1}-{\mathbf{\text{f}}}_i)\hfill \end{array}.$$

(19)

where $\mathbf{\text{D}}=\mathbf{\text{M}}+\gamma (\mathrm{\Delta }t){\mathbf{\text{C}}}_0+\beta {(\mathrm{\Delta }t)}^2{\mathbf{\text{K}}}_0$. Note that ${\mathbf{\text{C}}}_0$ and ${\mathbf{\text{K}}}_0$ are in correspondence to the initial damping coefficient matrix and the initial stiffness matrix. Clearly, ${\mathbf{\text{C}}}_0$ and ${\mathbf{\text{K}}}_0$ may be different from ${\mathbf{\text{C}}}_{i+1}$ and ${\mathbf{\text{K}}}_{i+1}$ for a nonlinear system. The displacement difference equation is a function of the initial structural properties, such as M, ${\mathbf{\text{C}}}_0$ and ${\mathbf{\text{K}}}_0$, and the step size $\mathrm{\Delta }t$. Hence, CFM is a SDIM. In addition, the displacement difference equation can be completely determined from the previous step data and thus CFM is an explicit integration method.

4 Dual Implementations for CFM

An explicit formulation for a SDIM is addressed in the development of structure-dependent methods since it has a better computational efficiency in contrast to an implicit formulation. However, it is worth noting that a SDIM can also have an implicit formulation. In general, structure-dependent coefficients are determined from the initial structural properties for a SDIM and then it can be explicitly implemented. In contrast, if initial structural properties are replaced by current structural properties to derive structure-dependent coefficients for a SDIM, an implicit formulation of the SDIM can be obtained and it can be implicitly implemented. Both implementations for CFM are shown next.

4.1 An Explicit Implementation

After conducting the time integration of the i-th time step, the displacement vector for the next time step can be determined from the second line of Equation (18) and is numerically equivalent to solve the following equation:

${\mathbf{\text{Dd}}}_{i+1}$	$={\mathbf{\text{Dd}}}_i+[\mathbf{\text{M}}+\gamma (\mathrm{\Delta }t){\mathbf{\text{C}}}_0](\mathrm{\Delta }t){\mathbf{\text{v}}}_i$
	$\mathrm{\hspace{1em}}+\left[{\displaystyle \frac{1}{2}}\mathbf{\text{M}}-\left(\beta -{\displaystyle \frac{1}{2}}\gamma \right)(\mathrm{\Delta }t){\mathbf{\text{C}}}_0\right]{(\mathrm{\Delta }t)}^2{\mathbf{\text{a}}}_i.$
	$\mathrm{\hspace{1em}}+\beta {(\mathrm{\Delta }t)}^2({\mathbf{\text{f}}}_{i+1}-{\mathbf{\text{f}}}_i)$	(20)

After obtaining ${\mathbf{\text{d}}}_{i+1}$, an assumed force-displacement relation can be used to determine the restoring force vector ${\mathbf{\text{r}}}_{i+1}={\mathbf{\text{K}}}_{i+1}{\mathbf{\text{d}}}_{i+1}$. Next, the velocity vector can be determined from the first and third lines of Equation (18) after eliminating the acceleration vector and is numerically equivalent to solve:

$[\mathbf{\text{M}}+\gamma (\mathrm{\Delta }t){\mathbf{\text{C}}}_{i+1}]{\mathbf{\text{v}}}_{i+1}$

$=\mathbf{\text{M}}[{\mathbf{\text{v}}}_i+(1-\gamma )(\mathrm{\Delta }t){\mathbf{\text{a}}}_i]+\gamma (\mathrm{\Delta }t)({\mathbf{\text{f}}}_{i+1}-{\mathbf{\text{r}}}_{i+1}).$

Similarly, the damping force vector can be alternatively denoted by ${\mathbf{\text{s}}}_{i+1}={\mathbf{\text{C}}}_{i+1}{\mathbf{\text{v}}}_{i+1}$ for an assumed nonlinear damping type. Finally, the acceleration vector can be calculated by:

$${\mathbf{\text{Ma}}}_{i+1}={\mathbf{\text{f}}}_{i+1}-{\mathbf{\text{s}}}_{i+1}-{\mathbf{\text{r}}}_{i+1}.$$

(22)

Since the matrices of D and M are invariant, a direct elimination method can be applied to solve Equation (4.1) and (22). Notice that Equation (4.1) must be iteratively solved if a nonlinear damping matrix is adopted. However, a constant damping matrix is generally used in nonlinear dynamic analysis, i.e., ${\mathbf{\text{C}}}_{i+1}={\mathbf{\text{C}}}_0$ is taken. Hence, a direct elimination method can be also used to solve Equation (4.1). An elimination method consists of a triangulation and a substitution; and a triangulation consumes most computational efforts. Only a triangulation is needed if D, $\mathbf{\text{M}}+\gamma (\mathrm{\Delta }t){\mathbf{\text{C}}}_0$ and M are invariant.

4.2 An Implicit Implementation

It is seen that Equation (4.1) will become a nonlinear equation if current structural properties are used to replace the initial structural properties. Thus, $\mathbf{\text{D}}=\mathbf{\text{M}}+\gamma (\mathrm{\Delta }t){\mathbf{\text{C}}}_0+\beta {(\mathrm{\Delta }t)}^2{\mathbf{\text{K}}}_0$ becomes $\mathbf{\text{D}}=\mathbf{\text{M}}+\gamma (\mathrm{\Delta }t){\mathbf{\text{C}}}_{i+1}+\beta {(\mathrm{\Delta }t)}^2{\mathbf{\text{K}}}_{i+1}$. As a result, an iteration procedure must be adopted to solve Equation (4.1) since ${\mathbf{\text{C}}}_{i+1}$ is generally a function of the current step data. Next, the calculations of the velocity and acceleration vectors are the same as those of the explicit implementation and will not be elaborated again.

It is important to note that an explicit implementation or an implicit implementation of CFM can have a comparable accuracy. This is because that CFM meets the major goal to successfully develop an EDIM, where low frequency modes can be accurately integrated while there exists no numerical instability in high frequency modes. This also indicates that the numerical accuracy of CFM is almost unaffected by using $\gamma (\mathrm{\Delta }t){\mathbf{\text{C}}}_0+\beta {(\mathrm{\Delta }t)}^2{\mathbf{\text{K}}}_0$ or $\gamma (\mathrm{\Delta }t){\mathbf{\text{C}}}_{i+1}+\beta {(\mathrm{\Delta }t)}^2{\mathbf{\text{K}}}_{i+1}$ to compute D while both choices can lead to a stable solution for high frequency modes. Therefore, an explicit implementation is preferred over an implicit implementation for CFM due to the saving of computational efforts.

Notice that an implicit implementation of CFM must be adopted for solving systems with nonlinear damping forces that are velocity dependent. This is because that an implicit difference equation is adopted by CFM for calculating the next step velocity, i.e., the third line of Equation (4) or (18). In this case, the implicit difference equation cannot be directly used to compute the next step velocity. As a result, an iterative procedure must be involved.

5 Comparisons of Accuracy

It is shown in the previous section that CFM can be both explicitly and implicitly implemented for time integration and an unconditional stability can be achieved for both implementations for solving general structural dynamics problems. Thus, it is of interest to explore the accuracy of NFM and the two implementations of CFM for both linear elastic and nonlinear systems. For this purpose, a parameter, which is named instantaneous degree of nonlinearity ${\delta }_i$, has been defined to monitor the stiffness change for a nonlinear system [2, 5]. In general, for the i-th time step, it can be written as ${\delta }_i=k_i/k_0$, where $k_i$ is the stiffness of the i-th time step while $k_0$ is the initial stiffness.

The general formulation of the well-known Newmark family of integration methods for a single degree of freedom system can be written as:

$$\begin{array}{c}m{a}_{i+1}+c{v}_{i+1}+k{d}_{i+1}=f_{i+1}\hfill \\ d_{i+1}=d_i+(\mathrm{\Delta }t)v_i+{(\mathrm{\Delta }t)}^2[(1-\beta )a_i+\beta a_{i+1}]\hfill \\ v_{i+1}=v_i+(\mathrm{\Delta }t)[(1-\gamma )a_i+\gamma a_{i+1}]\hfill \end{array}.$$

(23)

Besides, its characteristic equation for solving a general nonlinear system is found to be:

$$\lambda \left\{{\lambda }^2-\left[2-\frac{\left(\gamma +\frac{1}{2}\right){\mathrm{\Omega }}_{i+1}^2}{1+\beta {\mathrm{\Omega }}_{i+1}^2}\right]\lambda +\left[1-\frac{\left(\gamma -\frac{1}{2}\right){\mathrm{\Omega }}_{i+1}^2}{1+\beta {\mathrm{\Omega }}_{i+1}^2}\right]\right\}=0.$$

(24)

In the subsequent study, an Explicit implementation of CFM is referred as ECFM while its Implicit implementation is referred as ICFM. The characteristic equation of CFM is also derived for comparison. It is interesting to find that the characteristic equation of ICFM is the same as that of NFM. Whereas, that of ECFM has a different equation and is:

$$\lambda \left\{{\lambda }^2-\left[2-\frac{\left(\gamma +\frac{1}{2}\right){\mathrm{\Omega }}_{i+1}^2}{1+\beta {\mathrm{\Omega }}_0^2}\right]\lambda +\left[1-\frac{\left(\gamma -\frac{1}{2}\right){\mathrm{\Omega }}_{i+1}^2}{1+\beta {\mathrm{\Omega }}_0^2}\right]\right\}=0.$$

(25)

After comparing Equation (24) with (25), it is evident that the only difference between NFM and ECFM is the denominator. Notice that both characteristic equations will become the same for a linear elastic system.

Figure 1 Variations of relative period errors with $\mathrm{\Delta }t/T_0$.

To assess the numerical accuracy of NFM and ECFM for solving nonlinear systems, variations of the relative period errors versus $\mathrm{\Delta }t/T_0$, where $T_0=2\pi /{\omega }_0$, are calculated and plotted in Figure 1. The results of $\beta =1/4$ and $\gamma =1/2$ are shown in Figure 1(a) while those of $\beta =\gamma =1/2$ are plotted in Figure 1(b). It is seen in Figure 1(a) that the relative period error is small as $\mathrm{\Delta }t/T_0\le 1/20$ for ${\delta }_{i+1}=0.5$, 1.0 and 2.0 for both NFM and ECFM. Thus, it is evident that ECFM can have a comparable accuracy with NFM. It is also found that ECFM has the same period distortion as that of NFM for ${\delta }_{i+1}=1.0$. Whereas, the relative period error of ECFM for ${\delta }_{i+1}=0.5$ is larger than that of NFM while it is less than that of NFM for ${\delta }_{i+1}=2.0$. This implies that ECFM can have a better performance than for NFM for stiffness hardening systems while NFM performs better for stiffness softening systems. Similar phenomena are also found in Figure 1(b) except that the member of $\beta =\gamma =1/2$ exhibits a slightly larger period distortion in contrast to the member of $\beta =1/4$ and $\gamma =1/2$.

To corroborate the analytical results as shown in Figure 1, a simplified Duffing equation is solved and it can be expressed as:

$$\ddot{u}+u+\alpha u^3=0.$$

(26)

In general, $\alpha >0$ can simulate a hardening spring while a softening spring is mimicked by $\alpha <0$. Notice that $\alpha =0$ implies a linear elastic spring. Thus, the cases of $\alpha =-0.5$, 0 and 0.5 are used to simulate a softening spring, a linear elastic spring and a hardening spring correspondingly. The member of $\beta =1/4$ and $\gamma =1/2$ for NFM, ECFM and ICFM is chosen to calculate the free vibration responses of $u(0)=1$ and $\dot{u}(0)=0$. A time step of $\mathrm{\Delta }t=0.5\sec$ is used for time integration and the numerical results are plotted in Figures 2 to 4. Notice that a relative tolerance of ${10}^{-6}$ is used for NFM and ICFM in this paper. Figure 2(b) reveals that a softening spring is simulated since ${\delta }_{i+1}$ varies in the interval of $0.5\le {\delta }_{i+1}\le 1$. Whereas, a hardening spring is revealed by Figure 4(b) since ${\delta }_{i+1}$ varies in the interval of $1\le {\delta }_{i+1}\le 1.5$. Clearly, Figure 3(b) implies a linear elastic spring since ${\delta }_{i+1}=1$ is found. It is seen in Figures 2(a) to 4(a) that the results obtained from ICFM are overlapped with those obtained from NFM. This implies that the performance of ICFM is the same as that of NFM for both linear elastic and nonlinear systems. It is also found in Figure 2(a) that the result obtained from NFM is more accurate than that obtained from ECFM for the softening spring while ECFM can give a more accurate solution than that obtained from NFM for the hardening spring in Figure 4(a). The calculated results of NFM and ECFM coincide together for linear elastic spring in Figure 3(a). Clearly, these numerical results are consistent with the analytical predictions in Figure 1, where ECFM has more period distortion for softening systems and less period distortion for hardening systems in contrast NFM.

Figure 2 Free vibration responses to a softening spring.

Figure 3 Free vibration responses to a linear elastic spring.

Figure 4 Free vibration responses to a hardening spring.

The results for plotting Figures 2 to 4 are calculated from $\mathrm{\Delta }t=0.5\sec$, which is in correspondence to $\mathrm{\Delta }t/T_0=0.08$. As a consequence, these results considerably deviate from the reference solutions. More accurate solutions can be obtained for NFM, ICFM and ECFM if a smaller time step is adopted, such as a time step of $\mathrm{\Delta }t=0.2\sec$. This time step corresponds to $\mathrm{\Delta }t/T_0=0.03$. The results obtained from $\mathrm{\Delta }t=0.2\sec$ are plotted in Figure 5. Since $\mathrm{\Delta }t/T_0\le 1/20$ can be considered as a rule of thumb to yield an accurate solution, it is anticipated that reliable solutions can be achieved for using a step size of $\mathrm{\Delta }t=0.2\sec$. It is evident from Figure 5 that ECFM can have almost the same performance as that of NFM in the step-by-step solution of either linear elastic or nonlinear systems.

Figure 5 Free vibration responses to Duffing equations.

6 Applications

6.1 A Non-classical Damping System

To affirm that problem-dependent methods can be used to solve systems with non-classical damping, a second order equation of motion with a non-classical damping is solved and it can be expressed as:

$$\ddot{u}+\frac{1}{1+|\dot{u}|}\dot{u}+25u=0.$$

(27)

where the coefficient of the velocity term $\dot{u}$ is $1/(1+|\dot{u}|)$. It is evident that the damping force is not only non-classical but also nonlinear. Therefore, this equation must be solved iteratively by using ICFM and ECFM is inappropriate due to the adoption of an implicit velocity difference equation.

Figure 6 Numerical solution of single degree of freedom system with a non-classical damping force.

An initial natural frequency of the system is found to be 5 rad/sec. The initial condition of $u(0)=0$ and $\dot{u}(0)=10$ is considered in the subsequent calculations. Evidently, the nonlinearity of the system comes from the change of $\dot{u}$. The reference solution is computed from AAM with a step size of $\mathrm{\Delta }t=0.001$ sec for comparison. On the other hand, both AAM and ICFM are also applied to solve the same problem with the step sizes of $\mathrm{\Delta }t=0.05$ sec and 0.03 sec. All calculated results are plotted in Figure 6. It is seen in Figure 6(a) that the solution obtained from ICFM almost coincides together with that calculated by AAM with the time step of $\mathrm{\Delta }t=0.05$ sec although considerable errors are found for both integration methods when compared to the reference solution. On the other hand, an accurate solution is yielded for both methods if a smaller time step of $\mathrm{\Delta }t=0.03$ is adopted to compute the results as shown in Figure 6(b). Consequently, it is affirmed that ICFM can have a comparable accuracy with AAM in the solution of the second order ODE with a non-classical type of nonlinear damping force.

Figure 7 A multistory building of 10 stories.

6.2 A Multi-story Building of 10 Stories

To show that CFM can be used to solve a system with a non-classical damping type for a multiple degree of freedom system. A multi-story building of 10 stories is considered as shown in Figure 7. For simplicity, the building is mathematically modelled as a shear building. This simplified model is very useful to yield dynamic data for the preliminary design. The damping force of each story is nonlinear velocity dependent and is a non-classical type. Meanwhile, the spring force of each story is nonlinear displacement dependent. The lumped mass of each story is taken to be $m_i={10}^5$ kg. The damping coefficient of each story is assumed as:

$c_1$	$={10}^6(1-0.1\sqrt{|{\dot{u}}_1|})N-s/m$
$c_i$	$={10}^6\left[1-0.1\left({\displaystyle \frac{u_i-u_{i-1}}{u_i+u_{i-1}}}\right)\sqrt{|{\dot{u}}_i-{\dot{u}}_{i-1}|}\right]N-s/m.$	(28)

where $i=2,3,\mathrm{\dots },10$. It is evidently that a nonlinear damping force is not only velocity dependent but also displacement dependent. Notice that some buildings were equipped with a variety of nonlinear viscous and/or viscoelastic dampers to largely absorb seismic energy and thus to suppress dynamic responses [28, 29]. Such dampers are often classified as velocity-dependent dampers since the damping forces are mainly a function of the relative velocity or the frequency of motion [30]. On the other hand, the stiffness of each story is:

$k_1$	$={10}^8(1-2\sqrt{|u_1|})N/m$
$k_i$	$={10}^8(1-2\sqrt{|u_i-u_{i-1}|})N/m.$	(29)

where $i=2,3,\mathrm{\dots },10$. In general, the story stiffness will be degraded if $u_i-u_{i-1}\ne 0$.

The equation of motion for the multi-story building as shown in Figure 7 can be written as shown in Equation (1), where the mass matrix is a diagonal matrix and each diagonal term is $m_i={10}^4$ kg, for $i=1,2,\mathrm{\dots },10$; the damping matrix C and stiffness matrix K can be explicitly written as:

C	$=\left[\begin{array}{cccc}\hfill c_1+c_2\hfill & \hfill -c_2\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -c_2\hfill & \hfill c_2+c_3\hfill & \hfill -c_3\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -c_3\hfill & \hfill \mathrm{\cdots }\hfill & \hfill -c_n\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -c_n\hfill & \hfill c_n\hfill \end{array}\right],$
K	$=\left[\begin{array}{cccc}\hfill k_1+k_2\hfill & \hfill -k_2\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -k_2\hfill & \hfill k_2+k_3\hfill & \hfill -k_3\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -k_3\hfill & \hfill \mathrm{\cdots }\hfill & \hfill -k_n\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -k_n\hfill & \hfill k_n\hfill \end{array}\right].$	(30)

where C is a nonlinear function of displacement as well as velocity and K is a nonlinear function of displacement. The smallest and largest initial natural frequencies of the building are found to be 4.73 and 62.54 rad/sec, respectively. This building is excited by a ground acceleration record of Chi-Chi earthquake in 1999 in Taiwan. TCU129 is chosen as the seismic input and it has a peak ground acceleration of 989.2 cm/sec${}^2$. In this study, the peak ground acceleration of TCU129 is scaled to 1g.

Figure 8 Seismic responses of a multi-story building of 10 stories.

The seismic responses calculated from AAM with the time step of $\mathrm{\Delta }t=0.005$ sec is treated as a reference solution. Besides, both AAM and ICFM with $\mathrm{\Delta }t=0.05$ sec are also used to compute the lateral displacement responses. All the numerical results of the 1st, 5th and 10th floors are shown in Figure 8. It is revealed by the three plots that ICFM gives almost the same solutions as those calculated from AAM. This implies that both ICFM and AAM can faithfully seize the nonlinear variations of the damping and restoring forces due to a complicated ground shaking arising from an earthquake. Hence, it is substantiated that CFM can be also used to solve multiple degree of freedom systems with a damping force that is non-classical and nonlinear.

Figure 9 A 1000-DOF spring-mass system.

6.3 A Series of Large Spring-Mass Systems

A 1000-degree-of-freedom system is particularly chosen for examining the performance of ECFM in the solution of nonlinear multiple degree of freedom systems. The system is shown in Figure 9, where the lumped mass $m_l$ and the spring constant of $k^{(l)}$ are correspondent to the l-th degree of freedom. Clearly, the spring constant $k^{(l)}$ will decrease for a nonzero story drift, i.e., $|u^{(l)}-u^{(l-1)}|\ne 0$. It is found that the initial lowest and highest natural frequencies are found to be 4.965 and 6325 rad/sec, respectively. This spring-mass system is subjected to a sinusoidal acceleration $10\sin (0.5t)$ at its base. For simplicity, zero viscous damping is assumed in all computations.

Figure 10 Forced vibration responses to a spring-mass system.

Both ECFM and NFM with $\beta =1/4$ and $\gamma =1/2$ are adopted to calculate the responses by using a time step of $\mathrm{\Delta }t=0.05$ and 0.1 sec. Besides, the solution calculated from NFM with $\mathrm{\Delta }t=0.01\sec$ is treated as a reference solution for comparison. Numerical results are plotted in Figure 10. It is seen in Figure 10(a) that the result obtained from ECFM almost coincides together with that calculated from NFM for using $\mathrm{\Delta }t=0.05$ sec. A very similar phenomenon is also found in Figure 7(b) for $\mathrm{\Delta }t=0.1\sec$ although these results have more significant errors in contrast to the reference solution. Since the value of ${\mathrm{\Omega }}_0={\omega }_0^{(1000)}(\mathrm{\Delta }t)$ is as large as 632.5 as $\mathrm{\Delta }t=0.1$ sec, an unconditional stability of ECFM is implied. Besides, it is affirmed that the low frequency modes can be accurately integrated while there is no excessive amplitude growth or even instability for high frequency modes. As a consequence, a reliable solution can be yielded.

The CPU time consumed for each nonlinear dynamic analysis is also recorded. In fact, NFM consumes a CPU time of 7590 and 4192 sec for using $\mathrm{\Delta }t=0.05$ and 0.1 sec, respectively. In contrast, ECFM only consume a CPU time of 26 and 13, correspondingly. Hence, many computational efforts can be saved for ECFM since it involves no nonlinear iterations for each step. Meanwhile, an iteration procedure is required for NFM for each step. In fact, an average iteration number is about 2.8 for using $\mathrm{\Delta }t=0.05$ sec while for $\mathrm{\Delta }t=0.1\sec$ it is about 3.2.

7 Conclusions

Although structure-dependent integration methods have been successfully developed for time integration in the literature, it is still undisclosed for the derivation details and there does not exist a solid fundamental base to prove or support its feasibility for time integration. In this work, a typical procedure to develop structure-dependent integration methods is constructed from an eigen-based theory and then it seems that it can provide a fundamental basis for validating the feasibility of structure-dependent integration methods for conducting dynamic analysis of both linear elastic and nonlinear systems. For verification purpose, the first family of non-dissipative structure-dependent integration methods is developed herein in details for illustrating the developing procedure and its feasibility for time integration can be evident from the eigen-based theory. Because the numerical properties of this family structure-dependent integration methods have been explored for both linear elastic and nonlinear systems and its performance have been explored, they will not be elaborated again.

Acknowledgements

The author is grateful to acknowledge that this study is financially supported by the Ministry of Science and Technology, Taiwan, R.O.C., under Grant No. MOST-110-2221-E-027-020.

Conflict of Interest

The author declares that he has no conflict of interest.

References

[1] Chang, S.Y. (2002), “Explicit pseudodynamic algorithm with unconditional stability”, Journal of Engineering Mechanics, ASCE, 128 (9), 935–947. https://doi.org/10.1061/(ASCE)0733-9399(2002)128:9(935).

[2] Chang, S.Y. (2009), “An explicit method with improved stability property”, International Journal for Numerical Method in Engineering, 77 (8), 1100–1120. https://doi.org/10.1002/nme.2452.

[3] Chang, S.Y. (2010), “A new family of explicit method for linear structural dynamics”, Computers & Structures, 88 (11–12), 755–772. https://doi.org/10.1016/j.compstruc.2010.03.002.

[4] Chang, S.Y. (2015), “Dissipative, non-iterative integration algorithms with unconditional stability for mildly nonlinear structural dynamics”, Nonlinear Dynamics, 79 (2), 1625–1649. https://doi.org/10.1007/s11071-014-1765-7.

[5] Chang, S.Y. (2019) “A dual family of dissipative structure-dependent integration methods”, Nonlinear Dynamics, 98(1), 703–734. https://doi.org/10.1007/s11071-019-05223-y.

[6] Chang, S.Y. (2020) “Non-iterative methods for dynamic analysis of nonlinear velocity-dependent problems.” Nonlinear Dynamics, 101, 1473–1500. https://doi.org/10.1007/s11071-020-05836-8.

[7] Newmark, N.M. (1959), “A method of computation for structural dynamics”, Journal of Engineering Mechanics Division, 85 (3), 67–94.

[8] Hilber, H.M., Hughes, T.J.R., Taylor, R.L. (1977), “Improved numerical dissipation for time integration algorithms in structural dynamics”, Earthquake Engineering and Structural Dynamics, 5, 283–292. https://doi.org/10.1002/eqe.4290050306.

[9] Wood, W.L., Bossak, M., Zienkiewicz, O.C. An alpha modification of Newmark’s method. International Journal for Numerical Methods in Engineering 15 (1981) 1562–1566. https://doi.org/10.1002/nme.1620151011.

[10] Chung, J., Hulbert, G.M. A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-$\alpha$ method. Journal of Applied Mechanics 60(6) (1993) 371–375. https://doi.org/10.1115/1.2900803.

[11] Noels, L., Stainier, L., Ponthot, J.P., Bonini, J. Combined implicit-explicit algorithms for non-linear structural dynamics, European Journal of Computational Mechanics 11(5) (2002) 565–591. https://journals.riverpublishers.com/index.php/EJCM/article/view/2551.

[12] Papathanasiou, T.K. A Linearised $\theta$ Numerical Scheme for the Vibrations of Inextensible Beams. European Journal of Computational Mechanics 30(1) (2021) 121–144. https://doi.org/10.13052/ejcm1779-7179.3015.

[13] Zhou, X, Tamma, K.K. (2004), “Design, analysis, and synthesis of generalized single step single solve and optimal algorithms for structural dynamics.” International Journal for Numerical Methods in Engineering, 59, 597–668. https://doi.org/10.1002/nme.873.

[14] Krenk, S. (2006), “Energy conservation in Newmark based time integration algorithms”, Computer Methods in Applied Mechanics and Engineering, 195, 6110–6124. https://doi.org/10.1016/j.cma.2005.12.001.

[15] Tang, Y., Ren, D., Qin, H., Luo, C. (2021), “New family of explicit structure-dependent integration algorithms with controllable numerical dispersion”, Journal of Engineering Mechanics, 147(3), 04021001. http://10.1061/(ASCE)EM.1943-7889.0001901.

[16] Li, J., Yu, K., Li, X. (2018), “A generalized structure-dependent semi-explicit method for Structural dynamics”, ASME. Journal of Computational and Nonlinear Dynamics, 13(11), 111008. https://doi.org/10.1115/1.4041239.

[17] Chang, S.Y. An explicit structure-dependent algorithm for pseudo-dynamic testing, Engineering Structures, 46 (2013) 511–525. http://dx.doi.org/10.1016\%2Fj.engstruct.2012.08.009.

[18] Namadchi, A.H., Jandaghi, E., Alamatian, J. (2020) “A new model-dependent time integration scheme with effective numerical damping for dynamic analysis”, Engineering with Computers, 37(4) 2543–2558. https://doi.org/10.1007/s00366-020-00960-w.

[19] Qing, L.C., Zhong, J.L. (2018) “Unconditionally stable explicit displacement method for analyzing nonlinear structural dynamics problems”, Journal of Engineering Mechanics, 144(9), 04018079. https://doi.org/10.1061/(ASCE)EM.1943-7889.0001454.

[20] Dahlquist, G. Convergence and stability in the numerical integration of ordinary differential equations, Mathematica Scandinavica 4 (1956) 33–53. https://doi.org/10.7146/math.scand.a-10454.

[21] Clough, R.W., Penzien, J., Dynamics of Structures, 2nd Edition, McGraw-Hill, 1993. https://doi.org/10.12989/sem.2012.43.1.001.

[22] Belytschko, T. and Hughes, T.J.R. (1983), Computational Methods for Transient Analysis, Elsevier Science Publishers B.V., North-Holland, Amsterdam.

[23] Chang, S.Y. An unusual amplitude growth property and its remedy for structure-dependent integration methods, Computer Methods in Applied Mechanics and Engineering, 330 (2018) 498–521. https://doi.org/10.1016/j.cma.2017.11.012.

[24] Lax, P.D., Richtmyer, R.D. (1956), “Survey of the stability of linear finite difference equations”, Communications on Pure and Applied Mathematics, 9, 267–293. https://doi.org/10.1002/cpa.3160090206.

[25] Dahlquist, G. (1963), “A special stability problem for linear multistep methods”, BIT, 3, 27–43. https://doi.org/10.1007/BF01963532.

[26] Chang, S.Y. (2018) “Weak instability of CR explicit method for structural dynamics,” Journal of Computational and Nonlinear Dynamics, 13(5), 054501 https://doi.org/10.1115/1.4039379.

[27] Strang, G. (2020), Linear Algebra and Its Applications, 4th Editions, Cengage Learning, Inc.

[28] Oua, J.P., Long, X., Li, Q.S. (2007), “Seismic response analysis of structures with velocity-dependent dampers. Journal of Constructional Steel Research, 63, 628–638. https://doi.org/10.1016/j.jcsr.2006.06.034.

[29] Sadek, F., Mohraz, B., Riley, M.A. (2000), “Linear procedures for structures with velocity-dependent dampers”, Journal of Structural Engineering, 126(8), 887–895. https://doi.Org/10.1061/(ASCE)0733-9445(2000)126:8(887).

[30] Building Seismic Safety Council (BSSC). NEHRP guidelines for the seismic rehabilitation of buildings. FEMA 273, 1997, Washington, D.C.

Biography

Shuenn-Yih Chang received a Ph.D. degree from the University of Illinois, Urbana-Champaign in 1995. Prior to joining National Taipei University of Technology, Taipei, Taiwan, in 2002, he was an associate research fellow of National Centre for Research on Earthquake Engineering (NCREE). His research focuses on the structural dynamics and earthquake engineering. His work addresses on the dynamic testing of large-scale earthquake resistant structures, including the developments of novel pseudo-dynamic techniques. He is also of great interest in the design of accurate and efficient computational methods for dynamic problems of contemporary engineering interest. He has published more than one hundred and fifty refereed journal articles since 1988.

European Journal of Computational Mechanics, Vol. 31_5–6, 555–582.
doi: 10.13052/ejcm2642-2085.31561
© 2023 River Publishers

Abstract

1 Introduction

2 Concept of Eigen-based Theory

3 Development of CFM

3.1 Eigen-Decomposition of Equation of Motion

3.2 An Eigen-Dependent Integration Method

3.3 A Structure-Dependent Integration Method

4 Dual Implementations for CFM

4.1 An Explicit Implementation

4.2 An Implicit Implementation

5 Comparisons of Accuracy

6 Applications

6.1 A Non-classical Damping System

6.2 A Multi-story Building of 10 Stories

6.3 A Series of Large Spring-Mass Systems

7 Conclusions

Acknowledgements

Conflict of Interest

References

Biography