質量靜矩和慣性矩對水翼流致振動及噪聲影響的數值研究

劉 龍， 張懷新， 姚慧嵐

(上海交通大學 海洋工程國家重點實驗室 高新船舶與深海開發裝備協同創新中心，上海 200240)

在船舶行業中，有很多水翼[1-3]結構，像減搖鰭、舵、水翼艇的滑行水翼、潛艇的圍殼舵和尾翼等都是典型的水翼結構。水翼結構在水中前行時會產生升力和力矩，如果設計得當可以用翼型的受力特點來改善船舶的水動力性能，應該注意到的是，翼型結構自身可以看成是一個彈性體，在流場力的作用下會產生振動和變形，當流體作用力較小或者說翼型結構強度較大時，這種振動會隨時間趨于收斂，但是，當流場作用力達到一種臨界狀態時，翼型會產生發散或者顫振，發散是指速度大于一定值時，翼型的位移會不斷增大，不圍繞平衡點做往復振動的情況；顫振是指速度增大到一定值時，振幅不衰減，并有增大趨勢。現在，船舶正向高速化發展，當翼型的內部結構(主要是質量靜矩和質量慣性矩)設計不合理時，可能會影響水翼的性能和結構安全，另外，潛艇上翼型結構的振動產生的噪聲也會威脅其隱身性能。因此，研究水翼的振動和噪聲具有重要意義。

流致振動[4-7]是十分常見的問題，在翼型結構的流致振動問題的研究中，人們更多關注的是空氣中機翼的流致振動，Theodorsen[8]提出的非定常氣動理論可以快速求解顫振問題，從而它在機翼的氣動力問題上得到了廣泛應用，但是，由于它采用了簡諧振動的假定，所以它不能求解發散問題。目前國內外對水翼的研究還比較少，水翼的研究主要集中在空泡[9-11]問題上。相對而言，彈性水翼流致振動問題的研究很少，試驗研究主要集中在上個世紀。Woolston等[12]對NACA16-010翼型的穩定性做了試驗，研究了相對質量比(質量比的平方根)和臨界速度(剛好產生發散或顫振時的速度)的關系，結果表明，相對質量比小于3時，會產生發散；相對質量比大于3時，會發生顫振。Besch等[13]針對水翼，試驗研究了NACA16-012翼型在四種不同質量比狀態下的臨界速度，結果表明，質量比為0.963的模型的顫振速度為24.7 kn，而另外三種質量比相對較小的模型沒有發生顫振，只是出現了靜態失效(靜態發散)，靜態失效速度為36 kn。Ducoin等[14]試驗研究了NACA66312水翼由于邊界層轉捩而產生的振動。結果表明，邊界層轉捩會產生較大的振動，而振動的頻率和幅度取決于渦脫落的頻率。在數值研究上，劉曉宙等[15]進行了機翼的渦激振動和聲輻射的研究，計算側重于中、低雷諾數。計算結果表明，翼型在大攻角，當系統的頻率恰當時，會出現共振；流體繞彈性翼引起的聲輻射大于流體繞固定翼引起的聲輻射，當彈性翼發生共振時，聲輻射達到最大值。Chae等[16]根據Woolston等的試驗數據，采用k-ωSST湍流模型數值計算了二維翼型NACA16-010的流致振動問題，數值計算結果和試驗結果相吻合。Akcabay等[17]研究了兩個自由度二維彈性水翼在空泡情況下的振動和受力特征。劉胡濤等[18]采用大渦模擬的方法數值計算了具有兩個自由度的二維彈性水翼NACA0012的振動和聲輻射問題，他在文章中指出，隨著流速的增加，兩個自由度振動的頻率會趨于接近，最后會發生顫振，水翼振動主要影響低頻噪聲。對于水翼噪聲問題，通常包括流噪聲和流激噪聲，流噪聲由水翼壁面偶極子噪聲和流場中的湍流噪聲(四極子噪聲)組成，在湍流不強烈的情況下通常不考慮四極子噪聲。Kim等[19]采用NACA0018翼型計算了該翼型的渦脫落特性和噪聲問題，計算結果表明，渦脫落頻率可以在聲壓頻譜圖上捕捉到；偶極子噪聲明顯大于四極子噪聲。

