Using USim to Solve a Magnetic Nozzle Problem

In this tutorial we show how to use USim to solve a problem with both an induced and imposed magnetic field using the gas dynamic form of the mhd equations.

Solving Problems Using The Gas Dynamic Form of the MHD Equations

This tutorial follows many of the same concepts as the previous tutorial Using USim to Solve MHD with General Equation of State. In this case we are solving the MHD system written in gas dynamic form. This form is significant since it is not conservative, but simplifies some considerations when there are both imposed and induces magnetic fields.

First of all we want to split the field into an imposed field generated by external coils and and induced field generated by the plasma motion. To do this we have 3 data structures, the first, backgroundB stores the field generated by a magnetic field coil

<DataStruct backgroundB>
   kind = nodalArray
   onGrid = domain
   numComponents = 3
 </DataStruct>

The second variable stores the conserved variables along with the total magnetic field, i.e., the induced field + the background field

<DataStruct qModified>
  kind = nodalArray
  onGrid = domain
  numComponents = NUMCOMP
</DataStruct>

The 3rd variable stores the conserved variables and only the induced magnetic field

<DataStruct q>
  kind = nodalArray
  onGrid = domain
  numComponents = NUMCOMP
</DataStruct>

The imposed magnetic field is generated at startup using a wire source wireFieldEqn and a Hyperbolic Equations updater. The Hyperbolic Equations updater works with all USim Algebraic Equations blocks. Since the problem is not axisymmetric we use a wire, but could use a coilFieldEqn to model a current carrying loop. The wireFieldEqn initialization of the background magnetic field is given below

<Updater initMagField>
  kind = equation2d
  onGrid = domain
  in = []
  out = [backgroundB]

  <Equation coil>
    kind = wireFieldEqn

    outRange = [0, 1, 2]

    point = [$COILRAD$, -0.001, 0.0]
    mu0 = MU0
    normal = [0.0, 0.0, 1.0]
    current = $COILCURRENT$
  </Equation>
</Updater>

Now the fields are split into induces q and imposed fields backgroundB. The imposed field backgroundB does not evolve in time, but does affect the flow of fluid through the \(J\times B\) force which is treated as a source in this example. Taking this into account we apply the hyperbolic update to q (which ignores the imposed field) not qModified which includes the imposed field. The classicMusclUpdater (1d, 2d, 3d) updater is used with the input file seen below

<Updater hyper>
  kind = classicMuscl2d
  timeIntegrationScheme = none
  gradientType = leastSquares

  numericalFlux = hlleFlux
  variableForm = conservative
  preservePositivity = true

  limiter = [muscl, none, muscl, none]

  onGrid = domain
  in = [q, J, E, Z]
  out = [qnew]

  cfl = CFL

  equations = [mhd]

  <Equation mhd>
    kind = gasDynamicMhdEqn
    gasGamma = ION_GAMMA
    mu0 = MU0
    correctionSpeed = 0.0
    ionMass = ION_MASS
    chargeState = ZRATIO
    fundamentalCharge = FUNDAMENTAL_CHARGE
  </Equation>
</Updater>

Since this equation system does not use the general equation of state the inputs are somewhat simpler than in the case of the twoTemperatureMhdEosEqn. The inputs are the conserved variables q, the total current J, the electric field E and the charge state Z. Also, note that we have not included the mhdSrc which was included in the case of twoTemperatureMhdEosEqn. Instead we add the mhdSrc in a separate step. First we compute qModified which includes the conserved variables and both the induced and imposed field (we use a combiner (1d, 2d, 3d))

<Updater computeQMod>
   kind = combiner2d
   onGrid = domain

   in = [q, backgroundB]
   out = [qModified]

   indVars_q = ["rho", "mx", "my", "mz", "en", "bx", "by", "bz", "phi", "ee"]
   indVars_backgroundB = ["B0x","B0y","B0z"]

   exprs = ["rho", "mx", "my", "mz", "en", "bx+B0x", "by+B0y", "bz+B0z", "phi", "ee"]
 </Updater>

The current density is computed from the perturbed magnetic field only \(J=\frac{1}{\mu_{0}}\nabla \times B\) since it is already known that that the curl of the imposed field is zero everywhere inside the domain

<Updater computeJ>
    kind = vector2d
   onGrid = domain
   derivative = curl
   numScalars = 1
   orderAccuracy = 2
   coefficient = 1.0
   numberOfInterpolationPoints = 8

   in = [B]
   out = [J]
 </Updater>

Now that we’ve computed qModified we use this to compute the mhdSrc through the use of a Hyperbolic Equations so that both the induced and imposed fields provide a force to the fluids. The source term is stored in src.

<Updater computeSource>
   kind = equation2d
   in = [qModified, J, E, Z]
   out = [src]
   onGrid = domain

   <Equation gdSrc>
      kind = mhdSrc
      model = gasDynamicMhdEqn
      gasGamma = ION_GAMMA
      electronGasGamma = ELECTRON_GAMMA
      ionMass = ION_MASS
      chargeState = ZRATIO
      fundamentalCharge = FUNDAMENTAL_CHARGE
      correctionSpeed = 0.0
      mu0 = MU0
   </Equation>
 </Updater>

The source term is then added to the update q, qnew through the use of a uniformCombiner (1d, 2d, 3d). The uniformCombiner (1d, 2d, 3d) acts just like a combiner (1d, 2d, 3d) except that it assumes all input variables and output variables are of the same size and that the same operation is being applied to all components. In the block below we have added src generated by the computeSource Updater above

