船舶槳-軸系統(tǒng)雙向流固耦合動(dòng)力學(xué)建模方法研究

鄒冬林，焦春曉，徐江海，塔 娜，饒柱石

（上海交通大學(xué)a.振動(dòng)、沖擊、噪聲研究所;b.機(jī)械系統(tǒng)與振動(dòng)國家重點(diǎn)實(shí)驗(yàn)室，上海200240）

0 引 言

船舶螺旋槳在水中工作時(shí)，由于船體艉部繞流、軸系或船體振動(dòng)、海洋湍流等因素導(dǎo)致螺旋槳伴流場(chǎng)的不均勻性，從而不可避免地在槳葉上產(chǎn)生激勵(lì)力。螺旋槳激勵(lì)力通過軸系各軸承傳遞至船體進(jìn)而引起船體振動(dòng)并輻射低頻噪聲，從而降低船舶的舒適性。眾所周知，低頻噪聲衰減慢，傳播距離遠(yuǎn)，且難以控制。而由螺旋槳激勵(lì)力激起的槳-軸系-艇體耦合振動(dòng)產(chǎn)生的低頻噪聲除具備通常低頻噪聲的特點(diǎn)外，還具有其獨(dú)特性——線譜和窄帶譜結(jié)構(gòu)特征極其明顯（主要由軸頻、葉頻及其倍頻激勵(lì)槳葉、軸系及船體等結(jié)構(gòu)共振頻率激勵(lì)導(dǎo)致）。針對(duì)這類低頻噪聲，降噪的關(guān)鍵是掌握其頻譜結(jié)構(gòu)特征（即線譜和窄帶譜結(jié)構(gòu)特征）的成因[1]，從而有針對(duì)性地對(duì)不同的頻譜特征采用不同的降噪方法。因?yàn)闃?軸系-船體系統(tǒng)與艉部不均勻伴流場(chǎng)構(gòu)成一個(gè)復(fù)雜的流-固-聲耦合系統(tǒng)，致使目前很難認(rèn)識(shí)這類低頻噪聲的產(chǎn)生機(jī)理，使得當(dāng)前的振動(dòng)噪聲治理技術(shù)難以有效應(yīng)用，減振降噪的效果不明顯。因此，掌握螺旋槳激勵(lì)導(dǎo)致的槳-軸系-船體耦合系統(tǒng)產(chǎn)生的低頻振動(dòng)噪聲機(jī)理進(jìn)而加以控制，已經(jīng)成為當(dāng)前我國船舶減振降噪急需解決的關(guān)鍵問題之一。

從螺旋槳激勵(lì)力引起的低頻噪聲的產(chǎn)生原因（前面已闡述）可以知道，這類低頻噪聲的頻譜結(jié)構(gòu)特征與螺旋槳激勵(lì)力特性以及槳-軸系統(tǒng)動(dòng)態(tài)特性均密切相關(guān)。因而對(duì)螺旋槳激勵(lì)力特性以及槳-軸系統(tǒng)振動(dòng)特性的研究有助于我們掌握這類低頻噪聲頻譜結(jié)構(gòu)特征的成因。盡管螺旋槳激勵(lì)力特性以及槳-軸系統(tǒng)耦合振動(dòng)特性是兩個(gè)不同的研究?jī)?nèi)容，但是他們之間又是一個(gè)不可分割的整體。一方面，掌握螺旋槳激勵(lì)力特性可為槳-軸系統(tǒng)振動(dòng)響應(yīng)計(jì)算提供準(zhǔn)確的輸入條件；另一方面，槳-軸系統(tǒng)耦合振動(dòng)又將對(duì)螺旋槳激勵(lì)力產(chǎn)生反饋影響。因此不應(yīng)將螺旋槳激勵(lì)力特性以及槳-軸系統(tǒng)動(dòng)態(tài)特性割裂開來單獨(dú)研究，而是需要從流體、螺旋槳及軸系三者間的耦合關(guān)系入手，從中闡明這些復(fù)雜耦合因素對(duì)槳-軸系統(tǒng)振動(dòng)特性及螺旋槳激勵(lì)力特性的影響規(guī)律，并最終揭示這類低頻噪聲中線譜和窄帶譜結(jié)構(gòu)特征的成因?？傊铌P(guān)鍵的問題是要建立一個(gè)流體-槳-軸系統(tǒng)雙向流固耦合動(dòng)力學(xué)模型。

