Abstract

Within a continuum framework, flows featuring shock waves can be modelled by means of either shock capturing or shock fitting. Shock-capturing codes are algorithmically simple, but are plagued by a number of numerical troubles, particularly evident when shocks are strong and the grids unstructured. On the other hand, shock-fitting algorithms on structured grids allow to accurately compute solutions on coarse meshes, but tend to be algorithmically complex. We show how recent advances in computational mesh generation allow to relieve some of the difficulties encountered by shock capturing and contribute towards making shock fitting on unstructured meshes a versatile technique.

1. Introduction

A shock wave is a surface in a flow-field across which various thermodynamic and kinematic properties of the fluid appear, from a macroscopic viewpoint, to change discontinuously. In reality, shock fronts have a finite thickness, of the order of a few mean free paths, across which the physical properties change continuously. Since the shock thickness is much smaller than any macroscopic length scale, such as the radius of curvature of a curved shock, the assumption that the discontinuity thickness is negligible compared with any macroscopic dimension is a fundamental one in most engineering analyses of flows of practical interest.

From the mathematical modeling viewpoint, although it might be expected that a continuum approach based upon the Navier-Stokes equations fails to be valid at the microscopic length scale that characterises the width of a shock wave, Navier-Stokes solutions of the shock layer (in this context, by the term “shock layer’’ we refer to the microscopic region across which the fluid properties change from their upstream to their downstream values) agree reasonably well with experimental data, at least when the shocks are not too strong (in this paper, except where otherwise noted, when we mention the strength of a shock wave we refer to the magnitude of the jump in fluid properties that occurs through the shock; this might be defined in several ways, for example, by the ratio between the downstream and upstream pressures). Moreover, since changes across a shock wave occur only in the shock-normal direction, the problem is of one-dimensional nature so that closed-form or semianalytical solutions of the 1D Navier-Stokes equations can be obtained and provide the distribution of the various flow properties within the shock layer, see, for example, the book by Hayes [1] or the more recent one by Salas [2, Chapter 3.7].

In the vanishing viscosity limit, the shock width shrinks to zero, the fluid properties change discontinuously through the shock, which becomes a true discontinuity (in the mathematical sense), and the Rankine-Hugoniot jump relations (algebraic equations that connect the upstream and downstream flow properties with the shock speed) are obtained from the 1D Navier-Stokes equations (see [3, page 136] and also the related discussion in [2, Chapter 3.7.4] for details).

Shifting from the theoretical to the numerical models actually used in the computation of flows with shocks, the dualism between viscous shocks, having a finite though microscopic thickness, and inviscid shocks, being a true discontinuity, shows up again in the two numerical approaches: shock capturing and shock fitting, that have been traditionally used in Computational Fluid Dynamics (CFD) to simulate shock waves.

Since shock thicknesses are of the order of a few mean free paths and the length of a mean free path is of the order of 107 m for air at standard conditions, it is easily seen that fully resolved Navier-Stokes simulations of a shock wave would require a large number of grid points within a microscopic region surrounding the shock front, in addition to those grid points that are needed to properly resolve the macroscopic features of the flow field. Although the tremendous advances in computer power “may bring computational fluid dynamics to a molecular scale, where shocks will cease to be discontinuities and have a real thickness, full of viscous and thermodynamical effects …’’ [4], even nowadays, the inner structure of the shock waves is necessarily underresolved in a realistic Navier-Stokes simulation. This is precisely the situation that occurs when shock-capturing schemes, which represent the numerical model most commonly used nowadays, are used to solve flows with shocks. Shock-capturing discretizations, which date back to the work by von Neumann and Richtmyer [5], lay their foundations in the mathematical theory of weak solutions, which allows to compute all type of flows, including those with shocks, using the same discretization of the conservation law form of the governing equations at all grid cells. This has obvious consequences in terms of coding simplicity, since the same set of operations is repeated within all control volumes of the mesh, no matter how complicated the flow might be. Coding simplicity is the primary reason for the widespread use of shock-capturing schemes in the simulation of compressible flows. Coding simplicity comes not for free, however, and shock capturing solutions of flows featuring strong shock waves are sometimes characterised by large numerical errors and by the appearance of bizarre anomalies. The sources of the numerical errors are various. The most evident cause is the numerical thickness of the captured shock, which depends on the local mesh size and, for the reasons explained above, is always much larger than the physical shock thickness. Two other possible causes are the generation of spurious numerical oscillations [68] along the shock front and the reduction of the order of accuracy of the computed solution to a mere first order in the entire region downstream of a captured shock [6, 911]. Spurious oscillations arise when the faces of the control volumes of the mesh are not aligned with the shock front; loss of accuracy is probably due to the fact that the flow state of those gridpoints that fall inside a captured shock bear little resemblance with the physical shock structure and are rather to be seen as numerical artefacts. Among the various anomalies exhibited by shock-capturing discretizations, the most evident ones are the stagnation point anomalies [1216] and the carbuncle phenomenon [1719]. A large body of the literature has emerged over the last decades addressing possible cures to the aforementioned problems. Trying to give even a brief account of these attempts is well beyond the scope of the present paper, but it seems clear to the authors that none of these efforts has yet given birth to the “ultimate’’ shock-capturing scheme. Given up hope of finding such an ideal scheme, another chance to improve the quality of shock-capturing calculations consists in coupling the flow solver with a posteriori mesh generation/adaptation/refinement procedures. Indeed, the generation of a new mesh, tailored on the shock pattern and finer in the shock-normal direction, or the refinement/adaptation of an existing mesh not only allows to bring the captured shock thickness closer to (but still far from) its physical value, but also allows to improve the solution quality by aligning the faces of the control volumes with the shock front. Of course, improving the mesh quality alleviates the symptoms, but does not cure the disease. Nonetheless, it is true that good meshes help considerably in reducing the numerical errors induced by captured shock waves, thus significantly enhancing the quality of the computed solution in the entire region downstream of a shock. A discussion on the possible strategies to construct meshes that are better suited to capture shock waves is one of the goals of the present paper.

