諶勝 劉春川 李鳳明?

(1．哈爾濱工業(yè)大學(xué)航天學(xué)院，哈爾濱 1 50001)(2．哈爾濱工業(yè)大學(xué)機電學(xué)院，哈爾濱 1 50001)

基于CFD/CSD方法的亞音速平板結(jié)構(gòu)氣動彈性分析*

諶勝1劉春川2李鳳明1?

(1．哈爾濱工業(yè)大學(xué)航天學(xué)院，哈爾濱 1 50001)(2．哈爾濱工業(yè)大學(xué)機電學(xué)院，哈爾濱 1 50001)

采用CFD/CSD雙向流固耦合算法研究平板結(jié)構(gòu)的氣動彈性耦合特性．首先，采用CFD/CSD算法計算平板結(jié)構(gòu)的顫振臨界速度，并與已有文獻中的實驗結(jié)果進行比較驗證．然后，分別對簡支和固支邊界條件的三維平板結(jié)構(gòu)進行氣動彈性特性分析，計算不同約束情況下流場分布的變化和平板結(jié)構(gòu)的位移響應(yīng)．同時還考慮加肋和結(jié)構(gòu)材質(zhì)對平板結(jié)構(gòu)氣動彈性特性的影響．

平板結(jié)構(gòu)， 亞音速氣流， 氣動彈性耦合特性， CFD/CSD算法， 時域響應(yīng)

引言

隨著航空工業(yè)的迅猛發(fā)展，結(jié)構(gòu)的氣動彈性問題受到了研究人員的廣泛關(guān)注．目前，針對結(jié)構(gòu)氣動彈性方面的研究，已經(jīng)獲得了很多的成果，但主要集中在跨音速與超音速的范圍內(nèi)．1940年，美國Tacoma懸索橋在非常低風(fēng)速(18m/s)氣流的作用下發(fā)生倒塌［1］，使得人們認識到了氣動彈性問題有可能在亞音速氣流作用下明顯的表現(xiàn)出來．另外，由于工程的要求，需要在復(fù)雜的外界條件下減輕結(jié)構(gòu)的質(zhì)量，這樣就導(dǎo)致氣動彈性問題越來越突顯．

關(guān)于平板顫振的問題，國內(nèi)外已有很多研究成果，其中以超聲速的研究居多［2－4］．在亞音速平板顫振研究方面，美國學(xué)者Dowell對前人的研究成果進行了總結(jié)［5－7］．給出了線性和非線性壁板顫振的計算模型，以及平板結(jié)構(gòu)亞音速和跨音速氣動力的計算方法．之后國內(nèi)外有很多學(xué)者采用Galerkin法、微分求積法等方法對亞音速平板結(jié)構(gòu)的氣動彈性問題進行了研究［8－11］．

氣動彈性問題除了上述研究中使用的經(jīng)典方法以外，還發(fā)展出了一些數(shù)值計算方法．其中運用最為廣泛的就是應(yīng)用于飛行器顫振分析的偶極子網(wǎng)格法(DLM)，這是一種計算氣動力的簡化數(shù)值方法，MSC．NASTRAN有限元分析軟件的氣動彈性模塊就是使用的偶極子網(wǎng)格法．但這種算法對于固體場和流體場都有很強的限制，例如該方法只適用于線性、勢流的氣動力模型，無法考慮機翼結(jié)構(gòu)厚度的迎角，對于其它復(fù)雜結(jié)構(gòu)的模型也無法處理［12］．近年來，隨著計算機技術(shù)的快速發(fā)展，計算流體動力學(xué)(CFD)和計算固體動力學(xué)(CSD)都得到了廣泛的應(yīng)用．但是，由于流體域和固體域的基本方程有很大差別，在數(shù)學(xué)特性上也完全不一樣．因此，建立整體的流固耦合方程并耦合求解只能適用于簡單問題的科學(xué)研究，無法進行復(fù)雜問題的大規(guī)模計算．此后，基于CFD/CSD的雙向迭代算法逐漸運用于求解耦合場的問題．

