Plasma sphere interaction - vonkarmaninstitute/pantera-pic-dsmc GitHub Wiki

This tutorial case is intended to show the user how to setup a plasma simulation using hybrid PIC approach, i.e. simulating electrons as equilibrium fluid through Boltzmann relation. All the setup files of this test case can be found in the tutorials folder 02_charged_drag_on_sphere. For more tutorials, visit the tutorials page.

Introduction

This example is inspired by the paper of C. Capon et al., 2019 [1] that did a numerical reconstruction of original experiment on charged drag by E. D. Knechtel and W. C. Pitts, 1964 [2] where a conductive sphere is exposed to an energetic plasma beam. Charged drag comes from the interaction of charged particles and charged objects, sphere in this case, through electrostatic fields. The aim of this example is to reproduce the same result as in the paper by building the Pantera simulation step-by-step. This setup consists of hybrid Particle-in-Cell method, i.e. simulating electrons as fluid and other species (mainly ions) as particles while solving for the Poisson equation. This also includes the computation of the drag forces (either from direct collisions or from the electrostatic interaction) on any solid body using boundary analysis through separate script.

Physical considerations

The accelerated beam consists of energetic mercury ions and electrons. The paper from the original experiment [2] states following parameters of the beam

Boundary (beam) conditions
$n_\infty$ 1014 m-3
$M_\infty$ 20
$U_\infty$ 10 km/s
$T_\infty$ 3000 K

Within a chamber, they placed a conductive sphere of a given radius which was electrically biased with respect to the plasma. In our case, we select one value of the possible voltages the sphere acquires accordingly

Sphere conditions
$r_{1/2}$ 9.5 / 12.5 mm
$T_\mathrm{s}$ 300 K
$\phi_\mathrm{s}$ - 100 V

The electron part of the beam is assumed to have following parameters

Electron parameters
$n_\mathrm{e}$ 1014 m-3
$T_\mathrm{e}$ 32810 K
$\phi_0$ 0 V

where $\phi_0$ is the reference potential, which is assumed to be 0 V for undisturbed plasma in infinity.

In this case, to assure our simulation is stable, we have to consider all the numerical constraints. Due to the fast movements and oscillations of electrons (due to their much lower mass then heavy ions), one has to resolve spatially the electron Debye length

$ \large \lambda_\mathrm{D} = \left(\frac{\varepsilon_0 k T_\mathrm{e}}{n_\mathrm{e} e^2}\right)^{\frac{1}{2}} \approx 1.2 \ \mathrm{mm}$

The stability condition requires that the size of the cell within domain should be generally $\Delta x < \large \lambda_{\mathrm{D}}$, but in practice, $\Delta x = \large \lambda_{\mathrm{D}} / 4$ is usually preferred.

Moreover, the time stepping needs to be fine enough to resolve the electron plasma oscillations given by

$ \large \omega_p = \left( \frac{n_\mathrm{e} e^2}{\varepsilon_0 m_\mathrm{e}} \right)^{\frac{1}{2}} \quad \rightarrow \quad \omega_p^{-1} \approx 1.8 \ \mathrm{ns}$

It is usually recommended to introduce temporal stability condition $\omega_p \Delta t = 0.2$. However, these constraints for both cases, when simulating plasma dynamics and interactions, are quite heavy for numerical simulations and might require more computational resources. One of the ways how to overcome these constraints is to model electrons as equilibrium fluid through Boltzmann relation as follows

$n_{\mathrm{e}} \left(\phi\right) = n_0 \exp\left(\frac{e (\phi-\phi_0)}{k T_0} \right)$

where $n_0 = 10^{14} \ \mathrm{m}^3$ is the reference density, $T_0 = 32810 \ \mathrm{K}$ is the reference temperature and $\phi_0 = 0 \ \mathrm{V}$ is the reference potential in infinity. This relation essentially assumes thermal equilibrium within the electric field for which is solved self-consistently through the Poisson equation.

