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

(海軍工程大學艦船綜合電力技術(shù)國防科技重點實驗室，武漢 430033)

0 引言

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

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

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

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

1 理論計算方程

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

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

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

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

1.1 電磁場控制方程

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

(1)

(2)

J=σE

(3)

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

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

(4)

電樞和導軌所受洛倫茲力的計算方程為：

(5)

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

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

(6)

(7)

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

1.2 溫度場控制方程

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

(8)

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

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

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

(9)

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

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

εT=αΔT

(10)

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

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

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

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

(11)

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

w=Nae

(12)

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

(13)

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

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

Kae=0

(14)

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

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

2.1 計算模型

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

圖3 電樞和導軌的幾何模型

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

通入導軌和電樞的激勵電流如圖4所示，峰值電流為200 kA，峰值電流時刻為2 ms，激勵時間為3 ms。

圖4 激勵電流曲線

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

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

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

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

(15)

(16)

Z為z軸上的坐標值，單位為m。

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

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

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

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

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

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

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

圖8為1.5 ms時洛倫茲力密度在電樞和導軌上的分布情況，洛倫茲力密度分布與電樞和導軌的結(jié)構(gòu)有關(guān)。對長方體導軌，在寬度方向上洛倫茲力密度均勻分布；在長度方向上，越靠近樞軌接觸面后側(cè)，洛倫茲力密度越大，這是因為該位置處電流密度集中。對于C型電樞結(jié)構(gòu)，洛倫茲力密度垂直于電樞臂分布，既有推動電樞前進的分力，又有使電樞外擴以保持樞軌良好接觸的分力。

圖8 1.5 ms時洛倫茲力密度分布

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

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

圖9 不同時刻電樞溫度分布圖

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

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

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

電樞和導軌上應(yīng)力分布如圖11所示，不同時刻電樞上的應(yīng)力主要集中在電樞喉部區(qū)域，這是因為電樞喉部電流密度最大。峰值電流時刻即2 ms時，喉部區(qū)域最高應(yīng)力為422 MPa。對比圖9可以看到，在不同時刻電樞喉部區(qū)域的溫升也最高，造成電樞喉部區(qū)域抗拉強度和屈服強度的下降。如果電樞喉部區(qū)域所受到的平均應(yīng)力高于該位置平均溫度下的屈服強度，則會導致電樞的變形甚至是斷裂。

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

4 結(jié)論

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