二代高温超导REBCO(稀土-钡-铜-氧)带材具有良好的载流能力及较高的机械强度，在核磁共振成像和回旋加速器等领域得到了广泛应用[1-2]．2010年，Hahn等[3]提出了高温超导线圈的无绝缘绕制技术，匝间绝缘层的移除使得电流可以沿线圈径向分流，降低了线圈失超期间的温升，进而提高线圈的热稳定性．这些优点使得由REBCO带材绕制的无绝缘线圈成为制造极高场磁体的有力候选者．2019年，美国国家高磁场实验室将REBCO带材绕制的无绝缘磁体内插于电阻式磁体中，产生了45.5 T的磁场[4]；2020年，中国科学院电工所使用内插无绝缘REBCO线圈和低温超导线圈实现了32.35 T的全超导磁体[5]．以往研究很少考虑力学应变对超导线圈电磁性能的影响，使得高背景磁场下的数值模拟结果和实验测量结果存在误差．Kolb-Bond等[6]将带材受电磁体力后产生的偏转角引入了线圈的电磁场分析；Yan等[7-8]基于此并考虑了应变对临界电流密度的影响，通过有限元T-A方法分析了超导线圈的磁化过程，数值模拟结果与实验结果符合较好．通常模拟无绝缘线圈的充电过程须基于基尔霍夫定律构建等效电路模型，利用得到的环向电流求解电磁场．Mataira等[9]提出一种径向施加电流的方式，可直接通过有限元方法模拟无绝缘线圈的充电过程．本研究基于该方法构建REBCO无绝缘线圈的二维轴对称电-磁-力耦合模型，并与常规的电路模型结果进行对比验证结果的可靠性；最后模拟无绝缘线圈在较高背景磁场下磁化和充电过程中的电磁及力学行为，分析线圈变形对临界电流密度退化的影响．1 电磁场控制方程1.1　棱边单元不同于常见的有限元节点型单元，棱边单元将自由度赋予单元的各边[10]．如图1所示的矩形棱边单元e，x和y方向的长度分别为ae和be，(xce,yce)是单元中心的坐标，Qke(k=1,2,3,4)是矢量场Q沿单元各边的切向分量，于是单元内部的矢量场Qe为Qe=∑k=14QkeNke , (1)式中Nke为插值形函数，N1e=(yce+be/2-y)ex /be;N2e=(y-yce+be/2)ex/be;N3e=(xce+ae/2-x)ey /ae;N4e=(x-xce+ae/2)ey /ae, (2)其中ex和ey分别为x和y方向的基矢量，可以验证插值形函数Nke有如下性质∇∙Nke=0．(3)式(3)表明：用棱边单元表达无源矢量场时，场函数散度为零将自动满足，无须添加额外的约束方程．10.13245/j.hust.239302.F001图1矩形棱边单元e1.2　H方法超导电磁场计算需要的麦克斯韦方程组如下∇×Ε=-∂B/∂t; ∇×H=J;∇∙B=0 , (4)式中t为时间．超导的电磁本构方程为E=ρJ;B=μ0μrH, (5)式中：E，B，H和J分别为电场强度、磁感应强度、磁场强度和电流密度；ρ为电阻率；μ0和μr分别为真空磁导率和相对磁导率，通常取μr=1．联立式(4)和(5)可得H方法的控制方程为[11]∇×(ρ ∇×H)+μ0∂H/∂t=0．(6)设权函数为f，由伽辽金加权残值法及分部积分可得式(6)在每个单元e上的等效积分弱形式为      ∫Ωeρ (∇×H)∙(∇×f)+μ0∂H∂t∙f dΩe+∮Γeρ ne× (∇×H)∙f  dΓe=0  , (7)式中Ωe，Γe和ne分别为单元e的面元矢量、边界线元矢量和边界单位外法线矢量．在伽辽金有限元方法中，权函数即为单元插值基函数，若使用棱边单元并考虑矢量混合积的性质，则有      ∑lHle∫Ωeρ (∇×Nke)∙(∇×Nle)  dΩe+∑l∂Hle∂t∫Ωeμ0 Nke∙Nle  dΩe-∮ΓeNke×E∙ne  dΓe=0  . (8)在每个单元e上，令Kkle=∫Ωeρ (∇×Nke)∙(∇×Nle)  dΩe ;Lkle=∫Ωeμ0 Nke∙Nle  dΩe;Pkle=∮ΓeNke×E∙ne  dΓe, (9)以Ke,Le和Pe分别表示单元e上由Kkle,Lkle和Pkle得到的单元矩阵，进而组装整体矩阵，即K=∑eKe;L=∑eLe;P=∑ePe, (10)可得到求解域离散后的方程组为KH̃+L∂H̃∂t-P=0. (11)式中H̃为所有单元各边的切向磁场强度组成的列向量．磁通量守恒可由棱边单元形函数的性质式(3)保证．2 二维轴对称电-磁-力耦合模型2.1　二维轴对称电磁模型2.1.1　各向异性电阻率真实的超导线圈为图2(a)所示的螺旋结构．以带材的切线方向(eτ)、表面法线方向(en)和高度方向(ez)建立局部坐标系，线圈的电阻率为如下二阶张量ρlocal=ρn000ρτ000ρz, (12)式中ρn，ρτ和ρz分别为径向、切向和轴向电阻率，则局部坐标系中的本构方程为Elocal=ρlocalJlocal, (13)式中Elocal和Jlocal分别为局部坐标系下的电场强度和电流密度．10.13245/j.hust.239302.F002图2线圈结构、轴对称模型及局部/全局柱坐标系变换在轴对称模型中，线圈被视作多匝同心圆环，建立全局柱坐标系，如图2(b)所示．局部坐标系和全局柱坐标系的变换矩阵为g，即[12]g=cosα-sinα0sinαcosα0001 , (14)同一位置处eτ和eθ的夹角α=±d/(2πr)，其中：d为带材厚度；r为该位置的半径．若线圈的缠绕方向如图2(b)所示，则取正号，否则取负号．联立式(13)和(14)可得全局柱坐标系中的本构方程为[9]Eglobal=gρlocalJlocal=gρlocalg-1Jglobal , (15)式中：Eglobal为全局柱坐标系下的电场强度；定义gρlocalg-1为全局柱坐标系下的电阻率张量．对式(12)中各分量做出如下说明：ρτ和ρz均等于超导电阻率ρsc，即[13]ρτ=ρz=ρsc=EcJc(B)JJc(B)n-1, (16)式中：临界电场Ec=1×10-4 V/m[14]；J为电流密度的模；与磁场有关的临界电流密度Jc(B)是平行于和垂直于带材的磁感应强度分量B∥和B⊥的函数[15]JcB=Ic0wd11+k2B∥2+B⊥2/Bch , (17)其中，Ic0为带材的自场临界电流[15]；w为带材宽度．ρn=ρct/d为等效体接触电阻率[9]，ρct=70 μΩ∙cm2[16-17]．2.1.2　二维有限元法计算无绝缘线圈的电磁场高温超导无绝缘线圈二维轴对称电磁模型的几何结构如图3所示．电流由正极流出，沿径向从线圈最内匝进入，在最外匝流出，沿电路返回电源负极[9]．方程(6)在全局柱坐标系下求解，电磁场矢量的各分量均参与计算．由于模型是轴对称的，因此各分量的环向导数为零．10.13245/j.hust.239302.F003图3二维轴对称线圈电磁模型的几何结构 (黄色箭头为rz平面内电流方向)2.2　二维轴对称力学模型2.2.1　REBCO带材的均匀化REBCO带材由铜、银、YBCO(钇-钡-铜-氧)及Hastelloy基底等组分组成．为求得宏观上的平均力学响应，须将结构不均匀的REBCO带材等效为一种均匀的材料，因此要构建代表性体积单元(RVE)[18]．设1，2和3分别代表空间的三个方向，不失一般性，将RVE构建成几何尺寸a1×a2×a3的长方体，内部结构按照REBCO带材各组分的真实顺序及厚度的占比排布，代表性体积单元RVE(见图4)．10.13245/j.hust.239302.F004图4REBCO带材的代表性体积单元RVE令σxy和εxyx,y=1,2,3分别表示应力分量和应变分量．在有限元计算中，设RVE的6个应变分量为1，同时其余分量均为0[19]，可得到每种情形下的应力分量，将结果代入广义胡克定律σ=Cε(σ和ε分别为应力和应变张量)，可以获得刚度矩阵C．而柔度矩阵为刚度矩阵的逆，对于REBCO带材这类横观各向同性材料，表达式为[18]S=C-1=1Y1-υ12Y2-υ13Y3000-υ21Y11Y2-υ23Y3000-υ31Y1-υ32Y21Y30000001G230000001G310000001G12,式中：Yk为弹性模量；υkl为泊松比；Gkl为剪切模量．由此可获得REBCO带材均匀化后的Yk，υkl和Gkl．在轴对称模型中，1，2和3分别对应r，θ和z方向．2.2.2　力学控制方程假设线圈的力学响应为准静态，因此不考虑惯性项的平衡方程为∇∙σ+F=0,式中F为线圈受到的电磁体力，在轴对称模型中为F=FrFz=JθBz-JzBθJrBθ-JθBr,式中：Jr，Jθ和Jz分别为径向、环向和轴向电流密度；Br，Bθ和Bz分别为径向、环向和轴向磁感应强度．边界条件为uzbottom=0，其中uzbottom为线圈及芯轴下沿的轴向位移．力学模型考虑了带材之间的接触特性，不考虑线圈中的匝间摩擦[20]．2.3　力学变形对临界电流密度的影响2.3.1　偏转角的影响式(17)的磁感应强度垂直分量B⊥和平行分量B∥即为径向分量Br和轴向分量Bz．当背景磁场较大时，带材将发生明显的变形，带材变形后磁感应强度的垂直分量与平行分量如图5所示．此时B⊥和B∥满足B⊥B∥=cos β-sin βsin βcos βBrBz ,式中：β=arctan(∂ur/∂z)为带材的偏转角；ur为径向位移[6-7]．10.13245/j.hust.239302.F005图5带材变形后磁感应强度的垂直分量与平行分量2.3.2　应变导致临界电流退化的影响临界电流密度也与应变相关，此时Jc=p(ε)Jc(B)，p(ε)为与应变有关的归一化临界电流密度．Barth等[21]对4.2 K下p(ε)进行了测量，其中0.67%是不可逆应变，拟合得到p(ε)的表达式为[22]p(ε)=1-6.28×10-8exp(2174.35ε)  (ε≤0.67%);0.867-582.466(ε-0.0067)  (ε0.67%).3 模型验证3.1　和等效电路模型的对比图6为等效轴对称电路模型的示意图[23]．图中：Mul为第u匝和第l匝之间的互感，当u=l时表示第u匝的自感；第u匝中超导层电阻Rscu和其他层的电阻Rnormu并联组成第u匝的环向电阻Rθu；Rctu为第u匝的接触面电阻；电源电流Is进入第u匝后分流为环向电流Iθu和径向电流Iru．10.13245/j.hust.239302.F006图6等效轴对称电路模型若线圈未发生失超，则环向电流均在超导层中流动，不计超导层电阻，因此线圈环向仅存在感应电动势．由基尔霍夫定律可得[24]Iθu+Iru=Is;∑u=1mMuldIθudt-RctuIru=0,式中：m为线圈总匝数；接触电阻Rctu的表达式为Rctu=ρctw(2πru)   (u=1,2,…,m-1),Rctu-1   (u=m),其中ru为第u匝的半径．通过分析一个60匝的无绝缘REBCO线圈的充电过程，对本文方法和常规的等效电路法进行比较．本文结果由有限元软件COMSOL计算得到．带材和线圈的参数如下[14,15,25-26]：带材宽度为4 mm，厚度为0.1 mm；铜层厚度为20 μm，银层厚度为20 μm，YBCO层厚度为1 μm，基底厚度为50 μm；式(16)中n=25，式(17)中k=0.038 13，h=0.712 2，Bc=631 mT；4.2 K带材自场临界电流为720 A；Yr=154 GPa，Yθ=164 GPa，Yz=164 GPa；Grθ=57 GPa，Gθz=62 GPa，Grz=57 GPa；υrθ=0.333，υθz=0.328，υrz=0.333；线圈匝数为60，内半径为50 mm，外半径为56 mm，带材总长度为19.98 m，线圈自感系数为0.78 mH，自场线圈磁场常数为0.71 mT/A，4.2 K线圈自场临界电流为419.1 A．电源充电速率为2 A/s，目标电流为120 A．未施加背景磁场，此时力学变形对电磁场的影响可以忽略．图7为两种模型线圈充电过程中轴向中心磁场的变化；图8为电路模型得到的环向电流密度Jθ和H方法得到的切向电流密度Jτ．可以看出：H方法的计算结果与电路模型的结果基本一致．10.13245/j.hust.239302.F007图7两种模型线圈充电过程中轴向中心磁场的变化10.13245/j.hust.239302.F008图8电路模型得到的环向电流密度Jθ和H方法得到的切向电流密度Jτ (色标单位：109 A/m2)3.2　和实验数据的对比图9为实验测量与H方法得到的线圈中心磁场，实验测量结果由文献[27]给出．实验中，电源充电速率为5 A/s，目标电流为30 A．由图9可见，模拟结果与实验测量数据符合较好．10.13245/j.hust.239302.F009图9实验测量与H方法得到的线圈中心磁场综合对比H方法与电路模型以及实验测量得到的结果，可以证明H方法准确可靠．4 无绝缘线圈的磁化与充电将前文所述线圈放置在沿轴向的均匀背景磁场Bbg中，线圈在Bbg=15 T下的临界电流为371 A．芯轴的弹性模量和泊松比分别为35.4 GPa和0.21[28]．由于背景磁场较高，因此力学变形对临界电流密度的影响必须考虑．4.1　磁化与充电过程中线圈的电磁及力学行为第300，301和450 s分别是磁化结束、充电前和充电结束的3个时刻．图10为不同时刻线圈的切向电流密度(Jτ)、偏转角(β)、环向应变(εθ)和环向应力(σθ)的分布．在磁化过程中，顺时针的切向电流密度由线圈的上下两端往内部穿透，下端向外偏转而上端环向受压向内偏转，线圈整体环向受压，压缩应变最大点位于靠外匝的高度中心．当背景磁场开始稳定而充电尚未开始时，逆时针的切向电流密度开始在线圈上下两端出现，这是为了保证在一定时间后，线圈内部的净电流为零．然而，由于超导体内电阻率非常低，因此电流密度会长期存在，这使得线圈在充电开始前依然受到电磁体力，并且偏转角、应变和应力均不为零．充电开始后，线圈上下两端的逆时针切向电流密度继续增大并向内部穿透．充电结束时，最大偏转角为1.38°，最大环向应力达到324 MPa，最大环向应变为0.196%．10.13245/j.hust.239302.F010图10不同时刻切向电流密度Jτ、偏转角β、环向应变εθ和环向应力σθ的分布4.2　与顺序求解结果的对比图11是不考虑力学变形影响的模型450 s时的偏转角和环向应变．由图11可见：最大偏转角为1.29°，最大环向应变为0.193%．这意味着在考虑的外场及假设条件下，若仅将应变和偏转角等作为线圈在电磁载荷下的力学响应，则顺序求解与双向耦合得到的结果差异并不大．然而，考虑力学变形与否对电磁场的计算结果有较为显著的影响．10.13245/j.hust.239302.F011图11不考虑力学变形影响的模型450 s时的偏转角和环向应变图12为充电阶段不考虑和考虑力学变形影响的临界电流密度的分布．不考虑力学变形采用先算电磁场，再算电磁体力作用下力学响应的顺序求解方式；考虑力学变形则意味着在电磁场计算中加入变形的影响，实现电磁场-力学的双向耦合．与顺序求解相比，双向耦合的临界电流密度依然具有上下两个临界电流密度较小的区域，但不具有如图12(a)中的对称特征．由于双向耦合中环向应变较小，即归一化临界电流密度p(ε)几乎等于1，因此造成图12中(a)和(b)差异的主要原因是偏转角不同．10.13245/j.hust.239302.F012图12充电阶段不考虑和考虑力学变形影响的临界电流密度分布 (色标单位：109A/m2)在双向耦合模型中，由于电磁力的作用，线圈每匝带材的上端向磁感线的方向偏转，使得线圈最顶端的临界电流密度受垂直场B⊥的影响较小，衰减不太明显．双向耦合模型上端临界电流密度较小区域的面积大于顺序求解模型，而最小值较之于顺序求解模型更小．因此，线圈在高背景磁场下充电时，靠近上端的部位有较大的失超风险，而顺序求解模型低估了该位置的失超风险．可以通过在线圈外侧加固不锈钢套等手段抑制带材的偏转角，进而延缓线圈上端临界电流密度的衰减，保障线圈在高背景磁场下运行的安全性．另外，两种情况得到的屏蔽磁场Bsc也有区别.图13是双向耦合和顺序求解下的屏蔽磁场，充电开始时两种结果非常接近，但在到达峰值后，双向耦合结果快速衰减，而顺序求解结果的衰减速率小于前者，即不考虑力学变形影响会高估屏蔽磁场，这说明在线圈的电磁场分析中考虑力学影响尤为重要．10.13245/j.hust.239302.F013图13双向耦合和顺序求解下的屏蔽磁场Bsc5 结语给出了计算超导线圈电磁场分布的H方法控制方程及其有限元离散格式．基于该方法，构建了高温超导无绝缘线圈二维电-磁-力耦合模型，该模型从带材角度偏转和临界电流退化两个方面考虑了线圈受电磁体力后变形对临界电流密度的影响．将线圈零背景磁场充电的模拟结果与常规的电路模型结果进行了对比，验证了模型的可靠性．分析了高温超导无绝缘线圈在均匀高背景磁场下的磁化和充电过程．结果表明：磁化和充电过程中电流密度、偏转角、环向应变和环向应力的方向均相反；背景磁场稳定而充电尚未开始之前线圈内已出现和磁化电流反向的感应电流．通过与不考虑力学变形影响的计算结果对比可知：在高背景磁场下充电时，考虑变形影响的临界电流密度分布不再具有上下对称性，且靠近上端的区域临界电流密度较小，有较大的失超风险；不考虑变形影响会高估线圈的屏蔽磁场．接下来，须综合考虑双向耦合条件下力学变形对平衡方程的影响，为实际工程应用提供参考．