A shock-fitting technique for 2D unstructured grids

Computers & Fluids 38 (2009) 715–726

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

A shock-fitting technique for 2D unstructured grids Renato Paciorri a, Aldo Bonfiglioli b,* a b

Dipartimento di Meccanica e Aeronautica, Università di Roma, La Sapienza, Via Eudossiana, 18 00184 Roma, Italy Dipartimento di Ingegneria e Fisica dell’Ambiente, Università della Basilicata, Viale dell’Ateneo Lucano 10, 85100 Potenza, Italy

a r t i c l e

i n f o

Article history: Received 23 August 2007 Received in revised form 10 June 2008 Accepted 25 July 2008 Available online 9 September 2008

a b s t r a c t A new floating shock-fitting technique featuring the explicit computation of shocks by means of the Rankine–Hugoniot relations has been implemented on two-dimensional unstructured grids. This paper illustrates the algorithmic features of this original technique and the results obtained in the computation of the hypersonic flow past a circular cylinder and a steady Mach reflection. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction Shock waves occur very frequently in nature and technological applications. Their presence characterizes compressible flows not only in aeronautics and aerospace, but also in other areas of theoretical and applied physics. The extracorporeal shock wave lithotripsy, used for breaking kidney stones, and sonoluminescence, which is the emission of short bursts of light from imploding bubbles in a liquid excited by sound, are two such examples arising in different fields. In all flows where shock waves occur, these play an important role that affects the overall flow behavior. We could mention numerous examples testifying the importance of shock waves and their effects. Among these are: the flow around space vehicles re-entering the atmosphere at hypersonic speeds, the buffeting phenomenon induced by oscillating shock waves on transonic wings and turbomachinery blades, shock-induced noise and shock-induced boundary-layer separation that occur inside highly over-expanded nozzles. It is therefore not surprising that computing shocks correctly has always been one of the most important issues in the numerical simulation of compressible flows. The numerical models used in the simulation of compressible flows can be cast into two main categories: shock-capturing [1] and shock-fitting [2]. The technique for compressible flow computations known as shock-fitting was proposed and developed by Gino Moretti [2,3] in the 60s. It consists in identifying the shock as a singular line and computing the shock motion and upstream and downstream states according to the Rankine–Hugoniot (R–H) equations. Besides Moretti, other researchers applied this technique in the 70s and 80s (see, for instance, Refs. [4–7]) and in the early 90s [8,9]. In those years, the shock-fitting technique, used in conjunction with * Corresponding author. Tel.: +39 0971205203; fax: +39 0971205160. E-mail addresses: [email protected] (R. Paciorri), aldo.bonfi[email protected] (A. Bonfiglioli). 0045-7930/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2008.07.007

numerical schemes based on the quasi-linear equations, allowed to accurately compute flows with strong discontinuities using the modest computational resources available at that time. In the 90s, however, the development of modern shock-capturing schemes based on the integral conservation equations along with the strong growth in computing resources diminished the interest in shock-fitting schemes, which were considered too complex and not general enough as compared to shock-capturing schemes. The development and the application of the shock-fitting technique was pursued with good results by only one research group [10,11]. Despite the widespread use of shock-capturing codes and the continuous efforts made over the last 20 years, shock solutions obtained by means of shock-capturing schemes are plagued by a number of problems pertaining to accuracy [12–17], stability [18,19] and, more in general, solution quality. All these shortcomings have not been completely overcome and appear to be particulary severe when unstructured-grids are used for hypersonic calculations [20–24]. The shock-fitting technique, on the contrary, does not suffer from these numerical problems but, as it was said before, it has been abandoned over the last decades. We believe there exist today good reasons for reconsidering the shock-fitting technique especially in the context of unstructured grids. During their historical development, computational fluid dynamics (CFD) techniques have first been developed to use structured grids. Compared to unstructured-grid methods, the former offer a number of advantages, most notably a reduced algorithmic complexity which results in improved computational efficiency. Nevertheless, starting in the early 90s, the CFD community has shown an increasing interest in unstructured-grid techniques [25]. This is mainly due to the following features, which are peculiar to unstructured grids: (i) the ability to describe complex geometries and (ii) the possibility of locally adapting the grid to follow the flow features. This latter capability makes them particularly well suited to simulate compressible flows, which may develop discontinuities such as shock waves and contact discontinuities.

716

R. Paciorri, A. Bonfiglioli / Computers & Fluids 38 (2009) 715–726

Moretti’s shock-fitting technique was born at a time when CFD practitioners used solvers based only on structured meshes. The algorithmic complexity and the difficulties encountered in the implementation of the shock-fitting methodology within structured-grid solvers have been the main reasons for the loss of interest in shock-fitting. The traditional implementations adopted in structured-grid solvers were based on either the boundary-shock or the floatingshock concepts. Using the boundary-shock implementation, the structured gasdynamic solver considers the shocks as boundaries which split the domain into various sub-domains; each of these is then meshed using a structured-grid block. Therefore, in the boundary-shock approach the shock points are boundary mesh points whose nodal values are stored and updated by the gasdynamic solver. In the floating-shock implementation, on the contrary, the shocks move over and independently of a fixed background mesh. In this case, the shock points do not coincide with the grid points where the dependent variables are stored and, therefore, the additional problem arises of how to transfer correctly the informations from the shock points to the mesh points and vice versa. From a coding point of view, these two types of implementations differ considerably. Indeed, the boundary-shock approach is easy to implement because it does not introduce any modifications into the computational kernel of the code. On the other hand, the floating-shock approach is more complex and its implementation obliges to introduce corrections within the computational kernel of the CFD code. These difficulties can be somewhat mitigated by adopting special gasdynamic solvers, such as those based on the so-called k-scheme [26]. Despite the major simplicity of the shock-boundary implementation, problems arise since shock points may appear or disappear when shock waves interact or weaken into Mach lines or when compression waves coalesce into a new shock. In these cases, the definition of a mesh topology suitable to follow the evolution of all shock fronts quickly becomes a ‘‘topological nightmare” [27]. For these reasons, the few attempts to develop general-purpose shock-fitting codes [4,28,8,10] were always based on the floating-shock approach, whereas the shock-boundary concept was applied to many different codes mainly to solve blunt-body problems (see, for instance, Refs. [29–32]). To our knowledge, only a few attempts have been made to incorporate shock-fitting ideas within shock-capturing, unstructured-grid solvers. Van Rosendale [33] and Trepanier and co-workers [35] used explicitly the term ‘‘shock-fitting” in their papers’ titles, but none of these contributions can be considered as a genuine shock-fitting technique. The strategies adopted in Refs. [33,35] are similar: rather than explicitly solving the R–H equations to compute the shock, use is made of the property of Roe’s approximate Riemann solver [37] to return the exact jump relations when the discontinuity is aligned with the mesh. Therefore, both methods try to locally align the edges of the triangulation with the discontinuities by means of grid motion and/or adaptation. Following this, a shock-capturing scheme is used to compute the shocks on the adapted grid. For these reasons these schemes should be addressed as front-tracking rather than shock-fitting methods. Similar approaches were proposed by Parpia and Parikh [34] and, very recently, by Mavriplis [38]. A special attention has to be devoted to the work of Gloth et al. [36]. In this paper, the authors illustrate an original method for unstructured meshes that they classify as a front-tracking method. Despite of its label, this method incorporates many shock-fitting ideas although it differs from Moretti’s technique for at least one important aspect.

