基于VBA 的螺旋槳流固耦合數(shù)值模擬＊

賀偉 李廷秋 李子如

（武漢理工大學(xué)交通學(xué)院 武漢 430063）

0 引 言

對(duì)于采用鎳鋁青銅或鎳錳青銅等銅合金材料的船用螺旋槳，材料彈性模量大，槳葉運(yùn)轉(zhuǎn)過(guò)程中變形量小，因而在槳葉設(shè)計(jì)／分析中常假設(shè)螺旋槳為剛性，忽略槳葉結(jié)構(gòu)變形與水動(dòng)力之間的相互影響.但對(duì)于遠(yuǎn)離設(shè)計(jì)工況（如系柱或倒車）的重負(fù)荷螺旋槳、大側(cè)斜螺旋槳，以及復(fù)合材料或塑料螺旋槳，葉片在水動(dòng)力和離心力作用下將會(huì)產(chǎn)生明顯的變形，葉片變形反過(guò)來(lái)又會(huì)對(duì)周圍的流場(chǎng)產(chǎn)生影響，這是一個(gè)典型的流固耦合問(wèn)題［1］.尤其是近年來(lái)被反復(fù)提及的復(fù)合材料螺旋槳水彈性自適應(yīng)性［2-5］，其機(jī)理就在于利用非各向同性復(fù)合材料的彎曲-扭轉(zhuǎn)耦合特性，螺旋槳螺距能根據(jù)負(fù)荷的輕重產(chǎn)生自適應(yīng)行為，使得復(fù)合材料螺旋槳在偏離設(shè)計(jì)點(diǎn)或船后非均勻流場(chǎng)中具有優(yōu)于金屬螺旋槳的水動(dòng)力表現(xiàn)，此時(shí)常規(guī)基于剛性假設(shè)將水動(dòng)力載荷和結(jié)構(gòu)響應(yīng)分開(kāi)考慮的研究思路顯然已無(wú)法勝任，而必須在統(tǒng)一的框架下考慮水動(dòng)力載荷和葉片變形的相互影響，即開(kāi)展螺旋槳水彈性研究.

目前，船用螺旋槳水彈性數(shù)值分析多基于弱耦合模型展開(kāi)，求解過(guò)程中流體和結(jié)構(gòu)計(jì)算相互獨(dú)立，通過(guò)迭代來(lái)考慮流體和結(jié)構(gòu)的相互影響，計(jì)算模塊的完整性得以保持.其中，結(jié)構(gòu)計(jì)算幾乎都采用有限元分析方法，不同之處在于水動(dòng)力分析模型是采用勢(shì)流模型［6-8］還是粘性模型［9-10］.考慮到勢(shì)流理論在螺旋槳水動(dòng)力分析方面的成熟性及計(jì)算效率，本文采用面元法進(jìn)行螺旋槳水動(dòng)力分析，基于Mechanical APDL（ANSYS）進(jìn)行槳葉結(jié)構(gòu)分析，借助Excel VBA 進(jìn)行流場(chǎng)與結(jié)構(gòu)計(jì)算的迭代，實(shí)現(xiàn)了螺旋槳穩(wěn)態(tài)性能的雙向流固耦合數(shù)值模擬.

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

1.1 螺旋槳水動(dòng)力分析模型

螺旋槳水動(dòng)力分析采用基于擾動(dòng)速度勢(shì)的面元法進(jìn)行.在固定于螺旋槳上的運(yùn)動(dòng)坐標(biāo)系中，物面上擾動(dòng)速度勢(shì)φ滿足如下積分方程：

式中：SB為包含槳葉表面和槳轂表面的物面；SW為尾渦面；Δφ為尾流面上偶極子分布；r為場(chǎng)點(diǎn)到源點(diǎn)的距離；n 為物面或尾渦面的法向；?為Cauchy主值積分.

式（1）的定解條件包括物面不可穿透條件和隨邊等壓Kutta條件：

式中：Vn為法向速度；▽in為來(lái)流速度；Δp為槳葉葉背和葉面壓力差.

槳葉和槳轂表面離散為四邊形雙曲面元，螺旋槳尾渦面采用經(jīng)驗(yàn)?zāi)Ｐ头譃榻擦鲄^(qū)和遠(yuǎn)尾流區(qū)，也采用四邊形雙曲面元離散.物面和尾渦面面元上均布置等強(qiáng)度偶極子，其強(qiáng)度通過(guò)求解式（1）～（3）可惟一確定.對(duì)擾動(dòng)速度勢(shì)微分可得擾動(dòng)速度，結(jié)合伯努利方程可得表面壓力，粘性貢獻(xiàn)則通過(guò)經(jīng)驗(yàn)?zāi)Ｐ痛_定［11］.

