Fusion Res.
Res. SERIES,SERIES,Vol.
Vol.6 (2004)
6 (2004)395–398
000–000
J. Plasma Fusion

Development of the Robust MHD Code Using CIP Scheme
NISHINO Nobuhiro
Hiroshima University, Hiroshima 739-8527, Japan(Received: 9 December 2003 / Accepted: 17 June 2004)

Abstract
A robust MHD code was developed to study the behavior of hot plasma with steep density and/or pressure gradient. MHD equations included in thermal conductivity, viscosity and resistivity which are strongly dependent on plasma
parameters had to be solved; using splitting fractional method an advection (convection) equations were solved by theConstrained Interpolation Profile (explicit-CIP), and Poisson equations in the non-advection terms were solved by
Alternating Direction Implicit scheme (ADI). Total time step was marching by Runge-Kutta method. The code wasexamined for simple MHD shock tube problem. The differences of the magnetic compression and the density ratio
ahead and behind the shock front from the calculation results were within a few percents when the density and/or pressure ratios were less than 103 and 105, respectively.
Keywords:
steep density gradient, steep pressure gradient, robustness, CIP scheme, MHD code

problem with wide temperature range and it had the limitation to use.
In this paper, we have developed new code to treat thevariable plasma properties. Using splitting fractional method
we introduce a linear interpolation at the shock front to solve
the advection equation stable and we apply semi-implicit
method to solve the pressure equation, also full implicit
method to solve the viscosity and thermal conduction equations because of numerical stability.

1. Introduction
Recently pellet injection and/or super-sonic gas jet areused as fuelling methods of hot plasma in large experimentaldevice. These fuelling may change the characteristics of plasma confinement. Therefore, it is interesting to study thebehavior of hot plasma under these fuelling conditions. To
solve these phenomena numerically the code has to include
in the complex atomic processes and it can treat non-thermalequilibrium.
As already known, fluids model such as MHD only
treats quasi-thermal equilibrium states. Therefore, MHD code
has the limitation to solve these drastic states, and it cannot
treat the transition phenomena like plasma just after thickfuelling. However, in general the duration of the transition
processes are shorter than that of interested macroscopic
plasma behavior e.g. a “snake” after pellet injection.
Therefore, the macroscopic density diffusion process inwhole plasma after transition processes can be treated as
MHD phenomena. For this purpose a MHD code has to havethe robustness in a large density and/or pressure gradient(ratio) within a few calculation cells.
We have already developed a single fluid MHD code [1]using normal explicit constrained interpolation profile
scheme (explicit-CIP) [2]. CIP has robustness under severe
conditions such as density ratio of 103 within a few cells.
This code could treat ideal MHD, however, thermal conductivity, resistivity and viscosity were constant in space andtime, and viscosity and thermal conduction were solved by an
explicit method. Therefore, it could not be applied to the2. Governing equations for a single fluidmodel
The governing equations of a single fluid for MHD
model are as follows,∂ρ+ ∇· ( ρ u) = 0,∂t∂uρ+ ρ u ·∇u = −∇p + j × B + ∇⋅τ ij ,∂t∂p= −∇· ( ρ u) − (γ − 1)[p∇· u + η j · j + φ ] − γ ∇· q,∂t⎞⎛⎛
1 ∂u ⎞
1
τ ij = 2µ p ⎜e ij − δ ij k⎟ , φ = 2µ p ⎜e ij e ij − (∇ · u)2⎟ ,3 ∂ xk⎠
3⎠⎝⎝

e ij =1 ⎛ ∂ ui ∂ uj⎞+⎜⎟ , p = ρ Rp T , q = − λ p∇T ,2 ⎝ ∂ xj ∂ xi ⎠∂B= − ∇ × E, E = − u × B + η j, µ 0 j = ∇ × B,∂t©2004 by The Japan Society of Plasma
Science and Nuclear Fusion ResearchCorresponding author’s e-mail: nishino@hiroshima-u.ac.jp1395

