Mathematical Problems in Engineering

Volume 2018, Article ID 9646017, 12 pages

https://doi.org/10.1155/2018/9646017

## A Modified Precise Integration Method for Long-Time Duration Dynamic Analysis

School of Engineering, Sun Yat-Sen University, No. 135 West Xingang Road, Guangzhou, China

Correspondence should be addressed to Minghui Fu; nc.ude.usys.liam@hmfsts

Received 10 February 2018; Revised 24 May 2018; Accepted 30 May 2018; Published 14 June 2018

Academic Editor: Giovanni Falsone

Copyright © 2018 Ce Huang and Minghui Fu. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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