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

賀偉 李廷秋 李子如

（武漢理工大學交通學院 武漢 430063）

0 引 言

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

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

1 數學模型

1.1 螺旋槳水動力分析模型

螺旋槳水動力分析采用基于擾動速度勢的面元法進行.在固定于螺旋槳上的運動坐標系中，物面上擾動速度勢φ滿足如下積分方程：

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

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

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

槳葉和槳轂表面離散為四邊形雙曲面元，螺旋槳尾渦面采用經驗模型分為近尾流區和遠尾流區，也采用四邊形雙曲面元離散.物面和尾渦面面元上均布置等強度偶極子，其強度通過求解式（1）～（3）可惟一確定.對擾動速度勢微分可得擾動速度，結合伯努利方程可得表面壓力，粘性貢獻則通過經驗模型確定［11］.

1.2 槳葉結構分析模型

螺旋槳槳葉考慮為根部約束的線彈性結構.同樣在固定于螺旋槳上的運動坐標系中，槳葉結構運動方程為

式中：M，C和K分別為結構的質量矩陣、阻尼矩陣和剛度矩陣；和u分別為節點加速度、速度和位移矢量；Fce，Fco和Fh分別為離心力，科氏力和水動力載荷.對于瞬態問題，水動力載荷同時包含剛性槳運動和槳葉彈性變形的貢獻.對于穩態問題，外載荷與結構內力相平衡且不隨時間變化，結構分析中僅需考慮結構變形量的影響，而不用考慮結構變形速度的影響，因此結構運動方程中可省去所有與加速度和速度有關的項，槳葉彈性變形對水動力載荷的貢獻也可省去，槳葉結構運動方程退化為

式中：Fr為剛性槳運動引起的水動力載荷，通過前面介紹的水動力模型確定.在研究的初始階段，本文首先關注穩態問題.

式（5）利用有限元軟件ANSYS求解，將材料屬性、單元類型、有限元節點、水動力載荷等信息按ANSYS提供的APDL 語言格式寫入數據文件，并通過ANSYS批處理讀入該數據文件可自動進行結構分析，從而避免了復雜的交互式建模過程，為流固耦合迭代計算自動化的實現提供了可能.

1.3 螺旋槳流固耦合實現

在Excel VBA 環境中，借助于Shell函數分別調用BEM 模塊和FEM 模塊，可實現螺旋槳的流固耦合分析，調用方式如下.

其中：BEM.exe為自行開發的水動力分析可執行文件；〈InputFile〉為前文提及按APDL 語言格式生成的結構計算輸入數據文件，〈OutputFile〉為結構分析結果輸出文件.

圖1給出了螺旋槳穩態流固耦合分析流程圖.選用推力系數kt（表征水動力分析）和最大變形umax（表征結構分析）這2個參數來判斷迭代計算是否收斂，收斂準則定義為

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

圖1 流固耦合分析流程圖

2 實例分析

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

表1 螺旋槳材料參數

2.1 非耦合計算結果

首先選用鋁青銅材料螺旋槳針對設計進速系數工況J=0.833進行了非耦合計算（或稱為單向耦合），以初步驗證各個模塊的計算精度以及網格合理性.

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

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

式中：p，p0，ρ分別為槳葉表面壓力，參考壓力和流體密度.與通常基于剖面相對來流速度VR（r）進行無量綱相比，前者可更直觀的反映槳葉不同位置處壓力相對大小關系.可以看出，在遠離槳轂的外半徑0.7R和0.9R處計算值與試驗值吻合較好，而在靠近槳轂的0.3R處誤差稍大，但壓力差基本相當，且內半徑載荷的力臂相對短，認為該誤差對結構彎曲變形和扭轉變形均影響不大.

圖2 水動力分析模型

圖3 結構分析模型

圖4 不同徑向位置剖面壓力分布計算值與試驗值比較

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

槳葉的變形特征及應力分布特征在下一部分與耦合計算結果統一分析比較.

2.2 流固耦合計算結果

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

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

圖6 槳葉壓力分布對比

圖7 槳葉變形輪廓對比（線框為變形前，實體為變形后）

圖8 槳葉軸向變形對比

圖9 槳葉等效應力分布對比

另外，由于變形對水動力影響有限，流固耦合基本上迭代3次即達到式（6）所定義的收斂準則.

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

圖10 螺旋槳敞水性能計算值與試驗值比較

圖11 槳葉最大變形隨進速系數變化關系

圖12 槳葉最大等效應力隨進速系數變化關系

3 結束語

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

［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］張建國，岳 金，宋春生，等.碳纖維復合材料螺旋槳鋪層角度研究［J］.武漢理工大學學報：信息與管理工程版，2014（2）：207-210.

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

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

［9］張 帥，朱 錫，侯海量.船舶螺旋槳流固耦合穩態求解算法［J］.哈爾濱工程大學學報，2012，33（5）：615-621.

［10］洪 毅.高性能復合材料螺旋槳的結構設計及水彈性優化［D］.哈爾濱：哈爾濱工業大學，2011.

［11］譚廷壽.非均勻流場中螺旋槳性能預報和理論設計研究［D］.武漢：武漢理工大學，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.