A Fast Frequency Domain Method for Steady-State Solution of Forced Vibration of System with Complex Damping

The conventional frequency domain method (CFDM) and dual-force-based time domain method (DTDM) are often used to solve the steady-state response of system with complex damping under an arbitrary force. However, the calculation efficiency of the DTDM is low due to the straightforward summation operation of series even if the solution of the DTDM is the exact real part of the solution. In addition, since the CFDM only can obtain the real part of solution not the complete solution, it gives misleading information that the solution does not have an imaginary part. In this paper, a fast frequency domain method (FFDM) is proposed to calculate the complete response of complex damping system including the imaginary part with a higher accuracy in a much faster manner. The new FFDM uses half of the Fourier series of the discrete Fourier transform of the actual arbitrary force to construct the Fourier series of the dual force, followed by calculating the time history response using the inverse fast Fourier transform. The new developed method is validated through three numerical examples with harmonic and seismic excitations. The numerical results show that the accuracy of the new FFDM is compatible to the DTDM but with much higher computational efficiency.


Introduction
In structural dynamics, the damping force of a complex damping system [1][2][3] which is also known as hysteretic damping, rate-independent damping, structural damping, f d , can be calculated as: in which i = √ −1 is a pure imaginary number, and k is the stiffness of the structure, η is the complex damping factor, and u is the displacement. From Equation (1), it can be seen that the energy dissipation of the complex damping is independent of the excitation frequency, which is satisfied with the energy dissipation characteristics of many materials and has been widely used. Taking soil as an example, overserved from numerous dynamic triaxial tests, the dissipated energy per cycle of many kinds of soil is almost independent of the loading frequency [4,5]. Laboratory shaking table test results of a single-layer reticular dome [6] showed that comparing with the identified damping ratios of the first ten vibration modes, complex damping model has smaller errors than Rayleigh damping model. The generalized complex damping is used to uncouple the coupled damping matrix caused by the offshore and added hydrodynamic damping [7]. In addition, the damping of sandwich-like plate and beam with viscoelastic cores is usually also expressed as complex damping [8,9].
For the forced vibration of system with complex damping, two methods are commonly used to obtain the steady-state responses of such a system, including the conventional frequency domain method (CFDM) and the dual-force-based time domain method (DTDM). The seismic response of soil layer and soil-structure interaction [10][11][12] are often calculated by the CFDM due to the complex damping of soil deposit. For example, the CFDM has been used to compute the nonlinear seismic response of one-dimensional soil layers [13,14], two-dimensional irregular-shaped basins with horizontal soil layers [15], and a hill in a layered elastic half space [16]. The CFDM also was used to analyze complex damping system under forced vibrations of axial and torsional forces [17].
In theory, the differential equation with complex coefficients can be equivalent to two differential equations with real coefficients which use the real part and imaginary part of solution as variable. The real and imaginary parts of the solution not only affect each other but also have a dual relationship. The real parts of displacement response and force are observed and measured. Although the imaginary part is not the actual process, it is the complete condition for the equation with the complex constitutive. In the process of establishing and solving the equation, the influence of the imaginary part should be considered, although in most of the current engineering applications, the imaginary part may not be necessary to be applied due to complexity in solving the imaginary part. Therefore, methods using frequency domain which can only obtain the real part of the solution satisfy the engineering demand. To obtain both the real and imaginary parts of the solution, the DTDM can be applied. Previous research demonstrated that the solution of differential equation with complex coefficients in the time domain analysis under the real-valued force is inaccurate, and the imaginary part of the complex force corresponding to real-valued force should be established according to the dual principle [18]. Thus, the calculation method of dual force was proposed to be improved [19]. However, since the eigenvalues of the complex damping system are complex eigenpairs that are opposite numbers to each other [20], it is inevitable that the real part of one eigenvalue is positive. Since the real part of eigenvalue is a positive quantity, the complementary solution would be divergent, resulting in contrary to the physical rule. For an analytical solution, the divergent term can be manually deleted to make the calculation result stable. The direct integration method of the complex damping system converges to the sum of particular and complementary solutions, which still contains the divergent terms and will induce unstable results [21,22]. Therefore, in the DTDM, the complex force is expressed as the summation of harmonic functions from discrete Fourier transform (DFT), and the complex responses for each harmonic force element are calculated. The response of the complex system can then be obtained by summing all of the complex solutions from all harmonic force elements. The procedure of the DTDM indicates one of its main drawbacks of low computational efficiency due to the direct summation operation of series. To improve the computational efficiency, dual force has been suggested for both the CFDM and DTDM [23], however, not yet fully investigated.
In this paper, a fast frequency domain method (FFDM) without constructing dual force is developed, which can simultaneously compute the imaginary part of the solution to comprehensively realize the dual relation in the complex system. In addition, the comparison of the differences between the CFDM and DTDM shows that the dual force is not necessary for the CFDM to calculate the real part response. The comparison of numerical examples and theoretical solutions on the steady-state response of system under harmonic and seismic excitations demonstrated that the developed FFDM has a high accuracy and fast computational efficiency compared with the CFDM and DTDM.

