基于懸臂梁法和面元法耦合的槳葉應(yīng)力分布預(yù)報(bào)＊

葉禮裕 王 超 孫 帥 孫文林

（哈爾濱工程大學(xué)船舶工程學(xué)院 哈爾濱 150001）

0 引 言

隨著現(xiàn)代艦船逐步向大型化、高速化發(fā)展，以及大功率主機(jī)的應(yīng)用，引起螺旋槳槳葉表面負(fù)荷增大；由于減振降噪方面有特殊的優(yōu)勢，大側(cè)斜螺旋槳在船舶領(lǐng)域得到廣泛的應(yīng)用［1］，但特殊的幾何外形使其強(qiáng)度問題較為突出.另外，為了使螺旋槳實(shí)際應(yīng)用中更加節(jié)能化和環(huán)保化，設(shè)計(jì)槳應(yīng)當(dāng)符合輕量化的要求.因此，開展螺旋槳強(qiáng)度的研究不僅在于槳葉強(qiáng)度符合規(guī)范要求，更應(yīng)當(dāng)能夠準(zhǔn)確地預(yù)報(bào)槳葉應(yīng)力的分布.

目前，螺旋槳槳葉應(yīng)力的預(yù)報(bào)通常采用有限元法和模型試驗(yàn)［2］.對于有限元法，開展較多的是采用CFD計(jì)算和有限元分析軟件結(jié)合的方法預(yù)報(bào)槳葉應(yīng)力分布［3－5］，也有學(xué)者基于面元法和有限元法耦合的方法預(yù)報(bào)槳葉應(yīng)力分布［6］.在模型試驗(yàn)方面，Boswell［7］對大側(cè)斜螺旋槳葉片進(jìn)行了靜態(tài)應(yīng)力測量試驗(yàn)；趙波［8］開展了大側(cè)斜螺旋槳的靜態(tài)應(yīng)力試驗(yàn)和動態(tài)應(yīng)力試驗(yàn)，并與理論計(jì)算結(jié)果進(jìn)行對比；楊向暉等［9］通過在大側(cè)斜槳模表面粘貼應(yīng)變片的方法研究不同工況下的槳葉應(yīng)變和應(yīng)力分布.有限元法能夠較為準(zhǔn)確預(yù)報(bào)螺旋槳槳葉應(yīng)力分析，但是需要復(fù)雜的模型建立、網(wǎng)格劃分過程，不利于開展螺旋槳的設(shè)計(jì).模型試驗(yàn)方法成本較高、實(shí)驗(yàn)難度大且需要花費(fèi)大量時間.螺旋槳運(yùn)轉(zhuǎn)時主要受到周向力、軸向力、離心力，以及偶然外力作用，從而產(chǎn)生彎曲、轉(zhuǎn)矩和拉伸變形，可將槳葉簡化為變截面扭曲的懸臂梁［10］.懸臂梁法也是一種比較方便且可行的槳葉應(yīng)力預(yù)報(bào)方法.以往學(xué)者開展懸臂梁預(yù)報(bào)螺旋槳應(yīng)力，因螺旋槳幾何外形復(fù)雜且橫截面為非對稱的，計(jì)算截面所受彎矩、離心力，以及面積、抗彎模量等的過程中運(yùn)用了一些經(jīng)驗(yàn)公式，造成只能預(yù)報(bào)槳葉有限個點(diǎn)的應(yīng)力或者應(yīng)力預(yù)報(bào)結(jié)果，比較保守.

為準(zhǔn)確預(yù)報(bào)槳葉各半徑切面所受的彎矩和離心力，將懸臂梁法理論中的由螺旋槳推力、旋轉(zhuǎn)阻力產(chǎn)生的彎矩、離心力及其產(chǎn)生的彎矩的積分計(jì)算公式進(jìn)行離散處理.槳葉切面的形心和慣性矩可借助于工程力學(xué)中截面靜矩和慣性矩等計(jì)算公式來求解，以準(zhǔn)確求解槳葉各半徑切面上各個點(diǎn)的抗彎剖面模數(shù).本文采用FORTRAN語言編譯懸臂梁法槳葉應(yīng)力預(yù)報(bào)程序，將其與螺旋槳的定常面元法性能預(yù)報(bào)程序?qū)?以無側(cè)斜槳和大側(cè)斜螺旋槳為例，通過本文計(jì)算結(jié)果和有限元計(jì)算結(jié)果對比分析來驗(yàn)證本文所提方法的有效性.

