MHD现代算法介绍 张嫚 2011.03.24主要内容1 MHD算法简介2 MHD算法存在的主要问题3 MHD在太阳风模拟中应用1 MHD算法简介理想MHD方程 迎风差分格式Lax单步差分格式 迎风差分格式Lax单步差分格式是一阶精度，耗散效应很强，激波会被拉宽或抹平，不存在色散效应，因此不会出现振荡Classical schemesLax-Wendroff schemeMacCormack schemeLax-Wendroff 和MacCormack都是二阶精度，耗散效应很弱，存在色散效应，激波附近会产生振荡Others like Warming-Beam, Rusanov,etc.Riemann Solvers 从流动的物理特性出发，采用间断分解概念求解流动方程组，在相邻的网格单元之间形成一个局部Riemann问题，转化为求解黎曼问题，在此基础上发展了一系列Riemann solvers流动量在每个网格内以常数分布求解过程物理概念清晰，能较好反应流动特性，计算效果较好，但是该算法只有一阶精度。Godunov Riemann solvers发展高阶精度的Godunov 是必要的，MUSCUL和PPM差分算法是这方面的典型 要获得黎曼精确解是非常复杂的，并且很大程度上也没有必要，因此数值方法上很少涉及到精确解，广泛应用的是近似的黎曼解，下面介绍几种常用的近似黎曼解Roe (linearized) type solvers 假定间断面两侧流动量常数分布，构造一个合理的线性近似系数矩阵 来代替非线性Jiacobian系数矩阵 把非线性矢量方程转化为线性矢量方程，求出近似解HLL (non linear) type solvers These non-linear solvers are based on approximatingthe wave structure of the full Riemann problem by a simplified wave structure.HLL Riemann solver 这种方法可以保持正压，保持正压在low-beta是非常重要的，但是具有较强的耗散性，在此基础上发展的HLLC和HLLE克服这一缺点HLLC Riemann solver The HLLC flux based on the approximate Riemann solver is a modification to account for the shortcoming of the HLL flux offset the influence of intermediate wavesHLLD Riemann solverHigh resolution comparable to Roe schemeHigh resolution comparable to Roe schemeHLLE, HLLEM .etc2 MHD算法存在的主要问题保正性MHD模拟的困难之一就是保持压力为正，由于MHD总能量包含三个分量：内能，磁能和动能，在low-beta区域，磁能相比于动能很大，因此，在计算总能量，动能和磁能时产生的离散误差很大程度上很导致负压的产生，出现非物理解。为了保持正压性，我们在局部区域用熵密度方程或者压强方程来代替能量方程求解MHD方程，为此我们需要用开关函数来判断可能出现负压的地方参考（Dinshaw S. Balsara* and Daniel Spicer. Jcp)给出开关函数It consists of 3 sub-switches which includes or excludes some flow conditions. and it is switched off in other situations,Switch 1 (SW1 for short) is designed to indicate whether the pressure is very low. So SW1 is on when Switch 2 the second switch SW2 is intended to exclude the regions of strong magnetosonic shocks. This is accomplished by ensuring that the local minima and maxima pressure are not too far apart. It is on whenThe third switch SW3 is used to exclude the regions where mildly compressive motions(also strong compressive motions) may take place. It is on when can be adjusted in different situations. Typically, we use and the pressure equation is used only when all of the three switches are on保持磁场散度为零MHD不同于Euler方程的就是存在磁场效应，如何处理磁场使之散度为零也是求解MHD的一个重要方面，有很多种方法可以用来处理磁场 Projection method， 8-wave formulation8-wave formulation， Constrained Transport (CT) methodConstrained Transport (CT) method， Central differenceCentral difference Diffusive controlDiffusive control，GLM-MHDGLM-MHD， Vector potential Vector potential Hyperbolic divergence cleaning Hyperbolic divergence cleaning Splitting-based schemesSplitting-based schemesProjection method在每个计算一个时间步之后对磁场进行更新8-wave formulation8-wave formulation在方程右边加入Powell源项The eight-wave scheme is computationally inexpensive, easy to add an existing code, and quite robust,Constrained Transport (CT) methodConstrained Transport (CT) method采用交错网格方法，计算面上的磁场，之后将其投影到网格中心处 这种方法是严格守恒的，可以通过不同的构造电场磁场方法达到高阶精度，但是延伸到一般的网格上比较困难，不能直接应用于AMRCentral differenceCentral differenceCD算法是在简化CT算法基础上发展的，在时间和空在时间和空间上都是二阶精度间上都是二阶精度Diffusive controlDiffusive controlModifying the magnetic equation by adding a source term proportional to the gradient of divergence so that the divergence satisfies an advection-diffusion equationThis technique, referred as diffusive control of divergence , has the same advantages and disadvantages as eight-wave approach, but appears to keep the level of divergence lower than the eight-wave approachGLM-MHDGLM-MHD采用9个方程来控制divergencehyperbolic correctionHyperbolic divergence cleaning显式的抛物型方法中 的取值受到稳定性条件的限制，结果不是很满意，为了克服这一缺点，将抛物型和双曲型方法结合起来Vector potentialVector potentialuse vector potential A as the primary variable instead of Buse vector potential A as the primary variable instead of B，and the divergence free property is automatically and the divergence free property is automatically ensuredensured3 MHD在太阳风模拟中应用球坐标下太阳风模拟 目前太阳风模拟大多是在直角坐标系下进行，计算量比较大，直接应用于AMR比较困难，鉴于计算区域为球面网格划分，所以考虑在球坐标下实现太阳风模拟，达到稳定所需时间比较短，关于磁场处理部分，用一种可以直接应用于AMR的CT方法实现，不用求解Possion方程，大大减少了计算量方程形式直角坐标下 球坐标下考虑太阳自传，Powell 源项，加速加热，磁场分裂模拟区域 有限体积法，将流场和磁场分开处理，磁场部分用CT方法，采用一种新 central scheme 的方法，不需要求解黎曼问题，也不需要求解特征值和特征向量，该方法适应于任何坐标，半离散形式，时间上用两步 Runge-kutta 法达到二阶精度数值方法符号说明流场重构用线性重构达到二阶精度流场离散形式磁场重构 满足二阶精度，其余方向上的重构类似磁场离散形式CT方法在面上对磁场进行离散，左边均为面上的平均值，之后有面上的值可以得到体积平均的磁场计算区域考虑整个球面时，在极区离散方程会存在奇异性，需要特别处理流场流通量计算每个小单元被分割为27部分，以边和脚点为中心的矩形单元，以面为中心的非矩形单元，剩余的内部区域 分割依据：在 时间内，由边界上产生的波并没有影响到内部区域这样每个小单元的边界上都是光滑的，边界上流通量的计算直接求积方法每个小单元上面27个区域同时进行更新，之后把结果投影到单元中心上，得到新时刻体平均上的流体变量的值电场计算类似于流场，用周围四点的值投影到计算区域边沿速度的确定源项处理 不能简单用体心或者几何中心的值来近似体平，那样不能达到二阶精度，由于三角函数的存在，源项众不同部分取单元中不同位置的值，以此可以写出不同源项的表达式，这里不再一一列举 当源项中包含有Powell源项时，需要考虑divergence B在极区的奇异性，另外就是用CT处理磁场时，需要计算divergence B在面上的表示初始输入用太阳表面参数进行无量钢化，利用太阳风Parker解偶极场或者螺旋磁场作为初值偶极场螺旋磁场边界条件在半径方向上采用固定边界条件