Data Assimilation of Ideally Expanded Supersonic Jet Using RANS Simulation for High-Resolution PIV Data

Abstract

Data assimilation using particle image velocimetry (PIV) and Reynolds-averaged Navier–Stokes (RANS) simulation was performed for an ideally expanded supersonic jet flying at a Mach number of 2.0. The present study aims to efficiently reconstruct all the physical quantities in the aeroacoustic fields that match well with a realistic, experimentally obtained flow field. The two-dimensional, two-component PIV measurement was applied to the jet axis plane, and the time-averaged velocity field was obtained using single-pixel ensemble correlation. Two-dimensional axisymmetric RANS simulation using the Menter shear stress transport (SST) model was also performed, and the parameters of the SST model were optimized via data assimilation using the ensemble Kalman filter. The standard deviation of the observation noise $\sigma$, which is a parameter of the ensemble Kalman filter, is estimated by the previously proposed method (Nakamura et al., Low-Grid-Resolution-RANS-Based Data Assimilation of Time-Averaged Separated Flow Obtained by LES. Int. J. Comp. Fluid. Dyn., 2022), and its effectiveness was investigated for the first time. This method effectively estimated the magnitude of $\sigma$ at each generation without tuning the hyperparameters. The assimilated flow fields exhibited similar flow structures observed in PIV such as the potential core length or shear layer. Therefore, the present framework can be used to estimate time-averaged full flow fields that match well with experimentally observed flow fields, and has the potential to construct a database for the Navier-Stokes-based stability analysis that requires a full flow field.

1. Introduction

A supersonic jet, which is the exhaust flow from a propulsion device on a rocket or an airplane, generates strong acoustic waves due to the presence of complex flow fluctuations. These acoustic waves cause noise pollution and fatigue failure due to the acoustic loads. Therefore, the noise generation mechanism of supersonic jets has been studied in recent decades utilizing various experimental techniques and CFD/CAA (computational fluid dynamics/computational aeroacoustics) [1,2,3,4]. In light of the requirements for better understanding and modeling of jet noise, reduced-order models have been employed [5,6]. Among these models, a resolvent analysis has recently been attracting attention for applications in stability analysis in the aeroacoustic field [7,8,9,10,11]. Moreover, super-resolution measurements based on reduced-order models have recently proposed, and they reveal the unsteady behaviors of jets [12,13,14,15,16]. Super-resolution measurements rely on modal decomposition [17] such as proper orthogonal decomposition (POD) [18] or dynamic mode decomposition (DMD) [19,20], and thus, it is possible to replace them with resolvent modes. Here, resolvent analysis requires knowing the full physical quantities of flow fields because the resolvent is a linear operator derived from linearized Navier–Stokes equations. This feature forces the use of computational data for resolvent analysis, and thus, the resolvent analysis has not yet been applied to an experimentally measured supersonic jet flow.