As already mentioned, however, a different approach for modeling shock waves within the continuum framework is possible: it is called shock fitting and, from a historical perspective, it is even older than shock capturing, since its first appearance is credited [2, Chapter 1], [20] to Emmons [21, 22]. Shock fitting consists in locating and tracking the motion of the discontinuities which are treated as boundaries between regions where a smooth solution to the governing PDEs exists. The Rankine-Hugoniot jump relations are used to calculate the space-time evolution of the discontinuity and the governing equations are discretized only within the smooth regions of the flow field. Shock-fitting methods have enjoyed a remarkable popularity at the dawn of CFD thanks to their ability to accurately compute flows featuring strong shocks on coarse grids using the modest computational resources available at that time. For a given mesh resolution, shock-fitting discretizations are characterized by smaller numerical errors than shock-capturing solutions because they are immune from all those sources of error, previously described, that arise along a captured shock front. Moreover, shock-fitting schemes do not suffer from the other numerical anomalies described previously.

Nonetheless, the exponential growth of computing power, along with the algorithmic difficulties that hindered the development of general-purpose codes based on the shock-fitting approach in the framework of structured meshes, decreed the gradual dismissal of the technique over the last decades. Simply stated, the slogan “one code for all flows’’ that has so much contributed to the popularity of shock-capturing discretizations does not carry over to shock fitting. Indeed, the two variations of the basic technique that have been so far proposed in the structured grid setting, that is, “boundary’’ shock fitting and “floating’’ shock fitting, are either algorithmically simple, but limited to topologically simple configurations or more widely applicable, but algorithmically complex. Recently, a new shock-fitting technique for unstructured two- and three-dimensional meshes has been proposed and developed [2327] by two of the authors of this article. This new technique reduces significantly the algorithmic difficulties encountered in the previous implementations of the shock-fitting technique within the structured-grid context. In this newly developed technique, the fitted shock fronts are moving boundaries that are free to float inside a computational domain discretized using triangles, in two space dimensions, or tetrahedra in three. This is achieved by locally regenerating the computational mesh in the neighbourhood of the fitted shocks while ensuring that the shock lines (or the shock surfaces in three dimensions) are always part of the mesh during their motion. The development of this new shock-fitting technique has been greatly furthered by the availability of public domain software libraries, developed over the last years, implementing a variety of tools for the generation of unstructured meshes.

This brief outline of the two types of numerical techniques used in the computation of flows featuring strong shocks highlights the role played by mesh generation/adaptation/refinement techniques in providing high-quality solutions. On the one hand, the use of these techniques improves the quality of shock-capturing solutions; on the other hand, the development of modern unstructured mesh generation tools is opening up new perspectives in the development of shock fitting techniques. To support this assertion, three different applications of mesh generation/adaptation/refinement techniques in the context of shock-capturing computations will be discussed and analysed in Section 2. Section 3 is devoted to shock-fitting: a brief review of the traditional “boundary’’ and “floating’’ techniques is given first, followed by the newly developed technique for unstructured meshes. Applications of the unstructured shock-fitting methodology both in two and three space dimensions is illustrated in Sections 3.1 and 3.2.

2. Shock Capturing

In this section we shall describe three different mesh adaptation and/or refinement techniques that can help improving the resolution of strong shock waves computed by means of shock-capturing discretizations. The results to be presented in the following paragraphs have been obtained using different codes, both in house and commercial; the discretization techniques used in the various codes will not be described here, but the reader will be referred to the relevant references for details.

2.1. Structured Mesh Adaptation

