鄒冬林，張建波，田佳彬，塔 娜，饒柱石，4

（1.上海交通大學(xué) 振動(dòng)、沖擊、噪聲研究所，上海200240；2.上海交通大學(xué) 機(jī)械系統(tǒng)與振動(dòng)國家重點(diǎn)實(shí)驗(yàn)室，上海200240；3.武漢第二船舶設(shè)計(jì)研究所，武漢430205；4.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海200240）

0 引 言

由于海洋來流的不均勻，船舶螺旋槳工作中產(chǎn)生靜推力的同時(shí)也不可避免地產(chǎn)生脈動(dòng)力。該脈動(dòng)力傳遞到軸上引起軸系振動(dòng)，進(jìn)而引起船體振動(dòng)并輻射噪聲。過大的振動(dòng)會導(dǎo)致軸系、船體等結(jié)構(gòu)疲勞損壞，產(chǎn)生安全隱患。同時(shí)由振動(dòng)而產(chǎn)生的輻射噪聲也會影響乘員舒適性。因此在船舶推進(jìn)軸系設(shè)計(jì)階段，需要考慮如何減少螺旋槳脈動(dòng)力的產(chǎn)生以及優(yōu)化螺旋槳脈動(dòng)力到軸系的傳遞途徑等一系列問題。為解決這一問題，需要全面分析流體-螺旋槳-軸系間的耦合動(dòng)力學(xué)特性。流體-螺旋槳-軸系耦合系統(tǒng)動(dòng)力學(xué)分析本質(zhì)上是一個(gè)流固耦合問題，且非常復(fù)雜。這是因?yàn)橐环矫娌痪鶆蛄黧w引起的螺旋槳脈動(dòng)力會導(dǎo)致軸系振動(dòng)，而軸系振動(dòng)又會帶動(dòng)螺旋槳振動(dòng)，從而改變螺旋槳的脈動(dòng)力。二者之間相互耦合，形成反饋機(jī)制。另一方面，即使螺旋槳工作在均勻來流中，由于螺旋槳或轉(zhuǎn)軸不可避免地存在偏心，使得船舶推進(jìn)軸系在繞自身中心線旋轉(zhuǎn)的同時(shí)，又發(fā)生空間渦動(dòng)（又稱為進(jìn)動(dòng)），導(dǎo)致附在其上的螺旋槳也做復(fù)雜的空間旋轉(zhuǎn)運(yùn)動(dòng)，從而使螺旋槳的進(jìn)流速度每時(shí)每刻都在變化，進(jìn)一步導(dǎo)致脈動(dòng)力的產(chǎn)生?？傊?，螺旋槳脈動(dòng)力和軸系振動(dòng)是相互影響的，并且可能存在增強(qiáng)效應(yīng)。根據(jù)相關(guān)文獻(xiàn)資料，由于船體和軸系振動(dòng)誘發(fā)了螺旋槳和伴流場耦合面的流體振蕩，可使螺旋槳激勵(lì)力產(chǎn)生的低頻聲輻射增加10～15 dB，因此非常有必要考慮軸系振動(dòng)對螺旋槳水動(dòng)性能的影響。目前針對流體-螺旋槳-軸系的雙向流固耦合這一復(fù)雜問題的研究文獻(xiàn)很少。國內(nèi)外已有的文獻(xiàn)均是把這一問題割裂為兩個(gè)問題：流體-螺旋槳流固耦合動(dòng)力學(xué)分析及流體載荷激勵(lì)下的軸系振動(dòng)。通常是先研究流體-螺旋槳流固耦合動(dòng)力學(xué)特性，求解出螺旋槳對根部的激勵(lì)力。研究中認(rèn)為螺旋槳根部固定，忽略軸系振動(dòng)的影響。然后把激勵(lì)力作為集中載荷施加在軸系末端上（與螺旋槳連接的地方），進(jìn)而研究軸系的振動(dòng)響應(yīng)等等，此時(shí)螺旋槳通常作為一個(gè)集中質(zhì)量。該研究思路忽略了軸系振動(dòng)對流場的反饋?zhàn)饔?。?dāng)軸系振動(dòng)幅值很小時(shí)，由此導(dǎo)致的螺旋槳水動(dòng)力學(xué)性能變化可以忽略，這種分開處理的方法是滿足工程要求的。而對于大型柔性船舶軸系，由于其細(xì)長比通常很小，軸系振動(dòng)幅值不能忽略；同時(shí)由于大型船舶載重增加，通常選用大尺寸螺旋槳，由此導(dǎo)致的螺旋槳微小偏心都會產(chǎn)生很大的不平衡載荷，從而引起軸系劇烈振動(dòng)。因此把流體、螺旋槳和軸系作為一個(gè)統(tǒng)一的整體去研究其動(dòng)力學(xué)特性及其演化規(guī)律是必要的。

