7.3. Hypersonic Flow around the 70° Cone (DSMC) - 2D Mesh

A widespread validation case for rarefied gas flows is the wind tunnel test of the 70° blunted cone in a hypersonic, diatomic nitrogen flow [72], [73]. Before beginning with the tutorial, copy the dsmc-cone directory from the tutorial folder in the top level directory to a separate location

cp -r $PICLAS_PATH/tutorials/dsmc-cone-2D .
cd dsmc-cone-2D

The general information needed to setup a DSMC simulation is given in the previous tutorial Adiabatic Box/Reservoir (DSMC, Relaxation/Chemistry). The following focuses on case-specific differences.

7.3.1. Mesh Generation with HOPR (pre-processing)

Before the actual simulation is conducted, a mesh file in the HDF5 format has to be supplied. The mesh files used by piclas are created by supplying an input file hopr.ini with the required information for a mesh that has either been created by an external mesh generator or directly from block-structured information in the hopr.ini file itself. Here, a conversion from an external mesh generator is required. The mesh and the corresponding boundaries are depicted in Fig. 7.7.

../../../_images/dsmc-cone-mesh-bcs.jpg — Fig. 7.7 Mesh of the 70° Cone.

7.3.1.1. Mesh generated by HEXPRESS (CGNS format)

An example setup for such a conversion is given in hopr.ini. In difference to a direct generation the FileName of the original mesh has to be specified and the Mode has to be set to 3, to read a CGNS mesh. For unstructured ANSA CGNS meshes the BugFix_ANSA_CGNS needs to be active.

!=============================================================================== !
! MESH
!=============================================================================== !
FileName         = 70degCone_2D_mesh.cgns
Mode             = 3
BugFix_ANSA_CGNS = T

The BoundaryName needs to be the same as in the CGNS file.

!=============================================================================== !
! BOUNDARY CONDITIONS
!=============================================================================== !
BoundaryName = IN
BoundaryType = (/3,0,0,0/)
BoundaryName = OUT
BoundaryType = (/3,0,0,0/)
BoundaryName = WALL
BoundaryType = (/4,0,0,0/)
BoundaryName = SYMAXIS
BoundaryType = (/4,0,0,0/)
BoundaryName = ROTSYM
BoundaryType = (/4,0,0,0/)

To create the .h5 mesh file, you would simply run

hopr hopr_cgns.ini

This would create the mesh file 70degCone_2D_mesh.h5 in HDF5 format. For more information see directly at https://github.com/hopr-framework/hopr.

7.3.1.2. Mesh generated by Gmsh (msh format)

Instead of using the given .cgns mesh file, the mesh can be generated by using the open-source mesh generator Gmsh.

First, create a new file in gmsh: 70DegCone_2DSurf.geo. In general, the mesh can be generated using the GUI or by using the .geo script environment. In the GUI, the script can be edited via Edit script and loaded with Reload script. This tutorial focuses on the scripting approach.

After opening the .geo script file, select the OpenCASCADE CAD kernel and open the provided 70DegCone_2DSurf.step file with the following commands:

SetFactory("OpenCASCADE");
v() = ShapeFromFile("70degCone_2DSurf.step");

Two-dimensional and axisymmetric simulations require a mesh in the $xy$-plane, where the $x$-axis is the rotational axis and $y$ ranges from zero to a positive value. Additionally, the mesh shall be centered around zero in the $z$-direction with a single cell row. As the loaded file only consists of a surface, which is meshed with Gmsh and then extruded using Hopr, the surface needs to be translated:

Translate {0, 0, -1} {
    Surface{1};
}

Physical groups are used to define the boundary conditions at all curves:

Physical Curve("SYMAXIS") = {3};
Physical Curve("WALL") = {1, 2};
Physical Curve("IN") = {4};
Physical Curve("OUT") = {5, 6};

The mesh options can be set with the following commands:

Mesh.MeshSizeMin = 0;
Mesh.MeshSizeMax = 10;
Field[1] = MathEval;
Field[1].F = "0.2";
Field[2] = Restrict;
Field[2].EdgesList = {1, 2};
Background Field = 2;
Mesh.Algorithm = 1;
Mesh.RecombinationAlgorithm = 2;

