4.1. Super-temporal-resolution reconstruction with down-sampled observations

Data assimilation is first applied for STR reconstruction using the observations of down-sampled LES data containing many small-scale turbulence structures. The computational cases are listed as cases 1–6 in table 1. The length of the intermediate DA window $20\Delta t$ equals 100 LES time steps. For cases 4–6, 50 DA windows are computed within a total time T corresponding to 9 turnovers of the largest-scale mode (Hussain & Zaman Reference Hussain and Zaman2006) according to the Strouhal number St = 0.3. Figure 2 presents the iteration residuals for different step-length strategies in clear and noisy cases. In the initial field DA phase as shown in figure 2(a), the instantaneous flow at a random instant is used as the first guess. The adaptive step lengths are computed using (2.19a,b). The residual has an increasing descent rate with respect to the initial step length, with only 10 iterations required using $\lambda _{\boldsymbol{s}}^0 = 500$. Nevertheless, the residual grows rapidly after reaching a minimum when a constant step length is used; this growth is well suppressed by the adaptive formulation. In the model-error DA phase as shown in figure 2(b), there is an increasing descent rate with respect to the initial step length. When adopting the adaptive scheme, a gradual increase in the step length according to (2.20a) is effective in accelerating convergence whereas a large step length is prone to cause divergence. This problem can be solved using the adaptive step length with the combined use of (2.20a) and (2.20b). Although this adaptive scheme relies on trial and error to determine the optimal initial step length, the complementary computations can be performed only in the first time step or the first DA window, yielding an increase in the computational cost of less than 1 % in total in this study. In noisy cases, the residual gradient with respect to the iteration number is better suited to convergence evaluation as the residual level is strongly associated with the regularisation coefficient $\alpha $. Obviously, larger $\alpha $ results in better de-noising effects in the DA procedure and thus greater discrepancy between the DA results and observations. As shown in figure 2(c), the residual gradient has the highest convergence speed in the DA-Noisy5e-5 case, reducing to below 10⁻⁵ after 10 rounds of iteration. The number of iterations is fixed at 50 in the model-error DA phase of all noisy cases.

Figure 2. Residuals and their gradients in DA computations: (a) initial field DA, (b) model error DA in the case of DA-Clear30 and (c) residual gradients in the model error DA of noisy cases.

In the STR reconstruction, the flow at the middle instant between the two observations is undoubtedly subject to the largest reconstruction error. This suggests the need to inspect the reconstructed flow snapshot at the middle instant in the case of DA-Clear30 as shown in figure 3. The clear LES data and the POD optimal reconstruction (Appendix C) from the down-sampled observations with time interval $30\mathrm{\Delta }t$ are also presented for comparison. The DA result is similar to the LES data on both longitudinal and cross-section planes as shown in figure 3(a–d). There are also slight differences in the upstream region near the inflow boundary, where some plume-like signatures are missing from the jet mainstream in the DA results. This is largely due to the error in the assimilated inflow conditions. In addition, the temporal evolution of these reconstructed fields shows the advantage of DA in STR reconstruction compared with interpolation-based approaches, which basically rely on temporally resolving the turbulence events with a sampling rate up to that meeting the Nyquist law. A POD full mode reconstruction (reconstruction using all of the modes) using the existing observations disrupts the real evolution process of the convecting small-scale vortices and is thus applicable only to slowly evolving vortices. A POD optimal reconstruction, as shown in figure 3(e,f), selects the most energetic modes that are temporally correlated for interpolation. This correctly recovers the flow convective properties but preserves less than 90 % of the total kinetic energy.

Figure 3. Instantaneous flow fields at the middle instant in the DA window: (a,b) DA reconstruction, (c,d) clear LES data and (e,f) POD optimal reconstruction. (a,c,e) The longitudinal middle sections and (b,d,f) the cross-sections at the location marked by the dashed line.

Figure 4 presents the pointwise error defined as

(4.1)\begin{gather}{\varepsilon _t} = {{\sum\limits_{(x,y,z)} {||\boldsymbol{u}(t) - {\boldsymbol{u}^{LES}}(t)||} } / {\sum\limits_{(x,y,z)} {||{\boldsymbol{u}^{LES}}(t)||} }},\end{gather}

(4.2)\begin{gather}{\varepsilon _x} = {{\sum\limits_{(t,y,z)} {||\boldsymbol{u}(x) - {\boldsymbol{u}^{LES}}(x)||} } / {\sum\limits_{(t,y,z)} {||{\boldsymbol{u}^{LES}}(x)||} }},\end{gather}

for each reconstruction scheme in clear cases. The summations are applied for all possible coordinates $(x,y,z)$ and $(t,y,z)$, respectively. In this application, the terminal observation is fixed at $t/\Delta t = 30$ whereas the initial observations are selected at $t/\Delta t = 0$, 10 and 20 for the cases of DA-Clear30, DA-Clear20 and DA-Clear10, respectively. The direct simulation by solving the NS equation for the initial field at $t/\Delta t = 0$ yields an error that increases monotonically with respect to time, as shown in figure 4(a). This error is largely induced by the boundary conditions imposed on the compact numerical domain and is absorbed by the model error in this study. Data assimilation is performed from each selected initial field and yields the largest error at the middle time of the DA window before reaching a minimum at the end. The error is slightly higher than zero at each start and end time owing to the finite residual level at the final iterative loop. With assimilation of the model error $\boldsymbol{\xi }$, the error of the DA reconstruction is appreciably lower than that of the direct simulation and maintains a clear descending trend with respect to the decreasing DA window length $\Delta T$. Additionally, the POD full modes and optimal reconstructions are compared with the DA results for each corresponding observational time interval. As the full mode reconstruction at an observational time (start or end of the DA window) produces the exact fields, the error is reduced to zero. However, the reconstruction error at the intermediate instants remains much higher than that of the DA results for each window length. The POD optimal reconstruction reduces the maximum error for large observational time intervals but performs even worse than the POD full mode reconstruction for small time intervals. The above results indicate that although POD reconstruction plays important roles in the feature recovery of large-scale evolution, the lack of temporal information on small scales in the raw observations remains a barrier for the temporal resolution enhancement of turbulence details. In the DA reconstruction, the error is mainly concentrated in the upstream region, where the jet speed is high, as shown in figure 4(b); the error is manifested as the convection lag of the flow events, as shown by the instantaneous vertical velocity distributions. In addition, the inflow and outflow boundary conditions have notable local effects. Nevertheless, most of the error is below 1 % in the case of DA-Clear30, indicating the high accuracy of the proposed DA scheme for STR reconstruction even with a long time interval.