對這一復(fù)雜耦合系統(tǒng)的動(dòng)力學(xué)分析首先要研究以下三個(gè)方面的子問題，如圖1所示。第一個(gè)子問題是流體-螺旋槳雙向流固耦合問題的研究；第二個(gè)子問題是螺旋槳脈動(dòng)力引起軸系振動(dòng)這一單向流固耦合問題的研究；第三個(gè)子問題是軸系振動(dòng)導(dǎo)致的螺旋槳脈動(dòng)力變化的研究。只有將這三個(gè)關(guān)鍵子問題摸透了，才能將所有子問題串聯(lián)在一起，進(jìn)而建立一個(gè)復(fù)雜的流體-螺旋槳-軸系雙向流固耦合動(dòng)力學(xué)模型。

第一個(gè)子問題主要研究螺旋槳葉片變形、槳葉動(dòng)態(tài)應(yīng)力及水動(dòng)力預(yù)報(bào)等等。第二個(gè)子問題主要研究軸系振動(dòng)響應(yīng)、固有頻率及軸承力等等。目前對第一個(gè)和第二個(gè)子問題均有大量研究文獻(xiàn)[1-4]，而對第三個(gè)子問題的研究幾乎沒有。第三個(gè)子問題是第二個(gè)子問題的逆問題，其核心點(diǎn)是如何預(yù)報(bào)隨軸系做復(fù)雜空間運(yùn)動(dòng)的螺旋槳的脈動(dòng)力。而本文的研究重點(diǎn)集中在第三個(gè)子問題上，目的在于探索一種復(fù)雜空間運(yùn)動(dòng)的船舶螺旋槳水動(dòng)性能計(jì)算方法。

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

對于具有復(fù)雜運(yùn)動(dòng)的船舶螺旋槳水動(dòng)性能預(yù)報(bào)的文獻(xiàn)，國內(nèi)外研究均很少。目前有少量文獻(xiàn)研究具有簡單運(yùn)動(dòng)形式的螺旋槳的水動(dòng)性能。如陶堯森等[5]利用勢流理論近似計(jì)算了波浪導(dǎo)致螺旋槳做升沉運(yùn)動(dòng)時(shí)的水動(dòng)性能。喻欣[6]用Fluent軟件計(jì)算了螺旋槳做升沉運(yùn)動(dòng)時(shí)的水動(dòng)性能。黃苗苗等[7]也用Fluent軟件做了類似的工作。郭春雨等[8]利用試驗(yàn)研究了具有升沉運(yùn)動(dòng)的螺旋槳的水動(dòng)性能。Politis[9-10]利用不定常面元法預(yù)報(bào)了螺旋槳做升沉運(yùn)動(dòng)時(shí)的水動(dòng)性能。Kinnas[11-12]等也利用不定常面元法預(yù)報(bào)了螺旋槳做升沉或者縱向或者扭轉(zhuǎn)振動(dòng)時(shí)的水動(dòng)性能。在所有這些研究中，螺旋槳的運(yùn)動(dòng)形式均比較簡單，基本都限定在單一的平動(dòng)形式。而對于隨軸系振動(dòng)的螺旋槳而言，由于其有6個(gè)空間自由度，使得其運(yùn)動(dòng)軌跡復(fù)雜，不同自由度間還具有耦合現(xiàn)象。因此上面這些文獻(xiàn)中的方法不具有通用性。