Finally, we can drop most of the constraints tied to the electron motion and focus mostly on the energetic ion motion, although the Debye length still needs to be resolved in the vicinity of the sphere that will interact with the plasma flow. Later it will be shown in the geometric files, that we maintain a fine mesh only in the vicinity of the sphere while towards boundaries, we are allowed to mesh the domain with larger size cells. To sum it up, we come to the main constraints that need to be satisfied for the simulation of plasma-sphere interaction

Constraints
$\Delta x_\text{min} = 0.25 \lambda_\mathrm{D}$ 0.3 mm
$\Delta x_\text{max}$ 5 mm
$\Delta t U_{\infty} / x_\text{min}$ $\Delta t =$ 10 ns

Pantera input files

Pantera needs a main input file that contains the setup for the simulation, including numerical parameters, emitting surfaces and species properties and it also requires a .su2 file with generated mesh for the simulation.

The file is read line-by-line. Lines that start with "!" are considered comments and are therefore discarded. Generally, the order in which the lines appear doesn't matter, but there are some exceptions, so try to keep the ordering shown in this example.

We can now go through the file one section at the time.

Mesh generation

We show how to generate a mesh file in .su2 format using GMSH. We assume a 2D axisymmetric domain with a sphere of radius $r_1 = 9.5 \ \text{mm}$ (as given by the experimental setup, although $r_2 = 12.5 \ \text{mm}$ can be also used) as shown in the picture below. Here we assume the flow or mercury plasma beam from the left boundary, while top and right boundary are taken as open boundary.

Screenshot from 2026-04-23 09-35-55

The mesh is generated from the mesh.geo file with the following input (while taking in account the domain constraints from physical considerations):

// INPUT FOR MESH.GEO FILE

// Domain parameters //
upstream_distance = -0.045;
downstream_distance = 0.08;
domain_width = 0.035;

sphere_radius = 0.0095;

cell_size_sphere = 0.0003; // Cell size on the sphere
cell_size_boundary = 0.005; // Cell size extending to domain boundaries

// ----- MESH ----- //
Point(1) = {upstream_distance, 0., 0., cell_size_boundary};
Point(2) = {downstream_distance, 0., 0., cell_size_boundary};
Point(3) = {downstream_distance, domain_width, 0., cell_size_boundary};
Point(4) = {upstream_distance, domain_width, 0., cell_size_boundary};

Point(5) = {-sphere_radius, 0., 0., cell_size_sphere};
Point(6) = {0., sphere_radius, 0., cell_size_sphere};
Point(7) = {sphere_radius, 0., 0., cell_size_sphere};
Point(8) = {0., 0., 0., cell_size_sphere};

Line(1) = {1,5};
Circle(2) = {5, 8, 6};
Circle(3) = {6, 8, 7};
Line(4) = {7,2};
Line(5) = {2,3};
Line(6) = {3,4};
Line(7) = {4,1};

Curve Loop(8) = {1,2,3,4,5,6,7};
Plane Surface(9) = {8};

Physical Curve("Sphere", 14) = {2, 3};
Physical Curve("Inlet", 10) = {7};
Physical Curve("Outlet", 11) = {5};
Physical Curve("Top", 13) = {6};
Physical Curve("Symmetry", 15) = {1,4};
Physical Surface("Domain", 16) = {9};

where the physical groups will be used for boundary condition definition in the Pantera input file.

The generated mesh should be then exported in .msh format of version 2 in ASCII format, which is supported by Pantera. SU2 format can be used as well, but in this case the user needs to check the 'Save all elements' option. In terminal, for the former, this is easily done with

$ gmsh mesh.geo -2 -format msh2

Simulation settings for input file

We can now proceed to the actual input file. The input file will consist of several keywords to set up the proper domain, boundary conditions and species-specific parameters. Please refer to commands wiki page to get a full description of all the keywords mentioned below

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!                 Input file for PANTERA                                !!!!
!!!!                                                                       !!!!
!!!! The program reads is line by line, and when it finds a known line, it !!!!
!!!! reads and saves the value in the next one.                            !!!!
!!!! All other lines are ignored (you can comment etc).                    !!!!
!!!! You can put trailing comments on lines, starting with "!" symbol:     !!!!
!!!! the program will remove them before checking the content of a line.   !!!!
!!!!                                                                       !!!!
!!!! All units are SI [m, kg, s, ...]                                      !!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


