粘弹性结构动力学分析中的一种数值方法彭 凡 傅 衣 铭(湖南大学工程力学系, 长沙 410082, 中国)摘要：针对材料具积分型本构关系，及松弛模量为Prony 级数形式的粘弹性结构动力学问题，本文结合Newmark 方法与Taylor 方法，建立了计算该类问题的一种数值算法。

且以简支梁为例，应用该方法具体地分析了考虑线性与非线性粘弹性时梁的强迫振动响应。

关键词：粘弹性 动力学 数值方法 响应1 引言随着人们对结构材料物理与力学性质了解的不断深入，以及新型材料的广泛应用，粘弹性结构的动力学研究受得了愈来愈多的重视，数值计算已成为一种主要的分析手段。

文[1]基于Newmark 方法建立了粘弹性结构动力学响应的有限元法，但只涉及到线性问题，而且在每一计算步，卷积积分的计算量较大。

桂洪斌等[2]提出将粘弹性结构的动力学方程进行Laplace 变换，然后在相域中求解问题，显然这种处理方式同样只适应于线性情况。

当考虑几何，物理包括损伤等非线性因素时，粘弹性结构动力学的数值分析就变得十分复杂与困难了。

文[3,4]通过将微分－积分型非线性动力学方程化成高阶的微分方程，最终由Runge-Kutta 法来获得数值解，但只有当材料为标准线性固体或Prony 级数取较少项数时，这种方法才比较容易实现。

本文针对材料服从积分型本构关系，且松弛模量为Prony 级数形式的粘弹性结构动力学问题，建立了从时域内直接求解的数值算法，它是基于Newmark 方法与Taylor 方法而得出的。

其中Taylor 方法为卷积积分的递规算法，能使计算量显著降低[5]。

文中通过对粘弹性梁的受迫振动分析来说明方法的应用。

2 简支粘弹性梁受迫振动的动力学方程考虑一简支梁，其跨度为L ，高为h ，中点受横向周期激励t H θsin 。

