朱媛媛 李瑩 程昌鈞

(1．上海師范大學計算機科學與技術(shù)系，上海 2 00234)(2．上海大學力學系，上海 2 00444)

孔隙熱彈性梁的非線性氣彈性特性*

朱媛媛1?李瑩2程昌鈞2

(1．上海師范大學計算機科學與技術(shù)系，上海 2 00234)(2．上海大學力學系，上海 2 00444)

根據(jù)作者由Hamilton變分原理導出的一個孔隙熱彈性梁的非線性數(shù)學模型和氣彈性原理中的一階修正線性活塞理論，本文首先給出了位于高速或者超高速流動中兩端固定的平面孔隙熱彈性梁的控制微分方程和定解條件，其中基本未知量是梁的軸向和橫向位移以及孔隙百分比和溫度變化引起的“力矩”．為了考察孔隙熱彈性梁在橫向載荷和氣彈性載荷聯(lián)合作用下的非線性力學特性，采用微分求積方法對問題進行空間離散，得到一組關(guān)于時間的非線性常微分方程，然后在給定初始條件下采用變步長Runge－Kutta方法對方程組進行數(shù)值求解，由此研究了孔隙熱彈性梁的氣彈性特性，考察了參數(shù)的影響，得到了一些有益的結(jié)論．

孔隙熱彈性梁， 有限變形， 線性活塞理論， 微分求積方法， 氣彈性特性

引言

孔隙熱彈性材料兼具結(jié)構(gòu)和功能雙重用途，多見于天然多孔材料、人造多孔材料和生物工程材料等．這種材料具有相對密度低、比強度高、比表面積高、重量輕、隔音、隔熱、滲透性好等優(yōu)點，在航空、航天、化工、建材、冶金、原子能、石化、機械、醫(yī)藥和環(huán)保等領(lǐng)域具有廣泛的應用前景．同時，新型孔隙材料正逐漸地被用作絕緣、緩沖、吸收沖擊能量的材料，從而發(fā)揮了其由多孔結(jié)構(gòu)決定的獨特的綜合性能．因此，對孔隙熱彈性材料理論的研究引起了很多學者的興趣．主要研究內(nèi)容包括兩個方面，其一是孔隙熱彈性材料性質(zhì)的合理描述與本構(gòu)關(guān)系的表達，其二是孔隙熱彈性材料和結(jié)構(gòu)的相關(guān)問題的求解．經(jīng)過Goodman，Cowin和Iesan等眾多學者的發(fā)展，已經(jīng)形成了一個較為統(tǒng)一的孔隙熱彈性理論［1－4］．

作為最基礎(chǔ)的結(jié)構(gòu)元件，梁在建筑結(jié)構(gòu)、橋梁結(jié)構(gòu)、水工結(jié)構(gòu)、地下結(jié)構(gòu)、飛機結(jié)構(gòu)、船舶結(jié)構(gòu)、海洋工程結(jié)構(gòu)、空間工程以及微機電系統(tǒng)等許多領(lǐng)域具有廣泛應用．文獻［5］在有限變形的條件下，把Hamilton變分原理推廣到孔隙熱彈性梁中，由此給出了孔隙熱彈性梁的一個以中性軸的3個位移、孔隙百分比和溫度引起的“力矩”所表示的非線性數(shù)學模型，并考慮了軸力、中性層慣性和轉(zhuǎn)動慣性的影響．本文以此模型為基礎(chǔ)，研究了兩端固定的孔隙熱彈性平面梁在氣動力作用下的非線性力學行為，其中氣動力采用修正的一階線性活塞理論來模擬．由于問題的非線性性，難以得到其解析解或者半解析解，因此采用微分求積方法(Differential Quadrature Method，簡稱DQM)對數(shù)學模型進行空間離散，再采用Runge－Kutta方法對離散化的動力系統(tǒng)進行數(shù)值求解，在此基礎(chǔ)上研究了孔隙熱彈性梁的氣彈性特性，得到了一些有益的結(jié)論．

