3. Diffusion of magnetic field in an immovable plane layer of plasma with regard to joule heating and its effect on the diffusion and heat conduction coefficients

The problem of magnetic diffusion in a plane layer of material has many applications in practice [14]. In its detailed formulation, the problem was considered in paper [14] for mega gauss fields. Hydrodynamic motion, magnetic diffusion, heat conduction by electrons, and radiant heat exchange in the "back and forth" approximation were taken into account. Since finding an exact solution to such a problem causes difficulties, the original formulation needs to be simplified. Self-similar solutions to the problem obtained with simplifying assumptions were also presented in [15].

A model problem is considered with the following assumptions:


With such assumptions, the problem is reduced to solving equations

$$\frac{d\mathbf{H}}{dt} = -rot(\mathbf{v} \cdot rot\mathbf{H}), \ \rho \frac{dT}{dt} = (\gamma - 1)\nu(rot\mathbf{H} \cdot rot\mathbf{H}),\tag{4}$$

where <sup>ν</sup> <sup>¼</sup> <sup>c</sup><sup>2</sup>=4πσ, <sup>σ</sup> <sup>¼</sup> <sup>σ</sup>0ð Þ <sup>T</sup>=Ry <sup>α</sup> , α ¼ 3=2, γ � 1 ¼ R=CV, R ¼ 1, CV ¼ 1:5

Here, r is the density of plasma, γ is the heat capacity ratio (adiabatic index), σ<sup>0</sup> is the conductivity at T = Ry (it is expressed via atomic constants), and Ry is the Rydberg constant. Energy units have been chosen for temperature.

At initial time t=0, all quantities depend on one space coordinate. It is assumed that the magnetic field has only one component, H ¼ ð Þ 0; 0; Hz . The solution is considered for the problem with initial conditions Eq. (5) and boundary conditions Eq. (6):

$$H\_z(\mathbf{x}, t=0) = \begin{cases} 0 & \text{if} \quad \mathbf{x} < 0 \\ H\_0 & \text{if} \quad \mathbf{x} > 0 \end{cases}, T(\mathbf{x}, t=0) = T\_0 \ \rho(\mathbf{x}, t=0) = \rho\_0 \tag{5}$$

$$H\_z(\mathbf{x}\rightarrow -\mathbf{e}\bullet, t) = 0,\\ H\_z(\mathbf{x}\rightarrow \mathbf{e}\bullet, t) = H\_0,\\ T(\mathbf{x}\rightarrow \pm \mathbf{e}\bullet, t) = T\_0. \tag{6}$$

For dimensionless variables, hz ¼ Hz=H0, τ ¼ T=T<sup>0</sup> Eq. (4) are reduced to the form:

$$\frac{dh\_z}{dt} = \frac{\partial}{\partial \mathbf{x}} \nu(\tau) \frac{\partial h\_z}{\partial \mathbf{x}},\\\frac{d\tau}{dt} = \eta \nu(\tau) \left(\frac{\partial h\_z}{\partial \mathbf{x}}\right)^2,\\\nu(\tau) = \nu\_0 \tau^{-a},\\\nu\_0 = \frac{c^2}{4\pi\sigma\_0} \left(\frac{T\_0}{Ry}\right)^a,\\\eta = (\gamma - 1)\frac{H\_0^2}{\rho T\_0}.\tag{7}$$

In an infinite region (�∞ < x < ∞), the problem has a self-similar solution depending on the variable <sup>ξ</sup> <sup>¼</sup> <sup>x</sup><sup>=</sup> ffiffiffiffiffiffi <sup>ν</sup>0<sup>t</sup> <sup>p</sup> . The solution can be obtained by integrating the system of ordinary differential equations:

$$\frac{\xi}{2}\frac{dh\_z}{d\xi} + \frac{d}{d\xi}\left(\tau^{-\alpha}\frac{dh\_z}{d\xi}\right) = 0,\\ \frac{\xi}{2}\frac{d\tau}{d\xi} + \eta\tau^{-\alpha}\left(\frac{dh\_z}{d\xi}\right)^2 = 0. \tag{8}$$

with boundary conditions

LÞ ¼ 1 is imposed on the right boundary z = L. On the left boundary z=0, the field components take constant values according to Eq. (3). In simulations with 2D and 3D codes, boundary conditions ∂H=∂n ¼ 0 (n is a normal vector to a face) are imposed on lateral faces. By varying parameters ν and β, we can study the effect of the diffusion and Hall terms in Eq. (2) on the diffusion wave parameters. Consider an option with the Hall effect dominating over the diffusion effect: ν = 0, β = 1, Hx0 = Hy0 = Hz0 = 1. Profiles of magnetic field's components Hx, Hy at time t = 0.01, 0.1 are shown in Figures 1 and 2. With the grid refinement, convergence to the reference

3. Diffusion of magnetic field in an immovable plane layer of plasma with regard to joule heating and its effect on the diffusion and heat conduction

The problem of magnetic diffusion in a plane layer of material has many applications in practice [14]. In its detailed formulation, the problem was considered in paper [14] for mega gauss fields. Hydrodynamic motion, magnetic diffusion, heat conduction by electrons, and radiant heat exchange in the "back and forth" approximation were taken into account. Since finding an exact solution to such a problem causes difficulties, the original formulation needs to be simplified. Self-similar solutions to the problem obtained with simplifying assumptions

solution takes place.

coefficients

were also presented in [15].

• heat conduction is absent.

• plasma has Coulomb conductivity,

A model problem is considered with the following assumptions:

• plasma is immovable, it has a constant heat capacity,

Figure 2. Profiles of field components at time t = 0.1: (a) Hx and (b) Hy.

220 Numerical Simulations in Engineering and Science

$$h\_z(\xi \to -\infty) = 0, \quad h\_z(\xi \to \infty) = 1, \quad \tau(\xi \to \pm \infty) = 1. \tag{9}$$

Note that for a linear case, α = 0, the solution of Eqs. (8), (9) can be found in quadratures

$$h\_z(\xi) = 0.5(1 + \text{sign}(\xi)\text{erf}(\xi/2)), \tau(\xi) = 1 - \frac{\eta}{4\pi} \text{Ei}(-\xi^2/4). \tag{10}$$

Since Eið Þ¼ �<sup>x</sup> <sup>C</sup> <sup>þ</sup> ln <sup>x</sup> <sup>þ</sup> <sup>P</sup> i ð Þ �<sup>1</sup> <sup>i</sup> xi <sup>i</sup>�i! , temperature in the vicinity of interface <sup>ξ</sup><sup>2</sup> ~0 for the linear case <sup>α</sup> = 0 has the logarithmic profile τ ξð Þ��<sup>η</sup> ln <sup>ξ</sup><sup>2</sup> =4π.

In general, if α > 0, one does not manage to establish the asymptotic law for temperature in the vicinity of ξ<sup>2</sup> ~0, because of no integral curves of Eq. (8) satisfying the boundary conditions at infinity Eq. (9).

Now, let us build the reference solution to the problem with regard to heat conduction. In this case temperature near the interface takes a finite value. The diffusion equations and energy equation of magnetic field with regard to Joule heating and heat conduction are considered. As it was assumed earlier, all quantities depend on one space coordinate, and the magnetic field has only one component, H ¼ ð Þ 0; 0; Hz . For dimensionless variables, hz = Hz/H0 and τ = T/T0 equations are reduced to the forms

$$\frac{dh\_z}{dt} = \frac{\partial}{\partial \mathbf{x}} \nu\_0 \tau^{-\alpha} \frac{\partial h\_z}{\partial \mathbf{x}},\\\frac{d\tau}{dt} = \eta \nu\_0 \tau^{-\alpha} \left(\frac{\partial h\_z}{\partial \mathbf{x}}\right)^2 + \frac{\partial}{\partial \mathbf{x}} \kappa\_0 \tau^{\beta} \frac{\partial \tau}{\partial \mathbf{x}}.\tag{11}$$

Confine oneself to the consideration of case a = 0.5. Introduce new variables Ψð Þ¼ ξ Ψ ξð Þ

<sup>d</sup><sup>ξ</sup> ¼ � <sup>W</sup>τ�<sup>β</sup>

<sup>d</sup><sup>ξ</sup> <sup>¼</sup> <sup>ξ</sup><sup>W</sup> <sup>1</sup> � <sup>τ</sup>�<sup>β</sup> 2a

=2 with regard to boundary conditions Eq. (15). The

<sup>a</sup> exp � <sup>ξ</sup><sup>2</sup>

2 ,

Benchmarks for Non-Ideal Magnetohydrodynamics http://dx.doi.org/10.5772/intechopen.75713

1,

(16)

223

(17)

<sup>þ</sup> <sup>η</sup>Ψ<sup>2</sup> τα

exp ξ<sup>2</sup>

<sup>=</sup><sup>4</sup> , <sup>W</sup>ð Þ¼ <sup>ξ</sup> <sup>w</sup>ð Þ <sup>ξ</sup> exp �ξ<sup>2</sup>

dhz

dΨ

with boundary conditions:

the self-similar coordinate, x = ξ.

Figure 3. Profiles of self-similar functions: (a) W and (b) Ψ.

replacement of variables gives us the equation system:

<sup>d</sup><sup>ξ</sup> <sup>¼</sup> <sup>Ψ</sup>τ<sup>α</sup> exp � <sup>ξ</sup><sup>2</sup>

