ho-conv

ho-conv - the equation that is satisfied at convergence is...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: the equation that is satisfied at convergence is given by ¨¥ ©§¤ ¥ ¦¤ The consequence of this is that the final answer we obtain depends on the time-step 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 final 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 Higher-Order Schemes We have seen thus far that both the first order upwind and central difference schemes have severe limitations, the former due to artificial 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 higher-order interpolation. These higher-order schemes aim to obtain at least a second-order 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 profile of φ is essentially constant. That is, for Fe 0, Instead of a constant profile assumption for φ , we may use higher-order profile assumptions, such as linear or quadratic, to derive a set of upwind weighted higher-order 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 Second-Order Upwind Flux Expressions We may derive a second-order accurate expressions for the face flux by making a linear profile assumption. This is equivalent to retaining the first 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 second-order 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 second-order fully upwind expression for face flux  ¢ Schemes using this form of face flux are referred to in the literature as Beam-Warming schemes. 5.6.2 Third-Order Flux Expressions Third-order accurate expressions for the face flux can be derived by retaining both the first and second derivative in the Taylor series expansion, Equation 5.95: ¤ ¢ £ and using cell-centroid values to write the derivatives  and we may write  ¡ ¨ Re-arranging, 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 third-order 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 fluxes 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 flux expressions given by Equation 5.108 with the first order backward Euler differencing for the unsteady term that we employed earlier we will find 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 Lax-Wendroff scheme described in Section 5.5.5. It is possible to formulate a general procedure whereby any of the second order flux 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 Lax-Wendroff scheme, most of these schemes also result in time-step dependence in the steady state solution. As we noted earlier, this occurs because the spatial profile 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 Runge-Kutta method used for ODEs) or an implicit temporal discretization. However, the resulting methods are more difficult 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 flux expressions but present one scheme that results from the fully upwinded second order expression for face flux, Equation 5.101. The scheme is known 114 ¨ Here, κ 1 yields the Beam-Warming 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 Beam-Warming Exact Solution 1.5 Beam-Warming 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: Beam-Warming solutions for linear convection of (a) sine wave (b) square wave as the Beam-Warming second-order 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 profiles the second order upwind flux expression is clearly much less diffusive than the first 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 Beam-Warming solutions shown in Figure 5.10 are almost complimentary to the solutions obtained using the Lax-Wendroff scheme (Figure 5.9). For symmetric profiles, 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 five 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 flux limiters applicable to any second order scheme. To summarize, we have seen so far that the first order upwind scheme produces physically plausible results but is very diffusive. The second order central-difference 1 There is an implicit variant of the scheme which is also known as Beam-Warming scheme but for the purposes of this chapter we will only be considering the explicit version and refer to it as the Beam-Warming 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 von-Neumann sense when used with first order explicit time stepping but even when they can be made stable, they produce non-physical oscillations in the solution. Attempts at formulating upwind-biased (i.e, non-symmetric) higher order schemes by improving the interpolations used in the first 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 first 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 first one involves keeping the discretization as it is but damping out the oscillations produced by introducing artificial 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 sufficient 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 coefficients ε 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 ad-hoc 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 upwind-biased schemes. Such schemes do have some dissipative character but the higher order derivatives used in constructing the face fluxes 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 “second-order”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 non-oscillatory 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 reflect the properties of the original transport equation so that we always obtain physically plausible solutions, even on finite sized meshes. This led us to identify some criteria, such as positivity of the coefficients, that we could require our discretizations to satisfy. Unfortunately, for the convection equation the situation is not so straightforward and no definitive criteria can be identified. This is specially true in the case of non-linear 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 briefly 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 first 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 first 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 Beam-Warming 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 first coefficient is negative when 1 ν 2 and the third coefficient is negative when 0 ν 1. Indeed, it can be shown that all monotonic schemes have to be first 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 satisfied 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 first order accurate. For higher order accuracy therefore we need to devise non-linear 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, defined as ∂φ TV φ dx (5.114) ∂x does not increase in time. For a discrete solution, the total variation is defined 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 non-monotonic profile on a finite 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 amplified. 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 second-order 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 difficult 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 profile 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 satisfied 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 sufficient 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 non-oscillatory 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 non-linear equations as well. To see how the upwind range conditions can be imposed or verified 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 coefficients. 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 satisfied 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 sufficient 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 satisfied if a e is identically zero and ! 2 ! For both positive and negative δφ 0 , this equation is satisfied only if 0 a w 1. Thus, f Equation 5.121 is both a necessary and a sufficient 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 first 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 first order upwind scheme the result happens to be the same. The Lax-Wendroff 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 Beam-Warming scheme as being in the form of a general 3 point scheme with the quantity in brackets as the coefficient 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 Beam-Warming scheme which had a linear form (i.e., the coefficients were not functions of φ ) in the form of Equation 5.128 we have made it look non-linear. Of course, at this point Equation 5.109 and Equation 5.128 are identical but if we modify the coefficient in brackets in any way we will make the scheme actually non-linear. 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 Beam-Warming scheme, we start by observing that the violations are due to the higher order terms we added in the profile interpolation (Equation 5.96) since without them we would have the first 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 Beam-Warming scheme, Equation 5.127, in the form of a first order scheme plus higher order terms :  ¤ ¨ We can now write a limited form of the Beam-Warming scheme by introducing coefficients Ψ 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 Beam-Warming scheme, in a general five point stencil based scheme, the difference across the face e will also be involved and in that case it will have a coefficient Ψ e . In order to formulate a general limiting procedure, we need to be able to write the coefficients Ψ at any face as a function of some quantity defined at that face. Once we have decided what quantity the coefficients 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 coefficients 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 coefficients Ψ f in terms of these ratios so that ¤2 ¨ Ψf Ψf ¨ Ψ rf Ψ rf Ψ rf £ and so on. The limited form of the Beam-Warming 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 Ψ satisfies ( ¥ ¤ 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 sufficient conditions for the limiter function such that when it is applied to the Beam-Warming scheme we will obtain an upwind range satisfying second order accurate scheme. Figure 5.11 shows the region in which the Beam-Warming 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 Beam-Warming scheme corresponds to choosing Ψ r 1. The reader can show that by choosing Ψ r r in Equation 5.134 we obtain the original Lax-Wendroff scheme (Equation 5.91). Both of these choices for the Ψ function are also shown in Figure 5.11. We can see that the unlimited Beam-Warming scheme will not meet the upwind range condition when r 0 5 while the Lax-Wendroff 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 identifies an extremum in φ . At such locations we would like to limit the scheme to first 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 satisfied for all ν in the stability range of 0 ν 1 as long as the function satisfies ( 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 non-smooth regions the limiter function is required to adaptively reduce the influence of the higher order terms so that the scheme is upwind range satisfying and at extrema it should make the scheme fully first order. We have derived conditions for the limiter function by considering the BeamWarming scheme. We can apply the same procedure for the Lax-Wendroff scheme by first 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 Lax-Wendroff 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 Beam-Warming 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 Lax-Wendroff and Beam-Warming schemes. This property leads to the realization that universal limiter functions, applicable to any second-order 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 infinite 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 defined 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 Beam-Warming 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 Beam-Warming scheme for the two scalar convection problems we used earlier. Figure 5.14 shows the resulting profiles after 25 time steps, using a CFL number ν 0 5 for the minmod scheme. For the sine wave profile, the results are certainly less diffuse than the first order upwind results (Figure 5.8) but still more diffuse than the unlimited Beam-Warming 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 sine-wave profile. 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 first 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 profiles, but for the sine wave, Figure 5.15(a) shows an artificial sharpening of the smooth profiles 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 artificial viscosity. We saw at the beginning of the section that the upwind scheme has significant artificial diffusion while the central diffusion has none. We can improve the performance of the first 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 Beam-Warming 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 profiles. 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 Steady-State Applications If we are interested in only the steady state solution, in principle we only need to consider the spatial discretization. Consider the semi-discretized 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 coefficients 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 Beam-Warming scheme. We start with the scheme written in the form of a first order scheme plus higher order correction and then introduce a coefficient 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 coefficient 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 semi-discretized version of the limited explicit Beam-Warming scheme s ¤ 82 2 ¤2 P % % % The ratios rw and rww are as defined 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 final solution is desired and in that case the positivity condition is not sufficient to ensure non-oscillatory solutions. Most codes therefore still use limiter functions that lie within the region shown in Figure 5.12. An additional consideration for steady-state applications is the ability of the scheme to converge. Because the use of limiters makes the scheme non-linear, 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 flux 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 trade-off 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 non-linear 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 fixed 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 Multi-Dimensional Problems The general theoretical framework for multi-dimensional 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 infinite 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 multi-dimensional 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-fitted 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 Second-Order 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 define 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 defintion the limiting coefficients 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 satisfied by the nominally linear equations presented to the iterative solver. Consequently, typical implementations of higher-order 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

This note was uploaded on 12/29/2011 for the course ME 608 taught by Professor Na during the Fall '10 term at Purdue.

Ask a homework question - tutors are online