孫蘭鳳，張成波

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

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

孫蘭鳳1，張成波2

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

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

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

0 引言

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

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

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

1 翼型與水域模型

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

圖1 舵葉殼體翼型模型

圖2 水域模型

1.1 模型單元

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

圖3 翼型和水域單元

2 材料選擇

2.1 舵葉殼體的材料模型

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

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

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

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為震動；p為粒子；該曲線采用S1、S2、S3三個系數(shù)多項式擬合，是us-up曲線的斜率系數(shù)；C是粒子速度為0時的波速，即曲線在y軸的截距；γ0是Gruneisen參數(shù)；a是對γ0的修正系數(shù)；E是初始內(nèi)能，常溫下通常設(shè)為0；V是初始相對體積，即相對沒有任何變形的體積，初始無體積應(yīng)變時則設(shè)為1；根據(jù)相對體積V，壓縮量μ定義為：μ=(1/V)-1。

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

水介質(zhì)材料狀態(tài)方程不用計算偏應(yīng)力，也可以選擇定義粘度。在計算時，初步定義水密度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，相對體積V=1。

3 計算方法

3.1 計算方法的討論比較

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

ALE方法最初出現(xiàn)于數(shù)值模擬流體動力學(xué)問題的有限差分方法中。這種方法兼具Lagrange方法和Euler方法二者的特長，即首先在結(jié)構(gòu)邊界運(yùn)動的處理上引進(jìn)Lagrange方法的特點(diǎn)，因此能夠有效地跟蹤物質(zhì)結(jié)構(gòu)邊界的運(yùn)動；其次在內(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)重的畸變。這種方法在分析大變形問題時是非常有利的。使用這種方法時網(wǎng)格與網(wǎng)格之間物質(zhì)也是可以流動的。

3.2 分析時間及時間積分

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.2 提取計算結(jié)果

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

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

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

圖6 不同單元的壓力

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[12]余志興.粘性流場中的水彈性計算[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-），女，講師，研究生，主要研究方向為流體機(jī)械及工程。