An efficient time-history analysis method for linear complex structures based on optimized integral parameters

Abstract

The present invention discloses an efficient time-history analysis method for linear complex structures based on optimized integration parameters. The optimal integration parameters of the integration iteration gain matrix are obtained by considering the characteristics of the structure itself, making it have the characteristics of unconditional stability and controllable high-order numerical damping. The specific steps are as follows: According to the discrete time-domain formula of the structural dynamics equation, set the integration parameters, establish the general form of the integration algorithm, and further obtain the iteration gain matrix; obtain the initial parameters and initial conditions according to the characteristics of the structure itself, and at the same time set the integration time step; according to the initial parameters and initial conditions of the structure, use optimization algorithms such as genetic algorithms and particle swarm optimization algorithms to obtain the optimal integration parameters; substitute the obtained optimal integration parameters into the iteration gain matrix to obtain the response of the structure at the next step; update the initial conditions of the structure, substitute them into the discrete time-domain motion equation for solution until the preset time is reached, and output the response time-history curve of the structure.

Description

Technical Field

[0001] The present invention belongs to the field of solving discrete-time motion equations of structural dynamics in civil engineering, and particularly relates to an efficient time-history analysis method for linear complex structures based on optimized integration parameters. Background Art

[0002] In structural dynamics, accurate solution of partial differential equations is the basis for structural dynamics analysis of various systems. The dynamic response analysis of structural systems mainly adopts two methods: modal decomposition method and direct integration method. Due to its high-precision characteristics, the modal decomposition method is widely used in the dynamic analysis of simple linear structures. For the transient response analysis of large complex structures, the direct integration method has the advantages of avoiding high-order mathematical analysis and being easier to calculate, and has more obvious applicability. Therefore, the direct integration method is widely used in the solution of various complex civil engineering problems such as structural dynamics analysis, finite element analysis, and real-time substructure testing.

[0003] Direct integration algorithms can be divided into two categories: explicit integration algorithms and implicit integration algorithms according to their characteristics. The implicit integration algorithm has the advantages of unconditional stability and controllable numerical damping, and is widely used in linear systems with low modal responses; however, the implicit method requires repeated construction and decomposition of the structural system matrix until a specific convergence criterion is met, so the computational cost is large and it is not applicable to complex structures. The explicit algorithm avoids unnecessary matrix decomposition and convergence determination, and has the advantage of high computational efficiency, and is widely used in the dynamic analysis of various structural systems.

[0004] However, explicit integration algorithms, such as the Newmark method, are not unconditionally stable and require extremely small time steps to ensure their stability, which obviously increases the computational cost of structural dynamic analysis. To solve this problem, in recent years, many experts and scholars at home and abroad have achieved certain results in developing new explicit integration algorithms, among which the model-based integration algorithms have attracted much attention, such as the Chang algorithm and the CR algorithm. At present, although the model-based integration algorithms can ensure unconditional stability, they have no numerical damping in high-order modes and cannot guarantee their computational accuracy. Therefore, it is necessary to develop a new model-based integration algorithm that can ensure computational efficiency, computational accuracy, and applicability. Summary of the Invention

[0005] To solve the above problems, the present invention proposes an efficient time history analysis method for linear complex structures based on optimized integration parameters. By using an optimization algorithm and considering the structural characteristic parameters (mass and stiffness) of the structure, the optimal integration parameters of the integration iteration gain matrix are obtained through the optimization algorithm, thereby updating the iteration gain matrix in the integration process. The algorithm has the characteristics of unconditional stability and controllable high-order numerical damping, and has the advantages of high precision, high computational efficiency, and wide application range. It is of great significance for the structural dynamics analysis, finite element analysis, and real-time substructure test of linear complex systems.

[0006] The technical solution of the present invention is as follows:

[0007] An efficient time history analysis method for linear complex structures based on optimized integration parameters. First, according to the discrete time domain formula of the structural dynamics equation, integration parameters are set, the general form of the integration algorithm is established, and the iteration gain matrix is further obtained. Secondly, the initial parameters and initial conditions are obtained according to the structural characteristics, and the integration time step is set at the same time; then, according to the structural initial parameters and initial conditions, an optimization algorithm is used to obtain the optimal integration parameters; next, the obtained optimal integration parameters are substituted into the iteration gain matrix to obtain the response of the structure at the next step; then, the initial conditions of the structure are updated and substituted into the discrete time domain motion equation for solution until the preset time is reached, and the response time history curve of the structure is output.

[0008] The specific steps are as follows:

[0009] Step 1: Set the integration parameters to obtain the discrete time domain solution of the motion equation.

[0010] (1) Set the integration parameters according to the structural dynamic equilibrium equation and obtain the gain matrix.

[0011] For a single-degree-of-freedom system, according to the dynamic equilibrium equation:

[0012]

[0013] where k, c, and m represent the stiffness, viscous damping, and mass of the structure respectively; f i+1 = f(t i+1 ) represents the external force, x i+1 is the displacement at the (i + 1)-th step. Correspondingly, is the velocity, is the acceleration. Set the discrete time domain solution of the motion equation as:

[0014]

[0015]

[0016] where Δt is the time interval; α 1,i 、α 2,iand β i are the optimal integral parameters obtained through optimization. The matrix form of the motion equation for discrete time domain solution is derived as follows:

[0017]

[0018]

[0019] I = {0 0 Δt 2 / m} T

[0020]

[0021]

[0022] where Ω = ωΔt; ξ = c / 2mω; is the natural vibration frequency of the structure, related to mass and stiffness; A is the gain matrix; and I is the force guiding vector.

[0023] For a multi-degree-of-freedom system:

[0024]

[0025]

[0026]

[0027] where M, C, and K are the mass, damping, and stiffness matrices of the multi-degree-of-freedom system respectively, and

[0028] i can be expressed in the following form:

[0029]

[0030]

[0031]

[0032] where j = 1, 2, 3, ……, n is the number of degrees of freedom.

[0033] Step 2: Determine the initial parameters and initial conditions of the structure.

[0034] (2.1) Determine the initial structural characteristic parameters (mainly mass m, stiffness k, and damping c), and at the same time determine a specific time step Δt.

[0035] (2.2) Obtain the initial state of the structure, including the initial displacement x 0 , the initial velocity v 0 and the initial acceleration a0 。

[0036] Step 3: Determination of optimization objectives.

[0037] (3.1) Define the optimization objectives as the minimum errors between the displacement, velocity, and acceleration and their expected values under free vibration conditions, which are respectively:

[0038]

[0039]

[0040]

[0041] (3.2) Use a linear weighted model to process the three optimization objectives in (3.1), and set adjustable weight coefficients according to the emphasis on different optimization objectives to obtain a comprehensive optimization objective function:

[0042]

[0043] where a, b, and c are weight coefficients, a + b + c = 1, and the specific values can be set according to the importance requirements of different optimization objectives.

[0044] Step 4: Set the specific characteristic requirements of the algorithm as the constraint conditions in the optimization process, mainly targeting the stability of the algorithm, numerical damping, period distortion, and the error from the exact solution

[0045] (4.1) To ensure the stability of the algorithm, the spectral radius of the gain matrix must be kept below 1, and the free vibration response must maintain a sinusoidal form.

[0046] The characteristic equation of the gain matrix A is:

[0047] |A - λE| = λ 3 -2A 1 λ 2 +A 2 λ - A 3

[0048] where E is the identity matrix and λ is the eigenvalue of matrix A; A 1 = 1 - 1 / 4Ω 2 α 2,i - ξΩβ i is half of the trace of matrix A, A 2 = -1 / 2Ω 2 α 2,i + Ω 2 α 1,i β i - 2ξΩβ i represents the sum of the cofactors of the main terms of matrix A, A 3= 0 is the determinant of matrix A; the roots of the characteristic equation are λ 1,2 = P ± Qi and the spurious roots λ 3 = 0.

[0049] Therefore, the stability constraint condition is:

[0050] ρ(A) = max|λ 1,2 | ≤ 1

[0051] After rearrangement, we get:

[0052]