More specifically, the method proposed by Gloth et al. computes flows with shocks on a fixed grid which is not adapted during the computation. Contrary to Moretti’s original shock-fitting technique, in which floating shock points are moved according to the Lagrangian equation of motion, Gloth et al. use a level set method, which is an Eulerian approach. This method is based on the solution of a PDE (Hamilton–Jacobi type equation) for a scalar field function and the moving shock is a particular level of this scalar field. The shock front moves with a local speed which is computed by solving the R–H equations, with values on both sides of the discontinuity obtained by means of one-sided, zero-order projection from the neighbouring nodes of the fixed unstructured grid. A special treatment, which uses these states satisfying the R–H relation, is required to compute the flux integrals in those cells that are crossed by the shock front. According to the authors, the method described in Ref. [36] is affected by two drawbacks: it is first-order accurate in the neighbourhood of the shock front and it is not rigorously conservative. It is our opinion that these two drawbacks are closely related to the absence of floating shock points. In the approach proposed here, which borrows ideas from both the boundary-shock and floating-shock implementations of Moretti’s shock-fitting technique, the shock front is discretized as a polygonal curve and treated as an internal boundary by the CFD code: its local speed and downstream nodal values are computed by coupling the upstream moving wave predicted by the shockcapturing code with the R–H relations. The discontinuity is allowed to move over and independently of a background triangular grid which is used by the gasdynamic solver to compute the shockless regions. In the neighbourhood of the shock fronts, a local re-meshing technique is used to ensure that the edges of the shock are part of the computational grid. The present paper gives a detailed description of the proposed shock-fitting algorithm along with a brief discussion of the main features of the gasdynamic solver used in association with the proposed shock-fitting technique. Finally, results are presented for a normal shock moving through a constant section duct, the bow shock that is formed ahead of a circular cylinder moving at hypersonic speed and a Mach reflection. 2. Shock fitting algorithm The present shock-fitting algorithm is able to compute isolated shocks occurring within a two-dimensional flow field discretized by an unstructured mesh made of triangles. In its current version, the algorithm does not include any procedure for shock detection and shock interaction. Therefore, the fitted shocks must be present within the flow domain at the beginning of the simulation and no shock interactions between the fitted shocks are presently allowed. The lack of appropriate shock detection and shock interaction procedures reduces the generality of the present shock-fitting technique. Nevertheless, it is worth underlining that the use of this shock-fitting technique in conjunction with unstructured gas-dynamic solvers based on the conservative formulation reduces the negative impact of these restrictions. Indeed, those shocks which may eventually appear within the flow-field during the simulation would be captured by the conservative solver and the interaction between a fitted shock and a captured shock may also arise and would be properly handled by the current algorithm, as demonstrated in [39] and in the last example of Section 4. In order to illustrate the algorithmic features of the proposed methodology, let us consider a two-dimensional domain and a shock front crossing this domain at time level t (see Fig. 1). The shock front is shown in bold by a polyline whose endpoints are the shock points, marked by squares. A background triangular

R. Paciorri, A. Bonfiglioli / Computers & Fluids 38 (2009) 715–726

717

2.2. Local re-meshing around the shock front Shock

Upstream

Downstream

Fig. 1. Shock front moving over the background triangular mesh at time t.

mesh, whose nodes are denoted by circles, covers the entire computational domain. The background mesh is built using an hybrid frontal-Delaunay algorithm developed by Müller et al. [40,41]. It is worth observing that the position of the shock points is completely independent of the location of the nodes of the background grid. While each node of the background mesh is characterized by a single set of state variables, two sets of values, corresponding to the upstream and downstream states, are assigned to each shock point. We assume that at time t the solution is known at all grid and shock points. The computation of the subsequent time level t + Dt can be split into several steps that will be described in detail in the following sub-sections. 2.1. Cell removal around the shock front

2.3. Computation of the tangent and normal unit vectors

The first step consists in removing cells around the shock front. More precisely, all cells of the background mesh that are crossed by the shock line are removed along with the mesh points which are located too close to the shock front (see Fig. 2). Specifically, if the relative distance of a mesh point from the nearest shock edge satisfies the following inequality: min

d=l 6 rlsh

where l is the shock edge length, d the grid-point distance from the min shock edge and rlsh the minimum relative distance preset by the user, then the node is declared as a ‘‘phantom” node and it is removed from the mesh. Finally, all cells having at least one phantom node among their vertices are also removed from the background triangulation.

d

Phantom nodes

Following the cell removal step, the remaining part of the background mesh is cut by a hole containing the shock, as shown in Fig. 3a. The hole is then re-meshed using a constrained Delaunay triangulation (CDT): both the shock segments and the edges belonging to the boundary of the hole are constrained to be part of the final triangulation. No further mesh point is added during this triangulation process which is performed using the TRIANGLE [42] code developed by Shewchuk. It is important that this triangulation process generates a good quality mesh, meaning that the modified mesh around the shock front has to be composed by triangles which do not tend to degenerate into segments. To obtain a good quality triangulation two conditions have to be met: (i) the mesh points must not be located too close to the shock edges, and (ii) the shock edge length has to be approximately equal to the local length of the triangle edges. These conditions can be satisfied by means of the action of the first (2.1) and the last (2.9) steps of the shock-fitting algorithm. In the first step (cell removal around the shock front), the choice of an appropriate value of min rdsh (ranging between 0.2 and 0.5) allows to remove the gridpoints that are too close to the shock front. The length of the shock edge can then be set approximately equal to that of the triangle edges in step 2.9 (shock point redistribution) by assigning an appropriate value to the shock edge length. At this stage, the computational domain is discretized by a computational mesh in which the shock points and the shock edges are part of the triangulation. This mesh differs from the background mesh only in the neighbourhood of the shock front (compare Fig. 3a and b).

