BEM/FEM耦合螺旋槳靜強度計算方法

葉禮裕，王超，李鵬，王錫棟

哈爾濱工程大學船舶工程學院，黑龍江哈爾濱150001

BEM/FEM耦合螺旋槳靜強度計算方法

葉禮裕，王超，李鵬，王錫棟

哈爾濱工程大學船舶工程學院，黑龍江哈爾濱150001

［目的］螺旋槳強度的可靠性與船舶安全航行直接相關。為了快速且準確地計算螺旋槳強度，［方法］將邊界元（BEM）和有限元法（FEM）耦合開發螺旋槳強度預報程序。運用低階邊界元法程序對螺旋槳進行水動力性能計算，并使用普朗特—許力汀平板摩擦阻力公式進行粘性修正，然后將計算得到的槳葉表面壓力和粘性修正力作為有限元法結構計算的面力輸入。針對螺旋槳結構形狀的特殊性，發展實體螺旋槳有限元結構單元的自動劃分方法，運用有限元結構計算方程計算出螺旋槳在表面壓力和體積力作用下的應力與位移分布。以DTRC 4119模型槳為例，對提出的方法進行收斂性和網格無關性分析，并與文獻的有限元軟件計算結果進行比較，以驗證其有效性。［結果］計算結果表明，提出的方法能夠準確計算螺旋槳的應力和位移分布。［結論］該方法避免了人工建模及有限元網格劃分的過程，具有實施程序簡便、計算效率高等優點，可嵌入到螺旋槳的理論設計和優化設計過程中，形成快速計算螺旋槳強度的能力，提高螺旋槳設計的效率。

螺旋槳；強度預報；邊界元法；有限元法；應力分布

Abstract：［Objectives］The reliability of propeller stress has a great influence on the safe navigation of a ship.To predict propeller stress quickly and accurately,［Methods］a new numerical prediction model is developed by coupling the Boundary Element Method（BEM）with the Finite Element Method(FEM).The low order BEM is used to calculate the hydrodynamic load on the blades,and the Prandtl-Schlichting plate friction resistance formula is used to calculate the viscous load.Next,the calculated hydrodynamic load and viscous correction load are transmitted to the calculation of the Finite Element as surface loads.Considering the particularity of propeller geometry,a continuous contact detection algorithm is developed; an automatic method for generating the finite element mesh is developed for the propeller blade;a code based on the FEM is compiled for predicting blade stress and deformation;the DTRC 4119 propeller model is applied to validate the reliability of the method;and mesh independence is confirmed by comparing the calculated results with different sizes and types of mesh.［Results］The results show that the calculated blade stress and displacement distribution are reliable.This method avoids the process of artificial modeling and finite element mesh generation,and has the advantages of simple program implementation and high calculation efficiency.［Conclusions］The code can be embedded into the code of theoretical and optimized propeller designs,thereby helping to ensure the strength of designed propellers and improve the efficiency of propeller design.

Key words：propeller；strength prediction；Boundary Element Method（BEM）；Finite Element Method（FEM）；stress distribution

0 引 言

螺旋槳作為船舶航行過程中的動力來源，一直以來是船舶設計的重要關注點。螺旋槳結構的可靠性是保證船舶航行性能滿足要求的前提條件，故對航行安全有著重要意義。然而，隨著船舶朝大型化、高速化方向發展，大功率主機的應用導致螺旋槳槳葉的表面負荷增大，而大側斜螺旋槳的廣泛應用使螺旋槳的強度問題變得更加突出。當螺旋槳工作時，其產生的推力和旋轉阻力對槳葉具有彎曲和扭轉的作用，而且螺旋槳旋轉產生的離心力還會造成槳葉彎曲以及沿展向外拉伸。若螺旋槳強度不足，可能會使螺旋槳發生破損或斷裂，或因變形大而導致螺旋槳的水動力性能無法滿足設計要求。因此，為了提高螺旋槳設計效率和保證螺旋槳槳葉的強度，迫切需要開發一種能準確、快速預報螺旋槳槳葉強度的方法。

目前，螺旋槳槳葉強度預報可采用規范校核、數值預報和模型試驗等方法。國內外有關規范都規定了強度校核的方法和要求，但均是基于大量的使用經驗而提出，其預報結果比較保守。在模型試驗方面，Boswell［1］對大側斜螺旋槳槳葉（片）進行了靜態應力測量試驗；趙波［2］開展了大側斜螺旋槳的靜態應力試驗和動態應力試驗，并與理論計算結果進行了對比；楊向暉等［3］通過在大側斜槳模表面粘貼應變片的方法研究了不同工況下槳葉的應變和應力分布。從以上研究和試驗的結果來看，螺旋槳強度模型試驗成本較高，試驗難度大且耗時多，無法廣泛應用。

