Home A time-step-robust algorithm to compute particle trajectories in 3-D unstructured meshes for Lagrangian stochastic methods
Article
Licensed
Unlicensed Requires Authentication

A time-step-robust algorithm to compute particle trajectories in 3-D unstructured meshes for Lagrangian stochastic methods

  • Guilhem Balvet ORCID logo EMAIL logo , Jean-Pierre Minier ORCID logo , Christophe Henry ORCID logo , Yelva Roustan ORCID logo and Martin Ferrand ORCID logo
Published/Copyright: March 15, 2023

Abstract

The purpose of this paper is to propose a time-step-robust cell-to-cell integration of particle trajectories in 3-D unstructured meshes in particle/mesh Lagrangian stochastic methods. The main idea is to dynamically update the mean fields used in the time integration by splitting, for each particle, the time step into sub-steps such that each of these sub-steps corresponds to particle cell residence times. This reduces the spatial discretization error. Given the stochastic nature of the models, a key aspect is to derive estimations of the residence times that do not anticipate the future of the Wiener process. To that effect, the new algorithm relies on a virtual particle, attached to each stochastic one, whose mean conditional behavior provides free-of-statistical-bias predictions of residence times. After consistency checks, this new algorithm is validated on two representative test cases: particle dispersion in a statistically uniform flow and particle dynamics in a non-uniform flow.

MSC 2010: 76M35

Award Identifier / Grant number: 2020/1387

Funding statement: G. Balvet has received a financial support by ANRT through the EDF-CIFRE contract number 2020/1387. The authors acknowledge the infrastructures at EDF R&D and the CEREA laboratory for providing access to computational resources.

A Details of the particle tracking algorithm for 3-D unstructured mesh

This appendix presents the trajectory algorithm that is used in the present study. This algorithm is able to track the motion of particles even in 3-D fully unstructured meshes with warped faces. The tracking algorithm is first described in Section A.1. It includes the description of the original tracking algorithm with the detection of face-crossing events in Section A.1.1. Then the algorithm is extended in Section A.1.2 to compute the location and exit time when a particle crosses a face. Last, the algorithm is validated in Section A.2 by comparing the numerical results obtained using various 3-D unstructured meshes in a simple non-uniform flow.

A.1 Principle of the neighbor search algorithm

The algorithm is based on a successive neighbor search. This means that the cell inside which a particle currently resides is determined by browsing through the neighboring cells. Such algorithms require three pieces of information: (a) the origin particle location X O = X n , (b) the corresponding cell inside which it was initially and (c) the destination particle location X D = X n + 1 (as depicted in Figure 4). The principle is then to determine if the particle leaves the current cell assuming a free-flight motion between point X O and X D . This is performed by

  1. computing which faces of the current cell are intersected by the line ( X O X D ) ;

  2. checking if the intersection belongs to the straight-line vector X O X D = X D X O .

The key issue is then to have a robust method to detect the intersection between a displacement vector and any face. It is of prime importance to prevent any particle from being permanently lost in the computational domain. For that reason, the method uses Boolean elementary tests which are reproducible from one cell to another so that it can handle pathological cases such as when the vector X O X D crosses a face through one of its edges (to the machine precision).

A.1.1 Method to detect face-crossing events through warped-faces

The method to detect face-crossing events is based on the decomposition of each face into a set of triangular sub-faces (see also Figure 17). Each triangular sub-face is built using one of the oriented edges of the face (formed by two consecutive vertices X i and X j ) and the center of gravity of the face X f . This decomposition of faces ensures that each triangular sub-face treated is planar, hence making the method tractable even for warped meshes (i.e., with faces whose vertices do not belong to the same plane).

Once a face is decomposed into a set of triangular sub-faces, the intersection between a vector and each planar sub-face is detected using a kind of Möller–Trumblemore algorithm [26]. The principle is to project the vertices of each sub-face in the oriented plane orthogonal to the displacement vector X O X D passing through X O denoted ( X O , X O X D ) . Then, to compute if the line ( X O X D ) crosses a sub-face, we simply need to check if the point X O belongs to the triangle formed by the projection of the sub-face on this plane (as displayed in Figure 16, where the superscript ( ) corresponds to the projection in the plane perpendicular to the displacement). For that purpose, we resort here to a series of three elementary Boolean tests, each one allowing to verify if the point X O is located on the “proper side” of a projected edge.

Elementary Boolean tests

For each of the three edges forming a projected sub-face (namely X f X i , X f X j and X i X j ), we have to verify whether the point X O is on the proper side of the projected edge. To that extent, we resort here to simple logical tests. For the sake of clarity, let us consider the case of an oriented edge connecting two points X α X β (where 𝛼 and 𝛽 are the indexes of two vertices of a sub-face) on the projection plane. In that case, the logical test L α , β edge reads