! ========= Computational domain and grid ===========
Dimensions:
2
Axisymmetric:
.TRUE.
Number_of_cells:
1    1    1
Domain_limits:
-0.1    0.1    0.    0.5    -0.5   0.5 ! xmin[m] xmax[m] ymin[m] ymax[m] zmin[rad] zmax[rad]
Domain_periodicity:
.FALSE.   .FALSE.   .FALSE.
Domain_specular:
.FALSE.    .FALSE.    .TRUE.    .FALSE. ! Here ymin is set to .TRUE. for axisymmetric case
Domain_diffuse:
.FALSE.    .FALSE.    .FALSE.   .FALSE.

! PIC Settings, we choose 'hybrid' to include electron fluid model
PIC_type:
hybrid

Mesh_file:
'./mesh.msh'

Partition_style:
stripsx !strips in x direction

! specify all paths or just "All_output_path:" + path
All_output_path:
'./dumps/'

where PIC_type: chooses the solver used for simulating electrostatic interaction, the value hybrid essentially allows for electrons to be modelled as equilibrium fluid.

Mesh_file: requires a path to the generated mesh file. In case of .su2 mesh file, keyword 'Mesh_file_SU2:' needs to be used.

! ========= Numerical parameters ====================

Fnum:  ! ratio of real-to-simulated particles
4e3
Timestep:
1e-8
Number_of_timesteps:
3000
Stats_every: ! screen output every th timestep
10
Save_initial_timestep:
.TRUE.

RNG_seed:
15950

Binary_output:
.TRUE.


Dump_grid_start:
0
Dump_grid_avgevery: ! average flow field every th time step
1
Dump_grid_numavgs:  ! how many averages before saving
100

Dump_bound_start:
0
Dump_bound_avgevery: ! average values on boundaries every th time step
1
Dump_bound_numavgs:  ! how many averages before saving
100

Dump_part_start: ! start dumping particles
1000000
Dump_part_every: ! dumping particles every th timestep
1000000

Perform_checks: ! .csv files -- saves integral quantitites (moments, etc)
.FALSE.

Domain_type:
Domain fluid 1.0

where Domain_type: defines the main domain part in which we solve for particles. It is the identical physical group that we saved in our mesh.geo file.

! ========= Multispecies ==============================

Species_file:
'mercury_exp.species'

! First: name   second: specie     third: molecular fraction    (if more: +specie - fraction)
Def_mixture:
Ions Hg+ 1.0

Def_mixture:
Neutrals Hg 1.0

! ========= Collisions ==============================

Collision_type: ! [MCC / DSMC / BGK / NONE]
NONE

Wall_reactions_file:
mercury_exp.wall

Here Species_file requires the path for the species files, which consists of particle species that we wish to simulate. In this example, we have the mercury neutrals Hg and Ions Hg+ as follows

Hg+  200.00        3.32E-25  0   0.0   0   0.0    0.0    1.0      1.0
Hg   200.00        3.32E-25  0   0.0   0   0.0    0.0    1.0      0.0

Wall_reactions_file: requires a path for a file that specifies possible reactions on walls (defined by physical groups). In this case, we expect that ions hitting the sphere walls will neutralize (by accepting an electron from the wall) and then bounce from the walls as neutrals. In our case, the reaction file consists of one reaction

Hg+ --> Hg
1.0

where the number below the reaction specifies the probabilty for such reaction to occur.

! ========= Initial particle seed ===================

! mixture   density     3 components velocity   3 components temperature    rotational_temp     vibrational_temp
Initial_particles:
Ions   1e14    10000.    0.0    0.    3000.    3000.    3000.    0.    0.

! Settings for fluid electrons when choosing 'hybrid' PIC type
! n0   phi0    Te0
Fluid_electrons:
1e14    0.  32810.

'Fluid_electrons' keyword is the main part in which we set the following properties of the electron fluid. This is only accepted when PIC_type is set to hybrid