基于CFD/CSD的流固耦合算法早期主要用于計算固體和液體的耦合作用．如水下導(dǎo)彈的研究、船和水的耦合作用等，后來逐漸應(yīng)用于結(jié)構(gòu)的氣動彈性問題分析．Bhardwaj 和 Kapania［13，14］采 用 CFD/CSD的耦合算法對機翼結(jié)構(gòu)進行了靜氣動彈性計算，并與實驗進行對比，驗證了CFD/CSD算法對靜氣動彈性問題的計算準(zhǔn)確性．Biedron和Lee［15］采用非結(jié)構(gòu)化網(wǎng)格建立了UH－60A型機翼的模型，分析了其在高速運行時的流固耦合問題．另外，國內(nèi)外很多學(xué)者采用CFD/CSD算法對AGARD445．6型機翼的氣動彈性問題進行了計算，并與實驗結(jié)果以及理論定性分析結(jié)果進行對比驗證［12，16－19］．除了機翼結(jié)構(gòu)之外，還有學(xué)者使用CFD/CSD算法對渦輪葉片的顫振特性進行了研究［20，21］．總體來說，目前基于CFD/CSD方法的氣動彈性問題研究主要集中在飛機機翼結(jié)構(gòu)上，對于平板等其它結(jié)構(gòu)的研究較少．

本文采用CFD/CSD雙向迭代算法分析亞音速氣流作用下平板結(jié)構(gòu)氣動彈性耦合特性，計算了平板結(jié)構(gòu)的顫振臨界速度，并與已有的實驗結(jié)果進行對比，驗證了采用CFD/CSD算法分析結(jié)構(gòu)氣動彈性耦合特性的有效性．計算了加肋三維平板結(jié)構(gòu)氣動彈性瞬態(tài)響應(yīng)，并討論了不同邊界條件對平板結(jié)構(gòu)氣動彈性響應(yīng)的影響．本文提出的研究方法為研究復(fù)雜平板結(jié)構(gòu)及其組合結(jié)構(gòu)的氣動彈性耦合振動特性提供了一種有效的思路．

1 CFD/CSD算法的基本理論

1．1 基本思想

流體域和固體域分別采用有限單元法和有限體積法單獨求解，然后通過數(shù)據(jù)交換技術(shù)在流固交界面處耦合迭代，最后求得整個域的變量．流固耦合計算的過程與時間相關(guān)，因此是一個典型的瞬態(tài)動力學(xué)問題．

固體場的計算主要采用有限單元法進行結(jié)構(gòu)動力學(xué)的求解．采用有限單元法對固體結(jié)構(gòu)進行離散．在動力學(xué)分析中引入了時間坐標(biāo)，因此這是一個四維(x，y，z，t)問題．通過有限元的離散，連續(xù)性結(jié)構(gòu)的無限維問題就變成與節(jié)點位移有關(guān)的有限維數(shù)的問題．寫成矩陣表達式:

M¨a(t)+Ca(t)+Ka(t)=Q(t)(1)式中a為節(jié)點位移列向量，M為質(zhì)量矩陣，C為阻尼矩陣，K為剛度矩陣，Q為氣動壓力．式中的變量均與時間有關(guān)．式(1)是一個二階常微分方程組，采用Newmark數(shù)值積分方法可以對其進行求解．

流體場的計算則采用有限體積法．有限體積法能保證離散單元的守恒性，且離散方程物理意義清晰，因此在流體計算時有限體積法應(yīng)用最為廣泛．與有限單元法相似，有限體積法也是離散求解的思想．有限體積法將流體計算區(qū)域劃分成若干個控制體，每個控制體都由一個節(jié)點代表，通過對控制體進行積分從而達到離散化的目的．本文主要采用二階迎風(fēng)格式對流體場進行離散處理．

1．2 計算流程

圖1是CFD/CSD雙向流固耦合算法的流程圖．首先，由CFD方法計算出流體域的壓力分布，然后在流固交界面處將流體計算出的壓力施加到固體域上，對固體進行CSD瞬態(tài)動力學(xué)分析．固體產(chǎn)生變形以后，使得流體邊界發(fā)生變化，通過動網(wǎng)格技術(shù)產(chǎn)生新的流體域，計算新的流場分布，然后再作用到固體上，依次循環(huán)．當(dāng)前后兩次的流場參數(shù)分布以及固體場參數(shù)分布的差值小于某一設(shè)定數(shù)值時，則認為迭代收斂．之后進入下一個時間步計算，計算方法與上一時間步相同．最終計算完成以后可以得到固體場和流體場在每個時刻的結(jié)果，如壓力分布、速度、位移等．