1 問題的數(shù)學描述

由Iesan［4］提出的孔隙熱彈性材料的動力學理論，有如下孔隙的運動微分方程、熵均衡方程和各向同性孔隙熱彈性梁的本構(gòu)方程

這里，ρ為參考構(gòu)形的密度，φ為基體材料體積百分比的變化，χ和hi是平衡慣量和平衡應力向量，g和ρl是內(nèi)外平衡體積力，η和S分別是單位質(zhì)量的熵和熱源，qi是熱流矢量，σx為梁的應力分量，θ=T－T0表示溫度的變化，其中T是物體的絕對溫度，T0是參考構(gòu)形的絕對溫度．D0=E(1－v)/［(1+v)(1－2v)］，β=αtE/(1－2v)，E，v分別是彈性模量和泊松比，αt是線性熱膨脹系數(shù)，K是熱傳導系數(shù)，ce是比熱，同時 αv，bv，ξv，mv是與孔隙有關(guān)的材料參數(shù)，它們的物理意義可在文獻［4］中找到．

現(xiàn)在考察圖1所示長度為l，寬度為b，厚度為h的孔隙熱彈性平面梁．取坐標原點o與梁一端的截面的慣性主軸的交點重合，ox－軸為梁的中性軸，oy－軸和oz－軸為截面的慣性主軸，并組成右手坐標系．設(shè)中性層的位移分別為u，v，w，基于 Kirchhoff－Love假設(shè)，當撓度w與寬度b和厚度h同階，并且梁的長度遠大于寬度和厚度時，有非線性幾何關(guān)系

其中，εx為應變分量，εx0為中性軸的應變分量．

圖1 梁的模型Fig．1 The model of the beam

其中，ρh¨u和ρIy¨w，xx分別是中性層的慣性項和轉(zhuǎn)動慣性項，p是橫向的分布保守力，N0是端部軸力．另外，根據(jù)文獻［5］，qφ和qθ是兩個分別與Mφ和Mθ有關(guān)的量，被給定為

邊界條件:對于x=0，l兩端固定的不受軸力N0作用的孔隙熱彈性梁，有如下邊界條件:

初始條件:假設(shè)初始時刻t=0時，有初始條件

其中，u0，w0，，，都是x的已知函數(shù)，特別當初始時刻梁處于靜止時，這些函數(shù)為零．

2 孔隙熱彈性梁非線性氣彈性分析的控制微分方程

處于高速和超高速氣流中的彈性結(jié)構(gòu)往往會發(fā)生顫振，即所謂的氣動彈性問題，這種現(xiàn)象是由于空氣動力及其引起的彈性結(jié)構(gòu)反作用力之間的相互影響而產(chǎn)生的．飛行器以及處于自由風速流動中的高層結(jié)構(gòu)容易會發(fā)生這種問題．氣動彈性力學大致可分成兩大分支:一是飛行器氣動彈性力學，主要研究飛機、直升機、導彈、火箭直至航天飛機和衛(wèi)星等飛行器的氣動彈性問題，這方面的知識已成為設(shè)計和制造先進的飛行器所必備的基礎(chǔ);另一是建筑結(jié)構(gòu)氣動彈性力學，主要研究大型高層建筑、橋梁、高空電纜、龍門吊車以及細長的船舶桅桿、電視尖塔、煙囪等各種具有彈性性質(zhì)的建筑物的氣動彈性問題，是近年來蓬勃發(fā)展的一個學科分支．由于這些建筑結(jié)構(gòu)幾乎與國民經(jīng)濟的所有部門有關(guān)，因此保證它們的安全和工作可靠具有極其重要的意義［6－7］．

氣動力計算是氣動彈性分析的重要組成部分，活塞理論由于形式簡單，且具有較好的精度，被廣泛應用于計算分析中，本文中采用修正的一階線性活塞理論來模擬氣動壓載荷［7］

