基于socket通信运行方式的等离子体平衡反演计算

11超导托卡马克HT-7U的等离子体平衡反演技术原理等离子体平衡反演技术[1]主要是指通过对于圆形或非圆形等离子体放电外部的磁场和磁通的测量，来确定平衡磁面的技术。该算法的任务是计算R、Z平面上极向磁通的分布和设想的电流密度tJ分布，使得该数据与诊断数据以最小二乘法拟合，同时满足Grad-Shafranov方程),(0*RRJtp。总的极向磁通为：coilp，p是等离子体电流产生的磁通，coil是除了等离子体外的所有电流源产生的磁通。平衡解由覆盖整个真空室的矩形网格上的极向磁通值和等离子体电流密度tJ值组成。2超导托卡马克HT-7U的等离子体平衡反演系统算法超导托卡马克HT-7U的等离子体平衡反演系统主要完成两种模式的计算功能：平衡计算和反演计算。所谓固定边界平衡计算就是在给定一个初始等离子体边界的条件下，通过计算图1平衡计算流程图算流程图图2反演计算流程图开始计算等离子体电流密度计算极向场线圈电流计算等离子体网格区域的极向磁通值计算确定等离子体边界和磁轴位置及磁通值判断迭代前后计算的极向磁通值是否满足给定收敛条件结束是否开始计算等离子体电流密度计算极向场线圈电流计算等离子体网格区域的极向磁通值计算确定等离子体边界和磁轴位置及磁通值判断是否满足最小二乘意义上的拟合要求结束是否计算反演参数极向场线圈电流以及多项式系数迭代次数等于1？是否2极向场线圈电流和相应的等离子体电流密度分布，得到满足Grad-Shafranov平衡方程的极向磁通值，使得得到的等离子体边界与给定的等离子体边界相吻合。所谓反演计算就是在符合最小二乘意义上计算值与诊断数据的拟合计算，也就是使得计算得到的极向场电流和等离子体电流分布最大程度上满足磁测量值。其中平衡计算模式和反演计算模式算法流程分别如图1和图2所示。2.1计算等离子体电流密度等离子体平衡反演系统计算等离子体电流密度选取了多项式模型［2］。多项式模型如下：nmtFRPRJ,~1,~)~~(),~(10mimiimP，)~~(),~(10niniinF这里P是等离子体压力，FF是和极向电流有关的量，F表示对于极向磁通的偏导，m和n则是该模型中的自由参数。)/()(~axisbdyaxis是相对于放电边界磁通和磁轴磁通差的归一化磁通，axis是位于磁轴的磁通，bdy是位于最后一个闭合磁通面的磁通。2.2计算极向场线圈电流等离子体平衡反演系统计算极向场线圈的算法模型如下：),(),,,(),,,(),(1mmtmmmmmninininimZRJZRZRGdZdRIZRZRGZRM),(ZRM为给定或迭代计算得到的磁通值，),,,(ninimZRZRG表示电流源与等离子体格点的互感系数，),,,(mmmZRZRG表示等离子体之间的互感系数，),(mmtZRJ表示等离子体电流密度。此算法主要通过解矩阵方程AX=B得到极向场电流值。A对应对应算法模型中的),,,(1nininimZRZRG；B对应算法模型中的),(),,,(),(mmtmmmmmZRJZRZRGdZdRZRM。具体实现步骤如下：（1）分别计算极向场线圈与固定边界的互感系数和等离子体格点与固定边界点的感应系数（2）设置极向场线圈与固定边界点的互感系数作为A矩阵（3）设置单位矩阵ut（4）用单值分解法求A的逆矩阵A－1（5）设置给定或迭代计算得到的总极向磁通值与等离子体的极向磁通贡献部分之差作为B向量3（6）解得X=A-1B，即可得到极向场线圈电流2.3计算网域内的极向磁通超导托卡马克HT-7U的等离子体平衡反演系统是基于Grad-Shafranov平衡方程TJ0*的。Δ*算子可以写成以下形式：其中，),(1),(),;,(nnnnTFRPRRJ，P和F分别是等离子体压强和角向电流。这在计算等离子体电流密度算法中已经得到。Grad-Shafranov方程的右端与等离子体电流有关，若它是关于的线性函数或与无关，可精确求解方程。在给定了P和F的具体函数形式后，加上合适的边界条件，就可以从Grad-Shafranov方程解出),(ZR。对于一个矩形的计算区域，内部点的可以通过转化*运算符，利用fastbuneman方法和picard迭代方案，选择性地得到结果。具体实现步骤如下：（1）首先保存电流迭代更新前的总磁通值（2）计算网格区域边界处的磁通（3）计算网格区域中间处的磁通tRJ（4）通过buneman方法转化*，得到等离子体的磁通贡献。（5）通过格林函数法得到极向场线圈电流磁通贡献（6）最后，得到总的极向场磁通值2.4确定等离子体边界和磁轴此算法主要寻找等离子体边界、磁轴和确定等离子体边界上的x-point点，同时，计算归一化磁通函数~。2.4.1寻找等离子体边界和磁轴算法（1）调用双三样条插值函数通过插值得到格点之间的值。（2）给定一个初始的等离子体边界r坐标rad，rad=(rin+rout)/2。其中，rin=xctr(假定的边界中心r坐标)，rout=xlmin(限制器r最大或最小坐标)，得到对应的i，j，kk，计算边界磁通值。特别说明：边界psi值的计算过程如下。jiijijijijijiijijijiJrzrrr,21,,1,,1,12,1,,12)(221)(2x1x2y1y2Ψ(kk)Ψ(kk+nh)Ψ(kk+nh+1)Ψ(kk+1)Xta3a4●●21●●●●图3计算边界磁通值4从图3可以看出，若(xt,yt)在1点处时，即yt=y(j)时dxixxtkknhkkkk)()]()([)(若(xt,yt)在2点处时dxixxtkknhkkkk)()]1()1([)1(（3）根据i，j，kk得到四个相邻格点围成一个矩形区域，确定边界点相对这四个格点的位置。①如果该边界点在格点上，则取新的i，j，kk，得到新的四个格点围成新的矩形区域，然后检查区域内有没有限制器点。②如果该边界点不在格点上，直接检查区域内有没有限制器点。（4）根据限制器上的磁通值值psilx和边界点的磁通值psivl比较，分以下两种情况：①若psilxpsivl，则检查边界点(xt,yt)在不在限制器内，然后执行extrap。②若psilx=psivl，则直接执行extrap。（5）找到两个点(xt1,yt1)和(xt2,yt2)供选择作为下一个边界点。选取与前一个边界点距离较远的那个点作为下一个边界点。（6）取新的i，j，kk，判断当前边界点是否回到起点。若没有，则回到第（3）步继续找下一个边界点。确定等离子体边界的条件如下：①dxyyxxnn1.0)()(2121②0001.0bprevb③dyyyn5.01且dxxxn5.01特殊情况：①psilxpsivl且（xt,yt）不在限制器内。②psi误差大于给定的etol=1.E-04。这时，确定新的边界起点rad（通过二分法逐步找到真正的等离子体边界起点）和新的psi值，从第（2）步开始重新寻找边界点。无法确定等离子体边界的情况：①边界点数超过500还没有满足确定等离子体边界的条件。②dyyyn5.01或dxxxn5.01。在所有格点中寻找最大的psi值及对应的点即为磁轴点。2.4.2等离子体边界上x-point点的确定x-point点为等离子体边界磁通面上磁场强度BP＝0的点，具体寻找边界上xpoint点的步骤如下。（1）计算等离子体边界每点上的磁场强度drdrBPz1，22)()(1dzddrdrBP（2）找到所有边界点上磁场强度最小的点。（3）确定x-point点，即找到磁场强度BP＝0的点。5222222)(drdzddzddrd，drddzddzddrdzdx222，dzddrddrddrdzdy222修正点为：xs=xs+x,ys=ys+y，当100.122eyx，则找到的xpoint点收敛于边界，否则，认为找到的xpoint点不收敛。2.4.3计算归一化磁通函数~从以上计算可知，网格区域内所有点的极向磁通值、等离子体体边界上的磁通值b和磁轴上的磁通值m。所以相应地我们可以得到归一化磁通函数bmm~。2.5判断两次计算网域中磁通值的收敛性主要计算R-Z网格上磁通的变化值，如果每次迭代前后的最大相对误差值小于给定的误差值，则迭代过程结束。这时已经找到Grad-shafranov方程的平衡解。2.6反演矢量参数的计算在反演计算模式时，在线性最小二乘法的基础上，通过诊断数据，如磁场和磁通的测量来拟合等离子体电流剖面多项式模型中的多项式系数，从而得到等离子体的平衡位形。（1）同样是解矩阵方程AX=B得到极向场电流值和相应的多项式系数。A矩阵是与极向场线圈以及与P和F相关的互感系数矩阵，B向量是对应单匝环、小探针以及等离子体总电流的诊断数据，X向量对应于我们需要拟合得到的极向场线圈电流和等离子体电流剖面多项式模型的系数。（2）分别计算单匝环、磁探针和总的等离子体电流的计算值。（3）计算liiiCM222)(，主要涉及到单匝环、小探针和等离子体电流的测量值和计算值。（4）判断是否满足最小二乘法的误差要求。