孫蘭鳳，張成波

● （1.四川電力職業(yè)技術(shù)學(xué)院 動(dòng)力工程系，四川成都 610071；2.東方電氣股份有限公司，四川成都 610036）

旋轉(zhuǎn)翼型沖擊水介質(zhì)動(dòng)態(tài)模擬研究

孫蘭鳳1，張成波2

● （1.四川電力職業(yè)技術(shù)學(xué)院 動(dòng)力工程系，四川成都 610071；2.東方電氣股份有限公司，四川成都 610036）

簡要介紹了旋轉(zhuǎn)翼型聲輻射的研究現(xiàn)狀，建立了NACA系列翼型以及環(huán)形水域模型。對(duì)翼型采取剛性體，對(duì)水介質(zhì)采用 Gruneisen狀態(tài)方程來定義材料的壓力，對(duì)模型采用 8節(jié)點(diǎn)實(shí)體單元進(jìn)行網(wǎng)格劃分。采用ALE算法模擬了水域單元的聲壓值，尤其是翼型在同一起始和終止時(shí)間的速度—時(shí)間曲線(直線、正弦線、拋物線以及正切線)被研究了。對(duì)同一單元在不同旋轉(zhuǎn)方式下的聲壓值以及同一旋轉(zhuǎn)方式下不同單元的聲壓值進(jìn)行了比較。

翼型；沖擊；動(dòng)態(tài)模擬

0 引言

翼型是工程技術(shù)領(lǐng)域常用的典型形體，在水利水電工程、航空、動(dòng)力工程等許多技術(shù)部門都有廣泛的應(yīng)用[1]。

近年來，旋轉(zhuǎn)翼型的聲輻射研究受到了極大關(guān)注。許多學(xué)者在實(shí)驗(yàn)和理論方面首先對(duì)振動(dòng)翼型的流場進(jìn)行了測量。文獻(xiàn)[2]引用了Anderson論文中對(duì)振動(dòng)翼非定常流問題的研究，對(duì)翼型流場進(jìn)行了顯示實(shí)驗(yàn)。動(dòng)態(tài)失速問題研究也就是前緣渦分離問題研究以及在推進(jìn)方面對(duì)翼型最優(yōu)運(yùn)動(dòng)方式研究，同時(shí)完成了振動(dòng)翼型全面的推進(jìn)能力和效率的測量。文獻(xiàn)[3]運(yùn)用Rankine奇點(diǎn)一階面元法，以在無限水域中斜航運(yùn)動(dòng)的舵為研究對(duì)象，考慮線性自由表面邊界條件計(jì)算了穿透自由表面定常斜航運(yùn)動(dòng)物體的水動(dòng)力。文獻(xiàn)[4]主要討論了自由液面下三維水翼定常升力繞流的水動(dòng)力系數(shù)計(jì)算問題。文獻(xiàn)[5]運(yùn)用有限體積法結(jié)合雙時(shí)間推進(jìn)技術(shù)求解三維非定常歐拉方程，模擬了有限翼展機(jī)翼在同時(shí)具有上下拍動(dòng)和俯仰運(yùn)動(dòng)狀態(tài)下的非定常流場，計(jì)算了不同狀態(tài)下?lián)湟淼纳巴屏?，分析了各個(gè)影響因素對(duì)撲翼氣動(dòng)特性的影響。文獻(xiàn)[6]運(yùn)用有限體積法結(jié)合雙時(shí)間推進(jìn)技術(shù)數(shù)值求解三維非定常Euler方程，模擬了機(jī)翼在運(yùn)動(dòng)狀態(tài)下的非定常流場，采用雙時(shí)間推進(jìn)法求解非定常Euler方程，驗(yàn)證此方法可以極大地提高非定常流場計(jì)算效率。同時(shí)，劉曉宙[7-10]在余志興[11-12]研究粘性流場結(jié)果的基礎(chǔ)上，對(duì)流體通過渦激振動(dòng)機(jī)翼的聲輻射進(jìn)行了研究，首先考慮了機(jī)翼的強(qiáng)迫轉(zhuǎn)動(dòng)，只有一個(gè)自由度，進(jìn)而考慮了更為復(fù)雜機(jī)翼振動(dòng)，有兩個(gè)自由度，分別是上下振動(dòng)和繞彈性中心的轉(zhuǎn)動(dòng)，得出流體繞渦激振動(dòng)的機(jī)翼引起的聲輻射大于流體繞固定機(jī)翼引起的聲輻射，特別是當(dāng)渦脫頻率和機(jī)翼的固有的振動(dòng)頻率一致時(shí)，聲輻射達(dá)到最大，并且根據(jù)翼型的周期性渦發(fā)放，需要較大的攻角，而攻角的具體大小與翼型有關(guān)系，提出了可以利用改變翼型形狀控制渦激噪聲的可能性。

