船舶螺旋槳流固耦合穩(wěn)態(tài)求解算法

張帥，朱錫，侯海量

(海軍工程大學(xué) 船舶與動(dòng)力工程學(xué)院，湖北 武漢 430033)

水面艦艇螺旋槳大多采用錳-鎳-鋁-銅或鎳-鋁-青銅合金制成.盡管合金材料具有屈服強(qiáng)度高和可靠性好等優(yōu)勢(shì)，但加工成螺旋槳幾何形狀的難度較高，且金屬螺旋槳較差的聲學(xué)性能使得其極易因振動(dòng)而產(chǎn)生噪聲.而纖維增強(qiáng)復(fù)合材料具有比強(qiáng)度高、耐腐蝕性好、良好的阻尼特性以及材料可設(shè)計(jì)性強(qiáng)等優(yōu)點(diǎn).因此復(fù)合材料應(yīng)用于艦船螺旋槳和潛艇推進(jìn)系統(tǒng)中［1］將改善螺旋槳的綜合性能.但不同于金屬螺旋槳，復(fù)合材料螺旋槳在水動(dòng)力作用下會(huì)產(chǎn)生大的變形，其水動(dòng)力性能必然發(fā)生變化.因此要得到復(fù)合材料螺旋槳的水動(dòng)力性能，螺旋槳的流-固耦合算法是基礎(chǔ).傳統(tǒng)的螺旋槳理論設(shè)計(jì)與計(jì)算建立在勢(shì)流理論基礎(chǔ)之上，未能全面考慮粘性的影響且不考慮旋度，因而無法準(zhǔn)確預(yù)測(cè)槳葉邊界層、螺旋槳尾流場(chǎng)的結(jié)構(gòu)及槳葉梢渦的形成等真實(shí)情況下的流動(dòng)特性［2］.基于RANS方程的粘性流場(chǎng)計(jì)算螺旋槳的流場(chǎng)特性的方法日趨完善，黃勝等［3-5］等分析了螺旋槳在不同工作狀態(tài)下的水動(dòng)力性能.關(guān)于螺旋槳流-固耦合算法的研究，LIN H J等［6］采用升力面法和九節(jié)點(diǎn)退化殼單元耦合算法，實(shí)現(xiàn)了求解復(fù)合材料螺旋槳的水動(dòng)力性能的算法.Young Y L［7］研究了面元法和軟件ABAQUS耦合的螺旋槳流-固耦合計(jì)算方法，但這些方法均是基于勢(shì)流理論的螺旋槳水動(dòng)力計(jì)算.

本文首先通過不同湍流模型和試驗(yàn)，驗(yàn)證求解螺旋槳的敞水性能，再將粘性流場(chǎng)計(jì)算軟件和有限元耦合起來，推導(dǎo)變形后螺旋槳的幾何參數(shù)，實(shí)現(xiàn)螺旋槳流-固耦合的穩(wěn)態(tài)求解算法.本文計(jì)算采用實(shí)驗(yàn)測(cè)得的復(fù)合材料板拉伸模量和泊松比，將復(fù)合材料考慮為各向同性，沒有考慮復(fù)合材料的鋪層結(jié)構(gòu).

1 數(shù)學(xué)模型

1.1 流體控制方程

考慮螺旋槳在粘性湍流中旋轉(zhuǎn)，其連續(xù)性方程和動(dòng)量方程可表示為:

連續(xù)方程:

動(dòng)量方程:

RANS方程雖然不用求解流場(chǎng)中的瞬時(shí)量，但是方程卻引入了新的未知量-雷諾應(yīng)力，這時(shí)的方程解不封閉.要封閉求解方程，就必須引入新的方程，這些方程通過各種湍流模型來定義.

1.2 湍流模型的選取