Figure 4. Pointwise error of the reconstructed flow in each clear case: (a) domain integration at different instants and (b) cross-wise and temporal integrations at different downstream locations. Instantaneous vertical velocities are inset in (b) for comparison.

As an analogy to the measurement noise in a real experiment (e.g. PTV), the noisy observations are synthesised by adding white noise to the clear LES data according to (3.2) and box filtering using 3 × 3 × 3 elements on a DA grid. In these cases, DA is performed successively in 50 windows using the terminal fields of the previous window for initialisation. Detailed parameters are given in table 1. To assess the accuracy of the DA field throughout the domain, the cross-correlation coefficient for the streamwise velocity is calculated using the clear LES data as reference for each instant in the DA window:

(4.3) \begin{equation}{C_x}(t) = \frac{{\sum\limits_{\boldsymbol{x}} {[u(\boldsymbol{x},t)\boldsymbol{\cdot }{u^{LES}}(\boldsymbol{x},t)]} }}{{\sqrt {\sum\limits_{\boldsymbol{x}} {{{[u(\boldsymbol{x},t)]}^2}} \boldsymbol{\cdot }\sum\limits_{\boldsymbol{x}} {{{[{u^{LES}}(\boldsymbol{x},t)]}^2}} } }},\end{equation}

where $t \in [{t_0},{t_0} + \Delta T]$, and the summations are applied for each spatial location $\boldsymbol{x} = (x,y,z)$ on the DA grid. The cross-correlation coefficients for other velocity components are defined similar to (4.3). The final cross-correlation coefficients are averaged for each relative time with respect to the start instants in all the DA windows. The correlation coefficients for all the noisy cases, as well as the noisy observations, are shown in figure 5. The correlation coefficients of the noisy observations remain appreciably lower than those of all the DA cases owing to the noise contamination. This clearly indicates the strong capacity of DA to recover the turbulence properties even when taking noisy measurement data as observations. Comparatively, the regularisation term in the objective function not only eliminates noise but also plays an important role in improving the facticity of the reconstructed flow fields. The correlation coefficient is the highest at the middle instant of the DA window for the cases of DA-Noise5e-6 and DA-Noise5e-7 owing to the propagation of the governing equations. This trend is obviously different from that of the clear cases, where the error is largest in the middle of the window. When a smaller value of $\alpha $ is used, DA applies a larger weighting to the noisy observational data and thus introduces more error into the instantaneous field near the initial and terminal instants. We thus observe straighter correlation lines along time with larger $\alpha $, and the best performance with nearly constant correlations in the DA window is achieved in the case of DA-Noise5e-5. The optimal $\alpha $ is not readily obtained from figure 5 for different DA cases. The value of $\alpha $ is obviously case-dependent but can be cheaply determined adopting the DA procedure with only one window. In the remaining discussion of noisy cases, only DA-Noise5e-5 is taken as the representative example.

Figure 5. Cross-correlation coefficients of the reconstructed velocity fields with the referential LES data: (a) x component, (b) y component and (c) z component. The noisy LES data are averaged over all time. Different ordinate scales are used for clear illustration.

The assimilated $\boldsymbol{\xi }$ term acts as a compensator of error between the model prediction and observation. It also plays an important role in the turbulence evolution. The physical interpretability of the compensator has increasingly drawn attention for numerous data-driven techniques. The compensator not only directly reduces the predictive error between two datasets but also is a feature vital for the flow to evolve correctly. Term $\boldsymbol{\xi \; }$ has equivalent effects on flows with pressure gradients in the NS equations but is always divergence-free according to (2.17). Helmholtz decomposition casts $\boldsymbol{\xi \; }$ in two parts as

(4.4)\begin{equation}\boldsymbol{\xi } = \boldsymbol{\psi } - \boldsymbol{\nabla }\varphi ,\end{equation}

where $\boldsymbol{\psi }$ is stochastic forcing due to, for example, the subgrid stress. The reconstruction error is included in $\boldsymbol{\psi }$ in the present DA framework. Term $\boldsymbol{\nabla }\varphi $ is the irrotational part, which is lumped into the pressure gradient. The pressure and pressure gradient component in the x direction are first explored in figure 6. There is a visible difference in the pressure and its gradient distributions between the clear LES result and the other results. Simulation by solving the original NS equations directly does not reproduce the pressure and gradient correctly, as shown in figure 6(a–d). The quantitative comparison presented in figure 7 shows the model error effects on the confined computational domain. The results obtained by simulation and DA have high similarity, as shown in figure 6(c–f), with only slight discrepancy observed quantitatively in figure 7. The above results demonstrate the important contribution of the irrotational part $\boldsymbol{\nabla }\varphi $ in compensating for the pressure gradient discrepancies stemming from the compact computational domain. The results also indicate that further computation steps after DA are required to reconstruct the pressure fields correctly. An explicit calculation of $\boldsymbol{\nabla }\varphi $ and $\boldsymbol{\psi }$ can be readily made considering the difference in the pressure gradients between the LES and DA, i.e.

(4.5)\begin{gather}\boldsymbol{\nabla }\varphi = \boldsymbol{\nabla }{p^{LES}} - \boldsymbol{\nabla }{p^{DA}},\end{gather}

