复阻尼结构动力方程的增维精细积分法

吴泽玉,王东炜,李玉河

(1.华北水利水电大学 土木与交通学院，郑州 450045；2 郑州大学 建工学院，郑州 450000 )

复阻尼结构动力方程的增维精细积分法

吴泽玉1,王东炜2,李玉河1

(1.华北水利水电大学 土木与交通学院，郑州 450045；2 郑州大学 建工学院，郑州 450000 )

为了避开求解复阻尼结构强迫激励动力学方程的积分运算，引入增维精细积分方法。根据复阻尼系统复化对偶原则，将动力学方程和激励波对偶复化为实部和虚部的形式，推导出增维矩阵的精细积分求解过程。结果表明，由于不用求解迭代矩阵H的逆矩阵，避免了矩阵奇异带来的计算解的不稳定性。在计算矩阵仅增加一维的情况下，化积分运算为代数运算，扩大了精细积分法的应用范围。通过对比增维精细积分法和频域法计算结果，二者结果保持较高的一致性。

复阻尼；动力时程分析；复化对偶；材料损耗因子；增维精细积分法

阻尼是结构系统的固有特征参量，是计算结构动力响应的重要参数。在动力时程分析时，阻尼具有消减峰值、降低结构响应的作用；同时，阻尼也是一个十分复杂的问题。在众多的阻尼理论中，以黏滞阻尼模型和复阻尼模型为主。黏滞阻尼理论每周能量消耗与激励荷载频率相关，这与实验现象矛盾[1]。复阻尼模型不仅引入了材料损耗因子的概念，而且是在大量固体材料实验基础上总结得到，所以能更好的表征结构的阻尼性质[2-5]。在求解复阻尼动力学方程时应遵循对偶复化原则，将激励荷载和结构响应均表示成复数形式进行分析[6-9]。

随着动力时程分析理论发展和计算机技术的巨大进步，动力时程分析方法在工程中得到广泛应用。结构动力时程计算方法主要包括Wilson-θ和Newmark-β法等。近年来，因精细时程分析法具有计算精度高、计算步长大和计算结果稳定等优点在时程分析中得到广泛关注[10-11]。然而用精细积分法求动力时程分析强迫振动问题时需通过积分求得方程解，求解过程较为复杂。利用增维精细积分法求解复阻尼动力方程，化积分运算为代数运算，仅以增加一维计算矩阵求结构时程响应，尚未见相关文献阐述。

本文根据复阻尼结构系统复化对偶基本原则，将离散激励荷载表示成实部加虚部的形式；进而推导出多自由度体系增维矩阵精细积分状态方程。最后通过计算例题，计算结果同频域法(视为精确解)进行对比分析，验证增维精细积分法的计算精度和数值稳定性。

1 激励向量的复化对偶分析

复阻尼多自由度系统动力学方程式为：

(1)

对于具有N个偶数离散激励向量，数据可表示为：

aj=a(jΔt),j=0,1,2,…,N

(2)

此离散向量可用三角函数插值表示为：

AN/2/2cosθN/2t

(3)

其中参量θn、An和Bn计算式分别为：

θn=2πn/(NΔt),n=1,2,…,N/2

(4)

(5)

(6)

对于式(1)中f2(t)的具体表达形式，根据对偶复化要求，f1(t)中A0/2的对偶项为ηA0/2；Ancosθnt的对偶项为Ansinθnt；Bnsinθnt的对偶项为-Bncosθnt。故f2(t)的三角函数表达式为：

AN/2/2sinθN/2t

(7)

2 复阻尼结构动力方程增维积分公式推导

在复荷载作用下结构动力响应也由实部和虚部组成，式(1)结构位移响应可表示为：

X=X1+iX2

(8)

将式(8)代入到式(1)中，可得到两个n维实数动力学方程式为：

(9)

将式(9)表示成矩阵形式：

(10)

令：

(11)

将式(11)写成状态方程形式：

(12)

式中：I为主对角元素为1的方阵。

式(12)为4n维状态空间方程，分别由2n维实数向量和2n维虚数向量组成。当方程为非齐次方程时求解需做积分运算，求解过程较为复杂。为了简化计算，引入常数等式a(t)=1，化积分运算为代数运算。式(12)可表示为：

(13)

令：

则式(12)可表示为：

(14)

式(14)的解为：

Y(tk+1)=T(τ)Y(tk)

(15)