流體-槳-軸系統(tǒng)流固耦合動(dòng)力學(xué)涉及流體動(dòng)力學(xué)、螺旋槳?jiǎng)恿W(xué)及軸系動(dòng)力學(xué)等三方面的理論知識(shí)，研究過程非常復(fù)雜，因此，國內(nèi)外直接對(duì)流體-槳-軸系統(tǒng)流固耦合振動(dòng)研究的文獻(xiàn)非常少。已有研究中，大多數(shù)把這一問題割裂為兩個(gè)問題單獨(dú)研究：流體-螺旋槳流固耦合動(dòng)力學(xué)分析及槳-軸系統(tǒng)振動(dòng)分析。對(duì)于第一個(gè)問題，通常將螺旋槳從槳轂處斷開，忽略軸系的影響。國內(nèi)外對(duì)這個(gè)問題作了大量研究：Jang 等[2]利用升力面理論結(jié)合有限元方法提出了考慮流固耦合效應(yīng)的螺旋槳設(shè)計(jì)方法；Lee等[3]利用邊界元理論（BEM）耦合有限元法（FEM）分析了復(fù)合材料螺旋槳在非均勻來流下的流固耦合特性，并將結(jié)果與基于商業(yè)軟件的CFD-FEM（Star-CCM+/Abaqus）方法進(jìn)行比較[4]；類似地，Li等[5]也利用邊界元耦合有限元方法計(jì)算了螺旋槳在非均勻來流下的流固耦合性能，重點(diǎn)討論了螺旋槳濕模態(tài)的影響因素；Ghassemi 等[6]基于商業(yè)軟件的CFD-FEM（ANSYS/CFX）方法計(jì)算了復(fù)合材料螺旋槳在均勻來流下的流固耦合力學(xué)性能；同樣，Hong 等[7]也基于商業(yè)軟件的CFD-FEM（ANSYS/CFX）方法計(jì)算了復(fù)合材料螺旋槳在均勻來流下的流固耦合力學(xué)性能；He 等[8]利用ANSYS/CFX 軟件對(duì)比研究了復(fù)合材料螺旋槳與金屬螺旋槳的流固耦合動(dòng)力學(xué)特性；曾志波等[9]基于邊界元理論與有限元法開展了復(fù)合材料螺旋槳流固耦合數(shù)值方法研究，重點(diǎn)探討了流固耦合中載荷及幾何變形量的傳遞問題，并與已有文獻(xiàn)比較，驗(yàn)證了算法的可行性；熊鷹等[10]采用瞬態(tài)計(jì)算方法研究了非均勻流場(chǎng)下螺旋槳的流固耦合性能，隨后他們又在ANSYS Workbench 平臺(tái)下，利用ACP 模塊和CFX 模塊研究了均勻流場(chǎng)下復(fù)合材料螺旋槳的流固耦合力學(xué)性能，并與實(shí)驗(yàn)結(jié)果比較，驗(yàn)證了ANSYS ACP 模塊分析復(fù)合材料螺旋槳流固耦合性能的可行性[11]；安邦等[12]基于ANSYS/CFX 軟件，利用單向流固耦合方法研究了螺旋槳在水動(dòng)載荷作用下的結(jié)構(gòu)強(qiáng)度。

對(duì)于第二個(gè)問題，通常忽略了流體與槳之間的流固耦合效應(yīng)。同樣，國內(nèi)外對(duì)這個(gè)問題也做了大量研究：Murawski[13]介紹了使用簡(jiǎn)單的集總參數(shù)模型快速評(píng)估船舶推進(jìn)軸系縱向振動(dòng)與扭轉(zhuǎn)振動(dòng)固有頻率的方法；Polic 等[14]研究了在冰雪中工作的推進(jìn)軸系扭轉(zhuǎn)振動(dòng)響應(yīng)；Zhang 等[15]利用傳遞矩陣法研究了推進(jìn)軸系的縱向振動(dòng)，探討了推力軸承剛度及位置對(duì)縱向固有頻率的影響。在上述研究中，螺旋槳均被簡(jiǎn)化為一個(gè)集中質(zhì)量單元附加在軸系末端，從而忽略了槳的彈性效應(yīng)。同時(shí)該假設(shè)也忽略了槳與軸系間的彈性耦合效應(yīng)。該處理方法在螺旋槳?jiǎng)偠群艽髸r(shí)，有很高的工程精度，這也正是國內(nèi)外學(xué)者廣泛采用的原因。目前隨著螺旋槳大型化，以及復(fù)合材料螺旋槳的出現(xiàn)，螺旋槳變得越來越“柔軟”，因此該處理方法將會(huì)引入較大的誤差。比如李小軍等[16]研究了彈性槳-軸耦合系統(tǒng)的固有頻率特性，其研究表明槳、軸間的彈性耦合效應(yīng)對(duì)系統(tǒng)固有頻率有較大影響；同樣，樓京俊等[17]將螺旋槳槳葉簡(jiǎn)化為梁-彈簧振子模型附在軸系末端，專門研究了槳葉彈性效應(yīng)對(duì)槳-軸系統(tǒng)縱振動(dòng)力學(xué)特性的影響；熊晨熙等[18]研究了螺旋槳在寬帶激勵(lì)下引起的槳-軸-船體耦合振動(dòng)問題。他們采用質(zhì)量槳-軸系和彈性槳-軸系兩種模型對(duì)槳-軸系的縱向振動(dòng)特性做了細(xì)致的對(duì)比研究。其研究結(jié)果表明，槳、軸系間的彈性耦合效應(yīng)使得槳-軸系統(tǒng)間出現(xiàn)耦合固有頻率。因此，將螺旋槳作為質(zhì)量單元考慮時(shí)，會(huì)產(chǎn)生較大誤差。