從前文的敘述可知，Theodorsen非定常理論在翼型顫振問題上運用廣泛，它可以快速在頻域里提供一個理論解，對水翼問題有一定的參考意義，但是它也有很多局限性，如不能求解發散問題和沒有考慮黏性。基于Theodorsen 非定常理論的不足之處，本文將采用流固耦合的方法計算水翼的流致振動和聲輻射問題，流場在Fluent里計算，運動方程在UDF(User Defined Function)里用四階龍格庫塔法求解，在每個時間步里，用UDF控制流場力和水翼位移的相互傳遞，在時域里模擬水翼的振動。另一方面，質量比會影響翼型振動的臨界狀態，而實際問題中，翼型的質量比往往是一定的，在這種情況下，要想提高翼型的臨界速度就得從翼型的內部結構入手。翼型的內部結構關系到翼型的剛性位置和重心位置，更具體的說法是翼型的內部結構會影響翼型的質量靜矩和慣性矩。由于以往研究者對水翼的質量靜矩和慣性矩研究較少，所以，本文將采用NACA16-010翼型數值計算具有兩個自由度的三維剛性水翼的流致振動和聲輻射問題，重點研究質量靜矩、質量慣性矩對水翼振動以及臨界狀態的影響，并分析不同振動情況下的聲輻射問題。

1 數值計算

1.1 物理模型

本文將研究兩自由度三維剛性水翼(只考慮結構運動，不考慮結構變形)在不可壓縮流體中的流致振動問題，流體介質為水，三維水翼具有兩個自由度，分別對應垂向運動和繞彈性軸的扭轉運動。簡化模型如圖1所示。

圖1 兩自由度三維水翼模型Fig.1 The three-dimensional hydrofoil with 2DOF

水翼弦長c=2b，展長為L，來流速度為v，初始攻角為α0。E.A為三維水翼彈性中心，位于翼弦中點后ab處，a為彈性中心到水翼弦長中點距離的無量綱量，定義彈性中心靠近隨邊，a為正弦長，b為半弦長。翼型在彈性軸處由一個線彈簧和一個扭轉彈簧支撐，垂向位移h(t)與扭轉角θ(t)均以彈性軸為參考點測量的，h(t)向上為正，θ(t)逆時針為正，在模態上分別對應水翼的第一彎曲模態和第一扭轉模態。H.C為水動力中心，是流體作用力合力的作用點，位于彈性中心前緣eb處。重心C.G位于彈性中心后緣xθb處，xθ為重心到彈性軸的無量綱距離。

Woolston等試驗研究了質量比對翼型穩定性的影響，翼型為NACA 16-010，材料為桃木，弦長c=0.305 m，展長L=1.219 m。試驗的流體介質是空氣和氟利昂的混合物，通過改變空氣和氟利昂的比例來改變流體的密度，質量比μ=m/πρb2L,m為機翼的質量，ρ為流體介質的密度,在本次計算中，流體介質為水，水的密度不變，質量比是通過改變水翼的質量來實現的，模型的具體參數，如表1所示。

表1 Woolston和本文的模型參數Tab.1 Physical parameters of Woolston and this paper

在本文的驗證工作中，參數選取與Woolston等試驗模型的參數基本一致，在后面計算質量靜矩和慣性矩對振動的影響時，主要計算質量靜矩半徑xθ和慣性矩半徑rθ對振動以及噪聲的影響。

1.2 兩個自由度剛體的運動方程

(1)

式中：{Ffluid}=[LfluidMfluid]T，Lfluid，Mfluid分別為水翼受到的升力和繞彈性軸的力矩；{X}=[hθ]T,h為水翼的垂向運動位移，定義向上為正；θ為水翼繞彈性軸的扭轉角，定義逆時針為正。

(2)