[0053] Meanwhile, to ensure that the free vibration response of the structure is in the form of a sine wave, the roots λ of the characteristic equation 1,2 need to be complex numbers, and we get:

[0054]

[0055] (4.2) Period Error (PE) and Amplitude Decay (AD) have always been important indicators for measuring time-discretized integration algorithms. By setting the period error and numerical damping as constraints and setting different limits, the accuracy requirements of different structural dynamic calculations can be met.

[0056] AD and PE can be defined respectively as:

[0057]

[0058] where T = 2πΔt / Ω are the calculated period and the true period respectively; is the calculated frequency. Therefore, PE can be expressed as:

[0059]

[0060] For AD, it is commonly represented by the equivalent damping ratio ξ. After rearrangement, the constraint conditions related to period error and numerical damping are:

[0061]

[0062]

[0063] where P and Q represent the real part and the imaginary part of λ 1,2 respectively. It should be noted that the limit here can be set according to the accuracy requirements of structural dynamic calculations.

[0064] (4.3) Set the limits of the errors (including amplitude error and phase error) between the algorithm responses (displacement, velocity, and acceleration) and the exact solutions as the constraint conditions.

[0065] For a free vibration system, the algorithm displacement and the exact solution displacement can be expressed as follows:

[0066]

[0067]

[0068] The algorithm velocity and the exact solution velocity can be expressed as follows:

[0069]

[0070]

[0071] The algorithm acceleration and the exact solution acceleration can be expressed as follows:

[0072]

[0073]

[0074] Where,

[0075]

[0076]

[0077] In the above formula, t n = nΔt, x 1 and x 2 can be obtained by solving the motion equation through time-domain discretization. Assume that A and A e represent the algorithm amplitude and the exact amplitude respectively, and θ and θ e represent the algorithm phase and the exact phase respectively. The constraint conditions can be sorted out as follows:

[0078] A / A e -1 ≤ limit

[0079] θ - θ e ≤ limit

[0080] Where,

[0081] Step 5: Use an optimization algorithm to find the minimum value of the optimization objective under the constraint conditions.

[0082] (5.1) Adopt optimization algorithms such as Genetic Algorithm (GA), Particle Swarm Optimization (PSO), etc., and select the number of iterations and population size in the optimization algorithm according to the balance between calculation accuracy and calculation efficiency; at the same time, refer to other model-based integral algorithms to set the upper and lower limits of the optimization results and the search speed.

[0083] (5.2) Use the optimization algorithm to perform an optimization to obtain the optimal integral parameters α 1 、α 2 and β, and then obtain the corresponding gain matrix A.

[0084] Step Six: Perform iterative solution of the discrete-time domain motion equation to obtain the structural response.

[0085] (6.1) Substitute the gain matrix A obtained in (5.2) into the matrix form of the discrete-time domain motion equation, solve for the structural response, and obtain the next response.

[0086] (6.2) Define the structural response in (6.1) as the new initial state of the structure, update it, substitute it into the discrete-time domain motion equation for the next solution, until the pre-set time condition is met, terminate the iteration and output the structural response time history curve.

[0087] According to the above entire process, the optimal integral parameters based on the model can be obtained, forming an efficient time history analysis method system for linear complex structures based on optimized integral parameters.

[0088] The beneficial effects of the present invention: This method is a completely explicit integral algorithm with the characteristic of unconditional stability. Since the constraint conditions in the optimization process can be set by itself, the numerical damping of the algorithm in high-order modes is controllable; on the other hand, the accuracy of the algorithm can be changed by adjusting the limits in the constraint conditions. Therefore, compared with the traditional explicit integral algorithm, it has a wider scope of application; furthermore, the determined integral parameters and the diagonal form of the gain matrix make it have higher calculation efficiency in complex linear structures. In summary, an efficient time history analysis method for linear complex structures based on optimized integral parameters proposed by the present invention can demonstrate its advantages in systematic structural dynamics analysis, finite element analysis, and real-time substructure testing. Brief Description of the Drawings

[0089] Figure 1 is the basic process of the method proposed by the present invention.

