An Unusual Damped Stability Property and its Remedy for an Integration Method

Shuenn-Yih Chang

Department of Civil Engineering, National Taipei University of Technology, NTUT Box 2653, No. 1, Section 3, Jungshiau East Road, Taipei 10608, Taiwan, Republic of China
E-mail: changsy@ntut.edu.tw

Received 16 December 2021; Accepted 17 March 2022; Publication 28 April 2022

Abstract

An unusual stability property is found for a structure-dependent integration method since it exhibits a different nonlinearity interval of unconditional stability for zero and nonzero damping. Although it is unconditionally stable for the systems of stiffness softening and invariant as well as most systems of stiffness hardening, an unstable solution that is unexpected is obtained as it is applied to solve damped stiffness hardening systems. It is found herein that a nonlinearity interval of unconditional stability for a structure-dependent method may be drastically shrunk for nonzero damping when compared to zero damping. In fact, it will become conditionally stable for any damped stiffness hardening systems. This might significantly restrict its applications. An effective scheme is proposed to surmount this difficulty by introducing a stability factor into the structure-dependent coefficients of the integration method. This factor can effectively amplify the nonlinearity intervals of unconditional stability for structure-dependent methods. A large stability factor will result in a large nonlinearity interval of unconditional stability. However, it also introduces more period distortion. Consequently, a stability factor must be appropriately selected for accurate integration. After choosing a proper stability factor, a structure-dependent method can be widely and easily applied to solve general structural dynamic problems.

Keywords: Damping, stiffness hardening, stability, accuracy, structure-dependent method.

1 Introduction

Structure-dependent methods are a different type of integration methods when compared to conventional integration methods, whose coefficients are scalar constants [1–9]. This type of integration methods is characterized by structure dependency since its coefficients of the two difference equations can be functions of the product of the step size and initial structural properties for defining the problem under analysis [10–16]. It has been verified that there exists no explicit method that can be unconditionally stable in the linear multistep methods [17]. This Dahlquist barrier is deeply rooted in people’s minds in the realm of computational methods. To resolve this barrier, the fundamental basis of the structure-dependent methods was disclosed recently [18], where an eigen-based theory or a concept of eigenmodes can provide a solid foundation for the development and feasibility of this type of integration methods. Matrix coefficients for structure-dependent methods instead of scalar coefficients for conventional methods is the key issue to combine unconditional stability and explicit formulation together. Hence, the Dahlquist barrier can be surmounted.

Since an equation of motion is a second order differential equation and it has a natural frequency that is always positive, its stability can be defined based on the product of the step size $\mathrm{\Delta }t$ and natural frequency $\omega$, i.e., $\mathrm{\Omega }=\omega (\mathrm{\Delta }t)$. As a consequence, an unconditional stability implies that there is no limitation on any positive value of $\mathrm{\Omega }$, i.e., $0<\mathrm{\Omega }\le \mathrm{\infty }$. Whereas, a conditional stability implies that the positive value of $\mathrm{\Omega }$ is constrained by an inequality of $0<\mathrm{\Omega }\le {\mathrm{\Omega }}^{(u)}$, where ${\mathrm{\Omega }}^{(u)}$ is known as an upper stability limit. Hence, a stable computation can be yielded only if $0<\mathrm{\Omega }\le {\mathrm{\Omega }}^{(u)}$ is satisfied. This definition of stability is different from that for a general first order differential equation, whose eigenvalue can be a complex number with a negative real part and then A-stability is defined. The stability property of structure-dependent methods is also affected by the change of natural frequency or stiffness for a constant mass. To monitor the stiffness change, an instantaneous degree of nonlinearity $\delta$ is introduced [11]. In fact, it is defined as the instantaneous stiffness $k$ at a specific time over the initial stiffness $k_0$. As a result, $\delta =k/k_0$ is defined. Based on this definition, $0<\delta <1$, $\delta =1$ and $1<\delta <\mathrm{\infty }$ represent that a structure experiences the case of stiffness softening, invariant and hardening at the specific time, respectively.

The first structure-dependent method was successfully developed for pseudo-dynamic tests in 2002 [10], which will be referred as SD1, and it was found that it can only have an unconditional stability for stiffness softening and invariant systems, i.e., $0<\delta \le 1$ while it will become conditionally stable for the systems of stiffness hardening, i.e., $1<\delta <\mathrm{\infty }$ [11]. As a consequence, its applications to solve structural dynamics or earthquake engineering problems will be largely limited or inconvenient since the type of nonlinear behaviors of the problems under analysis might not be known in advance. To circumvent this difficulty, an improved structure-dependent method has been proposed later [12, 13], which is referred as SD2. Although it has an unconditional stability property for most stiffness hardening systems in addition to stiffness softening and invariant systems for zero viscous damping, i.e., $0<\delta \le 2$, an unusual instability property was recently experienced in the solution of a damped stiffness hardening system. However, there exist no reports for this problem in the literature. Hence, it must be further explored what causes this numerical instability.

It will be shown that SD2 can have an unconditional stability for the systems of stiffness softening and invariant as well as most systems of stiffness hardening for the case of zero viscous damping. Whereas, for nonzero viscous damping, it can only have an unconditional stability for stiffness softening and invariant systems and it becomes conditionally stable for stiffness hardening systems. In addition to revealing a conditional stability for damped stiffness hardening structures, it is also of interest to devise a technique to surmount this adverse stability property. A stability factor will be introduced into the initial stiffness of the denominator of structure-dependent coefficients to enhance its stability property. This modification of structure-dependent coefficients not only improves stability but affects accuracy. Hence, it is required to thoroughly investigate its stability and accuracy property so that an appropriate stability factor can be chosen. Besides, some numerical examples will be adopted to substantiate the feasibility of this modified formulation of the structure-dependent coefficients for both SD1 and SD2.

2 Formulation

Both SD1 and SD2 are classified as structure-dependent methods and their formulation can be generally expressed as:

$$m{a}_{i+1}+c{v}_{i+1}+k{d}_{i+1}=f_{i+1}$$
$$d_{i+1}=d_i+{\beta }_1(\mathrm{\Delta }t)v_i+{\beta }_2{(\mathrm{\Delta }t)}^2{a}_i+p_{i+1}$$	(1)
$$v_{i+1}=v_i+\frac{1}{2}(\mathrm{\Delta }t)(a_i+a_{i+1}).$$

where $m$, $c$, $k$ and $f$ are the mass, viscous damping coefficient, stiffness, and external force, respectively; $d_i$, $v_i$ and $a_i$ are the nodal displacement, velocity, and acceleration at the i-th step. The structure-dependent coefficients ${\beta }_1$, ${\beta }_2$ and $p_{i+1}$ [19] for SD1 and SD2 are found to be:

$${\beta }_1=\frac{1}{D_1}\left(1+\xi {\mathrm{\Omega }}_0\right),{\beta }_2=\frac{1}{D_1}\frac{1}{2},$$
$$p_{i+1}=\frac{1}{D_1}\frac{1}{4m}{(\mathrm{\Delta }t)}^2(f_{i+1}-f_i)\mathit{\hspace{1em}}\mathrm{SD1}$$
$${\beta }_1=\frac{1}{D_2}(1+\xi {\mathrm{\Omega }}_0),{\beta }_2=\frac{1}{D_2}\frac{1}{2}(1+\xi {\mathrm{\Omega }}_0),$$
$$p_{i+1}=\frac{1}{D_2}\frac{1}{2m}{(\mathrm{\Delta }t)}^2(f_{i+1}-f_i)\mathit{\hspace{1em}}\mathrm{SD2}.$$	(2)

where $D_1=1+\xi {\mathrm{\Omega }}_0+\frac{1}{4}{\mathrm{\Omega }}_0^2$; $D_2=1+\xi {\mathrm{\Omega }}_0+\frac{1}{2}{\mathrm{\Omega }}_0^2$; $\xi$ is often known as a viscous damping ratio and ${\mathrm{\Omega }}_0={\omega }_0(\mathrm{\Delta }t)$; ${\omega }_0=\sqrt{k_0/m}$ is the initial natural frequency and it is calculated from the initial stiffness of $k_0$. Since ${\beta }_1$, ${\beta }_2$ and $p_{i+1}$ are functions of initial structural properties and step size, SD1 and SD2 are classified as structure-dependent methods and are different from conventional methods.

In the subsequent study, the instantaneous degree of nonlinearity at the end of the i-th step is ${\delta }_i=k_i/k_0$, where $k_i$ is the stiffness at the end of the i-th step. Although the definition of ${\delta }_i$ is defined for a single degree of freedom system, it can be used to describe the corresponding modal stiffness change of each mode for a multiple degree of freedom system. There is a way to determine the instantaneous degree of nonlinearity for each mode at a specified time step by calculating the natural frequencies in correspondence to the updated structural properties. Consequently, after obtaining the natural frequency of the j-th mode at the i-th time step, the instantaneous degree of nonlinearity at the i-th step for this mode can be computed by:

$${\delta }_i^{(j)}={\left[\frac{{\omega }_i^{(j)}}{{\omega }_0^{(j)}}\right]}^2.$$

(3)

where ${\omega }_0^{(j)}$ and ${\omega }_i^{(j)}$ are the natural frequencies of the j-th mode corresponding to the initial stiffness and the stiffness at the end of the i-th step.

It has been verified that SD1 and SD2 can have unconditional stability in the interval of $0<{\delta }_i\le 1$ and $0<{\delta }_i\le 2$ for zero viscous damping, respectively. In fact, the stability property of SD1 for zero viscous damping has been derived in the reference [11] and is:

$$\begin{array}{cc}0<{\mathrm{\Omega }}_0\le \mathrm{\infty }\mathit{\hspace{1em}}\hfill & \mathrm{if}{\delta }_i\le 1\hfill \\ 0<{\mathrm{\Omega }}_0\le {\mathrm{\Omega }}_0^{(u)}=\frac{2}{\sqrt{{\delta }_i-1}}\mathit{\hspace{1em}}\hfill & \mathrm{if}{\delta }_i>1\hfill \end{array}.$$

(4)

while that for SD2 has been derived in the reference [12] and is found to be:

$$\begin{array}{cc}0<{\mathrm{\Omega }}_0\le \mathrm{\infty }\mathit{\hspace{1em}}\hfill & \mathrm{if}{\delta }_i\le 2\hfill \\ 0<{\mathrm{\Omega }}_0\le {\mathrm{\Omega }}_0^{(u)}=\frac{2}{\sqrt{{\delta }_i-2}}\mathit{\hspace{1em}}\hfill & \mathrm{if}{\delta }_i>2\hfill \end{array}.$$

(5)

For brevity, the interval of ${\delta }_i$ will be referred as a nonlinearity interval since it can reveal an interval for the variation of the instantaneous degree of nonlinearity. As a result, it is revealed by Equation (4) that SD1 will become conditionally stable in the nonlinearity interval of ${\delta }_i>1$. Similarly, Equation (5) reveals that SD2 has a conditional stability in the nonlinearity interval of ${\delta }_i>2$. Since there is rare for a real structure that its stiffness becomes twice of the initial stiffness (${\delta }_i>2$), a nonlinearity interval of unconditional stability of $0<{\delta }_i\le 2$ seems large enough for practical applications for a structure-dependent method.

One may take it for granted that both SD1 and SD2 can inherit the stability property of zero viscous damping. In other words, it is anticipated that both integration methods can have the same nonlinearity intervals of stability as shown in Equations (4) and (5) for nonzero viscous damping. However, the following example can be used to reveal that the prediction is incorrect. In fact, only SD1 can have the same nonlinearity intervals of stability for both zero and nonzero damping cases; whereas different nonlinearity intervals of stability are found for SD2. This phenomenon is first found in this work and therefore its root cause is thoroughly explored next.

