1. Introduction
Shock wave phenomena constitute fundamental characteristics of supersonic flows, with their systematic investigation forming the cornerstone of gas dynamics. Since the initial paper in the 1890s (Mach Reference Mach1893), shock wave research has evolved through distinct phases: early linear analyses of normal, oblique and conical shocks preceding the more complex curved shock investigations. The systematic investigation of curved shocks originated with the foundational work by Crocco (Reference Crocco1937) and Thomas (Reference Thomas1947) on zero-flow curvature and zero pressure gradient conditions, and was subsequently extended by Thomas (Reference Thomas1949) through detailed post-shock streamline curvature profiles. Gerber & Bartos (Reference Gerber and Bartos1960) established the coefficients for curved shock equations under steady, irrotational and uniform inflow assumptions. A seminal derivation by Emanuel & Liu (Reference Emanuel and Liu1988) of post-shock property gradients laid the groundwork for coordinate-transformation-based vector analyses (Emanuel Reference Emanuel2012). A significant breakthrough came with the derivation by Mölder (Reference Mölder2016) of curved shock equations for compound shock curvatures at any point along curved shocks, providing a simplified framework for analysing various shock wave dynamics phenomena. However, these equations remain constrained to scenarios where the flow plane aligns with the osculating plane, such as in conical flow with zero angle of attack. Mölder’s incorporation of vorticity relations, albeit limited to two dimensions, marked another important development. These theoretical advancements have enabled researchers to extend their investigations beyond two-dimensional curved shocks, addressing the more complex and practical scenarios of doubly curved shocks in non-uniform and rotational flows, which is a crucial condition for practical applications in aerospace engineering.
The transition to three-dimensional curved shock analysis began with Kanwal (Reference Kanwal1959), who pioneered the study of streamwise gradients of flow parameters following three-dimensional shocks. Serrin (Reference Serrin1959) subsequently developed a method for solving the gradients of flow parameters behind three-dimensional curved shocks, albeit with the limitation of requiring zero pressure gradient on the vertical flow surface. Tembharey & Sachdeva (Reference Tembharey and Sachdeva1973) expanded this research by establishing conditions for the existence of complex laminar flows and Beltrami flows following three-dimensional stationary curved shocks. Recent advances in this field come from Emanuel & Mölder (Reference Emanuel and Mölder2022). Benefiting from their previous work, they derived a three-dimensional curved shock theory for uniformly stable incoming flow and applied it to the analysis of gradient parameter distributions, vorticity and reflection coefficients behind blunt body shock waves. However, the existing theoretical framework exhibits inherent limitations in its mathematical formulation, particularly regarding its applicability to asymmetric three-dimensional flow. This fundamental constraint was recently overcome through the work of Zhang et al. (Reference Zhang, Yang, Shi, Zhu, Zheng, You and Wang2025). They established a shock-surface-based coordinate system by the method of moving frames, thereby achieving the solutions for arbitrary shock geometries.
While these theoretical developments have enhanced our understanding of three-dimensional curved shocks, they primarily reveal flowfield structures only near the shocks. Comprehensive analysis of three-dimensional shock-flowfields requires computationally efficient analytical tools. Although computational fluid dynamics (CFD) provides Euler equation solutions for three-dimensional inviscid flows (Wang Reference Wang2007; Laurence et al. Reference Laurence, Parziale and Deiterding2012), its computational demands and limited resolution hinder detailed flow structure analysis and inverse design. Consequently, algebraic methods, particularly the method of characteristics (MOC), have gained prominence in these applications (Harrop Reference Harrop1953). Zucrow & Hoffman (Reference Zucrow and Hoffman1976) conducted a detailed numerical study of MOC accuracy for axisymmetric, supersonic flows, encompassing both irrotational and rotational flows. The fundamental principle of MOC involves transforming quasi-linear partial differential equations into total differential equations, solved through numerical methods like the Euler predictor-corrector difference scheme. Chu & Nadir (Reference Chu and Nadir1976) extended MOC to unsteady flows using generalised finite-difference approximations, enabling calculations of unsteady two-dimensional supersonic flowfields. The advantages of MOC make it play a great advantage in flowfield analysis and inverse design. Building on supersonic flowfield analysis, You & Liang (Reference You and Liang2009) developed a design concept for three-dimensional section-controllable internal waverider hypersonic inlets, while Zhang (Reference Zhang2020) applied MOC to the inverse design of curved shock wave compression surfaces. The field witnessed significant progress with Shi et al. (Reference Shi, Zhu, You and Zhu2021) who developed the method of curved-shock characteristics (MOCC), which improved accuracy and efficiency through high-order parameters. Zhang et al. (Reference Zhang, Xu, Shi, Zhu and You2023) subsequently integrated Mölder’s curved shock theory with MOCC to construct a reflection model for planar flowfields, successfully predicting shock reflection structures and Mach stem.
Despite MOC’s extensive development and application in two-dimensional analyses, the inherent three-dimensional nature of real flows necessitates corresponding three-dimensional methods. Three-dimensional algebraic approaches primarily fall into two categories: (i) direct extensions of MOC and (ii) osculating methods. The former, exemplified by Thornhill (Reference Thornhill1952) in a hexahedral bi-characteristic network and by Sauerwein (Reference Sauerwein1964) in a stabilised implementation, achieves high accuracy through global iteration. Recent applications by Huang et al. (Reference Huang, Zuo and Qiao2018) have demonstrated the three-dimensional MOC’s utility in curved shock flowfield design, albeit limited to shock dependence domains. However, three-dimensional MOC suffers from computational inefficiency due to space intersection calculations and multiple iterations. Alternative approaches like the cross-stream marching method (Jones et al. Reference Jones, Sobieczky, Seebass and Dougherty1995) trade spatial complexity for directional derivative approximations, with the cross-stream marching plus method (Qiao et al. Reference Qiao, Yu and Wang2019) variant improving robustness through adaptive solution optimisation. Unfortunately, it is difficult to guarantee accuracy when dealing with complex cases.
Osculating methods provide an effective approach to overcome the spatial complexity inherent in three-dimensional MOC. These methods simplify the three-dimensional flowfield by discretising the flowfield into superimposed axisymmetric flowfield slices (Sobieczky et al. Reference Sobieczky, Dougherty and Jones1990). This discretisation strategy relies on a fundamental assumption: the flow remains confined within individual osculating planes, thereby neglecting transverse flow components and lateral pressure gradients between adjacent planes. While this assumption enables efficient flowfield solutions through independent planar computations, it inherently limits accuracy in scenarios where lateral flow effects become significant – particularly when analysing shockwaves with pronounced lateral intensity variations. To mitigate these limitations, subsequent studies have introduced pressure gradient corrections. Chauffour & Lewis (Reference Chauffour and Lewis2003, Reference Chauffour and Lewis2004) and Lewis & Chauffour (Reference Lewis and Chauffour2005) systematically validated the azimuthal pressure gradient assumption under specific Mach number regimes and implemented lateral pressure gradient compensations. However, their corrections failed to fully account for three-dimensional flow interactions. Building upon the conventional osculating cones methodology, Zheng et al. (Reference Zheng, Li, Zhu and You2020b ) developed a multiple osculating cones approach that partially addressed lateral pressure gradient effects. Nevertheless, this method retained the restrictive assumption of flow confinement within fixed osculating planes, leading to persistent errors in cases involving large lateral gradients. The subsequent local-turning osculating cones method (LTOCs) by Zheng et al. (Reference Zheng, Hu, Li, Zhu, You and Song2020a ) attempted to resolve three-dimensional flow characteristics through flow deflection approximations in osculating planes. While LTOCs reduced wall pressure distribution errors to approximately 2 %, it introduced computational inefficiencies. Based on this idea, Sun et al. (Reference Sun, Zheng and You2024) developed the stream-surface iteration-based flowfield calculation (SIFC) method, which calculates the flowfield according to the wall pressure distribution. More recently, Zhou et al. (Reference Zhou, Jin, Zhuang and Shi2022) proposed the space-streamline-based MOC (SMOC) method, which incorporates pressure gradient corrections during streamline calculations. By determining osculating planes through discrete pressure gradient analysis and solving them via axisymmetric MOC, SMOC demonstrates improved three-dimensional flowfields calculation. While this approach marks an algorithmic advancement, its validation remains limited to a single test case, insufficient for comprehensive error quantification or robust assessment of its applicability to complex configurations.
These methodological developments collectively highlight two critical unresolved challenges: (i) the need to trace how the shock geometry shapes the downstream flow physics, beyond merely computing the post-shock state, and (ii) the persistent trade-off between computational efficiency and physical fidelity in three-dimensional flowfield reconstruction. In this paper, the governing equations of the method of three-dimensional curved characteristics (MOTC) are derived from the steady Euler equations. Built upon the three-dimensional curved shock theory (3D-CST-boos), which provides the post-shock pressure gradients and streamline curvatures at the shock surface, MOTC advances these quantities along streamlines and characteristics to construct the complete downstream flowfield. The method retains the lateral pressure gradient and cross-flow information that conventional osculating methods discard, achieving errors below 0.5 %, while maintaining the computational efficiency required for parametric design studies. More importantly, by preserving the physical link between shock geometry and downstream flow evolution, MOTC enables a shock-based analysis of three-dimensional flowfield physics: the influence of shock curvature and axial ratio on compression performance, lateral flow redistribution and reflected shock structures can be systematically examined. This provides a rational basis for selecting three-dimensional shock shapes in forebody or inlet inverse design, where the designer must understand, a priori, how a chosen shock geometry will determine the downstream flow behaviour. Numerical validation across multiple benchmark cases confirms the methodology’s precision, particularly in resolving strong reflected shocks and lateral pressure gradients that are essential to three-dimensional shock-flowfield physics.
The remainder of this paper is organised as follows. Section 2 presents the theoretical foundation and computational framework of MOTC. Section 2.1.1 introduces the three-dimensional curved shock theory (3D-CST-boos), which supplies the post-shock pressure gradients and streamline curvatures at the shock surface, and § 2.1.2 derives the MOTC governing equations from the steady Euler equations in a streamline-oriented coordinate system. Section 2.2 describes the computational procedure, including the shock surface discretisation, the interior unit process and the shock unit process for reflected shock calculations. Section 3 validates the method against the Taylor–Maccoll equation for a conical shock and against inviscid CFD for external and internal elliptic curved shock flowfields. Section 4 applies MOTC to a parametric study of curved shock flowfields, analysing how the shock geometry – characterised by the transverse shock curvature and the axial ratio – governs the streamwise compression, the lateral pressure gradient distribution, and the reflected shock structure in both external and internal flows. Conclusions are drawn in § 5.
2. Method of three-dimensional curved characteristics
2.1. Introduction to 3D-CST-boos and derivation of MOTC equations
Osculating methods simply segment a three-dimensional flowfield into planes, discarding the lateral flow between slices from the starting point. When this lateral flow starts from a shock wave, the streamline follows a curved surface, as shown in figure 1. The three-dimensional curved shock theory (Zhang et al. Reference Zhang, Yang, Shi, Zhu, Zheng, You and Wang2025) establishes the relationship between the geometric parameters of the three-dimensional shock wave and the aerodynamic parameters before and after the shock with their gradients so that the lateral gradient parameters on the shock can be obtained. Based on these lateral gradient parameters, MOTC provides a flowfield that is more consistent with the actual flow than traditional osculating methods.
Schematic diagram of three-dimensional flow.
The three-dimensional curved shock theory provides the post-shock parameters and their gradients at the shock surface. To propagate these quantities into the downstream flowfield, governing equations along streamlines and characteristics are required. The following subsections present these two components in turn: § 2.1.1 summarises the 3D-CST-boos theory that supplies the shock-surface data and § 2.1.2 derives the MOTC equations that advance the solution downstream.
2.1.1. Three-dimensional curved shock theory based on orthogonal frame of shock surface
The 3D-CST-boos constructs a local orthogonal frame
$(\boldsymbol{e}_n, \boldsymbol{e}_m, \boldsymbol{e}_l)$
attached to the shock surface, as shown in figure 2. When a free incoming flow with velocity
$\boldsymbol{V}_1$
passes through a three-dimensional curved shock, the downstream velocity is
$\boldsymbol{V}_2$
. The upstream velocity is decomposed into a component
$V_{1m}$
tangent to the shock surface and a component
$V_{1n}$
normal to it. The unit vector
$\boldsymbol{e}_n$
is aligned with the shock normal,
$\boldsymbol{e}_m$
is taken along the direction of
$V_{1m}$
and
$\boldsymbol{e}_l$
completes the right-handed system, lying perpendicular to the
$\boldsymbol{e}_n$
–
$\boldsymbol{e}_m$
plane.
Schematic diagram of 3D-CST-boos (Zhang et al. Reference Zhang, Yang, Shi, Zhu, Zheng, You and Wang2025).
In this orthogonal frame, the steady Euler equations take the form
\begin{align} \left . \begin{aligned} \rho \frac {\partial V_{m}}{\partial m}+\rho \frac {\partial V_{n}}{\partial n}+V_{m}\frac {\partial \rho }{\partial m}+V_{n}\frac {\partial \rho }{\partial n}-\rho V_{n}k_{nm}-\rho V_{n}k_{nl}-\rho V_{m}k_{gl}=0, \\[4pt]\rho V_{n}\frac {\partial V_{l}}{\partial n}+\rho V_{m}^{2}k_{gm}-\rho V_{m}V_{n}\tau _{gm}+\frac {\partial p}{\partial l}=0, \\[4pt]\rho V_{m}\frac {\partial V_{m}}{\partial m}+\rho V_{n}\frac {\partial V_{m}}{\partial n}-\rho V_{n}V_{m}k_{nm}+\frac {\partial p}{\partial m}=0, \\[4pt]\rho V_{n}\frac {\partial V_{n}}{\partial n}+\rho V_{m}\frac {\partial V_{n}}{\partial m}+\rho V_{m}^{2}k_{nm}+\frac {\partial p}{\partial n}=0, \\[4pt]\frac {\gamma V_{m}}{\gamma -1}\left (\frac {\partial p}{\partial m}-\frac {p}{\rho }\frac {\partial \rho }{\partial m}\right )+\rho V_{m}^{2}\frac {\partial V_{m}}{\partial m}+\rho V_{m}V_{n}\frac {\partial V_{n}}{\partial m} \\[4pt]+\frac {\gamma V_{n}}{\gamma -1}\left (\frac {\partial p}{\partial n}-\frac {p}{\rho }\frac {\partial \rho }{\partial n}\right )+\rho V_{n}V_{m}\frac {\partial V_{m}}{\partial n}+\rho V_{n}^{2}\frac {\partial V_{n}}{\partial n}=0. \\\end{aligned} \right \} \end{align}
Here,
$\rho$
and
$p$
denote the local density and pressure,
$V_m$
,
$V_n$
and
$V_l$
are the velocity components along
$\boldsymbol{e}_m$
,
$\boldsymbol{e}_n$
and
$\boldsymbol{e}_l$
, and
$\gamma$
is the specific heat ratio. The quantities
$k_{nm}$
and
$k_{nl}$
denote the normal curvatures of the shock surface in the
$\boldsymbol{e}_m$
and
$\boldsymbol{e}_l$
directions, respectively, while
$k_{gm}$
and
$k_{gl}$
are the corresponding geodesic curvatures. The quantity
$\tau _{gm}$
arises from the motion of the frame along the surface. The subscripts
$m$
,
$n$
and
$l$
refer to components in the
$\boldsymbol{e}_m$
,
$\boldsymbol{e}_n$
and
$\boldsymbol{e}_l$
directions, respectively.
To close the system at the shock, the Euler equations must be supplemented by the shock-jump conditions. For a three-dimensional curved shock, not only the Rankine–Hugoniot relations for the zero-order quantities but also their directional derivatives along the shock surface are required. In the following, the subscripts 1 and 2 denote quantities immediately upstream and downstream of the shock, respectively. Differentiating the Rankine–Hugoniot relations in the
$\boldsymbol{e}_l$
and
$\boldsymbol{e}_m$
directions, and expressing them in the shock-based frame yields the following differential shock conditions:
\begin{align} \left . \begin{aligned} V_{1n} \frac {\partial \rho _1}{\partial l} + \frac {\partial V_{1n}}{\partial l} \rho _1 &= V_{2n} \frac {\partial \rho _2}{\partial l} + \frac {\partial V_{2n}}{\partial l} \rho _2, \\[4pt]V_{1n} \frac {\partial \rho _1}{\partial m} + \frac {\partial V_{1n}}{\partial m} \rho _1 &= V_{2n} \frac {\partial \rho _2}{\partial m} + \frac {\partial V_{2n}}{\partial m} \rho _2, \\[4pt]\frac {\partial p_1}{\partial l} + V_{1n}^2 \frac {\partial \rho _1}{\partial l} + 2\rho _1 V_{1n} \frac {\partial V_{1n}}{\partial l} &= \frac {\partial p_2}{\partial l} + V_{2n}^2 \frac {\partial \rho _2}{\partial l} + 2\rho _2 V_{2n} \frac {\partial V_{2n}}{\partial l}, \\[4pt]\frac {\partial p_1}{\partial m} + V_{1n}^2 \frac {\partial \rho _1}{\partial m} + 2\rho _1 V_{1n} \frac {\partial V_{1n}}{\partial m} &= \frac {\partial p_2}{\partial m} + V_{2n}^2 \frac {\partial \rho _2}{\partial m} + 2\rho _2 V_{2n} \frac {\partial V_{2n}}{\partial m}, \\[4pt]\frac {\partial V_{1m}}{\partial l} &= \frac {\partial V_{2m}}{\partial l}, \\[4pt]\frac {\partial V_{1m}}{\partial m} &= \frac {\partial V_{2m}}{\partial m}, \\[4pt]\frac {\gamma }{\gamma - 1} \left ( \frac {1}{\rho _1} \frac {\partial p_1}{\partial l} - \frac {p_1}{\rho _1^2} \frac {\partial \rho _1}{\partial l} \right ) + V_{1n} \frac {\partial V_{1n}}{\partial l} &= \frac {\gamma }{\gamma - 1} \left ( \frac {1}{\rho _2} \frac {\partial p_2}{\partial l} - \frac {p_2}{\rho _2^2} \frac {\partial \rho _2}{\partial l} \right ) + V_{2n} \frac {\partial V_{2n}}{\partial l}, \\[4pt]\frac {\gamma }{\gamma - 1} \left ( \frac {1}{\rho _1} \frac {\partial p_1}{\partial m} - \frac {p_1}{\rho _1^2} \frac {\partial \rho _1}{\partial m} \right ) + V_{1n} \frac {\partial V_{1n}}{\partial m} &= \frac {\gamma }{\gamma - 1} \left ( \frac {1}{\rho _2} \frac {\partial p_2}{\partial m} - \frac {p_2}{\rho _2^2} \frac {\partial \rho _2}{\partial m} \right ) + V_{2n} \frac {\partial V_{2n}}{\partial m}. \end{aligned} \right \} \end{align}
Equations (2.1) and (2.2) together form the complete system that governs the flow variables and their gradients across a three-dimensional curved shock.
To extract the post-shock gradients, the dimensionless pressure, density and velocity gradients are defined as
\begin{align} \left . \begin{aligned} & \frac {\partial p}{p \partial l} = P_{L}, \quad \frac {\partial p}{p \partial m} = P_{M}, \quad \frac {\partial p}{p \partial n} = P_{N},\\[4pt]& \frac {\partial \rho }{\rho \partial l} = R_{L}, \quad \frac {\partial \rho }{\rho \partial m} = R_{M}, \quad \frac {\partial \rho }{\rho \partial n} = R_{N}, \\[4pt]& \frac {\partial V_{m}}{V_{m}\partial l}=V_{ML}, \quad \frac {\partial V_{m}}{V_{m}\partial m}=V_{\textit{MM}},\quad \frac {\partial V_{n}}{V_{n}\partial l}=V_{NL}, \quad \frac {\partial V_{n}}{V_{n}\partial m}=V_{\textit{NM}}. \\ \end{aligned} \right \} \end{align}
Combining the Euler equations (2.1) with the differential shock conditions (2.2) and introducing the dimensionless variables (2.3) yields the 3D-CST-boos equations
\begin{align} \left . \begin{aligned} P_{2L} & = J_{1}R_{1L} + J_{2}V_{1NL} + J_{3}P_{1L}, \\[4pt]P_{2M} & = J_{1}R_{1M} + J_{2}V_{1NM} + J_{3}P_{1M}, \\[4pt]P_{2N} & = (L_{1}J_{1} + L_{2}K_{1})R_{1M} + (L_{1}J_{2} + L_{2}K_{3})V_{1NM} + (L_{1}J_{3} + L_{2}K_{2})P_{1M} \\[4pt]&\quad+ L_{2}V_{1MM} - L_{2}k_{gl} - L_{3}(k_{nl} + k_{nm}) - L_{4}k_{nm}. \end{aligned} \right \} \end{align}
The coefficients can be found in Appendix A. The detailed derivation of (2.1)–(2.4) can be found from Zhang et al. (Reference Zhang, Yang, Shi, Zhu, Zheng, You and Wang2025).
Equations (2.4) supply the post-shock pressure gradients
$P_{2L}$
,
$P_{2M}$
,
$P_{2N}$
in the shock-based frame as explicit functions of the shock curvature and the upstream flow non-uniformity.
However, the aerodynamic parameters and their gradients before or after the shock that the three-dimensional curved shock theory gives are only tenable near the shock. As the flow moves away from the shock, the flow parameters will also change accordingly. Therefore, the MOTC equations still need to be derived from the Euler equations applicable to the supersonic flow.
2.1.2. Derivation of MOTC equations
The continuity and momentum equations of Euler equations in differential form are
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{V}) = 0, \\[-28pt] \nonumber \end{align}
\begin{align} \rho (\boldsymbol{V} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{\cdot }\boldsymbol{V} + \boldsymbol{\nabla }p = 0, \\[0pt] \nonumber \end{align}
where
$p$
,
$\rho$
and
$\boldsymbol{V}$
represent the local static pressure, density and velocity vector, respectively. By expanding them along
$s,n,l$
(
$s$
is the direction of the streamline,
$n$
and
$l$
are any two directions perpendicular to the direction of
$s$
), equations can be obtained along the streamline and in the directions normal to the streamline (Appendix B)
\begin{align} \frac {\partial (\rho V)}{\partial s} + \rho V \frac {\partial \delta }{\partial n} + \rho V \frac {\partial \eta }{\partial l} = 0, \\[-28pt] \nonumber \end{align}
\begin{align} \rho V \frac {\partial V}{\partial s} + \frac {\partial p}{\partial s} = 0, \\[-28pt] \nonumber \end{align}
\begin{align} \rho V^2 \left ( \frac {\partial \delta }{\partial s} \boldsymbol{e}_n + \frac {\partial \eta }{\partial s} \boldsymbol{e}_l \right ) + \frac {\partial p}{\partial n} \boldsymbol{e}_n + \frac {\partial p}{\partial l} \boldsymbol{e}_l = 0 , \\[-12pt] \nonumber \end{align}
where
$\delta$
and
$\eta$
represent the deflection angles of the streamline in the
$n$
-direction and
$l$
-direction, respectively. Here,
$\boldsymbol{e}_s$
,
$\boldsymbol{e}_n$
and
$\boldsymbol{e}_l$
are the unit vectors in the respective directions. The energy equation along the streamline can be written as
\begin{align} \frac {\partial p}{\partial s}-\frac {\gamma p}{\rho }\frac {\partial \rho }{\partial s}=0. \end{align}
Then, the pressure gradients and streamline curvatures in the three directions are projected to the direction of the characteristics (represented by
$c$
)
\begin{align} \frac {\partial }{\partial c}=\frac {\partial }{\partial s}\frac {\partial s}{\partial c}+\frac {\partial }{\partial n}\frac {\partial n}{\partial c}+\frac {\partial }{\partial l}\frac {\partial n}{\partial c}. \end{align}
For the pressure and flow direction angles, the projection equations are
\begin{align} \frac {\partial p}{\partial c}=\frac {\partial p}{\partial s}\frac {\partial s}{\partial c}+\frac {\partial p}{\partial n}\frac {\partial n}{\partial c}+\frac {\partial p}{\partial l}\frac {\partial n}{\partial c}, \\[-28pt] \nonumber \end{align}
\begin{align} \frac {\partial \delta }{\partial c}=\frac {\partial \delta }{\partial s}\frac {\partial s}{\partial c}+\frac {\partial \delta }{\partial n}\frac {\partial n}{\partial c}+\frac {\partial \delta }{\partial l}\frac {\partial l}{\partial c}, \\[-28pt] \nonumber \end{align}
\begin{align} \frac {\partial \eta }{\partial c}=\frac {\partial \eta }{\partial s}\frac {\partial s}{\partial c}+\frac {\partial \eta }{\partial n}\frac {\partial n}{\partial c}+\frac {\partial \eta }{\partial l}\frac {\partial l}{\partial c}. \\[-12pt] \nonumber \end{align}
Define the dimensionless pressure gradients in each direction as
$P, D, L$
:
\begin{align} P=\frac {\partial p}{\rho V^{2}\partial s}, \\[-28pt] \nonumber \end{align}
\begin{align} D=\frac {\partial p}{\rho V^{2}\partial n}, \\[-28pt] \nonumber \end{align}
\begin{align} L=\frac {\partial p}{\rho V^{2}\partial l}. \\[-12pt] \nonumber \end{align}
At the shock surface, the initial values of
$P$
,
$D$
,
$L$
are obtained from the 3D-CST-boos gradients
$P_{2L}$
,
$P_{2M}$
,
$P_{2N}$
derived in § 2.1.1. In the
$(s, n, l)$
system,
$s$
is aligned with the post-shock velocity
$\boldsymbol{V}_2$
,
$n$
is perpendicular to
$s$
within the
$\boldsymbol{e}_m$
–
$\boldsymbol{e}_n$
plane and
$l$
coincides with
$\boldsymbol{e}_l$
. Let
$\delta$
denote the flow deflection angle behind the shock, i.e. the angle between
$\boldsymbol{V}_2$
and
$\boldsymbol{e}_m$
. The directional derivatives in the two frames are related by the chain rule:
\begin{align} \frac {\partial }{\partial s} = \cos \delta \,\frac {\partial }{\partial m} + \sin \delta \,\frac {\partial }{\partial n}, \quad \frac {\partial }{\partial n} = -\sin \delta \,\frac {\partial }{\partial m} + \cos \delta \,\frac {\partial }{\partial n}, \quad \frac {\partial }{\partial l} = \frac {\partial }{\partial l}. \end{align}
Applying these relations to the pressure gradient and converting the dimensionless forms yields
\begin{align} \left . \begin{aligned} P &= \frac {1}{\rho V^{2}}\frac {\partial p}{\partial s} = \frac {p}{\rho V^{2}}\big (\cos \delta \,P_{2M} + \sin \delta \,P_{2N}\big ), \\[4pt]D &= \frac {1}{\rho V^{2}}\frac {\partial p}{\partial n} = \frac {p}{\rho V^{2}}\big (-\sin \delta \,P_{2M} + \cos \delta \,P_{2N}\big ), \\[4pt]L &= \frac {1}{\rho V^{2}}\frac {\partial p}{\partial l} = \frac {p}{\rho V^{2}}\,P_{2L}, \end{aligned} \right \} \end{align}
where
$p$
,
$\rho$
and
$V$
are the post-shock values, and
$p/(\rho V^{2}) = 1/(\gamma M^{2})$
, with
$M$
the post-shock Mach number.
According to (2.7), (2.8) and (2.10),
\begin{align} \frac {\partial \delta }{\partial n} + \frac {\partial \eta }{\partial l} = - \frac {\partial (\rho V)}{\rho V \partial s} = \big(1 - M^2 \big) P. \end{align}
Replace
${\partial \eta }/{\partial l}$
in (2.14) to reduce the unknowns:
\begin{align} \frac {\partial \eta }{\partial c}=-L\frac {\partial s}{\partial c}+\frac {\partial \eta }{\partial n}\frac {\partial n}{\partial c}+ \left[ \big(1-M^{2} \big)P-\frac {\partial \delta }{\partial n} \right]\frac {\partial l}{\partial c}. \end{align}
After sorting out the above-mentioned equations, the equations of partial differential equations of aerodynamic parameter
$p$
, and flow direction angles
$\delta$
and
$\eta$
along the streamline and characteristics can be obtained:
\begin{align} & \frac {\partial \delta }{\partial c}=-D\frac {\partial s}{\partial c}+\frac {\partial \delta }{\partial n}\frac {\partial n}{\partial c}+\frac {\partial \delta }{\partial l}\frac {\partial l}{\partial c}, \nonumber \\[4pt]& \frac {\partial \eta }{\partial c}=-L\frac {\partial s}{\partial c}+\frac {\partial \eta }{\partial n}\frac {\partial n}{\partial c}+ \left[\big(1-M^{2} \big)P-\frac {\partial \delta }{\partial n} \right]\frac {\partial l}{\partial c}, \nonumber \\[4pt]& \frac {\partial p}{\partial c}=\rho V^2\left (P \frac {\partial s}{\partial c}+D \frac {\partial n}{\partial c}+L \frac {\partial l}{\partial c}\right ) \!, \nonumber \\[4pt]& \frac {\partial \delta }{\partial s}=-D, \nonumber \\[4pt]& \frac {\partial \eta }{\partial s}=-L, \nonumber \\[4pt]& \frac {\partial p}{\partial s}=\rho V\frac {\partial V}{\partial s}, \nonumber \\[4pt]& \frac {\partial p}{\partial s}=\rho V^{2}P, \nonumber \\[4pt]& \frac {\partial p}{\partial s}=\frac {\gamma p}{\rho }\frac {\partial \rho }{\partial s}. \end{align}
The partial differential equations (2.22) of aerodynamic parameters along the streamline and characteristics are transformed into the total differential equations along the streamline or characteristics, and finally the governing equations suitable for the three-dimensional curved characteristics are obtained:
\begin{align} & \frac {\mathrm{d}\delta }{\mathrm{d}c} = -D(\boldsymbol{e}_{s} \boldsymbol{\cdot }\boldsymbol{e}_{c}) + \frac {\partial \delta }{\partial n}(\boldsymbol{e}_{n} \boldsymbol{\cdot }\boldsymbol{e}_{c}) + \frac {\partial \delta }{\partial l}(\boldsymbol{e}_{l} \boldsymbol{\cdot }\boldsymbol{e}_{c}), \nonumber \\[4pt]& \frac {\mathrm{d}\eta }{\mathrm{d}c} = -L(\boldsymbol{e}_{s} \boldsymbol{\cdot }\boldsymbol{e}_{c}) + \frac {\partial \eta }{\partial n}(\boldsymbol{e}_{n} \boldsymbol{\cdot }\boldsymbol{e}_{c}) + \left[ \big(1 - M^{2} \big)P - \frac {\partial \delta }{\partial n} \right](\boldsymbol{e}_{l} \boldsymbol{\cdot }\boldsymbol{e}_{c}) ,\nonumber \\[4pt]& \frac {\mathrm{d} p}{\mathrm{d}c} = \rho V^{2} \big[P(\boldsymbol{e}_{s} \boldsymbol{\cdot }\boldsymbol{e}_{c}) + D(\boldsymbol{e}_{n} \boldsymbol{\cdot }\boldsymbol{e}_{c}) + L(\boldsymbol{e}_{l} \boldsymbol{\cdot }\boldsymbol{e}_{c}) \big] ,\nonumber \\[4pt]& \frac {\mathrm{d}\delta }{\mathrm{d}s} = -D ,\nonumber \\[4pt]& \frac {\mathrm{d}\eta }{\mathrm{d}s} = -L, \nonumber \\[4pt]& \frac {\mathrm{d }p}{\mathrm{d}s} = -\rho V \frac {\mathrm{d} V}{\mathrm{d} s}, \nonumber \\[4pt]& \frac {\mathrm{d} p}{\mathrm{d}s} = \rho V^{2} P, \nonumber \\[4pt]& \frac {\mathrm{d} p}{\mathrm{d}s} = \frac {\gamma p}{\rho } \frac {\mathrm{d}\rho }{\mathrm{d} s}. \end{align}
From the ordinary differential equations (2.23), it can be seen that in a unit process, there are eleven unknowns, including
$p, \rho , V, P, D, L, \delta , \eta , ({\partial \delta }/{\partial n}), ({\partial \delta }/{\partial l})$
and
$({\partial \eta }/{\partial n})$
. There are five equations along the streamlines and three along the characteristics. Thus, at least two curved characteristics are required to solve a unit process, which means two sets of equations for curved characteristics. MOTC employs two sets of characteristics equations to facilitate computation: one for the curved characteristics within the flow surface and the other obtained by subtracting the characteristics equations on both sides of the flow surface. This extends the application of the three-dimensional curved shock theory by MOTC from the vicinity of the shock wave to the post-shock flowfield. Moreover, unlike the 3-D-MOC, which requires four characteristics and one streamline for iteration, solving the MOTC equations needs at most three characteristics and one streamline. There is no need to iteratively solve for the coordinates and parameters of the fifth line, as MOTC can directly calculate the solution point from the information of four nodes. Additionally, MOTC does away with the need for complex spatial surface intersections. These advantages make MOTC a simple and time-saving method for flowfield computation. Furthermore, MOTC establishes a relationship between three-dimensional pressure gradients, streamline curvatures and flow parameters. After computation, pressure and flow deflection angles, as well as their gradients, can be obtained across the entire flowfield. Pressure gradients
$P$
,
$D$
,
$L$
are important variables in flow analysis, describing the variations in flow variables. Specifically, according to (2.20) and (2.9),
$P$
represents the convergence or divergence of the flow, and
$D$
and
$L$
represent the curvature of the flow. With information representing these patterns of change, a flowfield can be solved with fewer grids while maintaining high accuracy.
In the solution process, because the relationship between aerodynamic parameters and pressure gradients is established, the pressure gradient parameters
$P, D, L$
in the flowfield can be solved synchronously for more specific flow analysis.
At the same time, the two streamline curvatures are also non-negligible in three-dimensional flows. In (2.20),
$({\partial \delta }/{\partial n})$
and
$({\partial \eta }/{\partial l})$
represent the convergence and divergence of the flow in two directions. Since the vorticity equation in streamline coordinates is defined as
\begin{align} \omega _{s} = V \frac {\partial \eta }{\partial n} - V \frac {\partial \delta }{\partial l}, \\[-28pt] \nonumber \end{align}
\begin{align} \omega _{n} = \frac {\partial V}{\partial l} - V \frac {\partial \eta }{\partial s}, \\[-28pt] \nonumber \end{align}
\begin{align} \omega _{l} = V \frac {\partial \delta }{\partial s} - \frac {\partial V}{\partial n}. \\[-12pt] \nonumber \end{align}
Therefore,
$({\partial \delta }/{\partial l})$
and
$({\partial \eta }/{\partial n})$
represent the rotation of the streamlines. Hereafter,
$\varGamma$
will be used to denote the normalised vorticity:
\begin{align} \varGamma _{s} = \frac {\partial \delta }{\partial l} - \frac {\partial \eta }{\partial n}, \\[-28pt] \nonumber \end{align}
\begin{align} \varGamma _{n} = \frac {\partial V}{V\partial l} - L, \\[-28pt] \nonumber \end{align}
\begin{align} \varGamma _{l} = D - \frac {\partial V}{V\partial n}, \\[-12pt] \nonumber \end{align}
where
$({\partial V}/{V\partial l})$
and
$({\partial V}/{V\partial n})$
are the gradients of the velocity field perpendicular to the direction of the streamlines. The above-mentioned difference equation does not calculate this term. If the vorticity of the entire flowfield needs to be obtained, it needs to be calculated through a separate formula:
\begin{align} \begin{aligned} \frac {{\rm d}V}{{\rm d}c}=V\left [\frac {\partial V}{V\partial s}(\boldsymbol{e}_{s}\boldsymbol{\cdot }\boldsymbol{e}_{c})+\frac {\partial V}{V\partial n}(\boldsymbol{e}_{n}\boldsymbol{\cdot }\boldsymbol{e}_{c})+\frac {\partial V}{V\partial l}(\boldsymbol{e}_{l}\boldsymbol{\cdot }\boldsymbol{e}_{c})\right ] \\[5pt]=V\left [-P(\boldsymbol{e}_{s}\boldsymbol{\cdot }\boldsymbol{e}_{c})+\frac {\partial V}{V\partial n}(\boldsymbol{e}_{n}\boldsymbol{\cdot }\boldsymbol{e}_{c})+\frac {\partial V}{V\partial l}(\boldsymbol{e}_{l}\boldsymbol{\cdot }\boldsymbol{e}_{c})\right ]\! .\end{aligned} \end{align}
By solving (2.30) along two characteristic lines, the velocity gradients required for calculating vorticity can be obtained, and thus, the vorticity of the entire flowfield can be determined. These gradients can also be used to study physical phenomena such as shock wave reflection.
2.2. Algorithms
With the post-shock parameters and their gradients determined by 3D-CST-boos along each discrete shock curve, the downstream flowfield is constructed by marching along streamlines and characteristics. The computation proceeds through two types of unit processes: the interior process, which advances the solution within the smooth region of the flowfield, and the shock process, which incorporates the reflected shock in internal flow calculations.
Before the marching procedure, the shock surface is discretised into a family of stream surfaces following the method of Zheng et al. (Reference Zheng, Hu, Li, Zhu, You and Song2020a
), as illustrated in figure 3. The local osculating plane is determined by the free stream vector and the surface normal vector. Starting from an initial shock point
$A_{1}$
, the intersection of the shock surface with its local osculating plane
$P_{1}$
yields the next discrete shock point
$A_{2}$
. Repeating this procedure produces a family of discrete shock curves
$A_{n}{-}A_{3}$
that serve as the computational boundaries for MOTC.
Determination of the shock curve in the stream surface (Zheng et al. Reference Zheng, Hu, Li, Zhu, You and Song2020a ).
Along each discrete shock curve, the shock angle
$\theta$
is obtained from the dot product of the incoming flow velocity vector
$\boldsymbol{V}$
and the shock surface unit normal vector
$\boldsymbol{n}$
:
\begin{align} \theta = \arcsin \left ( \frac {\boldsymbol{V} \boldsymbol{\cdot }\boldsymbol{n}}{|\boldsymbol{V}|} \right ) \!. \end{align}
The zero-order post-shock parameters then follow from the standard Rankine–Hugoniot relations. The first-order gradients (
$P$
,
$D$
,
$L$
and the streamline curvatures) are supplied by 3D-CST-boos, as described in § 2.1.1. With these initial data established along each discrete shock curve, the interior and shock unit processes are applied as follows.
2.2.1. Interior unit process
The interior process computes the flow parameters at an unknown point from the known data on the shock curve and the preceding streamlines. Figure 4(a) shows the geometric configuration, where the red, blue and black curves represent the shock curve, streamline and characteristics, respectively. The upstream free stream velocity
$\boldsymbol{V}_{\infty }$
and the post-shock velocity
$\boldsymbol{V}_2$
at point 2 are indicated to clarify the flow direction. The step-by-step procedure, summarised in the flowchart of figure 4(b), is as follows.
Interior process and flowchart of MOTC. (a) Unit interior process of MOTC. (b) Flowchart of MOTC interior process.
Step 1: Obtain post-shock parameters at the shock boundary. For two adjacent shock points 1 and 2 (figure 4 a), the post-shock pressure, density, velocity and flow direction are determined from the Rankine–Hugoniot relations using the local shock angle given by (2.31).
Step 2: Obtain gradients of post-shock parameters. The dimensionless pressure gradients
$P$
,
$D$
,
$L$
and the streamline curvatures
$({\partial \delta }/{\partial n})$
,
$ ({\partial \delta }/{\partial l})$
,
$ ({\partial \eta }/{\partial n})$
at points 1 and 2 are supplied by 3D-CST-boos, which relates these gradients to the shock curvature and the upstream flow conditions.
Step 3: Determine the position of the unknown point. The unknown point 3 lies at the intersection of the streamline emanating from point 2 and the characteristic running from point 1. Its position is obtained from the geometric vector equations:
\begin{align} t_{1}\boldsymbol{e}_{c1} - t_{2}\boldsymbol{e}_{s2} = \boldsymbol{P}_{2} - \boldsymbol{P}_{1}, \\[-28pt] \nonumber \end{align}
\begin{align} \boldsymbol{P}_{3} = \boldsymbol{P}_{2} + t_{2}\boldsymbol{e}_{s2}, \\[0pt] \nonumber \end{align}
where
$\boldsymbol{P}$
denotes the coordinate vector,
$\boldsymbol{e}_{c}$
is the characteristic direction,
$\boldsymbol{e}_{s}$
is the streamline direction, and
$t_{1}$
,
$t_{2}$
are the length coefficients along the respective directions.
Step 4: Provide initial estimates of the parameters at point 3. Using the gradients known at point 2, a Taylor expansion along the streamline gives the initial values of the flow deflection angles and pressure at point 3:
\begin{align} \delta _{3} = \delta _{2} + (s_{3} - s_{2})\,D_{2}, \\[-28pt] \nonumber \end{align}
\begin{align} \eta _{3} = \eta _{2} + (s_{3} - s_{2})\,L_{2}, \\[-28pt] \nonumber \end{align}
\begin{align} p_{3} = p_{2} + (s_{3} - s_{2})\,\rho _{2}V_{2}^{2}\,P_{2}. \\[-12pt] \nonumber \end{align}
Step 5: Determine the auxiliary points on the initial-value line. The intersections of the curved inverse Mach cone from point 3 with the initial-value line connecting points 4, 2 and 5 define the positions and parameters of the auxiliary points 4 and 5 (figure 4 a).
Step 6: Solve for the parameters and gradients at point 3. The flow parameters at point 3 are updated by integrating the MOTC governing equations (2.23) along the characteristics 4–3, 5–3, 1–3 and the streamline 2–3. This yields the pressure, flow deflection angles and the gradients
$P$
,
$D$
,
$L$
at point 3.
Step 7: Check convergence and iterate if necessary. If the difference between the updated and previous values of the aerodynamic parameters exceeds a prescribed tolerance, the position of point 3 is refreshed using the new parameter values, and Steps 5–6 are repeated. The positions and parameters of points 4 and 5 are updated accordingly in each iteration. Once convergence is achieved, point 3 provides the final parameters and their gradients, which serve as input data for the next downstream marching step.
Compared with the conventional interior MOC process, the present method offers two advantages. First, the Taylor expansions in Step 4 incorporate gradient information that is absent in MOC, yielding initial estimates closer to the converged values and reducing the number of iterations. Second, MOTC eliminates the need for complex surface intersection or bivariate interpolation.
2.2.2. Shock unit process
In internal flow inverse design, the flowfield contains a reflected shock, which requires a different unit process. Figure 5(a) shows the configuration, where the reflected shock, streamlines and characteristics are drawn in red, blue and cyan, respectively; the initial-value lines (black) are the terminal characteristics of the shock-dependent domain. The velocity
$\boldsymbol{V}_2$
at point 2 (upstream of the reflected shock) and the post-shock velocity
$\boldsymbol{V}_{\textit{post}}$
at point 3 are indicated to clarify the flow direction. The procedure, summarised in the flowchart of figure 5(b), is as follows.
Shock process and flowchart of MOTC. (a) Unit shock process of MOTC. (b) Flowchart of MOTC shock process.
Step 1: Initialise the known points. Points 1 and 2 lie on the incident shock or on the preceding characteristics; their parameters and gradients are known from the upstream interior computation. To close the shock process, additional shock conditions are required, which may include the post-shock flow deflection angle, the reflected shock angle or their variation rates along the shock surface, depending on the specific inverse design problem.
Step 2: Determine the position and curvature of the reflected-shock point. The position of point 3 on the reflected shock and the local shock curvature at that point are obtained using the secant method. The geometric relation between points 1, 2 and 3 is
\begin{align} t_{1}\boldsymbol{e}_{s1} - t_{2}\boldsymbol{e}_{s2} = \boldsymbol{P}_{2} - \boldsymbol{P}_{1}, \\[-28pt] \nonumber \end{align}
\begin{align} \boldsymbol{P}_{3} = \boldsymbol{P}_{2} + t_{2}\boldsymbol{e}_{s2}, \\[0pt] \nonumber \end{align}
where
$\boldsymbol{e}_{s1}$
is the direction of the reflected shock curve,
$\boldsymbol{e}_{s2}$
is the streamline direction at point 2,
$\boldsymbol{P}$
denotes the coordinate vector, and
$t_{1}$
,
$t_{2}$
are the length coefficients.
Step 3: Provide initial estimates of the parameters at point 3. Using the post-shock conditions prescribed in Step 1, the R-H relations are applied to obtain the initial estimates of the pre-shock parameters immediately upstream of the reflected shock at point 3.
Step 4: Solve for the parameters and gradients at point 3. The parameters and gradients at point 3 are updated by applying the MOTC governing equations (2.23) along the shock curve, incorporating the known data from points 1 and 2, and the initial estimates from Step 3.
Step 5: Check convergence and iterate if necessary. If the difference between the updated and previous values of the aerodynamic parameters exceeds a prescribed tolerance, the position and curvature of point 3 are refreshed using the secant method, and Steps 3–4 are repeated until convergence is achieved.
Step 6: Recover the remaining gradients. After convergence, the streamline curvatures
$ ({\partial \delta }/{\partial n})$
,
$ ({\partial \delta }/{\partial l})$
and
$ ({\partial \eta }/{\partial n})$
at point 3, which are not directly obtained in the shock process, are recovered from the inverse characteristics equations in (2.23) on both sides of point 3. This yields the complete set of parameters and gradients at point 3, which serve as the starting data for the subsequent interior marching.
For all internal flow calculations presented in this paper, the shock conditions are specified as follows: the flow deflection angle behind the reflected shock is set to zero, aligning the post-shock flow with the centre-body surface, and the variation rate of the reflected shock angle along the shock surface is set to zero.
The MOTC shock process does not involve lateral information during iteration, so it is essentially a curved two-dimensional calculation involving only simple geometric operations. The use of gradient information provides initial values for point 3 that are closer to the converged solution, reducing the number of iterations required. The remaining lateral gradient information is retrieved only after convergence in Step 6.
3. Verification
In this paper, three cases are selected to verify the accuracy of the MOTC: the conical shock flowfield and the external and internal elliptic curved shock flowfield. The flowfield behind the conical shock is verified by the Taylor–Maccoll equation, and the external and internal elliptic curved shock is verified by inviscid CFD simulation.
3.1. T-M flow
In the first case, the uniform incoming flow is given, the Mach number is 6, the half-cone angle of the conical shock is 26.57
$^\circ$
, the cone vertex position is (–0.1,0,0) and the axis is the
$x$
-axis. The numerical solution of the conical flow can be obtained through Taylor–Maccoll equations, so this case is used to verify the accuracy of the method under axisymmetric conditions.
Figure 6 shows the comparison diagram of pressure ratio distribution calculated by the two methods under the same conditions. The coloured areas and red lines are the results of the T-M equation, and the black dash-dot-dotted lines are the result of MOTC’s solution. The two have a strong consistency, indicating that MOTC can calculate axisymmetric flowfields with high precision. Specifically, the maximum error between the two is no more than 0.014 %.
Flowfields calculated using T-M equations and MOTC in uniform conical flow.
In all cases presented in this paper, spatial coordinates are non-dimensionalised by a reference length
$S$
, defined as the total length of the prescribed shock surface in the
$x$
-direction. Let
$x$
,
$y$
,
$z$
denote the physical coordinates; the dimensionless coordinates are then
$x/S$
,
$y/S$
,
$z/S$
. In each case,
$S$
is set to unity (
$S = 1$
), so the numerical values of the dimensionless coordinates coincide with those of the physical coordinates. Unless otherwise stated, all coordinates appearing in the subsequent figures refer to the dimensionless forms
$x/S$
,
$y/S$
,
$z/S$
.
3.2. External flow
The second case has the same incoming flow condition as Case 1. The shock surface function is an elliptic curved cone, and the expression is in the form of a parametric equation:
\begin{align} r(u, v) = \begin{cases} x = v ,\\[4pt] y = 0.2(v+1)^{2}\cos (u), \\[4pt] z = 0.15(v+1)^{2}\sin (u). \end{cases} \end{align}
Figure 7 shows the external elliptic curved shock flowfield calculated according to this equation. It can be seen from the figure that in addition to the change of lateral curvature, it also includes the change of longitudinal curvature, which has a stronger three-dimensional effect. The axial ratio
$a/b$
is 1.33. Therefore, the flowfield behind the shock exhibits a large pressure ratio variation in both the lateral and longitudinal directions. When traditional osculating methods or even improved LTOCs are used, there will be a large error.
Schematic of external elliptic curved shock flowfield.
The shock generator is simulated by inviscid CFD and the flowfield is compared with the MOTC. Figure 8(a) shows the comparison of contours of pressure ratios obtained by MOTC and inviscid CFD in the symmetric plane. The coloured areas and red solid lines represent the inviscid CFD flowfield. The black dashed curves in figure 8(a) represent the flowfield obtained by MOTC. It can be seen from the figure that the shock curves of the two are completely consistent and the isobaric lines inside the flowfield are also completely coincident. Due to the shock area being limited by the grid, CFD cannot accurately simulate and generate some numerical jitter. The
$x/S = 0.4$
and
$x/S = 0.7$
sections are sliced for comparison. The comparison results of the pressure ratio contours are shown in figure 8(b), where the coloured areas and red solid lines are the result of CFD simulation, and the black dashed curves are the results of MOTC calculation. It can be seen from the figure that both the upstream and downstream of the flowfield, the shock curve and the isobaric line all coincide, and there is no systematic deviation.
Comparison of CFD and MOTC non-dimensional pressure contours of external elliptic curved shock flowfield. (a) Symmetric plane. (b)
$x/S=0.4$
and
$x/S=0.7$
plane.
Detailed pressure distributions of the wall along the
$x$
-direction and
$y$
-direction are important to double-check the accuracy of the MOTC. Figure 9(a) shows the non-dimensional wall pressure distribution on three
$x$
-direction curves (
$y/S = 0$
,
$y/S = 0.15$
and
$y/S = 0.3$
). In the
$x$
-direction, the pressure obtained by the MOTC agrees well with the pressure gained by the inviscid CFD simulation. The pressure ratio distribution curves of the two have no systematic deviation and the maximum error is not more than 0.5 %. Figure 9(b) demonstrates the non-dimensional wall pressure distribution on three
$y$
-direction curves (
$x/S = 0.2$
,
$x/S = 0.4$
and
$x/S = 0.6$
). In the
$y$
-direction curves, the pressure obtained by the MOTC is nearly identical to the pressure gained by the inviscid CFD simulation and the maximum error is not more than 0.5
Comparison of CFD and MOTC non-dimensional wall pressure of external flow: (a)
$x$
direction; (b)
$y$
direction.
The flowfield structure obtained by MOTC is not significantly different from that obtained by inviscid CFD simulation regardless of longitudinal or lateral comparison, indicating that MOTC can accurately determine the flowfield information after the elliptic curved shocks with both longitudinal and lateral curvature changes. Therefore, MOTC can meet the accuracy requirements of supersonic external flowfield calculation and aircraft aerodynamic profile inverse design.
3.3. Internal flow
Due to the three-dimensional reflected shock and the strong lateral pressure gradient, the internal flowfield is more complex than the external flowfield. This is also why the osculating methods find it difficult to calculate the three-dimensional internal flow. This case uses the same flow conditions as the above-mentioned cases to verify the applicability of the MOTC method in calculating the internal flowfield. The equations of the incident shock surface are
\begin{align} r(u, v) = \begin{cases} x = v, \\[4pt]y = 1.2 \big(-0.1v^{2} - 0.3v + 1 \big)\cos (u), \\[4pt]z = \big(-0.1v^{2} - 0.3v + 1 \big)\sin (u) .\end{cases} \end{align}
Figure 10 shows the internal elliptic curved shock flowfield calculated by MOTC according to this equation. It can be seen from the figure that in addition to the change of lateral curvature, it also includes the change of longitudinal curvature, which has a stronger three-dimensional effect. The axial ratio
$a/b$
reaches 1.2. Therefore, the flowfield behind the shock exhibits a large pressure ratio variation in both the lateral and longitudinal directions. To avoid the Mach reflection, a centre-body is placed in this flowfield.
Schematic of internal elliptic curved shock flowfield.
The shock generator is simulated by inviscid CFD and the flowfield is compared with the MOTC. Figure 11(a) shows the comparison of contours of pressure ratios obtained by MOTC and CFD in the symmetric plane. The coloured areas and red solid lines represent the CFD flowfield. The black dashed curves represent the flowfield obtained using MOTC. From the figure, the curves of the two basically coincide. The
$x/S=1.0$
section is captured for comparison. The comparison results of the pressure ratio contours are shown in figure 11(b), where the coloured areas and red solid lines represent the CFD flowfield. The black dashed curves are the result of MOTC calculation. It can be seen from the figure that the pressure variation trend of the two is essentially the same and the resulting shock wave shapes are nearly identical. Due to the shock wave area being limited by the grid, CFD cannot accurately simulate and generate some numerical jitter. This phenomenon is particularly obvious at the reflected point.
Comparison of CFD and MOTC non-dimensional pressure contours of internal elliptic curved shock flowfield: (a) symmetric plane; (b)
$x/S=1.0$
plane.
Detailed pressure distributions of the wall along the
$x$
-direction and
$y$
-direction are important to double-check the accuracy of the MOTC. Figure 12(a) shows the non-dimensional wall pressure distribution on three
$x$
-direction curves (
$y/S = 0$
,
$y/S = 0.4$
and
$y/S = 0.8$
). In the
$x$
-direction, the pressure obtained by the MOTC agrees well with the pressure gained by the inviscid CFD simulation. The pressure ratio distribution curves of the two have no systematic deviation and the maximum error is not more than 0.5 %. Figure 12(b) demonstrates the non-dimensional wall pressure distribution on three
$y$
-direction curves (
$x/S = 0.4$
,
$x/S = 0.8$
and
$x/S = 1.1$
). It should be noted that
$x/S = 0.4$
is in the shock dependence domain, and
$x/S = 0.8$
and
$x/S = 1.1$
are in the isentropic compression region. In the
$y$
-direction curves, the pressure obtained by the MOTC is nearly identical to the pressure gained by the inviscid CFD simulation and the maximum error is not more than 0.5 %. Therefore, MOTC can accurately calculate the flowfield in the shock-dependent region and isentropic compression region.
Comparison of CFD and MOTC non-dimensional wall pressure of internal flow: (a)
$x$
direction; (b)
$y$
direction.
The above-mentioned comparison shows that MOTC can accurately calculate the flowfield, including the reflected shock, even in internal flows where the lateral pressure has a significant influence. MOTC can play an important role in internal flowfield analysis and aerodynamic design of internal compression basic flowfields.
4. Flowfield analysis with first-order parameters
Three-dimensional flows introduce significant spanwise variations compared with two-dimensional flows. When the flow around a streamline is not uniform, it can impact the overall flow. The aerodynamic design of hypersonic vehicles often employs non-axisymmetric shock waves and the effects of these shock waves must be considered when they are used. This paper alters the geometry of the basic shock surface by varying the coefficients of a parametric curved shock, resulting in different basic flowfields for both internal and external flows. The first-order parameters of the flowfield can be analysed. When using shocks for flowfield inverse design, it is best to know the features of the chosen shock waves and what advantages and disadvantages they will bring. This is of great significance in guiding the design. The benefit of the precision and efficiency of the MOTC method is that the flow features of the selected shock can be solved quickly and accurately. The shock that best meets the design requirements can then be selected based on the results obtained. This article selects only some parameters for analysis.
4.1. Three-dimensional features of supersonic external flow
The present analysis considers a family of elliptic curved shocks defined by the parametric equation
\begin{align} r(u, v) = \begin{cases} x = v, \\[4pt]y = \big(0.1v^{2} + \alpha v + 0.1 \big)\cos (u), \\[4pt]z = \big(0.1v^{2} + \beta v + 0.1 \big)\sin (u), \end{cases} \end{align}
in which the coefficients
$\alpha$
and
$\beta$
govern the expansion rate of the shock surface in the
$y$
- and
$z$
-directions, respectively. At the initial station (
$v=0$
), the shock cross-section is a circle of radius 0.1; downstream, the difference between
$\alpha$
and
$\beta$
causes the cross-section to evolve from circular into elliptical cross-section shapes.
According to curved shock theory (Mölder Reference Mölder2016), the transverse shock curvature
$S_b$
quantifies the spanwise bending of the shock surface and plays a central role in determining the post-shock pressure gradients. For the parametric surface defined by (4.1), the elliptic cross-section at any station
$v$
has semi-axes
$a(v)=0.1v^{2}+\alpha v+0.1$
in the
$y$
-direction and
$b(v)=0.1v^{2}+\beta v+0.1$
in the
$z$
-direction. In the symmetry plane
$y/S=0$
, the transverse shock curvature is given by
$S_b = (b/a^{2})(-\cos \theta )$
, where
$\theta$
is the local shock angle. Since
$\beta$
is the same for all subsequent cases considered, the shock angle
$\theta$
in this symmetry plane is identical and the variation of
$S_b$
across different configurations is determined solely by the geometric factor
$b/a^{2}$
, which depends on
$\alpha$
.
To illustrate the three-dimensional features captured by MOTC, three external curved shocks with distinct lateral shapes are selected, as shown in figure 13. In all three cases,
$\beta$
is fixed at 0.3, so the shock profiles in the symmetry plane
$y/S=0$
are identical. The differences among the three shocks arise solely from variations in
$\alpha$
, which alter the spanwise expansion rate, the cross-sectional ellipticity and, consequently, the transverse shock curvature
$S_b$
. The black line corresponds to the axisymmetric shock with
$\alpha =\beta =0.3$
. Increasing
$\alpha$
to 0.4 (orange line) produces an elliptic shock whose major axis aligns with the
$y$
-direction, while reducing
$\alpha$
to 0.2 (blue line) yields an elliptic shock with its minor axis along the
$y$
-direction.
External shock configurations with different
$\alpha$
values (
$\beta =0.3$
for all cases), front view along the
$x$
-direction.
Figure 14 compares the pressure and the pressure gradient along the streamline for the three shocks at
$y/S = 0$
, where the solid line represents the streamwise pressure gradient
$(P)$
and the dashed line represents the pressure ratio
$p/p_{0}$
. Because all three shocks share the same
$\beta$
, both the local shock angle
$\theta$
and the shock curvature in the
$y/S=0$
plane are identical; consequently, the initial pressure and the initial streamwise pressure gradient immediately behind the shock at
$x/S = 0$
are the same for all three cases. The differences that develop downstream are introduced solely through
$\alpha$
, which determines the spanwise expansion rate and thereby governs the streamwise evolution of the transverse shock curvature
$S_{b}$
.
Pressure and pressure gradient along the streamline at
$y/S = 0$
of external flow.
In the near-shock region, the flow is dominated by shock-induced compression. The transverse shock curvature in the symmetry plane
$y/S=0$
is
$S_b = (b/a^{2})(-\cos \theta )$
. For the external flow,
$\theta \in [0,\pi /2]$
and
$S_b$
is negative. Since
$\beta$
is fixed,
$\theta$
is the same for all three cases, so the variation of
$S_b$
is determined by the geometric factor
$b/a^{2}$
. A larger
$\alpha$
produces a larger semi-axis
$a$
, hence, a smaller
$b/a^{2}$
; because
$-\cos \theta \lt 0$
, a smaller
$b/a^{2}$
corresponds to a larger algebraic value of
$S_b$
. Consequently, the algebraic value of
$S_b$
is largest for
$\alpha =0.4$
and smallest for
$\alpha =0.2$
. According to curved shock theory, a larger algebraic value of
$S_b$
leads to a more rapid decay of the streamwise pressure gradient
$P$
. As a result,
$P_{\alpha =0.4} \lt P_{\alpha =0.3} \lt P_{\alpha =0.2}$
in the first half of the domain and the pressure ratio
$p/p_{0}$
follows the same ordering.
Further downstream, the influence of the shock weakens and the transverse flow becomes dominant. For
$\alpha =0.4$
, the higher pressure on both sides of the symmetry plane drives the flow to converge towards
$y/S=0$
, which enhances the compression and causes
$P$
to recover; for
$\alpha =0.2$
, the lower lateral pressure leads to flow divergence away from the symmetry plane, further weakening the compression. This interplay between shock-curvature-induced and transverse-flow-induced compression causes the ordering of both
$P$
and
$p/p_{0}$
to cross-over. The dominant mechanism shifts from the former to the latter along the streamline. In the rear half of the domain,
$P_{\alpha =0.2} \lt P_{\alpha =0.3} \lt P_{\alpha =0.4}$
. The
$\alpha =0.4$
case ultimately attains the highest pressure level.
The reversal of the streamwise pressure gradient observed in figure 14 is driven by lateral flow effects, which are characterised by the lateral pressure gradient
$L$
. To understand how different shock geometries generate distinct lateral flow patterns, we revisit the three curved shocks defined by (4.1) and use MOTC to compute the distribution of
$L$
in the post-shock flowfield. Figure 15 presents the wall-surface contours of
$L$
from a front-view perspective, together with the corresponding streamline patterns.
Comparison of streamline and lateral pressure gradient contours on the wall of external flowfield with different shock axial ratio.
The three shocks share the same initial cross-section – a circle of radius 0.1 at
$v=0$
– but develop different three-dimensional features due to the choices of
$\alpha$
and
$\beta$
. The case
$\alpha =\beta =0.3$
(black lines) is axisymmetric: the wall contour is circular, the streamlines appear as straight radial lines in the front-view and the lateral pressure gradient is everywhere zero. This serves as the baseline for identifying streamline deflection in the non-axisymmetric cases.
For
$\alpha =0.2,\;\beta =0.3$
(blue lines, left panel), the shock expands more slowly in the
$y$
-direction, leading to a pronounced lateral pressure gradient that peaks near the symmetry plane
$y/S=0$
and decays downstream. The corresponding streamlines deflect substantially towards the major axis, with the strongest deflection occurring in regions of elevated
$L$
. In contrast, for
$\alpha =0.4,\;\beta =0.3$
(orange lines, right panel), the shock expands more rapidly in the
$y$
-direction, producing a more moderate lateral pressure gradient whose maximum value shifts towards the major axis and diminishes further downstream. The streamline deflection is weaker and, in some downstream regions, the sign of
$L$
reverses, indicating a change in the direction of the lateral force.
These differences are directly related to the axial ratio
$a/b$
: a larger axial ratio (major axis in the
$y$
-direction) generates a stronger and more localised lateral pressure gradient near the major axis, while a smaller axial ratio produces a weaker and more widely distributed gradient. Consequently, when selecting a basic shock shape for flowfield design, the lateral pressure gradient distribution must be matched to the design objectives. For applications requiring nearly uniform outlet conditions, an elliptical shock with a modest axial ratio is preferable and the region near the major axis should be avoided. Conversely, when lateral flow is intentionally used – for instance, to sweep the boundary layer – a larger axial ratio combined with a cut near the major axis can enhance the desired transverse motion.
The origin of the lateral pressure gradient can be understood from the circumferential variation of the shock angle. For a shock with given coefficients
$\alpha$
and
$\beta$
, the shock angle varies along the circumference according to the local surface inclination, which is governed by the ratio of the expansion rates in the
$y$
- and
$z$
-directions. When the shock angle is circumferentially uniform, the post-shock pressure is constant and the lateral pressure gradient adjacent to the shock vanishes. A non-uniform shock angle, in contrast, produces circumferentially varying post-shock pressures and thus drives lateral flow between adjacent osculating planes. The strength of this lateral pressure gradient depends not only on the pressure difference, but also on the circumferential distance over which it occurs: a given pressure difference concentrated over a shorter arc length produces a stronger gradient. For an elliptical shock, the shock angle changes most rapidly near the major axis, which is why the largest lateral pressure gradient is consistently found in that region.
As illustrated in figure 15, the initial shock cross-section is circular at
$v=0$
. On the right side (
$\alpha =0.4,\;\beta =0.3$
), as the shock evolves downstream, the rate of change of the shock angle is initially greatest midway between the symmetry planes. As the cross-section gradually becomes elliptical, the lateral pressure gradient induced by this shock-angle variation increasingly dominates, and its peak shifts towards the major axis. Once the shock shape becomes fully elliptical, the maximum of the wall lateral pressure gradient remains near the major axis throughout the remainder of the flowfield. This understanding of the relationship between shock geometry and lateral pressure gradient distribution provides a basis for tailoring basic flowfields to meet specific design requirements.
It is worth noting that although the shock profiles at
$y/S=0$
are identical for the three cases in figure 15, the corresponding wall-surface streamlines differ significantly in length. When the
$y/S=0$
plane coincides with the major axis of the elliptical shock (
$\alpha =0.4$
), the wall streamline extends further downstream than in the axisymmetric case; when it coincides with the minor axis (
$\alpha =0.2$
), the streamline is shortened. This difference arises because the local curvature of the shock in the symmetry plane, and consequently the streamline curvature
$D$
within that plane, varies with the axial ratio, thereby altering the size of the shock-determination domain.
The increase in vorticity following a shock wave is one of the salient characteristics accompanying shock waves and this increase in vorticity can also have an impact on aerospace applications. Figure 16 presents the distribution of streamwise vorticity on the wall of flowfields as described by (2.27). On the left side of the figure,
$\alpha =0.35$
and
$\beta =0.3$
, while on the right side,
$\alpha =0.4$
and
$\beta =0.3$
. For shocks with different axial ratios, the streamwise vorticity increases with the increase of shock axial ratios. For the same shock, in addition to the symmetry plane, the streamwise vorticity progressively intensifies from upstream to downstream, which may have an effect on the viscous flow in the boundary layer.
Comparison of streamwise vorticity contours on the wall of flowfields with different shock axial ratios.
4.2. Three-dimensional features of supersonic internal flow
The internal flowfields considered in this section are generated by a family of curved shocks defined by the parametric equation
\begin{align} r(u, v) = \begin{cases} x = v, \\ y = (-0.1v^{2} - \alpha v + 1)\cos (u), \\ z = (-0.1v^{2} - \beta v + 1)\sin (u), \end{cases} \end{align}
where the coefficients
$\alpha$
and
$\beta$
govern the contraction rate of the shock surface in the
$y$
- and
$z$
-directions, respectively. At the initial station (
$v=0$
), the shock cross-section is a circle of radius 1; downstream, the difference between
$\alpha$
and
$\beta$
causes the cross-section to evolve from circular into elliptical contours.
As in the internal flow case, the transverse shock curvature
$S_b$
in the symmetry plane
$y/S=0$
is given by
$S_b = (b/a^{2})(-\cos \theta )$
, where
$\theta$
is the local shock angle and
$a(v)=-0.1v^{2}-\alpha v+1$
,
$b(v)=-0.1v^{2}-\beta v+1$
are the semi-axes of the contracting elliptic cross-section. With
$\beta$
fixed at 0.3, the shock angle
$\theta$
in this symmetry plane is identical for all three configurations, so the variation of
$S_b$
is governed by
$\alpha$
through the geometric factor
$b/a^{2}$
.
To examine the influence of the axial ratio on the internal flow structure, three configurations are selected with a common
$\beta =0.3$
and varying
$\alpha$
, as shown in figure 17.
Internal shock configurations with different
$\alpha$
values (
$\beta =0.3$
), front view along the
$x$
-direction.
Setting
$\alpha =\beta =0.3$
yields an axisymmetric shock (black line). Increasing
$\alpha$
to 0.4 produces an elliptic shock whose major axis lies in the
$y$
-direction (orange line); conversely, reducing
$\alpha$
to 0.2 results in an elliptic shock with its minor axis along the
$y$
-direction (blue line). Because
$\beta$
is held constant, the shock profiles in the symmetry plane
$y/S=0$
are identical for all three cases; the differences in the post-shock flowfield therefore arise from the variation in
$\alpha$
, which alters the axial ratio, the spanwise contraction and, consequently, the transverse shock curvature
$S_b$
. This configuration allows a direct assessment of how the axial ratio influences the three-dimensional flow structure, including the pressure distribution and the reflected shock pattern discussed subsequently.
Figure 18 compares the pressure and the pressure gradient along the streamline for the three shocks at
$y/S = 0$
, where the solid line represents the streamwise pressure gradient
$(P)$
and the dashed line represents the pressure ratio
$p/p_{0}$
. Because all three shocks share the same
$\beta$
, both the local shock angle
$\theta$
and the shock curvature in the
$y/S=0$
plane are identical; consequently, the initial pressure and the initial streamwise pressure gradient immediately behind the shock at
$x/S = 0$
are the same for all three cases. The differences that develop downstream are introduced solely through
$\alpha$
, which determines the spanwise contraction rate and thereby governs the streamwise evolution of the transverse shock curvature
$S_{b}$
.
Pressure and pressure gradient along the streamline at
$y/S=0$
of internal flow.
In the shock-dependent region, the flow is dominated by shock-induced compression. For the internal flow,
$\theta \in [\pi /2,\pi ]$
and
$S_b = (b/a^{2})(-\cos \theta )$
is positive. Since
$\beta$
is fixed,
$\theta$
is the same for all three cases, so the variation of
$S_b$
is determined by
$b/a^{2}$
. A larger
$\alpha$
accelerates the contraction in the
$y$
-direction, producing a smaller semi-axis
$a$
and therefore a larger
$b/a^{2}$
; with
$-\cos \theta \gt 0$
, this means a larger positive
$S_b$
. Hence, the algebraic value of
$S_b$
is largest for
$\alpha =0.4$
and smallest for
$\alpha =0.2$
. A larger
$S_b$
produces a more rapid decay of
$P$
, so that
$P_{\alpha =0.2} \gt P_{\alpha =0.3} \gt P_{\alpha =0.4}$
throughout this region and the pressure ratio follows the same ordering.
Downstream of the shock-dependent region, the influence of the shock weakens and the transverse pressure gradient becomes dominant. For
$\alpha =0.4$
, the smaller shock angle in the
$z$
-direction produces higher pressure on the lateral sides, which drives the flow to converge towards the symmetry plane
$y/S=0$
; this convergence enhances the compression and causes
$P$
to rise immediately upon entering the isentropic compression region. For
$\alpha =0.2$
, in contrast, the lower lateral pressure causes the flow to diverge away from
$y/S=0$
, leading to an abrupt drop in
$P$
. This interplay between shock-curvature-induced and transverse-flow-induced compression causes the ordering of
$P$
to cross-over to
$P_{\alpha =0.4} \gt P_{\alpha =0.3} \gt P_{\alpha =0.2}$
in the isentropic compression region. The dominant mechanism shifts from the former to the latter downstream of the shock-dependent domain. The pressure ratio
$p/p_{0}$
exhibits a corresponding cross-over: the
$\alpha =0.4$
case, which initially had the lowest pressure, experiences the most rapid pressure rise downstream and ultimately attains the highest pressure level, while the
$\alpha =0.2$
case shows a decelerated pressure increase.
The reversal of the streamwise pressure gradient observed in figure 18 is driven by lateral flow effects, which are characterised by the lateral pressure gradient
$L$
. Figure 19 presents the wall-surface contours of
$L$
from a front-view perspective, together with the corresponding streamline patterns, for the three internal shock configurations.
Comparison of streamline and lateral pressure gradient contours on the wall of internal flowfield with different shock axial ratios.
The case
$\alpha =\beta =0.3$
(black lines) is axisymmetric: the wall contour is circular, the streamlines appear as straight radial lines in the front-view and the lateral pressure gradient is everywhere zero. This serves as the baseline for identifying streamline deflection in the non-axisymmetric cases.
For
$\alpha =0.2,\;\beta =0.3$
(blue lines, left panel), the shock contracts more rapidly in the
$y$
-direction, resulting in an elliptic cross-section whose major axis lies in the
$z$
-direction. The larger axial ratio produces a stronger lateral pressure gradient, with the maximum value located in the vicinity of the major axis. The corresponding streamlines deflect substantially towards the major axis, with the most pronounced deflection occurring in regions of elevated
$L$
. In the isentropic compression region downstream of the reflected shock, the lateral pressure gradient becomes discontinuous due to the first-order continuity of the wall surface, but remains elevated, sustaining the lateral deflection.
In contrast, for
$\alpha =0.4,\;\beta =0.3$
(orange lines, right panel), the shock contracts more slowly in the
$y$
-direction, yielding an elliptic cross-section whose major axis lies in the
$y$
-direction. The smaller axial ratio produces a more moderate lateral pressure gradient that is distributed more uniformly along the circumference and the streamline deflection is correspondingly weaker.
Taken together, the results in figures 18 and 19 reveal a consistent two-regime behaviour for internal flows. In the upstream region, the flow is dominated by the shock-induced compression, where the axial ratio governs the transverse shock curvature
$S_{b}$
and thus the decay rate of the streamwise pressure gradient. Further downstream, or beyond the shock-dependent domain, the transverse flow becomes dominant: the lateral pressure gradient drives flow convergence or divergence, which in turn determines the recovery or further decline of the streamwise compression. This transition from shock-dominated to transverse-flow-dominated behaviour is a characteristic feature of three-dimensional internal flowfields and provides a systematic basis for tailoring the compression process through the choice of shock geometry.
After that, a reflection analysis was conducted on the shock wave given in figure 10. Elliptical shocks exhibit different reflection structures at different locations. Since the reflection structure directly affects the flowfield performance, it is essential to study the reflection structure during the design of non-axisymmetric flowfields. Typically, in basic flowfield design, a centre-body is used to avoid Mach reflection. This paper plots shock polar in the reflection point area to analyse the reflection structure of non-axisymmetric flowfields and the design issues of the centre-body. Figure 20 shows a series of shock polar and polar streamline diagrams for the symmetrical plane of the above-mentioned elliptical flowfield at different centre-body heights, solved using 3D-CST-boos and MOTC. The horizontal axis represents the flow deflection angle and the vertical axis represents the pressure ratio. The red solid line indicates the incident shock, which has the same shock polar when the incoming flow conditions are given. The red dashed lines represent the reflected shock at different centre-body heights. The blue solid line represents the polar streamlines. The flow structure around the reflection point can be represented by points 1, 2 and 3, where point 1 represents the incoming flow state and is deflected to point 2 by the shock. Due to the limitation of the outlet, the flow after the reflected shock wave needs to be deflected to point 3. Figure 20(a) shows the shock polar in the
$z/S = 0$
section, where the left dashed line is the reflected shock in the flowfield of figure 10, and the right dashed line is the reflected shock wave at the minimum centre-body height determined by the detachment criterion. When performing internal flow calculations, characteristics intersect near the axis, making the subsequent flowfield unsolvable. The position of the reflected shock at this singular line is indicated by the middle dashed line as determined by MOTC. Figure 20(b) shows the shock polar in the
$y/S = 0$
section, where the left dashed line is the reflected shock in the flowfield of figure 10, the right dashed line is the reflected shock wave at the cone point and the middle dashed line is the reflected shock wave at the singular line position. Due to the presence of the singular line, the reflected shock wave on the right theoretically does not exist. Therefore, for a given internal flow shock wave, the minimum centre-body height should be determined by this singular line. For elliptical flowfields, the flow at the major and minor axes has different centres of curvature, which also leads to differences in reflection patterns. Comparing figures 20(a) and 20(b), it can be found that the minimum centre-body height in the
$z/S = 0$
section is
$h = 0.48$
, while in the
$y/S = 0$
section, it is
$h = 0.21$
, meaning that a higher centre-body is needed at the major axis to ensure that the shock shape matches the pre-designed shock. This is due to the smaller radius of curvature at the major axis. This also makes it harder for the reflected shock at the major axis to leave the double solution region, easily leading to Mach reflection. Therefore, for the inverse design of internal flow shock waves, selecting a basic flowfield near the minor axis can rely on a shorter centre-body to avoid Mach reflection.
Shock polar diagrams and polar streamlines for the internal elliptic shock flowfield: (a)
$z/S = 0$
plane; (b)
$y/S = 0$
plane.
In figure 20, the blue line represents the polar streamline, whose expression is
$d(\ln (p/p_0))/\mathrm{d}\delta = (\gamma M^2 P)/D$
, which can be used to solve for the aerodynamic parameters behind the Mach stem. It can be observed from the figure that due to three-dimensional flow effects, as the centre-body height decreases, the slope of the ‘polar streamline’ behind the incident shock first changes from negative to positive and then back to negative, with the singular line appearing between these two transitions, that is, between the Crocco point and the Thomas point. This may imply that in this form of internal flow shock wave, it is not possible to generate a Thomas point while preserving the prescribed shock shape.
5. Conclusions
A method of three-dimensional curved characteristics (MOTC) has been developed to uncover the downstream flow physics that originate at a three-dimensional curved shock, but cannot be inferred from the shock surface alone. The approach combines two components: the three-dimensional curved shock theory (3D-CST-boos), which supplies the post-shock pressure gradients and streamline curvatures at the shock surface, and the MOTC marching procedure, which propagates these initial quantities along streamlines and characteristics to construct the complete downstream flowfield. Validated against inviscid CFD with errors below 0.5 % for both external and internal supersonic flows, MOTC reveals the following shock-initiated flow phenomena that are inaccessible to conventional planar osculating methods:
-
(i) Cross-over of the streamwise pressure gradient. In the near-shock region, a larger transverse shock curvature
$S_b$
produces a more rapid decay of the streamwise pressure gradient
$P$
, so the ordering of
$P$
is governed by the shock geometry. Downstream, the lateral pressure gradient
$L$
becomes dominant, driving flow convergence or divergence that reverses the ordering of
$P$
. The shock with initially weaker compression may thus attain the highest terminal pressure. -
(ii) Axial ratio control of lateral compression. The lateral pressure gradient
$L$
is governed by the axial ratio of the elliptic shock cross-section. A larger axial ratio concentrates a stronger
$L$
near the major axis, producing pronounced streamline deflection and lateral convergence, which provides a guideline for shock shape selection in forebody or inlet design. -
(iii) Shock-geometry-dependent reflection asymmetry. The reflected shock structure inherits the asymmetry of the incident shock. The major axis plane, with its smaller radius of curvature, requires a higher centre-body to avoid Mach reflection, while the minor axis admits a shorter centre-body. The minimum admissible centre-body height is constrained by a singular line originating from the three-dimensional shock curvature.
These findings demonstrate that the combination of 3D-CST-boos and MOTC establishes a complete analytical chain from the shock surface to the far field: the shock theory defines the initial gradients that seed the downstream physics, while MOTC advances these gradients to reveal flow phenomena that are not evident from the shock alone. This shock-to-flowfield framework provides a systematic basis for both the physical understanding and the inverse design of three-dimensional supersonic flows where transverse shock curvature and lateral pressure gradients play a critical role. The gradients that describe three-dimensional flow trends make MOTC a good candidate for three-dimensional flow analysis and inverse design of three-dimensional supersonic flowfields.
Funding
The authors would like to acknowledge the support of the National Natural Science Foundation of China (grant numbers 12372295, 12302389 and 12472337) and Natural Science Foundation of Fujian Province (grant number 2025J010011).
Declaration of interests
The authors report no conflict of interest.
Appendix A.
The coefficients in (2.4) can be expressed as
\begin{align} \left . \begin{aligned} J_{1} & =\left ( -\frac {B_{2}}{A_{2}} + \frac {-B_{2}A_{2}C_{2} - B_{2}^{2}C_{2} + B_{1}B_{2}C_{2}}{A_{2}^{2}D_{2} + A_{2}^{2}C_{2} - A_{2}B_{2}C_{2}} + \frac {B_{1}}{A_{2}} + \frac {B_{2}C_{1}}{A_{2}D_{2} + A_{2}C_{2} - B_{2}C_{2}} \right )\!, \\[4pt]J_{2} & =\left ( -\frac {B_{2}}{A_{2}} + \frac {-B_{2}A_{2}C_{2} - B_{2}^{2}C_{2}}{A_{2}^{2}D_{2} + A_{2}^{2}C_{2} - A_{2}B_{2}C_{2}} + \frac {2B_{1}}{A_{2}} + \frac {2B_{1}B_{2}C_{2} - B_{2}A_{2}D_{1}}{A_{2}^{2}D_{2} + A_{2}^{2}C_{2} - A_{2}B_{2}C_{2}} \right )\!, \\[4pt]J_{3} & =\left ( \frac {A_{1}}{A_{2}} + \frac {A_{1}B_{2}C_{2} - A_{2}B_{2}C_{1}}{A_{2}^{2}D_{2} + A_{2}^{2}C_{2} - A_{2}B_{2}C_{2}} \right )\!, \\[4pt]K_{1} & =\left ( \frac {-A_{2}C_{2} - B_{2}C_{2}}{A_{2}D_{2} + A_{2}C_{2} - B_{2}C_{2}} + \frac {B_{1}C_{2} + A_{2}C_{1}}{A_{2}D_{2} + A_{2}C_{2} - B_{2}C_{2}} \right )\!, \\[4pt]K_{2} & =\left ( \frac {B_{1}C_{2} + A_{2}C_{1}}{A_{2}D_{2} + A_{2}C_{2} - B_{2}C_{2}} \right )\!, \\[4pt]K_{3} & =\left ( \frac {-A_{2}C_{2} - B_{2}C_{2} + 2B_{1}C_{2} - A_{2}D_{1}}{A_{2}D_{2} + A_{2}C_{2} - B_{2}C_{2}} \right )\!, \\[4pt]A_{1} & = p_{1}, \\[4pt]B_{1} & = V_{1n}^{2}\rho _{1}, \\[4pt]A_{2} & = p_{2}, \\[4pt]B_{2} & = V_{2n}^{2}\rho _{2}, \\[4pt]C_{1} & = \frac {\gamma }{\gamma - 1}\frac {p_{1}}{\rho _{1}}, \\[4pt]D_{1} & = V_{1n}^{2}, \\[4pt]C_{2} & = \frac {\gamma }{\gamma - 1}\frac {p_{2}}{\rho _{2}}, \\[4pt]D_{2} & = V_{2n}^{2}, \\[4pt]L_{1} & = \frac {V_{2m}}{V_{2n}} / \left ( \frac {\gamma p_2}{\rho _2 V_{n2}^{2}} - 1 \right )\!, \\[4pt]L_{2} & = \frac {\gamma V_{2m}}{V_{2n}} / \left ( \frac {\gamma p_2}{\rho _2 V_{2n}^{2}} - 1 \right )\!, \\[4pt]L_{3} & = \gamma / \left ( \frac {\gamma p_2}{\rho _2 V_{2n}^{2}} - 1 \right ) \!, \\[4pt]L_{4} & = \frac {\gamma V_{2m}^{2}}{V_{2n}^{2}} / \left ( \frac {\gamma p_2}{\rho _2 V_{2n}^{2}} - 1 \right )\!. \end{aligned} \right \} \end{align}
Appendix B.
The
$s$
-
$n$
-
$l$
forms an orthogonal system consisting of the flow direction and the other two perpendicular directions, and their unit vectors satisfy the following relationships:
\begin{align} \left |\boldsymbol{e}_i\right |=1,\ \ i=s,n,l, \\[-28pt] \nonumber \end{align}
\begin{align} \boldsymbol{e}_i\boldsymbol{\cdot }\boldsymbol{e}_j= \begin{cases} 0,\quad i\neq j \\[4pt] 1,\quad i=j \end{cases} ,\quad i,j=s,n,l. \\[0pt] \nonumber \end{align}
In curvilinear coordinates, the partial derivatives of the unit vectors are expressed as follows:
\begin{align} \frac {\partial e_{i}}{\partial i}=-\sum _{j=1, j\neq i}^{3}\frac {1}{H_{j}}\frac {\partial H_{i}}{\partial j}e_{j}, \quad i=s, n, l, \\[-28pt] \nonumber \end{align}
\begin{align} \frac {\partial e_{i}}{\partial j}=\frac {1}{H_{i}}\frac {\partial H_{j}}{\partial i}\, e_{j}, \quad i, j=s, n, l, \quad i\neq j. \\[0pt] \nonumber \end{align}
The coefficients are expressed as the partial derivatives of the deflection angles
$\delta$
and
$\eta$
, as follows:
\begin{align} \frac {\partial \delta }{\partial n}=\frac {\partial \boldsymbol{e}_{s}}{\partial n}\boldsymbol{e}_{n}=-\frac {\partial \boldsymbol{e}_{n}}{\partial n}\boldsymbol{e}_{s}=\frac {1}{H_{s}}\frac {\partial H_{n}}{\partial s} , \\[-28pt] \nonumber \end{align}
\begin{align} \frac {\partial \eta }{\partial l}=\frac {\partial \boldsymbol{e}_{s}}{\partial l}\boldsymbol{e}_{l}=-\frac {\partial \boldsymbol{e}_{l}}{\partial l}\boldsymbol{e}_{s}=\frac {1}{H_{s}}\frac {\partial H_{l}}{\partial s}, \\[-28pt] \nonumber \end{align}
\begin{align} \frac {\partial \delta }{\partial s}=\frac {\partial \boldsymbol{e}_{s}}{\partial s}\boldsymbol{e}_{n}=-\frac {\partial \boldsymbol{e}_{n}}{\partial s}\boldsymbol{e}_{s}=-\frac {1}{H_{n}}\frac {\partial H_{s}}{\partial n}, \\[-28pt] \nonumber \end{align}
\begin{align} \frac {\partial \eta }{\partial s}=\frac {\partial \boldsymbol{e}_{s}}{\partial s}\boldsymbol{e}_{l}=-\frac {\partial \boldsymbol{e}_{l}}{\partial s}\boldsymbol{e}_{s}=-\frac {1}{H_{l}}\frac {\partial H_{s}}{\partial l}. \\[0pt] \nonumber \end{align}
Mass conservation equation:
\begin{eqnarray} \begin{aligned} \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{V}) & = \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho V \boldsymbol{e}_s) \\[4pt]& = \boldsymbol{e}_s \frac {\partial (\rho V \boldsymbol{e}_s)}{\partial s} + \boldsymbol{e}_n \frac {\partial (\rho V \boldsymbol{e}_s)}{\partial n} + \boldsymbol{e}_l \frac {\partial (\rho V \boldsymbol{e}_s)}{\partial l} \\[4pt]& = \frac {\partial (\rho V)}{\partial s} + \rho V \boldsymbol{e}_s \frac {\partial \boldsymbol{e}_s}{\partial s} + \rho V \boldsymbol{e}_n \frac {\partial \boldsymbol{e}_s}{\partial n} + \rho V \boldsymbol{e}_l \frac {\partial \boldsymbol{e}_s}{\partial l} \\[4pt]& = \frac {\partial (\rho V)}{\partial s} + 0 + \rho V \boldsymbol{e}_n \frac {\partial \boldsymbol{e}_s}{\partial n} + \rho V \boldsymbol{e}_l \frac {\partial \boldsymbol{e}_s}{\partial l} \\[4pt]& = \frac {\partial (\rho V)}{\partial s} + \rho V \frac {\partial \delta }{\partial n} + \rho V \frac {\partial \eta }{\partial l} = 0. \end{aligned} \end{eqnarray}
Equation (2.7) can be obtained.
Momentum equation:
\begin{eqnarray} \begin{aligned} & \rho (\boldsymbol{V} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{\cdot }\boldsymbol{V} + \boldsymbol{\nabla }p \\[4pt]& \quad = \rho V \boldsymbol{e}_s \boldsymbol{\cdot }\left [ \boldsymbol{e}_s \frac {\partial (V \boldsymbol{e}_s)}{\partial s} + \boldsymbol{e}_n \frac {\partial (V \boldsymbol{e}_s)}{\partial n} + \boldsymbol{e}_l \frac {\partial (V \boldsymbol{e}_s)}{\partial l} \right ] \\[4pt]&\quad\quad + \frac {\partial p}{\partial s} \boldsymbol{e}_s + \frac {\partial p}{\partial n} \boldsymbol{e}_n + \frac {\partial p}{\partial l} \boldsymbol{e}_l \\[4pt]& \quad = \rho V \frac {\partial (V \boldsymbol{e}_s)}{\partial s} + \frac {\partial p}{\partial s} \boldsymbol{e}_s + \frac {\partial p}{\partial n} \boldsymbol{e}_n + \frac {\partial p}{\partial l} \boldsymbol{e}_l \\[4pt]& \quad = \rho V \frac {\partial V}{\partial s} \boldsymbol{e}_s + \frac {\partial p}{\partial s} \boldsymbol{e}_s + \rho V^2 \frac {\partial \boldsymbol{e}_s}{\partial s} + \frac {\partial p}{\partial n} \boldsymbol{e}_n + \frac {\partial p}{\partial l} \boldsymbol{e}_l \\[4pt]& \quad = \left (\rho V \frac {\partial V}{\partial s} + \frac {\partial p}{\partial s} \right ) \boldsymbol{e}_s + \left (\rho V^2 \frac {\partial \delta }{\partial s} + \frac {\partial p}{\partial n} \right ) \boldsymbol{e}_n + \left (\rho V^2 \frac {\partial \eta }{\partial s} + \frac {\partial p}{\partial l}\right ) \boldsymbol{e}_l = 0. \end{aligned} \end{eqnarray}
x/S=0.4
x/S=0.7
x
y
x/S=1.0
x
y
α
β=0.3
x
y/S=0
α
β=0.3
x
y/S=0
z/S=0
y/S=0