The commands Mesh.MeshSizeMin and Mesh.MeshSizeMax define the minimum and maximum mesh element sizes. With the prescribed Field options, the size of the mesh can be specified using an explicit mathematical function using MathEval and restriced to specific surfaces with Restrict. In this tutorial, a mesh refinement at the frontal wall of the cone is enabled with this. Background Field = 2 sets Field[2] as background field. Different meshing algorithms for creating the 2D and 3D meshes can be chosen within Gmsh.

Next, the 2D mesh is created and all cells are recombined to quads with the following commands:

Mesh.RecombineAll = 1;
Mesh 2;

The following commands are required to save all elements even if they are not part of a physical group and to use the ASCII format, before saving the mesh as 70degCone_2D.msh:

Mesh.SaveAll = 1;
Mesh.Binary = 0;
Mesh.MshFileVersion = 4.1;
Save "70degCone_2D.msh";

The mesh can be created by simplying executing Gmsh from the terminal:

gmsh 70degCone_2DSurf.geo

The resulting mesh shall consist of quad elements and not triangles. Finally, it has to be converted to the file format used by piclas using HOPR by supplying an input file hopr.ini using the corresponding mode:

Mode = 5

To extrude the 2D quadrangular mesh to a 3D hexahedral mesh with only one element in the third direction, the following commands next to be set in the hopr.ini:

MeshDim      = 2
zLength      = 1
nElemsZ      = 1
lowerZ_BC    = (/3,0,0,0/)
upperZ_BC    = (/3,0,0,0/)

While zLength specifies the length of the mesh in $z$-direction, nElemsZ ensures a single cell row. lowerZ_BC and upperZ_BC are the boundary conditions on the surfaces parallel to the $xy$-plane, that need to be defined as stated.

Again, to create the .h5 mesh file, you would simply run

hopr hopr_gmsh.ini

This would create the mesh file 70degCone_2D_mesh.h5 in HDF5 format.

7.3.2. Flow simulation with DSMC

Install piclas by compiling the source code as described in Section Installation and make sure to set the correct compile flags. For this setup, we are utilizing the regular Direct Simulation Monte Carlo (DSMC) method

PICLAS_TIMEDISCMETHOD = DSMC

or simply run the following command from inside the build directory

cmake ../ -DPICLAS_TIMEDISCMETHOD=DSMC

to configure the build process and run make afterwards to build the executable. It is recommended to either utilize a separate build folder (e.g. build_DSMC/) or to delete the contents of the folder beforehand to avoid conflicts between different compiler options (e.g. the setting PICLAS_EQNSYSNAME = poisson from the plasma wave tutorial is in conflict with the DSMC method). An overview over the available solver and discretization options is given in Section Solver settings. The values of the general physical properties are listed in Table 7.4.

Table 7.4 Physical properties at the simulation start
Property	Value (initial and surfaceflux)
Species	$\text{N}_2$
Molecule mass $m_{\text{N}_2}$	$\pu{4.65e-26 kg}$
Number density $n$	$\pu{3.715e+20 m^{-3}}$
Translational temperature $T_{\text{trans}}$	$\pu{13.3 K}$
Rotational temperature $T_{\text{rot}}$	$\pu{13.3 K}$
Vibrational temperature $T_{\text{vib}}$	$\pu{13.3 K}$
Velocity $v_{\text{x}}$	$\pu{1502.57 \frac{m}{s}}$

To define an incoming, continuous flow, the procedure is similar to that for initialization sets. For each species, the number of inflows is specified via Part-Species[$]-nSurfaceFluxBCs. Subsequently, the boundary from which it should flow in is selected via Part-Species[$]-Surfaceflux[$]-BC. The rest is identical to the initialization explained in Section Particle solver.

Part-Species1-nSurfaceFluxBCs = 1

Part-Species1-Surfaceflux1-BC                   = 1
Part-Species1-Surfaceflux1-velocityDistribution = maxwell_lpn
Part-Species1-Surfaceflux1-MWTemperatureIC      = 13.3
Part-Species1-Surfaceflux1-TempVib              = 13.3
Part-Species1-Surfaceflux1-TempRot              = 13.3
Part-Species1-Surfaceflux1-PartDensity          = 3.715E+20
Part-Species1-Surfaceflux1-VeloIC               = 1502.57
Part-Species1-Surfaceflux1-VeloVecIC            = (/1.,0.,0./)