引入如下的無量綱變量和參數(shù):

將(9)代入方程(5)～(8)中，則有如下無量綱的控制微分方程和定解條件:

其中，θ+，θ－是上下表面的溫度．Δp為氣動力的無量綱化，即

邊界條件無量綱化為

假設(shè)初始時刻靜止，則有無量綱初始條件為

3 求解方法與數(shù)值解的驗證

現(xiàn)在，采用微分求積方法(Differential Quadrature Method簡稱DQM)對控制方程(10)～(11)和定解條件(12)～(13)進行空間離散．DQM是由Bellman等人［8－9］于 1971 ～1972 年提出的一種數(shù)值求解方法，現(xiàn)已廣泛應用于各個領(lǐng)域［10－12］．其基本概念是將函數(shù)對某方向的自變量的偏導數(shù)近似表達為沿該自變量方向各離散點(節(jié)點)處相應函數(shù)值的加權(quán)和，其加權(quán)系數(shù)取決于試函數(shù)和節(jié)點分布，而與具體問題無關(guān)．因此，利用這些系數(shù)，任何偏微分方程都能轉(zhuǎn)化為相應的常微分方程或者代數(shù)方程，然后再求解．根據(jù) DQM的概念和公式［10－11］，可得無量綱控制微分方程(10)的 DQ 離散化形式為:其中，Δpi為無量綱氣動力(11)的DQ離散化，被給定為

在(14)～(15)式中，i=1，2，…，N，N為所布置的結(jié)點數(shù)表示對x的j階偏導數(shù)的權(quán)系數(shù)，并且上面具有圓點的量表示對時間τ的導數(shù)．本文采用多項式作為試函數(shù)來構(gòu)造權(quán)系數(shù)，并利用Chebyshev－Lobatto多項式的零點的坐標來布置節(jié)點［10－11］．

邊界條件(12)的DQ離散化形式為

同時，初始條件(13)的DQ離散化

這樣，對于孔隙熱彈性平面梁的氣動力學問題，就是在給定邊界條件(16)和初始條件(17)條件下，求解非線性常微分方程組(14～(15)．

圖2 DQ數(shù)值解與解析解的對比(a)=10－5;(b)=10－4Fig．2 Comparison between the numerical solutions and the analytical solutions(a)=10－5;(b)=10－4

為了說明本文方法的正確性，考慮受均布力作用兩端簡支彈性梁的彎曲問題．在小變形的情況下可得撓度的解析解．圖2(a)示出了均布力p=10－5，布點數(shù)為N=7時，所得DQ數(shù)值解和解析解的比較，兩者非常吻合．圖2(b)示出了均布力p=10－4，布點數(shù)為N=7，9時，DQ 數(shù) 值解和解析解的比較．可以看到，當均布力增大時，線彈性小變形理論不再適用，由于非線性的影響，DQ數(shù)值解略小于解析解．由圖2(b)還看到，DQ數(shù)值解的收斂性很好．

4 孔隙熱彈性梁的DQ數(shù)值結(jié)果與氣彈性特性分析

現(xiàn)在考察橫向動載荷作用和氣動力同時作用下孔隙熱彈性梁的動力學特性，按照Chebyshev零點在梁的中性軸上布置N=15節(jié)點．并取材料參數(shù)為［13－14］:β1=0．1，h=0．1，a3=0．32，a4=12．0，a5=0．75，a6=0．6，a7=0，a8=0．5，a9=0．0804，a10=1．0 ×10－8，a11=4．0 ×106．

4．1 梁的動力學特性

首先考察兩種不同類型的簡單動載荷，即(1)Heaviside載荷:(a)(τ)= － 1．0 ×10－4，τ≥0;(b)(τ)= －1．0 ×10－3，τ≥0．(2)正弦載荷:(a)(τ)= －1．0 ×10－3× s in(2πτ)，τ≥0;(b)(τ)= －1．2 ×10－2×sin(2πτ)，τ≥0．圖3～圖4和圖5～圖6分別示出了在Heaviside載荷和正弦載荷作用下梁中點處變量的時程曲線和相圖．可見，動載荷的類型對梁的力學行為的影響是明顯的．

計算中，我們也研究了四種材料梁(即孔隙熱彈性梁(TEVB)、無孔隙熱彈性梁(TEB)、孔隙彈性梁(EVB)和無孔隙彈性梁(EB))的動力學特性和(7)中參數(shù)的影響，限于篇幅不再一一列出．總的結(jié)論是:孔隙的存在使梁的撓度增加．特別是位移－孔隙耦合參數(shù)主要引起撓度和孔隙力矩的變化，對溫度力矩的影響不太明顯;位移－熱耦合參數(shù)主要引起撓度和溫度力矩的變化，對孔隙力矩也有一定的影響;孔隙－熱耦合參數(shù)主要引起孔隙力矩和溫度力矩的變化，對撓度影響較?。旱膿隙入S參數(shù)a1，a3，a5的增加而增加，但隨參數(shù)a2，a4，a8的增加而減小;其他參數(shù)對梁的撓度影響均比較小．在動力學情形下，由于外力的功部分轉(zhuǎn)化為熱能被耗散，溫度對撓度和孔隙力矩的影響不像靜力學情形那樣明顯．

圖3 梁中點處變量的時程曲線和相圖()= －1．0 ×10 －4，τ≥0)(a)撓度的時程曲線;(b)撓度的相圖;(c)孔隙百分比力矩時程曲線;(d)溫度力矩時程曲線Fig．3 The time－h(huán)istory curves and the phase portrait at the midpoint()= －1．0 ×10 －4，τ≥0)(a)The time－h(huán)istory curves of the deflection;(b)The phase portrait of the deflection;(c)The time－h(huán)istory curves of the void moment;(d)The time－h(huán)istory curves of the temperature moment