1 懸臂梁離散方法

懸臂梁法預(yù)報(bào)螺旋槳強(qiáng)度是將螺旋槳葉梢視為自由端，而葉根固定于槳轂的一根變截面的懸臂梁，一般按正車最大航速工況預(yù)報(bào)螺旋槳靜載荷，將作用在螺旋槳上的外力載荷視為不變.以平切面作為校核面，該校核面為垂直于槳葉參考線的面且與槳葉相交所截的平面，對大側(cè)斜槳有更好的適用性［11］.

槳葉切面各點(diǎn)抗彎彈性模量計(jì)算過程如下：首先，對該切面進(jìn)行分割，計(jì)算出其面積和靜矩，進(jìn)而求出該剖面的形心坐標(biāo)；然后，計(jì)算出該剖面的慣性矩；最后，由該剖面的形心坐標(biāo)和慣性矩，可求處其各點(diǎn)的抗彎彈性模量.具體的切面形心、靜矩、慣性矩和抗彎彈性模量求解方法可參考文獻(xiàn)［12］.

離散法可準(zhǔn)確地計(jì)算槳葉各半徑切面的彎矩.根據(jù)離散法原理，需對求解的螺旋槳幾何模型進(jìn)行分割.為此，將槳葉沿展向分割成有限個區(qū)域，半徑r處取一葉元體，厚度為dr.為更好對接面元法預(yù)報(bào)程序，懸臂梁法的槳葉展向劃分方式與面元法相同，即采用余弦分割，節(jié)點(diǎn)為

1.1 推力和旋轉(zhuǎn)阻力引起的彎矩積分公式離散方法

對于常規(guī)螺旋槳，推力和旋轉(zhuǎn)阻力是產(chǎn)生彎矩的主要成分.設(shè)該葉元體所受的推力為dT和旋轉(zhuǎn)阻力dF，則槳葉推力和旋轉(zhuǎn)阻力對半徑處的彎矩分布為

為獲得上式的數(shù)值解，將該式離散處理，沿展向分割成Nr份，則

為計(jì)算推力和旋轉(zhuǎn)阻力產(chǎn)生的彎矩，必須計(jì)算出ΔTj和ΔFj大小，而ΔTj和ΔFj可通過面元法計(jì)算得到.面元法計(jì)算螺旋槳總的推力和轉(zhuǎn)矩的離散公式為［13－14］

面元法是基于勢流理論的分析方法，計(jì)算中沒有考慮流體粘性的影響.為使結(jié)果更符合實(shí)際情況，計(jì)算水動力性能時，采用表面摩擦阻力的方法近似計(jì)入粘性的影響，用以下公式求解推力和扭矩的粘性修正量.

只要對面元法計(jì)算螺旋槳推力和扭矩的離散公式稍加處理，即可用于求解每個葉元體的推力和旋轉(zhuǎn)阻力，即

1.2 離心力及其產(chǎn)生的彎矩積分公式離散方法

螺旋槳高速運(yùn)轉(zhuǎn)時，槳葉各徑向剖面將受到離心力的作用.離心力引起的彎矩相對于推力和旋轉(zhuǎn)阻力要小，但有側(cè)斜和縱傾的螺旋槳由于離心力作用線已不是徑向剖面幾何中心的連線，此時產(chǎn)生的離心力彎矩較大.

假設(shè)螺旋槳的轉(zhuǎn)速為n（r／s），材料的重量密度為γ（10kN／m3），計(jì)算半徑處的切面面積為S（m2），重力加速度為g（9.8m／s2），則取一小微元dr，它對半徑rP處的離心力為

與推力和旋轉(zhuǎn)阻力產(chǎn)生彎矩?cái)?shù)值計(jì)算方法類似，將上述進(jìn)行離散處理，得到

式中：Sj為切面的面積.

對于有縱傾的螺旋槳，將因槳葉縱傾產(chǎn)生離心力彎矩MC.假設(shè)在不同半徑處槳葉縱傾分布為h（r），下式縱傾引起的槳葉rP處切面彎矩MC及其離散結(jié)果