L α , β edge = { true if ( X α X β X α X O ) X O X D > 0 , false otherwise .

As displayed in Figure 14, this elementary Boolean test provides information on whether a point X O is in the half plane on the left (true) or on the right (false) of the projected line ( X α X β ) . It is worth noticing that, the inequality in this test being strict, the behavior on the oriented line is asymmetric. Indeed, if point X O belongs to this line at the machine precision, test L α , β edge returns false. In other words, we arbitrarily consider that the line ( X α X β ) belongs to the closed right half plane and not to the open left one. As a result, the Boolean test L α , β edge not ( L α , β edge ) is a partition of the domain. This means that the computation of L α , β edge is reproducible for a given displacement vector X O X D for two faces sharing the same edge[1] provided that the ordering of the vertices 𝛼 and 𝛽 is fixed. In the present work, we have imposed to sort the vertices in the following order: first X f , then X i and finally X j (with i < j ).

Figure 14 
                           Sketch illustrating how the elementary boolean test 
                                 
                                    
                                       
                                          L
                                          
                                             α
                                             ,
                                             β
                                          
                                          edge
                                       
                                    
                                    
                                    \mathcal{L}^{\mathrm{edge}}_{\alpha,\beta}
                                 
                               works: the point 
                                 
                                    
                                       
                                          (
                                          
                                             X
                                             0
                                          
                                       
                                    
                                    
                                    (\mathbf{X}_{0}
                                 
                               can either be located on the left-hand side or on the right-hand side of the oriented line 
                                 
                                    
                                       
                                          (
                                          
                                             
                                                X
                                                α
                                                †
                                             
                                             ⁢
                                             
                                                X
                                                β
                                                †
                                             
                                          
                                          )
                                       
                                    
                                    
                                    (\mathbf{X}_{\alpha}^{\dagger}\mathbf{X}_{\beta}^{\dagger})
                                 
                              .
Here, the figure is seen from above the plane 
                                 
                                    
                                       
                                          (
                                          
                                             X
                                             0
                                          
                                          ,
                                          
                                             
                                                X
                                                O
                                             
                                             ⁢
                                             
                                                X
                                                D
                                                ⟂
                                             
                                          
                                          )
                                       
                                    
                                    
                                    (\mathbf{X}_{0},\mathbf{X}_{O}\mathbf{X}_{D}^{\perp})
                                 
                              , meaning that the elementary Boolean test is true on the left-hand side of the figure.
Note that the line belongs to the closed half-plane (i.e., false in red color).
Figure 14

Sketch illustrating how the elementary boolean test L α , β edge works: the point ( X 0 can either be located on the left-hand side or on the right-hand side of the oriented line ( X α X β ) . Here, the figure is seen from above the plane ( X 0 , X O X D ) , meaning that the elementary Boolean test is true on the left-hand side of the figure. Note that the line belongs to the closed half-plane (i.e., false in red color).

Figure 14 also shows that this elementary Boolean test is not enough to determine if a point X O belongs to the projected triangle sub-face. In fact, a point can be located either on the right side or on the left side of the line depending on the orientation of the vector X i X j but also on whether the point is above or below the plane. While the orientation of the vector is now fixed thanks to the sorted vertices introduced in the previous paragraph, we have to introduce the notion of plane orientation to deal with the second issue. For that purpose, we rely on the orientation of the sub-face, and we introduce an additional logical test to know if the displacement vector X O X D is aligned with the oriented sub-face. It returns true if they are aligned and false if they are not. This condition reads

A f , i , j face = { true if ( X f X i X f X j ) X O X D > 0 , false otherwise .

By combining both logical tests ( L i , j edge A f , i , j face ) (where ≡ means that both Boolean tests have the same value), we are able to determine on which side of the line X i X j the particle lies with respect to the face orientation.

Combined elementary Boolean tests

To determine if the intersection point X i is within the projected face, we have then to combine these elementary logical tests. As displayed in Figure 15, the intersection point belongs to the projected face if three conditions are met. These conditions depend on the orientation of the sub-face with respect to the displacement vector.

  • When the displacement vector is aligned with the sub-face orientation (i.e., A f , i , j face = true ), the point is located to the left of each of the projected edges provided that we follow the sorted vertices (namely X f X i , X i X j and X j X f ). This means that the following three conditions have to be met: first, L f , i edge = true ; second, L f , j edge = false (due to its reverse orientation); third, L i , j edge = true . This case is displayed in the LHS of Figure 15.

  • When the displacement vector is not aligned with the sub-face orientation (i.e., A f , i , j face = false ), the point is located to the right of each of the projected edges provided that we follow the sorted vertices (namely X f X i , X i X j and X j X f ). This means that the following three conditions have to be met: first, L f , i edge = false ; second, L f , j edge = true (due to its reverse orientation); third, L i , j edge = false . This case is displayed in the RHS of Figure 15.

Figure 15 
                           Sketch of the Boolean test for alignment 
                                 
                                    
                                       
                                          A
                                          
                                             f
                                             ,
                                             i
                                             ,
                                             j
                                          
                                          face
                                       
                                    
                                    
                                    \mathcal{A}^{\mathrm{face}}_{f,i,j}
                                 
                              .
It determines if the displacement vector 
                                 
                                    
                                       
                                          
                                             X
                                             O
                                          
                                          ⁢
                                          
                                             X
                                             D
                                          
                                       
                                    
                                    
                                    \mathbf{X}_{O}\mathbf{X}_{D}
                                 
                               is aligned with the oriented sub-face 
                                 
                                    
                                       
                                          
                                             X
                                             f
                                          
                                          ⁢
                                          
                                             X
                                             i
                                          
                                          ⁢
                                          
                                             X
                                             j
                                          
                                       
                                    
                                    
                                    \mathbf{X}_{f}\mathbf{X}_{i}\mathbf{X}_{j}
                                 
                               (and returns true in that case).
The sketch shows the two possible cases: on the left-hand side, the displacement vector is aligned with the sub-face orientation; on the right-hand side, they are not aligned.
In each case, the point 
                                 
                                    
                                       
                                          X
                                          O
                                       
                                    
                                    
                                    \mathbf{X}_{O}
                                 
                               lies within the projected sub-face if it is located on the proper side of all oriented edges (namely 
                                 
                                    
                                       
                                          
                                             X
                                             f
                                          
                                          ⁢
                                          
                                             X
                                             i
                                          
                                       
                                    
                                    
                                    \mathbf{X}_{f}\mathbf{X}_{i}
                                 
                              , 
                                 
                                    
                                       
                                          
                                             X
                                             f
                                          
                                          ⁢
                                          
                                             X
                                             j
                                          
                                       
                                    
                                    
                                    \mathbf{X}_{f}\mathbf{X}_{j}
                                 
                               and 
                                 
                                    
                                       
                                          
                                             X
                                             i
                                          
                                          ⁢
                                          
                                             X
                                             j
                                          
                                       
                                    
                                    
                                    \mathbf{X}_{i}\mathbf{X}_{j}
                                 
                              , with 
                                 
                                    
                                       
                                          i
                                          <
                                          j
                                       
                                    
                                    
                                    i<j
                                 
                              ).
The logical tests 
                                 
                                    
                                       
                                          
                                             L
                                             
                                                α
                                                ,
                                                β
                                             
                                             edge
                                          
                                          =
                                          true
                                       
                                    
                                    
                                    \mathcal{L}^{\mathrm{edge}}_{\alpha,\beta}=\mathrm{true}
                                 
                               are displayed according to their result (
                                 
                                    
                                       
                                          green
                                          =
                                          true
                                       
                                    
                                    
                                    \mathrm{green}=\mathrm{true}
                                 
                               and 
                                 
                                    
                                       
                                          red
                                          =
                                          false
                                       
                                    
                                    
                                    \mathrm{red}=\mathrm{false}
                                 
                              ).
Figure 15

Sketch of the Boolean test for alignment A f , i , j face . It determines if the displacement vector X O X D is aligned with the oriented sub-face X f X i X j (and returns true in that case). The sketch shows the two possible cases: on the left-hand side, the displacement vector is aligned with the sub-face orientation; on the right-hand side, they are not aligned. In each case, the point X O lies within the projected sub-face if it is located on the proper side of all oriented edges (namely X f X i , X f X j and X i X j , with i < j ). The logical tests L α , β edge = true are displayed according to their result ( green = true and red = false ).

Figure 16 
                           Sketch illustrating how the tracking algorithm determines if a particle displacement from 
                                 
                                    
                                       
                                          X
                                          O
                                       
                                    
                                    
                                    \mathbf{X}_{O}
                                 
                               to 
                                 
                                    
                                       
                                          X
                                          D
                                       
                                    
                                    
                                    \mathbf{X}_{D}
                                 
                               crosses a face.
Here, the particle exits the cell through the face on the right.
As a result, when the edges of this face are projected on the plane normal to the displacement vector (here with superscripts 
                                 
                                    
                                       †
                                    
                                    
                                    \dagger
                                 
                              ), the various logical tests 
                                 
                                    
                                       
                                          L
                                          
                                             α
                                             ,
                                             β
                                          
                                          edge
                                       
                                    
                                    
                                    \mathcal{L}^{\mathrm{edge}}_{\alpha,\beta}
                                 
                               and 
                                 
                                    
                                       
                                          A
                                          
                                             f
                                             ,
                                             i
                                             ,
                                             j
                                          
                                          face
                                       
                                    
                                    
                                    \mathcal{A}^{\mathrm{face}}_{f,i,j}
                                 
                               confirm that one of the triangular sub-face is detected as an exit face (green color) while the three other sub-faces are not (red color).
Figure 16

Sketch illustrating how the tracking algorithm determines if a particle displacement from X O to X D crosses a face. Here, the particle exits the cell through the face on the right. As a result, when the edges of this face are projected on the plane normal to the displacement vector (here with superscripts ), the various logical tests L α , β edge and A f , i , j face confirm that one of the triangular sub-face is detected as an exit face (green color) while the three other sub-faces are not (red color).

Figure 17 
                           Sketch showing a particle displacement from 
                                 
                                    
                                       
                                          X
                                          O
                                       
                                    
                                    
                                    \mathbf{X}_{O}
                                 
                               to 
                                 
                                    
                                       
                                          X
                                          D
                                       
                                    
                                    
                                    \mathbf{X}_{D}
                                 
                               going through a warped face: it can be seen that the particle exits the cell through 
                                 
                                    
                                       
                                          X
                                          
                                             I
                                             1
                                          
                                       
                                    
                                    
                                    \mathbf{X}_{I_{1}}
                                 
                               (which belongs to the sub-face 
                                 
                                    
                                       
                                          X
                                          f
                                       
                                    
                                    
                                    \mathbf{X}_{f}
                                 
                              , 
                                 
                                    
                                       
                                          X
                                          
                                             v
                                             1
                                          
                                       
                                    
                                    
                                    \mathbf{X}_{v_{1}}
                                 
                              , 
                                 
                                    
                                       
                                          X
                                          
                                             v
                                             4
                                          
                                       
                                    
                                    
                                    \mathbf{X}_{v_{4}}
                                 
                              ) and re-enters the cell through 
                                 
                                    
                                       
                                          X
                                          
                                             I
                                             2
                                          
                                       
                                    
                                    
                                    \mathbf{X}_{I_{2}}
                                 
                               (which belongs to the sub-face 
                                 
                                    
                                       
                                          X
                                          f
                                       
                                    
                                    
                                    \mathbf{X}_{f}
                                 
                              , 
                                 
                                    
                                       
                                          X
                                          
                                             v
                                             2
                                          
                                       
                                    
                                    
                                    \mathbf{X}_{v_{2}}
                                 
                              , 
                                 
                                    
                                       
                                          X
                                          
                                             v
                                             3
                                          
                                       
                                    
                                    
                                    \mathbf{X}_{v_{3}}
                                 
                              ).
This is naturally handled by the present algorithm which browse through all sub-faces and then counts the number of times a sub-face is crossed: a pair number means that it stays in the current cell while an odd number means that it exits the current cell.
Figure 17

Sketch showing a particle displacement from X O to X D going through a warped face: it can be seen that the particle exits the cell through X I 1 (which belongs to the sub-face X f , X v 1 , X v 4 ) and re-enters the cell through X I 2 (which belongs to the sub-face X f , X v 2 , X v 3 ). This is naturally handled by the present algorithm which browse through all sub-faces and then counts the number of times a sub-face is crossed: a pair number means that it stays in the current cell while an odd number means that it exits the current cell.

To sum it up, the intersection point is inside a given sub-face if the following condition is respected:

( L i , j edge A f , i , j face ) and ( L f , i edge A f , i , j face ) and ( not ( L f , j edge ) A f , i , j face ) = true ,

where the symbol ≡ corresponds to the operator “equivalent” (i.e., it is true if the two Boolean variables have the same value). This test is made for each sub-face. The algorithm is illustrated in Figure 16 (where the intersection point is on the right-hand face).

In the case of warped faces, it is possible for a line ( X O X D ) to cross the same face several times. This is depicted in Figure 17. In that case, the algorithm monitors the number of times that the line ( X O X D ) crossed the face. If this number is even, it means that the particle remains in the current cell (it has left and reentered the cell). If this number is odd, it means that the particle leaves the cell through this face.

A.1.2 Method to estimate the intersection time and position

As mentioned in Section 4.2.1, the new algorithm not only requires information on the cell containing the particle but also on the intersection time and location. This means that the trajectory algorithm described previously has to be extended to provide this information.

Having determined that the line ( X O X D ) does cross a sub-face, the relative time θ = t cross / Δ t necessary to reach this sub-face can be estimated using the free-flight assumption. It gives

(A.1) θ = X O X f ( X f X i X f X j ) X O X D ( X f X i X f X j ) .

This equation is well-posed since we apply it only when we have previously determined that the line actually crosses the face. In fact, the value of 𝜃 actually provides additional information. When 𝜃 is negative, it means that the intersection point is an entrance point for the oriented axis ( X O X D ) . When 𝜃 is positive, it means that the intersection point is an exit point. The number of sub-faces through which the oriented axis ( X O X D ) enters ( n in ) and leaves ( n out ) is then counted. We then check if the particle is indeed in the correct cell thanks to these numbers. The particle is in the correct cell if the number of sub-faces through which the line ( X O X D ) enters in the cell is equal to the number of sub-faces through which it exits the cell and if both of these numbers are not zero (i.e., n in = n out > 0 ).

If the value of 𝜃 defined in equation (A.1) is in the interval θ [ 0 , 1 [ , the particle does actually cross a face during the time step. The position of the particle at the intersection is then simply given using a simple linear interpolation (this linear interpolation is coherent with the assumption of free-flight motion),

X I = X O + θ × X O X D .

When a particle leaves a cell after crossing one face several times (as in wrapped faces), the exit time is considered equal to the largest value in the range [ 0 , 1 [ (i.e., the last exit time).

A.2 Validation on 3-D unstructured meshes

The trajectory algorithm has been tested using various meshes obtained from the FVCA6 benchmark test cases [7]. These meshes were selected to be representative of a range of different meshes, going from a regular Cartesian mesh to a highly distorted mesh with different refinements. The four meshes used here are displayed in Figure 18.

Figure 18

Type of mesh used in the present case, which span a range of possible configurations (from simple Cartesian mesh without wrapped faces to highly distorted meshes with wrapped faces).

(a) 
                        Hexa mesh.
(a)

Hexa mesh.

(b) 
                        Tetra mesh.
(b)

Tetra mesh.

(c) 
                        PrT mesh.
(c)

PrT mesh.

(d) 
                        PrG mesh.
(d)

PrG mesh.

Figure 19 
                     Maximum dimensionless distance between the particles and their center of gravity using the Hexa mesh (
                           
                              
                                 ■
                              
                              
                              \blacksquare
                           
                        ), the Tetra mesh (
                           
                              
                                 ▲
                              
                              
                              \blacktriangle
                           
                        ), the PrT mesh (
                           
                              
                                 ▼
                              
                              
                              \blacktriangledown
                           
                        ) and the PrG mesh (∙).
(All maximum dimensionless distances are smaller than or equal to 1, meaning that the particles stay within the physical domain.)
Figure 19

Maximum dimensionless distance between the particles and their center of gravity using the Hexa mesh ( ), the Tetra mesh ( ), the PrT mesh ( ) and the PrG mesh (∙). (All maximum dimensionless distances are smaller than or equal to 1, meaning that the particles stay within the physical domain.)

Figure 20 
                     Evolution of 
                           
                              
                                 
                                    ⟨
                                    
                                       X
                                       ⁢
                                       X
                                    
                                    ⟩
                                 
                              
                              
                              \langle XX\rangle
                           
                        , 
                           
                              
                                 
                                    ⟨
                                    
                                       U
                                       ⁢
                                       U
                                    
                                    ⟩
                                 
                              
                              
                              \langle UU\rangle
                           
                        , 
                           
                              
                                 
                                    ⟨
                                    
                                       X
                                       ⁢
                                       U
                                    
                                    ⟩
                                 
                              
                              
                              \langle XU\rangle
                           
                         for the point source dispersion in the ballistic limit case.
Comparison between the analytical solution (black line) and numerical results obtained with the new algorithm for the various meshes.
The numerical results are all in agreement with the analytical solution regardless of the mesh used.
Figure 20 
                     Evolution of 
                           
                              
                                 
                                    ⟨
                                    
                                       X
                                       ⁢
                                       X
                                    
                                    ⟩
                                 
                              
                              
                              \langle XX\rangle
                           
                        , 
                           
                              
                                 
                                    ⟨
                                    
                                       U
                                       ⁢
                                       U
                                    
                                    ⟩
                                 
                              
                              
                              \langle UU\rangle
                           
                        , 
                           
                              
                                 
                                    ⟨
                                    
                                       X
                                       ⁢
                                       U
                                    
                                    ⟩
                                 
                              
                              
                              \langle XU\rangle
                           
                         for the point source dispersion in the ballistic limit case.
Comparison between the analytical solution (black line) and numerical results obtained with the new algorithm for the various meshes.
The numerical results are all in agreement with the analytical solution regardless of the mesh used.
Figure 20 
                     Evolution of 
                           
                              
                                 
                                    ⟨
                                    
                                       X
                                       ⁢
                                       X
                                    
                                    ⟩
                                 
                              
                              
                              \langle XX\rangle
                           
                        , 
                           
                              
                                 
                                    ⟨
                                    
                                       U
                                       ⁢
                                       U
                                    
                                    ⟩
                                 
                              
                              
                              \langle UU\rangle
                           
                        , 
                           
                              
                                 
                                    ⟨
                                    
                                       X
                                       ⁢
                                       U
                                    
                                    ⟩
                                 
                              
                              
                              \langle XU\rangle
                           
                         for the point source dispersion in the ballistic limit case.
Comparison between the analytical solution (black line) and numerical results obtained with the new algorithm for the various meshes.
The numerical results are all in agreement with the analytical solution regardless of the mesh used.
Figure 20

Evolution of X X , U U , X U for the point source dispersion in the ballistic limit case. Comparison between the analytical solution (black line) and numerical results obtained with the new algorithm for the various meshes. The numerical results are all in agreement with the analytical solution regardless of the mesh used.

The case considered for validation actually corresponds to the uniform flow described in Section 5.1.1. It consists in a point source dispersion within homogeneous isotropic turbulence. To ensure that no particle is lost, the distance between the particles and the cell center (point source) is tracked. In order to have a representative quantity independent of the mesh, the quantity followed is the dimensionless distance d * . It is defined as d * = X c X p max X Ω c X c X with X c the center of gravity of the cell, X p the position of the particle and Ω c the domain defined by the cell 𝑐. When particles are properly tracked, this distance is always smaller than 1. However, if one of the face crossed by a particle would be missed, the particle would be permanently lost since it could continue its motion without bound.[2] In such cases, the distance from the cell center could diverge and become much greater than 1.

Results obtained with various meshes are compared using a time scale made dimensionless using the Lagrangian time scale t * = t / T L (which is constant for all meshes considered). The results are displayed in Figure 19: we can see that, with 100 000 particles dispersed initially from the point source, the maximum distance to the cell center converges toward 1 but remains always smaller than unity. This proves that the current algorithm is tractable even for 3-D unstructured meshes.

As in Section 5.1, we can also analyze the results obtained for the different second-order moments. We focus here on verifying that the results are consistent regardless of the mesh used, using a given time step (the role of the time step has been detailed in Section 5.1). The results are displayed in Figure 20: it can be seen that all numerical results match the analytical values regardless of the mesh used in such cases. This confirms the accuracy of the algorithm even when 3-D unstructured meshes are used. At this stage, it is also worth noting that 1-D simulations provide the same results as 3-D ones since each of the three directions can be treated independently of the other ones (this is actually a typical characteristic of homogeneous isotropic turbulence).

References

[1] F. Archambeau, N. Méchitoua and M. Sakiz, Code Saturne: A finite volume code for the computation of turbulent incompressible flows—industrial applications, Int. J. Finite Vol. 1 (2004), no. 1, 1–62. Search in Google Scholar

[2] M. L. Bahlali, C. Henry and B. Carissimo, On the well-mixed condition and consistency issues in hybrid Eulerian/Lagrangian stochastic models of dispersion, Boundary-Layer Meteorol. 174 (2020), no. 2, 275–296. 10.1007/s10546-019-00486-9Search in Google Scholar

[3] C. E. Brennen and C. E. Brennen, Fundamentals of Multiphase Flow, Cambridge University, Cambridge, 2005. 10.1017/CBO9780511807169Search in Google Scholar

[4] S. Chibbaro and J.-P. Minier, A note on the consistency of hybrid Eulerian/Lagrangian approach to multiphase flows, Int. J. Multiph. Flow 37 (2011), no. 3, 293–297. 10.1016/j.ijmultiphaseflow.2010.10.010Search in Google Scholar

[5] R. Clift, J. R. Grace and M. E. Weber, Bubbles, Drops, and Particles, Dover Civil Mech. Eng. Ser., Dover, New York, 2005. Search in Google Scholar

[6] T. D. Dreeben and S. B. Pope, Wall-function treatment in pdf methods for turbulent flows, Phys. Fluids 9 (1997), no. 9, 2692–2703. 10.1063/1.869381Search in Google Scholar

[7] R. Eymard, G. Henry, R. Herbin, F. Hubert, R. Klöfkorn and G. Manzini, 3D benchmark on discretization schemes for anisotropic diffusion problems on general grids, Finite Volumes for Complex Applications VI Problems & Perspectives, Springer, Berlin (2011), 895–930. 10.1007/978-3-642-20671-9_89Search in Google Scholar

[8] C. W. Gardiner, Handbook of Stochastic Methods. Vol. 3, Springer, Berlin, 1985. Search in Google Scholar

[9] D. C. Haworth, Progress in probability density function methods for turbulent reacting flows, Progr. Energy Combust. Sci. 36 (2010), no. 2, 168–259. 10.1016/j.pecs.2009.09.003Search in Google Scholar

[10] D. C. Haworth and S. B. Pope, A generalized Langevin model for turbulent flows, Phys. Fluids 29 (1986), no. 2, 387–405. 10.1063/1.865723Search in Google Scholar

[11] D. C. Haworth and S. B. Pope, A pdf modeling study of self-similar turbulent free shear flows, Phys. Fluids 30 (1987), no. 4, 1026–1044. 10.1063/1.866301Search in Google Scholar

[12] R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles, CRC Press, Boca Raton, 1988. 10.1201/9781439822050Search in Google Scholar

[13] A. Innocenti, R. O. Fox and S. Chibbaro, A Lagrangian probability-density-function model for turbulent particle-laden channel flow in the dense regime, Phys. Fluids 33 (2021), no. 5, Article ID 053308. 10.1063/5.0045690Search in Google Scholar

[14] A. Innocenti, N. Mordant, N. Stelzenmuller and S. Chibbaro, Lagrangian stochastic modelling of acceleration in turbulent wall-bounded flows, J. Fluid Mech. 892 (2020), Paper No. A38. 10.1017/jfm.2020.203Search in Google Scholar

[15] R. Löhner and J. Ambrosiano, A vectorized particle tracer for unstructured grids, J. Comput. Phys. 91 (1990), no. 1, 22–31. 10.1016/0021-9991(90)90002-ISearch in Google Scholar

[16] J.-P. Minier, On Lagrangian stochastic methods for turbulent polydisperse two-phase reactive flows, Progr. Energy Combust. Sci. 50 (2015), 1–62. 10.1016/j.pecs.2015.02.003Search in Google Scholar

[17] J.-P. Minier, Statistical descriptions of polydisperse turbulent two-phase flows, Phys. Rep. 665 (2016), 1–122. 10.1016/j.physrep.2016.10.007Search in Google Scholar

[18] J.-P. Minier, A methodology to devise consistent probability density function models for particle dynamics in turbulent dispersed two-phase flows, Phys. Fluids 33 (2021), no. 2, Article ID 023312. 10.1063/5.0039249Search in Google Scholar

[19] J.-P. Minier, R. Cao and S. B. Pope, Comment on the article “An effective particle tracing scheme on structured/unstructured grids in hybrid finite volume/PDF Monte Carlo methods” by Li and Modest, J. Comput. Phys. 186 (2003), no. 1, 356–358. 10.1016/S0021-9991(03)00006-8Search in Google Scholar

[20] J.-P. Minier, S. Chibbaro and S. B. Pope, Guidelines for the formulation of Lagrangian stochastic models for particle simulations of single-phase and dispersed two-phase turbulent flows, Phys. Fluids 26 (2014), no. 11, Article ID 113303. 10.1063/1.4901315Search in Google Scholar

[21] J.-P. Minier and E. Peirano, The pdf approach to turbulent polydispersed two-phase flows, Phys. Rep. 352 (2001), no. 1–3, 1–214. 10.1016/S0370-1573(01)00011-4Search in Google Scholar

[22] J.-P. Minier, E. Peirano and S. Chibbaro, Weak first- and second-order numerical schemes for stochastic differential equations appearing in Lagrangian two-phase flow modeling, Monte Carlo Methods Appl. 9 (2003), no. 2, 93–133. 10.1515/156939603322663312Search in Google Scholar

[23] J.-P. Minier, E. Peirano and S. Chibbaro, PDF model based on Langevin equation for polydispersed two-phase flows applied to a bluff-body gas-solid flow, Phys. Fluids 16 (2004), no. 7, 2419–2431. 10.1063/1.1718972Search in Google Scholar

[24] J.-P. Minier and J. Pozorski, Derivation of a PDF model for turbulent flows based on principles from statistical physics, Phys. Fluids 9 (1997), no. 6, 1748–1753. 10.1063/1.869291Search in Google Scholar

[25] J.-P. Minier and J. Pozorski, Wall-boundary conditions in probability density function methods and application to a turbulent channel flow, Phys. Fluids 11 (1999), no. 9, 2632–2644. 10.1063/1.870125Search in Google Scholar

[26] T. Möller and B. Trumbore, Fast, minimum storage ray-triangle intersection, J. Graphics Tools 2 (1997), no. 1, 21–28. 10.1080/10867651.1997.10487468Search in Google Scholar

[27] M. Muradoglu and S. B. Pope, Local time-stepping algorithm for solving probability density function turbulence model equations, AIAA J. 40 (2002), no. 9, 1755–1763. 10.2514/2.1880Search in Google Scholar

[28] H. C. Öttinger, Stochastic Processes in Polymeric Fluids: Tools and Examples for Developing Simulation Algorithms, Springer, Berlin, 1996. 10.1007/978-3-642-58290-5Search in Google Scholar

[29] E. Peirano, S. Chibbaro, J. Pozorski and J.-P. Minier, Mean-field/PDF numerical approach for polydispersed turbulent two-phase flows, Progr. Energy Combust. Sci. 32 (2006), no. 3, 315–371. 10.1016/j.pecs.2005.07.002Search in Google Scholar

[30] E. Peirano and J.-P. Minier, Probabilistic formalism and hierarchy of models for polydispersed turbulent two-phase flows, Phys. Rev. E 65 (2002), no. 4, Article ID 046301. 10.1103/PhysRevE.65.046301Search in Google Scholar PubMed

[31] S. B. Pope, Application of the velocity-dissipation probability density function model to inhomogeneous turbulent flows, Phys. Fluids A 3 (1991), no. 8, 1947–1957. 10.1063/1.857925Search in Google Scholar

[32] S. B. Pope, Lagrangian PDF methods for turbulent flows, Annu. Rev. Fluid Mech. 26 (1994), no. 1, 23–63. 10.1146/annurev.fl.26.010194.000323Search in Google Scholar

[33] S. B. Pope, On the relationship between stochastic Lagrangian models of turbulence and second-moment closures, Phys. Fluids 6 (1994), no. 2, 973–985. 10.1063/1.868329Search in Google Scholar

[34] S. B. Pope, Turbulent Flows, Cambridge University, Cambridge, 2000. 10.1017/CBO9780511840531Search in Google Scholar

[35] S. B. Pope and Y. L. Chen, The velocity-dissipation probability density function model for turbulent flows, Phys. Fluids A 2 (1990), no. 8, 1437–1449. 10.1063/1.857592Search in Google Scholar

[36] P. P. Popov, R. McDermott and S. B. Pope, An accurate time advancement algorithm for particle tracking, J. Comput. Phys. 227 (2008), no. 20, 8792–8806. 10.1016/j.jcp.2008.06.021Search in Google Scholar

[37] J. Pozorski and J.-P. Minier, On the Lagrangian turbulent dispersion models based on the Langevin equation, Int. J. Multiph. Flow 24 (1998), no. 6, 913–945. 10.1016/S0301-9322(98)00016-0Search in Google Scholar

[38] K. K. Sabelfeld, Random Fields and Stochastic Lagrangian Models: Analysis and Applications in Turbulence and Porous Media, Walter de Gruyter, Berlin, 2012. 10.1515/9783110296815Search in Google Scholar

[39] D. P. Schmidt and C. J. Rutland, A new droplet collision algorithm, J. Comput. Phys. 164 (2000), no. 1, 62–80. 10.1006/jcph.2000.6568Search in Google Scholar

[40] H. Sigurgeirsson, A. Stuart and W.-L. Wan, Algorithms for particle-field simulations with collisions, J. Comput. Phys. 172 (2001), no. 2, 766–807. 10.1006/jcph.2001.6858Search in Google Scholar

[41] S. Subramaniam and D. C. Haworth, A probability density function method for turbulent mixing and combustion on three-dimensional unstructured deforming meshes, Int. J. Engine Res. 1 (2000), no. 2, 171–190. 10.1243/1468087001545128Search in Google Scholar

[42] J. Xu and S. B. Pope, Assessment of numerical accuracy of PDF/Monte Carlo methods for turbulent reacting flows, J. Comput. Phys. 152 (1999), no. 1, 192–230. 10.1006/jcph.1999.6241Search in Google Scholar

Received: 2022-10-25
Revised: 2023-02-12
Accepted: 2023-02-16
Published Online: 2023-03-15
Published in Print: 2023-06-01

© 2023 Walter de Gruyter GmbH, Berlin/Boston

Downloaded on 24.9.2025 from https://www.degruyterbrill.com/document/doi/10.1515/mcma-2023-2002/html
Scroll to top button