2.1. Experiment and measurements

The experimental flow configuration is shown in figure 1(a). The Mach $M_\infty =6.14$ flow developed as part of the dynamics of a Ludwieg tube, whereby heated ($T_0=450~\text {K}$) and pressurized gas was contained in the charge tube upstream of the converging–diverging nozzle. Once the fast-acting valve was opened, the gas accelerated through the nozzle to calculated free-stream conditions of $U_\infty =904\ \text {m}~\text {s}^{-1}$, $T_\infty =54~\text {K}$ and $\rho _\infty =0.0274\ \text {kg}~\text {m}^{-3}$, yielding a unit free-stream Reynolds number $Re_\infty /L=7.11\times 10^6~\text {m}^{-1}$. The total test time was approximately 0.2 s, which was separated into two approximately time-stationary periods, lasting 100 ms. Data from the second period are analysed for this work (for detailed characteristics of the AFRL Ludwieg Tube, see Kimmel et al. Reference Kimmel, Borg, Jewell, Lam, Bowersox, Srinivasan, Fuchs and Mooney2017; Kennedy et al. Reference Kennedy, Jewell, Paredes and Laurence2022).

Figure 1. (a) Schematic of the flow configuration. (b) Extended time series of the PCB pressure data and (c) a detail for a 80 $\mathrm {\mu }\text {s}$ window. Signals are offset by 2 kPa for clarity.

The test article was a $\theta =7^\circ$ half-angle circular cone with a sharp nose (tip radius $r_n=0.508$ mm) and a total length of 414 mm. The model was installed at zero incidence to the streamwise direction; it was at room temperature ($T_w=300$ K) prior to the start of the experiments, and changes to the surface temperature during the brief test time are small enough that they can be ignored. The cone was instrumented with six PCB model 132A piezoelectric pressure sensors, positioned at four streamwise positions $\{s_1,s_2,s_3,s_4\}=\{215, 241, 266, 291, 316\}~\text {mm}$, measured downstream of the nose. In addition, the $s_4$ row has two probes offset from the primary ray of sensors by $\pm 8.5^\circ$ in order to assess azimuthal dependence. Simulations and test data suggest that the effective sensing area has a diameter of approximately $d_s=0.97$ mm (Ort & Dosch Reference Ort and Dosch2019). For the case considered, time-resolved pressure data are available for six probes, acquired at a frequency of 5 MHz. A sample of the measurements is shown in figure 1(b,c). In figure 1(c), the pressure signatures show the presence of second-mode wave packets amplifying and advecting downstream. For any instance in figure 1(b,c), the signals do not appear to be chaotic with a broad range of time scales, suggesting the flow has not transitioned to turbulence at the sensor locations. High-speed schlieren measurements were also acquired, which provide additional points of validation above the wall. However, our focus is on wall data, which are considered to be the primary modality taken during flight.

An analysis of the experimental data can aid in the choice of inflow frequencies and wavenumbers adopted in the data-assimilation simulations. The time-resolved pressure measurements are Fourier-analysed using the Welch method, similar to previous efforts (Casper et al. Reference Casper, Beresh, Henfling, Spillers, Pruett and Schneider2016; Kennedy et al. Reference Kennedy, Jewell, Paredes and Laurence2022). The window size of 80 $\mathrm {\mu }\text {s}$ (figure 1c) furnishes $10^3$ non-overlapping Hann windows, which yields converged spectral amplitudes. The spectral resolution was halved by combining every two frequency bins, while conserving energy, in order to reduce the dimensionality of the data-assimilation problem. Although the flow is instantaneously 3-D (figure 1c), it is statistically axisymmetric under nominal conditions, i.e. zero incidence. Since the experimental set-up was nominally axisymmetric, we assume homogeneity in the azimuthal direction and average the spectra of the three probes at $s_4$. The post-processed measurements are then concatenated into two vectors,