(4.6)\begin{gather}\boldsymbol{\psi } = \boldsymbol{\xi } + \boldsymbol{\nabla }{p^{LES}} - \boldsymbol{\nabla }{p^{DA}}.\end{gather}

Probability density functions (PDFs) of ${\xi _x}$, ${\psi _x}$ and ${\nabla _x}\varphi $ are plotted in figure 8, where ${\nabla _x}$ is the streamwise component of the gradient. All quantities are normalised by $U_0^2/D$. Herein, PDFs are separately calculated for all DA windows at the second instant immediately after the initial time (figure 8a), at the middle instant (figure 8b) and at the terminal instant (figure 8c) using data in the region $r/D < 1$ (where $r = \sqrt {{y^2} + {z^2}} $). The results of averaging across all instants are shown in figure 8(d). The curves for the negative part of the abscissa are mirrored to the positive part to demonstrate the asymmetry with respect to the ordinate axis. The PDFs of ${\xi _x}$ have a sharp peak at the origin and rapidly decay in both positive and negative directions of the abscissa. There are clear differences between the positive and negative bunches, especially at the first instant, as shown in figure 8(a). This indicates that ${\xi _x}$ has a pronounced bias towards the positive direction and complements the pressure gradient from the perspective of the whole region of interest. Recalling that direct simulation yields a lower speed of convection in the downstream direction as depicted in figure 4(b), ${\xi _x}$ is physically interpreted as the driving force for the model prediction to keep up with the observation evolution. However, the importance of the driving force varies in time and largely prevails in the first half of the DA window. This depends on ${\nabla _x}\varphi $, which is biased against the driving force. The PDFs of ${\xi _x}$ and ${\nabla _x}\varphi $ averaged across all the instants and shown in figure 8(d) clearly demonstrate the above mentioned complementation mechanism during the flow evolution in the DA window, although there is slightly anomalous behaviour at the terminal instant in figure 8(c), probably due to computational error. Supplemental evidence for this supposition is provided by the PDFs of ${\psi _x}$, which show good agreement for the positive and negative bunches at each instant in the figure. This result is consistent with our previous conjecture as the remaining part of ${\xi _x}$, when subtracting the irrotational part $\boldsymbol{\nabla }\varphi $, has unbiased stochastic behaviours that are responsible for the production and dissipation of turbulence.

Figure 6. Instantaneous pressure (a,c,e) and streamwise-component pressure gradient (b,d,f) fields in (a,b) clear LES, (c,d) simulation and (e,f) the case of DA-Noise5e-5.

Figure 7. Quantitative comparison of (a) the pressure and (b) its streamwise-component gradient along a streamwise line at $y/D = 0.5$.

Figure 8. The PDFs of the quantities ${\xi _x}$, ${\psi _x}$ and ${\nabla _x}\varphi $ at (a) the second instant, (b) the middle instant, (c) the terminal instant and (d) all instants. The data on the negative abscissa are mirrored to be positive. Here $\varTheta $ denotes ${\xi _x}$, ${\psi _x}$ or ${\nabla _x}\varphi $ as indicated by the arrows.

A direct view of the STR reconstruction is acquired from figure 9, where the temporal variations of the original and reconstructed velocities are plotted. The streamwise velocities at $x/D = 3.5D$ on the jet centreline are quantitatively compared in figure 9(a). The rapid variation of the clear LES results represents the small-scale coherent structures in the jet turbulence. The down-sampling of the observations (‘Noisy LES obs.’) only captures the large-scale events even when adopting spline interpolation, and the observational noise induces appreciable deviation from the clear LES data. However, the DA results recover the small-scale fluctuations well, even though slight discrepancies are still observed. A salient feature of DA is that the STR results do not strictly follow the variation of the observations, manifesting remarkable de-noising and correcting effects when erroneous observations are used. The spatiotemporal plots of the streamwise velocities at the cross-section $x/D = 3.5D$ are present in figure 9(b–d). The clear LES data in figure 9(b) reveal appreciable temporal fluctuations along with the small-scale variation showing apparent coherent tubular structures. The detailed small vortex organisations are seen in the vorticity contours. In the down-sampled observations shown in figure 9(c), although a smooth isosurface is obtained, the loss of details between successive snapshots is an issue. This situation is equivalent to direct temporal interpolation, which makes sense only when the original sampling frequency meets the requirement of the Nyquist sampling law. The DA reconstruction shows obvious advantages in the recovery of flow details even when using the noisy down-sampled observations. The DA results with the sampling rate equalling that of the clear LES data are shown in figure 9(d). The contours of the streamwise vorticity clearly demonstrate the reconstructed flow details. The red circles mark the moderate-scale vortices that exist in the clear LES data and are successfully reconstructed by DA but are lost in the down-sampled observations.

Figure 9. Temporal plots of the velocity and vortical structures (Q = 500) from the original and DA reconstructed flow fields at $x/D = 3.5D$: (a) streamwise velocity variation on the jet centreline, (b) clear LES data with the time interval $\Delta t = 0.001\ \textrm{s}$, (c) noisy LES observations with the time interval $\Delta T = 0.02\ \textrm{s}$, (d) STR reconstructions with the time interval $\Delta t = 0.001\ \textrm{s}$. The green dots and the dashed line in (a) denote the observational data and the spline fitting, respectively. The black dashed lines in (b–d) indicate the observational instants.

