4.1. General features of flow fields

The main focus of the present study is the relatively small region of the incident oblique shock wave reflection from the plane of symmetry. Figure 3 presents the numerical results of the distributions of dimensionless temperature in this region (see (4.1)) obtained by the NSF, R13 equations and DSMC. The point with the coordinates $(0,0)$ is the reflection point obtained analytically from the inviscid solution (see figure 2). The black lines present the shock-wave discontinuities in the inviscid solution.

Figure 3. Temperature distribution fields for NSF (a), R13 (b) and DSMC (c).

All three considered viscous flow models predict non-uniform temperature behind the RS that shows essentially viscous and non-continuum effects since the inviscid continuum solution given by the Euler equations predicts the uniform flow behind the RS. In particular, they all predict the appearance of a 2-D maximum (overshoot) to the right of the point $(0,0)$. The zone with the temperature maximum behind the RS is called the non-RH zone (Sternberg Reference Sternberg1959; Ivanov et al. Reference Ivanov, Bondar, Khotyanovsky, Kudryavtsev and Shoev2010; Shoev et al. Reference Shoev, Kudryavtsev, Khotyanovsky and Bondar2023) (this term was originally coined for the flow near the triple point in the irregular shock reflection). Indeed, the temperature value behind two oblique shocks given by RH conditions correspond to the orange colour in the legend. It is quite clear that all three approaches predict values of temperature that exceed the RH values.

The presence of a temperature extremum for the DSMC and R13 solutions might be expected in connection with similar results in the problem of the 1-D structure of the shock wave (see the § 4.2). The presence of an overshoot in the planar shock wave structure problem at large Mach numbers has been shown in many papers using both kinetic description (Elliott & Baganoff Reference Elliott and Baganoff1974; Erofeev & Friedlander Reference Erofeev and Friedlander2002; Dodulad & Tcheremissine Reference Dodulad and Tcheremissine2013; Malkov et al. Reference Malkov, Bondar, Kokhanchik, Poleshkin and Ivanov2015) and extended gas dynamics methods (Jin, Pareschi & Slemrod Reference Jin, Pareschi and Slemrod2002; Torrilhon & Struchtrup Reference Torrilhon and Struchtrup2004; Erofeev & Friedlander Reference Erofeev and Friedlander2007; Ivanov et al. Reference Ivanov, Kryukov, Timokhin, Bondar, Kokhanchik and Ivanov2012; Timokhin et al. Reference Timokhin, Bondar, Kokhanchik, Ivanov, Ivanov and Kryukov2015). In addition, it was shown that this maximum is observed at all types of molecular interaction potentials with the Mach number $Ma \geq 3.9$ (Yen Reference Yen1984; Erofeev & Friedlander Reference Erofeev and Friedlander2002). At the same time, it can be shown (Struchtrup Reference Struchtrup2005) that it is impossible to obtain a non-monotonic temperature distribution in a planar shock wave with the NSF equations. It is interesting to obtain a temperature maximum in the NSF solution of this 2-D problem. The same effect was also demonstrated with the NSF equations by Khotyanovsky et al. (Reference Khotyanovsky, Bondar, Kudryavtsev, Shoev and Ivanov2009). These results suggest that the appearance of such a maximum is not due to non-equilibrium effects as in the planar shock wave case, but rather can be explained by 2-D effects (see § 4.4). Another difference between these two cases is that the planar shock overshoot is located within the shock front while the present 2-D results clearly indicate that the maximum is separated from the RS front (see figure 3). For the sake of clarity, in what follows the term ‘temperature overshoot’ is used for the planar shock wave maximum (which is also observed in oblique shock waves) and the term ‘temperature maximum’ is applied to the maximum which is observed on the symmetry plane in the RR problem.

As it can be seen from the comparison of temperature flow fields in figure 3, the numerical results by the R13 equations turn out to be not only qualitatively but quantitatively very close to the DSMC kinetic solution. Figure 4 presents a comparison of temperature isolines for all three numerical approaches (figure 4a gives NSF versus R13, and figure 4b gives DSMC versus R13). The NSF numerical solution is far from the R13 numerical solution – R13 demonstrates greater thickness for both IS and RS waves. This disagreement is expected since the NSF equations fail to predict the internal structure of a shock wave accurately if the Mach number of the free stream is higher than $Ma=2.0$ (Pham-Van-Diep, Erwin & Muntz Reference Pham-Van-Diep, Erwin and Muntz1991). The R13 results agree with the DSMC benchmark solution quite well. The quantitative agreement between the two methods improves with the decreasing of the local Mach number (as it moves downstream). This R13 equations’ behaviour is similar to that in the 1-D shock wave structure problem, where the fast upstream part of a shock wave is predicted much more poorly than the downstream slower flow part (Timokhin et al. Reference Timokhin, Struchtrup, Kokhanchik and Bondar2017).

Figure 4. Temperature isolines for NSF–R13 (a) and R13–DSMC (b).

Figure 5 presents the comparison of the density and temperature profiles for all three models over horizontal cross-sections at different distances from the symmetry plane. The curves at finite distances from the symmetry plane consist of distinct incident oblique and RS-wave profiles and areas of nearly uniform flow (zones (1), (2), (3) in figure 1) while the curves for the $y=0$ cross-section resemble ordinary planar 1-D shock-wave profiles (IS and RS fronts are merged into one). The profiles along the symmetry plane for all three models also clearly reveal a temperature maximum which was observed in the flow fields. The maximum exceeds 5 per cent for the DSMC and R13 and is much smaller for the NS. The density profile along the symmetry plane especially for the R13 and DSMC case lies below the RH value of density behind IS and RS even a hundred mean free paths downstream of the reflection point. It is consistent with the long wake observed in the temperature flow fields. It is interesting that at the distance of five and 10 mean free paths from the symmetry plane all the density profiles reach the RH value rapidly as in a typical shock wave, although tiny overshoots are observed in the DSMC and R13 density profiles.

Figure 5. Profiles of density (a) and temperature (b) for NSF, R13 and DSMC along $x$ at $y=0.0$, $y=5\lambda _\infty$, $y=10\lambda _\infty$ and $y=15\lambda _\infty$. Dashed lines show RH values behind incident (bottom) and reflected (top) shocks.