綜上所述，目前國內(nèi)外直接針對(duì)流體-彈性槳-軸系統(tǒng)流固耦合振動(dòng)特性的研究很少，大部分是將其割裂為流體-螺旋槳流固耦合動(dòng)力學(xué)分析以及螺旋槳-軸系統(tǒng)振動(dòng)特性分析兩個(gè)問題單獨(dú)研究。將這一復(fù)雜問題割裂開來處理時(shí)，將會(huì)直接忽略槳-軸系統(tǒng)振動(dòng)與螺旋槳激勵(lì)力間的耦合效應(yīng)。這是因?yàn)橐环矫嬗捎诼菪龢蜣D(zhuǎn)軸不可避免地存在偏心，使得船舶推進(jìn)軸系在繞自身中心線旋轉(zhuǎn)的同時(shí)，又發(fā)生空間渦動(dòng)（又稱為進(jìn)動(dòng)），導(dǎo)致附在其上的螺旋槳也做復(fù)雜的空間運(yùn)動(dòng)，從而使螺旋槳的進(jìn)流速度每時(shí)每刻都在變化，從而進(jìn)一步導(dǎo)致流場(chǎng)的不均勻性；另一方面，即使螺旋槳工作在均勻來流中，由于軸系振動(dòng)帶動(dòng)的螺旋槳運(yùn)動(dòng)，同樣會(huì)誘發(fā)螺旋槳和伴流場(chǎng)耦合面間的流體振蕩，使螺旋槳進(jìn)一步產(chǎn)生激勵(lì)力。

通常來說，當(dāng)軸系振動(dòng)幅值非常小且頻率很低時(shí)，由軸系振動(dòng)導(dǎo)致的螺旋槳水動(dòng)力學(xué)性能變化可以忽略，此時(shí)這種分開處理的方法滿足工程要求。而對(duì)于大型柔性船舶軸系，由于其細(xì)長(zhǎng)比通常比較小，軸系振動(dòng)幅值不能忽略；同時(shí)由于大型船舶載重增加，通常選用大尺寸螺旋槳，由此導(dǎo)致的螺旋槳微小偏心都會(huì)產(chǎn)生很大的不平衡載荷，從而引起軸系劇烈振動(dòng)。因此，把流體、螺旋槳和軸系作為一個(gè)統(tǒng)一的整體去研究其動(dòng)力學(xué)特性及其演化規(guī)律很有必要。

本文的研究目的正是在已有研究的基礎(chǔ)上，從流體、螺旋槳及軸系三者間的耦合關(guān)系入手，將兩個(gè)被割裂的問題統(tǒng)一起來，建立一個(gè)完整的流體-槳-軸系統(tǒng)雙向流固耦合動(dòng)力學(xué)模型。另外，目前已有的針對(duì)流體-螺旋槳流固耦合動(dòng)力學(xué)分析的研究文獻(xiàn)中，很多是基于商業(yè)軟件的方法。采用商業(yè)軟件（比如ANSYS+CFX）計(jì)算雙向流固耦合問題時(shí)，需要大量的計(jì)算時(shí)間，對(duì)計(jì)算機(jī)性能要求（包括CPU、內(nèi)存以及存儲(chǔ)空間等）非常高。特別是如果在已有的螺旋槳流固耦合模型中更進(jìn)一步考慮軸系振動(dòng)影響后，用商業(yè)軟件的計(jì)算時(shí)間將進(jìn)一步急劇增加。因此，本文所建立的流體-槳-軸系統(tǒng)雙向流固耦合動(dòng)力學(xué)模型需要具備計(jì)算速度快，求解精度高等優(yōu)點(diǎn)。

1 動(dòng)力學(xué)建模

一個(gè)完整的流體-槳-軸系統(tǒng)流固耦合動(dòng)力學(xué)模型涉及螺旋槳流體動(dòng)力學(xué)、螺旋槳結(jié)構(gòu)動(dòng)力學(xué)及轉(zhuǎn)子動(dòng)力學(xué)等多方面理論知識(shí)，異常復(fù)雜。而對(duì)這一復(fù)雜耦合系統(tǒng)的動(dòng)力學(xué)建模首先要解決三個(gè)子問題：（1）流體-螺旋槳雙向流固耦合動(dòng)力學(xué)建模；（2）螺旋槳-軸系動(dòng)力學(xué)建模；（3）軸系振動(dòng)對(duì)螺旋槳激勵(lì)力反饋?zhàn)饔媒?，如圖1所示。只有解決了這三個(gè)關(guān)鍵子問題，才能在此基礎(chǔ)上建立一個(gè)完整的流體-槳-軸系統(tǒng)雙向流固耦合動(dòng)力學(xué)模型。下面簡(jiǎn)單介紹這三個(gè)子問題，更詳細(xì)的描述參考文獻(xiàn)[19]。

圖1 流體-螺旋槳-軸系雙向流固耦合建模的三個(gè)子問題Fig.1 Three key problems for fluid structure interaction of propeller-shaft system

1.1 槳-軸系統(tǒng)動(dòng)力學(xué)模型

圖2 為某船舶推進(jìn)軸系示意圖，其通常由螺旋槳、轉(zhuǎn)軸、后艉軸承、前艉軸承、中間軸承及推力軸承等組成。本文利用有限元方法(FEM)對(duì)推進(jìn)軸系進(jìn)行動(dòng)力學(xué)建模。其中葉片使用20節(jié)點(diǎn)六面體單元；轉(zhuǎn)軸使用Timoshenko 梁?jiǎn)卧?；軸承簡(jiǎn)化為彈簧阻尼單元。由于槳轂與轉(zhuǎn)軸屬于過盈配合，因此建模時(shí)將其當(dāng)成轉(zhuǎn)軸的一部分，也用梁?jiǎn)卧M（大量算例表明這種處理帶來的誤差很小[19]）。

圖2 典型的船舶推進(jìn)軸系示意圖Fig.2 Diagram of a typical marine propeller-shafting-bearings system