The experimental techniques in the aeroacoustic field have evolved alongside improvements in the specifications of measurement devices and the development of superior analytic algorithms. For example, schlieren, particle image velocimetry (PIV), and pressure-sensitive paint (PSP) methods can visualize density gradient fields [21,22,23,24], velocity fields [25,26,27,28,29], and surface pressure fields [30], respectively. These experimental techniques are productive and enable us to easily conduct parametric studies. However, their measurement spaces and resolutions are dependent on the image sensor of the camera and the optics. In addition, the number of measurable physical quantities is often limited to one or two. These points make conducting resolvent analysis impossible using experimental data alone. On the other hand, CFD/CAA methods allow us to obtain all of the physical quantities in the aeroacoustic field with sufficiently high spatial resolution. In recent years, a method for predicting supersonic jet noise at a relatively high Reynolds number ($Re\sim {10}^6$) with an error of 2 dB using a large-eddy simulation (LES) was proposed [31]. This kind of high-fidelity computation of a supersonic jet generally requires higher-order schemes and a high grid resolution, resulting in huge computational costs [32,33,34]. Reynolds-averaged Navier–Stokes (RANS) simulations are favorable in terms of their computational costs and applicability to engineering fields, although they can only be used to calculate the averaged flow field, which is required for resolvent analysis. Semlitsch and Mihăescu [35] computed a supersonic jet using the Menter shear stress transport (SST) model [36] and compared it with the PIV results under the same conditions. The computed flow field exhibited similar shock structures, although the computation overestimated the turbulent viscosity, resulting in the rapid growth of the shear layer. Overestimations of the turbulent kinetic energy in the shear layer have also been observed in similar computations reported by Mishra and Iaccarino [37], and the uncertainty estimation for the turbulence model has been performed by [38]. Chauhan and Massa [39] performed a RANS analysis of a supersonic jet and a comparative analysis of the vorticity thickness for different turbulence models. They showed that the SST model with a compressibility correction is able to estimate the experimental data well. However, the results of those computations still highly depend on the inflow boundary condition or the turbulence model [40]. Therefore, obtaining aeroacoustic fields with sufficient prediction accuracy requires tuning the parameters, which often takes a lot of time.

In consideration of the research background described above, the present study aims to efficiently reconstruct all the physical quantities in an aeroacoustic field that match well with a realistic supersonic jet. This is achieved by integrating the experimental and computational results using a data assimilation technique. Data assimilation, which incorporates observed values into a model to improve the prediction accuracy, has recently been introduced to CFD [41]. Kato et al. [42] applied data assimilation using an ensemble Kalman filter to a RANS simulation around an airfoil. The experimentally measured pressure coefficients were used as the observations, and the parameters in a turbulence model were optimized to improve the prediction accuracy. Nakamura et al. [43] employed a similar framework to that by Kato et al. [42] for a low-grid-resolution RANS simulation of a square cylinder wake and obtained similar flow fields to those found in an LES analysis. Although previous studies successfully improved the prediction accuracy of wake flows from low-speed to transonic regimes, the method’s applicability to a highly compressible flow such as that from a supersonic jet has not been evaluated. In addition, there are several studies involving the optimization of a turbulence model using assimilated PIV data as observations [44,45]. Those studies employed subsonic jets, and their PIV data were processed by spatial correlation, which tends to decrease the spatial resolution. However, this is not suitable for high-Reynolds-number jets with thin shear layers, and the spatial resolution of the PIV data used for observations should be sufficiently high.

Therefore, the present study reconstructs the flow field of a supersonic jet by integrating PIV and RANS simulation using data assimilation. Single-pixel ensemble PIV was performed for an axisymmetric supersonic jet of a Mach number of 2.0, and a high-spatial-resolution mean velocity field, which resolves the thin shear layers, was obtained. Then, data assimilation based on the two-dimensional axisymmetric RANS simulation using the Menter SST model [36] was performed. Similar flow fields to the experimental ones were successfully reproduced by an ensemble Kalman filter after optimizing the SST model parameters. Here, the estimation method of the variable hyperparameter, previously proposed by [43], is introduced into the ensemble Kalman filter for robust data assimilation, and it is evaluated by comparing the fixed hyperparameter computations for the first time. Finally, the effectiveness of the present framework was discussed by illustrating the data-assimilated flow fields.

2. Jet Conditions

An ideally expanded supersonic jet at the Mach number of 2.0 was employed in the present study. The jet conditions were controlled by the nozzle pressure ratio and the designed Mach number of the nozzle. Figure 1 illustrates the cross-sectional geometry of an axisymmetric convergent–divergent nozzle for which the designed Mach number is 2.0. The nozzle’s exit diameter D is 10 mm, and the nozzle’s contour was designed by the method of characteristics. The nozzle pressure ratio is defined as the pressure ratio between the stagnation chamber and the ambient and can be calculated using the following equation:

$$\frac{p_c}{p_{\infty }}={\left[1+\frac{1}{2}(\gamma -1)M_j^2\right]}^{\frac{\gamma }{(\gamma -1)}}.$$

(1)

Here, $M_j$, $p_c$, $p_{\infty }$, and $\gamma$ are the jet’s Mach number at the nozzle’s exit, the stagnation pressure, the ambient pressure, and the specific heat ratio, respectively. The ideally expanded conditions were reproduced at the nozzle pressure ratio of 7.82, corresponding to a Mach number of 2.0. The Reynolds number at the nozzle’s exit was $1.0\times {10}^6$, and the temperature ratio between the stagnation chamber and ambient conditions was 1.0.

3. Methods

3.1. Experimental Apparatus

The two-dimensional, two-component (2D2C) PIV measurement was performed, and the time-averaged velocity field of the supersonic jet was obtained for the data assimilation. The PIV measurement was performed in an anechoic room equipped with a jet-generating device at Tohoku University. Readers can refer to Ozawa et al. [27,46] for more details about the experimental facilities. Figure 2 illustrates the experimental setup of the PIV measurement. The PIV measurement system is composed of a double-pulsed laser (LDY-303, Litron, Tokyo, Japan) and a high-speed camera (Phantom V611, Vision Research, Sydney, Austrilia). The hardware specifications and measurement conditions are summarized in

Table 1. Hardware specifications and measurement conditions.

Specifications of Laser:
Laser system	LDY-300PIV (Litron)
Laser type	Nd:YLF
Laser wavelength	527 nm
Pulse energy	2 mJ at 10 kHz
Laser sheet width	approx. 0.8 mm
Specifications of Camera:
High-speed camera	Phantom V611 (Vision Research)
Image sensor	$1280\times 800$ pixels
Pixel pitch	20 µm
Camera lens	Nikkor 80–200 mm f/2.8
Measurement conditions:
Measurement area	$160\times 50$ mm
Pixel resolution	$1280\times 400$ pixels
Spatial discretization	12.5 µm/pix
Time between laser pulses	1.2 µs
Sampling rate	1 kHz
Number of snapshots	20,000 pairs

. The seeding particles are generated using Laskin nozzles and a 50% aqueous solution. Both the jet flow and ambient air were filled with seeding particles. The short time intervals for the particle image pairs were set to be 1.2 µs, and the image pairs were acquired at the sampling rate of 1 kHz. The obtained particle images were analyzed using single-pixel ensemble correlation [46,47], and the pixel-by-pixel velocity fields were obtained. Single-pixel ensemble correlation calculates the cross-correlation using the change in time series intensity in a single pixel instead of a spatial interrogation window. Therefore, only the mean velocity field can be obtained, but the spatial resolution is high. This approach is suitable for data assimilation because the thin shear layer responsible for the acoustic generation can be properly resolved. The spatial resolution of the calculated velocity vectors was $1256\times 392$ vectors (12.5 µm/vector).

3.2. Numerical Apparatus

Two-dimensional axisymmetric RANS simulations were performed using the open source software SU2CFD (v8.0.1) [48] developed by Stanford University. The unstructured grid was generated by Gmsh software (4.12.2), and the boundary conditions are shown in Figure 3. The calculation area, defined as $0<x/D<30$ and $0<r/D<10$ for the outside of the nozzle and the inside of the nozzle, was also calculated, as shown in Figure 1. A prism mesh of 30 layers was applied inside the nozzle, and the minimum grid spacing inside the nozzle was 0.0001 D. The total grid number is approximately 40,000 points. The inflow conditions inside the nozzle were defined as $p_c=765$ kPa, $T_c=291$ K. The far-field conditions were defined as freestream. The freestream pressure and velocity were set to atmospheric pressure and 0.01 % of the jet velocity at the nozzle’s exit, respectively.

