基于進(jìn)化規(guī)劃算法的導(dǎo)波頻散曲線計(jì)算

,,

(1. 92330部隊(duì), 青島 266103;2.海軍駐武漢438廠軍代室, 武漢 4300613.海軍工程大學(xué) 動(dòng)力工程學(xué)院, 武漢 430033)

作為一種新興的無損檢測(cè)技術(shù)，基于磁致伸縮效應(yīng)的導(dǎo)波無損檢測(cè)方法可以實(shí)現(xiàn)非接觸、大范圍及長(zhǎng)距離的快速檢測(cè)，是近年來國內(nèi)外無損檢測(cè)界研究的熱點(diǎn)技術(shù)之一[1-2]。該技術(shù)的理論基礎(chǔ)是導(dǎo)波傳播理論及導(dǎo)波特性，其中導(dǎo)波最重要的兩個(gè)特性就是頻散和多模態(tài)特性：頻散特性即波的傳播速度隨著頻率的變化而變化，多模態(tài)特性就是一定頻率下可以存在多個(gè)模態(tài)的導(dǎo)波，而最直觀反映導(dǎo)波這兩種特征的就是頻散曲線，其是研究分析導(dǎo)波特性和利用導(dǎo)波實(shí)施無損檢測(cè)時(shí)選擇激勵(lì)頻率的重要工具。由于不考慮導(dǎo)波能量的泄漏，相對(duì)而言，求自由管道頻散方程的實(shí)數(shù)解比較簡(jiǎn)單，并且計(jì)算方程實(shí)數(shù)根的算法很多，如二分法、牛頓迭代法、弦線法等。相對(duì)于自由系統(tǒng)中導(dǎo)波頻散方程的求解，泄漏系統(tǒng)中的方程求解要困難得多，而導(dǎo)波在泄漏系統(tǒng)中的衰減特性通常使用復(fù)頻率或復(fù)波數(shù)來描述[3-6]。ARISTéGUI等[7]討論了復(fù)頻率和復(fù)波數(shù)的差別，發(fā)現(xiàn)復(fù)波數(shù)能更準(zhǔn)確地描述導(dǎo)波的衰減。筆者提出了一種基于進(jìn)化規(guī)劃的加載流體管道中頻散方程復(fù)波數(shù)求解的方法，充分利用了進(jìn)化規(guī)劃智能搜索的特點(diǎn)，大大簡(jiǎn)化了粗略和精確搜索反復(fù)求解的過程。

1 加載流體管道中導(dǎo)波頻散方程的建立

假設(shè)管道材料特性是均勻、各向同性的線彈性體，管道為軸對(duì)稱且無限長(zhǎng)。對(duì)于各向同性的彈性固體介質(zhì)，若不考慮體力的影響，其一般的彈性動(dòng)力學(xué)柱坐標(biāo)下的運(yùn)動(dòng)方程(Navier-Stokes方程)為[8]