3 Illustrated Example

To demonstrate that the actual performance of SD1 and SD2 in the solution of a damped and an undamped stiffness hardening system, a 11-story regular building is considered and it can be mathematically modelled as a shear building. In general, there is no rotation of a horizontal section at each floor level for a shear building due to a very large rigidity of the beam and floor systems. As a result, the second order differential equation of motion for the shear building can be expressed as:

$\left[\begin{array}{ccccc}\hfill m_1\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill m_2\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill \mathrm{\cdots }\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill m_{10}\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill m_{11}\hfill \end{array}\right]\left\{\begin{array}{c}\hfill {\ddot{u}}_1\hfill \\ \hfill {\ddot{u}}_2\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\ddot{u}}_{10}\hfill \\ \hfill {\ddot{u}}_{11}\hfill \end{array}\right\}+\mathbf{\text{C}}\left\{\begin{array}{c}\hfill {\dot{u}}_1\hfill \\ \hfill {\dot{u}}_2\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\dot{u}}_{10}\hfill \\ \hfill {\dot{u}}_{11}\hfill \end{array}\right\}$
$\mathrm{\hspace{1em}\hspace{1em}}\mathit{\hspace{1em}}+\left[\begin{array}{ccccc}\hfill k_1+k_2\mathit{\hspace{1em}}\hfill & \hfill -k_2\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill -k_2\mathit{\hspace{1em}}\hfill & \hfill k_2+k_3\mathit{\hspace{1em}}\hfill & \hfill -k_3\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill -k_3\mathit{\hspace{1em}}\hfill & \hfill \mathrm{\cdots }\mathit{\hspace{1em}}\hfill & \hfill -k_{10}\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill -k_{10}\mathit{\hspace{1em}}\hfill & \hfill k_{10}+k_{11}\mathit{\hspace{1em}}\hfill & \hfill -k_{11}\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill -k_{11}\mathit{\hspace{1em}}\hfill & \hfill k_{11}\hfill \end{array}\right]\left\{\begin{array}{c}\hfill u_1\hfill \\ \hfill u_2\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill u_{10}\hfill \\ \hfill u_{11}\hfill \end{array}\right\}$
$\mathrm{\hspace{1em}\hspace{1em}}=-\left\{\begin{array}{c}\hfill m_1\hfill \\ \hfill m_2\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill m_{10}\hfill \\ \hfill m_{11}\hfill \end{array}\right\}a.$	(6)

where $u_1$ to $u_{11}$ are the lateral displacement at each floor from the bottom to top story. In addition, C is a damping coefficient matrix and it will be defined later. The lumped mass concentrated at each floor is specified as the same value of $m_i=2\times {10}^{4}kg$ for $i=1,2,\mathrm{\dots },11$; the input ground acceleration is assumed as $a=5\times {10}^2\sin (\pi t)m/{\sec }^2$; and a total stiffness for each story is given as:

$k_1$	$={10}^8(1+1.2\sqrt{|u_1|})$
$k_i$	$={10}^8(1+1.2\sqrt{|u_i-u_{i-1}|}),\mathrm{for}i=2,3,\mathrm{\dots },11.$	(7)

It is evident that the story stiffness will become greater than the initial stiffness after the building deforms. Hence, a stiffness hardening system can be mimicked. In Equation (3), the damping coefficient matrix is assumed to be a linear combination of the mass and stiffness matrices, i.e., $\mathbf{\text{C}}=a_0\mathbf{\text{M}}+a_1\mathbf{\text{K}}$, which is known as a Rayleigh damping. Besides, the same viscous damping ratio of $\xi$ is assumed for the first two modes. As a result, the coefficients $a_0$ and $a_1$ can be determined from the following equation [20]:

$$a_0=\frac{2{\omega }_0^{(1)}{\omega }_0^{(2)}}{{\omega }_0^{(1)}+{\omega }_0^{(2)}}\xi ,a_1=\frac{2}{{\omega }_0^{(1)}+{\omega }_0^{(2)}}\xi .$$

(8)

The initial natural frequencies of the 11-story building are found to be ${\omega }_0^{(1)}=9.65$ and ${\omega }_0^{(2)}=28.77$ $rad/\sec$ in correspondence to the first and second modes. In the subsequent calculations, $\xi =0.1$ is taken and it will lead to $a_0=1.45\times {10}^{-1}$ and $a_1=5.21\times {10}^{-3}$.

Equation (3) with zero and nonzero viscous damping is solved by using SD1 and SD2. It is worth noting that structure-dependent integration methods generally can have an explicit implementation and hence they are non-iterative for each step in conducting step-by-step integration. The implementation details for a multiple degree of freedom system for SD1 were presented in [11] while those for SD2 were shown in [12]. In addition, the well-known constant Average Acceleration Method (AAM) is also used to carry out computations for comparison. The result obtained from AAM with $\mathrm{\Delta }t=0.001$ sec is treated as a reference solution for both the damped and undamped cases. At first, AAM, SD1 and SD2 are applied to solve Equation (3) with $\mathrm{\Delta }t=0.02$ sec for zero viscous damping. The displacement responses for the 1${}^{st}$, 6${}^{th}$ and 11${}^{th}$ stories and the corresponding time histories of the instantaneous degree of nonlinearity are plotted in Figure 1. It is manifested from Figures 1(b), 1(d) and 1(f) that the shear building experiences stiffness hardening since $1\le {\delta }_i^{(j)}\le 2$ is generally found for each mode, i.e., for $j=1,2,\mathrm{\dots },11$. This is consistent with the story stiffness defined in Equation (3), where each story stiffness will increase after the building deforms. Besides, very significant fluctuations of ${\delta }_i^{(j)}$ are found in Figures 1(b), 1(d) and 1(f). It is revealed by Figures 1(a), 1(c) and 1(e) that the result obtained from SD2 almost coincide together with that calculated from AAM. Evidently, these results are reliable when compared to the reference solution.