In-depth analysis of the small scales recovery of DA is conducted using the velocity spectra. Power spectral densities (PSDs) are computed using Welch's method (Welch Reference Welch1967) with all 1000 snapshots for each case. Although the window size of 200 (corresponding to 1.8 turnovers of the largest-scale flow fluctuation) with an overlap of 100 in the PSD computation does not provide a good estimation of the spectra at low frequencies due to the data size, the results are not visibly different from those obtained using larger windows. Our present analysis focuses on high-frequency structures with Strouhal numbers larger than 0.2, which can be well captured by the present PSD algorithm. Figure 10 shows the PSDs of the streamwise velocity at various radial locations. The PSDs of the time-sparse noisy LES data are lower than the others due to the normalisation using the root mean square of the velocity signals, and it is only possible to resolve vortical structures up to $St = 0.1$ with a relatively flat distribution due to the noise. The POD optimal reconstruction recovers the large-scale vortices well but yields a rapid decay of the PSDs at $St > 0.5$ due to the lack of small-scale information. The DA reconstruction produces the best results up to $St = 5$, although with a plateau higher than that for the clear LES data due to the noise in the observations. An observable feature of the PSDs at different radial locations is the underestimation of the DA reconstruction at $r/D = 0$ but a slight overestimation at $r/D = 1.5$. This feature is reasonable in the present computations due to the high noise ratio at far radial locations where the flow speed is low.

Figure 10. PSDs on the streamwise velocity at the downstream position $x/D = 3.5D$ at various radial locations: (a) $r/D = \textrm{0}$, (b) $r/D = 1$ and (c) $r/D = 1.5$.

4.2. Super-temporal-resolution reconstruction for tomo-PIV with small scales injection

Although tomo-PIV has become one of the most widely used volumetric measurement techniques in fluid mechanics, the compromise between the measurement domain size and the spatial resolution resulting from the cross-correlation (Lawson & Dawson Reference Lawson and Dawson2014; Naka et al. Reference Naka, Tomita, Shimura, Fukushima, Tanahashi and Miyauchi2016; Fiscaletti et al. Reference Fiscaletti, Ragni, Overmars, Westerweel and Elsinga2022) is the major issue that limits its application in certain circumstances. The present study evaluates the field correction for the low-sampling-rate tomo-PIV results through DA reconstruction, with much effort being given to the regeneration of the turbulence structures by small scales (SS) injection from the observations. It is worth noting that the small vortices cannot be injected from the initial condition as done by Li et al. (Reference Li, Zhang, Dong and Abdullah2019) because of the weak-constraint nature of the present DA algorithm and the time-sparse observations.

The case of DA-Tomo (see table 1) is first discussed on the basis of the snapshots selected at the terminal instant of the DA window. The comparison of the tomo-PIV data (observation) and the DA result presented in figure 11(a,b) shows excellent similarity in terms of the velocity contours and even the three-dimensional vortical structures (Q-criterion isosurface, Q = 500). However, a notable feature of the DA procedure is the prominent ability to remove the velocity divergence error, which is important for post-processing, such as vorticity computation and pressure determination, using the existing three-dimensional fields. A snapshot obtained by box filtering (denoted $\widetilde {\boldsymbol{u}^{LES}}$, filtered using 7 × 7 × 7 elements on the DA grid, which is close to the physical size of the tomo-PIV interrogation volume) from the clear LES data at a selected instant is shown in figure 11(c). The velocity contours barely show any difference between the different data, but the isosurface of the Q criterion in the box-filtered field is more coherent with many elongated vortex tubes distributed in the jet shear layer. This indicates that the cross-correlation process in tomo-PIV measurements involves far more than spatial filtering than what has been commonly understood before, and more complex effects should be considered.

Figure 11. Instantaneous velocity fields and vortical structures $(Q = 500)$ of (a) the synthetic tomo-PIV measurements, (b) DA reconstruction and (c) filtered LES data on the DA grid.

The turbulence fields of the reconstructed flow can be inspected in more depth by plotting the joint PDFs of the invariants ${Q^\ast }$ and ${R^\ast }$, which are defined as (Rishita & Girimaji Reference Rishita and Girimaji2022)

(4.7)\begin{gather}{Q^\ast } =- {\textstyle{1 \over 2}}{b_{ij}}{b_{ji}},\end{gather}

(4.8)\begin{gather}{R^\ast } =- {\textstyle{1 \over 3}}{b_{ij}}{b_{jk}}{b_{ki}},\end{gather}

where ${b_{ij}} = {A_{ij}}/||A||$. Here, the asterisk is added to the invariants to distinguish them from the vortex Q criterion. The velocity gradient tensor is ${A_{ij}} = \partial {u_i}/\partial {x_j}$, and $||\cdot ||$ denotes the Frobenius norm. The zero-discriminant lines are determined according to ${Q^{{\ast} 3}} + 27{R^{{\ast} 3}}/4 = 0$. It has been well established that the ${Q^\ast }$–${R^\ast }$ PDFs have a characteristic teardrop shape (i.e. a Vieillefosse tail) (Vieillefosse Reference Vieillefosse1984) in various turbulent flows, with a high probability of occurrence along the right discriminant line. This characteristic probability map represents the predominance of sheet-like structures in the fourth quadrant and that of vortex stretching in the second quadrant. In this study, ${Q^\ast }$–${R^\ast }$ maps are calculated using all the flow instantaneous fields in the region of $r/D < 0.75$. A map of the clear LES data is plotted using solid contours in each panel of figure 12 for comparison. Figure 12(a) shows that both the clear and filtered LES data have the well-known teardrop shape in the ${Q^\ast }$–${R^\ast }$ plane, whereas the filtered flow has a slightly slimmer tail indicating fewer sheet-like structures and less vortex compression effects. However, there is an appreciable difference in the tomo-PIV data, where the predominant occurrences of the high PDFs are located slightly above the discriminant lines, as shown in figure 12(b). Misalignment in tomo-PIV data has previously been observed (Khashehchi et al. Reference Khashehchi, Elsinga, Ooi, Soria and Marusic2010) but with high PDFs slightly below the discriminant lines. This misalignment is obviously case-dependent, whereas tomo-PIV invariably fails to capture the detailed turbulence dynamics. The DA reconstruction produces a ${Q^\ast }$–${R^\ast }$ map having a well-organised teardrop shape, which is shown in figure 12(c) and is more like the filtered LES results. This implies an improvement of the DA procedure in turbulence reconstruction even though the improved reconstruction can be hardly distinguished from the vortical structures, as shown in figure 11. Additionally, there is a slight difference, with a sharper Vieillefosse tail, relative to both the clear and filtered LES data. In contrast, Vieillefosse tails were observed to be shorter and blunter in the work of Li et al. (Reference Li, Zhang, Dong and Abdullah2019) owing to the use of different DA algorithms.

