3.1 Fourier transform method

When the test beam passes through a free jet with a variable refractive-index field, the background fringes are changed into the deformed fringe patterns because of the phase shift caused by the variations of the light speed as the beam passes through the test field. Typical profiles taken using a camera showing background and the deformed fringe patterns are illustrated in figure 3 as lines of constant phase. The parallel and equally spaced fringes shown as blue solid lines in figure 3(a) are also referred to as wedge fringes. The interval $b$ in figure 3(a) denotes the distance between two successive crests of the background fringes, and it is a function of the intersection angle between the reference and the test beams and the wavelength of the laser light used in the experiment. The red dashed line in figure 3(b) shows the intensity profile $g(y,z_{0})$ corresponding to the deformed fringe pattern at a particular axial position $z_{0}$, and it can be given by (Takeda et al. Reference Takeda, Ina and Kobayashi1982; Yagi et al. Reference Yagi, Inoue, Nakao, Ono and Miyazato2017)

(3.1)$$\begin{eqnarray}g(y,z_{0})=g_{0}(y,z_{0})+g_{1}(y,z_{0})\cos [k_{0}y-\unicode[STIX]{x0394}\unicode[STIX]{x1D711}(y,z_{0})].\end{eqnarray}$$

Here, the phase shift $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}(y,z_{0})$ contains the desired information on the density field in the free jet, and $g_{0}(y,z_{0})$ and $g_{1}(y,z_{0})$ represent unwanted irradiance variations arising from the non-uniform light reflection or transmission when the test beam passes through a free jet, and $k_{0}=2\unicode[STIX]{x03C0}/b$. The coordinates $y$ and $z$ form the vertical plane, which is perpendicular to the test beam propagation direction, and the $x$-axis is taken as the direction in which the test beam propagates after being reflected at mirror 2, as shown in figure 1.

Figure 3. Variation of background and deformed fringes by the finite fringe method.

Equation (3.1) can be rewritten in the following form:

(3.2)$$\begin{eqnarray}g(y,z_{0})=g_{0}(y,z_{0})+c(y,z_{0})\exp (\text{i}k_{0}y)+c^{\ast }(y,z_{0})\exp (-\text{i}k_{0}y),\end{eqnarray}$$

with

(3.3)$$\begin{eqnarray}c(y,z_{0})={\textstyle \frac{1}{2}}g_{1}(y,z_{0})\exp [-\text{i}\unicode[STIX]{x0394}\unicode[STIX]{x1D711}(y,z_{0})],\end{eqnarray}$$

where $\text{i}$ is the imaginary unit and the asterisk $^{\ast }$ denotes the complex conjugate.

The Fourier transform of (3.2) with respect to $y$ is given by

(3.4)$$\begin{eqnarray}G(k,z_{0})=G_{0}(k,z_{0})+C(k-k_{0},z_{0})+C^{\ast }(k+k_{0},z_{0}),\end{eqnarray}$$

where the capital letters denote the Fourier transforms of the respective primitive functions, and $k$ is the spatial wavenumber in the $y$ direction. Since the spatial variations of $g_{0}(y,z_{0})$, $g_{1}(y,z_{0})$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}(y,z_{0})$ are slow compared with the spatial frequency $k_{0}$ when the interval between fringes is sufficiently small, the Fourier spectra in (3.4) are separated by the wavenumber $k_{0}$ and have the three independent peaks as schematically shown in figure 4(a). We make use of either of these two spectra on the carrier, say $C(k-k_{0},z_{0})$, and translate it by $k_{0}$ on the wavenumber axis towards the origin to obtain $C(k,z_{0})$ as shown in figure 4(b). The unwanted background variation $G_{0}(k,z_{0})$ has been filtered out at this stage by the pertinent bandpass filter.

Figure 4. Fourier transform method for fringe-pattern analysis.

Applying the inverse Fourier transform of $C(k,z_{0})$ with respect to $k$ to obtain $c(y,z_{0})$ defined by (3.3) and taking the logarithm of (3.3) leads to

(3.5)$$\begin{eqnarray}\ln c(y,z_{0})=\ln \frac{g_{1}(y,z_{0})}{2}-\text{i}\unicode[STIX]{x0394}\unicode[STIX]{x1D711}(y,z_{0}).\end{eqnarray}$$

