LINEBURG

 << Пред. стр. страница 5(всего 23)ОГЛАВЛЕНИЕ След. стр. >>
where П† is the temperature and b is the thermal conductivity. The heat
stepping with the CFL condition in equation (4.7) as well.
equation is the most basic equation of the parabolic model.
The stringent O(( x)2 ) time-step restriction resulting from the forward
When П† is a signed distance function, bОє| П†| and b П† are identical, and
Euler time discretization can be alleviated by using an ODE solver with
either of these can be used to calculate the right-hand side of equation (4.5).
a larger stability region, e.g., an implicit method. For example, п¬Ѓrst-order
However, once this right-hand side is combined with a forward Euler time
4.3. Convection-Diп¬Ђusion Equations 45 46 4. Motion Involving Mean Curvature

accurate backward Euler time stepping applied to equation (4.6) gives and the two can be used interchangeably if one maintains a signed distance
approximation for П† oп¬Ђ the interface. These equations can be solved using
П†n+1 в€’ П†n the upwind methods from the last chapter on the V В· П† term and cen-
= b П†n+1 , (4.8)
t tral diп¬Ђerencing on the parabolic b П† or bОє| П†| term. A TVD RK time
discretization can be used with a time-step restriction of
which has no time step stability restriction on the size of t. This means
that t can be chosen for accuracy reasons alone, and one typically sets
|u| |v| |w| 2b 2b 2b
t = O( x). Note that setting t = O( x) as opposed to t = O(( x)2 ) t + + + + + <1 (4.12)
2 2 2
x y z ( x) ( y) ( z)
lowers the overall accuracy to O( x). This can be improved upon using
the trapezoidal rule satisп¬Ѓed at every grid point.
Suppose the O(1) size b term is replaced with an O( x) size term
П†n+1 в€’ П†n П†n + П†n+1
=b , (4.9) that vanishes as the mesh is reп¬Ѓned with x в†’ 0. Then equation (4.10)
t 2
becomes
which is O(( t)2 ) in time and thus O(( x)2 ) overall even when t =
П†t + V В· П†= П†, (4.13)
O( x). This combination of the trapezoidal rule with central diп¬Ђerencing
which asymptotically approaches equation (3.2) as в†’ 0. The addition of
of a parabolic spatial operator is generally referred to as the Crank-Nicolson
an artiп¬Ѓcial П† term to the right-hand side of equation (3.2) is called the
scheme.
artiп¬Ѓcial viscosity method. Artiп¬Ѓcial viscosity is used by many authors to
The price we pay for the larger time step achieved using either equa-
stabilize a central diп¬Ђerencing approximation to the convective П† term in
tion (4.8) or equation (4.9) is that a linear system of equations must
be solved at each time step to obtain П†n+1 . Luckily, this is not diп¬ѓcult equation (3.2). This arises in computational п¬‚uid dynamics, where terms
given the simple linear structure of П†n+1 . Unfortunately, an implicit dis- of the form П† are added to the right-hand side of convective equations
to pick out vanishing viscosity solutions valid in the limit as в†’ 0. This
cretization of equation (4.5) requires consideration of the more complicated
nonlinear Оєn+1 | П†n+1 | term. vanishing viscosity picks out the physically correct weak solution when no
classical solution exists, for example in the case of a discontinuous shock
We caution the reader that one cannot substitute equation (4.6) for equa-
tion (4.5) when using an implicit time discretization. Even if П†n is initially wave. It is interesting to note that the upwind discretizations discussed
a signed distance function, П†n+1 will generally not be a signed distance in the last chapter have numerical truncation errors that serve the same
function after the linear system has been solved. This means that П†n+1 purpose as the П† term. First-order accurate upwinding has an intrinsic
is not a good approximation to Оєn+1 | П†n+1 | even though П†n may be O( x) artiп¬Ѓcial viscosity, and the higher-order accurate upwind methods
have intrinsic artiп¬Ѓcial viscosities with magnitude O(( x)r ), where r is the
exactly equal to Оєn | П†n |. Although we stress (throughout the book) the
order of accuracy of the method.
conceptual simpliп¬Ѓcations and computational savings that can be obtained
In , Sethian suggested an entropy condition that required curves to
when П† is a signed distance function, e.g., replacing N with П†, Оє with
п¬‚ow into corners, and he provided numerical evidence to show that this
П†, etc., we caution the reader that there is a signiп¬Ѓcant and important
entropy condition produced the correct weak solution for self-interesting
diп¬Ђerence between the two in the case where П† is not a signed distance
curves. SethianвЂ™s entropy condition indicates that Оє| П†| is a better form
function.
for the vanishing viscosity than П† for dealing with the evolution of lower-
dimensional interfaces. This concept was rigorized by Osher and Sethian in
, where they pointed out that
4.3 Convection-Diп¬Ђusion Equations
П†t + V В· П† = Оє| П†| (4.14)
The convection-diп¬Ђusion equation is a more natural choice than equation (4.13) for dealing with level set
methods, although these two equations are interchangeable when П† is a
П†t + V В· П†=b П† (4.10)
signed distance function.
includes both the eп¬Ђects of an external velocity п¬Ѓeld and a diп¬Ђusive term.
The level set version of this is
П†t + V В· П† = bОє| П†|, (4.11)
48 5. Hamilton-Jacobi Equations

