Axisymmetric simulations - vonkarmaninstitute/pantera-pic-dsmc GitHub Wiki

Positions and velocities in Pantera are always stored in 3D cartesian coordinates. The only difference between 1D, 2D, and 3D simulations are the number of cells and boundary conditions in the additional dimensions. This means that for a 2D simulation, 1 cell will be used in the z-direction, and the zmin and zmax boundaries would be typically set to periodic. For 1D simulations, the same applies to the y-direction as well. 3D simulations are very expensive, since number of particles $\propto$ number of cells $\propto(D/\lambda)^3$, so when possible we resort to 2D axisymmetric simulations. This is made possible in PIC-DSMC with clever modifications to a series of procedures.

Fundamental idea

Axisymmetric_simulations_cylinder_with_particles.svg

The idea for axisymmetric simulations is to keep storing the velocities and positions in cartesian coordinates, such that the cell sorting and collision procedures remain valid. In Pantera then, we assume $x$ as the symmetry axis, $y$ as the radial coordinate (this means that nothing should be defined such that $y<0$), and $z$ as the azimuthal coordinate, although internally measured in linear rather than angular units. Now we have to be cautious because the following procedures require a special treatment:

  • Injection of particles from [faces of the simulation box] or from [line sources]
  • Movement of the particles and collision checking with [faces of the simulation box] or with [internal walls]
  • Solution of Poisson's equation
  • Weighting of particles in the domain
  • Averaging procedure in the domain

Particle injection

Consider filling the "pie slice" domain with particles at uniform density. You will realize that they are not uniformly distributed in the radial direction. In fact, since the cross-sectional area is proportional to $r$, the probability of finding a particle with a specific value of $r$ is also proportional to $r$ [^boyd] . Therefore, when we have a source for the injection of particles at a certain uniform density, we have to take this in consideration. Rather than sampling the position uniformly in $r$, we use the following expression:

$r = \sqrt{y_1^2 + R(y_2^2-y_1^2)}$

Where:

  • $R$ is a random fraction uniformly distributed in $[0,1]$
  • $y_1, y_2$, with $y_1<y_2$, are the radial coordinates between which we want the particle to be created

If the line source from which we want to inject is not vertical, in this case representing a conical section, me must also generate the $x$ coordinate for the particle. We do this after the radial position has been generated, by computing

$x = x_1 + \frac{x_2-x_1}{y_2 - y_1}(r - y_1)$

Notice that we have a problem if $y_2=y_1$, but in this case the line is horizontal, and there is no need to generate a radial position. So in this case the algorithm generates directly the axial coordinate.

Particle movement

Checking collisions

In Pantera, motion of the particles during a time step is rectilinear. While in cartesian simulations it is quite simple to check collisions with the internal walls and with the simulation box, this is more difficult in axisymmetric simulations, since we have to compute the intersection between the rectilinear trajectory of the particles and cone or cylinder sections. If a wall is characterized by the extremes $(x_1,y_1)$, $(x_2,y_2)$, it can be parametrized by

$\begin{cases} x = x_1 + s(x_2-x_1)\\ y = y_1 + s(y_2-y_1) \end{cases}$

Where $s \in [0,1]$. The particle trajectory can be expressed as

$\begin{cases} x = x_0 + v_x \Delta t\\ y = \sqrt{(y_0+v_y \Delta t)^2 + (v_z \Delta t)^2} \end{cases}$

Where $(x_0,y_0,z_0)$ is the position of the particle at the beginning of the time step, and $(v_x, v_y, v_z)$ is the particle velocity. We suppose that the $z_0=0$, a condition that will be enforced by the algorithm. The possible solutions can be found by solving the resulting second order equation for the time-to-collide $\Delta t$:

$(v_y^2+v_z^2-\alpha^2 v_x^2)(\Delta t)^2 + 2(y_0v_y-(y_1+\beta)\alpha v_x) \Delta t + y_0^2-(y_1+\beta)^2=0$

Where $\alpha = \dfrac{y_2-y_1}{x_2-x_1}$ and $\beta = (x_0-x_1)\alpha$.

Notice that solutions also have to satisfy the inequality:

$\beta+\alpha(v_x \Delta_t)+y_1 \geq 0$

Of course, since we have many possible solutions, the smallest positive solution still smaller than the time step should be accepted. If a solution is accepted, further checks are performed:

  • We check if the collision point is between the extremes of the wall segment, i.e. if $ 0 \leq s_\text{coll} \leq 1$.
  • We check if the velocity at collision time is directed towards the wall, i.e. if $\mathbf{v}_\text{coll} \cdot \hat{\mathbf{n}} < 0$

If both conditions are satisfied, we proceed to the next part of the algorithm: actually moving the particle.

Moving the particle

Once the time by which the particle should be moved has been determined (either the time to a collision or the full time step), we perform the movement. The final position will be

$\begin{cases} x = x_0 + v_x \Delta t\\ y = \sqrt{(y_0+v_y \Delta t)^2 + (v_z \Delta t)^2}\\ z = 0 \end{cases}$

Notice that, as anticipated, we bring back the particle to the initial plane. This requires rotating the velocity vector in the $y-z$ plane. We compute the rotation angle $\theta$ as

$\begin{cases} \sin \theta = \dfrac{v_z \Delta t}{y}\\ \cos \theta = \text{sign}(y_0 + v_y \Delta t) \sqrt{1-\sin^2 \theta} \end{cases}$

Where $\text{sign}(\cdot)$ is the signum function.

Then we perform the rotation, according to the classical rotation matrix

$\begin{bmatrix} v_z'\\ v_y' \end{bmatrix}= \begin{bmatrix} \cos\theta & -\sin\theta\\ \sin\theta & \cos\theta \end{bmatrix} \begin{bmatrix} v_z\\ v_y \end{bmatrix}$

After which the particle is ready for the next collision iteration or time step.

Poisson's equation in cylindrical coordinates

The discretization of Poisson's equation must be adapted to the now cylindrical geometry.

Sampling procedure

Not much has to be changed regarding the sampling procedure, but the code must be aware that the cells now have different volumes even for an uniform grid. Therefore the array CELL_VOLUMES must be correctly initialized and then used during the computation of collisions and during sampling of macroscopic properties.

References

[^boyd]: I. D. Boyd, T. E. Schwartzentruber, Nonequilibrium Gas Dynamics and Molecular Simulation, Cambridge University Press, 2017