非線性輸流曲管面內(nèi)振動的樣條小波有限元方法研究

曹建華， 黃穎青， 仝國軍， 趙年順

(1. 黃山學(xué)院 機(jī)電工程學(xué)院，安徽 黃山 245021; 2. 中國科學(xué)技術(shù)大學(xué) 工程科學(xué)學(xué)院 近代力學(xué)系，合肥 230027; 3. 河北水利電力學(xué)院 巖土工程安全與變形控制重點(diǎn)實(shí)驗(yàn)室，河北 滄州 061001)

由于小波具有多尺度和多分辨的優(yōu)越特性，基于小波的數(shù)值方法，尤其是小波有限元方法，在工程和科學(xué)計(jì)算的應(yīng)用研究吸引了眾多國內(nèi)外學(xué)者，并取得了巨大的進(jìn)步。

在小波有限元研究方面，Williams等[1-2]做了很多基礎(chǔ)性和開創(chuàng)性的研究工作，并研究第二代小波多分辨在有限元中的應(yīng)用。國內(nèi)眾多學(xué)者在小波有限元方面做了很多工作。Chen等[3]應(yīng)用二維多尺度區(qū)間樣條小波單元于自適應(yīng)有限元分析中，并用算例驗(yàn)證了其在奇異性問題上的有效性和可靠性。Wang等[4]又采用第二代小波構(gòu)建小波有限元，并應(yīng)用于偏微分的求解算例上。Xiang等[5]采用區(qū)間樣條小波有限元，求解了轉(zhuǎn)子動力學(xué)中轉(zhuǎn)子軸承系統(tǒng)的頻率問題。Xiang等[6]也用同樣的方法，推導(dǎo)了圓錐厚殼的有限元矩陣，計(jì)算其模態(tài)數(shù)據(jù)，并用小波分解和支持向量機(jī)進(jìn)行了探傷分析。

小波有限元與其他傳統(tǒng)方法、新技術(shù)的結(jié)合，能夠提高小波有限元方法的精度，縮短其計(jì)算時(shí)間。Shen等[7]結(jié)合小波有限單元和基于FFT(fast Fourier transform)的譜分析法，構(gòu)建了高精度動剛度矩陣，研究了有裂紋或者脫落的桿和梁的波動特性。Zhang等[8]將小波有限單元和拉普拉斯變換相結(jié)合，提出了一種既能減少單元數(shù)又不用減小時(shí)間間隔的新方法，并將其應(yīng)用于一維桿狀結(jié)構(gòu)中的超聲波模擬。Hao等[9]采用區(qū)間樣條小波構(gòu)建復(fù)合板單元，基于圖形處理單元(graphics processing unit,GPU)，實(shí)現(xiàn)一種小波板單元并行計(jì)算技術(shù)，應(yīng)用于復(fù)合材料層合板的結(jié)構(gòu)健康監(jiān)測系統(tǒng)中，并指出基于GPU的并行計(jì)算比基于CPU的計(jì)算快了140多倍。

由以上可知，小波有限元計(jì)算可靠，應(yīng)用廣泛，但文獻(xiàn)大多集中在線性計(jì)算方面，本文將以非線性輸流曲管為研究對象，應(yīng)用小波有限元求解其振動特性。輸流管道在工業(yè)中應(yīng)用廣泛，受到學(xué)者們的關(guān)注。Misra等[10-11]采用傳統(tǒng)有限元分別求解了基于軸線不可伸縮理論和軸線可伸縮理論的輸流曲管振動問題。Jung等[12]采用非線性應(yīng)變推導(dǎo)了軸線可伸長的輸流曲管面內(nèi)振動微分方程，并用傳統(tǒng)有限元求解了系統(tǒng)頻率。李寶輝等[13]采用波動法求解輸流曲管的面內(nèi)振動?？紤]穩(wěn)態(tài)組合力，Dai等[14]應(yīng)用傳遞矩陣法分析了空間輸流管道的流致振動。uczko等[15-16]采用伽遼金法分析了各類螺旋的三維輸流曲管流致振動問題，也分析了流速和曲率對平面輸流曲管的振動模態(tài)和頻率的影響。Qu等[17]結(jié)合有限元和隨機(jī)振動離散法，研究了受隨機(jī)激勵(lì)的液壓曲管的振動特性。Zhou等[18]應(yīng)用有限元研究了初始幾何缺陷對微曲懸臂輸流管道靜平衡的影響，且提出其臨界流速大小取決于靜平衡構(gòu)形。Ye等[19]采用伽遼金法、微分積分單元法和離散傅里葉法分析研究了微曲輸流管道的超臨界動力特性?；诔叨葹?、階數(shù)為6的區(qū)間樣條小波，文獻(xiàn)[20]構(gòu)建了一維小波曲管單元，計(jì)算了軸線不可伸長的輸流曲管頻率。