l

Cells enclosing phantom nodes

Cells crossed by shock

Fig. 2. Dashed lines mark the cells to be removed; dashed circles denote the phantom nodes.

The tangent and normal unit vectors along the shock front have to be calculated within each shock point since this is required to properly compute the jump relations. The computation of these vectors in a generic shock point is carried out through specific finite-difference formulas which involve the coordinates of the shock point itself and of its neighbouring shock points. Let us consider the shock point i whose position at time level t is denoted by Pti (see Fig. 4). Shock points i + 1 and i  1 are located on both sides of shock point i. Their position Ptiþ1 and Pti1 at time level t can be used to compute the tangent and normal unit vector in the shock point i. A preliminary test is necessary in order to verify if these adjacent shock nodes fall within the domain of dependence of shock point i. Let us fix the attention on node i + 1. If a disturbance originating from shock node i + 1 is able to reach shock point i, then point i + 1 belongs to the domain of dependence of point i. This is checked through the following simple test. Let us consider the tangent unit vector s of the shock edge i + 1/2 oriented from point i to i + 1. If:

utdiþ1  siþ1=2  atdiþ1 < 0 utdiþ1

ð1Þ

atdiþ1

and being the flow and acoustic velocity of the downstream state of shock point i + 1 at time level t, then the shock point i + 1 belongs to the domain of dependence of the shock point i. After having repeated this test in the shock point i  1, three different situations may occur: (1) both shock points i + 1 and i  1 are in the domain of dependence of the shock point i; (2) only the shock point i  1 is in the domain of dependence of the shock point i;

718

R. Paciorri, A. Bonfiglioli / Computers & Fluids 38 (2009) 715–726

a

b Shock

Hole boundary

Hole Hole

Fig. 3. (a) The background mesh is split in two parts by a hole which encloses the shock. (b) The triangulation around the shock has been rebuilt.

t i+1 ad

t

ud

.τ

1

i+

2

τ i+1/2

Pti+1

2

ð4Þ

2

while n is obtained from s as before. Eq. (4) involves the shock points i, i  1 and i  2, whereas the shock point i + 1 does not appear. Also in this case the accuracy is second order. Finally, the third case is specular to the second one. Indeed, Eq. (4) is replaced by a similar formula which involves the shock points i, i + 1 and i + 2.

1/

i+1

2

i+

udt

 ƒƒƒ! ƒƒƒ! vs ¼ P i1 Pi li1 þ li3  Pi2 P i li1

l i+1/2 t

Pi

2.4. Solution update using the capturing code

l i−1/2 t Pi−1

Upstream

Downstream l i−3/2 t Pi−2

Fig. 4. Test to check if the point Pi+1 belongs to the dependence domain of Pi.

(3) only the shock point i + 1 is in the domain of dependence of the shock point i. When the first situation occurs, the computation of the vectors can involve the adjacent shock points on both sides. Specifically, the tangent unit vector s is obtained by:

s¼

vs jvs j

ð2Þ Internal boundary

where vs is a vector defined as:

ƒƒƒ! ƒƒƒ! vs ¼ Pti P tiþ1 li1 þ Pti1 P ti liþ1 2

2

ð3Þ

ƒƒƒ! ƒƒƒ! Pti Ptiþ1 and Pti1 P ti being two vectors joining the shock points, liþ1 and li1 the length of the shock edges adjacent to the shock point 2

The solution is updated at time level t + Dt using the unstructured, gasdynamics code that will be described in the next section. The gasdynamics solver uses the computational mesh generated in step 2.2 as input and, therefore, the computational domain is cut in non-communicating parts by the shock which is treated by the code as if it were an internal boundary (see Fig. 5). This is achieved by replacing each shock point by two superimposed mesh points: one belonging to the downstream region and the other to the upstream region. The downstream and upstream states are correspondingly assigned to each pair of mesh nodes associated to each shock node. A single time step calculation is performed by the unstructured solver which provides updated nodal values within all grid and shock points at time level t + Dt. However, the downstream state of the shock points is not correctly computed except for the following combination of variables:

Downstream state

Shock point Upstream state

2

i. These relations guarantee second-order accuracy in the computation of s, even if the shock edge lengths are different. The normal unit vector n is obtained from s by imposing that n is orthogonal to s. The orientation is chosen such that n points from the shock downstream region to the upstream region. In the second situation, the shock point i + 1 must not be used in the computation of the tangent and normal vectors. The vector s is computed using Eq. (2), having replaced Eq. (3) with the following one-sided formula:

Fig. 5. The shock is treated as an internal boundary.

719

R. Paciorri, A. Bonfiglioli / Computers & Fluids 38 (2009) 715–726

RdtþDt ¼ adtþDt þ

c1 2

udtþDt  n

ð5Þ

where n is the shock normal, adtþDt and udtþDt are the values of the acoustic and flow velocity of the downstream state of the shock nodes computed by the unstructured solver. It is worth observing that RdtþDt is correct even if the quantities adtþDt and udtþDt appear to be wrong. Indeed, Rd is associated to the acoustic wave which moves from the downstream region towards the shock and the numerical scheme implemented in the unstructured code is able to compute correctly the evolution of this wave without requiring any informations upstream of the shock. 2.5. Shock calculation Each shock point is characterized by an upstream state, a downstream state and a shock speed (w). The upstream state has been updated at time level t + Dt in the previous step 2.4. Concerning the downstream state, only the quantity Rd has been correctly updated. The complete downstream state and the shock speed are computed through a system of five algebraic equations. Four equations are the Rankine–Hugoniot relations:

qdtþDt ðudtþDt  n  wÞ ¼ qutþDt ðuutþDt  n  wÞ qdtþDt ðudtþDt  n  wÞ2 þ pdtþDt ¼ qutþDt ðuutþDt  n  wÞ2 þ putþDt c pdtþDt ðudtþDt  n  wÞ2 c putþDt ðuutþDt  n  wÞ2 þ þ ¼ tþ D t c  1 qd c  1 qutþDt 2 2