固體場產(chǎn)生變形傳遞給流體場時涉及到動網(wǎng)格的技術(shù)．本文主要使用基于彈性變形的網(wǎng)格調(diào)整方法．基于彈性變形的網(wǎng)格調(diào)整方法將整個網(wǎng)格視為一個彈簧系統(tǒng)，通過邊界節(jié)點上的位移來調(diào)整整個流體域的節(jié)點位置，這樣的方法不改變網(wǎng)格之間的連貫性．將節(jié)點想象成一個質(zhì)點，節(jié)點間的連線想象成一個彈簧，邊界節(jié)點上的位移導(dǎo)致與節(jié)點相連的“彈簧”產(chǎn)生相應(yīng)的變形，從而通過整個網(wǎng)格傳播出去．每個節(jié)點周圍的“彈簧”對節(jié)點產(chǎn)生的力是平衡的．在邊界位移已知的情況下，可通過Jacobi矩陣對流域內(nèi)部所有節(jié)點位移進行求解，從而得到網(wǎng)格變形以后的節(jié)點位移．

圖1 CFD/CSD算法流程圖Fig．1 The flow diagram of CFD/CSD algorithm

2 CFD/CSD算法驗證

本文選取美國學(xué)者 Dugundji和Dowell［22］文獻中的結(jié)果進行對比，以驗證算法的有效性．文獻［22］中計算了兩邊簡支的彈性支撐平板結(jié)構(gòu)在亞音速氣流作用下的顫振特性．平板材料是鋁，其長2．64m，寬 0．61m，厚 0．46mm．彈性支撐的彈性系數(shù)K=18879．91N/m2．文獻中實驗測得的顫振臨界速度為42．47m/s．對于線性顫振問題，當(dāng)來流速度大于顫振臨界速度時，位移響應(yīng)發(fā)散;小于顫振臨界速度時位移響應(yīng)收斂;來流速度恰好等于顫振臨界速度時，位移響應(yīng)處于等幅運動的臨界狀態(tài)．本文基于以上原則判斷結(jié)構(gòu)的顫振臨界速度．

本文采用CFD/CSD算法對上述模型進行計算．圖2是通過CFD/CSD算法計算出的平板中點時域響應(yīng)曲線．如圖2(a)所示，初始時刻平板上作用有一個外激勵，平板中點位移初始有一些波動，位移大致趨于衰減．t=1s左右運動趨于穩(wěn)定，隨著時間的增長，位移持續(xù)發(fā)散．因此可以說明，40m/s是大于顫振臨界速度的．但此時位移發(fā)散的速率不大，此時比較接近臨界速度．圖2(b)中當(dāng)來流速度為37m/s時，位移初始衰減然后趨于小幅的等幅運動．因此，從圖2可以判斷出顫振點在37m/s附近，與文獻［22］中得到的顫振臨界速度42．47m/s是很接近的．這表明，采用CFD/CSD算法對平板結(jié)構(gòu)進行氣動彈性分析是可行的．

圖2 平板中點處時域響應(yīng)Fig．2 The time domain responses at the midpoint of plate

3 算例分析

在實際工程中很多平板結(jié)構(gòu)通常可以簡化成四邊簡支或四邊固支的情況，本節(jié)針對加肋板結(jié)構(gòu)的四邊分別為簡支和固支情況進行分析研究．

圖3 計算模型圖Fig．3 The calculating model

圖4 四邊簡支加肋平板壓力云圖Fig．4 The pressure cloud pictures of simply supported ribbed plate

圖3(a)是平板結(jié)構(gòu)和流體場的計算模型．平板長L=1m，寬b=1m，厚度h=3mm，平板結(jié)構(gòu)沿X方向加了三根肋條，其尺寸為高2cm，寬1cm．來流速度為140m/s，采用層流理論進行計算．左邊為流體入口，右邊為壓力出口，平板處為流固耦合交界面，其它部分為默認壁面設(shè)置．平板初始時刻作用有均布壓力p=1000Pa，作用時間0．01s．

