電磁軌道發(fā)射裝置的動(dòng)態(tài)多物理場耦合分析

楊 帆，翟小飛，張 曉，劉 華

(海軍工程大學(xué)艦船綜合電力技術(shù)國防科技重點(diǎn)實(shí)驗(yàn)室，武漢 430033)

0 引言

伴隨著電樞的高速運(yùn)動(dòng)，電磁軌道發(fā)射裝置內(nèi)電磁場、溫度場和結(jié)構(gòu)場相互耦合在一起，使裝置內(nèi)彈道工作環(huán)境十分惡劣[1]。為了設(shè)計(jì)出既能保證樞軌間有足夠的接觸力而不發(fā)生轉(zhuǎn)捩，又能保證不會因溫度過高、受力過大而產(chǎn)生破壞的電樞結(jié)構(gòu)，研究發(fā)射狀態(tài)下的內(nèi)彈道多物理場耦合模型十分必要。

目前對構(gòu)建電磁軌道發(fā)射裝置的多物理場耦合計(jì)算模型已取得了一定成果：在國外，一些計(jì)算程序如EMAP3D[2]，MEGA[3]和HERB[4]等被開發(fā)出來，用于模擬高速滑動(dòng)電接觸條件下的電磁-溫度-結(jié)構(gòu)耦合，并得到了一些基本規(guī)律。在國內(nèi)，譚賽等[5-6]采用自由度平移法建立了考慮摩擦熱的電磁-溫度-運(yùn)動(dòng)耦合場有限元數(shù)值計(jì)算模型，并與由EMAP3D所得的結(jié)果相吻合，但沒有對結(jié)構(gòu)場進(jìn)行耦合計(jì)算。Lin等[7]利用Ls-dyna軟件，基于有限元和邊界元耦合計(jì)算的方法，求解了電磁發(fā)射過程中電樞和導(dǎo)軌上電流密度、溫度和應(yīng)力的分布情況，但該模型只進(jìn)行了單向耦合分析，且沒有考慮材料的非線性問題。

COMSOL Multiphysics軟件可以用于多物理場耦合仿真，該軟件的突出優(yōu)點(diǎn)是各物理場可在一套相同的網(wǎng)格節(jié)點(diǎn)中進(jìn)行雙向耦合計(jì)算。田振國等[8]和鄭杜成等[9]均通過該軟件構(gòu)建多物理場耦合模型，以計(jì)算導(dǎo)軌上的應(yīng)力分布情況。但前者沒有考慮速度趨膚效應(yīng)帶來的影響；后者僅通過“移動(dòng)電導(dǎo)率”的方式處理電樞運(yùn)動(dòng)問題，忽略了電樞結(jié)構(gòu)帶來的影響。王世禮[10]、張占群[11]詳細(xì)介紹了動(dòng)網(wǎng)格移動(dòng)方法的種類及特點(diǎn)，采用基于Laplacian方程映射的動(dòng)網(wǎng)格移動(dòng)方法能夠有效保證結(jié)構(gòu)體邊界上網(wǎng)格的質(zhì)量，有利于提高模型的收斂性和準(zhǔn)確性。

針對上述模型中存在的不足，通過COMSOL Multiphysics軟件，采用基于Laplacian方程映射的動(dòng)網(wǎng)格移動(dòng)方法處理電樞運(yùn)動(dòng)問題。在考慮材料非線性的同時(shí)，構(gòu)建電磁軌道發(fā)射裝置電磁-溫度-結(jié)構(gòu)耦合場的三維瞬態(tài)有限元模型。最后經(jīng)模型計(jì)算，得到在電樞運(yùn)動(dòng)情況下的電流、溫度和應(yīng)力的分布特點(diǎn)。

1 理論計(jì)算方程

電磁軌道發(fā)射裝置的工作原理[12]如圖1所示，電流由上導(dǎo)軌流入，經(jīng)過電樞從下導(dǎo)軌流出，形成一閉合回路。電流經(jīng)過導(dǎo)軌時(shí)在膛內(nèi)產(chǎn)生磁場，與流經(jīng)電樞的電流相互作用產(chǎn)生洛倫茲力加速電樞運(yùn)動(dòng)。

圖1 電磁發(fā)射工作原理圖