consider that in the subsequent time steps the shock front might move sufficiently far away from its current position that some of the phantom nodes may be re-inserted into the computational mesh. The update of the phantom nodes is performed by first locating each phantom node within the computational mesh (a search operation which can be efficiently performed on a Delaunay mesh) and then using linear interpolation which is consistent with the linear elements adopted in the unstructured code. In Fig. 6 the background mesh (dashed lines) is superimposed to the computational mesh (continuous lines). This plot makes evident the position of the phantom nodes, the cells of the computational mesh which contain them and the mesh points involved in the interpolation. After updating the phantom nodes, the computational mesh used in the current time step has finished its task and it can be completely removed. 2.7. Shock displacement The new position of each shock point at time t + Dt is computed using the shock speed calculated in the previous step 2.5 by the following first-order accurate integration formula:

PtþDt ¼ Pt þ wtþDt nDt: ð6Þ

udtþDt  s ¼ uutþDt  s where p and q are, respectively, pressure and density, n and s are the normal and tangent unit vector computed in step 2.3, the subscript ‘d’ denotes the downstream state and ‘u’ the upstream state. Finally, the last equation of the system is the variable combination of Eq. (5) where the value Rd has been computed in the previous step 2.4. These equations constitute a system of five equations with respect the unknown quantities adtþDt , qdtþDt , uutþDt and wtþDt . The solution of this system of equations within each shock point provides the correct downstream state and the shock speed at time level t + Dt. 2.6. Interpolation of the phantom nodes In the previous steps, all nodes of the computational mesh have been updated at time level t + Dt. Nevertheless, the phantom nodes, which belong to the background mesh but had been excluded from the computational mesh, have not been updated. The need to also update these nodes should be clear if we

ð7Þ

The low order of accuracy of this formula does not affect the accuracy of steady state solutions which depends on the space-accuracy of the gasdynamics solver and of the tangent and normal vector computation. Therefore, Eq. (7) can be used along with higher order solvers for steady computations. On the contrary, the computation of unsteady flows would require an improvement of the accuracy of the integration scheme used to compute the shock motion: it should have at least the same order of time accuracy as that of the unstructured solver. The new position of the shock line is obtained (see the dashed line in Fig. 7) by means of the displacement of all shock points. 2.8. Interpolation of the jumped nodes During the shock displacement step (2.7), it may occur that the shock sweeps over several nodes of the background mesh. For instance, in Fig. 8 a mesh node whose position at time level t was located in the upstream region, is jumped by the shock, so that its position at time t + Dt is in the downstream region. The state of these ‘‘jumped” nodes has to be correctly re-computed. This is performed through linear interpolation. Specifically, the interpolation t

Pι+1

Pιt

t+ Δ t

P ι+1

t+ Δ t

t w ι nΔ

t+ Δ t

Pι

Phantom node

Cell containg the phantom node

Fig. 6. Interpolation of the phantom nodes.

Fig. 7. During shock motion some of the phantom nodes might be re-inserted into the mesh and some background nodes might be overcome by the shock front.

720

R. Paciorri, A. Bonfiglioli / Computers & Fluids 38 (2009) 715–726 t

Pι+1

t+ Δ t

P ι+1

jumped node

jumping edge

A t+ Δ t

Pι

t

Pι

during the temporal evolution, the overall number of shock points is also likely to change during the shock point re-distribution process. The states of the new distribution of shock points are computed from the previous one by using linear interpolation based on the curvilinear coordinate measured along the shock line. At this stage the numerical solution has correctly been updated at time level t + Dt, taking into account the shock front displacement. The next time level can be computed re-starting from the first step 2.1 of the algorithm. 3. Unstructured solver

Fig. 8. Mesh points jumped by the shock and their re-computation.

uses the states of the two shock points belonging to the shock edge that has swept the mesh point. In the case shown in Fig. 8, the swept node is re-computed by interpolating between the downstream states of the shock points marked in red. The distance between the projection (point A) of the jumped grid point onto the Dt tþDt and P iþ1 is used to compute shock edge and the shock points Ptþ i the weights used in the interpolation formula. 2.9. Shock point redistribution During its motion the shock front may shorten or lengthen and this can cause remarkable variations of the shock edge lengths. When this occurs, it is appropriate to redistribute the shock points along the shock line in order to make sure that the shock edge lengths are approximately equal to the edges of the neighbouring cells. As it was said in Section 2.2, this is one of the two conditions needed to obtain a good quality mesh during the re-meshing process around the shock. The shock point re-distribution along the shock line can be carried out imposing, for instance, that all shock edges have the same fixed length lsh (see Fig. 9). This length scale, which is preset by the user, can be determined from the characteristics of the background mesh. If the length of the shock line varies

lsh

As described in Section 2, the proposed shock-fitting algorithm provides a time dependent mesh in which the shocks are treated as internal boundaries of zero thickness. Communications with the CFD code are only required in step 2.4, when the modified grid containing the shock front is fed to the code which returns provisional values of the nodal variables at the subsequent time level. This reduced amount of communication allows treating the CFD code as a black box, thus making the proposed shock-fitting algorithm reusable in conjunction with virtually any vertex-centred unstructured CFD code. It is, however, worth describing the particular CFD solver used to produce the numerical results presented in Section 4. Specifically, the CFD solver is EULFS, a 2D/3D unstructured Euler/RANS code developed by Bonfiglioli [43]. In the test cases presented in Section 4, which are restricted to two spatial dimensions, the governing conservation equations are the Euler equations for a non-viscous, non-conducting, perfect gas. The computational domain is partitioned into a set of nonoverlapping triangular cells. Roe’s parameter vector [37] is chosen as dependent variable: it is stored in the vertices of the mesh and assumed to vary linearly and continuously in space, just as with linear finite elements (FEs). The control volumes over which the system of conservation equations is discretized are the median dual cells which are obtained by joining the centroids of gravity of all the cells surrounding a given mesh-point with the midpoints of all the edges that connect that mesh-point with its nearest neighbours. This is schematically shown in Fig. 10 where the solid lines denote the boundaries of the triangular cells and the dashed lines those of the control volumes. To enhance accuracy, the Euler equations are solved in preconditioned form, using the van Leer, Lee, Roe matrix [44]. This approach is often referred to as the hyperbolic–elliptic (HE) splitting [45,46] as it allows to achieve full or partial (depending on the flow regime) diagonalisation of the original unsteady Euler system. In supersonic flows, this is re-cast into a set of four scalar