由于梁?jiǎn)卧辛鶄€(gè)自由度，實(shí)體單元只有三個(gè)自由度，因此有限元模型中，葉片與槳轂間的連接需特殊處理，本文采用約束方程實(shí)現(xiàn)。對(duì)于平面問題，最簡(jiǎn)單的約束方程如圖3所示。

圖3 平面梁與實(shí)體間的約束方程Fig.3 Link of the beam element and solid element

建好后的槳-軸系統(tǒng)有限元模型如圖4所示。其動(dòng)力學(xué)方程可表示為

圖4 槳-軸-軸承系統(tǒng)有限元模型Fig.4 FEM model of propeller-shaft-bearings system

1.2 流體-螺旋槳流固耦合模型

流場(chǎng)模型用邊界元法（BEM，又稱為面元法，是基于勢(shì)流理論的方法）描述。盡管勢(shì)流理論忽略了流體的粘度與可壓縮性，但是大量算例均表明其能取得很好的工程精度，目前已被工程師們廣泛應(yīng)用[20]。假設(shè)螺旋槳處于不均勻流場(chǎng)中，如圖5 所示。OXYZ 為靜止慣性坐標(biāo)系，其中OX 指向船舶下游。oxyz 坐標(biāo)系固定在螺旋槳葉片上，隨葉片一起旋轉(zhuǎn)。為了避免使用動(dòng)網(wǎng)格技術(shù)，流場(chǎng)求解均在該坐標(biāo)系中進(jìn)行。

假設(shè)螺旋槳旋轉(zhuǎn)速度為ω，V0表示來流速度。螺旋槳引起的擾動(dòng)速度勢(shì)Φ(t)滿足Laplace方程：

在螺旋槳表面上滿足法向速度為零的運(yùn)動(dòng)邊界條件，即

圖5 不均勻流場(chǎng)中的螺旋槳及坐標(biāo)系示意圖Fig.5 Fixed and rotating coordinate system and an inflow wake of propellers

式中，Vin= V0+ ω × r 表示螺旋槳進(jìn)流速度，Vb為螺旋槳葉片的變形速度，其反映了葉片彈性效應(yīng)的影響。

當(dāng)Vb?Vin時(shí)，由非定常Bernoulli方程可知，葉片表面壓力分布可表示為

在物體變形不大時(shí)，物體與流場(chǎng)間的作用力可以線性分解為剛性物體與流場(chǎng)間作用力（流固耦合界面上不考慮物體的變形）和彈性物體與流場(chǎng)間作用力（流固耦合界面上不考慮物體的復(fù)雜運(yùn)動(dòng)）的疊加[19]，即

式中，φ由進(jìn)流速度Vin產(chǎn)生；φ由葉片彈性變形速度Vb產(chǎn)生。因此相應(yīng)的擾動(dòng)速度勢(shì)可以分解成兩個(gè)問題：

忽略高階小量，則式（4）可分解為

在劃分螺旋槳葉片邊界元網(wǎng)格時(shí)（采用雙曲四邊形面元，因此一個(gè)面元有四個(gè)節(jié)點(diǎn)），節(jié)點(diǎn)布置有限元節(jié)點(diǎn)完全一致，如圖6 所示，從而使得數(shù)據(jù)在流場(chǎng)和固體場(chǎng)間傳遞時(shí)更快速、更準(zhǔn)確。同時(shí)有限元節(jié)點(diǎn)與邊界元節(jié)點(diǎn)均采用余弦方式劃分[19]。這樣可以保證葉片根部、頂部、導(dǎo)邊及隨邊處的網(wǎng)格更密，從而提高計(jì)算精度。

圖6 有限元節(jié)點(diǎn)與邊界元節(jié)點(diǎn)示意圖Fig.6 Schematic diagram of finite element nodes and boundary element nodes

1.3 軸系振動(dòng)對(duì)流場(chǎng)反饋模型

由于存在外激勵(lì)（流體激勵(lì)、不平衡激勵(lì)、軸承摩擦激勵(lì)等），船舶推進(jìn)軸系同時(shí)發(fā)生彎曲振動(dòng)、扭轉(zhuǎn)振動(dòng)及縱向振動(dòng)等，從而帶動(dòng)螺旋槳一起做復(fù)雜的空間運(yùn)動(dòng)，使螺旋槳進(jìn)一步產(chǎn)生激勵(lì)力。如圖7（a）所示，軸系彎曲振動(dòng)（回旋振動(dòng)）在轉(zhuǎn)子動(dòng)力學(xué)中又稱之為渦動(dòng)或進(jìn)動(dòng)。本文使用轉(zhuǎn)子動(dòng)力學(xué)的理論來描述軸系的空間運(yùn)動(dòng)形式，引入三個(gè)坐標(biāo)系，如圖7（b）所示。OXYZ為慣性靜止坐標(biāo)系；ox1y1z1為隨體坐標(biāo)系附在葉片上，且始終與OXYZ保持平行；oxyz為隨體坐標(biāo)系，跟著葉片一起旋轉(zhuǎn)，其坐標(biāo)軸oy與葉片參考線始終重合。由此可知，ox1y1z1由OXYZ平移得到，設(shè)平移向量為R=(δx,δy,δz)T。oxyz由ox1y1z1轉(zhuǎn)動(dòng)得到，設(shè)轉(zhuǎn)動(dòng)向量為Θ=(θx,θy,θz)T。因此軸系末端振動(dòng)位移可表示為{ }X =[δx,δy,δz,θx,θy,θz]T，其振動(dòng)速度可表示為{ }X˙ =[δ˙x,δ˙y,δ˙z,θ˙x,θ˙y,θ˙z]T。由轉(zhuǎn)子動(dòng)力學(xué)理論可知：θy=?δy/?x，θz=?δz/?x。