對(duì)流體與翼型相互作用產(chǎn)生聲輻射的問題，大部分按照流噪聲的思路在研究。對(duì)于運(yùn)動(dòng)翼型，文獻(xiàn)[7-8，10-12]也只是考慮了翼型轉(zhuǎn)動(dòng)角度按照正弦曲線變化的情況，在潛艇實(shí)際航行過程中，輻射噪聲的極值點(diǎn)常出現(xiàn)在舵葉開始加速轉(zhuǎn)動(dòng)的過程中。為了研究旋轉(zhuǎn)翼型從靜止?fàn)顟B(tài)加速旋轉(zhuǎn)至一定角度的過程中聲輻射的特性，本文采用了翼型沖擊水介質(zhì)模型，分別選取適當(dāng)?shù)牟牧夏M了剛性的翼型在旋轉(zhuǎn)過程中對(duì)周圍靜止水域的沖擊作用，研究了旋轉(zhuǎn)方式對(duì)水域壓力的變化的影響，從而得到了旋轉(zhuǎn)翼型沖擊水介質(zhì)的特性以及聲輻射的規(guī)律。

1 翼型與水域模型

為了準(zhǔn)確的表達(dá)翼型，采取工程中常用的NACA系列翼型。針對(duì)翼型的旋轉(zhuǎn)運(yùn)動(dòng)方式，在水域模型上采取了環(huán)形水域，一方面便于計(jì)算時(shí)劃分網(wǎng)格，另一方面在察看結(jié)果時(shí)便于選取觀測點(diǎn)進(jìn)行對(duì)比分析。圖1、2分別給出了建立的翼型模型以及水域模型。

圖1 舵葉殼體翼型模型

圖2 水域模型

1.1 模型單元

翼型和水域模型由8節(jié)點(diǎn)單元?jiǎng)澐帧卧诿總€(gè)節(jié)點(diǎn)在x，y，z方向都具有位移、速度以及加速度。單元默認(rèn)一個(gè)節(jié)點(diǎn)進(jìn)行粘性沙漏控制來加快單元變形，這樣可以節(jié)省計(jì)算時(shí)間以及在大變形的時(shí)候不會(huì)產(chǎn)生單元失效。翼型和水域單元見圖3。

圖3 翼型和水域單元

2 材料選擇

2.1 舵葉殼體的材料模型

在材料選擇方面，因?yàn)檠芯康膫?cè)重點(diǎn)是翼型沖擊水介質(zhì)，目的是為了得到水域壓力的變化，因此翼型采用剛性體，假設(shè)翼型在旋轉(zhuǎn)過程中未產(chǎn)生自身的顫振，并且在與水介質(zhì)作用過程中，不發(fā)生變形。

用剛體來定義模型中較硬的部分能夠大大減少顯示動(dòng)力分析的計(jì)算時(shí)間，剛體將自由度耦合在質(zhì)心，因此無論有多少個(gè)節(jié)點(diǎn)，但各剛體只有6個(gè)自由度。質(zhì)量、質(zhì)心和慣性矩將根據(jù)剛體的體積與單元密度自動(dòng)計(jì)算。作用在剛體上的力與力矩在每個(gè)時(shí)間步由各節(jié)點(diǎn)值相加而成。剛體的運(yùn)動(dòng)首先在質(zhì)心處計(jì)算，然后轉(zhuǎn)換到各個(gè)節(jié)點(diǎn)上，并且剛體不需要網(wǎng)格連續(xù)。尤其是要計(jì)算接觸剛度，剛體材料參數(shù)值要用實(shí)際值。