Consequently, the phase shift $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}(y,z_{0})$ in the imaginary part of (3.5) can be completely separated from the unwanted amplitude variation $g_{1}(y,z_{0})$ in the real part. It is necessary to make the interval $b$ between parallel fringes as narrow as possible to avoid islands in the fringe pattern. Wide fringes move farther away from the original location than narrower fringes, and, as a result, they cover regions that are considerably different from the optical retardations. This results in the creation of islands (Winckler Reference Winckler1948).

3.2 Abel inversion method

A free jet with an axisymmetric refractive-index field is depicted in figure 5(a). The refractive index $n$ in the flow field depends only on the radial distance $r$ from the centreline O ($x=y=0$) for a fixed axial distance $z_{0}$ measured from the nozzle exit plane, i.e.

(3.6a,b)$$\begin{eqnarray}n=n(r,z_{0}),\quad \text{for }0\leqq r\leqq R,\quad n=n_{a},\quad \text{for }R<r,\end{eqnarray}$$

where $n_{a}$ is the refractive index for the ambient air and $R$ is the radial distance to the jet boundary.

Refractive-index values can be obtained directly from the measured fringe positions on a given cross-section normal to the jet axis, and each measured cross-section is independent of others. As illustrated in figure 5(a), the optical path difference $\unicode[STIX]{x1D6EC}(y,z_{0})$ of the test beam that passes through the field with and without the jet can be expressed by

(3.7)$$\begin{eqnarray}\unicode[STIX]{x1D6EC}(y,z_{0})=2\int _{y}^{R}\frac{[n(r,z_{0})-n_{a}]}{\sqrt{r^{2}-y^{2}}}r\,\text{d}r.\end{eqnarray}$$

The fringe shift $\unicode[STIX]{x0394}y$ at any location on a recording medium can be related to $\unicode[STIX]{x1D6EC}(y,z_{0})$ as follows:

(3.8)$$\begin{eqnarray}\unicode[STIX]{x0394}y=\frac{b}{\unicode[STIX]{x1D706}_{0}}\unicode[STIX]{x1D6EC}(y,z_{0}),\end{eqnarray}$$

where $\unicode[STIX]{x1D706}_{0}$ is the wavelength of light in vacuum for the laser used. The use of the Gladstone–Dale relation

(3.9)$$\begin{eqnarray}n(r,z)=1+K\unicode[STIX]{x1D70C}(r,z)\end{eqnarray}$$

with the combination of (3.7) and (3.8) gives

(3.10)$$\begin{eqnarray}\unicode[STIX]{x0394}y=\frac{2bK}{\unicode[STIX]{x1D706}_{0}}\int _{y}^{R}\frac{[\unicode[STIX]{x1D70C}(r,z_{0})-\unicode[STIX]{x1D70C}_{a}]}{\sqrt{r^{2}-y^{2}}}r\,\text{d}r,\end{eqnarray}$$

where $K$ is the Gladstone–Dale constant and $\unicode[STIX]{x1D70C}_{a}$ is the density of ambient air.

Figure 5. Light beam passing through axisymmetric refractive-index field.

Once the fringe shifts for a fixed streamwise location $z_{0}$ are obtained by the Fourier transform method, the density field is found by inverting (3.10) using the Abel inversion (Sneddon Reference Sneddon1972):

(3.11)$$\begin{eqnarray}\unicode[STIX]{x1D70C}(r,z_{0})=\unicode[STIX]{x1D70C}_{a}-\frac{\unicode[STIX]{x1D706}_{0}}{\unicode[STIX]{x03C0}bK}\int _{r}^{R}\frac{1}{\sqrt{y^{2}-r^{2}}}\frac{\text{d}}{\text{d}y}(\unicode[STIX]{x0394}y)\,\text{d}y.\end{eqnarray}$$

Measurements of the fringe shifts only provide a discrete set of values, and therefore (3.11) can only be solved numerically. Several algorithms for solving the Abel inversion equation have been proposed so far (Nestor & Olsen Reference Nestor and Olsen1960; Bradley Reference Bradley1968; Behjat, Tallents & Neely Reference Behjat, Tallents and Neely1997; Malka et al. Reference Malka, Coulaud, Geindre, Lopez, Najmudin, Neely and Amiranoff2000; Álvarez, Rodero & Quintero Reference Álvarez, Rodero and Quintero2002; Shimomiya, Ono & Miyazato Reference Shimomiya, Ono and Miyazato2013). Among them, the integral in (3.11) is evaluated by using the algorithm proposed by Nestor & Olsen (Reference Nestor and Olsen1960) as follows:

(3.12)$$\begin{eqnarray}\unicode[STIX]{x1D70C}(r_{k},z_{0})=\unicode[STIX]{x1D70C}_{a}-\frac{\unicode[STIX]{x1D706}_{0}}{\unicode[STIX]{x03C0}^{2}K\unicode[STIX]{x0394}w}\mathop{\sum }_{n=k}^{N}B_{k,n}\unicode[STIX]{x0394}\unicode[STIX]{x1D711}(y_{n},z_{0}),\end{eqnarray}$$

where $r_{k}=k\unicode[STIX]{x0394}w$ ($k=0,1,2,\ldots ,N$) is the radial distance from the jet centreline, $\unicode[STIX]{x0394}w$ is the sampling interval, which is the same as the distance between neighbouring pixels on the recording medium, $N$ is the total number of intervals in the medium, $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}(y_{n},z_{0})=2\unicode[STIX]{x03C0}y_{n}/b$, $y_{n}\leqq y<y_{n+1},$$y_{n}=n\unicode[STIX]{x0394}w$ and

(3.13a,b)$$\begin{eqnarray}B_{k,n}=-A_{k,n}\quad \text{if }n=k,\quad B_{k,n}=A_{k,n-1}-A_{k,n}\quad \text{if }n\geqq k+1,\end{eqnarray}$$

with

(3.14)$$\begin{eqnarray}A_{k,n}=\frac{\sqrt{(n+1)^{2}-k^{2}}-\sqrt{n^{2}-k^{2}}}{2n+1}.\end{eqnarray}$$

Density values can be captured at the spatial resolution determined by the sampling interval. Interferogram images of a free jet obtained by the Mach–Zehnder interferometer were recorded in 30 pictures at a sampling frequency of 10 kHz and formed onto a complementary metal oxide semiconductor (CMOS) camera (The Imaging Source, DFK72BUC02), which recorded a JPEG RGB image (8-bit each colour) at a resolution of 2592 $\times$ 1944 square pixels, and then were turned into an HSV image (8-bit each colour) according to the hue ($H$)-saturation ($S$)-value ($V$) colour model. Since the intensity distributions of interferograms recorded are represented by a single parameter, $V$, with an 8-bit greyscale image, the distributions of background and deformed fringes with 256 different possible intensities can be calculated from the Mach–Zehnder interferometer for the density field in the free jet. The plane of focus was located in the nozzle axis. Based on the physical and imaged dimensions of the test nozzle, the spatial resolution of the present imaging system is $4.1~\unicode[STIX]{x03BC}\text{m}\pm 130~\text{nm}$. The average and standard deviation for the experiment were obtained over the 30 samples recorded in the experiment.

The Fourier transform method utilized in the present study for the phase shift analysis can be seen as a sub-fringe measuring technique similar to heterodyne interferometer (Malacara Reference Malacara2008), in which the spatial resolution is smaller than one fringe. The spatial resolution ($4.1~\unicode[STIX]{x03BC}\text{m}\pm 130~\text{nm}~\text{per}~\text{pixel}$) of the present Mach–Zehnder interferometer system is identical to the distance, $\unicode[STIX]{x0394}w$, between two adjacent pixels on the recording medium after being focused by the imaging lens (see figure 1) and also corresponds to the minimum fringe distortion that can be reliably measured. Since the density field is reconstructed using equation (3.12), the density uncertainty is affected by the choices of $\unicode[STIX]{x0394}w$, $b$, $\unicode[STIX]{x1D706}_{0}$ and $K$ as well as the Abel inversion algorithms. Considering that $\unicode[STIX]{x0394}w$ and $b$ are simultaneously focused and that the phase shift is calculated by $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}=2\unicode[STIX]{x03C0}n\unicode[STIX]{x0394}w/b$ ($n$ is the natural number), it is possible to remove $b$ during the reconstruction of the density. Note that $b$ is set by the deliberate misalignment of the mirrors and the beam splitters before starting experiments. Therefore, $b$ can be altered arbitrarily regardless of the focal point of the imaging lens. However, the range of $b$ usually remains narrow since $K$ and $\unicode[STIX]{x1D706}_{0}$ have little effect on the density uncertainty relative to the phase measurement uncertainty (Mercer & Raman Reference Mercer and Raman2002). This leads to the conclusion that the effects of the measurement uncertainty on $\unicode[STIX]{x0394}w$ as well as numerical errors of the Abel inversion algorithms on the density uncertainty could also be insignificant. As a result, the density uncertainty of the Mach–Zehnder interferometer with the finite fringe method would depend significantly on the uncertainty in the phase shift measurement.