Figure 12. Joint PDFs of the velocity gradient invariables ${Q^\ast }$ and ${R^\ast }$. The solid contours show the clear LES data, while the coloured dashed lines show (a) filtered LES results, (b) tomo-PIV data and (c) the DA reconstruction. White dashed lines in the third and fourth quadrants represent the zero-discriminant lines ${Q^{{\ast} 3}} + 27{R^{{\ast} 3}}/4 = 0$. The PDFs are normalised by the mean values on the ${Q^\ast }$–${R^\ast }$ plane.

The filtered flow $\widetilde {\boldsymbol{u}^{LES}}$ is then subtracted from the clear LES fields to derive the small-scale structures, which can be added to the tomo-PIV observations for the SS injection. To remove the correlation between the injected small scales and the tomo-PIV snapshots, a time shift of 0.5 s (approximately 4.5 turnovers of the largest-scale jet oscillation) is applied for the observations in the cases of DA-SS and DA-SSN. This approach is analogous to real tomo-PIV work where small-scale information is injected from independent numerical data. Recalling that the SS injection from initial fields is subject to extremely expensive computation, taking 15 days to determine a single initial field in a long DA window as noted by Li et al. (Reference Li, Zhang, Dong and Abdullah2019), the observation injection in the present study requires less than 5 minutes to determine each optimal field on a desktop computer, with the number of computational grid elements being twice that used by Li et al. (Reference Li, Zhang, Dong and Abdullah2019). This SS injection is principally achieved by the adjoint source, which transfers the large-scale information from the observations to the primary flow quantities, with the system equations serving as the physical filter for the vortex correction and noise removal. A qualitative understanding of the DA reconstruction for small scales recovery is acquired from figure 13, which presents different observations and the DA results. It is noted that the Q value (3000) used to show the three-dimensional vortical structures is much larger than that used for figure 11. The tomo-PIV field and its DA result in figure 13(a,d) thus show few vortical structures and a low vorticity magnitude with much emphasis on the small scales in the cases of SS injection. Recalling that the regularisation coefficient $\alpha $ is selected to minimise the difference between the mid-time and terminal-time reconstructed fields and for de-noising in noisy cases, the DA reconstructions in figure 13(e,f) have a certain attenuation of the vorticity magnitude relative to the observations in figure 13(b,c) but are still considerably higher than those in the tomo-PIV field. These results are clear demonstrations of SS injection, with small-scale structures regenerated in the DA fields differing from those in the observations. The authenticity of these reproduced small-scale turbulence dynamics is thus evaluated subsequently.

Figure 13. Instantaneous vorticity magnitude contours and vortical structures (Q = 3000) of the observations (a–c) and DA reconstructions (d–f): (a) tomo-PIV observation, (b) tomo-PIV observation with SS injection, (c) tomo-PIV observation with SS injection and noise, (d) DA-Tomo, (e) DA-SS and (f) DA-SSN. The instant at the terminal of a DA window is taken.

The joint PDFs of the invariants ${Q^\ast }$ and ${R^\ast }$ of the observations and DA results are shown in figure 14. These invariant maps are computed according to (4.7) and (4.8). The invariant map of the clear LES data is plotted in each panel for comparison. It is observed that with SS injection, the tomo-PIV generates the teardrop shape of the invariant map much better than the original measurements (figure 14a). The Vieillefosse tails shift onto the right discriminant line but are still shorter than those of the clear LES data. The noise greatly deteriorates the map pattern such that the characteristic teardrop shape is not seen, as shown in figure 14(b). The induced error is prone to raise the high-PDF regions in the third and fourth quadrants above the discriminant lines, which is also seen in the original tomo-PIV data in figure 12(b). Despite the unphysical turbulence in the observations, DA is able to improve the dynamical feature substantially, resulting in the well-organised teardrop shapes of PDFs on the ${Q^\ast }$–${R^\ast }$ plane as shown in figure 14(c,d).

Figure 14. Joint PDFs of the velocity gradient invariables ${Q^\ast }$ and ${R^\ast }$. The solid contours show the clear LES data, while the coloured dashed lines show (a) tomo-PIV observations with SS injection, (b) tomo-PIV observations with SS injection and noise, (c) the case of DA-SS and (d) the case of DA-SSN. White dashed lines in the third and fourth quadrants represent the zero-discriminant lines ${Q^{{\ast} 3}} + 27{R^{{\ast} 3}}/4 = 0$. The PDFs are normalised by the mean value on the ${Q^\ast }$–${R^\ast }$ plane.