The Menter SST model was employed for the turbulence model. The SST model is a two-equation eddy viscosity model that combines the $k-\omega$ turbulence model and $k-ϵ$ turbulence model. The $k-\omega$ and $k-ϵ$ models are suitable for estimating the inner region of the boundary layer and the freestream region, respectively. The SST model switches those two models using a blending function. The closure coefficients for the $k-\omega$ model are $({\sigma }_{k1},{\sigma }_{\omega 1},{\beta }_1,a,{\beta }^{*},\kappa ,{\gamma }_1)$, and those for the $k-ϵ$ model are $({\sigma }_{k2},{\sigma }_{\omega 2},{\beta }_2,{\beta }^{*},\kappa ,{\gamma }_2)$.

Table 2. Default SST parameters.

${\mathit{\sigma }}_{\mathit{k}1}$	${\mathit{\sigma }}_{\mathit{k}2}$	${\mathit{\sigma }}_{\mathit{\omega }1}$	${\mathit{\sigma }}_{\mathit{\omega }2}$	${\mathit{\beta }}_1$	${\mathit{\beta }}_2$	${\mathit{\beta }}^{*}$	a	$\mathit{\kappa }$
0.85	1.0	0.5	0.856	0.075	0.0828	0.09	0.31	0.41

shows the default SST parameters. These default parameters cannot reproduce the experimentally obtained velocity fields well. Therefore, data assimilation was applied to the RANS simulation, and the SST parameters were optimized. Here, eight parameters of the SST model $({\sigma }_{k1},{\sigma }_{k2},{\sigma }_{\omega 1},{\sigma }_{\omega 2},{\beta }_1,{\beta }_2,{\beta }^{*},a)$ are the targets of the optimization in the present study. In addition, although several correction terms for compressible flows have been proposed for the SST model [49,50,51,52], the present study employed the basic SST model for a simple evaluation of the effectiveness of the data assimilation.The following section describes the details of the data assimilation technique used in the present study.

3.3. Data Assimilation

The present study employs the ensemble Kalman filter for data assimilation because this method is easy to incorporate into numerical simulations. Figure 4 illustrates the schematics of the data assimilation. The flow field calculated using the default SST parameters was used as the initial solution, and then randomly selected SST parameters were used for the first calculation. Then, the SST parameters and the updated flow fields were input into the ensemble Kalman filter, and the filtered SST parameters were calculated. In the next step, the filtered SST parameters were used for RANS simulation and the same process was repeated. In the present study, the ensemble number was set to be $N=16$, and the calculations were repeated 20 times.

The ensemble Kalman filter considers the state equation and observation equation of a discrete-time nonlinear system as follows:

$$\begin{array}{c}{\mathbf{x}}_t^{\left(l\right)}=f\left({\mathbf{x}}_{t-1}^{\left(l\right)}\right)\hfill \\ {\mathbf{y}}_t^{\left(l\right)}={\mathbf{H}}_t{\mathbf{x}}_t^{\left(l\right)}+{\mathbf{w}}_t^{\left(l\right)}\hfill \end{array}$$

(2)

Here, f represents a function that updates the state vector. This corresponds to the RANS calculation performed by SU2CFD in the present study. In addition, $\mathbf{y}$, $\mathbf{H}$, and $\mathbf{w}$ represent the observation vector, observation matrix, and observation noise, respectively. The index l is used for an ensemble member ($1<l<N$). The matrix $\mathbf{H}$ corresponds to the operator mapping the velocity field obtained by the RANS simulation onto the PIV measurement grid, as shown in Figure 5. Here, the state vector of the ensemble Kalman filter includes the SST parameters as well as the state and observation quantities of the flow field. The state vector ${\mathbf{x}}_t^{\left(l\right)}$ and the observation matrix $\mathbf{H}$ are defined as follows:

$$\begin{array}{ccc}\hfill {\mathbf{x}}_t^{\left(l\right)}& =& {\left[{\mathbf{x}}_{\mathrm{CFD}}^{\mathsf{T}},{\mathbf{x}}_{\mathrm{OBS}}^{\mathsf{T}},{\mathbf{x}}_{\mathrm{SST}}^{\mathsf{T}}\right]}^{\mathsf{T}}.\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill {\mathbf{H}}_t^{\left(l\right)}& =& \left[\mathbf{0},\mathbf{I},\mathbf{0}\right],\hfill \end{array}$$

(4)

Here, ${\mathbf{x}}_{\mathrm{CFD}}$ is a state vector that contains the density $\mathit{\rho }$; the momentum $\mathit{\rho }\mathit{u},\mathit{\rho }\mathit{v}$; the energy $\mathit{e}$; the turbulent energy $\mathit{k}$; and the dissipation rate $\mathit{\omega }$. ${\mathbf{x}}_{\mathrm{OBS}}$ has a velocity field that is mapped onto the PIV grid as the observation. Moreover, the state vector includes the SST parameters ${\mathbf{x}}_{\mathrm{SST}}$. Note that the SST parameters do not evolve over time, but alter when the filtering process is executed after obtaining the steady flow fields with the previous parameters.

$$\begin{array}{ccc}{\mathbf{x}}_{\mathrm{CFD}}\hfill & =\hfill & {\left[{\mathit{\rho }}^{\mathsf{T}},\mathit{\rho }{\mathit{u}}^{\mathsf{T}},\mathit{\rho }{\mathit{v}}^{\mathsf{T}},{\mathit{e}}^{\mathsf{T}},{\mathit{k}}^{\mathsf{T}},{\mathit{\omega }}^{\mathsf{T}}\right]}^{\mathsf{T}},\hfill \\ {\mathbf{x}}_{\mathrm{OBS}}\hfill & =\hfill & {\left[{\mathit{u}}_{\mathrm{PIV}}^{\mathsf{T}},{\mathit{v}}_{\mathrm{PIV}}^{\mathsf{T}}\right]}^{\mathsf{T}},\hfill \\ {\mathbf{x}}_{\mathrm{SST}}\hfill & =\hfill & {\left[{\sigma }_{k1},{\sigma }_{k2},{\sigma }_{\omega 1},{\sigma }_{\omega 2},{\beta }_1,{\beta }_2,{\beta }^{*},a\right]}^{\mathsf{T}}.\hfill \end{array}$$

(5)

In addition, $\mathbf{I}$ is the identity matrix of the size of the number of observations; the sizes of $\mathbf{0}$ inside ${\mathbf{H}}_t^{\left(l\right)}$ are appropriately determined so that ${\mathbf{H}}_t^{\left(l\right)}{\mathbf{x}}_t^{\left(l\right)}={\mathbf{x}}_{\mathrm{OBS}}$ is satisfied.

The filtering process of the ensemble Kalman filter is described as follows. The ensemble average of the state vector is calculated using the following equation:

$$\overline{{\mathbf{x}}_t}=\frac{1}{N}\sum _{l=1}^N{\mathbf{x}}_t^{\left(l\right)}.$$

(6)

The matrix $\delta {\mathbf{X}}_t$, which summarizes the differences of each ensemble from the ensemble average, and the variance–covariance matrix $\overline{{\mathbf{V}}_t}$ can be defined as follows:

$$\delta {\mathbf{X}}_t=\frac{1}{\sqrt{N-1}}\left[{\mathbf{x}}_t^{\left(1\right)}-\overline{{\mathbf{x}}_t},{\mathbf{x}}_t^{\left(2\right)}-\overline{{\mathbf{x}}_t},\cdots ,{\mathbf{x}}_t^{\left(N\right)}-\overline{{\mathbf{x}}_t}\right],$$

(7)

$$\overline{{\mathbf{V}}_t}=\delta {\mathbf{X}}_t\delta {\mathbf{X}}_t^{\mathsf{T}}.$$

