1无碰撞激波的MHD模拟方法

无碰撞激波的理想MHD模拟国家自然科学基金重点项目(10135020)和国家863国家自然科学基金重点项目(10135020)和国家863高技术惯性约束聚变主题资助陆启明，杨维纮

(中国科学技术大学近代物理系，合肥230027)［摘要］本文提出了通过求解理想MHD方程模拟空间无碰撞激波的方法，并且使用该方法模拟了垂直无碰撞激波与行星际反向磁场结构和高密度等离子体团的相互作用过程，同时与粒子模拟的结果进行比对。两者的结果非常类似，而且在MHD模拟中得到了一些粒子模拟中没有观察到的现象。模拟结果表明了理想MHD模拟是准确且可行的，同时相对于粒子模拟又有很好的计算效率，便于扩展至高维的情形。［关键词］无碰撞激波MHD方程数值模拟0.引言行星际空间和星际空间中充满着完全电离的稀薄等离子体，粒子平均自由程非常大，经典库仑碰撞效应往往可以略去。这些无碰撞的等离子体通常以超声速运动，形成太阳风和星风。当太阳风和星风遇到存在磁场的行星和恒星的阻挡时，在界面处将形成各种间断面，如地球的磁顶层和弓激波、太阳系的日球顶层等［1］。以弓激波为例，观察资料表明，弓激波是无碰撞的激波，上游是未扰动的超声速太阳风，而下游的等离子体以亚声速绕过地球的磁顶层［2］。无碰撞激波是等离子体物理、空间物理和天体物理学中的重要基础性课题，对它的深入研究有助于了解激波本身的产生、演化、耗散机制以及各种行星际结构与激波的相互作用问题。同时，随着计算机技术的不断发展，使得通过计算机来模拟无碰撞激波成为可能。对无碰撞激波的数值研究通常可分为两类：粒子模拟和MHD模拟。前者通过求解Maxwell方程和每个粒子的运动方程得到无碰撞激波的结构及伴随的物理过程，如文献［3］种介绍的方法。如果模拟足够多的粒子，这种方法可以较好的再现激波的结构。但限于当前计算机的处理能力，如果要在大尺度和高空间维数情况下使用粒子模拟比较困难。与之相对应的MHD模拟方法是通过求解MHD方程来模拟无碰撞激波的结构和各种物理过程,在高空间维数的模拟中有着较高的计算效率，容易用当前的计算机条件来达到。本文第一节中讨论了MHD模拟无碰撞激波的方法,第二节中使用该方法模拟了垂直无碰撞激波与行星际空间结构的相互作用并与粒子模拟的结果进行了对比。1.无碰撞激波的MHD模拟方法忽略了粘滞力和电阻的理想MHD方程组如下:^=-V・(pv)^=-V・(pv)atQ(pv)at( BB)=—V・pvv+p*-——0聲=Vx(vxB)(y1)ly-12丿(vxB)xB1）2）3）4）其中的p，p，v,B和£分别表示流体的密度、压强、速度、磁场和总能量，p国家自然科学基金重点项目（10135020）和国家863高技术惯性约束聚变主题资助国家自然科学基金重点项目（10135020）和国家863高技术惯性约束聚变主题资助2卩0B2 py为绝热系数，£=pv2+ + 。另外磁场同时必须满足散度为0：V-B=0。2卩y-10通常可用L（参考尺度），a（自由流离子声速），p （自由流密度）将方程（1）~884）无量纲化［4］。另外磁场可用4）无量纲化［4］。另外磁场可用归一化，此时原方程组中便不含％。本文以下讨论的是无量纲化的MHD方程。在一维情况下将上述方程展开可得：