<sup>d</sup><sup>ξ</sup> ¼ � <sup>ξ</sup><sup>Ψ</sup> <sup>1</sup> � τα ð Þ

4 , <sup>d</sup><sup>τ</sup>

<sup>2</sup> , dW

hzð Þ¼ ξ ¼ ξ<sup>1</sup> 0:5 1ð Þ þ erfð Þ ξ1=2 , τ ξð Þ¼ ¼ ξ<sup>1</sup> 1,

hzð Þ¼ ξ ¼ 0 0:5, wð Þ¼ ξ ¼ 0 0:

parameter ξ<sup>1</sup> from ξ<sup>1</sup> = 10 to ξ<sup>1</sup> = 1 does not affect simulation results.

<sup>Ψ</sup>ð Þ¼ <sup>ξ</sup> <sup>¼</sup> <sup>ξ</sup><sup>1</sup> <sup>c</sup>1, Wð Þ¼ <sup>ξ</sup> <sup>¼</sup> <sup>ξ</sup><sup>1</sup> <sup>c</sup><sup>2</sup> <sup>þ</sup> ηξ1c<sup>2</sup>

The set of Eqs. (16), (17) was solved numerically with the methods of automatically selecting an integration step. The following values of parameters were used in simulations: ξ<sup>1</sup> = 10, η = 20/3, α = 3/2, D<sup>0</sup> = 1, and k<sup>0</sup> = aD<sup>0</sup> = 1/2. The values of constants satisfying the boundary conditions Eq. (17) were obtained: C1 = 0.10231 and C2 = 1.79474. Since the behavior of functions near the right boundary corresponds to asymptotic laws Eq. (17), the reduction of

Results of simulations are illustrated in Figures 3 and 4. With the use of such regularity method (with regard to heat conduction), temperature at the central point of the computational domain takes its finite value. Note that with t = 1/ν<sup>0</sup> the space coordinate coincides with

A self-similar solution depending on the variable <sup>ξ</sup> <sup>¼</sup> <sup>x</sup><sup>=</sup> ffiffiffiffiffiffi <sup>ν</sup>0<sup>t</sup> <sup>p</sup> can be obtained by integrating the system of ordinary differential equations:

$$\frac{\xi}{2}\frac{dh\_z}{d\xi} + \frac{d}{d\xi}\left(\tau^{-a}\frac{dh\_z}{d\xi}\right) = 0, \quad \frac{\xi}{2}\frac{d\tau}{d\xi} + \eta\tau^{-a}\left(\frac{dh\_z}{d\xi}\right)^2 + a\frac{d}{d\xi}\tau^{\beta}\frac{d\tau}{d\xi} = 0, \quad a = \frac{\kappa\_0}{D\_0},\tag{12}$$

with boundary conditions:

$$h\_z(\xi \to -\infty) = 0, \quad h\_z(\xi \to \infty) = 1, \quad \pi(\xi \to \pm \infty) = 1. \tag{13}$$

To find the reference solution, it is convenient to use the first-order system with an increased number of unknowns instead of the second-order system Eq. (12). The first-order system relative to variables hz, <sup>τ</sup>, <sup>Ψ</sup> <sup>¼</sup> <sup>τ</sup>�<sup>α</sup>dhz=dξ, w ¼ �aτ<sup>β</sup>dτ=d<sup>ξ</sup> looks like

$$\frac{d\hbar\_z}{d\xi} = \Psi \tau^a, \quad \frac{d\tau}{d\xi} = -\frac{w\tau^{-\beta}}{a}, \quad \frac{d\Psi}{d\xi} = -\frac{\xi\Psi\tau^a}{2}, \quad \frac{dw}{d\xi} = -\frac{\xi w\tau^{-\beta}}{2a} + \eta\Psi^2\tau^a. \tag{14}$$

Consider the numerical solution of Eq. (14) for the right half plane (0 < ξ < ∞). The solution in the left half plane ð Þ �∞ < ξ < 0 follows from the symmetry conditions:

$$h\_z(-\xi) = 1 - h\_z(\xi), \quad \tau(-\xi) = \tau(\xi), \quad \Psi(-\xi) = \Psi(\xi), \quad w(-\xi) = -w(\xi).$$

Consider the numerical solution of Eq. (14) in a bounded domain 0ð Þ < ξ < ξ<sup>1</sup> . To formulate boundary conditions for this bounded domain, it is required to find the asymptotic behavior of functions with ξ ! ∞. Asymptotic laws can be formulated, if we assume a = 0.5. In this case boundary conditions have the form:

$$\begin{aligned} \eta\_z(\xi \to \infty) &= 0.5(1 + \text{erf}(\xi/2)), \quad \tau(\xi \to \infty) = 1, \\ \Psi(\xi \to \infty) &= c\_1 \exp\left(-\xi^2/4\right), \quad w(\xi \to \infty) = \left(c\_2 + \eta \xi c\_1^2\right) \exp\left(-\xi^2/2\right). \end{aligned} \tag{15}$$

Constants C1, C2 are taken so that the following conditions are satisfied on the left boundary of the computational domain:

$$h\_z(\xi=0) = 0.5, \quad w(\xi=0) = 0.5$$

Confine oneself to the consideration of case a = 0.5. Introduce new variables Ψð Þ¼ ξ Ψ ξð Þ exp ξ<sup>2</sup> <sup>=</sup><sup>4</sup> , <sup>W</sup>ð Þ¼ <sup>ξ</sup> <sup>w</sup>ð Þ <sup>ξ</sup> exp �ξ<sup>2</sup> =2 with regard to boundary conditions Eq. (15). The replacement of variables gives us the equation system:

$$\begin{aligned} \frac{dh\_z}{d\xi} &= \Psi \tau^a \exp\left(-\frac{\xi^2}{4}\right), & \frac{d\tau}{d\xi} &= -\frac{W\tau^{-\beta}}{a} \exp\left(-\frac{\xi^2}{2}\right), \\\frac{d\Psi}{d\xi} &= -\frac{\xi\Psi(1-\tau^a)}{2}, & \frac{dW}{d\xi} &= \frac{\xi W(1-\tau^{-\beta})}{2a} + \eta\Psi^2\tau^a \end{aligned} \tag{16}$$

with boundary conditions:

Now, let us build the reference solution to the problem with regard to heat conduction. In this case temperature near the interface takes a finite value. The diffusion equations and energy equation of magnetic field with regard to Joule heating and heat conduction are considered. As it was assumed earlier, all quantities depend on one space coordinate, and the magnetic field has only one component, H ¼ ð Þ 0; 0; Hz . For dimensionless variables, hz = Hz/H0 and τ = T/T0

dt <sup>¼</sup> ην0τ�<sup>α</sup> <sup>∂</sup>hz

<sup>d</sup><sup>ξ</sup> <sup>þ</sup> ητ�<sup>α</sup> dhz

To find the reference solution, it is convenient to use the first-order system with an increased number of unknowns instead of the second-order system Eq. (12). The first-order system

<sup>d</sup><sup>ξ</sup> ¼ � ξΨτ<sup>α</sup>

Consider the numerical solution of Eq. (14) for the right half plane (0 < ξ < ∞). The solution in

hzð Þ¼ �ξ 1 � hzð Þ ξ , τð Þ¼ �ξ τ ξð Þ, Ψð Þ¼ �ξ Ψð Þ ξ , wð Þ¼� �ξ wð Þ ξ :

Consider the numerical solution of Eq. (14) in a bounded domain 0ð Þ < ξ < ξ<sup>1</sup> . To formulate boundary conditions for this bounded domain, it is required to find the asymptotic behavior of functions with ξ ! ∞. Asymptotic laws can be formulated, if we assume a = 0.5. In this case

<sup>=</sup><sup>4</sup> � �, wð Þ¼ <sup>ξ</sup> ! <sup>∞</sup> <sup>c</sup><sup>2</sup> <sup>þ</sup> ηξc<sup>2</sup>

Constants C1, C2 are taken so that the following conditions are satisfied on the left boundary of

hzð Þ¼ ξ ¼ 0 0:5, wð Þ¼ ξ ¼ 0 0:

dξ � �<sup>2</sup>

∂x � �<sup>2</sup>

þ ∂ ∂x

þ a d <sup>d</sup><sup>ξ</sup> <sup>τ</sup><sup>β</sup> <sup>d</sup><sup>τ</sup>

hzð Þ¼ ξ ! �∞ 0, hzð Þ¼ ξ ! ∞ 1, τ ξð Þ¼ ! �∞ 1: (13)

<sup>2</sup> , dw

<sup>d</sup><sup>ξ</sup> ¼ � <sup>ξ</sup>wτ�<sup>β</sup> 2a

1 � � exp �ξ<sup>2</sup>

<sup>κ</sup>0τ<sup>β</sup> <sup>∂</sup><sup>τ</sup> ∂x

<sup>ν</sup>0<sup>t</sup> <sup>p</sup> can be obtained by integrating

<sup>d</sup><sup>ξ</sup> <sup>¼</sup> <sup>0</sup>, a <sup>¼</sup> <sup>κ</sup><sup>0</sup>

<sup>þ</sup> ηΨ<sup>2</sup>

: (11)

D<sup>0</sup>

, (12)

τ<sup>α</sup>: (14)

<sup>=</sup><sup>2</sup> � �: (15)

