Open Access is an initiative that aims to make scientific research freely available to all. To date our community has made over 100 million downloads. It’s based on principles of collaboration, unobstructed discovery, and, most importantly, scientific progression. As PhD students, we found it difficult to access the research we needed, so we decided to create a new Open Access publisher that levels the playing field for scientists across the world. How? By making research easy to access, and puts the academic needs of the researchers before the business interests of publishers.
We are a community of more than 103,000 authors and editors from 3,291 institutions spanning 160 countries, including Nobel Prize winners and some of the world’s most-cited researchers. Publishing on IntechOpen allows authors to earn citations and find new collaborators, meaning more people see your work not only from your own field of study, but from other related fields too.
To purchase hard copies of this book, please contact the representative in India:
CBS Publishers & Distributors Pvt. Ltd.
www.cbspd.com
|
customercare@cbspd.com
Cartesian Cell Registration (CCR) is an effective method to reduce computational time for particle search and positioning, for example, in the Direct Simulation Monte Carlo (DSMC) method for dilute fluid flows and multiphase flow calculations. In this chapter, an efficient FORTRAN algorithm for the CCR method is presented to further reduce both the computational time for registration and computer memory. With this algorithm, the computation time for searching the target cell where the target particle exists is reduced by as much as 44,234 times. Moreover, this algorithm was successfully applied to the DSMC simulation, resulting in a 369-fold reduction in computational time compared to the brute-force approach, that is, searching all target cells until the cell in which the target particle is present is found.
*Address all correspondence to: ogami@se.ritsumei.ac.jp
1. Introduction
Particle search and positioning are the main time-consuming processes in the direct simulation Monte Carlo (DSMC) method [1, 2] for dilute fluid flows and the calculation of multiphase flows [3, 4, 5]. When a Cartesian (rectangular) grid system is used for these simulations, particle search and positioning are quite straightforward because the identification number (ID number) associated with the target cell in which a target particle is located can be calculated using a simple algebraic expression. However, structured and unstructured grid systems require a time-consuming process for searching a large number of target cells to position the target particles, although these grid systems can simulate flows around complex geometries with better quality of body fitting and better accuracy of object boundary calculation than the Cartesian grid system [6].
Taking advantage of the Cartesian grid system to efficiently search for target cells, Liang et al. [7] proposed a DSMC method that combines two levels of Cartesian grid systems with an unstructured triangular grid system. However, as Wang et al. [8] point out that although this method improves the calculation efficiency, multiple calculations must be performed separately at the junction of the two types of grids, which can be relatively tedious. Wang C et al. [9] divided the computational domain into rectangular regions for Cartesian grid systems by using a small number of triangular elements. The trajectory of the target particle is first tracked by searching the rectangular regions, followed by the triangular elements. This method also improves the calculation efficiency. However, the target grid cell system is limited to a triangular grid system, and searching for target grids in small areas can be time-consuming.
Wang Z et al. [8] proposed a Background Cartesian Grid Positioning (BCP) technique for two-dimensional computations. In this technique, a structured grid system is superposed on a layer of a background Cartesian grid system, and the geometric relationship between the structured and background grids is determined in advance to limit the number of target cells to be searched later. Although the results showed a significant reduction in the computational time required for particle positioning, the BCP was limited to structured grid systems as target grid systems. In addition, applying BCP to three-dimensional grids can be challenging.
To overcome these limitations, Ogami [10] introduced two methods to limit the target cells to search for the positioning of a target particle in a structured or unstructured grid system. The first is the Surrounding Cell Registration (SCR) method. In this method, the ID numbers of the target cells surrounding the central target cell are registered in advance with the central target cell to limit the cells to be searched for later. It has been demonstrated that a program using this method runs faster than other methods, as SCR has the smallest number of cells to search among the methods studied. However, SCR has two disadvantages: the ID numbers of the target cells into which the particles are introduced must be known in advance or predetermined, and it cannot be applied to particles that move through more than two cells during one-time step.
These drawbacks are overcome by the second method, called Cartesian cell registration (CCR). CCR uses a Cartesian grid system as a background two-dimensional grid system, as in the method proposed by Wang et al. [8]. However, this method can be applied in a completely different way to both structured and unstructured grids as target grids. Using SCR and CCR for a two-dimensional potential flow around a circular cylinder as an application of the proposed methods, it is shown that computational programs using both methods are almost 200 times faster than a program using the brute-force approach (searching all target cells until the target particle is found).
The ideal condition for CCR is that only one target cell is registered in each cell of the Cartesian grid system. To achieve this, the cell size of the Cartesian grid system must be as small as the smallest cell of the target grid system, in which the flow is to be simulated. However, as the number of cells of the Cartesian grid system increases (while the cell size of the Cartesian grid system decreases), the computer time for registration and the use of computer memory increase immeasurably. To solve this problem, we propose an efficient computer algorithm for registration. Using this method, simulation results of the DSMC method for dilute flows around two-dimensional circular cylinders are presented. It was found that the program for searching target cells (finding the ID number of a target cell) in which target particles exist runs 44,234 times faster, and the entire program, including the search program and DSMC simulation, runs 448 times faster than the brute-force approach.
2. CCR method and the problem of previous FORTRAN algorithm
In this section, we explain the problem of the FORTRAN program for Cartesian cell registration (CCR) method [10]. Consider a grid system around a circular cylinder as the target grid system (shown in red in Figure 1), which is to be registered with a Cartesian grid system (shown in black in Figure 1). The target grid system can be structured, unstructured, or have a polygonal structure that is not a quadrangular structure. Also, the ID placed on the target cells were either sequential or arbitrary.
Figure 1.
Target grid cell system (red) and Cartesian grid cell system (black).
Consider a square cell S0 of width W drawn in bold black lines on the Cartesian grid system, as shown in Figure 2. The ID numbers of all red cells in the target grid system whose centers are within a distance R from the center of S0 are registered with the square cell S0. This registration process was performed for all Cartesian and target red cells. Some red target cells may have been registered more than once with different Cartesian grid cells. However, it is very important to adjust parameters W and R so that no red target cells remain unregistered. The number of registered target cells may be zero or excessive, depending on the W and R parameters and the respective cell size of the target grid system.
Figure 2.
Target red cells in the blue circle are registered with the black Cartesian grid cell S0.
The core components of the FORTRAN algorithm for CCR are shown in Appendix I (the author is familiar with the old FORTRAN).
Suppose the x and y coordinates of the target particle are xp and yp, respectively. The ID number of the target cell in which this particle is located can be determined using the program in Appendix II.
The above program can be used to determine the ID number of the target cell. Using this algorithm, the computation time of potential flows around a circular cylinder decreases as the width of the Cartesian cell decreases (i.e., as the average number of target cells registered in the Cartesian cells decreases), resulting in a maximum speed up ratio of 200 [10]. The average number of fastest computations was 14.859, and if this number can be reduced to one, the speed up ratio may be a limiting factor. However, to reduce this average number, the computation time for registering the above algorithm increases immeasurably, as does the computer memory. Moreover, the number of elements in the dimensional array iSR(200,200,400) of the above program increases, even though much of the memory in the coarse cell regions of the target grid system is unused.
In the next section, a new algorithm is presented that overcomes these drawbacks.
The most effective change to the FORTRAN algorithm is to change the order of the triple-do loops. In the previous algorithm, the do loop for the target cells (do k = 1, ic) was the innermost loop. In the modified algorithm, this loop is placed as the outermost loop, as shown in the program in Appendix III.
See how the above algorithm works. Figure 3 shows an example of the target grid system of a circular cylinder (red) and a Cartesian grid system (black). The Cartesian cell (i = 9, j = 7) is shown enlarged on the right side of the figure. This cell covers or partially touches the target cell with ID numbers 13–16, 99–102, and 326–329.
Figure 3.
Target grid cell system (red) and Cartesian grid cell system (black). The Cartesian cell 9–7 is magnified on the right side.
Running the above algorithm gives the elements of the array iSR, as shown in Table 1 (only the elements of Cartesian cell 9–7 are shown). The elements of parameters ir and k are in ascending order. However, the parameters i and j are not. To use this data more efficiently to search for the ID number of the target cell in which a target particle exists, the elements in the array iSR are rearranged using the algorithm in Appendix IV.
ir
…
205
206
207
208
209
210
…
iSR(ir,1) = k
…
13
14
15
16
16
17
…
iSR(ir,2) = i
…
9
9
9
9
10
10
…
iSR(ir,3) = j
…
7
7
7
7
7
7
…
400
401
402
403
404
…
99
100
101
102
102
…
9
9
9
9
10
…
7
7
7
7
7
…
Table 1.
Elements of array iSR.
Table 2 lists the elements in the array iRec from i = 101 to i = 112, that were determined using the above algorithm. These elements are the ID numbers of the target cells (red cells in Figure 3 on the right) in the Cartesian cell 9–7 (black square cell in Figure 3 on the right), where the numbers i for the first and last cells are given as index (9,7,1) = 101 and index (9,7,2) = 112, respectively (Table 3).
i
…
101
102
103
104
105
iRec(i)
…
13
14
15
16
99
106
107
108
109
110
111
112
…
100
101
102
326
327
328
329
…
Table 2.
Elements of the array iRec (i) from i = 101 to 121.
index(9, 7, 1)
101
index(9, 7, 2)
112
Table 3.
Elements of the array index for the Cartesian cell 9–7.
Only the elements of the Cartesian cell 9–7 are shown.
Only, the elements of the Cartesian cell 9–7 are shown.
If the x and y coordinates of the target particle are xp and yp, respectively, the ID number of the target cell in which this particle exists can be easily determined using the algorithm in Appendix V.
The subroutines “INorOUT” and “intersect” are shown in Appendices VI and VII. Whether a line AB and a line CD intersect and whether a point (x, y, x) exists in a square ABCD can be easily judged by comparing areas, as shown in the Appendices.
In a previous study [10], the CCR method was not used for the DSMC simulation, but for simple potential flows. In this section, dilute flows around a circular cylinder are simulated with DSMC using the new FORTRAN algorithm shown in Section 3. We consider two-dimensional flows of Argon around a circular cylinder with a radius 0.05 m and an inlet velocity of 2634.1 m/s. The parameters of the target and Cartesian grid cell systems are listed in Tables 4 and 5, respectively. The parameters for the dilute gas and DSMC simulations are listed in Tables 6 and 7, respectively. To compare the results, the same parameters as in [11] were used for the gas and DSMC simulations. The full FORTRAN programs were written with double precision, compiled using the commercial compiler PGI 19.10–0 with an optimization level of four, and executed on a computer with an Intel Core i7-4770K processor with a clock speed of 3.50 GHz.
Radius of the circular cylinder
0.05 m
Number of divisions in the radial direction
260
Number of divisions in the circumferential direction
462
Total grid cell count
120,120
Table 4.
Parameters of target grid cell system.
Maximum x-coordinate
0.654 m
Minimum x-coordinate
−0.396 m
Maximum y-coordinate
0.525 m
Minimum y-coordinate
−0.525 m
Number of cells in the x-direction
1 ∼ 2000
Number of cells in the y-direction
1 ∼ 2000
Table 5.
Parameters of Cartesian grid cell system.
Temperature of the cylinder surface
500 K
Stream temperature
200 K
Inlet velocity
2634.1 m/s
Number density
4.247×1020
Molecular diameter
3.595 × 10−10 m
Molecular mass of Ar
6.63×10−26kg
Table 6.
Parameters of dilute gas.
Time step
0.5×10−6s
Total updates
40,000
Table 7.
Parameters for the DSMC simulation.
Figures 4 and 5 show the converged results for Mach number and temperature. Because the solutions were vertically symmetric, only the upper half is shown. These solutions agree well with the results in [11] (pp. 219 and 220).
Figure 4.
Mach number. Only the upper half is displayed.
Figure 5.
Temperature. Only the upper half is displayed.
Figure 6 shows the CPU time for registration (red) and that for DSMC simulations (black) with different numbers of N (N × N Cartesian grid cell systems). For example, if a 1 × 1 Cartesian grid cell system, i.e., the brute-force approach, is used, the CPU time for registration is 0.42 s and that for DSMC simulation for 100 time steps is 165282.07 s. When the 2000 × 2000 Cartesian grid cell system was used, the CPU time for registration was 13535.78 s and that for DSMC simulation was 447.79 s. The DSMC speed up ratio was calculated to be 165282.07/447.79≈369.11. Table 8 lists the CPU times and speed up ratios for N = 1–2000.
Figure 6.
Registration time and DSMC time for 100 time steps.
N
Registration CPU time
DSMC CPU time
Speed up ratio
1
0.42
165282.07
1.00
2
0.57
46141.01
3.58
4
0.62
14822.28
11.15
10
0.51
3015.89
54.80
50
0.91
602.53
274.31
100
1.91
490.39
337.04
500
124.77
457.20
361.51
1000
1194.68
457.33
361.41
2000
13535.78
447.79
369.11
Table 8.
Registration time, DSMC time, and speed up ratio of DSMC calculation.
To study the speed up ratio for searching for target cells where the target particles exist, it is impractical to measure the CPU time for the search process because this process is mixed in the DSMC simulation. In addition, measuring the CPU time for all particles would affect total CPU time.
Instead, the average number of cells searched for one particle during one-time step was measured, and is summarized in Table 9 and Figures 7 and 8. For example, if a 1 × 1 Cartesian grid cell system is used (the brute-force approach), this number is calculated to be 58208.83, which is nearly half of the total grid cell number of 120,120. Notably, this implies that the probability of any cell being searched is 50%.
N
Averaged searched cell number for one particle during one-time step
Speed up ratio
1
58208.83
1.00
10
893.87
65.11
100
15.86
3670.44
300
3.95
14740.89
500
2.53
22989.85
800
1.87
31082.32
1000
1.64
35445.57
2000
1.32
44234.29
Table 9.
Averaged searched cell number for one particle during one-time step, and speed up ratio.
Figure 7.
Averaged searched cell number for one particle during one-time step.
Figure 8.
Speed up ratio when searching cells.
As the number of cells N of the Cartesian grid cell system increased, the number of searched cells decreases dramatically (Table 9 and Figure 7), and the speed up ratio increased exponentially (Table 9 and Figure 8). Most importantly, the CPU time for the DSMS reached a plateau at nearly N = 100 (Figure 6), indicating that N = 100 is sufficient for our simulation, although the speed up ratio for searching the target cells in Figure 8 did not reach a plateau. This means that when N is more than 100, the CPU time for searching the target cells is negligible compared to the CPU time for the entire DSMC simulation. Consequently, the efficiency of our fast particle search algorithm reaches a maximum.
Moreover, as explained in Section 2, for the previous method, the average number of target cells registered in Cartesian cells for the fastest simulation was 14.859, but the average number of modified algorithms for N = 2000 was 1.59, excluding Cartesian cells with no registered target cells. This means that the new algorithm can reduce this average number (i.e., the size of the Cartesian cells can be reduced, and the Cartesian cell number can be increased) without sacrificing computational time or computer memory.
5. Application for the measurement of multiple physical quantities
It was mentioned in the Introduction that particle searching and positioning are the most time-consuming processes in the direct simulation Monte Carlo (DSMC) method [1, 2] for dilute fluid flows, and in multiphase flow calculations [3, 4, 5]. As shown in Section 4, Cartesian cell registration (CCR) is an effective method for reducing the computational time required for DSMC calculations. In this section, we briefly explain how the CCR method can be used in fields of computational calculations other than DSMC.
We have worked on computational simulations for measuring multiple physical quantities using a single-motion sensor system [12]. Figure 9a shows a simple model of a motion sensor system consisting of two temperature sensors and a heater placed in between. When this motion sensor was accelerated in the right direction, the top of the temperature distribution was inclined in the right direction owing to the buoyancy effect, as shown in Figure 9b. Analysis of a sufficient amount of data obtained by computer simulations showed that the magnitude of the acceleration is related to the temperature and temperature variations detected by the temperature sensors. Therefore, a motion sensor equipped with this data can measure acceleration due to temperature and the variation in temperature detected by the temperature sensors.
Figure 9.
Motion sensor system consisting of two temperature sensors and a heater placed in between (a) and temperature distributions with and without acceleration (b).
When this motion sensor is accelerated and rotated at the same time, the signal detected by the temperature sensors contains mixed information about both acceleration and rotation speed. The problem is to separate this signal into acceleration and rotation speeds. Using computer simulations, the two temperature outputs on the two temperature sensors can be calculated from the two inputs for acceleration and rotational speed. For example, when the acceleration is 10 m/s2 and the rotational speed is 250°/s (black point A in Figure 10a), two temperature sensors in this motion sensor system detect 477.9 and 407.1 K (black point A’ in (b)). Similarly, the rectangular region ABCD in (a) is converted to quadrilateral A’B’C’D′ in (b). The remaining rectangular regions in (a) have been transformed into quadrilaterals, as shown in (b).
Figure 10.
Inputs of acceleration and rotation (a) and outputs of temperature (b). The region of input ABCD in (a) is transformed into the region of output A’B’C’D′ in (b).
When a real motion sensor system installed in cars, drones, phones, game consoles, and virtual reality systems, for example, detects temperatures of 490 and 382 K with the simulation data in Figure 9, as indicated by the red circle in Figure 10b, acceleration and rotation speed can be estimated by inverting this red point in (b) to a red point in (a). The CCR method can be used to determine the quadrilateral A’B’C’D′, in which the red point in (b) lies. That is, the red quadrilaterals in (b) are registered with black Cartesian cells also in (b). Then, the black Cartesian cell in which the red point exists can be calculated algebraically, and the quadrilateral A’B’C’D′ is determined by the CCR method. Finally, even though the signal detected by the motion sensor contains mixed information about acceleration and rotation speed, as mentioned above, the corresponding rectangle ABCD in (a), where the inversely transformed red point is located, can be determined, and the acceleration and rotational speed can be easily calculated by our fast search algorithm.
In this study, an efficient FORTRAN algorithm for the Cartesian cell registration (CCR) method is presented to further reduce both the computation time for registration and computer memory. Using this algorithm, the computational time for searching the target cell where the target particle is located was reduced by as much as 44,234 times. In addition, this algorithm was successfully applied to the DSMC simulation, resulting in a 369-fold reduction in computational time compared to the brute-force approach (searching all target cells until the target particle is present).
Future work will apply the new CCR method to multiphase flow calculations [3, 4, 5], the detection of static particle positions [13], and measurement of multiple physical quantities [12], as presented in Section 5. When measuring multiple physical quantities, the set of inputs and outputs must be a one-to-one correspondence. Removing this limitation will be our main task. We will also increase the number of simultaneously measured physical quantities from two to more and write our algorithm in Python, one of the most popular languages in engineering today, to achieve even faster processing.
! The following arrays are x and y coordinates of four corners of target cell with ID number k. ! cell(k,1,1), cell(k,1,2) ! cell(k,2,1), cell(k,2,2) ! cell(k,3,1), cell(k,3,2) ! cell(k,4,1), cell(k,4,2) ! The maximum cell number 40000 is tentative. ! The order of the four corners is clockwise or counterclockwise.
! The following arrays are x and y coordinates of the centers of Cartesian cell with ID number ! i and j. ! recent(i,j,1), recent(i,j,2) ! The maximum ID number 200 is tentative. ! i is the ID number of the Cartesian cell in x direction ! j is the ID number of the Cartesian cell in y direction
! iSR(i,j,1) is the number of registered target cells in the Cartesian cell ij. ! iSR(i,j,2 ∼ 400) are ID numbers of the target cells which are registered in the Cartesian cell ij. ! 400 (tentative) is the maximum number of the target cells minus 1 which Cartesian cell ij ! can memorize. ! Initialize iSR(i,j,1); the number of registered target cells in the Cartesian cell ij. do i = 1,nRecx! nRecx is the maximum number of i do j = 1,nRecy! nRecy is the maximum number of j iSR(i,j,1) = 0 enddo enddo ! Start registration do i = 1,nRecx do j = 1,nRecy do k = 1,ic ! do loop for the target cells ! ic is the total number of the target cells ! The center of the target cell xjc = (cell(k,1,1) + cell(k,2,1) + cell(k,3,1) + cell(k,4,1))/4 yjc = (cell(k,1,2) + cell(k,2,2) + cell(k,3,2) + cell(k,4,2))/4
! The center of the rectangle cell xc = recent (i,j,1) yc = recent (i,j,2) ! The distance between the centers al = dsqrt((xc-xjc)**2 + (yc-yjc)**2) if(al.le.ak*dx) then! ak is arbitrary number. dx is width of Cartesian cell iii = iSR(i,j,1) iii = iii + 1 iSR(i,j,1) = iii! iii is the total number of registered target cells in the Cartesian cell ij iSR(i,j,iii + 1) = k! k is the ID number of the target cell endif enddo enddo enddo
! ID numbers of the Cartesian cell i and j i = (xp-xRecmin)/dx + 1! xRecmin is the minimum x coordinate of the Cartesian grid system j = (yp-yRecmin)/dx + 1! yRecmin is the minimum y coordinate of the Cartesian grid system
! The total number of target cells registered in Cartesian cell ij iscn = iSR(i,j,1)
! Find out the target cell where the target particle exists do k = 1,iscn ii = iSR(i,j,k + 1)! candidate ID number of the target cell x1 = cell(ii,1,1)! x coordinate of the first corner on the target cell y1 = cell(ii,1,2)! y coordinate of the first corner on the target cell x2 = cell(ii,2,1)! x coordinate of the second corner on the target cell y2 = cell(ii,2,2)! y coordinate of the second corner on the target cell x3 = cell(ii,3,1)! x coordinate of the third corner on the target cell y3 = cell(ii,3,2)! y coordinate of the third corner on the target cell x4 = cell(ii,4,1)! x coordinate of the fourth corner on the target cell y4 = cell(ii,4,2)! y coordinate of the fourth corner on the target cell
call INorOUT(x1,y1,x2,y2,x3,y3,x4,y4,xp,yp,iflag)! subroutine to investigate ! whether the target particle exists in this target cell. ! This is shown in Appendix VI. if(iflag.eq.0) then! the target particle exists in this target cell goto 999 endif enddo
999 continue! the ID number (the element of parameter k) has been found.
Dimension iSR(100,000,3), rec(100,100,4,2), cell(100,000,4,2) Dimension index(100,100,2), iRec(100000) ! The element numbers 100,000 and 100 are tentative. These depend on the total cell number ! of the target grid system and the Cartesian grid system
! Initialize iSR
do i = 1,100,000 iSR(i,1) = 0 iSR(i,2) = 0 iSR(i,3) = 0 enddo
ir = 0! ir is the total number of registered set of k, i and j
! The following arrays are x and y coordinates of four corners of target cell with ID number k. ! cell(k,1,1), cell(k,1,2)! the first corner ! cell(k,2,1), cell(k,2,2)! the second corner ! cell(k,3,1), cell(k,3,2)! the third corner ! cell(k,4,1), cell(k,4,2)! the forth corner
! do loop for the target cells is the outermost loop do k = 1,isc0 ! do loop for the target cells ! isc0 is the total number of the target cells
! The maximum and minimum values of the four corners of the target cells xmax2 = −1.0d10 if(xmax2.lt.cell(k,1,1)) xmax2 = cell(k,1,1) if(xmax2.lt.cell(k,2,1)) xmax2 = cell(k,2,1) if(xmax2.lt.cell(k,3,1)) xmax2 = cell(k,3,1) if(xmax2.lt.cell(k,4,1)) xmax2 = cell(k,4,1)
! Find out in which Cartesian cell these corners exist
is = (xmin2-xRecmin)/dx + 1! xRecmin is the minimum x coordinate of the Cartesian grid ie = (xmax2-xRecmin)/dx + 1
js = (ymin2-yRecmin)/dx + 1! yRecmin is the minimum y coordinate of the Cartesian grid je = (ymax2-yRecmin)/dx + 1 ! Start registration ! do loops for the Cartesian cells are inner loops do i = is,ie do j = js,je
do ii = 1,4 Ax = cell(k,ii,1) Ay = cell(k,ii,2) iii = ii + 1 if(ii.eq.4) iii = 1 Bx = cell(k,iii,1) By = cell(k,iii,2) ! The following arrays are x and y coordinates of the four corners of Cartesian cell ij. ! rec(i,j,1,1), rec(i,j,1,2)! the first corner ! rec(i,j,2,1), rec(i,j,2,2)! the second corner ! rec(i,j,3,1), rec(i,j,3,2)! the third corner ! rec(i,j,4,1), rec(i,j,4,2)! the fourth corner
do jj = 1,4 Cx = rec(i,j,jj,1) Cy = rec(i,j,jj,2) jjj = jj + 1 if(jj.eq.4) jjj = 1 Ex = rec(i,j,jjj,1) Ey = rec(i,j,jjj,2)
! See if the side of the target cell (AB) and the side of the Cartesian cell (CE) intersect
call intersect(Ax,Ay,Bx,By,Cx,Cy,Ex,Ey,iflag,err)! this subroutine is shown in Appendix VII.
if(iflag.eq.0) goto 555! if they intersect, go to 555 ! and register the target cell k with the Cartesian cell ij enddo enddo
! if no sides of the target cell (AB) and of the Cartesian cell (CD) intersect, ! then see if the whole region of target cell k is inside the Cartesian cell ij.
ifs = 0 do lll = 1,4 x = Scell(k,lll,1) y = Scell(k,lll,2) call INorOUT(xA,yA,xB,yB,xC,yC,xD,yD,x,y,iflag) ifs = ifs + iflag enddo if(ifs.eq.0) goto 555! if the whole region of target cell k is inside the Cartesian cell ij, go to 555 ! and register the target cell k with the Cartesian cell ij ! Or inversely see if the whole region of the Cartesian cell ij is inside the target cell k
Ax = cell(k,1,1) Ay = cell(k,1,2) Bx = cell(k,2,1) By = cell(k,2,2) Cx = cell(k,3,1) Cy = cell(k,3,2) Ex = cell(k,4,1) Ey = cell(k,4,2)
ifs = 0 do lll = 1,4 x = rec(i,j,lll,1) y = rec(i,j,lll,2) call INorOUT(Ax,Ay,Bx,By,Cx,Cy,Ex,Ey,x,y,iflag)
ifs = ifs + iflag
enddo
if(ifs.eq.0) goto 555! if the whole region of the Cartesian cell ij is inside the target cell k, ! go to 555 and register the target cell k with the Cartesian cell ij
goto 666! if there is nothing to register, go to 666.
555 continue
ir = ir + 1
iSR(ir,1) = k! id number of the target cell iSR(ir,2) = i! id number of the Cartesian cell in x direction iSR(ir,3) = j! id number of the Cartesian cell in y direction
666 continue
enddo enddo enddo
! end of registration
Appendix IV
! Rearrange the ID number of target cells in the order of Cartesian cell i and j in = 0
do i = 1,nRecx! nRecx is the maximum number of i do j = 1,nRecy! nRecy is the maximum number of j is = in do k = 1,ir! ir is the total number of registered sets of k, i and j
if(iSR(k,2).eq.i.and.iSR(k,3).eq.j) then in = in+1 iRec(in) = iSR(k,1)! the ID number of the target cell k registered in Cartesian cell ij is ! put into array iRec of ID number in endif enddo ! ID numbers of target cells registered in Cartesian cell ij starts from index(i,j,1)-th element ! to index(i,j,2)-th element in arrary iRec(in). Namely, in = index(i,j,1) ∼ index(i,j,2)
index(i,j,1) = is+1 index(i,j,2) = in
enddo enddo ! End of rearrangement
Appendix V
! ID numbers of the Cartesian cell i and j i = (xp-xRecmin)/dx + 1! xRecmin is the minimum x coordinate of the Cartesian grid j = (yp-yRecmin)/dx + 1! yRecmin is the minimum y coordinate of the Cartesian grid
is = index(i,j,1) ie = index(i,j,2)
! Find out the target cell where the target particle exists do k = is, ie ii = iRec(k)! ID number of target cell registered in Cartesian cell ij. x1 = cell(ii,1,1)! x coordinate of the first corner on the target cell y1 = cell(ii,1,2)! y coordinate of the first corner on the target cell x2 = cell(ii,2,1)! x coordinate of the second corner on the target cell y2 = cell(ii,2,2)! y coordinate of the second corner on the target cell x3 = cell(ii,3,1)! x coordinate of the third corner on the target cell y3 = cell(ii,3,2)! y coordinate of the third corner on the target cell x4 = cell(ii,4,1)! x coordinate of the fourth corner on the target cell y4 = cell(ii,4,2)! y coordinate of the fourth corner on the target cell
call INorOUT(x1,y1,x2,y2,x3,y3,x4,y4,xp,yp,iflag)! subroutine to investigate ! whether the target particle exists in this target cell if(iflag.eq.0) then! the target particle exists in this target cell goto 999 endif enddo
999 continue! the ID number is the element of parameter k
Appendix VI
********************************************************** subroutine INorOUT(xA,yA,xB,yB,xC,yC,xD,yD,x,y,iflag) ********************************************************** ! Check if a point (x,y,z) exists in a square ABCD ! The order of the four corners ABCD is clockwise or counterclockwise. ! Just comparing areas shows if a point (x,y,z) exists in a square ABCD or not.
iflag = 1! false e1 = dabs(er/area) if(e1.lt.2.5d-9) iflag = 0! The threshold value should be adjusted depending on the cell size if(area.eq.0.0d0) iflag = 1
return end ********************************************************** subroutine rectangle(xA,yA,xB,yB,xC,yC,xD,yD,area) ********************************************************** implicit real*8 (a-h,o-z)
***************************************************************** subroutine intersect(Ax,Ay,Bx,By,Cx,Cy,Dx,Dy,iflag,err) ***************************************************************** * Check if line AB and line CD intersect. * Just comparing areas shows if line AB and line CD intersect or not. * implicit real*8 (a-h,o-z) iflag = 1 call triangle(Ax,Ay,Cx,Cy,Dx,Dy,ACD) call triangle(Bx,By,Cx,Cy,Dx,Dy,BCD) call triangle(Cx,Cy,Ax,Ay,Bx,By,CAB) call triangle(Dx,Dy,Ax,Ay,Bx,By,DAB)
if(err.lt.1.0d-10) iflag = 0! The threshold value should be adjusted depending on the cell size return end ***************************************************************** subroutine triangle(xA,yA,xB,yB,xC,yC,area) ***************************************************************** implicit real*8 (a-h,o-z)
! Area of triangle ABC area = dabs((xA-xB)*(yC-yA)-(xC-xA)*(yA-yB))/2
return end
References
1.Bird GA. Molecular Gas Dynamics and the Direct Simulation of Gas Flows. Oxford: Oxford University Press; 1994
2.Bird GA. Direct simulation and the Boltzmann equation. Physics of Fluids. 1970;13:2676-2681
3.Shojaee S, Hosseini SH, Razavi BS. Computational fluid dynamics simulation of multiphase flow in structured packings. Journal of Applied Mathematics. 2012;2012:1-17
4.Florice NM, Andrei K. Modelling and simulation of multiphase flow applicable to processes in oil and gas industry. Chemical Product and Process Modeling. 2019;20170066:1-16. DOI: 10.1515/cppm-2017-0066
5.Parsi M, Kara M, Agrawal M, Kesana N, Jatale A, et al. CFD simulation of sand particle erosion under multiphase flow conditions. Wear. 2017;376-377(B):1176-1184
6.Blazek J. Principles of grid generation. In: Computational Fluid Dynamics: Principles and Applications. 3rd ed. Amsterdam: Elsevier; 2015
7.Liang J, Yan C, Du BQ. An algorithm study of three dimensional DSMC simulation based on two-level Cartesian coordinates grid structure. Acta Aerodynamica Sinica. 2010;28:466-471
8.Wang Z, Li L, Zhang B, Liu HBCP. Particle positioning techniques for DSMC method. Journal of Aeronautics, Astronautics and Aviation. 2019;51:225-236
9.Wang C, Cheng J, Ji L, Lu Y, Sun Y, et al. 2-D DSMC algorithm based on Delaunay triangles. Journal of Tsinghua University (Science and Technology). 2015;10:1079-1086. DOI: 10.16511/j.cnki.qhdxxb.2015.22.010
10.Ogami Y. Fast algorithms for particle searching and positioning by cell registration and area comparison. Trends in Computer Science and Information Technology. 2021;6:7-16
11.Bird GA. The DSMC Method. Version 1.2. 2013. Available from: http://www.gab.com.au/index.html.
12.Siddique K, Ogami Y. Computational study on thermal motion sensors that can measure acceleration and rotation simultaneously. Sensors. 2022;22(18):6744. DOI: 10.3390/s22186744 [Accessed: April 17, 2023]
13.Thapa S, Lukut N, Selhuber-Unkel C, Cherstvy AG, Metzler R. Transient superdiffusion of polydisperse vacuoles in highly motile amoeboid cells. The Journal of Chemical Physics. 2019;150:144901-1-144901-18
Written By
Yoshifumi Ogami and Siddique Kamran
Submitted: 26 May 2023Reviewed: 02 June 2023Published: 20 June 2023
Open access peer-reviewed chapter