Theoretical Background
The governing equation of motion for a single-degree-of-freedom (SDOF) system with complex damping in forced vibration can be expressed as [20]: in which m and k are the mass and stiffness, respectively,ü and u are the complex acceleration and displacement of the SDOF system, respectively, and f (t) is the real-valued force whose duration is T p . To solve Equation (2) in the frequency domain, the duration, T p , is divided into N equal intervals, ∆t, where N equals to 2 γ and γ is an integer from fast Fourier transform (FFT). The discrete time, t n , can then be calculated as t n = n∆t (n = 0,1,2, . . . , N). The response of displacement at a discrete time instant t n is u n . The applied force f (t) at a discrete time instant t n is f n . Thus, the discrete Fourier transform (DFT) of force f n , F j = F(θ j ), can be given as: in which, θ j is the discrete frequency which can be obtained as θ j = j∆θ with ∆θ = 2π N∆t . In the DFT, the highest effective frequency is the Nyquist frequency θ N/2 . The amplitude F j of j > N/2 is the conjugate of that of j < N/2 on the Nyquist frequency, that is, in which, the superscript "-" denotes complex conjugation. Under harmonic excitation of F j e iθ j t n ( j = 0, 1, 2, . . . , N/2), the displacement amplitude of steady-state response can be expressed as: in which H j is the discrete form of the complex frequency response function H(θ j ). For a single degree-of-freedom (SDOF) system with complex damping in forced harmonic vibration as shown in Equation (2), the frequency response function H j can be given as [24]: With an excitation frequency of θ N/2+ j ( j = 1, 2, . . . , N/2−1), the amplitude of the steady-state response can be directly calculated by the conjugate method according to the Nyquist criterion, that is, Then based on the conventional frequency domain method, the inverse Fourier transform should be used to evaluate the dynamic response in the discrete form, which can be expressed as: in which F 0 = 1 N N−1 n=0 f n , Re denotes the real part of the dynamic displacement response. Specifically, the real and imaginary parts (u Rn and u In ) of dynamic displacement response can be expressed as: It can be seen from Equations (9)-(11) that the complex dynamic displacement response from Equation (11) is very small. Thus, for the conventional frequency domain method (CFDM), only the real part of the dynamic displacement response is considered using Equation (8).

Development of Fast Frequency Domain Method
To overcome the limitation of the CFDM for ignoring the complex response of the system and obtain both the real and imaginary parts of the dynamic response, the dual imaginary part of applied force can be brought back to Equation (2). In this case, Equation (2) with the governing equation of motion under complex force excitations can be rewritten as: in which g(t) is dual force which is constructed to match the measured real-valued force f (t). The force g(t) sampled at the discrete time t n is g n . Using the discrete Fourier transform, f n can then be expressed as the summation of simple harmonic functions, that is: where, According to the dual criterion, g n can be expressed as: Appl. Sci. 2020, 10, 3442