The NSF equations provide results that are significantly different from the benchmark DSMC solution. As can be expected, disagreement is observed inside both IS and RS: the NSF model significantly underpredicts thickness of the shock front and does not predict a temperature overshoot inside the IS at various distances from the symmetry plane (a purely 1-D effect predicted by the R13 and DSMC for high Mach numbers). This 1-D overshoot is shown in the zoomed zone of the figure 5(b). The other difference is that the effects mentioned above (the temperature maximum on the symmetry plane and long downstream density tail lying below the RH value) are not prominent in the NSF solution. On the contrary, the distributions of the macroparameters of the R13 moment solution are quantitatively consistent with the DSMC kinetic solution in all cross-sections. The exceptions include the overestimated 1-D temperature overshoot in the IS and temperature distributions in the region of formation of the leading front of shock waves. Both of these effects are well known (see the results of Timokhin et al. (Reference Timokhin, Bondar, Kokhanchik, Ivanov, Ivanov and Kryukov2015) and Timokhin et al. (Reference Timokhin, Struchtrup, Kokhanchik and Bondar2017)). Another small difference is slight overestimation of the temperature maximum on the symmetry plane by the R13 equations.

Figure 6 presents the results for temperature and density in the cross-sections perpendicular to the flow direction in the free stream: along $x=-20.6\lambda _\infty$, $x=0$ and $x=20.6\lambda _\infty$ lines. Two velocity components in the same cross-sections are shown in figure 7. Solid black lines present the inviscid solution on both sides of the shock waves. The dashed vertical lines show the positions of the IS and RS. The R13 results inside the fronts of oblique shock waves (parts of the profiles for $x = -20.6\lambda _\infty$ and $x =20.6\lambda _\infty$) turn out to be much closer to the reference DSMC solution than the NSF equations solution. The qualitative behaviour of the solutions by all three approaches is similar for the main flow macroparameters.

Figure 6. Profiles of density and temperature for NSF, R13 and DSMC along $y$ at $x=-20.6\lambda _\infty$, $x=0$ and $x=20.6\lambda _\infty$.

Figure 7. Profiles of velocity components for NSF, R13 and DSMC along $y$ at $x=-20.6\lambda _\infty$, $x=0$ and $x=20.6\lambda _\infty$.

Note, that the inviscid solution for temperature, density and flow velocity is constant for $x=0$ and equal to their values behind the IS. For all three models of the viscous flow one observes clearly non-constant solutions in this cross-section. In particular, temperature exhibits more than a 50 per cent maximum in comparison with its value behind the IS. The NSF equations underpredict the temperature in this zone while the R13 equations overpredict the maximum value of the temperature in the vicinity of the origin. There are also visible non-uniformities in the R13 and NS density profiles and an $x$-velocity minimum in the origin. The value of this minimum is overpredicted by the R13 model while the NSF model predicts it accurately.

The quantitative difference between the three models under consideration for all the variables become smaller with increasing $x$-coordinate. Most significant differences are observed for the upstream cross-section ($x=-20.6\lambda _\infty$), while the downstream cross-section ($x=20.6\lambda _\infty$) reveals the smallest difference between all three models. In this case, the results of the R13 equations almost repeat the DSMC data profiles. This trend correlates with the decrease in the local Mach number and increase of density of the flow downstream. It is interesting to notice that in the downstream cross-section near the symmetry plane the DSMC and R13 models predict temperature higher and $x$-velocity lower than that given by the RH conditions. It agrees with the figures shown above which confirm the existence of the large non-RH-zone behind the point of the shock-wave reflection.

4.2. Comparison of structures of the planar shock wave and the zone of the RR of an oblique shock wave

The flow fields and profiles of macroparameters along the symmetry plane presented in the previous subsection demonstrates that the area in the vicinity of the reflection point where IS and RS merge resembles an ordinary normal shock (such as Mach stem in the case of Mach reflection, see e.g. Ben-Dor Reference Ben-Dor2007) with macroparameter isolines normal to the symmetry plane and a steep rise of density and temperature along the streamline. The present section is devoted to the analysis of the similarities and differences in the macroparameters’ distributions of the NSF, R13 and DSMC solutions in the 1-D normal shock and along the symmetry plane in the 2-D RR flow. Recall, that we consider the normal shock wave with Mach number $Ma= 8.0$ equal to the normal Mach number for the IS in the RR case, $Ma_n=8.0$ to ensure comparable degree of thermal non-equilibrium inside the shocks in both flows. Below, macroparameters’ comparison is presented starting with zeroth- and first-order moments of the distribution function (density and velocity) and moving on to higher-order moments such as temperature, viscous stress tensor and heat flux. The origin for the normal shock wave case is located in the shock centre – a point with the density equal to the mean value of the densities upstream and downstream of the shock.

The figures 8, 9 and 10 present density, velocity and temperature profiles for both flows obtained with all three gas flow models. The black lines here illustrate the inviscid solutions. As it can be seen from the profiles, the thickness of the 1-D shock wave (characterized by maximal slope of the profile) is approximately 10 mean free paths, while the thickness of the 2-D structure along the symmetry plane in the RR case is more than two times greater. The reference DSMC solution in the 2-D RR case reveals that a very long tail is observed after the steep part of the profile with all the macroparameters reaching the RH values rather slowly in comparison with the normal shock case. While it takes approximately five mean free paths from the origin to reach the RH values in the normal shock case, in the 2-D RR case density, velocity and temperature are substantially different from the RH values 35 mean free paths downstream of the origin, or in other words the non-RH-zone is observed. The velocity profiles for the two cases considered have a qualitative difference in the region behind the shock wave front. A minimum of velocity is observed in the problem of RR in the non-RH-zone (figure 9b). Its position coincides with the position of the temperature maximum (figure 10b).

Figure 8. The NSF, R13 and DSMC density profiles for the 1-D shock structure problem (a) and 2-D RR problem (b).

Figure 9. The NSF, R13 and DSMC velocity profiles for the 1-D shock structure problem (a) and 2-D RR problem (b).

Figure 10. The NSF, R13 and DSMC temperature profiles for the 1-D shock structure problem (a) and 2-D RR problem (b).

The NSF equations greatly underestimate the thickness of the shock structures in both cases (the slope of the profiles is almost two times steeper than for the DSMC solution). The R13 solution agrees well with the DSMC one for density. A significant difference is observed in velocity and temperature at the leading edge of the shock wave which is manifested in a longer upstream tail of the DSMC profile (figures 9a and 10a). A detailed discussion of the behaviour of shock-wave solutions of various variants of the R13 system is given by Timokhin et al. (Reference Timokhin, Struchtrup, Kokhanchik and Bondar2017).

As it was mentioned above, both the R13 equations and the DSMC method predict temperature overshoot in the planar shock wave (see figure 10a) while the NSF solution is monotonic. For $Ma_\infty =8$ the size of the overshoot in the DSMC solution is approximately $1\,\%$ with respect to the value behind the shock front. The overshoot in the R13 solution is much bigger (approximately $3\,\%$), and the NSF temperature profile does not have this peak at all. For the flow along the plane of symmetry (see figure 10b), all models considered predict a temperature maximum, and its value in the DSMC and R13 solutions is approximately 5 % with respect to the downstream temperature, which is substantially greater than the maximum predicted by the NSF solution (approximately 2 %).