equations are reduced to the forms

222 Numerical Simulations in Engineering and Science

ξ 2 dhz dξ þ

with boundary conditions:

dhz

<sup>d</sup><sup>ξ</sup> <sup>¼</sup> Ψτα, <sup>d</sup><sup>τ</sup>

boundary conditions have the form:

the computational domain:

Ψ ξð Þ¼ ! <sup>∞</sup> <sup>c</sup><sup>1</sup> exp �ξ<sup>2</sup>

dhz dt <sup>¼</sup> <sup>∂</sup> ∂x

the system of ordinary differential equations:

d <sup>d</sup><sup>ξ</sup> <sup>τ</sup>�<sup>α</sup> dhz dξ � �

<sup>ν</sup>0τ�<sup>α</sup> <sup>∂</sup>hz ∂x , dτ

> <sup>¼</sup> <sup>0</sup>, <sup>ξ</sup> 2 dτ

relative to variables hz, <sup>τ</sup>, <sup>Ψ</sup> <sup>¼</sup> <sup>τ</sup>�<sup>α</sup>dhz=dξ, w ¼ �aτ<sup>β</sup>dτ=d<sup>ξ</sup> looks like

the left half plane ð Þ �∞ < ξ < 0 follows from the symmetry conditions:

hzð Þ¼ ξ ! ∞ 0:5 1ð Þ þ erfð Þ ξ=2 , τ ξð Þ¼ ! ∞ 1,

<sup>a</sup> , <sup>d</sup><sup>Ψ</sup>

<sup>d</sup><sup>ξ</sup> ¼ � <sup>w</sup>τ�<sup>β</sup>

A self-similar solution depending on the variable <sup>ξ</sup> <sup>¼</sup> <sup>x</sup><sup>=</sup> ffiffiffiffiffiffi

$$\begin{aligned} h\_{\mathfrak{z}}(\xi = \xi\_1) &= 0.5(1 + \text{erf}(\xi\_1/2)), & \tau(\xi = \xi\_1) &= 1, \\ \Psi(\xi = \xi\_1) &= \mathbf{c}\_1, & W(\xi = \xi\_1) &= \mathbf{c}\_2 + \eta \xi\_1 \mathbf{c}\_1^2, \\ h\_{\mathfrak{z}}(\xi = 0) &= 0.5, & w(\xi = 0) &= 0. \end{aligned} \tag{17}$$

The set of Eqs. (16), (17) was solved numerically with the methods of automatically selecting an integration step. The following values of parameters were used in simulations: ξ<sup>1</sup> = 10, η = 20/3, α = 3/2, D<sup>0</sup> = 1, and k<sup>0</sup> = aD<sup>0</sup> = 1/2. The values of constants satisfying the boundary conditions Eq. (17) were obtained: C1 = 0.10231 and C2 = 1.79474. Since the behavior of functions near the right boundary corresponds to asymptotic laws Eq. (17), the reduction of parameter ξ<sup>1</sup> from ξ<sup>1</sup> = 10 to ξ<sup>1</sup> = 1 does not affect simulation results.

Results of simulations are illustrated in Figures 3 and 4. With the use of such regularity method (with regard to heat conduction), temperature at the central point of the computational domain takes its finite value. Note that with t = 1/ν<sup>0</sup> the space coordinate coincides with the self-similar coordinate, x = ξ.

Figure 3. Profiles of self-similar functions: (a) W and (b) Ψ.

hrð Þ¼ r; t

problem [17].

Since <sup>r</sup> r0 ¼ r2 <sup>0</sup>∂r<sup>0</sup>

<sup>1</sup> � <sup>β</sup> rFð Þ<sup>t</sup> r � �<sup>3</sup>

8 >>><

>>>:

netic field in internal domain.

the solution of Eq. (19) looks like

Hxð Þ¼ x; y; z; t

hrð Þ¼ r; t

hϑð Þ¼ r; t

hφð Þ¼ r; t

dx

RF = 13.467.

Hzð Þ¼ x; y; z; t Hz<sup>0</sup>

depending on one space coordinate:

Hr Hz<sup>0</sup> cos <sup>ϑ</sup> <sup>¼</sup> <sup>1</sup>

H<sup>ϑ</sup> Hz<sup>0</sup> sin <sup>ϑ</sup> ¼ � <sup>1</sup>

H<sup>φ</sup> Hz<sup>0</sup> sin <sup>ϑ</sup> <sup>¼</sup> <sup>1</sup>

dz <sup>¼</sup> xz hð Þ <sup>r</sup>ð Þþ <sup>r</sup>; <sup>t</sup> <sup>h</sup>ϑð Þ <sup>r</sup>; <sup>t</sup>

div<sup>H</sup> <sup>¼</sup> Hz<sup>0</sup> cos <sup>ϑ</sup> <sup>1</sup>

<sup>1</sup> � <sup>β</sup> � � <sup>r</sup>0ð Þ <sup>r</sup>;<sup>t</sup> r � �<sup>2</sup>

r2

<sup>r</sup>2∂<sup>r</sup> , then 1 � <sup>β</sup> <sup>¼</sup> <sup>1</sup> <sup>þ</sup> <sup>β</sup>

Hz0xz

z2

Hz<sup>0</sup>

Hz<sup>0</sup>

�r<sup>2</sup>hϑð Þþ <sup>r</sup>; <sup>t</sup> <sup>z</sup><sup>2</sup>ð Þ hrð Þþ <sup>r</sup>; <sup>t</sup> <sup>h</sup>ϑð Þ <sup>r</sup>; <sup>t</sup> ,

Hz<sup>0</sup>

The magnetic field lines can be obtained by integrating equations

0 @

∂r<sup>2</sup>hr ∂r þ 2 hϑ r

� �

2 � � <sup>γ</sup>�<sup>1</sup>

, r > rFð Þt

, r ≤ rFð Þt

, hϑð Þ¼ r; t

Here, the functions r0ð Þ r; t , rð Þ r; t are defined from the self-similar solution to the point blast

Unknown constant β can be found from the condition of the solenoidal distribution of mag-

<sup>γ</sup>þ<sup>1</sup>. That is why <sup>β</sup> <sup>¼</sup> <sup>4</sup>

Note that in external domain this condition is satisfied automatically. In Cartesian coordinates,

It is convenient to compare the numerical and exact solutions using the field components

�yHxð Þþ x; y; z; t xHyð Þ x; y; z; t r � �

dy

We consider the process stage, at which the numerical simulation becomes self-similar. In this case, the shock wave is considerably far (compared to the energy release region) from the blast center. For example, at the final time t = 3, the wave front is located at a distance of

¼ 2 Hz<sup>0</sup> r0

<sup>r</sup><sup>2</sup> ð Þ hrð Þþ <sup>r</sup>; <sup>t</sup> <sup>h</sup>ϑð Þ <sup>r</sup>; <sup>t</sup> , Hyð Þ¼ <sup>x</sup>; <sup>y</sup>; <sup>z</sup>; <sup>t</sup>

<sup>r</sup><sup>2</sup> <sup>ð</sup>hrð Þþ <sup>r</sup>; <sup>t</sup> <sup>h</sup>ϑð Þ <sup>r</sup>; <sup>t</sup> Þ � <sup>h</sup>ϑð Þ <sup>r</sup>; <sup>t</sup> � �:

Hzð Þþ x; y; z; t

Hzð Þ� x; y; z; t

�<sup>1</sup> � <sup>β</sup> 2

8 >>><

>>>:

cos <sup>ϑ</sup> <sup>1</sup> � <sup>β</sup> � � <sup>r</sup><sup>2</sup>

� <sup>1</sup> <sup>þ</sup> <sup>β</sup> 2

rFð Þt r � �<sup>3</sup>

<sup>0</sup>∂r<sup>0</sup>

Hz0yz

xHxð Þþ x; y; z; t yHyð Þ x; y; z; t z

z xHxð Þþ x; y; z; t yHyð Þ x; y; z; t � �

r<sup>2</sup> � z<sup>2</sup>

dz <sup>¼</sup> yz hð Þ <sup>r</sup>ð Þþ <sup>r</sup>; <sup>t</sup> <sup>h</sup>ϑð Þ <sup>r</sup>; <sup>t</sup>

¼ 0:

�r<sup>2</sup>hϑð Þþ <sup>r</sup>; <sup>t</sup> <sup>z</sup><sup>2</sup>ð Þ hrð Þþ <sup>r</sup>; <sup>t</sup> <sup>h</sup>ϑð Þ <sup>r</sup>; <sup>t</sup> : (22)

!

<sup>3</sup>γþ<sup>1</sup> :

<sup>r</sup><sup>2</sup>∂<sup>r</sup> � <sup>1</sup> <sup>þ</sup> <sup>β</sup>

� �

2 � � r

<sup>r</sup><sup>2</sup> ð Þ hrð Þþ <sup>r</sup>; <sup>t</sup> <sup>h</sup>ϑð Þ <sup>r</sup>; <sup>t</sup> ,

,

1 A

, (21)

r0

γ � 1 γ þ 1

¼ 0:

(20)

� � <sup>r</sup>ð Þ <sup>r</sup>; <sup>t</sup> <sup>r</sup>ð Þ <sup>γ</sup> � <sup>1</sup>

, r > rFð Þt

:

225