5.1 Comparison of experimental results with simulations

A comparison between the experimental measurement of the density contour plot at the cross-section including the jet centreline and the corresponding numerical result is depicted in figure 7. The contour levels with an interval of $0.1~\text{kg}~\text{m}^{-3}$ are shown in the colour bar at the top, and the spatial resolution in the experimental density map is around $4~\unicode[STIX]{x03BC}\text{m}$. Figure 7 illustrates the various flow features of the shock-cell structures quantitatively, such as the shape and the size of the Mach disk in the first shock-cell as well as the expansion and compression regions, the shock-cell intervals, the jet boundaries, the slip streams produced from the triple points of a Mach disk and the shear layers outside of the jet boundaries. Although the unsteady features are ignored here and a relatively coarse mesh is used, there is very favourable overall agreement between the experiment and the simulation on the shock-cell structures. However, the simulation displays more distinct slip streams when compared with the experiment.

Figure 7. Density contours of underexpanded microjets.

The experimental data also show that the edge location of the shear layer remains almost the same ($r=0.5~\text{mm}$) until the end of the domain, but the prediction does not. This discrepancy might be attributed to two issues: first, axial symmetry is assumed, and thus the nature of the fully three-dimensional effects could be overlooked; and second, the unsteady features are not captured by this RANS simulation. These issues should be investigated in the future. A Mach disk with a diameter of around $300~\unicode[STIX]{x03BC}\text{m}$ can be clearly recognized in the contour plots of the experiment and simulation of figure 7 in this study. This is consistent with the measurement done by André, Castelain & Bailly (Reference André, Castelain and Bailly2014), where the presence of a Mach disk with a small diameter in the first shock-cell is observed by the measurement using particle image velocimetry (PIV). The experiments by André et al. (Reference André, Castelain and Bailly2014) were performed under a nozzle pressure ratio of $\text{NPR}=3.67$ using an axisymmetric convergent nozzle with an inner diameter of 38.7 mm at the nozzle exit.

Figure 7(b) shows two distinct oblique shocks in the second shock-cell. Both of them originate from near the jet centreline and are towards the jet free boundaries. It is noteworthy that a second Mach disk can be observed indistinctly at the positions where the oblique shocks are produced. The slip streams produced from the intersections between the second Mach disk and the oblique shocks extend downstream keeping the shape almost parallel to the jet centreline. The existence of the third slip stream can be found in the third shock-cell. The presence of a second Mach disk was displayed in a velocity contour map obtained from PIV measurements by Edgington-Mitchell, Honnery & Soria (Reference Edgington-Mitchell, Honnery and Soria2014), who demonstrated that the shear layer originating from the triple point of the first Mach disk persists across multiple shock-cells downstream. Their experiment was performed at a nozzle pressure ratio of $\text{NPR}=4.2$, which is slightly higher compared with our experimental condition ($\text{NPR}=4.0$). However, the second Mach disk cannot be observed in figure 7(a). A more sophisticated instrumentation to investigate such a local phenomenon needs to be done.

It should be noted that the conventional schlieren and shadowgraph pictures sometimes show the incident shock or barrel shock just upstream of the Mach disk, but such a shock cannot be clearly seen in the contour plot representation at the jet cross-section for both the experiment and simulation in this study. Therefore, what was considered to be the incident shock wave or the barrel shock wave seen in the conventional pictures so far could be considered to be a compression wave, where there is a noticeable difference between density distributions integrated in the line-of-sight direction and those at the jet cross-section shown in figure 7. As described later, there is some evidence to support this conclusion.

Figure 8. Streamwise density profiles of underexpanded microjets.

A comparison between the experiment and the simulation of the streamwise density profiles at radial distances of $r/D_{e}=0$, 0.25 and 0.5 is shown in figure 8(a–c), where the experimental results are represented with precision error bars. These radial locations are shown as the three horizontal dotted lines in the contour plots of figure 7. The density value at the nozzle exit plane, which is estimated based on the assumption of one-dimensional isentropic flow from the nozzle inlet to the exit, is shown as a leftward arrow on a vertical axis as a reference. The two-way vertical arrow indicates the ideal normal shock jump, which is estimated using the flow properties just upstream of a Mach disk with the assumption of an isentropic flow from the nozzle inlet to a Mach disk.