圖7 軸系振動(dòng)與三個(gè)坐標(biāo)系示意圖Fig.7 Schematic of shaft vibration and three coordinate systems

坐標(biāo)系oxyz 與坐標(biāo)系ox1y1z1的關(guān)系可以用投影角法表示[21]，如圖8 所示。設(shè)螺旋槳的自旋軸ox 在x1y1和x1z1平面上的投影線與ox1軸的夾角為θy和θz。θy和θz稱之為投影角。坐標(biāo)系oxyz與坐標(biāo)系OXYZ的關(guān)系可表示為

圖8 投影角示意圖Fig.8 Schematic of projected angle

式中：φ = ωt + θx，ω為螺旋槳的自轉(zhuǎn)角速度；[T ]為旋轉(zhuǎn)矩陣。

由此可見，只要軸系末端位移向量{ }X =[δx,δy,δz,θx,θy,θz]T完全確定，則慣性坐標(biāo)系OXYZ 與隨葉片轉(zhuǎn)動(dòng)的坐標(biāo)系oxyz的關(guān)系亦完全確定。而軸系的位移只要外激勵(lì)確定，可以通過1.1節(jié)中的有限元方法計(jì)算得到。

由式（11）可知，在坐標(biāo)系oxyz中，螺旋槳葉片表面的速度可以表示為

將式（12）代入式（6），即可考慮軸系振動(dòng)對(duì)螺旋槳進(jìn)流速度的影響。

對(duì)于螺旋槳葉片這種升力體，從葉片隨邊泄露出的尾渦會(huì)影響葉片表面的環(huán)量，因此需要考慮尾渦的影響。對(duì)尾渦的建模主要考慮尾渦的強(qiáng)度和尾渦的形狀。對(duì)于傳統(tǒng)的根部固定不隨軸系運(yùn)動(dòng)的螺旋槳，尾渦的形狀通常假定為螺旋槳面，泄露的強(qiáng)度通常按Morino 庫塔條件或壓力庫塔條件處理[20]。而對(duì)于隨軸系做復(fù)雜空間運(yùn)動(dòng)的螺旋槳，由于其運(yùn)動(dòng)軌跡復(fù)雜，使得不同時(shí)刻隨邊的位置也不一樣，從而從隨邊泄露的第一個(gè)尾渦位置也不斷變化。因此需要對(duì)尾渦進(jìn)行合理建模，以便能考慮不同時(shí)刻螺旋槳位置及速度變化對(duì)尾渦幾何形狀的影響。

在本文中，假設(shè)尾渦的泄露是一個(gè)按時(shí)間變化的過程，如圖9 所示。在初始時(shí)刻，假設(shè)螺旋槳靜止，此時(shí)沒有尾渦泄露。在Δt時(shí)刻，螺旋槳往前移動(dòng)一個(gè)距離，此時(shí)泄出第一個(gè)尾渦，其強(qiáng)度用簡(jiǎn)單的Morino 庫塔條件表示，即

式中，Δφ(rT)表示葉片半徑為rT的隨邊處泄露的尾渦速度勢(shì)，φ+(rT)為葉背（吸力面）隨邊處的速度勢(shì)，φ-(rT)為葉面（壓力面）隨邊處的速度勢(shì)。

圖9 尾渦泄露過程Fig.9 Process of vortex-shedding

在2Δt 時(shí)刻，螺旋槳繼續(xù)往前移動(dòng)一個(gè)距離，此時(shí)泄出第二個(gè)尾渦，強(qiáng)度仍然按式（13）確定。而Δt 時(shí)刻泄露的尾渦其強(qiáng)度保存不變，在原地運(yùn)動(dòng)并發(fā)生收縮、卷曲等變形。以此時(shí)間類推，尾渦的泄露是一個(gè)連續(xù)的過程。

泄露的尾渦是不受力的，由庫塔-茹科夫斯基定理可知，泄露尾渦的速度必定與當(dāng)?shù)氐牧鲌?chǎng)速度平行。也就是說泄露的尾渦片必定按當(dāng)?shù)氐牧骶€運(yùn)動(dòng)。在OXYZ 坐標(biāo)系中，尾渦的運(yùn)動(dòng)速度為VW=?φ。因此泄露的尾渦每個(gè)時(shí)間步運(yùn)動(dòng)的距離為

假設(shè)螺旋槳共有Z 個(gè)葉片，將一個(gè)葉片及相應(yīng)輪轂劃分成Np個(gè)四邊形面元（沿弦向面元數(shù)為ns，沿展向面元數(shù)為nr，輪轂面元數(shù)為Nh，則Np= ns?nr+ Nh）。泄露尾渦面元的展向數(shù)目為nr，弦向數(shù)目由時(shí)間總步數(shù)nt決定。比如在it個(gè)時(shí)間步，弦向泄露共it個(gè)尾渦面元，其中只有緊靠葉片隨邊的第一列面元強(qiáng)度未知，其它尾渦面元強(qiáng)度均已知，如圖10所示。