lsh j

lsh

M ij M

jk

G

lsh Shock nodes before the redistribution

k

i

M ik

Shock nodes after the redistribution

Fig. 9. Shock point redistribution.

Fig. 10. Control volumes in the unstructured grid code.

721

R. Paciorri, A. Bonfiglioli / Computers & Fluids 38 (2009) 715–726

equations describing convection of entropy and total enthalpy along the streamlines and convection of the acoustic variables along the Mach lines. In subsonic flows, the acoustic component is described by a coupled 2  2 Cauchy–Riemann system, while entropy and total enthalpy are convected along the streamlines, just as in supersonic flows. The spatial discretization relies upon the so-called residual distribution (RD) (or fluctuation splitting) schemes, see Refs. [47,48], which share common features with both FE and FV techniques. RD schemes operate by evaluating the flux integrals (or cell residuals) over each triangular cell and distributing them among the forming vertices. The nodal residual, which is driven to zero as steady state is approached, is obtained by assembling the fractions of cell residuals scattered from the elements surrounding the node. A distinctive feature adopted in many RD codes is that the flux balance over a triangular cell is not computed using numerical quadrature, but using the quasi-linear form of the equations and a multi-dimensional generalisation of Roe’s linearisation [49]. Contrary to most FV schemes, no Riemann problems are involved in RD schemes since the functional representation of the dependent variables is continuous across the elements. The discretization of the viscous fluxes can be shown to be equivalent to a standard FE discretization, while the inviscid fluxes are distributed in an upwind-biased fashion. These upwinding directions are based upon the eigen structure of the quasi-linear form of the equations, as described in details in [43]. Concerning time integration, a simple Euler explicit scheme has been used and mass lumping applied to the time derivative term.

is placed at x = 0.2 and the shock line is made of 50 shock edges of equal length. Fig. 12 plots the initial density distribution (t = 0) along the line y = 0.5 and that computed by the EULFS code using the shock-fitting technique (shock-fitting mode) at t = 0.08 after 10 integration time steps. The diagram also reports the exact value of the density in the upstream and downstream regions and the exact trajectory of the shock. For comparison, the solution obtained by EULFS without shock-fitting (shock-capturing mode) over the background mesh (dashed line) is also plotted. It is evident that for this special case the ‘‘shock-fitted” solution is exactly superimposed to the analytical one. This is due to the problem peculiarities. Indeed, the unstructured solver provides exact values when the computed regions are uniform and the shock-fitting technique becomes exact when the shocks are straight, the computation of the normal and tangent unit vectors being errorless. Fig. 13 shows an enlargement of the background mesh and the computational ones used at two different time levels. It is worth to observe that when the shock arrives at a specific position, the mesh in the proximity of the shock is modified to take into account the shock points and edges. Once the shock has moved away, the mesh resumes its original configuration. 4.2. Blunt body problem The second test case deals with the high speed flow past a circular cylinder at free-stream Mach number, M1 = 20. This test is more complex from a fluid dynamics point of view than the previous one and it represents a significant test for the proposed shockfitting technique. Also the present calculations have been performed using the EULFS code both in shock-capturing and shockfitting mode. The computational domain surrounds a half circular cylinder having radius R = 1 (see Fig. 14). Two grids of increasing spatial resolution (see Table 1) have been created using the DELAUNDO [41] mesh generator by specifying the distribution of the boundary nodes. These are evenly spaced and the distance between two adjacent boundary nodes of the coarse mesh equals 0.08. The fine mesh is obtained by halving the distance between two adjacent boundary nodes and re-meshing the interior of the computational domain. For this reason the number of cells and grid points is approximately four times that of the coarse mesh, whereas the number of shock points is nearly twice.

4. Applications 4.1. Shock moving through a constant section duct The first test case considers a normal shock wave moving through a constant section duct. The shock relative Mach number is 1.89. This test case is very simple from a fluid dynamic point of view, but it is useful to explain the working principles of the present shock fitting technique. Fig. 11a reports the sketch of the test case and the non-dimensional upstream and downstream states. The computational domain is a square of non-dimensional unit length (see Fig. 11b) and it is discretized by a background mesh made of 1969 nodes and 3772 triangles. Initially, the normal shock

a

b

1 Initial shock position

0.8

Shock

w = 1.89 p2 = 4.001 a 2 = 1.265 u 2 = 1.134

Y

0.6

Upstream

Downstream

0.4

p1 = 1 a 1= 1 u1= 0

0.2

0

0

0.2

0.4

X

0.6

0.8

Fig. 11. Test case 1: (a) test case sketch and (b) computational domain and background mesh.

1

722

R. Paciorri, A. Bonfiglioli / Computers & Fluids 38 (2009) 715–726

Density distribution at y=0.5 t=0.08

t=0

1

4 3.5 3.5

ρ = 3.5

t=0.08

ρ = 3.5

t=0

0.8

3 3

0.6

2.5

2

X=

2

w

t+

0.

0.4

2

t

ρ

ρ

2.5

0.2

1.5

ρ = 1.4

1.5

ρ = 1.4

0

1 1 -0.4

-0.2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

X Fig. 12. Density distribution at y = 0.5.

0.9

3

Computational domain

Background mesh

0.8

0.7

2

0.1

0.2

0.3

0.4

1

0.5

Y

0.9

Modified mesh at t1=0

0.8

0.7

M ∞=20

Circular cylinder

0

-1 0.1

0.2

0.3

0.4

0.5

-2

0.9

final shock postition

Initial shock postition x=-0.242415 y2 + 1.467

-3

Modified mesh at t2 > t1 0.8

-2

-1

0

X

1

2

3

4

Fig. 14. Computational domain: initial and final shock position.

0.7

0.1

0.2

0.3

0.4

0.5

Fig. 13. Moving shock wave: background and computational meshes at different time levels.

Table 1 Nodes and cells of the background and computational meshes Level