5)6)7)U二5)6)7)tx其中：U=Qpv,pv,pv,B,B,B,JxyzxyzF=(pv,pv2+p*,pvv一BB,pvv一BB,0,vB一vB,xxxyxyxzxzxyyxvB-vB,v1+p*丿-B(vB+vB+vB))txzzxxxxxyyzz用数值方法求解以上方程组时采用本征波分解的方法[4],[5],[6],。理想磁流体方程包含了三对特征波，分别是快波、中间波邙阿尔芬波)、慢波，它们的特征速度分别为C、c、c。fas另外再加上信息波(entropywave)和磁流波(magnetic-fluxwave)，特征速度v。这样理x想磁流体方程可以分解为8个波组成的本征系统，其本征值分别为：九=v+c,九=v+c,九=v+c,九=v (8)1,8xf2,7xa3,6xs4,5x同时对应每一个本征值有一个左本征矢量L和右本征矢量R，且满足L.R.=<5,s s iji,j使用波分解方法求解MHD方程，首先根据U.和计算x.1/2处的雅可比矩阵TOC\o"1-5"\h\zii+1 j+1/2A =dF 壬U (9)j+1/2 j+1/2； j+1/2并得到左本征矢量L.1/2和右本征矢量R.1/2。然后用左本征矢量将网格点上的F.和U.投射j1/2 j1/2 j j到左本征矢量空间中的F.s和U.s：jjFs=L -F Us=L-U (10-11)j j+1/2 j j j+1/2j其中a=maxsFs+=其中a=maxsFs+=j=max6sj+1/2 j+1/21(Fs±aUs)2j sj)。12)此时可使用迎风格式[4],[7]、ENO格式[8]或者改进的WENO格式[9]等能捕捉激波结构的高精度格式分别计算正、负流半格点处的值，然后再次合并：Fs =Fs++Fs- (13)j+1/2 j+1/2 j+1/2用右本征矢量将合并后的值还原至原空间：TOC\o"1-5"\h\zF Fs・Rs （14）j+1/2 j+1/2 j+1/2s=1时间推进可采用四阶Runge-Kutta方法［10］。（5）式的方程中定义空间差分算符L（U）=-丄G -F）Axj+1/2 j-1/2K=At・L（U）00K=At・L（U）=At・LfU+-K11（020丿K=At・L（U）=At・LfU+-K22（021丿K=At・L（U）=At・L（U+K）3302Ut+1= -（K+2K+2K+K）j 0 6 0 1 2 3无碰撞垂直激波的两个数值实验本节采用以上描述的理想MHD方法模拟了一维空间无碰撞激波的产生、演化以及行星际空间结构与无碰撞垂直激波的相互作用过程。计算时采用WENO(WeightedENO)五阶高精度格式［9］捕捉激波的空间结构，时间上采用四阶Runge-Kutta方法。另外计算中的时间步长根据空间步长实时计算得到：15）0.8Ax15）At= /—maxIv+cx在这里max为沿所有的空间格点取最大值，Cf为快波的速度，CFL值取0.8。空间格点取均匀分布，共400个格点。2.1无碰撞激波我们设定如下的初始条件以产生激波：上游(xvO)： p=0.25 v=1.0 v=0.0 v=0.0 B=0.1 B=0.0 p=0.01x y z y z下游(x>0)： p=0.25 v=0.0 v=0.0 v=0.0 B=0.3 B=0.0 p=0.01x y z y z首先取Bx=0.1,另外丫=2。该条件下上游的离子声速C=0.28,上游流的离子声马赫数大约为Mc=3.5。当高速粒子流到达静止的等离子体时，两者接触的间断面上即产生压缩等离子体激波。在该初始条件下,由于磁压与粒子流压力平衡,激波的激波面保持在初始的间断面处，便于观察激波的发展情况。图1是以上初始条件下经过t=0.8的MHD方程的解。图中可以看到一系列MHD系统的特征波：一个右行的快波FR(对应波速为v+c)；xf两个慢波SR1和SR2(对应的波速分别为v土c)，由于该条件下有c沁c，因此SR1和xsasSR2中同时复合了阿尔芬波；接触间断面C,为初始间断面随着磁流体以v运动，因此介x于SR1和SR2之间；以及我们关心的激波SF,该激波为快模式的斜激波，上游激波角8=n/4。在激波参考系下，上游流速超声速、超快波波速，经过激波面后速度降为亚快波波速，同时下游处的磁场强度、等离子体热压强以及等离子体密度均增加。激波面的厚度大约为3~4个空间格点的距离，与文献［7］、［9］中的值接近。以下在模拟行星际空间结构与激波的相互作用中，为使结果相对简单且易于观察，我们取B=0，此时C=C=0，因此该系统中只剩下右行的快波FR,切向间断面C以及快模式的x as垂直激波SF。