[0090] Figures 2(a), 2(b), 2(c), 2(d), 2(e) and 2(f) are respectively the error comparisons of the response results of four integral algorithms for a single-degree-of-freedom undamped free vibration system under different Ω in the first embodiment; among them, Ω = 0.0548 in Figure 2(a), Ω = 0.1 in Figure 2(b), Ω = 0.2 in Figure 2(c), Ω = 0.5 in Figure 2(d), Ω = 0.8 in Figure 2(e), and Ω = 1 in Figure 2(f).

[0091] Figures 3(a) and 3(b) are respectively the displacement time history curves of the free vibration response of a linear single-degree-of-freedom system considering the influence of damping in the first embodiment; among them, Ω = 0.0548 in Figure 3(a) and Ω = 0.1 in Figure 3(b).

[0092] Figures 4(a) and 4(b) are respectively the error comparisons of different integral algorithms for a linear single-degree-of-freedom system considering the influence of damping in the first embodiment; among them, Ω = 0.0548 in Figure 4(a) and Ω = 0.1 in Figure 4(b).

[0093] Figure 5 is the transmission tower model described in the second embodiment.

[0094] Figures 6(a) and 6(b) are respectively the comparison of displacement time history curves obtained by using different integral algorithms for the transmission tower structure in the second embodiment; among them, PGA = 0.2g in Figure 6(a) and PGA = 0.8g in Figure 6(b). Detailed implementation manners

[0095] The following further illustrates the specific implementation manners of the integral algorithm in the present invention in combination with the accompanying drawings and technical solutions.

[0096] The main steps of the present invention are as Figure 1 shown. In the first embodiment, a single-degree-of-freedom system is selected to compare and analyze the accuracy of a high-efficiency time history analysis method for linear complex structures based on optimized integral parameters and other typical discrete-time integral algorithms (explicit Newmark method with β = 0 and γ = 0.5, Chang method, and CR method). In the second embodiment, taking a certain real transmission tower as an example, the linear time history analysis of the structure under seismic action is carried out to compare the calculation accuracy and calculation efficiency of different algorithms.

[0097] The specific implementation manner of the undamped system is as follows:

[0098] (1) Define the initial condition of the structure as x(0) = 1.25m, The parameters of the undamped single-degree-of-freedom system are m = 300 kg, k = 9000 N / m, and the damping ξ of the single-degree-of-freedom damped system is 0.1217.

[0099] (2) By changing the integration time step to select different Ω values, a high-efficiency time history analysis method for linear complex structures based on optimized integration parameters and three other typical explicit discrete-time integration algorithms are respectively used to obtain the structural responses of this system under free vibration.

[0100] (3) The response results obtained by different methods are compared with the exact solutions, and then the accuracies of different algorithms are compared. The absolute error results of the single-degree-of-freedom undamped system are shown in Figures 2(a), 2(b), 2(c), 2(d), 2(e) and 2(f); the response results of the single-degree-of-freedom damped system are shown in Figures 3(a) and 3(b), and the comparison of the absolute errors of different algorithms is shown in Figures 4(a) and 4(b).

[0101] (4) It can be seen from the results of Figures 2(a), 2(b), 2(c), 2(d), 2(e) and 2(f) that in the free vibration case of the single-degree-of-freedom undamped system, the calculation error of the explicit Newmark method is significantly greater than the algorithm proposed in the present invention. In addition, the CR method and the Chang method have similar numerical error characteristics, which are slightly larger than the algorithm proposed in the present invention in the low-order modes, and the gap is more obvious in the high-order modes. For the algorithm proposed in the present invention, although there are still errors in the high-order modes, they are much smaller than the other three typical methods. These results provide strong evidence for the superiority of the algorithm proposed in the present invention in terms of accuracy.

[0102] (5) It can be seen from Figures 3(a) and 3(b) that when Ω = 0.0548, the displacement time history curves of the four integration methods are basically the same. When Ω = 0.1, there are differences in the response amplitudes calculated by different integration methods; and Figures 4(a) and 4(b) show that when considering the influence of damping, the error of the algorithm proposed in the present invention is significantly smaller than that of other methods, further proving that the algorithm proposed in the present invention has the characteristics of high accuracy.

