復阻尼結構動力方程的增維精細積分法

吳澤玉,王東煒,李玉河

(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