因此在第it個(gè)時(shí)間步，式（6）可離散成

式中：Δφm,1為各葉片上剛泄露出的尾渦強(qiáng)度，Δφm,l為各葉片上之前時(shí)間步泄露出的尾渦強(qiáng)度；δi,j為Kronecker 函數(shù)，l和為影響系數(shù)，定義如下：

圖10 螺旋槳及泄露尾渦面元分布Fig.10 Panel arrangement of propeller and wake

一旦求得φ，即可通過式（8）與式（9）求得水動(dòng)力載荷{Fφ} 。

1.4 求解流程

至此，圖1 中的三個(gè)子問題全部解決，將其有機(jī)整合在一起，即可建立流體-槳-軸系統(tǒng)雙向流固耦合動(dòng)力學(xué)模型。完整的流體-槳-軸系統(tǒng)雙向流固耦合動(dòng)力學(xué)方程可寫為

圖11 雙向流固耦合求解流程圖Fig.11 Flowchart of bidirectional fluid structure interaction

虛線框圖中螺旋槳有限元?jiǎng)恿W(xué)模型只需在初始時(shí)刻計(jì)算一次，這是因?yàn)楸M管在流固耦合計(jì)算過程中，每個(gè)時(shí)間步的葉片幾何形狀均不一樣，理論上結(jié)構(gòu)質(zhì)量矩陣及結(jié)構(gòu)剛度矩陣均會(huì)變化，但在工程實(shí)際中，由于結(jié)構(gòu)變形很小，通常都只需計(jì)算初始幾何形狀的結(jié)構(gòu)質(zhì)量與剛度矩陣，忽略幾何形狀導(dǎo)致的幾何附加剛度矩陣。只有在變形很大時(shí)，才考慮這種幾何非線性效應(yīng)的影響。同樣地，本文將這種處理方式也應(yīng)用于附加質(zhì)量矩陣[ Ma]和附加阻尼矩陣[ Ca]。當(dāng)螺旋槳葉片變形不大時(shí)，忽略由葉片變形導(dǎo)致的流體附加質(zhì)量矩陣及附加阻尼矩陣的變化，因此本文按螺旋槳葉片的初始幾何（沒有變形）形狀計(jì)算槳葉的附加質(zhì)量矩陣及附加阻尼矩陣，且只需計(jì)算一次。

剛性葉片及軸系振動(dòng)在非均勻流場(chǎng)中產(chǎn)生的水動(dòng)力{Fφ} 的求解流程在虛線框圖外。由于本文考慮了軸系振動(dòng)對(duì)螺旋槳激勵(lì)力的影響，因此在每個(gè)時(shí)間步計(jì)算{Fφ} 時(shí)，都重新更新了軸系末端振動(dòng)速度Vs和ωs。同時(shí)，本文也考慮了葉片變形對(duì){Fφ} 的影響，并實(shí)時(shí)更新了葉片的幾何數(shù)據(jù)。

在每個(gè)時(shí)間步下得到水動(dòng)力載荷向量{Fφ} 后，通過流固耦合動(dòng)力學(xué)方程就可以求解葉片及軸系的動(dòng)態(tài)響應(yīng)，本文采用Newmark-β方法求解。得到軸系動(dòng)態(tài)響應(yīng)后，可以進(jìn)一步求解葉片的動(dòng)態(tài)應(yīng)力及軸承支反力等。

2 數(shù)值仿真驗(yàn)證

本章只介紹有限幾個(gè)數(shù)值仿真例子，驗(yàn)證本文所建模型及編制程序的正確性（更多算例詳見文獻(xiàn)[19]）。首先通過兩個(gè)算例驗(yàn)證BEM 模型及程序的正確性;計(jì)算時(shí)，假設(shè)螺旋槳從槳轂處斷開，且槳葉為剛性，目的是驗(yàn)證模型中{ }Fφ的求解精度;然后通過兩個(gè)算例驗(yàn)證FEM 模型及程序的正確性。在此計(jì)算時(shí)，不考慮水動(dòng)載荷{ }Fφ的影響。最后通過一個(gè)算例，驗(yàn)證了耦合的FEM-BEM程序的正確性。

2.1 BEM程序驗(yàn)證

本小節(jié)首先計(jì)算P4119螺旋槳在均勻流場(chǎng)中葉片壓力分布及敞水性能，并與文獻(xiàn)[22]的實(shí)驗(yàn)結(jié)果比較。計(jì)算中，雖然來流是均勻的，但是由于考慮尾渦的時(shí)變過程，因此仍然按不均勻來流的方式處理，即按時(shí)間步求解。圖12 為槳葉各半徑處的壓力系數(shù)計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果比較。從圖中可以看出，在葉片導(dǎo)邊和隨邊附近誤差較大，這是因?yàn)榱黧w在葉片邊緣處的流動(dòng)非常復(fù)雜，比如導(dǎo)邊分離等，而這些現(xiàn)象在目前的模型中均不考慮。另外在隨邊處，葉面上下壓力差并不相等，這是由于本文采用Morrino庫塔條件，而不是壓力庫塔條件導(dǎo)致的。

圖12 P4119螺旋槳壓力系數(shù)比較Fig.12 Pressure distribution comparison of P4119(J=0.833)