為使RANS方程封閉，將雷諾應(yīng)力用低階或時(shí)均量表達(dá)，即湍流模型.封閉RANS方程主要有REYNOLDS應(yīng)力方程模型(RSM)與渦粘模型2種.渦粘模型以湍流各向同性為基礎(chǔ)，認(rèn)為雷諾應(yīng)力和時(shí)均速度呈線性關(guān)系，該類模型求解簡(jiǎn)單，計(jì)算容易收斂.渦粘模型應(yīng)用較廣泛的有κ-ε，RANG κ-ε模型和SST κ-ω模型.而RSM模型以湍流各向異性為基礎(chǔ)，考慮雷諾應(yīng)力的對(duì)流和擴(kuò)散作用，直接尋找雷諾應(yīng)力的輸運(yùn)方程.本文通過比較分析幾種常用湍流控模型，通過求解時(shí)間和求解精度方面的對(duì)比以及考慮到流固耦合本身對(duì)計(jì)算機(jī)性能的要求等方面進(jìn)行綜合考慮.螺旋槳水動(dòng)力性能計(jì)算模型選擇SST κ-ω湍流模型.該模型是利用混合函數(shù)將κ-ε和κ-ω2方程相結(jié)合而構(gòu)建的湍流模型，在近壁區(qū)采用κ-ω方程，其他區(qū)域則采用κ-ε方程以獲得湍流粘性作用，考慮了κ-ω方程近壁區(qū)模擬時(shí)的有效性及遠(yuǎn)場(chǎng)區(qū)無法準(zhǔn)確模擬的不足［8］.最終選擇SST κω模型求解螺旋槳的流固耦合特性.

1.3 計(jì)算方法

計(jì)算螺旋槳流場(chǎng)的控制方程是一系列非線性偏微分?jǐn)?shù)學(xué)物理方程，需借助數(shù)值方法對(duì)其進(jìn)行求解.本文利用Ansys/cfx軟件完成該數(shù)值計(jì)算.收斂判據(jù)設(shè)定為0.000 1.計(jì)算區(qū)域分為內(nèi)、外2個(gè)流域，在內(nèi)流場(chǎng)建立一個(gè)固定于螺旋槳的旋轉(zhuǎn)坐標(biāo)系，采用MRF坐標(biāo)系模型對(duì)螺旋槳周圍的旋轉(zhuǎn)流場(chǎng)進(jìn)行計(jì)算，外圍流場(chǎng)則在絕對(duì)靜止坐標(biāo)系下進(jìn)行求解.將進(jìn)口邊界設(shè)置為速度入口條件，給定均勻來流的速度分量;出口邊界給定靜壓分布，外邊界設(shè)為開放界面;考慮粘性影響，螺旋槳表面定義為不可滑移壁面.

2 流體計(jì)算模型

文中采用DTMB4119螺旋槳，它是一種無側(cè)斜無后傾分布的三葉螺旋槳，被ITTC選為驗(yàn)證數(shù)值方法預(yù)報(bào)精度的標(biāo)準(zhǔn)模型.槳葉直徑為0.304 8 m，盤面比為0.6，槳葉剖面為NACA-66(mod)型，轂徑比為0.2，螺距比(0.7R)為1.084.文中采用Excel編制計(jì)算過程文件，計(jì)算出螺旋槳葉面和葉背各個(gè)半徑處的型值點(diǎn)，將計(jì)算出的型值點(diǎn)轉(zhuǎn)換為文本文件，然后導(dǎo)入SolidWorks軟件中，建立三維實(shí)體螺旋槳模型.螺旋槳幾何模型和坐標(biāo)如圖1所示.坐標(biāo)軸的定義為:x軸與螺旋槳的旋轉(zhuǎn)軸一致，指向下游;y軸與槳葉參考線一致;z軸符合右手法則.計(jì)算采用全尺寸模型，轉(zhuǎn)速為n=10 r/s和20 r/s，通過改變來流速度來實(shí)現(xiàn)不同的進(jìn)速系數(shù).

網(wǎng)格質(zhì)量直接決定計(jì)算的收斂性、效率和精度，因此，應(yīng)根據(jù)流場(chǎng)中各物理量的分布特點(diǎn)對(duì)計(jì)算域進(jìn)行合理的網(wǎng)格劃分［1］.

圖1 4119槳幾何模型Fig.1 Geometry of propeller 4119