針對上述文獻(xiàn)中的不足，本文構(gòu)建了螺旋槳在流場中做復(fù)雜空間運(yùn)動(dòng)時(shí)的水動(dòng)力預(yù)報(bào)數(shù)學(xué)模型，并利用不定常面元法求解了這一問題。利用本文的方法可以方便地研究螺旋槳隨軸系振動(dòng)而做復(fù)雜空間運(yùn)動(dòng)時(shí)的水動(dòng)性能，分析軸系振動(dòng)對螺旋槳水動(dòng)力的影響，從而解決了上文提到的流體-螺旋槳-軸系流固耦合問題中的第三個(gè)關(guān)鍵子問題，為研究流體-螺旋槳-軸系雙向流固耦合動(dòng)力學(xué)特性做基礎(chǔ)。同時(shí)本文的方法也可以用來預(yù)報(bào)船舶轉(zhuǎn)彎、船舶升沉及縱搖振蕩、螺旋槳啟停及加速等復(fù)雜工況下的螺旋槳水動(dòng)性能。

1 數(shù)學(xué)模型

目前已有多種成熟方法預(yù)報(bào)螺旋槳水動(dòng)力，如升力線理論、升力面理論、面元法（又稱邊界元方法BEM）、CFD法等。升力線理論含有大量假設(shè)，因而只適用于輕載螺旋槳。升力面理論雖然可以較為準(zhǔn)確地預(yù)報(bào)總推力和總扭矩，但由于沒有考慮螺旋槳葉片的厚度，因而其預(yù)報(bào)的槳葉面壓力分布不夠準(zhǔn)確。CFD方法盡管能準(zhǔn)確計(jì)算槳葉表面壓力分布，但是由于其計(jì)算量大，耗費(fèi)時(shí)間長，因而應(yīng)用并不廣泛。相對升力面理論來說，面元法考慮了葉片厚度的影響，因此不僅可以準(zhǔn)確計(jì)算總推力與總扭矩，而且可以較準(zhǔn)確地計(jì)算槳葉表面的壓力分布；相對于CFD方法來說，面元法計(jì)算量小，耗時(shí)短，同時(shí)相比于早期來說，面元法有了很大的改進(jìn)，應(yīng)用范圍越來越廣，比如可以考慮空化現(xiàn)象，結(jié)合邊界層理論還可以考慮流動(dòng)分離等復(fù)雜流動(dòng)現(xiàn)象。因此本文采用基于擾動(dòng)速度勢的面元法來預(yù)報(bào)隨軸系做復(fù)雜空間運(yùn)動(dòng)的螺旋槳的水動(dòng)性能。

要建立隨軸系振動(dòng)的螺旋槳水動(dòng)力預(yù)報(bào)數(shù)學(xué)模型，需要解決三個(gè)問題：一是建立軸系復(fù)雜空間運(yùn)動(dòng)的數(shù)學(xué)模型；二是建立具有復(fù)雜空間運(yùn)動(dòng)固體的流體邊界條件；三是建立時(shí)變的尾渦數(shù)學(xué)模型。

1.1 軸系復(fù)雜空間運(yùn)動(dòng)的數(shù)學(xué)模型

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

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

圖3 投影角示意圖Fig.3 The schematic of projected angle

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

從而坐標(biāo)系OXYZ與坐標(biāo)系oxyz的關(guān)系可表示為：

1.2 流體邊界數(shù)學(xué)模型

式中：Vin表示螺旋槳的進(jìn)流速度。當(dāng)Vs與ωs為零［T］為單位陣時(shí)，其與傳統(tǒng)的不考慮軸系振動(dòng)的螺旋槳進(jìn)流速度一致。r表示螺旋槳表面某點(diǎn)的位置向量，在oxyz中表示。上述式子假設(shè)螺旋槳葉片剛性，不發(fā)生彈性變形，如果葉片有彈性變形，則還需要考慮其變形速度。

假設(shè)：（1）流體無粘、無旋且不可壓縮；（2）螺旋槳浸水足夠深，即不考慮自由液面影響，同時(shí)流體域延伸到無限遠(yuǎn)；（3）不考慮空化影響；（4）軸系振動(dòng)為小幅振動(dòng)，葉片上沒有渦的分離等復(fù)雜流動(dòng)現(xiàn)象。取一足夠大的外部控制面將其封閉在內(nèi)。如圖4所示。

流域的邊界面由物面SB，尾渦面SW和外邊界面S∞組成。在該流場中可用擾動(dòng)速度勢φ來表示螺旋槳的擾動(dòng)。在oxyz坐標(biāo)系中，Φ滿足Laplace方程：

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

式中：nQ是邊界面上的單位法向量，在oxyz坐標(biāo)系中表示。

圖4 螺旋槳及周圍流場示意圖Fig.4 Propeller and fluid around it

1.3 時(shí)變尾渦數(shù)學(xué)模型

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

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

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