<sup>r</sup>0r0ð Þ <sup>r</sup>; <sup>t</sup> ð Þ <sup>γ</sup> <sup>þ</sup> <sup>1</sup> , r <sup>≤</sup> rFð Þ<sup>t</sup>

Benchmarks for Non-Ideal Magnetohydrodynamics http://dx.doi.org/10.5772/intechopen.75713

Figure 4. Profiles: (a) magnetic field and (b) temperature.

### 4. A point explosion in a perfectly non-conducting atmosphere

Let us consider the problem of a point blast in the presence of a uniform magnetic field (for definiteness) along the z axis (Hz = Hz0 = 0.01) in a perfectly non-conducting atmosphere. Initial data are chosen in such a way that the field has no effect on the motion of matter ε0ð Þ r2=r<sup>1</sup> <sup>3</sup> >> H<sup>2</sup> <sup>z</sup>0=r0. It is assumed that behind the shock front, the medium becomes perfectly conducting. If a self-similar solution to the problem of a point explosion is known, then one can calculate the magnetic field components at some time t > 0. In an external domain ð Þ r > rFð Þt , the magnetic field's vector potential is the solution to stationary equation rot rot ð Þ¼ Ψ 0. With regard to conditions at infinity, this solution takes the form [16]:

$$\Psi\_r(r,\theta,t) = 0, \quad \Psi\_\psi(r,\theta,t) = 0.5H\_{z0}r(1-\mathcal{C}(t)/r^3)\sin\theta, \quad \Psi\_\delta(r,\theta,t) = 0. \tag{18}$$

Here, unknown constant C tðÞ¼ br<sup>3</sup> <sup>F</sup>ð Þt can be found from the condition of coupling with the solution in an internal domain ð Þ r < rFð Þt . Write components of magnetic field H ¼ rotΨ:

$$H\_r(r, \\$, t) = H\_{\text{z0}} \left( 1 - \mathbb{C}(t) / r^3 \right) \sin \theta, \quad H\_{\text{\overline{\theta}}}(r, \\$, t) = 0, \quad H\_{\text{\overline{\theta}}}(r, \\$, t) = -H\_{\text{z0}} \left( 1 + 0.5 \mathbb{C}(t) / r^3 \right) \cos \theta.$$

The solution in the internal domain (r<rF(t)) is found from the freezing-in condition of the magnetic field:

$$\frac{d}{dt}\left(\frac{H\_r}{\rho}\right) = \frac{H\_r}{\rho}\frac{\partial u\_r}{\partial r}, \quad \frac{d}{dt}\left(\frac{H\_\vartheta}{\rho}\right) = \frac{H\_\vartheta}{\rho}\frac{u\_r}{r}, \quad \frac{d}{dt}\left(\frac{H\_\vartheta}{\rho}\right) = \frac{H\_\vartheta}{\rho}\frac{u\_r}{r}.$$

The integration of these equations with regard to the solution in external domain Eq. (18) gives us

$$H\_r(r, \theta, t) = H\_{z0} h\_r(r, t) \cos \theta, \; H\_{\theta}(r, \theta, t) = H\_{z0} h\_{\theta}(r, t) \sin \theta, \; H\_{\theta}(r, \theta, t) = 0,\tag{19}$$

where

Benchmarks for Non-Ideal Magnetohydrodynamics http://dx.doi.org/10.5772/intechopen.75713 225

:

$$h\_{\mathbf{r}}(r,t) = \begin{cases} 1 - \beta \left(\frac{r\_{\mathbf{f}}(t)}{r}\right)^3, & r > r\_{\mathbf{f}}(t) \\ \left(1 - \beta\right) \left(\frac{r\_{\mathbf{f}}(r,t)}{r}\right)^2, & r \le r\_{\mathbf{f}}(t) \end{cases},\\ h\_{\vartheta}(r,t) = \begin{cases} -1 - \frac{\beta}{2} \left(\frac{r\_{\mathbf{f}}(t)}{r}\right)^3, & r > r\_{\mathbf{f}}(t) \\ -\left(1 + \frac{\beta}{2}\right) \frac{\rho(r,t)r(\mathbf{\check{y}} - 1)}{\rho\_0 r\_0(r,t)(\mathbf{\check{y}} + 1)}, & r \le r\_{\mathbf{f}}(t) \end{cases}$$

Here, the functions r0ð Þ r; t , rð Þ r; t are defined from the self-similar solution to the point blast problem [17].

Unknown constant β can be found from the condition of the solenoidal distribution of magnetic field in internal domain.

$$\text{div}\mathbf{H} = H\_{z0}\cos\theta\left(\frac{1}{r^2}\frac{\partial r^2 h\_r}{\partial r} + 2\frac{h\_\theta}{r}\right) = 2\frac{H\_{z0}}{r\_0}\cos\theta\left((1-\beta)\frac{r\_0^2 \partial r\_0}{r^2 \partial r} - \left(1 + \frac{\beta}{2}\right)\frac{\rho}{\rho\_0}\frac{\gamma - 1}{\gamma + 1}\right) = 0.$$

$$\text{Since }\frac{\rho}{\rho\_0} = \frac{r\_0^2 \partial r\_0}{r^2 \partial r}\text{ then }1 - \beta = \left(1 + \frac{\beta}{2}\right)\frac{\gamma - 1}{\gamma + 1}. \text{ That is why }\beta = \frac{4}{3\gamma + 1}.$$

Note that in external domain this condition is satisfied automatically. In Cartesian coordinates, the solution of Eq. (19) looks like

$$\begin{split}H\_{\mathbf{x}}(\mathbf{x},y,z,t) &= \frac{H\_{\mathbf{z}0}\mathbf{x}z}{r^2} (h\_r(r,t) + h\_{\vartheta}(r,t)), \quad H\_y(\mathbf{x},y,z,t) = \frac{H\_{\mathbf{z}0}yz}{r^2} (h\_r(r,t) + h\_{\vartheta}(r,t)),\\H\_{\mathbf{z}}(\mathbf{x},y,z,t) &= H\_{\mathbf{z}0} \left(\frac{z^2}{r^2} (h\_r(r,t) + h\_{\vartheta}(r,t)) - h\_{\vartheta}(r,t)\right). \end{split} \tag{20}$$

It is convenient to compare the numerical and exact solutions using the field components depending on one space coordinate:

$$\begin{split} h\_{r}(r,t) &= \frac{H\_{r}}{H\_{z0}\cos\mathfrak{S}} = \frac{1}{H\_{z0}} \Biggl( H\_{z}(x,y,z,t) + \frac{xH\_{x}(x,y,z,t) + yH\_{y}(x,y,z,t)}{z} \Biggr), \\ h\_{\theta}(r,t) &= \frac{H\_{\theta}}{H\_{z0}\sin\mathfrak{S}} = -\frac{1}{H\_{z0}} \Biggl( H\_{z}(x,y,z,t) - \frac{z\Big(xH\_{x}(x,y,z,t) + yH\_{y}(x,y,z,t)\Big)}{r^{2} - z^{2}} \Biggr), \\ h\_{\psi}(r,t) &= \frac{H\_{\psi}}{H\_{z0}\sin\mathfrak{S}} = \frac{1}{H\_{z0}} \Biggl( \frac{-yH\_{x}(x,y,z,t) + xH\_{y}(x,y,z,t)}{r} \Biggr) = 0. \end{split} \tag{21}$$

The magnetic field lines can be obtained by integrating equations

4. A point explosion in a perfectly non-conducting atmosphere

regard to conditions at infinity, this solution takes the form [16]:

Ψrð Þ¼ r; ϑ; t 0, Ψφð Þ¼ r; ϑ; t 0:5Hz0r 1 � C tð Þ=r

ε0ð Þ r2=r<sup>1</sup>

<sup>3</sup> >> H<sup>2</sup>

Here, unknown constant C tðÞ¼ br<sup>3</sup>

Figure 4. Profiles: (a) magnetic field and (b) temperature.

224 Numerical Simulations in Engineering and Science

d dt

Hr r  <sup>¼</sup> Hr r ∂ur <sup>∂</sup><sup>r</sup> , <sup>d</sup> dt

Hrð Þ¼ r; ϑ; t Hz<sup>0</sup> 1 � C tð Þ=r

magnetic field:

where

Let us consider the problem of a point blast in the presence of a uniform magnetic field (for definiteness) along the z axis (Hz = Hz0 = 0.01) in a perfectly non-conducting atmosphere. Initial data are chosen in such a way that the field has no effect on the motion of matter

conducting. If a self-similar solution to the problem of a point explosion is known, then one can calculate the magnetic field components at some time t > 0. In an external domain ð Þ r > rFð Þt , the magnetic field's vector potential is the solution to stationary equation rot rot ð Þ¼ Ψ 0. With

solution in an internal domain ð Þ r < rFð Þt . Write components of magnetic field H ¼ rotΨ:

H<sup>ϑ</sup> r 

The solution in the internal domain (r<rF(t)) is found from the freezing-in condition of the

The integration of these equations with regard to the solution in external domain Eq. (18) gives us

<sup>3</sup> sin <sup>ϑ</sup>, Hφð Þ¼ <sup>r</sup>; <sup>ϑ</sup>; <sup>t</sup> <sup>0</sup>, Hϑð Þ¼� <sup>r</sup>; <sup>ϑ</sup>; <sup>t</sup> Hz<sup>0</sup> <sup>1</sup> <sup>þ</sup> <sup>0</sup>:5C tð Þ=<sup>r</sup>

