Using USim to Solve Multi-Fluid Problems with Collisions¶

Required DataStructs¶

In this example we have 3 different fluids. The first thing we need to define is a variable to contain the collision matrix. Since there are 3 fluids the collision matrix will have 3X3=9 components

<DataStruct collisionMatrix>
  kind = nodalArray
  onGrid = domain
  numComponents = 9
  writeOut = 1
</DataStruct>

So that we can look at the total density, energy and momentum of the system we define qTotal for convenience (this variable is not necessary for this simulation)

<DataStruct qTotal>
  kind = nodalArray
  onGrid = domain
  numComponents = 5
</DataStruct>

For each of the 3 fluids we define a vector to contain the mass density, momentum density and energy density. The first fluid DataStruct is defined as q1 along with its updated value qnew1

<DataStruct q1>
  kind = nodalArray
  onGrid = domain
  numComponents = 5
</DataStruct>

<DataStruct qnew1>
  kind = nodalArray
  onGrid = domain
  numComponents = 5
  writeOut = false
</DataStruct>

Also, we need a variable to hold the momentum and energy collisional source terms for each of the fluids. The length of the vector is 5 to match the number of conserved variables even though the collisional source for the first variable (density) is 0

 <DataStruct q1Mom>
  kind = nodalArray
  onGrid = domain
  numComponents = 5
  writeOut = false
</DataStruct>

In addition data for storing collisions due to thermal energy exchange for each species are also created. The length of this DataStruct is 1

<DataStruct tempRelax1>
  kind = nodalArray
  onGrid = domain
  numComponents = 1
  writeOut = false
</DataStruct>

Data for the species number density is required for each species

<DataStruct N1>
  kind = nodalArray
  onGrid = domain
  numComponents = 1
  writeOut = 1
</DataStruct>

Data for the temperature of each species is also required

<DataStruct T1>
  kind = nodalArray
  onGrid = domain
  numComponents = 1
</DataStruct>

Data is also required for the velocity of each species

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

Computing the collision tensor using Collision Frequency¶

In general the collisionFrequency could be filled up manually using a combiner (1d, 2d, 3d), but it is much easier to use collisionFrequency. The in variable takes temperature number density and velocity for each species (repeated as below for any number of species). In the case where plasma collision are modeled type=ionized charge states Z are also required. In the example below type=neutrals is used so collision cross sections are computed from species diameters given by speciesDia. Notice that inverse=false which means collision frequencies are returned in collisionMatrix. If inverse=true then collision times are returned

<Updater collisionFrequency>
  kind = equation1d
  onGrid = domain

  in =  [T1, N1, V1, T2, N2, V2, T3, N3, V3]

  out = [collisionMatrix]

  <Equation momentum>
    kind = collisionFrequency
    inverse = false
    type = neutrals
    speciesMass = [MI1, MI2, MI3]
    speciesDia = [DI1, DI2, DI3]
  </Equation>
</Updater>

Computing the momentum and energy exchange terms using Momentum Energy Exchange¶

Once collisionMatrix is computed it can be used to compute the momentum and energy exchange term can be computed. Note that this term ignores thermal collisions which are computed using a different operator. momentumEnergyExchange takes in the conserved variables for each of the fluids along with the collisionMatrix and produces momentum and energy exchange terms that can then be added to the fluid equations

<Updater momentumSource>
  kind = equation1d
  onGrid = domain
  in =  [q1, q2, q3, collisionMatrix]
  out = [q1Mom, q2Mom, q3Mom]

  <Equation momentum>
    kind = momentumEnergyExchange
    speciesMass = [MI1, MI2, MI3]
  </Equation>
</Updater>

Computing the temperature relaxation terms using Temperature Relaxation¶

collisionMatrix is also used to compute the temperature relaxation between fluids. temperatureRelaxation takes mass densities or number densities (depending on whether isNumberDensity=1 (takes number densities) or isNumberDensity=0 (takes mass densities). q1,`q2` and q3 are conserved variable vectors where the first component is mass density, the remaining 4 components of the vector are ignored by temperatureRelaxation. The resulting relaxation terms are stored in the output variables and then are added to energy term of the respective fluid equations to compute the effect of temperature relaxation

<Updater energySource>
  kind = equation1d
  onGrid = domain
  in =  [q1, q2, q3, T1, T2, T3, collisionMatrix]
  out = [tempRelax1, tempRelax2, tempRelax3]

  <Equation momentum>
    kind = temperatureRelaxation
    isNumberDensity = 0
    speciesMass = [MI1, MI2, MI3]
  </Equation>
</Updater>

Adding the momentum/energy sources to the fluid equations¶

Now that the momentum and energy exchange terms have been computed along with the temperature relaxation term have been computed these are added to the solution after the hyperbolic part is computed. The source terms are added using a combiner. Notice that the energy term contains contributions from both the momentum exchange and the thermal relaxation term. The updater for the first fluid is given below, the updaters for the remaining fluids will be similar

<Updater addThermalRelaxation1>
  kind = combiner1d
  onGrid = domain

  in = [qnew1, tempRelax1, q1Mom]
  out = [qnew1]

  indVars_qnew1 = ["rho","mx","my","mz","en"]
  indVars_q1Mom = ["dRho","dMx","dMy","dMz","dEn"]
  indVars_tempRelax1 = ["dT"]

  exprs = ["rho","mx+dMx","my+dMy", "mz+dMz","en+dT+dEn"]

</Updater>

Computing the time step for collisions using Frequency¶

Collisions add new time scales that need to be applied when an explicit approach is used. A special updater is used to compute the smallest time scale introduced by collisionMatrix. We use Time Step Restriction which takes in collisionMatrix and then used the restriction kind given by kind=frequency. The frequency restriction takes components which tells USim how many of the components in collisionMatrix should be used. In this case all components are used. This restriction will ensure that the explicit solution is stable to the collisions

<Updater timestepRestriction>
  kind = timeStepRestrictionUpdater1d
  in = [collisionMatrix]
  onGrid = domain
  restrictions = [inverseTime]

  <TimeStepRestriction inverseTime>
    kind = frequency
    components = 9
    cfl = 0.5
  </TimeStepRestriction>
</Updater>

An Example Simulation¶

The input file for the problem Multi-Fluids with collisions in the USimHEDP package demonstrates each of the concepts described above.