1.2 槳葉結(jié)構(gòu)分析模型

螺旋槳槳葉考慮為根部約束的線彈性結(jié)構(gòu).同樣在固定于螺旋槳上的運(yùn)動(dòng)坐標(biāo)系中，槳葉結(jié)構(gòu)運(yùn)動(dòng)方程為

式中：M，C和K分別為結(jié)構(gòu)的質(zhì)量矩陣、阻尼矩陣和剛度矩陣；和u分別為節(jié)點(diǎn)加速度、速度和位移矢量；Fce，F(xiàn)co和Fh分別為離心力，科氏力和水動(dòng)力載荷.對(duì)于瞬態(tài)問(wèn)題，水動(dòng)力載荷同時(shí)包含剛性槳運(yùn)動(dòng)和槳葉彈性變形的貢獻(xiàn).對(duì)于穩(wěn)態(tài)問(wèn)題，外載荷與結(jié)構(gòu)內(nèi)力相平衡且不隨時(shí)間變化，結(jié)構(gòu)分析中僅需考慮結(jié)構(gòu)變形量的影響，而不用考慮結(jié)構(gòu)變形速度的影響，因此結(jié)構(gòu)運(yùn)動(dòng)方程中可省去所有與加速度和速度有關(guān)的項(xiàng)，槳葉彈性變形對(duì)水動(dòng)力載荷的貢獻(xiàn)也可省去，槳葉結(jié)構(gòu)運(yùn)動(dòng)方程退化為

式中：Fr為剛性槳運(yùn)動(dòng)引起的水動(dòng)力載荷，通過(guò)前面介紹的水動(dòng)力模型確定.在研究的初始階段，本文首先關(guān)注穩(wěn)態(tài)問(wèn)題.

式（5）利用有限元軟件ANSYS求解，將材料屬性、單元類型、有限元節(jié)點(diǎn)、水動(dòng)力載荷等信息按ANSYS提供的APDL 語(yǔ)言格式寫(xiě)入數(shù)據(jù)文件，并通過(guò)ANSYS批處理讀入該數(shù)據(jù)文件可自動(dòng)進(jìn)行結(jié)構(gòu)分析，從而避免了復(fù)雜的交互式建模過(guò)程，為流固耦合迭代計(jì)算自動(dòng)化的實(shí)現(xiàn)提供了可能.

1.3 螺旋槳流固耦合實(shí)現(xiàn)

在Excel VBA 環(huán)境中，借助于Shell函數(shù)分別調(diào)用BEM 模塊和FEM 模塊，可實(shí)現(xiàn)螺旋槳的流固耦合分析，調(diào)用方式如下.

其中：BEM.exe為自行開(kāi)發(fā)的水動(dòng)力分析可執(zhí)行文件；〈InputFile〉為前文提及按APDL 語(yǔ)言格式生成的結(jié)構(gòu)計(jì)算輸入數(shù)據(jù)文件，〈OutputFile〉為結(jié)構(gòu)分析結(jié)果輸出文件.

圖1給出了螺旋槳穩(wěn)態(tài)流固耦合分析流程圖.選用推力系數(shù)kt（表征水動(dòng)力分析）和最大變形umax（表征結(jié)構(gòu)分析）這2個(gè)參數(shù)來(lái)判斷迭代計(jì)算是否收斂，收斂準(zhǔn)則定義為

式中：ε為小量，在本文研究中取為1×10-3.

圖1 流固耦合分析流程圖

2 實(shí)例分析

選用DTMB4119槳對(duì)本文提出的螺旋槳流固耦合數(shù)值模擬平臺(tái)進(jìn)行數(shù)值測(cè)試，該槳為3葉無(wú)縱傾無(wú)側(cè)斜槳，盤(pán)面比為Ae／Ao=0.6，直徑為D=0.3048m，槳葉具體型值見(jiàn)文獻(xiàn)［12］.轉(zhuǎn)速設(shè)為n=10r／s，不同進(jìn)速系數(shù)通過(guò)改變來(lái)流速度實(shí)現(xiàn).螺旋槳材料分別選用“剛性”鋁青銅和“彈性”S玻璃纖維，均為各項(xiàng)同性材料，材料參數(shù)見(jiàn)表1.

表1 螺旋槳材料參數(shù)

2.1 非耦合計(jì)算結(jié)果

首先選用鋁青銅材料螺旋槳針對(duì)設(shè)計(jì)進(jìn)速系數(shù)工況J=0.833進(jìn)行了非耦合計(jì)算（或稱為單向耦合），以初步驗(yàn)證各個(gè)模塊的計(jì)算精度以及網(wǎng)格合理性.