The first technique to be described consists in an in-house developed, mesh adaptation code, applied within a structured multiblock context. The selected test-case consists in the hypersonic flow past the ESA IXV capsule [28, 29], a lifting reentry body, see Figure 1(a), equipped with two variable deflection flaps (𝛿𝐿, 𝛿𝑅) that allow an aerodynamically controlled reentry. The free-stream flow conditions considered in the present case are those typical of a high speed and high altitude point along the reentry trajectory of the space capsule and they are summarised in Table 1. For the given flight conditions, the flow past the capsule is completely laminar and chemically reacting.

The numerical solutions have been computed by means of the commercial code CFD++, developed by Metacomp [30]. The Navier-Stokes equations are supplemented with the 5-species nonequilibrium model by Park [31] and a diffusion model based upon the Gupta-Yos model [32]. The boundary conditions at the solid walls employ the fully catalytic model and radiative equilibrium.

As mentioned previously, the most important source of numerical errors in the shock-capturing calculation of an hypersonic flow is the shock wave itself, which affects the solution quality of the entire flow-field downstream of the discontinuity. An effective mesh adaptation procedure must therefore be tailored for the shock.

Figure 1(b) shows the computational domain around half of the capsule. The computational domain has been split in ten blocks, whose topology has been chosen so as to allow an easy management of the cell clustering in the critical regions: the windward side, the bow shock, and the flaps.

Starting from this block-topology, a baseline mesh is generated using a commercial mesh generator. This mesh is characterised by cell clustering near the wall in order to capture the boundary layers, but no attempt has been made to refine the mesh within the shock regions. Figure 2(a) shows some details of the ten-block mesh; the pressure field computed on the this mesh is displayed in Figure 2(b). The baseline mesh has been modified by changing the spatial distribution of the mesh nodes belonging to the coordinate lines that are perpendicular to the body surface, on the basis of the pressure gradient distribution.

More specifically, the sequence of N mesh nodes (𝑥baseline𝑖) belonging to each coordinate line is replaced by a new sequence of N nodes (𝑥adapted𝑖), whose mesh spacing Δ𝑠𝑖=𝑥𝑖+1𝑥𝑖 is computed as follows:Δ𝑠adapted𝑖=1𝑤adapted𝑖𝑖=1,𝑁1,(1) where 𝑤adapted𝑖 is the weighting function that drives the cell clustering. Using (1), the mesh spacing will be finer in regions where the corresponding weighting function is large and vice versa. The weighting function 𝑤adapted𝑖 associated with the adapted mesh, defined as1𝑤adapted𝑖=𝑐𝑤baseline𝑖+(1𝑐)𝑤𝑝𝑖(2) is a weighted average of the corresponding weighting function, 𝑤baseline𝑖, which is computed from the known mesh spacing of the baseline mesh:𝑤baseline𝑖=1Δ𝑠baseline𝑖𝑖=1,𝑁1,(3) and a weighting function 𝑤𝑝𝑖 based upon the pressure gradient distribution computed on the baseline mesh and normalised in such a way that:𝑁1𝑖=11𝑤𝑝𝑖=𝐿,(4) where 𝐿=𝑁1𝑖=1Δ𝑠baseline𝑖. The constant 𝑐 in (2) is set equal to 0.5.

The comparison between Figures 2(a) and 2(c) shows how the baseline mesh has been modified using the aforementioned pressure-driven adaptation strategy. Note, in particular, that the procedure has been activated only in those blocks where strong shocks are present and it has the effect of clustering the grid points around the shock front and aligning the cell faces with the shock front. Moreover, the procedure does not significantly alter the cell distribution of the baseline mesh inside the boundary layer, because (2) combines the solution-dependent weighting function 𝑤𝑝𝑖 with 𝑤baseline𝑖 that keeps track of the original distribution of mesh points within the near-wall region. As a result, only grid points located near the shocks, either on their upstream or downstream side, are moved towards the shock fronts.

The positive effects of the mesh adaptation procedure on the quality of the numerical solution can be observed by comparing Figures 2(b) and 2(d). In particular, the shock thickness computed on the adapted mesh is thinner than the one computed on the baseline mesh and the shock-shock interaction more clearly outlined. Moreover, the beneficial effects of the mesh adaptation procedure are not limited to the shock region but can also be observed on the surface of the capsule. Figure 3 compares the pressure distribution along the windward side of the vehicle obtained on the baseline mesh with that computed on the adapted mesh. The positive effects of mesh adaptation around the shocks are highlighted by the significant increase of the pressure levels that can be observed both in the nose region and on the surface of the flaps. Improving the prediction of the pressure distribution on the capsule’s surface significantly enhances the reliability of the numerical estimate of the aerodynamic coefficients; this is particularly true for the moment coefficients that strongly depend upon the pressure distribution on the flaps.

A more comprehensive report of the results obtained by the application of this mesh adaptation technique to the computation of the flow field past the IXV vehicle can be found in [33].

