2.1 IHVDT analysis framework

The best analysis framework for IHVDT is not clear (see e.g. the review in Brocchini & Peregrine Reference Brocchini and Peregrine2001). When one fluid is dispersed throughout the other as is the case for the immiscible flows studied here, conditional phase-weighted averaging is helpful (Drew Reference Drew1983; Aliod & Dopazo Reference Aliod and Dopazo1990; Brocchini & Peregrine Reference Brocchini and Peregrine2001). The framework takes heterogeneous mixed flow and replaces it with a mixture of coexisting equivalent fluids. However, the technique tends to smooth out continuous phase fluctuations caused by the discrete phase (Lance & Bataille Reference Lance and Bataille1991). To assess the nature of IHVDT in its most basic form, we choose an unconditioned temporal averaging with Reynolds decomposition. As such, our notation is the instantaneous velocity field $U_{i}(\boldsymbol{x},t)=\overline{U}_{i}(\boldsymbol{x})+u_{i}^{\prime }(\boldsymbol{x},t)$ for $i=1\ldots 3$ , with $\overline{U}_{i}$ the temporal mean value and $u_{i}^{\prime }$ the random centred fluctuation such that $\overline{u_{i}^{\prime }\,}=0$ . The other flow quantities are total pressure $P=\overline{P}+p^{\prime }$ , density $\unicode[STIX]{x1D70C}=\overline{\unicode[STIX]{x1D70C}}+\unicode[STIX]{x1D70C}^{\prime }$ and viscosity $\unicode[STIX]{x1D707}=\overline{\unicode[STIX]{x1D707}}+\unicode[STIX]{x1D707}^{\prime }$ . With these definitions, the governing equations for the mean continuity (2.1), momentum (2.2) and kinetic energy (2.3) are

(2.1) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}t}}+{\displaystyle \frac{\unicode[STIX]{x2202}\left(\overline{\unicode[STIX]{x1D70C}}\overline{U}_{k}\right)}{\unicode[STIX]{x2202}x_{k}}}=-{\displaystyle \frac{\unicode[STIX]{x2202}\,}{\unicode[STIX]{x2202}x_{k}}}\left(\overline{\unicode[STIX]{x1D70C}u_{k}^{\prime }}\right), & & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}\,\overline{U}_{i}}{\unicode[STIX]{x2202}t}}+{\displaystyle \frac{\unicode[STIX]{x2202}\left(\overline{\unicode[STIX]{x1D70C}}\,\overline{U}_{i}\,\overline{U}_{j}\right)}{\unicode[STIX]{x2202}x_{j}}}={\displaystyle \frac{\overline{\unicode[STIX]{x1D70C}}}{Fr^{2}}}g_{i}-{\displaystyle \frac{\unicode[STIX]{x2202}\overline{P}}{\unicode[STIX]{x2202}x_{i}}}-{\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }u_{j}^{\prime }}}{\unicode[STIX]{x2202}x_{j}}}+{\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70F}}_{ij}}{\unicode[STIX]{x2202}x_{j}}}-A_{i},\nonumber\\ \displaystyle & & \displaystyle \qquad \text{where}\quad A_{i}={\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }}}{\unicode[STIX]{x2202}t}}+{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}}\left(\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }}\,\overline{U}_{j}+\overline{\unicode[STIX]{x1D70C}u_{j}^{\prime }}\,\overline{U}_{i}\right)\end{eqnarray}$$

and

(2.3) $$\begin{eqnarray}\displaystyle & & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\left(\frac{1}{2}\overline{\unicode[STIX]{x1D70C}}\,\overline{U}_{i}\,\overline{U}_{i}\right)}{\unicode[STIX]{x2202}t}}+{\displaystyle \frac{\unicode[STIX]{x2202}\left[\left(\frac{1}{2}\overline{\unicode[STIX]{x1D70C}}\,\overline{U}_{i}\,\overline{U}_{i}\right)\overline{U}_{j}\right]}{\unicode[STIX]{x2202}x_{j}}}={\displaystyle \frac{\overline{\unicode[STIX]{x1D70C}}\,\overline{U}_{i}}{Fr^{2}}}g_{i}-{\displaystyle \frac{\unicode[STIX]{x2202}\overline{P}\,\overline{U}_{i}}{\unicode[STIX]{x2202}x_{i}}}-A_{i}\overline{U}_{i}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,{\displaystyle \frac{\unicode[STIX]{x2202}\left(\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }u_{j}^{\prime }}\,\overline{U}_{i}\right)}{\unicode[STIX]{x2202}x_{j}}}+\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }u_{j}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{i}}{\unicode[STIX]{x2202}x_{j}}}+{\displaystyle \frac{\unicode[STIX]{x2202}\left(\overline{\unicode[STIX]{x1D70F}}_{ij}\,\overline{U}_{i}\right)}{\unicode[STIX]{x2202}x_{j}}}-\overline{\unicode[STIX]{x1D70F}}_{ij}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{i}}{\unicode[STIX]{x2202}x_{j}}}.\end{eqnarray}$$

In the mean governing equations, $\unicode[STIX]{x1D70F}_{ij}=\unicode[STIX]{x1D707}\frac{1}{2}(U_{i,j}+U_{j,i})/\mathit{Re}$ is the instantaneous viscous stress tensor and $g_{i}=-\unicode[STIX]{x1D6FF}_{i3}$ is the gravitational force. When written in this form, the different contributions of the turbulent mass flux ( $\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }}$ ) are explicit: (i) the source term on the right-hand side of (2.1); (ii) the $A_{i}$ term in (2.2); and (iii) the $A_{i}\,\overline{U}_{i}$ term in (2.3). In addition to understanding the overall turbulence characteristics in the mixed-phase region of the transom stern wake, we focus on investigating the relative importance of these TMF terms for IHVDT closure modelling. For completeness, we briefly discuss the use and implications of mass-weighted (or Favre) averaging in appendix C.