设材料具非线性粘弹性，可由Leaderman [6]本构关系描叙，则有00()()(())(())()tE t t E g t g d t τσεετττ∂-=+∂-⎰(1)式中)0(0E E =，)(t E 为松弛函数，)(εg 为应变ε的非线性函数：23()g εεβεγε=++ (2) 其中β与γ为常数。

在小挠度情况下，梁的受迫振动方程为：()3234522024223345224220(,)(,)(,) 1280()(,)(,) sin ()1280tw x t h w x t h w x t A E t x x x E t h w x h w x d H x L t t x x x ργτττγτδθτ⎡⎤⎛⎫∂∂∂∂⎢⎥++ ⎪∂∂∂∂⎢⎥⎝⎭⎣⎦⎡⎤⎛⎫∂-∂∂∂⎢⎥++=- ⎪∂-∂∂∂⎢⎥⎝⎭⎣⎦⎰ （3）式中A ,ρ分别为梁的质量密度及横截面面积，δ为Dirac 函数，满足两个简支端条件，即,,(0,)(0,)(,)(,)0xx xx w t w t w L t w L t ====的挠度),(t x w 取为1(,)()sink k k xw x t f t Lπ∞==∑ （4） 为说明问题起见，式(4)中只考虑1=k 的项，且令)()(1t f t f =。

将式(4)代入式(3)后，作Galerkin 积分，并记)()(E t E t D =，且取无量纲位移h t f t q )()(=，经运算和整理后得到：t H dq D dq D qθκωsin 032=*+*+ （5） 式中430212⎪⎭⎫⎝⎛=L A h E πρω，842809⎪⎭⎫⎝⎛=L h γπωκ，AhL HH ρ20=；⎰-∂-∂+=*td q t t D t q dq D 0)()()()(ττττ， ⎰-∂-∂+=*td q t t D t q dq D 0333 )()()()( ττττ当材料为线粘弹性时，简支梁的受迫振动方程为t H dq D qθωsin 02=*+ （6） 式(5)与式(6)中的)(t D 可表为Prony 级数的形式∑=-+=Kk t c k k e X X t D 10)( （7）其中材料参量0,,0>k k c X X ，且110=+∑=Kk k X X 。

3 数值算法对式(5)进行数值求解，设时间增量步为t ∆，基于Newmark [7] 方法，t t ∆+时刻的运动方程为)sin()()(032t t H dq D dq D qt t t t t t ∆+=*+*+∆+∆+∆+κω （8） 又()t t t t t t t qq t q q t q⎪⎭⎫ ⎝⎛--∆--∆=∆+∆+12111ααα （9）式中α和λ是按积分精度及稳定性要求而决定的参数。

再由Taylor [8]的卷积积分数值递规算法有())1(11)(+∆++∆++-=*n t t t n t t q q dq D ψμ （10）())3(133113])()[(+++∆++-=*n t n n tt q q dq D ψμ （11）式中 ∑=+++=Kk n k k n hX X 1)1(01μ，其中 )()1(1n k k t c n kh tc e hk =∆-=-+；(1)(1)1011 K n n t k k X q ψξ++==+∑(3)3(1)1031() Kn n tk k X q ψξ++==+∑ 其中 ()(1)()()11kc t n n n k k k k t t t X h q q e ξξ-∆+-∆⎡⎤=+-⎣⎦，()(1)()()3333()()kc t n n n k k k k t t t X h q q eξξ-∆+-∆⎡⎤=+-⎣⎦ 而t t q ∆-为t t ∆-时刻q 的值。

且当1=n ，即0=t 时，0)0(3)0(1==k k ξξ。

将式(9)~(12)代入式(8)后得到0)(3=++∆+∆+F Bq q A t t t t （12） 式中1+=n A κμ；120++=n B μωγ；232(1)(3)00231111sin()()t t t n t n t n n F H t t q q q q q γγγωμκμωψκψ++++⎡⎤=-+∆+++++--⎣⎦在计算步t ，A ,B 及F 已知，通过Newton 法解非线性方程(12)来求出t t ∆+时刻q 的值。

然后再反过来计算tt ∆+时刻q及q 的值 t t t t t t t q q q q q320)(γγγ---=∆+∆+，t t t t t t q q q q ∆+∆+++= 76γγ （13）这里 21t ∆=αγ，t ∆=αγ12，1213-=αγ，14-=αλγ，⎪⎭⎫ ⎝⎛-∆=225αλγt ，()λγ-∆=16t ，t ∆=λγ7。

对线粘弹性情形B F q t t '=∆+ (16)式中 )1(12123200)sin(++-++++∆+='n t n t t t q q q q t t H F ψωμωγγγ 。

以上算法的具体计算步骤可参考文献[7]，这里为节省篇幅，不再叙述。

4 算例及分析首先取材料为标准线性固体t c e X X t D 110)(-+=，此时，通过对式(6)微分及消去积分项的运算后，得到)cos sin ()1()1(310311112221t t c H q X c q X c q q q q c qθθθκωκω+=-+-++++ （17） 为对比结果，本文同时对式(18)用二级三阶Runge-Kutta 法进行计算。

算例中，取05.0 ,5.0110===c X X ，1=θ，25.00=H ，0)0(=q ，(0)0q =,05.0=∆t ，243293.0ωκ-=，45.0=α，65.0=λ，图1为8.0=ω时非线性粘弹性梁的受迫振动响应，图2为8.0=ω时线粘弹性梁的受迫振动响应，可以看出依本文方法得出的数值结果与依Runge-Kutta 法得出的结果基本吻合，对比表明无量纲振动幅值的最大相对差值约为6%。

图18.0=ω时非线性粘弹性梁的响应(K=1) 图2 8.0=ω时线性粘弹性梁的响应(K=1)图38.0=ω时非线性粘弹性梁的响应(K=5) 图4 8.0=ω时线性粘弹性梁的响应(K=5)取K=5的Prony 级数来表述松弛模量，其参数如表1所示表1 Prony级数中的各项参数分别考虑非线性与线性粘弹性时，梁的强迫振动响应计算示于图3与图4。

值得指出的是，本文的方法能被推广到含多种非线性因素的问题，以及多自由度系统的求解。

由于该方法是以Newmark方法为基础，具有较好的数值稳定性。

参考文献1吴琪泰. 粘弹性结构动力学响应的一个数值方法. 计算结构力学与应用, 1989, 6(4):21～272桂洪斌, 赵德有, 金咸定. 关于粘弹性结构和复合结构有限元动力学方程的探讨. 非线性动力学学报,2002, 9(1~2):1~53陈立群, 程昌均. 非线性粘弹性柱的稳定性和混沌运动. 应用数学与力学, 2000, 21(9):890~8964陈立群, 程昌均. 非线性粘弹性梁的动力学行为. 应用数学与力学, 2000, 21(9):897~9015Bradshaw RD, Brinson LC. Mechanical response of linear viscoelastic composite laminates incorporating nonisothermal physical aging effects. Comp. Sci. Tech., 1999, 59(8):1411~14276Leaderman H. Large longitudinal retareded elastic deformation of rubberlike network polymers. Polymer Trans Soc Rheol,1962,6(4):361~3827王瑁成, 劭敏. 有限单元法基本原理和数值方法. 北京: 清华大学出版社, 19978Taylor RL,Pister KS and Goudreau GL. Thermomechanical analysis of viscoelastic solids. Int. J. Numerical Method in Eng., 1970,2(1):45~51A NUMERICAL METHOD OF DYNAMIC ANALYSIS FORVISCOELASTIC STRUCTURESPeng Fan Fu Yiming(Dept. of Engineering Mechanics, Hunan University, Changsha, 410082, China) Abstract:Based on Newmark method for dynamic problems and Taylor algorithm for convolution integrals, a numerical method is constructed to treat the dynamic analysis of viscoelastic structures with hereditary constitutive relationship, and the relaxation modulus expressed in Prony series as well. To scheme the procedure of the method, the calculation of the responses of forced vibration for simply supported beams, which are modeled by Leaderman constitutive equation, is presented, and the validity of the numerical approach is verified by comparing the calculation results with ones obtained by applying Runge-Kutta method.Key Words: viscoelasticity; structural dynamics; numerical method; response。