螺旋槳強度數值分析方法主要采用懸臂梁法和有限元法。懸臂梁法是一種比較方便且可行的槳葉應力預報方法，但該方法將槳葉簡化為變截面扭曲的懸臂梁，這一缺陷使得不能精準預報螺旋槳的強度［4］。對于有限元法，開展較多的是采用CFD計算與有限元分析軟件結合的方法來預報槳葉應力分布［5-6］，也有學者將邊界元法和有限元分析軟件對接預報槳葉應力分布［7-8］。這種方法雖然能準確預報螺旋槳槳葉的應力分布，但需要復雜的建模、網格劃分過程，不利于螺旋槳的設計，而且將邊界元程序與有限元軟件對接存在接口穩定性不足的問題。還有學者自編有限元程序開展螺旋槳的強度計算，如胡壽根等［9］將螺旋槳槳葉視為厚殼單元，將槳葉拆分為12個單元，編制了相應的應力分析程序；王玉華［10］專門針對大側斜槳開發了有限元強度計算程序HPROP，將升力面法計算得到的壓力值輸入到有限元程序中進行計算；劉竹青等［11］將邊界元和有限元計算程序HPROP相結合進行了螺旋槳水動力載荷作用下的強度計算。對于采用自編的有限元程序進行的強度計算，上述文獻并未詳細介紹具體的實施過程。從計算圖可知，它們對實體螺旋槳的有限元結構單元劃分存在一定的局限性，其結構單元數量較少，可能會帶來計算精度不足的問題。

本文將研究BEM/FEM耦合螺旋槳強度計算方法，提出一種準確和快速預報螺旋槳強度的方法，為設計階段的螺旋槳強度評估提供一種較好的思路。總體上，該方法是將邊界元法計算得到的螺旋槳表面壓力傳遞到螺旋槳有限元結構計算中，這里需要解決的問題是如何在螺旋槳表面網格劃分固定的情況下建立一種有限元結構單元的自動劃分方法，以實現水動力載荷在兩種方法之間的傳遞。為此，本文將詳細介紹有限元法計算螺旋槳強度的有關理論及具體的數值計算過程，采用Fortran語言編譯有限元法預報螺旋槳強度的程序，并將其與螺旋槳定常邊界元法的性能預報程序對接。最后，以某螺旋槳槳葉強度預報為例，驗證本文所提方法的有效性。

1 邊界元法理論

螺旋槳邊界元法不對物體形狀作出任何假設，其在真實物面上滿足邊界條件，能比較精確地預報螺旋槳的水動力性能，故近年來得到了廣泛應用。

圖1所示為固定于槳葉的直角坐標系O-XYZ和柱坐標系O-XRθ。圖中，R，θ分別為徑向坐標向量和角度坐標向量。假設螺旋槳在來流速度V0的情況下以角速度ω轉動，利用格林公式將用于描述不可壓、無粘、無旋的勢流問題的拉普拉斯方程轉化為物體邊界上的積分方程，使求解物體繞流問題轉化為求解任意物面上未知節點的強度問題［12］。

圖1 螺旋槳的坐標系Fig.1 Ccoordinates of propeller

在螺旋漿表面SB上,有

式中：SW為尾渦面；?為擾動速度勢；RPQ，RPQ1分別為場點P到螺旋槳表面點Q以及到尾渦面點Q1的距離；nQ，nQ1分別為螺旋槳表面上點Q和尾渦面點Q1的單位法向量；Δ?為通過尾渦面SW的速度勢跳躍，在SW上可表示為

式中：上標“+”和“-”分別表示在尾渦面上、下表面的值。

對于螺旋槳的定常問題，尾渦面速度勢跳躍Δ?在同一半徑處為常量。由法向偶極子分布與渦環分布的等價關系可知，Δ?即為尾渦強度，可在升力體的尾緣處滿足庫塔條件來確定該尾渦強度。庫塔條件有多種形式，這里采用壓力庫塔條件，它要求在升力體的尾緣處上、下表面的壓力差(Δp)TE為0，即

式中，下標TE表示螺旋槳隨邊。