The distribution of temperatures associated with thermal motion of molecules in different directions can provide additional information on the structure of the flows under consideration. Gas temperature can be represented as follows:

(4.2)\begin{equation} T= \tfrac{1}{3} (T_x + T_y + T_z), \end{equation}

where $T_x$, $T_y$ and $T_z$ are kinetic temperatures defined by mean energy of thermal motion of molecules in the $x$, $y$ and $z$ directions (temperature components or $x$-, $y$- and $z$-temperature, respectively). The following expressions relate these temperatures to the diagonal components of stress tensor in the non-dimensional form (4.1):

(4.3)$$\begin{gather} \theta_x=T_x=2(p+\sigma_{xx})/\rho, \end{gather}$$

(4.4)$$\begin{gather}\theta_y=T_y=2(p+\sigma_{yy})/\rho, \end{gather}$$

(4.5)$$\begin{gather}\theta_z=T_z=2(p+\sigma_{zz})/\rho. \end{gather}$$

In the problem of the structure of a planar shock wave, so called transverse temperatures $T_y$ and $T_z$ are equal to each other. So an overall gas temperature (4.2) in the 1-D case can be written as $T= \frac {1}{3} ( T_x + 2T_y )$.

Figures 11 and 12 show the comparisons of the temperature components for both problems obtained using all three approaches. In figure 11 NSF and DSMC results are presented. As for other macroparameters the NSF solution predicts much steeper slope and much less gradual onset of the profiles for both problems. Except for these differences the NSF and DSMC profiles are qualitatively similar for the planar shock wave case: streamwise temperature $T_x$ has a significant maximum while transverse temperature $T_y$ increases monotonically and the RH postshock value is reached within five free stream mean free paths by both temperatures. Moreover, the value of the streamwise temperature maximum in the planar shock wave is similar for both gas flow models, which is predicted by Yen's solution (Yen Reference Yen1966). The profiles of temperature components for the 2-D RR profile reveal the similar differences between the NSF and DSMC results from the planar shock wave cases witnessed above: the zone of steep gradients is approximately two times thicker and a very long non-RH-zone is observed downstream.

Figure 11. Different components of temperature predicted by the NSF and DSMC. The 1-D shock-wave structure (a) and internal structure of RR along the symmetry plane (b).

Figure 12. Different components of temperature predicted by the R13 and DSMC. The 1-D shock-wave structure (a) and internal structure of RR along the symmetry plane (b).

All temperature distributions provided by the R13 are much closer to the DSMC data for both problems (see figure 12) as it can be expected from the profiles of the macroparameters analysed above. However, longer upstream tails in the DSMC distributions are observed similar to other macroparameters. Note, that in the RR $x$ and $y$ temperature components are overpredicted by the R13 near the origin. These overestimated values of these two temperature components lead to an overestimated value of the total temperature in this region (see figure 5b for $y=0$).

The most significant qualitative difference in distributions of directional temperatures between the two problems under consideration, is that in the RR the mentioned streamwise temperature maximum is not observed at all. Instead all three approaches predict prominent maxima of $T_y$ temperature near the origin, while $T_x$ and $T_z$ temperature are very close to each other and steadily increase up to the point of equalization of all three directional temperatures. Farther downstream all the temperatures decrease down to the RH value, so, strictly speaking, all three temperatures exhibit a non-monotonic behaviour. This $T_y$ maximum can be explained by a sharp flow deceleration along the $y$-direction to a zero value of velocity $y$-component on the plane of symmetry in the vicinity of the point of origin (0,0). This part of the kinetic energy goes into thermal molecular motion. Obviously, the translational $y$-temperature must get more energy. At the same time, the molecules need time to undergo a sufficient number of collisions to dissipate thermal motion over all translational molecular degrees of freedom. This local thermal non-equilibrium leads to such behaviour of directional temperatures observed in figures 11(b) and 12(b).

Figures 13 and 14 present a comparison of the components of the viscous stress tensor and the longitudinal heat flux component in the same format as above. The heat flux is the highest moment of the local distribution function among considered. Therefore, one can expect a greater difference in the results of continuum approaches from the kinetic solution. Figures 13(a) and 14(a) compare the results of the NSF, R13 and DSMC for the normal shock. As it can be seen, the results of the R13 equations are in good agreement with the low-speed part of the DSMC shock wave front. The extrema for both macroparameters are predicted fairly accurately by the R13 equations. At the same time, significant differences are observed in the high-speed part of the structure. The NSF equations provide ${\mathsf{\sigma}} _{xx}$ and $q_x$ distributions which is significantly different from the DSMC ones: the value of the viscous stress maximum is overpredicted by 20 % while the value of the heat flux minimum is underpredicted by the same degree.

Figure 13. The $\sigma$-components predicted by NSF, R13 and DSMC. The 1-D shock-wave structure (a) and internal structure of RR along the symmetry plane (b).

Figure 14. The $q_{x}$ predicted by NSF, R13 and DSMC. The 1-D shock-wave structure (a) and internal structure of RR along the symmetry plane (b).

Due to the presence of a large transverse gradient of the $y$-velocity component (see § 4.4) near the RR point, the behaviour of the components of the viscous stress tensor differs qualitatively from the results of the 1-D problem. Figure 13(b) shows the distributions of the components ${\mathsf{\sigma}} _{xx}$ and ${\mathsf{\sigma}} _{yy}$ obtained by all the proposed methods. Figure 14(b) presents a similar comparison of the longitudinal heat flux $q_x$. As it can be seen from the both figures, the difference between the R13 and the DSMC solution appears to be significant in the entire considered region. The absolute values of the extrema are underestimated by the R13 equations (the differences reach approximately 20 %). Moreover, the NSF equations overestimate the absolute values of extrema by up to 30 %. The qualitative behaviour of the viscous stress tensor diagonal components is in agreement with the distributions of directional temperatures (see figures 11 and 12) which are related to them by (4.3)–(4.5).

4.3. The 2-D effects in the confluence of the incident and reflected shock wave

The present subsection aims at explaining the effects in the flow along the symmetry plane streamline in the RR discussed in the previous subsection. Let us first start with the flow at some distance from the symmetry plane where the IS and RS fronts can be separated from one another.

As mentioned above, the considered free stream Mach number $Ma=20$ and the wedge angle $\theta _w=17.061$ correspond to the normal Mach number $Ma_n=8$ of the incident oblique shock wave. The normal Mach number for the oblique RS wave in the considered case is $Ma_n=2.1$. It is well known that the RH conditions across any oblique shock are similar to the normal shock if one formulates them in terms of the normal Mach number and takes into account that the tangential flow velocity component is equal on both sides of the shock front (see e.g. Landau & Lifshitz Reference Landau and Lifshitz1987). Moreover, the problem of internal structure of an oblique shock wave can be reduced to a problem of internal structure of a normal shock by considering it in the frame of reference which is moving with the component of the free stream velocity tangential to the shock front. In this frame of reference tangential gas velocity becomes equal to zero and the problem formulation fully coincides with that of the 1-D normal shock. Hence the structures of the normal and oblique shocks differ only in the presence of the constant tangential component of flow velocity in the latter case.