We average each case over a quasi-steady time period of $t=UT/D=40$ with a sampling rate of $\text{d}t=0.01$ , or 4000 samples. We estimate the statistical error in the mean velocity, density, r.m.s. velocity and density fluctuations, turbulent mass flux $\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }}$ and Reynolds stresses $\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }u_{j}^{\prime }}$ , using standard deviations of the entire sample time and multiple sub-time blocks of the simulation, to obtain 95 % confidence in these mean quantities (Friedberg & Cameron Reference Friedberg and Cameron1970; Morales, Nuevo & Rull Reference Morales, Nuevo and Rull1990). The maximum statistical error for all quantities considered is $\pm 5.3$  %.

2.2 Definition of the mixed-phase region ${\mathcal{R}}$ and conditioning averaged quantities

We define the mixed-phase region ${\mathcal{R}}$ to be the variable density region where the average density satisfies $\unicode[STIX]{x0394}\overline{\unicode[STIX]{x1D70C}}\leqslant \overline{\unicode[STIX]{x1D70C}}\leqslant 1-\unicode[STIX]{x0394}\overline{\unicode[STIX]{x1D70C}}$ , where the boundary width $\unicode[STIX]{x0394}\overline{\unicode[STIX]{x1D70C}}$ is an unknown parameter. We choose the value of $\unicode[STIX]{x0394}\overline{\unicode[STIX]{x1D70C}}$ such that ${\mathcal{R}}$ : (i) contains much of the relevant quantities such as turbulent kinetic energy and turbulent mass flux; (ii) is well behaved; and (iii) does not contain noise generated by individual entrainment or spray events (caused by the 1000:1 water-to-air density ratio). Note that when defined in this manner and combined with the unconditioned time average, ${\mathcal{R}}$ includes surface fluctuations as well as entrained air and spray (see figure 1).

As the higher density water motion can have a significant influence on $\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }}$ , we perform a careful statistical analysis to determine $\unicode[STIX]{x0394}\overline{\unicode[STIX]{x1D70C}}$ . The full analysis in appendix A shows that $\unicode[STIX]{x0394}\overline{\unicode[STIX]{x1D70C}}=0.05$ allows for comparable contributions around the mean interface. Appendix A also contains the definitions of the conditional moments and conditioned averages for the flow. Unless specified, quantities plotted in the forthcoming sections are depth-conditioned averages (A 2) at a location $(x,y)$ or area-conditioned (A 3) at a location $x$ .

2.3 Mean flow properties within the mixed-phase region in the wake

Figure 3(a) shows transverse ( $y$ – $z$ plane) cuts of $\overline{\unicode[STIX]{x1D70C}}$ and ${\mathcal{R}}$ in the wake at successive streamwise $x$ locations. The pockets of entrained air are visible at the locations of quasi-steady wave breaking in the converging-corner-wave and diverging-wave regions. The surface fluctuations and spray areas above the mean interface $\overline{\unicode[STIX]{x1D70C}}=0.5$ are also evident. Generally, at each wake location the $z$ extent of ${\mathcal{R}}$ is directly related to the breaking and spray formation. When associated only with surface fluctuations, the $z$ extent of ${\mathcal{R}}$ is small. Figure 3(b) shows $S^{{\mathcal{R}}}$ (see (A 3) for notation) along the wake for all three cases. Cases A and B are very similar and case C contains the largest mixed-phase region. The scaling of the mixed-phase region along the wake with $Fr$ is consistent with the convective nature of the flow at these ship speeds.

Figure 3.

Transverse cuts of (a) $\overline{\unicode[STIX]{x1D70C}}$ for case C and (b) conditioned transverse area $S^{{\mathcal{R}}}$ for

 case A;

 case B; and

 case C along the wake.

Figure 4 shows the conditioned mean velocity for all three cases at two $\hat{x}$ locations (transverse cuts) and along the wake (streamwise). In the rooster-tail region $\hat{x}<\hat{x}_{B}$ (figure 4 a,d,g), the velocity profiles collapse with beam scaling, consistent with the ballistic behaviour identified in Part 1. The wake deficit magnitude $U-\overline{U}_{1}$ is proportional to $B$ . The magnitudes of the transverse and vertical velocities exhibit no distinct scaling with geometry at this location. We note here that there exists some asymmetry with $y$ in the mean velocity profiles (most notably case A e.g. figure 4 h) and r.m.s. fluctuations in figure 5, which we hypothesize are a function of a buffer domain boundary rather than physical phenomenon or a function of statistical convergence.

Figure 4.

Values of $\langle {\overline{U}_{i}}^{{\mathcal{R}}}\rangle _{z}$ and $\langle {\overline{U}_{i}}^{{\mathcal{R}}}\rangle _{yz}$ for

 case A;

 case B; and

 case C. (a,d,g) transverse cut $\hat{x}=1.5$ ; (b,e,h) transverse cut $\hat{x}=4.0$ ; and (c,f,i) along the wake.