結合庫塔條件，可以迭代求解得到線性方程組的數值解，即物體表面的擾動速度勢?。可以利用柳澤發展的方法由物體表面上的速度勢確定物體表面上的速度，然后通過伯努利方程計算螺旋槳表面上的壓力。

2 有限元法計算螺旋槳的強度

2.1 螺旋槳有限元網格自動劃分方法

有限元法是將復雜的連續介質求解區域離散為一組有限多個且按一定方式相互連接在一起的單元組合體。由此可知，網格劃分是有限元分析的關鍵技術之一，也是有限元數據準備過程中耗時最多、工作量最大的工作，所以在有限元技術發展過程中一直倍受關注。本節著重介紹螺旋槳有限元網格的自動劃分方法。通過程序編譯后，使用者只需輸入螺旋槳的型值表，即可實現螺旋槳有限元結構的體單元在計算機上自動完成。該方法使用簡單、性能可靠，生成的單元質量較好。

圖2 槳葉剖面的表達Fig.2 Expression of blade section

由于螺旋槳的幾何形狀比較復雜，為生成結構單元的節點坐標，首先需要對槳葉表面的坐標進行幾何表達。在圖2所示的柱坐標系中，s1為槳葉剖面上點到導邊的弦向距離，c1為槳葉剖面上導邊至母線的距離，xr為槳葉剖面處的縱傾，θs為剖面側斜角，β為螺旋槳幾何螺距角，yb和yf分別為槳葉葉背、葉面上的點到弦線的距離，下標b，f分別表示螺旋槳的葉背和葉面。則在柱坐標系O-XRθ下，螺旋槳半徑R處葉剖面上點的坐標可表示為

在坐標系O-XYZ下的相應坐標為

在邊界元法計算中，將螺旋槳表面劃分為一系列的小單元，并用雙曲面元代替每個單元，如圖3所示。這里在弦向和展向分別采用余弦分割，展向節點rj表示為

弦向節點si表示為

式中：R0，rh分別為螺旋槳半徑和槳轂半徑；bj為rj處的葉剖面弦長；Nr，Nc分別為展向和弦向的網格數；βrj，βci分別為展向節點角度和弦向節點角度，并表示如下：

考慮到螺旋槳結構的特殊性，本文對螺旋槳的實體結構沿展向、弦向以及厚度方向進行網格劃分，形成8節點的六面體單元，如圖3所示。為了更好地對接邊界元法預報程序，以便將水動力載荷傳遞到有限元結構計算中，槳葉展向和弦向的劃分方式與邊界元法相同，故螺旋槳有限元實體結構外層的節點坐標與邊界元法面網格節點坐標是重合的。這里弦向和展向采用余弦分割，主要目的是對導邊/隨邊、葉根/葉梢進行加密以反映這些區域的幾何特點。槳葉厚度方向采用平均分割。

圖3 槳葉有限元模型Fig.3 FE mesh model of propeller blade

與四面體相比，六面體單元具有更好的收斂性，達到同樣精度所需的六面體單元和節點數要遠小于四面體單元［13］。六面體單元劃分的實體不僅分析結果比四面體單元的好，而且離散后的單元數也遠小于四面體單元的單元數［14］。此外，六面體單元具有從幾何形態上更加容易辨認的優點。因此，分析人員很愿意采用六面體單元來完成三維實體的有限元分析。通過對螺旋槳實體結構的展向、弦向以及厚度方向的網格劃分，除導邊和隨邊外，其他部位被劃分為8節點的六面體，而導邊和隨邊處則被劃分為五面體單元。本文將五面體單元中導邊上的線條看成是空間四邊形退化成一條直線段的情況，在有限元計算中依然可以將其看作是8節點的六面體，如圖4所示。

圖4 六面體結構單元的生成Fig.4 Generation of the hexahedral element

2.2 總體剛度矩陣和力的平衡方程

由上節可知，本文考慮到螺旋槳結構的特殊性，采用8節點的六面體單元來劃分螺旋槳的實體結構。這里以8節點的六面體單元為例來介紹有限元法相關理論。

對于在旋轉坐標系中水動力載荷作用下的螺旋槳，其總體有限元結構動力學方程可以表示為

式中：M，C和K分別為總體附加慣性力矩陣、附加阻尼力矩陣、剛度矩陣；u?，u?，u分別為節點的加速度、速度和位移；Fce，Fco和Fr分別為離心力、科氏力和水動力載荷。