<sup>¼</sup> <sup>H</sup><sup>ϑ</sup> r ur <sup>r</sup> , <sup>d</sup> dt

Hrð Þ¼ r; ϑ; t Hz0hrð Þ r; t cos ϑ, Hϑð Þ¼ r; ϑ; t Hz0hϑð Þ r; t sin ϑ, Hφð Þ¼ r; ϑ; t 0, (19)

<sup>z</sup>0=r0. It is assumed that behind the shock front, the medium becomes perfectly

<sup>3</sup> sin <sup>ϑ</sup>, <sup>Ψ</sup>ϑð Þ¼ <sup>r</sup>; <sup>ϑ</sup>; <sup>t</sup> <sup>0</sup>: (18)

<sup>¼</sup> <sup>H</sup><sup>φ</sup> r ur r :

<sup>3</sup> cos ϑ:

<sup>F</sup>ð Þt can be found from the condition of coupling with the

H<sup>φ</sup> r 

$$\frac{dx}{dz} = \frac{xz(h\_r(r,t) + h\_\vartheta(r,t))}{-r^2h\_\vartheta(r,t) + z^2(h\_r(r,t) + h\_\vartheta(r,t))},\\\frac{dy}{dz} = \frac{yz(h\_r(r,t) + h\_\vartheta(r,t))}{-r^2h\_\vartheta(r,t) + z^2(h\_r(r,t) + h\_\vartheta(r,t))}.\tag{22}$$

We consider the process stage, at which the numerical simulation becomes self-similar. In this case, the shock wave is considerably far (compared to the energy release region) from the blast center. For example, at the final time t = 3, the wave front is located at a distance of RF = 13.467.

Figure 5. The magnetic field lines at time t = 3 in plane y = 0. A dashed line shows the shock front. Strong blast in a perfectly nonconducting atmosphere (a) and in a uniform conducting atmosphere (b) [11].

figures. Such a comparison of reference versus numerical solution indicates whether the spherical

Figure 7. Profiles of the non-dimensionalized field components at time t = 3: (a) radial hr and (b) angular hθ.

Simulation setup: The energy release region is a sphere of radius r1 = 0.1, in which the initial specific energy per unit mass is set to ε = ε<sup>0</sup> = 107 and the initial density is set to r = r<sup>0</sup> = 1. In the spherical layer r1 <R<r2 = 15, the specific energy and the density are equal to ε = 0, r = r<sup>0</sup> = 1, respectively. The EOS is P = r(γ�1)ε, γ = 1.4. The problem domain in the three-dimensional

This problem requires taking into account the magnetic field diffusion in external domain (outside the shock front). It is assumed that behind the shock front, the medium becomes perfectly conducting due to ionization effects. The magnetic viscosity is approximated by the

> , ε ≤ 0, ν1ð Þ 1 � ε=ε<sup>1</sup> , 0 < ε < ε<sup>1</sup> ¼ 1,

> > , t ≤ tk ¼ 3 is satisfied. Accounting of diffusion

Benchmarks for Non-Ideal Magnetohydrodynamics http://dx.doi.org/10.5772/intechopen.75713 227

0, ε<sup>1</sup> ≤ ε:

Parameter ε<sup>1</sup> is chosen to provide that the magnetic viscosity behind the shock front is always

in external domain leads to the necessity of increasing the size of computational domain (rF(tk)/ L))3 < <1) in comparison with the ideal MHD problems to be able to set boundary conditions corresponding to the initial undisturbed state. Thus, the size of computational domain must be

The problem formulation and its analytical solution have been taken from [16]. In contrast to

<sup>ν</sup><sup>1</sup> <sup>¼</sup> 105

setup is a cube L � L� L. All boundary faces of the domain are rigid walls.

8 ><

>:

5. Diffusion of magnetic field into a spherical plasma cloud

this paper, consider the diffusion problem (the plasma cloud motion is neglected):

symmetry is preserved during numerical simulations.

ν εð Þ¼ <sup>c</sup><sup>2</sup>

zero, i.e., the condition <sup>ε</sup><sup>1</sup> <sup>&</sup>lt; <sup>ε</sup>Fð Þ<sup>t</sup> <sup>≃</sup>ε0ð Þ <sup>r</sup>1=rFð Þ<sup>t</sup> <sup>3</sup>

almost five times larger, L ≈ 75.

4πσ ¼

following dependence:

The flow parameters in this problem depend on a single spatial variable, r, and the field components, on two variables, r and θ. One can restrict the consideration to any plane passing through the z axis. Magnetic field lines in the plane y = 0 at t = 3 are shown in Figure 5. It follows from this figure that these lines of force stretch along axis z and, thereby, prevent the spread of gas in the direction orthogonal to this axis. This effect is small in the given problem because of the field smallness. With an increased strength of the field, the pressure zone gets out of its spherical shape due to occurrence of the singled out direction.

Profiles of the field components Hy and Hz along the line x=z in this plane are presented in Figure 6. The self-similar functions of fluid parameters and radial and angular field components hr and h<sup>θ</sup> depend on one spatial variable r. These functions are shown in Figure 7. In testing numerical methods, the values of grid functions for all cells in the domain can be mapped onto the

Figure 6. Profiles of field components along line x = z and y = 0 at time t = 3: (a) Hx and (b) Hz.

Figure 7. Profiles of the non-dimensionalized field components at time t = 3: (a) radial hr and (b) angular hθ.

The flow parameters in this problem depend on a single spatial variable, r, and the field components, on two variables, r and θ. One can restrict the consideration to any plane passing through the z axis. Magnetic field lines in the plane y = 0 at t = 3 are shown in Figure 5. It follows from this figure that these lines of force stretch along axis z and, thereby, prevent the spread of gas in the direction orthogonal to this axis. This effect is small in the given problem because of the field smallness. With an increased strength of the field, the pressure zone gets

Figure 5. The magnetic field lines at time t = 3 in plane y = 0. A dashed line shows the shock front. Strong blast in a

Profiles of the field components Hy and Hz along the line x=z in this plane are presented in Figure 6. The self-similar functions of fluid parameters and radial and angular field components hr and h<sup>θ</sup> depend on one spatial variable r. These functions are shown in Figure 7. In testing numerical methods, the values of grid functions for all cells in the domain can be mapped onto the

out of its spherical shape due to occurrence of the singled out direction.

Figure 6. Profiles of field components along line x = z and y = 0 at time t = 3: (a) Hx and (b) Hz.

perfectly nonconducting atmosphere (a) and in a uniform conducting atmosphere (b) [11].

226 Numerical Simulations in Engineering and Science

figures. Such a comparison of reference versus numerical solution indicates whether the spherical symmetry is preserved during numerical simulations.

Simulation setup: The energy release region is a sphere of radius r1 = 0.1, in which the initial specific energy per unit mass is set to ε = ε<sup>0</sup> = 107 and the initial density is set to r = r<sup>0</sup> = 1. In the spherical layer r1 <R<r2 = 15, the specific energy and the density are equal to ε = 0, r = r<sup>0</sup> = 1, respectively. The EOS is P = r(γ�1)ε, γ = 1.4. The problem domain in the three-dimensional setup is a cube L � L� L. All boundary faces of the domain are rigid walls.

This problem requires taking into account the magnetic field diffusion in external domain (outside the shock front). It is assumed that behind the shock front, the medium becomes perfectly conducting due to ionization effects. The magnetic viscosity is approximated by the following dependence:

$$\nu(\varepsilon) = \frac{\varepsilon^2}{4\pi\sigma} = \begin{cases} \nu\_1 = 10^5 \varkappa & \varepsilon \le 0, \\ \nu\_1 (1 - \varepsilon/\varepsilon\_1) \searrow & 0 < \varepsilon < \varepsilon\_1 = 1, \\ 0, & \varepsilon\_1 \le \varepsilon. \end{cases}$$

Parameter ε<sup>1</sup> is chosen to provide that the magnetic viscosity behind the shock front is always zero, i.e., the condition <sup>ε</sup><sup>1</sup> <sup>&</sup>lt; <sup>ε</sup>Fð Þ<sup>t</sup> <sup>≃</sup>ε0ð Þ <sup>r</sup>1=rFð Þ<sup>t</sup> <sup>3</sup> , t ≤ tk ¼ 3 is satisfied. Accounting of diffusion in external domain leads to the necessity of increasing the size of computational domain (rF(tk)/ L))3 < <1) in comparison with the ideal MHD problems to be able to set boundary conditions corresponding to the initial undisturbed state. Thus, the size of computational domain must be almost five times larger, L ≈ 75.

### 5. Diffusion of magnetic field into a spherical plasma cloud

The problem formulation and its analytical solution have been taken from [16]. In contrast to this paper, consider the diffusion problem (the plasma cloud motion is neglected):

$$\frac{\partial H}{\partial t} = -rot(\nu \cdot rot \mathbf{H} + b[\mathbf{H} \times rot \mathbf{H}]).\tag{23}$$

It is assumed that the magnetic field at infinity is uniform and directed along axis z: Hð0; 0; H<sup>0</sup> <sup>¼</sup> ffiffiffi 2 <sup>p</sup> <sup>Þ</sup> (see Figure 8). The magnetic viscosity coefficient is constant inside and outside the cloud: νð Þ¼ r ν1, r < r<sup>0</sup> ¼ 1 ν2, r > r<sup>0</sup> � .

