邊光滑有限元-邊界元耦合法計(jì)算二維瞬態(tài)渦流場(chǎng)

王洋洋，蔣興良

(輸配電裝備及系統(tǒng)安全與新技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室(重慶大學(xué))，重慶 400044)

在瞬態(tài)渦流場(chǎng)研究中，常用的數(shù)值計(jì)算方法有積分法[1-2]、有限元法[3]、邊界元法以及有限元邊界元耦合法[4-7]. 雖然傳統(tǒng)有限元法(FEM)[8]作為最成功的數(shù)值方法之一，已被廣泛用于解決工程和科學(xué)的各種實(shí)際問題，但是因?yàn)槠溆泄逃械娜秉c(diǎn)，如剖分網(wǎng)格密度以及計(jì)算階數(shù)的復(fù)雜程度都對(duì)有限元計(jì)算精度影響非常大. 單元數(shù)量的增加，及其插值函數(shù)精度的提高，才能得到更精確的數(shù)值解. 特別是在某些復(fù)雜的設(shè)計(jì)優(yōu)化實(shí)例中，有限元法需要重復(fù)多次進(jìn)行網(wǎng)格剖分以及重要部分網(wǎng)格加密操作，這不僅增加了設(shè)計(jì)的難度，而且還降低了設(shè)計(jì)計(jì)算的效率. 因此國(guó)內(nèi)外學(xué)者一直致力于研究如何克服傳統(tǒng)有限元法的缺點(diǎn).

近年來(lái)，在努力開發(fā)先進(jìn)數(shù)值方法時(shí)，文獻(xiàn)[9]將無(wú)網(wǎng)格中的應(yīng)變光滑技術(shù)與傳統(tǒng)有限元法相結(jié)合，引入梯度光滑的操作概念，提出一系列廣義梯度光滑有限元法. 其包括用梯度光滑技術(shù)分別對(duì)構(gòu)造的節(jié)點(diǎn)光滑域、邊光滑域以及面光滑域進(jìn)行光滑操作而形成的點(diǎn)光滑有限元法(NS-FEM)[10]、邊光滑有限元法(ES-FEM)[11]以及面光滑有限元法(FS-FEM)[12]. 3種方法各有優(yōu)缺點(diǎn). 相對(duì)于傳統(tǒng)有限元法，NS-FEM在精度和效率方面有所提高. 但是NS-FEM在空間上是穩(wěn)定的在時(shí)間上可能是不穩(wěn)定的，并且由于“過軟”特性而不能直接應(yīng)用于解決動(dòng)態(tài)問題[13]. 然而，ES-FEM可以克服時(shí)間不穩(wěn)定性[14]，由于本文研究的是時(shí)變渦流場(chǎng)，因此克服時(shí)間不穩(wěn)定性對(duì)于時(shí)變渦流場(chǎng)分析尤其重要. 此外，對(duì)于固體力學(xué)問題，具有線性積分的ES-FEM比線性有限元法更加準(zhǔn)確和有效[15]. 至此，ES-FEM近年來(lái)得到了廣泛的應(yīng)用. 對(duì)于FS-FEM更適用于三維問題.

本文提出一種用ES-FEM與邊界元耦合法計(jì)算二維瞬態(tài)渦流場(chǎng). 該耦合法同時(shí)具有ES-FEM計(jì)算精度高以及邊界元法占用計(jì)算機(jī)內(nèi)存少等優(yōu)點(diǎn). 首先將問題域離散化為一組三角形單元，然后進(jìn)一步構(gòu)造與三角形單元邊相關(guān)的積分域，采用伽遼金加權(quán)余量法，推導(dǎo)出基于邊光滑有限元-邊界元耦合算法的離散公式，在導(dǎo)體區(qū)域內(nèi)采用ES-FEM，非導(dǎo)體區(qū)域內(nèi)采用邊界元法. 研究導(dǎo)體區(qū)域的磁場(chǎng)分布，并與測(cè)量結(jié)果和有限元-邊界元耦合法計(jì)算結(jié)果進(jìn)行對(duì)比. 計(jì)算結(jié)果表明，在相同網(wǎng)格密度條件下，基于ES-FEM邊界元耦合法計(jì)算精度更高，該方法具有更好的應(yīng)用前景.

1 電磁場(chǎng)的分析模型

電磁分析的問題實(shí)際上是在給定邊界條件下求解一組麥克斯韋方程組的問題. 如圖1所示，為開域渦流場(chǎng)求解區(qū)域示意圖. 整個(gè)電磁場(chǎng)求解域分為渦流區(qū)V1和非渦流區(qū)V2，Γ1是V1和V2的內(nèi)部邊界，非渦流區(qū)V2包括無(wú)源電流的非渦流區(qū)和有源電流的非渦流區(qū).V1區(qū)域內(nèi)采用ES-FEM，V2區(qū)域內(nèi)采用BEM, 在內(nèi)部邊界上滿足的邊界條件為