7.3.2.1. Exploiting symmetry

In axially symmetrical cases, the simulation effort can be greatly reduced. For this, 2D must first be activated via Particles-Symmetry-Order = 2. Particles-Symmetry2DAxisymmetric = T enables axisymmetric simulations.

! Symmetry
Particles-Symmetry-Order         = 2
Particles-Symmetry2DAxisymmetric = T

First of all, certain requirements are placed on the grid. The $y$-axis acts as the symmetry axis, while the $x$-axis defines the radial direction. Therefore grid lies in the $xy$-plane and should have an extension of one cell in the $z$-direction, the extent in $z$-direction is irrelevant whilst it is centered on $z=0$. In addition, the boundary at $y = 0$ must be provided with the condition symmetric_axis and the boundaries parallel to the $xy$-plane with the condition symmetric.

Part-Boundary4-SourceName  = SYMAXIS
Part-Boundary4-Condition   = symmetric_axis

For the .cgns mesh, the following commands need to be enabled:

Part-nBounds               = 5
Part-Boundary5-SourceName  = ROTSYM
Part-Boundary5-Condition   = symmetric

For the .msh mesh instead, the following commands need to be enabled:

Part-nBounds               = 6
Part-Boundary5-SourceName  = lowerZ_BC
Part-Boundary5-Condition   = symmetric
Part-Boundary6-SourceName  = upperZ_BC
Part-Boundary6-Condition   = symmetric

To fully exploit rotational symmetry, a radial weighting can be enabled via Particles-RadialWeighting = T, which will linearly increase the weighting factor towards $y_{\text{max}}$, depending on the current $y$-position of the particle. Thereby the Particles-RadialWeighting-PartScaleFactor multiplied by the MacroParticleFactor is the weighting factor at $y_{\text{max}}$. Since this position based weighting requires an adaptive weighting factor, particle deletion and cloning is necessary. Particles-RadialWeighting-CloneDelay defines the number of iterations in which the information of the particles to be cloned are stored and Particles-RadialWeighting-CloneMode = 2 ensures that the particles from this list are inserted randomly after the delay.

! Radial Weighting
Particles-RadialWeighting                 = T
Particles-RadialWeighting-PartScaleFactor = 60
Particles-RadialWeighting-CloneMode       = 2
Particles-RadialWeighting-CloneDelay      = 5

For further information see 2D/Axisymmetric Simulations.

7.3.2.2. Octree

By default, a conventional statistical pairing algorithm randomly pairs particles within a cell. Here, the mesh should resolve the mean free path to avoid numerical diffusion. To circumvent this requirement, an octree-based sorting and cell refinement can be enabled by Particles-DSMC-UseOctree = T. The resulting grid is defined by the maximum number Particles-OctreePartNumNode and minimum number Particles-OctreePartNumNodeMin of particles in each subcell. Furthermore, the search for the nearest neighbour can be enabled by Particles-DSMC-UseNearestNeighbour = T.

! Octree
Particles-DSMC-UseOctree           = T
Particles-DSMC-UseNearestNeighbour = T
Particles-OctreePartNumNode        = 40
Particles-OctreePartNumNodeMin     = 28

For further information see Pairing & Collision Modelling.

7.3.2.3. Sampling

The outputs of the Projectname_DSMCState_Timestamp.h5 (data in domain) and Projectname_DSMCSurfState_Timestamp.h5 (surface data) files are based on sampling over several time steps. There are two different approaches that can not be used at the same time. The first method is based on specifying the sampling duration via Part-TimeFracForSampling as the fraction of the simulation end time (as defined by TEnd) between 0 and 1, where 0 means that no sampling occurs and 1 that sampling starts directly at the beginning of the simulation

$t_\text{samp,start} = T_\text{end} \cdot \left(1 - f_\text{samp}\right)$

Particles-NumberForDSMCOutputs then indicates how many samples are written in this period. The data is not discarded after the respective output and the sampling is continued. In other words, the last output contains the data for the entire sampling period.