(2.1a,b)\begin{align} \boldsymbol{m}_S=[\hat{p}\hat{p}^\star(s_1,f_1), \ldots, \hat{p}\hat{p}^\star(s_4,f_{N_f})]^{\rm T}, \quad \boldsymbol{m}_I=\left[{\sum}_i\hat{p}\hat{p}^\star(s_1,f_i), \ldots,{\sum}_i\hat{p}\hat{p}^\star(s_4,f_i)\right]^{\rm T}, \end{align}

representing the individual spectral amplitudes and overall intensity of the signal.

The experimental spectra and signal intensity, used in the data assimilation, are shown with symbols in figure 2. At position $s_1$, the spectrum is dominated by a peak between $f=225$ and $275$ kHz, which is expected due to the presence of unstable second modes on straight cones (Kennedy et al. Reference Kennedy, Jewell, Paredes and Laurence2022). As the boundary layer thickens downstream, the frequency for the peak energy decreases. By $s_4$, the peak amplitude is near 200 kHz, and higher harmonics, near 400 kHz, also appear, possibly due to nonlinearity from these high-amplitude waves. The total intensity of the signal amplifies between $s_1$ and $s_3$ and appears to saturate between $s_3$ and $s_4$, also signalling the presence of nonlinear effects. The goal of this work, then, is to identify the unknown 2-D and 3-D instability waves upstream of $s_1$ that reproduce these measurements, and to study the associated flow field beyond the scope of the sensor data.

Figure 2. Schematic of the ensemble-variational (EnVar) data assimilation framework.

2.2. Simulation configuration

The simulation domain is shown in figure 1. The flow behind the cone-generated shock is considered in order to interpret the measured wall-pressure spectra in terms of post-shock boundary-layer disturbances. The base flow is axisymmetric and obeys the Taylor–Maccoll approximation above the boundary layer, which is a solution to the Blasius equations in coordinates parallel ($\xi$) and normal ($\eta$) to the cone surface using the Mangler–Levy–Lees transformation. Prescribing the initial base state, $\boldsymbol {q}_B(\xi,\eta )$, in this manner is a well-established approach (Sivasubramanian & Fasel Reference Sivasubramanian and Fasel2015). The inflow Reynolds number, based on the post-shock boundary-layer-edge conditions ($U_e=880\ \text {m}~\text {s}^{-1}$, $T_e=65~\text {K}$ and $\rho _\infty =0.0458\ \text {kg}~\text {m}^{-3}$) is $Re_o\equiv \rho _e U_e L_o/\mu _e =1345$ using the length scale $L_o = \sqrt {\nu _e \xi _o /U_e}$. At this streamwise position, the inflow condition is expressed as a superposition of the base state and instability waves,

(2.2)\begin{equation} \boldsymbol{q}_o=\boldsymbol{q}_B(\xi_0,\eta) + Re\left(\sum_m\sum_n c_{{{n,m}}} \hat{\boldsymbol{q}}_{{n,m}}(\eta) \exp[\mathrm{i}k_m \phi-\mathrm{i}\omega_n t + \mathrm{i} \varphi_{{{n,m}}}] \right), \end{equation}

where $c_{{{n,m}}}$ is the amplitude, $\hat {\boldsymbol {q}}_{{{n,m}}}$ is the discrete slow-mode profile, and $k_m$, $\omega _n$ and $\varphi _{n,m}$ are the azimuthal wavenumber, frequency and relative phase, respectively, for each $(n,m)$ instability wave.

The data assimilation attempts to identify the vector $\boldsymbol {c}=[\ldots,c_{{{n,m}}},\ldots ]^{\rm T}$ of amplitudes of these waves at the inlet. We assume the average spectrum is independent of relative modal phase ($\varphi _{n,m}$), which is reasonable given the random nature of free-stream tunnel forcing and the relatively long-time acquisition of the spectra. To reflect this assumption and lack of information especially for 3-D waves, each phase is independently selected from a uniform random distribution $\varphi _{n,m}=[-{\rm \pi}, {\rm \pi}]$ and remains fixed during the course of the data assimilation. In contrast, if the objective were to assimilate the instantaneous wall-pressure signals, the relative phases of the inflow modes would need to be included in the control vector and accurately estimated (Buchta & Zaki Reference Buchta and Zaki2021).

