蘆 逸 云

(洛陽理工學(xué)院 電氣工程與自動化學(xué)院，河南 洛陽 471023)

三維模型下永磁軌道上超導(dǎo)體臨界電流對磁懸浮力影響研究

蘆 逸 云

(洛陽理工學(xué)院 電氣工程與自動化學(xué)院，河南 洛陽 471023)

采用數(shù)值計算方法對高溫超導(dǎo)(HTSC)塊材電磁場進(jìn)行模擬，計算了永磁軌道(PMR)上方超導(dǎo)塊材所受到的磁懸浮力。文中給出了一種3D數(shù)值計算。模型采用了電磁場磁矢量法(H-方法)；基于Fortran語言開發(fā)平臺編寫了有限元求解代碼；使用E-J冪指數(shù)模型來描述高溫超導(dǎo)塊材的電場強(qiáng)度E和電流密度J的關(guān)系；采用Kim模型來描述高溫超導(dǎo)塊材的臨界電流密度Jc；使用了兩塊虛擬超導(dǎo)塊材疊加的方法處理高溫超導(dǎo)材料臨界電流密度各向異性的特性。采用文中描述的數(shù)值計算方法，成功地對由一個永磁軌道和一塊高溫超導(dǎo)塊材組成的磁懸浮系統(tǒng)進(jìn)行了計算研究，研究了磁懸浮系統(tǒng)中超導(dǎo)塊材的臨界電流密度對磁懸浮力的影響。

超導(dǎo); HTSC；三維模型；數(shù)值模擬

高溫超導(dǎo)塊材能夠很穩(wěn)定地懸浮(或者懸掛)在永磁體上方(或者下方)。這種內(nèi)在的自穩(wěn)定懸浮(或者懸掛)特性使得無摩擦軸承、儲能飛輪和自穩(wěn)定磁懸浮傳輸系統(tǒng)的潛在運用成為可能[1-2]。為了優(yōu)化設(shè)計這些設(shè)備的性能，很有必要對高溫超導(dǎo)塊材的電磁特性進(jìn)行模擬計算。由于YBCO塊材內(nèi)部存在仔晶邊界效應(yīng)，其內(nèi)部臨界電流密度沿a-b面的大小是電流密度沿c-軸的兩倍。對于其E-J本構(gòu)關(guān)系，與常導(dǎo)體的描述完全不相同。E-J冪指數(shù)假定電場強(qiáng)度E是與電流密度J的冪指數(shù)緊密相關(guān)的。Kim 模型假定高溫超導(dǎo)塊材的臨界電流密度Jc不僅僅與該位置的電密度J相關(guān)，還和該位置的磁場強(qiáng)度大小相關(guān)， 這些是高溫超體的E-J非線性特性和高溫超導(dǎo)材料非線性特性。由于高溫超導(dǎo)體的這些電磁強(qiáng)非線性特性，很難用數(shù)值的方法來模擬高溫超導(dǎo)體的電磁特性。

在最近這些年，針對此類電磁問題，出現(xiàn)了一些數(shù)值計算方法和解析計算方法。對于2D模型，局限于與其相關(guān)的簡單幾何模型以及在各向同性的假設(shè)下和軸對稱外磁場的假設(shè)下，研究者提出了幾種解析計算方法[3-5]。對于3D模型，由于很難處理高溫超導(dǎo)體的電磁非線性特性和受限于電腦運算能力，很少有相關(guān)報到見諸。

本文給出了一個基于有限元方法的3D模型，運用此數(shù)值模型對高溫超導(dǎo)電磁特性進(jìn)行了模擬計算。模型使用了兩塊虛擬超導(dǎo)塊材疊加的處理方法來解決非線性材料問題。模型中使用了一種迭代的方法用于迭代求解E-J非線性特性，還使用了向前差分法用于處理時間空間問題。有限元求解代碼是在FORTRAN語言平臺下編寫的。

運用文中給出的3D模型研究了由高溫超導(dǎo)塊材與永磁軌道組成的磁懸浮系統(tǒng)中超導(dǎo)塊材的臨界電流密度對磁懸浮力的影響。仿真計算的結(jié)果顯示了超導(dǎo)塊材的臨界電流密度的增加能夠增加最大磁懸浮力，并且最大磁懸浮力的增加幅度受到超導(dǎo)塊材所受到的外磁場大小的限制。

1 基本方程

