双撕裂模不稳定性的伽辽金谱方法模式分析

Mode analysis on double tearing mode instability using Galerkin spectral method

摘要

本研究旨在探究双撕裂模（Double Tearing Mode， DTM）不稳定性非线性演化过程复杂三波相互作用的特征。采用伽辽金谱方法数值求解平板位形下的约化磁流体动力学两场方程，该方法利用全局正交基函数（x向正弦级数，y向傅里叶级数）进行空间离散，其模式滤波优势适用于非线性模耦合分析。首次揭示了准线性阶段模结构的对称性：基模本征结构为偶对称；其经奇数次自相互作用产生的谐波模表现为奇对称，而经偶数次作用产生的谐波模则表现为偶对称。此对称性仅由谐波生成路径所涉基模自相互作用次数的奇偶决定，与电流片间距无关。通过包含多初始扰动的模拟，得到DTM谐波模的空间分布是各基模独立自相互作用结果的线性叠加，且各独立结果的对称性符合各自作用路径的奇偶特征。本文利用伽辽金谱方法有效地解析复杂非线性模耦合过程，阐明了DTM三波相互作用中模结构的对称性与可叠加性规律。

Abstract

This study investigates three-wave interactions in the nonlinear evolution of the double tearing mode （DTM） instability. We solve the reduced magnetohydrodynamics two-field equations in slab geometry using the Galerkin spectral method. This method employs global orthogonal bases （x-sine， y-Fourier series） for spatial discretization， leveraging its mode-filtering strength for nonlinear coupling analysis. A novel symmetry rule governing the mode structure in the quasi-linear phase is revealed： the fundamental mode eigenstructure exhibits even symmetry； harmonics generated via an odd number of self-interactions of the fundamental mode possess odd symmetry， while those generated via an even number of self-interactions exhibit even symmetry. This symmetry depends only on the interaction path parity （odd/even count）， not on current sheet separation. Simulations with multiple initial perturbations show the DTM harmonic distribution is a linear superposition of outcomes from each fundamental mode's independent self-interactions. Each outcome's symmetry strictly follows its path's parity characteristic. This work effectively utilizes the Galerkin spectral method to decipher complex nonlinear mode coupling， elucidating the symmetry and linear superposition governing DTM three-wave mode structures.

Graphical abstract

关键词

双撕裂模 / 磁重联 / 磁流体动力学 / 托卡马克

Key words

double tearing mode / magnetic reconnection / magnetohydrodynamics / Tokamak

磁重联作为磁场拓扑结构的动态重构过程，自Giovanelli首次提出解释太阳耀斑机制^［1-2］以来，已成为天体物理^［3-4］、空间环境^［5-7］及实验室等离子体^［8-10］研究的重要内容。该过程通过磁力线断开与重联，将磁能迅速转化为等离子体的动能与热能，在太阳耀斑能量释放^［11-12］、日冕加热^［4］、日冕物质抛射（Coronal Mass Ejection， CME）^［13］、磁层亚暴^［14］以及托卡马克破裂事件^［15］中起着关键作用。非稳态磁场重联主要由各种撕裂模不稳定性所触发，引起电流片内贮存磁能的释放，并在有限电阻效应下形成磁岛结构。

当两个电流片相邻时，撕裂模之间相互驱动使得重联过程更加剧烈，发生双撕裂模（Double Tearing Mode， DTM）不稳定性的演化。该过程在具有负磁剪切特性的先进托卡马克装置中尤为显著，其引发的离轴锯齿振荡直接制约等离子体约束性能。Pritchett等^［16］基于非常数磁通假设解析电阻磁流体动力学（Magnetohydrodynamics， MHD）方程，首次建立DTM线性增长率与电阻率的定标律，发现DTM解耦条件。Ishii等^［17］通过数值模拟给出非线性过程各个扰动的演化，捕捉到DTM扰动能量爆发性增长特征；Wang等^［18］采用平板模型明晰了DTM非线性演化的多阶段定标关系。DTM的非线性演化，其本质为各个扰动相互作用并交换能量的过程。扰动模通过级联（cascade）将能量传递给高次谐波，形成大量三波相互作用。DTM的相互驱动，以及扰动模之间的相互作用，使得准非线性理论无法完整描述其非线性演化过程。