Figure 1 Forced vibration responses to 3-story building with zero viscous damping.

Unlike the solutions obtained from AAM and SD2, an unstable result is obtained from SD1 and it becomes unstable very early. This is because SD1 can only have a conditional stability for the systems of stiffness hardening and the violation of the upper stability limit might be responsible for the instability. In fact, a maximum value of ${\delta }_i^{(11)}=1.85$ is found in Figure 1(f) for the 11th mode, whose initial natural frequency is found to be ${\omega }_0^{(11)}=140.10$ rad/sec. Hence, ${\mathrm{\Omega }}_0^{(11)}={\omega }_0^{(11)}\times \mathrm{\Delta }t=2.80$ is found and is greater than the upper stability limit ${\mathrm{\Omega }}_0^{(u)}=2.17$, which can be calculated from Equation (4) for the case of zero viscous damping. As a result, the violation of the upper stability limit of ${\mathrm{\Omega }}_0^{(11)}=2.80>{\mathrm{\Omega }}_0^{(u)}=2.17$ is the root cause of numerical explosions in the results calculated from SD1.

Numerical results for the case of nonzero viscous damping are plotted in Figure 2. In these calculations, it is found that an accurate solution can be still achieved by using AAM with a step size of $\mathrm{\Delta }t=0.06$ sec, which is much larger than that for using AAM with $\mathrm{\Delta }t=0.02$ sec to solve the zero viscous damping system. This might be due to the fact that the high frequency responses are suppressed or removed by the viscous damping. Thus, a large step size can be applied to yield an accurate solution. Again, it is also found that the variation of the instantaneous degree of nonlinearity for each mode and each step is in the nonlinearity interval of $1\le {\delta }_i^{(j)}\le 2$ as shown in Figures 2(b), 2(d) and 2(f). Both SD1 and SD2 exhibit numerical explosions as shown in Figures 2(a), 2(c) and 2(e). Apparently, the cause of instability for SD1 is the same as the zero viscous damping case. However, an unstable solution for SD2 seems to be totally inconsistent with the analytical prediction that it has a nonlinearity interval of unconditional stability of $0<{\delta }_i\le 2$ for zero viscous damping. This instability strongly indicates that SD2 might have a different stability property for nonzero viscous damping in contrast to zero viscous damping. Therefore, the root cause of this instability is thoroughly investigated next.

Figure 2 Forced vibration responses to 3-story building with nonzero viscous damping.

4 Stability Property

The illustrated example reveals that SD1 and SD2 might have different stability properties for damped and undamped stiffness hardening systems. Hence, the cause of this difference in stability property is explored next. At first, the variation of the upper stability limit with instantaneous degree of nonlinearity for each step, which will be defined later, will be numerically calculated and plotted for different viscous damping ratios for both SD1 and SD2. Subsequently, a stability factor will be introduced to improve the stability property of both methods [21]. Finally, the accuracy of SD1 and SD2 affected by this stability factor will be explored.

4.1 Damped Stability Property

Equation (2) can be used to solve a nonlinear system with nonzero damping of $\xi \ne 0$ and zero external force of $f_{i+1}=0$. The time integration of each step can be alternatively expressed as a recursive matrix form [8–16] and thus the solution procedure at the i-th step can be rewritten as:

$${\mathbf{\text{x}}}_i={\mathbf{\text{A}}}_i{\mathbf{\text{x}}}_{i-1},{\mathbf{\text{x}}}_i={[d_i,(\mathrm{\Delta }t)v_i,{(\mathrm{\Delta }t)}^2{a}_i]}^{\mathrm{T}}.$$

(9)

where ${\mathbf{\text{A}}}_i$ is an amplification matrix in correspondence to the i-th time step and it might vary for each step for a nonlinear system. The characteristic equation of this amplification matrix can be determined from solving an eigenvalue problem of $|{\mathbf{\text{A}}}_i-\lambda \mathbf{\text{I}}|=0$ and it can be applied to determine the numerical properties of the method at the i-th time step. As a result, it can be written as:

$$\lambda ({\lambda }^2-A_1\lambda +A_2)=0.$$

(10)

where

$A_1$	$=2-{\displaystyle \frac{1}{G}}2\xi {\mathrm{\Omega }}_i-{\displaystyle \frac{1}{G}}\left({\displaystyle \frac{1}{2}}{\beta }_1+{\beta }_2\right){\mathrm{\Omega }}_i^2,$
$A_2$	$=1-{\displaystyle \frac{1}{G}}2\xi {\mathrm{\Omega }}_i+{\displaystyle \frac{1}{G}}\left({\displaystyle \frac{1}{2}}{\beta }_1-{\beta }_2\right){\mathrm{\Omega }}_i^2.$	(11)

where $G=1+\xi {\mathrm{\Omega }}_i$ is defined.

It is manifested from Equation (9) that the recursive matrix form is generally different from that for a linear elastic system. The amplification matrix is a constant matrix for a linear elastic system while for a nonlinear system it may vary for each step. Hence, evaluations of an integration method are conducted for a specific time step but not for a complete step-by-step integration procedure. However, it is still indicative for a complete step-by-step integration procedure since a complete integration procedure consists of each step. Thus, the numerical properties for SD1 or SD2 at the i-th step can be determined from Equation (10).

It is evident that there is a zero eigenvalue of ${\lambda }_3=0$ in Equation (10). In addition, if it also has two principal complex conjugate eigenvalues of ${\lambda }_{1,2}$ in addition to $|{\lambda }_{1,2}|\le 1$, a bounded oscillatory solution can be yielded. Thus, the two principal roots of Equation (10) at the end of the i-th time step can be expressed in an exponential form as:

$${\lambda }_{1,2}=\frac{1}{2}{A}_1\pm j\sqrt{A_2-\frac{1}{4}{A}_1^2}=e^{-\overline{\xi }{\overline{\mathrm{\Omega }}}_i\pm j{\overline{\mathrm{\Omega }}}_i^D}.$$