The solution computed by the unstructured code in shock-capturing mode is used to initialize the flow-field and to determine the initial position of the shock front. In the shock-fitting simulation, the initial upstream state in the shock points is set equal to the free stream conditions, while the initial downstream state is computed from the upstream state and the shock slope assuming zero shock speed (w = 0). Successively, the flow-field and the shock position are integrated in time until steady state is reached. Fig. 14 shows the shock displacement which occurs between the initial position and the one reached at steady state. During the computation the number of nodes of the computational mesh varies with respect to that of the background mesh. For instance, the coarse background mesh is made of 351 nodes and 29 shock points are preset

Coarse Fine

Background

Initial

Final computational

Mesh

Shock nodes

Mesh

Cells

Nodes

610 2544

351 1365

29 61

Shock nodes

Cells

Nodes

624 2580

358 1383

29 57

at the beginning of the simulation, whereas the computational mesh used in the last time iteration is made of 358 grid points, 29 of which are shock points. Fig. 15 displays the coarse background mesh and the computational mesh used during the last iteration. The box on the right, which shows an enlargement of the near shock region, allows comparing the two meshes. It can be observed that these are superimposed everywhere except in

723

R. Paciorri, A. Bonfiglioli / Computers & Fluids 38 (2009) 715–726 Enlargement Background mesh

Computational mesh of the last time-iteration

3

3 Enlargement box

2

2

1

1

0

0

Y

Y

Enlargement box

-1

-1

-2

-2 Shock

-3 0

1

-3 0

2

X

1

2

X

Fig. 15. Comparison between the coarse background and the computational mesh at the last time-iteration.

2D flow around a circular cylinder at M ∞=20 pressure isolines (Coarse mesh)

3

2

2D flow around a circular cylinder at M ∞=20 pressure isolines (Fine mesh)

3

2 fitted shock

1

Y

the region adjacent to the shock, where the differences between the background (dashed lines) and the computational mesh (continuous lines) are due to the addition of the shock points and edges. Fig. 15 clearly shows that the re-meshing technique does not significantly increase the number of points and cells with respect to the background mesh, since the addition of the shock nodes in the computational mesh is partly balanced by the removal of the phantom nodes. Figs. 16 and 17 show the comparisons between the solutions computed by the EULFS code working in shock-capturing and shock-fitting mode on the coarse and fine meshes, respectively. Figs. 16 and 17 also include a detailed view of the stagnation point region showing the modifications introduced by the shock-fitting algorithm to the background mesh in order to take into account the shock line. The improvements in solution quality due to the

0

-1 captured shock

-2

0

1

2

3

X

4

5

6

Fig. 17. Fine grid solutions: comparison between the shock-capturing and shockfitting solutions.

fitted shock

Y

1

0

-1 captured shock

-2

0

1

2

3

4

5

6

X Fig. 16. Coarse grid solutions: comparison between the shock-capturing and shockfitting solutions.

shock-fitting technique are evident in the aforementioned figures; particularly when comparing the coarse grid solutions (see Fig. 16). More specifically, the coarse grid solution computed in shock capturing mode is characterized by a very large shock thickness and by strong spurious disturbances affecting the shock layer region. These disturbances originate from the shock and are caused by the mis-alignment between the mesh and the shock. On the contrary, the shock-fitting solution shows a very neat field. Moving on to the analysis of the fine grid solutions, Fig. 17, the one computed in shock-fitting mode appears to be very similar to that obtained on the corse mesh, whereas the shock-capturing solution still shows the presence of spurious oscillation even if the shock thickness has been halved. Also in this case, the fitted solution shows a significant gain in terms of solution quality with respect the shock-capturing solution. Of course, the resolution of the

R. Paciorri, A. Bonfiglioli / Computers & Fluids 38 (2009) 715–726

2

shock 45° 1

300

r 0

200

1

2

1

1.2

1.4

1.6

1.8

r Fig. 20. Normalized pressure along the 45° line: comparison among the computed solutions and the reference solution.

Y

4.3. Steady Mach reflection 1

1

0

0.8

Θ

d sh 0

1

X

2

3

0.6 0.4 Coarse Fine Ref. [50]

0.2 0

30

Θ

60

90

Fig. 18. Shock-wall distance: comparison among the computed solutions and the reference solution.

500

The present test case demonstrates that the proposed methodology allows also a hybrid mode of operation in which some shocks are fitted whereas all others are captured. The chosen test case consists in a steady Mach reflection, as schematically shown in Fig. 21. A uniform, supersonic (M1 = 2) stream of air undergoes a H = 14° deflection through a shock. Regular reflection of the incident shock is impossible for the chosen set of parameters (M1, H) so that a steady triple point arises which joins the incident and reflected shocks, the Mach stem and the contact discontinuity. The simulation has been carried out both in fully capturing and hybrid modes. In this latter case, the incident shock has been captured, whereas the reflected shock and the Mach stem have been fitted. In both calculations the slip line has been captured. These two solution are compared in Fig. 22. Fig. 22 displays an overall view of the Mach iso-contours, a zoom of the region surrounding the triple point and the Mach number distribution along a vertical line located at x = 1.5 within the subsonic region which originates downstream of the Mach

400

Computational domain

p/p∞

300

200

M∞=2.0 Coarse Fine Ref [50]

100

0

X

Ref. [50] S-C coarse S-C fine S-F coarse S-F fine

100

0

body 0

2

1.2

dsh

3

shock

On the contrary, the differences between the shock fitting solutions and the reference one are very small and the fine grid solution is nearly superimposed to the reference one.

1.4

0

body

400 Y

shock-capturing solutions could be significantly improved if local grid-refinement was used in the shock region, but this typically involves a significant increase in the number of cells and nodes. A quantitative analysis of the shock-fitting solutions was carried out by comparing the estimates of the shock position and the pressure distribution computed on both grid levels with the reference solution computed by Lyubimov and Rusanov [50]. In Figs. 18 and 19 the distance (dsh) between the shock and the cylinder wall and the normalized pressure at the wall are plotted against the azimuthal angle (H). This comparison shows that the solutions computed by the proposed shock-fitting methodology are not only in good agreement with the reference solution, but nearly superimposed. Moreover, this analysis clearly proves that the coarse grid solution is in practice grid-independent in spite of the extremely limited number of cells enclosed in the shock layer. A comparison between the shock-fitting and the shock-capturing solutions is displayed in Fig. 20 where the normalized pressure distributions along the line at 45° is plotted. The shock-capturing solutions are characterized by a finite shock thickness that significantly affects the distribution even on the fine grid. Moreover, visible differences between the shock-capturing solutions and the reference one are also visible in regions far from the shock and close to the body.