[0103] The specific implementation manner of the transmission tower structure is as follows:

[0104] (1) Use MATLAB R2010a software to establish a transmission tower model (as Figure 5 shown), this transmission tower is a lattice steel pipe structure, with a total height of 53.9m and a root opening of 8.32m, mainly composed of Q345 and Q420 steel pipes, and its first-order mode is calculated to be 1.70Hz; select the Northridge Beverly Hills-14145 Mulhol ground motion, adjust its intensity to 200gal and 800gal respectively, and input the ground motion along the Figure 5 shown X direction; Computer performance: Intel Core I i5-6402P CPU, 24GB memory.

[0105] (2) Without considering the nonlinearity of the structure, the top displacement responses of this transmission tower under seismic action are obtained by using an efficient time-history analysis method for linear complex structures based on optimized integration parameters and three other typical explicit discrete-time integration algorithms respectively, as shown in Figs. 6(a) and 6(b). Meanwhile, the computational time for the analysis and calculation is listed in Table 1.

[0106] (3) As shown in Figs. 6(a) and 6(b), the displacement time-history curves obtained by using different algorithms are basically the same, indicating the applicability of the algorithm proposed in the present invention in complex structures; and it can be seen from the computational time in Table 1 that the computational time of the algorithm proposed in the present invention is half of that of the Newmark method and much less than that of the Chang method and the CR method, indicating that the algorithm has high computational efficiency.

[0107] It can be seen from the comparative analysis that the algorithm proposed in the present invention can be applied to complex linear elastic structures and has higher computational accuracy and computational efficiency compared with other typical algorithms. Therefore, an efficient time-history analysis method for linear complex structures based on optimized integration parameters proposed in the present invention can demonstrate its advantages in the structural dynamics analysis, finite element analysis and real-time substructure test of amplitude linear systems. Table 1 shows the comparison of the computational time for the transmission tower structure described in Example 2 calculated by different integration algorithms.

[0108] Table 1 Comparison of the computational time for the transmission tower calculated by different algorithms

[0109]

Claims