水翼的運動方程組是一個質量耦合方程組，水翼兩個自由度的方程實際上是靠xθ和rθ耦合在一起，因此，xθ和rθ會對水翼的振動產生較大影響。求解水翼的振動必須先求解水翼的受力，水翼的受力通常可以用理論公式或CFD(Computational Fluid Dynamics)計算，理論方法一般用Theodorsen非定常理論在頻域里求解，而CFD可以考慮黏性的影響，是在時域里采用流固耦合的方式求解水翼振動的時間歷程信息。

1.3 流固耦合方法求解水翼的運動方程

本文采用時域法，用流固耦合的方式模擬水翼的流致振動，流體介質是水，流固耦合方式是在Fluent里用UDF實現的，流場在Fluent里求解，水翼的運動方程用UDF求解，并用UDF控制流場力和水翼位移的相互反饋。流場的計算是求解非定常的URANS方程，對于當前的計算，使用k-ωSST湍流模型可以達到一個較高的計算精度。

振動水翼在流場中的計算是一個相互影響的問題，由于本文的計算將水翼作剛性體考慮，不考慮結構變形，只考慮水翼位移，因此，在計算的過程中只需求解流體方程和剛體的運動方程，流體方程和剛體運動方程是交替求解的。耦合計算過程是在Fluent里采用UDF以動網格的形式完成，耦合計算過程如下：首先，Fluent迭代流場至穩定狀態，將計算出的流場力傳遞給UDF，然后UDF根據水翼網格節點受力情況求出總的升力Lfluid和相對于彈性軸的力矩Mfluid，并求解兩個自由度的運動方程，得到水翼的位移、速度和加速度并將位移傳回流場中，流場網格更新，重新迭代至穩定狀態，重復這個過程直到計算時間結束。兩個自由度運動方程是采用四階龍格庫塔法求解的，對于流固耦合問題，流場和固體場要相互反饋信息，時間步長對計算精度有較大影響，在計算中將采取小步長和多迭代次數來保證計算精度。具體的計算過程如圖2所示。

圖2 計算流程圖Fig.2 Coupled algorithm flow chart

流固耦合方法是在時間范圍內模擬水翼振動的，能得到水翼的振動曲線，可求解水翼的一般振動和發散問題。臨界速度是指剛好出現顫振或發散時的來流速度，判斷依據是通過觀察水翼振動的時間歷程信息得到，當速度增加到一定值時，剛好出現位移和力的振動曲線的振幅隨時間不發生變化，此時的速度為顫振的臨界速度；當剛好出現位移和力隨時間不斷增大，不產生簡諧運動，此時對應的速度為發散臨界速度。

2 計算結果分析

2.1 網格無關性與時間步長的驗證

流場的計算網格如圖3所示。采用結構網格和非結構網格結合的形式，邊界層采用結構網格以保證計算精度，遠場網格采用結構網格以減少網格數量，邊界層與遠場之間采用非結構網格以實現網格的變形。邊界層的y+接近1，并對遠場的尾流區進行了加密。

圖3 流場計算網格Fig.3 Mesh of hydrofoil

圖4 網格無關性的驗證Fig.4 The test of mesh

對于流固耦合的數值計算，時間步長對計算結果影響很大，如果時間步長偏大，會帶來較大的計算誤差，甚至有可能出現完全不同的計算結果，因此，有必要驗證時間步長的合理性。在時間步長的驗證中，選取三個時間步長計算，三個時間步長分別為dt=0.001，滿足T/dt>50；dt=0.000 5，滿足T/dt>100；dt=0.000 2，滿足T/dt>200。計算結果如圖5所示。從圖5可知，相對于網格數量，時間步長的影響更大，特別是對扭轉角的計算。理論上講，小的時間步長可以得到更高的計算精度，dt=0.000 5與dt=0.000 2的計算結果比較接近，考慮到計算耗時問題，在后面的計算中，一般都選取dt=0.000 5為計算時間步長，但對于流速較大的計算，采用dt=0.000 2來計算，以保證計算精度。

圖5 時間步長的驗證Fig.5 The test of time step

2.2 臨界速度與質量比的關系

圖的振動曲線(收斂狀態)Fig.6 Vibration curves of