p/p∞

724

0

30

Θ

60

90

Fig. 19. Normalized pressure at wall: comparison among the computed solutions and the reference solution.

Θ=14° Fig. 21. Steady Mach reflection: flow configuration.

725

R. Paciorri, A. Bonfiglioli / Computers & Fluids 38 (2009) 715–726 Mach reflection M=2, Θ = 14° Mach number conturlines

Enlargement around the triple point

0.2

Mach number distibution at x=1.5

0.2

0.1

0.1

0

0

0

Y

1

Y

Y

Captured solution

-0.1

-1

-0.1

hybrid solution

0

0.5

Captured solution

hybrid solution

1

1.5

2

-0.2 0.6

X

0.7

0.8

-0.2 0.7 0.8 0.9

1

1.1

M

X

Fig. 22. Steady Mach reflection: Mach iso-contours and detail of the triple point region.

stem. Remark that the hybrid solution produces a smoother Mach number distribution downstream of the reflected shock wave and the Mach stem. Moreover, observe the meandering shape taken by the captured Mach stem on the isotropic triangular grid (generated with DELAUNDO [41]) made of almost equilateral triangles. It can be observed that the meandering shape of the captured Mach stem gives rise to a strongly non-uniform Mach number distribution which is preserved downstream, as can be seen by the saw-tooth Mach number profile displayed on the right of Fig. 22. On the contrary, in the hybrid solution the fitted Mach stem produces the expected smooth Mach distribution which is also displayed in Fig. 22. As mentioned previously, this result proves the code capability to operate in hybrid mode. Nevertheless, it is important to underline that this is not the unique possible approach to compute shock interactions. Indeed, the proposed shock-fitting algorithm can be extended by fitting shock interactions (such as those occurring at triple points) as well as contact discontinuities. Work is currently underway in this direction and interesting results will be published in forthcoming papers. 5. Conclusions The shock-fitting technique illustrated in this paper has been successfully implemented in conjunction with an unstructured cell-vertex solver. The numerical simulation of two simple test cases characterized by the presence of an isolated shock wave has demonstrated that the proposed shock-fitting technique provides numerical solutions characterized by very small numerical errors (in the special case of a moving normal shock it returns the analytical solution) even when the computational mesh is very coarse. These results closely follow those obtained in the past when the shock-fitting technique had been used in conjunction with structured-grid solvers. This observation is not surprising because of the strong similarities between our shock-fitting technique and those used in the past. It is our opinion, however, that the proposed shock-fitting technique for unstructured solvers has a greater potential for development with respect to that displayed in the past by the shock-fitting techniques used in conjunction with structured meshes. On the one hand, the present approach borrows several features from both the shock-boundary and shock-floating implementations adopted in the past in structured-grid solvers based on the shock-fitting technique. On the other hand, it introduces several original features that greatly simplify its algorithmic implementation.

Indeed, the shock is treated as an internal boundary by the unstructured gasdynamic solver so that no modifications in the computational kernel of the CFD code are required. Nonetheless, the shock points can freely move over the mesh points of the background grid, as it was done in the floating-shock fitting approaches proposed in the past. To achieve this, the computational mesh is locally re-meshed at each time-step using the background grid and the shock edges and it differs from the background mesh only in the neighbourhood of the shock front, thus keeping the overall number of grid points to a minimum. It is also worth to underline that the use of a conservative gasdynamic solver makes available an hybrid computational technique in which some shocks are fitted whereas the remaining ones are captured by the gasdynamic solver. This has been demonstrated through the computation of a steady Mach reflection in which the incident shock is captured whereas the reflected shock and the Mach stem are fitted. It is our opinion that the original features proposed within the present shock-fitting approach could considerably ease the extension of this technique to multiple-shock cases. This is because in an unstructured framework the background mesh can be split and remeshed as much as one likes without encountering any major difficulty in terms of mesh topology and grid generation. References [1] Von Neumann J, Richtmyer RD. A method for the numerical calculation of hydrodynamic shocks. J Appl Phys 1950;21:232–7. [2] Moretti G, Abbett M. A time-dependent computational method for blunt body flows. AIAA J 1966;4:2136–41. [3] Moretti G. Computation of flows with shocks. Annu Rev Fluid Mech 1987;19:313–7. [4] Salas MD. Shock-fitting method for complicated two-dimensional supersonic flows. AIAA J 1976;14:583–8. [5] Yamamoto O, Anderson DA, Salas MD. Numerical calculations of complex Mach reflection. In: AIAA fluid dynamics, plasma dynamics, and lasers conference, Snowmass, CO, AIAA-1984-1679; 1984. [6] Marconi F, Salas M. Computation of three-dimensional flows about aircraft configurations. Comput Fluids 1973;1:185–95. [7] Yamamoto O, Anderson DA, Salas M. Numerical calculations of complex Mach reflection, AIAA Paper 84-1679; 1984. [8] Hartwich PM. Fresh look at floating shock fitting. AIAA J 1990;29:1084–91. [9] Dadone A, Fortunato B. Three-dimensional flow computations with shock fitting. Comput Fluids 1994;23:539–50. [10] Nasuti F, Onofri M. Analysis of unsteady supersonic viscous flows by a shockfitting technique. AIAA J 1996;34:1428–34. [11] Nasuti F, Onofri M. Viscous and inviscid vortex generation during startup of rocket nozzles. AIAA J 1998;36:809–15. [12] Lee TK, Zhong X. Spurious numerical oscillations in simulation of supersonic flows using shock-capturing schemes. AIAA J 1999;37(3): 313–9.

726

R. Paciorri, A. Bonfiglioli / Computers & Fluids 38 (2009) 715–726