In the diverging-wave region $\hat{x}>\hat{x}_{B}$ (figure 4 b,e,h), the wake deficit magnitude peaks on the outer boundaries of ${\mathcal{R}}$ , decreases and then increases again on the centreline $y/B=0$ . This is in direct contrast with the general transverse spreading of the wake deficit observed in bluff-body wakes (Pope Reference Pope2000) due to the presence of the breaking waves on the outer boundaries of ${\mathcal{R}}$ . When we consider the aggregate plot along the wake (see figure 4 c,f,i), we see that the streamwise and transverse velocities show no significant dependence on geometry. However, the vertical velocity shows an inverse relationship with $B$ , in particular near $\hat{x}=1$ . This is due to the presence of the converging corner waves as identified in Part 1. This location represents (for cases A and B) where the converging jets collide on the centreline, generating a large vertical velocity. For case C, the converging corner waves break near $\hat{x}\approx 0.5$ . Thus, the converging jets are weaker at $\hat{x}\approx 1$ generating a decreased vertical velocity.

3.1 Root mean square fluctuations

Figure 5 plots transverse cuts of the r.m.s. velocity $u_{i}^{rms}=\sqrt{\overline{u_{i}^{^{\prime }2}}}$ and r.m.s. density fluctuations $\unicode[STIX]{x1D70C}^{rms}$ . Generally, the transverse cuts across the rooster tail (figure 5 a) contain a near-homogeneous behaviour within $|y/B|\leqslant 1$ for velocity components and peak at the centreline for the density components. Transversely across the diverging wave (figure 5 b), the fluctuations have peaks at $|y/B|\sim 3$ (caused by the quasi-steady breaking waves) that frame a near-homogeneous behaviour within $|y/B|<2$ . Along the wake (see figure 5 c), all of the fluctuations peak near $\hat{x}=1$ . The rooster-tail region $1<\hat{x}<\hat{x}_{B}$ shows near constant density fluctuations with a decrease in the velocity fluctuations. In the diverging-wave region $\hat{x}>\hat{x}_{B}$ both density and velocity fluctuations increase due to the presence of breaking waves. There exists a subtle difference in the magnitudes of the velocity fluctuations (both transversely and along the entire wake) implying a degree of anisotropy.

Figure 5.

Root mean square velocity and density fluctuations. (a) Transverse cut at $\hat{x}=1.5$ ; (b) transverse cut at $\hat{x}=4.0$ ; (c) along the wake.

  $u_{1}^{rms}$ ;

  $u_{2}^{rms}$ ; 

  $u_{3}^{rms}$ ; and ○  $\unicode[STIX]{x1D70C}^{rms}$ . Vertical 

 represents from left to right $\hat{x}_{c}$ and $\hat{x}_{B}$ , respectively. Data are case A.

3.2 IHVDT anisotropy

Figure 6 shows the IHVDT Reynolds stress $\unicode[STIX]{x1D70E}_{ij}=\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }u_{j}^{\prime }}$ and the anisotropy tensor components $\unicode[STIX]{x1D623}_{ij}=\unicode[STIX]{x1D622}_{ij}/2k=\unicode[STIX]{x1D70E}_{ij}/2k-1/3\unicode[STIX]{x1D6FF}_{ij}$ along the wake (here $\unicode[STIX]{x1D622}_{ij}=\unicode[STIX]{x1D70E}_{ij}-2/3k\unicode[STIX]{x1D6FF}_{ij};k=1/2\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }u_{i}^{\prime }}$ ). The diagonal components of the Reynolds stresses are negative and peak near $\hat{x}_{c}$ . The magnitude decreases in the rooster tail and then increases again in the diverging-breaking-wave region $\hat{x}>\hat{x}_{B}$ . The diagonal components of the anisotropy tensor (see figure 6 c) show that only the transverse and vertical components of the velocity fluctuations contribute to the turbulent kinetic energy (TKE) immediately off of the stern ( $\unicode[STIX]{x1D623}_{22}>0$ and $\unicode[STIX]{x1D623}_{33}>0$ at $\hat{x}\gtrsim 0$ ). The streamwise component quickly becomes dominant $0.25<\hat{x}<\hat{x}_{c}$ and the energy is in both the streamwise and transverse directions near $\hat{x}_{c}$ . The transverse and vertical components exchange relative importance in the rooster tail ( $\unicode[STIX]{x1D623}_{11}>\unicode[STIX]{x1D623}_{33}\sim 0$ and $\unicode[STIX]{x1D623}_{22}<0$ ). This reverses itself further downstream in the diverging-wave region ( $\unicode[STIX]{x1D623}_{11}>\unicode[STIX]{x1D623}_{22}\sim 0$ and $\unicode[STIX]{x1D623}_{33}<0$ ). The off-diagonal terms of the Reynolds-stress tensor (figure 6 b) are an order of magnitude smaller than their diagonal counterparts. Figure 6(d) shows $\unicode[STIX]{x1D623}_{13}$ to be the most significant anisotropy tensor term along the entire wake.

Figure 6.

