A Truncation Error Analysis of Third-Order MUSCL Scheme for Nonlinear Conservation Laws
aa r X i v : . [ phy s i c s . c o m p - ph ] J un A Truncation Error Analysis of Third-Order MUSCLScheme for Nonlinear Conservation Laws
Hiroaki Nishikawa ∗ National Institute of Aerospace, Hampton, VA 23666, USA
Abstract
This paper is a rebuttal to the claim found in the literature that the MUSCL scheme cannot be third-orderaccurate for nonlinear conservation laws. We provide a rigorous proof for third-order accuracy of the MUSCLscheme based on a careful and detailed truncation error analysis. Throughout the analysis, the distinctionbetween the cell average and the point value will be strictly made for the numerical solution as well as for thetarget operator. It is shown that the average of the solutions reconstructed at a face by Van Leer’s κ -schemerecovers a cubic solution exactly with κ = 1 /
3, the same is true for the average of the nonlinear fluxes evaluatedby the reconstructed solutions, and a dissipation term is already sufficiently small with a third-order truncationerror. Finally, noting that the target spatial operator is a cell-averaged flux derivative, we prove that the leadingtruncation error of the MUSCL finite-volume scheme is third-order with κ = 1 /
3. The importance of the diffusionscheme is also discussed: third-order accuracy will be lost when the third-order MUSLC scheme is used witha wrong fourth-order diffusion scheme for convection-diffusion problems. Third-order accuracy is verified bythorough numerical experiments for both steady and unsteady problems. This paper is intended to serve as areference to clarify confusions about third-order accuracy of the MUSCL scheme, as a guide to correctly analyzeand verify the MUSCL scheme for nonlinear equations, and eventually as the basis for clarifying third-orderunstructured-grid schemes in a subsequent paper.
Keywords : Advection-diffusion, Finite volume, Finite difference, Convection, Compressible flow, Viscous flows
The MUSCL (Monotonic Upstream-centered Scheme for Conservation Laws) approach developed by Van Leerin a series of papers [1, 2], has become one of the most successful approaches to achieving higher-order accuracyin finite-volume methods. MUSCL refers to a comprehensive approach to achieving second- and higher-orderaccuracy in a finite-volume method with a monotonicity preserving mechanism designed for accurately capturingdiscontinuous solutions. In the MUSCL approach, a finite-volume discretization is constructed with cell-averagedsolutions stored at cells as numerical solutions, where a higher order accurate flux is achieved at the cell face utilizinghigher order accurate solutions reconstructed at the face from left and right cells using these numerical cell-averagedsolutions. For the solution reconstruction, Van Leer proposed a one-parameter-family of solution reconstructionscheme, the κ -reconstruction scheme, where κ is a free parameter, and demonstrated that the finite-volume schemecan achieve third-order accuracy with κ = 1 / κ = 1 / κ = 1 / κ = 1 / ∗ Associate Research Fellow ([email protected]) κ = 1 / κ = 1 / κ = 1 / κ = 1 / κ = 1 / κ = 0 [30, 31]). The clarification is important for truly achieving third-order accuracy and takingfull advantage of some economical third-order unstructured-grid schemes. In seeking the clarification, we havefound that the confusion is largely rooted in the ever-present confusion over third-order convection schemes: theMUSCL scheme and the QUICK scheme [32]. This paper is the first paper in a series. The second paper will focuson the QUICK scheme.The paper is organized as follows. In Section 2, we discuss the integral form of a nonlinear conservation law asa basis for a finite-volume discretization and point out the difference between the cell averaged and the point valuesolutions and operators. In Section 3, we describe the MUSCL scheme with the κ -reconstruction scheme appliedto a nonlinear conservation law, and make it clear that the MUSCL scheme is based on cell-averaged solutions. InSection 4, we derive the truncation error of the MUSCL scheme and prove that third-order accuracy is achievedonly with κ = 1 /
3. In Section 5, we derive a family of diffusion schemes compatible with the third-order MUSCLscheme. In Section 6, we present a series of numerical test results to verify the analysis and demonstrate third-orderaccuracy in the cell-averaged solution, not in the point-valued solution unless it is accurately recovered from thecell-averaged solution. In Section 7, we conclude the paper with final remarks.2 i − x i − x i u i x i +1 x i +2 u ( x ) xu (a) Cell-averaged solutions and control volumes on a uniform grid. x i xu u i u i u ( x ) (b) The cell-averaged solution and thepoint-valued solution within the cell i . Figure 1: Cell-averaged solutions on a uniform grid and the comparison of the cell average and the point value ina cell i . Consider a one-dimensional conservation law in the differential form: u t + f x = s ( x ) , (1)where u is a solution variable, f is a flux and s ( x ) is a forcing term, and the subscripts t and x denote the partialderivatives with respect to time and space, respectively. To derive an integral form, we consider a one-dimensionalgrid with a uniform spacing h : x i +1 − x i = h , i = 1 , , , · · · , n , where n is an integer. See Figure 1(a). In thispaper, we focus on the interior of a domain and do not discuss any boundary effects. Then, we integrate thedifferential form (1) over a control volume around x = x i , x ∈ [ x i − / , x i +1 / ] = [ x i − h/ , x i + h/ Z x i + h/ x i − h/ u t dx + Z x i + h/ x i − h/ f x dx = Z x i + h/ x i − h/ s ( x ) dx, (2)and obtain du i dt + 1 h [ f ( u i +1 / ) − f ( u i − / )] = s i , (3)where u and s are the cell-averaged solution and forcing term: u i = 1 h Z x i + h/ x i − h/ u dx, s i = 1 h Z x i + h/ x i − h/ s ( x ) dx, (4)and u i − / and u i +1 / are the point-valued solutions at the left and right faces, respectively.Note that the integral form is exact at this point and no approximation has yet been made. It thus representsthe exact evolution equation for the cell-averaged solution u i . Then, a natural discretization approach is to store thecell-averaged solution at a cell as a numerical solution, compute the left and right fluxes with point-valued solutionsreconstructed from the cell-averages, and then update the cell-averaged solution in time. This is the basis of thefinite-volume method and the essence of the MUSCL approach of Van Leer for achieving second- and higher-orderaccuracy [1, 2]. Note also that it is absolutely clear already at this point from the exactness of the integral formthat the MUSCL finite-volume method is arbitrarily high-order accurate with an arbitrarily high-order of solutionreconstruction scheme. Therefore, the claim that the MUSCL scheme cannot be third-order [10] is false.However, a confusion can easily arise when one attempts to prove third-order accuracy by a truncation erroranalysis for nonlinear equations. One of the main reasons is the lack of distinction of the cell-averaged solution u i and the point-valued solution u i = u ( x i ), which differ by a second-order error (see Figure 1(b)): u i = 1 h Z x i + h/ x i − h/ u ( x ) dx = u i + 124 ( u xx ) h + O ( h ) , (5)3 i − x i − x i q i ( x ) L R x i +1 x i +2 xu Figure 2: Quadratically reconstructed point-valued solution within each control volume on a uniform grid.where u ( x ) is a smooth solution that can be expanded around x = x i as u ( x ) = u i + ( u x )( x − x i ) + 12 ( u xx )( x − x i ) + 16 ( ∂ xxx u )( x − x i ) + 124 ( ∂ xxxx u )( x − x i ) + O ( h ) . (6)Another important point that can be easily missed is that it is not the flux derivative f x but the cell-averaged fluxderivative f x that the flux difference [ f ( u i +1 / ) − f ( u i − / )] /h represents: du i dt + f x = s i , (7)where f x = 1 h Z x i + h/ x i − h/ f x dx = 1 h [ f ( u i +1 / ) − f ( u i − / )] . (8)Therefore, the truncation error must be derived for the cell-averaged flux derivative f x , not the flux derivative f x at a point. To clarify the confusion, we provide a detailed derivation of the third-order truncation error of theMUSCL scheme for a general nonlinear conservation law. On the one-dimensional grid defined in the previous section, we store the cell-averaged solution u i at a cell i asa numerical solution, and define a finite-volume scheme as the direct application of the integral form: du i dt + 1 h [ F ( u i +1 / ,L , u i +1 / ,R ) − F ( u i − / ,L , u i − / ,R )] = s i , (9)where F denotes a numerical flux computed from two point-valued solutions reconstructed at a face, e.g., u i +1 / ,L and u i +1 / ,R at the right face. In this paper, we consider an upwind numerical flux in the form: F ( u L , u R ) = 12 [ f ( u L ) + f ( u R )] − D ( u R − u L ) , (10)where D is a dissipation coefficient D = | ∂f /∂u | , and u L and u R are reconstructed solutions at a face from the leftand right cells, respectively. See Figure 2 for an example of piecewise quadratic polynomials reconstructed fromthe cell-averaged solutions. Here, we assume that the forcing term will be integrated exactly or by a sufficientlyaccurate quadrature formula. Then, the above discretization is exact if the reconstructed solutions are computedexactly; of course, they cannot be exact, and thus the solution reconstruction accuracy determines the order ofaccuracy of the finite-volume discretization. For the solution reconstruction, we consider Van Leer’s κ -reconstruction scheme, which first appeared in Ref.[3]and later was formulated specifically for a finite-volume discretization in Ref.[4]: for the right face at i + 1 / u L = u i + ( u x ) i (cid:18) h (cid:19) + 3 κ u xx ) i "(cid:18) h (cid:19) − h , (11) u R = u i +1 + ( u x ) i +1 (cid:18) − h (cid:19) + 3 κ u xx ) i +1 "(cid:18) − h (cid:19) − h , (12)where ( u x ) i = u i +1 − u i − h , ( u xx ) i = u i +1 − u i + u i − h , (13)( u x ) i +1 = u i +2 − u i h , ( u xx ) i +1 = u i +2 − u i +1 + u i h , (14)or they can be simplified as u L = 12 ( u i + u i +1 ) − − κ u i +1 − u i + u i − ) , (15) u R = 12 ( u i +1 + u i ) − − κ u i +2 − u i +1 + u i ) . (16)Note that we have denoted u i − / ,L and u i − / ,R simply as u L and u R for brevity; we consider only the right faceacross cells i and i + 1 from now on since the left face is treated in a similar manner by shifting the index. Noticethat the κ -reconstruction scheme has been generated from a local quadratic polynomial defined in each cell. Forexample, in the cell i shown in Figure 2, the quadratic polynomial q i ( x ) is given by q i ( x ) = u i + ( u x ) i ( x − x i ) + 12 ( u xx ) i (cid:20) ( x − x i ) − h (cid:21) , (17)where ( u x ) i = u i +1 − u i − h , ( u xx ) i = u i +1 − u i + u i − h , (18)which gives u L at x = x i + h/ u i :1 h Z x i + h/ x i − h/ q i ( x ) dx = u i . (19)Note also that the derivatives ( u x ) i and ( u xx ) i are point values as can be easily proved as dq i ( x ) dx (cid:12)(cid:12)(cid:12)(cid:12) x = x i = ( u x ) i , d q i ( x ) dx (cid:12)(cid:12)(cid:12)(cid:12) x = x i = ( u xx ) i . (20)As pointed out in Ref.[4], the κ -reconstruction scheme corresponds to Fromm’s scheme with κ = 0, a quadraticpoint-valued solution reconstruction (from cell averages) with κ = 1 /
3, a quadratic interpolation scheme (i.e., theQUICK interpolation scheme [32]) with κ = 1 /
2, and the central scheme with κ = 1. As we will show, third-orderaccuracy is achieved for the finite-volume scheme only with κ = 1 / Remark : Readers who are not familiar with the κ -reconstruction scheme written as in Equations (11) and (12)may find the following list of equivalent forms useful: Equations (11) and (12) are equivalent to u L = u i + 14 (cid:8) (1 − κ )∆ − i + (1 + κ )∆ + i (cid:9) , u R = u i +1 − (cid:8) (1 − κ )∆ + i +1 + (1 + κ )∆ − i +1 (cid:9) , (21)where ∆ − i = u i − u i − , ∆ + i = u i +1 − u i , (22)5hich may be a convenient form to incorporate classical limiters [33, 34], or u L = u i + κ u i +1 − u i ) + (1 − κ )( u x ) i h , u R = u i +1 − κ u i +1 − u i ) − (1 − κ )( u x ) i +1 h , (23)which has been found useful for extending κ -reconstruction scheme to unstructured grids [11, 29], or u L = κ u i + u i +1 − κ ) (cid:20) u j + 12 ( u x ) i h (cid:21) , u R = κ u i +1 + u i − κ ) (cid:20) u i +1 −
12 ( u x ) i +1 h (cid:21) . (24)which has been found useful for explaining how an unstructured-grid version of κ -reconstruction scheme can losethe linear exactness on unstructured grids [35]. Note also that the reconstruction scheme proposed in Ref.[11],which they call the β -scheme with a parameter β , is equivalent to the κ -reconstruction scheme with κ = 1 − β ;thus it is a re-arrangement of the κ -reconstruction scheme. The scheme proposed in Ref.[36] is not equivalent butclosely related to the κ -reconstruction scheme. With an additional term added, it is claimed to achieve third-orderaccuracy with κ = − / To prove third-order accuracy of the MUSCL scheme with κ = 1 /
3, we first derive the so-called modifiedequation for the MUSCL scheme by expanding the numerical solution in a Taylor series [37]. Then, we identify thetruncation error as the leading term in the difference between the modified equation and the target conservationlaw. For convenience, we will use the Taylor expansion of a point-valued solution and then take its cell-average toexpand the cell-averaged solution.In the rest of the section, we proceed as follows. First, we examine the accuracy of the solution reconstructionand show that the averaged solution at a face is exact for a cubic solution with κ = 1 /
3. Then, we prove thatthe same is true for the averaged flux term in the numerical flux (10); this is an important point for nonlinearequations. Furthermore, we will prove that the dissipation term in the numerical flux (10) is a fourth-order quantityat a face and therefore will only generate a third-order truncation error. Finally, we derive a modified equation ofthe MUSCL scheme by converting the resulting point-valued flux derivative to a cell-averaged derivative. Below,we will present as much detail as possible to leave no room for misunderstanding.
To examine the accuracy of the κ -reconstruction scheme, we consider a Taylor expansion of a smooth point-valued solution u ( x ) around a cell center x = x i : u ( x ) = u i + ( u x )( x − x i ) + 12 ( u xx )( x − x i ) + 16 ( u xxx )( x − x i ) + 124 ( u xxxx )( x − x i ) + O ( h ) , (25)where u i is a point value at x = x i . It is very important to note here that the derivatives such as ( u x ) and ( u xx )are also point values defined at x = x i , which can be easily verified by differentiating u ( x ) and evaluating the resultat x = x i . Here, we do not use the subscript i for the derivatives to distinguish them from the finite-differenceapproximations in Equations (13) and (14). It is possible and actually more convenient for our purpose to usea Taylor expansion of a cell-averaged solution with cell-averaged derivatives [19], but we do not employ such anapproach here to illustrate and emphasize the importance of recognizing the difference between the cell averageand the point value.The exact reconstruction of a smooth solution at the right face is given by u exacti +1 / = u ( x i + h/
2) = u i + 12 ( u x ) h + 18 ( u xx ) h + 148 ( u xxx ) h + 1384 ( u xxxx ) h + O ( h ) . (26)This is the target expression we wish to approximate as accurately as possible with the κ -reconstruction scheme.The κ -reconstruction scheme is expected to be exact up to the quadratic term since it is derived from a quadraticpolynomial as discussed in the previous section, which is indeed true. To see this, we express (26) in terms of thecell-average stored at the cell i : take a cell-average of Equation (25), u i = 1 h Z x i + h/ x i − h/ u ( x ) dx = u i + 124 ( u xx ) h + 11920 ( u xxxx ) h + O ( h ) , (27)6olve it for u i , and substitute it into Equation (25) to get u ( x ) = u i + ( u x )( x − x i ) + 12 ( u xx ) (cid:20) ( x − x i ) − h (cid:21) + 16 ( u xxx )( x − x i ) + O ( h ) , (28)which gives u exacti +1 / = u i + ( u x ) h u xx ) "(cid:18) h (cid:19) − h + 16 ( u xxx ) (cid:18) h (cid:19) + O ( h ) (29)= u i + 12 ( u x ) h + 112 ( u xx ) h + 148 ( u xxx ) h + O ( h ) . (30)Comparing the κ -reconstruction scheme (11) with this exact expression, we find κ = 1 / u L and u R are averaged. This fact is pointed out in Ref.[38] by quoting Ref.[39], where a third-order finite-differencescheme is proposed, which is equivalent to the MUSCL scheme with κ = 1 / Remark : The cubic exactness is a special property of a quadratically exact algorithm on a regular grid; thequadratic exactness is sufficient to design a third-order scheme on arbitrary grids. The resulting one-order highertruncation error is the reason that the truncation error order matches the discretization error (i.e., solution error) or-der on regular grids; the truncation error order is typically one-order lower on irregular grids. See Refs.[40, 41, 42, 31]and references therein for further details about accuracy on irregular grids.To prove the cubic exactness, we expand the κ -reconstruction schemes (11) and (12) by substituting the expansionsof the cell-averages: u i − = 1 h Z x i − h/ x i − h/ u ( x ) dx = u i − ( ∂ x u ) h + 1324 ( u xx ) h −
524 ( u xxx ) h + O ( h ) , (31) u i = 1 h Z x i + h/ x i − h/ u ( x ) dx = u i + 124 ( u xx ) h + 11920 ( ∂ xxxx u ) h + O ( h ) , (32) u i +1 = 1 h Z x i +3 h/ x i + h/ u ( x ) dx = u i + ( ∂ x u ) h + 1324 ( u xx ) h + 524 ( u xxx ) h + O ( h ) , (33) u i +2 = 1 h Z x i +5 h/ x i +3 h/ u ( x ) dx = u i + 2( ∂ x u ) h + 4924 ( u xx ) h + 1712 ( u xxx ) h + O ( h ) , (34)and obtain u L = u i + 12 ( u x ) h + 14 ( u xx ) (cid:18) κ + 16 (cid:19) h + 548 ( u xxxx ) h + 60 κ + 11920 ( u xxxxx ) h + O ( h ) ,u R = u i + 12 ( u x ) h + 14 ( u xx ) (cid:18) κ + 16 (cid:19) h − (cid:18) − κ (cid:19) ( u xxxx ) h − (cid:18) − κ (cid:19) ( u xxxxx ) h + O ( h ) , thus u L + u R u i + 12 ( u x ) h + 14 ( u xx ) (cid:18) κ + 16 (cid:19) h − (cid:18) − κ (cid:19) ( u xxxx ) h − (cid:18) − κ (cid:19) ( u xxxxx ) h + O ( h ) , (35)which matches the exact reconstruction (26) up to the cubic term if κ + 16 = 12 , − κ = − , (36)both of which are satisfied with κ = 1 /
3, thus resulting in u L + u R u i + 12 ( u x ) h + 18 ( u xx ) h + 148 ( u xxxx ) h − u xxxxx ) h + O ( h ) . (37)7herefore, the κ -reconstruction scheme leads to the average of the reconstructed solutions exact for a cubic functionson a uniform grid. For a linear problem with f = au , where a is a constant, this means that the flux is reconstructedexactly for a cubic flux and thus it leads to a third-order MUSCL scheme for linear equations. However, fornonlinear equations, we must prove that the average of the nonlinear fluxes evaluated with the reconstructedsolutions [ f ( u L ) + f ( u R )] / For nonlinear equations, the flux is a nonlinear function of u and thus it is not immediately clear if the averagedflux [ f ( u L ) + f ( u R )] / f ( u L ): f ( u L ) = f u i + ( u x ) i (cid:18) h (cid:19) + 3 κ u xx ) i "(cid:18) h (cid:19) − h , (38)which one might simply expand with respect to u i , f ( u L ) = f ( u i ) + O (( u L − u i ) ) . (39)However, we need a point-wise expansion of the flux around the cell center based on f ( u i ), not f ( u i ), which differsby O ( h ) and does not make sense. This is an important point because the flux derivatives arising from the fluxexpansion and to be used to represent the target equation in the modified equation should be point values, notthose evaluated by the cell-averaged solution. The distinction is not made and it seems to be one of the majorproblems in the proof presented in Ref.[5]. To derive a correct truncation error, we write f ( u L ) = f ( u i + du ) , du L = u L − u i , (40)and expand it around the point value f ( u i ) as f ( u L ) = f ( u i ) + ( f u ) du L + 12 ( f uu ) du L + 16 ( f uuu ) du L + O ( du ) , (41)and similarly for f ( u R ), f ( u R ) = f ( u i ) + ( f u ) du R + 12 ( f uu ) du R + 16 ( f uuu ) du R + O ( du ) , (42)where du R = u R − u i . Note that all the derivatives, e.g., ( f u ) and ( f uu ), are point values at the cell center x = x i .Then, we take the average to get f ( u L ) + f ( u R )2 = f ( u i ) + 12 ( f u )( u x ) h + 18 (cid:20) κ + 13 ( f u )( u xx ) + ( f uu )( u x ) (cid:21) h + 148 (cid:2) ( f uuu )( u x ) + (6 κ + 1) ( f uu )( ∂ x u )( u xx ) + (6 κ −
1) ( f u )( ∂ xxx u ) (cid:3) h + O ( h ) . (43)Comparing it with the exact flux reconstruction expressed in the Taylor series, f exactj +1 / = f ( u i ) + 12 ( f x ) h + 18 ( f xx ) h + 148 ( f xxx ) h + O ( h ) , = f ( u i ) + 12 ( f uu )( u x ) h + 18 (cid:2) ( f u )( u xx ) + ( f uu )( u x ) (cid:3) h + 148 (cid:2) ( f uuu )( u x ) + 3( f uu )( u x )( u xx ) + ( f u )( u xxx ) (cid:3) h + O ( h ) , (44)we find that the averaged flux matches the exact flux up to the cubic term if6 κ + 13 = 1 , κ + 1 = 3 , κ − , (45)8hich are all satisfied with κ = 1 /
3, thus giving f ( u L ) + f ( u R )2 = f ( u i ) + 12 ( f u )( u x ) h + 18 (cid:2) ( f u )( u xx ) + ( f uu )( u x ) (cid:3) h + 148 (cid:2) ( f uuu )( u x ) + 3( f uu )( u x )( u xx ) + ( f u )( u xxx ) (cid:3) h + O ( h ) . (46)Therefore, the average of the left and right fluxes evaluated with the solutions reconstructed with the κ -reconstructionscheme is exact for a cubic flux with κ = 1 / In this section, we show that the dissipation term is of O ( h ) and does not contribute to the second-ordertruncation error. Consider the solution jump at the right face, which can be expanded by substituting the expandedcell-averages (31)-(34): u R − u L = ( κ − (cid:20)
14 ( u xxx ) h + 18 ( u xxxxx ) h + O ( h ) (cid:21) . (47)Observe that the jump vanishes for κ = 1 and it is consistent with the fact that κ = 1 corresponds to the centralscheme, which has no dissipation. For a linear equation, the dissipation coefficient D in the upwind flux (10) is aglobal constant and the leading third-order term will cancel with the same term arising from the expansion of thejump from the left face. However, for nonlinear equations, the dissipation coefficient is generally a function of thesolution. It is typically evaluated with the averaged solution: D = D (cid:18) u L + u R (cid:19) . (48)Note that other averages can be used but all differ from the arithmetic average by O ( u R − u L ) = O ( h ), whichis negligibly small for the purpose of our analysis. Hence, it suffices to consider the arithmetic average as in theabove equation. This coefficient needs to be expanded as well in the truncation error analysis. In doing so, weassume that the dissipation coefficient is differentiable with a smoothing technique (e.g., see Ref.[43]) applied tothe absolute value of the eigenvalue in constructing the absolute Jacobian D = | ∂f /∂u | . Then, we expand thedissipation coefficient (48) by using Equation (37), and find D = D ( u i ) + 12 D u ( u x ) h + O ( h ) . (49)where D u is the derivative of D with respect to u at x = x i . Therefore, the dissipation term at the right face isexpanded as D ( u R − u L ) = (cid:20) D ( u i ) + 12 D u ( u x ) h + O ( h ) (cid:21) ( κ − (cid:20)
14 ( u xxx ) h + 18 ( u xxxxx ) h + O ( h ) (cid:21) (50)= D ( u i ) κ −
14 ( u xxx ) h + κ −
18 [ D u ( u x )( u xxx ) + ( u xxxx )] h + O ( h ) . (51)By performing a similar expansion for the dissipation term at the left face, we obtain exactly the same third-orderleading term (and the same fourth-order term with the opposite sign). Therefore, the third-order contribution fromthe dissipation term will cancel out in the flux difference in Equation (9). Then, as the flux difference is divided by h , the leading error coming from the dissipation term will be O ( h ). Hence, the truncation error of the dissipationin the MUSCL scheme is already sufficiently small for third-order accuracy. For this reason, we do not considerthe dissipation term in the analysis below. We now derive the spatial truncation error of the MUSCL scheme by substituting the Taylor expansions of thefluxes and solutions as derived in the previous sections. Substituting the expansion of the averaged flux (43) at theright face and a similar expression for that at the left face, we obtain du i dt + f x + 124 (cid:2) f uuu ( u x ) + (6 κ + 1) f uu u x u xx + (6 κ − f u u xxx (cid:3) h + O ( h ) = s i , (52)9here we have used f u u x = f x , which is true since all derivatives are point values at x = x i . There is a second-ordererror term, but the derivation is not completed yet. As mentioned earlier and repeatedly pointed out by Leonard[9, 28], the spatial operator that the finite-volume discretization is trying to approximate is not the pointwise f x but the cell average f x . To derive a correct truncation error, we replace f x by f x . Consider the cell-average of f x . f x = 1 h Z x i + h/ x i − h/ f x dx = f x + 124 f xx h + 11920 f xxxxx h + O ( h ) , (53)from which we find f x = f x − f xxx h − f xxxxx h − O ( h ) , (54)and substituting this into Equation (52), we obtaiin du i dt + f x + 124 (cid:2) f uuu ( u x ) + (6 κ + 1) f uu u x u xx + (6 κ − f u u xxx − f xxx (cid:3) h + O ( h ) = s i , (55)which becomes by f xxx = f uuu ( u x ) + 3 f uu u x u xx + f u u xxx du i dt + f x + 3 κ −
112 [ f uu u x u xx + f u u xxx ] h + O ( h ) = s i . (56)Then, the second-order term vanishes for κ = 1 / du i dt + f x + O ( h ) = s i . (57)This is the correct modified equation for the third-order MUSCL scheme. By comparing with the target integralform (7), we find that the truncation error is O ( h ). Therefore, the MUSCL scheme is third-order accurate with κ = 1 / κ = 1 /
2. However, itdoes not disprove third-order accuracy of the QUICK scheme of Leonard [32] because the QUICK scheme is notequivalent to the k = 1 / κ = 1 / Before we move on to numerical experiments, we consider a MUSCL scheme for diffusion, which is required tosolve a convection-diffusion equation: f x − νu xx = 0 , (58)where ν is a positive constant. This is an important subject but rarely discussed in the literature. To show itssignificance, we consider the following fourth-order diffusion scheme: f x − ν − u i − + 12 u i − − u i + 12 u i +1 − u i +2 h = 0 . (59)which is indeed fourth-order accurate as a finite-difference scheme, f x − νu xx + 37 u xxxxx h + O ( h ) = 0 , (60)but only second-order accurate as a finite-volume scheme approximating the integral form,1 h Z x i + h/ x i − h/ f x dx − h Z x i + h/ x i − h/ ( νu xx ) dx + ν
24 ( u xxxx ) h + O ( h ) = 0 . (61)Therefore, if the third-order MUSCL scheme is used with this fourth-order diffusion scheme, the resulting schemewill be second-order accurate unless the diffusion term is negligibly small.10o construct a third-order MUSCL scheme for the convection-diffusion equation, we must construct a fourth-order finite-volume diffusion scheme. Consider the finite-volume discretization of the integral form of the diffusionterm: − h Z x i + h/ x i − h/ ( νu xx ) dx ≈ − h h F di +1 / − F di − / i , (62)where F d is the alpha-damping diffusive flux [44, 45], F d ( u L , u R ) = −
12 [ ν ( u x ) L + ν ( u x ) R ] + να h ( u R − u L ) , (63) α is a constant damping coefficient to be determined later, and u L , ( u x ) L , u R , and ( u x ) R are reconstructed point-valued solutions and derivatives at a face from the left and right cells, respectively. The left and right reconstructedsolutions are evaluated by the κ -reconstruction scheme: Equations (15) and (16). As mentioned earlier, thesesolution values are evaluated by reconstructed quadratic polynomials, e.g., q i ( x ) given in Equation (17) for u L .Then, it would be reasonable to compute ( u x ) L by differentiating and evaluating q i ( x ):( u x ) L = dq i ( x ) dx (cid:12)(cid:12)(cid:12)(cid:12) x = x i + h/ = ( u x ) i = u i +1 − u i − + 3 κ ( u i +1 − u i + u i − )2 h , (64)which becomes for κ = 1 /
3, ( u x ) L = u i +1 − u i h . (65)Similarly, we find( u x ) R = dq i +1 ( x ) dx (cid:12)(cid:12)(cid:12)(cid:12) x = x i +1 − h/ = ( u x ) i +1 = u i +2 − u i + 3 κ ( u i +2 − u i +1 + u i )2 h . (66)which again becomes for κ = 1 /
3, ( u x ) R = u i +1 − u i h . (67)Therefore, the derivative of the quadratic reconstruction is continuous across the face. Then, the diffusion schemeis given by − νh h F di +1 / − F di − / i = − ν C i − u i − + C i − u i − + C i u i + C i +1 u i +1 + C i +2 u i +2 h , (68)where C i − = C i +2 = 18 [ α ( κ − − κ − , C i − = C i +1 = 18 [24 κ − α ( κ − , C i = 18 [6 α ( κ − − κ + 1)] , (69)which can be expanded as − νh h F di +1 / − F di − / i = − νu xx −
18 [ α ( κ − − κ − νu xxxx ) h + O ( h ) , (70)thus leading to − νh h F di +1 / − F di − / i = 1 h Z x i + h/ x i − h/ ( νu xx ) dx − (cid:20) α ( κ − − κ − − (cid:21) ( νu xxxx ) h + O ( h ) . (71)Therefore, to achieve fourth-order accuracy, we must set α ( κ − − κ − −
13 = 0 , (72)or α = 2(4 − κ )3(1 − κ ) , (73)11 (a) Initial solution (cell averages). (b) Final solution (cell averages). Figure 3: Initial and final solutions for the unsteady test case.giving for κ = 1 / α = 1 (74)Note that it gives α = 8 / κ = 0, which is consistent with the value derived for fourth-order accuracy with κ = 0 in Refs.[44, 45]. Note also that the diffusion scheme (59) can be obtained with κ = 1 / α = 3 /
2, whichdoes not satisfy the condition (73). It is interesting to note that the fourth-order diffusion scheme compatible withthe MUSCL scheme is unique as the substitution of Equation (73) into the scheme (69) gives − h [ F di +1 / − F di − / ] = − ν − u i − + 16 u i − − u i + 16 u i +1 − u i +2 h . (75)This is the diffusion scheme required to achieve third-order accuracy for the convection-diffusion equation (58). We consider a time-dependent problem for Burgers’s equation: u t + f x = 0 , (76)where f = u / x ∈ [0 ,
1] with the initial solution, u ( x ) = sin(2 πx ) , (77)The domain is taken to be periodic (i.e., there is no boundary in this problem), and the solution is computed atthe final time t = t f = 0 . D is computed by the arithmetic average of the left and right values: D = ( | u L | + | u R | )2, which corresponds to theRoe linearization [46] with a smoothing technique of Harten and Hyman [43] applied to D . The time integrationis performed by the three-stage SSP Runge-Kutta scheme [47] for the total of 840 time steps with a constant timestep ∆ t = 0 . κ = 0, 1 /
2, and 1 / i , u i = 1 h Z x i + h/ x i − h/ sin(2 πx ) dx = 12 πh [cos ( π ( h − x i )) − cos ( π ( h + 2 x i ))] = 1 πh sin(2 πx i ) sin ( πh ) , (78)12nd then begin the time integration. The initial cell-averaged solution is plotted in Figure 3(a). Note that if theinitial solution is given by the point value u i = sin(2 πx i ), then a second-order error will be introduced immediately,even before we begin the time integration. Later, we will demonstrate this numerically.To emphasize the difference between the cell average and the point value, we measure the discretization errorin three different ways, locally at a cell i , E p ( x i ) = | u i − u exact ( x i , t f ) | , E c ( x i ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u i − h Z x i + h/ x i − h/ u exact ( x, t f ) dx (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , ˆ E p ( x i ) = | ˆ u i − u exact ( x i , t f ) | , (79)where u exact ( x, t f ) is the exact point-valued solution at the final time, u exact ( x, t f ) = sin(2 π ( x − u exact t )) , (80)which is defined implicitly and thus solved iteratively by the fix-point iteration (see Figure 3(b) for the finalsolution plotted on the coarsest grid), and ˆ u i is a point-valued solution at the cell center x = x i recovered from thecell-averaged solution u i as ˆ u i = u i − (cid:18) u i − − u i + u i +1 h (cid:19) h , (81)which is an accurate representation of the point value at x = x i up to fourth-order and thus sufficiently accurateto verify third-order accuracy of the numerical solution. Finally, we define the discretization error by the L norm: L ( E p ) = n X i =1 E p ( x i ) n , L ( E c ) = n X i =1 E c ( x i ) n , L ( ˆ E p ) = n X i =1 ˆ E p ( x i ) n , (82)where n is the number of cells in a grid. -3.5 -3 -2.5 -2 -1.5 -1-10-9-8-7-6-5-4-3-2 (a) Error convergence for L ( E c ). -3.5 -3 -2.5 -2 -1.5 -1-10-9-8-7-6-5-4-3-2 (b) Error convergence for L ( ˆ E p ). -3.5 -3 -2.5 -2 -1.5 -1-10-9-8-7-6-5-4-3-2 (c) Error convergence for L ( E p ). Figure 4: Error convergence results for the unsteady test case: cell-averaged initial solutions.Error convergence results are shown in Figure 4 for the case of cell-averaged initial solutions. As can be seenclearly in Figure 4(a), the numerical solution is third-order accurate as a cell-averaged solution with κ = 1 / κ = 1 /
2, which corresponds to the QUICK interpolation scheme, leads to second-orderaccuracy although slightly more accurate than Fromm’s scheme ( κ = 0). Figure 4(b) shows the result for thepoint-valued solution recovered from the cell-averaged solution as in Equation (81). Clearly, the recovered point-value solution is third-order accurate for κ = 1 /
3, but remains second-order accurate for the other choices. Finally,Figure 4(c) shows the error convergence for the pointwise error norm L ( E p ). As expected, the numerical solutionis second-order accurate as a point value for any κ .To illustrate the importance of using cell-averaged initial solutions, we performed the same numerical experi-ments with point-valued initial solutions: u i = sin(2 πx i ) at t = 0 for all cells, i = 1 , , , · · · , n . Error convergenceresults are shown in Figure 5. As expected, the error is second-order for all error norms, indicating that the erroris dominated by the initial second-order error committed in the initial solution. If we wish to avoid computing the13 (a) Error convergence for L ( E c ). -3.5 -3 -2.5 -2 -1.5 -1-10-9-8-7-6-5-4-3-2 (b) Error convergence for L ( ˆ E p ). -3.5 -3 -2.5 -2 -1.5 -1-10-9-8-7-6-5-4-3-2 (c) Error convergence for L ( E p ). Figure 5: Error convergence results for the unsteady test case: point-valued initial solution -3.5 -3 -2.5 -2 -1.5 -1-10-9-8-7-6-5-4-3-2 (a) Error convergence for L ( E c ). -3.5 -3 -2.5 -2 -1.5 -1-10-9-8-7-6-5-4-3-2 (b) Error convergence for L ( ˆ E p ). -3.5 -3 -2.5 -2 -1.5 -1-10-9-8-7-6-5-4-3-2 (c) Error convergence for L ( E p ). Figure 6: Error convergence results for the unsteady test case: corrected point-valued initial solution.cell-averaged initial solution, we may compute a cell-averaged initial solution from a point-valued initial solutionwith sufficient accuracy by u i = u i + 124 (cid:18) u i − − u i + u i +1 h (cid:19) h . (83)In this way, the initial solution does not been to be integrated over each cell. Results obtained with this correctionare shown in Figure 6. As can be seen in Figures 6(a) and 6(b), third-order accuracy is achieved for the cell-averagedsolution and the recovered point-value solution. To demonstrate third-order accuracy for a convection equation with a forcing term, we consider a steady problemfor Burgers’s equation in x ∈ [0 , f x = s ( x ) , (84)where f = u /
2, with the forcing term, s ( x ) = 2 sin(2 x ) cos(2 x ) , (85)so that the exact solution is given by u ( x ) = sin(2 x ) . (86)14or this problem, we drop the time-derivative term and define the residual at a cell i as Res i = 1 h [ F ci +1 / − F ci − / ] − s i , (87)where the forcing term is exactly cell-averaged over each cell: s i = 1 h Z x i + h/ x i − h/ s ( x ) dx = 12 h (cid:2) cos ( h − x i ) − cos ( h + 2 x i ) (cid:3) . (88)Then, we solve the system of nonlinear residual equations for the numerical solution: u , u , · · · , u n − , u n − . Notethat we will provide the exact cell-averaged solutions at the left two cells i = 1 and 2, and at the right cells i = n − n , in order to exclude boundary effects, which are beyond the scope of the present study. An implicit solverbased on the exact Jacobian of the first-order scheme is used to solve the residual equations. See Ref.[48], forexample, for further details of the implicit solver for a one-dimensional finite-volume scheme. To verify the orderof accuracy, we solve the steady problem over a series of grids with 15, 31, 63, 127 cells. As in the unsteady case,we consider two discretization error norms: L ( E p ) = 1 n − n − X i =3 | u i − u exacti | , L ( E c ) = 1 n − n − X i =3 | u i − u exacti | , (89)where u exacti is the exact cell-averaged solution: u exacti = 1 h Z x i + h/ x i − h/ sin(2 x ) dx = 12 h [cos ( h − x i ) − cos ( h + 2 x i )] = 12 h sin(2 x i ) sin ( h ) . (90)We also verify the order of truncation error numerically by substituting the exact solution into the residual (87),and taking the L norm over the cells. We consider both the point-valued and cell-averaged solutions, and definethe following two truncation error norms: L ( T p ) = 1 n − n − X i =3 Res i ( { u exacti } ) , L ( T c ) = 1 n − n − X i =3 Res i ( { u exacti } ) , (91)where Res i ( { u exacti } ) is the residual with the point-valued exact solution substituted, and Res i ( { u exacti } ) is theresidual with the cell-averaged exact solution substituted,Error convergence results are shown in Figure 7. We first consider our target case, where the errors are measuredwith the exact cell-averaged solution. Figure 7(a) shows the truncation error convergence for L ( T c ). As expected,it is third-order only with κ = 1 /
3. The discretization error is also third-order in both the cell-averaged solutionand the recovered point-value solution. These results confirm our truncation error analysis.On the other hand, if we compute the truncation error with the point-valued exact solution, we see third-orderaccuracy achieved with κ = 1 / κ = 1 / κ = 1 / κ = 1 / κ = 1 / To demonstrate the importance of the diffusion scheme, we consider a steady problem for the viscous Burgersequation with a forcing term: f x = νu xx + s ( x ) , (92)15 (a) L ( T c ). -2.2 -2 -1.8 -1.6 -1.4 -1.2 -1-7-6.5-6-5.5-5-4.5-4-3.5-3-2.5 (b) L ( E c ). -2.2 -2 -1.8 -1.6 -1.4 -1.2 -1-7-6.5-6-5.5-5-4.5-4-3.5-3-2.5 (c) L ( ˆ E p ). -2.2 -2 -1.8 -1.6 -1.4 -1.2 -1-7-6-5-4-3-2 (d) L ( T p ). -2.2 -2 -1.8 -1.6 -1.4 -1.2 -1-7-6.5-6-5.5-5-4.5-4-3.5-3-2.5 (e) L ( E p ). Figure 7: Truncation and discretization error convergence results for the case of the steady Burgers equation.where f = u / s ( x ) = 2 sin(2 x ) cos(2 x ) + 4 ν sin(2 x ) , (93)so that the exact solution is given by u ( x ) = sin(2 x ) . (94)Again, we integrate the forcing term is integrated exactly over each cell: s i = 1 h Z x i + h/ x i − h/ s ( x ) dx = 12 h (cid:2) cos ( h − x i ) − cos ( h + 2 x i ) (cid:3) − νh [cos( h − x i ) − cos( h + 2 x i )] . (95)In particular, we consider the case ν = 1, for which the convective and diffusion terms are equally important. Asbefore, we solve the steady problem in x ∈ [0 ,
1] by the implicit solver with κ = 0, 1 /
2, and 1 /
3, for a series of gridswith 15, 31, 63, 127 cells. The damping coefficient α is determined by Equation (73) for a given κ and the same u L and u R (i.e., the same κ ) are used for both convective and diffusive fluxes. Note that u L and u R are used tocompute the damping term in the diffusive flux. The analysis predicts that third-order accuracy is obtained onlyfor κ = 1 /
3. To verify the analysis, we computed the truncation errors as well, and the discretization errors forboth the cell-averaged and point-value exact solutions as described in the previous section.Figure 8 shows the results in terms of the point-value solution. As shown in Figure 8(a), the truncation error L ( T p ) is second-order for all choices of κ . Figure 8(b) confirms that the discretization errors are also second-orderfor all choices of κ . Therefore, unlike the previous case, third-order accuracy is not achieved with κ = 1 /
2. This is16 (a) L ( T p ). -2.2 -2 -1.8 -1.6 -1.4 -1.2 -1-8-7-6-5-4-3-2 (b) L ( E p ). -2.2 -2 -1.8 -1.6 -1.4 -1.2 -1-7-6-5-4-3-2-1 (c) L ( T c ). -2.2 -2 -1.8 -1.6 -1.4 -1.2 -1-8-7-6-5-4-3-2 (d) L ( E c ). -2.2 -2 -1.8 -1.6 -1.4 -1.2 -1-8-7-6-5-4-3-2 (e) L ( ˆ E p ). Figure 8: Truncation and discretization error convergence results for the case of the steady viscous Burgers equation.The fourth-order finite-volume diffusion scheme with α = − κ )3(1 − κ ) . -2.2 -2 -1.8 -1.6 -1.4 -1.2 -1-7-6-5-4-3-2-1 (a) L ( T c ). -2.2 -2 -1.8 -1.6 -1.4 -1.2 -1-8-7-6-5-4-3-2 (b) L ( E c ). -2.2 -2 -1.8 -1.6 -1.4 -1.2 -1-8-7-6-5-4-3-2 (c) L ( ˆ E p ). Figure 9: Truncation and discretization error convergence results for the case of the steady viscous Burgers equation.The fourth-order finite-difference diffusion scheme (59) implemented as the alpha-damping scheme with α = 3 / O ( h ) error and loses fourth-order accuracy. Then, the resulting convection-diffusion scheme is second-order accurate at best unless the diffusion term is negligibly small compared with the convective term. Therefore,the choice of the diffusion scheme is critically important to achieve third-order accuracy.On the other hand, the scheme is third-order with the cell-averaged solution. Figure 8(c) shows the truncationerror convergence for L ( T c ). As expected, the scheme achieves third-order accuracy with κ = 1 /
3. Also, third-orderaccuracy is confirmed in the discretization errors. See Figures 8(d) and 8(e).Finally, we performed the computations with α = 3 /
2, which corresponds to the incompatible diffusion scheme(59) when κ = 1 /
3. Results are shown in Figure 9. Clearly, as predicted in Section 5, third-order accuracy isnever achieved. Therefore, it is very important to choose the right diffusion scheme in order to achieve third-orderaccuracy: the diffusion scheme must be a fourth-order finite-volume scheme with cell-averaged solution.
We have provided a detailed truncation error analysis and proved that the MUSCL scheme is third-orderaccurate with κ = 1 / κ -reconstruction scheme is exact for a cubic function on a uniform grid, when averaged across a face (i.e.,averaging the reconstructed solutions at the left and right of the face), and the same is true for the average of thefluxes evaluated with the reconstructed solutions. The importance of the diffusion scheme is also discussed: third-order accuracy is lost with an incompatible four-order diffusion scheme. Third-order accuracy has been verified bynumerical experiments for both unsteady and steady problems with Burgers’ equation. For all the cases, it hasbeen demonstrated that the MUSCL scheme gives third-order accurate cell-averaged solutions as well as third-orderaccurate point values at cell centers if the point values are accurately recovered from the third-order accurate cellaverages. In a subsequent paper, we will discuss third-order accuracy of the QUICK scheme and clarify the reasonthat third-order accuracy was observed in the point value solution with κ = 1 / Acknowledgments
The author is grateful to Emeritus Professor Bram van Leer for illuminating discussions and helpful suggestions.The author also would like to thank Jeffery A. White (NASA Langley Research Center) for valuable comments anddiscussions. The author gratefully acknowledges support from Software CRADLE, part of Hexagon, the U.S. ArmyResearch Office under the contract/grant number W911NF-19-1-0429 with Dr. Matthew Munson as the programmanager, and the Hypersonic Technology Project, through the Hypersonic Airbreathing Propulsion Branch of theNASA Langley Research Center, under Contract No. 80LARC17C0004.
References [1] van Leer, B., “Towards the Ultimate Conservative Difference Scheme. IV. A New Approach to NumericalConvection,”
J. Comput. Phys. , Vol. 23, 1977, pp. 276–299.[2] van Leer, B., “Towards the Ultimate Conservative Difference Scheme. V. A Second-Order Sequel to Godunov’sMethod,”
J. Comput. Phys. , Vol. 32, 1979, pp. 101–136.[3] van Leer, B., “Towards the Ultimate Conservative Difference Scheme. III. Upstream-centered Finite DifferenceSchemes for Ideal Compressible Flow,”
J. Comput. Phys. , Vol. 23, 1977, pp. 263–275.[4] van Leer, B., “Upwind-Difference Methods for Aerodynamic Problems Governed by the Euler Equations,”
Fluid Mechanics, Lectures in Applied Mathematics , Vol. 22, 1985, pp. 327–335.[5] Liu, M., “Computational Study of Convective-Diffusive Mixing in a Microchannel Mixer,”
Chemical Engineer-ing Science , Vol. 66, 2011.[6] Kresta, S. M., III, A. W. E., Dickey, D. S., and Atiemo-Obeng, V. A., editors,
Advances in Industrial Mixing:A Companion to the Handbook of Industrial Mixing , John Wiley & Sons, 2015.187] Song, B. and Amano, R. S., “A High-Order Bounded Discretization Scheme,”
Computational Fluid Dynamicsand Heat Transfer: Emerging Topics , edited by R. S. Amano and B. Sud´en, WIT Press, 2011, pp. 3–5.[8] Zhang, D., Jiang, C., Liang, D., and Cheng, L., “A Review on TVD Schemes and a Refined Flux-Limiter forSteady-State Calculations,”
J. Comput. Phys. , Vol. 302, 2015.[9] Leonard, B. P. and Mokhtari, S., “Beyond First-Order Upwinding: The Ultra-Sharp Alternative for Non-Oscillatory Steady-State Simulation opf Convection,”
Int. J. Numer. Meth. Fluids , Vol. 30, 1990, pp. 729–766.[10] Wu, H., Wand, L., and Sun, G., “Non-Existence of Third Order MUSCL Schemes Unified Construction ofENO Schemes and a New Discontinuity Sharpenning Technique - Stiff Source Term Approach,”
ComputationalFluid Dynamics Review 1998 , 1998, pp. 300–317.[11] Debiez, C. and Dervieux, A., “Mixed-element-volume MUSCL methods with weak viscosity for steady andunsteady flow calculations,”
Comput. Fluids , Vol. 29, 2000, pp. 89–118.[12] Abalakin, I. and A. Dervieux, T. A. K., “A vertex centered high order MUSCL scheme applying to linearisedEuler acoustics,” INRIA Report, N 4459, April 2002.[13] Camarri, S., Salvetti, M. V., Koobus, B., and Dervieux, A., “A Low-Diffusion MUSCL Scheme for LES onUnstructured Grids,”
Comput. Fluids , Vol. 33, 2004, pp. 1101–1129.[14] Abalakin, I., Bakhvalov, P., and Kozubskaya, T., “Edge-based reconstruction schemes for unstructured tetra-hedral meshes,”
Int. J. Numer. Meth. Fluids , Vol. 81, 2015, pp. 331–356.[15] Abalakin, I., Bakhvalov, P., and Kozubskaya, T., “Edge-Based Reconstruction Schemes for Prediction of NearField Flow Region in Complex Aeroacoustics Problems,”
Int. J. Aeroacoust. , Vol. 13, 2014, pp. 207–233.[16] Shu, C.-W. and Osher, S. J., “Efficient Implementation of Essentially Non-Oscillatory Shock-CapturingSchemes, II,”
J. Comput. Phys. , Vol. 83, 1989, pp. 32–78.[17] Buchm¨uller, P. and Helzel, C., “Improved Accuracy of High-Order WENO Finite Volume Methods on CartesianGrids,”
J. Sci. Comput. , Vol. 61, 2014, pp. 343–368.[18] Buchm¨uller, P., Drehe, J., and Helzel, C., “Finite Volume WENO Methods for Hyperbolic Conservation Lawson Cartesian Grids with Adaptive Mesh Refinement,”
Applied Mathematics and Computation , Vol. 272, 2016,pp. 460–478.[19] Tamaki, Y. and Imamura, T., “Efficient Dimension-by-Dimension Higher Order Finite-Volume Methods for aCartesian grid with Cell-Based Refinement,”
Comput. Fluids , Vol. 144, 2017, pp. 74–85.[20] Barth, T. J. and Frederickson, P. O., “Higher Order Solution of the Euler Equations on Unstructured Gridsusing Quadratic Reconstruction,” AIAA Paper 90-0013, 1990.[21] Luo, H., Luo, L., Nourgaliev, R., Mousseau, V. A., and Dinhb, N., “A reconstructed discontinuous Galerkinmethod for the compressibleNavier–Stokes equations on arbitrary grids,”
J. Comput. Phys. , Vol. 229, 2010,pp. 6961–6978.[22] Jalali, A. and Ollivier-Gooch, C., “Higher-order unstructured finite volume RANS solution of turbulent com-pressible flows,”
Comput. Fluids , Vol. 143, 2017, pp. 32–47.[23] Koren, B., “A Robust Upwind Discretization Method for Advection,” Report NM-R9308, April 1993.[24] Margolin, L. G. and Rider, W. J., “Numerical Regularization: The Numerical Analysis of Implicit SubgridModels,”
Implicit Large Eddy Simulations: Computing Turbulent Fluid Dynamics , edited by F. F. Grinstein,L. G. Margolin, and W. J. Rider, chap. 5, Cambridge University Press, Rhode-Saint-Genese, Belgium, 2006.[25] Johnson, R. W. and Mackinon, R. J., “Equivalent Versions of the QUICK Scheme for Finite-Difference andFinite-Volume Numerical Methods,”
Communications in Applied Numerical Methods , Vol. 8, 1992, pp. 841–847.[26] Chen, Y. and Falconer, R. A., “Modified forms of the third-order convection, second-order diffusion schemefor the advection-diffusion equation,”
Advances in Water Resources , Vol. 17, 1994, pp. 147–170.1927] Leonard, B. P., “Comparison of Truncation Error of Finite-Difference and Finite-Volume Formulations ofConvection Terms,”
Appl. Math. Modell. , Vol. 18, 1994, pp. 46–50.[28] Leonard, B. P., “Order of Accuracy of QUICK and Related Convection-Diffusion Schemes,”
Appl. Math.Modell. , Vol. 19, 1995, pp. 644–653.[29] Burg, C. O. E., “Higher Order Variable Extrapolation for Unstructured Finite Volume RANS Flow Solvers,”AIAA Paper 2005-4999, 2005.[30] Katz, A. and Work, D., “High-Order Flux Correction/Finite Difference Schemes for Strand Grids,”
J. Comput.Phys. , Vol. 282, 2015, pp. 360–380.[31] Nishikawa, H. and Liu, Y., “Accuracy-Preserving Source Term Quadrature for Third-Order Edge-Based Dis-cretization,”
J. Comput. Phys. , Vol. 344, 2017, pp. 595–622.[32] Leonard, B. P., “A Stable and Accurate Convective Modelling Procedure based on Quadratic Upstream Inter-polation,”
Computer Methods in Applied Mechanics and Engineering , Vol. 19, 1979, pp. 59–98.[33] Lynn, J. F.,
Multigrid Solution of the Euler Equations with local preconditioning , Ph.D. thesis, University ofMichigan, Ann Arbor, Michigan, 1995.[34] Hirsch, C.,
Numerical Computation of Internal and External Flows , Vol. 2, A Wiley - Interscience Publications,1990.[35] Nishikawa, H., “On the Loss and Recovery of Second-Order Accuracy with U-MUSCL,”
J. Comput. Phys. ,Vol. 417, 2020, pp. 109600.[36] Yang, H. Q. and Harris, R. E., “Development of Vertex-Centered High-Order Schemes and Implementation inFUN3D,”
AIAA J. , Vol. 54, 2016, pp. 3742–3760.[37] Hirt, C. W., “Heuristic Stability Theory for Finite-Difference Equations,”
J. Comput. Phys. , Vol. 2, 1968,pp. 339–355.[38] Waterson, N. P. and Deconinck, H., “Design principles for bounded higher-order convection schemes - a unifiedapproach,”
J. Comput. Phys. , Vol. 224, 2007, pp. 182–207.[39] Agarwal, R. K., “A Third-Order-Accurate Upwind Scheme for Navier-Stokes Solutions at High ReynoldsNumbers,” Aiaa paper 81-0112, 1981.[40] Diskin, B. and Thomas, J. L., “Accuracy Analysis for Mixed-Element Finite-Volume Discretization Schemes,”
NIA Report No. 2007-08 , 2007.[41] Diskin, B. and Thomas, J. L., “Notes on Accuracy of Finite-Volume Discretization Schemes on IrregularGrids,”
Appl. Numer. Math. , Vol. 60, 2010, pp. 224–226.[42] Katz, A. and Sankaran, V., “Mesh Quality Effects on the Accuracy of CFD Solutions on Unstructured Meshes,”
J. Comput. Phys. , Vol. 230, 2011, pp. 7670–7686.[43] Harten, A. and Hyman, J. M., “Self-Adjusting Grid Methods for One-Dimensional Hyperbolic ConservationLaws,”
J. Comput. Phys. , Vol. 50, 1983, pp. 235–269.[44] Nishikawa, H., “Beyond Interface Gradient: A General Principle for Constructing Diffusion Schemes,”
Proc.of 40th AIAA Fluid Dynamics Conference and Exhibit , AIAA Paper 2010-5093, Chicago, 2010.[45] Nishikawa, H., “Robust and Accurate Viscous Discretization via Upwind Scheme - I: Basic Principle,”
Comput.Fluids , Vol. 49, No. 1, October 2011, pp. 62–86.[46] Roe, P. L., “Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes,”
J. Comput. Phys. ,Vol. 43, 1981, pp. 357–372.[47] Gottlieb, S., Shu, C.-W., and Tadmor, E., “Strong Stability-Preserving High-Order Time Discretization Meth-ods,”
SIAM Rev. , Vol. 43, No. 1, 2001, pp. 89–112.[48] Nishikawa, H. and Liu, Y., “Hyperbolic Advection-Diffusion Schemes for High-Reynolds-Number Boundary-Layer Problems,”