(1)

(2)

ν

(3)

圖1 開域渦流場(chǎng)求解示意圖

2 邊光滑有限元-邊界元耦合法介紹

2.1 邊光滑有限元法

邊光滑有限元法是在求解有限元的基礎(chǔ)上對(duì)其求解域進(jìn)一步光滑，形成多個(gè)光滑子單元，然后在新形成的光滑子單元內(nèi)引入梯度光滑操作，梯度光滑操作是指在新形成的光滑域上對(duì)有限元求解的系數(shù)矩陣進(jìn)行光滑平均運(yùn)算，其數(shù)學(xué)表達(dá)式為

其中Ax為光滑域的面積. 圖2為邊光滑有限元法的計(jì)算流程圖.

圖2 計(jì)算流程圖

對(duì)于二維求解域，先用劃分軟件Hypermesh對(duì)其進(jìn)行三角形劃分，其中生成N個(gè)節(jié)點(diǎn)和p條邊. 傳統(tǒng)有限元的離散矩陣[16]為

(4)

(5)

式中AL為光滑域的面積，在二維求解域單元內(nèi)，AL的表達(dá)式可根據(jù)參考文獻(xiàn)[9]計(jì)算得到.

圖3 基于邊的光滑域

對(duì)式(7)進(jìn)行格林公式換算，那么對(duì)光滑域內(nèi)位函數(shù)的面積分可以轉(zhuǎn)化為光滑域邊界上的線積分，即

(6)

將式(6)代入式(5)中，可得

(7)

其中：

(8)

(9)

對(duì)光滑域ΩL采用高斯積分，式(12)的光滑域邊界積分形式為

(10)

式中：NA為電磁場(chǎng)求解域邊界L的個(gè)數(shù)，NB為高斯點(diǎn)的個(gè)數(shù)，δq為高斯點(diǎn)積分的權(quán)函數(shù).

由式(10)可知，dij是由光滑域各個(gè)邊中點(diǎn)處的形函數(shù)求得，而各個(gè)中點(diǎn)處的形函數(shù)可通過場(chǎng)節(jié)點(diǎn)插值得到，如圖4所示，點(diǎn)1、2、3、4分別為場(chǎng)節(jié)點(diǎn)，A、B分別為中心節(jié)點(diǎn)，g1、g2、g3、g4分別為高斯積分點(diǎn)[17].

基于邊L的光滑有限元系統(tǒng)的局部光滑域的系數(shù)矩陣可以寫成：

(11)

對(duì)上述局部系數(shù)矩陣進(jìn)行整體組裝，得到域內(nèi)整體系數(shù)矩陣為

(12)

圖4 ES-FEM形函數(shù)分布

2.2 邊界元法

對(duì)于非渦流區(qū)應(yīng)用格林公式，可得對(duì)應(yīng)邊界積分方程[5]為

(13)

采用伽遼金余量法對(duì)式(13)離散，得

(14)

2.3 基于邊光滑有限元邊界元耦合算法

將式(4)和式(14)聯(lián)立，可得

(15)

對(duì)式(15)中的導(dǎo)數(shù)項(xiàng)，采用向后差分，即

(16)

(17)

(18)

3 計(jì)算算例與分析

3.1 脈沖線圈-鋁板表面磁場(chǎng)分布

脈沖線圈-鋁板的結(jié)構(gòu)參數(shù)：矩形截面為0.6 mm×4.8 mm的銅絲纏繞成內(nèi)徑、外徑分別為3.2、25.5 mm,匝數(shù)為30的圓柱形線圈[19]. 脈沖線圈流過的瞬時(shí)電流波形，如圖5所示. 在線圈正上方放置0.8 mm厚的鋁板，電導(dǎo)率δ=3.48×107s/m. 線圈與鋁板的間隙為2 mm. 采用本文所提的基于邊有限元邊界元耦合法求得鋁板表面的磁感應(yīng)強(qiáng)度，其結(jié)果分別與測(cè)量值以及有限元-邊界元法的計(jì)算結(jié)果進(jìn)行對(duì)比，如圖6～8所示.

圖5 脈沖電流隨時(shí)間變化

圖6 磁矢量位A線的分布云圖

圖7 切向磁感應(yīng)強(qiáng)度隨時(shí)間的變化(r=10 mm)

由圖6可知，離線圈越近，磁矢量位A值越大，與實(shí)際情況一致，說(shuō)明本文的計(jì)算方法合理. 為了驗(yàn)證本文所提方法的計(jì)算精度，設(shè)相對(duì)誤差Bi為

(19)

式中：B0為測(cè)量值，B1是計(jì)算值，i=x、y分別代表磁場(chǎng)沿切向、垂直方向的磁場(chǎng)分量.