式中，T(τ)的表达式为：

T(τ)=exp(Hτ)

(16)

其中τ为离散时间步长。

式(13)比式(12)增加了一维向量，即(4n+1)维，即可化积分运算为代数运算，简化了求解过程。同时，方程求解是通过不同时刻矩阵H和位移与速度确定，避开了对矩阵H的求逆运算造成的计算不稳定性问题，拓宽了精细积分法求解动力方程的应用范围。

3 指数矩阵T(τ)的求解

利用增维精细积分法求解式(15)的关键是T(τ)矩阵的计算。文献[12]提出了十九种指数矩阵计算方法，作者认为都不能取得可靠解。文献[10-11]创造性提出采用加法定理计算指数矩阵，计算精度可媲美精确解，计算过程如下：

T(τ)=eHτ=[eHτ/m]m

(17)

其中可选用m=2N，通常令N=20，则m=1 048 576。一般情况下荷载步τ值比较小，故时间段Δt=τ/m将是非常小数值，则在Δt时间段内有

eHΔt≈I+HΔT+(HΔT)2/2=I+Ta

(18)

从而

T=(I+Ta)2N=[(I+Ta)(I+Ta)]2N-1=

[(I+2Ta+Ta×Ta)]2N-1

(19)

式(19)相当于计算机程序语句：