5.2 Connection with Conservation Laws
Consider the one-dimensional scalar conservation law

5 ut + F (u)x = 0; (5.3)

where u is the conserved quantity and F (u) is the п¬‚ux function. A well-
Hamilton-Jacobi Equations known conservation law is the continuity equation

ПЃt + (ПЃu)x = 0 (5.4)

for conservation of mass, where ПЃ is the density of the material. In compu-
tational п¬‚uid dynamics (CFD), the continuity equation is combined with
equations for conservation of momentum and conservation of energy to ob-
tain the compressible Navier-Stokes equations. When viscous eп¬Ђects are
ignored, the Navier-Stokes equations reduce to the compressible inviscid
Euler equations.
The presence of discontinuities in the Euler equations forces one to con-
sider weak solutions where the derivatives of solution variables, e.g., ПЃx ,
can fail to exist. Examples include linear contact discontinuities and non-
linear shock waves. The nonlinear nature of shock waves allows them to
develop as the solution progresses forward in time even if the data are ini-
5.1 Introduction tially smooth. The Euler equations may not always have unique solutions,
and an entropy condition is used to pick out the physically correct solu-
In this chapter we discuss numerical methods for the solution of general tion. This is the vanishing viscosity solution discussed in the last chapter.
Hamilton-Jacobi equations of the form For example, vanishing viscosity admits physically consistent rarefaction
waves, ruling out physically inadmissible expansion shocks.
BurgersвЂ™ equation
П†t + H( П†) = 0, (5.1)
u2
ut + =0 (5.5)
where H can be a function of both space and time. In three spatial 2 x
dimensions, we can write
is a scalar conservation law that possesses many of the interesting non-
linear properties contained in the more complex Euler equations. BurgersвЂ™
П†t + H(П†x , П†y , П†z ) = 0 (5.2)
equation develops discontinuous shock waves from smooth initial data and
exhibits nonphysical expansion shocks if the vanishing viscosity solution is
as an expanded version of equation (5.1). Convection in an externally gen- not used to force these to become smooth rarefaction waves. Many of the
erated velocity п¬Ѓeld (equation (3.2)) is an example of a Hamilton-Jacobi numerical methods developed to solve BurgersвЂ™ equation can be extended to
equation where H( П†) = V В· П†. The level set equation (equation (4.4)) treat both the one-dimensional and the multidimensional Euler equations
is another example of a Hamilton-Jacobi equation with H( П†) = Vn | П†|. of gas dynamics.
Here Vn can depend on x, t, or even П†/| П†|. Consider the one-dimensional Hamilton-Jacobi equation
The equation for motion by mean curvature (equation (4.5)) is not a
П†t + H(П†x ) = 0, (5.6)
Hamilton-Jacobi-type equation, since the front speed depends on the sec-
ond derivatives of П†. Hamilton-Jacobi equations depend on (at most) the
which becomes
п¬Ѓrst derivatives of П†, and these equations are hyperbolic. The equation for
motion by mean curvature depends on the second derivatives of П† and is (П†x )t + H(П†x )x = 0 (5.7)
parabolic.
5.3. Numerical Discretization 49 50 5. Hamilton-Jacobi Equations