對于在均勻流中的螺旋槳，當以固定轉速旋轉時，其受到的水動力載荷為定常載荷，節點的加速度u?、速度u?以及科氏力Fco為0，則式（10）可以簡化為

由于有限元法是將實體結構劃分為單元，故將所有單元的剛度集成和疊加可得到總體剛度矩陣K，總體節點力列陣F=Fce+Fr則由所有單元的等效節點力集成和疊加而成。式（11）可離散為一個大型線性方程組，結合已知的位移邊界條件和力邊界條件，可求解得到位置的節點位移和節點力。

式中：σx，σy，σz分別為x，y，z方向上的正應力；τxy為法向面x方向且平行于y軸的剪切應力；τyz為法向面y方向且平行于z軸的剪切應力；τzx為法向面z方向且平行于x軸的剪切應力。

2.3 單元基本方程和剛度矩陣

應用虛位移原理和虛功方程，可以推導得到空間問題中單元的基本方程為［15］

式中：Ke為空間單元的剛度矩陣，上標e代表單元；ue為單元節點的位移列陣；Fe為單元等效節點力列陣。計算方法將在第2.4節介紹。

空間單元剛度矩陣Ke可表示為

式中：De為單元彈性矩陣；Be為單元應變矩陣，寫出分塊矩陣的形式為

其中，

式中，Ni為插值函數；分別為

式中，a為單元邊長；ξ，η，ζ為等參元局部坐標系。

單元彈性矩陣De為由彈性模量E與泊松比μ決定的常數矩陣，由式（18）計算得到。

式（18）是在單元結構比較規則的條件下推導得出。對于8節點的六面體單元，式（18）只適合于計算正六面體單元。由于螺旋槳結構復雜，實際上對螺旋槳實體結構劃分后得到的六面體是不規則的，需要引入等參元進行坐標轉換。通過等參元坐標轉換后，得到如下在局部坐標系中的單元剛度矩陣通式。

2.4 非節點載荷移置

結構離散后，各單元通過節點連接，結構位移近似地由所有節點的位移表示，所有外載荷都要等效地移置到單元節點上，以用于單元特性分析。根據彈性力學的虛位移原理，可將外載荷移置到單元節點上。

式中，N為形函數矩陣。螺旋槳不受到集中力作用，這個力的大小設置為0。

螺旋槳因旋轉效應將引起離心力，可以將其處理成體積力，其計算表達式為

式中：ρ為螺旋槳材料密度；ω為螺旋槳的旋轉角速度；X為節點坐標向量。如果單元e的某一界面上分布有面力，將微分面dA上的力（A為單元面積）作為集中力，可得到移置后的等效節點力列陣為

對于螺旋槳而言，主要是螺旋槳旋轉過程中槳葉受到水動力載荷和粘性阻力的作用。基于本文建立的螺旋槳有限元網格自動劃分方法，可使螺旋槳有限元實體結構外層面元與邊界元法面元重合。因此，邊界元法計算得到的螺旋槳面元中心處的壓力分布，可作為面力施加到有限元結構中，并通過式（26）等效地移置到單元節點上。

從而得到單元基本方程Keue=Fe中的單元等效節點力列陣Fe。

3 計算流程

本文采用雙向流—固耦合方法，即邊界元法預報的計算結果和有限元法預報的結構計算結果是相互傳遞的，兩種方法在迭代時由于邊界條件的改變而開始了新的計算。信息在流體和結構模塊中往復傳遞，直至獲得滿足收斂條件的解為止。圖5所示為本文計算流程圖，具體的迭代求解過程如下：

1）基于速度勢的低階邊界元法，對每個槳葉布置的雙曲面元求解式（1）獲得擾動速度勢?，利用伯努利方程求解得到面元控制點的壓力分布和水動力性能。

2）將螺旋槳表面壓力分布和結構單元旋轉引起的離心力施加到有限元體單元中，通過整個結構的總體力平衡方程式（11），求解螺旋槳的應力和位移分布。

3）將有限元法計算得到的螺旋槳表面節點位移與邊界元法點坐標相加，求解式（1）獲得擾動速度勢?，利用伯努利方程求取面元控制點的壓力分布和水動力性能。

4）重復步驟2）和步驟3），直到滿足最大位移收斂條件為止。

圖5 螺旋槳的流—固耦合求解過程Fig.5 Calculation process of fluid-solid interaction for propeller

4 計算模型和參數設置