! ========= Boundary emit conditions ================
Boundary_emit:
    Inlet   Ions    1e14  10000.  0.  0.  3000.   0. 0.
Boundary_emit:
    Top   Ions    1e14  10000.  0.  0.  3000.   0. 0.


! ========= Boundary conditions =====================

! BC for electrostatic potential
Boundary_condition:
    Inlet   dirichlet 0.0

Boundary_condition:
    Outlet  neumann 0.0

Boundary_condition:
    Top neumann 0.0

Boundary_condition:
    Symmetry neumann 0.0

! BC for sphere with particle BC
Boundary_condition: ! Potential of the sphere with respect to the beam
    Sphere dirichlet -100.

Boundary_condition: ! Particles diffusively scatter on the sphere with a temperature of 300 K
    Sphere diffuse 300.

Boundary_condition: ! Particle reactions on wall using 'Wall_reactions_file'
    Sphere react

Boundary_condition: sets boundary condition for either electrostatic potential or particles. For the elecrostatic potential, we follow what is depicted in the previous picture. Thus, grounding the inlet to zero while the axis and outlets are set to Neumann BC. The sphere maintains a single potential using the Dirichlet BC.

Since we have also defined the wall reactions, these need to be specified by the keyvalue react on the Sphere physical group.

Running the simulation

You should have a version of MPI installed on your machine and compiled Pantera.

In the same folder where you have pantera.exe, you should also have a folder named "dumps", and the three input files "input", "ar.species", and "arlofthouse.vss" and the generated ".su2" mesh file. You can then run the simulation with

mpirun -np [number of processes] ./pantera.exe

For hybrid simulations, it is highly recommended to add a following flag -ksp_type cg. This essentailly tells PETSc to use a conjugate gradient method, which heavily boosts the calculation speed for symmetrical matrices, which is our case.

Results

The flow field quantities are written directly in the VTK format, and can be visualized using Paraview.

Apart from the simulation files, a python script plot_force.py together with experimental and numerical reference datas are provided to the user

###########################################################################################################
# Script intended for a tutorial test case 02_charged_drag_on_sphere in Pantera PIC-DSMC code.            #
# --------------------------------------------------------------------------------------------------------#
# This script reads the simulation data from a VTK file, processes it to calculate the force              #
# on a sphere due to charged drag, and then plots the results alongside experimental data.                #
#                                                                                                         #
# NOTE: Make sure to run this script with Paraview's Python environment due to the paraview package usage #
###########################################################################################################

# Import necessary libraries
import numpy as np
import scipy.constants as sc
import matplotlib.pyplot as plt
from scipy.special import erf
from paraview.simple import *

# Path for VTK files generated by Pantera simulation with time step 2000, adjust if necessary
folder_path = 'dumps'
timestep = 3000


# Set up species parameters from the simulation
# We use mercury gas and assume boundary conditions in default input file for Pantera
particle_mass        = sc.m_u*200
particle_density     = 1e14
particle_velocity    = 10000
particle_temperature = 3000

sphere_radius      = 9.5e-3
sphere_temperature = 300
sphere_potential   = -100


##################################################
##################################################
##################################################

# ----------------------------
# ----- SOURCE FUNCTIONS -----
# ----------------------------

# Function to calculate the drag coefficient (C_d) on the sphere from analytical solution
def get_Cd(S,T,Tw):
    return np.exp(-S**2)/(np.sqrt(np.pi)*S**3)*(2*S**2+1) + erf(S)*(4*S**4+4*S**2-1)/(2*S**4) + 2*np.sqrt(np.pi)*np.sqrt(Tw/T)/3/S