after one takes a spatial derivative of the entire equation. Setting u = П†x or HJ WENO schemes. For brevity, we discuss the two-dimensional numer-
in equation 5.7 results in ical approximation to H(П†x , П†y ), noting that the extension to three spatial
dimensions is straightforward. An important class of schemes is that of
ut + H(u)x = 0, (5.8)
monotone schemes. A scheme is monotone when П†n+1 as deп¬Ѓned in equa-
which is a scalar conservation law; see equation (5.3). Thus, in one spatial tion (5.9) is a nondecreasing function of all the П†n . Crandall and Lions
dimension we can draw a direct correspondence between Hamilton-Jacobi proved that these schemes converge to the correct solution, although they
equations and conservation laws. The solution u to a conservation law is are only п¬Ѓrst-order accurate. The numerical Hamiltonians associated with
the derivative of a solution П† to a Hamilton-Jacobi equation. Conversely, monotone schemes are important, and examples will be given below.
the solution П† to a Hamilton-Jacobi equation is the integral of a solution u The forward Euler time discretization (equation (5.9)) can be extended to
to a conservation law. This allows us to point out a number of useful facts. higher-order TVD Runge Kutta in a straightforward manner, as discussed
For example, since the integral of a discontinuity is a kink, or discontinuity in Chapter (3). The CFL condition for equation 5.9 is
in the п¬Ѓrst derivative, solutions to Hamilton-Jacobi equations can develop
|H1 | |H2 | |H3 |
kinks in the solution even if the data are initially smooth. In addition, t max + + < 1, (5.10)
x y z
solutions to Hamilton-Jacobi equations cannot generally develop a discon-
tinuity unless the corresponding conservation law develops a delta function. where H1 , H2 , and H3 are the partial derivatives of H with respect to П†x ,
Thus, solutions П† to equation (5.2) are typically continuous. Furthermore, П†y , and П†z , respectively. For example, in equation (3.2), where H( П†) =
since conservation laws can have nonunique solutions, entropy conditions V В· П†, the partial derivatives of H are H1 = u, H2 = v, and H3 = w. In
are needed to pick out вЂњphysicallyвЂќ relevant solutions to equation (5.2) as this case equation (5.10) reduces to equation (3.10). As another example,
well. consider the level set equation (4.4) with H( П†) = Vn | П†|. Here the partial
Viscosity solutions for Hamilton-Jacobi equations were п¬Ѓrst proposed by derivatives are slightly more complicated, with H1 = VN П†x /| П†|, H2 =
Crandall and Lions . Monotone п¬Ѓrst-order accurate numerical meth- VN П†y /| П†|, and H3 = VN П†z /| П†|, assuming that VN does not depend on
ods were п¬Ѓrst presented by Crandall and Lions  as well. Later, Osher П†x , П†y or П†z . Otherwise, the partial derivatives can be substantially more
and Sethian  used the connection between conservation laws and complicated.
Hamilton-Jacobi equations to construct higher-order accurate вЂњartifact-
freeвЂќ numerical methods. Even though the analogy between conservation
5.3.1 Lax-Friedrichs Schemes
laws and Hamilton-Jacobi equations fails in multiple spatial dimensions,
many Hamilton-Jacobi equations can be discretized in a dimension by di- Л†
The п¬Ѓrst approximation to H that we consider is the Lax-Friedrichs (LF)
mension fashion. This culminated in , where Osher and Shu proposed scheme from  given by
a general framework for the numerical solution of Hamilton-Jacobi equa-
П†в€’ + П†+ П†в€’ + П†+ П†+ в€’ П†в€’
П†+ в€’ П†в€’
tions using successful methods from the theory of conservation laws. We y y y y
Л† x x x x
в€’О±x в€’О±y
H=H , , (5.11)
2 2 2 2
where О±x and О±y are dissipation coeп¬ѓcients that control the amount of
numerical viscosity. These dissipation coeп¬ѓcients
5.3 Numerical Discretization О±x = max |H1 (П†x , П†y )|, О±y = max |H2 (П†x , П†y )| (5.12)
are chosen based on the partial derivatives of H.
A forward Euler time discretization of a Hamilton-Jacobi equation can be
The choice of the dissipation coeп¬ѓcients in equation (5.12) can be rather
written as
subtle. In the traditional implementation of the LF scheme, the maximum is
П†n+1 в€’ П†n Л†
+ H n (П†в€’ , П†+ ; П†в€’ , П†+ ; П†в€’ , П†+ ) = 0, (5.9) chosen over the entire computational domain. First, the maximum and min-
x x y y z z
t
imum values of П†x are identiп¬Ѓed by considering all the values of П†в€’ and П†+
x x
Л†xxyyzz
where H(П†в€’ , П†+ ; П†в€’ , П†+ ; П†в€’ , П†+ ) is a numerical approximation of H(П†x , П†y , x min max
on the Cartesian mesh. Then one can identify the interval I = [П†x , П†x ].
Л† A similar procedure is used to deп¬Ѓne I y = [П†y , П†max ]. The coeп¬ѓcients
min
П†z ). The function H is called a numerical Hamiltonian, and it is required y
Л† О±x and О±y are set to the maximum possible values of |H1 (П†x , П†y )| and
to be consistent in the sense that H(П†x , П†x ; П†y , П†y ; П†z , П†z ) = H(П†x , П†y , П†z ).
Recall that spatial derivatives such as П†в€’ are discretized with either п¬Ѓrst- |H2 (П†x , П†y )|, respectively, with П†x в€€ I x and П†y в€€ I y . Although it is oc-
x
casionally diп¬ѓcult to evaluate the maximum values of |H1 | and |H2 |, it is
order accurate one-sided diп¬Ђerencing or the higher-order accurate HJ ENO
5.3. Numerical Discretization 51 52 5. Hamilton-Jacobi Equations