電磁軌道發(fā)射裝置內(nèi)多物理場耦合關(guān)系如圖2所示。電磁場中產(chǎn)生的電磁力和焦耳熱在電樞和導(dǎo)軌上的分布，分別影響著結(jié)構(gòu)場中的應(yīng)力和溫度場中的溫度在電樞和導(dǎo)軌上的分布情況；溫度場中產(chǎn)生的溫升會降低電樞和導(dǎo)軌材料的電導(dǎo)率和抗拉/屈服強(qiáng)度，進(jìn)而分別影響電磁場中產(chǎn)生的焦耳熱和結(jié)構(gòu)場中的應(yīng)力應(yīng)變；結(jié)構(gòu)場中電樞高速滑動(dòng)產(chǎn)生的摩擦熱是溫度場中溫升的熱量來源之一，樞軌接觸面上的接觸力分布影響著電磁場中電流密度在樞軌接觸面上的分布情況。

圖2 多物理場耦合關(guān)系圖

1.1 電磁場控制方程

建立由電磁感應(yīng)定律、安培環(huán)路定律和歐姆定律為基礎(chǔ)的電磁學(xué)控制方程，其微分形式分別為：

(1)

(2)

J=σE

(3)

式中：E,B,J分別為電場強(qiáng)度,磁感應(yīng)強(qiáng)度,電流密度；μ為材料的磁導(dǎo)率；σ為材料的電導(dǎo)率；v為電樞速度，只有z軸方向上有分量ν。

將式(1)～式(3)聯(lián)立，可以得到用于計(jì)算磁感應(yīng)強(qiáng)度和電流密度分布與擴(kuò)散的方程：

(4)

電樞和導(dǎo)軌所受洛倫茲力的計(jì)算方程為：

(5)

其中V為結(jié)構(gòu)體的體積。

在不考慮電樞受到的摩擦力和空氣阻力的情況下，電樞的速度與位移方程為：

(6)

(7)

式中：Fz為電樞受到的電磁推力，即洛倫茲力在z軸方向上的分力；ma為電樞質(zhì)量；t為時(shí)間；wa為電樞位移;w0為電樞在導(dǎo)軌上的初始位置。

1.2 溫度場控制方程

如果只考慮電流焦耳熱作用，忽略接觸電阻焦耳熱和摩擦熱，導(dǎo)軌和電樞內(nèi)的熱傳導(dǎo)方程可寫為：

(8)

式中：T為溫度；j為電流密度模；ρ,c,κ分別為材料的密度,比熱容和熱傳導(dǎo)系數(shù)。

1.3 結(jié)構(gòu)場控制方程

基于動(dòng)量守恒的結(jié)構(gòu)場控制方程以張量的形式表達(dá)為：

(9)

s,f，w分別表示應(yīng)力張量，單位體積力和位移。

當(dāng)電樞和導(dǎo)軌產(chǎn)生溫升或溫度分布不均時(shí)，結(jié)構(gòu)體將產(chǎn)生熱應(yīng)變：

εT=αΔT

(10)

式中：α為材料的熱膨脹系數(shù)；ΔT為溫升。

1.4 基于Laplacian方程映射的動(dòng)網(wǎng)格移動(dòng)

基于Laplacian方程映射的動(dòng)網(wǎng)格移動(dòng)方法，是一種基于有限元法的網(wǎng)格移動(dòng)方法。其核心是將結(jié)構(gòu)體視為一彈性固體，對結(jié)構(gòu)體進(jìn)行網(wǎng)格劃分，當(dāng)結(jié)構(gòu)體的邊界移動(dòng)時(shí)，結(jié)構(gòu)體內(nèi)部網(wǎng)格節(jié)點(diǎn)的位移滿足Laplacian方程。即當(dāng)結(jié)構(gòu)體邊界移動(dòng)時(shí)，內(nèi)部網(wǎng)格節(jié)點(diǎn)會相應(yīng)的產(chǎn)生位移從而使網(wǎng)格移動(dòng)。

當(dāng)結(jié)構(gòu)體邊界移動(dòng)時(shí)，以求解內(nèi)部網(wǎng)格節(jié)點(diǎn)在z軸方向上的位移w為例，其滿足的泛函形式方程為：

(11)

在有限元計(jì)算方法中，網(wǎng)格節(jié)點(diǎn)位移以矩陣形式表示為：

w=Nae

(12)

式中:N為插值函數(shù)矩陣;ae為節(jié)點(diǎn)位移矢量。式(12)兩端分別對x，y，z取偏導(dǎo)可得：

