Abstract
This paper presents a modified Precise Integration Method (PIM) for long-time duration dynamic analysis. The fundamental solution and loading operator matrices in the first time substep are numerically computed employing a known unconditionally stable method (referred to as original method in this paper). By using efficient recursive algorithms to evaluate these matrices in the first time-step, the same numerical results as those using the original method with small time-step are obtained. The proposed method avoids the need of matrix inversion and numerical quadrature formulae, while the particular solution obtained has the same accuracy as that of the homogeneous solution. Through setting a proper value of the time substep, satisfactory accuracy and numerical dissipation can be achieved.
1. Introduction
In structural dynamic analysis, direct time integration algorithms are widely used to solve the motion equation [1]. Many researchers [2–4] have focused on the improvement of numerical properties of time integration algorithms, such as numerical stability, accuracy, numerical dissipation, and computational efficiency.
An alternative way to evaluate the dynamic response is to use the Precise Integration Method (PIM-L) [5], where a numerical exponential matrix function is accurately computed to approximate the fundamental solution matrix and the applied loading is simulated by using piecewise linear functions. For linear dynamic homogeneous systems, the PIM-L can obtain precise numerical results approaching the exact solution. However, its accuracy is limited by the simulation of the applied loading in nonhomogeneous systems. Lin et al. [6] applied Fourier series expansions to simulate the applied loading and proposed the PIM-F method. Both the PIM-L and PIM-F require matrix inversion to obtain the particular solution. However, if the objective matrix is singular or nearly singular, its inverse matrix either does not exist or is hard to obtain steadily.
In order to avoid matrix inversion, Gu et al. [7] converted nonhomogeneous equations into homogeneous ones at the expense of dimensional expansion, which significantly increased storage and computational cost. Afterwards, Wang et al. [8] proposed a PIM based on Lagrange piecewise polynomials and improved the dimensional expansion approach. Another practical option to compute the particular solution without matrix inversion is to use the numerical quadrature formulae. Wang et al. [9, 10] applied Gauss quadrature, Simpson quadrature, and Romberg integration formulae to integrate the particular solution. The accuracy of particular solution using these methods depends crucially on the number of integration points used in each time-step and the size of time-step. Basically, the accuracy of the particular solution obtained by these aforementioned methods may not achieve the same accuracy as that of the homogeneous solution, especially when a large time-step is used.
In conventional PIMs, the fundamental solution matrix in the first time substep is approximated by the Taylor series expansion of the exponential matrix function. Fung [11] halved the computational effort required using step-response and impulsive-response matrices and considering the symmetry of these matrices. Later on, the Krylov subspace method [12] and the Ritz vectors [13] were also applied to improve the computational efficiency. Gao et al. [14] utilized the special feature of the dynamic system to reduce the computational effort involved in establishing the exponential matrix function. Wang et al. [9] investigated the numerical stability and computational accuracy of the conventional PIMs in detail. Subsequently, Wang et al. [10] proposed the PIM based on Padé approximation, where the generalized Padé approximation was used to replace the Taylor series expansion to achieve unconditional stability and controllable numerical dissipation.
This paper aims to modify and improve the existing PIM for solving long-time duration dynamics problems. The numerical fundamental solution and loading operator matrices are obtained by applying an unconditionally stable time integration method and recursive formulae. Using this proposed approach, the same accuracy of particular solution as that of the homogeneous solution can be achieved even for a large time-step, which is considered as a major improvement of PIMs. Besides, the proposed method can obtain approximately the same numerical results as those using the original method with a small time-step; thus it can also be considered as an accelerator of conventional time integration methods.
2. The Precise Integration Methods
The motion equation of a discretized linear structural model with n-degrees-of-freedom is given bywhere , , and are the mass, damping, and stiffness matrices, respectively; is the vector of applied external loading; , , and are the displacement, velocity, and acceleration vectors of the system, respectively.
Assuming , , (1) can be transformed into the following first-order ordinary differential equation:where . An analytical solution of (2) can be expressed as follows [15]: Applying time discretization with a time-step and assuming the analysis starts at the instant , can then be represented by : where , , and . can be written aswhere , is a given integer. The exponential matrix can be approximated by using the Taylor series expansion or the generalized Padé approximation. Taking the Taylor series expansion as an example, can be written aswherewhere is the truncation order of the Taylor series expansion. Substituting (6) into (5) gives can be evaluated by applying the following recursive formula:
3. A Modified PIM for Long-Time Duration Dynamic Analysis
3.1. The Numerical Fundamental Solution Matrix
In the standard finite-element semidiscrete systems, only a few low-frequency modes can represent the original continuous system. The quality of time integration algorithms depends crucially on their dissipative property. It has been found that the PIM is conditionally stable [9]. In order to satisfy the stability conditions, the time substep has to be smaller than the highest modal period of the particular discretized structure. In this case, the dissipative property would be eliminated. Wang et al. [10] proposed the PIM based on Padé approximations to achieve unconditional stability and controllable dissipation, in which the Taylor series expansion of the matrix exponential function is replaced with the generalized Padé approximations. In this work, the fundamental solution matrix is approximated using numerical Green’s function matrix to achieve desirable numerical properties. An analytical expression of , can be written in terms of , [16–19]:where , represent Green’s function matrix of the model and its derivative, respectively. Thus, the analytical expression of is given as follows: and can be calculated numerically by solving the following equation:
In order to satisfy the unconditional stability, a known unconditionally stable time integration method (referred to as original method in the following discussion) with time-step is used to solve (12). By selecting a proper original method and setting the value of time substep, one can obtain satisfactory accuracy, numerical dissipation, and other numerical properties.
Substituting into the following equation, can be obtained:
Substituting into (8) and (9), can be evaluated. On the other hand, there is another method to evaluate , which does not need the inverse of the mass matrix. The column of can be obtained through solving the following equation: where is the entry of the unity vector . is zero if and 1 if (, being the Kronecker delta).
It can be seen that the method incorporates the PIM with Green’s function matrix based method [16–19]. For this method, the sparsity of the fundamental matrix would be weakened as the number of time substeps increases. As a result, this method requires more storage usage compared to conventional time integration methods. However, much larger time-step can be used, which causes less computational cost. Users can find a balance between storage usage and computation cost by setting proper . The improvement of the storage efficiency is feasible, which needs further investigation.
3.2. The Loading Operator Matrix
For homogeneous systems, the PIMs can give a precise numerical solution. However, for nonhomogeneous dynamic systems, the accuracy of numerical results depends essentially on the accuracy of the particular solution. The accuracy of the PIM-L and PIM-F methods is mainly limited by the error in the matrix inversion and the applied loading simulation. For the PIM proposed by Wang et al. [10], the main error was caused by the process of the numerical quadrature formulae. In order to obtain a high accuracy, a large number of integration points were used, which required extra computational effort. When a larger time-step is adopted, the particular solution obtained may not attain the same accuracy as that of the homogeneous solution. In this study, a novel method is proposed to evaluate the particular solution, which avoids the matrix inversion. It can retain the same accuracy of the particular solution as that of the homogeneous solution even though the time-step is very large.
In practical engineering problems, abounding loadings can be expressed in a periodical form, which can be approximated by using the Fourier series expansion where is the truncation order of the Fourier series expansion. According to Euler’s formula, the trigonometric functions in (15) can be transformed as an exponential function. An exponential function loading can be expressed as follows in the time interval ,where is a constant vector and is a real or complex number. Substituting (17) into (4), the following equation can be obtained:Assuming , (18) becomeswhere is the loading operator matrix. It can be seen that is not associated with the integration interval , thus it is calculated only once during the whole step-by-step integration procedure. As (19) is the analytical expression of the particular solution, once is properly computed, the particular solution would become extremely accurate. In order to ensure the same accuracy of the particular solution as that of the homogeneous solution, in (20) can be approximated by , which is generated in Section 3.1. can be written as It can be seen that is the sum of terms. Let denote the sum of the first terms of (21). The following recursive formula can be used to evaluate from : where . is obtained in Section 3.1. The matrices and are computed using the following looping algorithm:For conventional PIMs, the algorithm (23) can also be used to acquire . Thus the particular solution would possess the same accuracy as the homogeneous solution. In this work, is obtained numerically by solving the following equation using the original method with step-size :
During the recursive process(22), only matrix-vector multiplications are required. Therefore, it can be summarized that only times of matrix-matrix multiplication are needed to obtain the numerical fundamental solution and loading operator matrices (as the computational effort of matrix-vector multiplication is negligible compared to that of matrix-matrix multiplication, only matrix-matrix multiplication is considered).
3.3. Properties Analysis
As the recurrence relations ((9) and (22)) are analytical expressions, the accuracy, stability, dissipation, and other numerical properties of the proposed method are determined only by the original method and the size of time substep. In this section, the stability of the proposed method is investigated. For the analysis of other properties, one can refer to [16] for detailed information.
Considering the following undamped unforced single degree of freedom (SDOF) problem:where is the natural frequency of the model and the natural period is . Assuming the time-step is , time substep is (). The following recursive relationship can be obtained using the original method:where stands for the amplification matrix of the original method. According to (5), the amplification matrix of the proposed method is obtained as . The stability condition is , where is the spectral radius of and are the eigenvalues of .
In order to investigate the stability of the proposed method, the Newmark (, ), the Bathe [2] and the fourth-order Runge-Kutta (RK-4) methods are selected as the original method, respectively.
Figure 1 shows the spectral radius of the proposed method, where MPIM-Newmark, MPIM-Bathe, and MPIM-RK denote the results obtained using the Newmark, Bathe, and RK-4 method as the original method, respectively. For all curves the value of spectral radius is 1 until is equal to a certain value and thereafter it rapidly diminishes. As increases, this value increases. The MPIM-Newmark has lower limit values than the Newmark method. Compared to RK-4, the stability condition of MPIM-RK4 is relaxed.
(a)
(b)
(c)
4. Numerical Examples
In this section, two numerical examples are investigated to evaluate the validity of the proposed method. The proposed method solution is compared with those obtained by the SPIM (, , ), the ExGA-RK method [16], and the original method with a small time-step. The Newmark (, ), the Bathe [2], and the fourth-order Runge-Kutta (RK-4) methods are selected as the original method, respectively. The RK-4 method is selected just for the sake of the comparison with the ExGA-RK method. An unconditionally stable method can be used when solving particular dynamic problems.
In Figures 2–11, MPIM-RK, MPIM-Newmark, and MPIM-Bathe denote the solution obtained by the proposed modified PIM using the RK-4, Newmark, and the Bathe method as the original method, respectively.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
4.1. Two-DOF Model
The first example considers the following two-DOF system:
The natural periods of the system are (some values are rounded). In the following cases, the solution obtained by the mode superposition method is referred to as the reference solution. In Case 1, a constant loading is used, while the applied loading is a sine function in Case 2.
Case 1. The initial conditions are set as , , , and with . In this case, the time-step is selected as , the recursive parameter , and the time substep .
Figures 2 and 3 show the simulation results obtained using various methods. The error is defined as follows:where is the response obtained using numerical methods and is the reference solution. The ExGA-RK method employs the RK-4 to calculate the fundamental solution matrix and Simpson’s 1/3 rule to obtain the particular solution. The SPIM applies the generalized Padé approximation (, ) to calculate the fundamental solution matrix and the repeated Simpson integration formulae (using 8 integration points) to obtain the particular solution. In this case, the applied loading is a constant vector and the integrand of the particular solution remains unchanged. As a result, all methods perform very well. The error follows the sequence of MPIM- - Bathe > SPIM >ExGA-RK=MPIM-RK.
The accuracy of numerical fundamental solution matrices obtained by different methods is also investigated by fixing the time-step and choosing different recursive parameters . The error of is defined as follows:where is the numerical fundamental solution matrices obtained using different methods and is the reference fundamental solution matrix. Figure 4 depicts the error versus the recursive parameters . It can be seen that the error decreases with the increasing of . With the same , the error follows the sequence of MPIM- Newmark >MPIM- Bathe >SPIM>MPIM-RK.
Case 2. The initial conditions are set as , , , and and the load . In this case, the time substep is set as . The recursive parameter is chosen as . Thus and , where is the period of the applied loading.
Figures 5–10 show the simulation results using different . It can be seen that the errors of the ExGA-RK, SPIM solution increase as the time-step increases, while the solutions obtained by the proposed methods (MPIM-Bathe, MPIM-RK) is always the same as the reference solution. For instance, for , the loading is zero at the start and end points of each time-step; thus it would be considered as zero within each time-step. Although the SPIM uses 8 integration points to calculate the particular, the numerical results have large errors when .
In this case, the numerical convergence property of the proposed method is also investigated by calculating the time order of convergence rate. The time order of convergence rate is calculated by applying the following formula:where , are different time-steps, and , are the maximum error of corresponding numerical result. Figure 11 shows the log of the versus the log of time-step for displacement and velocity of node 2. It can be seen that the slopes of SPIM-RK and SPIM-Bathe curves remain unchanged as increases. The slopes of SPIM-RK and SPIM-Bathe curves are about 3.99 and 1.97, respectively.
4.2. Regular Diagonally Positioned Pyramid
The second example considers the regular diagonally positioned pyramid shown in Figure 12, which is composed of aluminium tubes. The section dimensions and the material properties are as follows: the outer diameter , the wall thickness , the density , the Poisson ratio , and Young’s modulus . The top level of the pyramid is fixed along four sides, while node 1 is subjected to a vertical load . By using spatial discretization (spatial bar element), a structural model with 543 degrees of freedom is obtained. A consistent mass matrix is used in the analysis and the material is assumed to be linearly elastic. The damping matrix is assumed to be (Rayleigh damping), where the coefficients and [3]. The maximum natural period of the model is .
The time-step is set as , the recursive parameter , and the time substep . Thus , , where is the period of the loading. The solution obtained by the Bathe method with time-step is set as the reference solution. The model was solved for over 10 seconds. Table 1 shows the CPU time occupied by different methods (using a desktop PC with Inter(R) Core(TM) i5-5200U CPU @ 2.2 GHz, Ram@16G). The time-step for the RK-4 and Bathe method is . It can be seen that the proposed method requires much less CPU time than the original method with the small time-step.
Figure 13 shows the simulation results of the vertical displacement of node 1 and node 2, respectively. Similar conclusions can be obtained for other nodes. It demonstrates that the MPIM-Bathe and MPIM-RK solution exactly follows the reference solution. As the time-step is set as nearly half the period of the load, larger error occurs in the ExGA-RK and SPIM solutions.
(a)
(b)
From the foregoing discussion, it can be concluded that the accuracy of the proposed method is determined by the size of time substep and the properties of the original method. Even with a large time-step, the particular solution has the same accuracy as the homogeneous solution. On the premise of obtaining almost the same numerical result, the proposed method can greatly improve the calculative efficiency of the original method.
5. Discussion
In recent years, many integration methods have been applied to calculate Green’s function matrix, resulting in different versions of time marching method, such as the Newmark family methods [17], the generalized- method [18], the Runge-Kutta, and the central difference schemes [16]. Using a substepping procedure, these aforementioned approaches can obtain accurate homogeneous solution. However, the technologies used to evaluate the particular solution are basically based on some approximations that are only suitable for a small time-step. For instance, in [16] the external loading vector is assumed to be linear functions in each time-step. In [19] the integrand of the particular solution is assumed to be a linear or parabolic function within each time-step. When a larger time-step is considered, these approximations are oversimplified and results in unexpected errors. For instance, the time-step is set as an integral multiple of the period of the periodical loading; the loading would be considered to be unchanged within each integral interval. As a result, the accuracy of the particular solution obtained would be inadequate compared to that of the homogeneous solution.
It is widely known that the computational effort required is inversely proportional to the size of time-step in the structural dynamic analysis. Compared to these aforementioned approaches, the proposed method can use a large time-step and can preserve the same accuracy of particular solution as that of the homogeneous solution; thus it is more suitable for the long-time duration dynamic analysis.
6. Conclusions
An improved version of PIM is presented in this paper. This approach replaces the Taylor series expansion or the generalized Padé approximation of the matrix exponential function in the conventional PIMs with numerical fundamental solution matrix and employs loading operator matrix to obtain particular solutions. Efficient recurrence expressions are developed to evaluate the loading operator matrices in the first time-step. The highlights of the approach can be summarized as follows,(1)The proposed method avoids matrix inversion operations and numerical quadrature formulae; thus it expands the scope of application. Meanwhile, the accuracy of particular solution obtained can be maintained the same as that of the homogeneous solution even for a large time-step.(2)As the intermediate variables required for evaluation of the loading operator matrix have already been generated in calculating the fundamental solution matrix, the proposed method improves the computational efficiency of conventional PIMs.(3)By selecting an appropriate original method and time substep, the satisfactory numerical properties, such as accuracy, numerical dissipation, and unconditional stability, can be obtained.(4)Although the proposed method uses a large time-step, it can obtain almost the same numerical results as the original method using a small time-step. Thus, the proposed method is actually an accelerator of conventional time integration methods.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant nos. 11672338 and 11202247).