An important feature of velocity gradients in turbulence flows is the alignment of the vorticity along the eigenvectors of the strain-rate tensor ${S_{ij}}$ (Ashurst et al. Reference Ashurst, Kerstein, Kerr and Gibson1987). This alignment largely determines the vortex stretching term ${S_{ij}}{\omega _i}{\omega _j}$, which is the main source of the enstrophy growth. By examining the PDFs of the cosines of angles between the vorticity vector and the eigenvectors of ${S_{ij}}$ (i.e. $|{\hat{\boldsymbol{e}}_1}\boldsymbol{\cdot }\hat{\boldsymbol{\omega }}|$, $|{\hat{\boldsymbol{e}}_2}\boldsymbol{\cdot }\hat{\boldsymbol{\omega }}|$ and $|{\hat{\boldsymbol{e}}_3}\boldsymbol{\cdot }\hat{\boldsymbol{\omega }}|$, where the eigenvectors ${\hat{\boldsymbol{e}}_1}$, ${\hat{\boldsymbol{e}}_2}$ and ${\hat{\boldsymbol{e}}_3}$ are arranged in descending order of the corresponding eigenvalues ${\lambda _1} > {\lambda _2} > {\lambda _3}$, and $\hat{\boldsymbol{\omega }}$ denotes the normalised vorticity vector), the preferential alignment of $\hat{\boldsymbol{\omega }}$ with ${\hat{\boldsymbol{e}}_2}$ is observed. Checking the vorticity alignments in the different fields provides a useful means for assessing the effectiveness of the DA scheme in recovering the turbulence dynamic features. Figure 15 shows that the alignment of the vorticity vector with each eigenvector in the clear LES data has a tendency similar to that in a turbulent jet, as shown by Buxton & Ganapathisubramani (Reference Buxton and Ganapathisubramani2010), and to that in spatially developing flow, as shown by Gomes-Fernandes, Ganapathisubramani & Vassilicos (Reference Gomes-Fernandes, Ganapathisubramani and Vassilicos2014), with the preferential alignment of the vorticity being along the second eigenvector. For other data in figure 15, there are certain discrepancies of the alignments compared with the clear LES. The vorticity alignments in the case of DA-tomo-PIV have a tendency almost identical to that of the original tomo-PIV results except for the third eigenvector, where there is a slight improvement, as shown in figure 15(c). This suggests that the DA procedure hardly changes the vortex orientation when using the tomo-PIV data as observations due to an absence of small-scale vortices. The SS injection improves the vorticity alignment with the first eigenvector, as shown in figure 15(a), whereas the case of DA-SS further augments the alignments with the second and third eigenvectors, as shown in figure 15(b,c). We see here evidence for the dynamic feature enhancement of small-scale turbulence structures; this is more conspicuous when there is noise in the observations, as the alignment is improved substantially in the case of DA-SSN from a rather poor tendency in the noisy observations.

Figure 15. Alignment of the vorticity vectors with the eigenvectors of the strain-rate tensor. (a) Alignment with the first (extensive) eigenvector, (b) alignment with the second (intermediate) eigenvector and (c) alignment with the third (compressive) eigenvector. The ordinate scales differ across panels.

Another important feature of turbulence flows is the modulation of the small-scale vortices by the large-scale gradients. Many studies have shown that large-scale structures modulate small scales both in amplitude and frequency (Mathis, Hutchins & Marusic Reference Mathis, Hutchins and Marusic2009; Ganapathisubramani et al. Reference Ganapathisubramani, Hutchins, Monty, Chung and Marusic2012; Fiscaletti, Ganapathisubramani & Elsinga Reference Fiscaletti, Ganapathisubramani and Elsinga2015) rather than the scales being independent of each other as had long been considered. In jet flows, the small-scale signal has a stronger amplitude for positive fluctuations of the large-scale signal independently of the radial location (Fiscaletti et al. Reference Fiscaletti, Ganapathisubramani and Elsinga2015). There is also preferential alignment between the vorticities on the small and large scales, indicating that the vorticity direction does not vary appreciably across the scales (Fiscaletti et al. Reference Fiscaletti, Attili, Bisetti and Elsinga2016). This body of research suggests that the anisotropy of the large scales is partially preserved at the small scale, which seems to be in contrast with the local isotropy hypothesis (Pope Reference Pope2000). The small scales have been defined to be smaller than the Taylor length scale ${\lambda _T}$. Fiscaletti et al. (Reference Fiscaletti, Attili, Bisetti and Elsinga2016) reported that the scale modulation was weakened at larger scales but remained appreciable at scales up to $3{\lambda _T}$. Inspection of the scale modulation of the DA results is thus important in evaluating the DA scheme on the accurate recovery of the scale organisation in the cascade mechanism.

Following Fiscaletti et al. (Reference Fiscaletti, Attili, Bisetti and Elsinga2016), the amplitude modulation is assessed by defining the local vorticity root mean square $A_{gL}^\ast $ and the large-scale velocity gradient $g_L^\ast $ as

(4.9)\begin{gather}A_{gL}^\ast (\boldsymbol{x}) = \sqrt {\frac{1}{N}\sum\limits_{i = 1}^N {(||{\boldsymbol{\omega }_i} - \tilde{\boldsymbol{\omega }}||)} } ,\end{gather}

(4.10)\begin{gather}g_L^\ast (\boldsymbol{x}) = \frac{1}{N}\mathop \sum \limits_{i = 1}^N \sqrt {\left. {\frac{{\partial {{\tilde{u}}_x}}}{{\partial y}}} \right|_i^2 + \left. {\frac{{\partial {{\tilde{u}}_x}}}{{\partial z}}} \right|_i^2 + \left. {\frac{{\partial {{\tilde{u}}_y}}}{{\partial x}}} \right|_i^2 + \left. {\frac{{\partial {{\tilde{u}}_y}}}{{\partial z}}} \right|_i^2 + \left. {\frac{{\partial {{\tilde{u}}_z}}}{{\partial x}}} \right|_i^2 + \left. {\frac{{\partial {{\tilde{u}}_z}}}{{\partial y}}} \right|_i^2} .\end{gather}