(λ+μ)(

(1)

式中：u為位移矢量;ρ為材料密度;λ和μ為L(zhǎng)amé常數(shù)。

時(shí)間諧振位移矢量u通過Helmholtz分解表示為可壓縮標(biāo)量勢(shì)φ的梯度和零散度矢量φ的旋度，即

(2)

將式(2)代入Navier 運(yùn)動(dòng)方程得

(λ+μ)·(φ+×φ)+

μ2(φ+×φ)=

(3)

文章以管道中充滿水為例進(jìn)行分析，由于非黏性流體不支持剪應(yīng)力，因而固-液和固-氣交界處的邊界條件不同于固-固交界面的邊界條件。在管道內(nèi)壁與流體的交界處，只有徑向位移ur，正應(yīng)力σrr及壓縮應(yīng)力σrz連續(xù)，而管道外表面上應(yīng)力自由，則邊界條件為

σrr=0,σrz=0，r=b

(4)

將管道中導(dǎo)波位移和應(yīng)力表達(dá)式以及水中導(dǎo)波的位移和應(yīng)力表達(dá)式代入邊界條件，同真空中管道一樣，產(chǎn)生一組特征方程，為以幅度A、B、C1、D1、E的矩陣形式

[Eij]·[H]=0

(5)

式中：H=[ABC1D1E]T;Eij為系數(shù)矩陣。

要使上述方程有非零解，其系數(shù)矩陣行列式必須為零，即

|Eij|=0，i,j=1,2,…,5

(6)

式(6)即為加載流體管道中的導(dǎo)波頻散方程。

2 導(dǎo)波頻散曲線計(jì)算

2.1 進(jìn)化規(guī)劃及變異算子

進(jìn)化規(guī)劃(Evolutionary Programming,EP)是三種典型進(jìn)化算法(Evolutionary Algorithms,EA)之一，最早由美國FOGEL L J等[9]于20世紀(jì)60年代提出，后經(jīng)FOGEL D B等[10]完善而成，其主要的應(yīng)用介于數(shù)值分析和人工智能之間。同其他進(jìn)化算法比較，EP主要有以下特點(diǎn)：① 對(duì)生物進(jìn)化過程的模擬主要著眼于物種的進(jìn)化過程，沒有交叉算子，變異是唯一的操作方法；② 常用q選擇運(yùn)算，著重于群體中個(gè)體的競(jìng)爭(zhēng)選擇；③ 不必對(duì)變量進(jìn)行編碼，直接利用實(shí)數(shù)編碼，以問題的可行解作為個(gè)體的表現(xiàn)形式，表述問題更加自然，搜索過程比較平穩(wěn)，便于應(yīng)用；④ 以n維實(shí)數(shù)空間上的優(yōu)化問題為主要處理對(duì)象。

進(jìn)化規(guī)劃采用十進(jìn)制編碼，沒有重組或交換算子，但有選擇，進(jìn)化主要依賴突變，操作控制簡(jiǎn)單。

在標(biāo)準(zhǔn)進(jìn)化規(guī)劃中，個(gè)體的表達(dá)形式為

i=xi+ζ·Fi(0,1)

(7)

式中：xi為父代個(gè)體；i為子代個(gè)體，由父代個(gè)體產(chǎn)生；ζ為一個(gè)變異尺度；Fi(0,1)為服從某一分布的隨機(jī)數(shù)，也就是下面所討論的變異算子。

式(7)表明，新個(gè)體是在舊個(gè)體的基礎(chǔ)上添加一個(gè)隨機(jī)數(shù)形成的，添加的隨機(jī)數(shù)的數(shù)值與個(gè)體的適應(yīng)度有關(guān)，適應(yīng)度大的個(gè)體添加值也大，反之亦然。

根據(jù)表達(dá)方式，進(jìn)化規(guī)劃首先產(chǎn)生μ個(gè)初始個(gè)體，對(duì)其添加變異。然后從μ個(gè)舊個(gè)體及μ個(gè)新個(gè)體，共2μ個(gè)個(gè)體中根據(jù)適應(yīng)度挑選出μ個(gè)個(gè)體組成新群體。如此反復(fù)迭代，直至得到滿意結(jié)果。進(jìn)化規(guī)劃的工作流程類似于其他進(jìn)化算法的流程，同樣經(jīng)歷“產(chǎn)生初始群體-突變-計(jì)算個(gè)體適應(yīng)度-選擇-組成新群體”的過程，然后反復(fù)迭代，一代一代地進(jìn)化，直至達(dá)到最優(yōu)解。

通過進(jìn)化算法的個(gè)體表現(xiàn)形式，即從式(7)可以看出，變異算子是該算法中主要的遺傳算子，對(duì)算法本身的效率有著重要的影響。高斯算子的局部開發(fā)能力較好，柯西算子的全局探索能力較強(qiáng)，兩類算子各有優(yōu)點(diǎn)和不足，同時(shí)兩者的變異性能具有一定的互補(bǔ)性。t分布的特性使得t算子能整合這兩類算子的長(zhǎng)處，從而更便利地找到最優(yōu)解[11]。鑒于t算子的優(yōu)良特性，文章將采用其作為進(jìn)化規(guī)劃算法的變異算子。

2.2 計(jì)算步驟

基于進(jìn)化規(guī)劃算法的頻散方程求解的流程如圖1所示，可具體描述為：

圖1 利用進(jìn)化規(guī)劃算法計(jì)算頻散曲線的流程圖

(1) 參數(shù)初始化。輸入管道內(nèi)徑ra、外徑rb，材料密度ρ，彈性模量E,泊松比σ，流體密度ρf，流體中縱波速度cf；

(2) 建立頻散方程。根據(jù)邊界條件建立流體管道的頻散方程|Eij|=0，設(shè)定初始計(jì)算頻率f，頻率步長(zhǎng)fstep，計(jì)算頻率最大值fmax；

圖2 頻散方程行列式與實(shí)波數(shù)之間的變化關(guān)系

(4) 基于進(jìn)化規(guī)劃的頻散方程求解。選擇合適的t算子，在各個(gè)模態(tài)波數(shù)的定義域內(nèi)求解計(jì)算頻率下對(duì)應(yīng)模態(tài)頻散方程的復(fù)波數(shù)解k；

(5) 計(jì)算相速度cp。求出波數(shù)的虛部kimag為衰減系數(shù)catten，由實(shí)部kreal計(jì)算各個(gè)模態(tài)的相速度cp=2πf/kreal；

(6) 若f>fmax，轉(zhuǎn)入步驟(7)，否則f=f+fstep，轉(zhuǎn)入步驟(3)；

(7) 計(jì)算群速度cg。根據(jù)cg=dω/dk計(jì)算各模態(tài)的群速度。

2.3 t算子的選擇

如上文所述，進(jìn)化規(guī)劃算法中主要通過變異控制整個(gè)進(jìn)化過程，故變異算子的選擇會(huì)直接影響算法性能的優(yōu)劣。一般來說，當(dāng)種群中目標(biāo)函數(shù)值相差較小的個(gè)體可能位于最優(yōu)點(diǎn)附近時(shí)，算法應(yīng)該以局部開發(fā)為主；與目標(biāo)函數(shù)值相差比較大的個(gè)體應(yīng)加大變異尺度，在更廣的空間范圍內(nèi)進(jìn)行全局搜索，盡快逃離不利區(qū)域進(jìn)入目標(biāo)區(qū)域，算法的這種功能就只有通過合理控制進(jìn)化變異算子來實(shí)現(xiàn)。t算子變異尺度由變異步長(zhǎng)ζ和自由度n控制，文章將采用保持自由度n，變異步長(zhǎng)ζ線性遞減的t算子控制變異尺度優(yōu)化進(jìn)化過程。圖3所示為使用不同變異步長(zhǎng)算子時(shí)，基于tEP求解導(dǎo)波頻散方程的進(jìn)化過程。當(dāng)變異步長(zhǎng)ζ=0.1時(shí)，在經(jīng)過很長(zhǎng)一段時(shí)間才開始有向最優(yōu)點(diǎn)進(jìn)化的趨勢(shì)，進(jìn)化速率有點(diǎn)緩慢。ζ=2時(shí)，開始時(shí)間段內(nèi)目標(biāo)函數(shù)迅速下降，但經(jīng)過二三十代后趨于平穩(wěn)，此時(shí)收斂精度不高。采用變異步長(zhǎng)ζ線性遞減的t算子時(shí)，開始階段ζ相對(duì)較大，可以看出目標(biāo)函數(shù)值迅速收斂進(jìn)入目標(biāo)區(qū)域，ζ隨著進(jìn)化代數(shù)的增加也逐漸變小，算法開始以局部搜索為主來尋找全局最優(yōu)值。變異步長(zhǎng)ζ線性遞減的t算子融合了大步長(zhǎng)和小步長(zhǎng)的優(yōu)點(diǎn)，故采用這種變步長(zhǎng)的t算子來求解頻散方程。

圖3 不同算子時(shí)基于t EP的頻散方程求解進(jìn)化過程

2.4 試驗(yàn)驗(yàn)證

使用自主開發(fā)的磁致伸縮導(dǎo)波檢測(cè)系統(tǒng)(GWNDT-II)對(duì)某充水管道進(jìn)行檢測(cè)。管道的相關(guān)參數(shù)為：內(nèi)半徑，27 mm;壁厚,4 mm;密度,7 932 kg·m-3;泊松比σ,0.29;彈性模量E,210 GPa;管長(zhǎng)6 m。水的屬性參數(shù)為:密度,1 000 kg·m-3;縱波速度cf,1 500 m·s-1。

利用文章提出的算法計(jì)算得到該充水管道的導(dǎo)波頻散曲線，如圖4所示。然后，依據(jù)頻散曲線選擇不同的導(dǎo)波激勵(lì)頻率，得到不同激勵(lì)頻率時(shí)的導(dǎo)波檢測(cè)信號(hào)如圖5所示。

將由頻散曲線得出的充水管道中不同激勵(lì)頻率導(dǎo)波的群速度理論值與試驗(yàn)檢測(cè)值進(jìn)行比較，如表1所示，可見各頻率下試驗(yàn)值與理論值基本一致，誤差很小。圖6為兩者的對(duì)比曲線，可更直觀地表現(xiàn)出兩者間的關(guān)系，從圖中可看出試驗(yàn)值與計(jì)算曲線相吻合。由此證明了基于進(jìn)化規(guī)劃算法的頻散方程求解方法的正確性。

圖4 充水管道的頻散曲線

圖5 不同激勵(lì)頻率時(shí)導(dǎo)波的檢測(cè)信號(hào)

m·s-1

圖6 群速度試驗(yàn)結(jié)果與理論計(jì)算值的比較

3 結(jié)論

提出了一種求解導(dǎo)波頻散曲線的新方法，將進(jìn)化規(guī)劃算法應(yīng)用到頻散方程的求解中。相對(duì)于其他進(jìn)化算法，進(jìn)化規(guī)劃算法的模擬主要著眼于物種的進(jìn)化過程，沒有交叉，突變是其唯一的控制方法，操作控制比較簡(jiǎn)單。頻散曲線是反映導(dǎo)波在管道中的傳播特征的最直觀的手段之一，鑒于頻散方程求解的復(fù)雜性和進(jìn)化規(guī)劃算法解決復(fù)雜問題的能力，將進(jìn)化規(guī)劃算法應(yīng)用到流體管道的頻散曲線求解中。通過比較，其計(jì)算結(jié)果和試驗(yàn)值擬合較好，且計(jì)算所得的頻散曲線反應(yīng)出的導(dǎo)波特性與試驗(yàn)反應(yīng)出的導(dǎo)波特性是一致的。通過試驗(yàn)驗(yàn)證，充分證實(shí)該計(jì)算方法的可行性，且較其他方法便捷。文章以充水管道為研究對(duì)象，而同樣的方法可推廣到充滿黏性液體管道、埋地管道等頻散曲線的求解中，只是頻散方程和邊界條件不同，求解思路和過程是一致的。