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

哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱150001

BEM/FEM耦合螺旋槳靜強(qiáng)度計算方法

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

哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱150001

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

螺旋槳；強(qiáng)度預(yù)報；邊界元法；有限元法；應(yīng)力分布

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 引 言

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

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

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

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

1 邊界元法理論

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

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

圖1 螺旋槳的坐標(biāo)系Fig.1 Ccoordinates of propeller

在螺旋漿表面SB上,有

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

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

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

式中，下標(biāo)TE表示螺旋槳隨邊。

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

2 有限元法計算螺旋槳的強(qiáng)度

2.1 螺旋槳有限元網(wǎng)格自動劃分方法

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

圖2 槳葉剖面的表達(dá)Fig.2 Expression of blade section

由于螺旋槳的幾何形狀比較復(fù)雜，為生成結(jié)構(gòu)單元的節(jié)點(diǎn)坐標(biāo)，首先需要對槳葉表面的坐標(biāo)進(jìn)行幾何表達(dá)。在圖2所示的柱坐標(biāo)系中，s1為槳葉剖面上點(diǎn)到導(dǎo)邊的弦向距離，c1為槳葉剖面上導(dǎo)邊至母線的距離，xr為槳葉剖面處的縱傾，θs為剖面?zhèn)刃苯?，β為螺旋槳幾何螺距角，yb和yf分別為槳葉葉背、葉面上的點(diǎn)到弦線的距離，下標(biāo)b，f分別表示螺旋槳的葉背和葉面。則在柱坐標(biāo)系O-XRθ下，螺旋槳半徑R處葉剖面上點(diǎn)的坐標(biāo)可表示為

在坐標(biāo)系O-XYZ下的相應(yīng)坐標(biāo)為

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

弦向節(jié)點(diǎn)si表示為

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

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

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

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

圖4 六面體結(jié)構(gòu)單元的生成Fig.4 Generation of the hexahedral element

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

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

對于在旋轉(zhuǎn)坐標(biāo)系中水動力載荷作用下的螺旋槳，其總體有限元結(jié)構(gòu)動力學(xué)方程可以表示為

式中：M，C和K分別為總體附加慣性力矩陣、附加阻尼力矩陣、剛度矩陣；u?，u?，u分別為節(jié)點(diǎn)的加速度、速度和位移；Fce，F(xiàn)co和Fr分別為離心力、科氏力和水動力載荷。

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

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

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

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

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

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

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

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

其中，

式中，Ni為插值函數(shù)；分別為

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

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

式（18）是在單元結(jié)構(gòu)比較規(guī)則的條件下推導(dǎo)得出。對于8節(jié)點(diǎn)的六面體單元，式（18）只適合于計算正六面體單元。由于螺旋槳結(jié)構(gòu)復(fù)雜，實際上對螺旋槳實體結(jié)構(gòu)劃分后得到的六面體是不規(guī)則的，需要引入等參元進(jìn)行坐標(biāo)轉(zhuǎn)換。通過等參元坐標(biāo)轉(zhuǎn)換后，得到如下在局部坐標(biāo)系中的單元剛度矩陣通式。

2.4 非節(jié)點(diǎn)載荷移置

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

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

螺旋槳因旋轉(zhuǎn)效應(yīng)將引起離心力，可以將其處理成體積力，其計算表達(dá)式為

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

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

從而得到單元基本方程Keue=Fe中的單元等效節(jié)點(diǎn)力列陣Fe。

3 計算流程

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

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

2）將螺旋槳表面壓力分布和結(jié)構(gòu)單元旋轉(zhuǎn)引起的離心力施加到有限元體單元中，通過整個結(jié)構(gòu)的總體力平衡方程式（11），求解螺旋槳的應(yīng)力和位移分布。

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

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

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

4 計算模型和參數(shù)設(shè)置

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

5 網(wǎng)格數(shù)和收斂性分析

根據(jù)以往在有限元方面的實踐和研究發(fā)現(xiàn)，實體結(jié)構(gòu)網(wǎng)格尺度的大小將直接影響計算結(jié)果的精度。本節(jié)借鑒相關(guān)性研究方法，預(yù)報螺旋槳在不同網(wǎng)格數(shù)情況下對計算結(jié)果的影響，并充分挖掘計算數(shù)據(jù)，以掌握計算結(jié)果與變量之間的相互關(guān)系，用于指導(dǎo)選取合適的螺旋槳網(wǎng)格數(shù)，使計算結(jié)果更準(zhǔn)確，同時也不影響計算速度。

