- Email: [email protected]
MAG - two-dimensional resistive MHD code using an arbitrary moving coordinate system
Computer Physics Communications ELSEVIER
Computer Physics Communications 106 (1997) 76-94
MAG - two-dimensional resistive MHD code using an arbitrary moving coordinate system Oleg V. Diyankov, Igor V. Glazyrin, Serge V. Koshelev Russian Federal Nuclear Center. All-Russian Research Institute of Technical Physics (RFNC - VNliTF), PO. Box 245. Snezhin.~k~ Chelyabinsk Region 456770, Russian Federation
Received 25 July 1996; revised 25 July 1997
Abstract The MAG code used for two-dimensional magnetohydrodynamics flows modelling is described. MAG algorithms are formulated for an arbitrary moving coordinate system. MAG can calculate flows with large deformations contained inside a region whose weakly deformed boundaries maintain the correct description of boundary conditions. Numerical experiments with two-dimensional MHD flows calculations show that the proposed numerical schemes are quite robust and accurate. Sample calculations illustrating the properties of the algorithms are presented. (~) 1997 Elsevier Science B.V. Keywo~s: 2 ½D MHD; Arbitrary moving coordinate system; TVD scheme
1. Introduction
The motivation for the creation of the MAG code was the necessity of modelling hot dense plasmas in a magnetic field. These plasmas are encountered in inertial fusion experiments such as laser-produced plasmas, tokamak plasmas, Z-pinch, etc. Typically, such plasmas have ," large ,,,~:,.,,,u~ D. . . . . la,. number and a large Lundquist number. Strong magnetic fields cause the transport of heat to be several orders smaller than without the fietd. These facts make the problem of code creation more difficult. There is a great variety of approaches in code constructing, which are used now for plasma flows in Megagauss magnetic fields simulation. Among them one can find Lagrangian, Eulerian, and Lagrange-Eulerian codes. Lagrangian codes have the advantage over the Eulerian ones in the simplicity of boundary condition formulation but ordinary Lagrangian codes do not permit to simulate long-term evolution of fluid instabilities because of numerical grid mixing. On the other hand, Eulerian schemes allow one to simulate flow motion with large deformations but there are difficulties with boundary conditions determination. To resolve boundary conditions one should introduce the high magnitude of electrical resistivity in the background non-ionized gas, disused in the region between the boundary and plasma: In this paper the code MAG [ 1,2] for plasma modelling in an arbitrary moving coordinate system at its present stage of development is presented. The use of an arbitrary moving coordinate system allows one to simulate flows with large deformations inside the flow region, conserving the correct description of conditions on its weakly deformated boundaries. This has proved to be highly successfid for strongly shocked flows. 0010-4655/97/517.001~) 1997 Elsevier Science B.V. All fights reserved. Pll SO0 ! 0-465 5 ( 97 )0008 3-0
O.V. Diyankov et al./Computer Physics Communications 106 (1997) 76-94
77
Arbitrary moving computational mesh permits the use of detailed resolution mesh in all main stages of plasma flow formation and reduces numerical diffusion of the magnetic field so that flows with a magnetic Reynolds number of the order of hundreds may be calculated. The code has been used for simulating three types of problems. The first one is the plate ablation under '~aser radiation. It is a well-known fact that during this process a spontaneous magnetic field may appear, and the estimation of its influence is an actual problem. The second problem is the simulation of tokamak plasma-wall interaction. The third one is modelling of Z-pinch implosion. 2. The model The system of equations used in the MAG code is delermined by the Braginskii [3] model for the onetemperature case (electron temperature is equal m ion temperature). The dynamic evolution of the system is described by a solution of the following equations: - continuity:
30 =.r_+ V. (pu) = 0 , at - momentum:
p.
+(u.V)u
(I)
)
+VP=-4~
[B×rotB],
(2)
- magnetic field diffusion: riB=rot O"t-
C2 [uxB]-~rotB
47ro"
c[(,)
+- V e
xVPe
~
)
]
-
c rot ( [rOtziniB x B] ) ~e
(3) '
- energy:
(a£ ) m p ( j £ ) + PVu + Vq p . - ~ - + (u. V)£ ANe (l+Zi) Zi _ _ -
-
j2 +
c
Per"
f
rot-B'~
+ p"
(Qext - Qrad)
(4)
heat flux: q = -~e V'T,
(5)
- Maxwell's equations: ~TB=0,
(6)
j = ~ rotB, (7) 4~ where t is time, u is the mass velocity, p is the mass density, ni =- p/(Aimp ) (with Ai for the atomic weight, and mt, for the mass of the proton), P is pressure, T,£ are the temperature and specific internal energy, B is the magnetic field, j is the current density, q is the heat conductivity flow, o', ~e are coefficients of electrical conductivity and heat condtJctivity, Qext is the energy released, Qrad is the radiation loss (bremsstrahlung
78
O, V. Diyankov et aL /Computer Physics Communications 106 (1997) 76-94
model at the present time), and Zi is the average charge of an ion (at present taken in the form Zi = Z ( p , T ) from paper [4]). Other variables are generally accepted. Assuming an independence of unknown variables (p, T, P, u, B) from the angle ~, in the case of cylindrical geomet~, and from z in the case of plane geometry, one obtains a system of equations, in which all unknown ~ a b l ~ de~nd on 2 spatial v~ables but the mass velocity vector and magnetic field vector have 3 components. This system is solved in MAC; code. Therefore, MAG code is said to be a 2½D MHD code.
2.1. Equation of state and transport coefficients 2.1.1. Equation of state for ideal gas Pressure and specific internal energy of ideal gas are obtained from the equations R
P f ( Zi + I ) "~i PT,
£ =
(Z~ +
1)
R
T
AiT-I
+
5Joniz,
Zi P~fP l+z~
(8)
where R is the universal gas constant, 7 = 5/3, and ~'ioniz is the energy for matter ionization.
2.1.2. Transport coeOfcients Transport coefficients are obtained in two ways: - magnetized transport coefficients for ideal gas from paper [ 3 ]; - widerange coefficients for average charge, heat and electrical conductivity from paper [4]. To account for magnetized effects all coefficients are multiplied by the correspondent factors taken from [ 3 ]. Heat conduction and magnetic diffusion terms are then in sca!~ f o x .
3. B o u n d a D " c o n d i t i o n s
A wide range of boundary conditions is implemented in the code. Three types of boundary conditions are considered.
3.1. Gas dynamics conditions They are: - applied pressure: PIb = f ( t ) , where PIb means the pressure at the corresponding boundary, and f ( t ) is a given time dependent function; - rigid wall: Unlb = 0, where un is a normal to boundary component of mass velocity; - piston: Unjb- f ( t ) ; - incoming flow (only for Eulerian variant of the code)" Ulb = fu(t); £1b = fc(t); Plb -- fp(t), where all variables at the correspondent boundaries are a given function of time.
3.2. Conditions for heat conductivity They are: given temperature, given heat flux, and heat flux as a function of temperature. These conditions may be written in the form
c~VnTjb+ I~l'jb =/(Tjb,
t) •
O.V. Diyankovet aL/Computer Physics Communications106 (1'997) 76-94
79
Here, TIb is the temperature at the corresponding boundary, a,/3 are numerical parameters ( a = 0 and ,8 = 1 for a given temperature and a = ! and/3 = 0 for the heat flux), f(T[b,t) is a given function• VnTIb is a normal to the boundary component of the temperature gradient.
3.3. Conditions for magnetic field These conditions may be written in the form: - symmetry: ~Bz/dn = O, OB~,/an = 0, B~ = 0; - conducting wall: Bz = 0, B~, = 0, Bn = 0; - axis: Bz = O, B~ = 0, B, = O, where B~, is used for the cylindrical symmetry and Bz for the plane one; - given current: Bz = 2¢rj/c - plane symmetry; B~, = 21/rc - cylindrical symmetry, where Bn is the normal to the boundary component of the magnetic field, j is the current density, 1 is the whole current in the z-direction, and r is the upper radius.
4.
Details
of
numerics
Below, the technique of MAG difference scheme derivation will be described• We will omit the detailed formulas of the difference scheme because of their complexity• Let us assign indices x and y in the case of the cylindrical geometry to the r- and z- components of vector variables and r- and z- components of independent spatial variables correspondingly, and index z to the ~ocomponent of vector variables,, Then we obtain a basic system of equations, which depends on the parameter of geometry u ( ~, = 0 for plane geometry, and u = 1 when the geometry is cylindrical). Because of Eq. (6), one can, substitute an equation t0r the z-component of the magnetic vector potemiai A instead of two equations for the x- and y- components of the magnetic field,
Bx = x -v maA 3y
By = - x - v --.3A Ox
(9)
Eqs. ( 1 ) - ( 5 ) have been written in an arbitrary moving coordinate system, and then splitting into two systems has been performed. The diffusive terms have been split into a separate system of equations, and the remained terms produce a quasilinear hyperbolic system. The first system is a hyperbolic one (compare with ( I ) - ( 4 ) ): - Equation of continuity: Opx"v at
lOOx"v/g(uk--vk) k
=0.
(iO)
- Euler equations for different components of velocity: • x component: ap x
u., +
Ot • y
cgPX~'v~(Uk -- Vk)Ux O~~
+
~"aPf ~xvrg O~k
=
i { 3x'v/~B'Bx 4¢r O~k
_
u B 2 1 + ~,pu2 ," ~ "
(11)
component: ap X~'V~ ltv- +
Ot
~pXUv~(Uk--Vk)U, ,
x~OPIT~yv~
1{OxVv~BkB~' 1
(12)
O. V. Diyankov et aL/ Computer Physics Communications 106 (1997) 76-94
80
• z component: Op x 2~'vfg u z
+ af'x2% (uk vk)u
i f Px2"v/'~BkBz 1
-
(13)
at -
Equation for A:
3A + Ot -
8A
( u k - v k) ~
=0.
(14)
Equation for z component of magnetic field:
av~Bz + av'gBz(u k - v k ) ) = aV~U~B~ at a~k o~k "
(15)
- Equation for total energy:
a£/ ' x~v/gat + a~-fi{£Ix~v~(uk - vk) + PYx~j;°uk} = PX u x/rg (Qext - - Qrad) 4- ~ -
1 {aX~v/gB'(Bu)} 3~e
.
(16)
Equation for the square root of the metric tensor:
at
(17)
a~k
The sum is taken over the same upper and lower indices. Here u k is a projection of the mass velocity vector to the kth vector of a local basis, u I' = u. T t, where ~ is a covariant local basis vector, v k, B k are determined in the same way, v is the coordinate system velocity, g is a determinant of metric tensor with covariant component go, ~k are mathematical coordinates, the total energy Ef = p . (£ + ½. [ul 2) + ( ! / 8 7 r ) . [B[2, and the total pressure P / = P(p,£) + (1/87r). [B[2. The second system of equations is a diffusive one, p
:
equation of continuity, - equation for determinant of metric tensor, - Euler equation for velocity.
coast.
-
ffi c o a s t . u
-
-
coast.
(18)
Equation for A:
aAvrg = x ~ 41rtr "-'7 a~ ~kv~x-P " ~ at
(19)
"
- Equation for Bz:
0 { c2 Ot
:Ox~Bz)}
cmrAi e Z
{ 0 (I)3P "~ a~J
- Equation for energy: O~'xJ'Pv~
! XVv~ fgiJx-2V axVBz ax~'Bz
(
a
oP
(20)
O. V. Diyankov et al. / Computer Physics Communicatimas 106 (1997) 76-94
I (_~ ..
d-- "
gtj x - p V ~
g
O~j)2}
trip c9 ( "t- n . Ai .
e
cmraix/"g(OBzx v OP "~41re Zi p d~i O~J
~
8l
Jj.~_,
)
( l + Z i ) . Z i "vc~x"
OBzx" OP) tgscj d~i •
(21)
The flow region is a geometric figure, topologically equivalent to a rectangle, and it is initially covered by regular rectangular mesh. All dependent variables are specified at the center of the numerical mesh cell, besides coordinates which are vertex-centred. The following approach is used for the hyperbolic system solving. The process of the solution of the equation on each time step is broken into two stages. At the first stage the coordinate system velocity is considered equal to the mass velocity (and this stage is called a Lagrangian stage), and the corresponding hyperbolic system is solved, using difference scheme, which was derived from the TVD approach [5], first suggested by Harten in paper [ 6]. At the second stage the grid moving is performed and fluxes between cells are calculated. This approach allows one to simulate flows with rather large deformations inside a flow region on regular grids. The second system (of diffusive type) is solved in the following manner. The equations of magnetic diffusion are written for the z-component of magnetic field and z-component of magnetic vector potential, and the resulting system consists of three nontrivial diffusive equations: equation of heat conductivity; equation of z-component of magnetic field diffusion; equation of z-component of magnetic vector potential diffusion. The main terms of these equations are moved to the left sides and a difference scheme [7], analogous to -
-
-
t k a lkr--rel~olu e~homo [.~1 ~...,~,_
I.~_T~
O-=-='~--.-
~ = .
. . . . .
.
v
~.,
ie , t ~ r l ft~r ~ l u ; n ~ .
.
.
.
. .
.
.
.
.
,,,.v..
_-~,,,~,
.=.~k ,~f' tll~ae~ onn,~ti~ne ~
._
"--""
'~''--"-"--"~
•
~ ; m n l o i t p r n | i ~ n ie neoti I n ro~:nlv~ tho_
",---'.'"~'-P-"
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
nonlinearity of the right sides and the ICCG method [ 9] is used to solve linear algebraic systems. Five algorithms are used for coordinate system moving. They are Lagrangian, when the coordinate system velocity equals the mass velocity; Eulerian, when the coordinate system velocity equals to zero; local, when the coordinate system velocity at a mesh vertex is determined using mass velocities in the 8 neighbour sites; - algebraic, when the coordinate system velocity at a mesh vertex is determined using mass velocities on the boundaries of the calculated area; Poisson-type, when the coordinate system velocity field is determined by solving a Poisson equation, the boundary conditions for which are determined by the mass velocities on the calculated area boundaries. -
-
-
-
4.1. The scheme for the hyperbolic system of equations The hyperbolic system has been written in the form of hyperbolic conservation laws,
Ow ot +
OF
~G +
(22)
= R
Then the first-order TVD scheme [6] may be written in the form •
__ Wi * .4.
.
-- ~7(W)i_l/2.k
-~-
•
G(W)i.k+l/2 -- G ( w ) i . , _ i / 2
where ~" is the time step, and h E, h,7 are spatial differt~nce steps. The numerical fluxes have the following tbrm:
=0,
(23)
82
O.V.Diyanlu~v et ai.IComputer Physics Communications 106 (1997) 76-94 ½. (~'(.,,.~) + F(.,,.f...,,))
.~,),"+,/~.~ =
--2' " ~ q=l
Q ( ~ ' q ) " l~..q i+l/2.k .
(
L (i+i/2.k .q
. [Wi+ I .k -- wi.nk] n
)}
a(.,);..,,+,/; ='. (G(_,?,.~)+c,,w,~+,)) tl
__!. 2
Q (~ v l ' q) • R vl'q • i,k+ i/2
(
L Tl'q • [wi.nk+ - w n i,k+ !/2 i i,k
•
(24)
.q and L ~,q n ~), corresponding Here R~/+1/2.~ i+ I/2.k are the right and left eigenveetors of the matrix (0F/dw)(wi+l/2. to the eigenvalue A¢'q. R~.k+t/2 ~'q •'q n and Li.k+l/2 are the right and left eigenvectors of the matrix (cgG/Jw)(wi.k+i/2), corresponding to the eigenvalue An'q. Q(A) is a function, which should satisfy the following conditions: 1" The
Q(A) > lal.
• Q(A) _< I,
i terms V~+I/2, k and V2k+i/2are
equal to
'q n Q( A¢'q) " R~/+t/2.k Li+,/2.kE'q " [W i+,.k -- W,~]
V~+,/2.k = ~
(
q=l
6 • R;.k+t/2.
Li.k+l/2
q=!
t,,+ w.
_will
)I )}
.
(25)
and are described in the appendix in detail. They have the meaning of artificial viscosity. The viscosity for contravariant local basis vectors T t and T 2 and the square root of the metric tensor determinant v/~ are used for coordinate viscosity determination, thus the difference equation for coordinates may be written in the form r "+Z = r " + r .
(u " + V,.).
The term Vc is also described in the appendix. This scheme is stable under the following Courant-Fridrichs-Levi condition: 1-a I
Ta 2
h---~-+ T ; < t, where
a' =IT'I"
, [,n,2 ,/f, i2 L~ +C~+v,~ +c0 -!T'I~
~.
i
,
,
and •
F-Ir=l.
m
I ! I [B[2 C2 ~//Is[2 ~.. L4.p+
+
~+c2
/2
C2
(B2)2]
- IT=I2" 4~rp "
This result is a consequence of the theorem proved in [ t0].
O.V. Diyankov et al./Computer Physics Communications 106 (1997) 76-94
83
At the second stage the new mesh is constructed. The algorithms for mesh construction mentioned above are used for this: (1) Lagrangian and Eulerian algorithms are trivial: in the Lagrange case the mesh is moving together with the flow; in the Eulerian case the mesh is motionless; (2) Local algorithm. In this case coordinates are determined using the coordinates of the current edge and 4 neighbours, r i+ llLm&' m ( ! 12,j+ 112 --
I - o') • "i+112,j+112 - ola _ add ÷ IT. ri+i/2,j+l/2 ,
radd i+112.j+i12---- " ~l " { ( d ; q - d ; ) .
[r°ld [ i-- !/20j+
. old ! 15 ÷ ui+312,j+i
]
12
+ (d~ + d~ )- [r old + . old 1 L i+t12,j-112 "i+l12,j+3/2J j ' d = ( d'~ + d ; + d ; + d ; ) , d f = Ir~+,/~.s+~/2 - r ~ - , / 2 j + , / 2 1 .
d(+ = [ r i + l / 2 j + l / 2 d~ Iri+ll2,S+ii2 =
-
ri+3/2j+l/2 [ , r~+l/2.S-llZl ,
d~ = lri+n/z,s+,/2- ri+,;2,s+./2l' • ~
(26)
,
where o" is a numerical parameter. It is used for decreasing the difference between coordinate and mass ve|oclties, thus making the mesh reconstruction process more stable. (3) Algebraic algorithm. In this c~e special algorithms are u~ed fer boundary grid rccons!r,scHr~n l,,~er edges coordinates are calculated according to the formula l
r n~v
[-
old
i+J/2.s+,/2 = - ~ • l f t/2.S+l/2" (N~ - i) +
r old
N~+I/2.S+~/2" i
)
! t' old .old .j) +-~n " ~ri+!/2'J/2" ( N" - j ) + *i+l/2,N,+i/2 ! (N¢-i)
N~
1 (N¢_i)
_old -_~I_ ( N~ - j ) . "1/~.!/2 1 • ..old
i
I ( N~ - j ) . r°ld N~+l/2,1/2
~ i~
I
I . _ola
(4) Algorithm, based on Poisson equation solving. This algorithm is described in details in paper [ I I ]. After the stage of mesh reconstruction, the coordinate velocities are calculated by new
Ui+ 112,j+ 112 ~" Ill+ ! 12,j+ i 12 -I-
ri+i/2,3+ 1/2
. old -- " i + 112,3+
112
7"
After that the numerical fluxes across the cell edge are calculated. For the cell with index (i, j) the flux across the edge (i + l / 2 , j ) is determined as follows: one should calculate 8Si+l/2.j for the plane geometry as t ! I the square of the rectangular [ri+l/2.j-l/2, ri+l/2.j_l/2, ri+l/24÷l/2, ri+l/2.j+i/2] (see Fig. 1) which equals to 8S~+~/2.s = [ ~ . ( u - v ) • T I • v ~ "
=1.
h~li+~/2j
{ ( X i + l / 2 , j - I ] 2 - Xi+t/2,j+l]2) ' , " (Yi+l/2,j-I/2-
Yi-t-l/2.j+l/2)
t --(Yi+i/2,j-I]2 - - . Vti+i/2,j+l[2 ) " ( Xi+l/2,j-I/2 -- Xi+l/2,j+i/2) }.
(27)
O.V. Diyankov et al.lCompuler Physics Communications 106 (1997) 76-94
84
r i+ l12,j+ ! / 2 ~
r t
Y i.
r t i+ l 1 2 , j , i1 2 !/2 r i.!/2.j-i/2 Fig. I. Fluxes between cells.
'l
" .......
j
"
. . . . . . . . . . . . . . . . . . . . . .
J
i
°f
.-- 1
e 1 10sL
0
2J_.../ r 0
I
2
3
i 4
5
6
7
8
9
10
======================
!, i . . . . . . . . . . -1 [ ................................................................................................... J 0 1 2 3 4 5 6 7 8 9 10
xlcml
X[~l Fig. 2. ID MHD test from paper [ 12]. The time is equal to 0.7.
If 8Si+tl2.j > 0, then the flux is directed out of the cell and the matter fraction 8 = 8 S / ~ h~h~ leaves this "-'-h ~ i+i h~" comes from cell cell and comes into cell (i + 1, j). In the other case the matter fraction 8 = 8S"/~/gij ( i + l , j ) into cell (i,j).
5. Tests A series of common gas dynamics ID test has been performed. The results were satisfying. Here the ID MHD test presented is taken fr~om paper [ 12] The results are shown in Fig. 2. They are very close to the results presented in the original paper. The initial conditions for the test are shown in Fig. 3.
6. Some examples of calculations 6. I. Laser beam interaction with thin aluminium plate
The modelling was perforrrled in cylindrical symmetry. Initially a plate of 5/~m thickness and 340/zm length with density of 2.7 g/cm 3 was heated by a laser beam, focused onto 50/zm. The energy with an amplitude of
O. V. Diyankov et al./Computer Physics Communications 106 (1997) 76-94
p=l
p=4
rigid wall T = 2 T=2 c)Hz '..~l,i '4 ~ 3---if'=0 0 Hz=l 14 Hz=2 10 t at
I ~
~Hz wall ~---ff'=
g•igid i
i
i
85 0
i
. . . .
X
Fig. 3. Geometry for ID MHD test calculation.
0.2 J with a triangle form for its time dependence (0.5 ns front duration and 0.5 ns decay time) was absorbed in the density range from 3 × 10--4 g/cm 3 up to 3 × 10-3 g/cm 3 (see Fig. 4, top). The plate was expanding, and in the middle of the expanded region the hot plasma has been produced. The problem has been modelled taking into account the spontaneous magnetic field (the term proportional to I / Z . [V ( I / p ) × ~'P] in Eq. (3) is responsible tbr this field). See the magnetic field distribution in Fig. 4, bottom. The maximum value of the magnetic field is up to 150 kGs. The magnetized electrical conductivity and heat conductivity coefficients were taken into account. From the temperature distributions, presented from Fig. 5, one can easily see that the maximum temperature is 10% larger when the spontaneous magnetic field is taken into account. The main reason for its increase is the process of magnetizing heat conductivity coefficients. In the calculation the numerical mesh was 120×80. The algorithm used for coordinate system moving was "local".
6.2. Plasma-wall interaction in a tokamak
Fig. 6 shows the initial geometry for this case. A particle beam is falling at the angle a,, = 4 ° on the wall. At the initial time particles of the beam transfer their energy to the wall. The heat wave begins to propagate inside the wail. As a result some portion of the wall is evaporated and expanded. Then the energy of the particles is transferred to the wall as well as to the expanding matter. Thus the matter obtains impulse of the particles too. The power flow of about I kJ/cm 2 is released during 100/zsec. The matter moves into the region where a magnetic field of the magnitude of H0 = 3 × 104 Gs is disposed. The magnetic field has also the direction with the angle of cz,, to the wall and influences the matter motion. The Euler variant of MAG code has been used. In the region ABCF (see Fig. 6), a background gas with the density of p = 10-t2g/cm 3 and temperature T = 2.5 × 10 - 7 keV was disposed. On the boundary GH the inflow condition specially developed for this problem was set. Results are presented in Figs. 7 and 8. The numerical mesh was 110×210.
6.3. Development of the m = 0 type instability in imploding plasma liners
This problem was previously modelled by Vikhrev et al. [ 131. Initially the plasma column was in Bennet equilibrium. In our case the initial data were: current 1 = ! MA, pl~sma density p,, = 0.2 g/cm 3, temperature 100 eV, initial radius r,, = 10 -2 cm. The perturbation region was taken in the center of the column and equals to 15% of the radius in the unperturbed region. Fig. 9 (top) shows the initial mesh. The other figures show pressure at the time of 0.6, 0.74, 0.9 ns. The evolution of plasma column is shown in Figs. 9 (bottom) and 10. These graphs illustrate the evolution of neck development up to the moment when cavity reaches the axis, reflects from axis, produces two shock waves propagating in opposite directions alone the axis. These shock waves meet at another place and a new hot spot appears. Then the process repeats again. In the calculation the numerical grid was roughly 30x40 and the algorithm for mesh reconstruction was "algebraic".
86
O. V. Diyankov et aL / Computer Physics Communications 106 (1997) 76-94
"RO" MAG t: 5.0016e-03 st: 885 N: 7 max: 1.779e+02
~I0e-2
4.03
3.28
2.52
1,76
1.01
0.25
-0.50
-1.26
-2.01 rain: 1.313o-05
-2.77 -3.40
-2.64
-1.89
-0.38
-1.13
0.38
1.13
1.89
2.64
x10e-2 3.40
"HZ" MAG t: 5.0016e-03 st: 885 N: 7 4.03 "1
max: 1.173O-02
xl0e-2
• 8.75e-03 8.00e-03 7,25e-03 6.50e.03 ' 5.75e-03 5.00e-03 4.25e-03 • 3.50e.03 2.75e-03 ' 2.00e-03 • 1.25e-03 5.00e-04 -2.50e-04 .I.00e.03 -1.73O-03
3.28
2.52
1.76
1.01
0.25
-0.50 -5,50e-03 -6.2~-93 -7.0~.03
-1.26
-2.01
i -9.25e-03 , .I .00e-02 ~
-2.77 -3.40
-2.64
-1.89
-1.13
mln: -1.161e-02 -0.38
0.38
1.13
1.89
2.64
x10e-2 3.40
Fig. 4. Density (top) and magnetic field (bottom) for laser plasma.
O. V Diyankov et aL/ Computer Physics Communications 106 (1997) 76-94 "T" M A G t: 5.0025e-03 st: 8 7 6 N: 7 4.05 -1
xl0e-2
max: 4.208e-01
' 3.75O-01
I
' 3.60e-01 ' 3.45O-01 "3.30e-01 • 3.15O-01
3.29
2.54
"3.00e-01
• 2.85O-01 2.70e.01 2.55O-01 2.40e-01 2.25o<)1 2.10e-01 t+95o-01
1.78
1.03
1.8oe.ol 1.65o.Ol
0.27
+1.50e-01 ,1.35O-01 ' 1.20e-01
-0.48
1,05o.01 9.ooe.o2 ,7.5oe.o2 ! 8,ooe-o2 ~4.5oe-o2 ; 3.00e-02 1.50e.02 • 1.00e-I0
-1.24
-1.99
-2.75
"-"'1 -3.40
-2.64
-1.89
-1.13
-0.38
0.38
1.13
1.89
2.64
rain: 1.000e-lO
xl0e-2 3.40
"T" M A G t: 5.0016e-03 st: 8 8 5 N: 7 max: 4.664e-01
4o3-1x l 0 e - 2
• 4.60e-01
3.28
4.40e.01 4,20e-01
2.52
' 3.80e.01
4.00e-01 3.60e-01
3.40e.01 3,20e-G1 ' 3.00e~1 ' 2.80e-01 +2.60e-01 ' 2.40e-01 • 2.20e-01
1.76
1.01
' 2.00e-01
0.25
,1.80e-01
1.60e-01 1.40e-01 1.20e-01 1.00e-01 8.00e-02 ' 6.00e-02 : 4+00e-02 • 2.00e-02
-0.50
!iZi+i i¸: -1.26
-2.01
'
"I
-2.77 " 1 - - - - - ' - 1 - - - - - - ' I - - " -3.40
-2.64
-139
-1.13
-0.38
0.38
1.!3
1.89
2.64
1.00e-10
rain: 1.000e-10
xl0e-2 3.40
Fig. 5. Tcmlmmturcs without (top) and with (bottom) magn:tic field.
87
O.V. Diyankov et aL / Computer Physics Communications 106 (1997) 76-.94
88
F
Partciel~beam
E
,D
1f----:-I-I
A
C
B Fig. 6. Initial geometry for plasma-wall interaction in a tokamak.
7. Conclusion
The full nonlinear MHD equations used for high power plasmas modelling are quite complex because of the presence of different spatial and temporal scales. Therefore, appropriate numerical methods and algorithms to model shock-like solutions are required for accurate simulations. Applications of MAG code to high-speed compressible flow simulations illustrate the effectiveness of the chosen numerical algorithms. The code can be applied to many other problems as well. In future, the MAG code will be developed in the following directions: - the average charge of ion will be calculated from the equations of dynamics of matter ionization and recombination; - the radiation transfer will be calculated in the close escape probabilities model from [ 14] ; - a two-temperature MHD model (electron temperature is not equal to ion temperature) will be developed.
Acknowledgements
The authors express their acknowledgment to Dr. Serge A, Terekhoff, whose h!itiative gave a strong impulse te the work on the MAG code development, and who has applied many efforts to move this work forward. The authors also want to thank Dr. V.A. Lykov for his friendly attention to this work and for multiple useful comments and discussions. We are glad to thank Igor Krasnogorov, Natali Savina and Rouslan Kotov for their help in preparing this work. This work is supported partially by ISTC, Projects # 107 and #009.
Appendix A. The flows for the TVD scheme
Here B ! GO ~
~
IB" T212 + IBa 12.
4~p
IT' 12
89
O.V. Diyankov et al./Computer Physics Communications 106 (1997) 76-94
RO g/cm~,MAGT'une: 70 mcsec xl0e2 2.00
max: 9.34%-06
7
1.00 0.00 .1.00 .2.00 .3.00 .4.00 .5.00
) -6.00 ~3 -7.00
XAUe~
7.00
6.00
8.00
9.00
10.00 11.00 12.00 13.00 14.00 15.00
RO g/cra~,MAG,Time: 70 mcsec
max: 3.524e-06
xlOe2 2.00
-1
1.00 0.00 -1.00 -2.00 -3.00 -4.00 -5.00 ) -6.00 13 -ZOO
,1%,A v ~ , f ~
6.00
7.00
8.00
9.00
10.00 ~.i.00 12.00 13.00 14.00 15.00
Fig. 7. E~n~ities without (top) and with (buttom) magnetic field.
90
O. V. Do,ankov et aL / Computer Physics Communications 106 (1997) 76-94
9
x 10 "=
/
~
%
2'0
4
'so
~'o
~b
.'o
.o
Fig. 8. Evaporation mass vs. angle of beam incidence.
I
Irll 2
+:)
-
47rpJ'
z; j'
where Ti are covariant local basis vectors, and BW,+l/2.k = Wi+t., - Wc, is the variation of value W. The artificial velo.cities for numerical fluxes in the ¢ direction are determined by Vu '
_---{ '~..~. C2" IT' 12 + I-~T~ " _f
a3 + a-~
ao . T}
Llaol
a2_L+ c. IT'I" laol
.
B2B ~
ITII 2 47rPv~la±l
+ iTil 2
v,,,=
t
~
T~. ao
+
'
T~! -
! } " A3 + As "~ul
{ B z • t~u2 - B2 . 8u z
a~ aoBz.
F.iT'I 2 + a . +2c . I T l l . l a o l A3 + A 5
'T'l
B2B ! 47rp-~-~
c. IT'l+laol TJ ao. l a - I } ,t3 + ,~5 + IT ~i----7",~3 + ,~s
( 1/g) Be" 8u2 + IT 112B: • 8Uz
~'lail
7'~
,I
4,rrPv,~
'T~ B2 B' ir~l~ 4~rpv~
I . 8u I a3 + a5
B2B I
c. IT'l+laol r~ ao.laill iao[ " IT~l'---~'4~rpv~lail " a' +': - [¢ I~" a' ¥ ~ ( l/g)B2.8u2 + IT1 ]2Bz • 8uz
4~-~. laii
+ a~
4~rpv ~
., '
O. V. Diyankov el aL / Computer Physics Communications 106 (1997) 76-94 ....
91
Method: MAG Time: O.O00000e+O0 ( Step 0 1.00e-02
6.00e.03
~~
-"
X,,
,,
,,q
2.00e-03
-2.00e-03
-
"
...
"
o,i
#
i
-6.00e-03
~'
i
" /
,,2 I
1.00e-05
O.O00e+O0
4,01e.03
8,0h-03
1.20e.02,
I
1.60e-02 2.00e-02 :~:~ O.O00e+O0
1.00e-02
6.00e-03
2.00e-03
-2.00e-03
-6.00e-03
-1.00e-02
I 1.60e-02
1.35%-0i
I
2.00e-02 5.495e-0i
Fig. 9. initial mesh for Z-pinch simulation (top) and pressure at the time t = 0.6 ns (bottom).
92
O.V. Diyankov et aL/Computer Physics Communications 106 (1997) 76-94
1.00e.02
2.00~, ' ~
-2.00e-03
.6.00e-03
-1.00e-02
,,
L2Oe-O: 1.224e.01
..... ~: Q!,ii,; ~,~ii~711:ii!:~!~! !!j~17i~3~!~I !?i,~
,,
1.60e-02
|
2.00e-02
i ¸
1.00e-02
6.01e-0~
2.01e,03
-2.00e.03
,6.00e-03
-1.00e-02 ,~.~,
1.300e-0|
....
......1 ~ ~i;:~! ~~8~ 1.60e-02 2.00e.02 l ....................................................................................................................... 1.897e+00 ......
Fig. I0. m = 0 type instability development. Pressure at the time t = 0.74 ns (top), time t = 0.9 ns (bottom)
O. V. Diyankov et aL / Computer Physics Communications 106 (1997) 76-94
Vu =
aoB2 {B:'Su2-B2"Su:)
"
-v/ga~
"
~
{ Bz -
93
ao }
~ ' A ' - - - -+' ' A ~
.13u'
}.(1/g)Bu.Su~,+I'I:"I2B:.Su,l a o l + c ' l r ' l
_(laol B~ laJ_l 4 ~
4v/T~. lall
"
B2" (Bz " 8B2 - B2" 8Bz) ( IT' 12 VB: = --ao" 4ztpga~ + Bz " i. A3 + As
a-~ + as
'
8P P
+ a~_ + laol" IT'i .c + laol 2 (l/g)n2.8B2 + IT'I 2 B: .sn~ ], A3 + A5
4~rpa~
V~=P {IT"" IT'I'c+ao. SPL p"
A3 + A 5
1
J '
(1/g)B:,.SB2+IT:]2Bz.6Bz}
pc + A 3 + A5
4~rp
Vv~=_V~.(iT, I. IT'I.c+ao. SP A3 + A 5
I pc + A3 + A 5 "
'
(I/g)B2.SB2+lril2Bz.SBz} 47rp
"
= 0, VT.I
VT? = 0,
B:-(B~.SB2-B2.SB~)
c. Ir'i+iaol a3 + a5
+
T~ IT'I
,
,
B I B2
8P
41rpv/g" laol
j.-~i 2 ,~ + A5 + signao. 4~rpv/~a ~
c. IT'l--t-laol h.3 + A5
T~
IT I I2 j,
( 1/g) B2 ¢~B2 +" IT I 12 Bz ••Bz 4~rp
VT: = --sign ao" f
c. ITll+laol T~
I
A3 + A 5
Vc.,=%, -'-- ~
+
lr'l
i~ii2 (l/g)g2"
cv
. ,7'2.
a3 + , :
[ -
4~'Pv/"~ a ~ . 4V/'4V/'~
I . 8P
B IB2 4 ~ p ~ . laol
B2 signa°'4~.Pv~a2 •
A3 + A5
"I'~
8B2 q-ITII 2 B z • 8Bz
4"trp
VT~
rE, = p ~ v~ +/,v~.
(,~-v,. +
1
,:.ov~, + ,:. v,:) + ~ ~ . (B~ . v,,.+ 8:.. v,,,+ 8: . v,,.) .
94
O. E Diyankov et ai./Computer Physics Communications 106 (1997) 76-94
References I 11 O.V Diyankov and S.A. Terekhoff, in: D-~as¢ Z-pinches, M. Haines and A. Knight, eds., AIP Conference Proc. 299 (ALP, New York, 1994) p. 121. 121 I.V. Giazyrin, O.V. Diyankov, R.A. Kotov and S.V. Koshelev, High School News 12 (Tomsk Univ., Physics, 1995) p. 23 l in Russian i. [31 S.L. Braginskii, Review of Plasma Physics, Vol. 1, M.A. Leontovich, ed. (Consultant Bureau, New York, 1965) p. 205. [41 V.V. Ermakov and N.N. Kalitkin, Tables of elect_~ical conductivity and electron heat conductivity coefficients for plasmas of ! I th matters (Appl. Math. Inst., Moscow, 1978), preprint lin Russiani. [51 O.V. Diyankov, TVD approach utilizing for Lagrangian gas dynamics conservation laws, RFNC-VNilTF preprint N 96 (1996) l in Russian I. 161 Ami Harten, High resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 49 (1983) 357-393. 171 O.V. Diyankov and S.A. Terekhoff, Difference scheme for thermal conductivity equation, RFNC-VNIITF preprint N 23 (1993) l in Russian !. [ 81 David S. Kershaw, Differencing of the ditthsion equation in Lagrangian hydrol)y~amic codes, J. Comput. Phys. 39 ( 1981 ) 375-395. [91 David S. Kershaw, The incomplete Cholesky-conjugate gradient method for the iterative solution of system of linear equations, .I. Comput. Phys. 26 (1978) 43-65. !1Ol O.V. Diyankov, Investigation of numerical scheme TVD-type stability for multidimensional linear hyperbolic system of equations, RFNC-VNIITF preprint N 97 (1993) [in Russian]. I I I I Stevhen A. Jordan ar,d Malcolm L. Spaulding, J. Comput. Phys. 104 (1993) ! 18-128. 1121 R.C. Sdvastava, K.G. Roesner and D. Leutbolf, Astrophys. & Space Sci. 135 (1987) 399-407. I13! V.V. Vikhrev, V.V. Ivanov wad G.A. Rozanova, Fizika Plasmy I Plasma Physicsi 15 (1989) 78 l in Russian I. 1141 J.P. Apruzese, J. Davis, D. Duston and K.G. Whitney, J. Quant. Spectrosc. Radiat. Transfer 23 (1980) 479.
Computer Physics Communications 106 (1997) 76-94
MAG - two-dimensional resistive MHD code using an arbitrary moving coordinate system Oleg V. Diyankov, Igor V. Glazyrin, Serge V. Koshelev Russian Federal Nuclear Center. All-Russian Research Institute of Technical Physics (RFNC - VNliTF), PO. Box 245. Snezhin.~k~ Chelyabinsk Region 456770, Russian Federation
Received 25 July 1996; revised 25 July 1997
Abstract The MAG code used for two-dimensional magnetohydrodynamics flows modelling is described. MAG algorithms are formulated for an arbitrary moving coordinate system. MAG can calculate flows with large deformations contained inside a region whose weakly deformed boundaries maintain the correct description of boundary conditions. Numerical experiments with two-dimensional MHD flows calculations show that the proposed numerical schemes are quite robust and accurate. Sample calculations illustrating the properties of the algorithms are presented. (~) 1997 Elsevier Science B.V. Keywo~s: 2 ½D MHD; Arbitrary moving coordinate system; TVD scheme
1. Introduction
The motivation for the creation of the MAG code was the necessity of modelling hot dense plasmas in a magnetic field. These plasmas are encountered in inertial fusion experiments such as laser-produced plasmas, tokamak plasmas, Z-pinch, etc. Typically, such plasmas have ," large ,,,~:,.,,,u~ D. . . . . la,. number and a large Lundquist number. Strong magnetic fields cause the transport of heat to be several orders smaller than without the fietd. These facts make the problem of code creation more difficult. There is a great variety of approaches in code constructing, which are used now for plasma flows in Megagauss magnetic fields simulation. Among them one can find Lagrangian, Eulerian, and Lagrange-Eulerian codes. Lagrangian codes have the advantage over the Eulerian ones in the simplicity of boundary condition formulation but ordinary Lagrangian codes do not permit to simulate long-term evolution of fluid instabilities because of numerical grid mixing. On the other hand, Eulerian schemes allow one to simulate flow motion with large deformations but there are difficulties with boundary conditions determination. To resolve boundary conditions one should introduce the high magnitude of electrical resistivity in the background non-ionized gas, disused in the region between the boundary and plasma: In this paper the code MAG [ 1,2] for plasma modelling in an arbitrary moving coordinate system at its present stage of development is presented. The use of an arbitrary moving coordinate system allows one to simulate flows with large deformations inside the flow region, conserving the correct description of conditions on its weakly deformated boundaries. This has proved to be highly successfid for strongly shocked flows. 0010-4655/97/517.001~) 1997 Elsevier Science B.V. All fights reserved. Pll SO0 ! 0-465 5 ( 97 )0008 3-0
O.V. Diyankov et al./Computer Physics Communications 106 (1997) 76-94
77
Arbitrary moving computational mesh permits the use of detailed resolution mesh in all main stages of plasma flow formation and reduces numerical diffusion of the magnetic field so that flows with a magnetic Reynolds number of the order of hundreds may be calculated. The code has been used for simulating three types of problems. The first one is the plate ablation under '~aser radiation. It is a well-known fact that during this process a spontaneous magnetic field may appear, and the estimation of its influence is an actual problem. The second problem is the simulation of tokamak plasma-wall interaction. The third one is modelling of Z-pinch implosion. 2. The model The system of equations used in the MAG code is delermined by the Braginskii [3] model for the onetemperature case (electron temperature is equal m ion temperature). The dynamic evolution of the system is described by a solution of the following equations: - continuity:
30 =.r_+ V. (pu) = 0 , at - momentum:
p.
+(u.V)u
(I)
)
+VP=-4~
[B×rotB],
(2)
- magnetic field diffusion: riB=rot O"t-
C2 [uxB]-~rotB
47ro"
c[(,)
+- V e
xVPe
~
)
]
-
c rot ( [rOtziniB x B] ) ~e
(3) '
- energy:
(a£ ) m p ( j £ ) + PVu + Vq p . - ~ - + (u. V)£ ANe (l+Zi) Zi _ _ -
-
j2 +
c
Per"
f
rot-B'~
+ p"
(Qext - Qrad)
(4)
heat flux: q = -~e V'T,
(5)
- Maxwell's equations: ~TB=0,
(6)
j = ~ rotB, (7) 4~ where t is time, u is the mass velocity, p is the mass density, ni =- p/(Aimp ) (with Ai for the atomic weight, and mt, for the mass of the proton), P is pressure, T,£ are the temperature and specific internal energy, B is the magnetic field, j is the current density, q is the heat conductivity flow, o', ~e are coefficients of electrical conductivity and heat condtJctivity, Qext is the energy released, Qrad is the radiation loss (bremsstrahlung
78
O, V. Diyankov et aL /Computer Physics Communications 106 (1997) 76-94
model at the present time), and Zi is the average charge of an ion (at present taken in the form Zi = Z ( p , T ) from paper [4]). Other variables are generally accepted. Assuming an independence of unknown variables (p, T, P, u, B) from the angle ~, in the case of cylindrical geomet~, and from z in the case of plane geometry, one obtains a system of equations, in which all unknown ~ a b l ~ de~nd on 2 spatial v~ables but the mass velocity vector and magnetic field vector have 3 components. This system is solved in MAC; code. Therefore, MAG code is said to be a 2½D MHD code.
2.1. Equation of state and transport coefficients 2.1.1. Equation of state for ideal gas Pressure and specific internal energy of ideal gas are obtained from the equations R
P f ( Zi + I ) "~i PT,
£ =
(Z~ +
1)
R
T
AiT-I
+
5Joniz,
Zi P~fP l+z~
(8)
where R is the universal gas constant, 7 = 5/3, and ~'ioniz is the energy for matter ionization.
2.1.2. Transport coeOfcients Transport coefficients are obtained in two ways: - magnetized transport coefficients for ideal gas from paper [ 3 ]; - widerange coefficients for average charge, heat and electrical conductivity from paper [4]. To account for magnetized effects all coefficients are multiplied by the correspondent factors taken from [ 3 ]. Heat conduction and magnetic diffusion terms are then in sca!~ f o x .
3. B o u n d a D " c o n d i t i o n s
A wide range of boundary conditions is implemented in the code. Three types of boundary conditions are considered.
3.1. Gas dynamics conditions They are: - applied pressure: PIb = f ( t ) , where PIb means the pressure at the corresponding boundary, and f ( t ) is a given time dependent function; - rigid wall: Unlb = 0, where un is a normal to boundary component of mass velocity; - piston: Unjb- f ( t ) ; - incoming flow (only for Eulerian variant of the code)" Ulb = fu(t); £1b = fc(t); Plb -- fp(t), where all variables at the correspondent boundaries are a given function of time.
3.2. Conditions for heat conductivity They are: given temperature, given heat flux, and heat flux as a function of temperature. These conditions may be written in the form
c~VnTjb+ I~l'jb =/(Tjb,
t) •
O.V. Diyankovet aL/Computer Physics Communications106 (1'997) 76-94
79
Here, TIb is the temperature at the corresponding boundary, a,/3 are numerical parameters ( a = 0 and ,8 = 1 for a given temperature and a = ! and/3 = 0 for the heat flux), f(T[b,t) is a given function• VnTIb is a normal to the boundary component of the temperature gradient.
3.3. Conditions for magnetic field These conditions may be written in the form: - symmetry: ~Bz/dn = O, OB~,/an = 0, B~ = 0; - conducting wall: Bz = 0, B~, = 0, Bn = 0; - axis: Bz = O, B~ = 0, B, = O, where B~, is used for the cylindrical symmetry and Bz for the plane one; - given current: Bz = 2¢rj/c - plane symmetry; B~, = 21/rc - cylindrical symmetry, where Bn is the normal to the boundary component of the magnetic field, j is the current density, 1 is the whole current in the z-direction, and r is the upper radius.
4.
Details
of
numerics
Below, the technique of MAG difference scheme derivation will be described• We will omit the detailed formulas of the difference scheme because of their complexity• Let us assign indices x and y in the case of the cylindrical geometry to the r- and z- components of vector variables and r- and z- components of independent spatial variables correspondingly, and index z to the ~ocomponent of vector variables,, Then we obtain a basic system of equations, which depends on the parameter of geometry u ( ~, = 0 for plane geometry, and u = 1 when the geometry is cylindrical). Because of Eq. (6), one can, substitute an equation t0r the z-component of the magnetic vector potemiai A instead of two equations for the x- and y- components of the magnetic field,
Bx = x -v maA 3y
By = - x - v --.3A Ox
(9)
Eqs. ( 1 ) - ( 5 ) have been written in an arbitrary moving coordinate system, and then splitting into two systems has been performed. The diffusive terms have been split into a separate system of equations, and the remained terms produce a quasilinear hyperbolic system. The first system is a hyperbolic one (compare with ( I ) - ( 4 ) ): - Equation of continuity: Opx"v at
lOOx"v/g(uk--vk) k
=0.
(iO)
- Euler equations for different components of velocity: • x component: ap x
u., +
Ot • y
cgPX~'v~(Uk -- Vk)Ux O~~
+
~"aPf ~xvrg O~k
=
i { 3x'v/~B'Bx 4¢r O~k
_
u B 2 1 + ~,pu2 ," ~ "
(11)
component: ap X~'V~ ltv- +
Ot
~pXUv~(Uk--Vk)U, ,
x~OPIT~yv~
1{OxVv~BkB~' 1
(12)
O. V. Diyankov et aL/ Computer Physics Communications 106 (1997) 76-94
80
• z component: Op x 2~'vfg u z
+ af'x2% (uk vk)u
i f Px2"v/'~BkBz 1
-
(13)
at -
Equation for A:
3A + Ot -
8A
( u k - v k) ~
=0.
(14)
Equation for z component of magnetic field:
av~Bz + av'gBz(u k - v k ) ) = aV~U~B~ at a~k o~k "
(15)
- Equation for total energy:
a£/ ' x~v/gat + a~-fi{£Ix~v~(uk - vk) + PYx~j;°uk} = PX u x/rg (Qext - - Qrad) 4- ~ -
1 {aX~v/gB'(Bu)} 3~e
.
(16)
Equation for the square root of the metric tensor:
at
(17)
a~k
The sum is taken over the same upper and lower indices. Here u k is a projection of the mass velocity vector to the kth vector of a local basis, u I' = u. T t, where ~ is a covariant local basis vector, v k, B k are determined in the same way, v is the coordinate system velocity, g is a determinant of metric tensor with covariant component go, ~k are mathematical coordinates, the total energy Ef = p . (£ + ½. [ul 2) + ( ! / 8 7 r ) . [B[2, and the total pressure P / = P(p,£) + (1/87r). [B[2. The second system of equations is a diffusive one, p
:
equation of continuity, - equation for determinant of metric tensor, - Euler equation for velocity.
coast.
-
ffi c o a s t . u
-
-
coast.
(18)
Equation for A:
aAvrg = x ~ 41rtr "-'7 a~ ~kv~x-P " ~ at
(19)
"
- Equation for Bz:
0 { c2 Ot
:Ox~Bz)}
cmrAi e Z
{ 0 (I)3P "~ a~J
- Equation for energy: O~'xJ'Pv~
! XVv~ fgiJx-2V axVBz ax~'Bz
(
a
oP
(20)
O. V. Diyankov et al. / Computer Physics Communicatimas 106 (1997) 76-94
I (_~ ..
d-- "
gtj x - p V ~
g
O~j)2}
trip c9 ( "t- n . Ai .
e
cmraix/"g(OBzx v OP "~41re Zi p d~i O~J
~
8l
Jj.~_,
)
( l + Z i ) . Z i "vc~x"
OBzx" OP) tgscj d~i •
(21)
The flow region is a geometric figure, topologically equivalent to a rectangle, and it is initially covered by regular rectangular mesh. All dependent variables are specified at the center of the numerical mesh cell, besides coordinates which are vertex-centred. The following approach is used for the hyperbolic system solving. The process of the solution of the equation on each time step is broken into two stages. At the first stage the coordinate system velocity is considered equal to the mass velocity (and this stage is called a Lagrangian stage), and the corresponding hyperbolic system is solved, using difference scheme, which was derived from the TVD approach [5], first suggested by Harten in paper [ 6]. At the second stage the grid moving is performed and fluxes between cells are calculated. This approach allows one to simulate flows with rather large deformations inside a flow region on regular grids. The second system (of diffusive type) is solved in the following manner. The equations of magnetic diffusion are written for the z-component of magnetic field and z-component of magnetic vector potential, and the resulting system consists of three nontrivial diffusive equations: equation of heat conductivity; equation of z-component of magnetic field diffusion; equation of z-component of magnetic vector potential diffusion. The main terms of these equations are moved to the left sides and a difference scheme [7], analogous to -
-
-
t k a lkr--rel~olu e~homo [.~1 ~...,~,_
I.~_T~
O-=-='~--.-
~ = .
. . . . .
.
v
~.,
ie , t ~ r l ft~r ~ l u ; n ~ .
.
.
.
. .
.
.
.
.
,,,.v..
_-~,,,~,
.=.~k ,~f' tll~ae~ onn,~ti~ne ~
._
"--""
'~''--"-"--"~
•
~ ; m n l o i t p r n | i ~ n ie neoti I n ro~:nlv~ tho_
",---'.'"~'-P-"
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
nonlinearity of the right sides and the ICCG method [ 9] is used to solve linear algebraic systems. Five algorithms are used for coordinate system moving. They are Lagrangian, when the coordinate system velocity equals the mass velocity; Eulerian, when the coordinate system velocity equals to zero; local, when the coordinate system velocity at a mesh vertex is determined using mass velocities in the 8 neighbour sites; - algebraic, when the coordinate system velocity at a mesh vertex is determined using mass velocities on the boundaries of the calculated area; Poisson-type, when the coordinate system velocity field is determined by solving a Poisson equation, the boundary conditions for which are determined by the mass velocities on the calculated area boundaries. -
-
-
-
4.1. The scheme for the hyperbolic system of equations The hyperbolic system has been written in the form of hyperbolic conservation laws,
Ow ot +
OF
~G +
(22)
= R
Then the first-order TVD scheme [6] may be written in the form •
__ Wi * .4.
.
-- ~7(W)i_l/2.k
-~-
•
G(W)i.k+l/2 -- G ( w ) i . , _ i / 2
where ~" is the time step, and h E, h,7 are spatial differt~nce steps. The numerical fluxes have the following tbrm:
=0,
(23)
82
O.V.Diyanlu~v et ai.IComputer Physics Communications 106 (1997) 76-94 ½. (~'(.,,.~) + F(.,,.f...,,))
.~,),"+,/~.~ =
--2' " ~ q=l
Q ( ~ ' q ) " l~..q i+l/2.k .
(
L (i+i/2.k .q
. [Wi+ I .k -- wi.nk] n
)}
a(.,);..,,+,/; ='. (G(_,?,.~)+c,,w,~+,)) tl
__!. 2
Q (~ v l ' q) • R vl'q • i,k+ i/2
(
L Tl'q • [wi.nk+ - w n i,k+ !/2 i i,k
•
(24)
.q and L ~,q n ~), corresponding Here R~/+1/2.~ i+ I/2.k are the right and left eigenveetors of the matrix (0F/dw)(wi+l/2. to the eigenvalue A¢'q. R~.k+t/2 ~'q •'q n and Li.k+l/2 are the right and left eigenvectors of the matrix (cgG/Jw)(wi.k+i/2), corresponding to the eigenvalue An'q. Q(A) is a function, which should satisfy the following conditions: 1" The
Q(A) > lal.
• Q(A) _< I,
i terms V~+I/2, k and V2k+i/2are
equal to
'q n Q( A¢'q) " R~/+t/2.k Li+,/2.kE'q " [W i+,.k -- W,~]
V~+,/2.k = ~
(
q=l
6 • R;.k+t/2.
Li.k+l/2
q=!
t,,+ w.
_will
)I )}
.
(25)
and are described in the appendix in detail. They have the meaning of artificial viscosity. The viscosity for contravariant local basis vectors T t and T 2 and the square root of the metric tensor determinant v/~ are used for coordinate viscosity determination, thus the difference equation for coordinates may be written in the form r "+Z = r " + r .
(u " + V,.).
The term Vc is also described in the appendix. This scheme is stable under the following Courant-Fridrichs-Levi condition: 1-a I
Ta 2
h---~-+ T ; < t, where
a' =IT'I"
, [,n,2 ,/f, i2 L~ +C~+v,~ +c0 -!T'I~
~.
i
,
,
and •
F-Ir=l.
m
I ! I [B[2 C2 ~//Is[2 ~.. L4.p+
+
~+c2
/2
C2
(B2)2]
- IT=I2" 4~rp "
This result is a consequence of the theorem proved in [ t0].
O.V. Diyankov et al./Computer Physics Communications 106 (1997) 76-94
83
At the second stage the new mesh is constructed. The algorithms for mesh construction mentioned above are used for this: (1) Lagrangian and Eulerian algorithms are trivial: in the Lagrange case the mesh is moving together with the flow; in the Eulerian case the mesh is motionless; (2) Local algorithm. In this case coordinates are determined using the coordinates of the current edge and 4 neighbours, r i+ llLm&' m ( ! 12,j+ 112 --
I - o') • "i+112,j+112 - ola _ add ÷ IT. ri+i/2,j+l/2 ,
radd i+112.j+i12---- " ~l " { ( d ; q - d ; ) .
[r°ld [ i-- !/20j+
. old ! 15 ÷ ui+312,j+i
]
12
+ (d~ + d~ )- [r old + . old 1 L i+t12,j-112 "i+l12,j+3/2J j ' d = ( d'~ + d ; + d ; + d ; ) , d f = Ir~+,/~.s+~/2 - r ~ - , / 2 j + , / 2 1 .
d(+ = [ r i + l / 2 j + l / 2 d~ Iri+ll2,S+ii2 =
-
ri+3/2j+l/2 [ , r~+l/2.S-llZl ,
d~ = lri+n/z,s+,/2- ri+,;2,s+./2l' • ~
(26)
,
where o" is a numerical parameter. It is used for decreasing the difference between coordinate and mass ve|oclties, thus making the mesh reconstruction process more stable. (3) Algebraic algorithm. In this c~e special algorithms are u~ed fer boundary grid rccons!r,scHr~n l,,~er edges coordinates are calculated according to the formula l
r n~v
[-
old
i+J/2.s+,/2 = - ~ • l f t/2.S+l/2" (N~ - i) +
r old
N~+I/2.S+~/2" i
)
! t' old .old .j) +-~n " ~ri+!/2'J/2" ( N" - j ) + *i+l/2,N,+i/2 ! (N¢-i)
N~
1 (N¢_i)
_old -_~I_ ( N~ - j ) . "1/~.!/2 1 • ..old
i
I ( N~ - j ) . r°ld N~+l/2,1/2
~ i~
I
I . _ola
(4) Algorithm, based on Poisson equation solving. This algorithm is described in details in paper [ I I ]. After the stage of mesh reconstruction, the coordinate velocities are calculated by new
Ui+ 112,j+ 112 ~" Ill+ ! 12,j+ i 12 -I-
ri+i/2,3+ 1/2
. old -- " i + 112,3+
112
7"
After that the numerical fluxes across the cell edge are calculated. For the cell with index (i, j) the flux across the edge (i + l / 2 , j ) is determined as follows: one should calculate 8Si+l/2.j for the plane geometry as t ! I the square of the rectangular [ri+l/2.j-l/2, ri+l/2.j_l/2, ri+l/24÷l/2, ri+l/2.j+i/2] (see Fig. 1) which equals to 8S~+~/2.s = [ ~ . ( u - v ) • T I • v ~ "
=1.
h~li+~/2j
{ ( X i + l / 2 , j - I ] 2 - Xi+t/2,j+l]2) ' , " (Yi+l/2,j-I/2-
Yi-t-l/2.j+l/2)
t --(Yi+i/2,j-I]2 - - . Vti+i/2,j+l[2 ) " ( Xi+l/2,j-I/2 -- Xi+l/2,j+i/2) }.
(27)
O.V. Diyankov et al.lCompuler Physics Communications 106 (1997) 76-94
84
r i+ l12,j+ ! / 2 ~
r t
Y i.
r t i+ l 1 2 , j , i1 2 !/2 r i.!/2.j-i/2 Fig. I. Fluxes between cells.
'l
" .......
j
"
. . . . . . . . . . . . . . . . . . . . . .
J
i
°f
.-- 1
e 1 10sL
0
2J_.../ r 0
I
2
3
i 4
5
6
7
8
9
10
======================
!, i . . . . . . . . . . -1 [ ................................................................................................... J 0 1 2 3 4 5 6 7 8 9 10
xlcml
X[~l Fig. 2. ID MHD test from paper [ 12]. The time is equal to 0.7.
If 8Si+tl2.j > 0, then the flux is directed out of the cell and the matter fraction 8 = 8 S / ~ h~h~ leaves this "-'-h ~ i+i h~" comes from cell cell and comes into cell (i + 1, j). In the other case the matter fraction 8 = 8S"/~/gij ( i + l , j ) into cell (i,j).
5. Tests A series of common gas dynamics ID test has been performed. The results were satisfying. Here the ID MHD test presented is taken fr~om paper [ 12] The results are shown in Fig. 2. They are very close to the results presented in the original paper. The initial conditions for the test are shown in Fig. 3.
6. Some examples of calculations 6. I. Laser beam interaction with thin aluminium plate
The modelling was perforrrled in cylindrical symmetry. Initially a plate of 5/~m thickness and 340/zm length with density of 2.7 g/cm 3 was heated by a laser beam, focused onto 50/zm. The energy with an amplitude of
O. V. Diyankov et al./Computer Physics Communications 106 (1997) 76-94
p=l
p=4
rigid wall T = 2 T=2 c)Hz '..~l,i '4 ~ 3---if'=0 0 Hz=l 14 Hz=2 10 t at
I ~
~Hz wall ~---ff'=
g•igid i
i
i
85 0
i
. . . .
X
Fig. 3. Geometry for ID MHD test calculation.
0.2 J with a triangle form for its time dependence (0.5 ns front duration and 0.5 ns decay time) was absorbed in the density range from 3 × 10--4 g/cm 3 up to 3 × 10-3 g/cm 3 (see Fig. 4, top). The plate was expanding, and in the middle of the expanded region the hot plasma has been produced. The problem has been modelled taking into account the spontaneous magnetic field (the term proportional to I / Z . [V ( I / p ) × ~'P] in Eq. (3) is responsible tbr this field). See the magnetic field distribution in Fig. 4, bottom. The maximum value of the magnetic field is up to 150 kGs. The magnetized electrical conductivity and heat conductivity coefficients were taken into account. From the temperature distributions, presented from Fig. 5, one can easily see that the maximum temperature is 10% larger when the spontaneous magnetic field is taken into account. The main reason for its increase is the process of magnetizing heat conductivity coefficients. In the calculation the numerical mesh was 120×80. The algorithm used for coordinate system moving was "local".
6.2. Plasma-wall interaction in a tokamak
Fig. 6 shows the initial geometry for this case. A particle beam is falling at the angle a,, = 4 ° on the wall. At the initial time particles of the beam transfer their energy to the wall. The heat wave begins to propagate inside the wail. As a result some portion of the wall is evaporated and expanded. Then the energy of the particles is transferred to the wall as well as to the expanding matter. Thus the matter obtains impulse of the particles too. The power flow of about I kJ/cm 2 is released during 100/zsec. The matter moves into the region where a magnetic field of the magnitude of H0 = 3 × 104 Gs is disposed. The magnetic field has also the direction with the angle of cz,, to the wall and influences the matter motion. The Euler variant of MAG code has been used. In the region ABCF (see Fig. 6), a background gas with the density of p = 10-t2g/cm 3 and temperature T = 2.5 × 10 - 7 keV was disposed. On the boundary GH the inflow condition specially developed for this problem was set. Results are presented in Figs. 7 and 8. The numerical mesh was 110×210.
6.3. Development of the m = 0 type instability in imploding plasma liners
This problem was previously modelled by Vikhrev et al. [ 131. Initially the plasma column was in Bennet equilibrium. In our case the initial data were: current 1 = ! MA, pl~sma density p,, = 0.2 g/cm 3, temperature 100 eV, initial radius r,, = 10 -2 cm. The perturbation region was taken in the center of the column and equals to 15% of the radius in the unperturbed region. Fig. 9 (top) shows the initial mesh. The other figures show pressure at the time of 0.6, 0.74, 0.9 ns. The evolution of plasma column is shown in Figs. 9 (bottom) and 10. These graphs illustrate the evolution of neck development up to the moment when cavity reaches the axis, reflects from axis, produces two shock waves propagating in opposite directions alone the axis. These shock waves meet at another place and a new hot spot appears. Then the process repeats again. In the calculation the numerical grid was roughly 30x40 and the algorithm for mesh reconstruction was "algebraic".
86
O. V. Diyankov et aL / Computer Physics Communications 106 (1997) 76-94
"RO" MAG t: 5.0016e-03 st: 885 N: 7 max: 1.779e+02
~I0e-2
4.03
3.28
2.52
1,76
1.01
0.25
-0.50
-1.26
-2.01 rain: 1.313o-05
-2.77 -3.40
-2.64
-1.89
-0.38
-1.13
0.38
1.13
1.89
2.64
x10e-2 3.40
"HZ" MAG t: 5.0016e-03 st: 885 N: 7 4.03 "1
max: 1.173O-02
xl0e-2
• 8.75e-03 8.00e-03 7,25e-03 6.50e.03 ' 5.75e-03 5.00e-03 4.25e-03 • 3.50e.03 2.75e-03 ' 2.00e-03 • 1.25e-03 5.00e-04 -2.50e-04 .I.00e.03 -1.73O-03
3.28
2.52
1.76
1.01
0.25
-0.50 -5,50e-03 -6.2~-93 -7.0~.03
-1.26
-2.01
i -9.25e-03 , .I .00e-02 ~
-2.77 -3.40
-2.64
-1.89
-1.13
mln: -1.161e-02 -0.38
0.38
1.13
1.89
2.64
x10e-2 3.40
Fig. 4. Density (top) and magnetic field (bottom) for laser plasma.
O. V Diyankov et aL/ Computer Physics Communications 106 (1997) 76-94 "T" M A G t: 5.0025e-03 st: 8 7 6 N: 7 4.05 -1
xl0e-2
max: 4.208e-01
' 3.75O-01
I
' 3.60e-01 ' 3.45O-01 "3.30e-01 • 3.15O-01
3.29
2.54
"3.00e-01
• 2.85O-01 2.70e.01 2.55O-01 2.40e-01 2.25o<)1 2.10e-01 t+95o-01
1.78
1.03
1.8oe.ol 1.65o.Ol
0.27
+1.50e-01 ,1.35O-01 ' 1.20e-01
-0.48
1,05o.01 9.ooe.o2 ,7.5oe.o2 ! 8,ooe-o2 ~4.5oe-o2 ; 3.00e-02 1.50e.02 • 1.00e-I0
-1.24
-1.99
-2.75
"-"'1 -3.40
-2.64
-1.89
-1.13
-0.38
0.38
1.13
1.89
2.64
rain: 1.000e-lO
xl0e-2 3.40
"T" M A G t: 5.0016e-03 st: 8 8 5 N: 7 max: 4.664e-01
4o3-1x l 0 e - 2
• 4.60e-01
3.28
4.40e.01 4,20e-01
2.52
' 3.80e.01
4.00e-01 3.60e-01
3.40e.01 3,20e-G1 ' 3.00e~1 ' 2.80e-01 +2.60e-01 ' 2.40e-01 • 2.20e-01
1.76
1.01
' 2.00e-01
0.25
,1.80e-01
1.60e-01 1.40e-01 1.20e-01 1.00e-01 8.00e-02 ' 6.00e-02 : 4+00e-02 • 2.00e-02
-0.50
!iZi+i i¸: -1.26
-2.01
'
"I
-2.77 " 1 - - - - - ' - 1 - - - - - - ' I - - " -3.40
-2.64
-139
-1.13
-0.38
0.38
1.!3
1.89
2.64
1.00e-10
rain: 1.000e-10
xl0e-2 3.40
Fig. 5. Tcmlmmturcs without (top) and with (bottom) magn:tic field.
87
O.V. Diyankov et aL / Computer Physics Communications 106 (1997) 76-.94
88
F
Partciel~beam
E
,D
1f----:-I-I
A
C
B Fig. 6. Initial geometry for plasma-wall interaction in a tokamak.
7. Conclusion
The full nonlinear MHD equations used for high power plasmas modelling are quite complex because of the presence of different spatial and temporal scales. Therefore, appropriate numerical methods and algorithms to model shock-like solutions are required for accurate simulations. Applications of MAG code to high-speed compressible flow simulations illustrate the effectiveness of the chosen numerical algorithms. The code can be applied to many other problems as well. In future, the MAG code will be developed in the following directions: - the average charge of ion will be calculated from the equations of dynamics of matter ionization and recombination; - the radiation transfer will be calculated in the close escape probabilities model from [ 14] ; - a two-temperature MHD model (electron temperature is not equal to ion temperature) will be developed.
Acknowledgements
The authors express their acknowledgment to Dr. Serge A, Terekhoff, whose h!itiative gave a strong impulse te the work on the MAG code development, and who has applied many efforts to move this work forward. The authors also want to thank Dr. V.A. Lykov for his friendly attention to this work and for multiple useful comments and discussions. We are glad to thank Igor Krasnogorov, Natali Savina and Rouslan Kotov for their help in preparing this work. This work is supported partially by ISTC, Projects # 107 and #009.
Appendix A. The flows for the TVD scheme
Here B ! GO ~
~
IB" T212 + IBa 12.
4~p
IT' 12
89
O.V. Diyankov et al./Computer Physics Communications 106 (1997) 76-94
RO g/cm~,MAGT'une: 70 mcsec xl0e2 2.00
max: 9.34%-06
7
1.00 0.00 .1.00 .2.00 .3.00 .4.00 .5.00
) -6.00 ~3 -7.00
XAUe~
7.00
6.00
8.00
9.00
10.00 11.00 12.00 13.00 14.00 15.00
RO g/cra~,MAG,Time: 70 mcsec
max: 3.524e-06
xlOe2 2.00
-1
1.00 0.00 -1.00 -2.00 -3.00 -4.00 -5.00 ) -6.00 13 -ZOO
,1%,A v ~ , f ~
6.00
7.00
8.00
9.00
10.00 ~.i.00 12.00 13.00 14.00 15.00
Fig. 7. E~n~ities without (top) and with (buttom) magnetic field.
90
O. V. Do,ankov et aL / Computer Physics Communications 106 (1997) 76-94
9
x 10 "=
/
~
%
2'0
4
'so
~'o
~b
.'o
.o
Fig. 8. Evaporation mass vs. angle of beam incidence.
I
Irll 2
+:)
-
47rpJ'
z; j'
where Ti are covariant local basis vectors, and BW,+l/2.k = Wi+t., - Wc, is the variation of value W. The artificial velo.cities for numerical fluxes in the ¢ direction are determined by Vu '
_---{ '~..~. C2" IT' 12 + I-~T~ " _f
a3 + a-~
ao . T}
Llaol
a2_L+ c. IT'I" laol
.
B2B ~
ITII 2 47rPv~la±l
+ iTil 2
v,,,=
t
~
T~. ao
+
'
T~! -
! } " A3 + As "~ul
{ B z • t~u2 - B2 . 8u z
a~ aoBz.
F.iT'I 2 + a . +2c . I T l l . l a o l A3 + A 5
'T'l
B2B ! 47rp-~-~
c. IT'l+laol TJ ao. l a - I } ,t3 + ,~5 + IT ~i----7",~3 + ,~s
( 1/g) Be" 8u2 + IT 112B: • 8Uz
~'lail
7'~
,I
4,rrPv,~
'T~ B2 B' ir~l~ 4~rpv~
I . 8u I a3 + a5
B2B I
c. IT'l+laol r~ ao.laill iao[ " IT~l'---~'4~rpv~lail " a' +': - [¢ I~" a' ¥ ~ ( l/g)B2.8u2 + IT1 ]2Bz • 8uz
4~-~. laii
+ a~
4~rpv ~
., '
O. V. Diyankov el aL / Computer Physics Communications 106 (1997) 76-94 ....
91
Method: MAG Time: O.O00000e+O0 ( Step 0 1.00e-02
6.00e.03
~~
-"
X,,
,,
,,q
2.00e-03
-2.00e-03
-
"
...
"
o,i
#
i
-6.00e-03
~'
i
" /
,,2 I
1.00e-05
O.O00e+O0
4,01e.03
8,0h-03
1.20e.02,
I
1.60e-02 2.00e-02 :~:~ O.O00e+O0
1.00e-02
6.00e-03
2.00e-03
-2.00e-03
-6.00e-03
-1.00e-02
I 1.60e-02
1.35%-0i
I
2.00e-02 5.495e-0i
Fig. 9. initial mesh for Z-pinch simulation (top) and pressure at the time t = 0.6 ns (bottom).
92
O.V. Diyankov et aL/Computer Physics Communications 106 (1997) 76-94
1.00e.02
2.00~, ' ~
-2.00e-03
.6.00e-03
-1.00e-02
,,
L2Oe-O: 1.224e.01
..... ~: Q!,ii,; ~,~ii~711:ii!:~!~! !!j~17i~3~!~I !?i,~
,,
1.60e-02
|
2.00e-02
i ¸
1.00e-02
6.01e-0~
2.01e,03
-2.00e.03
,6.00e-03
-1.00e-02 ,~.~,
1.300e-0|
....
......1 ~ ~i;:~! ~~8~ 1.60e-02 2.00e.02 l ....................................................................................................................... 1.897e+00 ......
Fig. I0. m = 0 type instability development. Pressure at the time t = 0.74 ns (top), time t = 0.9 ns (bottom)
O. V. Diyankov et aL / Computer Physics Communications 106 (1997) 76-94
Vu =
aoB2 {B:'Su2-B2"Su:)
"
-v/ga~
"
~
{ Bz -
93
ao }
~ ' A ' - - - -+' ' A ~
.13u'
}.(1/g)Bu.Su~,+I'I:"I2B:.Su,l a o l + c ' l r ' l
_(laol B~ laJ_l 4 ~
4v/T~. lall
"
B2" (Bz " 8B2 - B2" 8Bz) ( IT' 12 VB: = --ao" 4ztpga~ + Bz " i. A3 + As
a-~ + as
'
8P P
+ a~_ + laol" IT'i .c + laol 2 (l/g)n2.8B2 + IT'I 2 B: .sn~ ], A3 + A5
4~rpa~
V~=P {IT"" IT'I'c+ao. SPL p"
A3 + A 5
1
J '
(1/g)B:,.SB2+IT:]2Bz.6Bz}
pc + A 3 + A5
4~rp
Vv~=_V~.(iT, I. IT'I.c+ao. SP A3 + A 5
I pc + A3 + A 5 "
'
(I/g)B2.SB2+lril2Bz.SBz} 47rp
"
= 0, VT.I
VT? = 0,
B:-(B~.SB2-B2.SB~)
c. Ir'i+iaol a3 + a5
+
T~ IT'I
,
,
B I B2
8P
41rpv/g" laol
j.-~i 2 ,~ + A5 + signao. 4~rpv/~a ~
c. IT'l--t-laol h.3 + A5
T~
IT I I2 j,
( 1/g) B2 ¢~B2 +" IT I 12 Bz ••Bz 4~rp
VT: = --sign ao" f
c. ITll+laol T~
I
A3 + A 5
Vc.,=%, -'-- ~
+
lr'l
i~ii2 (l/g)g2"
cv
. ,7'2.
a3 + , :
[ -
4~'Pv/"~ a ~ . 4V/'4V/'~
I . 8P
B IB2 4 ~ p ~ . laol
B2 signa°'4~.Pv~a2 •
A3 + A5
"I'~
8B2 q-ITII 2 B z • 8Bz
4"trp
VT~
rE, = p ~ v~ +/,v~.
(,~-v,. +
1
,:.ov~, + ,:. v,:) + ~ ~ . (B~ . v,,.+ 8:.. v,,,+ 8: . v,,.) .
94
O. E Diyankov et ai./Computer Physics Communications 106 (1997) 76-94
References I 11 O.V Diyankov and S.A. Terekhoff, in: D-~as¢ Z-pinches, M. Haines and A. Knight, eds., AIP Conference Proc. 299 (ALP, New York, 1994) p. 121. 121 I.V. Giazyrin, O.V. Diyankov, R.A. Kotov and S.V. Koshelev, High School News 12 (Tomsk Univ., Physics, 1995) p. 23 l in Russian i. [31 S.L. Braginskii, Review of Plasma Physics, Vol. 1, M.A. Leontovich, ed. (Consultant Bureau, New York, 1965) p. 205. [41 V.V. Ermakov and N.N. Kalitkin, Tables of elect_~ical conductivity and electron heat conductivity coefficients for plasmas of ! I th matters (Appl. Math. Inst., Moscow, 1978), preprint lin Russiani. [51 O.V. Diyankov, TVD approach utilizing for Lagrangian gas dynamics conservation laws, RFNC-VNilTF preprint N 96 (1996) l in Russian I. 161 Ami Harten, High resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 49 (1983) 357-393. 171 O.V. Diyankov and S.A. Terekhoff, Difference scheme for thermal conductivity equation, RFNC-VNIITF preprint N 23 (1993) l in Russian !. [ 81 David S. Kershaw, Differencing of the ditthsion equation in Lagrangian hydrol)y~amic codes, J. Comput. Phys. 39 ( 1981 ) 375-395. [91 David S. Kershaw, The incomplete Cholesky-conjugate gradient method for the iterative solution of system of linear equations, .I. Comput. Phys. 26 (1978) 43-65. !1Ol O.V. Diyankov, Investigation of numerical scheme TVD-type stability for multidimensional linear hyperbolic system of equations, RFNC-VNIITF preprint N 97 (1993) [in Russian]. I I I I Stevhen A. Jordan ar,d Malcolm L. Spaulding, J. Comput. Phys. 104 (1993) ! 18-128. 1121 R.C. Sdvastava, K.G. Roesner and D. Leutbolf, Astrophys. & Space Sci. 135 (1987) 399-407. I13! V.V. Vikhrev, V.V. Ivanov wad G.A. Rozanova, Fizika Plasmy I Plasma Physicsi 15 (1989) 78 l in Russian I. 1141 J.P. Apruzese, J. Davis, D. Duston and K.G. Whitney, J. Quant. Spectrosc. Radiat. Transfer 23 (1980) 479.