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./e CQ-2xmethod is unconditionally stable regardless of stiffness hardening or stiffness weakening, while the KR-α method is only conditionally stable in case of stiffness hardening. /e CQ-2xmethod needs only one solver within one time step, while the KR-αmethod needs two solvers within one time step, whichmakes the CQ-2xmethod show higher efficiency. Numerical examples are presented to demonstrate the potential of the proposed method.


Introduction
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]. e well-known Houbolt [3], Newmark [4], Wislon [5], HHT-α [6], WBZ α [7], generalized-α [8], and Bath [9][10][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. e 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. e 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][14][15][16], the CR [17] method, and the KR-α [18] method are examples of unconditionally stable noniterative algorithms. e first unconditionally stable noniterative structure-dependent integration method was developed by Chang [13]. An assessment of such structuredependent unconditionally stable noniterative and semiexplicit (SE) algorithms for application to dynamic problems has been presented by Kolay and Ricles [19]. ese 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 [18] 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. e 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 [22]. e 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 [23] 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. us, the development of a family of two-step unconditionally stable NDAs with controllable numerical damping is the main task of this paper. e article begins with an objective function for a twostep NDA. en, 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 b 5 , b 7 , and b 8 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.

Preliminary Form of NDA
For a single-DOF system, the linear structural dynamics equation can be written as (1) In the equation above, m, c, and k represent the mass, damping, and stiffness, respectively; a, v, and u represent the acceleration, velocity, and displacement of the mass, respectively; f is the external force; and t is the time.
When no nonlinear damping is considered, the nonlinear structural dynamics equation can be written as follows [18]: where f ext (t) is the external force. e 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 as where the nonlinear damping c(t) and the nonlinear stiffness k(t) 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 twostep NDA is to find u t+2Δt when u t+Δt , u t , f t+Δt , m, c t+Δt , 2 Shock and Vibration k t+Δt , and Δt are given; this can be expressed by the objective function given in equation (5) where u t+2Δt is regarded as a dependent variable and u t+Δt , u t , f t+Δt , m, c t+Δt , k t+Δt , and Δt are regarded as independent variables: where Δt 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 b i are undetermined coefficients. e detailed derivation of equation (6) based on dimensional analysis is presented in Appendix A. b i can be determined from the NDA's consistency and stability requirements.

Consistency Analysis to Determine b i
e 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 ε u 2+3Δt between the exact displacement u(t + 2Δt) and the computed displacement u t+2Δt , 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.
Substituting n � 2 into equation (15) yields the condition that must be met for the NDA to achieve a nominal secondorder time accuracy of displacement: Equation (16) is satisfied only when e 0 � e 1 � e 2 � 0, and the conditions for which can be obtained by simultaneously solving equations (B.2a), (B.2b), and (B.2c) (the detailed derivation can be seen in Appendix B): Equation (17) ensures that the NDA will achieve a nominal second-order time accuracy of displacement.
To confirm the nominal second-order time accuracy of displacement achieved with equation (17), we can substitute equation (17) into equation (12) and obtain the following equation: where e substitution of equation (17) into equation (6) yields an NDA with 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). e 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 [24]: Equation (21) has a nominal time accuracy of velocity of at least first order, and the local algorithmic velocity error 4 Shock and Vibration After the noniterative calculation of u t+2Δt and v t+2Δt , a t+2Δt can be calculated from the nonlinear dynamic equilibrium (as given in equation (3)) at step t + 2Δt: 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 ε a t+2Δt (� a(t + 2Δt) − a t+2Δt ) can be written as where e accuracy expressed above is called the nominal accuracy because the derivation of the displacement consistency assumes that the values of u at time steps t and t + Δt 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.