將求解域分成旋轉(zhuǎn)區(qū)和靜止區(qū)2個(gè)區(qū)，采用非結(jié)構(gòu)網(wǎng)格劃分方法，2個(gè)區(qū)采用CFX的GGI方式連接.即首先在螺旋槳表面生成三角形網(wǎng)格，再通過值要控制在30～300，通過多次試算確定表面網(wǎng)格大小和槳葉邊界層的過渡方式和層數(shù).螺旋槳表面網(wǎng)格和壁面棱柱過渡層網(wǎng)格如圖2所示.總網(wǎng)格數(shù)為120×104.

3 結(jié)構(gòu)模型

通過流體計(jì)算軟件求出螺旋槳的水動(dòng)力載荷以后通過Ansys-cfx軟件指定流固耦合界面將流體壓力通過表面效應(yīng)單元的方式傳遞給有限元單元.因流體計(jì)算和結(jié)構(gòu)計(jì)算采用的是非同種單元類型，為保證求解精度和數(shù)據(jù)傳遞的準(zhǔn)確性，在劃分有限元網(wǎng)格時(shí)保證導(dǎo)邊和隨邊以及葉梢附近的網(wǎng)格要密一些.經(jīng)多次試算后的有限元網(wǎng)格如圖3所示.考慮到螺旋槳自身質(zhì)量和運(yùn)轉(zhuǎn)工況，在ANSYS中設(shè)置旋轉(zhuǎn)軸和施加旋轉(zhuǎn)速度即可施加離心力和重力作用，槳葉根部邊界為固支端.計(jì)算模型采用的材料參數(shù)如表1所示，其中S玻璃纖維為實(shí)驗(yàn)所測(cè)板拉結(jié)果，計(jì)算出的鋁青銅材料槳葉的變形比S玻纖槳葉變形低一個(gè)數(shù)量級(jí)，因此主要考慮高強(qiáng)S玻纖的槳葉結(jié)構(gòu).計(jì)算用復(fù)合材料槳幾何和金屬槳幾何相同.

圖3 槳葉有限元計(jì)算模型Fig.3 Finite element model of the blade

表1 槳葉材料參數(shù)Table 1 Properties of the materials

4 螺旋槳敞水性能計(jì)算和驗(yàn)證

對(duì)于流體動(dòng)力載荷，由于槳葉工作于復(fù)雜的流場(chǎng)中，葉面和葉背受到分布載荷，這種載荷既不均勻，也不滿足一些簡(jiǎn)單規(guī)律，因此如何盡可能真實(shí)地模擬槳葉的載荷分布是流固耦合分析的關(guān)鍵所在.

4.1 螺旋槳敞水性能計(jì)算和對(duì)比

為了實(shí)現(xiàn)結(jié)構(gòu)載荷計(jì)算的準(zhǔn)確性，首先采用幾種常用的湍流模型(κ-ε，SST κ-ω和RSM)對(duì)比求解螺旋槳的水動(dòng)力敞水性能并和JESSUP S D［9］實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比.

對(duì)比結(jié)果如圖4、5所示，螺旋槳敞水動(dòng)力參數(shù)的計(jì)算值和實(shí)驗(yàn)結(jié)果吻合較好，計(jì)算所得的KT和10KQ與實(shí)驗(yàn)結(jié)果的最大誤差除了κ-ε模型為12.5%以外，SSTκ-ω為9.6%和RSM為8%，而對(duì)于螺旋槳效率的計(jì)算值僅有κ-ε模型超過了11.5%，而其他情況誤差均在5%.需要指出的是，最大誤差值均出現(xiàn)于最大進(jìn)速情況.原因在于相同轉(zhuǎn)速和相同直徑下，進(jìn)速越大，相應(yīng)的推力和扭矩越小，任何一個(gè)干擾就會(huì)導(dǎo)致預(yù)測(cè)值和實(shí)驗(yàn)值的誤差增加很大.整體來看SSTκ-ω模型和RSM模型均能得到精度較高的計(jì)算結(jié)果.圖6為在同臺(tái)計(jì)算機(jī)上不同湍流模型在不同進(jìn)速下的計(jì)算時(shí)間對(duì)比圖.由圖6可知，相同進(jìn)速下RSM模型的求解時(shí)間最長(zhǎng)，在低進(jìn)速時(shí)尤為明顯，κ-ε求解時(shí)間最短，而SSTκ-ω模型求解時(shí)間接近κ-ε模型.

