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

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

0 引言

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

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

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

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

1 理論計算方程

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

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

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

圖2 多物理場耦合關系圖

1.1 電磁場控制方程

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

(1)

(2)

J=σE

(3)

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

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

(4)

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

(5)

其中V為結構體的體積。

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

(6)

(7)

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

1.2 溫度場控制方程

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

(8)

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

1.3 結構場控制方程

基于動量守恒的結構場控制方程以張量的形式表達為：

(9)

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

當電樞和導軌產生溫升或溫度分布不均時，結構體將產生熱應變：

εT=αΔT

(10)

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

1.4 基于Laplacian方程映射的動網格移動

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

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

(11)

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

w=Nae

(12)

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

(13)

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

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

Kae=0

(14)

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

2 計算模型及仿真設置

2.1 計算模型

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

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

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

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

圖4 激勵電流曲線

2.2 仿真條件及設置

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

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

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

(15)

(16)

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

3 仿真結果及分析

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

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

3.1 電磁場仿真結果

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

圖6 電樞和導軌內的電流密度分布

圖7 電樞和導軌內的磁感應強度分布

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

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

3.2 溫度場仿真結果

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

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

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

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

3.3 結構場仿真結果

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

圖11 電樞應力分布

4 結論

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