The profiles across the IS and RS along the free stream direction ($x$ axis, symmetry plane) at different $y$-coordinate were compared with the planar shock wave profiles computed at $Ma=8.0$ and $Ma=2.1$. For the sake of such comparison the normal shock profiles were linearly transformed from the shock-front reference frame to the laboratory reference frame. The scheme of such a transformation is shown in figure 15. A planar shock wave lies schematically between two blue lines. The red segment AB is the shock wave profile along the normal to the shock front, and the green segment AC is the result of the transformation (non-orthogonal projection in the direction parallel to the shock front onto the $x$ axis). In order to obtain the required shock wave profile along the AC segment the coordinate across the original 1-D normal shock profile (AB segment) is divided by the sine of the oblique shock angle.

Figure 15. Scheme of comparison between 1-D and 2-D numerical results.

The comparison series of the NSF results of transformed temperature distributions for two 1-D shock waves and the results for the 2-D RR problem is presented in figure 16 in decreasing order of distance to the symmetry plane (for $y=\{15\lambda _\infty$, $8\lambda _\infty$, $3\lambda _\infty$ and $0\lambda _\infty \}$). The red and green curves represent the projections of the 1-D profiles of the IS and RS, respectively. Symbols indicate the results of the regular reflection computation. Far from the reflection point, the results of two 1-D shock waves exactly repeat the 2-D solution (see the first two plots in figure 16 for $y = 15\lambda _\infty$ and $y=8\lambda _\infty$). This indicates the absence of 2-D effects of the investigated flow at such a distance. As one approaches the symmetry plane, at a distance of the order of a mean free path, 2-D effects start playing a significant role in the distribution of macroparameters in the vicinity of the reflection point: two shock waves merge in the RR problem (see the other two plots in figure 16 for $y = 3\lambda _\infty$ and $y = 0\lambda _\infty$). In addition to the difference in the region around the point $x = 0$, the temperature maximum is observed in the 2-D problem when it cannot be predicted by 1-D computations.

Figure 16. Comparison between transformed 1-D temperature distributions and temperature distributions in the RR problem at $y=$ const. (NSF).

Similar results obtained with the DSMC method are presented in figure 17. At a distance of 15 mean free paths from the symmetry plane the IS and RS do not affect each other (2-D effects are negligible) and the RR profile coincides with the 1-D solutions as in the NSF case. At a distance of eight mean free paths, however, some minor 2-D effects are seen in the profile of the RS: temperature immediately behind the shock is slightly smaller than in the 1-D case and approaches the RH value farther downstream. This result is consistent with temperature isolines presented in figure 4. In the other two plots for $y = 3\lambda _\infty$ and $y = 0\lambda _\infty$ one can observe clear 2-D merging of two shocks with a much more prominent temperature maximum in the DSMC solution than in the NSF one. The R13 model provides qualitatively similar results to the DSMC (not presented in the plots).

Figure 17. Comparison between transformed 1-D temperature distributions and temperature distributions in the RR problem at $y=$ const. (DSMC).

Profiles of directional temperatures $T_x$, $T_y$ and $T_z$ for various distances from the symmetry plane are presented in figures 18 and 19 for all three gas flow models. Despite the clearly demonstrated fact that internal structure of an oblique shock coincides with that of a normal shock, these profiles exhibit unexpected pronounced overshoots of $T_y$ inside the fronts while $T_x$ and $T_z$ are monotonic. Recall, that an overshoot of $T_x$ is typical for a normal shock. This seeming contradiction is explained by the choice of the axis direction in these two problems. Indeed, the $T_x$ temperature in the RR problem is not associated with thermal motion of molecules in the direction normal to the front but rather in the direction of the free stream. Due to significant inclination of both IS and RS with respect to the free stream (more than 45 degrees) the thermal motion of the molecules in the direction normal to the shocks mainly contribute to the $T_y$ temperature which explains its overshoots inside the shock fronts. With decreasing distance from the symmetry plane, $T_y$ overshoots start to merge and eventually form the sole overshoot at the symmetry plane for all three gas flow models under consideration. This fact clearly explains the presence of the $T_y$ overshoot discussed in the previous subsection as a result of the confluence of two shocks each of which has its own $T_y$ overshoot.

Figure 18. Temperatures $T_x$, $T_y$ and $T_z$ at $y=$ const. for the RR problem. The NSF and DSMC results: $y=10\lambda _\infty$ (a), $y=5\lambda _\infty$ (b), $y=3\lambda _\infty$ (c), $y=0\lambda _\infty$ (d).

Figure 19. Temperatures $T_x$, $T_y$ and $T_z$ at $y=$ const. for the RR problem. The R13 and DSMC results: $y=10\lambda _\infty$ (a), $y=5\lambda _\infty$ (b), $y=3\lambda _\infty$ (c), $y=0\lambda _\infty$ (d).

4.4. Analyses of the flow with conservation equations

The results on RR presented above clearly demonstrate deviation from the RH values downstream of the reflection point, e.g. the existence of a non-RH zone. Recall, that the RH relations can be considered as a corollary of the 1-D conservation laws. Deviation from the RH relations can be explained by 2-D effects in the region of the shock reflections. In order to clarify which processes are important in this regard an analysis of mass, momentum and energy conservation equations is performed for the symmetry plane streamline. In addition, the normal shock problem is also considered. Note, that the conservation equations are presented in the general form. All of them directly follow from the Boltzmann equation and must hold for any dilute gas flow model, including the three models employed in the present paper. As the R13 equations provide results which are qualitatively consistent with the DSMC ones, only the NSF and R13 models were used in the analyses.

The mass and momentum conservation equations for a steady 2-D flow along the plane of symmetry streamline can be written in terms of the substantial derivatives

(4.6)\begin{equation} \frac{{\rm D}}{{\rm D} t} = \frac{\partial}{\partial t} + v_x \frac{\partial}{\partial x} +\underline{v_y \frac{\partial}{\partial y}} = v_x \frac{\partial}{\partial x} +\underline{v_y \frac{\partial}{\partial y}} \end{equation}

as follows:

(4.7)$$\begin{gather} \frac{{\rm D} \rho}{{\rm D} t} ={-} \rho \left(\frac{\partial v_x}{\partial x} + \underline{\frac{\partial v_y}{\partial y}}\right), \end{gather}$$