Reynolds-stress tensor $\unicode[STIX]{x1D70E}_{ij}=\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }u_{j}^{\prime }}$ and anisotropy tensor $\unicode[STIX]{x1D623}_{ij}$ along the wake. (a) Diagonal components $\unicode[STIX]{x1D70E}_{ii}$ : 

  $\unicode[STIX]{x1D70E}_{11}$ ; 

  $\unicode[STIX]{x1D70E}_{22}$ ; and

  $\unicode[STIX]{x1D70E}_{33}$ . (b) Off-diagonal components $\unicode[STIX]{x1D70E}_{ij\,i\neq j}$ : 

  $\unicode[STIX]{x1D70E}_{12}$ ; 

  $\unicode[STIX]{x1D70E}_{13}$ ; and

  $\unicode[STIX]{x1D70E}_{23}$ . (c) Diagonal components $\unicode[STIX]{x1D623}_{ii}$ : 

  $\unicode[STIX]{x1D623}_{11}$ ; 

  $\unicode[STIX]{x1D623}_{22}$ ; and

  $\unicode[STIX]{x1D623}_{33}$ . (b) Off-diagonal components $\unicode[STIX]{x1D623}_{ij\,i\neq j}$ : 

  $\unicode[STIX]{x1D623}_{12}$ ; 

  $\unicode[STIX]{x1D623}_{13}$ ; and

  $\unicode[STIX]{x1D623}_{23}$ . Vertical 

 represents from left to right $\hat{x}_{c}$ and $\hat{x}_{B}$ , respectively. Data are case A.

A measure of the flow isotropy from the invariants of the anisotropy tensor ( $I=\unicode[STIX]{x1D623}_{ii}$ ; $II=\unicode[STIX]{x1D623}_{ij}\unicode[STIX]{x1D623}_{ji}$ ; $III=\unicode[STIX]{x1D623}_{ij}\unicode[STIX]{x1D623}_{jk}\unicode[STIX]{x1D623}_{ki}$ ) is the two-dimensionality or anisotropy parameter ${\mathcal{J}}=1{-}9(1/2II{-}III)$ , where ${\mathcal{J}}=0$ is highly anisotropic and ${\mathcal{J}}=1$ is isotropic (Jovanovic Reference Jovanovic2004). In the rooster tail (figure 7 a), the flow is moderately isotropic on the wake centreline, reaching a value of 0.6–0.8. However, it is highly anisotropic at $|y/B|\approx 1$ with a secondary peak value of 0.4 for $y/B>1$ . In the diverging-wave region (figure 7 b), the area of moderately isotropic behaviour widens with the wake width. However, the area of increased anisotropy (while not as strong) still persists at $|y/B|\sim 1$ and $|y/B|\sim 4$ . Figure 7(c) shows the aggregate anisotropy parameter plotted along the wake for all three cases. The flow generally changes from highly anisotropic to moderately isotropic downstream with a general collapse of all cases in the far wake to a similar value of ${\mathcal{J}}\approx 0.6$ . The flow never achieves fully isotropic behaviour within $O(10)D$ behind the stern.

Figure 7.

Anisotropy parameter ${\mathcal{J}}$ . (a) Case A transverse cut $\hat{x}=1.5$ ; (b) case A transverse cut $\hat{x}=4.0$ ; and (c) along the wake for

 case A;

 case B; and

 case C.

The method of componentality contours (Emory & Iaccarino Reference Emory and Iaccarino2014) provides a technique of visualizing the anisotropy of complex flows similar to the anisotropy invariant map (Lumley & Newman Reference Lumley and Newman1977) and barycentric map (Banerjee et al. Reference Banerjee, Krahl, Durst and Zenger2007) based on the invariants and their eigenvalues $\unicode[STIX]{x1D706}_{i}$ . The componentality contours construct a colour map from the barycentric coordinates

(3.1a,b ) $$\begin{eqnarray}\displaystyle x_{B}=C_{1c}+{\textstyle \frac{1}{2}}C_{3c},\quad y_{B}=C_{3c};\quad \text{with}\quad C_{1c}=\unicode[STIX]{x1D706}_{1}-\unicode[STIX]{x1D706}_{2}\quad \text{and}\quad C_{3c}=3\unicode[STIX]{x1D706}_{3}+1,\qquad & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D706}_{i}$ is defined such that $\unicode[STIX]{x1D706}_{1}\geqslant \unicode[STIX]{x1D706}_{2}\geqslant \unicode[STIX]{x1D706}_{3}$ . Note here that the Euclidean domain mapping chosen by Banerjee et al. (Reference Banerjee, Krahl, Durst and Zenger2007) simplifies such that $C_{2c}=2(\unicode[STIX]{x1D706}_{2}-\unicode[STIX]{x1D706}_{3})$ is not necessary for these definitions. The colour map translates $C_{n}$ , $n=(1c,2c,3c)$ to RGB (red green blue) colour components for visualization

(3.2) $$\begin{eqnarray}\displaystyle [R~G~B]^{\text{T}}=C_{1c}[1~0~0]^{\text{T}}+C_{2c}[0~1~0]^{\text{T}}+C_{3c}[0~0~1]^{\text{T}} & & \displaystyle\end{eqnarray}$$

(for further colour map details see Emory & Iaccarino Reference Emory and Iaccarino2014). Overlaying the colour map onto flows of geometric complexity provides a powerful tool for identifying local regions of anisotropy.

Figures 8 and 9 show the componentality contours in the wake of the stern. As shown in the legend in figure 9(a), red represents 1-component turbulence $x_{1C}$ , green 2-component turbulence $x_{2C}$ and blue 3-component turbulence $x_{3C}$ . Outside of the mixed-phase region (see figure 9), the mainly 1-component bulk water turbulence becomes 2-component and isotropic in the diverging-wave region. The spreading wake of the geometry in the bulk air turbulence is also visible. Within the mixed-phase region, the converging corner waves inject 1-component turbulent flow into the rooster-tail region (see figures 8 and 9 b). The rooster-tail centreline at $\hat{x}_{c}$ (figure 9 b) as well as the scars formed by the breaking waves in the diverging-wave region (figure 9 c) are sources of 2-component turbulence in ${\mathcal{R}}$ . Also shown in figure 9(c), is that the wave breaking draws 1-component bulk water turbulence up into ${\mathcal{R}}$ .