2.2. Structured Mesh Refinement/Adaptation Using Overset Meshes

In this section we describe a more innovative mesh refinement/adaptation technique, suited for structured, multiblock, shock-capturing solvers, which is based upon the so-called “Chimera’’ or overset mesh approach.

Some early examples of application of the overset grid technology can be found in the pioneering papers by Volkov [34, 35], whereas some more recent applications are reported, for instance, in [36, 37]. This approach is extremely useful for the simulation of flows in domains with complex boundaries (as it happens to be the case in most engineering problems of practical interest), because it allows to mesh complex geometries with relative ease while retaining most of the computational advantages of algorithms based on block structured grids: simple data structure, easy application of iterative schemes based on geometrical multigrid, approximate factorisation techniques for implicit time integration, and so on. The flexibility and convenience of the Chimera approach lies in the process of grid generation, that can be carried out by producing individual pieces of grid around each portion of the frontier; the various blocks are then assembled in a global grid, letting them overlap. Of course, proper algorithms must be developed to compute grid connections at the block boundaries and possible “hole’’ cutting in the inner portion of each block.

Besides making grid generation more flexible, the overset grid capability can be exploited also for mesh refinement in regions where the solution is known to exhibit strong gradients, like wakes, strong expansions, and shocks. This capability will be illustrated with the following example, which consists in the calculation of the supersonic flow (free-stream Mach number equal to 3) past the space launcher VEGA, shown in Figure 4, during its ascent trajectory. The model equations are the Reynolds Averaged Navier-Stokes equations; the one-equation Spalart-Allmaras model [38] was used for turbulence closure. The CFD solver is an in-house code [39] that exploits the overset grid capability. The discretization is based on a finite volume scheme with second-order accurate ENO [40] approximation of the inviscid fluxes and centred approximation of the viscous terms.

The simulation was aimed at assessing the thermal flux on critical sections of the launcher surface (raceways, retro-rockets, junctions, etc.) and of the appendages (such as antennas, cameras, etc.); therefore, any single device protruding off the vehicle’s surface has been carefully discretized, as shown in the detailed view of Figure 4(b). In addition to the geometrical description of the finest details of the launcher, overset meshes have also been exploited for local refinement in critical regions of the flow field.

A preliminary computation on a grid which is not refined in the interior of the domain provided the approximate location of the areas featuring strong gradients. Then, on the basis of this preliminary solution, new blocks were designed with a commercial CAD code and then discretized using a commercial grid generator. The overlapping blocks used for mesh refinement in various critical flow regions can be clearly seen in the left half of both frames of Figure 5. The right half of the two frames of Figure 5 shows the density contours computed on the adapted mesh. It can be seen that the overlapping grid technology can be an effective tool to improve the resolution of specific flow features, such as the bow shock shown in the right frame of Figure 5. The quality of the aforementioned mesh refinement technique might be further improved, particularly for shock waves, if coupled with a shock-detection algorithm capable of providing a reasonable approximation to the shape of the shock surface. The overlapping blocks might then be built with one set of coordinate lines perpendicular to the guessed shock surface. This kind of semiautomatic procedure for constructing approximate shock shapes has been used in conjunction with shock fitting and shall be described in Section  3.2.

2.3. Unstructured Grid: Anisotropic Mesh Adaptation

The third methodology to be described consists in the application of an Anisotropic Mesh Adaptation (AMA) technique to unstructured, triangular meshes. AMA is more suitable than the techniques discussed in the previous sections to simulate flows characterised by complex shock topologies. Indeed, unstructured grids clearly have the potential of making mesh adaptation much simpler than it is in the context of structured grids and allow to easily refine the computational mesh even when multiple shock interactions occur. In this section we report the authors’ experience in using angener [41], a public domain anisotropic mesh generator. The AMA technique available in angener consists in building a mesh made of isotropic, almost equilateral triangles in a computational plane which are then transformed into the physical plane using a mapping governed by the Hessian of an appropriate scalar quantity, available from a previous calculation made on a different grid. In the calculations presented herein we have chosen the first component of the set of dependent variables: 𝜌, which is adequate to resolve both shocks and contact discontinuities. Further details concerning the AMA technique implemented in the angener code can be found in [42]. The flow solver is eulfs [43, 44] an in-house, shock-capturing, vertex-centred code which uses second-order accurate, multidimensional upwind, Fluctuation Splitting schemes for the spatial discretization.