As a global feature, the experimental density profile in figure 8(a) exhibits a representative property appearing in an underexpanded free jet, i.e. the flow expansion and compression quasi-periodically repeated downstream are responsible for the shock-cell structures. In addition, the density rapidly decreases below the ambient level by expansion waves originating from the nozzle lip. The simulated density profile exhibits a similar trend as the experiment over the entire region except that the simulation shows a smoother density variation at the entire streamwise locations compared with the experiment.

Let us focus on the spatial variation of the centreline density profile behind a Mach disk. The experimental density profile in figure 8(a) shows two distinctive bumps (at $z/D_{e}\simeq 1.2$ and 1.7) just downstream of a Mach disk, which are caused by the wavy behaviour with the divergence/convergence of the slip streams immediately after a Mach disk as clearly seen in figure 7(a). The first divergence of the slip streams leads to the rapid density rise after the density jump by a Mach disk because of the aerodynamic effect of the subsonic flow behind a Mach disk. The waveform with the two significant bumps behind a Mach disk also appears in the density profile obtained from the rainbow schlieren deflectometry by Takano et al. (Reference Takano, Kamikihara, Ono, Nakao, Yamamoto and Miyazato2016) as illustrated later in figure 12. Additionally, the similar waveform just downstream of a Mach disk can be observed in the velocity profile along the jet centreline obtained by PIV of André et al. (Reference André, Castelain and Bailly2014). However, such bumps cannot be seen in the simulated density profile. This may be due to the overlooked three-dimensional effects, unsteady flow features, and possible inadequacy of the current simulation model (e.g. coarse mesh). However, as can be seen in figure 8(a), the overall agreement between the measured and the simulated density profiles is favourable from the nozzle exit ($z/D_{e}=0$) to the exit of the simulation domain ($z/D_{e}=5$). Also, both the predicted shock-cell lengths and the density amplitudes agree well with the measurement quantitatively.

Figure 8(b) shows the measured and the simulated density profiles at the radial location of $r/D_{e}=0.25$. Both of the density profiles capture the reflected shock wave, and the overall agreement between the experimental and the simulated results is satisfactory. Figure 8(c) shows the streamwise density profile along the jet lipline ($r/D_{e}=0.5$). The experimental and the simulated density profiles show significant small bumps quasi-periodically appearing, which are marked as small black circles in the experimental density contour plot of figure 7(a). These bumps are responsible for the reflection as expansion waves of shock waves or compression waves at the jet boundaries. Again, good quantitative agreement for both of the density profiles is reached between experiment and simulation.

Comparisons between the experiment and the simulation for the radial density profiles along the four parallel vertical lines ($z/D_{e}=0.6$, 1.0, 1.5 and 1.8) in the density contour plots in figure 7 are shown in figure 9(a–d). The dashed line parallel to the vertical axis indicates the density value estimated based on the assumption that the flow from the nozzle is isentropically expanded to the back pressure, which is referred to as the fully expanded jet density $\unicode[STIX]{x1D70C}_{j}$ and is estimated to be $1.75~\text{kg}~\text{m}^{-3}$ ($\unicode[STIX]{x1D70C}_{j}/\unicode[STIX]{x1D70C}_{b}=1.48$) in this experiment. The solid and open circles shown on the experimental and the numerical density profiles in figures 9(a) and 9(b) are the radial positions of the local maximum and minimum values, respectively. The maximum density values in the simulated density profiles in both figures 9(a) and 9(b) are in quantitatively good agreement with the fully expanded jet density $\unicode[STIX]{x1D70C}_{j}$, and the corresponding radial positions nominally indicate the inviscid jet boundary (Panda & Seasholtz Reference Panda and Seasholtz1999). The density at the jet boundary approaches the ambient density $\unicode[STIX]{x1D70C}_{b}$ ($=1.18~\text{kg}~\text{m}^{-3}$) through the shear layer radially outside of the boundary. It is noted that $\unicode[STIX]{x1D70C}_{j}$ is different from the ambient density $\unicode[STIX]{x1D70C}_{b}$ even though the fully expanded jet pressure $p_{j}$ has to coincide exactly with the back pressure $p_{b}$. Therefore, we conclude that the position of the local maximum in the density profile for the experiment seems to correspond to that of the inviscid jet boundary. As mentioned before, figure 9(b) demonstrates that the radial density profile just upstream of a Mach disk shows a smooth radial variation over from the local minimum value to the local maximum value. It means that no incident shock or barrel shock is present in front of a Mach disk, but instead compression waves appear. The experimental data of Panda & Seasholtz (Reference Panda and Seasholtz1999) also show the same trend.

Figure 9. Radial density profiles of underexpanded microjets.