Here, $A_{gL}^\ast $ is a measure of small-scale vortical structures whereas $g_L^\ast (\boldsymbol{x})$ only includes the shear components of the large-scale velocity gradients. Both quantities are normalised by the corresponding mean values. The tilde denotes the large-scale quantities in each dataset, and N is the number of nodes in the filtering window where the large scales are computed. Herein, the definition of large scales is not straightforward, as simply imposing filters does not give a reasonable and clear separation of the different scales in the present jet. In the clear LES, noisy LES and DA-Noisy5e-5 results, the filtered LES results obtained using a filter size of 7 × 7 × 7 elements on the DA grid (as shown in figure 11c) are taken as the large scales. The tomo-PIV original data are used as large scales for observations of SS injection and SS injection with noise. In the DA-SS and DA-SSN cases, the DA-Tomo results are used as large scales. With these definitions, the filter size is estimated to be approximately $3{\lambda _T}$ using the dissipation formulation of (3.1), and the scale modulation should still be appreciable according to the above discussion. Statistics of the activity at small scales conditioned on the large-scale gradients are shown in figure 16(a). Computations are performed in the region of $r/D \le 0.75$. The results show that the clear LES results have a nearly linear correlation as obtained by Fiscaletti et al. (Reference Fiscaletti, Attili, Bisetti and Elsinga2016) for $0.5 < g_L^\ast < 2.5$. This confirms that large-scale gradients strongly modulate small-scale structures. Even with a high level of noise in observations with attenuated correlation (see ‘Noisy LES’), this modulation is substantially improved in the case of DA-Noisy5e-5. The accuracy of the modulation recovery by DA is indeed region-dependent, as much better recovery can be observed (not shown here) by checking only the jet shear layer region at $r/D \approx 0.5$. It is seen that in the case of DA-SS, there is notable improvement relative to the SS-injected observations even with temporally sparse observational data. A loss of the scale modulation in the SS injection with and without noise is expected as the large- and small-scale structures are fully uncorrelated. Data assimilation improves the modulation properties regardless of the observational noise. This modulation recovery is an important criterion with which to assess the DA scheme in the reconstruction of real turbulence even when the small-scale structures are synthetic or contaminated with noise in measurement applications.

Figure 16. Modulation and alignment of vortices at different scales. (a) Local vorticity root mean square $A_{gL}^\ast $ conditioned on the fluctuations of the large-scale gradients $g_L^\ast $. (b) Alignment of large-scale vorticity vectors with small-scale vorticity vectors.

The modulation is also reflected in the vortex alignment between the small- and large-scale structures. The PDFs of the cosine of the angle between large and small scales (expressed as $|{\boldsymbol{\omega }_L}\boldsymbol{\cdot }{\boldsymbol{\omega }_S}|$) are plotted in figure 16(b). It is seen that in the clear LES data, the small-scale vortices have a preferential alignment with the large-scale vortices, as demonstrated by the rapid increase in the PDFs as $|{\boldsymbol{\omega }_L}\boldsymbol{\cdot }{\boldsymbol{\omega }_S}|$ approaches 1. We also see that this alignment is almost fully recovered even with noisy observations in the case of DA-Noisy5e-5. There is nearly no alignment correlation of different scales in the observations of SS injection and that with noise, with PDFs being almost horizontal in the range of $0 \le |{\boldsymbol{\omega }_L}\boldsymbol{\cdot }{\boldsymbol{\omega }_S}|\le 1$. This is further evidence that the small- and large-scale structures are independent in the original SS-injected observations. However, the vortex alignment is clearly improved by DA with the PDFs increasing at large values of $|{\boldsymbol{\omega }_L}\boldsymbol{\cdot }{\boldsymbol{\omega }_S}|$.

4.3. Lagrangian particle tracking through STR reconstruction

One of the motivations of the present study on STR reconstruction is to establish an auxiliary role for the successive prediction of particle positions in 4D-PTV. Four-dimensional PTV relies on the successive tracking of densely distributed particles at different instants to form long trajectories. This greatly increases the spatial resolution of the velocity vectors and improves the measurement accuracy by eliminating ghost particles that exist in tomo-PIV measurements, when using the same instrument configuration. The shake-the-box (STB) technique (Schanz, Gesemann & Schrder Reference Schanz, Gesemann and Schrder2016) is a representative 4D-PTV technology that predicts the new particle positions after the trajectory initialisation, followed by position perturbation for the iterative correction. Physically constrained interpolation schemes are adopted to reconstruct velocity fields on an Euler grid (Gesemann et al. Reference Gesemann, Huhn, Schanz and Schroder2016; Schneiders & Scarano Reference Schneiders and Scarano2016). However, a high frequency of camera imaging (i.e. small displacements of the particles at two successive sampling instants) is required to identify the same particles in the image sequence (i.e. particle pairing) (Khojasteh et al. Reference Khojasteh, Yang, Heitz and Laizet2021). This requirement limits the use of the STB technique to low-speed flows for currently available measurement devices. Lagrangian particle tracking largely depends on the particle prediction for the next instant before the position correction (shaking). A Wiener filter (Wiener Reference Wiener1949) is used in the prediction phase of the STB technique, with the accuracy of prediction deteriorating drastically for short trajectories and large particle displacement (Schanz et al. Reference Schanz, Gesemann and Schrder2016). Although a multi-pulse strategy using more cameras and lasers overcomes this difficulty by equivalently increasing the sampling rate and thus reducing the particle displacement (Manovski et al. Reference Manovski, Novara, Mohan, Geisler and Schrder2021), we prefer less expensive methods based on STR reconstruction for artificially increasing the imaging frequency. This relies on the accurate recovery of the detailed turbulence between the two given snapshots, which can be tentatively determined by either tomo-PIV or conventional PTV using the double-exposure mode of the cameras. Although the practice of this measurement strategy is left as a topic of future study, Lagrangian particle tracking using the existing flow fields is evaluated here as a preliminary demonstration.

Once the velocity fields are obtained, the new particle position ${\boldsymbol{x}_t}$ is computed from the previous position ${\boldsymbol{x}_{t - 1}}$ and previous local velocity ${\boldsymbol{u}_{t - 1}}$ as

(4.11)\begin{equation}{\boldsymbol{x}_t} = {\boldsymbol{x}_{t - 1}} + {\boldsymbol{u}_{t - 1}}\boldsymbol{\cdot }(\Delta t\ \textrm{or}\ \Delta T).\end{equation}