Nishino N., Development of the Robust MHD Code Using CIP Scheme(ρ, p: density and pressure, u: velocity [vector], τij :stress tensor, i, j and k are indices of space, µp: viscosity ofplasma, Rp: gas constant for plasma = 2Rgas, γ : ratio of specific heats η : resistivity, q: thermal flux [vector], λp: thermal
conductivity, Φ : dissipation function due to stress in a fluid).
These equations are derived from the mass, momentum
and energy conservation laws, and Maxwell equations. Using
splitting fractional method the advection equations are solved
first by CIP method [1], and Poisson equations in the nonadvection terms are solved by ADI method. ADI is unconditional stable in two dimensions [3]. The time evolution ofone full step is marching by Runge-Kutta method. The orderof time accuracy is chosen 1 to 4 as an input parameter; usually we choose second or fourth order. Thermal conductivity
depends on the magnetic field direction (anisotropy), otherproperties are assumed to be only a function of plasma
parameters (isotropy). The plasma properties used in the codeare described in ref. [4].

3. Numerical test
The CIP and ADI routine were tested independently,
and after that we combined them to make full-MHD code.
The completed code was two-dimensional. For numericalexample, we choose a Riemann problem of ideal MHD shock
tube, which was introduced by G. Sod [5] firstly. The solution of this problem consisted various MHD shocks and it
was already analyzed numerically by some groups [6-8]. The
initial value of this problem is as follows,
The dimensionless initial left and right states, which aredenoted by l and r areρl = 1, ul = 0, pl = 1, Bx = 0.75, By = 1, γ = 2, andρr = 0.125, ur = 0, pr = 0.1, Bx = 0.75, By = −1, γ = 2,respectively.
Initial discontinuity is located in the middle of computational domain at the direction of X, and all states are uniform
in the direction of Y. Using this 1D initial states, numericalsolution, which were obtained for Nx = 800 and Ny = 5 gridFig. 1(a)Fig. 1(b)Fig. 1(c)Fig. 1(d)Fig. 1 Calculational results of a coplanar MHD Riemannproblem, (a) density, (b) pressure, (c) velocity normal
to the shock front, (d) velocity parallel to the shock
front and (e) magnetic field parallel to the shock front.
Grid numbers of calculation domain are Nx = 800, Ny= 5, and horizontal axis shows the grid of X-direction.
The results of 1000 time steps at CFL = 0.4 are shown.

Fig. 1(e)2396

Nishino N., Development of the Robust MHD Code Using CIP Scheme