Figure 8.

Componentality contours on $\overline{\unicode[STIX]{x1D70C}}$ isosurfaces. (a) $\overline{\unicode[STIX]{x1D70C}}=0.95$ ; (b) $\overline{\unicode[STIX]{x1D70C}}=0.75$ ; (c) $\overline{\unicode[STIX]{x1D70C}}=0.25$ ; (d) $\overline{\unicode[STIX]{x1D70C}}=0.05$ . Data are case A. See figure 9(a) for legend.

Figure 9.

Componentality contours at transverse locations. Data are case A. (a) Legend; (b) $\hat{x}=1.4$ ; (c) $\hat{x}=3.9$ . In each image, upper white line is $\overline{\unicode[STIX]{x1D70C}}=0.05$ , lower white line $\overline{\unicode[STIX]{x1D70C}}=0.95$ .

3.3 Turbulent kinetic energy

Figure 10(a) shows the mean turbulent kinetic energy $k=\frac{1}{2}\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }u_{i}^{\prime }}$ along the wake for all three cases. The initial peak of the TKE occurs at $\hat{x}\sim 1$ , as expected. For case C, there are two peaks in the CCW, the first at $\hat{x}\sim 0.5$ and then at $\hat{x}_{c}$ . This is due to the fact that the converging corner waves break prior to the collision point for this case. The TKE decreases in the rooster tail and then increases again in the diverging-wave region. Cases B and C have an additional local peak at $\hat{x}\sim 2.75$ .

The mean turbulent kinetic energy obeys the exact transport equation in variable density flows (Chassaing et al. Reference Chassaing, Antonia, Anselmet, Joly and Sarkar2002)

(3.3) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}k}{\unicode[STIX]{x2202}t}}+{\displaystyle \frac{\unicode[STIX]{x2202}(k\overline{U}_{j})}{\unicode[STIX]{x2202}x_{j}}} & = & \displaystyle -{\displaystyle \frac{\unicode[STIX]{x2202}\left(\frac{1}{2}\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }u_{i}^{\prime }u_{j}^{\prime }}\right)}{\unicode[STIX]{x2202}x_{j}}}-\underbrace{\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }u_{j}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{i}}{\unicode[STIX]{x2202}x_{j}}}}_{\text{a}}+\underbrace{{\displaystyle \frac{\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }}}{Fr^{2}}}g_{i}}_{\text{b}}\nonumber\\ \displaystyle & & \displaystyle -\,{\displaystyle \frac{\unicode[STIX]{x2202}\overline{p^{\prime }u_{i}^{\prime }}}{\unicode[STIX]{x2202}x_{i}}}+\overline{p^{\prime }{\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }}{\unicode[STIX]{x2202}x_{i}}}}+{\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70F}_{ij}^{\prime }u_{i}^{\prime }}}{\unicode[STIX]{x2202}x_{j}}}-\overline{\unicode[STIX]{x1D70F}_{ij}^{\prime }{\displaystyle \frac{\unicode[STIX]{x2202}u_{i}^{\prime }}{\unicode[STIX]{x2202}x_{j}}}}-\underbrace{\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }}\left({\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{i}}{\unicode[STIX]{x2202}t}}+\overline{U}_{j}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{i}}{\unicode[STIX]{x2202}x_{j}}}\right)}_{\text{c}}.\qquad\end{eqnarray}$$

Equation (3.3) does not include the effect of surface tension which is not included in the iLES. Dodd & Ferrante (Reference Dodd and Ferrante2016) established that the power of surface tension can act either as a source or sink of TKE, and that for turbulent $We<5$ , the change of the surface energy on the droplet due to droplet oscillation and coalescence can be significant and any increase in TKE energy is irreversibly lost to viscous dissipation. Thus, bubbles with small turbulent Weber number will impact turbulent dissipation of the carrier fluid and should be incorporated in any explicit sub-grid-scale turbulence model. However, as the smallest bubbles we resolve have $We$ several orders of magnitude greater than this critical value, the TKE analyses presented here do not account for such effects.

Figure 10(b) shows the terms of (3.3) for case A. To determine a best estimate of the turbulent dissipation for this iLES, we calculate $\overline{\unicode[STIX]{x1D70F}_{ij}^{\prime }\unicode[STIX]{x2202}u_{i}^{\prime }/\unicode[STIX]{x2202}x_{j}}$ and scale it by $\unicode[STIX]{x1D708}_{e}^{-1}$ value from table 1. We apply the same scaling to the viscous transport term $\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70F}_{ij}^{\prime }u_{i}^{\prime }}/\unicode[STIX]{x2202}x_{j}$ . The production by mean shear and turbulent dissipation peaks near $\hat{x}_{c}$ . The transport terms are also significant in the converging-corner-wave region, as expected. Further downstream, the production by mean shear and turbulent dissipation are the dominant terms. However, the convection and transport terms continue to make small contributions in this region.

Figure 10.

(a) Mean turbulent kinetic energy $k$ along the wake for

 case A;

 case B; and

 case C. (b) Steady-state terms of (3.3) for case A

 term a;

 dissipation;

 convection;

 all transport terms;

 term b;

 term c. Vertical

 represents from left to right $\hat{x}_{c}$ and $\hat{x}_{B}$ , respectively.

