周艷秋 余志強
(1.大連科技學院電氣工程學院,116052,大連;2.石家莊鐵道大學電氣與電子工程學院,050043,石家莊//第一作者,副教授)
近年來,隨著新型超導材料的成功研制,高溫超導技術得到了快速發展。其中,高溫超導磁懸浮具有無源自穩定的優良懸浮特性,其相關應用已經成為國內外研究的熱點之一[1-3]。高溫超導軸承是其中的重要代表(以下簡稱“超導軸承”)之一[4],由超導定子和永磁轉子兩部分組成,如圖1所示。在設計和優化超導軸承結構時,采用數字方法分析其懸浮特性是非常必要的。
本文采用H法(磁場強度法)在二維軸對稱空間建立了徑向型超導軸承的有限元數學模型。超導體E-J關系(電場強度與電流密度的關系)采用冪指數模型。通過試驗測量值與數值計算結果的比較,驗證了本模型的正確性。在此基礎上,應用本模型分析了徑向型超導軸承的電磁特性。本文提出的數學模型不僅能夠分析大功率負載的徑向型超導軸承的電磁特性,而且對于內轉子和外轉子兩種類型的超導軸承均適用。
本文試驗用的徑向型超導軸承裝置2D(二維)軸對稱模型原理圖如圖1左圖所示。原理圖中參數含義及其數值見表1。
如圖1 a)所示,徑向型超導軸承裝置的超導定子由16個超導塊材組成的2層超導環組成,永磁轉子由3個永磁環和4個聚磁鐵環組成。所建立的2D軸對稱模型結構和尺寸與試驗裝置完全相同。

a) 結構圖

參數參數含義參數數值/mmROP永磁轉子外徑18.0RIP永磁轉子內徑10.0ROS超導定子外徑32.0RIS超導定子內徑21.0WHTS超導塊材寬度10.0hHTS超導塊材高度16.0WPM永磁環寬度8.0hPM永磁環高度8.0hF永磁環軸向間高度2.0g永磁轉子與超導間氣隙3.0g1,g2定子之間氣隙1.5gb 超導塊材軸向氣隙1.0
如圖1 b)所示,在二維軸對稱空間里,向量H具有兩個自由度:Hr和Hz。電流密度J和電場強度E僅在φ方向具有非零分量:Jφ和Eφ。根據Ampère定律,電流密度J可表示為:
(1)
式中:
z——二維坐標縱軸坐標值;
r——二維坐標橫軸坐標值;
超導體E-J關系采用冪指數模型形式:
(2)
式中:
Esc——超導塊電場;
E0——常量,取值為1×10-4V/m;
Jc——臨界電流密度;
Jsc——感應電流密度;
ni——定義為U0/kT,其中U0和k分別為超導體的釘扎勢能和Boltzmann常量。
采用Kim模型[5]描述外場強度對臨界電流的影響。引入超導體的等效電阻率(effective resistivity)ρsc計算超導體的內部磁場,即其定義為:
(3)
將本構關系E=ρJ和Ampère定律(1)式代入到法拉第電磁感應定律▽×E=-μ?H/?t中,得到關于向量H的拋物型偏微分方程:
(4)
式中:
μ——磁導率;
ρ——電阻率;
▽——梯度算子;
t——時間。
(4)式即為本模型的電磁場控制方程。該方程較為簡潔,體現了H法的優越性。
本節采用伽遼金法和格林公式對上面得到的電磁場控制方程(4)式進行推導,得到其弱形式為:
(5)
式中:
S——整個求解域;
Г——S的邊界;
?H/?n——邊界上法向微分,n為S表面法向單位矢量;
W——向量H的權函數。
對計算域進行網格劃分,選擇一階三角形單元為試探函數,則在一個網格單元里,向量H的表達式為H=Hj,e(t)Nj,j=1,2,3,e表示單元;根據伽遼金法,權函數與試探函數相同即W=Ni,i=1,2,3;式中i,j是單元節點的編號。這樣,在網格單元的一個節點上的弱形式方程可以寫為:
(6)
在時域中,采用后向Euler法對時間導數進行離散化處理,將每個單元質量矩陣和單元剛度矩陣進行擴充,并把擴充后的矩陣中具有相同下標的元素進行累加,將單元列向量擴充為含所有節點向量的列向量。這樣,就得到了該非線性系統的最終的有限元方程,其形式如下:
([M(μ)]+[S(ρ)]·Δt){Hc}=
[M(μ))]{Hc-1}
(7)
式中:
Δt——連續時間步的時間間隔;
c、c-1——分別代表當前時間步和上一個時間步。
對應于單元矩陣,[M(μ)]和[S(ρ)]分別為非線性系統的質量矩陣和剛度矩陣。列向量{Hc-1}在本時間步為已知,放在方程的右側。
徑向型超導軸承的懸浮力F,是永磁轉子的磁場作用到超導定子中超導塊材的感應電流所產生的電磁力在軸向的分量。因為超導塊材中的感應電流密度Jsc為φ方向,所以產生懸浮力的外磁場方向為徑向,即磁場的徑向分量Br。根據Lorentz方程和有限元法的特性,懸浮力F的公式可表示為:
(8)
式中:
Br——磁場的徑向分量;
ΔSsc,e——超導域網格單元的面積均值;
Ng——網格單元的數量;