П†В± and П†В± at xi,j using the HJ WENO scheme. This type of scheme has
straightforward to do so in many instances. For example, in equation (3.2), x y
both H1 = u and H2 = v are independent of П†x and П†y , so О±x and О±y can been referred to as a Stencil Lax-Friedrichs (SLF) scheme, since it deter-
be set to the maximum values of |u| and |v| on the Cartesian mesh. mines the dissipation coeп¬ѓcient using only the neighboring grid points that
Consider evaluating О±x and О±y for equation (4.4) where H1 = VN П†x /| П†| are part of the stencil used to determine П†x and П†y . An alternative to the
and H2 = VN П†y /| П†|, recalling that these are the partial derivatives only if dimension-by-dimension neighborhoods is to use the 49 grid points in the
VN is independent П†x and П†y . It is somewhat more complicated to evaluate rectangle with diagonal corners at xiв€’3,jв€’3 and xi+3,j+3 to determine О±.
О±x and О±y in this case. When П† is a signed distance function with | П†| = 1 This idea of searching only locally to determine the dissipation coeп¬ѓ-
(or в‰€ 1 numerically), we can simplify to H1 = VN П†x and H2 = VN П†y . cients can be taken a step further. The Local Lax-Friedrichs (LLF) scheme
These functions can still be somewhat tricky to work with if VN is spatially proposed for conservation laws by Shu and Osher  does not look at
varying. But in the special case that VN is spatially constant, the maximum any neighboring grid points when calculating the dissipation coeп¬ѓcients in
values of |H1 | and |H2 | can be determined by considering only the endpoints a given direction. In , Osher and Shu interpreted this to mean that
О±x is determined at each grid point using only the values of П†в€’ and П†+
of Ix and Iy , respectively. This is true because H1 and H2 are monotone x x
at that speciп¬Ѓc grid point to determine the interval I x . The interval I y is
functions of П†x and П†y , respectively. In fact, when VN is spatially constant,
H1 = VN П†x /| П†| and H2 = VN П†y /| П†| are straightforward to work with as still determined in the LF or SLF manner (in the SLF case we rename LLF
as SLLF). Similarly, О±y uses an interval I y , deп¬Ѓned using only the values
well. The function H1 achieves a maximum when |П†x | is as large as possible
of П†в€’ and П†+ at the grid point in question while I x is still determined in
and |П†y | is as small as possible. Thus, only the endpoints of I x and I y need y y
be considered; note that we use П†y = 0 when the endpoints of I y diп¬Ђer the LF or SLF fashion. Osher and Shu  also proposed the Local Lo-
in sign. Similar reasoning can be used to п¬Ѓnd the maximum value of |H2 |. cal Lax-Friedrichs (LLLF) scheme with even less numerical dissipation. At
each grid point I x is determined using the values of П†в€’ and П†+ at that grid
One way to treat a spatially varying VN is to make some estimates. For x x
point; I y is determined using the values of П†в€’ and П†+ at that grid point;
example, since |П†x |/| П†| в‰¤ 1 for all П†x and П†y , we can bound |H1 | в‰¤ |VN |. y y
A similar bound of |H2 | в‰¤ |VN | holds for |H2 |. Thus, both О±x and О±y can and then these intervals are used to determine both О± and О±y . When H is
x