圖2～3分別給出了水動(dòng)力分析模型和結(jié)構(gòu)分析模型，槳葉表面網(wǎng)格采用展向均勻弦向余弦劃分方式，2個(gè)模型中坐標(biāo)系均定義為x軸指向螺旋槳后方為正，y軸垂直向上為正，z軸完成右手定則.

圖4給出了不同徑向位置剖面壓力分布計(jì)算結(jié)果與試驗(yàn)值［13］比較.其中，壓力系數(shù)選用直徑和轉(zhuǎn)速進(jìn)行量綱一的量化處理，即

式中：p，p0，ρ分別為槳葉表面壓力，參考?jí)毫土黧w密度.與通?；谄拭嫦鄬?duì)來(lái)流速度VR（r）進(jìn)行無(wú)量綱相比，前者可更直觀的反映槳葉不同位置處壓力相對(duì)大小關(guān)系.可以看出，在遠(yuǎn)離槳轂的外半徑0.7R和0.9R處計(jì)算值與試驗(yàn)值吻合較好，而在靠近槳轂的0.3R處誤差稍大，但壓力差基本相當(dāng)，且內(nèi)半徑載荷的力臂相對(duì)短，認(rèn)為該誤差對(duì)結(jié)構(gòu)彎曲變形和扭轉(zhuǎn)變形均影響不大.

圖2 水動(dòng)力分析模型

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

圖4 不同徑向位置剖面壓力分布計(jì)算值與試驗(yàn)值比較

圖5給出了槳葉葉背和葉面的壓力分布云圖，總體而言，槳葉外半徑承擔(dān)更多的負(fù)荷，且導(dǎo)邊20%弦長(zhǎng)部分壓差明顯要高于隨邊20%弦長(zhǎng)部分壓差.

槳葉的變形特征及應(yīng)力分布特征在下一部分與耦合計(jì)算結(jié)果統(tǒng)一分析比較.

2.2 流固耦合計(jì)算結(jié)果

圖6～9依次給出了設(shè)計(jì)進(jìn)速系數(shù)下分別采用鋁青銅和S玻纖材料流固耦合與非耦合壓力分布、槳葉變形輪廓、槳葉軸向變形量和槳葉等效應(yīng)力分布對(duì)比.其中，變形量放大倍數(shù)為200倍.可以看出：（1）對(duì)于鋁青銅材料，槳葉變形量非常小，葉稍最大變形量約為0.0018%D，幾乎可近似為剛性螺旋槳，耦合與非耦合計(jì)算在水動(dòng)力、槳葉變形及應(yīng)力分布方面基本沒(méi)有體現(xiàn)出差別；（2）對(duì)于S玻纖材料，槳葉整體的變形特征與鋁青銅保持一致，但變形量明顯高出一個(gè)數(shù)量級(jí)，最大變形量約為0.0105%D；（3）同一徑向位置導(dǎo)邊的軸向變形量要大于隨邊的軸向變形量（見(jiàn)圖8所示），這與圖5所反映的槳葉水動(dòng)力載荷特征吻合，但由于變形量量級(jí)相對(duì)于螺旋槳直徑而言均太小，槳葉的螺距分布變化非常有限，槳葉變形主要體現(xiàn)為軸向和周向的彎曲特征；（4）在槳葉等效應(yīng)力方面，最大等效應(yīng)力均出現(xiàn)在槳葉葉根的中部，兩種材料螺旋槳分布規(guī)律基本一致，從數(shù)值上看S玻纖材料略低.

圖5 槳葉壓力分布云圖（左：葉面；右：葉背）

圖6 槳葉壓力分布對(duì)比

圖7 槳葉變形輪廓對(duì)比（線框?yàn)樽冃吻?，?shí)體為變形后）

圖8 槳葉軸向變形對(duì)比

圖9 槳葉等效應(yīng)力分布對(duì)比

另外，由于變形對(duì)水動(dòng)力影響有限，流固耦合基本上迭代3次即達(dá)到式（6）所定義的收斂準(zhǔn)則.

圖10給出了不同進(jìn)速系數(shù)下敞水性能計(jì)算值與試驗(yàn)值比較，不難看出，針對(duì)所選定的槳型和工況，兩種材料螺旋槳的結(jié)構(gòu)變形均沒(méi)有體現(xiàn)出對(duì)水動(dòng)力性能的影響.圖11～12則給出了槳葉最大變形和最大等效應(yīng)力隨進(jìn)速系數(shù)的變化關(guān)系，可以看出“彈性材料”槳葉變形量更大，在各工況下其最大變形量約為“剛性材料”的6倍，近似為兩種材料彈性模量的反比關(guān)系.但即便是變形量相對(duì)較大的“彈性槳”，槳葉變形對(duì)水動(dòng)力影響也極其有限，槳葉總載荷變化很小，因而槳葉根部最大主值應(yīng)力也基本保持一致.