圖4 梁中點處變量的時程曲線和相圖)= －1．0 ×10 －3，τ≥0)(a)撓度的時程曲線;(b)撓度的相圖;(c)孔隙百分比力矩時程曲線;(d)溫度力矩時程曲線Fig．4 The time－h(huán)istory curves and the phase portrait at the midpoint= －1．0 ×10 －3，τ≥0)(a)The time－h(huán)istory curves of the deflection;(b)The phase portrait of the deflection;(c)The time－h(huán)istory curves of the void moment;(d)The time－h(huán)istory curves of the temperature moment

圖5 梁中點處變量的時程曲線和相圖= －1．0 ×10－3×sin(2πτ)，τ≥0)(a)時程曲線;(b)相圖Fig．5 The time－h(huán)istory curves and the phase portrait at the midpoint(p(τ)= －1．0 ×10－3×sin(2πτ)，τ≥0)(a)The time－h(huán)istory curves;(b)The phase portrait

圖6 梁中點處變量的時程曲線和相圖= －1．2 ×10－2×sin(2πτ)，τ≥0)(a)時程曲線;(b)相圖Fig．6 The time－h(huán)istory curves and the phase portrait at the midpoint()= －1．2 ×10－2×sin(2πτ)，τ≥0)(a)The time－h(huán)istory curves;(b)The phase portrait

4．2 梁的氣彈性特性

圖7 撓度的動力學特性(Q∞=2．365×10－3)(a)時程曲線;(b)功率譜;(c)相圖;(d)Poincare截面Fig．7 Aerodynamic characteristics of the deflection(Q∞ =2．365×10－3)(a)Time－h(huán)istory curve;(b)Power spectrum;(c)Phase portrait;(d)Poincar＇e map

圖8 撓度的動力學特性(Q∞=4．063×10－3)(a)時程曲線;(b)功率譜;(c)相圖;(d)Poincare截面Fig．8 Aerodynamic characteristics of the deflection(Q∞ =4．063×10－3)(a)Time－h(huán)istory curve;(b)Power spectrum;(c)Phase portrait;(d)Poincar＇e map