p-0.5-1.0-0.5X-1.0-0.50.00.5X0.10.0-0.1Xr-1.01.0-1.0-0.50.00.5X1.0p-0.5-1.0-0.5X-1.0-0.50.00.5X0.10.0-0.1Xr-1.01.0-1.0-0.50.00.5X1.0图1t=0.8时MHD方程的解2.2垂直激波与反向磁场结构的相互作用行星际磁场扇形边界两侧磁场呈相反方向，而行星际磁云除磁场强度异常增强外，磁场矢量相对于黄道面的倾角呈单调的旋转，当他们扫过地球时，局域行星际磁场将发生由南到北或由北到南的转变［3］。以下模拟了这种反向磁场的结构与垂直激波的相互作用。从t=0时刻开始，激波的上游边界处磁场强度逐渐变小、反向，直至t>2tc时，磁场大c小为其初值而方向相反。［（1一t/t）B,0<t<2tB=<cy c一b, tn2tyc上游的反向磁场结构随着上游等离子体向激波方向运动，各个不同时刻的磁场以及等离子体密度分布如图1、图2所示：图中激波面始终维持在x=0附近。可以看到t=0.5时上游已经形成了一个完整的反向磁场结构，并继续向下游运动。在t=0.9时开始穿越激波面，t=1.25时反向的磁场已经到达下

游，并开始出现一个负的尖峰，然后逐渐扩大至整个下游区域，表示反向的磁场结构可以顺利穿越激波面到达下游。磁场的模拟与文献［3］中的结果非常类似。图2磁场分布图0.80.2-1.0-0.50.00.51.00.60.4-0.50.00.51.0-1.0-1.0-0.50.00.51.0-1.0-0.50.00.51.00.80.2-1.0-0.50.00.51.00.60.4-0.50.00.51.0-1.0-1.0-0.50.00.51.0-1.0-0.50.00.51.0图3密度分布图该初始条件下形成的激波为快模式激波［2］，特点是经过激波面后密度增加、压强增加和磁场强度增加。负的磁场结构经过激波面后同样会增加磁场的强度，而方向不变，因此经过激波面后磁场跳变，形成了一个负的尖峰，并逐渐扩大。而由于上游区域磁场随时间变化产生了小的扰动，由方程（3）可知流体速度也会因此有一个扰动，并因此带动了上游区域密度的变化。如图2所示，上游的等离子体密度由于磁场随时间的扰动而形成了一个小的密度峰，并伴随着反向的磁场结构继续向下游运动。并在t=0.9开始穿越激波面。由于激波的压缩作用，使得上游的密度扰动在下游发展成一个较高的密度峰，并且继续随着反向的磁场结构向下游运动，直至移出计算区域。文献［3］的模拟结果同样可以看到下游产生一个较高密度的峰值，但是上游的密度扰动淹没于背景粒子随机扰动中，无法观察出来。另外文献［3］的结果中可观察到经过反向磁场结构后下游的剧烈湍动现象，而本文中的MHD模拟则没有这种现象。原因可能是等离子体经过激波后，由于部分质子被反射，导致下游的等离子体偏离麦氏分布，造成等离子体动力学上的不稳定。而MHD方程默认粒子、电子服从麦氏分布，无法描述该情况下的湍动现象。2.3垂直激波与行星际高密度等离子团的相互作用太阳等离子体输出的另一种形式是行星际高密度结构，它可能具有环状磁场位形，有时还伴有高速流。这里我们模拟这些高密度等离子体团与垂直无碰撞激波的相互作用。类似于文献［3］,设t=0时上游有一个高密度的突起，峰值密度为上游背景值的2倍，其余磁场、等离子体热压强等都不变。这个密度峰随着上游流向激波面运动，在t=0.7时到达激波面，然后穿过激波面继续向下游运动。由于激波的压缩作用，该密度峰经过激波面后变陡，密度峰值上升，但对其他的参数如磁场强度、压强、速度等都没有太大的影响。因此高密度等离子