Part-TimeFracForSampling          = 0.5
Particles-NumberForDSMCOutputs    = 2

The second method is activated via Part-WriteMacroValues = T. In this approach, Part-IterationForMacroVal defines the number of iterations that are used for one sample. After the first sample has been written, the data is discarded and the next sampling process is started.

Part-WriteMacroValues             = T
Part-IterationForMacroVal         = 1250

For further information see Particle Flow and Surface Sampling.

7.3.3. Run the simulation

Finally, you can start the simulation using the Message Passing Interface (MPI) on multiple cores

mpirun -np 8 piclas parameter.ini DSMC.ini | tee std.out

To continue a simulation after a successful run, you have to provide the state file (Projectname_State_Timestamp.h5) you want to restart from

mpirun -np 8 piclas parameter.ini DSMC.ini Projectname_State_Timestamp.h5 | tee std.out

The restart also redistributes the computational load and can thus significantly reduce the time to solution. In the following, additional automatic load balancing during the run-time is described.

7.3.3.1. MPI and Load Balancing

When using several cores, piclas divides the computing load by distributing the computing cells to the various cores. If a particle leaves the boundaries of a core, it is necessary that the surrounding grid is also known. This region is defined by the Particles-HaloEpsVelo and the time step. In general, it can be said that this velocity should be greater than or equal to the maximum velocity of any particle in the simulation. This prevents a particle from completely flying through this halo region during a time step.

Particles-HaloEpsVelo = 8.0E+4

If the conditions change, it could make sense to redistribute the computing load. An example is the build-up of a bow shock during the simulation time: While all cells have the same particle density during initialization, an imbalance will develop after a short time. The cores with cells in the area of the bow shock have significantly more computational effort, since the particle density is significantly higher. As mentioned at the beginning, piclas redistributes the computing load each time it is started.

The parameter Particles-MPIWeight indicates whether the distribution should be oriented more towards a uniform distribution of the cells (values less than 1) or a uniform distribution of the particles (values greater than 1). There are options in piclas to automate this process by defining load balancing steps during a single program call. For this, load balancing must have been activated when compiling piclas (which is the default). To activate load balancing based on the number of particles, DoLoadBalance = T and PartWeightLoadBalance = T must be set. piclas then decides after each Analyze_dt whether a redistribution is required. This is done using the definable Load DeviationThreshold. Should the maximum relative deviation of the calculation load be greater than this value, a load balancing step is carried out. If DoInitialAutoRestart = T and InitialAutoRestart-PartWeightLoadBalance = T are set, a restart is carried out after the first Analyze_dt regardless of the calculated imbalance. To restrict the number of restarts, LoadBalanceMaxSteps limits the number of all load balancing steps to the given number.

! Load Balancing
Particles-MPIWeight                      = 1000
DoLoadBalance                            = T
PartWeightLoadBalance                    = T
DoInitialAutoRestart                     = T
InitialAutoRestart-PartWeightLoadBalance = T
LoadBalanceMaxSteps                      = 2

Information about the imbalance are shown in the std.out and the ElemTimeStatistics.csv file. The default load balancing scheme will exchange the required data internally, but there is also the possibility to perform the re-balancing step via HDF5, which will create a state file and restart from this file by activating

UseH5IOLoadBalance = T ! default is False

7.3.4. Visualization (post-processing)

7.3.4.1. Ensuring physical simulation results

After running a simulation, especially if done for the first time it is strongly recommended to ensure the quality of the results. For this purpose, the Particles-DSMC-CalcQualityFactors = T should be set, to enable the calculation of quality factors such as the maximum collision probability and the mean collision separation distance over the mean free path. All needed datasets can be found in the *_DSMCState_*.h5 or the converted *_visuDSMC_*.vtu.