本文詳述基于曲管軸線可伸長的非線性輸流曲管面內(nèi)振動微分方程的推導(dǎo)過程，構(gòu)建小波曲管單元，并應(yīng)用小波有限元離散其微分方程，分析輸流曲管的頻率、模態(tài)和動力時(shí)間響應(yīng)。

1 非線性輸流曲管面內(nèi)振動微分方程的推導(dǎo)

1.1 管道橫截面上一點(diǎn)的應(yīng)變公式推導(dǎo)

如圖1所示，截取一段軸線弧長為ds、曲率半徑為R的ds的曲管微元，在截面上放置一個(gè)局部坐標(biāo)系，坐標(biāo)軸為x和y軸，單位矢量分別為i和j，曲管微元離軸線為x處的弧長為L0。變形后，軸線弧長增大εEds，曲率半徑為r，εE為軸線應(yīng)變，離軸線為x處曲管微元的弧長為L1。則曲管微元離軸線為x處一點(diǎn)的應(yīng)變εM[21]為

圖1 曲管微元變形前和變形后示意圖Fig.1 The diagram of the undeformed element and the deformed element of curved pipe

(1)

對于細(xì)長輸流曲管管，x=R,x=r，因此，由式(1)可得管道橫截面上一點(diǎn)的應(yīng)變?yōu)?/p>

εM=εE+xΔκ

(2)

其中

(3)

圖2 曲管軸線變形前和變形后示意圖Fig.2 The diagram of the undeformed and deformed axial line of curved pipe

變形后的矢徑r′表示為

r′=r+u(s)n+v(s)t

(4)

式中，u(s)，v(s)分別為軸線上一點(diǎn)沿切線單位矢量t和負(fù)法線單位矢量n方向的位移量，由圖2可知，軸線的拉伸應(yīng)變εE為

(5)