(13)

Ni表示節(jié)點(diǎn)i處的插值函數(shù)，假如用四面體網(wǎng)格劃分結(jié)構(gòu)體，則每一網(wǎng)格上有4個(gè)節(jié)點(diǎn)，所以i分別取1,2,3,4。

對式(11)中的泛函取變分，并將式(13)代入，則有：

Kae=0

(14)

式中：K為n×n的矩陣，n為網(wǎng)格節(jié)點(diǎn)的數(shù)量。通過把結(jié)構(gòu)體的邊界位移代入到式(14)中，便能求解出結(jié)構(gòu)體內(nèi)部網(wǎng)格節(jié)點(diǎn)在z方向上的位移。依此求解過程，可以計(jì)算出節(jié)點(diǎn)在x和y方向上的位移。

2 計(jì)算模型及仿真設(shè)置

2.1 計(jì)算模型

如圖3所示，為電樞和導(dǎo)軌的幾何模型及網(wǎng)格劃分。導(dǎo)軌為長方體結(jié)構(gòu)，長0.9 m，寬40 mm，高20 mm，間距30 mm，材料為銅合金。電樞為C型結(jié)構(gòu)，材料為鋁合金，電樞質(zhì)量ma=81.46 g。初始時(shí)刻，電樞末端距離導(dǎo)軌末段為0.12 m，樞軌接觸面長度為36 mm，因此初始時(shí)刻樞軌接觸面前端距離導(dǎo)軌末端為0.156 m，并將此長度范圍內(nèi)的導(dǎo)軌部分記為導(dǎo)軌1。在電樞運(yùn)動(dòng)過程中，即將與電樞接觸的導(dǎo)軌部分記為導(dǎo)軌2。材料的電導(dǎo)率和屈服強(qiáng)度會隨著溫度的升高而降低，相關(guān)材料隨溫度變化的性能參數(shù)如表1所示。

圖3 電樞和導(dǎo)軌的幾何模型

表1 材料性能參數(shù)

通入導(dǎo)軌和電樞的激勵(lì)電流如圖4所示，峰值電流為200 kA，峰值電流時(shí)刻為2 ms，激勵(lì)時(shí)間為3 ms。

圖4 激勵(lì)電流曲線

2.2 仿真條件及設(shè)置

為簡化計(jì)算，仿真過程中作如下假設(shè)：不考慮電樞磨損帶來的影響；不考慮接觸電阻焦耳熱和摩擦熱的影響；不考慮摩擦力和空氣阻力對電樞運(yùn)動(dòng)帶來的影響；初始時(shí)刻電樞與導(dǎo)軌接觸均勻。

在COMSOL Multiphysics中添加“磁場”、“電路”、“固體傳熱”、“固體力學(xué)”、“電磁熱”和“熱膨脹”等物理場接口，同時(shí)設(shè)置“動(dòng)網(wǎng)格(ALE)”接口用于計(jì)算電樞動(dòng)態(tài)發(fā)射過程中的多物理場耦合問題。

如圖3所示，用四面體網(wǎng)格劃分幾何體。在“動(dòng)網(wǎng)格(ALE)”接口設(shè)置中，電樞的邊界位移wa由式(7)決定，其中電樞受到的電磁推力Fz由“磁場”接口計(jì)算得到。導(dǎo)軌1和導(dǎo)軌2在電樞運(yùn)動(dòng)過程中所設(shè)置的邊界位移w1和w2分別為：

(15)

(16)

Z為z軸上的坐標(biāo)值，單位為m。

3 仿真結(jié)果及分析

仿真計(jì)算得到電樞在所受電磁推力作用下的速度、位移曲線如圖5所示，電樞最大位移為0.75 m，出口最大速度為450 m/s。

圖5 電樞的位移和速度曲線

3.1 電磁場仿真結(jié)果

發(fā)射過程中，電流密度和磁感應(yīng)強(qiáng)度在電樞和導(dǎo)軌上的分布如圖6～圖7所示。受速度趨膚效應(yīng)的影響，電樞和導(dǎo)軌內(nèi)的電流密度和磁感應(yīng)強(qiáng)度呈現(xiàn)出相同的分布規(guī)律：電流密度與磁感應(yīng)強(qiáng)度主要分布在導(dǎo)軌內(nèi)側(cè)與電樞后側(cè)的表層部分，并分別向?qū)к壨鈧?cè)和電樞前側(cè)擴(kuò)散；電流密度和磁感應(yīng)強(qiáng)度更集中于樞軌接觸面的后側(cè)；對于C型電樞結(jié)構(gòu)，電流密度和磁感應(yīng)強(qiáng)度集中在電樞臂轉(zhuǎn)角處的喉部分布，該位置處電流密度和磁感應(yīng)強(qiáng)度的值最高。