圖5 尾渦泄露過程Fig.5 The process of vortex-shedding

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

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

2 求解過程

本文采用基于擾動(dòng)速度勢的面元法求解（4）式與（5）式組成的定解問題。根據(jù)Green定理，當(dāng)場點(diǎn)P（x,y,z）在物面上時(shí)，擾動(dòng)速度勢可以表示為（結(jié)合邊界條件）

式中，Δφ（Q1）為尾渦面的上速度勢跳躍，按（6）式確定。

對于復(fù)雜形狀結(jié)構(gòu)，（8）式很難求解析解，因此采用數(shù)值面元法求解。假設(shè)螺旋槳共有Z個(gè)葉片，將一個(gè)葉片及相應(yīng)輪轂劃分成Np個(gè)四邊形面元（沿弦向面元數(shù)為N，沿展向面元數(shù)為M，輪轂面元數(shù)為Nh，則Np=N·M+Nh）。泄露尾渦面元的展向數(shù)目為M，弦向數(shù)目由時(shí)間總步數(shù)Nt決定。比如在it個(gè)時(shí)間步，弦向泄露共it個(gè)尾渦面元，其中只有緊靠葉片隨邊的第一列面元強(qiáng)度未知，其它尾渦面元強(qiáng)度均已知。如圖6所示。

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

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

由此可知，（9）式右邊各項(xiàng)均為已知量，方程可以求解。進(jìn)一步，（9）式寫成矩陣形式：

式中：“±”表示隨邊處上下面元。

在oxyz坐標(biāo)系中，不定常Bernoulli方程可表示為：

式中：P0為參考點(diǎn)處流體壓力；ρ為流體密度；V＝Vin+▽φ，為總擾動(dòng)速度。

求出葉片上速度勢分布后，即可由（13）式求出葉片表面的壓力分布，結(jié)合結(jié)構(gòu)動(dòng)力學(xué)，進(jìn)而可以求解螺旋槳葉片變形、軸承力等等。整個(gè)求解流程如圖7所示。

圖7 求解流程圖Fig.7 The process of solution

3 數(shù)值算例

本節(jié)針對實(shí)際使用中的螺旋槳給出兩個(gè)算例，以驗(yàn)證算法的準(zhǔn)確性，同時(shí)預(yù)報(bào)軸系振動(dòng)對螺旋槳激勵(lì)力的影響。從求解流程圖7可以看出首先必須確定進(jìn)流速度Vin，而由（3）式可知，需要知道軸系振動(dòng)速度Vs與ωs，才能求出進(jìn)流速度Vin。由于實(shí)際中軸的振動(dòng)形式比較復(fù)雜，一方面受多種激勵(lì)源影響（比如不平衡激勵(lì)，彎曲激勵(lì)，流體激勵(lì)等等）；另一方面激勵(lì)力與振動(dòng)之間還存在耦合效應(yīng)，比如流體激勵(lì)力引起軸系振動(dòng)，而軸系振動(dòng)又對流體激勵(lì)力有反饋?zhàn)饔?。因此需要建立一個(gè)完整的流體-槳-軸系流固耦合動(dòng)力學(xué)模型，才能較為準(zhǔn)確地預(yù)報(bào)軸系振動(dòng)，這已經(jīng)超出本文的研究范圍。為了研究方便，本文假設(shè)軸系的振動(dòng)形式已知，為簡單的簡諧振動(dòng)。本文以4381螺旋槳為研究對象，其為5葉片無側(cè)斜螺旋槳，設(shè)其直徑為5 m，其余詳細(xì)幾何參數(shù)在參考文獻(xiàn)[16]中給出。

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

第二個(gè)算例是研究螺旋槳在縱向與橫向（垂向）聯(lián)合振動(dòng)下的脈動(dòng)力變化。設(shè)軸系末端的振動(dòng)速度為（單位為m/s）：

式中ω既是軸的振動(dòng)頻率，也是螺旋槳的旋轉(zhuǎn)頻率。因此可以認(rèn)為軸系按轉(zhuǎn)頻振動(dòng)。計(jì)算中設(shè)ω為3 Hz。假設(shè)軸系縱向振動(dòng)幅值為1 mm，根據(jù)簡諧振動(dòng)理論，軸系縱向振動(dòng)速度大約為0.018 m/s。