本節(jié)介紹描述高溫超導(dǎo)塊材電磁特性的基本方程?？紤]一個高溫超導(dǎo)三維幾何塊材模型Ω1 和非超導(dǎo)區(qū)域Ω2 ，如圖1所示。

圖1 三維超溫超導(dǎo)模型圖

Maxwell 安培定律和法拉第電磁感應(yīng)定律:

×H=J,。

(1)

B=μH。

(2)

對于高溫超導(dǎo)體，其E-J本構(gòu)關(guān)系可描述為：

(3)

其中：Esc、Jsc和Hsc分別是高溫超導(dǎo)塊材的電場強(qiáng)度、電流密度和磁場強(qiáng)度。 Ec0是描述材料特性的一個溫度相關(guān)常數(shù)， Jc是高溫超導(dǎo)塊材的臨界電密度。

臨界電流密度是磁場強(qiáng)度相關(guān)的，可以由Kim模型來描述:

(4)

其中： Jc0和Hc0分別是與材料特性和溫度相關(guān)的常數(shù)。

為了簡化微分方程的推演過程，我們使用了虛擬歐姆定律：

(5)

其中：σsc是高溫超導(dǎo)塊材的有效電導(dǎo)率：

(6)

對于高溫超導(dǎo)三維模型的非線性材料問題，研究認(rèn)為沿c-軸的臨界電流密度Jc是沿a-b面臨界電流密度的三分之一。我們假定沿c-軸流動的電流密度所受到的電阻率三倍于沿a-b平面流動的電流密度所受到的電阻率。方程(5)又可以寫成：

(7)

將方程 (1) 和方程 (2) 代入方程 (7)，可以得到：

(8)

高溫導(dǎo)塊材數(shù)學(xué)上由兩部分組成，一部分是由各向同性超導(dǎo)材料組成，另一部分是由僅由沿c-軸方向電導(dǎo)率不為零的超導(dǎo)塊材組成。

著眼于所有計算區(qū)域的超導(dǎo)及非超導(dǎo)區(qū)域，將方程(1) 代入方程 (8) ，并且考慮到Jscy=?Hx/?z-?Hz/?x,于是方程(8)可以寫成如下方程：

(9)

其中：σ=σair,λ=0是指子區(qū)域Ω2；σ=σsc,λ=1 是指子區(qū)域Ω1 (見圖1)； 并且定義矢量Q為Q=[0 2(?Hx/?z-?Hz/?x) 0]T。

式(9)是一個含時間變化參量的電磁偏微分方程，它是由法拉第電磁感應(yīng)定律，安培定律和E-J冪指數(shù)關(guān)系式推演而出的。磁場H的分量Hx、Hy、Hz是待求解量。

對于描述某具體物理現(xiàn)象的偏微分方程，要想得到穩(wěn)定的定解，必須給出合適的邊界條件。對于控制方程(9), Dirichlet 和 Neumann 邊界條件用于描述這里的電磁邊界條件。在超導(dǎo)區(qū)Ω1與非超導(dǎo)抗磁區(qū)Ω2交界處的邊界應(yīng)該滿足連續(xù)邊界條件：

μ1H1n=μ2H2n。

(10)

計算區(qū)域的最外邊界Ω2 屬于動態(tài)邊界條件。這里使用一個時間相關(guān)的函數(shù)來描述：

HΩ2(r,t)=fΩ2(r,t)。

(11)

函數(shù)fΩ2(r,t)描述了高溫超導(dǎo)塊材外界的非均勻磁場在邊界Ω2處隨時間的變化而變化。

2 數(shù)值求解過程

本文中使用有限元方法(FEM)和有限差分方法來建立數(shù)值求解代碼，進(jìn)而對此高溫超導(dǎo)塊材電磁非線性強(qiáng)耦合問題進(jìn)行求解。

模型中使用了直接迭代方法來處理非線性偏微分方程(9)，如下所示：

(12)

其中：上標(biāo) i(i=1,2,…)暗示了迭代步驟； 下標(biāo) n 暗示了時間空間推演步驟。Hn-1是時間推演的第n步驟時的穩(wěn)定解；時間變分△t指明了時間推演步每次迭代時的時間步大小。當(dāng)前計算區(qū)域位于超導(dǎo)區(qū)域時，有效電導(dǎo)率σ和參數(shù)λ滿足：σ=σsc和λ=1。

方程(12)又可以變換為：

(13)