圖3(b)是流場網(wǎng)格劃分情況．流場是規(guī)則的矩形，可以方便地進行結(jié)構(gòu)化網(wǎng)格劃分，將流場分割成九個矩形以便于對網(wǎng)格尺寸進行控制．距離平板近的流域網(wǎng)格細分以較好地捕捉流場變化，距離平板較遠的地方由于流場變化較小，采用較粗的網(wǎng)格即可滿足計算要求．這樣的網(wǎng)格劃分既能滿足計算精度要求，又能合理控制單元數(shù)，提高計算效率．該模型流場網(wǎng)格共有約32萬個網(wǎng)格．

首先考慮四邊簡支的情況，計算中四邊簡支的加肋平板結(jié)構(gòu)材料分別取為鋼和鋁合金．

圖4是四邊簡支加肋鋼板的計算結(jié)果，圖中顯示了不同時刻流體域中間剖面和平板上的壓力分布．其中為了便于觀察，平板的變形經(jīng)過了放大處理．圖4(a)中由于初始時刻的壓力作用，變形較大，流場的變化也比較大．平板中間處于正壓區(qū)，邊界處有少數(shù)區(qū)域?qū)儆谪搲簠^(qū)．Y=0．5m的中間剖面上，平板附近的壓力變化比較明顯，平板以外的流域壓力變化較小．當(dāng)平板位移運動到Z軸正方向時，壓力分布則正好相反．對于圖4(b)，當(dāng)時間t=0．25s時，平板的位移已經(jīng)很小，流場的壓強也相對較?。?/p>

圖4(c)是t=0．25s時平板上的壓力分布圖．由于整個模型以及邊界條件、初始條件都是關(guān)于X方向?qū)ΨQ的，因此流場分布基本上也是關(guān)于X方向?qū)ΨQ．這一點可以從圖4(c)和圖4(d)中看出．圖4(d)中遠離平板的剖面上壓力變化較小，壓力值也比較?。f明平板的運動變形對于遠處的流場影響較小，在建模時將流場的區(qū)域取足夠大時就能夠得到更接近于真實的結(jié)果．從圖4中還可以明顯看到，流場沿Y方向是變化的，并且變化比較明顯．因此，在對平板結(jié)構(gòu)進行氣動彈性分析時，如果平板結(jié)構(gòu)是三維模型，采用二維流場模型進行流固耦合計算是不夠的，將會產(chǎn)生較大誤差．

本文對四邊簡支的加肋鋁合金板進行了同樣的計算，其流場分布規(guī)律與鋼板類似．圖5是四邊簡支加肋平板結(jié)構(gòu)中點處的響應(yīng)曲線，其中圖5(a)材料為鋼，圖5(b)材料為鋁合金．由于邊界條件為簡支，其整體結(jié)構(gòu)剛度比較大，兩種材料的平板結(jié)構(gòu)位移響應(yīng)幅值都比較小，加肋鋼板的最大位移幅值在1mm左右，鋁合金板的最大位移幅值在2．5mm左右．而且兩者運動趨勢均為收斂，這表明在來流速度為140m/s時，兩種材質(zhì)的平板結(jié)構(gòu)均未發(fā)生顫振．

圖5 四邊簡支加肋平板中點的時域響應(yīng)Fig．5 The time domain responses at the midpoint of simply supported ribbed plate

下面對四邊固支邊界條件下的鋼板和鋁合金板進行計算．如圖6為四邊固支加肋鋼板流固耦合計算的壓力分布云圖．分布狀況與四邊簡支板基本類似，但由于固支板的剛度較大，因此平板結(jié)構(gòu)位移響應(yīng)非常小，由此導(dǎo)致流場壓強也非常小．當(dāng)時間t=0．1s時，平板運動產(chǎn)生的壓強只有2Pa左右，而到t=0．25s時位移已經(jīng)接近于零，此時流場的變化非常小，壓強數(shù)值也趨近于零．

圖6 四邊固支加肋平板壓力分布云圖Fig．6 The pressure cloud pictures of clamped ribbed plate