(12)

where ${\overline{\mathrm{\Omega }}}_i={\overline{\omega }}_i(\mathrm{\Delta }t)$, ${\overline{\mathrm{\Omega }}}_i^D={\overline{\mathrm{\Omega }}}_i\sqrt{1-{\xi }^2}$ and $j=\sqrt{-1}$. Hence, a phase shift of ${\overline{\mathrm{\Omega }}}_i^D$ and a numerical damping ratio of ${\overline{\xi }}_i$ can be determined by:

$${\overline{\mathrm{\Omega }}}_i^D={\tan }^{-1}\sqrt{\frac{4A_2}{A_1^2}-1},{\overline{\xi }}_i=-\frac{\ln (A_2)}{2{\overline{\mathrm{\Omega }}}_i}.$$

(13)

This phase shift can be further applied to define a relative period error and is a measure of period distortion. In general, it is defined as:

$$E_p=\frac{{\overline{T}}_i-T_i}{T_i}=\frac{{\omega }_i}{{\overline{\omega }}_i}-1,{\overline{T}}_i=\frac{2\pi }{{\overline{\omega }}_i},T_i=\frac{2\pi }{{\omega }_i}.$$

(14)

where $T_i$ and ${\overline{T}}_i$ are used denote the true and calculated periods of the system, respectively. It is evident that ${\overline{\xi }}_i$ and ${\overline{\omega }}_i$ can be treated as the quantities corresponding to ${\xi }_i$ and ${\omega }_i$ in a numerical procedure.

The stability property of SD1 and SD2 has been studied for nonlinear systems [11–13]. However, only the case of zero viscous damping was explored for simplicity as shown in Equations (4) and (5). This is because that it is very complicated to analytically conduct a stability analysis for nonzero viscous damping. Hence, a numerical method is adopted as an alternative for conducting a stability analysis. For given $\xi$ and ${\delta }_i$, the three eigenvalues ${\lambda }_i$ of the characteristic equation as shown in Equation (10) can be calculated for a given ${\mathrm{\Omega }}_0$. Notice that ${\mathrm{\Omega }}_0$ varies from a small value, i.e., close to zero, to a large value for determining the upper stability limit ${\mathrm{\Omega }}_0^{(u)}$ in the subsequent calculations. The upper stability limit is the maximum value of ${\mathrm{\Omega }}_0$, where $|{\lambda }_i|$ for $i=1,2,3$ must be always less than or equal to 1.

Figure 3 shows the variation of upper stability limit of ${\mathrm{\Omega }}_0^{(u)}$ with ${\delta }_i$ for SD1 and SD2 for different values of $\xi =0$, 0.1 and 0.3. It is seen in Figure 3 that SD1 can have almost the same nonlinearity interval of unconditional stability for zero and nonzero viscous damping cases and the upper stability limit in the nonlinearity interval of conditional stability is only slightly affected by viscous damping. However, drastic differences are found for SD2. In fact, it has a nonlinearity interval of unconditional stability of $0<{\delta }_i\le 2$ for zero viscous damping whereas this interval will be shrunk to be $0<{\delta }_i\le 1$ for nonzero viscous damping. These calculated results of Figure 3 can be applied to thoroughly explain why the unstable solutions were found in Figure 2 for SD1 and SD2.

Figure 3 Variation of upper stability limit with ${\delta }_i$ for SD1 and SD2.

The maximum value of ${\delta }_i^{(11)}=1.83$ can be found from Figure 2(f). On the other hand, ${\mathrm{\Omega }}_0^{(11)}=4.20$ is found due to ${\omega }_0^{(11)}=140.10$ rad/sec as well as $\mathrm{\Delta }t=0.03$ sec. Hence, one can find an upper stability limit of ${\mathrm{\Omega }}_0^{(u)}=2.44$ from Figure 3 for the case of ${\delta }_i^{(11)}=1.83$ and $\xi =0.1$ for SD1. Similarly, it is found to be ${\mathrm{\Omega }}_0^{(u)}=1.67$ for SD2. Clearly, ${\mathrm{\Omega }}_0^{(11)}=4.20>{\mathrm{\Omega }}_0^{(u)}=2.44$ is found for SD1 and ${\mathrm{\Omega }}_0^{(11)}=4.20>{\mathrm{\Omega }}_0^{(u)}=1.67$ for SD2, the upper stability limit is violated for both SD1 and SD2. Thus, instability occurs.

4.2 Resolving Instability

It is beneficial for structure-dependent methods to resolve this difficulty since a conditional stability for any damped stiffness hardening systems will largely limit their applications or lead to a great inconvenience. A scheme can be used to magnify the nonlinearity interval of unconditional stability for both SD1 and SD2. Since the structure-dependent coefficients are functions of the initial stiffness $k_0$ and their stability properties are closely related to ${\delta }_i$ (or $k_0$), there is a great motive to enlarge the nonlinearity interval of unconditional stability for these integration methods. Based on the definition of ${\delta }_i$, a nonlinearity interval of $0<{\delta }_i\le 1$ implies $k\le k_0$, which means that an unconditional stability can be yielded if the instantaneous stiffness of $k$ is no more than the initial stiffness of $k_0$. Hence, there exists a great idea to virtually amplify the initial stiffness $k_0$ by a stability factor of $\sigma$. As a result, one can have $k\le \sigma k_0$, which implies a nonlinearity of $0<{\delta }_i\le \sigma$. This implies that a nonlinearity interval of unconditional stability will be effectively amplified from $0<{\delta }_i\le 1$ to $0<{\delta }_i\le \sigma$ if $\sigma$ is chosen to be greater than 1.

It is apparent that the slight modification from $k_0$ to $\sigma k_0$ will also alter the formulation of the structure-dependent methods. Thus, the denominator and loading-correction terms must be modified to be:

$D_1$	$=1+\xi {\mathrm{\Omega }}_0+{\displaystyle \frac{1}{4}}\sigma {\mathrm{\Omega }}_0^2,p_{i+1}={\displaystyle \frac{1}{D_1}}{\displaystyle \frac{1}{4}}\sigma {(\mathrm{\Delta }t)}^2(f_{i+1}-f_i)\mathit{\hspace{1em}}\mathrm{SD1}$
$D_2$	$=1+\xi {\mathrm{\Omega }}_0+{\displaystyle \frac{1}{2}}\sigma {\mathrm{\Omega }}_0^2,p_{i+1}={\displaystyle \frac{1}{D_2}}{\displaystyle \frac{1}{2}}\sigma {(\mathrm{\Delta }t)}^2(f_{i+1}-f_i)\mathit{\hspace{1em}}\mathrm{SD2}.$	(15)

The only difference between SD1 and SD2 is the coefficient between $\frac{1}{4}\sigma$ and $\frac{1}{2}\sigma$, which is consistent with the coefficients as shown in Equation (2). This modification is to modify the initial stiffness $k_0$, which is adopted to determine the structure-dependent coefficients of ${\beta }_1$, ${\beta }_2$ and $p_{i+1}$, to $\sigma k_0$. Consequently, all the terms related to the initial stiffness of $k_0$ or its equivalent form of ${\mathrm{\Omega }}_0^2={\omega }_0^2{(\mathrm{\Delta }t)}^2=(k_0/m){(\mathrm{\Delta }t)}^2$ must be replaced by $\sigma k_0$. Consequently, Equation (5.1) can be obtained. On the other hand, it can be also derived from the procedure for developing a structure-dependent integration method that is constructed in the reference [18].

To affirm the effectiveness of the stability factor for improving the stability property of SD1 and SD2, the variations of ${\mathrm{\Omega }}_0^{(u)}$ with ${\delta }_i$ for different values of $\sigma =1$, 2 and 3 as well as different values of $\xi =0$, 0.1 and 0.3 are plotted in Figure 4. The procedure to calculate the results for Figure 3 is also applied to obtain the calculated results for Figure 4. Figure 4(a) reveals that a nonlinearity interval of unconditional stability can be amplified from $0<{\delta }_i\le 1$ to $0<{\delta }_i\le \sigma$ for SD1 after introducing a stability factor $\sigma$ into the structure-dependent coefficients for both zero and nonzero viscous damping. On the other hand, for SD2, a nonlinearity interval of unconditional stability changes from $0<{\delta }_i\le 2$ to $0<{\delta }_i\le 2\sigma$ for zero viscous damping while for the case of nonzero viscous damping it changes from $0<{\delta }_i\le 1$ to $0<{\delta }_i\le \sigma$. Consequently, it can be concluded that the choice of the stability factor $\sigma >1$ can effectively magnify a nonlinearity interval of unconditional stability for both SD1 and SD2.

Figure 4 Variation of upper stability limit with ${\delta }_i$ for SD1 and SD2 with different $\sigma$.

4.3 Accuracy Affected by Stability Factor

Stability and accuracy are the most important properties of numerical methods [11–13]. In addition to the investigation of stability property, it is also of interest to explore whether period distortion is significantly affected by stability factor. The relative period error for a specified value of $\mathrm{\Delta }t/T_0$, where $T_0=2\pi /{\omega }_0$, for SD1 and SD2 can be calculated from Equation (14). As a consequence, the variation of relative period error with $\mathrm{\Delta }t/T_0$ for both SD1 and SD2 for different values of $\sigma =1$, 2 and 3 is shown in Figure 5. In general, SD2 will have larger period distortion than for SD1 for a given value of $\sigma =1$. In fact, the relative period error is found to be as large as 2.01%, 4.37% and 6.67% for SD2 in correspondence to $\sigma =1$, 2 and 3 as $\mathrm{\Delta }t/T_0=0.05$. Similarly, it is found to be 0.81%, 2.01% and 3.19% for SD1 correspondingly. It is evident that a large stability factor will result in a large period distortion.

To choose an appropriate stability factor for SD1 or SD2 for practical applications, some influencing factors must be considered. It is generally required to estimate the possible variation of the nonlinearity interval of ${\delta }_i$. In general, for a realistic structure, its stiffness is very rare to become twice of the initial stiffness and then the nonlinearity interval of unconditional stability $0<{\delta }_i\le 2$ is of general interest. This implies that a nonlinearity interval of unconditional stability of $0<{\delta }_i\le 2$ is required for an integration method so that there will be no constraint on step size. It is evident that both SD1 and SD2 with $\sigma =2$ and 3 can automatically satisfy this requirement. However, it is seen in Figure 5 that SD1 with $\sigma =3$ and SD2 with $\sigma =2$ and 3 will cause too much period distortion. As a consequence, it seems that the SD1 with $\sigma =2$ is the best one among these methods for practical applications. This integration method will be referred as MSD1 for brevity in the subsequent study.

Figure 5 Variation of relative period error with $\mathrm{\Delta }t/T_0$ for SD1 and SD2 with different $\sigma$.

5 Numerical Confirmation

After the analytical studies of stability and accuracy properties of SD1 and SD2 as well as their modified formulations, it is discovered that MSD1 can possess a nonlinearity interval of unconditional stability of $0<{\delta }_i\le 2$ for both zero and nonzero viscous damping as well as an acceptable period distortion. Therefore, it is of importance to further substantiate that MSD1 can be applied to solve a damped stiffness hardening system without experiencing any stability problem. For this purpose, two examples will be examined next.

5.1 A Mechanical System

To simulate mechanical systems that move along a straight line, the three basic elements of mass, spring and dashpot or damper are often adopted. A mechanical system, which is simulated by 3 lumped masses connected by 3 springs and the connection details can be found in Figure 6. The dynamic behaviors of the mechanical system are governed by a set of equations of motion as:

$\left[\begin{array}{ccc}\hfill m_1\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill m_2\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill 0\mathit{\hspace{1em}}\hfill & \hfill m_3\hfill \end{array}\right]\left\{\begin{array}{c}\hfill {\ddot{u}}_1\hfill \\ \hfill {\ddot{u}}_2\hfill \\ \hfill {\ddot{u}}_3\hfill \end{array}\right\}+\mathbf{\text{C}}\left\{\begin{array}{c}\hfill {\dot{u}}_1\hfill \\ \hfill {\dot{u}}_2\hfill \\ \hfill {\dot{u}}_3\hfill \end{array}\right\}$
$\mathrm{\hspace{1em}\hspace{1em}}+\left[\begin{array}{ccc}\hfill k_1+k_2\mathit{\hspace{1em}}\hfill & \hfill -k_2\mathit{\hspace{1em}}\hfill & \hfill 0\hfill \\ \hfill -k_2\mathit{\hspace{1em}}\hfill & \hfill k_2+k_3\mathit{\hspace{1em}}\hfill & \hfill -k_3\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill -k_3\mathit{\hspace{1em}}\hfill & \hfill k_3\hfill \end{array}\right]\left\{\begin{array}{c}\hfill u_1\hfill \\ \hfill u_2\hfill \\ \hfill u_3\hfill \end{array}\right\}$
$\mathrm{\hspace{1em}\hspace{1em}}=-\left\{\begin{array}{c}\hfill m_1\hfill \\ \hfill m_2\hfill \\ \hfill m_3\hfill \end{array}\right\}a.$	(16)

where

$$k_1={10}^7(1+\sqrt{|u_1|}),k_2={10}^7(1+\sqrt{|u_2-u_1|}),$$
$$k_3={10}^7(1+\sqrt{|u_3-u_2|}).$$	(17)

The lumped mass as shown in Figure 6 is specified as $m_1=m_3$ $={10}^{2}kg$ and $m_2={10}^{4}kg$. In addition, the input ground acceleration of $a={10}^2\sin (\pi t)m/{\sec }^2$ is adopted for the analysis. In Equation (5.1), $u_1$ to $u_3$ are introduced to represent the lateral displacement of each lumped mass as shown in Figure 6. The damping matrix of $C$ is also assumed to be a Rayleigh damping of $C=a_{0}M+a_{1}K$, where the scalar constants of $a_0$ and $a_1$ can be also calculated from Equation (8). The initial natural frequencies of this mechanical system are found to be ${\omega }_0^{(1)}=7.03$, ${\omega }_0^{(2)}=100.50$ and ${\omega }_0^{(3)}=141.60$ $rad/\sec$. In the subsequent calculations, the damping ratio of $\xi =0.1$ is taken and it leads to $a_0=1.31$ as well as $a_1=1.86\times {10}^{-3}$.

Figure 6 A mechanical system simulated by 3 masses and 3 springs.

Equation (7.1) is solved by MSD1 with $\sigma =2$ and $\mathrm{\Delta }t=0.02$ sec for zero viscous damping. In addition, the case of nonzero viscous damping is also solved except that a step size of $\mathrm{\Delta }t=0.05$ sec is taken. Numerical results of $u_3$ as well as the time histories of ${\delta }_i^{(j)}$ for each mode are shown in Figure 7. It is revealed by Figures 7(a) and 7(c) that MSD1 can have stable solutions for both zero and nonzero viscous damping. It seems that the results calculated from MSD1 almost coincide together with those calculated from AAM. Consequently, it is affirmed that MSD1 can have the same performance as that of AAM in the solution of both damped and undamped stiffness hardening systems. On the other hand, it is seen in Figures 7(b) and 7(d) that $1\le {\delta }_i^{(j)}\le 2$, for $j=1,2,3$, is found for both zero and nonzero viscous damping cases. The curve for ${\delta }_i^{(1)}$ is almost overlapped together with that of ${\delta }_i^{(3)}$ for both zero and nonzero viscous damping. On the other hand, $1\le {\delta }_i^{(2)}\le 1.15$ is generally found in Figures 7(b) and 7(d). Consequently, it is affirmed by this example that MSD1 can generally have a nonlinearity interval of unconditional stability of $0<{\delta }_i\le 2$ for zero and nonzero viscous damping in the solution of stiffness hardening systems.

Figure 7 Forced vibration response of the mechanical system.

5.2 A 11-Story Shear Building

In the illustrated example, it is shown that SD1 and SD2 result in numerical instability in the solution of the damped stiffness hardening system. To affirm the improved stability for MSD1, it is required to examine whether MSD1 can give an accurate solution for solving this damped stiffness hardening problem. In addition, this example is also solved by MSD1 for nonzero viscous damping.

The displacement responses of the 1st, 6th and 11th stories are presented in Figure 8. It is evident from Figures 8(b), 8(d) and 8(f) that MSD1 can provide very accurate solutions for the damped stiffness hardening system. In fact, the calculated results almost coincide together with those calculated from AAM with same step size and are very close to the reference solutions. Hence, the performance of MSD1 is drastically different from that of SD1 and SD2 in the solution of damped stiffness hardening system as shown in Figure 2. On the other hand, Figures 8(a), 8(c) and 8(e) also reveal that MSD1 can still give accurate solutions for the undamped stiffness hardening systems. This can be applied to corroborate that the slight modification of the formulation of SD1 to yield the new formulation of MSD1 does not affect its performance in the solution of undamped stiffness hardening systems.

Figure 8 Forced vibration responses of 11-story building calculated from MSD1.

5.3 Summary

Notice that SD1 results in an instability in the solution of an undamped stiffness hardening system as shown in Figures 1(a), 1(c) and 1(e). Besides, both SD1 and SD2 lead to unstable solutions in the solution of a damped stiffness hardening system as shown in Figures 2(a), 2(c) and 2(e). Hence, their applications to solve general structural dynamic problems will be severely limited or inconvenient since the nonlinear type of the structural problems may not be known in prior. It seems that the modified algorithm of MSD1 can have an improved stability property and then difficulty experienced by SD1 and SD2 can be surmounted.