模型中使用四面體單元對計算區(qū)域進(jìn)行離散剖分處理。這些屬于外邊界上的單元結(jié)點上的值由于屬于邊界條件方程(11)描述，由方程(11)計算得出。代數(shù)方程組可由下面矩陣表達(dá)式描述：

(14)

其中：[A]和[C]分別是方程(14)第二項和第三項的系數(shù)矩陣。[I]是方程(14)第一項相關(guān)的單位矩陣。[F]是與Hn-1相關(guān)的負(fù)載列向量。很明顯，矩陣[A]和[C]是與有效電導(dǎo)率σsc緊密相關(guān)的。

數(shù)值求解的主要過程可以很簡明地由下面所描述：

步驟1: 最開始，對方程(14)初始化初值如下所示：

t=0:H=H0,HΩ2(r,t=0)=fΩ2(r,t=0)。

(15)

賦于超導(dǎo)塊材內(nèi)電導(dǎo)率σsc非常大的初始值，這意味著當(dāng)高溫超導(dǎo)塊材在初始迭代步驟時屏蔽電流密度幾乎等于0。

步驟2:t=t+dt,n=t,i=1。

步驟3: 計算剛度矩陣[A]、[C]和[F],求解列向量矩陣[Hi+1] 的值。

(16)

δ是一個給定的值，用于控制/橫量迭代精度。如果超導(dǎo)塊材中所有的剖分節(jié)點都滿足表達(dá)式(16)，則執(zhí)行式(16),then[Hi]=[Hi+1] ,i=i+1,gotoStep3。

步驟5: 如果n

當(dāng)計算出來磁場強(qiáng)度和電流密度分布后，高溫超導(dǎo)塊材所受到的磁懸浮力可以由下式計算得出：

Fem=∫vJ×Bexdv。

(17)

其中： v代表高溫超導(dǎo)塊材的體積；J是電流密度；Bex是高溫超導(dǎo)所受到的外磁感應(yīng)強(qiáng)度。求解代碼的有效性已由實驗驗證[6]。

圖2 高溫超導(dǎo)塊材與永磁軌道組成的磁懸浮系統(tǒng)

3 結(jié)果與討論

基于文中描述的數(shù)值求解器，研究了永磁軌道上方高溫超導(dǎo)塊材內(nèi)臨界電流密度對磁懸浮力的影響。由永磁軌道和高溫超導(dǎo)塊材組成的磁懸浮系統(tǒng)簡圖如圖2所示。其中永磁軌道是由橫截面為40mm寬、40mm高的NdFeBs構(gòu)成。永磁軌道的兩個NdFeB的磁化方向向?qū)ε帕?。永磁軌道的橫截面平行于x-y平面，其尺寸為寬90mm、高40mm。由永磁軌道感生出的磁場在x-y平面內(nèi)呈現(xiàn)對稱分布。

圖3 磁感應(yīng)強(qiáng)度在永磁軌道上方5 mm處豎直方向和水平方向上的分布圖

模型中永磁軌道感生的磁感應(yīng)強(qiáng)度可以使用等效電流平面解析模型計算。永磁軌道上方5mm處磁感應(yīng)強(qiáng)度的分布圖如圖3所示。

磁懸浮系統(tǒng)中的高溫超導(dǎo)塊材尺寸是邊長為30mm、厚為15mm。通過以下這種方法進(jìn)行模擬計算出樣品的磁懸浮力-位移恢復(fù)曲線:超導(dǎo)樣品在永磁軌道正上方沿豎直方向由位置A以1mm/s的速度向下向位置B移動。當(dāng)超導(dǎo)樣品到達(dá)位置B后，又以同樣的速度返回向位置A處移動。當(dāng)高溫超導(dǎo)塊材的磁感應(yīng)強(qiáng)度分布采用本文描述的方法計算出來后，塊材樣品內(nèi)部的電流密度分布就可以通過數(shù)學(xué)公式(1)計算得出。最后可以依照公式(17)將超導(dǎo)樣品的磁懸浮力計算出來。

超導(dǎo)樣品在不同臨界電流密度Jc情況下磁懸浮力-位移曲線圖如圖4所示。從圖中可以看出，超導(dǎo)樣品的磁懸浮力-位移磁滯恢復(fù)曲線展示了一些典型的磁滯特性。這是因為高溫超導(dǎo)體材料的磁滯特性決定的。當(dāng)超導(dǎo)塊材沿豎直方向向下移動到永磁軌道表面處，然后又返回向上離開永磁軌道，不同的臨界電流密度對應(yīng)的磁懸浮力-位移磁滯恢復(fù)曲線所圍的面積也不相同。從圖4可以看出，磁懸浮力-位移磁滯恢復(fù)曲線所圍的面積隨著臨界電流密度Jc的增加而減少。這是因為，當(dāng)超導(dǎo)塊材的臨界電流密度較低時，樣品受到外磁場的作用，磁場更容易穿透到超導(dǎo)塊材內(nèi)部。