<Updater updateQ>
  kind = uniformCombiner2d
  onGrid = domain

  in = [qnew, src]
  out = [qnew]

  indVars_qnew = ["qn"]
  indVars_src = ["s"]
  exprs = ["qn+s"]
</Updater>

As in the previous tutorial the electric field is computed using a generalizedOhmsLaw (1d, 2d, 3d) updater. Notice that the conserved variables including the total magnetic field qModified is used in computing the electric field. In addition, we’ve included a resistivity term by defining resistivity=eta. eta is a nodalArray defined using a combiner (1d, 2d, 3d)

<Updater computeE>
  kind = generalizedOhmsLaw2d
  onGrid = domain

  in = [qModified, J, Z]
  out = [E]
  resistivity = eta

  hallTerm = false

  fundamentalCharge = FUNDAMENTAL_CHARGE
  ionMass = ION_MASS
  electronMass = ME
  boltzmannConstant = KB

  mu0 = MU0
</Updater>

Computing the \(\Delta t\) restriction for explicit resistive term

USim has a Time Step Restriction for explicit diffusion type terms. The restriction can be applied to viscous, thermal diffusion and resistive terms. In order to compute the restriction properly the diffusion coefficient needs to be computed properly. In the case of resistivity we want to compute the diffusion coefficient \(\gamma\) in the proper form

\[\frac{\partial B}{\partial t}=\gamma\nabla^{2}B\]

It turns out gamma=eta/mu_{0} in the case of magnetic field diffusion. In this simulation we first compute the resistivity

<Updater initEta>
   kind = initialize2d
   onGrid = domain
   out = [eta]

   <Function func>
     kind = exprFunc
     eta0 = RESISTIVITY
     exprs = ["eta0"]
   </Function>
 </Updater>

and then compute etaBymu0 = eta0/mu0 so that we can compute the explicit time step constraint

<Updater initEtaBymu0>
   kind = combiner2d
   onGrid = domain
   in = [eta]
   out = [etaBymu0]
   indVars_eta = ["eta0"]
   mu0 = $MU0$
   exprs = ["eta0/mu0"]
 </Updater>

etaByMu0 can then be used in the quadratic (1d, 2d, 3d) time step restriction updater

<Updater timeStepRestriction>
  kind = timeStepRestrictionUpdater2d
  in = [etaBymu0]
  onGrid = domain
  restrictions = [quadratic]

  <TimeStepRestriction quadratic>
    kind = quadratic
    cfl = 0.5
  </TimeStepRestriction>
</Updater>

Note that diffusion terms have explicit time step restriction given by \(\Delta t \approx \Delta x^{2}\) so that the time step reduces quadatically as the grid is refined. Currently we can get around this problem in USim by using super time stepping (see Advanced Time-Stepping Methods in USim).

Divergence Cleaning with an Imposed Field

Divergence cleaning with imposed fields is accomplished by cleaning the induced field only. It is already known that the imposed field is divergence free so that part can be ignored. Hyperbolic divergence cleaning in this case, but only applied to the perturbed field. A DataStructAlias is used to point just to the magnetic field and correction potential components in the conserved variable vector

<DataStructAlias qClean>
   kind = nodalArray
   target = q
   componentRange = [5,9]
   writeOut = 0
 </DataStructAlias>

This DataStructAlias is used in the hyperbolic cleaning step which consists of a classicMusclUpdater (1d, 2d, 3d) updater, a multiUpdater (1d, 2d, 3d) and associated boundary conditions.

<Updater hyperClean>
  kind = classicMuscl2d
  timeIntegrationScheme = none
  gradientType = leastSquares
  variableForm = conservative

  cfl = CFL

  numericalFlux = hlleFlux

  limiter = [muscl]

  onGrid = domain
  in = [qClean]
  out = [qCleanNew]

  equations = [clean]

  syncAfterSubStep = [qCleanNew]
  <Equation clean>
    kind = hyperbolicCleanEqn
    waveSpeed = CORRECTION_SPEED
  </Equation>
</Updater>

<Updater hyperClean2>
  kind = multiUpdater2d
  onGrid = domain
  timeIntegrationScheme = rk2
  updaters = [cleanCopyBc, cleanInflow, hyperClean]

  integrationVariablesIn = [qClean]
  integrationVariablesOut = [qCleanNew]
  dummyVariables_qClean = [dummyClean1,dummyClean2]

  syncAfterSubStep = [qCleanNew]
</Updater>

In this case the boundary conditions consist of copying out the magnetic field and reversing the sign of the correction potential. In can be the case that field build up can occur at the boundary and in these cases simply setting the perturbed magnetic field to zero resolves the issue. Boundary conditions used for the hyperbolic cleaning step with the perturbed magnetic field are shown below using a generalBc (1d, 2d, 3d)

<Updater cleanCopyBc>
  kind = generalBc2d
  onGrid = domain

  in = [qClean]
  dynVectors = []

  indVars_qClean = ["Bx","By","Bz","Phi"]

  exprs = ["Bx","By","Bz","-Phi"]
  out = [qClean]

  entity = ghost
</Updater>

<Updater cleanInflow>
  kind = generalBc2d
  onGrid = domain

  in = [qClean]
  dynVectors = []

  indVars_qClean = ["Bx","By","Bz","Phi"]

  exprs = ["Bx","By","Bz","-Phi"]
  out = [qClean]

  entity = ghost
</Updater>

An Example Simulation

The input file for the problem Magnetic Nozzle in the USimHEDP package demonstrates each of the concepts described above.