We consider 52 frequency and azimuthal wavenumber pairs that encompass Mack's first- and second-mode instabilities for the Reynolds-number range $1345 \leq \sqrt {Re_\xi } \leq 1872$. This size of the control vector was based on spectral analysis of the pressure-probe data to determine the dominant frequencies in the observations (see figure 3a) and preliminary testing. The frequencies and integer azimuthal wavenumbers are $f\in [50,75,\ldots,350]~\text {kHz}$ and $k\in [0,20,40,60]$, respectively. For each $(\,f,k)$ pair, we consider the most unstable slow modes (Fedorov Reference Fedorov2011) over the Reynolds-number range of interest, from the discrete Orr–Sommerfeld and Squire spectra. For the computation of the inflow instability profiles $\hat {\boldsymbol {q}}_{{{n,m}}}(\eta )$ only, the effect of curvature was neglected, i.e. we adopt the wavenumber $\beta _m=k_m / r_o$, where $r_o$ is the radius of the cone at the inflow, which is a reasonable assumption for this Reynolds number (Malik & Spall Reference Malik and Spall1991).

Figure 3. (a) Power spectra at measurement stations for the experiment (black circles), EnVar prediction (connected small blue squares) and linear prediction (red line). (b) Integral of power spectra. The confidence in each predicted observation is represented by $\pm 75\sigma _i$, where $\sigma _i$ is the ensemble standard deviation of each $i$th observation. (c) Cost and gradient during data assimilation.

Once the inflow condition (2.2) is prescribed, the downstream evolution in the DNS/LNS numerical simulations (see next paragraph for explanation of DNS and LNS) accounts for curvature fully. Since the experimental measurements lack 3-D information to guide the selection of relevant oblique modes, we consider a range of azimuthal wavenumbers that are linearly unstable ($k\leq 60$), and discretize this range in a manner that we can efficiently search and that enables the generation of important higher harmonics by nonlinear interactions. Beyond the inflow, the azimuthal discretization of the computations permits the formation of waves up to $k=540$, thus providing ample resolution for secondary instabilities that arise due to fundamental resonance between oblique and planar instability waves.

Most of the present simulations solve the compressible, nonlinear Navier–Stokes equations in curvilinear coordinates (termed DNS) for the simulation domain shown in figure 1. Linearized dynamics are acquired with exactly the same solver assuming a steady base state and neglecting nonlinear interactions only (termed LNS). LNS is used to provide an initial estimate of the unknown inflow spectra that will be refined to account for nonlinear effects; LNS will also be referenced in the discussion to contrast the disturbance field to the nonlinear evolution. The gas is assumed ideal with ratio of specific heats $\gamma = 1.4$ and temperature-dependent viscosity, following a power-law formula, $T^n$, where $n$ is determined from measured viscosity data between the boundary-layer edge and wall temperatures ($T_e=65$ K and $T_w=300$ K). The flow equations are solved using a standard fourth-order Runge–Kutta scheme and interior fourth-order finite differences. Near boundaries, stencils are biased and accuracy is reduced to second order. Second derivatives in the viscous terms are evaluated using repeated first derivatives, instead of discretizing the operator directly (Mattsson & Nordström Reference Mattsson and Nordström2004), which requires adding high-order numerical dissipation of short waves to stabilize the solution – for details, see Vishnampet, Bodony & Freund (Reference Vishnampet, Bodony and Freund2015). The computational domain size $(L_\xi, L_\eta,L_\phi )=(105~\text {mm}, 17.6~\text {mm}$, $36^\circ$) is discretized with $(N_s, N_y, N_\theta )=(751, 201, 108)$ grid points. The azimuthal size was chosen to extend beyond the peripheral probes at $s_4$ in order to reduce any artificial correlation by the periodic boundaries. In viscous units, the grid spacing along the wall at the inflow is ($\Delta \xi ^+$, $\Delta \eta ^+$, $r_o \Delta \phi ^+$) = (2.8, 0.1, 3.1). The azimuthal and streamwise grids are uniform, and the wall-normal grid spacing is stretched according to the transformation used in Pruett et al. (Reference Pruett, Zang, Chang and Carpenter1995).