本文以DTRC 4119模型槳為研究對象，考察網格數以及網格劃分對計算結果的影響，并評估本文所提方法的可行性。其中，螺旋槳直徑為0.305 m，轂徑比為0.2，無縱傾，無側斜，剖面翼型為NACA-66mod。螺旋槳材料選擇為密度為7 600 kg/m3的鎳—鋁青銅，材料的彈性模量E= 113 GPa，泊松比μ=0.34。鑒于螺旋槳槳葉通常與槳轂采取剛性連接，為便于計算，本文模型將槳葉葉根處節點進行了六自由度剛性約束。計算工況設定為：設計進速系數為0.833，轉速為600 r/min。

5 網格數和收斂性分析

根據以往在有限元方面的實踐和研究發現，實體結構網格尺度的大小將直接影響計算結果的精度。本節借鑒相關性研究方法，預報螺旋槳在不同網格數情況下對計算結果的影響，并充分挖掘計算數據，以掌握計算結果與變量之間的相互關系，用于指導選取合適的螺旋槳網格數，使計算結果更準確，同時也不影響計算速度。

5.1 展向和弦向網格數

基于上述方法，采用10種網格劃分方式，即弦向和展向網格數為10×10，12×12，14×14，16× 16，18×18，20×20，22×22，24×24，26×26和28×28，并將厚度方向的網格數固定為6，然后對計算結果進行分析。

圖6和圖7所示分別為將槳葉應力預報結果導入到Tecplot軟件后得到的槳葉應力與位移分布，即在計算工況下，弦向和展向網格數分別為12×12，16×16和20×20時槳葉葉面的應力和位移分布。預報時，為便于比較，設置了相同的等高線取值范圍。由圖6可知，隨著網格數的增加，槳葉應力不斷增大，槳葉應力分布更加均勻；網格數過少時，容易導致槳葉應力峰值集中于某一點上。由圖7可知，不同網格數得到的槳葉位移分布趨勢基本一致，但是位移量的區別較大。

圖6 弦向和展向不同網格數時計算得到的槳葉應力分布Fig.6 Equivalent stress distributions of blade with different chord-wise and span-wise mesh numbers

圖7 弦向和展向不同網格數時計算得到的槳葉位移分布Fig.7 Displacement distributions of blade with different chord-wise and span-wise mesh numbers

為了更好地分析弦向和展向網格數對計算結果收斂性的影響，圖8給出了在計算工況下弦向和展向不同網格數對應的最大等效應力與最大位移量。由圖8可知，隨著網格數的增加，最大等效應力和最大位移量均不斷增大，但增大的幅度有所減小；當網格數達到24時，雖然最大等效應力和最大位移量還有增長的趨勢，但增長的幅度已經很小。因此，弦向和展向網格數為24×24時基本可以認為計算結果收斂了。

圖8 網格收斂性分析Fig.8 Convergence process with different chord-wise and span-wise mesh numbers

5.2 厚度方向網格數目

基于上述方法，對厚度方向網格數為2，3，4，5，6，7，8的計算結果進行分析，將弦向和展向的網格數固定為24×24。厚度方向的網格數只對有限元網格單元的劃分有影響，不會影響邊界元法的網格劃分，故在未開始流—固耦合迭代前，不同厚度方向的網格數水動力性能的計算結果相同。

圖9和圖10所示為在計算工況下厚度方向網格數為2，4和6的槳葉葉面應力與位移分布。為便于比較，設置了相同的等高線取值范圍。由圖9可知，厚度方向網格數對槳葉應力計算結果影響很大，隨著網格數增加，紅色區域范圍越來越大，說明槳葉應力呈整體增大的趨勢。由圖10可知，不同網格數計算得到的槳葉位移分布趨勢基本一致，但隨著網格數的增加，紅色區域范圍變大，說明槳葉的位移也呈整體增大趨勢。由此可見，厚度方向網格數對螺旋槳強度的計算結果影響很大。

圖9 厚度方向不同網格數計算的槳葉應力分布Fig.9 Equivalent stress distributions of blade with different mesh numbers in thickness direction

圖10 厚度方向不同網格數計算的槳葉位移分布Fig.10 Displacement distributions of blade with different mesh numbers in thickness direction

為了更好地分析弦向和展向網格數對計算結果收斂性的影響，圖11給出了在計算工況下不同厚度方向網格數計算得到的最大等效應力和位移量。由圖11可知，隨著網格數的增加，最大等效應力和位移量均不斷增大，但增大的幅度迅速減小。因此當網格數增大到6以上時，基本可以認為計算結果收斂了。