As shown in figure 9(c), the radial density profile passing through the middle of the reflected shock exhibits a noticeable discrepancy between the experimental and the simulated results. The experimental density profile is not capable of capturing a reflected shock wave, although the simulated density profile shows that the density remains uniform until it crosses the reflected shock and suddenly drops to the density at the inviscid jet boundary just downstream of the reflected shock. Such a gradual density variation across the reflected shock is also observed in the experimental density profile of Panda & Seasholtz (Reference Panda and Seasholtz1999), but the primary reason remains unclear at the present stage. The experimental density profile in figure 9(c) shows a smooth variation across the slip stream and it seems to be responsible for a mixing caused by a significant spreading of the slip stream opposed to a somewhat sharp variation of the simulated density profile. The simulated density sharply drops to the $\unicode[STIX]{x1D70C}_{j}$ passing through the reflected shock and approaches the ambient condition via the shear layer after a region where the density remains constant at $\unicode[STIX]{x1D70C}_{j}$, while the experimental density profile shows a continuous variation between the reflected shock and the shear layer. In the radial density profiles downstream of the reflected shock shown in figure 9(d), relatively good agreement between the experimental and the simulated results is achieved, with the exception of the regions near the jet centreline.

Given the axisymmetric geometry, the density uncertainty may be affected by radial distance from the jet centreline. Behjat et al. (Reference Behjat, Tallents and Neely1997) demonstrated that the error in the densities deduced by the Abel inversion algorithm with the sixth-order polynomial fitting function for an axisymmetric jet increases towards the jet edges. Álvarez et al. (Reference Álvarez, Rodero and Quintero2002) also discussed the error introduced by the Abel inversion algorithm based on an exact solution including the ridge between the centre and outer edge, and further compared the exact solution with those analysed using different Abel inversion algorithms. It is shown that the Nestor–Olsen algorithm produces the highest error at the centre of the radial distribution and a few errors around the ridge, and that the discretization algorithm can reproduce the exact solution except for around the ridge. However, it should be noted that the studies of Behjat et al. (Reference Behjat, Tallents and Neely1997) and Álvarez et al. (Reference Álvarez, Rodero and Quintero2002) did not take into account a discontinuous variation caused by shocks in the distribution. Since it is not feasible to precisely estimate the density uncertainty attributed to only the Abel inversion algorithm for the case of the jets subject to unstable strong shocks, we have derived the precision errors based on the convolution back-projection (CBP) method (Decker Reference Decker, Lading, Wigley and Buchhave1994) and compared these errors with those from the Abel inversion method. Figure 10(a–d) shows the precision errors of jet radial densities at four different axial locations: $z/D_{e}=0.6$, 1.0, 1.5 and 1.8. It is found that the errors are slightly larger for the Abel inversion method than for the CBP method. However, the errors from both methods remain almost constant inside the jet boundary and then gradually decrease in the radial direction. Note that these errors (or fluctuations) are caused by the coupling effects, such as: strong density fluctuations due to an oscillating Mach disk, associated oblique shocks and slip streams (Panda Reference Panda1998; Emami et al. Reference Emami, Bussmann and Tran2010; Edgington-Mitchell et al. Reference Edgington-Mitchell, Honnery and Soria2014; Li et al. Reference Li, Zhou, Chen, Fan and Xie2018). These effects could be more dominant than the measurement uncertainty. However, it is impossible to separate the uncertainty error from the effects causing the density fluctuations. Clarifying the unsteady behaviour of shock-containing jets would be one of the main objectives in our future works.

Figure 10. Comparison of precision errors for radial densities between Abel inversion method and CBP method.

Figure 11. Radial static pressure profiles obtained from RANS simulation for $\text{NPR}=4.0$.

Figure 12. Comparison among centreline density profiles of underexpanded sonic jets obtained by various optical diagnostics for $\text{NPR}=4.0$.

Figure 11 shows the radial static pressure profiles along four parallel vertical lines ($z/D_{e}=0.6$, 1.0, 1.5 and 1.8) (see figure 7b). The horizontal axis is normalized by the back pressure $p_{b}$, which is the same as the fully expanded jet pressure $p_{j}$. As seen from a comparison between the simulated density and pressure profiles at the corresponding streamwise locations in figures 9 and 11, the static pressures coincide with the back pressure $p_{b}$ at the edge of the inviscid jet boundary indicated by the density profiles at $z/D_{e}=0.6$ and 1.0. For the simulated density profiles at $z/D_{e}=1.5$ (figure 9c) and 1.8 (figure 9d), for analogous reasons, the radial locations ($r/D_{e}\simeq 0.4$ and 0.5, respectively), in which the static pressures and densities are equal to $p_{b}$ and $\unicode[STIX]{x1D70C}_{j}$, seem to be the inviscid jet boundaries. The static pressure in the jet shear layer is equal to the back pressure at all axial locations. This contrasts with the density at the inviscid jet boundary, which gradually decreases towards the ambient density.