For all of the simulations, the flow develops for 1.75 flow-through times, based on the edge velocity $U_e$, before data acquisition commences. During 40 $\mathrm {\mu }\text {s}$, wall pressure is recorded at every time step; the time-step size is chosen to ensure the Courant–Friedrichs–Lewy number $\text {CFL}\approx 0.4$ throughout the simulated time horizon. The finite-sensing area of the PCB probes is modelled by averaging the wall pressure, $\tilde {p}_j = (1/N)\sum _{i=1}^N p_i$, for $N$ grid points that satisfy the inequality $(x_j-x_i)^2+(y_j-y_i)^2+(z_j-z_i)^2 \leq (d_s/2)^2$, where $d_s$ is the sensing diameter. To take advantage of statistical homogeneity in the simulation, this operation is performed at the streamwise sensor position for all azimuthal grid points. The simple averaging used to model the sensing area in DNS was shown to improve adherence to measurements of pressure spectra and intensity (Huang et al. Reference Huang, Duan, Casper, Wagnild and Bitter2020). The data from the simulated probe pressure $\tilde {p}$ are processed and concatenated into two vectors in the same manner as the experimental measurements,

(2.3a,b)\begin{align} \boldsymbol{h}_S = [\hat{p}\hat{p}^\star(s_1,f_1), \ldots, \hat{p}\hat{p}^\star(s_4,f_{N_f})]^{\rm T}, \quad \boldsymbol{h}_I = \left[{\sum}_i\hat{p}\hat{p}^\star(s_1,f_i), \ldots,{\sum}_i\hat{p}\hat{p}^\star(s_4,f_i)\right]^{\rm T}. \end{align}

These vectors are used as inputs for the following data assimilation framework.

2.3. Data assimilation procedure

The unknown amplitudes $\boldsymbol {c}$ of instability waves are sought using an ensemble-variational (EnVar) framework (Buchta & Zaki Reference Buchta and Zaki2021). To initiate the nonlinear optimization, an estimate for the control vector $\boldsymbol {c}_0$ is determined using linear dynamics of the flow. This choice provides a more physics-based approximation than adopting white noise across the spectrum and, as such, a more accurate starting estimate for the subsequent nonlinear optimization. For the linear estimate, the control vector is computed by direct differentiation of

(2.4)\begin{equation} \mathcal{J}_l=\tfrac{1}{2}\|\boldsymbol{m}-\boldsymbol{\mathsf{L}}\boldsymbol{v}\|^2 + \frac{\varrho}{2}\|\boldsymbol{v}\|^2 \end{equation}

with respect to the weights and identifying $\boldsymbol {v}$ at a stationary point $\boldsymbol {\nabla }\!\mathcal {J}_l = 0$. Each column of matrix $\boldsymbol{\mathsf{L}}=[\ldots |\boldsymbol {l}_{f,k}|\ldots ]$ represents the wall observations acquired from an independent evolution of the $i$th instability wave at ($f,k$) using the linearized flow equations and an axisymmetric Blasius base state. Since the performance of the linear method is predicated on the validity of a linear assumption, only measurements from the first sensor position $s_1$, where harmonics of the large-amplitude waves are absent, are included for the estimate of $\boldsymbol {c}_0$. The second term in (2.4) represents a penalty, which has two desired effects: (i) it mitigates against high-energy inflow disturbances, and (ii) it improves the conditioning of the inverse problem. The value of $\varrho$ is chosen such that $\mathcal {J}_l/\mathcal {J}_{\boldsymbol {v}=0}=10^{-4}$, which ensures that the measurements are accurately reproduced to within $1$% (see § 3). If the experimental conditions only trigger a linear boundary-layer response, then minimizing (2.4) will be sufficient to accurately assimilate the measurements and predict the inflow amplitudes.

