基于分層模型的功能梯度輸流管道耦合振動

朱竑禎， 王緯波， 殷學文， 高存法

(1.南京航空航天大學 航空學院 機械結構力學及控制國家重點實驗室，南京 210016；2.中國船舶科學研究中心 船舶振動噪聲重點實驗室，江蘇 無錫 214082)

功能梯度材料(Functionally Graded Material, FGM)的概念最早由日本科學家于1984年提出[1]，是指材料屬性隨空間位置呈連續梯度變化的新型復合材料。隨著材料制備工藝的精進，功能梯度材料目前已經在耐磨損[2]、緩解內部熱應力及抗斷裂[3]等應用中展現出傳統材料難以超越的獨有優勢。針對不同的使用需求，功能梯度材料具有很強的可設計性，因此在未來工業及工程領域具有極大的潛在應用前景。

輸流管道作為一種介質運輸方式，廣泛應用于水利、化工、航空及海洋工程等領域中。傳統管道一般采用金屬材料，例如潛艇中布置了大量的金屬管道承載各種低速液體的流動，包括海水以及各種機器運轉產生的廢水等。經過時間的積累，金屬極易受到這些液體的腐蝕，從而造成金屬管道內部的堵塞，因此常常需要斥巨資用于管道的防腐蝕處理和定期維護以保證其使用年限[4]。此外傳統的金屬材料聲學阻尼性能差，流體與管壁相互碰撞擠壓不僅引起耦合振動，增大噪聲，甚至還可能導致管道的疲勞破壞。新型復合材料具有優良的抗疲勞、阻尼特性以及很強的設計性，為工業設計引入新風，受到廣泛關注[5]。目前已有許多國外學者研究用新型復合材料制作水下螺旋槳，既耐腐蝕，且預期能比傳統材料有更好的聲學性能[6-8]。但一般復合材料層間剪切強度和拉伸強度低[9]，而且極易產生分層損壞從而導致復合材料的失效。而功能梯度材料則由于其材料特性連續變化，材料組份之間無顯著界面，因此與傳統材料和復合材料相比在優化應力分布方面有明顯的優勢。船舶管道系統的可靠性直接影響了其內部各儀器的正常運作與配合，因而本文考慮用功能梯度材料這種新型材料取代傳統金屬材料和復合材料用于輸流管道。

在工程中，相對截面尺寸而言，一般管道長度較大，因此多采用梁模型描述管道的振動。早期的理論中大多將內部流體視為附加質量，將流速和壓力視為常量，忽略流體與管道的耦合效果[10-12]，對于輸送液體的管道的描述是不精確的[13]。隨耦合理論的深入發展，William[14]和Walker等[15]考慮了流體在軸向的可壓縮性，Wiggert等[16]和Lesmez等[17]計入管道的泊松耦合效應影響，Lee等[18]分析了重力作用和流體的黏性，Gorman等[19]引入了管道的徑向變形和初始預應力，傳統材料輸流管道的振動理論已經發展得較為成熟。盡管對于功能梯度材料的理論研究不少，但大多僅限于純結構的研究，如功能梯度梁[20-22]、功能梯度板的數值分析等[23-25]，而目前對于功能梯度輸流管道方面的研究主要通過圓柱殼[26]、Euler梁[27-34]或Timoshenko梁[35]模型研究管內流速導致管道失穩的問題。圓柱殼的模型雖然從三維角度分析管道更詳細準確，但是難以模擬復雜的組裝管道。而Euler梁模型忽略了截面的剪切變形，只能適用于薄壁細長管道，難以滿足工程各類計算需求。此外在處理管道失穩問題中，一般只考慮管道的橫向方程，對于流體的處理是較為簡單的，未考慮流體在軸向的可壓縮性，甚至只將流體視為附加質量。而實際上對于管道這種鏈式結構來說，其軸向和流體的運動影響了管道上下游振動與聲的傳遞特性。為了分析功能梯度管道對于管道結構和流體的影響，將流速和壓力設為變量，計入流體在管道軸向的可壓縮性，基于Timoshenko梁模型，綜合考慮了功能梯度材料管壁與流體的相互作用，通過Hamilton原理及流體動量、連續方程得出功能梯度輸流管道的耦合振動方程。