The terms directly associated with the variable density turbulent flows are terms b and c – as they explicitly involve the turbulent mass flux $\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }}$ . We define the production of IHVDT turbulence ${\mathcal{P}}$ as the sum of that due to traditional production by the mean velocity gradients (term a) and the TMF terms b and c of (3.3). Thus,

(3.4) $$\begin{eqnarray}\displaystyle {\mathcal{P}}=\overline{\unicode[STIX]{x1D70C}}~\overline{u_{i}^{\prime }u_{j}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{i}}{\unicode[STIX]{x2202}x_{j}}}+\overline{\unicode[STIX]{x1D70C}^{\prime }u_{i}^{\prime }u_{j}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{i}}{\unicode[STIX]{x2202}x_{j}}}-\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }}\left({\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{i}}{\unicode[STIX]{x2202}t}}+\overline{U}_{j}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{i}}{\unicode[STIX]{x2202}x_{j}}}\right)+{\displaystyle \frac{\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }}}{Fr^{2}}}g_{i}. & & \displaystyle\end{eqnarray}$$

The first two terms arise from decomposing the density into its mean and average components in term a of (3.3).

Figure 11 shows the total production of TKE and its sub-components along the wake for case A and C, respectively. For case A, production due to mean shear and mean density ${\mathcal{P}}_{1}$ is responsible for most of the TKE production in the converging corner waves. The production from mean shear with variable density ${\mathcal{P}}_{2}$ and gravity effects ${\mathcal{P}}_{4}$ provides the rest of the production in this region, and ${\mathcal{P}}_{3}$ actually decreases the total TKE production moderately. For the case where the waves break within the converging-corner-wave region (case C), the gravity effects are as important (if not more) than the production by mean shear ${\mathcal{P}}_{1}+{\mathcal{P}}_{2}$ . For both case A and case C, the local peak of TKE at $\hat{x}\sim 2.75$ corresponds to a peak in the TMF specific term ${\mathcal{P}}_{3}$ and the mean shear term ${\mathcal{P}}_{1}$ .

Figure 11.

Value of ${\mathcal{P}}$ and its sub-components ${\mathcal{P}}_{i},i=1\ldots 4$ for (a) case A and (b) case C.

  ${\mathcal{P}}$ ; 

  ${\mathcal{P}}_{1}$ ; 

  ${\mathcal{P}}_{2}$ ; 

  ${\mathcal{P}}_{3}$ ; and

  ${\mathcal{P}}_{4}$ . Vertical

 represents $\hat{x}_{c}$ and $\hat{x}_{B}$ from left to right, respectively.

5.1 Model coefficients

Retaining the even-numbered terms of (5.1) creates a general 4 term model for a vector quantity, implying a minimum 12 potential values (assuming a global value will be adequate along the entire wake). Determining the model coefficients is a two-part process. The first part identifies the most relevant terms in the model by calculating a conditioned correlation coefficient at a wake location $x$ between the model term (in the absence of $C_{n}$ ) and a subset of the iLES data within ${\mathcal{R}}$ . The second part determines the coefficients using a linear-least squares fit between the relevant model terms and the iLES data.

With the available data, we select the subset of the iLES data to be case A and case B, leaving case C for a priori testing in § 5.2. Note that the estimate for $\unicode[STIX]{x1D716}$ is the same as that in § 3.3. Analysis of the conditioned correlation coefficients along the wake yields that a regional model focusing on the near wake and far wake has the best potential for success. For the near-wake region (defined here as including the converging-corner-wave and rooster-tail regions, $\hat{x}<\hat{x}_{B}$ ), the mean density gradient terms of the model (GTH) best correlate with the iLES data, namely $C_{2}$ for the $y$ - and $z$ -component and $C_{4}$ for the $x$ -component. Thus, we propose the near-wake explicit algebraic turbulent mass flux model (NW) as

(5.2) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\overline{\unicode[STIX]{x1D70C}u_{1}^{\prime }}=C_{4}{\displaystyle \frac{k^{2}}{\unicode[STIX]{x1D716}^{2}}}\left(\overline{u_{i}^{\prime }u_{k}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{j}}{\unicode[STIX]{x2202}x_{k}}}+\overline{u_{j}^{\prime }u_{k}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{i}}{\unicode[STIX]{x2202}x_{k}}}\right){\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}x_{j}}},\\[12.0pt] \overline{\unicode[STIX]{x1D70C}u_{2}^{\prime }}=C_{2}{\displaystyle \frac{k}{\unicode[STIX]{x1D716}}}\overline{u_{i}^{\prime }u_{j}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}x_{j}}},\\[12.0pt] \overline{\unicode[STIX]{x1D70C}u_{3}^{\prime }}=C_{2}{\displaystyle \frac{k}{\unicode[STIX]{x1D716}}}\overline{u_{i}^{\prime }u_{j}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}x_{j}}}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

All terms in this explicit model take the form: forcing $\times$ turbulent diffusivity. For this NW model, the forcing is the mean density gradient. The turbulent diffusivity for the streamwise component is the rate of production of stresses ${\sim}\overline{u_{i}^{\prime }u_{k}^{\prime }}(\unicode[STIX]{x2202}\overline{U}_{i}/\unicode[STIX]{x2202}x_{k})$ and for the transverse and vertical components the diffusivity is the actual stresses.