Stability Analysis to Determine b i
An algorithm is convergent if and only if it is consistent and stable. e 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). e NDA expressed in equation (20) can be rewritten as follows: where the amplification matrix A is A stability study similar to the one conducted by Hilber and Hughes [25] is performed (the detailed derivation can be seen in Appendix C). e stability conditions for a first-order time-accurate NDA can be written as follows: Shock and Vibration where ω t+Δt is the natural frequency and ξ t+Δt is the damping ratio; equations (30a) and (30b) at time t + Δt can be rewritten as follows: Substituting equations (31a) and (31b) into equations (29a) and (29b) yields the following: 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 b 5 � − 1, b 7 � 1/2, and b 8 � 0. 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 b 5 � − 1, b 7 � 1/2, and b 8 � 0 into equations (32a)-(32c). en, equations (32a)-(32c) take the following form: It can be seen that equations (33a)-(33c) yield the correct stability domain for the CDM, Ω 2 t+Δt ≤ 4. Since the stiffness ratio δ (where δ < 1 indicates stiffness softening and δ > 1 indicates stiffness hardening) does not appear in the amplification matrix A of CQ-2x but does appear in the amplification matrices of other algorithms (for example, Chang's algorithms [16]), 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 Ω ∈ [0, ∞) and ξ ∈ [0, ∞), which will be validated by Example 2.

CQ-2x Method with Controllable Numerical Dissipation
Up to this point, three coefficients have remained undetermined: b 5 , b 7 , and b 8 . ese coefficients can be determined from the algorithm's characteristics of unconditional stability and acceptable numerical dispersion and dissipation.

Determination of b 7 .
e unconditional stability of an NDA of the form given in equation (20) means that when Ω ∈ [0, +∞) and ξ ∈ [0, +∞), equations (32a)-(32c) are always satisfied, which requires that the following must hold: Shock and Vibration An intermediate variable x is introduced: Consequently, b 8 can be written as en, simplifying equation (34) leads to (37) An analysis shows that, for a damping structure, a larger b 7 leads to more satisfactory algorithmic characteristics; thus, b 7 is assigned its largest possible value of 0.5 according to equation (37): Substituting equations (36) and (38) into equation (20) yields the preliminary algorithmic form of the CQ-2x method: It can be seen from the above analysis that there are only two coefficients, b 5 and x, that are undetermined in the preliminary form of the CQ-2x method. ese coefficients will be determined by introducing controllable numerical dissipation into CQ-2x.

Determination of b 5 and x.
When ξ � 0, according to equation (36), equations (31a) and (31b) can be written as e effects of b 5 and x on the relative period error of the preliminary CQ-2x method have been analyzed. e results when Δt/T � 0.1, b 5 � − 0.5, − 0.3, 0, 0.3, 1, 2, and 4, and x ∈ [0, 4] are plotted in Figure 2. It can be seen that, as b 5 and x increase, the relative period error also increases. us, to achieve the best accuracy, x and b 5 should be assigned their lowest possible values of 0 and − 0.5, respectively. e analysis shows that when x � 0 and b 5 � − 0.5, 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 x � 0 nor b 5 � − 0.5 is appropriate.
It can be proven that when x � 0, the numerical damping ratio is also ξ � 0, as shown in the following equation: It can also be proven that when b 5 � − 0.5, 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, b 5 and x must be greater than their lowest values of 0 and − 0.5, respectively. e numerical damping ratio depends on the spectral radius of the amplification matrix. It may be preferable for the spectral radius of the amplification matrix A of the preliminary CQ-2x method to decrease with increasing Ω with no bifurcation. According to equations (D1a) and (D1b), no bifurcation occurs when A 2 ≥ A 2 1 , which requires the following:

Shock and Vibration
Equation (43) is satisfied when Since a lower value of b 5 results in a higher accuracy, as mentioned before, b 5 is assigned the lowest possible value according to equation (44): Substituting equation (45) Equation (46) is the final algorithmic form of displacement calculation in the CQ-2x method.
To determine x, while simultaneously confirming the unconditional stability of the CQ-2x method, consider the following two extreme cases: Ω ⟶ 0 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 Figures 3 and 4 show how the spectral radius of the CQ-2x method varies with ρ ∞ and Ω; the findings are fully consistent with the above analysis.

e 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

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 [17] can be used to explore an algorithm's stability for the cases of stiffness hardening and softening. e solution to equation (2) is also often written in the incremental form, as in the following equation: where k tan is the tangent stiffness at time step t + Δt [18].
When addressing structural dynamic nonlinearity in the form of equation (2), the CQ-2x method uses equations (46), (21), and (50). e stability of an integration algorithm is generally analyzed using one of two approaches, namely, recurrence relations using an amplification matrix or discrete control theory. e 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). e discrete z transform and its real translation property are defined as follows [18]: e stability of CQ-2x when solving nonlinear equation (2) can be analyzed by applying discrete control theory to equations (46), (21), and (50).
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).
is solution leads to the discrete transfer function of the algorithm, which is defined as the ratio of the z transforms of the output u i+1 and the input f i+1 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 a i+1 and the input f i+1 can be expressed in the following general form: where Since n 3 , n 2 , n 1 , and n 0 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 Δt, 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 k tan to k t+Δt only affects the value of the numerator of equation (54).
us, 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.

