Three-Dimensional Streamline Tracing Method over Tetrahedral Domains

: Getting a clear understanding of the fluid velocity field in underground porous media is critical to various engineering applications, such as oil/gas reservoirs, CO 2 sequestration, groundwater, etc. As an effective visualization tool and efficient transport behaviors solution algorithm, the streamline-based method was improved significantly by numerous studies conducted in the last couple of decades. However, the implementation of streamline simulation is still challenging while working with Finite Element Method (FEM) over 3D tetrahedral domains, where the mass conservation is not guaranteed. Considering the increased computational cost to enforce mass conservation in FEM and additional complexity, a new three-dimensional streamline tracing algorithm is presented that only relies on the velocity vector of a flow field on each vertex of a tetrahedron in a 3D unstructured mesh system. Owning to the shape functions and transformation equations between the master element and actual element, the exit coordinate leaving a tetrahedral element can be determined effectively. As a result, Time of Flight (TOF), the coordinate variable along each streamline, can be calculated accurately and efficiently because that the analytical solution depicting the trajectory in Master Element is deduced. The presented streamline-based method is tested under FEniCS, a programming framework for FEM, which eases the implementation and further development of the presented method.

Due to the application of numerical theories such as stream-tube concept, image-of-well, and the superposition principle [20][21][22], the streamline method has been developing since the 1950s in the petroleum industry, which can trace the fluid pathway. After the introduction of the new concept of "Time Of Flight" (TOF, denoted symbolically by τ) introduced in streamline simulation, the modern streamline simulation concept has been presented and studied by many researchers [7,23,24] since the 1990s. From then on, the streamline simulation is not just used to visualize a flow field tool but to solve transport problems in complex domains. One feature of streamline-based simulation is that a complex 3D fluid-flow simulation could be decoupled into a series of 1D solutions along each streamline, which makes the streamline-based simulation efficient when dealing with large, heterogeneous reservoir and complex enhanced oil recovery (EOR) processes. Some principles for most streamline simulations [25,26] are shown as follows: (1) tracing streamlines and constructing a 1D coordinate system according to the concept of TOF; (2) updating streamlines; (3) using TOF to express the mass equations; and (4) solving the 1D formulation of the transport equation along each streamline in terms of TOF.
The streamline tracing algorithm plays an important role in streamline-based simulation. One category of such algorithm is to use numerically methods such as Runge-Kutta-type solvers to trace streamlines. These methods are applicable if all we need is to generally visualize flow pathways. However, to solve the transport behaviors along these streamlines (known as the concept of modern streamline stimulation), Runge-Kutta-type methods are not acceptable, because these methods are not accurate or efficient enough for 3D large-scale domains and very difficult to construct the TOF grid system [27]. The other category is semi-analytical approaches, which are a good choice for tracing streamlines for the presented purpose. In early streamline-based simulation applications, Lin [22] adopted the image-well method to trace streamlines in a 2D reservoir, which is not robust to determine the positions and numbers through a trial-and-error procedure.
Then, Numbere [28] built a streamline system via the Boundary Element Method, which can address the problem of the image-well method, but it is difficult to be implemented in a 3D and heterogeneous computation domain. Fortunately, Pollock [29] has presented a very famous semi-analytical streamline tracing approach, assuming that the velocity in all directions is linearly distributed and does not affect each other. In addition, the concept of time of flight (TOF) has been introduced, which is efficient and widely applied to a lot of streamline-based simulation in 3D reservoirs [24,27,[30][31][32][33][34][35]. The Pollock method is widely used to establish streamlines over structured grids. Prevost et al. [30] extended the Pollock method to finite-volume flow solutions using distorted quadrilaterals (2-D) or hexahedra (3-D). Although the Pollock method has many advantages in tracing streamlines, the drawbacks are obvious. For instance, it relies on conservative flux at each facet of a cell, and it could not trace streamlines in a mesh system when the cells contain sink or source [26,29,36].
Secondly, all the formulations of tracing the streamlines should depend on a control volume. Thirdly, during the tracing streamline process, another numerical solver is needed to solve the Ordinary Differential Equations (ODEs). Therefore, the TOFs could not be computed efficiently, which limits the implementation of the streamline-based simulation.
Fourthly, this type of streamline tracing method is just a visualization tool for the fluid flow field and could not implement streamline-based simulation.
Fifth, although there are some presented methods for streamline tracing methods in 2D and 3D unstructured mesh system such as triangle mesh [31,39,44,47], Perpendicular Bisector (PEBI) cell [32], and fully unstructured 3D Voronoi grids [45,46], to the best knowledge of the authors, few research studies on tracing streamlines over tetrahedral domains (created by FEM) are based on the velocity vector on each vertex.
The objectives of this paper are, on one hand, to present a new three-dimensional streamline tracing method based on the flow field in unstructured tetrahedron mesh systems where the conservative flux is not necessary; on the other hand, the presented work is implemented and tested in FEniCS; when the flow solver runs in FEniCS, only the finite element mesh and the resultant velocity vector are needed. The advances and features of the new method are summarized as follows.
(1) During the process of streamline tracing, only the coordinate and velocity vector on the vertex are acquired. Neither flux on element facets nor the flux conservation within an element need to be considered. (2) For any given coordinate of an entrance point in an unstructured tetrahedron element, the corresponding exit coordinate can be determined. (3) The time of flight can be calculated accurately and efficiently because the analytical solution that could describe the pathway of a particle in an unstructured tetrahedron element is deduced.

Mathematics Fundamentals for Tracing Streamlines over Tetrahedron Grid System
To trace streamlines over the tetrahedron grid system created by the Finite Element Method (FEM), the mathematics fundamentals are summarized as following sections.

Transformation Formulation of Coordinates and Velocity Vectors between Actual and Mater Element
Arbitrary shape tetrahedron elements in Actual Element (x-y-z plane) are shown in Figure 1a, the geometrical is very complex for tracing a streamline, so it is very important to minimize the geometrical complexity. Therefore, tetrahedron elements in Actual Element need to transform into Master Element (ξ-η-ζ plane) where the geometrical complexity is minimized, as shown in Figure 1b, and the transform process schematic from Actual Element to Master Element is shown as Figure 1.
The first variable needed to transform between Master and Actual Elements is the coordinates. The coordinates in Actual Elements could transform to Master Elements easily by Equation (2), which only needs to take coordinates (x,y,z) into Equation (2) and replace U; the result is shown in Equation (5). The (x1,y1,z1), (x2,y2,z2), (x3,y3,z3) and (x4,y4,z4), are the coordinates of four vertexes of the tetrahedron in the Actual Element.  Then, transformation coordinates (x, y, z) in Actual Elements to coordinates (ξ, η, ζ) in Master Element is more complex. The transformation formulation is shown as Equation (6), which is derived from Equation (5).
Another variable needed to transform between Actual and Master Elements is the velocity vector. For any given velocity vector [vx, vy, vz] at a point (x, y, z) in an Actual Element, it needs to be transformed to the Master Element ([vξ, vη, vζ]). According to Equation (6), the transformation formula for velocity vectors from an Actual Element to a Master Element is expressed as shown in Equation (7).

The Analytical Solution for Describing the Trajectory in Master Element
In order to obtain the analytical solution for describing the trajectory in Master Element, it needs to know the motion equations at any point (ξ, η, ζ), that is to say, the velocity vector for any point is the key point. According to the shape function in the Master Element, the velocity vector of any point yields Equation (8) The initial conditions are shown in Equation (9).
Then, the matrix format of Equation (8) is described as shown in Equation (10): According to Gregory et al. [48], the general solution of Equation (11) is of the form where the particular function Φ1(t), Φ2(t) and Φ3(t), and coefficients E1, E2, E3, and C depend on the eigenvalues of A. There are nine separate cases that X(t) takes depending on different combinations eigenvalues. Here is Case 1 (Equation (12)), which means A has three nonzero real eigenvalues (0 ≠ r1 ≠ r2 ≠ r 3≠ 0), and the other cases are shown in Appendix A.
(1) Case 1: A has three nonzero real

Computing the TOF and Exit Coordinate
For any particle in the Mater Element whose start point is P0, it will leave tetrahedron elements at exit points Pe on one face of the tetrahedron element. The time that the particle spends to leave the tetrahedron element from the start point to the exit point is TOF (Time of Flight), which is denoted by ∆τ. The ∆τ and Pe could be computed by the following steps.
Step 3: Determination of the coordinate of Pe. Now that the real TOF is determined, the exit coordinate (ξe, ηe, ζe) could be determined by the following equations.

The Method to Trace Streamlines over Tetrahedron Grid created by FEM
A new method is presented to trace streamlines over the tetrahedron grid system created by FEM, and the flowchart is shown in Figure 2. Based on a given velocity field, the process of tracing one streamline is shown in the steps below: (1) Import a tetrahedron mesh system created by FEM and the velocity field.
(2) Give a start point P0 in the x-y-z plane (Actual Element); then, find the cell where the start point is located, and obtain these element properties such as the coordinates on each vertex of the cell and the corresponding velocity vector.

Tracing Streamline over a Tetrahedron Element
This subsection will elaborate how to trace a streamline over a tetrahedron element in an Actual Element, which includes Step 2 to Step 7.
In Step 2, for any given start point (P0(x0, y0, z0)) in the Actual Element, find the tetrahedron element where P0 is located and obtain the parameters of this tetrahedron element, such as the coordinates of each vertex (P1(x1, y1, [vξ4, vη4, vζ4] at ξ-η-ζ coordinate of (0,0,0), (1,0,0), (0,1,0), (0,0,1), respectively) in the Master Element via Equation (7). In Step 4, based on the start point coordinate (ξ0, η0, ζ0) and tetrahedron element properties (coordinates and velocity vectors at each vertex), the velocity vector [vξ0, vη0, vζ0] at the start point can be calculated by the shape function. Now, the start points properties (coordinate and velocity vectors) and tetrahedron element properties (coordinate and velocity vectors at each vertex of the tetrahedron element) have been converted into Master Element. In Step 5, now it is time to determine the pathway of a streamline over one tetrahedron element, noting that this is only a streamline segment. The equation for illuminating the trajectory of a streamline is presented in Section 2.3. When given a time, t, the position coordinate in the Master Element can be determined. Additionally, the TOF in one tetrahedron element could be calculated according the process expressed in Section 2.4. Then, the exit point coordinate could be calculated by taking the calculated TOF into the corresponding equation (such as Equation (17)).
The procedure of tracing a streamline segment in a tetrahedron element and the relative formulations in Section 2.3 and Section 2.4 are demonstrated, which are verified by an example shown in Figure 3. A tetrahedron is defined in Actual Element with four vertices, which are P1(1, 1, 1), P2(3, 2, 2), P3 (1,4,3), and P4(1, 5, 1) whose unit is m, respectively, as well as the corresponding velocity vector at four vertices are also given. When a particle starts from a start point which is located at one face of this tetrahedra element, then it will leave this tetrahedron element at an exit point. The procedure of tracing a streamline segment over a tetrahedron element in the Master Element is shown in Figure 3, which is referred to as Case 1, whose start point P0(1. , which yields Case 1, whose feature is that Equation (12) has three nonzero and real eigenvalues, as shown in Section 2.3. The velocities unit is m/s. The geometrical sketch and parameters of a tetrahedron element in Actual Element (such as the coordinates and velocity vectors) are shown in Figure 3a. Then, coordinates and velocity vectors are transformed into Master Element via Equations (6) and (7), respectively, and the result is shown in Figure 3b. Thereby, the streamline pathway, Time of Flight (TOF) and the exit point coordinate Pe * are determined by relative formulations in Section 2.3 and Section 2.4; the result is also shown in Figure 3b. Finally, the results of the streamline segment (streamline pathway and exit point) could be converted to Actual Element, which is shown in Figure 3c. This example verified the procedure of tracing a streamline segment over a tetrahedron element. The calculation results by relative formulations in Section 2.3 and Section 2.4 are verified by the RK4 method. The RK4 is the so-called fourth-order Runge-Kutta method, which is adopted to verify ODEs, since it is easily to implement during calculation. Equation (10) has nine cases according to the eigenvalue combination caused by the different velocity vectors at each vertex. Figure 4 shows nine verification cases for Case 1 (three nonzero real eigenvalues shown as Equation (12), Case 2 (one real and two equal-nonzero real eigenvalues shown as Equation (A1)), Case 3 (a single nonzero real eigenvalue shown as Equation (A2)), Case 4 (one nonzero real and two complex eigenvalues shown as Equation (A3)), Case 5 (one zero real and two nonzero real eigenvalues shown as Equation (A4)), Case 6 (one nonzero real and two zero eigenvalues shown as Equation (A5)), Case 7 (only zero eigenvalues shown as Equation (A6)), Case 8 (one zero and two equal-nonzero real eigenvalues shown as Equation (A7)), Case 9 (two complex and one zero eigenvalues shown as Equation (A8)).
The comparison results are shown in Figure 4 All the given velocities are in Actual Element, and all corresponding values are in Master Element, which is denoted as v1 * , v2 * , v3 * , and v4 * can be found in Figure 4. The streamline for different cases is calculated, and we compared the analytical solution and the 4th-order Runge-Kutta method, and the procedures of using the 4th-order Runge-Kutta method solve partial differential equation (such as Equation (10)), which is shown in Appendix B.
Despite RK4 having numerical errors, if reducing the value of Δt, the results will be convergent to an accurate solution. When Δt = 0.001 s, all the results of RK4 match perfectly with the result of the analytical solution, which is shown in Section 2.3 in this study. In addition, the TOFs (Δτ) calculated by RK4 give equal results in this study (value is shown in Figure 4). Although RK4 can also compute a streamlined pathway, the RK4 is not suited for streamline-based simulation due to its low efficiency and numerical errors. As is shown in Table 1, the error of TOF and the number of iterations by RK4 is up to large ∆t or small ∆t. Taking Case 1 as an example, when using small ∆t(∆t = 0.0001s), the RK4 could calculate an accurate TOF (error of TOF is 0.007%), but it needs to repeat 45645 iterations; while using large ∆t(∆t = 0.01s), the number of iterations decreased drastically (only 450 iterations), but the error of TOF is 1.42%, which is bigger. Considering the efficiency and the accuracy, 0.001 s is selected to calculate TOF by RK4, which was used to verify the analytical solution. Table 1 implies that although RK4 could determine the trajectory and TOF, there are some problems limiting its applications. (1) The computational accuracy and efficiency depend on the step(∆t) size. A small ∆t could determine an accurate TOF and coordinate, which decrease the computational efficiency, while a large ∆t could increase the computational efficiency but could not determine an accurate TOF and coordinate. (2) When tracing streamlines in the entire domain, only TOF and exit coordinates are needed to calculate each element, and it is not significant to know the detailed trajectory in one element. However, it requires repeating a large number of iterations to determine the TOF and exit coordinate by the RK4 method compared to only a few iterations in the analytical solution method; the process of determining the exit points and TOF by the analytical solution method is outlined in Section 2.4.

Tracing Streamline over the Whole Domain
This subsection demonstrates the process of tracing streamlines in the whole domain and building up a streamline system. After tracing a streamline over a tetrahedron element, it is important to make sure whether the exit point reaches the sink or not. If yes, this streamline is the end. However, if it does not, take the exit point as the start point in the next element and repeat the step of tracing a streamline over a tetrahedron element until the exit point reaches the sink. One streamline is formed by connecting the start point and exit point in each tetrahedron element. The process of tracing streamlines for the whole domain is verified by Case 1, which is shown in Section 4.1.

Results and Discussion
Some cases are studied to verify the process of tracing streamlines for the whole domain and illustrate the computational power of the method presented in this paper about tracing streamlines. The purpose of this paper is to propose a general tracing streamline method, so it does not involve any specific physical problems. The following case is based on the Laplace type equation, and the velocity and TOFs are calculated, but not in units. The flow field of the case in this paper is achieved in FEniCS; the details about FEniCS are shown in https://fenicsproject.org. The streamline tracing method in this paper is implemented via the FEniCS framework.

Case 1
The process of tracing streamlines for the whole domain is verified by a simple fluid flow problem (reference as Case 1), which is shown in Figure 5. Figure 5a is showing a simple fluid flow field of u in a unit cubic domain whose governing equation is Laplace's equation, ∇ 2 u = 0, where u, which is scalar, represents the potential field, and υ represents the velocity vector =  As the velocity field is a uniform velocity vector, the TOF of one streamline is linear distributed from 0 at the left face boundary (x = 0) to 1 at the right face (x = 1), which means the analytical solution of TOF of one streamline could be represented as a line (slope equals one) on the graph of the relationship between TOF and x. Therefore, the TOF of this study could be verified by an exact solution, which is shown in Figure 5b as the solid line. As is shown in Figure 6, the TOFs of four streamlines are in perfect agreement with the analytical solution of TOFs on the plot of TOF and x, which means the new streamline tracing method over the tetrahedron mesh system presented in this paper is accurate.

Case 2
Case 2 is used to study the computational power of tracing a streamline in a vertical direction and horizontal direction, respectively. The purpose of this paper is to propose a general tracing streamline method, so it does not involve any specific physical problems.
Case 2 is shown in Figure 7, the calculation domain is a cube with length equaling 1.0, the governing equation is a Laplace equation, whose boundary condition is that the six faces of cube under the no-flow Neumann boundary condition (∂u/∂n = 0), and the source and sink are under the Dirichlet boundary condition (u = 1 and u = 0, respectively). The source is noted as a cylinder whose top coordinate is (0.1, 0.1, 1.0), bottom coordinate is (0.1, 0.1, 1.0), and radius is 0.05, while the sink is noted as a cylinder whose top coordinate is (0.9, 0.9, 1.0), bottom coordinate is (0.9, 0.9, 0.0), and radius is 0.05. The discretization of the calculation domain is shown in Figure 7a. In the FEniCS framework, the resolution of a mesh system is determined via specifying a resolution number instead of the average mesh size, such as the '20′ in code "mesh = generate_mesh(geometry, 20)". The calculation domain was discretized into 8595 elements. Then, the Laplace equation is solved by FEniCS with boundary conditions in the calculation domain, and it extends the program to compute potential field (u), which is shown in Figure 7a. Figure 7b presents the result of tracing 100 streamlines whose start points are located at the surface of the source and linearly distributed along the direction of the cylinder centerline. The TOF color map in Figure 7b,c shows the traveling time of a particle spending from the start point to the sink. Near the start point, the color of the streamline is blue, while near the sink, the color of the streamline is red. In Figure 7b, all the streamlines could start from the start point and reach the surface of the sink, which implies that the new method can trace the streamline in the vertical direction successfully. Figure 7c shows that 100 streamlines were traced in the horizontal direction. The start points of 100 streamlines in Figure 7c are evenly distributed along the horizontal direction at the surface of the source, and all the z coordinates of start points equal 0.5. As shown in Figure 7c, 100 streamlines along the horizontal direction are traced successfully and are evenly distributed on a plane (whose plane equation is z = 0.5), which means that the method for tracing a streamline over the tetrahedron grid system presented in this study can trace all the streamlines in the horizontal direction.
According to the streamline tracing results in Figure 7b,c, the method to tracing a streamline over the tetrahedron grid system presented in this study is a robust method.

Case 3
The method for tracing a streamline over the tetrahedron grid system presented in this study could trace all the streamlines in both the vertical direction and horizontal direction, as implied from Case 2. However, the computational power of the algorithm is not studied when there is an obstacle in the computation domain. So, Case 3 is introduced, where an impermeable sphere is artificially introduced to add computational complexity. The center of the sphere is at (0.5, 0.5, 0.5), and the sphere radius is 0.4. The boundary condition at the surface sphere is the Neumann boundary condition ((∂u/∂n = 0), the boundary conditions are same as those in Case 2. Figure 8a shows the discretization of the calculation domain and the potential field (u), which represents that the fluids flow from the source and bypass the obstacles before finally flowing out from the sink. Figure 8b shows 100 streamlines (10 layers in the vertical direction, 10 streamlines launched in the horizontal direction) are traced successfully, and all the streamlines start from the source, then steer by obstacles and reach the sink.

Conclusions
In this paper, a new method to trace streamlines over unstructured tetrahedron elements based on vertex velocities is presented. The following conclusion are drawn: (1) The exit coordinate leaving a tetrahedral element and the Time of Flight (TOF) could be computed accurately and efficiently because the analytical solution that could depict the trajectory in Master Element is deduced. (2) Two verification cases are studied, whose results are agreeing to the new method to trace streamlines over unstructured tetrahedron elements based on vertex velocities. One of the verifications verifies the analytical solution for depicting the trajectory in Master Element by comparing it with the results of the RK4 method, while the other one verifies the process of tracing streamlines in the entire domain. (3) Two cases with different geometric complexity are studied to demonstrate the computational power of the method presented in this paper. The new method to trace streamlines over unstructured tetrahedron elements is a robust method that can handle a complex computational domain. (4) The presented method is implemented in FEniCS, which is a programming framework for Finite Element Method (FEM) that can be easily extend the method to other research fields. If the velocity field and tetrahedron mesh system are obtained during any fluid field solving in FEniCS, the new method can trace the streamlines of this fluid field, which can easily extend the method to other research fields.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
According to Gregory and II-Hong work (Gregory and II-Hong 1999), Equation (10) has nine solutions depending on the different combinations of eigenvalues. Case 1 is shown as Equation (12), which means A has three nonzero real eigenvalues (0≠ r1≠r2≠r3≠0), and the other cases are shown as follows.

Appendix B
The procedures of using the 4th-Order Runge-Kutta method to solve a partial differential equation (such as Equation (10)) are summarized as follows.    where h is step size (small incremental traveling time, ∆τ, as used in this study), ξn, ηn, and ζn represent the coordinate in the master element at current iterative step, and ξn+1, ηn+1, and ζn+1 are for the next iterative step. Via solving Equations A9-A13 sequentially and repeating the iterative procedure, the trajectory can be depicted in the master element, and the exit coordinate as well as the corresponding traveling time (TOF) can be determined accordingly.