of 13
Thus, the complex coefficient Equation (12) can be decomposed into a set of two differential equations with real coefficients: , (17) where, u R = Re(u) and u I = Im(u) represent the real and imaginary parts of the displacement u, respectively. In the DTDM, the real part and the imaginary part of response satisfy dual relationship and are the exact solution. In the DTDM, the steady-state solution of Equation (17) can be given as: Comparing the results of Equations (10) and (11) with Equations (18) and (19), it can be seen that the real part of the displacement obtained by the CFDM and the DTDM are same except the first term. The CFDM only obtain the real part solution of the j = 1, 2, . . . , N/2 − 1 term, resulting in two issues in the CFDM. The first issue of the CFDM is that the real part of the calculation result is reasonable, but the imaginary part solution of the response cannot be obtained. The second problem is that when the force is non-zero mean value, the CFDM would have calculation error because of the first term of Equation (10). For the DTDM, it has two problems as well, including a need to construct dual force in advance and the extremely intensive computation efforts needed to obtain the solution by direct use of Equations (18) and (19).
To address the problems of the two above mentioned algorithms, both the CFDM and the DTDM, a fast frequency domain method (FFDM) can be used to make the solution accurate and meanwhile reduce the computation efforts. Comparing Equation (3) with Equations (14) and (15), the F j , a j , and b j would satisfy the following relationships: Therefore, the complex load can be expressed as: In the FFDM, a new Fourier series D j is defined by the force Fourier series F j in Equation (3) as follows: Then the dynamic displacement response at time t n can be rewritten as: Appl. Sci. 2020, 10, 3442 6 of 13 Comparing Equation (23) to Equations (18) and (19), it can be seen that they are the same. In fact, the second term on the right of Equation (23) is usually a small quantity because the amplitude F N/2 of the highest effective frequency in the Fourier spectrum is far less than that of the prominent frequency.
by ignoring the second term on the right of Equation (23), an approximate solution can be obtained with negligible error, which is To express conveniently, Equation (23) is called the accurate FFDM and Equation (25) is called the approximate FFDM. The first term on the right of Equation (23) can be calculated by inverse fast Fourier transform to significantly improve the calculation efficiency. Therefore, the FFDM not only retains the efficiency of the CFDM, but also is as accurate as the DTDM. At the same time, there is no need to construct the dual load in the FFDM which improves the calculation efficiency.
For multi-degree-of-freedom (MDOF) system with complex damping, the governing equation of motion can be expressed as: in which [m] and [k] are the mass matrix and stiffness matrix, respectively. f (t) is the time variation of applied force, and {s} is spatial distribution vector. Given the natural frequencies ω n and associated mode shapes φ n (n = 1, 2, . . . , N), the displacement of the structure can be expressed into a linear combination of N modal coordinates. That is, Based on the mode superposition method, Equation (26) can be transformed into a series of generalized SDOF whose vibration equation for the q n (t) can be written as is the modal stiffness, γ n = − φ T n {s}/M n is the modal participation factor, q n is the modal coordinate. The computational efficiency of MDOF system is depended on that of the SDOF system by the mode superposition method. Therefore, the following is to mainly investigate the efficiency and accuracy of SDOF system.

The Algorithm Workflow
In order to make the procedure of the algorithm clearer, the workflow of the algorithm described above is shown in Figure 1, which is executed via three steps: Step 1: Given m, k, η, and f n , calculate the Fourier coefficients F j from Equation (3).

Numerical Analysis
To investigate the effectiveness and accuracy of the developed new FFDM, three examples with harmonic excitation with a zero and non-zero mean value, and seismic force are numerically analyzed and discussed, respectively, in this section.

Harmonic Force with a Zero Mean Value
The first example in this study uses a harmonic excitation with a zero mean on a single-degreeof-freedom (SDOF) system. The harmonic excitation with a zero mean can be expressed as a sine  The exact steady solution of such a dynamic problem can be expressed as:

Numerical Analysis
To investigate the effectiveness and accuracy of the developed new FFDM, three examples with harmonic excitation with a zero and non-zero mean value, and seismic force are numerically analyzed and discussed, respectively, in this section.

Harmonic Force with a Zero Mean Value
The first example in this study uses a harmonic excitation with a zero mean on a single-degree-of-freedom (SDOF) system. The harmonic excitation with a zero mean can be expressed as a sine force, f (t) = sin θt, with θ/ω = 1, η = 0.1, m = 1 kg, k = 4π 2 /1.

Numerical Analysis
To investigate the effectiveness and accuracy of the developed new FFDM, three examples with harmonic excitation with a zero and non-zero mean value, and seismic force are numerically analyzed and discussed, respectively, in this section.

Harmonic Force with a Zero Mean Value
The first example in this study uses a harmonic excitation with a zero mean on a single-degreeof-freedom (SDOF) system. The harmonic excitation with a zero mean can be expressed as a sine  The exact steady solution of such a dynamic problem can be expressed as: The exact steady solution of such a dynamic problem can be expressed as: However, if this problem is solved without dual forces, the governing equation of motion can be expressed as: To express conveniently, Equation (30) is called the direct time domain method (Direct TDM). The dynamic response calculated by the Direct TDM can be expressed as: When comparing Equation (31) with Equation (29), it can be seen that they are not equal to each other. Thus, to compare the difference between different methods, this given example was solved by using five different methods, including the Direct TDM, the CFDM, the DTDM, the new accurate FFDM and the new approximate FFDM. Figure 3 compares the dynamic displacement responses from the exact solution as in Equation (29)