6 方法驗證

對于采用邊界元法計算螺旋槳定常水動力性能的準確性方面，已有多篇文獻［16-18］進行過驗證，本文不再贅述，而重點關注本文提出的FEM/BEM耦合方法計算螺旋槳強度的準確性，即對螺旋槳強度計算的結果進行重點分析和驗證。

圖11 厚度方向網格收斂性分析Fig.11 Convergence process with different mesh numbers in thickness direction

為了驗證本文方法預報螺旋槳強度的準確性，將本文應力預報結果與文獻［19］進行了比較。根據上述網格數和收斂性分析，在計算螺旋槳的定常邊界元法水動力性能時，螺旋槳表面的弦向和展向均采用半余弦分割，弦向和展向網格數為26×26。在計算工況下，本文采用邊界元法計算程序得到的推力系數和扭矩系數分別為0.143 2和0.026 5，文獻［19］計算得到的分別為0.135 2和0.028 1，模型試驗測量得到的分別為0.141 2和0.027 8［20］。由此可見，與模型試驗測量得到的推力系數和轉矩系數相比，本文方法計算的結果與試驗測量結果較為接近；與文獻［19］的計算結果相比，本文計算得到的推力系數比文獻［19］的大，扭矩系數比文獻［19］的小。

本文對螺旋槳進行有限元靜強度分析時，螺旋槳實體結構的弦向和展向與邊界元法相同，厚度方向采用平均分割，網格數為8。

圖12所示為本文方法預報得到的DTRC 4119模型槳葉背和葉面的應力分布。通過與文獻［19］計算得到的分布結果比較，發現本文方法計算的槳葉應力分布趨勢與文獻［19］的結果基本一致。螺旋槳強度校核主要關注最大等效應力是否超過許用應力。本文計算的槳葉最大等效應力為1.31MPa，而文獻［19］預報的最大等效應力為1.18 MPa，可見本文方法預報的結果相較于文獻［19］的偏大。雖然兩種方法計算得到的最大等效應力有一定的偏差，但從量級來看還是比較合理的。導致計算結果出現偏差的原因很多，其中包括：邊界元法和CFD這2種方法預報的螺旋槳水動力載荷壓力分布不可能相同；兩種方法的有限元網格劃分類型與數量不同；兩種方法載荷施加和傳遞的方式不盡相同。另外應注意的是，本文方法計算得到的槳葉應力集中在弦向槳葉葉根的中部，與文獻［19］得到的結論一致。

圖12 槳葉應力分布Fig.12 Equivalent stress distributions of blade

圖13所示為采用本文方法預報得到的DTRC 4119模型槳葉背和葉面的位移分布。由圖13可知，葉面和葉背的位移量相同。本文方法計算得到的應力分布趨勢與文獻［19］計算的結果一致。槳葉位移主要影響螺旋槳的水動力性能，由于本文螺旋槳模型為剛性槳，位移量小，故不會引起水動力變化。本文計算得到的槳葉最大位移量為7.92×10-6mm，而文獻［19］預報的最大位移量為7.0×10-6mm，可見本文預報得到的值有些偏大。

考慮到目前國內對螺旋槳的強度校核大多采用懸臂梁法，這里將本文方法與懸臂梁法的校核結果進行了對比，以驗證本文方法的可信性。文獻［4］提出了懸臂梁法和邊界元法耦合的槳葉應力分布預報方法，由于懸臂梁的局限性，無法預報槳葉的位移。圖14所示為文獻［4］采取與本文相同工況時預報的DTRC 4119模型槳槳葉的應力分布。由圖12和圖14的對比可知，本文方法計算的槳葉應力分布趨勢與懸臂梁法計算的結果類似，槳葉最大應力均集中在葉根的中部，但槳葉應力最大值卻存在一定的差異。其中，懸臂梁方法預報的槳葉最大拉應力為1.03 MPa，相較于本文方法計算的結果，懸臂梁法預報的結果偏小。造成上述偏差的原因為：懸臂梁法將槳葉過度簡化；兩種方法采用的強度理論有所區別，懸臂梁法采用的是最大拉應力理論（第1強度理論），而本文方法采用的是形狀改變比能理論（第4強度理論）。

圖13 槳葉位移分布Fig.13 Displacement distributions of blade