圖9 撓度的動力學特性(Q∞=4．096×10－3)(a)時程曲線;(b)功率譜;(c)相圖;(d)Poincare截面Fig．9 Aerodynamic characteristics of the deflection(Q∞ =4．096×10－3)(a)Time－h(huán)istory curve;(b)Power spectrum;(c)Phase portrait;(d)Poincar＇e map

給定= －10－4× sin(2πτ)，軸力N0=10－4，圖7～圖10給出了不同氣動壓Q∞條件下?lián)隙鹊臅r程曲線、功率譜、相圖和Poincare映射．可見，當自由流的動壓力逐漸增強時，Poincare映射點集由落在一個點上演變成許多個無規(guī)則的點，而且數(shù)值越來越大，結(jié)合相應的相圖知，系統(tǒng)由周期運動轉(zhuǎn)化為混沌運動．當氣動壓超過4．096×10－3時，點集并不被吸引到一個有界區(qū)域，而是發(fā)散到無窮．

圖11給出了Q∞=4．007 ×10－3，梁中點的初始撓度分別為w0= －0．01，0．01時的時程曲線，可見，該狀態(tài)由于存在多個漸進穩(wěn)定的平衡點，即使在相同的系統(tǒng)參數(shù)下，當初值不同時，軌線趨于不同的穩(wěn)定的平衡點，這是由于幾何非線性所致．

圖10 撓度的動力學特性(Q∞=4．121×10－3)(a)時程曲線;(b)功率譜;(c)相圖;(d)Poincare截面Fig．10 Aerodynamic characteristics of the deflection(Q∞ =4．121×10－3)(a)Time－h(huán)istory curve;(b)Power spectrum;(c)Phase portrait;(d)Poincar＇e map

圖11 不同初始撓度下?lián)隙鹊臅r程曲線(a)w0=－0．01;(b)w0=0．01Fig．11 The time－h(huán)istory curves of the deflection for different initial values(a)w0= －0．01;(b)w0=0．01

5 結(jié)論

在本文中，首先根據(jù)文獻［5］中通過Hamilton變分原理建立的有限變形孔隙熱彈性梁的數(shù)學模型，并采用一階修正線性活塞理論給出的氣動力來模擬結(jié)構(gòu)在超高速流動中受到的氣彈性載荷，由此給出了孔隙熱彈性平面固定梁在動載荷和氣動力作用下的非線性數(shù)學模型．在此基礎(chǔ)上，采用DQM對變量進行空間離散，然后采用Runge－Kutta方法對離散化動力學問題進行數(shù)值求解．為了說明本文理論和方法的正確性，研究了兩端簡支彈性梁的小變形平面彎曲靜態(tài)彎曲問題，并與解析解進行了比較．可以看到，由DQM得到的數(shù)值解和解析解是非常吻合的，且DQ數(shù)值解的收斂性也非常好．詳細考察了氣動載荷作用下孔隙熱彈性固定梁的氣彈性行為．發(fā)現(xiàn)在氣動力的作用下，結(jié)構(gòu)具有周期的穩(wěn)定運動和有界混沌、無界發(fā)散兩種不穩(wěn)定的運動狀態(tài)．與無氣動力作用的孔隙熱彈性梁相比，在氣動力作用下孔隙熱彈性梁更容易發(fā)生無界不穩(wěn)定的現(xiàn)象．同時，增加自由流動壓Q∞會降低孔隙熱彈性梁的穩(wěn)定性．

1 Goodman M A，Cowin S C．A continuum theory for granular materials．Archive for Rational Mechanics and Analysis，1972，44(4):249 ～266

2 Nunziato J W，Cowin S C．Nonlinear theory of elastic materials with voids．Archive for Rational Mechanics and Analy－sis，1979，72(2):175 ～201

3 Cowin S C，Nunziato J W．Linear elastic materials with voids．Journal of Elasticity，1983，13(2):125 ～147