功能梯度管道的材料屬性沿厚度方向變化，本文采用離散模型來描述，將管壁均勻劃分若干層，每一層近似為均勻材料，以此模擬連續梯度變化。動剛度法是基于單位內控制方程的精確解，與有限元法相比具有較高的計算效率[36]，因此本文利用動剛度法進行數值求解，然后通過算例分析了梯度指數的變化對于結構模態、流體模態及時域響應的影響，為將來功能梯度材料的工程應用提供一定理論依據。

1 輸流管道振動方程

1.1 管道運動方程

本文計算的輸流管道為圓形截面，管道徑向坐標記為r，截面內半徑為r1，外半徑為r2。材料參數如楊氏模量、密度及泊松比均沿半徑呈冪律變化[37]。

式中：Q(r)表示沿厚度方向呈連續梯度變化的材料參數，包括了楊氏模量E，密度ρ和泊松比μ等。下標1和2分別表示內界面和外界面的材料參數，n為梯度指數，表征材料的體積分布規律。

本文采用層合近似模型[38-40]，將這種連續梯度變化的功能梯度材料視為由多層連續梯度材料疊加成的層合材料。如圖1所示，沿厚度方向把截面均勻劃分為K層，每一層視為均勻材料，以此近似模擬沿半徑方向變化的材料屬性。對于第m層，該層的材料參數取為中間位置的值，記為Qm，寫為如下的函數形式：

(1)

圖1 管道截面徑向分層模型Fig.1 Multi-layered model on the cross section of pipe

式中：rm1和rm2分別為第m層的內外半徑。當需要考慮阻尼的影響時，可用復彈性模量E(1+iη)代替原彈性模量，其中η=2ξ，ξ為阻尼比，i為虛數單位。

工程上一般采用梁模型來描述管道，在Timoshenko梁模型中需考慮梁的剪切變形。以x和z方向分別表示管道的軸向和橫向，對于梁上任意一點，軸向位移和橫向位移表示為：

(2)

式中：u和w表示梁的中心面(z=0)的位移，φ(x,t)為截面的法線轉角，t為時間。由幾何方程可得應變為：

(3)

式中：微分符號(′)=?/?x。

由本構方程，應力表示為：

(4)

將應力式(4)在面內積分即可得到橫截面上的軸力、彎矩和剪力：

(5)

在圖2中，管道與水平面的傾斜角為θ，則重力加速度在軸向和橫向的分量分別為gx=gsin(θ)和gz=gcos(θ)。管內流體對于管壁的單位長度垂直作用力和切向作用力分別為Nf和τ。

圖2 管道微段局部坐標系Fig.2 Local coordinate system of an infinitesimal pipe element.

根據Hamilton原理，對系統泛函取極值：

(6)

式中：T為動能，U為應變能，W為非保守力做功，其具體形式為：

(7)

將式(6)代入式(7)中可得管道運動方程為：

(8)

1.2 流體方程

管道內流體的動量方程與管壁結構無關，可參考已有的管道理論，表示為：

(9)

式中：p為流體壓力，Af為流體截面積(對于充液管道而言，即為內徑面積)，c為流速，mw為單位長度流體的質量，微分符號為(·)=?/?t，(′)=?/?x。

根據流體力學[41]，流體與管壁間的摩擦作用力為：

(10)

此外，管道內流體的連續方程為：

(11)

式中：ρw為管道內流體的密度，根據流體狀態方程有：

(12)

式中：Ev為流體的體積模量。

由于流體對于管壁的擠壓而造成面積的變化，此外由于泊松比的效應，軸向的應變也會引起橫向變形。因此，式中的第二項為[42]：

(13)

式中：ε2為橫向應變，ε1為軸向應變，μ為泊松比，由于泊松比數值較小，影響相對小，因此在該式中取為截面的平均泊松比。

式(13)中的軸向應變為ε1=u′。圖3為圓環截面的受力圖，其中F為內力，p為內流體對于管壁的壓力，根據受力平衡得：

圖3 圓環截面受力示意圖Fig.3 Force diagram of annular section.

因此連續方程化為：

(14)

2 輸流管道振動方程求解

為了對方程組進行簡化，將流速和流體壓力線性化，表示為：

c(x,t)=c0+cd(x,t),p(x,t)=p0+pd(x,t)

式中：c0和p0表示穩態時的定常流速和定常壓力，而cd和pd表示微小的擾動變量。