圖14 懸臂梁法預報的DTRC 4119模型槳槳葉應力分布Fig.14 Stress distributions of DTRC 4119 blade predicted by the cantilever beam method

為了驗證本文方法對大側斜螺旋槳的適用性，預報了側斜角度為36°的槳葉應力分布和位移分布，結果如圖15所示。由圖可知，槳葉最大應力發生在葉根且靠近隨邊的部位，與文獻［8］得到的結論一致。

圖15 本文方法預報的DTRC 4382模型槳槳葉應力和位移分布Fig.15 Stress and displacement distributions of DTRC 4382 predicted by current method

7 結 論

本文介紹了有限元法求解螺旋槳靜強度的相關理論，分別提出和研究了一種螺旋槳有限元結構單元劃分及BEM/FEM耦合預報螺旋槳靜強度的方法，探討了在槳葉弦向、展向及厚度方向采用不同網格數對計算結果和收斂性的影響，并將本文方法與相關文獻計算的結果進行了比較，驗證了本文方法的可行性。針對本文計算對象，分析得到如下結論：

1）對槳葉弦向和展向采用不同網格數對計算結果的影響分析表明：隨著網格數的增加，槳葉應力分布更均勻，最大等效應力和最大位移量均呈不斷增大的趨勢，但增大的幅度迅速減小，槳葉弦向和展向的網格數需在26×26以上才可收斂。

2）對槳葉厚度方向采用不同網格數對計算結果的影響分析表明：隨著網格數的增加，槳葉應力和位移呈整體增大的趨勢，但增大的幅度迅速減小，在厚度方向的網格數需在6以上才可收斂。

3）采用本文方法計算得到的槳葉應力和位移分布與相關文獻的計算結果吻合較好，說明本文方法簡便、快速，可確保計算精度。

應用本文方法便于對螺旋槳的靜強度進行分析，可快速計算螺旋槳槳葉的應力和位移分布，后續將在其他類型的螺旋槳中予以應用，以進一步對所提方法進行驗證，并將該程序嵌入到螺旋槳理論設計和優化設計過程中，提高螺旋槳在設計階段的強度評估和計算效率。

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

［2］趙波.大側斜螺旋槳強度研究［D］.武漢：華中科技大學，2003.

［3］楊向暉，程爾升.大側斜螺旋槳槳葉應力測試研究［J］.船海工程，2004（1）：6-9.YANG X H，CHENG E S.Study on stress measurement of the high skewed propeller blades［J］.Ship& Ocean Engineering，2004（1）：6-9（in Chinese）.

［4］葉禮裕，王超，孫帥，等.基于懸臂梁法和面元法耦合的槳葉應力分布預報［J］.武漢理工大學學報（交通科學與工程版），2015（5）：968-973.YE L Y，WANG C，SUN S，et al.Prediction of blade stress distribution based on cantilever beam method and panel method［J］.Journal of Wuhan University of Technology（Transportation Science&Engineering），2015（5）：968-973（in Chinese）.

［5］黃毅，許輝，姜治芳.大側斜螺旋槳強度校核探討［J］.中國艦船研究，2010，5（5）：44-48.HUANG Y，XU H，JIANG Z F.Strength analysis of highly-skewed propeller［J］.Chinese Journal of Ship Research，2010，5（5）：44-48（in Chinese）.

［6］柳章存.船用螺旋槳的數值模擬及流固分析［D］.鎮江：江蘇科技大學，2012：1-8.

［7］陳悅，朱錫，黃政，等.水動力載荷下復合材料螺旋槳強度評估［J］.中國艦船研究，2015，10（1）：19-26.CHEN Y，ZHU X，HUANG Z，et al.Strength evaluation of the composite propeller under hydrodynamic fluid load［J］.Chinese Journal of Ship Research，2015，10（1）：19-26（in Chinese）.

［8］孫海濤，熊鷹.考慮變形的螺旋槳水動力及變形特性研究［J］.哈爾濱工程大學學報，2013，34（9）：1108-1112，1118.SUN H T，XIONG Y.Study on hydrodynamic and deformation performance of propellers considering the blade deformation［J］.Journal of Harbin Engineering University，2013，34（9）：1108-1112，1118（in Chinese）.

［9］胡壽根，王國強，汪庠寶.螺旋槳強度的厚殼元分析［J］.上海交通大學學報，1988，22（2）：33-42.HU S G，WANG G Q，WANG X B.Thick shell element method of stress analysis of propeller blades［J］.Journal of Shanghai Jiaotong University，1988，22（2）：33-42（in Chinese）.