separable, i.e., H(П†x , П†y ) = H x (П†x ) + H y (П†y ), LLLF reduces to LLF, since
be set to the maximum value of |VN | on the Cartesian mesh. The price
О±x is independent of П†y , and О±y is independent of П†x . When H is not sep-
we pay for using bounds to choose О± larger than it should be is increased
numerical dissipation. That is, while the numerical method will be stable arable, LLF and LLLF are truly distinct schemes. In practice, LLF seems
and give an accurate solution as the mesh is reп¬Ѓned, some details of this to work better than any of the other options. LF and SLF are usually too
solution may be smeared out and lost on a coarser mesh. dissipative, while LLLF is usually not dissipative enough to overcome the
Since increasing О± increases the amount of artiп¬Ѓcial dissipation, decreas- problems introduced by using the centrally averaged approximation to П†x
ing the quality of the solution, it is beneп¬Ѓcial to chose О± as small as possible and П†y in evaluating H in equation (5.11). Note that LLF is a monotone
without inducing oscillations or other nonphysical phenomena into the so- scheme.
Л†
lution. In approximating Hi,j at a grid point xi,j on a Cartesian mesh,
it then makes little sense to do a global search to deп¬Ѓne the intervals I x
and I y . In particular, consider the simple convection equation (3.2) where 5.3.2 The Roe-Fix Scheme
О±x = max |u| and О±y = max |v|. Suppose that some region had relatively
As discussed above, choosing the appropriate amount of artiп¬Ѓcial dissipa-
small values of |u| and |v|, while another region had relatively large values.
tion to add to the centrally evaluated H in equation (5.11) can be tricky.
Since the LF method chooses О±x as the largest value of |u| and О±y as the
Therefore, it is often desirable to use upwind-based methods with built-in
largest value of |v|, the same values of О± will be used in the region where
artiп¬Ѓcial dissipation. For conservation laws, Shu and Osher  proposed
the velocities are small as is used in the region where the velocities are
using RoeвЂ™s upwind method along with an LLF entropy correction at sonic
large. In the region where the velocities are large, the large values of О± are
points where entropy-violating expansion shocks might form. The added
required to obtain a good solution. But in the region where the velocities
dissipation from the LLF entropy correction forces the expansion shocks
are small, these large values of О± produce too much numerical dissipation,
to develop into continuous rarefaction waves. The method was dubbed Roe
wiping out small features of the solution. Thus, it is advantageous to use
Fix (RF) and it can be written for Hamilton-Jacobi equations (see )
only the grid points suп¬ѓciently close to xi,j in determining О±. A rule of
as
thumb is to include the grid points from xiв€’3,j to xi+3,j in the x-direction
П†+ в€’ П†в€’
П†+ в€’ П†в€’
and from xi,jв€’3 to xi,j+3 in the y-direction in the local search neighborhood y y
Л† x x
H = H П†x , П†y в€’ О± x y
в€’О± , (5.13)
for determining О±. This includes all the grid nodes that are used to evaluate 2 2
5.3. Numerical Discretization 53 54 5. Hamilton-Jacobi Equations