從以上對(duì)比可知，SSTκ-ω模型可以在保證求解精度的情況下，求解時(shí)間減少很多.另考慮到流固耦合計(jì)算本身對(duì)計(jì)算機(jī)求解性能和求解時(shí)間的嚴(yán)格要求，綜合考慮，螺旋槳流體計(jì)算湍流模型采用SSTκω模型.

圖4 計(jì)算KT和10KQ與實(shí)驗(yàn)結(jié)果的對(duì)比Fig.4 Comparisons of calculated and experimental KT and 10 KQ

圖5 計(jì)算效率ETA和實(shí)驗(yàn)結(jié)果對(duì)比Fig.5 Comparisons of calculated and experimental ETA

4.2 螺旋槳表面壓力分布對(duì)比

圖7為螺旋槳在設(shè)計(jì)工況J=0.833時(shí)，r/R= 0.3、0.7和0.9半徑上JESSUP S D實(shí)驗(yàn)換算得到的弦向壓力分布和采用SSTκ-ω計(jì)算值的比較.

圖6 求解時(shí)間對(duì)比(n=10r/s)Fig.6 Comparisons of calculated durations(n=10 r/s)

圖7 葉切面壓力分布系數(shù)Fig.7 Comparisons of pressure coefficient

5 變形后槳葉幾何參數(shù)的推導(dǎo)

槳葉參考線即葉剖面鼻尾線中點(diǎn)的坐標(biāo):［10］

設(shè)導(dǎo)邊的坐標(biāo)為(xl，r，θl)或(xl，yl，zl)，隨邊的坐標(biāo)為(xt，r，θt)或(xt，yt，zt)，葉寬(即弦長(zhǎng))分布為C(r)，則槳葉的輪廓線可表示:

式中:C(r)為葉剖面弦長(zhǎng);φ(r)為葉剖面螺距角;下標(biāo)l、t分別表示導(dǎo)邊(取－)和隨邊(取+).

選擇無量綱弦長(zhǎng)s，導(dǎo)邊表示為0，隨邊表示為1，槳葉弦向中點(diǎn)為1/2，槳葉剖面的拱度和厚度分布分別表示為f(s)和t(s)，槳葉拱弧面的表達(dá)式為:

式中:下標(biāo)c表示為拱弧面;δk=2π(k－1)/K，k= 1，2…，K，為槳葉數(shù).

變形后各葉切面半徑為rref:

在推導(dǎo)式(9)的過程中，由于螺旋槳的變形量很小，所以認(rèn)為C(r)不變.將和代入式(4)就可得到變形后的縱傾(r).最后將所有變形后的節(jié)點(diǎn)坐標(biāo)代入到式(5)槳葉拱弧面的方程中，就可得到變形后的拱弧f＇(s)

6 螺旋槳流固耦合算法的實(shí)現(xiàn)

復(fù)合材料螺旋槳槳葉受到水動(dòng)力載荷和離心力的作用，會(huì)產(chǎn)生較大的變形，那么槳葉結(jié)構(gòu)的離散方程為

式中:K為槳葉剛度矩陣，u槳葉節(jié)點(diǎn)位移矢量矩陣，F(xiàn)為槳葉所受外載荷矩陣.而式(11)中位移矢量矩陣u需采用有限元軟件計(jì)算，而F求解需要流體計(jì)算軟件.兩者求解的平衡實(shí)現(xiàn)需借助流固耦合算法.因此設(shè)計(jì)了一種求解螺旋槳性能的流固耦合算法，求解流程如圖8所示.

圖8 流固耦合的算法求解流程圖Fig.8 Flow chart of computation algorithm of FSI

計(jì)算采用穩(wěn)態(tài)求解，不考慮瞬時(shí)效應(yīng).詳細(xì)求解過程為:

1)通過第2節(jié)中的方法建立螺旋槳的流體計(jì)算模型，計(jì)算出螺旋槳的敞水性能和壓力分布.

2)將1)計(jì)算的壓力載荷結(jié)果通過流固耦合界面的方法傳遞給第3節(jié)中建立的有限元模型，計(jì)算出螺旋槳的變形，將變形后的槳葉節(jié)點(diǎn)坐標(biāo)輸出到文本文件.