团同样能够穿越垂直无碰撞激波，而本身性质没有太大改变。相互作用的具体过程见图4、图5。其中的SF为激波面，FR为右行的快波，C为切向间断面。另外我们发现，在高密度等离子团与激波相互作用过程中，从激波面处激发出了一个新的波FR2,并快速的穿越接触间断面向下游运动，如图4所示。从该波传播的速度和特征来看FR2为一个右行的快波，但该波形成的具体机制有待于进一步的研究。对比粒子模拟的结果中未发现有这个波存在。同样在该模拟结果中没有观察到下游的湍动现象，因此也能间接证明这种湍动是由于离子偏离麦氏分布的动力学不稳定性引起的，无法在MHD模拟中看到类似的现象P1.21.0_t=0.71JI1.0-1.0-0.50.00.51.0P0.60.2X图4密度分布图P1.21.0_t=0.71JI1.0-1.0-0.50.00.51.0P0.60.2X图4密度分布图X图5磁场分布图0.60.6X图5磁场分布图结论相对于无碰撞激波的粒子模拟，本文使用了求解理想MHD方程的流体模拟方法。并应用该方法模拟了文献［3］中讨论的两种行星际结构与垂直无碰撞激波的相互作用过程来作为与粒子模拟的对比。两者的模拟结果十分相似，MHD模拟下同样能够观察到粒子模拟的大多数现象：反向磁场结构与高密度等离子体团可以穿过垂直无碰撞激波，而且本身性质没有太大的改变，同时激波也能保持原先的状态。两者的对比可以说明在对垂直无碰撞激波的模拟中，粒子模拟与MHD模拟都是可行的，能够得出相同的结果。同时，在MHD模拟的结果中还发现了一些粒子模拟没有观察到的现象，如反向磁场结构与激波相互作用中上游的密度扰动，高密度等离子团与激波作用过程中激发出一个新的快波等。另外，在相同的模拟尺度下，MHD模拟有着相对较高的精度（与计算用的格式有关），而且耗用的计算机时比较小，容易拓展至二维或者三维的情况。MHD模拟的缺点是无法反映粒子分布偏离麦氏分布时的情况，因此在本文的模拟结果中均没有出现经过行星际结构后下游的强烈湍动现象。参考文献1.王水,陆全明.空间无碰撞激波的数值研究.天文学进展,1997,Vol.15No.3:218-2302.MargaretGetal.《太空物理学导论》.科学出版社2001.83.陆全明,李毅,王水.行星际结构与垂直无碰撞激波的相互作用.空间科学学报,1997,Vol.17No.3:200-205KennethG.Powell,PhilipL.Roe,TimurJ.Lindeetal.ASolution-AdaptiveUpwindSchemeforIdealMagnetogydrodynamics.JournalofComputationalPhysics,1999,Vol.154:284-309NecdetAslan.MHD-A:AFluctuationSplittingWaveModelforPlanarMagnetohydrodynamics.JournalofComputationalPhysics,1999,Vol.153:437-466K.Sankaran,L.Martinelli,S.C.Jardinetal.Aflux-limitednumericalmethodforsolvingtheMHDequationstosimulatepropulsiveplasmaflows.Int.J.Numer.Meth.Engng2002,Vol.53:1415-1432M.Brio,CCWu.AnUpwindDifferencingSchemefortheEquationsofIdealMagnetohydrodynamics;JournalofComputationalPhysics.Vol75:400-422HartonA,EngqmistB,OsherSetal.Uniformalyhighorderaccurateessentiallynon-oscillataryschemeIII.JournalofComputationalPhysics.Vol49:231-303Guang-ShanJiang,Cheng-chinWu.AHigh-OrderWENOFiniteDifferenceSchemefortheEquationsofIdealMagnetohydrodynamics.JournalofComputationalPhysics,1999,Vol.150:561-94SigalGetal.StrongStability-PreservingHigh-OrderTimeDiscretizationMethods.SIAMReview,2001,Vol.43No.1:89-112SimulationofSpaceCollisionlessShockwithIdealMHDEquationLuQiming,YangWeihong

(Dept.ofModernPhysics,USTC,Hefei230027)AbstractHerewepresentamethodtosimulatespacecollisionlessshockbysolvingidealMHDequation.Andwiththismethod,wesimu