圖9 脈動(dòng)推力系數(shù)與扭矩系數(shù)Fig.9 The pulsant thrust and torque coefficients

圖10 脈動(dòng)推力與扭矩Fig.10 The pulsant thrust and torque

圖11 脈動(dòng)垂向力與橫向力Fig.11 The pulsant vertical force and lateral force

螺旋槳轉(zhuǎn)動(dòng)一圈的過程中計(jì)算360步，一共計(jì)算6圈。圖9為計(jì)算的推力系數(shù)與扭矩系數(shù)隨旋轉(zhuǎn)角度的變化過程。從圖中可以看出，由于螺旋槳從靜止啟動(dòng)，在啟動(dòng)瞬間，螺旋槳加速度無窮大，在槳葉表面瞬時(shí)產(chǎn)生一個(gè)啟動(dòng)渦，導(dǎo)致槳葉表面的脈動(dòng)壓力變化非常劇烈。在螺旋槳旋轉(zhuǎn)大約兩圈后，槳葉表面的壓力脈動(dòng)呈現(xiàn)周期性變化，表明此時(shí)計(jì)算已經(jīng)趨于穩(wěn)定。計(jì)算結(jié)果表明在（14）式的速度假設(shè)下，穩(wěn)定后的脈動(dòng)推力呈周期變化，變化頻率為軸頻，其推力脈動(dòng)分量大約為靜態(tài)分量的4.5/1 000，扭矩脈動(dòng)分量大約為靜態(tài)分量的4/1 000。對于大型船舶軸系來說，縱向振動(dòng)1 mm的量級是非常普遍的，由此可見其產(chǎn)生的脈動(dòng)力將對軸系振動(dòng)產(chǎn)生重要影響，因而不能忽略。

圖10為計(jì)算穩(wěn)定后，葉片旋轉(zhuǎn)一圈的脈動(dòng)推力與扭矩變化。圖11為脈動(dòng)的垂向力與側(cè)向力變化。此處，脈動(dòng)力均去掉了靜態(tài)分量，只考慮其脈動(dòng)部分。從圖中可以看出，單個(gè)葉片上旋轉(zhuǎn)一圈的載荷不是按簡諧變化的，表明此時(shí)有軸頻以外的其它頻率成分存在，進(jìn)一步研究后表明是軸頻的倍頻成分。此頻率的來源是（13）式中計(jì)算脈動(dòng)壓力時(shí)，對速度進(jìn)行平方而引入的。但是所有葉片的合力又是一個(gè)簡諧激勵(lì)力，表明不同葉片的倍頻分量相位差滿足嚴(yán)格的滯后關(guān)系，導(dǎo)致它們可以相互抵消。又由于所有葉片脈動(dòng)力的軸頻分量的相位均一致，使得總的合力有所增強(qiáng)。

4 結(jié) 語

本文構(gòu)建了螺旋槳在流場中隨軸系做復(fù)雜空間運(yùn)動(dòng)時(shí)的水動(dòng)力預(yù)報(bào)數(shù)學(xué)模型，并利用不定常面元法求解了這一問題。通過與已有文獻(xiàn)的研究結(jié)果做比較，證明了本文算法的有效性。利用本文的方法可以方便地研究螺旋槳隨軸系振動(dòng)而做復(fù)雜空間運(yùn)動(dòng)時(shí)的水動(dòng)性能，分析軸系振動(dòng)對螺旋槳水動(dòng)力的影響，從而解決了上文提到的流體-螺旋槳-軸系流固耦合問題中的第三個(gè)關(guān)鍵子問題，為研究流體-螺旋槳-軸系雙向流固耦合動(dòng)力學(xué)特性做基礎(chǔ)。同時(shí)該方法也可以用來預(yù)報(bào)船舶轉(zhuǎn)彎、船舶升沉及縱搖振蕩、螺旋槳啟停及加速等復(fù)雜工況下的螺旋槳水動(dòng)性能。最后預(yù)報(bào)了軸系縱向振動(dòng)幅值1 mm，振動(dòng)頻率為3 Hz并伴有小量回旋振動(dòng)下4381螺旋槳的脈動(dòng)力，結(jié)果表明推力脈動(dòng)分量大約為其靜態(tài)分量的4.5/1 000，扭矩脈動(dòng)分量大約為其靜態(tài)分量的4/1 000。因此軸系振動(dòng)對螺旋槳脈動(dòng)力的影響不能忽略。