# Function to read VTK file, extract data and calculate the charged force on the sphere
def get_force(folder,timestep,parameters):


    m, n, u, T, R, Tw, phi = parameters
    rho = m*n
    A = R**2*np.pi
    S = u/np.sqrt(2*sc.k*T/m)
    C_d = get_Cd(S,T,Tw)
    E = 0.5*m*u**2

    F_N = 0.5*rho*u**2*A*C_d

    alpha_arr = np.array([])
    force_arr = np.array([])

    # Load the VTK file
    data = OpenDataFile(folder+'/dsmc_boundary_'+str(timestep)+'.vtk')

    # Threshold to filter the specific PHYS_GROUP (adjusted for cell data)
    phys_group_id = 4  # Replace with the desired PHYS_GROUP ID (should be 4 from our settings)
    threshold = Threshold(Input=data)
    threshold.Scalars = ['CELLS', 'PHYSICAL_GROUP']
    threshold.LowerThreshold = phys_group_id
    threshold.UpperThreshold = phys_group_id

    geometry = ExtractSurface(Input=threshold)

    geometry.UpdatePipeline()
    data_info = servermanager.Fetch(geometry)

    # Sum variables (momentum + electromagnetic contribution)
    mom_sum = 0.0
    pem_sum = 0.0

    # Iterate over each line segment (cell)
    for j in range(data_info.GetNumberOfCells()):
        # Get cell length (assuming 2D line segments)
        cell = data_info.GetCell(j)
        p0 = cell.GetPoints().GetPoint(0)
        p1 = cell.GetPoints().GetPoint(1)
        length = ((p1[0] - p0[0])**2 + (p1[1] - p0[1])**2)**0.5 # Length of segment between points

        # Get the variable value for the cell (segment)
        pem_x = data_info.GetCellData().GetArray('pem_x').GetValue(j)
        mom_x = data_info.GetCellData().GetArray('mom_x_in_Hg+').GetValue(j)
        mom_x+= data_info.GetCellData().GetArray('mom_x_out_Hg+').GetValue(j)
        
        # Add contribution to total sum
        mom_sum += mom_x*length*(p0[1] + p1[1])/2*2*np.pi
        pem_sum += pem_x*length*(p0[1] + p1[1])/2*2*np.pi

    
    F_direct = mom_sum
    F_indirect = pem_sum

    alpha_arr = np.append(alpha_arr,sc.e*np.abs(phi)/E) # Normalized potential by the particle energy
    force_arr = np.append(force_arr,(F_direct+F_indirect)/F_N) # Normalized force by the neutral drag force

    return alpha_arr, force_arr

# Functionn to plot experimental and numerical results
def plot(folder,timestep,parameters):

    fig,axs = plt.subplots(1,1,figsize=(8,6))

    exp_data_r1 = np.loadtxt('reference_data/expdata_r9mm.csv',delimiter=',').T
    exp_data_r2 = np.loadtxt('reference_data/expdata_r12mm.csv',delimiter=',').T

    axs.plot(exp_data_r1[0],exp_data_r1[1],"^",
            label=r'Experiment (Knechtel et al, 1964) - $r_1 = 9.5 \: \mathrm{mm}$',
            color='darkred',
            markerfacecolor='none',
            ms=6)
    axs.plot(exp_data_r2[0],exp_data_r2[1],"^",
            label=r'Experiment (Knechtel et al, 1964) - $r_2 = 12.5 \: \mathrm{mm}$',
            color='darkblue',
            markerfacecolor='none',
            ms=6)

    capon_r9_n1 = np.loadtxt('reference_data/capon_r9_n1.csv',delimiter=',').T
    axs.plot(capon_r9_n1[0],capon_r9_n1[1],"s-",
            label=r'Numerical (Capon et al, 2019) - $r_1, n = 10^{14} \ \mathrm{m^{-3}}$',
            color='darkblue',
            ms=6)



    alpha_n1r1, force_n1r1 = get_force(folder,timestep,parameters)
    axs.plot(alpha_n1r1,force_n1r1,"o",
            label=r'Pantera simulation',
            color='red',
            mfc='none',
            ms=10,
            mew=3)


    axs.legend()
    fig.tight_layout()
    plt.show()


# Call the plotting function
plot(folder=folder_path, timestep=timestep,
     parameters = [particle_mass,
            particle_density,
            particle_velocity,
            particle_temperature,
            sphere_radius,
            sphere_temperature,
            sphere_potential])

By setting the appropriate time step of the VTK file, we can then compare the results with the reference data. In this case, the agreement should be satisfying.

Add a custom footer
⚠️ **GitHub.com Fallback** ⚠️