Harmonic Force with a Non-Zero Mean Value
The second example uses a non-zero mean harmonic excitation for a SDOF system. Specifically, the system is subjected to a cosine force, ( ) 1-cos f t t  = with /1  = , η = 0.1, m = 1kg, k = 4π 2 /1.28 2 N/m. The calculation duration is 20.48s and the time interval t  is 0.01s. Figure 4 shows the time history and Fourier spectrum of non-zero mean cosine force. The relative error e of the amplitude of displacements can be expressed as in which u and u* are the approximate and exact solutions, respectively. The relative errors of the amplitude of displacements are also shown in Table 1. It can be seen from Figure 3 and Table 1 that the real part and imaginary part of displacement obtained by the DTDM, the accurate FFDM and the approximate FFDM are identical with the exact solution. This is the reason that the DTDM and the accurate FFDM are identical and equal to the discrete values of exact solution. As for the approximate FFDM, the D N/2 is equal to zero for sine function, therefore, ignoring of the last term of Equation (23) would not cause errors.
As for the CFDM, the real part obtained by CFDM is identical with the exact solution. The Fourier series of U j is conjugated symmetric on the Nyquist frequency in Equation (8), therefore the imaginary part of displacement is approximately equal to zero after summation which is different with the exact solution. Fortunately, the real part of solution is the real response of system, and the imaginary part is a dual term for the real part and not used to decide the status of real system. Therefore, the CFDM can be used to solve the response of system with complex damping. For the direct time domain solution, the real part of the solution is equal to zero, which is obviously incorrect because the dual relationship of the load is ignored. It means that the direct time domain solution is not analyzed in the subsequent examples.

Harmonic Force with a Non-Zero Mean Value
The second example uses a non-zero mean harmonic excitation for a SDOF system. Specifically, the system is subjected to a cosine force, f (t) = 1 − cosθt with θ/ω = 1, η = 0.1, m = 1 kg, k = 4π 2 /1.28 2 N/m. The calculation duration is 20.48 s and the time interval ∆t is 0.01 s. Figure 4 shows the time history and Fourier spectrum of non-zero mean cosine force.

Harmonic Force with a Non-Zero Mean Value
The second example uses a non-zero mean harmonic excitation for a SDOF system. Specifically, the system is subjected to a cosine force, ( ) 1-cos f t t  = with /1  = , η = 0.1, m = 1kg, k = 4π 2 /1.28 2 N/m. The calculation duration is 20.48s and the time interval t  is 0.01s. Figure 4 shows the time history and Fourier spectrum of non-zero mean cosine force. The exact steady solution of this example can be expressed as: The exact steady solution of this example can be expressed as: This example problem was solved by using four different methods including the CFDM, the DTDM, the new accurate FFDM, and the new approximate FFDM. Figure 5 compares the exact solution in Equation (33) with the numerical results from the four different approximate methods. Table 2 lists the relative errors of the amplitude of displacements of various method. The real part and imaginary part of the dynamic displacement responses obtained by the DTDM, the accurate FFDM and the approximate FFDM are also identical with the exact solution. The real part obtained by the CFDM is only slightly different from the exact solution, which caused by the small difference between Re[1/k(1 + iη)] and 1/k. The imaginary part of the displacement response obtained by the CFDM is completely different with the exact solution, indicating that the DTDM and the FFDM are more accurate than the CFDM.
This example problem was solved by using four different methods including the CFDM, the DTDM, the new accurate FFDM, and the new approximate FFDM. Figure 5 compares the exact solution in Equation (33) with the numerical results from the four different approximate methods. Table 2 lists the relative errors of the amplitude of displacements of various method. The real part and imaginary part of the dynamic displacement responses obtained by the DTDM, the accurate FFDM and the approximate FFDM are also identical with the exact solution. The real part obtained by the CFDM is only slightly different from the exact solution, which caused by the small difference between  

Re 1/ (1 )
ki  + and 1/k. The imaginary part of the displacement response obtained by the CFDM is completely different with the exact solution, indicating that the DTDM and the FFDM are more accurate than the CFDM.

Seismic Excitation
The third example uses seismic excitation with ground acceleration () g ut on a SDOF with m = 1kg, k = 4π 2 /1.28 2 N/m, and η=0.1. The equivalent real-valued force f(t) is: In this example, the acceleration () g ut was the north-south component of the ground motion recorded at a site in El Centro, California during Imperial Valley earthquake of May 18, 1940. Figure  6 shows the ground motion and its Fourier spectrum. The time interval of seismic wave sampling is 0.02s and the end of the force is filled with zero after the duration becomes 40.96s.