對于有側(cè)斜的螺旋槳，將因槳葉側(cè)斜引起離心彎矩MS.通常，螺旋槳側(cè)斜的方向是與其旋轉(zhuǎn)方向相反的.因此，側(cè)斜引起的離心力彎矩MS與縱傾引起的彎矩MF的方向相反.假設(shè)各半徑側(cè)斜分布為p（r），則側(cè)斜引起的槳葉rp處切面彎矩MC及其離散結(jié)果為

1.3 切面上各點(diǎn)應(yīng)力求解

針對半徑rp處承受的彎矩情況進(jìn)行分析，見圖1.將推力、旋轉(zhuǎn)阻力和離心力引起的彎矩沿ξ－ξ軸和η－η軸方向分解，可以得到合成彎矩.

圖1 槳葉剖面承受彎矩情況

校核切面的面積及切面上任何一點(diǎn)的抗彎剖面模數(shù)的求解方法已在以上作簡要介紹.假設(shè)P（ξ，η）為切面上的一點(diǎn)，以Wξ（η）和Wη（ξ）分別表示P點(diǎn)對于ξ－ξ軸和η－η軸抗彎截面模量.那么，合成彎矩Mξ和Mη在P點(diǎn)處的應(yīng)力可表示為：

若P（ξ，η）點(diǎn)所在校核切面的面積為S，則離心力C在P點(diǎn)產(chǎn)生的應(yīng)力為

綜述所述，合成彎矩和離心力所產(chǎn)生的總應(yīng)力為

式中：σξ（η），ση（η）正負(fù)號取值需要根據(jù)具體情況而定，即ξ－ξ軸和η－η軸把校核剖面分割為4個區(qū)域，表1給出了每個區(qū)域受力狀態(tài)相同.若為拉應(yīng)力，則取正值.反之，取負(fù)值.

表1 槳葉剖面各區(qū)域的應(yīng)力性質(zhì)

2 算例分析

根據(jù)以上介紹的懸臂梁離散方法，筆者編譯了懸臂梁法預(yù)報(bào)程序，并與定常面元法水動力性能預(yù)報(bào)程序進(jìn)行對接，以便更好地把螺旋槳水動力性能預(yù)報(bào)傳入到懸臂梁法程序計(jì)算中.程序運(yùn)行計(jì)算后，可將槳葉的應(yīng)力預(yù)報(bào)結(jié)果傳入到Tecplot軟件中進(jìn)行后處理，得到槳葉應(yīng)力分布云圖.為了證明本文方法的可行性和有效性，將本文方法與有限元軟件的槳葉應(yīng)力的預(yù)報(bào)結(jié)果進(jìn)行比較.

這里需要強(qiáng)調(diào)的是，一般情況下，任意一點(diǎn)的主單元體的應(yīng)力狀態(tài)是三向的應(yīng)力.懸臂梁法是以最大拉壓應(yīng)力σ1為校核條件的，將其他2項(xiàng)應(yīng)力σ2和σ3忽略不計(jì)，而有限元法計(jì)算的以von Mises應(yīng)力（等效應(yīng)力）σε為校核條件的，見文獻(xiàn)［15］，σε計(jì)算方法見式（20）.根據(jù)螺旋槳的受力情況可知，槳葉上任意一點(diǎn)的最大拉壓應(yīng)力σ1遠(yuǎn)大于其他2項(xiàng)應(yīng)力σ2和σ3.因此，槳葉任意一點(diǎn)的最大拉壓應(yīng)力σ1的絕對值與von Mises應(yīng)力σs的大小相差不大.

本文選了2種類型的螺旋槳進(jìn)行算例分析，包括無縱傾無側(cè)斜的P4119、P4381槳，以及某大側(cè)斜調(diào)距槳.

2.1 P4119槳和P4381槳的模型算例分析

對于無縱傾無側(cè)斜槳的算例，本文選用22nd ITTC Propulsion Committee提供的P4119槳和P4381槳為算例，其幾何參數(shù)見表2.

表2 P4119槳和P4381的主要參數(shù)

為了驗(yàn)證本文方法預(yù)報(bào)螺旋槳強(qiáng)度的準(zhǔn)確性，將本文的應(yīng)力預(yù)報(bào)結(jié)果與文獻(xiàn)［16］進(jìn)行比較.預(yù)報(bào)P4119槳的應(yīng)力分布的工況為：進(jìn)速系數(shù)0.833，轉(zhuǎn)速600r／min.在該工況下.本文面元法程序和文獻(xiàn)［16］預(yù)報(bào)的P4119槳水動力性能見表3，可見2種的誤差較小.