### 5.1. Diffusion of magnetic field in the absence of the Hall effect

Assume that the Hall effect contribution is small, bH0=ν ≪ 1. Write the equation of diffusion relative to vector potential H ¼ rotΨ:

$$\frac{\partial \Psi}{\partial t} = -\nu \cdot rot \,\,rot \Psi. \tag{24}$$

For finite values of conductivity in external domain ν<sup>2</sup> > 0, the limit numerical solution of

Components of magnetic field H ¼ rotΨ are found by differentiating vector potential Eq. (29).

<sup>¼</sup> <sup>2</sup>f rð Þ ; <sup>t</sup> , hϑð Þ¼ <sup>r</sup>; <sup>t</sup> <sup>∂</sup>r<sup>2</sup><sup>f</sup> <sup>=</sup>r∂r, hφð Þ¼ <sup>r</sup>; <sup>t</sup> <sup>0</sup> (31)

<sup>r</sup><sup>2</sup> ð Þ hrð Þ� <sup>r</sup>; <sup>t</sup> <sup>h</sup>ϑð Þ <sup>r</sup>; <sup>t</sup> , Hyð Þ <sup>x</sup>; <sup>y</sup>; <sup>z</sup>; <sup>t</sup>

:

xHxð Þþ x; y; z; t yHyð Þ x; y; z; t

z xHxð Þþ x; y; z; t yHyð Þ x; y; z; t 

<sup>r</sup> <sup>¼</sup> <sup>0</sup>:

<sup>z</sup> ,

<sup>r</sup><sup>2</sup> � <sup>z</sup><sup>2</sup> ,

<sup>r</sup><sup>2</sup> ð Þ hrð Þ� <sup>r</sup>; <sup>t</sup> <sup>h</sup>ϑð Þ <sup>r</sup>; <sup>t</sup> ,

Benchmarks for Non-Ideal Magnetohydrodynamics http://dx.doi.org/10.5772/intechopen.75713

(32)

229

xz

Hrð Þ¼ r; ϑ; t hrð Þ r; t cos ϑ, Hϑð Þ¼� r; ϑ; t hϑð Þ r; t sin ϑ, Hφð Þ¼ r; ϑ; t h<sup>φ</sup> sin ϑ, hrð Þ r; t

z2

Since the problem is axially symmetric, any plane coming across axis z can be taken to calculate the magnetic field lines. For example, for plane y ¼ 0, the differential equation describing the

> dz <sup>¼</sup> xz hð Þ <sup>r</sup> � <sup>h</sup><sup>ϑ</sup> r<sup>2</sup>h<sup>ϑ</sup> þ z<sup>2</sup>ð Þ hr � h<sup>ϑ</sup>

The magnetic field lines for the reference solution at time t = 0.01 are shown in Figure 9.

Results of Eq. (32) are the formulas for the radial and angular components of the field depending

xHyð Þ� x; y; z; t yHxð Þ x; y; z; t

If function f rð Þ ; t is known, these components are calculated using formulas.

<sup>r</sup><sup>2</sup> ð Þ hrð Þ� <sup>r</sup>; <sup>t</sup> <sup>h</sup>ϑð Þ <sup>r</sup>; <sup>t</sup> , Hxð Þ¼ <sup>x</sup>; <sup>y</sup>; <sup>z</sup>; <sup>t</sup>

dx

hrð Þ¼ r; t Hzð Þþ x; y; z; t

hϑð Þ¼ r; t Hzð Þ� x; y; z; t

hφð Þ¼ r; t

In Cartesian coordinates the field components have the forms

Hzð Þ¼ x; y; z; t hϑð Þþ r; t

Eq. (27) has been taken for the reference solution.

Figure 8. The problem of diffusion into a plasma cloud [16].

¼ zy

slope of the magnetic field lines looks like

on a single space coordinate:

This is an axially symmetric problem, and, therefore, it is convenient to use the polar coordinate system ð Þ r; ϑ;φ . With no Hall effect and with regard to the conditions at infinity, the initial data for vector potential takes the form [16]

$$\Psi\_r(r, \theta, t = 0) = 0, \quad \Psi\_\theta(r, \theta, t = 0) = 0, \quad \Psi\_\phi(r, \theta, t = 0) = r f(r, t = 0) \sin \theta,\tag{25}$$

$$f(r, t = 0) = \frac{H\_0}{2} \begin{cases} 0 & r < r\_0 \\ \left(1 - (r\_0/r)^3\right) & r > r\_0 \end{cases} \tag{26}$$

The solution of Eq. (24) with initial data Eq. (25) is reduced to solving Eq. (27) with initial data Eq. (26) and boundary conditions Eq. (28):