Seismic Excitation
The third example uses seismic excitation with ground accelerationü g (t) on a SDOF with m = 1 kg, k = 4π 2 /1.28 2 N/m, and η = 0.1. The equivalent real-valued force f (t) is: In this example, the accelerationü g (t) was the north-south component of the ground motion recorded at a site in El Centro, California during Imperial Valley earthquake of 18 May 1940. Figure 6 shows the ground motion and its Fourier spectrum. The time interval of seismic wave sampling is 0.02 s and the end of the force is filled with zero after the duration becomes 40.96 s. Figure 7 shows the comparison of the dynamic displacement responses from the developed FFDM with the DTDM and the CFDM under El Centro ground motion excitation. Assuming that the exact displacements obtained by the DTDM, Table 3 lists the relative errors of the amplitude of displacements of various methods. Since D N/2 = 2.0819, max D j = 245.7599, that is D N/2 << max D j , ignoring the effect of highest effective frequency in Fourier series of force in approximate FFDM would induce very little error. The real and imaginary parts obtained by the accurate FFDM and the approximate FFDM are perfectly in agreement with those by using the DTDM. Therefore, it is reasonable to use the approximate FFDM to analyze the seismic response of system with complex damping. Since the El Centro wave time history applied is a zero-mean value, the real part obtained by the CFDM coincides with those by the DTDM, but the imaginary part is almost equal to zero. This finding further confirms that there exists dual imaginary part, although it has not physical meaning, it has theoretical meaning to make the solution complete.  The real and imaginary parts obtained by the accurate FFDM and the approximate FFDM are perfectly in agreement with those by using the DTDM. Therefore, it is reasonable to use the approximate FFDM to analyze the seismic response of system with complex damping. Since the El Centro wave time history applied is a zero-mean value, the real part obtained by the CFDM coincides with those by the DTDM, but the imaginary part is almost equal to zero. This finding further confirms that there exists dual imaginary part, although it has not physical meaning, it has theoretical meaning to make the solution complete.     Figure 7 shows the comparison of the dynamic displacement responses from the developed FFDM with the DTDM and the CFDM under El Centro ground motion excitation. Assuming that the exact displacements obtained by the DTDM, Table 3 Table 4 compares the computation time of all the four different numerical methods analyzed in all three examples. It can be seen that the time consumed by the developed FFDM is significantly shorter than the CFDM and the DTDM. The DTDM is the slowest due to the fact that a straightforward evaluation of the summation of Equations (18) and (19) would require lots of multiplications. The main computation time for the CFDM contains two parts, one for solving the FFT of force and the other for solving the IFFT of displacement. The procedure of the approximate FFDM is the same as the CFDM, but the IFFT series in the approximate FFDM is only half of those in the CFDM. In the accurate FFDM, only one extra term adding calculation is demanded on the last term of Equation (23) compared with approximate FFDM, which leads to the time consumed by the accurate FFDM is more than that by approximate FFDM, but less than by the CFDM. Therefore, the approximate FFDM outperforms the others when considering both the efficiency and accuracy.

Conclusions
In this paper, a fast frequency domain method is developed for the steady-state solution of the forced vibration of system with complex damping. From the results of theoretical analysis and simulation results for three examples, the following conclusions can be drawn: (1) For the complex damping system, the CFDM can obtain the exact solution of the real part of the response under zero mean load, but there is a small error in the non-zero mean case. The result obtained by the accurate FFDM is the same as that by the DTDM. Since the influence of the highest effective frequency of the Fourier spectrum is very small, the approximate FFDM still has high accuracy.
(2) The direct time domain method ignoring the effect of dual force would lead to inaccurate results. The CFDM can only obtain the real part solution of the response but not the imaginary part solution. Although the imaginary part may not have a physical meaning, it has theoretical meaning to make the solution complete. The FFDM can obtain the real part and the imaginary part of the response simultaneously.
(3) The FFDM has the same procedure as the CFDM, but only half of Fourier series are used in IFFT. Moreover, it can accurately calculate the real and imaginary parts of response, which is the more efficient algorithm for complex damping system, thus, recommended for practical applications.
In the future, more researches are needed to analyze the forced vibration for multi-degrees-of-freedom system with complex damping.