5.2 Comparison of centreline density profiles obtained by various optical diagnostics

To obtain density data on shock-containing free jets, various non-intrusive experimental techniques have been presented in the past, including rainbow schlieren deflectometry (Satti et al. Reference Satti, Kolhe, Olcmen and Agrawal2007; Kolhe & Agrawal Reference Kolhe and Agrawal2009; Takano et al. Reference Takano, Kamikihara, Ono, Nakao, Yamamoto and Miyazato2016), background oriented schlieren (Venkatakrishnan Reference Venkatakrishnan2005; Van Hinsberg & Rösgen Reference Van Hinsberg and Rösgen2014; Nicolas et al. Reference Nicolas, Donjat, Léon, Le Besnerais, Champagnat and Micheli2017), moiré schlieren (Tabei, Shirai & Takakusagi Reference Tabei, Shirai and Takakusagi1992), Mach–Zehnder interferometer (Nakamura & Iwamoto Reference Nakamura and Iwamoto1996; Shimomiya et al. Reference Shimomiya, Ono and Miyazato2013), filtered Rayleigh scattering (Forkey et al. Reference Forkey, Lempert and Miles1998; Panda & Seasholtz Reference Panda and Seasholtz1999; Gustavsson & Segal Reference Gustavsson and Segal2005), coherent anti-Stokes Raman scattering (Woodmansee et al. Reference Woodmansee, Iyer, Dutton and Lucht2004) and dual-hologram interferometer (Velásquez-Aguilar et al. Reference Velásquez-Aguilar, Toker, Zamudio-Lara and Arias-Estrada2007). However, to the best of our knowledge, comparative studies between density data obtained with the same nozzle shape and under the same experimental conditions are virtually non-existent. Figure 12 shows a comparison between the density profiles along the jet centreline obtained using the various quantitative flow visualization methods for underexpanded sonic jets emitted from axisymmetric convergent nozzles under the same nozzle operating condition ($\text{NPR}=4$) with different exit diameters. The density profiles normalized by the ambient density $\unicode[STIX]{x1D70C}_{b}$ from Tabei et al. (Reference Tabei, Shirai and Takakusagi1992), Nakamura & Iwamoto (Reference Nakamura and Iwamoto1996), Takano et al. (Reference Takano, Kamikihara, Ono, Nakao, Yamamoto and Miyazato2016) and Nicolas et al. (Reference Nicolas, Donjat, Léon, Le Besnerais, Champagnat and Micheli2017) are superimposed on our measurement (indicated by the blue line) using the Mach–Zehnder interferometer and the numerical result (indicated by the red line) from the RANS simulation with the SST $k$–$\unicode[STIX]{x1D714}$ turbulence model. In order to make it easier to view, the precision error bars are removed from our experimental result and the theoretical density rise by a normal shock wave is also included as a reference.

Before comparing the present experimental and simulated results with previous experimental data, it is important to note that each methodology relies on certain assumptions and restrictions. The data acquired from the Abel inversion method are based on the assumption of an axisymmetric jet for the Mach–Zehnder interferometers and for the rainbow schlieren deflectometry. For the moiré schlieren method, the data measured by computed tomography with multi-viewing schlieren pictures are obtained by rotating a test nozzle around its longitudinal axis in equal intervals. Finally, the instantaneous density field reconstructed three-dimensionally with the simultaneous schlieren pictures is taken by mounting 12 synchronized cameras around a test nozzle for the background oriented schlieren.

What needs to be emphasized here is that, for the density profiles after the first and second shock-cells, the density value obtained from the Mach–Zehnder interferometer developed by Nakamura & Iwamoto (Reference Nakamura and Iwamoto1996) and that from the moiré schlieren developed by Tabei et al. (Reference Tabei, Shirai and Takakusagi1992) have an additional density increment compared to those from the other measuring techniques. The primary reason could be attributed to the oscillation of the shock position (Di Rosa, Chang & Hanson Reference Di Rosa, Chang and Hanson1993; Panda Reference Panda1998; Edgington-Mitchell et al. Reference Edgington-Mitchell, Honnery and Soria2014). When a shock wave oscillates across the time-mean position, the flow Mach number relative to the upstream-travelling shock is higher than that just ahead of the stationary shock. This results in an excessive density increment in addition to the density jump produced by the Rankine–Hugoniot relations. Other possible reasons could be the accuracy in the measuring techniques or the analytical method for reconstructing the density field.