3)采用第5節(jié)中的方法確定變形后的螺旋槳幾何參數(shù)，建立變形后的螺旋槳幾何模型.

4)將變形后的幾何再次輸入到CFD軟件中求解螺旋槳的敞水性能和壓力分布.

5)判斷是否滿足平衡方程和收斂條件，如果滿足則輸出計(jì)算結(jié)果，包括結(jié)構(gòu)變形量、應(yīng)力場(chǎng)、推力、扭矩以及推進(jìn)效率等.如不滿足收斂條件，重復(fù)迭代2)、3)、4)步，直至結(jié)果收斂.

文中指定的收斂準(zhǔn)則為2個(gè)迭代步內(nèi)的推力系數(shù)和扭矩系數(shù)小于2%.計(jì)算過程一般迭代2次就可達(dá)到收斂.

7 應(yīng)力、變形計(jì)算及變形后的壓力分布

7.1 應(yīng)力水平和變形分布

圖9是不同轉(zhuǎn)速、不同進(jìn)速下的MISE應(yīng)力分布和變形對(duì)照?qǐng)D(變形均放大20倍).設(shè)計(jì)工況，n=10 r/s和 20 r/s時(shí)，最大等效應(yīng)力分別為2.3 MPa和8.62 MPa.由圖可以看出，相同轉(zhuǎn)速不同進(jìn)速情況下，低進(jìn)速情況下的MISE應(yīng)力比高進(jìn)速下的要大.不同轉(zhuǎn)速，相同進(jìn)速情況，轉(zhuǎn)速高的應(yīng)力要比轉(zhuǎn)速低的要大.最大MISE應(yīng)力分布均在葉根弦向中心位置.根據(jù)螺旋槳理論可知，相同轉(zhuǎn)速下，在低進(jìn)速時(shí)為螺旋槳的“重載”狀態(tài)，推力和扭矩均較大;高進(jìn)速時(shí)為“輕載”狀態(tài)，推力和扭矩較小，所以在低進(jìn)速時(shí)螺旋槳的變形要比高進(jìn)速的大，應(yīng)力分布也是一致的.相同進(jìn)速，不同轉(zhuǎn)速下，根據(jù)螺旋槳的無因次理論，推力和扭矩增大數(shù)為2的轉(zhuǎn)速倍數(shù)次方.彈性范圍應(yīng)力和變形也增加同樣倍數(shù).

圖10和圖11是不同轉(zhuǎn)速、進(jìn)速下的槳葉變形分布.從圖10可以看出，槳葉變形引起了槳葉側(cè)斜的改變.4119螺旋槳為無側(cè)斜槳，變形后側(cè)斜角為正值，即向隨邊傾斜.相同轉(zhuǎn)速時(shí)，低進(jìn)速的側(cè)斜改變比高進(jìn)速的大，這與前面分析的受力規(guī)律一致.

而圖11則顯示了在水動(dòng)力載荷的作用下，槳葉朝著船體方向變形，且在葉稍附近達(dá)到最大值，變形規(guī)律類似懸臂梁.變形后的槳葉剖面鼻尾線在螺旋槳軸向(x軸)方向的投影有零變?yōu)檎担磳?dǎo)致產(chǎn)生了縱傾的分布.另槳葉朝著船體方向變形，促使各半徑處的螺距角分布改變.即在隨邊附件變形大，導(dǎo)致螺旋槳槳葉螺距的變小.

圖9 不同轉(zhuǎn)速n、進(jìn)速J下的MISE應(yīng)力分布Fig.9 MISE stress distribution at different n and J

圖10 側(cè)斜的改變(線框?yàn)槲醋冃危瑢?shí)體為變形后)Fig.10 Changes of the blade skew

圖11 縱傾的改變(線框?yàn)槲醋冃危瑢?shí)體為變形后)Fig.11 Changes of the blade rake

7.2 變形前、后槳葉的壓力分布

圖12為設(shè)計(jì)進(jìn)速J=0.833時(shí)變形前后槳葉壓力面和吸力面的壓力分布對(duì)比圖.如圖12所示，變形的葉稍改變了葉稍附近壓力的分布.變形后正壓力峰值相應(yīng)減小，負(fù)壓力峰值相對(duì)升高.