[13] Carpenter MH, Casper JH. Accuracy of shock capturing in two spatial dimensions. AIAA J 1999;37(9):1072–9. [14] Casper J, Carpenter MH. Computational considerations for the simulation of shock-induced sound. SIAM J Sci Comput 1998;19(3):813–28. Available from: http://link.aip.org/link/?SCE/19/813/1. [15] Roy CJ. Grid convergence error analysis for mixed-order numerical schemes. AIAA J 2003;41(4):595–604. [16] Roy CJ, McWherter-Payne MA, Oberkampf WL. Verification and validation for laminar hypersonic flow fields. Part 1: verification. AIAA J 2003;41(10): 1934–43. [17] Suresh A. Interaction of a shock with a density disturbance via shock fitting. J Comput Phys 2005;206(1):6–15. [18] Robinet J-C, Gressier J, Casalis G, Moschetta J-M. Shock wave instability and the carbuncle phenomenon: same intrinsic origin? J Fluid Mech 2000;417:237–63. [19] Pandolfi M, D’Ambrosio D. Numerical instabilities in upwind methods: analysis and cures for the carbuncle phenomenon. J Comput Phys 2001;166(2):271–301. [20] Gnoffo P. Computational fluid dynamic technology for hypersonic applications. In: AIAA/ICAS international air and space symposium and exposition: the next 100 years, Dayton, OH, July 14–17; 2003. [21] Gnoffo P, White JA. Computational Aerothermodynamic simulation issues on unstructured grids, AIAA-CP-2004-2371; 2004. [22] Gnoffo P. Simulation of stagnation region heating in hypersonic flow on tetrahedral grids. In: 18th AIAA computational fluid dynamics conference, 25– 28 June, Miami, FL, Paper 2007-3960; 2007. [23] Luke E. On robust and accurate arbitrary polytope CFD solvers. In: 18th AIAA computational fluid dynamics conference, 25–28 June, Miami, FL, Paper 20073956; 2007. [24] Candler GV, Barnhardt MD, Drayna TW, Nompelis I, Peterson DM, Subareddy P. Unstructured grid approaches for accurate aeroheating simulations. In: 18th AIAA computational fluid dynamics conference, 25–28 June, Miami, FL, Paper 2007-3959; 2007. [25] Unstructured grid methods for advection dominated flows, Tech. Rep. AGARDR-787; 1992. [26] Moretti G. The k-scheme. Comput Fluids 1979;7:191–205. [27] Moretti G. Thirty-six years of shock fitting. Comput Fluids 2002;31(4– 7):719–23. [28] Moretti G. Efficient Euler solver with many applications. AIAA J 1988;26:655–60. [29] Kopriva DA, Zang TA, Hussaini MY. Spectral methods for the Euler equations: the blunt body problem revised. AIAA J 1991;29:1458–62. [30] Paciorri R, Favini B, Sabetta F. Fast Euler solver for non-equilibrium steady flows. Comput Fluid Dyn J 1994;3(3):333–42. [31] Zhong X. High-order finite-difference schemes for numerical simulation of hypersonic boundary-layer transition. J Comput Phys 1998;144:662–709. [32] Brooks GP, Powers JM. Standardized pseudo-spectral formulation of the inviscid supersonic blunt body problem. J Comput Phys 2004;197:58–85. [33] Rosendale JV. Floating shock fitting via lagrangian adaptive meshes, Tech. Rep., Institute for Computer Applications in Science and Engineering, NASA Langley

[34]

[35] [36] [37] [38]

[39]

[40] [41] [42]

[43]

[44]

[45]

[46]

[47]

[48]

[49]

[50]

Research Center, Hampton, VA 23681-0001, NASA Contractor Report 194997; ICASE Report No. 94-89; 1994. Parpia IH, Parikh P. A solution-adaptive mesh generation method with cellface orientation control. In: AIAA, aerospace sciences meeting and exhibit, 32nd, 10–13 January, Reno, NV, No. AIAA-1994-416; 1994. Trepanier J-Y, Paraschivoiu M, Reggio M, Camarero R. A conservative shock fitting method on unstructured grids. J Comput Phys 1996;126:421–33. Gloth O, Hänel D, Tran L, Vilsmeier R. A front tracking method on unstructured grids. Comput Fluids 2003;32:547–70. Roe PL. Approximate Riemann solvers, parameter vectors and difference schemes. J Comput Phys 1981;43:357–72. Mavriplis DJ. Unstructured mesh discretizations and solvers for computational aerodynamics. In: 18th AIAA computational fluid dynamics conference, 25–28 June, Miami, FL, Paper 2007-3955; 2007. Ivanov M, Khotyanovsky D, Paciorri R, Nasuti F, Bonfiglioli A. Numerical simulation of weak steady shock reflections. In: Proceedings of the 26th international symposium on shock waves, Göttingen, Germany; 2007. Müller J-D, Roe P, Deconinck H. A frontal approach for internal node generation in Delaunay triangulations. Int J Numer Meth Fluids 1993;17:241–56. Müller JD. DELAUNDO. Available from http://www.cerfacs.fr/~muller/ delaundo.html. Shewchuk JR. Triangle: engineering a 2D quality mesh generator and Delaunay triangulator. In: Lin MC, Manocha D, editors. Applied computational geometry: towards geometric engineering. Lecture notes in computer science, vol. 1148. Berlin: Springer; 1996. p. 203–22. from the first ACM workshop on applied computational geometry. Bonfiglioli A. Fluctuation splitting schemes for the compressible and incompressible Euler and Navier–Stokes equations. Int J Comput Fluid Dyn 2000;14:21–39. van Leer B, Lee W-T, Roe P. Characteristic time-stepping or local preconditioning of the Euler equations. In: AIAA 10th computational fluid dynamics conference, AIAA-91-1552-CP; 1991. Mesaros L, Roe P. Multidimensional fluctuation splitting schemes based on decomposition methods. In: 12th AIAA CFD conference, San Diego, Paper 951699; 1995. Paillère H, Deconinck H, Roe P. Conservative upwind residual-distribution schemes based on the steady characteristics of the Euler equations. In: 12th AIAA CFD conference, San Diego, Paper 95-1700; 1995. Deconinck H, Paillère H, Struijs R, Roe P. Multidimensional upwind schemes based on fluctuation-splitting for systems of conservation laws. Comput Mech 1993;11(5/6):323–40. van der Weide E, Deconinck H, Issman E, Degrez G. A parallel, implicit, multidimensional upwind, residual distribution method for the Navier–Stokes equations on unstructured grids. Comput Mech 1999;23:199–208. Deconinck H, Roe P, Struijs R. A multi-dimensional generalization of Roe’s flux difference splitter for the Euler equations. Comput Fluids 1993;22(2/ 3):215–22. Lyubimov AN, Rusanov VV. Gas flows past blunt bodies, part ii: table of the gasdynamic functions, NASA TT F-715; 1973.