Unformatted text preview: the equation that is satisﬁed at convergence is given by
¨¥ ©§¤ ¥ ¦¤ The consequence of this is that the ﬁnal answer we obtain depends on the timestep size ∆t ! Although the solution still has a spatial truncation error of O ∆x 2 , this dependence on the time step is clearly not physical. In this particular case, the resulting error in the ﬁnal solution maybe small since stability requirements restrict the time step size. We will see in later chapters that similar path dependence of the converged solution can occur even in iterative schemes if we are not careful. 5.6 HigherOrder Schemes
We have seen thus far that both the ﬁrst order upwind and central difference schemes have severe limitations, the former due to artiﬁcial diffusion, and the latter due to dispersion. Therefore there has been a great deal of research to improve the accuracy of the upwind scheme, by using higherorder interpolation. These higherorder schemes aim to obtain at least a secondorder truncation error, while controlling the severity of spatial oscillations. Thus far, we have assumed that, for the purposes of writing the face value φ e , the proﬁle of φ is essentially constant. That is, for Fe 0, Instead of a constant proﬁle assumption for φ , we may use higherorder proﬁle assumptions, such as linear or quadratic, to derive a set of upwind weighted higherorder schemes. If Fe 0, we write a Taylor series expansion for φ in the neighborhood of the upwind point P:
¤ ¡ ¤ ¢ £ ¡ ¤ ¡ ¨ ¤ ¨ φe φP φx
£ 5.6.1 SecondOrder Upwind Flux Expressions
We may derive a secondorder accurate expressions for the face ﬂux by making a linear proﬁle assumption. This is equivalent to retaining the ﬁrst two terms of the expansion. Evaluating Equation 5.95 at x e xP ∆x 2, we obtain
¡ ¨ φe φP ∆x ∂ φ 2 ∂x
¤ £ This assumption has a truncation error of O ∆x 2 . In order to write φ e in terms of cell centroid values, we must write ∂ φ in terms of cell centroid values. On our one∂x dimensional grid we can represent the derivative at P using either a forward, backward or central difference formula to give us three different secondorder schemes. If we write ∂ φ using ∂x φE φW ∂φ (5.97) ∂x 2∆x 112
¢ ¨ £ ¤ £ ¡ ¢ ¨ £ φP x xP ∂φ ∂x x xP 2! 2 ∂ 2φ ∂ x2 O ∆x 3 ¤ ¤ £ ¢ £ ¢ ¢ ¢ £ ¢ 2 ¡ ¡ φE φP u∆t φ 2∆x E φP φP φW
2 u∆t φ 2∆x P φW 0 (5.93) (5.94) (5.95) (5.96) we obtain
¡ ¨ 4 or, adding and subtracting φ P 4, we get
¤ ¢ £ ¢ This scheme is referred to as the Fromm scheme in the literature. If we write
¢ we obtain a secondorder fully upwind expression for face ﬂux
¢ Schemes using this form of face ﬂux are referred to in the literature as BeamWarming schemes. 5.6.2 ThirdOrder Flux Expressions
Thirdorder accurate expressions for the face ﬂux can be derived by retaining both the ﬁrst and second derivative in the Taylor series expansion, Equation 5.95:
¤ ¢ £ and using cellcentroid values to write the derivatives
and we may write
¡ ¨ Rearranging, we may write
¢ ¡ ¤ ¡ £ This scheme is called the QUICK scheme (Quadratic Upwind Interpolation for Convective Kinetics) [3]. This scheme may be viewed as a parabolic correction to linear 113 ¢ ¨ φe φE φP
2 ¡ 4 φE φW 8 2φP ¢ ¡ ¢ φe φP φE φW φE φW 8 2φP ¤ £ ¡ ¤ ∆x
£ 2 ¢ ¡ ¨ ∂ 2φ ∂ x2 φE φW 2φP O ∆x2 ¤ £ ¡ ¢ ¨ ∂φ ∂x φE φW 2∆x ¡ ¤ ¢ £ ¡ ¨ ¤ φx
£ φP x xP ¡ ¨ φe φP ¨ ∂φ ∂x φP φW ∆x φP 2 ∂φ ∂x ¡ ¡ ¨ φe φP φP φW
4 φE φW x O ∆x2 ¢ φe
φP φE φW (5.98) φP
4 (5.99)
∂φ ∂x using (5.100) (5.101) xP 2!
∂φ ∂x 2 ∂ 2φ ∂ x2
∂ 2φ . ∂ x2 (5.102) Using (5.103) and (5.104) (5.105) (5.106) interpolation for φ e . We can emphasize this by introducing a curvature factor C such that
¤ ¤ ¡ £ ¤ ¢ ¢ ¡ £ ¤ ¨ £ ¡ φe and C=1/8. The second and thirdorder expressions we have seen here may be combined into a single expression for φ e using
¤ ¡ 5.6.3 Second Order Upwind Schemes
So far we have only considered higher order expressions for the face ﬂuxes appearing in Equation 5.79. Before we can use these schemes in practical applications we also need to consider the solution advancement strategy. If we simply use any of the higher order ﬂux expressions given by Equation 5.108 with the ﬁrst order backward Euler differencing for the unsteady term that we employed earlier we will ﬁnd that the resulting scheme is unconditionally unstable. Historically, higher order explicit schemes were formulated by introducing additional terms that compensated for the instability caused by the higher order spatial discretizations. We have already encountered one such scheme, the LaxWendroff scheme described in Section 5.5.5. It is possible to formulate a general procedure whereby any of the second order ﬂux expressions discussed above can be used to create stable schemes that are second order accurate in both space and time. Many schemes have been derived over the years by researchers using such ideas. Like the LaxWendroff scheme, most of these schemes also result in timestep dependence in the steady state solution. As we noted earlier, this occurs because the spatial proﬁle assumption is tied to the temporal discretization. In modern CFD practice such schemes are generally not used and independent temporal and spatial discretizations are usually preferred. For higher order spatial discretizations this requires the use of either some sort of multilevel or multistage temporal discretization (such as the RungeKutta method used for ODEs) or an implicit temporal discretization. However, the resulting methods are more difﬁcult to implement or study analytically. We will study such techniques in a later chapter but for the present discussion we will continue to look at explicit schemes since they are easier to analyze. We will not discuss the general procedure for formulating explicit schemes using higher order ﬂux expressions but present one scheme that results from the fully upwinded second order expression for face ﬂux, Equation 5.101. The scheme is known 114 ¨ Here, κ 1 yields the BeamWarming scheme, κ 0 the Fromm scheme and κ 1 2 the QUICK scheme. For κ 1 we get the familiar central difference scheme.
¨ ¨ ¢ ¨ ¢ £ ¡¤ ¢ ¨ φe φP 1 4 κ φP φW 1 κ
4 ¢ £ 1 φ 2E
£ φP C φE φW 2φP (5.107) φE φP (5.108) 2 1 BeamWarming Exact Solution 1.5 BeamWarming Exact Solution 0.5
1 φ 0 φ
0.5 0.5
0 1 0 0.2 0.4 0.6 0.8 1
0.5 0 0.2 0.4 0.6 0.8 1 x/L x/L Figure 5.10: BeamWarming solutions for linear convection of (a) sine wave (b) square wave as the BeamWarming secondorder upwind scheme. 2
¢ £ ¢ £ ¡ ¢ £ ¡ 1 0 2φW ∆x
¤ £ Figure 5.10 shows the solution obtained using this scheme for the two convection problems we have used for studying the other schemes. We see that for smooth proﬁles the second order upwind ﬂux expression is clearly much less diffusive than the ﬁrst order upwind case. However, at slope discontinuities, we do get oscillations just as with central difference scheme. Although we have only shown it in an explicit scheme, this behaviour occurs for all types of second order schemes, whether in an implicit setting or indeed even when used in steady state discretizations. The astute reader will notice that the BeamWarming solutions shown in Figure 5.10 are almost complimentary to the solutions obtained using the LaxWendroff scheme (Figure 5.9). For symmetric proﬁles, the overshoots and undershoots produced by the two schemes mirror each other at each time step (assuming of course that both are run at the same CFL number). This immediately suggests that one might be able to create a better scheme by combining these two schemes. This was indeed attempted, and the resulting scheme is known as the Fromm scheme. But there is something even more remarkable about these two schemes. It is possible to show that any explicit, linear scheme on the ﬁve point stencil of Figure 5.3 can be written as a linear combination of these two schemes. We will make use of this fact in a later section to design ﬂux limiters applicable to any second order scheme. To summarize, we have seen so far that the ﬁrst order upwind scheme produces physically plausible results but is very diffusive. The second order centraldifference
1 There is an implicit variant of the scheme which is also known as BeamWarming scheme but for the purposes of this chapter we will only be considering the explicit version and refer to it as the BeamWarming scheme. 115 ¨ ∆t
¢ 2 ¤ ¡ ¤ ¤ φP 0 φP u 0 φP 0 φW ∆x u ∆x u∆t 0 φP 0 φWW 0 (5.109) scheme (as well as other symmetric schemes of higher order that we did not discuss) are unstable in the vonNeumann sense when used with ﬁrst order explicit time stepping but even when they can be made stable, they produce nonphysical oscillations in the solution. Attempts at formulating upwindbiased (i.e, nonsymmetric) higher order schemes by improving the interpolations used in the ﬁrst order upwind scheme meet with the same fate. Since the strongest of the oscillations seem to occur wherever there discontinuities in either the solution or its derivatives, we must somehow make the higher order scheme behave like the ﬁrst order upwind scheme while retaining second order accuracy elsewhere in the domain. 5.6.4 Added Dissipation Schemes
There are two broad approaches that can be followed to achieve physically plausible solutions with second order accuracy. The ﬁrst one involves keeping the discretization as it is but damping out the oscillations produced by introducing artiﬁcial viscosity explicitly in the domain. This is typically employed in conjunction with the space symmetric schemes, which have no natural dissipation. For the second order central 2 difference scheme 5.81, we cannot add a dissipative term of the form ∂ xφ since that is ∂2 O ∆x and would thus destroy the second order accuracy of our scheme. Therefore, a discrete term corresponding to the next higher derivative with a dissipative character, i.e., the fourth derivative, is added to the balance equation at each cell. For example, instead of Equation 5.80, the face value is written as
¤ ¡ ¡ This amounts to adding a term of type ∆x 3 ∂ xφ to the control volume balance for the ∂4 cell P. Since the additional term term is O ∆x 3 , it does not affect the formal accuracy of the scheme but it does help damp out any oscillations that might be generated during solution advancement (whether via time marching or iterations). At discontinuities, however, this is not sufﬁcient and stronger dissipation, in the form of a second derivative term, must be added. This leads to face value of the form
4 Of course to use such ideas successfully, we need to somehow specify the coefﬁcients ε 2 and ε 4 , for the second and fourth order dissipation terms respectively. We also require some strategy for detecting discontinuities in the domain so that in smooth regions ε 2 can be made small and near shocks etc. ε 4 can be made small. While the added dissipation procedure requires somewhat adhoc choices and numerical experimentation and does not guarantee that undershoots and overshoots will not generated, it has been quite successfully used by many researchers for some class of problems. We will not discuss this approach in detail but instead look at another option, which is used in conjunction with upwindbiased schemes. Such schemes do have some dissipative character but the higher order derivatives used in constructing the face ﬂuxes can lead to creation of new extrema. The main idea is still to limit 116
¤ ¢ ¡ ¢ £ ¡¤ ¢ £ ¢ 2 ¡ ¨ φe φP φE ε 2 φE φP εe 4 φEE 3φE ¢ ¤ ¢ £ £ 2 ¡ ¨ ¤ £ φe φP φE εe 4 φEE 3φE 3φP φW (5.110) 3φP φW (5.111) the amount of “secondorder”ness to avoid this behaviour. But instead of an explicit switch, we seek to achieve this by enforcing certain requirements on our discretization that can guarantee nonoscillatory solutions. These requirements are inspired by some well known mathematical properties of the scalar convection equation which we will discuss next. 5.6.5 Monotonicity and TVD Concepts
We saw in the discussion of the diffusion equation that it was important for the discrete equations to reﬂect the properties of the original transport equation so that we always obtain physically plausible solutions, even on ﬁnite sized meshes. This led us to identify some criteria, such as positivity of the coefﬁcients, that we could require our discretizations to satisfy. Unfortunately, for the convection equation the situation is not so straightforward and no deﬁnitive criteria can be identiﬁed. This is specially true in the case of nonlinear equations such as the momentum transport or Euler equations. Instead, there are several different conditions that one can impose on the discretizations and all of them have some merits and demerits. We will brieﬂy describe some of the most commonly used concepts in this section; detailed discussion of these ideas can be found in [?] and further references cited therein. Let us ﬁrst look at the monotonicity condition that we found useful for the diffusion equation. A linear discretization scheme, i.e, a scheme which can be expressed in the form 0 (5.112) φP Σanb φnb is said to be monotonic if all a nb are positive. It is possible to prove that monotonic schemes always result in physically plausible solutions. The reader can easily show that the ﬁrst order upwind scheme, Equation 5.84, is monotonic when it is stable. But all the higher order schemes we have discussed so far do not satisfy this condition. For example, the BeamWarming scheme, Equation 5.109, can be written as
¤ ¢ £ ¡ ¤ ¢ £ ¡ ¤ ¢ £ ¤ ¢ £ ¨ ¨ φP 1 ν 1 2 0 ν φP ν2 0 ν φW ν ν 2 0 1 φWW (5.113) and we see that the ﬁrst coefﬁcient is negative when 1 ν 2 and the third coefﬁcient is negative when 0 ν 1. Indeed, it can be shown that all monotonic schemes have to be ﬁrst order. For scalar convection equation monotonicity basically means that if we had two initial solutions φ0 and φ1 such that φ0 x 0 φ1 x 0 for all x, then at any subsequent time t , φ0 x t φ1 x t is also satisﬁed for all x. Therefore, unlike the diffusion equation, this condition by itself is therefore not very useful for discretizations of the convection term. A somewhat more useful but distinct concept is that of monotonicity preservation. This says that if the initial solution φ x 0 is a monotonically increasing function then the solution at any subsequent time is also monotonically increasing (and vice versa). There is a famous theorem by Godunov which shows that any linear monotonocity preserving scheme can be at most ﬁrst order accurate. For higher order accuracy therefore we need to devise nonlinear schemes even when the equation itself is linear. To get some guidance in how to go about devising such schemes we look at some other properties of the scalar convection equation. 117
¤ ¤ ! ! £ £ ¤ ! £ ¤! £ ¤! £ Both monotonicity and monotonicity preservation criteria described above do not address the case where the initial solution contains oscillations. To handle such cases we must introduce the concept of total variation. One of the properties of a scalar conservation law of the form 5.75 is that the total variation of any physically plausible solution, deﬁned as ∂φ TV φ dx (5.114) ∂x does not increase in time. For a discrete solution, the total variation is deﬁned as
' ¢ ' ¨ TV ∑ φP φW where the sum is over all the points of the domain. The total variation of a solution is a measure of the amount of “oscillation” in the solution. Numerical schemes which have the property that TV φ TV φ 0 (5.116) are said to be total variation diminishing(TVD). It can be shown that TVD schemes are also monotonicity preserving. Although TVD is a term one often encounters in the literature, most of the so called TVD schemes actually impose conditions that are more stringent than what is implied by Equation 5.116. This discussion of total variation presents an opportunity to bring up an important point about CFD which is that it is not always possible or desirable to use results and insights from the continuous solution on the discrete numerical solution. Thus while the continuous solutions of scalar equations are TVD, enforcing that on the numerical scheme may actually be detrimental in some cases. For example, in convecting any nonmonotonic proﬁle on a ﬁnite mesh the maxima (or minima) will not always lie at the center of a cell. The total variation should be allowed to increase and decrease as the maxima moves away and back to a cell center. Strict imposition of the TVD criterion can, however, clip the extrema at each time step. The above example demonstrates that TVD may be too strong a condition in some cases. On the other hand, since TVD only requires the total variation to diminish, it is possible, at least in principle, for a TVD scheme to allow creation of local oscillations of smaller magnitudes and so it may be a weak condition for some cases. A more appropriate condition to ensure physically plausible solutions is to require that the scheme create no new extrema and that existing extrema do not get ampliﬁed. This property has been termed range diminishing or local extremum diminishing (LED). Range diminishing schemes are also TVD. They also suffer from clipping errors. In fact, it is possible to show that because of the clipping errors, any range diminishing scheme is, at best, only secondorder accurate. All the conditions we have examined so far are based on physical properties of the scalar convection equation and are therefore suitable candidates for imposition on our numerical methods. However, it is difﬁcult to verify analytically whether a given numerical scheme meets any of these conditions. We mention one further condition, the upwind range condition which happens to be something we can directly enforce on a numerical scheme. The basic idea behind this condition is very simple and can be easily appreciated when we consider the transport of a scalar proﬁle such as the sine 118
¤ £ ( )¤ £ % % % % % &$"¤ #¨ % % % £ (5.115) or square wave shown in Figure 5.7. Suppose we have a uniform mesh and an explicit scheme such that the solution advances less than one cell width over a time step, i.e., a CFL number ν 1. We can see right away that for u 0 the solution in cell P at the current time step must be bounded by the solution over the cell W and P at the previous time step, min φ 0 x φP max φ 0 x (5.117) In terms of discrete values, this means that
¤ ! £ ( ( "¤ £ 0 0 £ xW x xP
0 0 xW x xP 0 0 min φW φP
! £ φP 0 0 max φW φP Of course, if the exact solution had an extrema between x W and xP at the previous time step, imposing Equation 5.118 will wind up clipping the extrema. But we can see that if our numerical scheme satisﬁed Equation 5.118, we can be assured that no new extrema can ever be generated. Before continuing we should point out the although strictly speaking the conditions we have derived (Equation 5.118), are based on the upwind range condition, in the literature they are often referred to as the TVD conditions. It can indeed be shown that schemes that satisfy the upwind range condition are TVD. But we should keep in mind that upwind range is a more stringent condition than required for TVD, i.e, it only represents a sufﬁcient condition and not a necessary one. We should also emphasize that these conditions are only meaningful for explicit schemes with ν 1 when applied to the linear advection equation. Strictly speaking the upwind range condition guarantee a nonoscillatory solution only for such cases. Nevertheless, this condition and schemes we shall derive based on it are often used for implicit and explicit schemes with ν 1 and for nonlinear equations as well. To see how the upwind range conditions can be imposed or veriﬁed on a given scheme, we start by noting that any explicit scheme can be written as
¡ ¨ φP 0 φP ΣaP δφ 0 f f
¢ ¨ where δφ f is the difference in φ across the face. For example, δφ e φE φP and δφw φP φW . The form of Equation 5.119 is required if the scheme is to identically preserve a constant solution. For a scheme on the 3 point stencil W P E , each δ f can potentially contribute to the discrete equations of two cells, viz., the ones on the left and right side of the face. Let a f and a f denote the corresponding coefﬁcients. With this notation, a three point scheme will have the form
2 ¡ 1 ¢ ¨ φP
0 φP
0 ae δφe 0 aw δφw It is easy to see that for u 0, the upwind range condition is satisﬁed by a scheme written in the form of 5.120 if at any face f ,
( 2 ( ! ¨ af
1 0 0
¨ 1 af 1 To demonstrate this we note that when a e
¢ ¨ 0, the scheme has the form
¤ 119 ¢ £ 2 φP 0 φP 0 aw φP 0 φW ¢ ¢ ¤ 2 ( ( )¤ (5.118) (5.119) 1 ¢ ¨ (5.120) (5.121) (5.122) 0 0 Now, when a w 0, φP φP and when a w 1, φP φW . For any other value of a w 0 and φ 0 . In other words, Equation 5.121 between zero and one, φ P must lie between φP W is a sufﬁcient condition for the discretization scheme to satisfy the upwind range condition. We can also show that this is a necessary condition by substituting Equation 5.120 into Equation 5.118 and rearranging to obtain
! 1 2 ! 0 0 We can see that for arbitrary δφ w and δφe this equation can only be satisﬁed if a e is identically zero and
! 2 ! For both positive and negative δφ 0 , this equation is satisﬁed only if 0 a w 1. Thus, f Equation 5.121 is both a necessary and a sufﬁcient condition for a three point scheme to satisfy the upwind range condition (Equation 5.118). We should also keep in mind that the upwind range condition also implies ν 1 and that Equation 5.121 was derived for u 0. The reader can show that for u 0 the conditions are reversed, i.e.,
¢ 3 1 3 ! ¨ af
2 0 0 af 1 If we express the ﬁrst order upwind scheme, Equation 5.84, in the form of Equa0 and a f ν . Thus, in its stability range of tion 5.120, we see that it implies a f 0 ν 1, the upwind scheme also meets the upwind range condition. It should be be emphasized that the stability and upwind range conditions come from different considerations, even though for the ﬁrst order upwind scheme the result happens to be the same. The LaxWendroff scheme, on the other hand, does not meet this condition even though it also is conditionally stable. This is easily seen by writing Equation 5.91 in the following form
6 ¤ ¡ ¡ ¤ ¨ We have seen that stability considerations require 0 ν 1 and the above equation tells us that for any value in that range the upwind range condition cannot be met. At this point we can begin to appreciate that as long as we have only a 3 point stencil we will not be able to formulate a second order accurate scheme that meets the upwind range condition and thus guarantees a physically plausible solution. We must therefore involve a wider stencil. The Beam Warming scheme, Equation 5.109, we saw earlier is an example of a second order scheme involving a wider stencil. Let us write it in terms of δφ f
6 ¤ ¤ ¨ This equation does not appear in the form of the general 3 point scheme, Equation 5.120, for which we have derived the upwind range conditions. However, we can rewrite it as follows 0 ν ν δφww 0 0 (5.128) φP φP 3ν 1ν δφw 0 2 2 δφw 120
¥ ¤ ¢ £ ¢ ¤ ¢ £ ¢ ¨ ¢ £ ¢ ¢ £ 4 ¢ φP 0 φP ν 3 2 0 ν δφw ν 1 2 0 ν δφww ( £ ¢ £ 4 5¢ φP 0 φP ν 1 2 0 ν δφe ν 1 2 0 ν δφw ( 2 ( ¤ ¢ £ ¨ 1 ( ¢ ¨ ( "¤ min 0
£ 0 δφw
¢ 0 aw δφw max 0 0 δφw (5.124) (5.125) (5.126) (5.127) 1 ¤ ¢ £ ( ¢ 2 ( )¤ min 0
£ 0 δφw
¢ 0 aw δφw 0 ae δφe max 0 0 δφw (5.123) 2 ¨ ¨ 2 ¨ ¨ 2 We can now think of the BeamWarming scheme as being in the form of a general 3 point scheme with the quantity in brackets as the coefﬁcient a w and with ae 0. We can also see that in general the Beam Warming scheme will not satisfy the upwind range conditions but this form now gives us a starting point to create schemes that do meet the condition. Clearly, we need to somehow ensure that the quantity in brackets is between 0 and 1. Before we look at exactly how to go about enforcing it, we note two important points. First, we observe that this quantity involves ratios of δφ f ’s across adjacent faces. We also note that in writing the original BeamWarming scheme which had a linear form (i.e., the coefﬁcients were not functions of φ ) in the form of Equation 5.128 we have made it look nonlinear. Of course, at this point Equation 5.109 and Equation 5.128 are identical but if we modify the coefﬁcient in brackets in any way we will make the scheme actually nonlinear. Although we have not proven it, Godunov’s theorem indeed necessitates the involvement of a stencil wider than 3 points in a nonlinear fashion if we want a higher order accurate scheme without oscillations. 5.6.6 Flux Limiters
In order to formulate an upwind range satisfying version of the BeamWarming scheme, we start by observing that the violations are due to the higher order terms we added in the proﬁle interpolation (Equation 5.96) since without them we would have the ﬁrst order scheme which is upwind range satisfying. It seems logical then to try and limit the contribution of these higher order terms. To this end, we write the BeamWarming scheme, Equation 5.127, in the form of a ﬁrst order scheme plus higher order terms :
¤ ¨ We can now write a limited form of the BeamWarming scheme by introducing coefﬁcients Ψ f that limit the contributions of the δφ ’s that appear in the higher order terms. This gives us
2 72 2 ¤ ¨ Here, once again we use the superscript when the difference across a face is involved in the equation for the cell downstream of the downstream cell and the superscript for the face difference when it appears in the equation for the immediate downstream cell. To clarify the nomenclature, it might be helpful to note that although it does not appear in the BeamWarming scheme, in a general ﬁve point stencil based scheme, the difference across the face e will also be involved and in that case it will have a coefﬁcient Ψ e . In order to formulate a general limiting procedure, we need to be able to write the coefﬁcients Ψ at any face as a function of some quantity deﬁned at that face. Once we have decided what quantity the coefﬁcients will be a function of, we can then determine the constraints on the function that are required for the scheme to be upwind range satisfying. We have seen before that the upwind range condition involves ratios of δφ ’s across the faces and it seems logical to have the limiter coefﬁcients be a function 121
2 2 72 1 ¢ ¢ £ ¢ ¢ φP 0 φP 0 ν δφw ν 1 2 ν 0 Ψw δφw ¢ ¢ £ ¢ ¢ φP 0 φP 0 ν δφw ν 1 2 ν 0 δφw 0 δφww 0 Ψww δφww ¨ 2 1 (5.129) (5.130) of these ratios. We identify the ratios of δφ in the forward and backward directions at each face as r and r such that at face w 2
¨ £ £ ¢ ( £ £ ¨ ¢ ¨ 2 1 2 rw
¨ £ ¤ £ 2 72 £ £ ¢ We now decide to write the coefﬁcients Ψ f in terms of these ratios so that
¤2 ¨ Ψf Ψf
¨ Ψ rf Ψ rf Ψ rf
£ and so on. The limited form of the BeamWarming scheme is then given by
¤ 82 2 ¤2 ¤ ¨ We now write this equation in the form of a 3 point scheme as
¤ 72 2 ¤2 To satisfy the upwind range condition, the quantity in brackets must be between 0 and 1. This will be the case if the function Ψ satisﬁes
( ¥ ¤ for all r and s. Considering the two inequalities in turn we can show that this condition implies the following two conditions
( 9¤ ¤ practice, a small number ε , is added to the denominators while computing r in order to prevent division by zero in the smooth regions where δφ 0.
2 In 122 ¤ £ @ ¢ ¤ Ψs ¢ Ψr r £ ¢ Ψr r Ψs ¢ £ 2 2 2 1 2 ν ν
(5.138) ¤ £ ¤ ¢ ¤ ¢ ¡ ( 0 ν1 1 ν Ψs 1 ν Ψr r ¥ ¤ 72 2 2 82 £ ¤ ¢ ¢ ¤2 £ ¤ £ ¡ ¢ ¨ φP 0 φP ν ν 1 2 2 82 Identifying the ratio 0 δφw 0 δφww as rww , we can write this as ν Ψ rw ν 1 2 ν Ψ rww rww 0 δφw 1 ¥ £ ¤ £ ¢ £ ¢ £ ¡ ¢ ¨ φP 0 φP ν ν 1 2 ν Ψ rw ν 1 2 ν Ψ rww £ ¢ £ ¢ ¢ ¢ φP 0 φP 0 ν δφw ν 1 2 ν 0 Ψ rw δφw ¤ 72 2 ¤1 Ψf
¨ 1 2 82 ¨ ¢ ¢ 2 ¨ 2 72 For the ww face, the ratio r is then given by rww φP φW φW φWW ¢ ¢ rw
1 φE φP φP φW φW φWW φP φW
¢ δφe δφw δφww δφw (5.131) δφw δφww (5.132) (5.133) 0 Ψ rww δφww (5.134) 0 δφww 0 δφw 0 δφw (5.135) (5.136) (5.137) ¨ ¤ £ Equations 5.139, 5.140 and 5.141 together give us sufﬁcient conditions for the limiter function such that when it is applied to the BeamWarming scheme we will obtain an upwind range satisfying second order accurate scheme. Figure 5.11 shows the region in which the BeamWarming scheme limiter function must lie for the upwind range condition to be met. Comparing Equations 5.129 and 5.130 we see that the original BeamWarming scheme corresponds to choosing Ψ r 1. The reader can show that by choosing Ψ r r in Equation 5.134 we obtain the original LaxWendroff scheme (Equation 5.91). Both of these choices for the Ψ function are also shown in Figure 5.11. We can see that the unlimited BeamWarming scheme will not meet the upwind range condition when r 0 5 while the LaxWendroff scheme will fail to do so when r 2 0. We also see that both schemes have Ψ 1 1. This is desirable for any limiter function as well since it implies full second order accuracy in smooth regions,
G G ¨ "¤ £ ¨ F¤ £ Since r is the ratio of successive differences of φ ’s, a negative value of r identiﬁes an extremum in φ . At such locations we would like to limit the scheme to ﬁrst order behaviour in order to avoid overshoots. This gives us one further condition for the limiter function Ψr 0 if 0 (5.139) With this condition imposed, the conditions Equation 5.138 are satisﬁed for all ν in the stability range of 0 ν 1 as long as the function satisﬁes ( Figure 5.11: Upwind range region associated with limiter functions for the BeamWarming function EEEEEEEEEEEEEEEEEEEEEEE ACBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBAD AC AC AC AC AC AC AC AC AC AC AC ACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE BABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABACB C C C C C C C C C C C CBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBDE ACBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBAD EEEEEEEEEEEEEEEEEEEEEEE DBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDE C CA CA CA CA CA CA CA CA CA CA CA CABCABCABCABCABCABCABCABCABCABCABCABCABCABCABCABCABCABCABCABCABCABCABAB ABBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBCBAD BABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABACB C C C C C C C C C C C CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE DE ACBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBAD EEEEEEEEEEEEEEEEEEEEEEE DBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDB CBC C C C C C C C C C C CBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBDE ABBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBAD ABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABACB EEEEEEEEEEEEEEEEEEEEEEE ACBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBAD AC AC AC AC AC AC AC AC AC AC AC DACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE BABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABACB C C C C C C C C C C C CBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBDE ACBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBAD EEEEEEEEEEEEEEEEEEEEEEE DBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDE ACBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBAD AC AC AC AC AC AC AC AC AC AC AC ACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACB EEEEEEEEEEEEEEEEEEEEEEE ACBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBAD AC AC AC AC AC AC AC AC AC AC AC ACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDBACDE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE BABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABACB C C C C C C C C C C C CBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBDE ACBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBAD EEEEEEEEEEEEEEEEEEEEEEE DBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDE C CA CA CA CA CA CA CA CA CA CA CA CABCABCABCABCABCABCABCABCABCABCABCABCABCABCABCABCABCABCABCABCABCABCABAB ABBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBCBAD BABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABACB C C C C C C C C C C C CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE CDE DE ACBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBAD EEEEEEEEEEEEEEEEEEEEEEE DBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDB CBC C C C C C C C C C C CBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBDE ABBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBAD ABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABACB EEEEEEEEEEEEEEEEEEEEEEE ACBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBAD AC AC AC AC AC AC AC AC AC AC AC DACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACBDACE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE BABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABABACB C C C C C C C C C C C CBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBDE ACBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBAD EEEEEEEEEEEEEEEEEEEEEEE DBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDE CBBBBBACBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBAD A AC AC AC AC AC AC AC AC AC AC ACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACB AC A AC AC AC A AC AC AC BCBBBBCBBBBACBBBBBBBBBBBBBBBBBBBBBBBBBBBAD AC AC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE CC CCC CC C C CBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBCBDE ABABABABABABABABABBABABABABABABABABABABABABABABABABABABABABABABABABABACB BBCBBBBCBBBACBBBBBBBBBBBBBBBBBBBBBBBBBBBAD EEEEEEEEEEEEEEEEEEEEEEE DBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDE AC AC AC A AC AC AC A AC A AC BBBCBBBBCBBCBBACBBBBBBBBBBBBBBBBBBBBBBBBBAD ACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACBACB CCCC CCC C ABABABABABABABABABABABABDE DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC DEAC BBBBCBBBBCBBCBCBACBBBBBBBBBBBBBBBBBBBBBBBBAD BDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDE EEEEEEEEEEEEEEEEEEEEEE DEBBBBBBBBBBBBBBBBBBBBBBBBD EEEEEEEEEEEEEEEEEEEEEEE BDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDE DBBBBBBBBBBBBBBBBBBBBBBBBD DEBBBBBBBBBBBBBBBBBBBBBBBBD DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE DE EEEEEEEEEEEEEEEEEEEEEEE DBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDBDE BBBBBBBBBBBBBBBBBBBBBBBBD 1 2 Ψ (r) 0.5 1.0 1.5 Ψr ( )¤ £ ( ( ¨ ¤ £ and Ψr 2.0 123
( )¤ £ 2r 2
2.5 Lax−Wendroff Ψ (r)=r 3.0 Beam−Warming Ψ (r)=1 r (5.140) (5.141) i.e., where the ratio r is nearly unity. In nonsmooth regions the limiter function is required to adaptively reduce the inﬂuence of the higher order terms so that the scheme is upwind range satisfying and at extrema it should make the scheme fully ﬁrst order. We have derived conditions for the limiter function by considering the BeamWarming scheme. We can apply the same procedure for the LaxWendroff scheme by ﬁrst writing it as a upwind scheme plus higher order terms and then applying limiter functions to the higher order terms to obtain the following limited form
¨ ¤ 1 ¤ 2 ¤ ! H ¤ £ ! ¨ The reader can show that for this scheme to meet the upwind range conditions, the constraints on the limiter function Ψ are the same as those we obtained for the BeamWarming scheme, viz., Equations 5.139, 5.140 and 5.141 which can be summarized as 0 Ψr min 2r 2 (5.143) It is interesting to note here that although the unlimited LaxWendroff scheme has a 3 point stencil, the ratios r w and re involved in the limiters do indeed have the effect of enlarging the stencil, as required by the Godunov theorem. Another condition that we can identify is that the limiter function should have a symmetry property such that 1 Ψr rΨ (5.144) r Although not apparent for the BeamWarming scheme (which is based on an upwinded derivative), this condition is necessary when one considers schemes based on central differenced derivative such as the Fromm scheme, Equation 5.99. It ensures that forward and backward differences get treated in the same way. We noted previously that any 5 point second order scheme could be written as a linear combination of the LaxWendroff and BeamWarming schemes. This property leads to the realization that universal limiter functions, applicable to any secondorder scheme, can be devised by further constraining them between the Ψ 1 and Ψ r lines. This region is shown in Figure 5.12. Of course, an inﬁnite number of functions can be devised that lie in the region. We look at a few examples next. The lower limit of the region is described by the following function
£ ¨ P¤ I ¨ ¤ £ 1 2 ¤ ! £ ( )¤ £ ( Ψr
¤ min r 1
¨ if r 0
( 0 (5.145) Ψr
£ 0 if r This function is one instance of the minmod function, which is deﬁned as a function that is equal to the smallest of its arguments if all the arguments have the same sign and is equal to zero otherwise. This limiter function is therefore known as the minmod limiter. The upper limit can be described by the function
R ¦¤ ! £ !¤ ! Q ¨ "¤ Ψr
£ max 0 min 2r 1 min r 2 (5.146) This limiter function is known as the Superbee limiter. 124 ¨ £ ¢ £ ¢ £ ¢ £ ¢ φP 0 φP 0 ν δφw ν 1 2 ν 0 Ψ rw δφw 0 Ψ re δφe (5.142) f f fffffffffffffffffffffffffffffffffffffffff `f``f``````````````````````````````````````````f iiiiiiiiii cccccccccc h`h`h`h`h`h`h`h`h`h` ``````````hi b`b`b`b`b`b`b`b`b`b` ``````````bc g``g`g`g````````````````````````````````````hi`````g ````````````````````````````````````````bc`````Ya g g g g g g g g g g g g g g g g g g g g g g g g g g hig hig hig hig hig hig hig hig hig hig g g g g g Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya Ya bcYa bcYa bcYa bcYa bcYa bcYa bcYa bcYa bcYa bcYa Ya Ya Ya Ya Ya hi hi hi hi hi hi hi hi hi hi ``````````hi bc bc bc bc bc bc bc bc bc bc ``````````bc hi hi hi hi hi hi hi hi hi hi ``````````hi bc bc bc bc bc bc bc bc bc bc ``````````bc iiiiiiiiii cccccccccc h`h`h`h`h`h`h`h`h`h` ``````````hi b`b`b`b`b`b`b`b`b`b` ``````````bc iiiiiiiiii cccccccccc h`h`h`h`h`h`h`h`h`h` ``````````hi b`b`b`b`b`b`b`b`b`b` ``````````bc hi hi hi hi hi hi hi hi hi hi ``````````hi bc bc bc bc bc bc bc bc bc bc ``````````bc hi hi hi hi hi hi hi hi hi hi ``````````hi bc bc bc bc bc bc bc bc bc bc ``````````bc hi hi hi hi hi hi hi hi hi hi ``````````hi bc bc bc bc bc bc bc bc bc bc ``````````bc hi hi hi hi hi hi hi hi hi hi ``````````hi bc bc bc bc bc bc bc bc bc bc ``````````bc pq`````````````````````````pq``````````hi d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d` b`b`b`b`b`b`b`b`b`b` `````````````````````````de ``````````bc `p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`h`h`h`h`h`h`h`h`h`h` qqqqqqqqqqqqqqqqqqqqqqqqiiiiiiiiii eeeeeeeeeeeeeeeeeeeeeeeeecccccccccc pq`````````````````````````pq d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d` `````````````````````````de `p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p` qqqqqqqqqqqqqqqqqqqqqqqq eeeeeeeeeeeeeeeeeeeeeeeee pq`````````````````````````pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq de de de de de de de de de de de de de de de de de de de de de de de de de `````````````````````````de pq`````````````````````````pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq de de de de de de de de de de de de de de de de de de de de de de de de de `````````````````````````de pq`````````````````````````pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq de de de de de de de de de de de de de de de de de de de de de de de de de `````````````````````````de pq`````````````````````````pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq de de de de de de de de de de de de de de de de de de de de de de de de de `````````````````````````de pq`````````````````````````pq d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d` `````````````````````````de `p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p` qqqqqqqqqqqqqqqqqqqqqqqq eeeeeeeeeeeeeeeeeeeeeeeee pq`````````````````````````pq d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d`d` `````````````````````````de `p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p`p` qqqqqqqqqqqqqqqqqqqqqqqq eeeeeeeeeeeeeeeeeeeeeeeee pq`````````````````````````pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq de de de de de de de de de de de de de de de de de de de de de de de de de `````````````````````````de pq`````````````````````````pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq de de de de de de de de de de de de de de de de de de de de de de de de de `````````````````````````de pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq pq `````````````````````````pq de de de de de de de de de de de de de de de de de de de de de de de de de `````````````````````````de
1 2 0.5 1.0 1.5 2.0 2.5 3.0 ST ST ST ST ST ST ST ST ST ST ST ST ST STBSTBSTBSTBSTBSTBSTBSTBSTBSTBSTB BBBBBBBBBBBBBBBBBBBBBBBBSU UV UV UV UV UV UV UV UV UV UV UV T T T T T T T T T T T T T SBSBSBSBSBSBSBSBSBSBST TTTTTTTTTT SBSBSBSBSBSBSBSBSBSBSBSBSBUVBUVBUVBUVBUVBUVBUVBUVBUVBUVBUVB BBBBBBBBBBBBBBBBBBBBBBBBUS ST ST ST ST ST ST ST ST ST ST ST ST ST STBSTBSTBSTBSTBSTBSTBSTBSTBSTBSTB BBBBBBBBBBBBBBBBBBBBBBBBSU UV UV UV UV UV UV UV UV UV UV UV ST ST ST ST ST ST ST ST ST ST ST ST ST STBSTBSTBSTBSTBSTBSTBSTBSTBSTBSTB BBBBBBBBBBBBBBBBBBBBBBBBSU UV UV UV UV UV UV UV UV UV UV UV T T T T T T T T T T T T T SBSBSBSBSBSBSBSBSBSBST TTTTTTTTTT SBSBSBSBSBSBSBSBSBSBSBSBSBUVBUVBUVBUVBUVBUVBUVBUVBUVBUVBUVB BBBBBBBBBBBBBBBBBBBBBBBBUS ST ST ST ST ST ST ST ST ST ST ST ST ST STBSTBSTBSTBSTBSTBSTBSTBSTBSTBSTB BBBBBBBBBBBBBBBBBBBBBBBBSU UV UV UV UV UV UV UV UV UV UV UV ST ST ST ST ST ST ST ST ST ST ST ST ST STBSTBSTBSTBSTBSTBSTBSTBSTBSTBSTB BBBBBBBBBBBBBBBBBBBBBBBBSU UV UV UV UV UV UV UV UV UV UV UV T T T T T T T T T T T T T SBSBSBSBSBSBSBSBSBSBST TTTTTTTTTT SBSBSBSBSBSBSBSBSBSBSBSBSBUVBUVBUVBUVBUVBUVBUVBUVBUVBUVBUVB BBBBBBBBBBBBBBBBBBBBBBBBUS ST ST ST ST ST ST ST ST ST ST ST ST ST STBSTBSTBSTBSTBSTBSTBSTBSTBSTBSTB BBBBBBBBBBBBBBBBBBBBBBBBSU UV UV UV UV UV UV UV UV UV UV UV T T T T T T T T T T T T T SBSBSBSBSBSBSBSBSBSBST TTTTTTTTTT SBSBSBSBSBSBSBSBSBSBSBSBSBUVBUVBUVBUVBUVBUVBUVBUVBUVBUVBUVB BBBBBBBBBBBBBBBBBBBBBBBBUS TTTTTTTTTTTTTT WXBBBBBBBBBBBBBBBBBBBBBBBBBBWX BBBBBBBBBBS WX WX WX WX WX WX WX WX WX WX WX WX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST BWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWB ST ST ST ST ST ST ST ST ST ST X X X X X X X X X X X X BX BX BX BX BX BX BX BX BX BX BX BX BX B WXBBBBBBBBBBBBBBBBBBBBBBBBBBWX BBBBBBBBBBS TTTTTTTTTTTTTT SBSBSBSBSBSBSBSBSBSBSBSBSBSBT T T T T T T T T WXBBBBBBBBBBBBBBBBBBBBBBBBBBWX BBBBBBBBBBS WX WX WX WX WX WX WX WX WX WX WX WX BWX BWX BWX BWX BWX BWX BWX BWX BWX BWX BWX BWX BWX BSBSBSBSBSBSBSBSBSBST TTTTTTTTTTTTTT WXBBBBBBBBBBBBBBBBBBBBBBBBBBWX BBBBBBBBBBS WX WX WX WX WX WX WX WX WX WX WX WX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST BWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWB ST ST ST ST ST ST ST ST ST ST X X X X X X X X X X X X BX BX BX BX BX BX BX BX BX BX BX BX BX B WXBBBBBBBBBBBBBBBBBBBBBBBBBBWX BBBBBBBBBBS TTTTTTTTTTTTTT SBSBSBSBSBSBSBSBSBSBSBSBSBSBT T T T T T T T T WX WX WX WX BBBBWXBBBBBBBBBBBBBBBBBBBBBBWX BBBBBBBBBBS WX WX WX WX WX WX WX WX BWX BWX BWX BWX BWX BWX BWX BWX BWX BWX BWX BWX BWX BSBSBSBSBSBSBSBSBSBST XXXX XXXX X X X ST X ST X ST X ST X ST X ST X ST X ST X ST X ST X ST X ST X ST X ST WBWBWBWBWBWBWBWBWBBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWB ST ST ST ST ST ST ST ST ST ST BBBBXBBBBBWXBBBBBBBBBBBBBBBBBWX BBBBBBBBBBS TTTTTTTTTTTTTT SBSBSBSBSBSBSBSBSBSBSBSBSBSB XXXX XXX X X X BX BX BX BX BX BX BX BX BX BX BX BX BX BT T T T T T T T T WXBBWBWBWBWXBBWBWBXBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWBWB SBSBSBSBSBSBSBSBSBST BWBBBBBWBBBWBBBBBBBBBBBBBBBBBWX BBBBBBBBBBS TTTTTTTTTTTTTT WX WX W WX WX WX WX W WX W WX BBXBBBBBXBBXBBWXBBBBBBBBBBBBBBBWX BBBBBBBBBBS WX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBWX SBST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST BX BX BX BX BX BX BX BX BX BX BX BX BX B XXX XXXX X WBWBWBWBWBWBWBWBWBWBWBWBBWBWBWBWBWBWBWBWBWBWBWBWBWB ST ST ST ST ST ST ST ST ST ST BBBXBBBBBXBBXBXBWXBBBBBBBBBBBBBBWX BBBBBBBBBBS BSBSBSBSBSBSBSBSBSBSBSBSBSBSBSBSBSBSBSBSBSBSBST TTTTTTTTTTTTTTTTTTTTTT STBBBBBBBBBBBBBBBBBBBBBBBBS STBBBBBBBBBBBBBBBBBBBBBBBBS ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST STBBBBBBBBBBBBBBBBBBBBBBBBS ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST ST TTTTTTTTTTTTTTTTTTTTTTT SBSBSBSBSBSBSBSBSBSBSBSBSBSBSBSBSBSBSBSBSBSBSBST BBBBBBBBBBBBBBBBBBBBBBBBS Ψ (r) Figure 5.12: Upwind range region for all second order upwind schemes 1 2 Ψ (r) 0.5 Figure 5.13: minmod and Superbee limiters 1.0 1.5 2.0 125
r
1 2 2.5 Lax−Wendroff Ψ (r) Ψ (r)=r 0.5 3.0 1.0 Beam−Warming Ψ (r)=1 1.5 r 2.0 2.5 3.0 r 1 2 minmod Exact Solution
0.5 1.5 minmod Exact Solution 1 φ 0 φ
0.5 0.5 0 1 0 0.2 0.4 0.6 0.8 1 0.5 0 0.2 0.4 0.6 0.8 1 x/L x/L Figure 5.14: Solutions using the BeamWarming scheme with the minmod limiter for linear convection of (a) sine wave (b) square wave To study the characteristics of these limiters we use them in the BeamWarming scheme for the two scalar convection problems we used earlier. Figure 5.14 shows the resulting proﬁles after 25 time steps, using a CFL number ν 0 5 for the minmod scheme. For the sine wave proﬁle, the results are certainly less diffuse than the ﬁrst order upwind results (Figure 5.8) but still more diffuse than the unlimited BeamWarming scheme in the smooth gradient regions. Of course, the main improvement is that the oscillations are now gone completely. Note also how the application of limiters, has “clipped” the peaks of the sinewave proﬁle. As we noted earlier, this is a direct consequence of enforcing the TVD conditions on our numerical scheme. This clipping is a different phenomenon from the diffusion introduced by the ﬁrst order upwind scheme which reduces the amplitude uniformally and not just at the peaks. We see some interesting results when we apply the Superbee limiter. Of course, no new extrema are created for either proﬁles, but for the sine wave, Figure 5.15(a) shows an artiﬁcial sharpening of the smooth proﬁles in some regions. This is clearly not desirable for smooth solutions but this overcompressive property does help the Superbee limiter capture φ discontinuities much better, as seen in Figure5.15. The usual way of comparing how well two schemes perform for discontinuities is to count the number of cells over which the discontinuity is resolved. We see that the Superbee limiter in this example contains it over within a four cell region while the minmod limiter causes it to spread over eight cells. We can understand the limiter behaviour by thinking in terms of artiﬁcial viscosity. We saw at the beginning of the section that the upwind scheme has signiﬁcant artiﬁcial diffusion while the central diffusion has none. We can improve the performance of the ﬁrst order upwind scheme by adding some negative diffusion but we cannot add too much of it or else we make the scheme unstable. The constraints on the limiter functions shown in Figure 5.12 basically tell us how much is acceptable. At the lower 126
G ¨ 1 2 Superbee Exact Solution
0.5 1.5 Superbee Exact Solution 1 φ 0 φ
0.5 0.5 0 1 0 0.2 0.4 0.6 0.8 1 0.5 0 0.2 0.4 0.6 0.8 1 x/L x/L Figure 5.15: Solutions using the BeamWarming scheme with the Superbee limiter for linear convection of (a) sine wave (b) square wave limit, the minmod limiter adds just enough so that the scheme is actually second order while the Superbee limiter adds the maximum amount that still keeps the scheme from producing overshoots. Thus, the minmod scheme ends up being somewhat diffusive and smoothens out the gradients while the negative viscosity implied by the Superbee limiter has the opposite effect of sharpening even what should have been smooth proﬁles. In practice, a limiter function that lies somewhere in between usually offers a good compromise for general applications. Several different limiters have been proposed by researchers during the last two decades. We will look at a few of them in the next section. 5.6.7 Considerations For SteadyState Applications
If we are interested in only the steady state solution, in principle we only need to consider the spatial discretization. Consider the semidiscretized version of the general 3 point scheme (Equation 5.120)
2 ¡ 3 3 1 ¨ ¢ ( 1 ( % % % % ∂φ ∂t 0 ae δφe 0 aw δφw (5.147) P Suppose the coefﬁcients at all faces satisy the following condition known as the positivity condition af 0 af 0 (5.148) Now if φP is a local maxima, δφe 0 and δφw 0. Therefore the left hand side side of Equation 5.147 is negative and thus d φ P 0. This means that in principle a local dt maxima will not increase and similarly a local minima will not increase. This condition is also therefore known as local extremum diminishing (LED) condition and sometimes even as the TVD condition. 127
' 2 2 1.5 Ψ 1 0.5 minmod Superbee van Leer van Albada
0 1 2 3 4 0 r Figure 5.16: van Leer and van Albada Limiters 2 1.5 Ψ 1 0.5 minmod Quadratic Cubic 0 0 1 r 2 3 4 Figure 5.17: Quadratic and cubic approximations to the mimnod limiter 128 Comparing the positivity condition, Equation 5.148, with the upwind range condition, Equation 5.121, we see that the former is less restrictive. Consequently the range of limiter functions permitted is larger. The general procedure for determining the allowed range is similar to the one we used to derive the conditions for the limiter functions for the explicit BeamWarming scheme. We start with the scheme written in the form of a ﬁrst order scheme plus higher order correction and then introduce a coefﬁcient that limits the higher order contribution. For example, the limited form of second order value for the face φ as given by Equation 5.96 is written as
% ¢ % % % ¢ ¡ ¨ r £ 1 ¨ ¢ ¡ φe φP Ψe ∆x ∂ φ 2 ∂x (5.149)
P As before, we choose the coefﬁcient as a function of the ratio of δφ ’s so that
% % % % ¤1 2 ¡ ¨ φe To get a fully upwinded scheme we use Equation 5.100 to evaluate the derivative ∂ φ P ∂x and write ∆x φe φP Ψ re φ φW (5.151) 2P With this expression for face values we get the semidiscretized version of the limited explicit BeamWarming scheme
s ¤ 82 2 ¤2 P
% % % The ratios rw and rww are as deﬁned by Equations 5.131 and 5.132. Note also that in deriving the above equation we have made use of the following identity implied by Equation 5.131 1 (5.153) re rw The reader can show that for the scheme to meet the positivity condition the limiter must satisfy
¨ t¤ 2 82 2 in addition to having the symmetry property, Equation 5.144. Thus, in principle, any symmetric function valued less than 2 can be used under the positivity condition for steady state applications. In practice, however, we need to use some form of time marching or iterative scheme even if only the ﬁnal solution is desired and in that case the positivity condition is not sufﬁcient to ensure nonoscillatory solutions. Most codes therefore still use limiter functions that lie within the region shown in Figure 5.12. An additional consideration for steadystate applications is the ability of the scheme to converge. Because the use of limiters makes the scheme nonlinear, the interaction between the limiters and the solution advancement scheme can cause the solution to 129 ( u¤ Ψr Ψr
£ £ 0 2 if r 0 (5.154) £ £ ¢ ¨ % ∆x ∂φ ∂t u δφw u Ψ rw δφw 2 Ψ rww δφww (5.152) ' ¤ £ ¤1 £ φP Ψ re ∆x ∂ φ 2 ∂x (5.150)
P ¨ get caught in limit cycles. Both the minmod and Superbee limiters we have seen have slope discontinuities which exacerbate these problems. For example, in the case of the minmod limiter the value of Ψ r is constant for all values of r greater than one but changes suddenly if r falls below 1. This could cause the solution to change in response such that r becomes larger than 1 in the next iteration and this can in turn cause solution to change back to r greater than one. Differentiable, i.e., smoother limiter functions can alleviate such problems. One such function is given by 2r Ψr (5.155) 1r Figure 5.16 shows that this function lies between the minmod and Superbee limiters. This function was effectively the limiter function in the van Leer scheme although the scheme was itself described in 1974, long before the general concepts of TVD and ﬂux limiters were actually formulated. An even more smoother function, also shown in Figure 5.16, is the van Albada limiter r2 r (5.156) 1 r2 When designing limiter functions for general purpose codes, one has to tradeoff accuracy against solution stability. We have seen that the lower limit, represented by the minmod function, is relatively more diffusive. More diffusivity generally results in better convergence properties. The minmod function is therefore often an appropriate conservative choice for limiter function when the governing equations are nonlinear and being solved on poor quality meshes. To get a differentiable limiter, we can approximate the minmod function by the following function
( ¡ ¨ ¤ Ψr
£ This function is shown in Figure 5.17. Note that the function does not pass through (1,1) indicating that we do not quite achieve full second order accuracy. A somewhat better choice is a cubic approximation, given by
( ¡ ¨ t¤ Ψr
£ 4r r3 4 r2 r3 1 r2
¡ ¡ This is the limiter function used in the Fluent code. Both the quadratic and cubic limiters have a slope discontinuity at r 2 but it is much milder than the one in minmod function at r 1. In the literature, one often encounters limiter functions expressed as a function of b two variables rather than the single variable r we have used. For example, with r a, for positive a and b the cubic limiter function can be written as
( ¨ v¤ Ψab
! £ 1
¨ if 2a b b b 2a b2 a b 2a
¤ ¤ ¡ ¡ £ £ ¡ (5.159) otherwise This form avoids possible division by zero in regions where φ is constant. 130 ¨ ¡ ¡ 2 1
¨ ¡ ¨ ¨ t¤ Ψr
£ 2r r2 r r2 r2 ¡ ¡ ¨ ¤ £ ¨ ¤ £ r 2 (5.157) r 2 (5.158) ¨ 5.6.8 Geometrical Interpretation of Limiting
To gain some physical insight into the ratios r and the action of the limiters let us consider the evalutaion of the face value φ e using Equation 5.151. In the upwind scheme if u 0, we decide to linearly extrapolate the value of φ e from the cell P but we have two choices for the slope we can use  either φ E φP or φP φW . We can now see that the ratio re is the ratio of these two possible slopes. In an unlimited scheme the choice of the slope is ﬁxed but the limiting process switches between these choices in an adaptive manner depending on the variation of φ in the vicinity of P. Figure 5.18 shows four possible scenarios for the variation of φ in the neighbourhood of cell P. The dots represent the values at the cell centers and the crosses the interpolated value at the face e. The heavy lines show the possible slopes and the light lines show the slopes used by the limiting process. If the variation between E , P and W is exactly linear as shown in Figure 5.6.8 then r e 1. In this case either choice of slope gives us the same value of φe which is also the exact result for a linear variation. Thus for r 1 we want the limiter function Ψ to have a value of 1. If the function φ is varying faster in the interval E P than in P W , then r 1. In this case, if we use the slope φP φW , we certainly won’t cause an overshoot in the face value as seen in Figure 5.6.8 so a value of ψ 1 is also safe. In fact, a higher value upto r 2 can be used since we know that the cell value at E is certanly larger than the extrapolated cell value we will get. However, if φ variation is slower in E P than in P W , i.e., r 1, then using the slope φ P φW will overestimate the face value φ e , as shown in Figure 5.6.8. The proper choice in this case is to use the slope φ E φP for the extrapolation. Looking at Equation 5.151, we see that choosing Ψ r does exactly that. Finally, if φP is a local maxima, then r 0. In this case, as Figure 5.6.8 show, the only safe choice to set φ e φP , i.e., Ψ 0. 5.6.9 Considerations for MultiDimensional Problems
The general theoretical framework for multidimensional convections is not well developed. For example, there is no easy extension of the TVD concept in two or three dimensions. There are additional diffuculties due to the fact that the waves can propagate in inﬁnite number of directions and therefore unlike the one dimensional case there is not always a cell in upwind direction. The situation is even more complicated for systems of equations as we will see in later chapters that deal with Euler and Navier Stokes equations. In practice, most codes deal with multidimensional convection by using only the 1D concepts we have discussed in this chapter. In case of Cartesian grids these are independently applied along each of the coordinate directions. For bodyﬁtted structured grids, the transformed coordinates provide the directions. For example, for the cell shown in Figure ??, a higher order value for the face e can be written just as it was for the 1D case, i.e., using Equation 5.151 (assuming of course that the face normal velocity is positive). Similarly the value at face n (again for the case of positive normal velocity) is given by ∆η φn φP Ψ rn φP φS (5.160) 2 131
¤ ¢ £ ¤1 £ ¡ ¨ ¤ ¢ ¢ ¨ £ ¨ ¤ ¢ ¢ £ ¤ ¢ ¢ ¨ £ ¢ ¤ ¨ £ ¢ ¨ £ 1 ¨ 1 ¨ ¢ φ Ψ=1 φ Ψ=2 Ψ=1 W w P e
1
@ E x W w P e E x φ Ψ=1 φ Ψ= r W w P e E x W w P e Figure 5.18: Geometrical interpretation of the limiters 132 x x (c) r 1 (d) r w (a) r (b) r 11 Ψ=0 E x
0 A f V f ∆r C1 0 C0 Figure 5.19: Schematic for SecondOrder Scheme where
¨ £ ¢ For unstructured grids we can also use the 1D ideas at each face in the local normal direction. The main difference is that we cannot use formulas like Equation 5.97 or 5.100 to evaluate the gradient ∇φ 0 . Instead we must use one of methods for the calculation of the cell centered gradient that we dicussed in the last section. Using such a gradient we can now write the face value at face f in Figure 5.19
y ' ¤ ¤1 ¡ ¨ Another difference for unstructured meshes is that we cannot deﬁne the ratio r f in terms of a cell value in the upwind direction since no unique cell exists. However, since we already have the cell gradient we can construct an imaginary point U as shown in Figure 5.19 which is located symmetrically in opposite direction from C0 to C1. The ratio r is then given by φ1 φ0 rf (5.163) ∇ φ 0 ∆r10 With this deﬁntion the limiting coefﬁcients can be evaluted using the same limiter functions that we used for 1D cases.
¤ ¢ £ ¨ 1 5.6.10 Implementation Issues
If iterative solvers are used to solve the resulting set of discrete equations, it is important to ensure that the Scarborough criterion is satisﬁed by the nominally linear equations presented to the iterative solver. Consequently, typical implementations of higherorder schemes use deferred correction strategies whereby the higher order terms are included 133 £ φf φ0 Ψ rf ∇ φ ¢ ¢ rn
1 φN φP φP φS (5.161) ¤ £ ∆r f 0 (5.162) ...
View
Full Document
 Fall '10
 NA
 Numerical Analysis, pq pq pq, st st st, MUSCL scheme

Click to edit the document details