(4.8)$$\begin{gather}\frac{{\rm D} v_x }{{\rm D} t} ={-} \frac{1}{\rho} \left(\frac{\partial p}{\partial x} + \frac{\partial {\mathsf{\sigma}}_{xx}}{\partial x} + \underline{\frac{\partial {\mathsf{\sigma}}_{xy}}{\partial y}}\right). \end{gather}$$

The underlined terms are caused by the two-dimensionality of the flow. The rest of the terms remain in the 1-D case of a planar shock wave structure problem.

The energy conservation law along the symmetry plane has the following form:

(4.9)\begin{equation} \frac{{\rm D} E }{{\rm D} t} ={-} \frac{1}{\rho} \left(\frac{\partial v_xp}{\partial x} + \frac{\partial v_yp}{\partial y}+ \frac{\partial v_x{\mathsf{\sigma}}_{xx}} {\partial x} + \frac{\partial v_x{\mathsf{\sigma}}_{xy}}{\partial y}+ \frac{\partial v_y{\mathsf{\sigma}}_{yy}}{\partial y} + \frac{\partial q_x}{\partial x}+ \frac{\partial q_y}{\partial y} \right), \end{equation}

where $E$ is the total internal energy per unit mass containing kinetic $E_K = {v^2}/{2}$ and internal $E_I = e$ parts,

(4.10)\begin{equation} E = E_K + E_I = \frac{v^2}{2} + e = \frac{v^2}{2} + \frac{3}{4} T. \end{equation}

This form of the conservation laws allows one to analyse an evolution of density, velocity and total internal energy of a ‘fluid element’ of a constant mass and variable volume as it moves along the streamline (see e.g. Anderson Reference Anderson2019). The rates of change of these three parameters are determined by the terms on the right-hand side of the equations, which describe various physical processes, to be discussed below in detail. The term ‘fluid element’ typical of conventional continuum fluid dynamics is used hereafter for convenience. Note, that it should not be understood literally because all three considered gas flow models presume exchange of molecules between such an imaginary element and the gas medium which surrounds it.

To illustrate the influence and physical meaning of the right-hand part terms of the energy conservation equation (4.9) let us rewrite it taking into account that the terms $v_y({\partial p}/{\partial y})$, $v_y ({\partial {\mathsf{\sigma}} _{yy}}/{\partial y})$ and ${\mathsf{\sigma}} _{xy} ({\partial v_x}/{\partial y})$ are equal to zero on the symmetry plane

(4.11)\begin{align} \frac{{\rm D} E }{{\rm D} t} &={-} \frac{1}{\rho} \left( v_x \frac{\partial p} {\partial x} + p \frac{\partial v_x} {\partial x}+ v_x \frac{\partial {\mathsf{\sigma}}_{xx}} {\partial x} + {\mathsf{\sigma}}_{xx} \frac{\partial v_x} {\partial x} + \frac{\partial q_x} {\partial x} \right. \nonumber\\ &\quad \left.+ \underline{v_x \frac{\partial {\mathsf{\sigma}}_{xy}} {\partial y} } + \underline{p \frac{\partial v_y} {\partial y} } + \underline{{\mathsf{\sigma}}_{yy} \frac{\partial v_y} {\partial y} } + \underline{\frac{\partial q_y} {\partial y} }\right), \end{align}

where the underlined terms are present only in the 2-D RR problem and hence describe 2-D effects on the symmetry plane. Equation (4.11) can be rewritten into two governing equations for kinetic and internal parts of $E$ (see e.g. Anderson Reference Anderson2019) as follows:

(4.12) $$\begin{gather} \frac{{\rm D}}{{\rm D} t} \left( \frac{v^2}{2} \right) = \underbrace{- \frac{1}{\rho} \left( v_x \frac{\partial p} {\partial x} + v_x \frac{\partial {\mathsf{\sigma}}_{xx}} {\partial x} \right)}_{{P_{xx}\ gradient\ power}\quad} \underbrace{-\underline{\frac{1}{\rho}\left(v_x \frac{\partial {\mathsf{\sigma}}_{xy}} {\partial y}\right)}}_{viscous\ friction\ power}, \end{gather}$$

(4.13)$$\begin{gather}\frac{{\rm D}}{{\rm D} t} \left( e \right) = \underbrace{ - \frac{1}{\rho} \left( p \frac{\partial v_x} {\partial x} + {\mathsf{\sigma}}_{xx} \frac{\partial v_x} {\partial x} + \underline{p \frac{\partial v_y} {\partial y} } + \underline{{\mathsf{\sigma}}_{yy} \frac{\partial v_y} {\partial y} } \right)}_{adiabatic\ compression\ power} \underbrace{- \frac{1}{\rho} \left( \frac{\partial q_x} {\partial x} + \underline{\frac{\partial q_y} {\partial y}}\right)}_{heat\ input\ power}. \end{gather}$$

Note, that the first kinetic energy equation is, strictly speaking, redundant, because it follows directly from the momentum conservation equation.

Figures 20–23 present the distribution of various terms on the right-hand sides of the conservation equations obtained from NSF and R13 numerical solutions of the planar shock and RR problems.

Figure 20. Mass conservation terms for 1-D shock-wave structure (a) and internal structure of RR along the symmetry plane (b) for the NSF equations. Panels (c,d) are similar results for the R13 equations.

Figure 21. Momentum conservation terms for 1-D shock-wave structure (a) and internal structure of RR along the symmetry plane (b) for the NSF equations. Panels (c,d) are similar results for the R13 equations.

Figure 22. Energy conservation terms (kinetic part) for 1-D shock-wave structure (a) and internal structure of RR along the symmetry plane (b) for the NSF equations. Panels (c,d) are similar results for the R13 equations.

Figure 23. Energy conservation terms (internal part) for 1-D shock-wave structure (a) and internal structure of RR along the symmetry plane (b) for the NSF equations. Panels (c,d) are similar results for the R13 equations.

The terms of the mass conservation equation (4.7) shown in figure 20 demonstrate the compression of the fluid element caused by the velocity gradient. The values of all terms are higher for the NSF solution than for the R13 one. This can be explained by the fact the R13 equations predict a thicker shock wave front. For the planar shock wave the net compression across the shock wave is similar for all the models and is given by the RH density value. This is why a thicker shock front must lead to the lower compression terms’ values. The compression terms caused by the $x$-velocity gradient are more than two times lower for the RR reflection than for the 1-D planar shock. The major contribution to the fluid element compression in the RR case is provided by the $y$-velocity gradient: it is one order of magnitude higher than that of the $x$-velocity gradient and is more than four times higher than the compression term for the planar shock problem. This fact demonstrates that even at the level of the mass conservation equation, the 2-D effects dominate the flow structure for the RR problem. The terms of the momentum conservation equation (4.8) shown in figure 21 differ similarly between the NSF and the R13 cases: on the one hand, the R13 shock front (or reflection region structure in the 2-D case) is approximately one and a half times thicker than the NSF one, and on the other hand the values of parameters are approximately one and a half times higher for the NSF model. There is also a noticeable feature of the R13 profiles which is also seen in some R13 profiles in figure 20: the profiles have many inflection points (extrema of the derivative) and therefore look less smooth than the NSF profiles.