圖10 螺旋槳敞水性能計(jì)算值與試驗(yàn)值比較

圖11 槳葉最大變形隨進(jìn)速系數(shù)變化關(guān)系

圖12 槳葉最大等效應(yīng)力隨進(jìn)速系數(shù)變化關(guān)系

3 結(jié)束語(yǔ)

針對(duì)螺旋槳的流固耦合問(wèn)題，結(jié)合自主開(kāi)發(fā)的螺旋槳水動(dòng)力分析代碼以及結(jié)構(gòu)商業(yè)軟件的二次開(kāi)發(fā)，在通用軟件Excel中基于VBA 技術(shù)實(shí)現(xiàn)了螺旋槳的流固耦合數(shù)值模擬.通過(guò)對(duì)“剛性”和“彈性”兩種材料P4119槳的流固耦合數(shù)值計(jì)算，驗(yàn)證了所開(kāi)發(fā)螺旋槳水彈性分析平臺(tái)的有效性.且分析表明，對(duì)于所選的無(wú)側(cè)斜無(wú)縱傾平衡式螺旋槳，采用各向同性材料時(shí)，即便是采用相對(duì)更軟、變形量更大的彈性材料，槳葉的變形主要表現(xiàn)出軸向和周向的彎曲形態(tài)，槳葉螺距的變化及其有限，因而其水動(dòng)力性能與完全剛性的螺旋槳幾乎一致.對(duì)于旨在利用螺旋槳流固耦合特性進(jìn)而優(yōu)化其水動(dòng)力性能的相關(guān)研究來(lái)說(shuō)，基于本文研究結(jié)論，建議采用非各向同性材料或者采用大側(cè)斜、大縱傾等非平衡式螺旋槳設(shè)計(jì)理念.

［1］PROPULSION COMMITTEE.Final report and recommendations to the 23rd ITTC［R］.Venice：ITTC，2002.

［2］YOUNG Y L.Fluid-structure interaction analysis of flexible composite marine propellers［J］.Journal of Fluids and Structures，2008，24（6）：799-818.

［3］MULCAHY N，PRUSTY B，GARDINER C.Hydroelastic tailoring of flexible composite propellers［J］.Ships and Offshore Structures，2010，5（4）：359-370.

［4］BLASQUES J P，BERGGREEN C，ANDERSEN P.Hydro-elastic analysis and optimization of a composite marine propeller［J］.Marine Structures，2010，23（1）：22-38.

［5］YOUNG Y L.Time-dependent hydroelastic analysis of cavitating propulsors［J］.Journal of Fluids and Structures，2007，23（2）：269-95.

［6］張建國(guó)，岳 金，宋春生，等.碳纖維復(fù)合材料螺旋槳鋪層角度研究［J］.武漢理工大學(xué)學(xué)報(bào)：信息與管理工程版，2014（2）：207-210.

［7］曾志波，姚志崇，王瑋波，等.復(fù)合材料螺旋槳流固耦合分析方法研究［J］.船舶力學(xué)，2012，16（5）：477-483.

［8］孫海濤，熊 鷹.考慮變形的螺旋槳水動(dòng)力及變形特性研究［J］.哈爾濱工程大學(xué)學(xué)報(bào)，2013，34（9）：1-5.

［9］張 帥，朱 錫，侯海量.船舶螺旋槳流固耦合穩(wěn)態(tài)求解算法［J］.哈爾濱工程大學(xué)學(xué)報(bào)，2012，33（5）：615-621.

［10］洪 毅.高性能復(fù)合材料螺旋槳的結(jié)構(gòu)設(shè)計(jì)及水彈性優(yōu)化［D］.哈爾濱：哈爾濱工業(yè)大學(xué)，2011.

［11］譚廷壽.非均勻流場(chǎng)中螺旋槳性能預(yù)報(bào)和理論設(shè)計(jì)研究［D］.武漢：武漢理工大學(xué)，2003.

［12］KOYAMA K.Comparative calculations of propellers by surface panel method-workshop organized by 20 th ITTC propulsor committee［J］.Ship Research Institute，Suppement，1993（15）：1-29.

［13］JESSUP S D.An experimental investigation of viscous aspects of propeller blade flow［D］.Washington：The Catholic University of America，1989.