4.2. The Method of Characteristics
89
In order to continue the above computation procedure for the next time
step, it is necessary to update concentrations of all the particles caused by
dispersion. Let
where ACp,k+1 results from dispersion. In general, this concentration variation is assumed to be equal to the concentration variation of the grid square
in which the particle P is located, i.e.,
ACp,k+1 = ACi,j,k+I'
(4.2.13)
During the practical calculation, to avoid the occurrence of negative
values of concentration, the following modification is often adopted. When
ACi,j,k+1 < 0, Eq. (4.2.13) will no longer be used to calculate ACp,k+l' but the
following proportional relationship will be used instead:
AC
- AC
Cp,k
p,k+l -
i,j,k+l ~
i,j,k+l
(4.2.14)
As long as IACi,j,k+11 < Ci~j,k+l' we will have I ACp,k+1I < C p . k , so that Cp.k+l
cannot be a negative value.
After the concentration values of all particles at time tk+l have been
obtained, we can repeat the above procedure to get the concentration distribution at time tk+2' This is the major procedure for solving advection-dispersion problems by the method of characteristics.
As stated above, Eq. (4.2.3), Eq. (4.2.4), and Eq. (4.2.6) are all solved by
the method of explicit finite difTerence. In order to maintain the stability of
the calculation procedure, the time step size must be strictly limited.
4.2.4 Treatment of Boundary Conditions and
SourcejSink Terms
When the method of characteristics is used to solve advection-dispersion
problems, various boundary conditions and source/sink terms require special
treatments.
When the particles placed on the inflow boundary move out of the grid
squares connected with the boundary, new particles should be added with
concentrations equal to the concentrations of the inflow fluid at that time.
On the outflow boundary, no special treatments are required as the particles
move out of the boundary. On a no-flow boundary, new locations of the
computed particles may cross the boundary. In this case we can use the
method of mirror images to move the particles to return to the region, see
Figure 4.9. It shows that the particles move along the no-flow boundary.
In grid squares having solute sources, new particles must be continuously
added according to the strength of the sources, while in the grid squares that
have solute sinks we have to remove a number of particles according to the
89
In order to continue the above computation procedure for the next time
step, it is necessary to update concentrations of all the particles caused by
dispersion. Let
where ACp,k+1 results from dispersion. In general, this concentration variation is assumed to be equal to the concentration variation of the grid square
in which the particle P is located, i.e.,
ACp,k+1 = ACi,j,k+I'
(4.2.13)
During the practical calculation, to avoid the occurrence of negative
values of concentration, the following modification is often adopted. When
ACi,j,k+1 < 0, Eq. (4.2.13) will no longer be used to calculate ACp,k+l' but the
following proportional relationship will be used instead:
AC
- AC
Cp,k
p,k+l -
i,j,k+l ~
i,j,k+l
(4.2.14)
As long as IACi,j,k+11 < Ci~j,k+l' we will have I ACp,k+1I < C p . k , so that Cp.k+l
cannot be a negative value.
After the concentration values of all particles at time tk+l have been
obtained, we can repeat the above procedure to get the concentration distribution at time tk+2' This is the major procedure for solving advection-dispersion problems by the method of characteristics.
As stated above, Eq. (4.2.3), Eq. (4.2.4), and Eq. (4.2.6) are all solved by
the method of explicit finite difTerence. In order to maintain the stability of
the calculation procedure, the time step size must be strictly limited.
4.2.4 Treatment of Boundary Conditions and
SourcejSink Terms
When the method of characteristics is used to solve advection-dispersion
problems, various boundary conditions and source/sink terms require special
treatments.
When the particles placed on the inflow boundary move out of the grid
squares connected with the boundary, new particles should be added with
concentrations equal to the concentrations of the inflow fluid at that time.
On the outflow boundary, no special treatments are required as the particles
move out of the boundary. On a no-flow boundary, new locations of the
computed particles may cross the boundary. In this case we can use the
method of mirror images to move the particles to return to the region, see
Figure 4.9. It shows that the particles move along the no-flow boundary.
In grid squares having solute sources, new particles must be continuously
added according to the strength of the sources, while in the grid squares that
have solute sinks we have to remove a number of particles according to the