Grad-Shafranov方程的有限元求解

1．Grad-Shafranov方程

Grad-Shafranov方程是轴对称，稳态条件下的磁性流体(MHD)方程。稳态状态下的MHD方程表述流体压力P和Lorentz力间的平衡。

J×B=∇P

在这里J和B为电流和磁场强度。

将该条件代入Maxwell方程 μ0J=∇×B ，并将磁场B用矢量势A的形式 B=∇×A 来表示可得稳态下的磁性流体方程

(1)∇×(∇×A)=μ0∇P

Grad-Shafranov方程主要用来核聚变发电装置Tokamak(图1)。（当然也可以应用于其他方面，如参考文献[2]）. Tokamak是一个对称的环形装置，因此我们需要把它改写为轴对称的形式。

Tokamak装置·[1]

(2)(∂2Φ∂r2 ∂2Φ∂z2)−1r∂Φ∂r=−μ0r2∂P∂Φ−g∂g∂Φ

该方程即为Grad-Shafranov方程。在这里磁束函数Φ(r,z)=rAϕ(r,z) 为矢量势A成分 Aϕ(r,z) 的函数, 环向场函数 g(r,z)=r(∂Ar∂z−∂Az∂r) 。该方程的右端 Jϕ=−μ0r2∂P∂Φ−g∂g∂Φ 表征被磁场约束在环形圈内的等离子体的环形方向的流动。

2. 方程的弱形式

在方程两侧乘上试函数V∈H1(Ω)并积分，得

(3)∫Ω∇Φ⋅∇VdΩ ∫∂Ω1r∂Φ∂rVdΩ=∫ΩJϕVdΩ ∫∂Ω∂Φ∂nVdS

由于我们一般不允许磁束向外流出。该方程的最后一项为零。该方程为有限元计算控制方程。

3. 验证计算

使用参考文献[4]的first example，这是一个有理论解(Soloviev solution)的问题，其理论解为

ψ(r,z)=f0r022(1−z2a2−(r−r0a (r−r0)22ar0)2)

右边源项 Jϕ=f0(r2 r02) . 当a=0.5, f0=1.0,r0=1.0 时开源软件

hillyuan/四元gitee.com/hillyuan_be3b/siyuan

的计算结果为

开源软件四元的计算结果

下面是参考文献[4]和[5](计算条件稍有不同)的计算结果，可见结果基本相符。

文献·[4]结果

文献[5]结果

参考文献

1）S. Li, H. Jiang, Z. Ren, C. Xu - S. Li et. al. "Optimal Tracking for a Divergent-Type Parabolic PDE System in Current Profile Control" , Abstract and Applied Analysis、Vol.2014

2）藤澤幸太郎: Grad-Shafranov方程式を用いた磁場星の計算

3）Dama Adams Daniel等: Derivation and Applications of Grad-Shafranov Equation In Magnetohydrodynamics (MHD)，Journal of Research in Applied Mathematics，Volume 7， Issue 4 (2021) pp: 34-38

4）Erwan Deriaz, Bruno Despres, Gloria Faccanoni, Kirill Pichon Gostaf, Lise-Marie Imbert-Gérard, Georges Sadaka, and Remy Sart. Magnetic equations with FreeFem : The Grad-Shafranov equation & the current hole. ESAIM: Proceedings, 32:76–94, November 2011

5）Linear VC Grad-Shafranov solver for Soloviev solution I (FEniCS)