4 Iesan D．A theory of thermoelastic materials with voids．Acta Mechanica，1986，60(1－2):67～89

5 Li Y，Cheng C J．A nonlinear model of thermoelastic beams with voids and its applications．Journal of Mechanics of Materials and Structures，2010，5(5):805 ～820

6 伏欣．氣動彈性力學原理．沈克揚譯．上海:上?？茖W技術(shù)文獻出版社，1982(Forsching H W．Grundlagen der aeroelastik． S pringer－Verlag， B erlin， H eidelberg， N ew York，1974)

7 Mei C，AbdeI－Motagaly K，Chen R．Review of nonlinear panel flutter at supersonic and hypersonic speeds．Applied Mechanics Reviews，1999，52(10):321 ～ 3 32

8 Bellmam R E，Casti J．Differential quadrature and long term integration．Journal of Mathematical Analysis and Applications，1971，34(2):235 ～238

9 Bellmam R E，Kashef B G，Casti J．Differential quadrature:a technique for the rapid solution of nonlinear partial differential equations．Journal of Computational Physics，1972，10(1):40～52

10 Bert C W，Jang S K，Striz A J．Two new approximate methods for analyzing free vibtation of structural components．American Institute of Aeronautics and Astronautics Journal，1988，26(5):612 ～618

11 Bert C W，Malik M．Differential quadrature method in computational mechanics:a review．Applied Mechanics Reviews，1996，49(1):1 ～27

12 李鵬，楊翊仁，魯麗．微分求積法分析二維亞音速壁板的失穩(wěn)問題．動力學與控制學報，2012，10(1):11～14(Li P，Yang Y R，Lu L．Instability analysis of two dimensional thin panels in subsonic flow with differential quadrature method．Journal of Dynamics and Control，2012，10(1):11～14(in Chinese))

13 Puri P，Cowin S C．Plane waves in linear elastic materials with voids．Journal of Elasticity，1985，15(2):167～183

14 Kumar R，Rani L．Interaction due to mechanical and thermal sources in thermoelastic half－space with voids．Journal of Vibration and Control，2005，11(4):499～517

* The project supported by Innovation Program of Shanghai Municipal Education Commission(12YZ074)

? Corresponding author E－mail:yuanyuan_zhu@hotmail．com

NONLINEAR AEROELASTIC CHARACTERISTIC OF A THERMOELASTIC BEAM WITH VOIDS*

Zhu Yuanyuan1?Li Ying2Cheng Changjun2
(1．Department of Computer Science and Technology，Shanghai Normal University，Shanghai200234，China)(2．Department of Mechanics，Shanghai University，Shanghai200072，China)

In this study，based on the complete nonlinear mathematical model obtained from the generalized Hamilton principle of isotropic thermoelastic beams with voids and the aerodynamic pressure loading presented by the first order modified piston theory，the governing differential equations and the deterministic conditions of solutions for a thermoelastic beam with fixed ends and located in a high－speed or an ultra－h(huán)igh speed flow are presented，which are expressed by the axial and lateral displacements and the two“moments”defined by changes of the volume fraction of voids and the temperature field．In order to consider the nonlinear aeroelastic characteristic of the thermoelastic beam with voids subjected to the transverse loading and aerodynamic pressure loading，the differential quadrature method is applied to discretize the governing differential equations on the spatial domains，and a system of ordinary differential equations with respect to time is yielded and solved through the fourth－order Runge－Kutta method．From this，the aeroelastic characteristic of the beam is studied and the effect of parameters is considered as well．

thermoelastic beam with voids， finite deformation， the first order modified piston， differential quadrature method，aeroelastic characteristic

10 July 2012，

5 August 2012．

10．6052/1672－6553－2013－047

2012－07－10 收到第 1 稿，2012－08－05 收到修改稿．

*上海市教委創(chuàng)新基金資助項目(12YZ074)

E－mail:yuanyuan_zhu@hotmail．com