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

（上海交通大學a.振動、沖擊、噪聲研究所;b.機械系統(tǒng)與振動國家重點實驗室，上海200240）

0 引 言

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

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

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

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

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

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

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

1 動力學建模

一個完整的流體-槳-軸系統(tǒng)流固耦合動力學模型涉及螺旋槳流體動力學、螺旋槳結構動力學及轉(zhuǎn)子動力學等多方面理論知識，異常復雜。而對這一復雜耦合系統(tǒng)的動力學建模首先要解決三個子問題：（1）流體-螺旋槳雙向流固耦合動力學建模；（2）螺旋槳-軸系動力學建模；（3）軸系振動對螺旋槳激勵力反饋作用建模，如圖1所示。只有解決了這三個關鍵子問題，才能在此基礎上建立一個完整的流體-槳-軸系統(tǒng)雙向流固耦合動力學模型。下面簡單介紹這三個子問題，更詳細的描述參考文獻[19]。

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

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

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

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

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

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

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

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

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

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

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

在螺旋槳表面上滿足法向速度為零的運動邊界條件，即

圖5 不均勻流場中的螺旋槳及坐標系示意圖Fig.5 Fixed and rotating coordinate system and an inflow wake of propellers

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

當Vb?Vin時，由非定常Bernoulli方程可知，葉片表面壓力分布可表示為

在物體變形不大時，物體與流場間的作用力可以線性分解為剛性物體與流場間作用力（流固耦合界面上不考慮物體的變形）和彈性物體與流場間作用力（流固耦合界面上不考慮物體的復雜運動）的疊加[19]，即

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

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

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

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

1.3 軸系振動對流場反饋模型

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

圖7 軸系振動與三個坐標系示意圖Fig.7 Schematic of shaft vibration and three coordinate systems

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

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

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

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

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

將式（12）代入式（6），即可考慮軸系振動對螺旋槳進流速度的影響。

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

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

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

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

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

泄露的尾渦是不受力的，由庫塔-茹科夫斯基定理可知，泄露尾渦的速度必定與當?shù)氐牧鲌鏊俣绕叫?。也就是說泄露的尾渦片必定按當?shù)氐牧骶€運動。在OXYZ 坐標系中，尾渦的運動速度為VW=?φ。因此泄露的尾渦每個時間步運動的距離為

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

因此在第it個時間步，式（6）可離散成

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

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

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

1.4 求解流程

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

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

虛線框圖中螺旋槳有限元動力學模型只需在初始時刻計算一次，這是因為盡管在流固耦合計算過程中，每個時間步的葉片幾何形狀均不一樣，理論上結構質(zhì)量矩陣及結構剛度矩陣均會變化，但在工程實際中，由于結構變形很小，通常都只需計算初始幾何形狀的結構質(zhì)量與剛度矩陣，忽略幾何形狀導致的幾何附加剛度矩陣。只有在變形很大時，才考慮這種幾何非線性效應的影響。同樣地，本文將這種處理方式也應用于附加質(zhì)量矩陣[ Ma]和附加阻尼矩陣[ Ca]。當螺旋槳葉片變形不大時，忽略由葉片變形導致的流體附加質(zhì)量矩陣及附加阻尼矩陣的變化，因此本文按螺旋槳葉片的初始幾何（沒有變形）形狀計算槳葉的附加質(zhì)量矩陣及附加阻尼矩陣，且只需計算一次。

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

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

2 數(shù)值仿真驗證

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

2.1 BEM程序驗證

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

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

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

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

2.2 FEM程序驗證

第一個算例是計算槳-軸系統(tǒng)在葉片均布載荷下的變形及軸承支反力。均布壓力只施加在第一個葉片的葉背（吸力面）上，這種不對稱的載荷分布使得各徑向軸承上存在支反力。均布壓力大小為500 Nm-2。將本文結果與ANSYS計算結果比較。

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

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

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

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

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

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

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

3 實驗驗證

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

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

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

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

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

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

圖20 為塑料槳各葉片在空氣中的加速度頻響函數(shù)實測及擬合結果（僅給出三個葉片的結果）。其中敲擊部位為槳葉1的中部，方向為縱向?？梢钥闯鍪軅鞲衅靼惭b位置及敲擊部位等因素的影響，各葉片被激起的模態(tài)數(shù)目也略有差異。同時，不同葉片同一階固有頻率在數(shù)值上有差異，實驗中發(fā)現(xiàn)這是受傳感器附加質(zhì)量的影響。因此表4中的理論計算結果考慮了傳感器的附加質(zhì)量影響。

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

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

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

4 結 語

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

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