针对上述非线性研究的挑战，谱方法展现出独特优势。谱方法起源于瑞兹-伽辽金法（Ritz–Galerkin Methods），主要包括伽辽金法、tau方法和配置法^［19］。与有限差分（Finite Difference）法或有限元（Finite Element）法的局部离散特性不同，谱方法通常采用全局无穷可微的正交多项式作为试函数，使任意节点的计算均蕴含全域信息^［19］。这一特性使得当解自身无穷光滑时，谱方法收敛速度超越1/N的任意有限次幂（N为展开项数）。近30年来，谱方法在计算流体力学领域的成功应用，验证了其在可压缩流动、不可压缩流动、边界层流动、平行流及湍流流动等多种问题中的卓越性能^［19-21］。近年来，在MHD研究中，全局代码（如NIMROD^［22］、XTOR^［23］与CLT^［24］）通过在环向或者极向等周期方向引入谱分解，显著提升导数计算精度以及简化计算量。这些代码的径向离散格式通常采用有限差分法或有限元法。

最近，在以伽辽金方法作为唯一空间离散形式来研究撕裂模线性特性的工作中，此方法展现出便于模式分析以及滤波的特点^［25］，这一优势为非线性模耦合的研究提供了有力支撑。本文通过伽辽金方法求解约化磁流体动力学（Reduced Magnetohydrodynamics， RMHD）两场方程，对DTM的非线性演化过程开展了数值研究，并通过模式分析探究了此过程中三波相互作用的特征。

1 物理模型与数值方法

1.1 RMHD两场模型

本文重点讨论非线性过程的模式耦合特性。为排除环几何位形下磁面位移引起的线性的环耦合等干扰因素，特采用平板位形下的RMHD两场模型开展研究。假设平衡磁场在x方向具有剪切特性，其可表示为磁通函数形式\mathit{B}=B_{\mathrm{T}}\widehat{z}+\nabla \psi \times \widehat{z}。其中，导向磁场设为B_{\mathrm{T}}=\mathrm{0}， \widehat{z}为垂直平面方向单位矢量，\psi \left(x,y\right)为磁通函数。在不可压缩近似条件下，等离子体密度\rho均匀且无扰动，且不考虑平衡流速（\mathit{v}=\mathrm{0}）。流速可以表示为流函数的形式\mathit{v}=\nabla \varphi \times \widehat{z} （其中，\varphi \left(x,y\right)为流函数），对流速取旋度可得涡量与流函数的关系：\mathit{U}=\nabla \times \mathit{v}=-{\nabla }_{\perp }^{\mathrm{2}}\varphi \widehat{z}=U\widehat{z}。基于感应方程与运动方程，可推导得到关于磁通函数和涡量的归一化两场方程^［26］，如式（1）和式（2）所示。

\frac{\partial \psi }{\partial t}=\left[\psi ,\varphi \right]+S^{-\mathrm{1}}{\nabla }_{\perp }^{\mathrm{2}}\psi(1)

\frac{\partial U}{\partial t}=\left[U,\varphi \right]+\left[j_z,\psi \right](2)

其中，泊松括号定义为\left[f,g\right]=\nabla f\times \nabla g=\frac{\partial f}{\partial x}\frac{\partial g}{\partial y}-\frac{\partial f}{\partial y}\frac{\partial g}{\partial x}，拉普拉斯算子为{\nabla }_{\perp }^{\mathrm{2}}=-(\frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}}{\partial y^{\mathrm{2}}})，电流密度为j_z=-{\nabla }_{\perp }^{\mathrm{2}}\psi。