The time step of the DA computation $\Delta t$ or the time interval of the observations (DA window length) $\Delta T$ is used in (4.11) depending on the data used in different cases. This particle tracking formulation neglects various forces exerted on the particles by the fluid or adjacent particles; i.e. particles are distributed at a certain distance from each other and fully follow the flow as required for particle seeding in 4D-PTV measurements. We assume that the terminal field has been iteratively corrected in real experiments but includes a certain level of error. This implies that the terminal field does not need to be strictly accurate in using the present STR reconstruction. The particle prediction using the field of POD optimal reconstruction (Appendix C) and the Wiener prediction (Schanz et al. Reference Schanz, Gesemann and Schrder2016; Tan et al. Reference Tan, Salibindla, Masuk and Ni2019) are evolved for comparison.

Prediction errors are evaluated at the terminal instants in different DA windows using the clear DA cases. This evaluation demonstrates the best prediction that the DA can achieve using the present STR reconstruction strategy. To compute the error, particles are initialised at different radial and streamwise locations before they are convected downstream until the terminal instant in the DA window. The errors are computed as the distance from the predicted particles to the true particles at the terminal instant, and they are normalised using the voxel size in the tomo-PIV measurement as mentioned in § 3.2. Figure 17 shows the azimuthally averaged prediction errors for different lengths of the DA window. The largest error is found to be 13 voxels in the upstream region when the window length $\Delta T = 30\Delta t$ as shown in figure 17(a). This clearly demonstrates that the error decreases as the location moves downstream and the window size decreases. For a window length $\Delta T = 20\Delta t$ as shown in figure 17(b), the error is smaller than 6 voxels at all the considered locations, and the error becomes less than 2 voxels in the downstream half-domain when the window length is reduced to $\Delta T = 10\Delta t$ as shown in figure 17(c).

Figure 17. Particle prediction errors in the clear DA cases at the terminal instant: (a) window length $\Delta T = 30\Delta t$, (b) window length $\Delta T = 20\Delta t$ and (c) window length $\Delta T = 10\Delta t$. Results are averaged in the azimuthal direction.

Particle prediction is performed for each DA window independently in noisy cases as shown in figure 18; i.e. the particles are moved to the true positions (black dots on the true track) at each initial instant and the errors are evaluated at the terminal instant in a DA window. Data assimilation needs only two flow snapshots to compute the track whereas the Wiener prediction requires at least four previous positions for track initialisation and many more historical positions for improved accuracy. At time $t = 0$, two particles are initialised at $r/D = 0$ and $r/D = 0.5$ on the inflow surface. Figure 18(a) shows the seventh to the tenth true positions of the centreline particle. The true trajectory is computed using the clear LES data with time interval $\Delta t$ in (4.11), but the time interval between two successive black dots shown in figure 18(a) is $\Delta T = 20\Delta t$. The fifth-order Wiener prediction is applied using all the previous true positions with a time interval $\Delta T$, as we assume that only the down-sampled observations are known. It is seen that the Wiener prediction using the previous eight true positions determines the ninth position far from the truth; this also occurs for the prediction of the tenth position. The small red dots denote the DA prediction in the case of DA-Noisy5e-5 with short-time marching realised through STR reconstruction. This prediction is thus more accurate, being much closer to the true positions. Figure 18(b) presents the error evolution with time. It is seen that the error in the DA prediction decays with time as the particle convects downstream. The largest error is approximately 30 voxels in the upstream region. The error decay in the DA prediction is mainly due to the increase in accuracy of the STR reconstruction in the downstream region. The Wiener prediction has an error greater than 100 voxels in the upstream region and decays as more historical information becomes available. However, the error in the Wiener prediction remains higher than that in the DA prediction even at the last several instants. Indeed, the Wiener prediction improves when there are more historical data, but it is not applicable as the particles exit the measurement domain on excessively long trajectories. The DA procedure is thus preferable over the Wiener filter for particle prediction with long time intervals.

Figure 18. Particle prediction and position errors at different times: (a) three-dimensional tracking process and (b) errors at each terminal instant. The true track is computed using the clear LES data whereas the DA prediction is based on the case of DA-Noisy5e-5.

Figure 19 shows the particle prediction error for conditions closer to those of a real experiment subject to various flow measurement errors. The computational procedures are similar to those used in obtaining the results in figure 17 but involve further averaging across all the DA windows. The noisy LES and tomo-PIV results are computed using the time-sparse flow data with the time interval $\Delta T = 20\Delta t$ as if the high-sampling-rate fields were not available, whereas all DA predictions use the time interval $\Delta t$. The case of DA-Clear20 is reproduced here as a lower bound for the present STR strategy, and the POD optimal reconstruction, which yields much higher particle error, is also shown for comparison. The particle prediction using the tomo-PIV fields, which has been widely used in conventional PTV algorithms, has errors larger than 15 voxels in most of the measurement domain. The DA predictions in the cases of DA-SS and DA-SSN have errors similar to the error in the tomo-PIV results, whereas the DA prediction in the case of DA-Tomo has slightly less error owing to the STR reconstruction with a much smaller time step. Among all the discussed DA approaches, the case of DA-Noisy5e-5 is the most promising and has a particle prediction error much smaller than that of tomo-PIV used in conventional PTV algorithms. This DA case, with particle error smaller than 10 voxels, is thus much more suitable for Lagrangian particle tracking in a fluid with high particle density than the Wiener prediction. Additionally, this suggests that the uncorrelated SS injection does not benefit the particle tracking. Thus, correct small-scale structures should be included in the observations to improve the particle prediction; this can be done iteratively and is left for future work.

Figure 19. Particle prediction errors in different cases at the terminal instants of DA windows with $\Delta T = 20\Delta t$: (a) $x/D = 1$, (b) $x/D = 3$ and (c) $x/D = 5$. Results are averaged in the azimuthal direction and in different DA windows.