表3 P4119槳的推力和轉(zhuǎn)矩系數(shù)對比

圖2為本文方法預(yù)報(bào)DTRC4119槳的葉背和葉面的應(yīng)力分布.文獻(xiàn)［16］是采用有限元計(jì)算軟件計(jì)算DTRC4119槳槳葉應(yīng)力分布，預(yù)報(bào)結(jié)果見圖3.由圖2～3可知，本文計(jì)算的槳葉的應(yīng)力分布趨勢與文獻(xiàn)［16］計(jì)算結(jié)果類似.本文方法預(yù)報(bào)槳葉最大拉應(yīng)力為，而文獻(xiàn)［16］預(yù)報(bào)最大von Mises應(yīng)力為，兩者的最大應(yīng)力預(yù)報(bào)結(jié)果偏差較小.另外，注意到文獻(xiàn)［16］提到其最大應(yīng)力發(fā)生在槳葉葉根弦向中部處，而本文預(yù)報(bào)結(jié)果最大應(yīng)力也是發(fā)生在葉根弦向中部處.

圖2 本文方法預(yù)報(bào)的P4119槳槳葉應(yīng)力分布

圖3 文獻(xiàn)［16］預(yù)報(bào)P4119槳槳葉應(yīng)力分布

P4381槳的槳葉應(yīng)力計(jì)算工況為進(jìn)速2.711 45m／s、轉(zhuǎn)速600r／min.P4381槳的有限元軟件應(yīng)力預(yù)報(bào)步驟如下：首先，在ICEM軟件中建立P4381槳的幾何模型并對其劃分網(wǎng)格，導(dǎo)入到Fluent軟件中進(jìn)行水動力計(jì)算；將P4381槳的幾何模型導(dǎo)入到ANSYS結(jié)構(gòu)力學(xué)模塊中進(jìn)行結(jié)構(gòu)網(wǎng)格劃分，將P4381槳的水動力載荷傳入到調(diào)距槳該槳的有限元模型中，計(jì)算出P4381槳葉面和葉背的應(yīng)力分布.

在計(jì)算工況下，本文的定常面元法程序和CFD方法預(yù)報(bào)得到的該槳推力系數(shù)和轉(zhuǎn)矩系數(shù)見表4，2種方法得到的該槳推力系數(shù)和轉(zhuǎn)矩系數(shù)誤差較小.

表4 P4381槳的推力和扭矩系數(shù)對比

運(yùn)行懸臂梁法預(yù)報(bào)程序，得到了P4381槳的葉背和葉面的應(yīng)力分布情況，見圖4.有限元法預(yù)報(bào)P4381槳的葉背和葉面的應(yīng)力分布見圖5.從葉背應(yīng)力預(yù)報(bào)結(jié)果分析，本文方法和有限元法預(yù)報(bào)的P4381槳葉背的分布趨勢基本是一致的；本文方法預(yù)報(bào)葉背最大壓應(yīng)力為3.21MPa，而有限元法預(yù)報(bào)葉背最大von Mises應(yīng)力為3.61MPa，2種方法計(jì)算得到葉背最大應(yīng)力值相近且均是發(fā)生于葉根弦向中部處.從葉面應(yīng)力預(yù)報(bào)結(jié)果分析，本文方法和有限元法預(yù)報(bào)的P4381槳葉面的分布趨勢有一定的差別，但是葉根處的分布趨勢比較相近；本文方法預(yù)報(bào)葉面最大拉應(yīng)力為3.16 MPa，而有限元法預(yù)報(bào)葉面最大von Mises應(yīng)力為3.61MPa，2種方法計(jì)算得到葉面最大應(yīng)力值相近且均是發(fā)生于葉根弦向中部處.

圖4 本文方法預(yù)報(bào)P4381槳槳葉應(yīng)力分布

圖5 有限元法預(yù)報(bào)P4381槳槳葉應(yīng)力分布

綜上所述，本文方法可用于預(yù)報(bào)常規(guī)槳的應(yīng)力分布趨勢，能夠準(zhǔn)確預(yù)報(bào)常規(guī)槳的最大應(yīng)力值及其發(fā)生的位置.但是，在預(yù)報(bào)常規(guī)槳葉面的外半徑處的應(yīng)力分布有所缺陷，主要原因是葉面外半徑處葉切面比較平坦，導(dǎo)致同一半徑處的葉切面處的應(yīng)力值趨向相同.