$$\frac{df}{dt} = \frac{\nu(r)}{r^4} \frac{\partial}{\partial r} r^4 \frac{\partial f}{\partial r'} \tag{27}$$

$$\frac{\partial f}{\partial r}(r=0,t) = 0,\\ \frac{\partial f}{\partial r}(r \to \ast, t) = 0. \tag{28}$$

The solution has the form

$$
\Psi\_r(r, \theta, t = 0) = 0, \quad \Psi\_\vartheta(r, \theta, t) = 0, \quad \Psi\_\vartheta(r, \theta, t) = r f(r, t) \sin \theta. \tag{29}
$$

Paper [16] considers the plasma cloud interaction with the magnetic field of vacuum, and, therefore, ν<sup>2</sup> ! ∞ is assumed. For this special case, the solution to Eq. (27) in quadratures has been obtained, and it has the forms:

$$f(r,t) = \frac{H\_0}{2} \begin{cases} 1 - \frac{6}{\zeta^2} \sum\_{n=1}^{\omega} \frac{(-1)^n}{\left(\pi n\right)^2} T\_n(t) \cdot \left(\cos \pi n \zeta - \frac{\sin \pi n \zeta}{\pi n \zeta}\right), & 0 < \zeta = \frac{r}{r\_0} < 1 \\ 1 - \frac{6}{\zeta^2} \sum\_{n=1}^{\omega} \frac{1}{\left(\pi n\right)^2} T\_n(t), & 1 < \zeta \end{cases}, \\ T\_n(t) = \exp\left(-\left(\pi n\right)^2 \frac{\nu\_1 t}{r\_0^2}\right). \tag{30}$$

Figure 8. The problem of diffusion into a plasma cloud [16].

∂H

.

5.1. Diffusion of magnetic field in the absence of the Hall effect

∂Ψ

f rð Þ¼ ; <sup>t</sup> <sup>¼</sup> <sup>0</sup> <sup>H</sup><sup>0</sup>

∂f ∂r 2

df dt <sup>¼</sup> <sup>ν</sup>ð Þ<sup>r</sup> r4 ∂ ∂r r <sup>4</sup> ∂f ∂r

ð Þ¼ r ¼ 0; t 0,

ð Þ <sup>π</sup><sup>n</sup> <sup>2</sup> Tnð Þ� <sup>t</sup> cosπn<sup>ζ</sup> � sinπn<sup>ζ</sup>

(

ν1, r < r<sup>0</sup> ¼ 1 ν2, r > r<sup>0</sup>

<sup>¼</sup> ffiffiffi 2

cloud: νð Þ¼ r

�

228 Numerical Simulations in Engineering and Science

relative to vector potential H ¼ rotΨ:

data for vector potential takes the form [16]

Eq. (26) and boundary conditions Eq. (28):

The solution has the form

f rð Þ¼ ;t

H<sup>0</sup> 2

8 >>>><

>>>>:

been obtained, and it has the forms:

ð Þ �<sup>1</sup> <sup>n</sup>

1

<sup>1</sup> � <sup>6</sup> ζ2 X∞ n¼1

<sup>1</sup> � <sup>6</sup> ζ2 X∞ n¼1 <sup>∂</sup><sup>t</sup> ¼ �rotð Þ <sup>ν</sup> � rot<sup>H</sup> <sup>þ</sup> <sup>b</sup>½ � <sup>H</sup> � rot<sup>H</sup> : (23)

<sup>∂</sup><sup>t</sup> ¼ �<sup>ν</sup> � rot rotΨ: (24)

r > r<sup>0</sup> :

, (27)

, TnðÞ¼ <sup>t</sup> exp �ð Þ <sup>π</sup><sup>n</sup> <sup>2</sup> <sup>ν</sup>1<sup>t</sup>

r2 0 :

(30)

� �

ð Þ¼ r ! ∞; t 0: (28)

(26)

It is assumed that the magnetic field at infinity is uniform and directed along axis z: Hð0; 0; H<sup>0</sup>

Assume that the Hall effect contribution is small, bH0=ν ≪ 1. Write the equation of diffusion

This is an axially symmetric problem, and, therefore, it is convenient to use the polar coordinate system ð Þ r; ϑ;φ . With no Hall effect and with regard to the conditions at infinity, the initial

The solution of Eq. (24) with initial data Eq. (25) is reduced to solving Eq. (27) with initial data

∂f ∂r

Paper [16] considers the plasma cloud interaction with the magnetic field of vacuum, and, therefore, ν<sup>2</sup> ! ∞ is assumed. For this special case, the solution to Eq. (27) in quadratures has

πnζ

� �

ð Þ <sup>π</sup><sup>n</sup> <sup>2</sup> Tnð Þ<sup>t</sup> , <sup>1</sup> <sup>&</sup>lt; <sup>ζ</sup>

Ψrð Þ¼ r; ϑ; t ¼ 0 0, Ψϑð Þ¼ r; ϑ; t ¼ 0 0, Ψφð Þ¼ r; ϑ; t ¼ 0 rf rð Þ ; t ¼ 0 sin ϑ, (25)

<sup>1</sup> � ð Þ <sup>r</sup>0=<sup>r</sup> <sup>3</sup> � �

0 r < r<sup>0</sup>

Ψrð Þ¼ r; ϑ; t ¼ 0 0, Ψϑð Þ¼ r; ϑ; t 0, Ψφð Þ¼ r; ϑ; t rf rð Þ ; t sin ϑ: (29)

, <sup>0</sup> <sup>&</sup>lt; <sup>ζ</sup> <sup>¼</sup> <sup>r</sup>

r0 < 1

<sup>p</sup> <sup>Þ</sup> (see Figure 8). The magnetic viscosity coefficient is constant inside and outside the

For finite values of conductivity in external domain ν<sup>2</sup> > 0, the limit numerical solution of Eq. (27) has been taken for the reference solution.

Components of magnetic field H ¼ rotΨ are found by differentiating vector potential Eq. (29). If function f rð Þ ; t is known, these components are calculated using formulas.

$$\begin{split} H\_r(r, \ $, t) &= h\_r(r, t) \cos \$ , H\_\theta(r, \ $, t) = -h\_\theta(r, t) \sin \$ , H\_\psi(r, \ $, t) = h\_\psi \sin \$ , h\_r(r, t) \\ &= 2f(r, t), h\_\theta(r, t) = \text{\∂}r^2 f(r \text{\'} \text{\'} r, h\_\theta(r, t) = 0) \end{split} \tag{31}$$

In Cartesian coordinates the field components have the forms

$$\begin{aligned} H\_z(\mathbf{x}, y, z, t) &= h\wp(r, t) + \frac{z^2}{r^2} (h\_r(r, t) - h\_\vartheta(r, t)), H\_\vartheta(\mathbf{x}, y, z, t) \\ &= \frac{zy}{r^2} (h\_r(r, t) - h\_\vartheta(r, t)), H\_\mathbf{x}(\mathbf{x}, y, z, t) = \frac{\chi\varpi}{r^2} (h\_r(r, t) - h\_\vartheta(r, t)), \end{aligned} \tag{32}$$

Since the problem is axially symmetric, any plane coming across axis z can be taken to calculate the magnetic field lines. For example, for plane y ¼ 0, the differential equation describing the slope of the magnetic field lines looks like

$$\frac{d\mathbf{x}}{dz} = \frac{\mathbf{x}z(h\_r - h\_\vartheta)}{r^2h\_\vartheta + z^2(h\_r - h\_\vartheta)}.$$

The magnetic field lines for the reference solution at time t = 0.01 are shown in Figure 9.

Results of Eq. (32) are the formulas for the radial and angular components of the field depending on a single space coordinate:

$$\begin{aligned} h\_r(r,t) &= H\_z(\mathbf{x},y,z,t) + \frac{xH\_x(\mathbf{x},y,z,t) + yH\_y(\mathbf{x},y,z,t)}{z}, \\ h\_\vartheta(r,t) &= H\_z(\mathbf{x},y,z,t) - \frac{z\left(xH\_x(\mathbf{x},y,z,t) + yH\_y(\mathbf{x},y,z,t)\right)}{r^2 - z^2}, \\ h\_\vartheta(r,t) &= \frac{xH\_\vartheta(\mathbf{x},y,z,t) - yH\_x(\mathbf{x},y,z,t)}{r} = 0. \end{aligned}$$

Figure 10. Profiles of the non-dimensionalized field components at time t = 0.1: (a) radial hr and (b) angular hθ.

Benchmarks for Non-Ideal Magnetohydrodynamics http://dx.doi.org/10.5772/intechopen.75713 231

Figure 11. Profiles of the field components along line x=z and y=0 at time t = 0.1: (a) Hx and (b) Hz.

Figure 12. Evolution of non-dimensionalized field components in the problem with parameters ν<sup>1</sup> = 1 and ν<sup>2</sup> ! ∞: (a)

radial hr and (b) angular hθ.

Figure 9. The magnetic field lines in the exact solution with parameters ν<sup>1</sup> = 1 and ν<sup>2</sup> ! ∞. A dashed line shows the plasma cloud position.

Figures 10 and 11 show profiles of field components for different magnetic viscosity values in external domain (r > r0). Note that there is a small difference between the profiles obtained with ν<sup>2</sup> = 50 and ν<sup>2</sup> ! ∞ corresponding to the simulation of the plasma cloud interaction with the magnetic field of vacuum ν<sup>2</sup> ! ∞. For the problem with ν<sup>2</sup> ! ∞, profiles of the nondimensionalized field components for the initial phase of diffusion, t<r0 3 /ν1, are given (see Figure 12).

Simulation setup: A computational domain ð Þ j j x ≤ 0:5L yj j ≤ 0:5L; j j z ≤ 0:5L is a cube with edges L = 10. Boundary conditions corresponding to the initial undisturbed state are imposed on its lateral faces for the components of field Hð Þj r; t <sup>r</sup><sup>∈</sup> <sup>Γ</sup> ¼ Hð Þ r; t ¼ 0 . The initial data can be set either for the magnetic field components Eq. (33) or the components of vector potential Eq. (34).

$$\begin{aligned} H\_z(\mathbf{x}, y, z, t = 0) &= 2f(r, t = 0) + z^2 \left( 1 + \eta f'(r, t = 0) \right) / r^2, \\ H\_y(\mathbf{x}, y, z, t = 0) &= \mathbf{x} \mathbf{y} f'(r, t = 0) / r, \quad H\_x(\mathbf{x}, y, z, t = 0) = \mathbf{x} \mathbf{z} f'(r, t = 0) / r, \end{aligned} \tag{33}$$

$$\begin{aligned} \Psi\_z(\mathbf{x}, y, z, t = 0) &= 0, & \Psi\_y(\mathbf{x}, y, z, t = 0) &= -zf(r, t = 0), \\ \Psi\_x(\mathbf{x}, y, z, t = 0) &= yf(r, t = 0). \end{aligned} \tag{34}$$

Figure 10. Profiles of the non-dimensionalized field components at time t = 0.1: (a) radial hr and (b) angular hθ.

Figure 11. Profiles of the field components along line x=z and y=0 at time t = 0.1: (a) Hx and (b) Hz.

Figures 10 and 11 show profiles of field components for different magnetic viscosity values in external domain (r > r0). Note that there is a small difference between the profiles obtained with ν<sup>2</sup> = 50 and ν<sup>2</sup> ! ∞ corresponding to the simulation of the plasma cloud interaction with the magnetic field of vacuum ν<sup>2</sup> ! ∞. For the problem with ν<sup>2</sup> ! ∞, profiles of the non-

Figure 9. The magnetic field lines in the exact solution with parameters ν<sup>1</sup> = 1 and ν<sup>2</sup> ! ∞. A dashed line shows the

Simulation setup: A computational domain ð Þ j j x ≤ 0:5L yj j ≤ 0:5L; j j z ≤ 0:5L is a cube with edges L = 10. Boundary conditions corresponding to the initial undisturbed state are imposed on its lateral faces for the components of field Hð Þj r; t <sup>r</sup><sup>∈</sup> <sup>Γ</sup> ¼ Hð Þ r; t ¼ 0 . The initial data can be set either for the magnetic field components Eq. (33) or the components of vector potential

Ψzð Þ¼ x; y; z; t ¼ 0 0, Ψyð Þ¼� x; y; z; t ¼ 0 zf rð Þ ; t ¼ 0 ,

ð Þ <sup>r</sup>; <sup>t</sup> <sup>¼</sup> <sup>0</sup> <sup>=</sup>r<sup>2</sup>,

<sup>Ψ</sup>xð Þ¼ <sup>x</sup>; <sup>y</sup>; <sup>z</sup>; <sup>t</sup> <sup>¼</sup> <sup>0</sup> yf rð Þ ; <sup>t</sup> <sup>¼</sup> <sup>0</sup> : (34)

ð Þ r; t ¼ 0 =r, Hxð Þ¼ x; y; z; t ¼ 0 xzf <sup>0</sup>

3

ð Þ r; t ¼ 0 =r,

/ν1, are given (see

(33)

dimensionalized field components for the initial phase of diffusion, t<r0

Hzð Þ¼ <sup>x</sup>; <sup>y</sup>; <sup>z</sup>; <sup>t</sup> <sup>¼</sup> <sup>0</sup> <sup>2</sup>f rð Þþ ; <sup>t</sup> <sup>¼</sup> <sup>0</sup> <sup>z</sup><sup>2</sup> <sup>1</sup> <sup>þ</sup> rf <sup>0</sup>

Hyð Þ¼ x; y; z; t ¼ 0 xyf <sup>0</sup>

Figure 12).

plasma cloud position.

230 Numerical Simulations in Engineering and Science

Eq. (34).

Figure 12. Evolution of non-dimensionalized field components in the problem with parameters ν<sup>1</sup> = 1 and ν<sup>2</sup> ! ∞: (a) radial hr and (b) angular hθ.

The EGIDA code uses a scheme preserving the field divergence at one step because difference operators DIV and ROT (div and rot) [18] determined at nodes and in cells of a grid, respectively, satisfy the vector analysis identities: DIVrot = 0 and ROTdiv = 0.

It has been found that for the first set of initial data the divergence norm depends on errors induced by the initial distribution of the H field components in the vicinity of sphere r=r0. Though the difference scheme does not change the magnetic field divergence, initial errors lead to a significantly distorted numerical solution (see Figure 13).

In the second case, the magnetic field components are determined using the operator numerically differentiating the vector potential, and, hence, the magnetic field divergence norm equals zero at initial time and at all later times. For this case, a good agreement between the calculated results and the exact solution has been achieved even on the coarsest grid (see Figure 14).

5.2. Diffusion of magnetic field with a low Hall effect

∂Ψ

plasma [16]:

Hy in section z = 0 (c).

required to take into account the Hall term in the diffusion equation:

Assume that the Hall effect contribution is small, bH0=ν ≪ 1, but finite. For this reason, it is

Figure 14. Calculation of the magnetic field diffusion into a spherical plasma cloud for the second set of initial data (34). Profiles of the field components at time t = 0.1: (a) radial hr and (b) angular hθ. Distribution of transverse field component

Benchmarks for Non-Ideal Magnetohydrodynamics http://dx.doi.org/10.5772/intechopen.75713 233

<sup>∂</sup><sup>t</sup> ¼ �<sup>ν</sup> � rot rot<sup>Ψ</sup> � b rot ½ � <sup>Ψ</sup> � rot rot<sup>Ψ</sup> :

The Hall effect leads to the occurrence of the azimuthal component, Hφ, of magnetic field in

Hrð Þ¼ r; ϑ; t hrð Þ r; t cos ϑ, Hϑð Þ¼� r; ϑ; t hϑð Þ r; t sin ϑ, Hφð Þ r; ϑ; t <sup>¼</sup> <sup>h</sup><sup>φ</sup> sin <sup>ϑ</sup>, hrð Þ¼ <sup>r</sup>; <sup>t</sup> <sup>2</sup>f rð Þ ; <sup>t</sup> , hϑð Þ¼ <sup>r</sup>; <sup>t</sup> <sup>∂</sup>r<sup>2</sup><sup>f</sup> <sup>=</sup>r∂r, hφð Þ¼ <sup>r</sup>; <sup>t</sup> <sup>Ψ</sup>ð Þ <sup>r</sup>; <sup>t</sup> <sup>=</sup>r:

Figure 13. Calculation of the magnetic field diffusion into a spherical plasma cloud for the first set of initial data (33). Profiles of the field components at time t = 0.1: (a) radial hr and (b) angular hθ. Distribution of transverse field component Hy in section z = 0 (c).

Figure 14. Calculation of the magnetic field diffusion into a spherical plasma cloud for the second set of initial data (34). Profiles of the field components at time t = 0.1: (a) radial hr and (b) angular hθ. Distribution of transverse field component Hy in section z = 0 (c).

### 5.2. Diffusion of magnetic field with a low Hall effect

The EGIDA code uses a scheme preserving the field divergence at one step because difference operators DIV and ROT (div and rot) [18] determined at nodes and in cells of a grid, respec-

It has been found that for the first set of initial data the divergence norm depends on errors induced by the initial distribution of the H field components in the vicinity of sphere r=r0. Though the difference scheme does not change the magnetic field divergence, initial errors

In the second case, the magnetic field components are determined using the operator numerically differentiating the vector potential, and, hence, the magnetic field divergence norm equals zero at initial time and at all later times. For this case, a good agreement between the calculated results and the exact solution has been achieved even on the coarsest grid (see

Figure 13. Calculation of the magnetic field diffusion into a spherical plasma cloud for the first set of initial data (33). Profiles of the field components at time t = 0.1: (a) radial hr and (b) angular hθ. Distribution of transverse field component

tively, satisfy the vector analysis identities: DIVrot = 0 and ROTdiv = 0.

232 Numerical Simulations in Engineering and Science

lead to a significantly distorted numerical solution (see Figure 13).

Figure 14).

Hy in section z = 0 (c).

Assume that the Hall effect contribution is small, bH0=ν ≪ 1, but finite. For this reason, it is required to take into account the Hall term in the diffusion equation:

$$\frac{\partial \Psi}{\partial t} = -\nu \cdot rot \; rot \Psi - b[rot \Psi \times rot \; rot \; \Psi].$$

The Hall effect leads to the occurrence of the azimuthal component, Hφ, of magnetic field in plasma [16]:

$$\begin{aligned} H\_r(r, \mathfrak{d}, t) &= h\_r(r, t) \cos \mathfrak{d}, H\_{\mathfrak{d}}(r, \mathfrak{d}, t) = -h\_{\mathfrak{d}}(r, t) \sin \mathfrak{d}, H\_{\mathfrak{q}}(r, \mathfrak{d}, t) \\ &= h\_{\mathfrak{q}} \sin \mathfrak{d}, h\_r(r, t) = 2f(r, t), h\_{\mathfrak{d}}(r, t) = \partial r^2 f(r \partial r, h\_{\mathfrak{q}}(r, t) = \Psi(r, t)/r. \end{aligned}$$

6. Conclusion

Acknowledgements

Author details

VNIIEF, Sarov, Russia

References

some benchmarks and references.

Sofronov Vasily, Zhmailo Vadim and Yanilkin Yury\* \*Address all correspondence to: n.yanilkina@mail.ru

Moscow: Binom; 2009 (in Russian)

tional Physics. 1999;149:270-292

of Computational Physics. 2000;161:605-652

An important property of difference schemes in multidimensional flow simulations is that they keep the magnetic field divergence-free in difference solutions. An adverse aspect of this defect

Benchmarks for Non-Ideal Magnetohydrodynamics http://dx.doi.org/10.5772/intechopen.75713 235

Note that the zero-divergence requirements to difference schemes get more stringent as applied to the solution of diffusion problems. A violation of this requirement results in the accumulation of errors and loss of solution structure, especially in problems with high conductivity gradients.

Some EGIDA calculations have been performed under contract no. 1239349 between Sandia National Laboratories and RFNC-VNIIEF. The authors thank sincerely M. Pokoleva for her assistance in simulations and J. Kamm and A. Robinson for their assistance in formulating

[1] Kulikovsky AG, Pogorelov NV, Semenov AY. Mathematical Issues of Numerical Solution

[2] Brushlinsky KV. Mathematical and Computational Problems of Magnetic Gas Dynamics.

[3] Balsara DS, Spicer DS. A staggered mesh algorithm using high Godunov fluxes to ensure solenoidal magnetic fields in magneto hydrodynamics simulations. Journal of Computa-

[4] Toth G. The ∇B=0 constraint in shock–capturing magneto hydrodynamics codes. Journal

[5] Gardiner TA, Stone JM. An unsplit Godunov method for ideal MHD via constrained

of Hyperbolic Equations. Moscow: Fizmatlit; 2001 (In Russian)

transport. Journal of Computational Physics. 2005;205:509-539

is the unphysical transport of matter orthogonal to the field H [2].

Figure 15. Profiles of the azimuthal component of non-dimensionalized field in the problem with parameters ν<sup>1</sup> ¼ 1, ν<sup>2</sup> ! ∞, b ¼ 0:01.

In view of the smallness of parameter bH0=ν1, we have Eqs. (27), (28) to calculate function f rð Þ ; t , while the calculation of small additive Ψ is reduced to solving the following boundary value problem:

$$\frac{d\Psi}{dt} = \frac{\partial}{\partial r}\nu \frac{\partial \Psi}{\partial r} - \frac{6\nu\Psi}{r^2} - 2f r^2 \frac{\partial}{\partial r} \left(\frac{b}{r^4} \frac{\partial}{\partial r} r^4 \frac{\partial f}{\partial r}\right), \\ \Psi(r, t = 0) = 0, \\ \Psi(r = 0, t) = 0, \\ \Psi(r \to \ast, t) = 0.$$

Figure 15 shows profiles of the non-dimensionalized azimuthal field component h<sup>φ</sup> at early times, t < r<sup>2</sup> <sup>0</sup>=ν1. With small values of parameter bH0=ν1, the rest two components—hr, hϑ remain unchanged, and they are shown in Figure 12. Note that, as it has been shown in [16], if the motion of plasma is accounted, the Hall effect may lead to the occurrence of azimuthal velocity, i.e., to the plasma cloud rotation.

The changeover to Cartesian coordinates is performed using formulas

$$\begin{aligned} H\_z(\mathbf{x}, y, z, t) &= h\_{\ $}(r, t) + z^2 (h\_r(r, t) - h\_{\$ }(r, t)) / r^2, \\ H\_y(\mathbf{x}, y, z, t) &= z y (h\_r(r, t) - h\_{\ $}(r, t)) / r^2 + x h\_{\$ }(r, t) / r, \\ H\_x(\mathbf{x}, y, z, t) &= \mathbf{x} z (h\_r(r, t) - h\_{\ $}(r, t)) / r^2 - y h\_{\$ }(r, t) / r. \end{aligned}$$

Simulation setup: The initial data is given in the previous section. This problem requires accounting the Hall effect. A local exchange parameter is b = 0.01. Figure 15 shows profiles of the azimuthal component of non-dimensionalized field.