在計(jì)算時(shí)，初步定義了翼型密度DENS為7800kg/m3，彈性模量為200GPa，泊松比為0.3，先不考慮翼型直線運(yùn)動(dòng)，對(duì)翼型各個(gè)坐標(biāo)方向的運(yùn)動(dòng)進(jìn)行約束；并且在翼型旋轉(zhuǎn)方面，只考慮翼型繞某一坐標(biāo)軸轉(zhuǎn)動(dòng)。

2.2 水介質(zhì)的材料選擇

在分析中，假設(shè)水介質(zhì)是可壓縮的，采用格林納森(Gruneisen)狀態(tài)方程定義壓力體積的關(guān)系，確定材料是壓縮的。

2.2.1 水介質(zhì)的材料狀態(tài)方程

具有立體撞擊速度—粒子的速度狀態(tài)Gruneisen方程定義壓縮材料的壓力為：

Gruneisen方程是一種絕熱熵增的狀態(tài)方程，通常利用測量沖擊波波速—波后質(zhì)點(diǎn)粒子速度曲線（又叫us-up曲線）來決定方程的參數(shù)，us-up曲線中s為震動(dòng)；p為粒子；該曲線采用S1、S2、S3三個(gè)系數(shù)多項(xiàng)式擬合，是us-up曲線的斜率系數(shù)；C是粒子速度為0時(shí)的波速，即曲線在y軸的截距；γ0是Gruneisen參數(shù)；a是對(duì)γ0的修正系數(shù)；E是初始內(nèi)能，常溫下通常設(shè)為0；V是初始相對(duì)體積，即相對(duì)沒有任何變形的體積，初始無體積應(yīng)變時(shí)則設(shè)為1；根據(jù)相對(duì)體積V，壓縮量μ定義為：μ=(1/V)-1。

2.2.2 水介質(zhì)材料模型設(shè)定

水介質(zhì)材料狀態(tài)方程不用計(jì)算偏應(yīng)力，也可以選擇定義粘度。在計(jì)算時(shí)，初步定義水密度DENS為1000kg/ m3，水中聲速C為1647m/s；us-up曲線的斜率系數(shù)S1=1.921，S2=-0.096，S3=0；Gruneisen參數(shù)γ0=0.35，并且初始內(nèi)能E=0，相對(duì)體積V=1。

3 計(jì)算方法

3.1 計(jì)算方法的討論比較

在翼型與水介質(zhì)作用過程中，由于翼型運(yùn)動(dòng)有空間相對(duì)較大的位移并且在實(shí)際情況中，水介質(zhì)和翼型間有物質(zhì)輸送。如果翼型周圍環(huán)繞流體介質(zhì)的單元基于 Lagrange型單元時(shí)，就會(huì)產(chǎn)生嚴(yán)重的變形，直接導(dǎo)致每次迭代計(jì)算時(shí)的時(shí)間步長變得很小，因此模型的計(jì)算時(shí)間變得很長。另外，由于計(jì)算網(wǎng)格的畸變，可能直接影響到數(shù)值計(jì)算的精度。Euler型的有限元模擬雖然提高了在固定網(wǎng)格上采用Navier-Stokes方程時(shí)間解的精度，同時(shí)也增加了計(jì)算的時(shí)間。所以，選擇ALE算法來進(jìn)行分析。

ALE方法最初出現(xiàn)于數(shù)值模擬流體動(dòng)力學(xué)問題的有限差分方法中。這種方法兼具Lagrange方法和Euler方法二者的特長，即首先在結(jié)構(gòu)邊界運(yùn)動(dòng)的處理上引進(jìn)Lagrange方法的特點(diǎn)，因此能夠有效地跟蹤物質(zhì)結(jié)構(gòu)邊界的運(yùn)動(dòng)；其次在內(nèi)部網(wǎng)格的劃分上，吸收了Euler的長處，也就是使內(nèi)部網(wǎng)格單元獨(dú)立于物質(zhì)實(shí)體而存在，但它又不完全和Euler網(wǎng)格相同，網(wǎng)格可以根據(jù)定義的參數(shù)在求解過程中適當(dāng)調(diào)整位置，使得網(wǎng)格不致出現(xiàn)嚴(yán)重的畸變。這種方法在分析大變形問題時(shí)是非常有利的。使用這種方法時(shí)網(wǎng)格與網(wǎng)格之間物質(zhì)也是可以流動(dòng)的。