圖的振動曲線(臨界狀態)Fig.7 Vibration curves of (critical state)

圖的振動曲線(顫振狀態)Fig.8 Vibration curves of

圖的振動曲線(收斂狀態)Fig.9 Vibration curves of

圖的振動曲線(發散狀態)Fig.10 Vibration curves of

圖11 臨界速度與相對質量比的關系Fig.11 The relationship between critical velocity and added mas rate

v-g法可以在頻域里快速的提供一個理論解，對于高質量比，v-g法提供的理論解具有參考價值，但它不適用于發散問題，且不能提供振動的時間歷程信息，對于水翼問題，一般質量比都較小，出現發散的概率大，同時，對于水動力問題，黏性也是不可忽略的因素。流固耦合的方法很好的彌補了v-g法的缺點，可以更真實的模擬水翼的振動，這也說明用流固耦合方法求解水翼振動問題是非常有意義且必要的。

2.3 質量靜矩和慣性矩對振動狀態的影響

在計算質量靜矩對水翼流致振動的影響時，慣性半徑保持不變，rθ=0.403，在U=1.4時，計算了三種狀態，分別對應的無量綱靜矩半徑xθ為0.068，0.1和0.2。質量靜矩對振動周期基本沒有影響，其主要影響振動的幅值，質量靜矩增大一定程度時，水翼的振幅會增加，當靜矩增加過大時，振動反而衰減的更快，這可能是質量分布不對稱，扭轉會衰減的更快，從而加速整體的衰減，如圖12所示。在計算xθ=0.1和xθ=0.2兩種情況的臨界狀態時發現，xθ=0.1對應的臨界速度Uf=1.53，其臨界狀態的振動曲線，如圖13所示。在來流速度U=1.53時，水翼振動剛好處于臨界狀態，振動不增大也不減小；xθ=0.2對應的臨界速度Uf=1.59，其臨界狀態的振動曲線，如圖14所示。前面計算的xθ=0.068所對應的臨界速度Uf=1.7，由此發現，質量靜矩會影響水翼的臨界速度。U=1.6時，對比xθ=0.1和xθ=0.2的振動曲線發現，在顫振狀態下，質量靜矩對垂向位移影響較小，而對扭轉位移影響較大，xθ=0.2的扭轉幅值明顯比xθ=0.1的大,如圖15所示。

在計算質量慣性矩對水翼流致振動的影響時，靜矩半徑保持不變，xθ=0.068，在U=1.7時，計算了三種狀態，分別對應rθ=0.35，rθ=0.403和rθ=0.45。對于垂向位移，慣性矩主要影響振動頻率，rθ增大，振動頻率會降低；對于扭轉角度，慣性矩主要影響振動頻率和振幅，rθ增大，振動頻率會減小，其減小趨勢比垂向振動頻率減小趨勢要大，當rθ過大時，振幅會變小，如圖16所示。由此可見，減小rθ會使fh與fθ的差距增大，因而要達到顫振臨界狀態(即fh=fθ)，就必須使速度增加。U=1.8時，水翼剛好發生顫振，說明rθ=0.35的臨界速度Uf=1.8，如圖17所示。在計算rθ=0.3的臨界速度時，發現在U=1.95時，剛好出現了發散而沒有出現顫振，這可能是由于fh和fθ差別過大，在兩者頻率相等之前，水翼的位移過大，以致發生破壞，如圖18所示。

圖12 質量靜矩對振動的影響,U=1.4Fig.12 The impact of mass moment to the flow induced vibration,U=1.4

圖13 xθ=0.1,Uf=1.53對應的振動曲線Fig.13 Vibration curves of xθ=0.1,Uf=1.53

圖14 xθ=0.2，Uf=1.59對應的振動曲線Fig.14 Vibration curves of xθ=0.2，Uf=1.59

圖15 xθ=0.1,xθ=0.2,U=1.6對應的振動曲線Fig.15 Vibration curves of xθ=0.1,xθ=0.2,U=1.6

圖16 質量慣性矩對振動的影響，U=1.7Fig.16 The impact of mass moment of inertia to the flow induced vibration U=1.7

