close

Вход

Забыли?

вход по аккаунту

?

853

код для вставкиСкачать
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.
Документ
Категория
Без категории
Просмотров
3
Размер файла
998 Кб
Теги
853
1/--страниц
Пожаловаться на содержимое документа