圖8 垂直磁感應(yīng)強(qiáng)度隨時(shí)間的變化(r=7.6 mm)

由圖7、8可知，基于邊光滑有限元邊界元耦合法與測(cè)量結(jié)果一致，因此可用本文所提的方法計(jì)算渦流場(chǎng). 在網(wǎng)格單元?jiǎng)澐置芏认嗤臈l件下，運(yùn)行兩種方法效率對(duì)比見表1，雖然本文所提方法的運(yùn)行時(shí)間稍長(zhǎng)于有限元邊界元耦合法，但其計(jì)算精度遠(yuǎn)遠(yuǎn)高于有限元邊界元耦合法，另外，對(duì)于目前計(jì)算機(jī)的計(jì)算能力而言，本文通過改變模型的大小，對(duì)比兩種方法的運(yùn)行時(shí)間. 由表2可知，兩種方法運(yùn)行時(shí)間的差別可以忽略. 由此可知，在計(jì)算渦流場(chǎng)時(shí)基于邊光滑有限元邊界元法能大大提高計(jì)算精度.

表1 效率對(duì)比

表2 運(yùn)行時(shí)間對(duì)比

綜上所述，在進(jìn)行渦流場(chǎng)分析時(shí)，應(yīng)優(yōu)先考慮基于邊光滑有限元邊界元耦合法.

3.2 鋁板表面磁場(chǎng)分布的影響因素

3.2.1 電流密度變化的影響

脈沖電流-鋁板模型參數(shù)如上所述，在進(jìn)行電流密度對(duì)其磁場(chǎng)影響分析時(shí)，保持模型的其他參數(shù)不變，改變流過脈沖線圈的電流密度Js從10×106A/m2到180×106A/m2變化，間隔為30×106A/m2. 則上述電流密度變化的情況下鋁板表面磁場(chǎng)分量分布如圖9所示.

圖9 磁場(chǎng)強(qiáng)度隨渦流密度的變化

如圖9所示，磁場(chǎng)強(qiáng)度隨渦流密度的增大而增大，這是因?yàn)闇u流密度Js與系數(shù)矩陣F成正相關(guān)，即當(dāng)Js增大時(shí)，F(xiàn)增大. 由式(20)可知,隨著F增大，所求磁矢量位A增大，因此磁場(chǎng)強(qiáng)度增大.

3.2.2 脈沖線圈-鋁板間隙變化的影響

脈沖電流-鋁板模型參數(shù)如上所述，在進(jìn)行脈沖線圈-鋁板間隙對(duì)其磁場(chǎng)影響分析時(shí)，保持模型的其他參數(shù)不變，改變間隙h從2 mm到8 mm之間變化，間隔為1 mm. 則上述間隙變化的情況下鋁板表面磁場(chǎng)分量分布如圖10所示.

圖10 磁場(chǎng)強(qiáng)度隨間隙的變化

由圖10可知，隨著脈沖線圈與鋁板之間間隙h的增大，磁場(chǎng)強(qiáng)度逐漸減小. 這說(shuō)明如果為了增大鋁板中的渦流作用，即產(chǎn)生大的脈沖力，那么脈沖線圈需盡可能近的接近鋁板.

4 結(jié) 論

1)提出了一種計(jì)算渦流場(chǎng)分布的基于邊光滑有限元邊界元耦合法，該方法即具有ES-FEM計(jì)算精度高和邊界元法占用內(nèi)存少的優(yōu)點(diǎn)，并詳細(xì)給出了該方法的推導(dǎo)過程. 同時(shí)，利用本文所提的基于邊光滑有限元邊界元耦合法計(jì)算分析脈沖線圈上方鋁板表面的磁場(chǎng)分布，與測(cè)量結(jié)果對(duì)比表明，該方法是正確合理的.

2)在單元網(wǎng)格密度劃分相同的條件下，將本文得到的計(jì)算結(jié)果與傳統(tǒng)有限元邊界元耦合法計(jì)算得到的結(jié)果對(duì)比，雖然本文所提方法的運(yùn)行時(shí)間稍長(zhǎng)于有限元邊界元耦合法，但其計(jì)算精度遠(yuǎn)遠(yuǎn)高于有限元邊界元耦合法，另外，對(duì)于目前計(jì)算機(jī)的計(jì)算能力而言，兩種方法運(yùn)行時(shí)間的差別可以忽略.

3)在計(jì)算模型的其他參數(shù)保持不變，只改變渦流密度的情況下，磁場(chǎng)強(qiáng)度隨渦流密度的增大而增大且曲線增大率幾乎保持不變，同理，磁感應(yīng)強(qiáng)度隨著脈沖線圈與鋁板之間間隙h的增大而減小，且其曲線降低率也幾乎保持不變.