［10］王玉華.大側斜螺旋槳強度研究［J］.船舶力學，1998，2（2）：44-51.WANG Y H.Comparative studies on highly-skewed propeller strength［J］.Journal of Ship Mechanics，1998，2（2）：44-51（in Chinese）.

［11］劉竹青，陳奕宏，姚志崇.基于面元法及有限元法耦合的螺旋槳強度計算［J］.中國造船，2012，53（增刊1）：25-30.LIU Z Q，CHEN Y H，YAO Z C.Static strength analysis of civil ship propellers［J］.Shipbuilding of China，2012，53（Supp 1）：25-30（in Chinese）.

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

［13］肖翀，覃文潔，左正興.曲軸應力場有限元計算單元類型和尺寸對計算精度和規模的影響［J］.柴油機，2004（S1）：176-178.XIAO C，QIN W J，ZUO Z X.Influence of elementform and element-size on FE calculation of crankshaft stress field［J］.Diesel Engine，2004（Supp 1）：176-178（in Chinese）.

［14］廖宏偉.基于迭代的六面體網格生成算法［D］.杭州：浙江大學，2013：20-34.LIAO H W.All-hexahedral mesh generation by iterative method［D］.Hangzhou：Zhejiang University，2013：20-34（in Chinese）.

［15］趙均海，汪夢甫.彈性力學及有限元［M］.武漢：武漢理工大學出版社，2003：18-20.

［16］蘇玉民，黃勝.用面元法預報船舶螺旋槳的水動力性能［J］.哈爾濱工程大學學報，2001，22（2）：1-5.SU Y M，HUANG S.Prediction of hydrodynamic performance of marine propellers by surface panel method［J］.JournalofHarbin Engineering University，2001，22（2）：1-5（in Chinese）.

［17］胡健.螺旋槳空泡性能及低噪聲螺旋槳設計研究［D］.哈爾濱：哈爾濱工程大學，2006.HU J.Research on propeller cavitation characteristics and low noise propeller design［D］.Harbin：Harbin Engineering University，2006（in Chinese）.

［18］王超.螺旋槳水動力性能、空泡及噪聲性能的數值預報研究［D］.哈爾濱：哈爾濱工程大學，2010.WANG C.The research on performance of propeller's hydrodynamics，cavitation and noise［D］.Harbin：Harbin Engineering University，2010（in Chinese）.

［19］顧鋮璋，鄭百林.船用螺旋槳敞水性能與槳葉應力的數值分析［J］.力學季刊，2011，32（3）：440-443.GU C Z，ZHENG B L.Numerical analysis of propeller open-water performance and stress distribution in the blade of ship propeller［J］.Chinese Quarterly of Mechanics，2011，32（3）：440-443（in Chinese）.

［20］譚廷壽.非均勻流場中螺旋槳性能預報和理論設計研究［D］.武漢：武漢理工大學，2003.TAN T S.Performance prediction and theoretical design research on propeller in non-uniform flow［D］.Wuhan：Wuhan University of Technology，2003（in Chinese）.

Calculation of marine propeller static strength based on coupled BEM/FEM

YE Liyu，WANG Chao，LI Peng，WANG Xidong
College of Shipbuilding Engineering，Harbin Engineering University，Harbin 150001，China

U664.33

A

10.3969/j.issn.1673-3185.2017.05.008

2017-03-28< class="emphasis_bold">網絡出版時間：

時間：2017-9-26 10:53

國家部委基金資助項目；國家自然科學基金資助項目（51679052）

葉禮裕，男，1989年生，博士生。研究方向：冰區船舶推進器設計及性能評估。E-mail：yeliyuxrxy@126.com

王超（通信作者），男，1981年生，博士，副教授。研究方向：艦船推進及減振降噪技術，冰區推進技術。E-mail：wangchao0104@hrbeu.edu.cn

http://kns.cnki.net/kcms/detail/42.1755.TJ.20170926.1053.022.html期刊網址：www.ship-research.com

葉禮裕，王超，李鵬，等.BEM/FEM耦合螺旋槳靜強度計算方法［J］.中國艦船研究，2017，12（5）：64-74，83.

YE L Y，WANG C，LI P，et al.Calculation of marine propeller static strength based on coupled BEM/FEM［J］.Chinese Journal of Ship Research，2017，12（5）：64-74，83.