INTERNATIONALJOURNALFOR NUMERICAL METHODS IN ENGINEERING, VOL. 39, 1115-1 136 (1996) A PARTICLE TRACKING TECHNIQUE FOR THE LAGRANGIAN-EULERIAN FINITE ELEMENT METHOD IN MULTI-DIMENSIONS HWAI-PING CHENG, JING-RU CHENG AND GOUR-TSYH YEH Department of Civil and Environmental Engineering, 212, Sackett Building, The Pennsylvania State University, University Park, P A 16802, U.S.A. SUMMARY This paper presents a multi-dimensional particle tracking technique for applying the Lagrangian-Eulerian finite element method to solve transport equations in transient-state simulations. In the LagrangianEulerian approach, the advection term is handled in the Lagrangian step so that the associated numerical errors can be considerably reduced. It is important to have an adequate particle tracking technique for computing advection accurately in the Lagrangian step. The particle tracking technique presented here is designed to trace fictitious particles in the real-world flow field where the flow velocity is either measured or computed at a limited number of locations. The technique, named ‘in-element’ particle tracking, traces fictitious particles on an element-by-element basis. Given a velocity field, a fictitious particle is traced one element by one element until either a boundary is encountered or the available time is completely consumed. For the tracking within an element, the element is divided into a desired number of subelements with the interpolated velocity computed at all nodes of the subelements. A fictitious particle, thus, is traced one subelement by one subelement within the element. The desired number of subelements can be determined based on the complexity of the flow field being considered. The more complicated the flow field is, the more subelementsare needed to achieve accurate particle tracking results. A single-velocity approach can be used to efficiently perform particle tracking in a smooth flow field, while an average-velocityapproach can be employed to increase the tracking accuracy for more complex flow fields. KEY WORDS particle tracking; Lagrangian-Eulerian finite element methods INTRODUCTION It has been determined that by using an appropriate Lagrangian-Eulerian (LE) numerical approach, one can greatly reduce most types of numerical error, including spurious oscillation, numerical spreading, grid orientation and peak clipping/valley elevating, that are associated with the advection term in solving advection-dispersion transport equations. As a result, extremely accurate results can be obtained.’.’ However, this was achieved under the assumption of using an exact particle tracking scheme in the numerical approach. In general, a backward particle tracking scheme is used to trace fictitious particles, along characteristics, from specific locations such that the Lagrangian concentration assodated with those locations can be obtained through interpolation. A forward particle tracking scheme is employed to (1) determine rough elements, in which the concentration distribution is not smooth, for local grid refinement and (2) capture peaks/valleys in the concentration profile, with respect to space, so that the numerical diffusion due to peak capturing/valley elevating can be greatly reduced.’ For most non-uniform or unsteady flows, analytical particle tracking is not available. Thus, the particle tracking technique CCC 0029-5981/96/071115-22 0 1996 by John Wiley & Sons, Ltd. Received 27 January 1995 Revised 4 May 1995 1116 H.-P. CHENG. J.-R. CHENG AND G.-T. YEH plays a key role in the accuracy of a LE approach. Most tracking algorithms are rather simplified probably due to the simplicity of the flow fields used in common test problems. Baptista employed a fourth-order Runge-Kutta scheme to accurately trace particles in a one-dimensional d ~ m a i n Unfortunately, .~ the scheme is computationally expensive if it is applied to multidimensions. On the other hand, it is unlikely to be able to analytically describe the flow field for the real-world system. Instead, the flow field is either measured in the field or computed through numerical simulation. Thus, the velocity is known at a limited number of locations when contaminant transport is computed. These locations are usually specified as the nodes of the associated discretized domain for numerical modeling. For locations without known velocities, interpolation techniques can be used to compute the velocities at those locations if needed. The ‘in-element’ particle tracking technique has been developed to accurately and efficiently trace fictitiousparticles in the velocity fields of real-world systems. The algorithm of the technique is described in the next section. Detailed computations involved in the technique are stated in Appendix I. ALGORITHM Figure 1 shows a flow chart for this technique and a two-dimensional example. The principle of ‘in-element’tracking in a three-dimensionalvolume is the same as that in a two-dimensional area. By following the arrows in the flow chart, the whole scope of the ‘in-element’ tracking technique can be described. For convenience, forward particle tracking is considered in the following description. Action 1. Start tracking a particle: No matter how many particles are considered, each particle can be traced independently. Action 2. Determine the first element, say element M , which the particle will pass through: By considering the geometric relationship, it can be determined if the particle will pass through the element being checked under a given flow field (Appendix 1.1). Unless the particle is located on the boundary and the velocity is outward on this boundary, the first element, which the particle will pass through, can always be found by checking those elements adjacent to the particle in order. In Figure 1, element M is the first element. Action 3. Regularly refine element M into NW subelements: In this refinement, NW working subelementsand NPW working nodes are generated. Both the coordinates and the velocities for working nodes are interpolated from those of the nodes of element M. The connection relationship between working subelementsand working nodes is also computed. NW and NPW can be computed for a two-dimensional element as follows: (1) For a quadrilateral element: N W = NXW x NYW NPW = (NXW + 1) x (NYW + 1) (2) For a triangular element: N W = NXW x NXW N P W = (NXW + 1) x (NXW + 2)/2 1117 A PARTICLE TRACKING TECHNIQUE 2 ?’ NW4 Global elements: I , 456 SUbelemCntS~, 34 At= At1 + At2 + At3 + AtC + AtS Path 1: Element 1, subelement 1 Path 2: Element 1, subelement 3 Path 3: Element 1, subelement 4 Path 4:Element 4, subelement 2 Path 5: Element 5 , subelement 1 -- I ~tarttracking a particIe.1 ’ I 1+ Determine the 1-st global element, say element M,1 which the particle will pass through. refine elemeit M into NW subelements.) I F & l a r l y The LOOP which the particle will pass through. I 1 Compute the particle kcking in subelement MW. 1 The I 11 Deternine the successive subelement, say MWI, which the particle will pass through. b -[LetMW=MWI. L I I ~ Determine the successive element, say M1, which the particle will pass through. Figure 1. The flow chart of the ‘in-element’ particle tracking technique and a two-dimensionalexample While for a three-dimensional element, NW and NPW are calculated as follows: (1) For a hexahedral element: NW = NXW xNYW xNZW NPW = (NXW + 1)x (NYW + 1)x (NZW + 1) (2) For a triangular prism element: NW = NXW x NXW x NZW NPW = (NXW + 1)x (NXW + 2) x (NZW + 1)/2 1118 H.-P. CHENG, J.-R. CHENG AND G.-T. YEH (3) For a tetrahedral element: NW = NXW x NXW x NXW NPW = (NXW + 1) x (NXW + 2) x (NXW + 3)/6 NXW, NYW and NZW, prescribed as desired, are the numbers of refinement in the X-,Y- and Z-directions, respectively. For the two-dimensional example in Figure 1, both NXW and NYW are assigned to be 2. Thus, NW is 4 and NPW is 9. Action 4. Determine thejrst subelement, say subelement M W , which the particle will pass through: Under the same consideration as mentioned in action 2, it can be determined which of its surrounding subelements the particle will pass through. The first subelement, which the particle will pass through, can always be detected (e.g. subelement 1 in element 1 for path 1 in Figure 1). Action 5. Compute the particle tracking in subelement M W To compute the particle tracking in a subelement, the following two steps are taken: Step 1. Determine the subelement side on which the particle will end up. This can be accomplished under the consideration similar to that in action 2 (see Appendix 1.1). Step 2. Locate the ending point on the ending side. The basic concept involved in this determination is that both the coordinates and velocities of the ending point can be interpolated from those of the nodes of the ending side. Thus, there are interpolation parameters to be computed, such that the ending point can be described through interpolation. The interpolation parameters are computed based on the linear velocity-displacement relationship. In this computation, either the single-velocity approach (taking the velocity of the starting point as the tracking velocity) or the average-velocity approach (taking the average velocity of the starting and the ending points as the tracking velocity) can be used. For both situations, the Newton-Raphson method is employed to compute the interpolation parameters (Appendix 1.2). Judgement 1. Is the time completely consumed? This can be examined by comparing the available tracking time, which is the time-step size before the first tracking, with the time consumed for the tracking just computed in subelement MW. If the available time is completely consumed, then a particle has reached the end of tracking (e.g. the end of path 5 in Figure 1). Otherwise, the available time is updated by subtracting the time just consumed from the prior available time, and the next tracking is to proceed (e.g. the ends of paths 1 through 4 in Figure 1). Judgement 2. Will the particle pass through any other part of element M? This can be examined by examining (1)whether or not the newly-determinedending point is on the boundary of element M and (2) whether or not the velocity at this ending point will move the particle away from element M if it is on the boundary of element M. If the particle will continue its trip in element M (e.g. the ending points of paths 1 and 2 in Figure I), then the subelement associated with the next tracking is to be determined. Otherwise, the next tracking will be in an element adjacent to element M (e.g. the ending point of path 3 in Figure 1). Action 6. Determine the successive subelement, say M W l , which the particle will pass though The newly-determined ending point becomes the starting point for the next tracking. The subelement for the successive tracking can be determined by examining the subelements connected with the new starting point based on the same consideration described in action 4 (see Appendix 1.1). Subelement 3 in element 1 is the successive subelement after finishing path 1 in Figure 1. A PARTICLE TRACKING TECHNIQUE 1119 Action 7. Let MW = MW1: Update the working subelement to continue the particle tracking. Action 8. Determine the successive element, say MI,which the particle will pass through The computation in this action is basically the same as that in action 6, except that the fictitious particle here is on the interface between elements and the search is for the successive element, . rather than subelement. In Figure 1, element 4 is the successive element after three particle tracking paths in element 1 in Figure 1. Action 9. Let M = M I : Update the working element to continue the particle tracking. Action 10. End tracking a particle: The end of tracking a particle is reached. 5 6 1 * Element 1: Triangular prism element, 1-2-56-7-1 0 * Element 2: Hexahedral element, 2-3-4-5-7-8-9-10 * Element 3: Tetrahedral element, 6-7-10-11 * Element 4: Hexahedral element, 7-8-9-10-13-14-1 5-1 6 Element 5: Triangular prism element, 7-10-11-13-16-12 * Startingpoint of particle tracking:Node 1 Ending point of particle tracking : Point Q * 10 tracking paths are included to complete the particle tracking: Path 1 is in Element 1, Paths 2 through 5 are in Element 3, Paths 6 through 8 are in Element 5, Paths 9 and 10 are in Element 4. Figure 2. A threedimensional example of particle tracking 1120 H.-P. CHENG, J.-R. CHENG AND G.-T. YEH Figure 3. The breakdownof the particle trackingprocess for a three-dimensional example:(a) Element 1, (b) Element 3, (c) Element 5 and (d) Element 4 As shown in Figure 1, two loops are included in this technique. The outer loop is called the element loop which processes the particle tracking through elements, while the inner loop is the subelement loop which handles the particle tracking through subelements. Figure 2 demonstrates an example of ‘in-element’ tracking in a three-dimensional domain which is composed of two hexahedral, two triangular prism, and one tetrahedral elements. In this example, the numbers of refinement, NXW, NYW and NZW are set to 2, 2 and 1, respectively. Therefore, a hexahedral element is divided into four subelements, a triangular prism element is divided into four subelements,and a tetrahedral element is divided into eight subelements. Figure 3 shows the breakdown of the particle tracking process for this example. From this figure, one can see how the ‘in-element’tracking technique is implemented through elements and subelements. In this particle tracking, the fictitious particle starts from node 1, passes through elements 1,3,5 and stops at point Q which is inside element 4. Ten particle tracking paths are included one in element 1 (Figure 3(a)),four in element 3 (Figure 3(b)),three in element 5 (Figure 3(c))and two in element 4 (Figure 3(d)). It is also observed in Figure 2 that element 2 is not involved in the tracking process, therefore, it is not necessary to refine this element. The interface is a triangular side between elements 1 and 3, another triangular side between elements 3 and 5 and a quadrilateral side between elements 5 and 4. Once elements are connected compatibly,all three types ofelements can be used to discretize the domain of interest and there is no difficulty in executing the ‘in-element’tracking. The only requirement in using this technique is that every side of a subelement must be a plane. This is because the technique is based on the geometric relationships among points, lines, and planes. There is no doubt that a triangular side is a plane. However, all the nodes of a four-node side have to be assured to sit on the same plane. In 1121 A PARTICLE TRACKING TECHNIQUE other words, when hexahedral and/or triangular prism elements are employed, their four-node sides must be like quadrilateral elements in a two-dimensional area. Otherwise, the ending point cannot be identified. DEMONSTRATION In this demonstration, two examples are included. All the quantities in these two examples are normalized, hence, they are dimensionless. Points P and Q denote the starting and the ending points, respectively, in the two examples. V,, V, and V , are the velocity components. Example 1. A two-dimensional particle tracking A two-dimensional backward particle tracking is performed under the following flow field -A Y v, = and 500 AX 500 v y =- A region of [ - 3000, 30001 x [ - 3000, 3000) is discretized with 36 elements and 49 nodes. Each element is rectangular with a size of 1000 x 1000. The fictitious particle to be backward tracked is originally at (0,2000).After a tracking-time period of 500, the fictitious particle will be at (0, - 2000), which can be analytically determined by using the following relationships. jr dX = is"" 0 V,dt i: and dY = is' 0 Vydt (2) where XQ = 0 and YQ = - 2000 after solving equation (2) with equation (1). Eleven cases, with varied NXW, NYW and IJUDGE, are included in this example to examine the performance of the 'in-element' tracking technique. IJUDGE is an indicator for the approach for tracing fictitious particles. IJUDGE equals 1if the average-velocity approach is used, whereas it is set to 2 if the single-velocity approach is considered. One might choose the single-velocity approach for smooth flow fields to save computational time, whereas one might use the average-velocity approach for complex flow fields to increase the tracking accuracy. Table I lists the setup of the eleven cases. Figures 4-7 illustrate the associated particle tracking processes of Table I. List of setups for the particle tracking of Example 1 Case 1A 1B 1c 1D 1E 1F 1G 1H 11 1J 1K NXW NYW IJUDGE 1 1 2 2 3 3 1 1 2 2 3 3 4 4 4 5 5 10 4 5 5 10 1 2 1 2 1 2 1 2 1 2 2 1122 H.-P. CHENG. LR. CHENG AND G.-T. YEH -3000 -2000 -1000 X 0 1000 2000 3000 Figure 4. The particle tracking process for case 1A (1B) 3000 2000 1000 Y O -1000 -2000 -3000 -3000 -2000 -1000 X 0 1000 2000 3000 Figure 5. The particle tracking process for case 1B -3000 -2000 -1000 X 0 1000 2000 3000 Figure 6 . The particle tracking process for case 1C 1123 A PARTICLE TRACKING TECHNIQUE (W 3000 2000 1000 Y O -lo00 -2000 -3000 -3000 -2000 -1000 X 0 1000 2000 3000 Figure 7. The particle tracking process for case 1K cases lA, lB, 1C and lK, respectively. The black dots in the figures represent either a starting point or an ending point of a particle tracking path with respect to a subelement. The dashed lines shown in those figures are used to indicate the numbers of refinement. In practice, only the elements which the fictitious particle passes through need to be refined in implementing the technique. As expected, if larger numbers of refinement (i.e. NXW and NYW) are used, a more accurate result can be obtained for both IJUDGE = 1 and IJUDGE = 2. Additionally, the results associated with the average-velocity approach are more accurate than those associated with the single-velocity approach. When NXW and NYW are greater than or equal to 2 with IJUDGE = 1 (cases lC, lE, lG, and lI), the associated numerical results are very close to the analytical solution (Figure 6). However, even though NXW = NYW = 10 is employed for the single-velocity approach (case lK), the numerical result is not as accurate as that using the average-velocity approach without any refinement (case 1A). This can be checked by comparing Figures 4 and 7. When the single-velocity approach is used without any refinement (case lB), the fictitious particle is backward tracked to the bottom boundary (point Q* in Figure 6) before the available tracking time is completely consumed. With the above observations, it becomes apparent that we should use the average-velocity approach to both accurately and efficiently trace fictitious particles under such a non-uniform flow field. Example 2. A three-dimensional particle tracking The following flow field is applied to a three-dimensional backward particle tracking. - nZ v, = v,= 500 ’ 2, nX v,=- 500 (3) A region of [ - 3000,3000] x [0,3000] x [ - 3000,3000] is discretized with 108 elements and 196 nodes. Each element is cubic with a size of 1000 x loo0 x 1000. The fictitious particle to be backward tracked is originally at (0,2000,2000). After a tracking-time period of 500, the fictitious particle will be at (0, 1000, - 2000) based on equation (3). Seven cases, with varied NXW, NYW, NZW and IJUDGE, are considered in this example (Table 11). The (a) part of Figures 8 through 12 illustrate the particle tracking results of cases 2A, 1124 H.-P. CHENG, J.-R.CHENG AND G.-T. YEH Table 11. List of setups for the particle tracking of Example 2 Case NXW NYW NZW ~ 2A 2B 1 1 3 3 3 5 5 2c 2D 2E 2F 2G IJUDGE ~ 1 1 3 3 1 5 5 ~ 1 1 1 3 3 3 2 1 2 2 5 1 5 2 3000 (a) 2000 1000 z o -1000 -3000 -2000 -1000 X 0 1000 2000 3000 2000 Y 1000 500 400 3 00 200 Available tracking time 100 0 Figure 8. The particle tracking process for case 2A on (a) the X - 2 plane and (b) the Y-Available tracking time plane 2B, 2D, 2E, and 2F on the projected X - 2 plane, while the (b) part shows the associated result on the projected Y-Available tracking time plane. Consistent conclusions with those made in Example 1 are obtained after analyzing the results. Because of the constant V y ,the displacement in the Y-direction is always a linear function of the available tracking time regardless of the use of the approach and the refinement in the 1125 A PARTICLE TRACKING TECHNIQUE 3000 (a) 2000 1000 z o -1000 -Zoo0 -3000 -3000 Y -2000 -1000 X 0 1000 ZOO0 3000 2oo 2 1000 500 400 300 200 100 0 Available tracking time Figure 9. The particle tracking process for case 2B on (a) the X-Z plane and (b) the Y-Available tracking time plane ‘in-element’ particle tracking. Thus, a straight line is observed in all the figures regarding the projected Y-Available tracking time plane. Accurate results of Y, always can be obtained except for case 2B. This is because the fictitious particle encounters the bottom boundary (Figure 9(a)) before the available tracking time is totally consumed. Figure 9(b) shows the incompletelyconsumed available tracking time for case 2B. Cases 2D and 2E have the same setup except that NYW is assigned to 3 in case 2D but 1in case 2E. This does not affect the displacement in the Y-direction due to the constant Vy.However, the displacementsin both the X-and Z-directions are influenced. In Figure l q a ) (for case 2D), we can observe one black dot inside the grid of [1666*7,2000]x [lW, 1333.33and one black dot inside the grid of [2000, 2333.31 x [ - 1333-3, - 10001, whereas no black dot is inside those two elements in Figure 1l(a) (for case 2E). The two extra black dots in Figure 10(a)are generated due to the refinement in the Y-direction. As a result, these two black dots make the particle tracking result in Figure 10(a)more accurate than that in Figure 1l(a) even though we might not be able to distinguish it visually. This is because finer subelements are involved in the particle tracking in case 2D. Based on the above, we can say that if the flow field has a constant component in a direction, then the tracking accuracy in that direction will not be affected by the numbers of refinement. However, if the velocity components in other directions are not constant, the number 1126 H.-P.CHENG, J.-R. CHENG AND G.-T. YEH 3000 (a) 2000 1000 z o -lo00 -2000 -3000 -3000 -2000 -1000 X 0 1000 2000 3000 3000 (b) ......................................................... "M L ................. 2 I ....................................... ................................... ...................................................... 1000 .................................................... t----t ......................................................... 0 1 500 : : 400 : : 300 : I : 200 ; 100 : I I 0 Available tracking time Figure 10. The particle tracking process for case 2D on (a) the X-Z plane and (b) the Y-Available tracking time plane of refinement in the direction of constant velocity will affect the tracking accuracy in those other directions to some degree. For a particle tracking within a subelement, if the average-velocity approach cannot be used due to the local complexity of the flow field, then the single-velocityapproach will be used to trace the particle for the time being. An example of this would be a converging flow in a subelement (Figure 13(a)). The average-velocity approach is not applicable because the ending point cannot be determined. In cases 2A, 2C and 2F, the average-velocity approach is employed to trace the fictitious particle (i.e. IJUDGE = 1). In Figure 8(a), seven particle tracking paths are included in the particle tracking process of case 2A. IJUDGE is assigned to 1 for this case. However, the sixth particle tracking path (i.e. from point A to point B) is achieved with IJUDGE = 2 because there is no appropriate ending point available with IJUDGE = 1 when point A is the starting point. In other words, the complexity of the flow field around point A and its surrounding subelements makes the average-velocity approach inapplicable. After using the single-velocity approach for the sixth particle tracking path in Figure 8(a),we resume IJUDGE to 1 for the seventh particle tracking path. As can be imagined, finer subelements are needed to make the flow field smoother over subelements, so that the average-velocity approach can be used. When the number of 1127 A PARTICLE TRACKING TECHNIQUE -3000 -2000 -1000 500 400 X 0 1000 300 200 Available tracking time 2000 100 3000 0 Figure 11. The particle tracking process for case 2E on (a) the X - 2 plane and (b) the Y-Available tracking time plane refinement is increased to 5 (case 2F), the numerical result is very close to the analytical solution (Figure 12). CONCLUSIONS This paper presents the ‘in-element’ particle tracking technique to help deal with the advection term in the Lagrangian step when the finite element method is used to solve transport equations. The ‘in-element’particle tracking technique is developed to provide accurate and efficient particle tracking results under practical flow fields where flow velocity is measured at a limited number of locations. Ten actions and two judgements are included in the algorithm for the technique. The demonstration results show that the single-velocity approach, which involves less computation than the average-velocity approach, can be used to enhance the computational efficiency if the flow field is smooth over global elements. More subelements are needed to increase the computational accuracy if the flow is more complicated. The average-velocity approach can be employed to handle rough flow fields, such that not only accuracy is greatly improved, but also less computational time is needed because less subelements are required to achieve accurate results. With an adequate combination of IJUDGE and the refinement numbers in the ‘inelement’ particle tracking technique, therefore, we can make the particle tracking both accurate and efficient. 1128 H.-P. CHENG, J.-R. CHENG AND G.-T. YEH 3000 2000 1000 z o - 1000 -2000 -3000 -3000 -2000 -1000 X 0 1000 3000 2000 3000 ......................................................... (b) :mm i wl ......................................................... t ......................................................... =.. 1 ..................................................... ...................................................... ..................................................... ...................................................... ...................................................... 1000 ......................................................... ......................................................... ......................................................... ......................................................... O 500 J : ; 400 : : 300 : : ; 200 ; I00 : I 0 Availabletracking time Figure 12. The particle tracking process for case 2F on (a) the X-Z plane and (b) the Y-Available tracking time plane ACKNOWLEDGEMENT This research is supported, in part, by the Office of Health and Environmental Research, Department of Energy, Grant No. DE-FG02-91ER61197and in part, by Robert Kerr Environmental Laboratory, U.S. EPA under cooperative agreement No. CR-18322-01-0. The technical revision by Karen M. Salvage is much appreciated. APPENDIX I For convenience, the particle tracking path mentioned in this Appendix represents only a tracking process within a subelement. In order to increase the accuracy of the particle tracking result, one can opt to use the average velocity of the starting and ending points of each particle tracking path to trace a particle in a subelement. Alternatively, one can opt to use the velocity of the starting point of a particle tracking path to perform the particle tracking to save computational time. We define IJUDGE = 1 to represent the first case (the average-velocity approach) and IJUDGE = 2 for the second case (the single-velocity approach). In practice, if a flow field has a drastic change around a certain location, one may not be able to use the average velocity to perform the particle tracking near that location (e.g. Figure 13(a)). This occurs when the subelement used for particle tracking cannot resolve the flow field appropriately. Under such A PARTICLE TRACKING TECHNIQUE 1129 Figure 13. An example demonstrating a converging flow over a quadrilateral subelement where a forward particle tracking with the average-velocity approach, starting from point P, cannot be achieved if one subelement is used for particle tracking (part (a)), whereas it can be achieved if four subelements are used (part (b)) a situation, IJUDGE will be temporarily set to 2 for computing the particle tracking around that location and be reset to 1 for the rest of the particle tracking. If the subelement for particle tracking is set finely enough to resolve the flow field, the average-velocity approach may be used. For example, by using four subelements in Figure 13(b), rather than one in Figure 13(a), for particle tracking in the same converging flow area, the velocity field can be resolved within each subelement such that the average-velocity approach can be used. Points P and Q denote the starting and ending points, respectively, of a forward particle tracking path within a subelement. VX,VY and V Z are the components of the tracking velocity for the tracking path. (XI, YI,Z , ) and (VXI, V Y , , V Z , ) represent the coordinate and velocity, respectively, of point I. Although we will focus on forward particle tracking in this Appendix, backward particle tracking can be achieved by exchanging the roles of points P and Q. Appendix 1.1. Determination of the direction that a fictitious particle will pass along With the velocity4isplacement relationship, the direction that a fictitious particle will pass along can be determined. If a line (for two-dimensional cases) or a plane (for three-dimensional cases) contains point P, then the geometrical relationship between point Q and the line or the plane can be determined with a vector analysis. Based on this, we can determine (1)if a fictitious 1130 H.-P. CHENG, J.-R.CHENG AND G.-T. YEH particle will pass through a subelement and (2) on which side of the 'pass-through' subelement point Q will be located. The following computations can be applied not only to subelements but also elements. I. 1. I. i%e determination in a two-dimensionaldomain. As shown in Figure 14, point P is on line L which is determined by points A and P. Point A represents a subelement node here. Given the locations and the velocities at both points A and P, the quantity, Det, shown in Figure 14 is computed as follows. Let AX = XA - X p (4) L\Y = YA - Yp (5) if IJUDGE = 2 V X , + VX,) if IJUDGE = 1 if IJUDGE = 2 (7) then Det = A X . V Y - AY - VX (8) By computing Det, the location of point Q can be determined as indicated in Figure 14. In a two-dimensional space, we may use quadrilateral and/or triangular elements to discretize the domain of interest when the finite element method is employed. Point P can be (1) a subelement node, (2) a point on a side of a subelement,or (3) a point inside a subelement during the process of the 'in-element' particle tracking. If the starting point is on a subelement side, the side will not be eligible to be the ending side. For instance, two subelement sides are not eligible to be the ending side of the quadrilateral subelement if point P is a node of the element, while there is one side not eligible for the case of a side point and all sides are eligible for the case of an interior point. As described above, we can examine the location of the ending point relative to a given line which contains the starting point. Thus, we also can determine whether or not an eligible subelement side is the ending side by computing the values of Det associated with the two lines determined by the starting point and the two end points of the subelement side. The two lines are treated the same as line L in Figure 14. On Line L: Det = 0 __..,....- ... Above Line L Det > 0 ,,,.... w: ...-p / ...,_,..... ,_./.... DetLine C 0 L: ,._." _...- ...'_.- ....._.. ._.. Line L Figure 14. The determination of the tracking direction of a fictitious particle in a two-dimensional space A PARTICLE TRACKING TECHNIQUE 1131 1.1.2. The determination in a three-dimensional domain. As shown in Figure 15, the starting point P is on plane I which is determined by points A, B and P. Points A and B denote the two end points of a subelement edge. Given the locations and the velocities at points A, B and P, the quantities, Detl and Det2, shown in Figure 15 is computed as follows. Let AX, = XA- x p (9) AY1 = YA - Yp (10) A z 1 = ZA- z p (11) AX2 = XB - Xp (12) AY2 = YB - Yp (13) A 2 2 = ZB - Zp (14) On Plane I: Detl = 0 and Det2 = 0 On the left side of Plane I: Detl<OmdDet2<Q n Plane I B Detl O; andDet2> 0 Figure 15. The determination of the tracking direction of a fictitious particle in a threedimensional space 1132 H.-P. CHENG, J.-R. CHENG AND G.-T. YEH then Detl = V1.(D1xD2) (21) Det2 = V2 .(D1x D2) (22) Dz = (AXz, A Y2, AZ2) (26) where By computing Detl and Det2, the location of point Q can be determined as described in Figure 15. One thing that should be noted is that we use plane P-A-B to represent plane I in Figure 15. The values of Detl and Det2 associated with plane P-A-B are different from those associated with plane P-EA by a negative sign. In applying the finite element method to solving a three-dimensional transport problem, we may use hexahedral, triangular prism, and/or tetrahedral elements to discretize the domain of interest. Point P can be (1) a subelement node, (2) a point on an edge of a subelement,(3) a point on a side of a subelement, or (4)a point inside a subelement in the ‘in-element’particle tracking process. Only when point P is not on a specific subelement side, is the side eligible to be the ending side. For a hexahedral subelement,for example, there are three sides eligible to be the ending side for the case of a subelement node, four sides eligible for the case of an edge point, five sides eligible for the case of a side point, and six sides eligible for the case of an interior point. By computing Detl and Det2, we can examine the location of the ending point relative to a given plane which contains the starting point. In the same way, we can determine whether or not an eligible subelement side is the ending side by computing the values of Detl and Det2 associated with the planes determined by the starting point and the two end points of each side line of the subelement side. One might notice that Detl and Det2 may not be both positive, both zero, or both negative when IJUDGE = 1 (i.e. using the average velocity) is considered. Such a situation indicates a drastic change in velocity within the subelement being considered. If there is no ending side detected with IJUDGE = 1, IJUDGE will be assigned to 2 which guarantees determination of the ending side. To summarize, the computations in Appendix 1.1 are used in determining (1) if a fictitious particle will pass through a subelement and (2) on which side of a ‘pass-through’subelement the ending point will be encountered. Because the same concept is employed to achieve both purposes, they are not differentiated above. Appendix 1.2. Determination of the ending point of a particle tracking As soon as the ending side is determined, the location of the ending point can be determined as follows. 1.2.1. The computation in a two-dimensional domain. For a two-dimension particle tracking path, the ending side is a line segment. Let points 1 and 2 be the end points of this line segment. 1133 A PARTICLE TRACKING TECHNIQUE Then, the following equations can be written by using linear interpolation: where ’X and V Y can be calculated as follows. 4 + VXp = (2VXp) &VXp VXQ)= [(VXp + VX2) + [(VXI - VX,)] if IJUDGE = 2 if IJUDGE = 1 (34) VYp = 4(2VYP) $(VYp VYQ) = *[(VYp if IJUDGE if IJUDGE (35) + + VY2) + l ( V Y 1 - VY,)] =2 = 1 To determine 5, equation (33) can be substituted into either (31) or (32) to eliminate At. If equation (33) is substituted into (31), then the following equation can be obtained. After rearrangement, equation (36) can be rewritten as + F = (XQ- xp)’[vx2 VY2]1’2 - [XQ - Xp)’ + (YQ - Yp)2]1’2’ Vx = 0 (37) In equation (37),Xp and Yp are known values, whereas XQ, YQ, VX, and V Y are functions of 5 according to equations (27H30), (34) and (35). In other words, F is a function of i.Therefore, 5 can be determined by using the Newton-Raphson method to solve equation (37). After computing 5 (5 E [0, 13) to satisfy both equations (31) and (32) both the location and velocity of point Q can be calculated using equations (27H30). The time consumed for the particle to move along this particle tracking path is then calculated using equation (33). If the value is smaller than the remaining available time, then the remaining available time is updated by using its prior value minus the consumed time and the particle tracking continues. If the value is larger than the remaining available time, the remaining available time is the real tracking time for this path and point Q is located as follows. XQ = xp + (XQ- 6t xp)’t (38) 1134 H.-P. CHENG, J.-R. CHENG AND G.-T. YEH where St is the real tracking time in this particle tracking path, and At is the tracking time computed by equation (33). In this case, the updated remaining a vailable time becomes zero and the particle tracking ends here. 1.2.2. The computation in a three-dimensional domain when the ending side is a quadrilateral side. For a three-dimension particle tracking path, the ending side can be either a quadrilateral side with four points, say points 1, 2, 3 and 4, as its plane points or a triangular side with three points, say points 1,2 and 3, as its end points. If it is a quadrilateral side, the following equations can be written by using the bilinear interpolation. (40) XQ=X,’Nl +XZ*N2+X3.N3+X4‘N4 + Y z * N Z + Y3”3 + Y4eN4 (41) ZQ=Zl.Ni +Zz*N2+Zj.N3+24*N4 (42) Y Q = YI*N1 + VXz.Nz + VX3.N3 + VX4*N4 V Y Q = VY1SN1 + V Y z . N z + V Y 3 . N 3 + V Y 4 . N 4 V Z , = V Z 1 . N I + V Z z - N z + V Z 3 . N 3 + VZ4mN4 V X Q = VX1.N1 (43) (44) (45) where N1,N 2 , N 3 and N4 are the shape functions in the local coordinate systems. N1 = d(1 - O(1 - rl) (46) + 5)(1 - V ) (47) + rl) N4 = d(1 - O(1 + rt) (48) Nz = N3 = ac1 + 5)(1 (49) < in which and q represent the two independent variables in the local coordinate system. Both variables must be located within [ - 1,1]. The following equations can be described based on the velocity-displacement relationship. x At=[..- X Q - Xp = A t . V X (50) Y , - Yp = A t . V Y (5 1) ZQ - z (52) p = At. V z + ( Y Q- Yp)’ + (ZQ - Zp)’ vx2+ V Y Z+ Y Z 2 )’ 1 1’2 (53) where V X , V Y and V Z are calculated as follows. vx i VXp =+(2VXP) = i ( V x p VXQ) =3[VXp N l . V X 1 + + + Nz * VX’ + N 3 . vx3 + N 4 . V X 4 ] VYp =f(2VYp) =j[VYp + N 1 . V Y I + N z ’ V Y 2 + N 3 - V Y 3 + Nq.VY4-J if IJUDGE = 2 if IJUDGE = 1 (54) if IJUDGE = 2 if IJUDGE = 1 (55) 1135 A PARTICLE TRACKING TECHNIQUE I/z = 1 vzp = +(2VZP) f(Vzp VZQ) = +[VZP + Nl.VZ1 + + N ” V Z ’ + N3-VZ3 + Nq.VZ41 if IJUDGE = 2 if IJUDGE = 1 (56) To determine 5 and q, equation (53)can be substituted into any two of the equations (50H52) to eliminate At. If equation (53) is substituted into equations (50) and (51), then the following equations can be obtained. After mathematical manipulation, equations (57) and (58) can be rewritten as + VY’ -!-Vz2]1’2 - [(XQ - Xp)’ + (YQ- Yp)’ + (ZQ - ZP)’] ”’.I/x = 0 G = (YQ - Yp).[Vx’ + I/Y2 + Vzz]1’2 F = (XQ- Xp)‘[Vx’ -[(xQ-Xp)’+(YQ- Yp)’ + ( z Q - z p ) z ] ’ / 2 . ~ Y = o (59) (60) In equations (59) and (60),Xp, Yp and Zpare known values, whereas XQ, YQ,ZQ,VX, V Y and VZ are functions of 4 and q based on equations (40H49) and (54H56). Therefore, F and G are functions 5 and q. 5 and q can be determined by using the Newton-Raphson method to solve equations (59) and (60).After computing and q to satisfy equations (50)-(52), the shape functions N1, N 2 , N 3 and N4can be calculated with equations (46H49). Both the location and the velocity of point Q, thus, can be obtained with equations (40H45). The time consumed along this particle tracking path is calculated with equation (53). The remaining available time is updated as described in Appendix 1.2.1. Point Q is located with the following equations if the particle tracking ends here. 1.2.3. T h e computation in a three-dimensional domain when the ending side is a triangular side. If the ending side is a triangular side, the following equations can be written with the twodimensional linear interpolation. X Q = X , * N l+X2’N2+X3’(1-N1-N2) (64) YQ = Y1*N1+ Y2.N’ (65) + Y,.(l - N1 - Nz) Z Q = Z I ‘ N 1+ Z z ’ N 2 + Z 3 ’ ( 1 - N 1 - N z ) (66) 1136 H.-P. CHENG, J.-R. CHENG AND G.-T. YEH V X o = V X 1 . N 1 + VXZ-NZ + V X , . ( l - N 1 - N z ) (67) VYl-N+ 1 V Y z - N 2+ VY3.(l - N1 - N,) (68) VY, = VZ, = VZ1.N1 + V Z z . N z + VZ3*(1- N 1 - Nz) (69) where N1 and N z represent the two independent variables in the two-dimensional natural coordinate system. Both variables must be located within [0, 11. By following the procedure displayed in Appendix 1.2.2, the location and the velocity of point Q as well as the remaining available time can be computed. REFERENCES 1. G. T. Yeh, J. R. Chang and T. E. Short, An exact peak capturing and oscillation-free scheme to solve advectiondispersion transport equations, Water Resources. Res., zs(1 l), 2937-2951 (1992). 2. G. T. Yeh, J. R. Chang, H. P. Cheng and C. H. Sung, An adaptive local grid refinement based on the exact peak capturing and oscillation free scheme to solve transport equations, Int. J . Comput. Fluids, 24(3), 293-332 (1995). 3. A. M. Baptista, ‘Solution of advection-dominated transport by Eulerian-Lagrangian methods using the backwards method of characteristics’, Ph.D. Dissertation,The Massachusetts Institute of Technology, 1987.