為了更加清楚地顯示出超導(dǎo)樣品最大磁懸浮力與樣品的臨界電流密度Jc的關(guān)系，將最大磁懸浮力與臨界電流密度曲線繪制到同一圖中，如圖5所示。從圖中可以看出，磁懸浮力最大峰值隨著臨界電流密度從1.0×108A/m2到2.0×108A/m2逐漸增加時，呈線性迅速增加。當(dāng)臨界電流密度大于2.0×108A/m2后,磁懸浮力最大峰值隨著臨界電流密度的增加開始變的緩慢增加。

圖4 不同電流密度下磁懸浮力-位移曲線 圖5 不同臨界電流密度下超導(dǎo)塊材最大懸浮力曲線

通過對YBCO塊材磁懸浮力數(shù)值計算結(jié)果可以得出：當(dāng)增強(qiáng)超導(dǎo)塊材材料特性(通常提高臨界電流密度)時，可以在一定范圍內(nèi)提高超導(dǎo)塊材的最大磁懸浮力值。最大懸浮力的增加幅度受到外界應(yīng)用磁場的限制。

4 結(jié) 語

基于磁矢量方法(H-方法)和有限元方法，成功地在FORTRAN平臺開發(fā)出了一個3D-模型數(shù)值計算代碼，此數(shù)值計算模型著眼于高溫超導(dǎo)電磁場的模擬計算。使用此數(shù)值計算系統(tǒng)，研究了高溫超導(dǎo)塊材不同臨界電流密度對磁懸浮力的影響。數(shù)值模擬計算顯示了增強(qiáng)超導(dǎo)塊材內(nèi)部臨界電流密度可以提高最大磁懸浮力，但是最大磁懸浮力提高的幅度受到外磁場的限制。

[1] John R H.Applications of high-temperature superconductors in power technology[J].Phys,2003(66):1865-1886.

[2] Murakami M.Progress in applications of bulk high temperature superconductors[J].Phys C,2000,341:2281- 2284.

[3] Brandt E H.Superconductors of finite thickness in a perpendicular magnetic field:strips and slabs[J].Phys Rev B,1996,54(6):4246-4264.

[4] Prigozhin L.The Bean model in superconductivity:variational formulation and numerical solution[J].Journal of Computational Physics,1996,129:190-200.

[5] Qin M J,Li G,Liu H K,et al.Brandt,Calculation of the hysteretic force between a superconductor and a magnet[J].Phys Rev B,2002,66:024516.

[6] Lu Y Y,Wang J S,Wang S Y,et al.3D-Modeling Numerical solutions of electromagnetic behavior of HTSC bulk above permanent magnetic guideway[J].J Supercond Nov Magn,2008,21:467-472.

Effect of Superconductor with Critical Current Above PMR on Levitating Force Based on 3D-modeling

LU Yiyun

(Luoyang Institute of Science and Technology, Luoyang 471023, China)

In this paper, the electromagnetic field simulation of the high-temperature superconductors (HTS) bulk was carried out and the magnetic force between the permanent magnet railway (PMR) was calculated. Based on finite element method, a 3D-modeling numerical calculation method is described. The models were formulated with the magnetic field vector (H-method). The resulting codes were developed with FORTRAN platform. The electric field intensity E and the current density J constitutive relation of HTS were described with E-J power law. Kim macro-model was used to describe critical current density Jc of HTS bulk. Two virtual HTS bulks were used to solve the critical current density Jc anisotropic properties of HTS materials. A superconducting levitation system composed of one HTS bulk and PMR was successfully investigated using the proposed method. By this method, the influence of critical current density on magnetic levitation force of the superconducting levitation system was mathematically studied.

superconductor; HTSC; 3D model; simulation

2017-01-06

蘆逸云(1972-)，男，新疆溫泉縣人，博士，副教授，主要從事超導(dǎo)電磁理論及應(yīng)用方面的研究。

國家青年自然基金項目(11205080).

10.3969/i.issn.1674-5403.2017.02.022

O441;TM12

A

1674-5403(2017)02-0083-05