For the experiment considered in § 2.1, signatures of nonlinearity are present in the measurements. As a result, the initial estimate $\boldsymbol {c}_0$ of the control vector based on (2.4) is not sufficient, and the data assimilation requires an iterative nonlinear optimization. We consider the following cost,

(2.5)\begin{equation} {\mathcal{J}}={\tfrac{1}{2}\|{\log_{10}\boldsymbol{m}_S} -\log_{10}\boldsymbol{h}_S\|_{\varSigma_m^{{-}1}}^2} +{\tfrac{1}{2}\|\boldsymbol{m}_I-\boldsymbol{h}_I\|_{\varSigma_n^{{-}1}}^2} +\tfrac{1}{2}\|\boldsymbol{c}-\boldsymbol{c}_i\|_{\varSigma_c^{{-}1}}^2, \end{equation}

which comprises three parts. The first term is the logarithm of the pressure spectra. Since the spectra span several orders of magnitude, the logarithm ensures that the cost function is not solely focused on the spectral peak but rather targets the entire range of frequencies. The second term captures the importance of the spectral peak in determining the overall pressure intensity. The final term represents the degree of trust in a prior control vector ($\boldsymbol {c}_i$) to avoid large steps from the previous control vector. The choice of the cost function can affect performance of the data assimilation, and our experience has shown that (2.5) furnishes better accuracy in fewer iterations relative to, for example, pressure spectra on a linear scale.

An EnVar technique is used to update the control vector in order to minimize the cost function (2.5). Figure 2 shows a schematic of the algorithm. An updated estimate of the control vector $\boldsymbol {c}$ is sought from its previous value $\boldsymbol {c}_i$ using a weighted superposition of ensemble members $\boldsymbol {c}^{(j)}$, specifically $\boldsymbol {c}_{i+1}=\boldsymbol {c}_i + \boldsymbol{\mathsf{P}}\boldsymbol {w}$, where $\boldsymbol{\mathsf{P}}=[\ldots | \boldsymbol {c}^{(j)}-\boldsymbol {c}_i|\ldots ]$ is a matrix of the perturbation vectors. Each control vector, $\boldsymbol {c}_i$ and $\boldsymbol {c}^{(i)}$, is evolved with the governing equations, in this case the Navier–Stokes equations, to produce the spatio-temporal representation of the state $\boldsymbol {q}_i$ and $\boldsymbol {q}^{(j)}$. Model observations are then extracted from the mean and ensemble states using $\boldsymbol {h}_i=\mathcal {M}(\boldsymbol {c}_i)$ and $\boldsymbol {h}^{(j)}=\mathcal {M}(\boldsymbol {c}^{(j)})$; they are then compared with the experimental data in the cost function. An estimate of the local gradient and Hessian furnish the optimal weights $\boldsymbol {w}$ by assuming a stationary point ($\boldsymbol {\nabla }\!\tilde {\mathcal {J}}=0$). Multiple updates of the control vector can be performed in the descent direction, until the cost stagnates or starts to increase. At this point, the entire EnVar process is repeated, iteratively and till convergence. The number of ensemble members in the main algorithm was $N_e=10$. Each EnVar iteration thus involves $N_e+1=11$ simulations to acquire the local gradient, where the additional DNS corresponds to the mean of the ensemble. Previous characterizations of the performance of this method are provided elsewhere (Jahanbakhshi & Zaki Reference Jahanbakhshi and Zaki2019; Buchta & Zaki Reference Buchta and Zaki2021; Mons, Du & Zaki Reference Mons, Du and Zaki2021). The results of the linear and nonlinear optimization are presented in the next section.