The momentum equation terms describe the contribution of various forces to the deceleration of the fluid element. In both problems considered the normal stress $P_{xx}$ gradient clearly contributes to the deceleration everywhere from the undisturbed flow upstream to the undisturbed flow downstream (see figure 21). However, if we consider the contribution of the gradients of pressure $p$ and normal viscous stress ${\mathsf{\sigma}} _{xx}$ (recall that ${\mathsf{P}}_{xx}=p+{\mathsf{\sigma}} _{xx}$) that the situation qualitatively differs between the two problems considered. The pressure gradient also contributes to the deceleration everywhere in both problems. As for the normal viscous stress ${\mathsf{\sigma}} _{xx}$, in the case of the planar shock (figure 21a,c) the term related to its gradient decelerates the fluid element in the upstream high-speed part of the shock front but accelerates it in the downstream part. On the contrary, in the RR problem (figure 21b,d) normal stress accelerates the flow in the upstream part of the structure and decelerates farther downstream. This fact is consistent with the sign of the ${\mathsf{\sigma}} _{xx}$ gradient which is evident from the profiles for both problems shown in figure 13.

Additionally to the normal stress in the RR problem there is an effect of the shear stress (dashed black curves in figure 21) which is more than two times higher and opposite to that of the normal stress. The effect of the shear stress is quite evident because accordingly to the 2-D flow structure ‘gas layers’ which are located farther from the symmetry plane start passing through the IS first and decelerate the ‘gas layer of the symmetry plane’ due to shear stress in the upstream part of the structure. On the other hand, farther downstream the symmetry-plane layer is accelerated via the shear stress by the neighbouring layers which did not yet pass through the RS.

The profiles of kinetic energy conservation equation (4.12) terms shown in figure 22 are fully consistent with the considered results for the momentum conservation equation and reflect power produced by all considered forces which result in the change of kinetic energy of the fluid element. Recall, that the kinetic energy conservation equation is a corollary of the momentum conservation equation, and on the symmetry plane, as it can be seen, all the right-side terms of the former are obtained by the multiplication of the corresponding terms of the later by the absolute value of local gas velocity $v_x$. Nevertheless, this result cannot be considered completely redundant because it allows us to provide comparison of the presented terms with the right-side terms of the internal energy conservation equation (4.13) given in figure 23. Also it is worth noting, that the values of power produced by normal stress are more than two times higher in the 2-D problem than in the 1-D shock, which can be explained by much higher non-dimensional velocity values (see figure 9).

In the planar shock problem internal energy of a fluid element is changing due to heat conduction and compression by pressure forces which can be divided into pressure and normal viscous stress terms. Their profiles are given in figure 23(a,c). The contribution of compression to the internal energy change is positive everywhere while the effect of heat conduction is positive in the upstream cold part of the front and negative in the downstream cold one. The differences between the NSF and R13 results reflect the familiar pattern: the R13 equations predict a thicker shock front and smaller local values of the terms than the NSF equations. The values of these terms in the RR problem (see figure 23b,d) are significantly smaller, moreover, the term which describes power of compression due to ${\mathsf{\sigma}} _{xx}$ viscous stress has a minimum instead of maximum (this difference is explained by the behaviour of ${\mathsf{\sigma}} _{xx}$ inside the front for these two problems shown in figure 13).

The most striking feature of the internal energy balance in the RR problem is that for considered flow parameters the terms related to 2-D effects (their profiles are presented by dashed curves in figure 23) are an order of magnitude higher than 1-D terms shown by solid curves which are present also in the planar shock problem. Again the difference between the NSF and R13 solution manifests itself mainly in the width of the structure which is higher in the R13 and absolute values of the terms which are lower in the R13. The green and blue dashed curves which are related to the work of normal stress in the $y$-direction resulting in compression of the fluid element are positive and have maxima. Their relatively large values are consistent with the already-discussed effect of high compression of the fluid element along the $y$ axis (see figure 20). The term related to the change of the internal energy of the fluid element due to heat flux along $y$ axis (black dashed curve) is positive in the upstream cold part of the shock structure and negative in the downstream hot part being much higher in the absolute value. The value of the term in the downstream part is comparable to the power of compression due to pressure forces (green curves).

Generally one can conclude from the conducted analysis of kinetic and internal energy equations that the terms related to the 2-D effect such as acceleration and deceleration of the flow due two shear stress, compression due to pressure and normal viscous stress in the $y$-direction and heat conduction in the $y$-direction contribute to the energy of a fluid element on the RR problem comparably to the $x$-direction terms which present also in the 1-D case. While the kinetic energy for the most part is governed by a 1-D term related to deceleration due to the pressure gradient, the internal energy is determined by virtually only 2-D terms related to compression and heat conduction in the direction normal to the symmetry plane.

In order to analyse the net effect of the processes represented by the terms on the right-hand sides of the conservation equation on the parameters of the fluid element, let us consider integration of the conservation equations along the trajectory of the element. Below, (4.7), (4.8), (4.12) and (4.13) are presented in the following form:

(4.14)$$\begin{gather} \frac{{\rm D} \rho}{{\rm D} t} = \sum_{i=1}^2 \dot{\rho}_i, \end{gather}$$

(4.15)$$\begin{gather}\frac{{\rm D} v_x}{{\rm D} t} = \sum_{i=1}^3 \dot{v}_{xi}, \end{gather}$$

(4.16)$$\begin{gather}\frac{{\rm D}}{{\rm D} t} (E_K) = \frac{{\rm D}}{{\rm D} t} \left(\frac{v^2}{2} \right) = \sum_{i=1}^3 \dot{E}_i,\end{gather}$$

(4.17)$$\begin{gather}\frac{{\rm D}}{{\rm D} t} (E_I) = \frac{{\rm D}}{{\rm D} t} (e) = \sum_{i=4}^9 \dot{E}_i, \end{gather}$$

where $\dot {E}_i$ is a power of mechanical force or heat input. The expressions for $\dot {\rho }_i$, $\dot {v}_{xi}$ and $\dot {E}_i$ are presented in table 1 with their physical meaning and the NSF expressions through transport coefficients and gradients of pressure, temperature and velocity.

Table 1. The terms and their physical meaning in the conservation equations of mass, momentum and energy across the normal shock front and along the symmetry plane in the RR problem.

Increments of the density, velocity, kinetic and internal energy of a fluid element with respect to their free stream values can be calculated for any position $x$ of the element on the symmetry plane by evaluating the following integrals:

(4.18)$$\begin{gather} \Delta \rho (x)=\int_{-\infty}^{t(x)} \sum_{i=1}^2 \dot{\rho}_i \,{\rm d} t = \int_{-\infty}^{x} \sum_{i=1}^2 \frac{\dot{\rho}_i}{v_x} {{\rm d}\kern0.7pt x},\end{gather}$$

(4.19)$$\begin{gather}\Delta v_x (x)=\int_{-\infty}^{t(x)} \sum_{i=1}^3 \dot{v}_{xi} \,{\rm d} t = \int_{-\infty}^{x} \sum_{i=1}^3 \frac{\dot{v}_{xi}}{v_x} {{\rm d}\kern0.7pt x}, \end{gather}$$

(4.20)$$\begin{gather}\Delta E_K (x)=\int_{-\infty}^{t(x)} \sum_{i=1}^3 \dot{E}_i \,{\rm d} t = \int_{-\infty}^{x} \sum_{i=1}^3 \frac{\dot{E}_i}{v_x} {{\rm d}\kern0.7pt x}, \end{gather}$$

(4.21)$$\begin{gather}\Delta E_I (x)=\int_{-\infty}^{t(x)} \sum_{i=4}^9 \dot{E}_i \,{\rm d} t = \int_{-\infty}^{x}\sum_{i=4}^9 \frac{\dot{E}_i}{v_x} {{\rm d}\kern0.7pt x}. \end{gather}$$

The profiles of these increments across the normal shock and along the symmetry plane streamline in the RR problem are given in figure 24. Qualitatively the behaviour of the increments is almost similar for the normal shock and the RR problem: density and internal energy are increasing while velocity and kinetic energy are decreasing with a lower rate for the R13 and with a higher rate for the NSF equations. Quantitatively the increments in the RR problem are substantially higher (e.g. more than two times higher for both energies). Some differences between the two problems are manifested in subtle features like minima of velocity and kinetic energy increments present in the RR results. There is a maximum of internal energy in the RR increment of the internal energy, however, it is also present in the normal shock problem but only in the R13 and not in the NSF results. This is not surprising, because the internal energy is proportional to temperature with a constant coefficient and the temperature profiles exhibit similar features (see figure 10). Note, that while in the normal shock problem all increments approach analytical RH values behind the shock, in the RR problem the downstream values depend on the gas flow model which is explained by the presence of a long non-RH wake behind the shock reflection zone more prominent in the R13 model.

Figure 24. The NSF and R13 results for $\Delta \rho (x)$, $\Delta v_x (x)$, $\Delta E_K (x)$ and $\Delta E_I (x)$ in 1-D shock wave structure (a) and internal structure of RR along the symmetry plane (b). Black dashed lines are the analytical downstream values.

The contribution of each term $\dot {E}_i$ to the increments of kinetic and internal energies for a 1-D planar shock wave and the RR problem along the symmetry line can be evaluated as follows:

(4.22)\begin{equation} \Delta E_i (x)=\int_{-\infty}^{t(x)} \dot{E}_i \,{\rm d} t = \int_{-\infty}^x \frac{\dot{E}_i}{v_x} {{\rm d}\kern0.7pt x}. \end{equation}

The profiles of the increments of the terms of the kinetic energy conservation equation are presented in figure 25. As could be expected from the already-discussed distribution of the right-side terms of the kinetic energy conservation equation $\Delta E_1$ increment which reflects the work of the pressure gradient force is decreasing for both problems and both models (reducing the kinetic energy of the fluid element) while $\Delta E_2$, the work of the normal stress gradient force, has a minimum in the normal shock and a maximum in the RR problem. A less evident fact is that the net contribution of the $E_2$ term to the total kinetic energy is negative for the normal shock and positive for the RR problem (see the downstream values). Another interesting fact is that the net contributions of $E_1$ and $E_2$ depend on the gas flow model. At the same time the sum of these two contributions for the normal shock $\Delta (E_1+E_2)$ (total work of the normal stress gradient force which results in deceleration of the fluid element) which equals $\Delta E_K$ in the 1-D case do not depend on the model for the normal shock and is related to the RH downstream velocity. This is not the case for the RR problem: total work of the normal viscous stress gradient force depends on the mathematical model of the gas flow. There is also a contribution of the shear stress gradient to the net kinetic energy change in the RR problem: consistently with already discussed ${\mathsf{\sigma}} _{xy}$ profiles, the $\Delta E_3$ increment has a minimum inside the shock structure and, more interestingly, its net contribution to the kinetic energy change is negative and depends on the gas flow model used (it is less for the R13 than for the NSF equations).

Figure 25. The $\Delta E_i$ contributions (kinetic part) for 1-D shock-wave structure (a) and internal structure of RR along the symmetry plane (b) for the NSF equations. Panels (c,d) are similar results for the R13 equations. Black dashed lines in (a,c) are the analytical downstream values.

The increments of the terms describing contribution of various processes to the change of internal energy is given in figure 26. It is clear that the term related to the heat flux input $\Delta E_8$ has a maximum inside the normal shock and its net contribution to the internal energy of a fluid element passing the shock is equal to zero for both models considered. Indeed, if we one takes into account that in the 1-D stationary case mass flux is constant $\rho _x v_x = \rho _{\infty } v_{\infty }$ and heat flux asymptotically tends to zero in the free stream and downstream, then the net increment $\Delta E_8$ across the shock front is

(4.23)\begin{equation} \Delta E_8^{total}=\int_{-\infty}^{\infty} \frac{\dot{E}_8}{v_x} {{\rm d}\kern0.7pt x}={-}\int_{-\infty}^{\infty} \frac{1 }{{\rho_x v_x}} \frac{\partial q_x} {\partial x} {{\rm d}\kern0.7pt x} ={-} \frac{1 } {\rho_{\infty} v_{\infty}} \int_{-\infty}^{\infty} \frac{\partial q_x} {\partial x} {{\rm d}\kern0.7pt x} = 0. \end{equation}

Figure 26. The $\Delta E_i$ contributions (internal part) for 1-D shock-wave structure (a) and internal structure of RR along the symmetry plane (b) for the NSF equations. Panels (c,d) are similar results for the R13 equations. Black dashed lines in (a,c) are the analytical downstream values.

The contribution of pressure and viscous stress forces to the internal energy is positive in the 1-D shock and hence increments $\Delta E_4$ and $\Delta E_5$ are increasing across the shock and approach positive asymptotic values which depend on the gas flow model. Similarly to the to the case of $E_1$ and $E_2$ the total increment of the sum of $E_4$ and $E_5$ which equals net internal energy increment $\Delta E_I^{total} = \Delta (E_4+E_5)^{total}$ is completely defined by the RH conditions on temperature and therefore is independent of the gas flow model employed. Note, that the total increment of $E_5$ has an opposite sign and is equal in absolute value to $E_2$, which means the net contribution of viscous stress to the energy of a fluid element in the 1-D case equals zero. Indeed, in the 1-D case taking into account that ${\mathsf{\sigma}} _{xx}$ tends to zero with increasing distance from the shock, one obtains