Accuracy Analysis for the Linear System.
e 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). e two principal roots (no spurious real root) can be used to analyze the numerical dispersion and energy dissipation characteristics. e first characteristic is generally expressed in terms of the relative period error, as defined in equation (D6). e 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 ρ ∞ � 0 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 ρ ∞ ≠ 1. 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 (ρ ∞ ∈ [0.8, 1]) 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 ρ ∞ � 1, the CQ-2x method exhibits no numerical energy dissipation, and the Pe value becomes identical to that of the KR-α (ρ ∞ � 1) method and Newmark's CAA method.

Self-Starting and Overshoot.
ere are several possible methods of self-starting, which is a term that refers to a parameterless unified starting procedure [26]. In this paper, the self-starting of CQ-2x relies on u − Δt , as calculated using the CDM, as shown in the following equation: When t � − Δt, equation (46) can be written as Substituting equation (56) into equation (57) yields a value for u Δt as follows: With known initial values of u 0 , v 0 , and a 0 , the value of u Δt can be calculated using equation (58).
us, a selfstarting 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.
For CQ-2x, as expressed in equation (46), the recursive relationship among u t+2Δt , u t+Δt , and u t can be expressed as follows: (61) 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). us, CQ-2x can be classified as a family of noniterative algorithms of the form

Calculation Procedure for MDOF Systems.
When applied in MDOF systems, the forms of the aforementioned equations have some changes. When applied in MDOF systems, m, c t+iΔt , and k t+iΔt appearing in a SDOF system should then be replaced by matrices M, C t+iΔt , and K t+iΔt ; u t+iΔt , v t+iΔt , a t+iΔt , and f t+iΔt appearing in a SDOF system should be replaced by vectors u t+iΔt , v t+iΔt , a t+iΔt , and f t+iΔt . e calculation procedure for MDOF systems is provided as follows: (1) Formulate the stiffness matrix K 0 , mass matrix M, and damping matrix C 0 of the systems.    (3) only needs solving linear decoupling equations whose coefficient matrix is the diagonal matrix M. us, 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.

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. e 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.

Example 1.
e 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 u(0) � v(0) � 1: (64) e exact solutions [27] to this problem are e circular self-frequency ω was assigned a value of 2π. us, the period is T � 1 s. e integration interval was chosen to be T Fin � 10T � 10 s. In the calculations, the step size Δt was varied from 10 − 1 s to 10 − 4 s. e 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. ese errors were computed as functions of Δt (on a logarithmic scale) using the CQ-2x, CAA, and KR-α methods. When ρ ∞ ≠ 1, the CQ-2x method shows first-order rates of convergence in displacement, velocity, and acceleration, identical to the KR-α method. When ρ ∞ � 1, CQ-2x shows second-order rates of convergence in displacement, velocity, and acceleration, identical to the CAA method.

Example 2.
e 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 [28] is considered. e nonlinear dynamics equation is where S 1 � 100 and S 2 − 10. e initial conditions are u 0 � 1.5 and v 0 � 0.

Shock and Vibration
Equation (66) represents a hardening elastic spring. is is an example of a conservative system that maintains a constant total energy E in its exact solutions. For equation (66), To assess the time integration scheme, the percentage error in terms of the energy is introduced [29]: where E 0 is the total energy at t � 0. Equation (66) has an exact solution expressed in terms of Jacobian elliptic functions [30]:  where e period T of this oscillator is equal to 0.15. e exact phase portrait is given in Figure 13. Also, shown in e maximum values of the percentage errors of the nine considered algorithms over this time duration of 100 T are summarized in Table 1 for various time step sizes. When Δt � (1/25)T, CQ-2x (ρ ∞ � 1) 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 Δt � (1/2)T over a time duration of 100 T 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 Δt � (1/2)T over a time duration of 100 T. It is meaningless to study an algorithm's accuracy when Δt � (1/2)T, but it is useful to explore its stability. Figures 23-26 show that even when Δt � (1/2)T, CQ-2x (ρ ∞ � 1, 0.8, 0.5, 0) is still stable. By contrast, the KRα method is robustly stable only when ρ ∞ � 0, as seen in Figure 30; when ρ ∞ � 1, 0.8, and 0.5, 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 10 − 6 ). e time consumptions of the different algorithms with different time steps over the time duration of 100 T are listed in Table 2. For the same time step size, the noniterative algorithms (CQ-2x and KR-α) are more efficient than the implicit algorithm (CAA). In addition, the CQ-2x method is more efficient than the KR-α method for the same time step size. For example, when Δt � (1/25)T, the time consumption of CQ-2x is at most 0.0015 seconds, whereas the KR-α method requires at least 0.0048 seconds.

