旋轉翼型沖擊水介質動態模擬研究

孫蘭鳳，張成波

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

旋轉翼型沖擊水介質動態模擬研究

孫蘭鳳1，張成波2

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

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

翼型；沖擊；動態模擬

0 引言

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

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

對流體與翼型相互作用產生聲輻射的問題，大部分按照流噪聲的思路在研究。對于運動翼型，文獻[7-8，10-12]也只是考慮了翼型轉動角度按照正弦曲線變化的情況，在潛艇實際航行過程中，輻射噪聲的極值點常出現在舵葉開始加速轉動的過程中。為了研究旋轉翼型從靜止狀態加速旋轉至一定角度的過程中聲輻射的特性，本文采用了翼型沖擊水介質模型，分別選取適當的材料模擬了剛性的翼型在旋轉過程中對周圍靜止水域的沖擊作用，研究了旋轉方式對水域壓力的變化的影響，從而得到了旋轉翼型沖擊水介質的特性以及聲輻射的規律。

1 翼型與水域模型

為了準確的表達翼型，采取工程中常用的NACA系列翼型。針對翼型的旋轉運動方式，在水域模型上采取了環形水域，一方面便于計算時劃分網格，另一方面在察看結果時便于選取觀測點進行對比分析。圖1、2分別給出了建立的翼型模型以及水域模型。

圖1 舵葉殼體翼型模型

圖2 水域模型

1.1 模型單元

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

圖3 翼型和水域單元

2 材料選擇

2.1 舵葉殼體的材料模型

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

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

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

2.2 水介質的材料選擇

在分析中，假設水介質是可壓縮的，采用格林納森(Gruneisen)狀態方程定義壓力體積的關系，確定材料是壓縮的。

2.2.1 水介質的材料狀態方程

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

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

2.2.2 水介質材料模型設定

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

3 計算方法

3.1 計算方法的討論比較

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

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

3.2 分析時間及時間積分

考慮阻尼影響的運動方程為：

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

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

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

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

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

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

4 對不同運動曲線的動態模擬

4.1 動態模擬速度曲線

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

首先，對于直線速度時間曲線，選取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)。

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

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

4.2 提取計算結果

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

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

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

圖6 不同單元的壓力

5 結果分析及討論

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

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

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

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

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

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

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

[6]周成剛.運動機翼流場的歐拉方程求解[D].西安:西北工業大學,2002.

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

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

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

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

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

[12]余志興.粘性流場中的水彈性計算[D].上海:上海交通大學,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-），女，講師，研究生，主要研究方向為流體機械及工程。