To assess the performances of the AMA technique, an Eulerian flow characterised by multiple interactions has been computed. The selected flow configuration is a type IV shock-shock interaction [45], as classified by Edney [46]. This type of interaction arises when a weak (with the term “weak’’ we here mean an oblique shock featuring a shock angle 𝜎 which is smaller than that corresponding to the maximum flow deflection for the given free-stream Mach number) oblique shock meets the bow shock produced by a supersonic stream impinging on a blunt nosed body. Figure 6 shows a sketch of the geometry and flow structure: in the present case the blunt body is a circular cylinder with radius 𝑅 and the origin of the reference frame is located in the centre of the circle. The free-stream Mach number is 𝑀=10 and the oblique shock, which forms an angle 𝜎=10 with respect to the undisturbed flow, crosses the 𝑥-axis at 𝑥𝑃/𝑅=1.667. A first triple point forms where the oblique shock impinges on the bow shock. A reflected shock and a contact discontinuity arise in this triple point and move downstream. The contact discontinuity separates the supersonic stream which has crossed the reflected oblique shock from the subsonic flow downstream of the bow shock. The reflected shock coming from the first triple point rejoins the bow shock in a second triple point where a new reflected shock and a contact discontinuity arise. The two contact discontinuities bound a supersonic jet which is directed towards the body surface, whereas the second reflected shock triggers a series of wave reflections inside the supersonic jet bounded by the two contact discontinuities. Figure 7 shows the density contours computed on an isotropic mesh which has been constructed using delaundo [47], a public domain frontal/Delaunay mesh generator [48]. Starting from the solution displayed in Figure 7(a), an adapted mesh has been constructed using angener and a solution has been computed on this new mesh using the eulfs code. The AMA/solution process has been repeated three times and the third mesh and the corresponding solution are displayed in Figure 7(b). The anisotropically adapted mesh features roughly twice as many nodes and triangles than the original mesh; notice that mesh refinement in the most “active’’ areas of the flow field is partially balanced by mesh coarsening in areas of smooth flow variations. It is clear that AMA allows to improve considerably the resolution of the interaction region; this allows not only to recognise the various flow features, which are indistinguishable in the solution computed on the uniform grid due to the large numerical thickness of the captured discontinuities in that region, but also to reduce the disturbances (even if these are not completely removed) that are generated along the bow shock: compare the density isolines computed within the shock layer (here and in the remainder of this paper, by the term “shock layer’’ we refer to the flow region bounded between the bow shock and the solid body), which are shown in the leftmost frames of Figures 7(a) and 7(b).

3. Shock Fitting

Nowadays, shock-fitting techniques are neither widely used nor even known to most CFD practitioners, so that it might be worth giving an historical background about the traditional shock-fitting techniques used in the structured grid setting and briefly explaining its new developments within the unstructured grid context.

The first shock-fitting calculation is credited to Emmons [21, 22], but the undisputed guru of this technique is certainly Moretti who, along with his students, has greatly contributed to the popularity of the technique. The recent book by Salas [2] describes the theoretical foundations of shock-fitting along with its practical implementation. Shock fitting, which until recently has only been used in the structured-grid setting, has historically evolved following two different algorithmic approaches: boundary shock fitting and floating shock fitting. In the boundary shock-fitting approach, the shock is made to coincide with one of the boundaries of the single- or multiblock mesh which is used to discretize the flow field. For example, Figure 8 describes the steady calculation of a bow shock about the fore-body a blunt nosed object: the fitted bow shock is the leftmost boundary of the computational mesh. Starting from the initial configuration shown in Figure 8(a), the grid is deformed, see Figure 8(b), due to the motion of the fitted shock front, until it settles to its final position, shown in Figure 8(c), which corresponds to vanishing shock speed. Figure 8(d) shows the computed flow field. Boundary shock fitting is an algorithmically simple technique, since the treatment of the algebraic relations that hold across the shock (the Rankine-Hugoniot relations) is confined to the boundary points and therefore no modification is required to the computational kernel of the code which is used to discretize the governing PDEs within the smooth regions of the flow field. This feature greatly simplifies the coding, but clearly makes the treatment of embedded and/or interacting shocks a “hard bone to chew’’ [49].

