磁性流体 Grad-Shafranov 方程
关键词： 轴对称 稳态 线性 磁性流体
Grad-Shafranov 方程是轴对称、稳态条件下的磁性流体(MHD)方程，主要用来分析核聚变发电装置 Tokamak。
Tokamak 是一种利用磁约束来实现受控核聚变的环形容器，在 20 世纪 50 年代由位于苏联莫斯科的库尔恰托夫研究所的阿齐莫维齐等人发明。Tokamak 的中央是一个环形的真空室，外面缠绕着线圈。在通电的时候 Tokamak 的内部会产生巨大的螺旋型磁场，将其中的等离子体加热到很高的温度，以达到核聚变的目的。
本节以 Grad-Shafranov 方程为例，讲解如何使用 FEtch 系统对柱坐标系下的线性椭圆型方程进行求解。
控制方程
稳态状态下的 MHD 方程表述了等离子体压强 $P$ 和 Lorentz 力间的平衡，即 \(J \times B=\nabla P\)
其中，$J$ 和 $B$ 为电流密度和磁感应强度。将该条件代入 Maxwell 方程 $\nabla \times B=\mu_0 J $，并将磁场 $B$ 用矢量势 $A$ 的旋度来表示，即 $B= \nabla \times A$ ，可得稳态下的磁性流体方程
\[\nabla \times( \nabla \times A)=\mu_0 \nabla P\]由于等离子体必须满足在径向和轴向的力学平衡，因此，上式可进一步化为 Grad-Shafranov 方程 \(\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial \Phi}{\partial r}\right)+\frac{\partial^2\Phi}{\partial z^2} =-\mu_0 r^2 \frac{\partial P}{\partial \Phi} -g \frac{\partial g}{\partial \Phi} \tag {1}\label{eq1}\)
其中，$\Phi$ 是极向磁通量，$g$ 是环向场函数，$r$ 和 $z$ 分别为径向和轴向坐标。 $\eqref{eq1}$ 式为 Grad-Shafranov 方程的一种通用形式。该方程的右端
\[J_\phi =-\mu_0 r^2 \frac{\partial P}{\partial \Phi} -g \frac{\partial g}{\partial \Phi}\]用于表征被磁场约束在环形圈内的等离子体的环形方向的流动。不同的右端项会有不同的解。为了方便我们的计算结果与现有的精确解进行对比，这里令 $J_\phi =f_0(r^2+r_0^2)$。
弱形式
根据变分原理，Grad-Shafranov 方程的弱形式为
\[\int_\Omega \left(\frac{\partial\Phi}{\partial r}\frac{\partial\delta\Phi}{\partial r}+\frac{\partial\Phi}{\partial z}\frac{\partial\delta\Phi}{\partial z}\right) r\mathrm{d}r\mathrm{d}z = \int_\Omega J_\phi \delta \Phi\ r\mathrm{d}r\mathrm{d}z + \int_{\partial \Omega} \frac{\partial \Phi}{\partial n} \delta \Phi\ r\mathrm{d}S \tag {2}\label{eq2}\]由于一般不允许磁束向外流出，方程的最后一项通常为零。上式即为有限元计算的控制方程。
如果与直角坐标系下的弱解形式进行对比，可以发现唯一的不同点在于柱坐标系下的积分项与径向坐标 $r$ 相关，因此，程序编制的关键是如何将径向坐标作为参数引入到单元矩阵的计算中。
对式 $\eqref{eq2}$ 进行空间离散化后，可进一步写成矩阵形式 \(SU=F\) 其中，$S$ 为刚度矩阵，$U$ 为磁通量向量，$F$ 为载荷向量。该式即为最终需要求解的线性方程组。
脚本编写
一个具体问题的有限元语言表述，通常需要包含 5 种类型的文件：PDE、FBC、MDI、GCN 和 NFE。由于方程 $\eqref{eq2}$ 不存在边界积分项，所以这里不需要写 FBC 文件。下面依次给出其他 4 个文件的详细介绍。
PDE 文件
在这个例子中，只有一个物理场，需要给出极向磁通量对应的 PDE 文件。根据有限元语言的语法规则，将方程 $\eqref{eq2}$ 的体积分项写入 PDE 文件如下。
| grad.pde | 说明 |
|---|---|
| defi disp u coor r,z shap %1 %2 gaus %3 mate f0 r0 1d0 1d0 stif $c6 fj=f0*(r*r+r0*r0) dist =+[u/r;u/r]*r+[u/z;u/z]*r load =+[u]*fj*r end |
系统关键字（信息定义） 定义未知量 定义坐标分量 定义采用的形函数 定义采用的积分方式 定义用到的材料参数和默认参数值 （空行） 系统关键字（单元刚度矩阵表达式定义） 刚度矩阵项，对应式 $\eqref{eq2}$ 左侧 （空行） 荷载向量项，对应式 $\eqref{eq2}$ 右侧第 1 项 （空行） 系统关键字（文件结束） |
这里将 shap 行和 gaus 行的参数设为待定，参数值 %1、%2 和 %3 将从 MDI 文件中取得。
MDI 文件
MDI 文件负责组织 PDE 和 FBC 文件，声明程序具体使用的坐标系统、单元类型和积分方法。这里文件命名为 gs.mdi，具体内容如下。
| gs.mdi | 说明 |
|---|---|
| 2drz #a 0 1 u pde grad t3g2 # |
在二维轴对称柱坐标系下求解 a场，0 个初值，1 个自由度，自由度名称为 u 调用 grad.pde，采用 3 节点三角形单元，2 点高斯积分 （文件结束符，可留空） |
MDI 文件声明的单元类型和积分方式将以参数的形式传递给 PDE 文件中的的 %1、%2 和 %3。如果后续想要使用其他的单元类型或者积分方法，只需修改 MDI 文件即可。
GCN 文件
GCN 文件给出计算流程信息，系统规定其文件名字要与 MDI 文件一致。作为线性椭圆方程，Grad-Shafranov 方程的计算流程非常简单，只需填写少量代码。文件的具体内容如下。
| gs.gcn | 说明 |
|---|---|
| defi a ell #startsin a #solvsin a |
算法引入段 a 场，采用的算法模板文件为 ell.nfe (空行) 开始命令流段 初始化 a 场 求解 a 场，使用直接法，对称矩阵求解器 |
NFE 文件
GCN 文件使用了 ell 算法，系统会基于模板文件自动生成以下的 gsa.nfe 文件。如果用户有其他需要，可以在此基础上进一步地修改。
| gsa.nfe | 说明 |
|---|---|
| defi stif s load f equation matrix = [s] forc = [f] solution u gidpost('u',1.0) u(1) end |
DEFI段 定义单元刚度矩阵s 定义荷载向量f （空行） EQUATION段 总刚矩阵公式 荷载向量公式 （空行） SOLUTION段 定义求解的场向量 将计算结果写入 GiD 后处理文件 （空行） 结束标记 |
算例
问题描述
计算区域为 \(\partial \Omega=\left\{(r, z) \mid r=r_{0} \sqrt{1+\frac{2 a \cos \alpha}{r_{0}}},z=a r_{0} \sin \alpha, \alpha=[0,2 \pi]\right\}\) 在边界上 $\Phi=0$。此时的理论解为 \(\psi\left(r,z\right) = \frac{f_{0}r_{0}^{2}a^2}{2}\left(1 - \frac{z^{2}}{a^{2}} - \left(\frac{r-r_{0}}{a} + \frac{\left(r-r_{0}\right)^{2}}{2ar_{0}}\right)^{2} \right)\)
这里取 $a=0.5$，$f_0=1.0$，$r_0=1.0$。试求该稳态的有限元解。
前处理
填写材料参数文件
系统自动生成了 gs.mat 文件，用于设置材料参数。参数初始值全部为 PDE 文件给出的默认值。可根据具体问题进行修改。对于本算例，修改完成后的文件内容如下。
| gs.mat | 说明 |
|---|---|
| 1 3 aet3g2 num f0 r0 1 1d0 1d0 # |
3 节点三角形单元 材料参数 参数值 分隔符 |
赋单元材料号
施加边界条件
其中：
u-I 为边界类型，-1 表示约束，1 表示自由
u-D 为相应的约束值或施加节点流量
网格剖分
根据所开发程序的计算要求，采用三角形线性单元进行网格剖分，共 2812 个节点，5422 个单元。
计算结果
计算结果
解析解
上图依次给出了计算结果和解析解的云图。可以看出，两者是一致的，充分证明了算法和程序的有效性。
参考文献
[1] Deriaz E, Despres B, Faccanoni G, et al. Magnetic equations with FreeFem++: the Grad-Shafranov equation & the current hole[C]//Esaim: Proceedings. EDP Sciences, 2011, 32: 76-94.
感谢您的支持，我们会继续努力的!
打开支付宝或微信扫一扫，即可进行扫码打赏哦