5.1 展向和弦向網(wǎng)格數(shù)

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

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

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

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

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

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

5.2 厚度方向網(wǎng)格數(shù)目

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

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

圖9 厚度方向不同網(wǎng)格數(shù)計算的槳葉應(yīng)力分布Fig.9 Equivalent stress distributions of blade with different mesh numbers in thickness direction

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

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

6 方法驗證

對于采用邊界元法計算螺旋槳定常水動力性能的準(zhǔn)確性方面，已有多篇文獻(xiàn)［16-18］進(jìn)行過驗證，本文不再贅述，而重點(diǎn)關(guān)注本文提出的FEM/BEM耦合方法計算螺旋槳強(qiáng)度的準(zhǔn)確性，即對螺旋槳強(qiáng)度計算的結(jié)果進(jìn)行重點(diǎn)分析和驗證。

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

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

本文對螺旋槳進(jìn)行有限元靜強(qiáng)度分析時，螺旋槳實體結(jié)構(gòu)的弦向和展向與邊界元法相同，厚度方向采用平均分割，網(wǎng)格數(shù)為8。

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

圖12 槳葉應(yīng)力分布Fig.12 Equivalent stress distributions of blade

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

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

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

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

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

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

7 結(jié) 論

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

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

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

3）采用本文方法計算得到的槳葉應(yīng)力和位移分布與相關(guān)文獻(xiàn)的計算結(jié)果吻合較好，說明本文方法簡便、快速，可確保計算精度。

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

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

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

［3］楊向暉，程爾升.大側(cè)斜螺旋槳槳葉應(yīng)力測試研究［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］葉禮裕，王超，孫帥，等.基于懸臂梁法和面元法耦合的槳葉應(yīng)力分布預(yù)報［J］.武漢理工大學(xué)學(xué)報（交通科學(xué)與工程版），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］黃毅，許輝，姜治芳.大側(cè)斜螺旋槳強(qiáng)度校核探討［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］柳章存.船用螺旋槳的數(shù)值模擬及流固分析［D］.鎮(zhèn)江：江蘇科技大學(xué)，2012：1-8.

［7］陳悅，朱錫，黃政，等.水動力載荷下復(fù)合材料螺旋槳強(qiáng)度評估［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］.哈爾濱工程大學(xué)學(xué)報，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］胡壽根，王國強(qiáng)，汪庠寶.螺旋槳強(qiáng)度的厚殼元分析［J］.上海交通大學(xué)學(xué)報，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］王玉華.大側(cè)斜螺旋槳強(qiáng)度研究［J］.船舶力學(xué)，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］劉竹青，陳奕宏，姚志崇.基于面元法及有限元法耦合的螺旋槳強(qiáng)度計算［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］.哈爾濱：哈爾濱工程大學(xué)出版社，2003.

［13］肖翀，覃文潔，左正興.曲軸應(yīng)力場有限元計算單元類型和尺寸對計算精度和規(guī)模的影響［J］.柴油機(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］廖宏偉.基于迭代的六面體網(wǎng)格生成算法［D］.杭州：浙江大學(xué)，2013：20-34.LIAO H W.All-hexahedral mesh generation by iterative method［D］.Hangzhou：Zhejiang University，2013：20-34（in Chinese）.

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

［16］蘇玉民，黃勝.用面元法預(yù)報船舶螺旋槳的水動力性能［J］.哈爾濱工程大學(xué)學(xué)報，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］胡健.螺旋槳空泡性能及低噪聲螺旋槳設(shè)計研究［D］.哈爾濱：哈爾濱工程大學(xué)，2006.HU J.Research on propeller cavitation characteristics and low noise propeller design［D］.Harbin：Harbin Engineering University，2006（in Chinese）.

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

［19］顧鋮璋，鄭百林.船用螺旋槳敞水性能與槳葉應(yīng)力的數(shù)值分析［J］.力學(xué)季刊，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］譚廷壽.非均勻流場中螺旋槳性能預(yù)報和理論設(shè)計研究［D］.武漢：武漢理工大學(xué)，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">網(wǎng)絡(luò)出版時間：

時間：2017-9-26 10:53

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

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

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

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

葉禮裕，王超，李鵬，等.BEM/FEM耦合螺旋槳靜強(qiáng)度計算方法［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.