當管內流體在高速高壓下方程中的非線性項對結構的動態響應才有明顯的影響[43]，在此我們主要考慮工程中的低速流動情況。因此為便于方程的求解，略去非線性項，由流體方程(9)得到Nfφ=mwgzφ。此外由于擾動變量遠小于常量即cd?c0，pd?0，舍去擾動項的高次項，從而消去管道方程式(8)中管壁與流體的相互作用力，結合連續方程式(14)可得功能梯度輸流管道的耦合振動方程組為：

(15)

將式前四個方程寫為矩陣形式：

MΦ=0

式中：Φ=[uwφcd]T，M為微分系數矩陣。

注意到|M|為關于x的八次微分方程，因此將位移未知量設為如下形式：

對于|M|=0求解，可得到八個根，記為λ1～λ8。假設橫向位移函數形式為：

式中：i為虛數單位(i2=-1)，ω為圓頻率。根據方程(15)，其他三個變量與橫向位移的關系為：

式中：C1～C8為常系數，fj，huj，hcj可由方程系數求得，在此不作贅述。然后將軸向位移和流速函數代入式(15)最后一個方程中，可求得壓力變量。

對于一根長為L的管道，將其首末兩端作為單元節點，位移向量寫為：

d=[u(0)w(0)φ(0)cd(0)u(L)w(L)φ(L)cd(L)]T

相應的節點力向量寫為：

f=[N(0)Q(0)W(0)pd(0)N(L)Q(L)W(L)pd(L)]T

顯然內力及壓力函數都可由位移表示，因此節點位移和節點力都可以寫成關于常系數C1～C8的矩陣形式：

d=K1C，f=K2C

(16)

式中：C=[C1C2C3C4C5C6C7C8]T為常數矩陣，K1和K2為與ω相關的系數矩陣。

由式可直接表示出力向量和位移向量的關系：

f=KDd

(17)

式中：KD=K2K1-1，即為動剛度矩陣。

由前文的推導過程可見，動剛度法的精度與劃分單元數量無關，因此對于沿軸向幾何形狀及材料不發生改變的情況下，任意長度的管道都僅需要用一個單元模擬，且對于計算頻域沒有特別限制。

3 算例分析與討論

由數據對比可見，本文的計算結果與文獻基本一致。不過本文計算的橫向模態數值略低，因為本文的理論中考慮了梁的剪切變形，而文獻中采用的是Euler梁模型。雖然本例計算的已為細長梁，但隨著模態階數升高和流速增大，誤差仍會略有增加。注意到在流體模態和軸向模態中，文獻結果的實數部分變化不明顯，而本文的計算結果中考慮到了摩擦因數隨流速的變化，隨流速增大，固有頻率實數部分減小，虛數部分增大。

表1 不同流速下均質輸流管道振動固有頻率Tab.1 Natural frequencies(Hz)for the vibration of homogenous pipes conveying fluid with different velocities Hz

(2)本例計算兩端簡支的功能梯度充液直管，長度為L=6 m，截面內徑為0.2 m，外徑為0.4 m，管壁由鋁-SiC陶瓷材料復合而成，內壁材料為SiC陶瓷，楊氏模量E1=427 GPa，泊松比μ1=0.17，密度ρ1=3 100 kg/m3，外壁材料為鋁，楊氏模量E2=70 GPa，泊松比μ2=0.33，密度ρ2=2 700 kg/m3，在本例中用二十層均勻離散單元。

圖4 在兩端簡支直管道上施加單位橫向載荷Fig.4 The unit transverse load applied on the simply supported straight pipe.

如圖4所示，在管道的L/12處施加單位橫向載荷，前文建立的方程中，軸向與橫向、流體均是耦合的，因此在結構上只施加了橫向載荷即可同時引起軸向及流體響應。圖5繪制了當梯度指數為n=0時點L/12處在頻域0～550 Hz的速度響應，響應曲線的峰值處的頻率對應了固有頻率。軸向及流體響應是由耦合項引起的，因此其振動的幅值小于橫向單位力引起的橫向振動。由圖5可見，流體響應曲線包含了所有模態，而軸向和流體響應曲線上比橫向響應曲線多出的波峰對應了軸向模態(A)及流體模態(F)，以此響應曲線即可區分結構模態和流體模態。