Evidently, the improved stability of MSD1 allows it easily achieve an accurate solution in the solution of damped and undamped stiffness hardening systems. This is thoroughly validated by the examples as shown in 7.1 and 7.2. An iteration procedure is required for AAM for each step in the solution of nonlinear systems while MSD1 is non-iterative due to a simultaneous combination of a nonlinearity interval of unconditional stability of $0<{\delta }_i\le 2$, explicit implementation and second order accuracy. Clearly, MSD1 can have a high computational efficiency for solving general structural dynamic problems when compared to conventional implicit methods. In the solution of the mechanical system as presented in 7.1, an average iteration number is about 3.10 for each step if using AAM is adopted to solve the undamped stiffness hardening system while for the damped stiffness hardening system it consumes about 3.47. A stop criterion of ${10}^{-6}$ is adopted for nonlinear iterations. On the other hand, for solving the responses of 11-story building, an average iteration number is about 7.24 and 8.71 for using AAM to solve the undamped and damped stiffness hardening systems. A stringent stop criterion of ${10}^{-10}$ is taken for this example during nonlinear iterations. Consequently, much more CPU demand is consumed for AAM in contrast to MSD1.

6 Conclusions

An unusual stability property is discovered for a structure-dependent method. In fact, it has different stability property for zero and nonzero damping. In general, SD2 can have a nonlinearity interval of unconditional stability of $0<{\delta }_i\le 2$ for zero damping while it will be shrunk to $0<{\delta }_i\le 1$ for nonzero damping. This implies that it can only have conditional stability for any damped stiffness hardening systems. Therefore, its applications are strictly constrained or inconvenient due to conditional stability. A scheme of using a stability factor to surmount this difficulty is proposed. This scheme is simple and effective to magnify a nonlinearity interval of unconditional stability. In fact, this factor $\sigma$ is applied to virtually enlarge the initial stiffness from $k_0$ to $\sigma k_0$ and consequently a nonlinearity interval of unconditional stability for SD2 can be amplified from $0<{\delta }_i\le 2$ to $0<{\delta }_i\le 2\sigma$ for zero damping and from $0<{\delta }_i\le 1$ to $0<{\delta }_i\le \sigma$ for nonzero damping. Similarly, it is shown that a nonlinearity interval of unconditional stability for SD1 can be enlarged from $0<{\delta }_i\le 1$ to $0<{\delta }_i\le \sigma$ for zero and nonzero damping. Although a large stability factor can result in a large nonlinearity interval of unconditional stability, it also leads to a large period distortion. Therefore, a stability factor must be appropriately selected for a structure-dependent method. In summary, MSD1, i.e., SD1 with $\sigma =2$ is strongly recommended for practical applications since it has a nonlinearity interval of unconditional stability of $0<{\delta }_i\le 2$ in addition to a slight period distortion. This is because that there is almost no actual structures, whose instantaneous degree of nonlinearity is greater than 2, i.e., ${\delta }_i>2$. It is important to improve the stability property for the structure-dependent method since it can combine the nonlinearity interval of unconditional stability of $0<{\delta }_i\le 2$, explicit formulation and second order accuracy. Therefore, it is very computationally efficient for solving general structural dynamic problems.

Acknowledgements

The author is grateful to acknowledge that this study is financially supported by Ministry of Science and Technology, Taiwan, ROC, under Grant No. MOST-109-2221-E-027-002.

Conflict of Interest

The author declares that he has no conflict of interest.

References

[1] Houbolt, J.C. A recurrence matrix solution for the dynamic response of elastic aircraft. Journal of the Aeronautical Sciences 17 (1950) 540–550. https://doi.org/10.2514/8.1722

[2] Newmark, N.M. A method of computation for structural dynamics. Journal of Engineering Mechanics Division, 85(3) (1959) 67–94. https://doi.org/10.1061/JMCEA3.0000098

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

[4] 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

[5] 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

[6] 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

[7] 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

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

[9] Hughes, T.J.R., The Finite Element Method, Prentice-Hall, Inc., Englewood Cliffs, N.J., U.S.A. (1987).

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

[11] Chang, S.Y. Improved explicit method for structural dynamics, Journal of Engineering Mechanics, ASCE 133(7) (2007) 748–760. https://doi.org/10.1061/(ASCE)0733-9399(2007)133:7(748)

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

[13] Chang, S.Y., Yang, Y.S., Hsu, C.W. A family of explicit algorithm for general pseudo-dynamic testing, Earthquake Engineering and Engineering Vibration, 10(1) (2011) 51–64. http://www.springerlink.com/content/1647t81504715251/

[14] 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

[15] Chang, S.Y. A dual family of dissipative structure-dependent integration methods, Nonlinear Dynamics, 98(1) (2019) 703–734.

[16] Chang, S.Y. Non-iterative methods for dynamic analysis of nonlinear velocity-dependent problems.” Nonlinear Dynamics, 101 (2020) 1473–1500.

[17] 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

[18] Chang, S.Y. An eigen-based theory for structure-dependent integration methods for nonlinear dynamic analysis, International Journal of Structural Stability and Dynamics 20(12) (2020) 2050130. https://doi.org/10.1142/S0219455420501308

[19] 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

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

[21] Chang, S.Y. A general technique to improve stability property for a structure-dependent integration method, International Journal for Numerical Methods in Engineering, 101(9) (2015) 653–669. https://doi.org/10.1002/nme.4806

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 Center 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_1, 41–64.
doi: 10.13052/ejcm2642-2085.3112
© 2022 River Publishers

Abstract

1 Introduction

2 Formulation

3 Illustrated Example

4 Stability Property

4.1 Damped Stability Property

4.2 Resolving Instability

4.3 Accuracy Affected by Stability Factor

5 Numerical Confirmation

5.1 A Mechanical System

5.2 A 11-Story Shear Building

5.3 Summary

6 Conclusions

Acknowledgements

Conflict of Interest

References

Biography