數值計算的計算流程如圖2所示。
本節應用所提出的數值計算方法,計算徑向型超導軸承模型的懸浮力并與測量結果進行比較,然后分析其電磁行為。在計算的過程中,超導定子保持靜止,冷卻溫度為77.3 K。永磁轉子的移動速度為1 mm/s,移動范圍為0~10 mm。在超導體的冪指數模型中,ni通過U0/kT計算,其值為15。
圖3所示為徑向型超導軸承懸浮力的計算結果,從圖中可以看出隨著永磁轉子的軸向移動,懸浮力呈現出非線性變化。懸浮力上升段曲線可以分成兩部分:線性上升段(對應的移動距離為0~5 mm)和飽和段(對應的移動距離為5~10 mm)。在線性上升段,懸浮力幾乎線性上升,臨界電流密度Jc高的曲線更接近測量曲線。在飽和段,懸浮力曲線的增加率逐漸變慢,在懸浮力達到最大值后變成負值。懸浮力的這些非線性特征主要是由于場冷條件下的磁通釘扎效應造成的。也就是在線性上升段,超導塊材的捕獲磁通的數量幾乎保持不變,感應電流較大,外磁場較強,所以懸浮力隨著永磁轉子的移動迅速增加。而在飽和段,隨著永磁轉子的繼續移動,超導塊材受到的外磁場變弱,導致了懸浮力增加緩慢進而減少。就懸浮力曲線的形狀而言,計算曲線反映了以上兩部分的特性。

圖2 計算流程圖

圖3 懸浮力數值計算結果
圖3亦清楚地表明了計算曲線的變化趨勢與測量曲線是一致的。隨著Jc的增加,計算曲線的最高點亦上升。由于試驗中的超導塊材為單晶熔融織構的YBCO材料,所以計算中Jc分別設置為6.5×107A/m2、8.0×107A/m2以及9.0×107A/m2,這些值均在材料可達到的正常范圍內。當Jc=6.5×107A/m2時,計算曲線在線性上升部分的最大誤差為9.2%,飽和部分的最大誤差僅為3.7%,與測量結果具有較好的一致性。
圖4為永磁轉子沿軸向移動的6個位移處的計算域中磁場的分布情形。由圖4可以看出,大量的磁力線集中分布在聚磁鐵環附近,一小部分磁通線滲入到超導塊材的邊緣然后穿出超導塊材返回到相鄰的磁極附近。由于超導體的強抗磁性,只有非常少的磁通線進入到超導塊材的內部。因此,大部分磁力線被限制在了2個小的氣隙中:永磁轉子和超導定子之間的氣隙(圖1b中標記為g)和2個超導環(圖1b中標記為超導塊材)之間的軸向氣隙(圖1b中標記為gb)。從圖5可以看出,超導塊材中感應電流主要分布在其邊緣部分,因此在超導塊材的內邊緣部分存在較大的電磁力。






圖4 不同位置時計算域的磁場分量Br分布






圖5 不同位置時超導塊材中感應電流密度Jsc的分布
為了定量說明超導塊材中的磁場和電流密度的變化情況,在圖1b)下方超導塊材的橫截面上選取了A—G 7個點,其位置如圖6示。

注:g3——A—G點各點到超導塊材邊緣的距離
下文定量地說明這7個點的磁場和感應電流密度。圖7~11分別為磁場和感應電流密度隨永磁轉子軸向移動時的數值變化圖。

圖7 永磁轉子軸向移動時A—G點磁場徑向分量Br

圖8 永磁轉子軸向移動時A—G點的磁場縱向分量Bz

圖9 永磁轉子軸向移動時A、B、C點感應電流密度

圖10 永磁轉子軸向移動時D點感應電流密度

圖11 永磁轉子軸向移動時E、F、G點感應電流密度
設置每個點到離它最近的超導塊材邊線的距離為1 mm。對于A、D、E點,它們磁場的幅值相對較高,如圖7~8所示。這是因為這些點距離永磁轉子最近,滲入進來的磁通線較多的緣故。盡管A點和E點距離超導塊材的兩條邊都是最近的,但是A點同時離氣隙gb很近,其磁場幅值是最大的。從數值上看,A點Br和Bz的幅值分別為0.107 6T和0.099 8T。圖7亦清楚的表明,對于B、C、E、F點,其磁場幅值隨永磁轉子的移動而減少。另一方面,由于磁通釘扎效應和磁場的連續性,在永磁轉子移動的過程中,每個點處磁場的變化曲線為光滑的波浪線,變化較為緩慢。
從圖9~11可看出:G點的感應電流密度的變化量是最小的,其范圍為-8.118 6×106~7.245 4×106A/m2,這是因為有較少的磁通線滲入到G點,磁場變化較小。注意到由于較多的磁通線滲入到A點和E點,使得它們的磁場幅值較大,但它們的感應電流密度并不是較大的。也就是說,較高幅值的磁場并不能產生較大的感應電流,即感應電流隨磁場變化的增大而增加。最大的電流密度出現在C點,其值為3.789×107A/m2。圖5中C點處的較深顏色也能定性地反映這個結論。
通過以上分析表明,本文提出的有限元數學模型能夠用于分析包含多層超導環和多層永磁環的徑向型超導軸承的電磁行為,計算其懸浮力。
本文提出的有限元數學模型能夠用于分析徑向型超導軸承的電磁行為,計算其懸浮力。該模型的優點在于它的應用性強,不僅可以分析包含多層超導環和多層永磁環的徑向型超導軸承,而且還能用于分析內轉子型和外轉子型兩種結構。應用超導軸承的飛輪儲能系統已經應用到地鐵站的電力調峰[3],必將對其產生深遠影響。