选取如下归一化参数：长度L、速度v、时间t和磁场B，分别以特征长度L、特征阿尔芬速度v_{\mathrm{A}}=B_L/{\left({\mu }_{\mathrm{0}}\rho \right)}^{\mathrm{1}/\mathrm{2}}、阿尔芬时间{\tau }_{\mathrm{A}}=L/v_{\mathrm{A}}及x=L处的特征磁场强度B_L进行归一化。S={\tau }_{\eta }/{\tau }_{\mathrm{A}}为磁雷诺数，其中，{\tau }_{\eta }={\mu }_{\mathrm{0}}{L}^{\mathrm{2}}/\eta为电阻扩散时间，{\mu }_{\mathrm{0}}和\eta则分别为真空磁导率和等离子体电阻率。特征电场E_L=v_{\mathrm{A}}{B}_L，特征电流密度J_{\mathrm{L}}=B_L/{\mu }_{\mathrm{0}}L，特征磁通量{\psi }_L=B_{L}L，特征流函数{\varphi }_L=v_{\mathrm{A}}L，特征涡量U_L=v_{\mathrm{A}}/L。

式（1）和式（2）中的各场变量f\left(x,y,t\right)均可表示为初始平衡量f_{\mathrm{00}}与随时间演化的扰动量f_{\mathrm{11}}之和，即f\left(x,y,t\right)=f_{\mathrm{00}}+f_{\mathrm{11}}(x,y,t)。那么式（1）和式（2）可写为式（3）和式（4）形式。

\frac{\partial {\psi }_{\mathrm{11}}}{\partial t}=\left[{\psi }_{\mathrm{00}}+{\psi }_{\mathrm{11}},{\varphi }_{\mathrm{11}}\right]+S^{-\mathrm{1}}{\nabla }_{\perp }^{\mathrm{2}}{\psi }_{\mathrm{11}}(3)

\begin{array}{l}\frac{\partial U_{\mathrm{11}}}{\partial t}=\left[U_{\mathrm{11}},{\varphi }_{\mathrm{11}}\right]+\left[-{\nabla }_{\perp }^{\mathrm{2}}{\psi }_{\mathrm{00}},{\psi }_{\mathrm{11}}\right]+\\ \left[-{\nabla }_{\perp }^{\mathrm{2}}{\psi }_{\mathrm{11}},{\psi }_{\mathrm{00}}+{\psi }_{\mathrm{11}}\right]\end{array}(4)

1.2 伽辽金谱方法

早前的工作采用伽辽金方法对式（3）和式（4）线性化后的方程组进行了数值求解^［25］。该方法在全空间采用谱分解，即通过选用满足相应边界条件的傅里叶基函数，将控制方程投影至谱空间进行求解。研究表明，该方法对单、双撕裂模以及m=1撕裂模线性结果的分析与撕裂模线性理论一致，包括模结构、增长率与电阻的定标等。特别地，验证了DTM解耦时电阻定标率从\gamma ~{\eta }^{\mathrm{1}/\mathrm{3}}到\gamma ~{\eta }^{\mathrm{3}/\mathrm{5}}的变化。同时，通过与伪谱法在收敛性测试等方面的对比，证实了伽辽金方法在模拟广域模（如DTM）时的优越性。

为研究共振面（或电流片）足够接近且距离模拟边界足够远条件下的DTM，可将式（3）与式（4）投影至满足相应边界条件的试探函数空间。x方向两侧（x=\mathrm{0},a）设置为理想导体壁，因此试探函数采用正弦级数\mathrm{s}\mathrm{i}\mathrm{n}\left(k_{x}x\right)；y方向两侧（y=\mathrm{0},b）为周期边界，故试探函数采取傅里叶级数{\mathrm{e}}^{\mathrm{i}{k}_{y}y}。平衡量f_{\mathrm{00}}和扰动量f_{\mathrm{11}}可分别表示为式（5）与式（6）形式。

f_{\mathrm{00}}\left(x,t\right)\approx {\sum }_{l=-N_x}^{N_x}{f}_l\mathrm{s}\mathrm{i}\mathrm{n}\left(k_{x}x\right)(5)

f_{\mathrm{11}}\left(x,y,t\right)\approx {\sum }_{m=-N_y}^{N_y}{\sum }_{l=-N_x}^{N_x}{f}_{l,m}\left(t\right)\mathrm{s}\mathrm{i}\mathrm{n}\left(k_{x}x\right){\mathrm{e}}^{\mathrm{i}{k}_{y}y}(6)