3.2 分析時(shí)間及時(shí)間積分

考慮阻尼影響的運(yùn)動(dòng)方程為：

時(shí)間積分采用顯式中心差分格式，其算式為：

由于采用集中質(zhì)量矩陣M，運(yùn)動(dòng)方程（3）的求解是非偶合的，不需要組集成總體積矩陣，因此大大節(jié)省了計(jì)算時(shí)間。但是顯式中心差分法是有條件穩(wěn)定的。模擬時(shí)采用變時(shí)不常增量解法。每一時(shí)刻的時(shí)步長由當(dāng)前構(gòu)形的穩(wěn)定性條件控制，其算法為：

線計(jì)算每一個(gè)單元的極限時(shí)步長Δtei，i=1，2… (顯式中心差分法穩(wěn)定性條件允許的最大時(shí)步長)則下一步長取其極小值，即：

式中，Δtei為第i個(gè)單元的極限時(shí)步長；m是單元數(shù)目。

對(duì)于舵葉殼體模型和水域模型所選擇的三維實(shí)體單元：

本文設(shè)定的翼型轉(zhuǎn)動(dòng)的最低速度為3o/s。為了滿足最低的要求，假設(shè)在1s的時(shí)間內(nèi)，舵葉旋轉(zhuǎn)的極值速度為0.0523598rad/s。

4 對(duì)不同運(yùn)動(dòng)曲線的動(dòng)態(tài)模擬

4.1 動(dòng)態(tài)模擬速度曲線

本文對(duì)不同的已知曲線進(jìn)行了模擬。在模擬曲線的選擇中，參考了典型的電梯曲線的組成，選擇了直線，正弦曲線，拋物線以及正切線進(jìn)行分析。

首先，對(duì)于直線速度時(shí)間曲線，選取V=0.0523598t(0≤t≤1s)。其次選取正弦速度時(shí)間曲線為：

V=0.0523598sin (πt/2) (0≤t≤1s)。

最后分別選取拋物線速度時(shí)間曲線及正切速度時(shí)間曲線為：

V= 0.0523598t2(0≤t≤1s)及V=0.0523598tan(πt/4)。

根據(jù)上述分析，選取相同的水域單元進(jìn)行對(duì)比，見圖4。

圖4 部分水域中的舵葉模型

4.2 提取計(jì)算結(jié)果

采取相同的模型和計(jì)算方法，僅改變翼型在相同的時(shí)間內(nèi)達(dá)到相同的極值速度時(shí)采用的速度曲線，分別得到不同的曲線下觀測點(diǎn)的壓力值，如圖5所示。

對(duì)某一水域單元進(jìn)行進(jìn)一步的對(duì)比分析，不同時(shí)刻，不同的舵葉殼體旋轉(zhuǎn)速度曲線在此單元產(chǎn)生的壓力值如圖6所示。

圖5 不同速度曲線的單元壓力

圖6 不同單元的壓力

5 結(jié)果分析及討論

在同一速度曲線的模擬過程下，尾緣處觀測點(diǎn)6502，6522處的聲壓明顯大于前緣處觀測點(diǎn)7579，7580處聲壓。對(duì)于直線速度時(shí)間曲線，尾緣觀測點(diǎn) 6502和 6522在t＜0.664287s時(shí)，6522處聲壓總是大于 6502處，而在t≥0.664287s時(shí)，6522處聲壓轉(zhuǎn)變?yōu)樾∮?502處聲壓，而對(duì)于正弦速度時(shí)間曲線，此分界時(shí)間點(diǎn)在t=0.626316s處，對(duì)于拋物線速度時(shí)間曲線，分界時(shí)間點(diǎn)在t=0.093956s處，正切線速度時(shí)間曲線，分界時(shí)間點(diǎn)在t=0.730527s。這表明在對(duì)翼型旋轉(zhuǎn)攪動(dòng)水介質(zhì)的聲壓研究中，可以直接選取尾緣觀測點(diǎn)，根據(jù)尾緣處聲壓值的比較和聲壓值變化的分界時(shí)間點(diǎn)，在分界時(shí)間的前段可選擇尾緣內(nèi)側(cè)水域作為分析對(duì)象（如水域單元 6522），在分界時(shí)間的后段可選擇外側(cè)水域作為分析對(duì)象（如水域單元6502）。