A step forward towards a technique capable of handling more complex flow configurations, as well as unsteady flows, was undertaken by Moretti with the development of the so-called floating shock fitting, see, for example, [50, 51]. In the floating version of the shock-fitting technique, the discontinuities are a collection of shock points, shown by open circles in Figure 9(a), which are allowed to move (float) freely over a fixed background structured grid. The shock front is described by its intersections with grid lines, which gives rise to 𝑥 and 𝑦 shock points, as sketched in Figure 9(b). The main features of the methodology are the same as those of the boundary fitting technique, except for the need of a special treatment for grid nodes neighbouring shocks. This is because approximating derivatives by differences between nodes located on opposite sides of a shock must be avoided and therefore, ad hoc finite difference formulae (see, e.g., [52, 53]) have to be used in this case. Moreover, since the 𝑥 and 𝑦 shock points are constrained to be located along the grid lines, see Figure 9(b), the technique is able to adequately describe isolated shocks, but it does not allow the precise treatment of the shock-shock interaction point, since the latter may be located in between grid-lines, as also sketched in Figure 9(b). The presence of an interaction point needs then to be modelled by appropriately modifying the calculation of the neighbouring shock points, see, for example, [54]. Floating shock-fitting codes have been used with success in the past to compute steady and unsteady two-dimensional flows involving shock reflections and shock interactions, see, for example, [50, 52, 55]. Very recently, the floating shock-fitting (denoted as front-tracking) technique has been reproposed by Rawat and Zhong [53] in the framework of high-order, structured-grid schemes. Nonetheless, floating shock fitting is an algorithmically complex technique, primarily because of the need to interface the inherently unstructured collection of points used to represent the shock fronts with the underlying structured mesh, which is described using the classical IJK data structure. A thorough account of these difficulties can be found in [53]. Taking advantage of the availability of public domain software libraries [5660] that accomplish unstructured, volumetric, and surface mesh generation, a novel unstructured shock-fitting technique that combines features of both the boundary and floating shock-fitting techniques has been recently developed by two of the authors [2327]. In this novel unstructured shock-fitting approach the shock front is discretized as a double-sided, triangulated surface (a polygonal curve in 2D) which is treated as a zero-thickness, internal boundary by the shock-capturing eulfs [43, 44] solver (whose main features have already been described in Section 2.3) which is used to discretize the governing PDEs within the computational domain. The two sides of the triangulated shock surface correspond to the upstream and downstream sides of the discontinuity. A set of dependent variables is assigned to each of the grid points which belong to both sides of the discontinuity. The local shock speed and nodal values on the high pressure side of the shock are computed by enforcing the Rankine-Hugoniot relations across the discontinuity. The shock is allowed to move over and independently of a background tetrahedral (triangular in 2D) grid which is locally adapted at each time step to ensure that the nodes and the triangles (edges in 2D) that make up the shock front are part of a constrained Delaunay tetrahedralization (triangulation in 2D). The coupling of the shock-fitting technique with a shock-capturing solver also enables a hybrid mode of operation in which some of the shocks are fitted, whereas all others are captured. The reader is referred to [2327] for a detailed description of the two-, respectively three-dimensional, version of this novel unstructured shock-fitting approach; here we report some examples demonstrating how the advances in unstructured mesh generation enabled us to bring shock fitting back to life.

3.1. Two-Dimensional, Unstructured Shock Fitting

The current capabilities of the two-dimensional version of the unstructured shock-fitting technique will be demonstrated by referring to the same test case already examined in Section 2.3. In two dimensions this technique is able not only to fit the shock waves, but also the contact discontinuities and the interaction points where various discontinuities meet. The shock fitting simulation uses the unadapted mesh shown in Figure 7(a) as the background triangulation. The fitted discontinuities are approximated by a set of linear splines (shown by the heavier solid lines in Figure 10) and their initial location is guessed using a preliminary shock-capturing calculation performed on the background grid. Starting from this initial solution, the shock-fitting calculation is advanced in time until steady state is reached: at each time step, the background triangulation is locally modified to accommodate the fitted discontinuities which are free to move at a speed which obeys the Ranking-Hugoniot jump relations. Local remeshing along the discontinuities requires a constrained Delaunay triangulation (CDT), a task which is carried out using the triangle code [56], a public domain software developed by Shewchuk [56] and based on Ruppert’s algorithm [61]. At steady state, the speed of the discontinuities vanishes and the shock-fitting grid does not change any longer. Note that the converged shock-fitting grid, Figure 10, and the background one, Figure 7(a), only differ in the neighbourhood of the fitted discontinuities. Figure 7(a) shows the density isocontours computed by shock capturing on the same grid also used as background triangulation in the shock-fitting simulation; the shock-fitting result is shown in Figure 10. The superiority of the shock-fitting solution over a shock-capturing solution computed on an unadapted mesh of comparable resolution is evident: the detailed view of the region of interaction, shown in the smaller frame of Figure 10, shows that the various discontinuities are clearly identifiable in the shock-fitting solution, whereas their features are completely lost in the shock-capturing calculation, shown in Figure 7(a). If we now compare the shock-fitting solution of Figure 10 with the shock-capturing solution on the anisotropically adapted mesh, shown in Figure 7(b), both calculations appear to be of comparable quality. In particular, both solutions feature smooth contours within the shock layer, even though small wavelength disturbances are visible in the proximity of the bow shock in the shock capturing solution. When comparing shock fitting with shock capturing computed on the anisotropically adapted mesh, it is worth underlining that the former technique requires roughly half the number of grid points and triangles than the latter. A more thorough comparison between shock fitting and shock capturing on anisotropically adapted meshes can be found in [24].

3.2. Three-Dimensional, Unstructured Shock Fitting