As seen in figure 12, there is good quantitative agreement between the density profile obtained from the present Mach–Zehnder interferometer and those from both the rainbow schlieren deflectometry and the background oriented schlieren. The streamwise locations of the first and second local minima showing the ends of the expansion regions in the first and second shock-cells, respectively, and that of the second local maximum showing the end of the compression region in the second shock-cell, and the density values at their corresponding local minima and maximum, agree extremely well. In addition, the density behind a Mach disk in the first shock-cell for each density profile is almost the same peak value. The density profiles obtained from both the rainbow schlieren and the background schlieren show good quantitative agreement with the simulation over the entire region. It should be noted that the nozzle used by Takano et al. (Reference Takano, Kamikihara, Ono, Nakao, Yamamoto and Miyazato2016) is geometrically similar to that shown in figure 2, while the wall contour of the nozzle utilized by Nicolas et al. (Reference Nicolas, Donjat, Léon, Le Besnerais, Champagnat and Micheli2017) is converged in the downstream direction at a constant half-angle of $30^{\circ }$ and followed by a short straight duct to provide a uniform sonic condition on the jet exit plane. From a comparison between the centreline density profile obtained by the present Mach–Zehnder interferometer and that by the rainbow schlieren deflectometry, it is important to note that there is almost no difference in effects of the nozzle size on the jet density field. Furthermore, the local minima in the density profiles for all of the experimental data as well as the simulation results gradually increase with increasing streamwise distance. This trend is also in good agreement with the results of Panda & Seasholtz (Reference Panda and Seasholtz1999) performed under almost the same experimental condition using filtered Rayleigh scattering.

Figure 13. Three-dimensional representation of shock surface of shock-containing microjet.

5.3 Flow topology of near-field shock systems

The shock dynamics of a shock-containing microjet can be quantitatively evaluated using the magnitude of the density gradient vector calculated by (Maeno et al. Reference Maeno, Kaneta, Morioka and Honma2005; Li et al. Reference Li, Zhou, Chen, Fan and Xie2018)

(5.1)$$\begin{eqnarray}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}|=\left[\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x}\right)^{2}+\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}y}\right)^{2}+\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}z}\right)^{2}\right]^{1/2}.\end{eqnarray}$$

This quantity seems to be capable of capturing inherent phenomena for a strongly underexpanded free jet including the interaction of shock waves with expansion waves. The density increases rapidly across a shock wave or the root of the Prandtl–Meyer expansion fan diverging from a sharp corner. Figure 13 shows a cutaway view of an isosurface of the three-dimensional density gradient field of the shock-containing microjet experimentally measured in this study. From figure 13, the spatial feature of the expansion fan emanating from the nozzle lip, the jet boundary, the short intercepting shock and the reflected shock can clearly be seen. Here, the Mach disk is around $300~\unicode[STIX]{x03BC}\text{m}$ in diameter. It is found that the slip stream just behind the Mach disk is successfully captured in this three-dimensional view.

Figure 14 displays a bird’s-eye view of the magnitude of the density gradient vector on the cross-section including the centreline of the microjet. The flow topology as well as the shock strength are clearly shown. The Mach disk consists of shock systems where the shock strength reaches a maximum value at the centre and then decreases drastically towards the triple line (i.e. a triple point in a two-dimensional representation), while the reflected shock has some local maxima in the shock strength distribution, forming the wavy pattern. The slip stream expands initially just in the downstream region of the Mach disk followed by a rapid contraction. These results were also observed experimentally by Edgington-Mitchell et al. (Reference Edgington-Mitchell, Honnery and Soria2014) and numerically by Li et al. (Reference Li, Zhou, Chen, Fan and Xie2018). Since the flow is subsonic just behind the Mach disk, it decelerates initially and accelerates soon. This contributes to the behaviour of the divergence, and then it contracts in the streamwise direction of the area surrounded by the slip stream. We may conclude that our measurement technique with a high spatial resolution is capable of capturing the micrometre-scale shock dynamics in three-dimensional matter.

Figure 14. Density gradient field of shock-containing microjet.