First of all, it should be ensured that a sufficient number of simulation particles were available for the averaging, which forms the basis of the shown data. The value *_SimPartNum indicates the average number of simulated particles in the respective cell. For a sufficient sampling size, it should be guaranteed that at least 10 particles are in each cell, however, this number is very case-specific. The value DSMC_MCSoverMFP is an other indicator for the quality of the particle discretization of the simulation area. A value above 1 indicates that the mean collision separation distance is greater than the mean free path, which is a signal for too few simulation particles. For 3D simulations it is sufficient to adjust the Part-Species[$]-MacroParticleFactor accordingly in parameter.ini. In 2D axisymmetric simulations, the associated scaling factors such as Particles-RadialWeighting-PartScaleFactor can also be optimized (see Section 2D/Axisymmetric Simulations).

Similarly, the values DSMC_MeanCollProb and DSMC_MaxCollProb should be below 1 in order to avoid nonphysical values. While the former indicates the averaged collision probability per timestep, the latter stores the maximum collision probability. If this limit is not met, more collisions should have ocurred within a time step than possible. A refinement of the time step ManualTimeStep in parameter.ini is therefore necessary. If a variable timestep is also used in the simulation, there are further options (see Section Variable Time Step).

Table 7.5 Target value to ensure physical results and a connected input parameter
Property	Target	Connected Input Parameter
*_SimPartNum	$\gt 10$	Part-Species[$]-MacroParticleFactor
DSMC_MCSoverMFP	$\lt 1$	Part-Species[$]-MacroParticleFactor
DSMC_MaxCollProb	$\lt 1$	ManualTimeStep

Finally, the time step and particle discretization choice is a trade-off between accuracy and computational time. For further information see Section Ensuring Physical Simulation Results.

7.3.4.2. Visualizing flow field variables (DSMCState)

To visualize the data which represents the properties in the domain (e.g. temperatures, velocities, …) the DSMCState-files are needed. They are converted using the program piclas2vtk into the VTK format suitable for ParaView, VisIt or many other visualisation tools. Run the command

./piclas2vtk dsmc_cone_DSMCState_000.00*

to generate the corresponding VTK files, which can then be loaded into your visualization tool. The resulting translational temperature and velocity in the domain are shown in Fig. 7.8. The visualized variables are Total_TempTransMean, which is mean translational temperature and the magnitude of the velocities Total_VeloX, Total_VeloX, Total_VeloX (which is automatically generated by ParaView). Since the data is stored on the original mesh (and not the internally refined octree mesh), the data initially looks as shown in the two upper halves. ParaView offers the possibility to interpolate this data using the CellDatatoPointData filter. The data visualized in this way can be seen in the lower half of the image.

../../../_images/dsmc-cone-visu.jpg — Fig. 7.8 Translational temperature and velocity in front of the 70° Cone, top: original data; bottom: interpolated data.

7.3.4.3. Visualizing surface variables (DSMCSurfState)

For postprocessing and visualization, the parameter TimeStampLength = 13 is set inparameter.ini . This limits the output filename length. This can be needed, as e.g. Paraview may sort the files incorrectly and display a faulty time solution. To visualize the data which represents the properties at closed boundaries (e.g. heat flux, force per area, etc. the DSMCSurfState-files are needed. They are converted using the program piclas2vtk into the VTK format suitable for ParaView, VisIt or many other visualization tools. Run the command

./piclas2vtk dsmc_cone_DSMCSurfState_000.00*

to generate the corresponding VTK files, which can then be loaded into your visualization tool. A comparison between experimental data by [72] and the simulation data stored in dsmc_cone_visuSurf_000.002000000.vtu is shown at Fig. 7.9. Further information about this comparison can be found at [74].

../../../_images/dsmc-cone-heatflux.svg — Fig. 7.9 Experimental heat flux data compared with simulation results from PIClas.

Property	Value (initial and surfaceflux)
Species	\(\text{N}_2\)
Molecule mass \(m_{\text{N}_2}\)	\(\pu{4.65e-26 kg}\)
Number density \(n\)	\(\pu{3.715e+20 m^{-3}}\)
Translational temperature \(T_{\text{trans}}\)	\(\pu{13.3 K}\)
Rotational temperature \(T_{\text{rot}}\)	\(\pu{13.3 K}\)
Vibrational temperature \(T_{\text{vib}}\)	\(\pu{13.3 K}\)
Velocity \(v_{\text{x}}\)	\(\pu{1502.57 \frac{m}{s}}\)