第二個(gè)算例是螺旋槳做縱向振動(dòng)，振動(dòng)速度為δ˙x= 0.2V0sin(ωt)。文獻(xiàn)[23]中分別用升力面(VLM)和面元法(BEM)對(duì)其推力系數(shù)與扭矩系數(shù)進(jìn)行了預(yù)報(bào)。將本文的計(jì)算結(jié)果與文獻(xiàn)[23]的預(yù)報(bào)結(jié)果做比較，以驗(yàn)證本文算法的正確性。將一個(gè)葉片弦向劃分50 個(gè)面元，展向劃分30 個(gè)面元。其余計(jì)算參數(shù)（比如計(jì)算步長(zhǎng)，旋轉(zhuǎn)角速度等）均按文獻(xiàn)[23]處理。圖13 是計(jì)算的脈動(dòng)推力系數(shù)與扭矩系數(shù)與文獻(xiàn)[23]的結(jié)果比較。從圖中可以看出，本文的計(jì)算結(jié)果與其吻合良好，求解精度在工程許可范圍以內(nèi)，從而驗(yàn)證了本文算法的正確性。

圖13 脈動(dòng)推力系數(shù)與扭矩系數(shù)比較Fig.13 Comparison of thrust and torque coefficients

2.2 FEM程序驗(yàn)證

第一個(gè)算例是計(jì)算槳-軸系統(tǒng)在葉片均布載荷下的變形及軸承支反力。均布?jí)毫χ皇┘釉诘谝粋€(gè)葉片的葉背（吸力面）上，這種不對(duì)稱的載荷分布使得各徑向軸承上存在支反力。均布?jí)毫Υ笮?00 Nm-2。將本文結(jié)果與ANSYS計(jì)算結(jié)果比較。

圖14 為5 個(gè)槳葉的變形云圖，從該圖可以看出ANSYS 計(jì)算結(jié)果與本文計(jì)算結(jié)果很吻合，最大變形約為2.8 mm。表1為所有軸承上的支反力比較，從表中可以看出相對(duì)誤差很小。因此，該算例證明了本文的槳-軸系統(tǒng)有限元模型及所編寫程序的正確性。

圖14 槳葉變形云圖比較（單位為m）Fig.14 Comparison of deformation contours for blades(the unit is‘m’)

表1 軸承支反力比較Tab.1 Comparison of bearing reaction forces

第二個(gè)算例是比較槳-軸系統(tǒng)在空氣中的固有頻率及模態(tài)振型。表2為軸系各階固有頻率比較。對(duì)于彎曲振動(dòng)，由于橫向與垂向的各階固有頻率相同，因此表格中只給出一個(gè)方向的值。從該表可以看出，本文模型計(jì)算結(jié)果與ANSYS 結(jié)果間相對(duì)誤差均很小。更進(jìn)一步地，圖15～17 為各階振型比較。從中可以看出，按模歸一化后的振型不論是最大值還是最大值發(fā)生的位置均高度一致。從該算例可以進(jìn)一步證明本文槳-軸系統(tǒng)有限元模型及所編程序的正確性。

表2 空氣中軸系固有頻率比較Tab.2 Comparison of natural frequencies of shaft system in air

圖15 第一階彎曲振型比較Fig.15 Comparison of the bend modal shape(the first)in air

圖16 第二階縱向振型比較Fig.16 Comparison of the longitudinal modal shape(the second)in air

圖17 第二階扭轉(zhuǎn)振型比較Fig.17 Comparison of the torsional modal shape(the second)in air

2.3 耦合的FEM-BEM 程序驗(yàn)證

流固耦合最典型的特征是流體會(huì)引起結(jié)構(gòu)附加質(zhì)量，從而降低結(jié)構(gòu)的固有頻率。因此計(jì)算結(jié)構(gòu)濕模態(tài)時(shí)，需要將BEM（求解流體引起的結(jié)構(gòu)附加質(zhì)量）與FEM（求結(jié)構(gòu)固有頻率）耦合在一起。本小節(jié)為了驗(yàn)證耦合的FEM-BEM 程序，計(jì)算了螺旋槳的濕模態(tài)，并與ANSYS 計(jì)算結(jié)果比較。選取P438X系列螺旋槳為研究對(duì)象，其幾何尺寸及材料屬性見文獻(xiàn)[19]。ANSYS 中，在螺旋槳周圍建立一個(gè)半徑非常大的流體域并用fluid30單元模擬流體產(chǎn)生的附加質(zhì)量效應(yīng)。由于本文計(jì)算附連水質(zhì)量時(shí)采用的勢(shì)流理論，即認(rèn)為水是不可壓縮的，因此為了能跟ANSYS 結(jié)果比較，在ANSYS 中將聲速設(shè)為無窮大，以模擬水的不可壓縮特性。

計(jì)算結(jié)果如表3 所示。從該表可以看出，本文模型計(jì)算的濕模態(tài)結(jié)果與ANSYS 濕模態(tài)結(jié)果差別很小，表明了本文模型以及編制程序的正確性。

表3 螺旋槳濕模態(tài)比較Tab.3 Comparison of natural frequencies of propeller in water

3 實(shí)驗(yàn)驗(yàn)證

圖18 重力式水洞Fig.18 Gravity tunnel

本章實(shí)驗(yàn)測(cè)量了槳-軸系統(tǒng)在空氣中及水中的固有頻率，并與本文模型計(jì)算結(jié)果比較，從而實(shí)驗(yàn)驗(yàn)證本文模型的正確性。本實(shí)驗(yàn)在重力式水洞槳-軸振動(dòng)實(shí)驗(yàn)平臺(tái)上進(jìn)行，如圖18所示。