points with ∆x = 1, ∆y = 1, CFL = 0.4, are shown in Fig. 1
after 1000 time steps; it is necessary to meet CourantFriedrihs-Lewy (CFL) condition due to explicit-CIP scheme.
The solution consists of following waves; the waves moving
to the left are a fast refraction wave (denoted by FR in thefigure), and a compound wave (SM), which consists of aintermediate shock and attached to it a slow refraction wave;
the waves moving to the right includes a discontinuity (C), aslow shock (SS), and a fast refraction wave (FR). The
obtained result is very similar to that of ref [6-8]. Using LaxWendroff and/or MacCormack without artificial viscosity orTVD we could not solve these waves precisely. SM and SS
are sharply resolved by a few points, however, the discontinuity is resolved by more than 10 points, and between SM
and C some numerical oscillations are found in this solution.
It seems that the numerical oscillations are due to the pressure term, which is solved by semi-implicit method this time.
To avoid these oscillations we have to remake this using afull implicit method such as CIP-combine and unified procedure (C-CUP) [9]. This is very near future work.
To simulate dense plasma diffusion such as MHD
behavior after pellet ablation, the robustness for high density/pressure gradient of this code has to be demonstrated. We
have solved high Mach number problems as above. Thisshock waveform is the same as the shock tube solution ofneutral gas. Choosing γ = 5/3, Bx = 0, and By = 1.0, we cancompare magnetic compression with the density ratio ahead
and behind the shock front; two rates should be the samevalue at the ideal perpendicular magnetic shock. Various initial states to create the shocks are shown in a Table 1. Also,
the density ratio and magnetic compression ahead and behindthe shock front from calculation results are shown in the rightcolumn. The magnetic compression is limited to the range 1< B2 / B1 < (γ + 1)/(γ − 1), where B1 and B2 are the magnetic
field strength ahead and behind the shock front, respectively,
i.e. B1, B2 = Bx2 + By2 . For γ = 5/3 the upper limit of thiscompression rate is 4. We tried to confirm the effect of theplasma properties on the shock formation. Calculating thermal conductivity, viscosity and resistivity, which depend on
the plasma parameters, we can also compare the density ratiowith the magnetic compression ahead and behind the shock
front. These cases are denoted by “yes” in each plasma properties of Table 1, while the other cases are denoted by “no”,
in which thermal conductivity, viscosity and resistivity areconstant in time and space.
Usually Mach number represents the shock strength inneutral gas, however the magnetic compression also shows
the shock strength in MHD case. Comparing the results ofno.3 and no.4 in the Table 1, it is found that the shock
depends on not only the upper stream plasma parameters but
also downstream plasma parameters because of the soundspeed of downstream region. In each case the magnetic compression is good agreement with the density ratio, when themagnetic compression is below than 1.5−2.0. These results
are denoted by “good” in the status column of Table 1.
However the discrepancy between the density ratio and themagnetic compression becomes large when the magnetic
compression is large. Numerical oscillations appear to follow
behind shock as the magnetic compression is approaching theupper limit 4; moreover the magnetic compression is over 4
in one case. Those results are denoted by “poor” and “bad” inthe status of Table 1, respectively. Although the calculation
results were broken in those cases, the numerical results haveno negative density and pressure value and the errors in theTable 1 are also within a few percent. The limitation of valid
results region of this code would be as follows,
With the adjacent grid the density ratio is below 103, andthe pressure ratio is than 105 corresponding to the density
ratio.
There is no need to concern above limitation as far as
we treat typical experimental plasma parameters of magnetic
confinement device. Therefore, it is expected that this code
can treat the behavior of plasma with high density locallysuch as after pellet ablation and/or supersonic gas injection.
There is a slight difference in error between the resultsof no.15 and no.16. Initial condition of the temperature is 100
eV, therefore the viscosity, which is proportional to the meanfree path is higher than other cases. In general the effect ofplasma properties on the shock formation was not foundclearly in all cases, this was good while it means the calcula-

Table 1 Various initial conditions for shock calculations.

3397

Nishino N., Development of the Robust MHD Code Using CIP Schemetion of thermal conduction, viscosity and resistivity does not
do anything bad. In reality the shock should disappear itself
with increasing time marching due to dissipation. For example, when high temperature plasma would be selected in theshock calculation, it might be possible to confirm the effectof thermal conduction and viscosity even though we needhuge computational domain and time steps.

necessary for this code. This is near future work.

Acknowledgements
This work is partially supported by JAERI collaboration
program.

References[1] N. Nishino, J. Plasma Fusion Res. Series 5, 541(2001).
[2] T. Yabe and T. Aoki, Comput. Phys. Comm. 66 219(1991).
[3] W.H. Press et al., Numerical recipes in C++(Cambridge 2002).
[4] S.I. Braginskii, Reviews of plasma physics 1, 205(Consultants Bureau, 1966).
[5] G. Sod, J. Comput. Phys. 27, 1 (1978).
[6] M. Brio and C.C. Wu, J. Comput. Phys. 75, 400 (1988).
[7] D. Ryu and T.W. Jones, Astrophys. J. 442, 228 (1995).
[8] G.S. J and C.C. Wu, J. Comput. Phys. 150, 561 (1999).
[9] T. Yabe and P.Y. Wang, J. Phys. Soc. Jpn. 60, 2105
(1991).

4. Conclusion
The robust MHD code was developed. Although thepressure term was solved by semi-implicit method, this code
could handle the large density and/or pressure gradient (ratio)
within a few calculation cells without any negative densityand pressure.
Test results of MHD shock tube calculation showed thatit could calculate MHD behavior under the initial density
ratio of 103 and pressure ratio of 105 within a few percent
error.
It was found that there are numerical oscillations
between severe states; therefore, the improvement would be4398