Research Article | Open Access
Changqing Li, Baoyi Sheng, Zhipeng Lai, Lizhong Jiang, Ping Xiang, "Two-Step Unconditionally Stable Noniterative Dissipative Displacement Method for Analysis of Nonlinear Structural Dynamics Problems", Shock and Vibration, vol. 2021, Article ID 4689090, 27 pages, 2021. https://doi.org/10.1155/2021/4689090
Two-Step Unconditionally Stable Noniterative Dissipative Displacement Method for Analysis of Nonlinear Structural Dynamics Problems
When solving structural dynamic problems, the displacement algorithm needs only calculating and storing structure’s displacements in the main calculation process, which makes the displacement algorithm have advantages over multivariable algorithms in calculation efficiency and storage requirements. By using a novel approach based on dimensional analysis firstly given by the first author, a one-parameter family of two-step unconditionally stable noniterative displacement algorithms, referred to as the CQ-2x method, is developed. Compared with other unconditionally stable noniterative multivariable algorithms such as the representative KR- method, the proposed method has advantages in several aspects. The CQ-2x method is unconditionally stable regardless of stiffness hardening or stiffness weakening, while the KR- method is only conditionally stable in case of stiffness hardening. The CQ-2x method needs only one solver within one time step, while the KR- method needs two solvers within one time step, which makes the CQ-2x method show higher efficiency. Numerical examples are presented to demonstrate the potential of the proposed method.
Direct integration methods are commonly used to solve structural dynamics equations. According to whether iteration is needed to solve nonlinear dynamic problems, these methods can be divided into two categories: implicit and noniterative. If the solution for the next time step depends on the responses in future time steps, then the algorithm is classified as an implicit algorithm [1, 2]. The well-known Houbolt , Newmark , Wislon , HHT- , WBZ , generalized- , and Bath [9–11] algorithms are examples of implicit algorithms. When the implicit methods are used to solve nonlinear structural dynamics equations, an iterative scheme (for example, the Newton–Raphson algorithm) is required in each time step which is often complicated and time consuming.
By contrast, equilibrium in noniterative algorithms is satisfied in the previous time step rather than the current step; consequently, such algorithms do not require the formation of the tangent stiffness matrix or any iterative calculations and are computationally more efficient than implicit algorithms in the nonlinear case when using same time step size.
The central difference method (CDM) [1, 2] is the most common and effective noniterative method since it not only needs no nonlinear iteration but also needs no solver for nonlinear problems. However, the CDM is only conditionally stable, which means extraordinary small time step size must be used in order to make the calculation stable. The conditionally stable CDM is preferred when the time step size required for accuracy is of the same order as the stability limit of the algorithm, e.g., in wave propagation problems [2, 12]. Nevertheless, unconditionally stable noniterative algorithms are more useful for solving nonlinear structural dynamics problems.
Chang’s algorithms [13–16], the CR  method, and the KR-  method are examples of unconditionally stable noniterative algorithms. The first unconditionally stable noniterative structure-dependent integration method was developed by Chang . An assessment of such structure-dependent unconditionally stable noniterative and semiexplicit (SE) algorithms for application to dynamic problems has been presented by Kolay and Ricles . These authors found that an SE algorithm may generate large damping forces when a realistic time step size is used for the nonlinear transient analysis of realistic structural systems subjected to seismic excitations, whereas the KR- method  may produce high-frequency overshoot. Furthermore, unconditionally stable noniterative model-dependent algorithms may become conditionally stable when severe stiffness hardening occurs, and most of them require the solution of linear equations whose coefficient matrix is nondiagonal twice within the one time step, which lowers their computational efficiency. When the number of DOFs of the structure is large, the majority of the computation time of an algorithm is consumed in solving linear equations whose coefficient matrices are nondiagonal. From the perspective of computational efficiency, it is preferable for an algorithm to solve a system of linear equations only once per time step.
In conjunction with the total energy framework and a generalized time-weighted residual approach, Har and Tamma [20, 21] have provided good guidance for developing new time integration algorithms and procedures for evaluating new and existing algorithms. However, the unconditionally stable noniterative structure-dependent integration methods can hardly be evaluated by those guidances. The algorithms including the structure-dependent algorithm may be easier to be designed by using a novel approach developed by the first author of the current study based on dimensional analysis. Most of the commonly used algorithms for structural dynamics analysis consist of one motion equation and two difference equations of displacement and velocity. By contrast, the novel dimensional analysis approach presented here relies on only a simple objective function or functions (the case of multiple functions will not be discussed in this paper), and the assumptions concerning displacement, velocity, and acceleration are implicitly included in their Taylor series expansions and the dynamic equilibrium equation.
When solving structural dynamic problems, according to the types of physical quantities that must be involved in the main calculation steps, the algorithm can be divided into the single variable method and multivariable method. As a single variable method, the displacement method needs only calculating and storing displacements in the calculation process, which makes the displacement method have advantages over other multivariable algorithms in calculation efficiency and storage requirements. So the novel approach was applied firstly to design the displacement method. By using this novel approach, a two-step unconditionally stable noniterative displacement algorithm (NDA) called the CQ-2x method (a two-step algorithm named after the first author) has been proposed . The CQ-2x method combines unconditional stability with noniterative formulations of both displacement and velocity, the use of one solver per time step, and no overshoot in either displacement or velocity. However, CQ-2x does not offer controllable numerical dissipation. To introduce controllable numerical damping, while simultaneously preserving the previously mentioned desirable numerical characteristics of the CQ-2x algorithm, a three-step displacement algorithm called the CQ-3 method  with controllable numerical dissipation has also been developed. However, the need for three steps is a disadvantage that may lower the efficiency of the algorithm. Thus, the development of a family of two-step unconditionally stable NDAs with controllable numerical damping is the main task of this paper.
The article begins with an objective function for a two-step NDA. Then, the preliminary form of the two-step NDA is obtained through dimensional analysis. Since there are 9 coefficients to be determined in the preliminary form of the NDA, equations for determining those coefficients from the algorithm’s consistency and stability are provided. Since the algorithm’s consistency and stability requirements cannot be used to determine all coefficients, the remaining free coefficients , , and are then determined from other characteristics of the algorithm, such as the relative period error and numerical damping ratio. A self-starting procedure is presented, and the overshoots in displacement and velocity are analyzed; subsequently, three numerical examples are provided to demonstrate the accuracy, stability, and efficiency of the proposed method.
2. Preliminary Form of NDA
For a single-DOF system, the linear structural dynamics equation can be written as
In the equation above, represent the mass, damping, and stiffness, respectively; , , and represent the acceleration, velocity, and displacement of the mass, respectively; is the external force; and is the time.
When no nonlinear damping is considered, the nonlinear structural dynamics equation can be written as follows :where is the external force.
The nonlinearity of a structure may be complex. In some cases, the structural nonlinearity can be expressed in terms of nonlinear damping and nonlinear stiffness. In these cases, the nonlinear structural dynamics equation can be rewritten aswhere the nonlinear damping and the nonlinear stiffness are given by
In order to obtain the numerical solution at the discrete time point, according to the definition of the “noniterative” displacement algorithm given above, the objective of a two-step NDA is to find when , , , , , , and are given; this can be expressed by the objective function given in equation (5) where is regarded as a dependent variable and , , , , , , and are regarded as independent variables:where is the time step size and the subscripts denote discrete time points. Because the displacement algorithm is to be designed, there are no velocity and acceleration occurring in equation (5). A displacement algorithm designed based on dimensional analysis to solve the objective function given in equation (5) can be written as follows:where are undetermined coefficients. The detailed derivation of equation (6) based on dimensional analysis is presented in Appendix A. can be determined from the NDA’s consistency and stability requirements.
3. Consistency Analysis to Determine
The consistency of an NDA of the form given in equation (6) can be analyzed by calculating the local algorithmic error. To calculate the local algorithmic error between the exact displacement and the computed displacement , it is assumed that, in equation (6),
where a quantity with brackets “()” denotes an exact quantity and a quantity without brackets “()” but with a subscript denotes a computed quantity.
Taylor series expansions of , , , and can be written as follows:
where a superscripted dot denotes differentiation with respect to time. The nonlinear dynamic equilibrium of equation (3) in step is written as
Rewriting equation (12) yields
The statement that the NDA has a nominal nth order time accuracy of displacement means that the following holds:
Substituting into equation (15) yields the condition that must be met for the NDA to achieve a nominal second-order time accuracy of displacement:
Equation (17) ensures that the NDA will achieve a nominal second-order time accuracy of displacement.
It can be seen from the derivation of equation (20) that the NDA is derived from only one objective function, equation (5). The possible assumed relationships among the displacement, velocity, and acceleration are all implied in the Taylor series expansions (as seen in equations (9a)–(9d)) and the nonlinear dynamic equilibrium equation (as seen in equation (10)). Newmark’s assumptions (with either constant coefficients or structure-dependent coefficients) concerning displacement and velocity are not needed for the NDA family derived in this paper although they are needed for most other algorithms; thus, this derivation is more comprehensive and allows no promising algorithms to be missed when the objective function is determined. It should be noted here that the NDA accuracy analysis applies to the case of nonlinear structural dynamics since nonlinear equation (10) is used in the derivation instead of linear equation (1), which is also a point of difference with respect to other algorithms.
When calculation of the velocity is necessary, it can be achieved with acceptable accuracy by means of any applicable displacement difference approximation. For example, one can use, but is not limited to, the displacement difference approximation used in the central-eccentric difference method :
Equation (21) has a nominal time accuracy of velocity of at least first order, and the local algorithmic velocity error when this approximation is applied in the proposed CQ-2x method is
After the noniterative calculation of and , can be calculated from the nonlinear dynamic equilibrium (as given in equation (3)) at step :
Equation (23) has a nominal time accuracy of acceleration of at least first order, and when this approximation is applied in the CQ-2x method, the local algorithmic acceleration error can be written aswhere
The accuracy expressed above is called the nominal accuracy because the derivation of the displacement consistency assumes that the values of at time steps and are known exactly, which cannot be ensured, meaning that this assumption leads to additional errors. However, it can be concluded that an NDA described by equations (20), (21), and (23) has a global time accuracy of at least first order, as will be validated by Example 1.
4. Stability Analysis to Determine
An algorithm is convergent if and only if it is consistent and stable. The consistency of the proposed NDA family has been analyzed above. In this section, the stability of an NDA with the form given in equation (20) will be studied. It should again be noted here that the NDA consistency analysis applies to the solution of nonlinear structural dynamics problems of the form given in equation (3) because equation (3), not equation (1), is used throughout the entire derivation. For the same reason, the following NDA stability analysis also applies to the solution of nonlinear structural dynamics problems of the form given in equation (3). The NDA expressed in equation (20) can be rewritten as follows:where the amplification matrix is
A stability study similar to the one conducted by Hilber and Hughes  is performed (the detailed derivation can be seen in Appendix C). The stability conditions for a first-order time-accurate NDA can be written as follows:
Equations (32a)–(32c) give the stability conditions for the NDA that apply when solving a nonlinear structural dynamics problem of the form given in equation (3) with at least first-order time accuracy. Equation (20) shows that the CDM is only one specific kind of NDA with , , and . To check the correctness of the NDA stability conditions given in equations (32a)–(32c), one can solve for the domain of stability of the CDM simply by substituting , , and into equations (32a)–(32c). Then, equations (32a)–(32c) take the following form:
It can be seen that equations (33a)–(33c) yield the correct stability domain for the CDM, . Since the stiffness ratio (where indicates stiffness softening and indicates stiffness hardening) does not appear in the amplification matrix of CQ-2x but does appear in the amplification matrices of other algorithms (for example, Chang’s algorithms ), which causes them to become unstable when stiffness hardening occurs, it may be concluded that when used to solve the nonlinear structural dynamics equation given in equation (3), the CQ-2x method is unconditionally stable regardless of whether stiffness hardening or softening occurs as long as equations (32a)–(32c) are always satisfied when and , which will be validated by Example 2.
5. CQ-2x Method with Controllable Numerical Dissipation
Up to this point, three coefficients have remained undetermined: , , and . These coefficients can be determined from the algorithm’s characteristics of unconditional stability and acceptable numerical dispersion and dissipation.
5.1. Determination of
An intermediate variable is introduced:
Consequently, can be written as
Then, simplifying equation (34) leads to
An analysis shows that, for a damping structure, a larger leads to more satisfactory algorithmic characteristics; thus, is assigned its largest possible value of 0.5 according to equation (37):
It can be seen from the above analysis that there are only two coefficients, and , that are undetermined in the preliminary form of the CQ-2x method. These coefficients will be determined by introducing controllable numerical dissipation into CQ-2x.
5.2. Determination of and
The effects of and on the relative period error of the preliminary CQ-2x method have been analyzed. The results when , , and are plotted in Figure 2. It can be seen that, as and increase, the relative period error also increases. Thus, to achieve the best accuracy, and should be assigned their lowest possible values of 0 and , respectively. The analysis shows that when and , the preliminary CQ-2x method has a nominal third-order time accuracy of displacement and a global second-order time accuracy; this will be validated by Example 1. However, if numerical dissipation is needed to filter out spurious vibrations in high-vibration modes, neither nor is appropriate.
It can be proven that when , the numerical damping ratio is also , as shown in the following equation:
It can also be proven that when , the numerical damping ratio of very high vibration modes is also zero, as shown in the following equation:
To introduce an appropriate level of numerical dissipation, and must be greater than their lowest values of 0 and , respectively. The numerical damping ratio depends on the spectral radius of the amplification matrix. It may be preferable for the spectral radius of the amplification matrix of the preliminary CQ-2x method to decrease with increasing with no bifurcation. According to equations (D1a) and (D1b), no bifurcation occurs when , which requires the following:
Equation (43) is satisfied when
Since a lower value of results in a higher accuracy, as mentioned before, is assigned the lowest possible value according to equation (44):
Equation (46) is the final algorithmic form of displacement calculation in the CQ-2x method.
To determine , while simultaneously confirming the unconditional stability of the CQ-2x method, consider the following two extreme cases: and . In accordance with equations (31a) and (31b), (36), (45), and (46) and considering the limits, it can be shown that equation (C1) reduces to the following:
It can be calculated from equation (47b) that
5.3. The Proposed Algorithm
Now that all coefficients in the preliminary NDA form presented in equation (6) have been determined, a family of unconditionally stable noniterative algorithms, collectively referred to as the CQ-2x method with one-parameter , can be obtained by combining equations (21), (23), (46), and (48), which are rewritten aswhere
5.4. Stability for Systems with Nonlinear Stiffness
It has been demonstrated that the CQ-2x method is unconditionally stable when it is used to solve the nonlinear structural dynamics equation given in equation (3). Here, the stability of CQ-2x when it is used to solve the nonlinear structural dynamics equation given in equation (2) will be analyzed. When a general nonlinear structural dynamics problem expressed in the form of equation (2) is to be solved, equation (2) is often written in the incremental form, and the z transform  can be used to explore an algorithm’s stability for the cases of stiffness hardening and softening. The solution to equation (2) is also often written in the incremental form, as in the following equation:where is the tangent stiffness at time step .
The stability of an integration algorithm is generally analyzed using one of two approaches, namely, recurrence relations using an amplification matrix or discrete control theory. The unconditional stability of CQ-2x when solving nonlinear equation (3) has been analyzed through the use of an amplification matrix. In this section, discrete control theory will be used to analyze the stability of CQ-2x when solving nonlinear equation (2). The discrete z transform and its real translation property are defined as follows :
First, the stability analysis of equation (46) using the z transform is considered. By applying the z transform and its real translation property to equation (46), this equation can be solved to determine X(z) in terms of the system properties, the time step, and F(z). This solution leads to the discrete transfer function of the algorithm, which is defined as the ratio of the z transforms of the output and the input and can be expressed in the following general form:where
It can be seen that the denominator of the transfer function in equation (52) has the same form as the characteristic equation given in equation (C1) for the amplification matrix of equation (46), whose unconditional stability has been proven.
Since the displacement calculation is unconditionally stable, the velocity calculation in equation (21), which depends solely on the displacement, is also naturally unconditionally stable.
Now, the stability of equation (50) will be analyzed. By combining equations (46), (21), and (50), the transfer function defined as the ratio of the z transforms of the output and the input can be expressed in the following general form:where
Since , , , and do not affect the algorithm’s stability, for brevity, they are not listed in this paper.
Although the denominators of equations (52) and (54) differ by a multiplicative factor , their representative characteristic equations are the same, which means that the acceleration calculation using equation (50) is also unconditionally stable. An analysis shows that the ratio of to only affects the value of the numerator of equation (54). Thus, it can be concluded that CQ-2x is unconditionally stable when it is used to solve the structural nonlinear dynamics equation given in equation (2); this will be validated by Example 2.
5.5. Accuracy Analysis for the Linear System
The accuracy of an integration algorithm generally depends on its numerical dispersion and energy dissipation characteristics. For the purposes of an accuracy analysis, the solutions to a free vibration problem obtained using CQ-2x can be written in the form of equation (D3). The two principal roots (no spurious real root) can be used to analyze the numerical dispersion and energy dissipation characteristics. The first characteristic is generally expressed in terms of the relative period error, as defined in equation (D6). The latter one can be expressed in terms of the numerical damping ratio , as defined in equation (D4).
Figures 5 and 6 show the relative period errors (Pe) and the equivalent damping ratios , respectively, of the proposed CQ-2x method and the KR- method for various values of . As shown in the figures, both Pe and increase with decreasing , and produces the maximum Pe and .
Figure 5 shows that the numerical dispersion of CQ-2x is similar to that of the KR- method for the same value of when . Figure 6 shows that the KR- method’s energy dissipation is better than that of CQ-2x. However, at small values of , both Pe and are small for CQ-2x regardless of the value of , which means that the lower-mode responses in an MDOF system are negligibly affected by the numerical damping. By contrast, increases with increasing , indicating that the undesired higher modes can be effectively damped out. Note that, for , the CQ-2x method exhibits no numerical energy dissipation, and the Pe value becomes identical to that of the KR- method and Newmark’s CAA method.
5.6. Self-Starting and Overshoot
There are several possible methods of self-starting, which is a term that refers to a parameterless unified starting procedure . In this paper, the self-starting of CQ-2x relies on , as calculated using the CDM, as shown in the following equation:
When , equation (46) can be written as
With known initial values of , , and , the value of can be calculated using equation (58). Thus, a self-starting displacement calculation procedure for CQ-2x can be obtained.
When the velocity is to be calculated, equation (59) can be used, with a second-order time accuracy of velocity:
Equation (59) represents the self-starting velocity calculation in CQ-2x.
When , equation (58) becomes
For CQ-2x, as expressed in equation (46), the recursive relationship among , , and can be expressed as follows:It can be seen from equations (60) and (61) that CQ-2x has zero-order displacement overshoot. Similarly, it is easy to find that CQ-2x has zero-order velocity overshoot by considering equations (21) and (59). Thus, CQ-2x can be classified as a family of noniterative algorithms of the form .
5.7. Calculation Procedure for MDOF Systems
When applied in MDOF systems, the forms of the aforementioned equations have some changes. When applied in MDOF systems, appearing in a SDOF system should then be replaced by matrices ; appearing in a SDOF system should be replaced by vectors . The calculation procedure for MDOF systems is provided as follows:(1)Formulate the stiffness matrix , mass matrix , and damping matrix of the systems.(2)Select appropriate coefficient and time increment . Determine by solving equation (49d).(3)Obtain initial values . Calculate by(4)Calculate by equations (58) and (59), respectively. Then, calculate nonlinear stiffness and damping by equation (4). Subsequently, calculate by the following equation:(5)For each time step (i = 2, 4, 5, …, n), calculate and by equations (49a) and (49b), respectively. Then, calculate nonlinear stiffness and damping by equation (4). Subsequently, calculate by equation (49c).
The main calculation efforts of CQ-2x are concentrated in Step 5 discussed above. The main calculation efforts in Step 5 are concentrated in equation (49a) because the coefficient matrix () in equation (49a) is nondiagonal and solving the coupling linear equations needs times of multiplication or division. The computational costs consumed by equation (49b) are , and those consumed by equations (4) and (3) are since solving equation (4) only needs vector calculations and solving equation (3) only needs solving linear decoupling equations whose coefficient matrix is the diagonal matrix . Thus, it is true that there is only one solver of equation (49a) within one time step for CQ-2x. Furthermore, when no velocity and acceleration need to be calculated, the entire calculation of CQ-2x can be completed with only the involvement of displacement. In this case, the time consumption of CQ-2x will become less, and the storage requirements of CQ-2x will become considerably lower than other unconditionally stable noniterative algorithms.
6. Numerical Examples
To evaluate the overall behavior and confirm the convergence rate, controllable numerical dissipation, unconditional stability in case of stiffness hardening, and computational efficiency of the derived algorithms, four examples are given. The CAA method, as an example of an implicit algorithm, and the KR- method, as an example of a family of unconditionally stable structure-dependent noniterative algorithms, are considered for comparison.
6.1. Example 1
The purpose of this example is to validate the convergence rate of the CQ-2x method. Consider a single-DOF linear resonance problem with the initial conditions :
The exact solutions  to this problem are
The circular self-frequency was assigned a value of . Thus, the period is s. The integration interval was chosen to be s. In the calculations, the step size was varied from s to s. The plots in Figures 7–12 show the curves of the maximum relative absolute errors and average relative absolute errors in displacement, velocity, and acceleration for the duration of the simulation. These errors were computed as functions of (on a logarithmic scale) using the CQ-2x, CAA, and KR- methods. When , the CQ-2x method shows first-order rates of convergence in displacement, velocity, and acceleration, identical to the KR- method. When , CQ-2x shows second-order rates of convergence in displacement, velocity, and acceleration, identical to the CAA method.
6.2. Example 2
The purpose of this example is to validate the controllable numerical dissipation and unconditional stability in case of stiffness hardening of the CQ-2x method. In this example, the well-known unforced and undamped Duffing oscillator  is considered. The nonlinear dynamics equation iswhere and . The initial conditions are and .
To assess the time integration scheme, the percentage error in terms of the energy is introduced :where is the total energy at .
The period of this oscillator is equal to 0.15. The exact phase portrait is given in Figure 13. Also, shown in Figures 14–22 are the phase portraits of numerical results obtained using different integration schemes with the same time step size of over a time duration of 100 . The maximum values of the percentage errors of the nine considered algorithms over this time duration of 100 are summarized in Table 1 for various time step sizes. When , CQ-2x has an error of 6.06%, whereas the KR- method has an error of at least 99.63%. It can be seen from Figures 14–17 that is smaller, the numerical damping is bigger, and the undamped Duffing oscillator decays faster. To validate the unconditional stability of CQ-2x in the stiffness hardening case, a time step size of over a time duration of 100 is considered. Figures 23–31 show the phase portraits of the numerical results obtained using the nine integration schemes with the same time step size of over a time duration of 100 . It is meaningless to study an algorithm’s accuracy when , but it is useful to explore its stability. Figures 23–26 show that even when , CQ-2x is still stable. By contrast, the KR- method is robustly stable only when , as seen in Figure 30; when , the KR- method becomes evidently unstable, as seen in Figures 27–29. e CAA method also exhibits some instability, as seen in Figure 31 (for the CAA method, the iteration termination criterion is that the absolute value of the error between the displacements in the previous and current iterations is less than ).