Although the three-dimensional version of the unstructured, shock-fitting methodology is conceptually the same as its two-dimensional version described in Section 3.1, the implementation of these ideas in 3D is very challenging and far from trivial. A three-dimensional shock-fitting calculation proceeds similarly to the two-dimensional case. A preliminary shock-capturing calculation on the background mesh provides the initial solution and the approximate shape and location of the fitted shock surface(s). The generation of the initial shock surface(s) proceeds as follows: a cloud of candidate shock points is extracted from the shock-capturing calculation using the algorithm by Ma et al. [62]; these shock points are then processed using various tools available in the meshlab [58] software and a preliminary, triangulated shock surface is obtained. The shock-fitting calculation is then advanced in time until steady state is reached. During its time evolution, the fitted shock surface is allowed to float over the background tetrahedral mesh using the shock speed computed within each of the grid points that belong to the shock surface. Volumetric remeshing in the neighbourhood of the shock surface requires a constrained Delaunay tetrahedralization (CDT). This complex operation is performed at each time step using TetGen, a public domain code recently developed by Si [57]. Differently from the two-dimensional case, the solution of a CDT problem in three dimensions may require the insertion of several additional mesh points (so-called Steiner points), a situation that needs to be properly handled by the shock-fitting algorithm. Moreover, since the shock surface may change its area and shape during its motion, it has to be periodically remeshed. Surface remeshing is required also to ensure that the local mesh size of the shock surface triangulation is comparable to that of the background tetrahedral mesh; the fulfilment of this criterion avoids the appearance of badly shaped tetrahedra in the neighbourhood of the shock front. Unfortunately, most of the available surface-meshing tools have been developed for computer graphics applications and are typically designed to optimise the distribution of triangles over the surface based upon geometrical features (e.g., the local curvature), rather than user-specified criteria. The solution to our specific surface meshing problem has therefore required a remarkable effort in the search and testing of appropriate software libraries. Among the various tools we have come across, Yams, a surface-meshing code developed by Frey [60], turned out to be the one that best suited our needs.

The three-dimensional version of the unstructured, shock-fitting technique has been applied to the calculation of hypersonic flows past geometrically simple bodies, such as the flow past a hemisphere, as well as rather complex ones, including the IXV vehicle described in Section 2.1. A detailed account of these efforts can be found in [26, 27]. The qualitative features of the three-dimensional version of the technique are here described by reference to the high speed flow past a blunt object which consists in a cylinder with a hemispherical nose and a conical flare forming an angle of 30 with respect to the cylinder’s axis. The free-stream conditions are Mach number 𝑀=4.04 and angle of attack equal to 20. This flow configuration has been studied both experimentally and numerically by Houtman et al. [63]; it is characterised by a bow shock and an embedded shock which originates at the cylinder-cone junction. Figure 11(a) shows an experimental visualisation of the flow field, which highlights the type VI shock-shock interaction, as classified by Edney [46], which occurs along the windward side of the body.

A preliminary shock-capturing calculation has been performed using a grid which also served as background grid for the shock-fitting calculations. Two different shock-fitting solutions have also been computed: in the first only the bow shock has been fitted, whereas in the second also the embedded shock has been fitted; Figure 11(b) shows both fitted shock surfaces at steady state. It is worth underlining that the current three-dimensional version of the shock-fitting algorithm is able to fit more than a single shock surface, but the fitted shock surfaces are not allowed to mutually interact. The cause of this limitation is not related to modeling issues (the two-dimensional interaction models described in [24] can be extended to the three-dimensional case without major problems), but it has to do with a specific mesh generation problem. Indeed, in the three-dimensional space, two shocks interact along a curve. This requires the capability to remesh different shock surfaces while prescribing the position of those shock points that are located along the curve where the shock surfaces meet. We are currently investigating this issue and we are confident that we shall be able to “fit’’ three-dimensional shock-shock interactions in the near future, as we already do in two space dimensions. For the time being, the interaction between different fitted shocks is handled by terminating the fitted surface of the embedded shock shortly before it reaches the bow shock and letting the shock-capturing code capture the shock-shock interaction. Although this is clearly sub-optimal than a fully fitted interaction, it however provides improved resolution compared to the situation in which only the bow shock is fitted and the embedded one is captured. This point will be supported with examples in the remainder of this section; for further details the reader is referred to [25, 27]. Figure 12, which shows the normalised pressure distribution on the body surface and on the symmetry plane for all three numerical solutions, allows to pinpoint some significant details of the flow structure predicted by the shock-fitting and shock-capturing simulations. The shock-capturing solution is characterised by a large shock thickness, which has nearly the same width as the shock layer thickness in the stagnation region. In the neighbourhood of the stagnation point, the shock-fitting solutions are seen to be better defined and predict a pressure peak which is much closer to the theoretical stagnation value and higher than that predicted by the shock-capturing calculation.