(8)

The present study assumed that the instances of observation noise $\mathbf{w}$ are uniform at all the observation points and do not interfere with each other. Then, the variance–covariance matrix $\mathbf{R}$ of the observation noise is expressed as follows:

$$\mathbf{R}={\sigma }^2\mathbf{I}.$$

(9)

Here, $\sigma$ is the standard deviation of the observation noise, which is basically unknown. This is generally not updated after determining an appropriate value before repeated calculations. When observation noise is large, optimization calculations take time and the number of iterations increases, leading to an increase in calculation costs. Therefore, the present study applied a method proposed by Nakamura et al. [43] that can estimate the standard deviation of the observation noise using a low-rank state. The estimated ${\sigma }_{\mathrm{est}}$ is defined as follows:

$$\mathbf{H}\delta {\mathbf{x}}_t^{\left(l\right)}=\mathbf{U}\mathbf{\Sigma }{\mathbf{W}}^{\mathsf{T}},$$

(10)

$${\sigma }_{\mathrm{est}}^2=\frac{1}{N(N-1)}\left[{∥{\mathbf{U}}^{\mathsf{T}}\left(\mathbf{H}{\mathbf{x}}_t^{\left(l\right)}-{\mathbf{Y}}_t\right)∥}_F^2-{\parallel \mathbf{\Sigma }\parallel }_F^2\right].$$

(11)

Equation (10) represents the singular value decomposition, and Equation (11) is used for the estimation at each step. Here, the first term inside the bracket on the right-hand side of Equation (11) is the estimation of the sum of squares of the observation noise and the projection of covariance of members into the observations; the second term is a square of the projection of covariance of members into the observations. This simple subtraction process yields the estimation of the square of the observation noises. Refer to Nakamura et al. [43] for another description of this estimation method. In this study, this estimation for $\sigma$ is used unless otherwise mentioned, and a constant $\sigma$ is partially used for comparison.

Finally, the filtering process can be conducted using the following equation:

$${\mathbf{x}}_t^{\left(l\right)}={\mathbf{x}}_t^{\left(l\right)}+{\mathbf{K}}_t\left({\mathbf{y}}_t-\mathbf{H}{\mathbf{x}}_t^{\left(l\right)}+{\mathbf{w}}_t^{\left(l\right)}\right).$$

(12)

where ${\mathbf{K}}_t$ is the Kalman gain, and it is defined as follows:

$$\begin{array}{ccc}\hfill {\mathbf{K}}_t& =& \overline{{\mathbf{V}}_t}{\mathbf{H}}_t^{\mathsf{T}}{\left(\mathbf{R}+{\mathbf{H}}_t\overline{{\mathbf{V}}_t}{\mathbf{H}}_t^{\mathsf{T}}\right)}^{-1}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em} \hfill \end{array}$$

(13)

$$\begin{array}{ccc}& =& \delta \mathbf{X}\delta {\mathbf{X}}^{\mathsf{T}}{\mathbf{H}}_t^{\mathsf{T}}({\sigma }^{-2}\mathbf{I}-{\sigma }^{-4}{\mathbf{H}}_t\delta {\mathbf{X}}_t(\mathbf{E}\left.{\left.+{\sigma }^{-2}\delta {\mathbf{X}}_t^{\mathsf{T}}{\mathbf{H}}_t^{\mathsf{T}}{\mathbf{H}}_t\delta {\mathbf{X}}_t\right)}^{-1}\delta {\mathbf{X}}_t^{\mathsf{T}}{\mathbf{H}}_t^{\mathsf{T}}\right),\hfill \end{array}$$

(14)

where $\mathbf{E}$ is an identity matrix of the size of the number of ensemble members. In the present study, the Kalman gain in Equation (13) is expanded by the inverse matrix lemma, as shown in Equation (14), and memory usage and computational load are reduced by the multiplication calculations expressed in Equations (12) and (14) in order to avoid making a large matrix.