實(shí)驗(yàn)螺旋槳使用P4381槳，該槳用ABS塑料制造，具體尺寸及材料屬性見文獻(xiàn)[19]。加工的螺旋槳以及傳感器布置如圖19所示。

采用力錘敲擊法測(cè)試系統(tǒng)固有頻率，測(cè)試結(jié)果如表4所示。從該表可以看出，總體上兩者的相對(duì)誤差很小（盡管軸系第一階縱向固有頻率相對(duì)誤差很大，但是其絕對(duì)誤差很小），從而表明了本文所建立的槳-軸系統(tǒng)動(dòng)力學(xué)模型及所編寫程序的正確性?？赡艿恼`差來源有以下幾個(gè)方面：（1）實(shí)驗(yàn)測(cè)試中的誤差；（2）理論計(jì)算時(shí)，對(duì)于軸系而言，由于推力軸承剛度值廠家并沒有給定，因此本文的剛度值依據(jù)實(shí)驗(yàn)結(jié)果推測(cè)而得，理論計(jì)算時(shí)本身也會(huì)產(chǎn)生誤差。

圖19 加工的螺旋槳及傳感器安裝圖Fig.19 Machined propeller and sensor installation

表4 槳-軸系統(tǒng)空氣中固有頻率理論計(jì)算與實(shí)驗(yàn)結(jié)果比較Tab.4 Comparison of natural frequencies between the experiment and calculation

圖20 為塑料槳各葉片在空氣中的加速度頻響函數(shù)實(shí)測(cè)及擬合結(jié)果（僅給出三個(gè)葉片的結(jié)果）。其中敲擊部位為槳葉1的中部，方向?yàn)榭v向。可以看出受傳感器安裝位置及敲擊部位等因素的影響，各葉片被激起的模態(tài)數(shù)目也略有差異。同時(shí)，不同葉片同一階固有頻率在數(shù)值上有差異，實(shí)驗(yàn)中發(fā)現(xiàn)這是受傳感器附加質(zhì)量的影響。因此表4中的理論計(jì)算結(jié)果考慮了傳感器的附加質(zhì)量影響。

圖20 槳-軸系統(tǒng)空氣中加速度頻響曲線結(jié)果（槳葉上）Fig.20 Results of acceleration FRF for propeller-shaft system in air(on the blade)

圖21為縱向敲擊槳三個(gè)不同葉片中部時(shí)，被敲擊葉片在水中的加速度頻響函數(shù)實(shí)驗(yàn)測(cè)試及擬合結(jié)果。實(shí)驗(yàn)測(cè)試時(shí)，槳葉第二階及第四階模態(tài)并沒有被測(cè)到，這是受傳感器安裝位置的限制（仿真分析時(shí)發(fā)現(xiàn)這兩階模態(tài)均是葉片邊緣的局部共振）。

圖21 槳-軸系統(tǒng)水中加速度頻響曲線結(jié)果（槳葉上）Fig.21 Results of acceleration FRF for propeller-shaft system in water (on the blade)

4 結(jié) 語

本文考慮流體、槳、軸系間的復(fù)雜耦合關(guān)系，利用有限元法（FEM）耦合邊界元法（BEM）建立了槳-軸系統(tǒng)雙向流固耦合動(dòng)力學(xué)模型。通過數(shù)值仿真分析與實(shí)驗(yàn)研究，驗(yàn)證了所建模型的正確性。本文所建立的模型可以為揭示螺旋槳激勵(lì)力產(chǎn)生、傳遞機(jī)理及槳-軸系統(tǒng)流固耦合振動(dòng)演化規(guī)律提供理論計(jì)算模型，從而能更進(jìn)一步掌握螺旋槳激勵(lì)力引起的低頻噪聲成因。

同時(shí)研究發(fā)現(xiàn)本文建立的雙向流固耦合動(dòng)力學(xué)模型相比于目前廣泛采用的商業(yè)軟件方法具有計(jì)算速度快、對(duì)計(jì)算機(jī)性能要求低等優(yōu)點(diǎn)，且計(jì)算精度滿足工程要求。這是因?yàn)椋海?）本文的模型中，流場(chǎng)采用BEM 建模，只需對(duì)物體的邊界劃分面網(wǎng)格，因而可大量節(jié)省數(shù)據(jù)存儲(chǔ)空間，減少了計(jì)算量；（2）BEM 網(wǎng)格是基于FEM 網(wǎng)格最外層“剝離”生成，因此可保證邊界元網(wǎng)格與有限元網(wǎng)格節(jié)點(diǎn)一一對(duì)應(yīng)，從而使兩種場(chǎng)之間的數(shù)據(jù)傳遞無需插值處理，使得數(shù)據(jù)傳遞更快、更準(zhǔn)確；（3）整個(gè)流場(chǎng)求解均在旋轉(zhuǎn)坐標(biāo)系中進(jìn)行，最后將結(jié)果轉(zhuǎn)換到全局坐標(biāo)系，可避免使用滑移網(wǎng)格或者動(dòng)網(wǎng)格技術(shù)。總之，本文建立的槳-軸系統(tǒng)雙向流固耦合動(dòng)力學(xué)模型是一個(gè)快速的、適合工程應(yīng)用的計(jì)算模型。