圖7是氣動力作用下四邊固支板中點的響應(yīng)曲線．由于固支加肋板的剛度比較大，因此位移幅值很小，最大值分別只有0．1mm和0．2mm．在初始沖擊載荷的作用下，響應(yīng)曲線中第一個周期位移響應(yīng)相對較大，隨著時間的增長，響應(yīng)逐漸衰減，并且很快就衰減到位移接近于零．此后，平板和來流的相互作用越來越小，逐漸衰減至最后靜止．此處研究結(jié)果表明，當(dāng)增強結(jié)構(gòu)的約束條件，也即增大結(jié)構(gòu)剛度時，對結(jié)構(gòu)的氣動彈性耦合效應(yīng)將有明顯的提高．

圖7 四邊固支加肋平板中點的時域響應(yīng)Fig．7 The time domain responses at the midpoint of clamped ribbed plates

4 結(jié)論

本文主要對CFD/CSD數(shù)值計算方法進行了對比驗證，并采用CFD/CSD算法對平板結(jié)構(gòu)的氣動彈性問題進行了研究．首先將CFD/CSD數(shù)值計算結(jié)果與已有文獻的實驗結(jié)果進行了比較，驗證了CFD/CSD算法對于計算氣動彈性問題的適用性．然后分別對三維簡支和固支約束條件的平板結(jié)構(gòu)進行計算，分析了不同情況下流場的變化和平板結(jié)構(gòu)的位移響應(yīng)．通過分析，得到以下結(jié)論:

(1)CFD/CSD雙向流固耦合算法適用于計算結(jié)構(gòu)的氣動彈性問題，計算結(jié)果較好．其優(yōu)點是受氣動力模型和結(jié)構(gòu)模型的限制較少，適用于工程問題中的大規(guī)模計算，并且其計算結(jié)果以圖形的方式直觀顯示，可以為理論分析提供參考．其不足之處是計算量大，尋找顫振點花費時間較多，算法還有待進一步改進提高．

(2)在對三維平板結(jié)構(gòu)進行氣動彈性分析時，如果流場采用二維模型是遠遠不夠的，將會產(chǎn)生較大誤差．

(3)當(dāng)增大結(jié)構(gòu)剛度時，對結(jié)構(gòu)的氣動彈性耦合效應(yīng)將有顯著影響．其中，鋼材比鋁合金材料的結(jié)構(gòu)位移小，四邊固支平板結(jié)構(gòu)的位移響應(yīng)較四邊簡支平板結(jié)構(gòu)衰減更快．

1 伏欣H W著，沈克揚 譯．氣動彈性力學(xué)原理．上海:上海科學(xué)技術(shù)文獻出版社，1982(Fxin H W(written)，Shen K Y(translation)．Principle of aero－elastic mechanics．Shanghai:Shanghai Science and Technology Literature Press，1982(in Chinese))

2 Lee I，Lee D M，OH I K．Supersonic flutter analysis of stiffened laminated plates subject to thermal load．Journal of Sound and Vibration，1999，224(1):49～67

3 Li F M．Active aeroelastic flutter suppression of a supersonic plate with piezoelectric material．International Journal of Engineering Science，2012，51:190～203

4 Song Z G，Li F M．Active aeroelastic flutter analysis and vibration control of supersonic composite laminated plate．Composite Structures，2012，94:702～713

5 Dowell E H．Panel flutter:a review of the aeroelastic stability of plates and shells．AIAA Journal，1970，8(3):385～399

6 Dowell E H．A modern course in aeroelasticity．Kluwer Academic Publishers，2004

7 Dowell E H．Aeroelasticity of plates and shells．Noordhoff International Publishing，1975

8 Tang D，Dowell E．H．Limit cycle oscillations of two－dimensional panels in low subsonic flow．International Journal of Non-Linear Mechanics，2002，37:1199～1209

9 Tang D M，Yamamoto H，Dowell E H．Flutter and limit cycle oscillations of two－dimensional panels in three－dimensional axial flow．Journal of Fluids and Structures，2003，17:225～242

10 Eloy C，Souilliez C，Schouveiler L．Flutter of a rectangular plate．Journal of Fluids and Structures，2007:23:904～919

11 李鵬，楊翊仁，魯麗．微分求積法分析二維亞音速壁板的失穩(wěn)問題．動力學(xué)與控制學(xué)報，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))