Example 3.
e purpose of this example is to validate the efficiency of the CQ-2x method compared with the implicit method and KR − α method. In this example, the forced vibration response of a n-DOF (n � 200, 500, 1000, 2000, and 4000) spring-mass system was considered, as shown in Figure 32. e structural properties of this system are assumed to be m i � 100 kg and k i � 10 7 [1 − (u i − u i− 1 ) 2 ] N/m, in which i � 1, 2, . . . , n. is system is subjected to a ground acceleration of 10 sin(πt). e highest natural frequency is 632.4 rad/s [31]. According to the highest natural frequency  Table 3. It can be seen that the accuracy of the CQ-2x method is sufficient. e maximum error ratio here is defined by the function of max (abs(calculated value-exact value)))/max(abs(exact value)), where max � maximum and abs � absolute. Importantly, without the need for iteration calculations, CQ-2x algorithms have approximately less than half the computational time consumption compared to that consumed by the implicit algorithm CAA. With solving the inverse matrix only once for each time step, CQ-2x algorithms are more efficient than the noniterative KR − α method in which solving the inverse matrix twice for each time step is needed, as listed in Table 4.

Example 4.
e purpose of this example is to validate the efficiency of the CQ-2x method in case of Vehicle-Bridge Coupled Vibration where the contact nonlinearity exists. In this example, a simple beam of span length � 25 m is subjected to a moving spring mass [32] (as shown in Figure 35). e following data are adopted: Young's modulus E � 2.87 Gpa, ρ � 0. 10 − 5 s is obtained for the CDM, and the calculation result is selected to be the exact solution, while with unconditional stability, the time step duration of CQ2x(ρ ∞ � 1), KRα(ρ ∞ � 1), and CAA is only determined by the requirements of the calculation accuracy; here, Δt � 0.001s. e midpoint vertical deflection of the beam and deflection of the spring mass versus time are plotted in Figures 36 and 37. e maximum displacement error ratios (defined as that in Example 3) of the midpoint vertical deflection of the beam and the deflection of the spring mass are listed in Table 5. It can be seen that the accuracy of the CQ-2x method is sufficient. e efficiency of the CDM (Δt � 0.00002 s), CQ2x(ρ ∞ � 1), KR-α(ρ ∞ � 1), and CAA (Δt � 0.001 s) is compared. e average time consumption of the four algorithms in which every algorithm had been performed five times is listed in Table 6. In this example, limited by its conditional stability, the CDM has to be performed with very small Δt � 0.00002 s, which leads to its efficiency much lower than the other three algorithms. Without the need for iteration calculations, CQ-2x and KR-α algorithms show higher efficiency than that of the implicit algorithm CAA although the maximum number of iteration in each time step of CAA is only 3. e CQ-2x algorithm are more efficient than the noniterative KR − α algorithm since the computation amount in each time step of CQ-2x is much less than that of KR-α.

Conclusions
A novel idea based on dimensional analysis for designing noniterative displacement algorithms to be applied in the analysis of nonlinear structural dynamics problems has been introduced and proven to be effective. Originating from only an objective function, a one-parameter family of two-step unconditionally stable noniterative displacement algorithms with controllable numerical energy dissipation has been developed. Analysis results show that the proposed CQ-2x method has unconditional stability in the cases of both stiffness hardening and softening. With the advantages of noniterative formulations of both displacement and velocity, controllable numerical dissipation, the use of one solver per time step, and no overshoot in either displacement or velocity, CQ-2x is expected to serve as a competitive displacement method for solving nonlinear structural dynamics problems.