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

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

0 引言

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

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

COMSOL Multiphysics軟件可以用于多物理場(chǎng)耦合仿真，該軟件的突出優(yōu)點(diǎn)是各物理場(chǎng)可在一套相同的網(wǎng)格節(jié)點(diǎn)中進(jìn)行雙向耦合計(jì)算。田振國等[8]和鄭杜成等[9]均通過該軟件構(gòu)建多物理場(chǎng)耦合模型，以計(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)確性。

針對(duì)上述模型中存在的不足，通過COMSOL Multiphysics軟件，采用基于Laplacian方程映射的動(dòng)網(wǎng)格移動(dòng)方法處理電樞運(yùn)動(dòng)問題。在考慮材料非線性的同時(shí)，構(gòu)建電磁軌道發(fā)射裝置電磁-溫度-結(jié)構(gòu)耦合場(chǎng)的三維瞬態(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)生磁場(chǎng)，與流經(jīng)電樞的電流相互作用產(chǎn)生洛倫茲力加速電樞運(yùn)動(dòng)。

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

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

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

1.1 電磁場(chǎng)控制方程

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

(1)

(2)

J=σE

(3)

式中：E,B,J分別為電場(chǎng)強(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 溫度場(chǎng)控制方程

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

(8)

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

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

基于動(dòng)量守恒的結(jié)構(gòu)場(chǎng)控制方程以張量的形式表達(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)體視為一彈性固體，對(duì)結(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)會(huì)相應(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)兩端分別對(duì)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。

對(duì)式(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)度會(huì)隨著溫度的升高而降低，相關(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è)置

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

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

如圖3所示，用四面體網(wǎng)格劃分幾何體。在“動(dòng)網(wǎng)格(ALE)”接口設(shè)置中，電樞的邊界位移wa由式(7)決定，其中電樞受到的電磁推力Fz由“磁場(chǎng)”接口計(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 電磁場(chǎng)仿真結(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è)；對(duì)于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)。對(duì)長方體導(dǎo)軌，在寬度方向上洛倫茲力密度均勻分布；在長度方向上，越靠近樞軌接觸面后側(cè)，洛倫茲力密度越大，這是因?yàn)樵撐恢锰庪娏髅芏燃?。?duì)于C型電樞結(jié)構(gòu)，洛倫茲力密度垂直于電樞臂分布，既有推動(dòng)電樞前進(jìn)的分力，又有使電樞外擴(kuò)以保持樞軌良好接觸的分力。

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

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

不同時(shí)刻電樞和導(dǎo)軌上的溫度分布如圖9所示，高溫區(qū)域主要集中于電樞的喉部和樞軌接觸面上。這是電流密度集中造成焦耳熱增大所導(dǎo)致的。電樞是載流元件，在焦耳熱的作用下，隨著時(shí)間的增加，電樞在發(fā)射過程中的溫升會(huì)越來越大。出口時(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)場(chǎng)仿真結(jié)果

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

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

4 結(jié)論

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