對(duì)不同的速度曲線下的同一觀測點(diǎn)，拋物線速度曲線在翼型初始啟動(dòng)時(shí)刻的聲壓最小，對(duì)于尾緣6502處，一直到t≥0.929066s后，才有正弦速度時(shí)間曲線下的此處聲壓最小。尾緣6522處，在t≥0.333159s后，正切速度時(shí)間曲線下的聲壓才小于拋物線速度時(shí)間曲線下的聲壓，并且正切曲線下的聲壓并未持續(xù)降低，在t≥0.408957s后，拋物線下的此點(diǎn)聲壓最小，最后在t≥0.777732s后，才有正弦速度時(shí)間曲線下的此處聲壓最小。這表明在翼型旋轉(zhuǎn)初始時(shí)刻，采用拋物線速度時(shí)間曲線可降低水域聲壓值，在翼型旋轉(zhuǎn)速度到達(dá)極值的末端時(shí)刻，采用正弦速度時(shí)間曲線可降低水域聲壓值，可供設(shè)計(jì)水下翼型旋轉(zhuǎn)的聲學(xué)設(shè)計(jì)做參考。

[1]夏國澤.船舶流體力學(xué)[M].武漢:華中科技大學(xué)出版社,2003.

[2]王肇.二維振動(dòng)機(jī)翼的流場測量及數(shù)值分析[D].長沙:中國科學(xué)技術(shù)大學(xué),2004.

[3]汪淳.穿透自由表面傾側(cè)舵的水動(dòng)力計(jì)算[D].武漢:武漢理工大學(xué),1999.

[4]劉華.三維水翼水動(dòng)力數(shù)值計(jì)算[D].武漢:華中科技大學(xué),2003.

[5]龔凱.有限翼展撲動(dòng)翼的歐拉方程數(shù)值模擬[D].西安:西北工業(yè)大學(xué),2003.

[6]周成剛.運(yùn)動(dòng)機(jī)翼流場的歐拉方程求解[D].西安:西北工業(yè)大學(xué),2002.

[7]劉曉宙,繆國平,余志興,等.流體通過二維振動(dòng)翼型的聲輻射研究[J].聲學(xué)學(xué)報(bào),2004,l29(2):122-130.

[8]劉曉宙,繆國平,余志興,等.流體通過渦激振動(dòng)機(jī)翼的聲輻射研究[J].聲學(xué)學(xué)報(bào),2005,30(1):56-62.

[9]劉曉宙,繆國平,唐春權(quán),等.粘性流體通過二維振蕩圓柱的聲輻射研究[J].船舶力學(xué),2002,6(4):75-88.

[10]劉曉宙.旋渦運(yùn)動(dòng)與渦激噪聲的研究[D].上海:上海交通大學(xué),2000.

[11]余志興,劉應(yīng)中,繆國平.二維機(jī)翼彈簧系統(tǒng)的渦激振動(dòng)[J].船舶力學(xué),2002,6(5):25-32.

[12]余志興.粘性流場中的水彈性計(jì)算[D].上海:上海交通大學(xué),1999.

Dynamic Simulation Study on Rotating Aerofoil Impacting Water Medium

SUN Lan-feng1,ZHANG Cheng-bo2
(1.Department of Power Engineering,Sichuan Electric Vocational and Technical College,Chengdu 610071,China; 2.Dongfang Electric Corporation,Chengdu 610036,China)

The present state of study on the sound radiation of rotating aerofoil is simply introduced.The models of NACA series aerofoil and annular fluid are built.The rigid material is adopted in aerofoil and the Gruneisen equation of state is used to define the pressure of the water medium,both models is meshed by eight nodes solid elements.The sound pressure of the water unit is simulated by using ALE algorithm.In particular,four different velocity-time curves with the same starting and ending time and velocity of aerofoil movement (straight line,sine curve,parabolic curve,tangent curve) are developed.The different kinds of sound pressure of the same element with different rotating style and the different elements with the same rotating style are compared with each other.

aerofoil; impact; dynamic simulation

TP391.9;TK733.3

A

孫蘭鳳（1978-），女，講師，研究生，主要研究方向?yàn)榱黧w機(jī)械及工程。