For the far-wake region (defined here as $\hat{x}>\hat{x}_{B}$ ), both the mean density gradients ( $C_{2}$ and $C_{4}$ ) and buoyancy terms ( $C_{6}$ and $C_{8}$ ) correlate well and the proposed far-wake model (FW) is

(5.3) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}\overline{\unicode[STIX]{x1D70C}u_{1}^{\prime }}=C_{4}^{1}{\displaystyle \frac{k^{2}}{\unicode[STIX]{x1D716}^{2}}}\left(\overline{u_{i}^{\prime }u_{k}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{j}}{\unicode[STIX]{x2202}x_{k}}}+\overline{u_{j}^{\prime }u_{k}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{i}}{\unicode[STIX]{x2202}x_{k}}}\right){\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}x_{j}}}+C_{6}{\displaystyle \frac{\overline{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x1D716}}}\overline{u_{i}^{\prime }u_{j}^{\prime }}g_{j}\\[12.0pt] \quad \qquad +\,C_{8}^{1}\overline{\unicode[STIX]{x1D70C}}{\displaystyle \frac{k}{\unicode[STIX]{x1D716}^{2}}}\left(\overline{u_{i}^{\prime }u_{k}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{j}}{\unicode[STIX]{x2202}x_{k}}}+\overline{u_{j}^{\prime }u_{k}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{i}}{\unicode[STIX]{x2202}x_{k}}}\right)g_{j},\\[12.0pt] \overline{\unicode[STIX]{x1D70C}u_{2}^{\prime }}=C_{2}^{2}{\displaystyle \frac{k}{\unicode[STIX]{x1D716}}}\overline{u_{i}^{\prime }u_{j}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}x_{j}}}+C_{4}^{2}{\displaystyle \frac{k^{2}}{\unicode[STIX]{x1D716}^{2}}}\left(\overline{u_{i}^{\prime }u_{k}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{j}}{\unicode[STIX]{x2202}x_{k}}}+\overline{u_{j}^{\prime }u_{k}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{i}}{\unicode[STIX]{x2202}x_{k}}}\right){\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}x_{j}}},\\[12.0pt] \overline{\unicode[STIX]{x1D70C}u_{3}^{\prime }}=C_{2}^{3}{\displaystyle \frac{k}{\unicode[STIX]{x1D716}}}\overline{u_{i}^{\prime }u_{j}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}x_{j}}}+C_{8}^{3}\overline{\unicode[STIX]{x1D70C}}{\displaystyle \frac{k}{\unicode[STIX]{x1D716}^{2}}}\left(\overline{u_{i}^{\prime }u_{k}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{j}}{\unicode[STIX]{x2202}x_{k}}}+\overline{u_{j}^{\prime }u_{k}^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{U}_{i}}{\unicode[STIX]{x2202}x_{k}}}\right)g_{j}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The forcing in the FW model still involves mean density gradients and then adds buoyancy effects to the streamwise and vertical components. The turbulent diffusivities of the streamwise component buoyancy effects are the Reynolds stresses and their rate of production. The vertical component turbulent diffusivity is only the rate of production of the stresses. For the transverse component, the FW model only depends on the mean density gradient forcing with turbulent diffusivity coming from both the stresses and their rate of production.

We determine the model coefficients for (5.2) and (5.3) using a linear-least squares fit at each $\hat{x}$ where all terms for that model correlate using the subset data cases A and B. We define the transition between the models based on $\hat{x}_{B}$ in table 1. We note here that this transition location is currently a value obtained from pre-existing simulation data and determining an appropriate scaling estimate of this location is ongoing work. However, knowing the scaling of the salient features in the wake with only the ship geometry and speed obtained in Part 1 and the convective nature of the flow established in §§ 3–4 suggests an estimate of this transition location to be $O(2.5{-}3)Fr$ behind the stern. Table 2 shows the resulting model coefficients.

5.2 Model performance

To assess the feasibility of the model to explicitly predict the TMF from resolved quantities, we perform a priori tests of the NW and FW model. Specifically, we use the iLES data of case C and model coefficients of table 2 to compute the model terms in (5.2) and (5.3). For notation purposes, we will name this the ‘MODEL’.

Table 2.

Corresponding model coefficients fit from the subset iLES data for (5.2) and (5.3).

Figures 15 and 16 show the MODEL prediction of the TMF to the actual iLES data TMF at two transverse locations. The highly localized nature of the TMF is evident in all three components (see figure 15 a,c,e). In general, MODEL does not predict some of the streamwise TMF at the core of the rooster tail (located $(y/B,z)\approx (\pm 0.5,0.4)$ in figure 15) and all three components near the quasi-steady wave breaking in the FW (located $(y/B,z)\approx (\pm 1.5,-0.2)$ in figure 16). This aside, the explicit model predicts the localization and magnitude extremely well at both wake locations. Table 3 contains the conditioned correlation coefficient ${\mathcal{C}}$ and conditioned normalized standard deviation $\unicode[STIX]{x1D70E}_{N}$ for the transverse locations in figures 15 and 16. In particular, the model captures the transverse and vertical TMF associated with the surface fluctuations and spray areas in the near-wake region and the streamwise TMF in the far-wake region remarkably well.

Figure 15.

Transverse cuts at $\hat{x}=2.0$ for TMF of (a,c,e) iLES data and (b,d,f) MODEL prediction. (a,b) $\overline{\unicode[STIX]{x1D70C}u_{1}^{\prime }}$ ; (c,d) $\overline{\unicode[STIX]{x1D70C}u_{2}^{\prime }}$ ; (e,f) $\overline{\unicode[STIX]{x1D70C}u_{3}^{\prime }}$ . Correlation coefficients and normalized standard deviations in table 3.