1. An efficient time - history analysis method for linear complex structures based on optimized integration parameters, characterized in that, the steps are as follows: Step 1: Set the integration parameters to obtain the discrete - time - domain solution of the motion equation; (1) Set the integration parameters according to the structural dynamic equilibrium equation and obtain the gain matrix; For a single - degree - of - freedom system, according to the dynamic equilibrium equation: where k, c, and m represent the stiffness, viscous damping, and mass of the structure, respectively; f i+1 = f(t i+1 ) represents the external force, x i+1 is the displacement at the (i + 1)-th step. Correspondingly, is the velocity, is the acceleration. The discrete-time domain solution of the motion equation is set as: where Δt is the time interval; α 1,i , α 2,i and β i are the optimal integral parameters obtained by optimization; the matrix form of the motion equation for discrete-time domain solution is obtained as follows: I = {0 0 Δt 2 / m} T x i+1 = Ax i + If i+1 where Ω = ωΔt; ξ = c / 2mω; is the natural vibration frequency of the structure, related to mass and stiffness; A is the gain matrix; and I is the guiding vector of the force; For a multi - degree - of - freedom system: where \(M\), \(C\), and \(K\) are the mass, damping, and stiffness matrices of the multi-degree-of-freedom system, respectively, and are expressed in the following form: where j = 1, 2, 3, ……, n is the number of degrees of freedom; Step 2: Determine the initial parameters and initial conditions of the structure; (2.1) Determine the initial structural characteristic parameters, including mass m, stiffness k, and damping c, and at the same time determine the time step Δt; (2.2) Obtain the initial state of the structure, including the initial displacement x 0 , the initial velocity v 0 and the initial acceleration a 0 ; Step 3: Determine the optimization objective; (3.1) Define the optimization objective as the errors between the displacement, velocity, and acceleration under free - vibration conditions and the expected values, which are respectively: (3.2) Use a linear - weighted model to process the three optimization objectives in (3.1), set adjustable weight coefficients according to the emphasis on different optimization objectives, and obtain the comprehensive optimization objective function: where a, b, and c are weight coefficients, a + b + c = 1, and the specific values are set according to the importance requirements of different optimization objectives; Step 4: Set the specific characteristic requirements of the algorithm as the constraint conditions in the optimization process, mainly aiming at the stability of the algorithm, numerical damping, period distortion, and the error from the exact solution (4.1) To ensure the stability of the algorithm, the spectral radius of the gain matrix must be kept below 1, and the free - vibration response must remain in a sinusoidal form; The characteristic equation of the gain matrix A is: |A - λE| = λ 3 -2A 1 λ 2 +A 2 λ - A 3 where \(E\) is the identity matrix and \(\lambda\) is the eigenvalue of matrix \(A\); \(A\) 1 = 1 - 1 / 4\(\Omega\) 2 \(\alpha\) 2,i -\(\xi\)\(\Omega\)\(\beta\) i is half of the trace of matrix \(A\), \(A\) 2 = -1 / 2\(\Omega\) 2 \(\alpha\) 2,i +\(\Omega\) 2 \(\alpha\) 1,i \(\beta\) i - 2\(\xi\)\(\Omega\)\(\beta\) i represents the sum of the cofactors of the main terms of matrix \(A\), \(A\) 3 = 0 is the determinant of matrix \(A\); the roots of the characteristic equation are \(\lambda\) 1,2 = \(P\pm Qi\) and the pseudo root \(\lambda\) 3 = 0; Therefore, the stability constraint condition is: ρ(A) = max|λ 1,2 | ≤ 1 After arrangement: Meanwhile, to ensure that the free vibration response of the structure is in the form of a sine wave, the roots λ of the characteristic equation 1,2 need to be complex numbers, resulting in: (4.2) Set the period error PE and numerical damping AD as constraint conditions, and meet the accuracy requirements of different structural dynamic calculations by setting different limit values; AD and PE are respectively defined as: where T = 2πΔt / Ω are the calculation period and the true period respectively; is the calculation frequency. Therefore, PE is expressed as: For AD, it is represented by the equivalent damping ratio and after arrangement, the constraint conditions related to the periodic error and numerical damping are obtained as follows: where P and Q represent the real and imaginary parts of λ 1,2 respectively; limit is set by oneself according to the accuracy requirements of structural dynamic calculation; (4.3) Set the error limit between the algorithm response result and the exact solution as a constraint condition; among them, the algorithm response includes displacement, velocity, and acceleration, and the error includes amplitude error and phase error; For a free - vibration system, the algorithm displacement and the exact - solution displacement are respectively expressed as: The algorithm velocity and the exact - solution velocity are respectively expressed as: The algorithm acceleration and the exact - solution acceleration are respectively expressed as: where, In the above formula, t n = nΔt, x 1 and x 2 are obtained by solving the motion equation discretely in the time domain; assuming that A and A e represent the algorithm amplitude and the exact amplitude respectively, and θ and θ e represent the algorithm phase and the exact phase respectively, and the constraint conditions are obtained by sorting: A / A e -1 ≤ limit θ - θ e ≤ limit Among them, Step 5: Use an optimization algorithm to find the minimum value of the optimization objective under the constraint conditions; (5.1) Use an optimization algorithm, select the number of iterations and the population size in the optimization algorithm according to the balance of calculation accuracy and calculation efficiency; at the same time, set the upper and lower limits of the optimization result and the search speed; (5.2) Use the optimization algorithm to perform the first optimization to obtain the optimal integration parameters α 1 , α 2 and β, and then obtain the corresponding gain matrix A; Step 6: Perform iterative solution of the discrete - time - domain motion equation to obtain the structural response; (6.1) Substitute the gain matrix A obtained in (5.2) into the matrix form of the discrete - time - domain motion equation to solve the structural response and obtain the next - step response; (6.2) Define the structural response in (6.1) as the new initial state of the structure to update it, substitute it into the discrete - time - domain motion equation for the next - step solution, and terminate the iteration and output the structural response time - history curve until the pre - set time condition is met.

Images