圖6 電樞和導(dǎo)軌內(nèi)的電流密度分布

圖7 電樞和導(dǎo)軌內(nèi)的磁感應(yīng)強(qiáng)度分布

圖8為1.5 ms時(shí)洛倫茲力密度在電樞和導(dǎo)軌上的分布情況，洛倫茲力密度分布與電樞和導(dǎo)軌的結(jié)構(gòu)有關(guān)。對長方體導(dǎo)軌，在寬度方向上洛倫茲力密度均勻分布；在長度方向上，越靠近樞軌接觸面后側(cè)，洛倫茲力密度越大，這是因?yàn)樵撐恢锰庪娏髅芏燃小τ贑型電樞結(jié)構(gòu)，洛倫茲力密度垂直于電樞臂分布，既有推動(dòng)電樞前進(jìn)的分力，又有使電樞外擴(kuò)以保持樞軌良好接觸的分力。

圖8 1.5 ms時(shí)洛倫茲力密度分布

3.2 溫度場仿真結(jié)果

不同時(shí)刻電樞和導(dǎo)軌上的溫度分布如圖9所示，高溫區(qū)域主要集中于電樞的喉部和樞軌接觸面上。這是電流密度集中造成焦耳熱增大所導(dǎo)致的。電樞是載流元件，在焦耳熱的作用下，隨著時(shí)間的增加，電樞在發(fā)射過程中的溫升會越來越大。出口時(shí)刻，電樞最高溫度已達(dá)386 ℃。

圖9 不同時(shí)刻電樞溫度分布圖

圖10為動(dòng)態(tài)發(fā)射過程中不同時(shí)刻熱功率密度在電樞上的分布情況。從圖中可以看出，熱功率密度分布首先出現(xiàn)在樞軌接觸面上的后側(cè)，隨著時(shí)間的推進(jìn)，逐漸沿著電樞長度方向向接觸面的前側(cè)集中。這是由速度趨膚效應(yīng)作用所產(chǎn)生的熔化波燒蝕現(xiàn)象[7,13]。

圖10 電樞上熱功率密度分布

3.3 結(jié)構(gòu)場仿真結(jié)果

電樞和導(dǎo)軌上應(yīng)力分布如圖11所示，不同時(shí)刻電樞上的應(yīng)力主要集中在電樞喉部區(qū)域，這是因?yàn)殡姌泻聿侩娏髅芏茸畲蟆７逯惦娏鲿r(shí)刻即2 ms時(shí)，喉部區(qū)域最高應(yīng)力為422 MPa。對比圖9可以看到，在不同時(shí)刻電樞喉部區(qū)域的溫升也最高，造成電樞喉部區(qū)域抗拉強(qiáng)度和屈服強(qiáng)度的下降。如果電樞喉部區(qū)域所受到的平均應(yīng)力高于該位置平均溫度下的屈服強(qiáng)度，則會導(dǎo)致電樞的變形甚至是斷裂。

圖11 電樞應(yīng)力分布

4 結(jié)論

通過構(gòu)建電磁軌道發(fā)射裝置的動(dòng)態(tài)多物理場耦合計(jì)算模型，仿真計(jì)算了電流密度、溫度、應(yīng)力等物理量在電樞和導(dǎo)軌上的分布情況。在所構(gòu)建的有限元計(jì)算模型中，各物理場通過一套相同的網(wǎng)格節(jié)點(diǎn)進(jìn)行直接耦合計(jì)算，能夠得到受速度趨膚效應(yīng)影響的電磁擴(kuò)散和熔化波燒蝕等在電磁軌道發(fā)射中出現(xiàn)的典型現(xiàn)象。仿真結(jié)果表明，構(gòu)建的仿真模型對描述和解決電磁軌道發(fā)射多物理場耦合問題是可行的，可以為電樞的優(yōu)化設(shè)計(jì)工作提供一定的指導(dǎo)。