其中，-N_x\le l\le N_x和-N_y\le m\le N_y为谱空间的模拟区域，f_l为平衡量的展开系数，f_{l,m}(t)为扰动量的展开系数（对应傅里叶分量(k_x,k_y)的振幅）。角标l和m分别代表x与y方向的模数，它们与波数分别满足关系k_x=\mathrm{\pi }l/a以及k_y=\mathrm{2}\mathrm{\pi }m/b。将式（5）和式（6）代入式（3）与式（4），并利用试探函数与检验函数的正交性，便可将原表达式转化为展开系数的演化方程，如式（7）和式（8）所示。

\frac{\partial {\psi }_{l,m}}{\partial t}=-S^{-\mathrm{1}}{k}^{\mathrm{2}}{\psi }_{l,m}-\sum _{l\text{'},l^{\text{'}\text{'}}}{P}_{\mathrm{1}}{l}^{\text{'}}m\mathrm{i}{\psi }_{l^{\text{'}}}{\varphi }_{l^{\text{'}\text{'}},m}+\sum _{m^{\text{'}},m^{\text{'}\text{'}}}\sum _{l^{\text{'}},l^{\text{'}\text{'}}}\left(P_{\mathrm{2}}{l}^{\text{'}\text{'}}{m}^{\text{'}}-P_{\mathrm{1}}{l}^{\text{'}}{m}^{\text{'}\text{'}}\right)\mathrm{i}{\psi }_{l^{\text{'}},m^{\text{'}}}{\varphi }_{l^{\text{'}\text{'}},m^{\text{'}\text{'}}}(7)

\frac{\partial U_{l,m}}{\partial t}=\sum _{l^{\text{'}},l^{\text{'}\text{'}}}-P_{\mathrm{1}}\left[l^{\text{'}}m\left({l^{\text{'}\text{'}}}^{\mathrm{2}}+m^{\mathrm{2}}\right)-{l^{\text{'}\text{'}}}^{\mathrm{3}}m\right]\mathrm{i}{\psi }_{l^{\text{'}}}{\psi }_{l^{\text{'}\text{'}},m}+\sum _{m^{\text{'}},m^{\text{'}\text{'}}}\sum _{l^{\text{'}},l^{\text{'}\text{'}}}\left(P_{\mathrm{2}}{l}^{\text{'}\text{'}}{m}^{\text{'}}-P_{\mathrm{1}}{l}^{\text{'}}{m}^{\text{'}\text{'}}\right)\mathrm{i}{U}_{l^{\text{'}},m^{\text{'}}}{\varphi }_{l^{\text{'}\text{'}},m^{\text{'}\text{'}}}+\left(P_{\mathrm{2}}{l}^{\text{'}\text{'}}{m}^{\text{'}}-P_{\mathrm{1}}{l}^{\text{'}}{m}^{\text{'}\text{'}}\right)\left({l^{\text{'}\text{'}}}^{\mathrm{2}}+{m^{\text{'}\text{'}}}^{\mathrm{2}}\right){\psi }_{l^{\text{'}},m^{\text{'}}}{\psi }_{l^{\text{'}\text{'}},m^{\text{'}\text{'}}}(8)

其中，y方向模数满足m=m^{\text{'}}+m^{\text{'}\text{'}}；P_{\mathrm{1}}与P_{\mathrm{2}}分别为x方向模数l、l^{\text{'}}、l^{\text{'}\text{'}}的函数，其表达式如式（9）和式（10）所示。

P_{\mathrm{1}}=\left\{\begin{array}{c}-\frac{\mathrm{1}}{\mathrm{2}},l^{\text{'}\text{'}}=l+l^{\text{'}}\\ \frac{\mathrm{1}}{\mathrm{2}},l^{\text{'}\text{'}}=l-l^{\text{'}}\mathrm{o}\mathrm{r}{l}^{\text{'}\text{'}}=l^{\text{'}}-l\end{array}\right.(9)

P_{\mathrm{2}}=\left\{\begin{array}{c}\frac{\mathrm{1}}{\mathrm{2}},l^{\text{'}\text{'}}=l+l^{\text{'}}\mathrm{o}\mathrm{r}{l}^{\text{'}\text{'}}=l-l^{\text{'}}\\ -\frac{\mathrm{1}}{\mathrm{2}},l^{\text{'}\text{'}}=l^{\text{'}}-l\end{array}\right.(10)

式（7）和式（8）中流函数与涡量的展开系数满足式（11）。

{\varphi }_{l,m}=\frac{U_{l,m}}{\left(l^{\mathrm{2}}+m^{\mathrm{2}}\right)}(11)

2 数值结果与分析

2.1 平衡与初始扰动

本文采用平板位形进行数值研究，模拟区域设置为\mathrm{0}\le x\le \mathrm{2},\mathrm{0}\le y\le \mathrm{2}。对于DTM，磁场平衡分布表达式如式（12）所示。

B_{\mathrm{00}}\left(x\right)=\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\frac{x-x_l}{d}-\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\frac{x-x_r}{d}+\mathrm{1}(12)

图1展示了电流片间距分别为x_r-x_l =0.2和x_r-x_l =1.0时各平衡量的分布。其中，2个哈里斯电流片（Harris Current Sheet）^［27］分别位于x_l与x_r处，其半宽为d=\mathrm{0.1}。

模拟区域的边界条件设定如下：x=0和x=2处为导体壁边界，y=0和y=2处为周期性边界，则谱方法的初始扰动设置如式（13）所示。

{\psi }_{|m|}\left(x,y\right)=\widetilde{{\psi }_{|m|}}\left\{\mathrm{e}\mathrm{x}\mathrm{p}\left[-{\left(\frac{x-x_l}{d}\right)}^{\mathrm{2}}\right]+\mathrm{e}\mathrm{x}\mathrm{p}\left[-{\left(\frac{x-x_r}{d}\right)}^{\mathrm{2}}\right]\right\}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\pi }my\right)(13)

其中，角标|m|代表正负模（撕裂模为驻波解），且在后文代表初始扰动（基模）模数，\widetilde{{\psi }_m}代表初始扰动的幅值。

其他相关参数如下：谱空间k_x 方向的离散网格数为N_x=\mathrm{512}，此参数在早前工作中已验证收敛性要求^［25］；伦德奎斯特数（Lundquist Number）设定为S={\mathrm{10}}^{\mathrm{6}}，以兼顾数值稳定性要求；时间积分选取四阶龙格-库塔（Runge-Kutta）格式。

2.2 磁通函数及磁岛演化

图2分别给出了电流片间距为x_r-x_l =0.2和x_r-x_l =1.0时撕裂模（m=\mathrm{0}\mathrm{、}\mathrm{1}\mathrm{、}\dots \mathrm{、}\mathrm{8}）的磁通函数与总磁能E_M随时间的演化。由图2b可见，当电流片间距为x_r-x_l =1.0时，两电流片相距过远，导致DTM解耦；且总磁能的转化显著发生在非线性阶段。

图3a~d对应图2a小共振面间距算例中t=22、24、28、40时刻磁通函数的等高图。2个单撕裂模的磁岛经过准线性阶段后，在t=\mathrm{22}时刻发展至足够宽度；由于电流片相距较近，两磁岛开始相互驱动。t=24时刻可见，两磁岛O点外侧分别形成X点，这进一步提升了总磁能E_M转化的速度。受磁张力影响，两磁岛发生相向位移（如t=\mathrm{28}时刻所示）。随后两磁岛大部分结构重联完毕，t=\mathrm{40}时刻可见残余部分。此演化过程与Wang等^［18］的研究成果一致。图3e~h对应图2b大共振面间距算例中t=22、26、30、50时刻磁通函数的等高图。由于αx_s =0.314（α²/S）^1/9=0.166（其中α=k_yd=0.1π），DTM发生解耦，表现为2个单撕裂模的发展，这与Pritchett等^［16］的理论描述相符。此时，2个单撕裂模的磁岛在非线性阶段增长相对缓慢，且m=\mathrm{0}模的增长改变了平衡状态，最终使其达到饱和，此过程中E_M的转化速度也相对较慢。

由于本文研究的三波相互作用涉及大范围模数的耦合，因此需对模数N_y进行测试。图4对比了线性模拟与不同截止模数N_y下的非线性模拟结果。图4a为DTM的线性模拟结果，图4b~d分别对应N_y=\mathrm{8}\mathrm{、}\mathrm{12}\mathrm{、}\mathrm{16}时DTM的非线性模拟结果。各非线性模拟中，m=1模的增长率与线性模拟近似一致，验证了结果的准确性；且不同N_y下m=\mathrm{1}\mathrm{、}\mathrm{2}\mathrm{、}\mathrm{3}\mathrm{、}\mathrm{4}模在准线性阶段的增长率也保持一致——即不同截止模数计算得到的低模数增长率收敛于同一数值。这表明，对于撕裂模这种长波不稳定的宏观不稳定性模式，m\mathrm{8}的更高模式对于低模数的级联可忽略不计。

2.3 模结构对称性

由图2可知，不论DTM解耦与否，磁通函数在准线性阶段均以固定增长率呈指数形式增长。Wang等^［28］通过准线性理论，解释了三波相互作用过程中各模数增长率之间的关系。本文从磁通函数的空间分布入手，探究其在准线性阶段进行三波相互作用的特征。为与撕裂模本征模结构相区分，称谐波准线性阶段的空间分布为非线性模结构。图5为不同电流片间距以及不同基模条件下，非线性模结构在x方向的分布，可见其空间分布的对称性随三波相互作用的复杂度而变化。

图5a对应图2a的算例，其中实线与虚线分别代表奇和偶对称的模结构。基模为\left|m\right|=\mathrm{1}，其本征模结构为偶对称；而经1次自相互作用得到的m=\mathrm{0}\mathrm{、}\mathrm{2}非线性模结构为奇对称；m=\mathrm{3}\mathrm{、}\mathrm{5}\mathrm{、}\mathrm{7}与m=\mathrm{4}\mathrm{、}\mathrm{6}\mathrm{、}\mathrm{8}非线性模结构则分别为偶对称与奇对称。图5b包含2个算例：当基模为\left|m\right|=\mathrm{2}时，m=\mathrm{0}\mathrm{、}\mathrm{4}与m=\mathrm{6}模分别呈奇对称与偶对称（此时非线性作用不产生奇数模）；当基模为\left|m\right|=\mathrm{3}时，m=\mathrm{0}\mathrm{、}\mathrm{6}与m=\mathrm{9}模分别呈奇对称与偶对称。同一模数的对称性时而为奇、时而为偶。这表明模结构的对称性与基模自相互作用次数的奇偶有关。如，m=\mathrm{6}模的作用过程为“6=1+1+1+1+1+1=2+2+2=3+3”。综上可知，基模本征结构为偶对称，基模经奇（或偶）数次自相互作用得到的非线性模结构为奇（或偶）对称。图5c、d为电流片间距x_r-x_l =1.0的对比算例。结果表明，不论DTM解耦与否，磁通函数对称性的规律一致。

2.4 模结构可叠加性

在实际等离子体环境中，扰动总是随机的，可通过傅里叶变换分解为若干模——即初始扰动并非单一模。因此，本节从包含2个初始扰动的模拟入手，探究三波相互作用较为复杂的情况下，上述规律是否仍然保持。

基模分别取\left|m\right|=\mathrm{2}\mathrm{、}\mathrm{3}、 \left|m\right|=\mathrm{2}和\left|m\right|=\mathrm{3}时，DTM磁通函数的时间演化与空间分布如图6所示。图6a为演化过程，不同颜色代表不同模数；其中m=\mathrm{1}\mathrm{、}\mathrm{5}\mathrm{、}\mathrm{7}模仅在\left|m\right|=\mathrm{2}\mathrm{、}\mathrm{3}算例中出现，且它们在准线性阶段的增长率恒定。而在\left|m\right|=\mathrm{2}\mathrm{、}\mathrm{3}算例中，m=\mathrm{0}\mathrm{、}\mathrm{6}模的增长率在准线性阶段发生变化：前期与\left|m\right|=\mathrm{3}算例的增长率一致，后期则与\left|m\right|=\mathrm{2}算例的增长率重合。如m=\mathrm{0}模划线先沿着点线发展，后沿着实线发展。图6b为t=\mathrm{5}时刻m=\mathrm{0}模的空间分布，可见\left|m\right|=\mathrm{2}（绿实线）与\left|m\right|=\mathrm{3}（蓝实线）算例中m=\mathrm{0}模结构的叠加结果（黑点线），与\left|m\right|=\mathrm{2}\mathrm{、}\mathrm{3}算例中的（黑划线）重合。图6c为t=\mathrm{10}时刻m=\mathrm{6}模的空间分布，其特征与m=\mathrm{0}模一致。上述时间演化与空间分布规律表明，多初始扰动时，各基模经自相互作用生成的同一模数，其模结构是线性叠加的。图6d为t=\mathrm{5}时刻m=\mathrm{1}\mathrm{、}\mathrm{5}\mathrm{、}\mathrm{7}模结构在任意单位下的对称性比较。它们分别由“1=3－2”、“5=2+3”和“7=2+2+3”作用过程生成，因此分别表现为奇对称、奇对称和偶对称。

图7取基模\left|m\right|=\mathrm{1}\mathrm{、}\mathrm{2}，此时除m=\mathrm{1}模外，其余模数均存在多条非线性作用途径。图7a给出了m=\mathrm{1}及部分低偶数模磁通函数的演化过程：其中m=\mathrm{2}\mathrm{、}\mathrm{4}\mathrm{、}\mathrm{6}\mathrm{、}\mathrm{8}模增长率发生1次或多次变化，其准线性阶段演化曲线（划线）的首、尾分别与\left|m\right|=\mathrm{2}（实线）、\left|m\right|=\mathrm{1}（点线）算例一致；而m=\mathrm{0}模的增长率恒定，这是因为\left|m\right|=\mathrm{2}算例的贡献（黑实线）远小于\left|m\right|=\mathrm{1}算例（黑点线），其影响可忽略不计。图7b为t=\mathrm{21}时刻m=\mathrm{2}模的空间分布，此时有“2=1+1=2” 2条路径，满足模结构可叠加性。t=\mathrm{15}时刻m=\mathrm{4}模的空间分布如图7c所示，其有“4=1+1+1+1=1+1+2=2+2” 3条路径，因此\left|m\right|=\mathrm{1}与\left|m\right|=\mathrm{2}两算例m=\mathrm{4}模结构（红、蓝实线）的叠加（橙划线），与\left|m\right|=\mathrm{1}\mathrm{、}\mathrm{2}算例（橙划线）不重合。根据对称性规律，路径“4=1+1+2”对应的非线性模结构为偶对称；由模结构可叠加性了，可得此模结构分布，即3个算例此模分布的相减（图7c中橙点线），符合规律。

3 结论

基于伽辽金谱方法数值求解RMHD两场方程，本文利用其便于模式分析优势，研究非线性DTM不稳定性扰动时间演化与空间分布，以及三波相互作用特征。结论如下： 1）揭示了模结构对称性规律——基模本征结构为偶对称，其经奇数次自相互作用生成的谐波模为奇对称，偶数次作用生成的谐波模为偶对称，此规律与电流片间距无关，仅由谐波模生成路径所涉及基模自相互作用次数的奇偶决定； 2）发现模结构可叠加性——通过包含多初始扰动的模拟，证明伽辽金方法能有效解析复杂非线性耦合，同一谐波模的空间分布是各基模独立自相互作用结果的线性叠加，且各独立结果的对称性符合其路径奇偶次数规律。本研究通过伽辽金谱方法，清晰揭示了DTM非线性演化中三波相互作用模结构的对称性与可叠加性特征。

参考文献

[1] Giovanelli R G.A theory of chromospheric flares ［J］.Nature， 1946， 158（4003）： 81-82.

[2] Giovanelli R G.Magnetic and electric phenomena in the Sun's atmosphere associated with sunspots ［J］. Mon Not R Astron Soc， 1947， 107（4）： 338-355.

[3] Wyper P F， Antiochos S K， DeVore C R.A universal model for solar eruptions ［J］.Nature， 2017， 544（7651）： 452-455.

[4] Tsuneta S.Structure and dynamics of magnetic reconnection in a solar flare ［J］.Astrophys J， 1996， 456： 840-849.

[5] Dungey J W.Interplanetary magnetic field and the auroral zones ［J］.Phys Rev Lett， 1961， 6（2）： 47-48.

[6] Boozer A H.Magnetic reconnection in space ［J］. Phys Plasmas， 2012， 19（9）： 092902.

[7] Angelopoulos V， Artemyev A， Phan T D， et al. Near-Earth magnetotail reconnection powers space storms ［J］.Nat Phys， 2020， 16（3）： 317-321.

[8] Furth H P， Rutherford P H， Selberg H.Tearing mode in the cylindrical tokamak ［J］.Phys Fluids， 1973， 16（7）： 1054-1063.

[9] Rutherford P H.Nonlinear growth of the tearing mode ［J］.Phys Fluids， 1973， 16（11）： 1903-1908.

[10] Wang J， Xiao C， Wang X， et al.Sawtooth-driven neoclassical tearing modes in HL-2A plasmas ［J］. Plasma Phys Control Fusion， 2012， 54（12）： 122001.

[11] Sturrock P A.Model of the high-energy phase of solar flares ［J］.Nature， 1966， 211（5050）： 695-697.

[12] Su Y， Veronig A M， Holman G D， et al.Imaging coronal magnetic-field reconnection in a solar flare ［J］.Nat Phys， 2013， 9（8）： 489-493.

[13] Dahlin J T， Antiochos S K， Devore C R. A model for energy buildup and eruption onset in coronal mass ejections ［J］.Astrophys J， 2019， 879（2）： 96.

[14] Angelopoulos V， McFadden J P， Larson D， et al.Tail reconnection triggering substorm onset ［J］.Science， 2008， 321（5891）： 931-935.

[15] Priest E R.The magnetohydrodynamics of current sheets ［J］.Rep Prog Phys， 1985， 48（7）： 955.

[16] Pritchett P L， Lee Y C， Drake J F. Linear analysis of the double‐tearing mode ［J］.Phys Fluids， 1980， 23（7）： 1368-1374.

[17] Ishii Y， Azumi M， Kurita G， et al.Nonlinear evolution of double tearing modes ［J］.Phys Plasmas， 2000， 7（11）： 4477-4491.

[18] Wang Z X， Wang X G， Dong J Q， et al. Fast resistive reconnection regime in the nonlinear evolution of double tearing modes ［J］. Phys Rev Lett， 2007， 99（18）： 185004.

[19] Canuto C G， Hussaini M Y， Quarteroni A， et al. Spectral methods in fluid dynamics ［M］.Cham： Springer Berlin Heidelberg， 1988.

[20] Canuto C G， Hussaini M Y， Quarteroni A， et al.Spectral methods： Evolution to complex geometries and applications to fluid dynamics ［M］.Cham： Springer Berlin Heidelberg， 2007.

[21] Maître O P L， Knio O M.Spectral methods for uncertainty quantification： with applications to computational fluid dynamics ［M］.Dordrecht： Springer Dordrecht， 2010.

[22] Glasser A H， Sovinec C R， Nebel R A， et al. The NIMROD code： A new approach to numerical plasma physics ［J］.Plasma Phys Control Fusion， 1999， 41（3A）： A747.

[23] Lütjens H， Luciani J F.The XTOR code for nonlinear 3D simulations of MHD instabilities in tokamak plasmas ［J］.J Comput Phys， 2008， 227（14）： 6944-6966.

[24] Wang S， Ma Z W.Influence of toroidal rotation on resistive tearing modes in tokamaks ［J］.Phys Plasmas， 2015， 22（12）： 122504.

[25] Sun W， Wang J， Wei L， et al.Application of Galerkin spectral method for tearing mode instability ［J］. Chin Phys B， 2022， 31（11）： 110203.

[26] Bi H L， Wei L， Fan D M， et al.Excitations of tearing mode and Kelvin-Helmholtz mode in rotating cylindrical plasmas ［J］.Acta Phys Sin， 2016， 65（22）： 225201.

[27] Harris E G.On a plasma sheath separating regions of oppositely directed magnetic field ［J］.Nuovo Cim， 1962， 23（1）： 115-121.

[28] Wang Z， Dai Z， Wang S.Nonlinear excitation of zonal flows by turbulent energy flux ［J］.Phys Rev E， 2022， 106： 035205.

基金资助

国家重点研发计划(2022YFE03040002)

四川省科技项目(2025ZNSFSC0067)