圖12 J=0.833螺旋槳變形前后壓力分布(n=10 r/s)Fig.12 Predicted pressure distribution of 4119 in the undefand def-configuration，J=0.833，n=10 r/s

8 結(jié)論

1)對(duì)螺旋槳進(jìn)行穩(wěn)態(tài)流固耦合分析，SST κ-ω湍流模型預(yù)報(bào)精度高且求解時(shí)間合理.

2)實(shí)現(xiàn)了螺旋槳的流固耦合算法，可以求解槳葉的變形和應(yīng)力分布及變形后的推進(jìn)性能.

3)螺旋槳在水動(dòng)力作用下產(chǎn)生的變形，改變了螺旋槳的初始幾何參數(shù)，因此改變了螺旋槳的壓力分布，從而改變了螺旋槳的水動(dòng)力性能.對(duì)于易于變形的復(fù)合材料螺旋槳設(shè)計(jì)、計(jì)算，要考慮槳葉變形的影響.

［1］MOURITZ A，GELLERT E，BURCHILL P，et al.Review of advanced composite structures for naval ships and submarines［J］.Composite Structures，2001，53:21-41.

［2］高富東，潘存云，蔡汶珊，等.基于CFD的螺旋槳敞水性能數(shù)值分析與驗(yàn)證［J］.機(jī)械工程學(xué)報(bào)，2010，46(8): 133-139.

GAO Fudong，PAN Cunyun，CAI Wenshan，et al.Numerical analysis and validation of propeller open-water performance based on CFD［J］.Journal of Mechanical Engineering，2010，46(8):133-139.

［3］黃勝，王超，王詩洋.不同湍流模型在螺旋槳水動(dòng)力性能計(jì)算中的應(yīng)用與比較［J］.哈爾濱工程大學(xué)學(xué)報(bào)，2009，30(5):484-485.

HUANG Sheng，WANG Chao，WANG Shiyang.Application and comparison of different turbulence models in the computation of a propeller’s hydrodynamic performance［J］Journal of Harbin Engineering University，2009，30(5): 484-485.

［4］沈海龍，蘇玉民.基于滑移網(wǎng)格技術(shù)的船槳相互干擾研究［J］.哈爾濱工程大學(xué)學(xué)報(bào)，2010，31(1):1-7.

SHEN Hailong，SU Yumin.Use of the sliding mesh technique to analyze interaction between ship hulls and propellers［J］.Journal of Harbin Engineering University，2010，31 (1):1-7.

［5］王超，黃勝，單鐵兵.基于多塊混合網(wǎng)格方法預(yù)報(bào)螺旋槳非正常工作狀態(tài)時(shí)的水動(dòng)力性能［J］.船舶力學(xué)，2010，14(1/2):51-55.

WANG Chao，HUANG Sheng，SHAN Tiebing.Computations of the propeller’s hydrodynamic performance in abnormal working condition based on multi-block meshes［J］.Journal of Ship Mechanics，2010，14(1/2):51-55.

［6］LIN H J，LIN J J.Nonlinear hydro-elastic behavior of propellers using a finite-element method and lifting surface theory［J］.Journal of Marine Science and Technology，1996，1:114-124.

［7］YOUNG Y L.Fluid-structure interaction analysis of flexible composite marine propellers［J］.Journal of Fluids and Structures，2008，24:799-818.

［8］李巍，王國(guó)強(qiáng)，汪蕾.螺旋槳粘流水動(dòng)力特性數(shù)值模擬［J］.上海交通大學(xué)學(xué)報(bào)，2007，41(7):1200-1203.

LI Wei，WANG Guoqiang，WANG Lei.The numerical simulation of hydrodynamics characteristic in propeller［J］.Journal of Shanghai Jiaotong University，2007，41(7):1200-1203.

［9］JESSUP S D.An experimental investigation of viscous aspects of propeller［D］.Washington DC:The Catholic University of America，1989:65-154.

［10］王國(guó)強(qiáng)，董世湯.船舶螺旋槳理論與應(yīng)用［M］.哈爾濱:哈爾濱工程大學(xué)出版社，2007:89-94，113-114.