12 盧學(xué)成，葉正寅，張陳安．基于ANSYS/CFX耦合的機翼顫振分析．計算機仿真，2010，27(9):88～91(Lu X C，Ye Z Y，Zhang C A．A coupled ANSYS/CFX method for the wing flutter calculation．Computer Simulation，2010，27(9):88～91(in Chinese))

13 Bhardwaj M K，Kapania R K．A CFD/CSD interaction methodology for aircraft wings．NASA Technical Reports，1998

14 Kapania R K，Bhardwaj M K．Aeroelastic analysis of modern complex wings．NASA Technical Reports，1996

15 Biedron B T，Lee R E M．Computation of UH－60A airloads using CFD/CSD coupling on unstructured meshes．67th American Helicopter Society International Annual Forum 2011，2011

16 Liu F，Cai J，Zhu Y，et al．Calculation of wing flutter by a coupled CFD－CSD method．38th Aerospace Sciences Meeting＆ Exhibit，2000

17 Liu F，Cai J，Zhu Y，et al．Calculation of wing flutter by a coupled fluid－structure method．Journal of Aircraft，2001，38(2):334～342

18 Xi R，Jia H．A strong coupled CFD－CSD method on computational aeroelastity．Mechanic Automation and Control Engineering，2011:5512～5515

19 曾憲昂，徐敏，安效民等．基于CFD/CSD耦合算法的機翼顫振分析．西北工業(yè)大學(xué)學(xué)報，2008，26(1):79～82(Zeng X Y，Xu M，An X M，et al．Wing flutter analysis using CFD/CSD algorithm．Journal of Northwestern Polytechnical University，2008，26(1):79 ～ 8 2(in Chinese))

20 Zhang J P，Pan L L．Three－dimensional modeling and aeroelastic coupling analysis for the wind turbine blade．World Non－Grid－Connected Wind Power and Energy Conference，2009:24 ～26

21 張興國，劉鋒，張陳安．基于CFD技術(shù)的葉片顫振分析．航空計算技術(shù)，2009，39(4):75～78(Zhang X G，Liu F，Zhang C A．Flutter computation of turbomachinery blades based on CFD．Aeronautical Computing Technique，2009，39(4):75～78(in Chinese))

22 Dugundji J，Dowell E H，Perkin B．Subsonic flutter of panels on continuous elastic foundations．AIAA Journal，1963，1(5):1146～1154

*The project supported by the National Basic Research Program of China(2011CB711100)and the National Natural Science Foundation of China(11172084，10672017)

? Corresponding author E－mail:fmli@hit．edu．cn

AEROELASTIC ANALYSIS OF PLATE STRUCTURES IN SUBSONIC AIR FLOW BASED ON CFD/CSD ALGORITHM*

Shen Sheng1Liu Chunchuan2Li Fengming1?
(1．School of Astronautics，Harbin Institute of Technology，Harbin150001，China)(2．School of Mechatronic Engineering，Harbin Institute of Technology，Harbin150001，China)

The two－way fluid－structure coupling algorithm，i．e．the CFD/CSD algorithm，is used to study the aeroelastic coupliing characteristics of plate structures in subsonic air flow．Firstly，the flutter critical velocity of the plate structure is computed employing the CFD/CSD algorithm．The results obtained by the CFD/CSD are verified by comparing with the experiment results in the open literature．Secondly，the aeroelastic characteristics for three－dimensional plate structures with simply supported and clamped edges are analyzed．The variation of flow field distribution and the displacement responses of the plate under different restricted conditions are claculated．Meanwhile，the influences of the stiffener and structural material parameters on the aeroelastic characteristices of the plate structure are also considered．

plate structure， subsonic air－flow， aeroelastic coupling characteristics， CFD/CSD algorithm，time domain response

5 July 2012，

15 November 2012．

10．6052/1672－6553－2013－037

2012－07－05 收到第 1 稿，2012－11－15 收到修改稿．

*國家重點基礎(chǔ)研究發(fā)展計劃(973計劃)(2011CB711100)和國家自然科學(xué)基金資助項目(11172084，10672017)

E－mail:fmli@hit．edu．cn