根據(jù)Hutchinson的課件(https://scholar.harv ard.edu/files/jiaweiyang/files/ curved_beams.pdf)，可以推出曲率的變化率Δκ為

(6)

根據(jù)拉格朗日應(yīng)變公式和式(5)，拉伸應(yīng)變ηE為

(7)

對于小應(yīng)變、適度轉(zhuǎn)動的曲管變形，εE?1，可得

(8)

則軸線應(yīng)變εE為

(9)

因此，由式(3)可知管道橫截面上一點(diǎn)的應(yīng)變公式εM為

即為[22]

(10)

式(10)與Jung等所引用的公式是一致的。

1.2 輸流曲管面內(nèi)流固耦合振動微分方程

本節(jié)基于式(10)，簡述Jung等推導(dǎo)輸流曲管面內(nèi)流固耦合振動微分方程的過程。如圖3所示兩端固定的細(xì)長輸流曲管，等截面且軸線為半圓，其軸線半徑為R。假設(shè)流經(jīng)管道的流體流速為U的柱塞流體，OXY坐標(biāo)系為固定在機(jī)架上的慣性坐標(biāo)系，θ為曲管截面與X軸正向的夾角。在截面上放置一個(gè)局部坐標(biāo)系，坐標(biāo)軸為x和y軸，單位矢量分別為i和j，則輸流曲管面內(nèi)徑向位移ur(x,θ,t) 和周向位移uθ(x,θ,t)表達(dá)式分別為

圖3 軸線為半圓的細(xì)長等截面輸流曲管Fig.3 The geometry of semi-circular curved pipe conveying fluid

ur(x,θ,t)=u(θ,t)uθ(x,θ,t)=v(θ,t)+xφ

(11)

式中：t為時(shí)間；u與v分別為管道軸線上的點(diǎn)沿徑向和周向位移；φ為管道振動時(shí)截面的面內(nèi)轉(zhuǎn)角，其表達(dá)式為式(8)。

管道中一點(diǎn)位移可以表示為

rp=(R+u+x)i+(v+xφ)j

(12)

管道一點(diǎn)的速度為

(13)

其管道中流體的速度為

(14)

其中t如圖3所示，含義與圖2一致。輸流曲管的動能為

(15)

式中，mp和mf分別為單位長度管道和流體的質(zhì)量。輸流曲管其勢能為

(16)

式中，εM采用非線性von Karman應(yīng)變，其表達(dá)式如式(10)所示。采用線性應(yīng)力表達(dá)式為

(17)

對于輸流曲管，哈密頓原理表述為如下

(18)

其中，非保守力做功為

(19)

經(jīng)過一系列變分，可得

(20)

(21)

式中，Q為軸力，其表達(dá)式為

(22)

兩端固定的細(xì)長輸流管道的邊界條件為

θ=0和π,u=u′=v=0

(23)

2 樣條小波有限元離散輸流曲管面內(nèi)振動微分方程

2.1 區(qū)間樣條小波尺度函數(shù)

(24)

式中：j為尺度；m為0和1的多重節(jié)點(diǎn)數(shù)。基于上述序列點(diǎn)，定義在ti點(diǎn)處j尺度，m階數(shù)的B-樣條函數(shù)表達(dá)式為

(25)

(26)

因此, 在區(qū)間[0,1]上樣條小波的尺度函數(shù)可寫成向量形式

(27)

式中：ξ∈[0,1]； 且2j0≥2m-1，j≥j0。

在小波數(shù)值計(jì)算中，進(jìn)行積分生成各類矩陣，主要有兩類方法：一是直接用數(shù)學(xué)軟件中函數(shù)生成樣條小波尺度函數(shù)，應(yīng)用積分函數(shù)生成單元矩陣，存儲起來，然后在程序中直接調(diào)用即可； 二是本文采用的方法，編程生成樣條函數(shù)，再構(gòu)建樣條小波尺度函數(shù)，利用高斯積分生成單元矩陣。在線性分析中矩陣不需要更新，兩種方法均能勝任，但在非線性分析中，迭代過程中需要更新單元矩陣，第一種方法需要不停地調(diào)用數(shù)學(xué)軟件的積分函數(shù)，花費(fèi)太多時(shí)間，而第二種方法是采用高斯方法積分，故能夠適用非線性計(jì)算。

2.2 樣條小波有限元離散輸流曲管面內(nèi)振動微分方程

本文采用樣條小波有限元離散輸流曲管的非線性微分方程。在樣條小波有限元中，本文采用尺度j1為3、階數(shù)m1為4區(qū)間樣條小條尺度函數(shù)(BSWI43)作為橫向位移插值函數(shù)，采用尺度j2為2、階數(shù)m2為3區(qū)間樣條小條尺度函數(shù)(BSWI32)作為軸向位移的插值函數(shù)，兩種區(qū)間樣條小條尺度函數(shù)分別如圖4所示。

圖4 區(qū)間樣條小波尺度函數(shù)Fig.4 The scaling function of BSWI43 and BSWI32 in the interval [0, 1]

輸流曲管的軸向位移、橫向位移表達(dá)式分別為

(28)

定義單元長度為le，單元徑向位移和單元周向位移分別為

(29)

其中，ξi=(i-1)/2j,i=1,2，…,n+1=2j+1。

將式(28)代入式(29)，可以得到

(30)

其中

將式(30)代入式(28)

(31)

式中，Nu(ξ)=Φu(Re,u)-1和Nv(ξ)=Φv(Re,v)-1分別為徑向位移和周向位移的樣條小波有限元的形函數(shù)。令

(32)

式中， “′”和“″”分別為對ξ的一階和二階導(dǎo)數(shù)。根據(jù)傳統(tǒng)有限元的離散過程，單元離散方程表示為如式(33)所示

(33)

其中

(34)

(35)

(36)

(37)

式中，“·”和“··”分別為對時(shí)間的一階和二階導(dǎo)數(shù)。將所有單元集合起來，可得全局離散常微分方程如下

(38)

式中，q為所有單元ue，ve位移的集合；Mg，Cg，Kg(q)分別為管道的全局質(zhì)量矩陣、全局阻尼矩陣和全局剛度矩陣。由于Kg(q)為q的函數(shù)，因此是一個(gè)非線性問題。

3 非線性問題的求解方法

3.1 曲管的自然頻率

為了求得靜平衡位置處曲管的自然頻率，需要將式(38)線性化，即為

(39)

式中，KTg為在靜平衡位置處的切線剛度矩陣。

(40)

式中，q0為靜平衡時(shí)的曲管變形，由式(41)求出

Kg(q0)q0=F

(41)

計(jì)算過程都是在施加邊界條件之后進(jìn)行的。求解式(41)非線性問題，主要采用直接迭代法和牛頓-拉斐森迭代法[24]，讀者可參考相關(guān)文獻(xiàn)，在此不再贅述。考慮q=QeΩt，式(39)改寫為

(42)

計(jì)算式(42)的特征值，即可求出自然頻率，定義無量綱頻率和流速分別為

(43)

3.2 切線剛度矩陣的求解

計(jì)算切線剛度矩陣KTg的數(shù)值方法有前向差分法、中心差分法、復(fù)數(shù)步進(jìn)法，差分法對步長要求有一定的要求，太小并不一定能提高精度，而復(fù)數(shù)步進(jìn)法卻沒有這方面的限制[25-26]。復(fù)數(shù)步進(jìn)法的原理是：將復(fù)變函數(shù)f(x+ih)在x處進(jìn)行泰勒展開，f(x+ih)的表達(dá)式如下

(44)

式中：x+ih中的i為虛數(shù)單位；h為沿虛軸的小擾動步長。忽略高階項(xiàng)，取其虛部系數(shù)，可得f(x)的導(dǎo)數(shù)為

(45)

根據(jù)上述原理，求解切線剛度矩陣的復(fù)數(shù)步長法表達(dá)式如下

(46)

式中：q+iheJ中的i為虛數(shù)單位； Im[·]為復(fù)數(shù)的虛部系數(shù)；Fint,I(U+iheJ)為以復(fù)數(shù)向量為變量的函數(shù)向量Fint的第I個(gè)分量。 由于沒有差分法中的加減，精度有所提高。

4 數(shù)值算例

根據(jù)第3章的分析，作者編寫了一套用于輸流管道流固耦合振動分析的小波有限元程序，其中包含：生成區(qū)間樣條小波尺度函數(shù)的基礎(chǔ)程序、用于計(jì)算樣條小波有限元單元矩陣的分段高斯積分程序和相關(guān)有限元的程序，應(yīng)用于本文中輸流曲管流固耦合面內(nèi)振動分析。本章將小波有限元計(jì)算的曲管頻率、模態(tài)和時(shí)間積分與文獻(xiàn)、傳統(tǒng)有限元進(jìn)行對比，驗(yàn)證精度。

表1 當(dāng)時(shí)輸流曲管面內(nèi)無量綱頻率對比結(jié)果Tab.1 The comparison of the dimensionless natural frequencies of fluid-conveying curved pipe with

圖5 輸流曲管前3階徑向和周向模態(tài)Fig.5 The first three modes of curved pipe conveying fluid

圖6 不同單元數(shù)計(jì)算下的輸流曲管前3階無量綱頻率實(shí)部隨流速變化曲線Fig.6 The first three dimensionless frequencies of a fluid-conveying semi-circular pipe as a function of flow velocity solved by different numbers of wavelet-based elements

圖7 輸流曲管前3階無量綱頻率實(shí)部隨流速變化曲線比較Fig.7 The comparison of the first three dimensionless frequencies of a fluid-conveying semi-circular pipe as a function of flow velocity

以上數(shù)值結(jié)果是采用復(fù)數(shù)步進(jìn)法生成切線剛度矩陣的，分別應(yīng)用復(fù)數(shù)步進(jìn)法、中心差分法、前向差分法和向后差分法計(jì)算輸流曲管頻率隨流速變化曲線，如圖8所示，曲線相差不大，中心差分法、前向差分法和向后差分法的計(jì)算結(jié)果保持一致。由局部放大圖可知，復(fù)數(shù)步進(jìn)法的計(jì)算結(jié)果稍微小一些。

圖8 不同計(jì)算切線剛度矩陣方法的輸流曲管前3階頻率 實(shí)部隨流速變化曲線比較Fig.8 The comparison of the first three dimensionless frequencies of a fluid-conveying semi-circular pipe as a function of flow velocity solved by different method of computing the tangent stiffness

以上計(jì)算均是關(guān)于輸流曲管的靜平衡位置的計(jì)算結(jié)果，接著計(jì)算不同流速對于曲管靜平衡位置的影響，如圖9所示。由圖9可知，隨著流速增大，輸流曲管靜平衡時(shí)徑向位移u和周向位移v的絕對值均隨之增大，這與式(41)相符：流度越大，力F越大，位移也就越大。

圖9 不同流速下輸流曲管的靜平衡位置Fig.9 The effects of different flow velocities on the static equilibrium of a fluid-conveying semi-circular pipe

圖10 非線性輸流曲管中點(diǎn)徑向位移的時(shí)域動力 響應(yīng)Fig.10 The time history of the midpoint’s radial displacement of curved pipe conveying fluid

圖11 不同流速下非線性輸流曲管中點(diǎn)徑向位移的 時(shí)域動力響應(yīng)Fig.11 The time history of the midpoint’s radial displacement of curved pipe conveying fluid with different flow velocities

圖12 非線性輸流曲管中點(diǎn)徑向位移的時(shí)域 動力響應(yīng)Fig.12 The time history of the midpoint’s radial displacement of curved pipe conveying fluid

表2 不同流速下時(shí)域動力響應(yīng)計(jì)算時(shí)間比較Tab.2 The comparison of the computational time of the time response solved by wavelet-based finite element method and traditional finite element method

5 結(jié) 論

本文在前人文獻(xiàn)的基礎(chǔ)上，詳細(xì)推導(dǎo)曲管截面上一點(diǎn)的應(yīng)變公式，并簡述了輸流曲管面內(nèi)振動非線性微分方程的推導(dǎo)過程，采用尺度為3、階數(shù)為4的區(qū)間樣條小條尺度函數(shù)作為徑向位移的插值函數(shù)，尺度為2、階數(shù)為3的區(qū)間樣條小條尺度函數(shù)作為周向位移的插值函數(shù)，推導(dǎo)了區(qū)間樣條小波曲管單元相關(guān)矩陣，并將其用于離散輸流曲管面內(nèi)振動微分方程，求解了輸流曲管的頻率、模態(tài)和時(shí)域響應(yīng)，并與文獻(xiàn)、傳統(tǒng)有限元對比，進(jìn)行驗(yàn)證。根據(jù)數(shù)值結(jié)果討論，主要結(jié)論如下：

(1) 中心差分法、前向差分法和向后差分法計(jì)算切線剛度矩陣的頻率結(jié)果保持一致，與復(fù)數(shù)步進(jìn)法相比，數(shù)值相差不大，說明了復(fù)數(shù)步進(jìn)法的精確性和可靠性。

(3) 與傳統(tǒng)有限元對比，采用小波有限元計(jì)算的時(shí)域動力響應(yīng)相差不大，且數(shù)值結(jié)果表明流速越大，位移時(shí)域響應(yīng)就越大。

(4) 無論是頻率計(jì)算還是時(shí)域動力響應(yīng)計(jì)算，小波有限元計(jì)算時(shí)間都比傳統(tǒng)有限元的少，收斂速度快。

綜上所述，與傳統(tǒng)有限元相比，小波有限元在非線性輸流曲管面內(nèi)振動計(jì)算結(jié)果準(zhǔn)確，計(jì)算時(shí)間少，有著一定的優(yōu)勢。