Abstract
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.
1. 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]. The well-known Houbolt [3], Newmark [4], Wislon [5], HHT- [6], WBZ [7], generalized- [8], 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 [17] method, and the KR- [18] method are examples of unconditionally stable noniterative algorithms. The first unconditionally stable noniterative structure-dependent integration method was developed by Chang [13]. 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 [19]. 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 [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. 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 [22]. 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 [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. 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 [18]: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.
Substituting equations (7a)–(7f) into equation (6) yields the following:
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
Substituting equations (9a), (9c), and (9d) into equation (10) yields
Substituting equations (9a), (9b), and (11) into equation (8) yieldswhere the expressions for , , , , , and are given in Appendix B.
Rewriting equation (12) yields
The statement that the NDA has a nominal nth order time accuracy of displacement means that the following holds:
The combination of equations (13) and (14) leads to
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 (16) is satisfied only when , 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
The 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). 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 [24]:
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 [25] 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:
where
By considering , , and , where is the natural frequency and is the damping ratio; equations (30a) and (30b) at time 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 , , 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 [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 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
The unconditional stability of an NDA of the form given in equation (20) means that when and , equations (32a)–(32c) are always satisfied, which requires that the following must hold:
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):
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, 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
When , according to equation (36), equations (31a) and (31b) can be written as
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):
Substituting equation (45) into equation (39) yields
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
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.
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 [17] 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 [18].
When addressing structural dynamic nonlinearity in the form of equation (2), the CQ-2x method uses equations (46), (21), and (50).
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 [18]:
The 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). 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 [26]. 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
Substituting equation (56) into equation (57) yields a value for as follows:
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 [20].
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 [27] 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 [28] is considered. The nonlinear dynamics equation iswhere and . The initial conditions are and .
Equation (66) represents a hardening elastic spring. This is an example of a conservative system that maintains a constant total energy 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 is the total energy at .
Equation (66) has an exact solution expressed in terms of Jacobian elliptic functions [30]:where
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 ).
The time consumptions of the different algorithms with different time steps over the time duration of 100 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 , the time consumption of CQ-2x is at most 0.0015 seconds, whereas the KR- method requires at least 0.0048 seconds.
6.3. Example 3
The purpose of this example is to validate the efficiency of the CQ-2x method compared with the implicit method and 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. The structural properties of this system are assumed to be and , in which . This system is subjected to a ground acceleration of . The highest natural frequency is 632.4 rad/s [31]. According to the highest natural frequency, a time step duration of is obtained from Newmark’s average acceleration method and is selected to be the exact solution. The responses of displacement and velocity versus time of the 200th DOF when N = 200 are plotted in Figures 33 and 34 . The analysis is accomplished by considering the time step duration equal to 0.02 s. The maximum displacement and velocity error ratios of the 200th DOF when n = 200 are listed in Table 3. It can be seen that the accuracy of the CQ-2x method is sufficient. The 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 method in which solving the inverse matrix twice for each time step is needed, as listed in Table 4.
6.4. Example 4
The 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). The following data are adopted: Young’s modulus E = 2.87 Gpa, , moment of inertia , mass per unit length , suspended mass , suspension stiffness , and speed . The frequency of vibration computed of the spring mass is . The frequencies of vibration computed of the bridge are . The maximum frequency of vibration computed of the bridge is . According to stability requirement of the CDM, the widest time step satisfying CDM’s stability is s. A time step duration of 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 , KR-, and CAA is only determined by the requirements of the calculation accuracy; here, s. The midpoint vertical deflection of the beam and deflection of the spring mass versus time are plotted in Figures 36 and 37. The 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. The efficiency of the CDM , , KR-, and CAA is compared. The 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 , 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. The CQ-2x algorithm are more efficient than the noniterative algorithm since the computation amount in each time step of CQ-2x is much less than that of KR-.
7. 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.
Appendix
A. Derivation of the Preliminary NDA Form
When , , and are chosen as the basic dimensional units, equation (5) can be expressed in terms of dimensionless quantities as follows:
To construct the function , a linear combination of , 1, , and can be used:where are coefficients. The effects of the dimensionless arguments and , which represent the structure’s material properties, on the computed displacements can be expressed by replacing the coefficients with linear combinations of , , and a constant term 1 as follows:where are coefficients. Substituting equations (A.3a)–(A.3c) into equation (A.2) yieldswhere . Equation (A.4) can be simplified by multiplying by as follows:
To reduce the number of undetermined coefficients , equation (A.5) can be further simplified by dividing both sides by :where are undetermined coefficients.
Equation (A.6) is the preliminary NDA form of the CQ-2x method.
B. Expressions for , , , , , and
The targets of computation may be either buildings or bridges, and the dynamic responses may be forced by wind or an earthquake; consequently, the values of , , and and the forms of , , , , and in the algorithm are arbitrary. Thus, by referring to equations (B.1a)–(B.1f), the following can be derived:
C. Derivation of the NDA Stability Requirements
The characteristic equation of the amplification matrix is given bywhere is the root of the characteristic equation of and
At the boundary line of the domain of NDA stability, the following holds:where is the spectral radius of the amplification matrix .
Equation (C3) being satisfied means that
For a complex-valued , equation (C5) is satisfied only whenwhere and is the argument of .
Substituting equation (C5) into equation (C1) yields
Considering , , and , equation (C6) can be rewritten as
Equation (C7) is satisfied only when
The boundary line of the NDA stability domain can be determined from equations (C8a) and (C8b) when . It can be derived that when equations (C8a) and (C8b)are satisfied when
When and , equations (C8a) and (C8b) are satisfied when
When , equations (C8a) and (C8b) are satisfied when
The boundary line of the NDA stability domain is thus composed of three lines described by equations (C9), (C10), and (C11) in the plane, as shown in Figure 1. Because is a continuous function of and , when and , the domain of NDA stability is the inner part of the triangle shown in Figure 1 and can be expressed as
D. Dissipation and Dispersion
The characteristic values of in equation (C1) are
The spectral radius of the amplification matrix is
When , equation (D1a) can be rewritten as [18]where the numerical circular frequency is given byand the numerical damping ratio can be written as
The relative period error is written as
Data Availability
All data included in this study are available from the first author [email protected] upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was jointly supported by the National Natural Science Foundation of China (Grant no. U1934207) and Innovation-Driven Project of Central South University (Grant no. 502501006).