圖17 rθ=0.35,Uf=1.8的振動曲線Fig.17 Vibration curves of rθ=0.35,Uf=1.8

圖18 rθ=0.3,U=1.95的振動曲線Fig.18 Vibration curves of rθ=0.3,U=1.95

3 振動水翼的噪聲分析

圖19 固定翼和振動翼的聲幅射對比圖，U=1.7Fig.19 The sound radiation characteristics of fixed and flexible hydrofoil,U=1.7

研究質量靜矩對噪聲的影響時，其他參數不變，計算了xθ=0.068和xθ=0.1在U=1.7時的噪聲。從前面的計算結果可知，xθ=0.068(其臨界顫振曲線見圖7)和xθ=0.1(其臨界顫振曲線見圖13)所對應的臨界速度分別是Uf=1.7和Uf=1.53。在U=1.7時，兩者都處于顫振狀態，從圖20(a)可知，在首尾兩點的聲壓值比較接近，但是，在其他位置，xθ=0.1的聲壓值都比xθ=0.068的聲壓值大。從圖20(b)可知，xθ=0.1的主峰值比xθ=0.068的主峰值稍微大一點，對于次波峰值，xθ=0.1的明顯比xθ=0.068的大。前面計算質量靜矩時發現，xθ=0.1對應的臨界速度低且振幅大，綜合噪聲的計算結果，可以說明，在相同條件下，振幅越大，其對應的噪聲也越大。

圖20 xθ=0.068,xθ=0.1臨界狀態下的聲輻射對比圖，U=1.7Fig.20 The sound radiation characteristics of the hydrofoil (xθ=0.068，xθ=0.1),U=1.7

研究慣性矩對噪聲的影響時，其他參數不變(見表1)，計算了水翼慣性半徑為rθ=0.35和rθ=0.405在速度U=1.7時的噪聲，從前面振動計算結果可知，rθ=0.35對應水翼的臨界速度Uf=1.8，其臨界顫振曲線，如圖17所示，則，U=1.7時，rθ=0.35對應水翼的振動為收斂；rθ=0.405對應的臨界速度Uf=1.7，臨界顫振曲線，見圖7。對比臨界狀態和收斂狀態的噪聲信息，如圖21所示。臨界狀態下的噪聲比收斂狀態下的噪聲大，收斂狀態下的聲壓頻譜圖上只有一個波峰，對應其振動頻率，而臨界狀態下的聲壓頻譜圖上有兩個波峰，主波峰的峰值比收斂狀態下的峰值高。

圖21 臨界狀態和收斂狀態的聲幅射對比圖，U=1.7Fig.21 The sound radiation characteristics of the hydrofoil,U=1.7

4 結 論

本文采用流固耦合的方法，用k-ωSST湍流模型數值模擬了兩自由度三維剛性水翼NACA16-010的流致振動和聲輻射問題。在一定范圍內，CFD法計算水翼的臨界速度與Theodorsen理論值和試驗值比較接近，在水翼發散臨界速度的計算上，Theodorsen理論會出現較大偏差，其原因是發散問題不滿足Theodorsen理論解中的簡諧運動假定。對于小質量比的水翼問題，出現發散的概率大，故采用流固耦合的方法求解水翼流致振動問題是非常有意義且必要的。

本文重點計算了質量靜矩和慣性矩對振動的影響，從計算結果上看，靜矩對振動頻率影響較小，但對振動幅值影響較大，特別是對扭轉角的幅值，同時臨界速度也會隨靜矩發生變化；質量慣性矩對振動的頻率和幅值都有影響，并且還會影響水翼的臨界速度，有的情況下，慣性矩的改變還會使顫振水翼變成發散水翼。原因是水翼的運動方程是質量耦合的，質量靜矩和慣性矩之間的關系對水翼的振動有重要影響。

最后，本文計算了水翼的聲輻射問題，計算結果表明，水翼振動會使偶極子噪聲顯著增加，并且振動越劇烈，噪聲增加的越多，在聲壓頻譜圖上會有一個與顫振頻率相對應的峰值。