圖5 n=0，點L/12處0-550Hz頻域內速度響應曲線Fig.5 n=0,velocity responses at the point L/12 in frequency domain 0～550 Hz

圖6 不同梯度指數下速度響應曲線Fig.6 Velocity responses with different gradient indexes.

圖6中分別為在不同梯度指數下橫向、軸向及流體的速度響應曲線。由圖線可見，隨n值的增大，結構速度響應曲線均向右移動，即橫向和軸向固有頻率均明顯提高，且隨著階數的升高，改變幅度也變大。而值得注意的是在流速響應曲線中，顯然有一些波峰的位置改變幅度很小，幾乎重合，這些峰值對應的是流體模態，在圖中以F表示。由此可見，管道材料的變化對于管道結構本身的影響較大，而對于其內部的流體影響很小。

(3)工程中常見的薄壁管道，流速增大情況下固有頻率降低，如算例1所示。當第一階固有頻率降低至零時，則管道發生失穩。而由算例2可見，將功能梯度材料應用于管道結構上，能顯著提升固有頻率，因此能有效地延緩因流速導致的管道失穩的發生。圖7展示了將算例2中的功能梯度材料應用于算例1的薄壁管道中，在流速增大時第一階固有頻率的變化情況。由圖7曲線對比可見，隨n值的增大，失穩臨界流速顯著提高，保證了管道在較低流速下的穩定性。

圖7 不同梯度指數下第一階固有頻率隨流速變化曲線Fig.7 Curves of the first natural frequency with the increase of flowing velocities under different gradient indexes.

當管內流速為10 m/s時，在管道一端1/10處施加簡諧激勵力sin(2πt)，對該點的頻域響應進行傅里葉逆變換可得到時域響應，管道材料取不同n值時在0～1 s內的橫向位移響應曲線如圖8所示。隨著n值的增大，管道的振動幅度也有了顯著的降低。此外值得注意的是，n值由0增加至0.5，位移降低得非常顯著，而隨著n值的進一步增大，尤其在n值增加到2以上，位移的變化幅度逐漸變小。從功能梯度的結構來分析，圖9展示了取不同梯度指數時沿半徑方向變化的楊氏模量變化曲線，其他材料參數也有類似的變化規律。當n增大時，楊氏模量呈增大趨勢，從而管道的有效剛度也增大，因此位移振幅減小。管壁的平均楊氏模量值顯然是與曲線下覆蓋的面積呈正比的，當梯度值從0增加至0.5時，楊氏模量的增長量大，而隨著n值的進一步提升，曲線的斜率逐漸降低，材料參數的增長逐漸緩慢。

圖8 不同梯度指數橫向位移時域響應曲線Fig.8 Transverse displacement responses in time domain with different gradient indexes.

從提高管道的穩定性和降低振動幅度的角度來說，這種功能梯度材料輸流管道的梯度指數n越大越好，但是顯然實際應用還要考慮到材料的成本、制作工藝等因素，此外以這種陶瓷-鋁功能梯度材料為例，n值越大，陶瓷的體積率就越高，而雖然陶瓷材料的剛度更大，但其耐沖擊性能低，極易發生脆性斷裂。以本文的材料參數變化規律來說，當n>2時，其對于管道動態特性的影響并不大，因此對于本例的功能梯度輸流管道的梯度指數取為2是比較合適的。未來在工程中若要在管道上應用功能梯度材料，在選取梯度指數時應綜合考慮其對結構的影響、材料參數的變化規律以及成本和材料優缺點等實際因素。

圖9 不同梯度指數的楊氏模量變化曲線Fig.9 Curves of Young’s modulus with different gradient indexes.

4 結 論

本文應用離散化的思想，將功能梯度材料沿厚度方向劃分為若干均勻層，考慮流體與管壁間的相互作用，建立了功能梯度輸流管道的耦合振動模型，并應用動剛度法求解了管道振動的固有頻率及頻域、時域響應。當材料退化為普通均質材料時，計算結果與文獻結果完全一致，說明了本文所建立公式及計算方法的正確性。此外功能梯度管道的算例反映了功能梯度材料梯度指數的變化在提高結構固有頻率、延緩管道失穩方面的積極意義。本文的理論對于功能梯度材料和層合復合材料具有普適性，為新型復合材料在管道上的應用提供了理論依據。