(4.24)\begin{align} \Delta (E_2+E_5)^{total} &= \int_{-\infty}^{\infty} \frac{\dot{E}_2+\dot{E}_5}{v_x} {{\rm d}\kern0.7pt x} ={-} \int_{-\infty}^{\infty} \frac{1 } {\rho_{x} v_{x}} \left(v_x \frac{\partial {\mathsf{\sigma}}_{xx}} {\partial x} + {\mathsf{\sigma}}_{xx} \frac{\partial v_x} {\partial x} \right ) {{\rm d}\kern0.7pt x} \nonumber\\ &={-} \frac{1 } {\rho_{\infty} v_{\infty}} \int_{-\infty}^{\infty} \left( v_x \frac{\partial {\mathsf{\sigma}}_{xx}} {\partial x} + {\mathsf{\sigma}}_{xx} \frac{\partial v_x} {\partial x} \right ) {{\rm d}\kern0.7pt x} \nonumber\\ &={-} \frac{1 } {\rho_{\infty} v_{\infty}} \int_{-\infty}^{\infty} \frac{\partial (v_x{\mathsf{\sigma}}_{xx})} {\partial x} {{\rm d}\kern0.7pt x} = 0. \end{align}

This means the viscosity inside the shock is responsible for the transfer of some portion of kinetic energy which depends on the gas flow model to internal energy (viscous dissipation) but it does not contribute to the total energy of the fluid element passing the shock. Due to this fact the total increment of total energy across the normal shock equals $\Delta E = \Delta (E_1+E_4)^{total}$. This is why the total increment of the sum of $E_1$ and $E_4$ (asymptotic values for the green curves in figures 25a,c and 26a,c) which represents the total work of the pressure forces does not depend on the model and is defined by the RH conditions.

In the RR problem the situation is much more complicated. First, the total $x$-heat-flux input increment $\Delta E_8^{total}$ is positive and depends on the mathematical model of the gas flow. Second, the total viscous stress $x$-compression increment $\Delta E_5^{total}$ is negative, while the pressure $x$-compression term increment $\Delta E_4^{total}$ is positive and the total increment of their sum in contrast to the 1-D case depends on the model. Third, the net contribution of $x$-viscous stress to the total energy $\Delta (E_2+E_5)^{total}$ is not equal to zero and depends on the model as well. Fourth, the 2-D internal energy terms’ increments (shown in dashed curves) are almost one order of magnitude higher than the increments of the 1-D terms presented by solid curves. It is not surprising that the increments of the terms related to the compression of the fluid element along the $y$ axis provide the most significant positive impact on the total internal energy rise in the shock structure. The increment of the term $E_9$ related to the heat conduction along the $y$ axis behaves non-monotonically: first it increases and approximately five mean free paths before the origin starts to decrease and finally provides a very strong negative contribution to the internal energy which is several times higher than that of all the 1-D terms combined. All these contributions of the 2-D terms strongly depend of the gas flow model: e.g. the contribution of the viscous stress $y$-compression and $y$-heat-conduction terms to the net increment of internal energy are approximately one and a half times higher for the NSF equations than for the R13 model.

It is clearly seen from the slopes of the curves for the RR problem, especially in the R13 results, that the $y$-compression and $y$-heat-conduction terms continue to significantly contribute to the change of internal energy of a fluid element even as far as 20–30 mean free paths downstream of the origin. These coordinates correspond to what was referred as ‘wake’ in the flow fields. Indeed, the $y$-profiles of $v_y$ in figure 7 indicate the visible $y$-gradient of $v_y$ in the wake area on the symmetry plane, which produce the mentioned compression. The temperature non-uniformity along the $y$ axis of this downstream wake region evident from figure 6 may be the reason for the $y$-heat-flux non-zeroth divergence which contributes to the change of internal energy.

Recall, that the net contributions of viscous stress and heat conduction to the total energy of a fluid element passing the shock wave are both equal to zero for all dilute gas flow models. It leads to the fact that the change of total energy across the shock is determined by the work of pressure forces and hence the change of total enthalpy across the shock is equal to zero (see e.g. Shoev, Timokhin & Bondar Reference Shoev, Timokhin and Bondar2020). This statement can also be applied to a fluid element moving along a streamline in the RR problem sufficiently far from the symmetry plane and passing two separate (incident and reflected) oblique shocks. This is clearly not the case for the plane-of-symmetry streamline in the RR problem. Indeed, it is clear from figures 25b,d and 26b,d that the total contribution of viscous stress to the fluid element energy (given by the total contribution of $E_2$, $E_3$, $E_5$ and $E_7$ terms) is positive and is somewhat higher for the NS equations than for the R13.

For the kinetic energy, positive normal viscous $x$-stress contribution (acceleration) and negative shear stress contribution (deceleration) nearly compensate each other, yet there is a small net negative effect more pronounced for the R13 equations. Probably, the presence of the decelerating shear stress leads to a velocity minimum with values predicted by the RH relations. For the internal energy, negative normal viscous $x$-stress contribution (expansion) is more than one order of magnitude smaller than the most significant of all four terms’ $y$-stress contribution (compression). The total positive contribution of the viscous stress terms is almost compensated by the negative contribution of heat conduction which consists of a small positive contribution of $x$-heat flux component (heating) and several times larger negative contribution of $y$-heat flux (cooling). It is also worth noting that for the plane-of-symmetry streamline the net work of pressure forces (in the $x$ direction $\Delta (E_1+E_4)^{total}$ or in both directions $\Delta (E_1+E_4+E_6)^{total}$) depends on the gas flow model hence are not determined by the RH conditions.

In general, one can conclude that in contrast to the flow across any number of normal or oblique shocks the total energy of the fluid element passing through the zone of RR along the symmetry plane is significantly affected by the total non-zeroth contribution of viscosity and heat conduction. The most important effects are essentially 2-D, namely, compression by the normal viscous stress perpendicular to the stream line and heat-conduction cooling in the same direction. One can conclude from the curves corresponding to these processes in figures 25 and 26 that first the fluid element is compressed by the normal viscous stress in the $y$ direction which leads to the increase of its internal energy (viscous dissipation) and temperature values different to those predicted by the RH relations on the IS and RS. This difference manifests itself in the presence of a non-RH zone with high temperature in the form of the wake downstream of the shock reflection point. The excessive temperature leads to intensive heat transfer in the $y$ direction from the wake towards fluid elements passing two separate shocks and hence having lower temperatures predicted by the RH relations.