For(iter=0;iter

T=I+Ta

(20)

4 增维精细积分求解过程

复阻尼结构动力学增维精细积分法求解过程如下：

(1)确定结构阻尼比或材料耗散因子，通常取η=2ζ，ζ表示结构阻尼比；然后将离散激励点表示成f1(t)和其对偶项f2(t)的形式。

(2)将复阻尼动力学式(1)表示成式(12)的形式，引入常数量a(t)=1进一步将式(12)表示为式(13)。

(3)计算每一步荷载作用下指数矩阵T(τ)的值，进而由式(15)计算每一荷载步的位移值和速度值。

(4)将位移值和速度值代入式(1)，可求出每一步荷载作用下结构的加速度值。

(5)位移和速度的实部与虚部对应关系为：

位移：实部第(i)项⟹虚部第(2n+i)项;速度：实部第(i+n)项⟹虚部第(3n+i)项且i≤n，n为动力学式(1)实部的维数。例如对于两个自由度体系结构，第1自由度实部位移和第3项虚部位移为对偶项；第1自由度实部速度和第7项虚部速度为对偶项。

5 计算实例

某三个自由度体系，质量矩阵为M=diag(2,1,1)，刚度矩阵为K=[6 -2 0；-2 4 -2；0 -2 2]。激励波选用EL-Centro地震波，峰值加速度调整为0.2g，经复化对偶实部和虚部波及原地震波如图1所示。不计材料损耗因子随应力的变化[1]，取为η=2ζ=0.10，增维精细积分法(记作MM)、高斯精细积分法(记作GM)和频域法(记作FM，视为精确解)计算结果进行对比。为了节省篇幅，这里仅给第一和第二自由度位移和速度对比图：第一、二自由度位移如图2所示；第一、二自由度速度如图3所示。

图1 EL-Centro波及复化对偶波Fig.1 EL-Centro wave and complex dual wave

图2 第一、二自由度位移对比图Fig.2 Displacement comparison of first and second freedom

图3 第一、二自由度速度对比图Fig.3 Velocity comparison of first and second freedom

从前两个自由度位移图和速度图可以看出，实部和虚部基本保持了对偶关系，二者振动幅值相等；相位角方面，实部峰值滞后虚部峰值大约π/4左右。增维精细积分法计算结果具有足够的稳定性，没有出现发散现象；计算精度和高斯精细积分法及精确解(即频域解)不论在实部解还是虚部解都保持较高的一致性。

6 结 论

本文研究了复阻尼系统动力学方程的时间历程求解。根据复阻尼系统荷载复化对偶原则，利用指数矩阵加法原则和精细积分方法，化积分运算为代数运算，推导出增维矩阵求解状态空间的新方法。计算结果表明，利用增维精细积分法求解动力时程问题，计算精度同精确解高度一致；由于不用求解迭代矩阵H的逆矩阵，避免了矩阵奇异带来的计算解的不稳定性；在计算矩阵仅增加一维数的情况下，化积分运算为代数运算，扩大了精细积分法的应用范围。

[1] LAZAN B J. Damping of material and members in structural mechanics [M]. London: Pergamon Press, 1968.

[2] GADE S，ZAVERI K， KONSTANTIN-HANSEN H，et al. Complex modulus and damping measurements using resonant and non-resonant methods[C]// SAE Technical Papers, 1995.

[3] CORDIOLI J A, BRATTI G, STUMPF C, et al. On the prediction of damping loss factor of fuselage panels with viscoelastic materials using periodic structure theory and finite element method[C]// Proceedings of ISMA,2010: 2257-2266.

[4] KOSMATKA J B, LIGUORE S L. Review of methods for analyzing constrained-layer damped structures[J]. Journal of Aerospace Engineering, 1993,6(3):268-283.

[5] 朱镜清. 频率相关黏性阻尼理论及有关问题的解[J]. 振动与冲击, 1992,21(4): 1-7. ZHU Jingqing. Frequency dependant viscous damping theory and some related problems[J]. Journal of Vibration and Shock, 1992,21(4): 1-7.

[6] 何钟怡. 复本构理论中的对偶原则[J]. 固体力学学报, 1994, 15(2): 177-180. HE Zhongyi. The dual principle in theory of complex constitutive equations[J]. Acta Mechnica Solid Sinica, 1994, 15(2):177-180.

[7] 王建平，刘宏昭. 复阻尼振动系统的对偶原则及其数值方法研究[J]. 振动工程学报，2004，17(1): 62-65. WANG Jianping, LIU Hongzhao. Dual principle of oscillation systems with complex damping and its numerical method[J]. Journal of Vibration Engineering, 2004, 17(1): 62-65.

[8] 张辉东，王元丰. 复阻尼模型结构地震时程响应研究[J]. 工程力学，2010, 27(1): 109-115. ZHANG Huidong, WANG Yuanfeng. Study on seismic time-history response of structures with complex damping[J]. Engineering Mechanics, 2010, 27(1): 109-115.

[9] 潘玉华，王元丰. 复阻尼结构动力方程的高斯精细时程积分法[J]. 工程力学，2012，29(2): 16-20. PAN Yuhua, WANG Yuanfeng. Gauss precise time integration of complex damping vibration systems[J]. Engineering Mechanics, 2012, 29(2): 16-20.

[10] 钟万勰. 结构动力方程的精细时程积分法[J]. 大连理工大学学报，1994, 34(2): 131-136. ZHONG Wanxie. On precise time integration method for structural dynamics[J]. Journal of Dalian University of Technology, 1994, 34(2): 131-136.

[11] 钟万勰. 暂态历程的精细计算方法[J]. 计算结构力学及其应用，1995, 12(1): 1-6. ZHONG Wanxie. Precise computation for transient analysis[J]. Computational Structural Mechanics and Applications, 1995, 12(1): 1-6.

[12] MOLER C B, VON LOAN C F. Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later[J]. Society for Industrial and Applied Mechanics,2003,45(1): 3-49.

Magnified dimension precise integration method for the dynamic equations of complex damped structures

WU Zeyu1, WANG Dongwei2, LI Yuhe1

(1 School of Civil Engineering and Communication, North China University of Water Resources and Electric Power, Zhengzhou 450045, China; 2 Civil Engineering Department, Zhengzhou University, Zhengzhou 450000, China)

In order to avoid the integral operation of the forced excitation dynamics equation of complex damped structures, a magnified precise integration method was introduced. According to the dual principle of complex damping system, the dual dynamic equations and the dual excitation waves were divided into real and imaginary parts, and the solving process by using the precise integration of the augmented matrix was derived. The results show that because the inverse matrix of the iteration matrixHdoesn’t need to be solved, the instability of the computational solution caused by the singularity of matrix is avoided. In the calculation of the matrix, only a one-dimensional integral calculation is increased and the integral operation is transformed into algebraic operations, so, the scope of application of the precise integration method can be expanded. The comparison between the results calculated by using the augmented precise integration method and the frequency domain method shows they are in good consistency.

complex damping; time history analysis; complex duality; material loss factor; magnified dimension precise integration method

国家自然科学基金(50978232);河南教育厅项目(14B560029);华北水利水电大学高层次人才基金(201246)

2016-03-14 修改稿收到日期：2016-05-05

吴泽玉 男，博士，讲师，1976年9月生 E-mail:wuzeyu1976@163.com

O175.3

A

10.13465/j.cnki.jvs.2017.02.017