Once we obtain the new ensemble members, including turbulent model parameters, the next computations using the new turbulent model parameters can be conducted.

4. Results and Discussion

The error of the data assimilation is first investigated at different observation noises. The error is defined as the difference between the two components in the velocity field between PIV and RANS:

$$E=\frac{{\sum }_{j=1}^M{∥u_j^{\mathrm{PIV}}-u_j^{\mathrm{RANS}}∥}_2^2+{\sum }_{j=1}^M{∥v_j^{\mathrm{PIV}}-v_j^{\mathrm{RANS}}∥}_2^2}{{\sum }_{j=1}^M{∥u_j^{\mathrm{PIV}}∥}_2^2+{\sum }_{j=1}^M{∥v_j^{\mathrm{PIV}}∥}_2^2}.$$

(15)

Here, u and v are the streamwise and radial velocity components in the PIV or RANS results, respectively, and M is the total number of all observation points. Figure 6 shows the error of the data assimilation at each instance of observation noise. Note that the error in each generation is the average of all ensembles, and the error bar corresponds to the deviation in the ensembles. The error decreases with the increase in the number of generations, regardless of the difference in the standard deviation of the observation noise. If $\sigma$ is set to be large for all iterations, the convergence speed becomes slow. Moreover, a large $\sigma$ tends to induce a relatively high error, even if the number of iterations increases.

It should be noted again that we tried a smaller fixed $\sigma$ value than 500, but the computations often failed using this number because of the substantial changes in the turbulent model parameters. This implies that a fixed $\sigma$ of 500 is the close to the best choice for the fast convergence of computing the ensemble Kalman filter with the fixed hyperparameters. Setting the best $\sigma$ value for the ensemble Kalman filter requires a trial-and-error process in which filtering computations are calculated and investigated several times. The estimation method of $\sigma$ skips this trial-and-error process and maintains a robust computation with reasonable convergence speed based on only a single run of the filtering computation.

Figure 7 shows the relationship between the estimated $\sigma$ and the error. The value of ${\sigma }_{\mathrm{est}}$ is extremely high in the first generation because the first calculation uses randomly selected SST parameters. However, the magnitude of ${\sigma }_{\mathrm{est}}$ decreases as the number of generations increases, along with the error. This result indicates that the estimation method of Nakamura et al. [43] is effective in the present problem, and the data assimilation can be carried out without manually tuning the hyperparameters. It should again be noted that the evaluation of the estimation method of Nakamura et al. [43] is shown here for the first time.

The flow fields obtained using the optimized SST parameters are compared with the experimentally measured velocity field, and the effectiveness of the data assimilation is discussed. Here, the flow field with the minimum error set to $\sigma ={\sigma }_{\mathrm{est}}$ was compared with the PIV results. Figure 8 depicts the streamwise velocity distributions obtained with PIV and RANS. The velocity field obtained using the default SST parameters exhibits a short potential core length and a rapidly growing shear layer compared with those of the PIV. However, the optimized SST parameters suppress the shear layer’s growth and make the potential core length longer. These features are obvious in Figure 9, which shows a comparison of the streamwise velocity distributions. The wavy behavior of the PIV results in Figure 9a is due to the existence of weak shock waves. The reason why these weak shock waves are observed even under the ideally expanded conditions is due to the slight change in the actual nozzle diameter. Since the ideally expanded nozzle used in the present study does not consider the thickness of the boundary layer inside the nozzle, the boundary layer makes the actual nozzle diameter smaller, resulting in the appearance of shock waves. Note that the effect of the boundary layer’s thickness inside the nozzle [27] is hard to completely removefrom the experiment.Although the weak shock waves observed in PIV are not reproduced in the optimized RANS result, velocity fields similar to those in the PIV results can be obtained through data assimilation.