2.2 大側(cè)斜槳算例分析

本文大側(cè)斜螺旋的計(jì)算模型選的是某大側(cè)斜調(diào)距槳，主要參數(shù)見表5.

表5 大側(cè)斜槳的主要參數(shù)

該槳的計(jì)算工況：進(jìn)流速度為11.8m／s，轉(zhuǎn)速為3.4r／s.有限元法計(jì)算該槳的步驟與P4119槳和P4381槳的模型算例一樣.

從敞水性能預(yù)報(bào)的準(zhǔn)確性分析，表6給出了計(jì)算工況下本文面元法和CFD方法計(jì)算得到的調(diào)距槳的敞水性能，兩者的誤差很小，即不會由于水動力性能計(jì)算結(jié)果的不同而引起懸臂梁法和有限元法預(yù)報(bào)該槳的槳葉應(yīng)力分布有較大的差別.

表6 大側(cè)斜槳的敞水性能對比

本文方法預(yù)報(bào)該槳的槳葉應(yīng)力分布云圖見圖6.由圖6可知，本文方法預(yù)報(bào)的槳葉最大壓應(yīng)力為83.6MPa，發(fā)生于葉背中部的弦向中部處，槳葉壓應(yīng)力也主要集中于葉背中部弦向中部處；最大拉應(yīng)力為92.7MPa，發(fā)生于葉面中部的隨邊處，槳葉的拉應(yīng)力主要集中與槳葉中部的導(dǎo)邊和隨邊處.有限元軟件預(yù)報(bào)該調(diào)距槳的槳葉應(yīng)力分布云圖，見圖7.由圖7可知，有限元法預(yù)報(bào)的葉背最大von Mises應(yīng)力為83.2MPa，發(fā)生于葉背中部的弦向中部靠近隨邊處，槳葉的拉應(yīng)力也集中于此處；葉面最大von Mises應(yīng)力為91.6 MPa，發(fā)生于葉面中部的弦向隨邊處.2種方法預(yù)報(bào)大側(cè)斜槳葉背和葉面的最大應(yīng)力大小差別小，且發(fā)生位置比較相近，驗(yàn)證了懸臂梁法在計(jì)算大側(cè)斜槳槳葉應(yīng)力分布的可行性和準(zhǔn)確性.但是，2種方法預(yù)報(bào)的槳葉應(yīng)力分布趨勢有所差別，即懸臂梁法預(yù)報(bào)拉應(yīng)力集中于葉面中部的導(dǎo)邊和隨邊處，而有限元法預(yù)報(bào)拉應(yīng)力只集中于葉面中部的隨邊處，主要原因是：由于懸臂梁法的彎矩計(jì)算原理的局限性，導(dǎo)致計(jì)算得到同一徑向處不同弦向位置的彎矩是相同的，而有限元法計(jì)算得到同一徑向處不同弦向位置的彎矩是不同的.

圖6 本文方法預(yù)報(bào)的大側(cè)斜槳槳葉的應(yīng)力分布

圖7 有限元法預(yù)報(bào)大側(cè)斜槳槳葉的應(yīng)力分布

4 結(jié) 論

由于螺旋槳的幾何外形復(fù)雜，懸臂梁法理論螺旋槳推力、旋轉(zhuǎn)阻力產(chǎn)生的彎矩、離心力及其產(chǎn)生的彎矩的積分公式很難用解析的方法求出.為準(zhǔn)確計(jì)算出槳葉切面的彎矩和離心力，本文建立了懸臂梁離散方法用于求解槳葉切面所受彎矩和離心力的.基于FORTRAN語言編譯懸臂梁法槳葉應(yīng)力預(yù)報(bào)程序，并與螺旋槳定常面元法性能預(yù)報(bào)程序?qū)?以無縱傾無側(cè)斜槳和大側(cè)斜螺旋槳為例，預(yù)報(bào)槳葉應(yīng)力分布，將計(jì)算結(jié)果與有限元法計(jì)算結(jié)果相比較，分析得到如下結(jié)論.

