DSMC and PICMC documentation
2024-09-20
5 Direct Simulation Monte Carlo (DSMC) - examples
5.1 Cylindrical tube flow
5.1.1 Geometric model
In this example, we analyze the gas flow through a cylindrical tube. The geometric model as shown in Fig. 5.1 contains the following elements:
- The cylindrical tube
- A buffer chamber at the left side of the tube for simulation of the gas inflow
- A buffer chamber at the right side of the tube for simulation of the gas outflow
- A gas inlet surface within the left buffer chamber
- A pumping surface within the right buffer chamber
- A box surrounding all above mentioned elements
5.1.2 General considerations about geometries for gas flow simulation
In order to create such a model, we use the open source meshing and postprocessing tool GMSH (Geuzaine and Remacle (2009)). There are two possible ways of creating geometric model with GMSH: Either the model can be created manually by GMSH's own interactive mode and scripting language, or an IGS file (or STEP or other formats) originating from a CAD system can be imported. For our example, we use the manual way in order to point out the most relevant parts of a geometric model definition. Before starting to create the model, the following points are to be considered:
- At first, it is important to mention that we will create a surface mesh rather than a volume mesh. Within the DSMC/PIC-MC code, the surface mesh will be used to track collisions between particles and chamber walls, while collisions between the particles will be treated within a cartesian volume grid, which is created independently from the surface mesh model.
- Second, each wall has to be represented by at least two surfaces (for front and back side). This is important for determining, whether a given point lies within a wall or within the free space. The algorithm which decides about that is pointed out in Fig. 5.2. From a given point in space, a virtual line is drawn along X axis1, and the number of intersections between this line and any mesh surface is counted. For an even number of intersections, the point is considered to be within the gas flow volume, while for an uneven number of intersections, the point is considered to be within a solid wall.
Depending on the position of the point with respect to the edges of the simulation volume, the line is either drawn in positive or negative X direction.↩︎
5.1.3 Definition of geometric points
Before creating the model, we need to think about some conventions about the coordinates, orientation, units etc. In this example, we stick to the following conventions:
- All length units are in mm
- The tube is oriented along X axis
- The coordinate origin (0,0,0) is located at the left end of the tube, centralized within the circular cross section
A geometric model always consists of points, lines, and surfaces. First, the geometric points have to be defined, second, the points are connected with lines and finally, surfaces are defined by line loops. The definition of point coordinates is done in GMSH’s scripting engine as shown in the following.
First, some global variables are defined in the file header
= 10; // characterictic mesh size
LC
= 400;
L_TUBE = 50;
R_TUBE
= 100;
L_BUFFER = 200;
R_BUFFER
= 10; WALL
Next is the definition of inlet and outlet points for the cylindrical tube
// Tube orientation along x axis
// Center point at left tube edge
(101) = {0, 0, 0, LC};
Point
// Four points on tube radius separated by 90 degree
(102) = {0, R_TUBE, 0, LC};
Point(103) = {0, 0, R_TUBE, LC};
Point(104) = {0, -R_TUBE, 0, LC};
Point(105) = {0, 0, -R_TUBE, LC};
Point
// The same is repeated for the right tube edge
(106) = {L_TUBE, 0, 0, LC};
Point(107) = {L_TUBE, R_TUBE, 0, LC};
Point(108) = {L_TUBE, 0, R_TUBE, LC};
Point(109) = {L_TUBE, -R_TUBE, 0, LC};
Point(110) = {L_TUBE, 0, -R_TUBE, LC}; Point
Adjacent to both tube ends are located two larger cylinders which act as buffer chambers for the gas inlet and for the pumping.
// Additionally, we need buffer chambers at
// left and right side of tube
// Left buffer chamber:
(201) = {-L_BUFFER, 0, 0, LC};
Point(202) = {-L_BUFFER, R_BUFFER, 0, LC};
Point(203) = {-L_BUFFER, 0, R_BUFFER, LC};
Point(204) = {-L_BUFFER, -R_BUFFER, 0, LC};
Point(205) = {-L_BUFFER, 0, -R_BUFFER, LC};
Point(206) = {0, R_BUFFER, 0, LC};
Point(207) = {0, 0, R_BUFFER, LC};
Point(208) = {0, -R_BUFFER, 0, LC};
Point(209) = {0, 0, -R_BUFFER, LC};
Point
// Right buffer chamber:
(302) = {L_TUBE, R_BUFFER, 0, LC};
Point(303) = {L_TUBE, 0, R_BUFFER, LC};
Point(304) = {L_TUBE, -R_BUFFER, 0, LC};
Point(305) = {L_TUBE, 0, -R_BUFFER, LC};
Point(301) = {L_TUBE+L_BUFFER, 0, 0, LC};
Point(306) = {L_TUBE+L_BUFFER, R_BUFFER, 0, LC};
Point(307) = {L_TUBE+L_BUFFER, 0, R_BUFFER, LC};
Point(308) = {L_TUBE+L_BUFFER, -R_BUFFER, 0, LC};
Point(309) = {L_TUBE+L_BUFFER, 0, -R_BUFFER, LC}; Point
Finally, the whole arrangement is surrounded by an outer wall, which is shaped as a rectangular box.
// Bounding box around the whole setup
= -L_BUFFER-WALL;
x0 = L_BUFFER+L_TUBE+WALL;
x1 = -R_BUFFER-WALL;
y0 = R_BUFFER+WALL;
y1 = -R_BUFFER-WALL;
z0 = R_BUFFER+WALL;
z1
(401) = {x0, y0, z0, LC*5};
Point(402) = {x1, y0, z0, LC*5};
Point(403) = {x0, y1, z0, LC*5};
Point(404) = {x1, y1, z0, LC*5};
Point(405) = {x0, y0, z1, LC*5};
Point(406) = {x1, y0, z1, LC*5};
Point(407) = {x0, y1, z1, LC*5};
Point(408) = {x1, y1, z1, LC*5}; Point
A few comments on this code are given in the following:
- In the first line, a variable
LC
is defined, which is subsequently used as fourth argument for eachPoint
. This is a characteristic length and defines the average mesh size around the points. - Further characteristics dimensions such as the tube length
(
L_TUBE
), its radius etc. are also defined as variables, which are subsequently being used in the point coordinates. This simplifies to create variations of the geometric model with different tube lengths etc. later on. - The syntax of a point definition in GMSH is
Point(number)={x, y, z, lc};
The point number must be unique within the whole GMSH model. It is not required that the set of numbers is continuous. It is convenient to assign the same major number to points which belong to a certain part of the geometry, e.- 1XX for the tube, 2XX for the left buffer chamber etc.
- Wherever a coordinate is required in GMSH it is possible to make use of a built-in formula parser.
- If you copy the above script into an ASCII file (named e. g. 'tube3d.geo') and open this file with GMSH, you will get a picture similar as in Fig. 5.3.
To do so, open the Tools->Options dialog, click on "Geometry" at the left side and chose the "Visibility" ribbon↩︎
5.1.4 Definition of straight lines and circle arcs
The next step is to connect all points with straight or circular lines. This can be either done by further script commands or by using GMSH interactively. We will follow the interactive way. First, use the GMSH menu sequence shown in Fig. 5.4 in order to select "circle arcs".
A circle arc is now created by clicking on three points subsequently: (i) A point on the circle as start point, (ii) the center of the circle and (iii) another point on the circle as end point. The left circular entrance of the tube is e. g. defined by clicking on points no. 102-->101-->103; 103-->101-->104; 104-->101-->105; 105-->101-->102. The right circular entrance of the tube is defined analogously (107-->106-->108; 108-->106-->109; 109-->106-->110; 110-->106-->107).
After having defined both circular surfaces, please click on "Straight line" in the GMSH menu and connect the following points: 102-->107; 103-->108; 104-->109; 105-->110. The interactive GMSH graph should now look as in Fig. 5.5.
While interactively defining lines between the points, GMSH
automatically appends accordant script lines to the geometry script. If
you look at the ASCII file tube2d.geo
, you will find the
following lines automatically appended at its end3.
(1) = {102, 101, 103};
Circle(2) = {103, 101, 104};
Circle(3) = {104, 101, 105};
Circle(4) = {105, 101, 102};
Circle(5) = {108, 106, 109};
Circle(6) = {109, 106, 110};
Circle(7) = {110, 106, 107};
Circle(8) = {107, 106, 108};
Circle(9) = {102, 107};
Line(10) = {103, 108};
Line(11) = {104, 109};
Line(12) = {105, 110}; Line
By continuing in the same manner, all other circle arcs and straight lines are to be defined. At the end, each point will be part of at least one line or circle arc, and the geometry will look like in Fig. 5.6. At this stage, please press the "q" button in order to leave the line selection mode. For your convenience, a zipped file of this stage (i. e. all points and lines but no surfaces) is provided for download.
If the file is still open in an editor, you may need to use a refresh function e.g. by pressing F5 in order to see the appended lines. It is not recommended to use the Windows notepad editor, instead a more functional editor such as Programmer's notepad is recommended, which will automatically notice any changes externally applied to the text file↩︎
5.1.5 Definition of ruled and plain surfaces
The next step will be to define the surfaces. In this case, we have to distinguish between
- Plane surfaces and
- Ruled surfaces.
Ruled surfaces are required for the cylindrical walls of the tube and buffer chambers. By the straight lines, the cyclindrical wall of the tube is divided into four parts. Each part has to be selected, individually.
In order to start building surfaces, please select the entry "Ruled Surface" within the GMSH menu shown in Fig. 5.4. Subsequently, select one of the four possible line loops on the cylindrical wall as shown in Fig. 5.7. This is accomplished by clicking onto lines rather than onto points. If you mis-clicked a line, pressing on "u" will undo the last selection. When all appropriate lines are highlighted as shown in Fig. 5.7, press "e" in order to confirm the selection.
This process has to be repeated for all four partial surfaces of the tube as well as of both buffer chambers. In total, twelfe "Ruled Surfaces" are required. After having defined all ruled surfaces, please press "2" in order to obtain a two-dimensional surface mesh, which will yield a picture similar as in Fig. 5.8.
Finally, the remaining plane surfaces should be defined. These are the plane surfaces of the surrounding box as well as front and back surfaces of the buffer chambers. Important: There should be no surfaces defined for both tube entrances, since otherwise the gas flow will be blocked. Instead, the adjacent surface of the buffer chamber is to be defined with the tube entrance as hole boundary. This is accomplished by first selecting the outer line loop of the large buffer chamber and then selecting the inner line loop of the tube entrance and finally pressing on "e". Fig. 9 shows a possible way of doing this. Please proceed with all remaining plain surfaces in order to come to the final step of defining physical surface groups. The attached ZIP file contains the geometry model with all surfaces but without any physical surface group.
5.1.6 Physical surface groups
The surfaces of our tube model are supposed to have different physical functions. While most of the surfaces are just acting as walls, where particles are being diffusively reflected, there is also an inlet and a pumping surface. Thus, we expect to define three different physical surface groups:
- Wall --> All surfaces except from inlet and outlet surfaces
- Inlet --> The left surface of the left buffer chamber
- Pump --> The right surface of the right buffer chamber
In order to select physical surface groups in GMSH, please use the menu sequence as shown in Fig. 5.10.
An appropriate strategy is to first select all surfaces in the model in order to form the physical surface group for "Wall". Afterwards, the "Inlet" and the "Pumping" surface should be individually selected. In order to select all surfaces, the following steps are required:
- With the center mousebutton, select an appropriate zoom in the GMSH window, where the whole geometry becomes visible.
- Hold down the "Ctrl" key and click once at a location in the upper left from the geometry.
- Release the "Ctrl" key and move the mouse cursor to the lower right. A dotted frame will appear, when this frame encloses the whole geometry, click on the mouse button the second time (see also Fig. 5.11 below).
- All surfaces will appear in red now, if this is the case, press "e" to confirm the selection.
Afterwards, please select the left plane surface of the left buffer chamber as well as the right surface of the right buffer chamber individually. After having completed the physical surface selection, GMSH has added three more lines to the geometry script, which look similar as follows:
(95) = {60, 92, 62, 52, 54, 72, 77, 86, 50, 58, 64, 56, 70, 75, 84, 88, 80, 66, 68, 82, 90, 94};
Physical Surface(96) = {77};
Physical Surface(97) = {82}; Physical Surface
Here, we have to apply some changes manually. First, we would rather use strings than numbers in order to identify the physical surfaces. This can be done by replacing the numbers in round brackets by quoted strings as follows:
("Wall") = {60, 92, 62, 52, 54, 72, 77, 86, 50, 58, 64, 56, 70, 75, 84, 88, 80, 66, 68, 82, 90, 94};
Physical Surface("Inlet") = {77};
Physical Surface("Pump") = {82}; Physical Surface
Since, we have selected all surfaces for the first physical surface group "Wall", the surfaces number 77 and 82 are now selected twice. These ambiguities must be removed; thus, the duplicate surface numbers (in our case 77 and 82) should be removed from the first physical surface group "Wall". As a final result, the last three lines of the GMSH file should now look as in the following:
("Wall") = {60, 92, 62, 52, 54, 72, 86, 50, 58, 64, 56, 70, 75, 84, 88, 80, 66, 68, 90, 94};
Physical Surface("Inlet") = {77};
Physical Surface("Pump") = {82}; Physical Surface
5.1.7 Final step: Saving the mesh file
As a final step, a 2D surface mesh is created by pressing the "2" button. Afterwards, the mesh file can be saved by clicking on the "save" button of the GMSH main dialog. A complete version of the tube model can be found as "tube.geo" within this ZIP file.
5.2 Setting up the parameter file
After having created the GMSH mesh file, which will be named
tube.msh
here, the next step is to copy this file to the
Linux cluster. In this example, it is copied into
/home/ap/sim/Tube
. In order to create an appropriate
parameter file, use the command initpicmc tube
. The
following will appear on the terminal:
ap@sim037:~/sim/Tube$ initpicmc tube
***** scr/init.r, 15642c0 *****
***** tube, 15642c0 *****
** Assuming NON-STEPFILE mode
** The mesh has 9978 nodes and 3 codecs
** The mesh dimensions are {-110, -210, -210} --> {510, 210, 210}
** Accordant to the mesh geomtry a parameter file tube.par has been created.
** Please edit the parameter file and invoke: rvmmpi -picmc [NumTasks] tube.par
Physical strings:
1 <--> 'Wall'
2 <--> 'Inlet'
3 <--> 'Pump'
***** tube, 15642c0 *****
ap@sim037:~/sim/Tube$ _
The output comprises the following informations:
- The initpicmc script recognized, that all mesh elements are contained within a box between (-110,-210,-210)-->(510,210,210). This volume will be set as simulation volume ("bounding box") within the parameter file.
- Three physical surface groups are identified.
- A parameter file "tube.par" has been created.
The parameter file is parsed by "RIG-VM", which is a script interpreter developed at Fraunhofer IST. This allows to use variables and formula for definition of numerical values. Comments can be placed within the parameter file by using the "#" character. The different sections of the parameter file are discussed in section (sec:par_file?).
5.2.1 Global control variables
The first section of the parameter file comprises some global variables for the simulation run and its runtime behavior. It should be edited like follows:
= 1e-5; # Time step width
DT = 1.0; # Total simulation time in seconds
TSIM = 1; # Time averaging mode (0 = off, 1 = on)
TAVG = 10; # Number of plots
NPLOTS = 1; # Number of particle data dumps
NDUMPS = 0; # Work mode (0 = DSCM gasflow simuation,
WORKMODE # 1 = PICMC plasma simulation,
# 2 = electric field solver)
= 300; # Initial particle and border temperature [K]
T0 = 1; # Aimed power dissipation in Watt (plasma simulation only)
PSOLL = 1e-7; # Sampling time for power dissipation (plasma simulation only)
PSAMPLE #VPMAX = 1000; # Maximum voltage for power control (default = 1000 volt)
5.2.2 Gas species and cell resolution
The next lines within the parameter file deal with the selection of gas species, their initial distribution and their statistical weight. In our test case, we will define a gas inlet operating with Argon at a pressure of 1.0 Pa and a pumping surface. This means, that the pressure will reach 1 Pa close to the gas inlet but will be significantly lower in the major part of the simulation volume. For Argon with a pressure of 1 Pa, the mean free path is approx. 6 mm, thus a cell dimension of about 0.5 mm would be required close to the gas inlet. Since we assume that the pressure in the tube will be significantly lower, we will start with a homogeneous pressure distribution of 0.5 Pa and use 1.0 cm as appropriate cell dimension.
In section ((sec:num_constraints_dsmc?)),
for Argon with a pressure of 1.0 Pa, a cell dimension of 5 mm and a
value of NREAL_Ar=3e12;
is suggested, which leads to an
average of 10 simulation particles per cell. To transfer this to the
case with a pressure 0.5 Pa with 10 mm cell dimension, we need to keep
in mind that NREAL_Ar
scales reciprocally with the square
of the pressure. Accordingly, the resulting value for
NREAL_Ar
is 1.2 x 1013. Thus, the lines within
the parameter file should be edited as follows:
= ["Ar"]; # List of species strings
SPECIES = 1.2e13; # Default scale factor for Ar
NREAL_Ar = 0.5; # Default pressure of initial Ar distribution [Pa] P0_Ar
In the next lines of the parameter file, the desired output of the simulation run is specified:
= ["Ar"];
NUMBER = ["Ar"];
PRESSURE = ["Ar"];
VELOCITY = ["PHI"]; FIELD
The next section of the parameter file deals with the
bounding box, i.e. the simulation volume and its cell
resolution. First, we have to define the length unit, which is used
within the mesh file. Our geometry model uses mm as
length unit, while the default setting in the parameter file is
m. Thus, it is required to change the length
UNIT
to 0.001.
Next, the bounding box coordinates are shown within the array
BB
. It is possible to edit the BB
array in
order to restrict the simulation volume to a fraction of the mesh
geometry, but we will use the full volume, here. The six numbers define
two points in 3D with lower and upper coordinates. In our case, the
bounding box spans between (-110,-210,-210)-->(510,210,210).
Accordingly, the length of the bounding box is 62 cm in X direction and
42 cm in Y and Z direction. With the cell dimension of 10 mm as
discussed above, we would like to have 62 cells along X direction and 42
cells in Y and Z direction, leading to a total of 109368 cells.
The cell resolution and segmentation of the simulation volume is
defined in the arrays NX, NY, NZ
. The segmentation is
important with respect to parallel processing: Each segment of the
simulation volume can be assigned to an indidividual CPU core. By this
way, the number of segments limits the maximal number of parallel
processes4. In the default settings, the
simulation volume is divided into two segments along each axis, and each
segment has a resolution of 10 grid cells. Given our model dimension,
this seems not appropriate. For our purpose, we should change these
lines to:
= 0.001; # Unit of mesh geometry (1.0 [m], 0.001 [mm], 0.0254 [inch], ...)
UNIT = [-110, -210, -210, 510, 210, 210]; # Bounding box dimensions
BB = [62]; # Cell resolution of subvolumes in x-direction
NX = [14, 14, 14]; # Cell resolution of subvolumes in y-direction
NY = [14, 14, 14]; # Cell resolution of subvolumes in z-direction NZ
By these settings, we divide the volume into three segments in Y and Z direction, while only one segment is used in X direction (see Fig. 5.12). The total number of segments is 3x3 = 9; thus, a maximal number of 9 parallel PIC-MC processes can be used for this simulation model. The numbers given for each segment denote the number of cells. In Y and Z direction we have have the 42 cells divided into three segments, i.e. 14 cells per segment, while in X direction we use 62 cells for one segment.
The selection of species, their statistical weights and the definition of the cell resolution is the most crucial part in the parameter file. The remaining steps are comparably easy.
However it is always possible to have less parallel processes than segments. In this case, some of the PIC-MC processes have to handle multiple segments↩︎
5.2.3 Simulation domain - bounding box
The simulation domain has the shape of a cuboid, which may be divided into segments and sub-cells. It is required to specify the boundary conditions in case, if a simulation particle reaches a cuboid surface. Possible boundary conditions are outlet, specular and periodic. We have to distinguish between the following cases:
- In 3D models, usually the gas volume is contained into a closed mesh
geometry. In this case, no particle should reach the cuboid boundary. If
the mesh geometry is not rectangular, there may be particles originating
from an initial homogeneous distribution, which are outside the
gas volume. In all such cases the most appropriate boundary condition is
outlet, i. e. a particle reaching the cuboid boundary
will be deleted. The respective lines in the parameter file should look
like
follows:
TYPE_x1 = "outlet"; TYPE_x2 = "outlet"; TYPE_y1 = "outlet"; TYPE_y2 = "outlet"; TYPE_z1 = "outlet"; TYPE_z2 = "outlet";
{.C#} - In some cases, there may be 3D models with a mirror symmetry plane.
In such cases it is possible to reduce computation time by factor of two
by placing one of the cuboid boundaries at the plane of symmetry. The
appropriate boundary condition for such plane is
specular, i. e. a particle hitting this surface will be
specularly reflected. An example, where the cuboid boundary at lower Y
coordinates defines a symmetry plane, is given as
follows:
TYPE_x1 = "outlet"; TYPE_x2 = "outlet"; TYPE_y1 = "specular"; TYPE_y2 = "outlet"; TYPE_z1 = "outlet"; TYPE_z2 = "outlet";
{.C#} - In 2D models, there are only mesh geometry elements for boundaries within e. g. the XY plane, while the geometry remains open in Z direction. In this case, both cuboid boundaries perpendicular to the Z axis should be either specular or periodic5.
- In 1D models, four of the six boundaries must be set on
periodic or specular. This may look
like
follows:
TYPE_x1 = "periodic"; TYPE_x2 = "periodic"; TYPE_y1 = "periodic"; TYPE_y2 = "periodic"; TYPE_z1 = "outlet"; TYPE_z2 = "outlet";
{.C#}
This makes no difference for DSMC gas flow simulation; for PIC-MC plasma simulation we recommend to use rather periodic↩︎
5.2.4 Surface functions
Finally, the geometric mesh file is to be specified, and the functions of the respective physical surfaces must be defined. The specification of the mesh file's filename follows just after the boundary conditions of the bounding box:
# ==================================================
# Specify mesh files and boundary conditions
# Available border types are: "wall" and "membrane"
# ==================================================
= "sample.msh"; # Name of mesh file for PICMC particle simulation
MESHFILE = ""; # Name of mesh file for BEM b-Field solver BEMMESH
Here, only for MESHFILE
a filename must be specified;
BEMMESH
is only used for computation of the magnetic field
in PVD plasma simulation setups and can be left empty, here.
In the following parameter file section, the functions of the
respective physical surfaces can be defined. A detailed description
about the possibilities here is given in section (sec:prep_wall_model?).
In our case we just need a particle source at surface Inlet
and an absorption coefficient at surface Pump
. The surface
Wall
is by default diffusely reflective at a temperature
T0
and can be left unchanged. With these two modifications,
the wall function definition block would look like follows:
block BORDER_TYPES: {= {
Border Wall = 1; # Codec number from meshfile
icodec type = "wall"; # Alternative type: "membrane"
= PAR.T0; # Wall temperature
T #epsr = 1.0; # Rel. dielectric permittivity (insulator only)
#vf = "VP*cos(2*PI*13.56e6*TIME)"; # Voltage function (conductor only)
#vf_coupling = ""; # Border name of counter electrode
reactions: {# Add reactions here - please check SIM-Wiki for available reaction schemes!
;
}
= {
Border Inlet = 2; # Codec number from meshfile
icodec type = "wall"; # Alternative type: "membrane"
= PAR.T0; # Wall temperature
T # epsr = 1.0; # Rel. dielectric permittivity (insulator only)
# vf = "VP*cos(2*PI*13.56e6*TIME)"; # Voltage function (conductor only)
# vf_coupling = ""; # Border name of counter electrode
reactions: {"Ar", 1.0, "Pa");
add_source(;
};
}
= {
Border Pump = 3; # Codec number from meshfile
icodec type = "wall"; # Alternative type: "membrane"
= PAR.T0; # Wall temperature
T # epsr = 1.0; # Rel. dielectric permittivity (insulator only)
# vf = "VP*cos(2*PI*13.56e6*TIME)"; # Voltage function (conductor only)
# vf_coupling = ""; # Border name of counter electrode
reactions: {"Ar", 0.25);
add_absorption(;
};
}; }
The completed simulation case including GEO and PAR file can be downloaded here.
5.3 Cylindrical tube flow - optimized domain decomposition
While the concept of embedding the whole mesh geometry into a single "bounding box" is simple and user-friendly, this method is not always efficient. In some cases, only a small fraction of the simulation volume will be actually used; thus the simulation model will consists of a majority of empty cells. In other cases, there may be big pressure differences within the simulation volume, which makes the definition of strongly varying cell dimensions necessary. With one single bounding box, this is not always possible in a flexible manner.
For this purpose, it is possible since picmc version 2.3.4 to use multiple joined rectangular bounding boxes as simulation volume:
- The bounding boxes can be defined in the parameter file with ''block MQUAD: { ... }; ''
- The connections between the bounding boxes are defined within ''block MJOIN: { ... }; ''
- If these blocks are omitted (or commented out), picmc will use the previous method of a large single bounding box.
- Due to limitations of the electric field solver with respect to multi-gridding schemes, the MQUAD / MJOIN feature is currently only available for DSMC simulation (WORKMODE=0 in the parameter file).
An example of the usage of multiple bounding boxes is shown in the following picture; it is based on the same tube geometry as in the DSMC example of section (sec:dsmc_example_1?).
Single box with wasted simulation volume | Use of multiple boxes |
---|---|
In the following, it will be shown how to realize the multiple boxes arrangement in the parameter file.
5.3.1 Tube example with four bounding boxes
In the tube example of the previous section, we have an inlet and an outlet chamber which are connected by a long tube. With one single bounding box, the simulation volume has a much larger cross section than the tube diameter resulting in a large number of wasted simulation cells. Thus, we would like to divide the simulation volume into the following bounding boxes,
- one for the inlet region,
- one for the left side of the tube with higher pressure,
- one for the right side of the tube with lower pressure,
- one for the outlet region.
The GEO and PAR file of the tube example with multiquad approach can be downloaded here as ZIP file.
5.3.1.1 Definition of multiple bounding boxes in block MQUAD
In order to define that in the parameter file, a
block MQUAD
must be defined. Within this block, all desired
bounding boxes have to be defined as shown in the following:
# ... geometry definition within parameter file
= 0.001; # Unit of mesh geometry
UNIT
block MQUAD: {
block Left: {= [23, 23, 23, 23, 23]; NX = [4, 4, 4, 4, 4];
DIV_X = [140, 140, 140]; NY = [28, 28, 28];
DIV_Y = [140, 140, 140]; NZ = [28, 28, 28];
DIV_Z = -110; Y0 = -210; Z0 = -210;
X0 =0.5;
NREAL_SCALE;
}
block Tube1: {= [65, 65, 65]; NX = [12, 10, 8];
DIV_X = [105]; NY = [20];
DIV_Y = [105]; NZ = [20];
DIV_Z = [0.5, 0.5, 1];
NREAL_SCALE_X ;
}
block Tube2: {= [195]; NX = [18];
DIV_X = [105]; NY = [10];
DIV_Y = [105]; NZ = [10];
DIV_Z ;
}
block Right: {=[12]; DIV_X=[115];
NX=[8, 8, 8]; DIV_Y=[140, 140, 140];
NY=[8, 8, 8]; DIV_Z=[140, 140, 140];
NZ;
}; }
The following information is to be specified:
For each bounding box, a block has to be specified with the following syntax:
block boundingbox_name: {# ... all parameters of the
# ... bounding box go here
; }
The geometric dimension in X, Y, Z direction is specified in the
arrays DIV_X
, DIV_Y
, DIV_Z
. If e.
g. a box extends over 195 mm in X direction and should have three
segments whereof each extends 65 mm, the specification reads:
= [65, 65, 65]; DIV_X
If only one segment in X direction is needed, the line reads:
= [195]; DIV_X
The number of cells for each sub-segment in X, Y and Z direction is
specified in the arrays NX
, NY
,
NZ
. The number of entries in e. g. NX
must be
the same as the number of entries in DIV_X
.
Further, it is now possible to specify a different NREAL
particle weighting factor in each bounding box. This is performed by
specifying:
= 0.5; NREAL_SCALE
which means, that the particle weighting factor within this bounding is 50% of the global particle weighting factor.
In case of a sub-segmentation It is also possible to specify different weighting factors for each segment along one of the directions X, Y, Z. If we have e. g. three sub-segments along X direction it is possible to specify three different weighting factors by:
= [0.25, 0.5, 1]; NREAL_SCALE_X
Exactly one of the bounding boxes needs absolute
coordinates, these are to be specifed in the variables X0
,
Y0
, Z0
. The line:
= -110; Y0 = -210; Z0 = -210; X0
means, that the lower coordinate edge of the bounding box is located
at (-110, -210, -210). The unit of the coordinates is given in the
UNIT
variable in the parameter file. The absolute
coordinates of the other bounding boxes will be automatically
determined, when all boxes are joined together (see
below).
5.3.1.2 Joining multiple bounding boxes together in block MJOIN
After specifying the dimension, segmentation and cell numbers of all
bounding boxes, they have to be joined together in order to form a
contiguous simulation volume. For this purpose, the block
MJOIN
is used in the parameter file as shown in the
following lines.
block MJOIN: {"x", -160, -160);
Left.join(Tube1, "x", 0, 0);
Tube1.join(Tube2, "x", 160, 160);
Tube2.join(Right, ; }
As can be seen, adjacent boxes are connected via the
join
command, whereof the syntax is
box1.join(box2, "side", da, db);
. The meaning of the
parameters is as follows:
box1, box2
= The names of the boxes which are to be joined.side
= The direction of the junction, can be either"x"
,"y"
or"z"
.da, db
= If the junction areas of both boxes are different, a coordinate offset is needed to specify the lateral coordinate offset between both boxes.
For the offsets, the coordinates a
and b
are the remaining axes which are vertical to the junction direction. The
offsets are obtained by drawing a vector from the lower edge of the
second box to the lower edge of the first box. Example: If the
coordinates of the lower edges are (x1, y1, z1) for the first box and
(x2, y2, z2) for the second box, we can obtain the offsets according to
Table 5.2.
Join direction | Direction of a | Direction of b | Offset da | Offset db |
---|---|---|---|---|
x | y | z | y1-y2 | z1-z2 |
y | x | z | x1-x2 | z1-z2 |
z | x | y | x1-x2 | y1-y2 |
Generally, if the lower edge of the second box is positioned within the surface of the first box, the offset is negative. Otherwise it is positive. Examples are given in the following figure:
Second box fits within first box surface | First box fits within second box surface |
---|---|
In case if all boxes are joined central-symmetrically, there is no
need to calculate the offsets manually. Instead, the command
box1.join_center(box2, "side");
can be used, where the
offset parameters are omitted. Our tube example has central-symmetric
junctions; thus, the MJOIN
block can be simplified to
block MJOIN: {"x");
Left.join_center(Tube1, "x");
Tube1.join_center(Tube2, "x");
Tube2.join_center(Right, ; }
5.3.1.3 Former join_rect function
In picmc versions prior 2.4.0.5 only the join_rect
function was available within the block MJOIN
. For
backwards compatibility reasons, this function can still be used;
however, for new simulation cases it is recommended to use the join command instead.
With the join_rect
command, our MJOIN
block
would look like follows:
block MJOIN: {
Tube1.join_rect(Tube2, "x", 0, 0, 105, 105, 0, 0);
Tube2.join_rect(Right, "x", 0, 0, 105, 105, 160, 160);
Left.join_rect(Tube1, "x", 160, 160, 260, 260, 0, 0);
};
Generally, the syntax of the join_rect
command between
two boxes - "FirstBox" and "SecondBox" - is as follows:
FirstBox.join_rect(SecondBox, "direction", a0, b0, a1, b1, delta_a, delta_b);
with the following parameters:
FirstBox
,SecondBox
= The names of the bounding boxes as defined in blockMQUAD
. The box with the lower coordinates (with respect to the direction of the junction) should always be specified first asFirstBox
, the other box with higher coordinates will be theSecondBox
.- For
"direction"
chose either"x", "y", "z"
. This is the direction of the junction. - The parameters a0,b0,a1,b1 specify the position of the rectangular orifice in relative coordinates with respect to the first box.
- The parameters delta_a, delta_b specify the offset of the rectangular orifice in relative coordinates with respect to the second box.
5.3.1.4 Examples of the join_rect function
The usage of the join_rect
command is demonstrated on
two examples. In the first example, a big bounding box is joined to a
smaller bounding box along X direction:
Here, the first box is joined to the second box via the hatched rectangular area. Since the hatched area is smaller than the full right side of the first box, the position of the rectangular area within the right side can be freely chosen. It spans between the coordinates (a0, b0)-->(a1, b1). In the given coordinate system, the connection is performed along X direction; then the a coordinate is parallel to Y, and the b coordinate is parallel to Z. The coordinate conventions are shown in the table in the previous section.
It has to be mentioned again, that the (a0, b0)-->(a1, b1) coordinates are relative with respect to the first box, i. e.
\(0 \le a0, a1 \le \Delta Y\textrm{ and }0 \le b0, b1 \le \Delta Z\)
Since the hatched junction surface is identical with the left surface of the smaller right box, the offset (delta_a, delta_b) is (0, 0) in this case.
The opposite situation occurs, if a smaller bounding box is joined to a bigger one:
Here, the hatched junction surface completely fills out the right surface of the first bounding box. Thus, the surface coordinates are
\((a0, b0) = (0,0)\)
\((a1, b1) = (\Delta Y, \Delta Z)\).
Now, the position of the junction surface within the left surface of the second box can be freely adjusted via the offset coordinates (delta_a, delta_b). Thus, the full command would read like
FirstBox.join_rect(SecondBox, "x", 0, 0, Delta_X, Delta_Y, delta_a, delta_b);