曹銀行,柳貢民,張 龍
(1.哈爾濱工程大學(xué) 動力與能源工程學(xué)院,哈爾濱 150001;2.海裝駐哈爾濱地區(qū)第二軍事代表室,哈爾濱 150001)
遺傳算法是一種基于自然選擇和遺傳機制,在計算機上模擬生物進化的搜索尋優(yōu)算法,被廣泛應(yīng)用于各個領(lǐng)域[1-4],如在振動主動控制領(lǐng)域,文獻(xiàn)[5-7]將遺傳算法用于柔性結(jié)構(gòu)傳感器/作動器的位置優(yōu)化,取得了較好的效果,證明了遺傳算法的高效性和較好的全局收斂性。分支管路作為一個復(fù)雜的結(jié)構(gòu)系統(tǒng),分支數(shù)量較多,設(shè)計變量較多,適于采用遺傳算法這種直接以目標(biāo)函數(shù)值作為搜索信息,且可以同時搜索多個點優(yōu)化算法。目前已有文章將遺傳算法用于分支管路路徑尋優(yōu)問題,但未涉及分支管路動力學(xué)優(yōu)化設(shè)計。
本文首先借鑒吸收傳遞矩陣法[8-9,20,22]的基本思想,給出了考慮內(nèi)部流體時直管場傳遞矩陣和任意管路分支元件點傳遞矩陣的推導(dǎo)、應(yīng)用和實驗證明過程(見附錄);之后基于選定的分支管路設(shè)計基礎(chǔ)模型,參考均勻設(shè)計原則[21]選取計算樣本,利用吸收傳遞矩陣法對樣本的固有頻率、某一點的振動響應(yīng)加速度級等動力學(xué)指標(biāo)進行快速計算,最后利用二階多項式響應(yīng)面函數(shù)逼近回歸模型并應(yīng)用遺傳算法對分支管路系統(tǒng)進行優(yōu)化設(shè)計。
輸流管路動力學(xué)研究中常采用基于一維梁理論及流體平面波理論的四方程模型[10]、八方程模型[11]及十四方程模型[12-14,22]等,其中十四方程模型是目前較為完善的輸流管路動力學(xué)模型。
輸流直管簡化模型如圖1所示,同時管路動力學(xué)研究滿足下列基本假設(shè):
圖1 直管簡化模型Fig.1 Simplified model of straight pipeline
① 管路為等截面、純彈性、均質(zhì)、各向同性的圓型管;
② 流體壓力和流速在同一管道截面內(nèi)是恒值;
③ 管道與流體之間的摩擦耦合可以忽略[15];
④ 管內(nèi)不存在充液管道的空穴現(xiàn)象[16];
⑤ 用一維運動描述管內(nèi)流體低頻壓力波是可行的[17,18];
⑥ 忽略管路徑向變形引起的流體徑向運動及流體繞管軸的旋轉(zhuǎn)運動;
⑦ 管路的橫向彎曲振動符合Timoshenko梁模型的運動形式[19],即彎曲力和位移滿足胡克定律描述的線性關(guān)系。
本文選取參考文獻(xiàn)[22]中的十四方程模型展開研究,十四方程模型如下,包括四個軸向方程,八個橫向方程及兩個扭轉(zhuǎn)方程:
(1) 軸向方程
(2)y-z平面方程
(3)x-z平面方程
(4) 扭轉(zhuǎn)方程
其中:
m=ρpAp+ρfAf,B=ρpIp+ρfIf,
引入包含各向運動變量的時域狀態(tài)向量[14,22]
則直管十四方程模型可寫成矩陣方程的形式
(1)
(2)
式(2)兩邊同時左乘A-1
(3)
求解矩陣A-1B的特征值矩陣T(s)和特征向量矩陣V(s),可以推出:
V-1A-1BV=T
T(s)=diag{λ1(s),λ2(s),…,λ14(s)}
令v(z,s)=V-1(s)Y(z,s),式(3)兩端同時左乘矩陣V-1(s)
(4)
式(4)的形式解為
v(z,s)=E(z,s)v0(s)
(5)
其中,
v(z,s)=[v1(z,s),v2(z,s),…,v14(z,s)]T
v0(s)為未定的積分常數(shù),其值取決于邊界條件,而
由式(5)及E(0,s)=I14×14(14階單位陣)可得
Y(0,s)=V(s)E-1(L,s)V-1(s)Y(L,s)
(6)
即:
Φs=V(s)E-1(L,s)V-1(s)Φe
Φs,Φe分別為長度為L的直管始末端頻域狀態(tài)向量,Φs=Y(0,s),Φe=Y(L,s)。
因此可得直管場傳遞矩陣U
U=V(s)E-1(L,s)V-1(s)
(7)
結(jié)合傳遞矩陣法的基本思想以及管路邊界條件矩陣,管路振動整體方程可用矩陣的形式表示為[14]
(8)
式中:Ds、De均為7行14列的矩陣,分別表示管路始末端邊界條件,具體形式見參考文獻(xiàn)[20];Fs、Fe分別為管路始末端頻域下的外部激勵列向量。
式(8)可簡寫為
DtotΦtot=Ftot
(9)
式中,Dtot、Φtot、Ftot分別表示整體傳遞矩陣、整體狀態(tài)向量以及整體邊界激勵向量。
變換式(9),可以得到
由Φtot=[ΦsΦe]T可知,Φtot的前部分為直管始端狀態(tài)向量,結(jié)合直管場傳遞矩陣(7),則可得管路任意一點的狀態(tài)向量Φ(z,s)。
當(dāng)研究的問題是系統(tǒng)固有特性時,則可令整體外部激勵向量Ftot=0。式(9)簡化為
DtotΦtot=0
(10)
式中,Φtot具有非零解的條件是系數(shù)矩陣Φtot的行列式等于0,由此可以求得管路的固有頻率。
分支元件是管路系統(tǒng)中非常常見的管路元件,分支間的夾角會根據(jù)實際工程需求而復(fù)雜多樣。
圖2為包含N個子分支的分支元件,第n個子分支與第一個子分支間的夾角(順時針)用αn表示,第n個子分支在分支點處的狀態(tài)向量用Φn表示。取主傳遞路徑為通過子分支1傳至子分支2。
圖2 包含N個子分支的分支管路Fig.2 Branch pipeline containing N sub-branches
參考文獻(xiàn)[9,22],考慮內(nèi)含流體的分支元件,得到流體連續(xù)方程及分支點處的力與力矩的平衡方程
{Fz+AfP}1+…+{Fz+AfP}NcosαN-
{Fy}NsinαN=0
{Fy}1+…+{Fy}NcosαN+{Fz+AfP}NsinαN=0
{Mx}1-{Mx}2…-{Mx}N=0
{Fx}1-{Fx}2-…-{Fx}N=0
{My}1+{My}2cosα2+{Mz}2sinα2+…+
{My}NcosαN+{Mz}NsinαN=0
{Mz}1+{Mz}2cosα2-{My}2sinα2+…+
{Mz}NcosαN-{My}NsinαN=0
(11)
子分支1子與分支2在分支點處存在如下連續(xù)及平衡條件[22]
{P}1={P}2
(12)
同時在分支點處,第2個子分支與第n(3≤n≤N)個子分支存在連續(xù)和平衡條件[22]
{P}n={P}2
{Fx}1-{Fx}2-…-{Fx}N=0
(13)
經(jīng)過拉普拉斯變換后,由式(11)和式(12)容易得到
P1Φ1=P2Φ2+P3Φ3+…+PNΦN
(14)
P1,P2,…,PN為與α2,α3,…,αn有關(guān)的14×14的系數(shù)矩陣。
式(13)可寫成如下矩陣方程形式
AnΦn=A2Φ2
(15)
其中矩陣An為7×14的系數(shù)矩陣。
第n(3≤n≤N)個子分支(即除主路徑外的所有子分支)自分支點至邊界端的狀態(tài)向量傳遞關(guān)系為
Φn=UnΦn,b
其中Φn,b為第n個子分支邊界端的狀態(tài)向量,同時任意邊界激勵條件下在第n個子分支滿足
DnΦn,b=Fn
因此可得
(16)
由式(15)及(16)可得
HnΦn=TΦ2+Fnn
其中,
容易驗證,矩陣H2與T均為14階方陣,故有:
(17)
將式(17)代入式(14),有
Up+Mp
(18)
對于包含t個分支元件的管路系統(tǒng),在求解管路系統(tǒng)動力學(xué)問題時,只需將Dtot中的U,Ftot中的0分別進行替換
在利用吸收傳遞矩陣方法計算分支管路動力學(xué)特性時,分支元件和直管段管路可以獨立地分塊計算,并且分支數(shù)目的增加或減少并不影響分支管路系統(tǒng)的總體傳遞矩陣維數(shù),因此吸收傳遞矩陣方法有利于編制通用程序,且計算量相對穩(wěn)定,計算效率也較高。
結(jié)合2.2和3.1可得吸收傳遞矩陣法計算分支管路動力學(xué)特性流程圖,如圖3所示。
圖3 吸收傳遞矩陣法動力學(xué)計算流程Fig.3 Flow chart of dynamic calculation by absorbing transfer matrix method
4.1.1 遺傳算法的基本原理和結(jié)構(gòu)
遺傳算法是通過提取生物進化主要的內(nèi)在因素和外在因素來模擬生物進化過程的一種數(shù)學(xué)模型,它的構(gòu)造遵循達(dá)爾文的進化論和現(xiàn)代遺傳學(xué)原理。標(biāo)準(zhǔn)遺傳算法提供了一種求解復(fù)雜系統(tǒng)優(yōu)化問題的通用框架。對于一個需要優(yōu)化的實際應(yīng)用問題,一般可按圖4基本步驟構(gòu)造遺傳算法。
圖4 簡單遺傳算法流程圖Fig.4 Flow chart of simple genetic algorithm
4.1.2 分支管路遺傳算法動力學(xué)設(shè)計基本結(jié)構(gòu)
本文擬基于選定的分支管路設(shè)計基礎(chǔ)模型,應(yīng)用吸收傳遞矩陣數(shù)值計算方法和遺傳算法對分支管路進行優(yōu)化設(shè)計。首先根據(jù)優(yōu)化變量的參數(shù)變化空間,參考均勻原則[21]選擇計算樣本;然后應(yīng)用吸收傳遞矩陣法進行樣本計算,并利用二階多項式響應(yīng)面函數(shù)逼近回歸模型;最后應(yīng)用遺傳算法尋優(yōu),獲得動力學(xué)性能優(yōu)越的分支管路設(shè)計方案。分支管路遺傳算法優(yōu)化設(shè)計方案流程圖,如圖5所示。
圖5 分支管路動力學(xué)優(yōu)化設(shè)計流程圖Fig.5 Flow chart for dynamic optimization design of branch pipeline
以往的研究經(jīng)驗表明,由于輸流泵額定轉(zhuǎn)速、最高或最低轉(zhuǎn)速等條件的限制,正常工作的管路系統(tǒng)內(nèi)流體速度都有一定的范圍,且在該范圍內(nèi)流體速度對分支管路動力學(xué)計算結(jié)果的影響很小,影響更大的是輸送流體的密度、體積模量等屬性參數(shù),因此在以下的兩個算例中都暫時假設(shè)管路中流體的為靜止?fàn)顟B(tài)。
4.2.1 雙分支管路分支點位置的優(yōu)化選擇
圖6為以管路中心線描述的雙分支管路系統(tǒng)設(shè)計基礎(chǔ)模型,管路設(shè)計的目的為在主路徑AB段和BC段分別選取分支點F和G建立分支管段DF和EG,組成雙分支管路,同時保證該雙分支管路系統(tǒng)的第一階固有頻率最低。
圖6 雙分支管路設(shè)計基礎(chǔ)模型Fig.6 Basic model of double branch pipeline design
各直管段均滿足內(nèi)半徑R=0.025 m;壁厚e=0.004 m;管路材料的泊松比ν=0.29;密度ρ=7 800 kg/m3,楊氏模量E=168 GPa,管內(nèi)流體為空氣。
D、E兩點到管段AC的垂直距離均為1 m;D點到A點,E點到B點的水平距離也均為1 m。
① 優(yōu)化模型的建立(基礎(chǔ)模型)
優(yōu)化模型的形式如下:
MinF=F(x1,x2)
其中,F(xiàn)為優(yōu)化目標(biāo)即雙分支管路的第一階固有頻率,x1,x2為自變量,本設(shè)計實例里選定為AF和BG段的長度??紤]管口、分支點和管路半徑的影響,設(shè)優(yōu)化變量的變化空間,即x1,x2滿足:
0.1≤x1,x2≤1.9
② 二階響應(yīng)面模型的建立(回歸模型)
利用實際可用的回歸模型或關(guān)系表達(dá)式來替代詳細(xì)的模擬分析,可以大大降低模擬計算的次數(shù)。本文采用二階多項式響應(yīng)面模型作為回歸模型來搜索最優(yōu)設(shè)計方案。二階響應(yīng)面模型的基本方程如下
其中,N為變量數(shù),(x1,x2,…,xN)為模型的輸入變量集合,a0,bi,ci,cij(i=1,…,N;j=1,…,N)為方程的多項式系數(shù)。
對于雙分支管路設(shè)計問題,其變量個數(shù)N=2。將上式展開,則回歸模型即為含有6個系數(shù)的二階響應(yīng)面函數(shù)
(19)
參考均勻設(shè)計原則[21]生成了19個計算樣本見表1,利用吸收傳遞矩陣法對表1中的計算樣本進行計算,得到每個計算樣本的第一階固有頻率見表2。
表1 雙分支管路計算樣本Tab.1 Computational sample of double branch pipeline
表2 吸收傳遞矩陣法計算結(jié)果 Tab.2 Calculation results of absorption transfer matrix method Hz
利用Matlab軟件非線性擬合程序擬合表2中的計算結(jié)果得到二階響應(yīng)面模型的具體表達(dá)式如下
③ 遺傳算法尋優(yōu)
根據(jù)二階響應(yīng)面函數(shù)的特點,應(yīng)用Matlab遺傳算法程序,以雙分支管一階固有頻率F最小為優(yōu)化目標(biāo)進行尋優(yōu)。
④ 優(yōu)化設(shè)計結(jié)果
遺傳算法到的最優(yōu)解如表3所示,之后將最優(yōu)解代入吸收傳遞矩陣法進行計算,并建立有限元模型進行仿真,得到的雙分支管路系統(tǒng)一階固有頻率分別為20.0 Hz和20.1 Hz,理論優(yōu)化目標(biāo)與吸收傳遞矩陣法計算結(jié)果相差6%,與仿真結(jié)果差6.48%。這主要是因為在遺傳算法中沒有加入損失模型模塊,造成數(shù)值模擬的結(jié)果與理論優(yōu)化目標(biāo)有一定偏差。同時可以看出此時管路系統(tǒng)為對稱布置,與分支管路對稱布置時固有頻率更低的研究經(jīng)驗相符。
表3 遺傳算法優(yōu)化結(jié)果Tab.3 Optimization result of genetic algorithm
4.2.2 分支管路支撐位置的優(yōu)化設(shè)計
圖7為以管路中心線描述的平面Y型分支水管路設(shè)計基礎(chǔ)模型,管路設(shè)計的目的為在直管BC段和兩分支的EF及HI段分別選取一點施加管路彈性支撐,使得在A點受到大小為80 kN,作用時間為0.002 s的x向的瞬態(tài)激勵力時,G點在某一頻率下的振動響應(yīng)加速度級最小(算例里選定為50 Hz)。
圖7 支撐位置設(shè)計基礎(chǔ)模型Fig.7 Basic model for design of support positions
由參考文獻(xiàn)[22]可知彈性支撐點左、右兩側(cè)狀態(tài)向量的傳遞關(guān)系滿足ΦL=PΦR,其中彈性支撐點傳遞矩陣為
假設(shè)彈性支撐各向線性和扭轉(zhuǎn)剛度為
kx=1×106;ky=1×106;kz=1×106;
kx,ky,kz表示線性剛度,單位:N/m
kax=1×106;kay=1×106;kaz=1×106;
kax,kay,kaz表示扭轉(zhuǎn)剛度,單位:N/m
圖7中,各段長度滿足AB=CD=DE=DH=FG=IJ=0.5 m;BC=2 m;EF=HI=1 m;管內(nèi)充滿壓力為0.2 MPa的靜止的水;管路系統(tǒng)A、G兩端自由、J端固支,其它參數(shù)與4.2.1雙分支管路系統(tǒng)相同。
① 優(yōu)化模型的建立
優(yōu)化模型的形式如下
MinF=F(x1,x2,x3)
其中,F(xiàn)為優(yōu)化目標(biāo)即G點在50 Hz下的振動加速度級,x1,x2,x3為自變量,該設(shè)計實例里選定為支撐1到B點、支撐2到E點、支撐3到H點的距離。
② 二階響應(yīng)面模型的建立。
參考均勻設(shè)計原則[21]生成21個計算樣本。利用Matlab軟件非線性擬合程序擬合21個計算樣本的吸收傳遞矩陣法計算結(jié)果,得到二階響應(yīng)面函數(shù)模型如下
③ 優(yōu)化設(shè)計結(jié)果
應(yīng)用自編程遺傳算法優(yōu)化程序得到最優(yōu)解見表4。同時將最優(yōu)解代入吸收傳遞矩陣法整體計算模型進行計算,得到瞬態(tài)激勵下G點在50 Hz處80 kN作用力下的振動加速度級為145.2 dB,理論優(yōu)化目標(biāo)與吸收傳遞矩陣法計算結(jié)果相2%。(參考加速度為10-6m/s2)。此時管路系統(tǒng)同樣為對稱布置。
表4 遺傳算法優(yōu)化結(jié)果Tab.4 Optimization result of genetic algorithm
本文借鑒吸收傳遞矩陣的基本思想給出了任意管路分支元件點傳遞矩陣的建立過程,實現(xiàn)了任意分支管路系統(tǒng)動力學(xué)問題的求解并通過與實驗結(jié)果的對比證明其正確性(見附錄);之后基于選定的分支管路設(shè)計基礎(chǔ)模型,利用均勻設(shè)計原則選取計算樣本并應(yīng)用傳遞矩陣法對樣本模型進行計算,最后利用二階多項式響應(yīng)面函數(shù)逼近回歸模型并應(yīng)用遺傳算法尋優(yōu),對分支管路進行優(yōu)化設(shè)計,通過吸收傳遞矩陣法、均勻設(shè)計原則和數(shù)學(xué)建模優(yōu)化的方法,更好地節(jié)省了計算成本和資源,對分支管路設(shè)計過程有一定的指導(dǎo)意義。
綜上所述,本文的特征在于:
(1) 給出了基于吸收傳遞矩陣法的分支管路系統(tǒng)動力學(xué)計算程式化方法,提高了待優(yōu)化分支管路系統(tǒng)樣本計算的效率;
(2) 參考均勻設(shè)計原則設(shè)置計算樣本,能在減少實驗次數(shù)的同時更好地反映各自變量對優(yōu)化目標(biāo)的影響;
(3) 選擇精度較高的二階響應(yīng)面模型作為回歸模型,直接代替數(shù)值計算進行優(yōu)化設(shè)計,大幅度地提高了優(yōu)化效率;
(4) 采用遺傳算法能很好地尋求到模型的最優(yōu)解。
附錄A 分支管路吸收傳遞矩陣法實驗證明
實驗管路安裝示意圖如圖A1所示。
圖A1 管路安裝示意圖Fig.A1 Pipeline installation diagram
吸收傳遞矩陣法計算結(jié)果與實驗測得 1、2 兩點的管壁振速傳函的比較結(jié)果分別如圖A2和圖A3所示。圖A2和圖A3的比較結(jié)果表明吸收傳遞矩陣法理論預(yù)測結(jié)果與實驗結(jié)果吻合較好,進而證明了其正確性。
圖A2 1點x方向的傳函Fig.A2 Transfer function of point 1 in x direction
圖A3 2點x方向的傳函Fig.A3 Transfer function of point 2 in x direction