The shock-shock interaction causes an abrupt change of the shock slope which is clearly visible in Figure 12. In all solutions, the shock-shock interaction is spread over a region of finite width. This is because in all three calculations the shock-shock interaction is captured by the shock-capturing solver, even when both the bow shock and the embedded shock are fitted. The pressure distribution on the body surface is characterised by a second peak which occurs in the proximity of the shock-shock interaction. This is an area where the three numerical simulations show relevant mutual differences. The pressure level which is reached at this second pressure peak in the shock-capturing solution is visibly lower than those predicted by the two shock-fitting solutions. Moreover, downstream of this second pressure peak, the shock-capturing calculation features smoother contour lines that those predicted by the shock-fitting calculations, which is an indication of the fact that the shock-shock interaction has been excessively smeared in the shock-capturing calculation. The relevant role played by the shock-shock interaction is clearly highlighted by the fact that also the two shock-fitting solutions display nonnegligible differences in the neighbourhood of the second pressure peak. This observation also points to the fact that the capability of fitting shock-shock interactions in three dimensions would be an important add-on to the current methodology.

4. Conclusions

Building upon the authors’ experience in the field of the numerical simulation and of the development of numerical techniques, this paper highlights the strong relationship that exists between the development and application of mesh adaptation/refinement/generation techniques and the computation of high-quality solutions of flows featuring strong shocks.

The examples provided show that mesh adaptation/refinement/generation techniques that adapt the computational mesh to the shock front are a key ingredient to achieve accurate solutions for this kind of flow fields.

The application of these techniques in the framework of the shock-capturing approach allows to significantly reduce the numerical errors induced by the capture of the shock waves. Three different examples of mesh adaptation/refinement/generation techniques, used in conjunction with three different shock-capturing solvers, have been discussed and analyzed. Specifically, the first two examples show how it is possible to improve the solution quality in the structured-grid framework by means of (i) a mesh adaptation technique that does not involve topological changes in the mesh, or (ii) an overlapping mesh technique whereby additional blocks are superimposed to a background computational mesh. Finally, the third example shows that the same target can be achieved, with greater flexibility, using an unstructured, anisotropic mesh adaptation technique.

The use of different shock-capturing codes, both commercial and in-house developed, structured, and unstructured, contributes towards making these observations of broad validity.

Nevertheless, more relevant and conclusive enhancements in terms of solution quality can be obtained by means of a shock-fitting technique for unstructured grids. Indeed, the use of this newly developed technique completely removes the capturing process that is the source of most of the problems that plague current shock-capturing schemes. The proposed shock-fitting technique heavily relies upon unstructured mesh adaptation/refinement/generation techniques and it is worth underlining that only the development of advanced algorithms and the availability of public domain tools for unstructured mesh generation have made the development of this novel shock-fitting technique possible.

All the analyses presented in this paper deal with steady flows and have been conducted using second-order accurate discretizations in space. We believe, however, that most of the conclusions that have been drawn are applicable to higher order discretizations and unsteady flows as well.

Very high-order methods on structured grids have been extensively used since almost two decades to perform LES or DNS of compressible turbulence, as reviewed in the two recent articles by Johnsen et al. [8] and Pirozzoli [64]. In this particular context, shock-capturing schemes are severely challenged by two conflicting requirements: on the one hand, the discretization scheme should be minimally dissipative in order to be able to adequately simulate all relevant turbulent scales; on the other hand, it must supply a sufficient amount of artificial dissipation around the shocks to avoid un-physical oscillations. This task can be accomplished using different discretization schemes, away and near the discontinuities, which are blended by means of a shock sensor. A number of variations of these hybrid shock-capturing techniques are described in the two aforementioned papers, along with their merits and limitations.

Although structured grid adaptation or refinement is hardly foreseeable in the context of LES or DNS, the advantages offered by shock-fitting, whenever the flow topology makes it applicable, are clearly acknowledged by the “compressible’’ LES/DNS community [8, 64]. It is also true, however, as demonstrated using the test-cases presented in [8], that not only is the accurate treatment of the shock wave important to achieve accurate results in the entire domain, but also the properties of the numerical scheme used in smooth regions of the flowfield. It follows that the coupling between high-order schemes and improved shock modeling (such as the high-order shock-fitting technique presented in [53]) appears fully justified, at least in the context of LES/DNS.

The last two decades have also seen a growing interest in the development of high-order methods for unstructured grids [65]. These are expected to be more efficient than second-order accurate schemes, at least for selected applications requiring high accuracy. High-order methods on unstructured grids naturally lend themselves to the coupling with anisotropic mesh adaptation [66]. From the viewpoint of mesh generation/adaptation, high-order methods are less demanding than second-order accurate ones as they require coarser meshes than the latter to achieve the same error level. This is certainly true within boundary layers and wakes, but not around shocks, since also high-order methods need to locally downgrade their accuracy to first order around captured discontinuities. In this respect, the coupling between our newly developed unstructured shock-fitting algorithm and high-order unstructured grid schemes looks particularly promising.