where О±x and О±y are usually set identically to zero in order to remove formulas and use these cheaper approximations to decide whether or not
the numerical dissipation terms. In the RF scheme, I x and I y are initially upwinding or LLF will be used. After making this decision, the higher-
determined using only the nodal values for П†В± and П†В± as in the LLLF order accurate HJ WENO (or HJ ENO) method can be used to compute
x y
the necessary values of П†В± and П†В± used in the numerical discretization,
scheme. In order to estimate the potential for upwinding, we look at the x y
partial derivatives H1 and H2 . If H1 (П†x , П†y ) has the same sign (either obtaining the usual high-order accuracy. Sonic points rarely occur in prac-
always positive or always negative) for all П†x в€€ I x and all П†y в€€ I y , we tice, and this strategy reduces the use of the costly HJ WENO method by
know which way information is п¬‚owing and can apply upwinding. Similarly, approximately a factor of two.
if H2 (П†x , П†y ) has the same sign for all П†x в€€ I x and П†y в€€ I y , we can
upwind this term as well. If both H1 and H2 do not change sign, we upwind
5.3.3 GodunovвЂ™s Scheme
completely, setting both О±x and О±y to zero. If H1 > 0, information is
п¬‚owing from left to right, and we set П†x = П†в€’ . Otherwise, H1 < 0, and we In , Godunov proposed a numerical method that gives the exact solu-
x
set П†x = П†+ . Similarly, H2 > 0 indicates П†y = П†в€’ , and H2 < 0 indicates tion to the Riemann problem for one-dimensional conservation laws with
x y
П†y = П†+ . piecewise constant initial data. The multidimensional Hamilton-Jacobi
y
If either H1 or H2 changes sign, we are in the vicinity of a sonic point formulation of this scheme can be written as
where the eigenvalue (in this case H1 or H2 ) is identically zero. This signiп¬Ѓes
Л†
H = extx exty H (П†x , П†y ) , (5.14)
a potential diп¬ѓculty with nonunique solutions, and artiп¬Ѓcial dissipation
is needed to pick out the physically correct vanishing viscosity solution. as was pointed out by Bardi and Osher . This is the canonical monotone
We switch from the RF scheme to the LLF scheme to obtain the needed scheme. Deп¬Ѓning our intervals I x and I y in the LLLF manner using only
artiп¬Ѓcial viscosity. If there is a sonic point in only one direction, i.e., x or y, the values of П†В± and П†В± at the grid node under consideration, we deп¬Ѓne
x y
it makes little sense to add damping in both directions. Therefore, we look extx and exty as follows. If П†в€’ < П†+ , then extx H takes on the minimum
x x
for sonic points in each direction and add damping only to the directions value of H for all П†x в€€ I x . If П†в€’ > П†+ , then extx H takes on the maximum
x x
that have sonic points. This is done using the I x and I y deп¬Ѓned as in the value of H for all П†x в€€ I x . Otherwise, if П†в€’ = П†+ , then extx H simply plugs
x x
LLF method. That is, we switch from the LLLF deп¬Ѓned intervals used П†в€’ (= П†+ ) into H for П†x . Similarly, if П†в€’ < П†+ , then exty H takes on the
x x y y
above to initially look for sonic points to the larger LLF intervals that are minimum value of H for all П†y в€€ I y . If П†в€’ > П†+ , then exty H takes on the
y y
even more likely to have sonic points. We proceed as follows. If H1 (П†x , П†y ) maximum value of H for all П†y в€€ I y . Otherwise, if П†в€’ = П†+ , then exty H
y y
y
x
does not change sign for all П†x в€€ ILLF and all П†y в€€ ILLF , we set П†x equal to simply plugs П†в€’ (= П†+ ) into H for П†y . In general, extx exty H = exty extx H,
y y
either П†в€’ or П†+ depending on the sign of H1 . In addition, we set О±x to zero so diп¬Ђerent versions of GodunovвЂ™s method are obtained depending on the
x x
to remove the artiп¬Ѓcial dissipation in the x-direction. At the same time, this order of operations. However, in many cases, including when H is separable,
means that a sonic point must have occurred in H2 , so we use an LLF-type extx exty H = exty extx H so this is not an issue.
method for the y-direction, setting П†y = (П†в€’ + П†+ )/2 and choosing О±y as Although GodunovвЂ™s method can sometimes be diп¬ѓcult to implement,
y y
dictated by the LLF scheme. A similar algorithm is executed if H2 (П†x , П†y ) there are times when it is straightforward. Consider equation (3.2) for mo-
y
x
does not change sign for all П†x в€€ ILLF and П†y в€€ ILLF . Then П†y is set to tion in an externally generated velocity п¬Ѓeld. Here, we can consider the
either П†в€’ or П†+ , depending on the sign of H2 ; О±y is set to zero; and an LLF x and y directions independently, since H is separable with extx exty H =
y y
method is used in the x-direction, setting П†x = (П†в€’ + П†+ )/2 while choosing extx (uП†x ) + exty (vП†y ). If П†в€’ < П†+ , we want the minimum value of uП†x .
x x x x
О±x as dictated by the LLF scheme. If both H1 and H2 change sign, we have Thus, if u > 0, we use П†в€’ , and if u < 0, we use П†+ . If u = 0, we obtain
x x
sonic points in both directions and proceed with the standard LLF scheme uП†x = 0 regardless of the choice of П†x . On the other hand, if П†в€’ > П†+ ,
x x
at that grid point. we want the maximum value of uП†x . Thus, if u > 0, we use П†в€’ , and ifx
With the RF scheme, upwinding in the x-direction dictates that either u < 0, we use П†+ . Again, u = 0 gives uП†x = 0. Finally, if П†в€’ = П†+ , then
x x x
в€’
П†x or П†+ be used, but not both. Similarly, upwinding in the y-direction uses uП†x is uniquely determined. This can be summarized as follows. If u > 0,
x
either П†в€’ or П†+ , but not both. Since evaluating П†В± and П†В± using high-order use П†в€’ ; if u < 0, use П†+ ; and if u = 0, set uП†x = 0. This is identical to
y y x y x x
accurate HJ ENO or HJ WENO schemes is rather costly, it seems wasteful the standard upwind diп¬Ђerencing method described in Chapter 3. That is,
to do twice as much work in these instances. Unfortunately, one cannot for motion in an externally generated velocity п¬Ѓeld, GodunovвЂ™s method is
determine whether upwinding can be used (as opposed to LLF) without identical to simple upwind diп¬Ђerencing.
computing П†В± and П†В± . In order to minimize CPU time, one can compute
x y
П†В± and П†В± using the п¬Ѓrst-order accurate forward and backward diп¬Ђerence
x y
56 6. Motion in the Normal Direction

