A Novel Finite Difference Method for Euler Equations in 2D Unstructured Meshes
AA Novel Finite Difference Method for Euler Equationsin 2D Unstructured Meshes
Meiyuan Zhen a , Kun Qu a, ∗ , Jinsheng Cai a a School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, Shaanxi,China.
Abstract
Finite difference method was extended to unstructured meshes to solve Eulerequations. The spatial discretization is made of two steps. First, numericalfluxes are computed at the middle point of each edge with high order accuracy.In this step, a one-dimensional curvilinear stencil is assembled for each edge toperform one-dimensional fast non-uniform WENO interpolation derived in thispaper, which is much easier and faster than multi-dimensional interpolation.The second step is to compute the divergence of fluxes at each vertex fromthe fluxes at nearby edges and vertices by least square approximation of multi-dimensional polynomials. The order of the WENO interpolation in the first stepand the degree of the polynomial in the second step determined the order ofaccuracy of the spatial scheme. After that, explicit RungeKutta time discretescheme is used to update conservative variables. Several canonical numericalcases were solved to test the accuracy, performance and the capability of shockcapturing of the developed method.
Keywords: finite difference method, high order method, unstructured mesh,non-uniform WENO, least square
1. Introduction
In decades, the research of high order methods is one of the most impor-tant work in CFD, with a great number of high order methods being pro- ∗ Corresponding author
Preprint submitted to Elsevier February 26, 2021 a r X i v : . [ phy s i c s . c o m p - ph ] F e b osed. Among them, finite difference method (FDM) [11, 7, 2, 9], discontinuousGalerkin (DG) [15, 4, 5], flux reconstruction (FR) [10, 21], spectral differencemethod (SDM) [14, 20], as well as high order finite volume method (FVM),[1, 8, 19] are studied most intensively. Today, high order finite difference schemesare so matured that they are able to resolve shear layers, vortices and captureshock waves and material interfaces in rectangular or curve-linear grids [22, 18].Besides, the performance of FDM is very high since only one-dimensional opera-tions are necessary. High order of accuracy in FDM can be easily implemented.Fifth or higher order schemes are widely used in the simulation with FDM.Thus high order FDM are adopted with DNS/LES in a wide range of researchessuch as flow instability and multi-phase flows. Although highly developed, theshortcoming of FDM is obvious that generating structural grids is very difficultfor complex geometries. This hinders its applications in engineering problems.Thus, high order methods on unstructured meshes, such as DG, FR and SD, at-tract more scientists. Different from FDM, these methods rely on elements withmultiple degrees of freedom. In each element, polynomials serve as shape func-tions, which naturally leads to the difficulty of capturing strong shock waves.Although some measures such as adding artificial viscosity and slope limitersor combined with FVM are proposed to overcome this problem, which increasethe complexity of the method and therefore still needs to be improved.At the same time, some researchers tried to develop generalized finite differ-ence method (GDM) in the context of meshless methods [12, 13]. In meshlessmethods, on a virtual edge between two closed points, numerical fluxes are cal-culated with a Riemann solver and utilized to discretize the divergence of theflux on each node. The weights of computing divergence are determined fromsome kind of local kernel function based on distance. This algorithm can cap-ture shock waves since it adopted Riemann solvers. Usually, the closest pointsare used and only the 2nd order accuracy. In the work of Li and Ren [13], morepoints are introduced into each stencil to achieve higher order accuracy.In this paper, the authors developed a FDM on unstructured meshes. Inthis method, similar to Ren’s work, the divergence of flux in the governing2quations describing a conservative law is calculated on each vertex directlyfrom the nearby vertices and edges by means of least square, without any kindof real or virtual control volume and integrals on faces or volumes. The twodifferences are: 1) a stencil is constructed based the topology of the mesh andable to adapt to the local non-isentropic characteristics. It also includes moreedges thus making the stencil more compact; 2) the numerical flux on each edgeis computed from a one dimensional curve stencil but not multi-dimensionalstencil, making our method more efficient. To the authors’ knowledge, it is thefirst attempt to extract one-dimensional interpolation stencil in the context ofunstructured meshes.The rest of this paper is organized as follows. Section 2 describes several keycomponents of the developed method, including flux divergence computing andnumerical flux computing on one-dimensional stencils. Section 3 gives severalnumerical test cases to verify the performance of the method. Conclusions andpossible future work are presented in section 4.
2. Numerical Method
The two-dimensional Euler equation describes inviscid compressible flow,which is given by ∂ w ∂t + ∂ f ( w ) ∂x + ∂ g ( w ) ∂y = 0 (1)where w = ρρuρvρE f ( w ) = ρuρuu + pρuv ( ρE + p ) u g ( w ) = ρvρvuρvv + p ( ρE + p ) v and ρ , u , v , p , E are density, velocity in x , y directions, pressure and totalenergy respectively, where E = e + u + v and p = ( γ − ρe . γ = 1 . f and g are used to denote a component inthe flux vector f and g respectively. 3 .2. Discretization of Flux Divergence From the governing equations, it can be seen that the task of spatial dis-cretization is to compute ∂ x f + ∂ y g which is the divergence of fluxes.For structured grids, the flux divergence at a node is calculated along 2(in 2Dcases) or 3(in 3D cases) grid lines with some FD schemes. For instance, the firstorder derivative of the flux computed in WCNS (Weighted Compact NonlinearScheme)[6] is approximated with high order central difference from the numericalfluxes at nearby middle points or both middle points and nodes. Artificialviscosity is not used explicitly as numerical dissipation is already introducedwhen computing the numerical flux at middle points, making the flux functionsmooth enough to perform central FDM.When it comes to unstructured meshes used herein, similar to the way thatWCNS used, flux divergence is computed based on the numerical flux on eachedge and analytical flux vector on each vertex, which is described in the followingsubsections. Herein, the two flux functions f and g around a vertex O are approximatedby two dimensional polynomials f ( x, y ) = f + a x + a y + a x + a xy + a y + · · · = f + p T ( x, y ) a (2) g ( x, y ) = g + b x + b y + b x + b xy + b y + · · · = g + p T ( x, y ) b (3)where f and g are the flux components at the vertex O which is regarded asthe origin point (0 , p ( x, y ) = (cid:2) x, y, x , xy, y , x , x y, xy , y , · · · (cid:3) T (4)contains M necessary terms of two-dimensional polynomials. Eq.(2) and Eq.(3)should be satisfied in a domain near the vertex O . For a vertex j nearby, p T ( x j , y j ) a = f j − f (5) p T ( x j , y j ) b = g j − g (6)4or an edge k , only the projection of the flux vector [ f, g ] T at the middle-edgepoint along the direction of the edge is provided. Thus, n x,k p T ( x k , y k ) a + n y,k p T ( x k , y k ) b = r k − ( n x,k f + n y,k g ) (7)where r k is the flux projection at the middle point of edge k and [ n x,k , n y,k ] T represents the unit directional vector of edge k . Picking enough vertices andedges nearby, an overdetermined linear system can be assembled as p T ( x j =1 , y j =1 ) OO p T ( x j =1 , y j =1 ) p T ( x j =2 , y j =2 ) OO p T ( x j =2 , y j =2 )... ... p T ( x j = J , y j = J ) OO p T ( x j = J , y j = J ) n x,k =1 p T ( x k =1 , y k =1 ) n y,k =1 p T ( x k =1 , y k =1 ) n x,k =2 p T ( x k =2 , y k =2 ) n y,k =2 p T ( x k =2 , y k =2 )... ... n x,k = K p T ( x k = K , y k = K ) n y,k = K p T ( x k = K , y k = K ) ab = R (8)where R = f j =1 − f g j =1 − g f j =2 − f g j =2 − g ... f j = J − f g j = J − g r k =1 − ( n x,k =1 f + n y,k =1 g ) r k =2 − ( n x,k =2 f + n y,k =2 g )... r k = K − ( n x,k = K f + n y,k = K g ) (9)5he linear system Eq.(8) can be simply put as A ab = R (10)which can be solved with normal matrix method or SVD, with later used herein, ab = A † R (11)where A † is the pesudo-inverse of A . It is noticed that by utilizing the polyno-mial approximation in Eq.(2) and Eq.(3), the flux divergence at Vertex O canexpressed as weighted sum of R components, ∂f∂x (cid:12)(cid:12)(cid:12)(cid:12) x =0 ,y =0 + ∂g∂y (cid:12)(cid:12)(cid:12)(cid:12) x =0 ,y =0 = a + b = ( r † a + r † b ) R (12)where r † a and r † b are the row vectors in A † corresponding to a and b respec-tively. The vector r † a + r † b is only related to the geometry of mesh and can becomputed only once at the initial time for static meshes.By building stencils with enough vertices and edges around a vertex, highorder accuracy of the flux divergence can be achieved easily. As described above, f and g are coupled in single linear system. In orderto obtain the two polynomials of p degree, p ( p + 3) conditions are required for2D problems. Each vertex gives two conditions, while each edge provides onecondition. In an unstructured mesh, for a vertex, it is easy to obtain all thedirect neighbor vertices around it to construct a one-level stencil (Fig. 1(a) and1(b)) which contains all the direct neighbor vertices and the edges (the thicklines in the figures) between them. Sometimes, artificial edges can be createdand included in the stencil if there are not enough edges in a one-level stencil.In Fig. 1(b) and 1(c), the dashed lines are artificial edges.A one-level stencil in a triangular mesh contains N neighbor vertices and 2 N edges where N = 5 or 6, providing 4 N conditions. Thus, the accuracy order of6 a) A one-level stencil ina triangular mesh. (b) A one-level stencilin a hybrid mesh.(c) A one-level stencil withonly 4 direct neighbor ver-tices. (d) A two-level stencil in a tri-angle mesh.Figure 1: The stencils for divergence computing at a vertex (the red circle). The blue circlesare the first level neighbor vertices, while the green circles are the second level neighborvertices. std::set in C++ programminglanguage. It contains 18 vertices and 42 edges, providing 78 conditions whichshould be enough for 2D polynomials higher than 6 degree. In order to test the order of accuracy of the flux divergence computed by theleast square method, a sequence of unstructured meshes similar to Fig. 2 withdifferent mesh sizes are used. Two-level divergence stencils similar to Fig. 1(d)are utilized to support polynomial approximation of five degree, ie, fifth orderaccuracy of flux divergence computing. Exact flux vector and flux projectionare assigned to mesh vertexes and mesh edges respectively, with which the fluxdivergence at vertexes are computed using the abovementioned least squaremethod. The analytical solution for the test is given by f ( x ) = sin ( πl x ) cos ( πl y ) g ( x ) = cos ( πl x ) sin ( πl y ) . (13)where l is the length scale of the computation domain. Four boundaries are allset to periodic boundary condition. By comparing numerical results and theanalytical divergence, three different kinds of error norms are evaluated. Alongwith the refinement of meshes, the decreasing tendency of all three error normsin Fig. 3 show fifth order accuracy, which meets with our design expectation. To achieve high order accuracy, interpolation or reconstruction based onhigh degree of polynomials is necessary. When using FDM in structured grids,interpolation or reconstruction is applied dimension by dimension. On the con-trary, in unstructured meshes, most methods apply interpolating/reconstructiondirectly in multidimensional space, which is much more complicated and ex-pensive. In addition, the operation should have WENO-like features so that8 Y
4 2 0 2 4 4 2024
Figure 2: Unstructuredmesh used to test divergencecomputing accuracy. e rr o r L1 error normL2 error normL error norm
Figure 3: Three kinds of er-ror norm of divergence com-puting accuracy test com-pared with the reference offifth order. discontinuities can be captured without oscillation. This means a number ofcandidate sub-stencils are considered and selected, making the operation evenmore complicated and expensive, especially in three-dimensional cases. This canexplain why FDM is one of the most efficient shock-capturing method among allavailable high order ones. Herein, the author tried to achieve one-dimensionalWENO interpolation in unstructured meshes so that the computing cost arelargely reduced compared to other methods.
In unstructured meshes, there is no explicit ordered grid lines as in structuralgrids. It is natural that people turn to multidimensional interpolation/reconstructiondirectly in each element (in DG, FR, and SD) or over several adjacent cells (inFVM). However, we can still assemble a string as smooth as possible by assem-bling several connected edges. Fig. 4 shows such a string assembled by 5 edges.Along such a string, one-dimensional interpolation can be readily applied. Flowstates w L and w R can be obtained by WENO5 interpolation at the middle point(the blue square) of the middle edge (the red edge). After that, numerical fluxat the middle point of the red edge can be computed with a Riemann solver.When performing WENO interpolation along the one-dimensional stencil,there are two choices. The first one is to interpolate in the generalized coordinate9 igure 4: A one-dimensional interpolation stencil for the red edge in a triangular mesh.WENO5 interpolation is applied along this stencil to obtain flow states w L and w R at themiddle point (the blue square) of the red edge.Figure 5: A curve stencil(the color curve) genrated from linked edges(the thin black lines). ξ . In the ξ domain, we assume that the length of the curve between two neighborvertices is the same as ∆ ξ = 1. However, it is not guaranteed that the midpointin ξ domain is the same midpoint as on the original edge. For nonuniformstencils constructed in unstructured meshes, the difference cannot be ignored.In order to make the interpolation right at the middle point, another choiceis interpolating along the curve length. As shown in Fig. 5, two curves areconstructed on each side of the middle edge. From each end of the middle edge,a curve starts with the same tangent direction as the middle edge and cross theother two vertices. Here, Catmull-Rom spline is adopted to construct the twocurves. The whole stencil is composed of the curves and a straight middle edge.Along this stencil, nonuniform WENO5 interpolation can be applied along thecurve length coordinate to obtain w L and w R .The one-dimensional interpolation used here is more efficient and can beapplied on characteristics variables, making the solution more stable compared10ith interpolating conservative variables. Besides, the interpolation and fluxcomputing methods used herein does not require a good-quality mesh as artificialedges can be used. Since interpolation stencils obtained from unstructured meshes are not equallyspaced, non-uniform WENO scheme of fifth order accuracy is derived. Detailsof the derivation is given in the appendix. The accuracy of the scheme is testedin this section. The same tests are conducted by using standard uniform WENOscheme for comparation. The exact solution is assumed to be a sinusoidal wavewhich reads f ( x ) = sin ( x ) . (14)A stencil with five equally spaced points is the base stencil of the test. Thesepoints can be moved randomly to left or to right by some percent of the originalinterval. If the points do not move, standard and non-uniform WENO showthe same order of accuracy as in Fig. 6(a), both presenting fifth order accuracy.By contrast, if points are moved randomly by 50 percent of original interval(which makes the stencil become non-uniform), non-uniform WENO scheme(NU-WENO5 in the figure) presents much better result compared with standardWENO scheme. Standard WENO5 scheme shows much larger interpolationerror and has around first order of accuracy, while non-uniform WENO5 keepsalmost the same small error as in the uniform-stencil case and preserves fifthorder of accuracy (See Fig. 6(b)). Therefore, it is necessary to use non-uniformWENO scheme if one intends to directly apply WENO schemes to non-uniformstencils without coordinate transformation. Another two accuracy test cases are conducted to check the performanceof standard WENO5 scheme and non-uniform WENO5 scheme along genuinecurves. It is already shown in section 2.3.2 that non-uniform WENO5 schemeperforms well in the case of one-dimensional stencil. In the following step, it is11 e rr o r NU WENO5U WENO55thOrder (a) Uniform stencils h e rr o r NU WENO5U WENO55thOrder1stOrder (b) Non-uniform stencilsFigure 6: L1 norm of the interpolation error of standard and non-uniform WENO schemestest with uniform and non-uniform stencils, where h is the interval between stencil nodes.NU-WENO5 stands for the non-uniform WENO scheme of fifth order of accuracy, while U-WENO5 means standard fifth order WENO scheme in reference [11]. needed to test whether the scheme preserves the accuracy in the case of two-dimensional curves, which are extracted from unstructured meshes similar tothe curve in Fig. 4. Standard and non-uniform WENO5 schems will be used tointerpolate the value at midpoint along these genuine curves.Each curve in Fig. 7(a) and Fig. 8(a) is essentially a part of a two-dimensionalfield. Therefore, the test function used herein is two-dimensional rather thanone-dimensional as in section 2.3.2, which reads f ( x, y ) = sin ( x ) ∗ cos ( y ) (15)Fig. 7 shows the curve used in the first test case and the result of the ac-curacy test. The curve in Fig. 7(a) is obtained by fifth-degree polynomial in-terpolation that connects six nodes with x = [ − . − . − . . . . , y =[3 . . . . − . − . O (∆ x ).Therefore, with the refining of the node interval ∆ x , the interpolation errordecrease in a fifth-order way. x y (a) Fifth-degree polynomial curve h e rr o r U WENO5NU WENO55thOrder (b) Interpolation error distributionFigure 7: The interpolation curve and WENO interpolation results for the ”refining” case.For the interpolation of the legend in this figure, the reader is referred to the figure legend inFig. 6.
In the second test case, it will be shown that WENO schemes no longerpreserves the high accuracy if the test case is conducted in a ”scaling” wayrather than in a ”refining” way.The curve in Fig. 8(a) is a two-dimensional curve constructed by the Catmull-Rom method exactly in the same way as in Fig. 5. It is naturally a part of atwo-dimensional field and is extracted from an unstructured mesh. Unlike in thefirst test case, a series of interpolation error is obtained by scaling the curve anddo the WENO interpolation, rather than refining the nodes on a curve. In thefirst test case, the function distribution along the whole curve will not change13uring the process of refinement. By contrast, the function value along the curvewill change when the curve is scaled as it is not at the same spatial position,which makes the interpolation error unpredictable. The result in Fig. 8(b) showsan accuracy of first order.The above two test cases in fact correspond to the numerical accuracy test forthe traditional structured grids, and for the method developed in this paper inthe context of unstructured meshes respectively. When one is doing the accuracytest for the structured grids, the grid line in a dimension basically remainsunchanged when the node distribution is refined, which exactly corresponds tothe abovementioned ”refining” case. However, for unstructured meshes, curvesextracted from meshes as in Fig. 5 are much more like scaled rather than refinedwhen the original unstructured meshes are refined. As a result, the accuracy testof the finite difference method developed in this paper for general unstructuredmeshes may have a similar performance as in Fig. 8, which shows only aroundfirst-order accuracy. In fact, this conjecture will be tested and verified in thefollowing paragraph. In section 3.1, the accuracy of spatial discretization erroris tested on unstructured meshes. It is found that the discretization accuracyis of around first order for general unstructured meshes. Further details will bedescribed in the corresponding section. x y (a) Catmull-Rom curves h e rr o r U WENO5NU WENO51stOrder (b) Interpolation error distributionFigure 8: The interpolation curve and WENO interpolation results for the ”scaling” case. . Numerical Examples Results of several inviscid problems are presented in this section to test theperformance of the developed method. These cases include the problem withsmooth solution and problems with discontinuities such as shock waves.In all of these cases, two layers of vertices and edges are used as the di-vergence stencil. One-dimensional six-point stencil is used for the fifth orderWENO interpolation to compute numerical fluxes at the midpoints of edges.Thus, the overall spatial discretization accuracy is of the order fifth. The thirdorder total variation diminishing Runge-Kutta scheme [16] (RK-TVD) is usedfor time integration.
The first test problem is from Yee et al. [23], where an isentropic vortexconvects in an inviscid free stream. The computation domain is of the size [-5,-5] × [5,5], in which an isentropic vortex is located at [0,0] at the initial time.The vortex is added to a mean flow with u ∞ = 1 and v ∞ = 1, the initial flowfield is given by ρ = (cid:20) − ( γ − β γπ e − r (cid:21) γ − , r = ¯ x + ¯ y , ( u, v ) =(1 ,
1) + β π e − r ( − ¯ y, ¯ x ) p = ρ γ where β is the vortex strength and the value of 5.0 is used. Here, (¯ x, ¯ y ) =( x − x c , y − y c ), where ( x c , y c ) is the center of the vortex at the initial time. Theentire flow field is required to be isentropic, thus for a perfect gas, p/ρ γ = 1.Periodic boundary conditions are used in both directions, so that the vortexconvects along the diagonal of the computational domain and reaches the initialposition after each period.We use uniform triangular mesh which is obtained by diagonalizing the struc-tured rectangular mesh of 50 ×
50, shown in Fig. 9.15 y Figure 9: Uniform triangular mesh used in isentropic vortex case.
Fig. 10 gives the results of calculation after four time periods. Densitycontours are presented in Fig. 10(a) and it’s distribution at the cut line y =0 is shown in Fig. 10(b) at the same time. It can be seen that after fourtime periods, the vortex arrives at the initial position and remains almost thesame distribution compared with the initial flow field, which indicates very lowdissipation.A sequence of triangular meshes similar to Fig. 9 with different resolutionsare used to test the spatial discretization accuracy of the present method. To re-move the error caused by time discretization, the computational error of the den-sity is recorded after just one small time step where CFL (Courant-Friedrichs-Lewy condition) number is set to 0.01. The result in Fig. 11, in terms of the L error norm of density, shows that the present method has the accuracy of fifthorder on regular triangular meshes as in Fig. 9, which agrees with the designexpectation.Another sequence of general triangular meshes generated by Delaunay methodare used to test the discretization accuracy. The unstructured meshes are simi-lar to the mesh in Fig. 2, but are refined in order while keeping the computationdomain unchanged as x ∈ [ − , , y ∈ [ − , a) x d e n s i t y
2 1 0 1 20.40.50.60.70.80.911.1
T=0T=4 (b)Figure 10: Isentropic vortex advection at T=4 periods: (a) 17 density contours from 0.5 to0.98; (b) density distribution at the line y = 0. h e rr o r L1 error norm5thOrder
Figure 11: Accuracy test of the isentropic vortex problem on regular triangular meshes. e rr o r L1 error normL2 error norm1stOrder
Figure 12: Accuracy test of the isentropic vortex problem on general triangular meshes. time step with CFL=0.01. The simulation results are compared to the exactsolution and L1 norm and L2 norm of density error are recorded. The resultsin Fig. 12 show the accuracy of around first order. The reason for this loworder accuracy lies in the way that unstructured meshes are refined during thetest. When those general unstructured meshes are refined, the interpolationstencils/curves on them are scaled rather than refined, which leads to the samesituation as in Fig. 8 in section 2.3.3 and is only first-order accurate. How toimprove the performance of the accuracy in the case of general unstructuredmesh is left to the future work. Herein the focus of the study lies in the ideaof constructing one-dimensional interpolation stencil on unstructured mesh. Itshows the potential to improve the accuracy and efficiency of the simulation withunstructured meshes and to be applied to cases with solution discontinuity.
This problem is an extension of the one-dimensional shock tube problemintroduced by Sod [17], where high pressure and high density are set a smallregion at the initial time. The computational domain has size x ∈ [ − , y ∈ [ − ,
1] and periodic boundary conditions are used in both directions. Frontadvancing method is used to generate the triangular mesh. The average meshsize is around 0.04 (with equally spaced 51 points on each boundary), and thereare 3058 solution points in total, see Fig. 13.18 y
1 0.5 0 0.5 1 1 0.500.51
Figure 13: The triangular mesh used for 2d Sod shock tube case.
The initial condition is given by( ρ, u, v, p ) = (1 , , , , | x | < . , | y | < . . , , , . , elsewhereThe ratio of specific heats γ is set to 1.4. With the time evolution, shockwaves, contact discontinuities and expansion fans are generated. Fig. 14 showsthe evolution of density at four different moments, specifically, from t = 0 to t = 0 . dt = 0 .
2. The solutions clearly capture shock wavesand contact discontinuities, but with a little overshoot, which may be causedby insufficient dissipation. To validate the solution, our results are comparedwith ones computed by WCNS scheme on structured grids, where midpoint-and-node-to-node difference (MND) scheme with sixth order accurate is used,which reads ∂F∂x (cid:12)(cid:12)(cid:12)(cid:12) i,j = 1∆ x (cid:20) (cid:16) ˜ F i +1 / ,j − ˜ F i − / ,j (cid:17) −
310 ( F i +1 ,j − F i − ,j )+ 130 (cid:16) ˜ F i +3 / ,j − ˜ F i − / ,j (cid:17)(cid:21) (16)where ˜ F i +1 / ,j are numerical fluxes at midpoints between nodes. These midpoint-fluxes are computed by fifth order accurate WENO interpolation combined withRoe flux solver, thus makes the overall accuracy fifth order. The Mach numberdistribution at time t = 0 . a) t=0.0 (b) t=0.2(c) t=0.4 (d) t=0.6Figure 14: Density contours from 0.07 to 0.25 of two-dimensional Sod shock tube problemfrom time t = 0 to t = 0 . a) (b) (c)Figure 15: 19 equally spaced Mach number contours from 0.1 to 1.0, at t = 0 .
4: (a) unstruc-tured mesh with 3058 solution points; (b) structured grid of the size 55 ×
55; (c) structuredgrid of the size 201 × using the method that we developed, while Fig. 15(b) is the solution obtainedon structured grid of 55 ×
55, which has a similar resolution to the one usedin Fig. 15(a). They present similar Mach number distribution, while unstruc-tured method shows less dissipation as the red arrow in Fig. 15(a) indicates,where four low Mach number regions are found. Those features are smeared inFig. 15(b) due to the numerical dissipation of the scheme. To further confirmthis conclusion, a finer structured grid with the size of 201 ×
201 is utilizedwith WCNS scheme. The solution is in Fig. 15(c). It can be found that theabovementioned low Mach number region becomes clearer than what is found inthe coarser grid, which verifies that the low-dissipation feature of the developedmethod is not a non-physical one. This case demonstrates that present methodhas the ability to capture strong discontinuities such as shock waves.
Shock-vortex interaction problem studied in [3] is numerically tested herein.This problem consists a stationary shock of Mach 1.2 and a strong isentropicvortex that convects downstream and interacts with the former. The initialconfiguration is shown in Fig. 16. The computation domain size is 80 ×
80 and x ∈ [ − , y ∈ [ − , ,
0) initially.The stationary shock wave is at the position of x = 0. Initial flow field is given21 D Post-shock Pre- shock vortex
X=0
Figure 16: Shock-vortex interaction problem configuration at the initial time. by ρ = (cid:16) − . ρ − M v e − r (cid:17) γ − p = 1 γ ρ γ ( δu, δv ) = M v e − r ( − ¯ y, ¯ x ) , ¯ x = x − x c , ¯ y = y − y c where M v = 1 .
0, ( x c , y c ) is the position of the vortex core, ie. (4,0). The ratioof specific heats is chosen to been γ = 1 . x = 0 . (a) (b)Figure 17: Unstructured triangular mesh used in shock-vortex interaction problem: (a) entiremesh; (b) enlarged view of the mesh near the shock-vortex interaction region. t = 16 is shown in Fig. 18(a). The shock isdistorted, multiple pressure waves and reflected shocks are produced. Thisproblem is also computed by WCNS method described in 3.2 on a uniformstructured grid for comparation. The structured grid has the resolution of 1024 × (a) (b)Figure 18: Continuous pressure contour from 0.9 to 1.15 at t = 16 of shock-vortex interaction:(a) computed by present method; (b) computed by WCNS on structured grid.
4. Conclusion
In this paper, we extend finite difference method to unstructured meshes.This method constructs one-dimensional stencils first to get numerical flux pro-jection along each edge, which is similar to the way WCNS does on struc-tured grids. It introduces the feature of highly efficient spatial discretizationof structured finite difference methods, which makes this method has the po-tential to largely reduce the computational cost in the context of unstructuredmeshes. Non-uniform WENO scheme is derived to deal with non-uniform one-dimensional stencils. Based on fluxes at edges and vertices, least square methodcomputes flux divergence and then Euler equations can be advanced readily.Owing to the use of fluxes both at edges and at vertices, the divergence stencilis quite compact. 23everal numerical tests are conducted and the results show that the presentmethod has the potential to accurately solve the flow field and to capture strongsolution discontinuities. In addition, the developed method can be easily imple-mented in existing unstructured flow solver frames. This method shows expectedaccuracy on meshes from diagonalizing uniform structured grids. However, thehigh order performance degenerates with arbitrary unstructured meshes, whichis tested and explain in detail. How to improve the interpolation accuracy willbe left to the future research.In this paper, the attention is focused on presenting a novel idea of finitedifference method in the context of unstructured meshes. The method con-structs one-dimensional stencils in unstructured meshes and uses least squaremethod to get the divergence of flux vectors, which form a compute simulationmethod for the inviscid compressible fluid problems. However, it does not giveabundant and detailed performance tests of the method. More detailed studyof the efficiency, accuracy and stability of this method will be conducted in fu-ture work. Besides, numerical examples studied in the present work all consistsimple boundary conditions such as periodic or Dirichlet boundary conditions.In our future work, problems in engineering level with walls or other generalboundary conditions will be studied.
Appendix
WENO interpolation on non-uniform stencils is described herein. Pointsin one-dimensional stencils assembled in unstructured meshes are not equallyspaced along the curve. Directly applying standard WENO interpolation [11] tothese stencils leads to a considerably large error. This is because the standardWENO is derived from a uniform stencil. To get rid of this problem, non-uniform WENO interpolation is derived in this section.In order to obtain a fifth order accurate WENO interpolation, a six-pointstencil is used to get the left and right states of an midpoint. Each state isnonlinearly weighted from the interpolates of three sub-stencils as shown in24 igure A1: Stencils for non-uniform WENO to approximate value u at x i +1 / . Fig. A1. The stencil is basically the same as the one in [11] except that non-uniform points are used herein. For simplicity, only the left state interpolationof midpoint x i +1 / is presented here. The right state can be obtained by a sim-ilar way using the mirrored stencil. The coordinate origin of points in Fig. A1is put at the midpoint x i +1 / i.e. ( x i + x i +1 ) /
2, so that x i +1 / = 0. All co-ordinates of points are normalized by the length scale of the stencil to reducethe error introduced by numerical operation. Values at nodes like x i et al. forinterpolation are written as u i et al. They are known before the upwind in-terpolation. These values can be primitive variables, conservative variables orother variables. Characteristic variables are used in this paper for the betterstability. Each sub-stencil approximates the solution with a quadratic polyno-mial using three known values. These polynomials are then evaluated at themidpoint x i +1 / = 0, which give˜ u ,i +1 / = a d u i − + a d u i − + a d u i , (A1)where a = x i − x i ( x i − x i − ) , a = − x i − x i ( x i − x i − ) , a = x i − x i − ( x i − − x i − ) d = ( x i − − x i − )( x i − x i − )( x i − x i − ) . ˜ u ,i +1 / = b d u i − + b d u i + b d u i +1 , (A2)where b = x i x i +1 ( x i +1 − x i ) , b = − x i − x i +1 ( x i +1 − x i − ) , b = x i x i − ( x i − x i − ) d = ( x i − x i − )( x i +1 − x i − )( x i +1 − x i ) . u ,i +1 / = c d u i + c d u i +1 + c d u i +2 , (A3)where c = x i +1 x i +2 ( x i +2 − x i +1 ) , c = − x i x i +2 ( x i +2 − x i ) , c = x i +1 x i ( x i +1 − x i ) d = ( x i +1 − x i )( x i +2 − x i )( x i +2 − x i +1 ) . It is noticed that the interpolated value at midpoint is the weighted summationof values at nodes, which has the same form as in the uniform WENO case,except that the coefficients of values now are calculated from coordinates ofpoints rather than constant numbers. The large stencil S of five points inFig. A1 utilizes a fourth-degree polynomial to approximate the solution. Fivecoefficients of the polynomial are determined by interpolation using values atfive points. The polynomial gives the value at x i +1 / as˜ u ,i +1 / = l d u + l d u + l d u + l d u + l d u , (A4)where l = x i − x i x i +1 x i +2 · ( x i − x i − )( x i +1 − x i − )( x i +1 − x i )( x i +2 − x i − )( x i +2 − x i )( x i +2 − x i +1 ) ,l = − x i − x i x i +1 x i +2 ( x i − x i − ) · ( x i +1 − x i − )( x i +1 − x i )( x i +2 − x i − )( x i +2 − x i )( x i +2 − x i +1 ) ,l = x i − x i − x i +1 x i +2 · ( x i − − x i − )( x i +1 − x i − )( x i +1 − x i − )( x i +2 − x i − )( x i +2 − x i − )( x i +2 − x i +1 ) ,l = − x i − x i − x i x i +2 · ( x i − − x i − )( x i − x i − )( x i − x i − )( x i +2 − x i − )( x i +2 − x i − )( x i +2 − x i ) ,l = x i − x i − x i x i +1 · ( x i − − x i − )( x i − x i − )( x i − x i − )( x i +1 − x i − )( x i +1 − x i − )( x i +1 − x i ) ,d =( x i − − x i − )( x i − x i − )( x i − x i − )( x i +1 − x i − )( x i +1 − x i − ) · ( x i +1 − x i )( x i +2 − x i − )( x i +2 − x i − )( x i +2 − x i )( x i +2 − x i +1 ) .
26t is noticed that coefficients in Eq. (A4) have a certain regularity and can besimply expressed as l k = ( − k i +2 (cid:89) p = i − ,p (cid:54) = i − k x p · i +2 (cid:89) m,n = i − ,m The authors acknowledge the support from the National Natural ScienceFoundation of China (Grant No. 91952203). References [1] T. Barth and P. Frederickson. Higher order solution of the Euler equationson unstructured grids using quadratic reconstruction. In , Reno,NV,U.S.A., Jan. 1990. American Institute of Aero-nautics and Astronautics.[2] R. Borges, M. Carmona, B. Costa, and W. S. Don. An improved weightedessentially non-oscillatory scheme for hyperbolic conservation laws. Journalof Computational Physics , 227(6):3191–3211, Mar. 2008.283] A. Chatterjee and S. Vijayaraj. Multiple Sound Generation in Interactionof Shock Wave with Strong Vortex. AIAA Journal , 46(10):2558–2567, Oct.2008.[4] B. Cockburn and C.-W. Shu. Tvb runge-kutta local projection discon-tinuous galerkin finite element method for conservation laws. ii. generalframework. Mathematics of computation , 52(186):411–435, 1989.[5] B. Cockburn and C.-W. Shu. Runge–kutta discontinuous galerkin meth-ods for convection-dominated problems. Journal of scientific computing ,16(3):173–261, 2001.[6] X. Deng and H. Zhang. Developing High-Order Weighted Compact Nonlin-ear Schemes. Journal of Computational Physics , 165(1):22–44, Nov. 2000.[7] A. K. Henrick, T. D. Aslam, and J. M. Powers. Mapped weighted essen-tially non-oscillatory schemes: Achieving optimal order near critical points. Journal of Computational Physics , 207(2):542–567, Aug. 2005.[8] C. Hu and C.-W. Shu. Weighted Essentially Non-oscillatory Schemes onTriangular Meshes. Journal of Computational Physics , 150(1):97–127, Mar.1999.[9] X. Hu, Q. Wang, and N. Adams. An adaptive central-upwind weightedessentially non-oscillatory scheme. Journal of Computational Physics ,229(23):8952–8965, Nov. 2010.[10] H. T. Huynh. A Flux Reconstruction Approach to High-Order SchemesIncluding Discontinuous Galerkin Methods. In , Miami, Florida, June 2007. American Instituteof Aeronautics and Astronautics.[11] G.-S. Jiang and C.-W. Shu. Efficient Implementation of Weighted ENOSchemes. Journal of Computational Physics , 126(1):202–228, June 1996.2912] G. Lacaze, B. Cuenot, T. Poinsot, and M. Oschwald. Large eddy simulationof laser ignition and compressible reacting flow in a rocket-like configura-tion. Combustion and Flame , 156(6):1166–1180, 2009.[13] X.-L. Li and Y.-X. Ren. High order compact generalized finite differencemethods for solving inviscid compressible flows. 82(1):18.[14] Y. Liu, M. Vinokur, and Z. J. Wang. Spectral difference method for un-structured grids i: basic formulation. Journal of Computational Physics ,216(2):780–801, 2006.[15] W. H. Reed and T. Hill. Triangular mesh methods for the neutron transportequation. Technical report, Los Alamos Scientific Lab., N. Mex.(USA),1973.[16] C.-W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shock-capturing schemes. Journal of Computational Physics ,77(2):439–471, Aug. 1988.[17] G. A. Sod. A survey of several finite difference methods for systems ofnonlinear hyperbolic conservation laws. Journal of computational physics ,27(1):1–31, 1978.[18] H. Su, J. Cai, K. Qu, and S. Pan. Numerical simulations of inert andreactive highly underexpanded jets. Physics of Fluids , 32(3):036104, 2020.[19] Q. Wang, Y.-X. Ren, and W. Li. Compact high order finite volume methodon unstructured grids II: Extension to two-dimensional Euler equations. Journal of Computational Physics , 314:883–908, June 2016.[20] Z. J. Wang, Y. Liu, G. May, and A. Jameson. Spectral difference method forunstructured grids ii: extension to the euler equations. Journal of ScientificComputing , 32(1):45–71, 2007.[21] F. Witherden, A. Farrington, and P. Vincent. PyFR: An open source frame-work for solving advection–diffusion type problems on streaming architec-30ures using the flux reconstruction approach. Computer Physics Commu-nications , 185(11):3028–3040, Nov. 2014.[22] M. L. Wong and S. K. Lele. High-order localized dissipation weighted com-pact nonlinear scheme for shock- and interface-capturing in compressibleflows. 339:179–209.[23] H. Yee, N. Sandham, and M. Djomehri. Low-Dissipative High-Order Shock-Capturing Methods Using Characteristic-Based Filters.