Figure 10 and Figure 11 depict the turbulence kinetic energy (TKE) obtained with PIV and RANS. Note that the TKE of the PIV data is not calculated using the single-pixel correlation. The single-pixel correlation is not able to estimate the turbulence properties because it relies on temporal information. Therefore, the TKE was alternatively estimated by the conventional spatial correlation using commercial PIV software (Dynamic Studio 6.7, Dantec Dynamics). Consequently, the spatial resolution of the TKE distribution of PIV is 0.1 mm, which is one-eighth of the distribution in the single-pixel ensemble correlation. The TKE in the PIV data is mainly distributed in the shear laye, r and its value is higher near the nozzle’s exit. This is due to the error vector calculated by the spatial correlation. Moreover, the turbulence intensity in the PIV data decreases at first and then increases downstream. This drop in turbulence intensity seems to be caused by the insufficient spatial resolution of the spatial correlation concerning the turbulence scale. Therefore, the turbulence intensity increases again downstream when the turbulence scale is sufficiently large for the used spatial resolution. Although the TKE in the PIV data includes errors and should thus be discussed carefully, the magnitude of the TKE on the downstream side ($10<x/D<15$) can be used as a reference. The RANS results with the default SST parameters overestimate the TKE in the shear layers, as reported in previous studies [35,37], compared with the TKE estimated in PIV data in the downstream side ($10<x/D<15$). The optimized SST parameters suppress the TKE in the shear layer and show a more accurate estimation, although the estimated TKE is still higher than that estimated in the PIV data, even in the downstream side ($10<x/D<15$). The optimized SST parameters are shown in

Table 3. Optimized SST parameters.

${\mathit{\sigma }}_{\mathit{k}1}$	${\mathit{\sigma }}_{\mathit{k}2}$	${\mathit{\sigma }}_{\mathit{\omega }1}$	${\mathit{\sigma }}_{\mathit{\omega }2}$	${\mathit{\beta }}_1$	${\mathit{\beta }}_2$	${\mathit{\beta }}^{*}$	a
0.474	0.574	0.393	0.634	0.070	0.102	0.080	0.434

. The optimized values of ${\sigma }_{k1}$, ${\sigma }_{k2}$, ${\sigma }_{\omega 1}$, and ${\sigma }_{\omega 2}$ are lower compared to the default ones. This indicates that the decrease in those parameters makes the estimated turbulent intensity lower.

Although the TKE is still high in the optimized RANS simulation, the present study indicates that the optimized RANS simulation enables us to obtain a time-averaged flow field similar to the observations. The time-averaged flow field is used as the so-called base flow in resolvent analysis, and data assimilation using PIV and RANS has the potential to construct a parametric database of base flows that reproduce the experimentally obtained flow fields.

5. Conclusions

The present study integrated PIV and RANS using data assimilation in order to obtain the full flow field of a supersonic jet that coincidences with the velocity field observed experimentally. The 2D2C PIV of a supersonic jet at a Mach number of 2.0 was performed, and time-averaged velocity fields were obtained. Then, the ensemble Kalman filter was employed for data assimilation, and the SST model parameters in the two-dimensional axisymmetric RANS simulation were optimized.

The method for estimating the hyperparameters of the ensemble Kalman filter is evaluated for the first time. The findings here illustrate that this technique enables us to conduct the ensemble Kalman filter’s computation without manually tuning the hyperparameters. We believe that this technique is very effective for the data assimilation of sensitive parameters such as the turbulent model parameters adopted in the present study.

Finally, the assimilated flow fields exhibited similar flow structures observed in the PIV model, such as the potential core length and shear layer. The presented results indicate that data assimilation using PIV and RANS can effectively estimate time-averaged full flow fields, which are the so-called base flows in resolvent analysis. Therefore, the present framework can be used to obtain a parametric database of the base flow, making resolvent analysis based on experimentally obtained flow fields possible.