Figure 16.

Transverse cuts at $\hat{x}=4.0$ for TMF of (a,c,e) iLES data and (b,d,f) MODEL prediction. (a,b) $\overline{\unicode[STIX]{x1D70C}u_{1}^{\prime }}$ ; (c,d) $\overline{\unicode[STIX]{x1D70C}u_{2}^{\prime }}$ ; (e,f) $\overline{\unicode[STIX]{x1D70C}u_{3}^{\prime }}$ . Correlation coefficients and normalized standard deviations in table 3.

Table 3.

Conditioned correlation coefficient ${\mathcal{C}}$ and conditioned normalized standard deviation $\unicode[STIX]{x1D70E}_{N}=\unicode[STIX]{x1D70E}_{MODEL}/\unicode[STIX]{x1D70E}_{iLES}$ for the transverse $\hat{x}$ locations in figures 15–18.

Figure 17.

Transverse cuts at $\hat{x}=1.5$ for depth-integrated $\langle \overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }}^{{\mathcal{R}}}\rangle _{z}$ : 

 iLES; ○ NW model prediction. (a) $\overline{\unicode[STIX]{x1D70C}u_{1}^{\prime }}$ ; (b) $\overline{\unicode[STIX]{x1D70C}u_{2}^{\prime }}$ ; (c) $\overline{\unicode[STIX]{x1D70C}u_{3}^{\prime }}$ . Every eighth point shown for clarity. Correlation coefficients and normalized standard deviations in table 3.

Figure 17 shows a comparison between the depth-integrated MODEL prediction of TMF and the actual iLES data TMF at another NW transverse location. MODEL predicts TMF magnitude and transverse shape at this location extremely well. It also captures the peaks of the TMF for the streamwise and vertical components. However, it consistently under predicts the streamwise component near the centreline and the general magnitude of the transverse component. Figure 18 shows the same comparison between the depth-integrated MODEL prediction and iLES data of the FW region. MODEL predicts the average behaviour of the TMF and the peaks at the outer edges of the wake $|y/B|\sim 1.5$ at this $\hat{x}$ . However, it fails to capture the strong peaks of transverse and vertical TMF associated with the scars caused by the breaking waves as observed in figure 16. The correlation coefficients and standard deviations for figures 17 and 18 in table 3 show an overall good performance of the explicit model at these locations.

Figure 18.

Transverse cuts at $\hat{x}=3.0$ for depth-integrated $\langle \overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }}^{{\mathcal{R}}}\rangle _{z}$ : 

 iLES; ○ MODEL prediction. (a) $\overline{\unicode[STIX]{x1D70C}u_{1}^{\prime }}$ ; (b) $\overline{\unicode[STIX]{x1D70C}u_{2}^{\prime }}$ ; (c) $\overline{\unicode[STIX]{x1D70C}u_{3}^{\prime }}$ . Every eighth point shown for clarity. Correlation coefficients and normalized standard deviations in table 3.

Figure 19 shows the model predictions along the wake in aggregate with the conditioned correlation coefficient ${\mathcal{C}}$ . Overall, MODEL predicts the TMF qualitatively and quantitatively well with a few exceptions. The first is the under-prediction of transverse and vertical TMF in the CCW region. This is due to the presence of the fully overturning waves in case C, showing that GTH term in (5.2) fit using datasets without overturning is not sufficient. The second is the under-prediction of the streamwise TMF in the rooster-tail region $2\lesssim \hat{x}\lesssim 3$ . This is also due to the poor prediction near the centreline in figure 15(a). In a few regions with poor performance (namely $\overline{\unicode[STIX]{x1D70C}u_{1}^{\prime }}$ and $\overline{\unicode[STIX]{x1D70C}u_{2}^{\prime }}$ for $\hat{x}>\hat{x}_{B}$ and $\overline{\unicode[STIX]{x1D70C}u_{3}^{\prime }}$ for $\hat{x}<\hat{x}_{B}$ ) the conditioned correlation coefficient is high, implying that improving model coefficient in this region will improve the overall performance.

Figure 19.

Comparison of MODEL predictions of area-integrated $\langle \overline{\unicode[STIX]{x1D70C}u_{i}^{\prime }}^{{\mathcal{R}}}\rangle _{yz}$ along the wake. Lower plot: —— iLES; ○ NW Model; $\times$ FW Model; and upper plot: —— conditioned correlation coefficient ${\mathcal{C}}$ . (a) $\overline{\unicode[STIX]{x1D70C}u_{1}^{\prime }}$ ; (b) $\overline{\unicode[STIX]{x1D70C}u_{2}^{\prime }}$ ; (c) $\overline{\unicode[STIX]{x1D70C}u_{3}^{\prime }}$ . Every sixteenth point shown for clarity.

 represents $\hat{x}_{c}$ and $\hat{x}_{B}$ from left to right, respectively.

At this point, one may suggest further refinement of the model by dividing it into additional streamwise regions or using additional iLES datasets to improve the form and estimate of the model coefficients. The main objective of the present modelling process is to show whether or not it is feasible to explicitly predict the TMF in the wake from resolved/known quantities. Overall, our model shows that two main physical effects are directly responsible for much of the TMF in the wake: density gradients and buoyancy effects. With further refinement, explicit algebraic models should have a high potential for success in this type of scalar flux modelling.