1）本文方法能較為準(zhǔn)確的預(yù)報(bào)無縱傾無側(cè)斜槳的葉面和葉背最大應(yīng)力值，且對于該類型的槳最大應(yīng)力發(fā)生于槳葉葉根的弦向中部.除了預(yù)報(bào)常規(guī)槳葉面的外半徑處的應(yīng)力分布有所缺陷，本文方法預(yù)報(bào)常規(guī)槳的槳葉的應(yīng)力分布趨勢與有限元法預(yù)報(bào)趨勢是相同.

2）本文方法可較為準(zhǔn)確預(yù)報(bào)大側(cè)斜槳的葉面和葉背最大應(yīng)力值.對于本算例的平衡式側(cè)斜槳，最大拉應(yīng)力一般發(fā)生于槳葉中部的隨邊處，最大壓應(yīng)力發(fā)生于槳葉中部的弦向中部處.懸臂梁法預(yù)報(bào)槳葉應(yīng)力分布趨勢有缺陷，即懸臂梁法預(yù)報(bào)螺旋槳的拉應(yīng)力在葉面中部的導(dǎo)邊和隨邊處有集中現(xiàn)象，而實(shí)際上對于該槳拉應(yīng)力只集中于葉面中部的隨邊處.

總而言之，應(yīng)用本文方法可快速得找到在計(jì)算工況下的螺旋槳槳葉最大應(yīng)力大小及其發(fā)生位置，便于對螺旋槳的強(qiáng)度校核.雖然本文方法預(yù)報(bào)在槳葉應(yīng)力分布趨勢與有限元法還存在某些缺陷，但是總體上還是可行的.在后續(xù)的工作中，作者將開展將該方法的程序嵌入到螺旋槳的理論設(shè)計(jì)和優(yōu)化設(shè)計(jì)過程中，以保證設(shè)計(jì)槳有足夠的強(qiáng)度.

［1］CUMMING R A，MORGAN W B，BOSWELL R J.Highly skewed propellers［J］.SNAME Transactions，1984（1）：90：98.

［2］S?NTVEDT T.Propeller blade stresses，application of finite element methods［J］.Computers ＆ Structures，1974，4（1）：193－204.

［3］文學(xué)棟，王永生，李堅(jiān)波.船用螺旋槳槳葉應(yīng)力數(shù)值計(jì)算［J］.船海工程，2010，39（1）：27－30.

［4］柳章存.船用螺旋槳的數(shù)值模擬及流固分析［D］.鎮(zhèn)江：江蘇科技大學(xué)，2012.

［5］任 弘，李范春，杜 玲.流固耦合作用對螺旋槳強(qiáng)度影響的數(shù)值計(jì)算［J］.武漢理工大學(xué)學(xué)報(bào)：交通科學(xué)與工程版，2015，29（1）：144－147.

［6］劉竹青，陳奕宏，姚志崇.基于面元法及有限元法耦合的螺旋槳強(qiáng)度計(jì)算［J］.中國造船，2012，53（1）：25－30.

［7］BOSWELL R J.Static stress measurements on a highly－skewed propeller blade［R］.Naval Ship Research and Development Center Washington DC，1969.

［8］趙 波.大側(cè)斜螺旋槳強(qiáng)度研究［D］.武漢：華中科技大學(xué)，2003.

［9］楊向暉，程爾升.大側(cè)斜螺旋槳槳葉應(yīng)力測試研究［J］.船海工程，2004（1）：6－9.

［10］盛振邦，劉應(yīng)中.船舶原理：下冊［M］.上海：上海交通大學(xué)出版社，2004.

［11］高同兵，程爾升.大側(cè)斜螺旋槳設(shè)計(jì)中的考慮［J］.船海工程，2004（2）：1－3.

［12］楊在林.工程力學(xué)［M］.哈爾濱：哈爾濱工程大學(xué)出版社，2010.

［13］蘇玉民，黃 勝.船舶螺旋槳理論［M］.哈爾濱：哈爾濱工程大學(xué)出版社，2003.

［14］楊在林，宋天舒，楊 勇.材料力學(xué)：下冊［M］.哈爾濱：哈爾濱工程大學(xué)出版社，2011.

［15］KOYAMA K.Comparative calculations of propellers by surface panel method－workshop organized by 20th ITTC propulsion committee［J］.Papers of Ship Research Institute，1993，16（S）：57－66.

［16］顧鋮璋，鄭百林.船用螺旋槳敞水性能與槳葉應(yīng)力的數(shù)值分析［J］.力學(xué)季刊，2011，32（3）：440－443.