1 1 1

0.5 0.5 0.5

0 0 0

6 в€’0.5 в€’0.5 в€’0.5

в€’1 в€’1 в€’1
в€’1 0 1 в€’1 0 1 в€’1 0 1
Motion in the Normal Direction 1 1 1

0.5 0.5 0.5

0 0 0

в€’0.5 в€’0.5 в€’0.5

в€’1 в€’1 в€’1
в€’1 0 1 в€’1 0 1 в€’1 0 1

1 1 1

0.5 0.5 0.5

0 0 0

в€’0.5 в€’0.5 в€’0.5

в€’1 в€’1 в€’1
в€’1 0 1 в€’1 0 1 в€’1 0 1
6.1 The Basic Equation
Figure 6.1. Evolution of a star-shaped interface as it moves normal to itself in
In this chapter we discuss the motion of an interface under an internally the outward direction.
generated velocity п¬Ѓeld for constant motion in the normal direction. This
velocity п¬Ѓeld is deп¬Ѓned by V = aN or Vn = a, where a is a constant. The
corresponding level set equation (i.e., equation (4.4)) is Thus, if П†n is initially a signed distance function (with | П†n | = 1), it stays
a distance function (with | П†| = 1) for all time.
П†t + a| П†| = 0, (6.1)
When the initial data constitute a signed distance function, forward Euler
time stepping reduces to solving the ordinary diп¬Ђerential equation П†t = в€’a
where a can be of either sign. When a > 0 the interface moves in the
independently at every grid point. Since a is a constant, this forward Euler
normal direction, and when a < 0 the interface moves opposite the normal
time stepping gives the exact solution up to round-oп¬Ђ error (i.e., there is no
direction. When a = 0 this equation reduces to the trivial П†t = 0, where
truncation error). For example, consider a point where П† = П†o > 0, which is
 << Пред. стр. страница 5(всего 23)ОГЛАВЛЕНИЕ След. стр. >>