2.1.1. Kinematic assumptions

The kinematics of the BMHV leaflets (figure 1b) consist of a rapid opening (${\rm \Delta} t\approx 20\text {--}50$ ms) at the onset of systole (Vennemann et al. Reference Vennemann, Rösgen, Heinisch and Obrist2018), a longer phase (${\rm \Delta} t\approx 350$ ms for a typical heart rate at rest) where the leaflets remain in the open position and a rapid closing (${\rm \Delta} t\approx 25$ ms) (Dasi et al. Reference Dasi, Ge, Simon, Sotiropoulos and Yoganathan2007; Borazjani et al. Reference Borazjani, Ge and Sotiropoulos2008; Yun et al. Reference Yun, Dasi, Aidun and Yoganathan2014a). The valve leaflets remain fully open at a low angle of incidence ($5^\circ$) during a large part of systole, which provides a suitable configuration for the creation of ILEV instabilities, as for larger angles of incidence the flow features will differ markedly from ILEV (Deniz & Staubli Reference Deniz and Staubli1997). Starting from the time instant that the leaflets are fully open, we focus on the systolic acceleration phase, which spans approximately 200 ms. Similar to Bellofiore et al. (Reference Bellofiore, Donohue and Quinlan2011) and following additional justifications presented in Zolfaghari & Obrist (Reference Zolfaghari and Obrist2019), we assume that the leaflets are fixed at an angle of $\theta =5^\circ$.

2.1.2. Three-dimensional aortic root and valve model

We define our 3-D model, as an extension of the 2-D model used in Zolfaghari & Obrist (Reference Zolfaghari and Obrist2019). We place the leaflets of a bileaflet mechanical heart valve in the fully open position (position 1 in black colour in figure 2a). The sinuses of Valsalva are included in the 3-D model as three spherical cavities with a radius of $r_s$ (figure 2b, bottom). The centres of these cavities are located on the $x=0$ plane and fall on the sides of an equilateral triangle, to which the aortic root's cross-section is circumscribed. Leaflets of the BMHV are modelled as blunt plates with a triangular leading- and trailing-edge geometry (figure 2a). The spanwise view of the leaflets is provided on the top panel of figure 2(b). The model geometry resembles the semi-elliptical trailing edge of a Regent BMHV. The leaflets’ thickness faces at leading and trailing edges follow the same geometry as this prosthesis. The main differences between the 3-D model geometry and a Regent BMHV are the absence of the housing ring and the valve ear/hinge recess in the 3-D model. However, these differences are located far from the centreline, hence they are not expected to notably affect the bulk flow downstream of the central orifice which is of interest in this paper.

Figure 2.

(a) A 3-D model of the BHMV in the aortic root. The 2-D cross-sections of leaflets in their plane of symmetry are shown in (1) open and (2) closed positions. (b, bottom) The cross-section of the root model at $x=0$ is shown. The three grey solid points indicate the centres of the spherical cavities representing the sinuses of the Valsalva. These cavities are positioned on the centres of the sides of an equilateral triangle, to which the aortic root cross-section is circumscribed ($r_s={r_r\sqrt {3}}/{2}$). (b, top) The geometry of the leaflet in the spanwise direction. The spanwise profile of the leading edge of the leaflet is modelled as a semicircle whose centre is located on the domain's axis of symmetry and on the leading edge (grey cross mark), and has a radius of $l_l$. All geometrical parameters of the model except those related to the cavities are given in table 1.

2.1.3. Boundary conditions and flow forcing

The flow is smoothly accelerated from zero to the mean velocity of $\mathscr {U}_{in}=2\mathscr {U}_0$ at $t=200\,{\rm ms}$ (figure 3). No-slip boundary conditions are imposed via an immersed boundary method on the rigid valve leaflets and aortic root boundaries. Periodic boundary conditions are specified in the streamwise direction $x$, where the systolic waveform is forced using a fringe region technique (Nordström, Nordin & Henningson Reference Nordström, Nordin and Henningson1999) upstream of the valve

(2.3)\begin{equation} \tilde{\boldsymbol{\mathsf{f}}}=\boldsymbol{\lambda(x)}(\tilde{\boldsymbol{\mathsf{u}}}- {\tilde{\boldsymbol{\mathsf{u}}_0}}(t)), \end{equation}

where $\lambda (x)$ is the fringe function and $\tilde {\boldsymbol {u}}_0(t)$ is a uniform flow profile. The amplitude of $\tilde {\boldsymbol {u}}_0(t)$ and the fringe function $\lambda (x)$ are tuned $\textit {ad~hoc}$ to provide the desired systolic flow acceleration (figure 3). The fringe region forces the appropriate inflow profile, and simultaneously damps out the outflow disturbances re-entering the domain at the domain's inlet due to periodic boundary conditions.

Table 1.

Geometrical parameters of the 3-D model.

Figure 3.

Acceleration part of the waveform used in the DNS. The horizontal dashed line marks the peak flow rate, corresponding to a mean flow velocity $\mathscr {U}_{in}=2\mathscr {U}_0$.

A highly resolved numerical simulation using $2561\times 513\times 257= 337\,644\,801$ degrees of freedom in a cuboid domain of size $3r_r\times 3r_r\times 15r_r$ is performed. Such a high spatial resolution is set to deal with the non-conforming geometry of the leaflets, and also to resolve the spatio-temporal instability waves and their interactions (refer to the effect of grid resolution in Zolfaghari & Obrist Reference Zolfaghari and Obrist2021). Grid stretching is applied in all three directions, so that a grid resolution similar to the one reported by Zolfaghari & Obrist (Reference Zolfaghari and Obrist2019) is achieved near the leaflets (31 grid points are placed along each leaflet's thickness length $\delta _l$). A dimensionless time-step size of $d{t}=10^{-4}$ was set, and a total of 29.1K time steps were integrated (up to the physical time of $t\approx 200\,{\rm ms}$). To the best of our knowledge, this is the highest resolution that has been used for a numerical simulation of BMHV flow using the incompressible Navier–Stokes equations. This simulation was completed using our novel hybrid multi-GPU task and data parallel Navier–Stokes solver on 20 GPUs in three days. Because our GPU-based implementation is approximately 150 times faster than the CPU-based MPI-parallel flow solver (using the same node configuration, i.e. 20 nodes of Cray XC40/50, Piz Daint), the present simulation ported onto an equivalent number of CPU cores ($20\times 16=320$ cores with hyper-threading) would have taken approximately 1.5 years to complete. The size of each data output of the velocity field amounted to $\approx$8.1 GB. Further details on this DNS and the numerical and parallelisation methodologies can be found in Zolfaghari & Obrist (Reference Zolfaghari and Obrist2021).

3.1.1. Governing equations

Decomposing the flow field at the leading edge ${\tilde {\boldsymbol{\mathsf{q}}}}=({\tilde {\boldsymbol{\mathsf{u}}}},\tilde {p})$ into a mean flow ${{\tilde {\boldsymbol{\mathsf{Q}}}}}=({\tilde {\boldsymbol{\mathsf{u}}}},\tilde {P})$ and a perturbation ${\tilde {\boldsymbol{\mathsf{Q}}}^\prime }$ (${\tilde {\boldsymbol{\mathsf{Q}}}}={\tilde {\boldsymbol{\mathsf{Q}}}} +{\tilde {\boldsymbol{\mathsf{Q}}}^\prime }$; $\|{\tilde {\boldsymbol{\mathsf{q}}}^\prime } \|\ll \|{\tilde {\boldsymbol{\mathsf{q}}}}\|$) and substituting into the Navier–Stokes equations leads to the linearised Navier–Stokes equations which read

(3.1a,b)\begin{equation} \partial_t{\tilde{\boldsymbol{\mathsf{u}}}^\prime}+{\tilde{\boldsymbol{\mathsf{u}}}^\prime}\boldsymbol{\cdot} \boldsymbol{\nabla}{\tilde{\boldsymbol{\mathsf{u}}}}+{\tilde{\boldsymbol{\mathsf{u}}}}\boldsymbol{\cdot} \boldsymbol{\nabla}{\tilde{\boldsymbol{\mathsf{u}}}^\prime}={-}\boldsymbol{\nabla}\tilde{p}^\prime+\frac{1}{Re_*}\nabla^2{\tilde{\boldsymbol{\mathsf{u}}}^\prime};\quad \boldsymbol{\nabla}\boldsymbol{\cdot}{\tilde{\boldsymbol{\mathsf{u}}}^\prime}=0. \end{equation}

Here, $\tilde {\boldsymbol{\mathsf{U}}}$ denotes the 2-D time-averaged flow field obtained as

(3.2)\begin{equation} {{\tilde{\boldsymbol{\mathsf{U}}}}}(\xi,\eta) =\frac{1}{\boldsymbol{{\tilde{T}}}} \int^{t_0+\tilde{\boldsymbol{T}}}_{t_0} \frac{\tilde{\boldsymbol{\mathsf{u}}}(\xi,\eta,t)}{U_{{max},\xi=0}(t)}{\rm d} t. \end{equation}

The Reynolds number $Re_*$ is based on the maximum streamwise velocity at the leading edge ($\overline {U_{max,\xi =0}}$) as the reference velocity and the distance between the leaflets at the leading edge ($\eta ^*_0=\eta ^*(\xi =0)$) as the reference length. We have

(3.3)\begin{equation} Re_* = \frac{{\overline{U_{max,\xi=0}}}\eta^*_0}{\nu}, \end{equation}

where

(3.4)\begin{equation} \overline{U_{max,\xi=0}} = \frac{1}{\tilde{\boldsymbol{{T}}}} \int^{t_0+\tilde{\boldsymbol{T}}}_{t_0}{U_{\text{max},\xi=0}(t)}\,{\rm d} t. \end{equation}

Introducing a global mode ansatz of the form $\tilde {\boldsymbol{\mathsf{q}}}^\prime (\boldsymbol{\mathsf{x}},t)=\hat{\boldsymbol{\mathsf{q}}} (\boldsymbol{\mathsf{x}})\,{\rm e}^{-{\rm i}\tilde {\beta }t}$, where $\boldsymbol{\mathsf{x}}$ and $t$ denote the 2-D spatial coordinates and time, into the LNSE (3.1a,b) yields a system of partial differential equations (PDEs)

(3.5)\begin{equation} \begin{bmatrix} \mathscr{L}\{\tilde{\boldsymbol{\mathsf{u}}};Re_*\}+\tilde{U}_\xi & \tilde{U}_\eta & \mathscr{D}_\xi \\ \tilde{V}_\xi & \mathscr{L}\{\tilde{\boldsymbol{\mathsf{u}}};Re_*\} & \mathscr{D}_\eta \\ \mathscr{D}_\xi & \mathscr{D}_\eta & 0 \end{bmatrix} \begin{bmatrix} \hat{u} \\ \hat{v} \\ \hat{p} \end{bmatrix}= \mathrm{i}\tilde{\beta} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} \hat{u} \\ \hat{v} \\ \hat{p} \end{bmatrix} \end{equation}

where $\mathscr {L}\{\boldsymbol {\tilde {U}};Re_*\}=\tilde {U}\mathscr {D}_\xi +\tilde {V}\mathscr {D}_\eta -1/Re_*(\mathscr {D}_{\xi \xi }+\mathscr {D}_{\eta \eta })$, and $\skew2\hat{{\boldsymbol{\mathsf{q}}}}=(\hat{{\boldsymbol{\mathsf{u}}}},\hat {{p}})$ with the velocity vector ${{\hat{{\boldsymbol{\mathsf{u}}}}}}=({\hat {u}},\hat {{v}})$. In the above expression $\mathscr {D}_\xi$ and $\mathscr {D}_{\xi \xi }$ denote the first and second derivatives with respect to the $\xi$-direction, respectively.

Temporal global modes with $\operatorname {Im}\{\tilde {\beta }\}>0$ are sought. To this end, we discretise (3.5) using second-order finite-differences with $n_p=123\,71$ grid points. This results in a complex generalised eigenvalue problem of the form

(3.6)\begin{equation} \boldsymbol{\mathsf{A}}_{3n_p\times 3n_p}(Re_*,{\tilde{\boldsymbol{\mathsf{u}}}}) \boldsymbol{\hat{\boldsymbol{\mathsf{q}}}}_{3n_p\times 1}=\mathrm{i}\tilde{\beta} \boldsymbol{\mathsf{B}}_{3n_p\times 3n_p}\boldsymbol{\hat{\boldsymbol{\mathsf{q}}}}_{3n_p\times 1}, \end{equation}

where $\boldsymbol{\hat{\boldsymbol{\mathsf{q}}}}$ is the discretised form of $\hat{{\boldsymbol{\mathsf{q}}}}$.

3.1.2. Integration of the immersed boundaries

The sharp-interface immersed boundary method (Mittal et al. Reference Mittal, Dong, Bozkurttas, Najjar, Vargas and von Loebbecke2008) is used to integrate the boundary conditions associated with the leaflet walls into the matrices $\boldsymbol{\mathsf{A}}$ and $\boldsymbol{\mathsf{B}}$ (see (3.6)). The workflow of labelling the grid points according to their position with respect to the boundaries is demonstrated in algorithm 1. Compatibility boundary conditions for pressure are imposed on the domain's inlet and outlet as well as on the solid walls. These conditions, in general, read

(3.7)\begin{equation} \boldsymbol{\nabla}\tilde{p}^\prime={-}{\tilde{\boldsymbol{\mathsf{u}}}} \boldsymbol{\cdot}\boldsymbol{\nabla}{\tilde{\boldsymbol{\mathsf{u}}}^\prime}+\frac{1}{Re_*}\nabla^2{\tilde{\boldsymbol{\mathsf{u}}}^\prime}. \end{equation}

On the leaflet walls, this condition simplifies since $\tilde {\boldsymbol{\mathsf{U}}}=\boldsymbol{\mathsf{0}}$, and (3.7) reduces to

(3.8)\begin{equation} \boldsymbol{\nabla}\tilde{p}^\prime=\frac{1}{Re_*}\nabla^2{{\tilde{\boldsymbol{\mathsf{u}}}^\prime}}. \end{equation}

The compatibility boundary conditions, together with zero-disturbance Dirichlet boundary conditions for the velocity vector given by

(3.9)\begin{equation} \tilde{\boldsymbol{\mathsf{u}}}^\prime(\boldsymbol{\mathsf{x}},t)=0,\quad \boldsymbol{\mathsf{x}}\in{\partial_{in}\cup \partial_{walls}\cup\partial_{out} \cup\partial_\infty}, \end{equation}

are discretised and integrated into the matrices $\boldsymbol{\mathsf{A}}$ and $\boldsymbol{\mathsf{B}}$. To save computational time, and owing to the symmetry of the velocity profiles, only half of the profile is considered, which is results in a symmetry boundary condition at $\tilde {\eta } = 0$.

Algorithm 1

Two-dimensional global instability solver with immersed boundaries

Figure 9 shows the real part of the streamwise velocity component of the two zero-frequency global eigenmodes with positive temporal growth rates. Mode A, which has a higher growth rate, has a similar spatial structure than the one found for X-junction flow (Lashgari et al. Reference Lashgari, Tammisola, Citro, Juniper and Brandt2014), stretching from the separation point to the trailing end of the wave-maker zone. This could be related to similarities in the two base flows: both include flow separation downstream of a sharp edge. Yet, our mode presents more waviness past the wave maker, which includes the area where the separated shear layer reattaches. This waviness may be attributed to the strong fluctuations that are observed over the reattachment zone of ILEV structures, which were also reported by Cherry, Hillier & Latour (Reference Cherry, Hillier and Latour1984). Mode B, on the other hand, is more concentrated near the leading edge of the wave-maker zone. It slowly emerges on the separated shear layer nearly two thickness lengths downstream of the leading edge, and vanishes towards the reattachment zone. Provided this mode is present where the shear layer has the highest reverse-to-forward flow ratio, it projects the dynamics of the shear layer seemingly uninfluenced by the presence of the wall. This is in contrast to mode A, where oscillatory effects can be observed at locations where this latter ratio is low. In line with this, mode A includes a wavy structure that is extended to the outflow boundary, before it is affected by the homogeneous Dirichlet boundary conditions there. We will show in Appendix A that the structure of this mode and mode B (as well as their associated eigenvalues) do not change when the outflow boundary is extended further downstream.

Figure 9.

Global modes (a) A and (b) B, visualised by the real part of the streamwise velocity component. Dashed lines mark the extent of the local absolute instability. Approximate separation streamlines for the mean flow are drawn as solid black lines. Dark red and dark blue colours indicate the colour map extrema.

3.2.1. Adjoint modal formulation

For investigating the sensitivity of the global modes characteristics (e.g. their growth rates) to the underlying flow, adjoints of the global modes are studied. For the pairs $(\tilde {\boldsymbol{\mathsf{u}}}^\prime,\tilde {p}^\prime )$ and $(\tilde {\boldsymbol{\mathsf{u}}}^{\prime,{\dagger} },\tilde {p}^{\prime,{\dagger} })$ the following Lagrange identity holds as:

(3.10)\begin{equation} {[\mathscr{P}\boldsymbol{\cdot} \tilde{\boldsymbol{\mathsf{u}}}^{\prime,{\dagger}} + \boldsymbol{\nabla}\boldsymbol{\cdot} \tilde{\boldsymbol{\mathsf{u}}}^\prime \tilde{p}^{\prime,{\dagger}}]} + {[\tilde{\boldsymbol{\mathsf{u}}}^\prime\boldsymbol{\cdot}\mathscr{P}^{\dagger}{+} \boldsymbol{\nabla} \boldsymbol{\cdot}\tilde{\boldsymbol{\mathsf{u}}}^{\prime,{\dagger}} {\tilde{p}}^\prime]}= \partial_t (\tilde{\boldsymbol{\mathsf{u}}}^\prime\boldsymbol{\cdot}\tilde{\boldsymbol{\mathsf{u}}}^{\prime,{\dagger}}) + \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{M}}^{\dagger}, \end{equation}

where

(3.11)\begin{equation} \mathscr{P} \{\tilde{\boldsymbol{\mathsf{u}}}^\prime,{\tilde{p}^{\prime}}, \tilde{\boldsymbol{\mathsf{u}}}, Re_*\} = \partial_t \tilde{\boldsymbol{\mathsf{u}}}^\prime+ \tilde{\boldsymbol{\mathsf{u}}} \boldsymbol{\cdot}\boldsymbol{\nabla}\tilde{\boldsymbol{\mathsf{u}}}^\prime + \tilde{\boldsymbol{\mathsf{u}}}^\prime\boldsymbol{\cdot}\boldsymbol{\nabla} \tilde{\boldsymbol{\mathsf{u}}}+ \boldsymbol{\nabla} {\tilde{p}}^\prime-\frac{1}{Re_*}\nabla^2 \tilde{\boldsymbol{\mathsf{u}}}^{\prime,{\dagger}},\end{equation}

and

(3.12)\begin{equation} \mathscr{P}^{\dagger}\{ \tilde{\boldsymbol{\mathsf{u}}}^{\prime,{\dagger}},{\tilde{p}^{\prime,{\dagger}}}, \tilde{\boldsymbol{\mathsf{u}}}, Re_*\}= \partial_t \tilde{\boldsymbol{\mathsf{u}}}^{\prime,{\dagger}}+ \tilde{\boldsymbol{\mathsf{u}}}\boldsymbol{\cdot}\boldsymbol{\nabla} \tilde{\boldsymbol{\mathsf{u}}}^{\prime,{\dagger}} - \tilde{\boldsymbol{\mathsf{u}}}^{\prime,{\dagger}} \boldsymbol{\cdot} (\boldsymbol{\nabla} \tilde{\boldsymbol{\mathsf{u}}})^{\rm T}+ \boldsymbol{\nabla} \tilde{p}^{\prime,{\dagger}} + \frac{1}{Re_*}\nabla^2 \tilde{\boldsymbol{\mathsf{u}}}^{\prime,{\dagger}}. \end{equation}

Here, $\boldsymbol{\mathsf{M}}^{\dagger}$ is the bilinear concomitant defined as

(3.13)\begin{equation} \boldsymbol{\mathsf{M}}^{\dagger}{=}\tilde{\boldsymbol{\mathsf{u}}} (\tilde{\boldsymbol{\mathsf{u}}}^\prime\boldsymbol{\cdot}\tilde{\boldsymbol{\mathsf{u}}}^{\prime,{\dagger}})+ \frac{1}{Re_*}(\boldsymbol{\nabla} \tilde{\boldsymbol{\mathsf{u}}}^{\prime,{\dagger}}\boldsymbol{\cdot} \tilde{\boldsymbol{\mathsf{u}}}^\prime- \boldsymbol{\nabla} \tilde{\boldsymbol{\mathsf{u}}}^\prime\boldsymbol{\cdot} \tilde{\boldsymbol{\mathsf{u}}}^{\prime,{\dagger}})+\tilde{p}^{\prime,{\dagger}} \tilde{\boldsymbol{\mathsf{u}}}^\prime+{\tilde{p}}^\prime \tilde{\boldsymbol{\mathsf{u}}}^{\prime,{\dagger}}. \end{equation}

Integrating the Lagrange identity over the entire domain and over the chosen time horizon and using the divergence theorem for the last term, we obtain the adjoint LNSE in the form

(3.14a,b)\begin{equation} \mathscr{P}^{\dagger}\{ \tilde{\boldsymbol{\mathsf{u}}}^{\prime,{\dagger}},{\tilde{p}^{\prime,{\dagger}}}, \tilde{\boldsymbol{\mathsf{u}}}, Re_*\} =0;\quad \boldsymbol{\nabla}\boldsymbol{\cdot} \tilde{\boldsymbol{\mathsf{U}}}^{\prime,{\dagger}}=0. \end{equation}

Fourier transformation of (3.10) using the ansatz ${\tilde {\boldsymbol{\mathsf{q}}}}^{\prime,{\dagger} } (\boldsymbol{\mathsf{x}},t)= \hat{{\boldsymbol{\mathsf{q}}}}^{\dagger} (\boldsymbol{\mathsf{x}}) \,{\rm e}^{-\mathrm {i}\tilde {\beta }^{\dagger} t}$, where $\boldsymbol{\mathsf{x}}$ and $t$ are the 2-D spatial coordinates and time, gives the coupled form of the adjoint LNSE

(3.15)\begin{align} \begin{bmatrix} \mathscr{L}^{\dagger}{\{\tilde{\boldsymbol{\mathsf{u}}};Re_*\}} - \partial_\xi\tilde{U} & - \partial_\xi\tilde{V} & \mathscr{D}_\xi \\ -\partial_\eta\tilde{U} & \mathscr{L}^{\dagger}{\{\tilde{\boldsymbol{\mathsf{u}}}; Re_*\}}-\partial_\eta\tilde{V} & \mathscr{D}_\eta \\ \mathscr{D}_\xi & \mathscr{D}_\eta & 0 \end{bmatrix} \begin{bmatrix} \hat{u}^{\dagger} \\ \hat{v}^{\dagger} \\ \hat{p}^{\dagger} \end{bmatrix} =\mathrm{i}\tilde{\beta}^{\dagger} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0\\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} \hat{u} ^{\dagger} \\ \hat{v} ^{\dagger} \\ \hat{p} ^{\dagger} \end{bmatrix}, \end{align}

where $\mathscr {L}^{\dagger} {\{\tilde {\boldsymbol{\mathsf{u}}};Re_*\}}= \tilde {U}\mathscr {D}_\xi +\tilde {V}\mathscr {D}_\eta +Re_*^{-1}(\mathscr {D}_{\xi \xi }+\mathscr {D}_{\eta \eta })$, ${\hat{ {\boldsymbol{\mathsf{q}}}}^{{\dagger} }}= (\hat{{\boldsymbol{\mathsf{u}}}}^{\dagger},\hat {{p}}^{{\dagger} })$ and ${\hat{{\boldsymbol{\mathsf{u}}}}^{{\dagger} }}=({\hat {u}}^{\dagger},\hat {{v}}^{{\dagger} })$. The operators $\mathscr {D}_{\xi }$ and $\mathscr {D}_{\xi \xi }$ denote the first and second derivatives with respect to $\xi$, respectively. Discretisation yields the generalised eigenvalue problem

(3.16)\begin{equation} \boldsymbol{\mathsf{A}}^{\dagger}_{3n_p\times 3n_p}({\tilde{\boldsymbol{\mathsf{U}}}},Re_*) \boldsymbol{\hat{\boldsymbol{\mathsf{q}}}}^{\dagger}_{3n_p\times 1}=\mathrm{i}{\tilde{\beta}}^{\dagger} \boldsymbol{\mathsf{B}}^{\dagger}_{3n_p\times 3n_p}({\tilde{\boldsymbol{\mathsf{U}}}},Re_*) \boldsymbol{\hat{\boldsymbol{\mathsf{q}}}}^{\dagger}_{3n_p\times 1}, \end{equation}

where $\boldsymbol{\hat{\boldsymbol{\mathsf{q}}}}^{\dagger}$ is the discretised form of adjoint state vector. Homogeneous Dirichlet boundary conditions for the adjoint velocity disturbances $\tilde {\boldsymbol{\mathsf{u}}}^{\prime,{\dagger} }$ are applied, together with compatibility boundary condition for the adjoint pressure disturbance $\tilde {p}^{\prime,{\dagger} }$. The compatibility boundary conditions for $\tilde {p}^{\prime,{\dagger} }$ prove to be the proper choice for the domain's inflow and outflow. To obtain the adjoint compatibility boundary conditions, we first substitute this condition for the direct problem (3.7) into the momentum component of the direct linearised Navier–Stokes equations (3.1a,b). This yields

(3.17)\begin{equation} \partial_t{{\tilde{\boldsymbol{\mathsf{u}}}^\prime}}+{{\tilde{\boldsymbol{\mathsf{u}}}^\prime}}\boldsymbol{\cdot} \boldsymbol{\nabla}{{\tilde{\boldsymbol{\mathsf{U}}}}}=0, \end{equation}

which, in the adjoint form, becomes

(3.18)\begin{equation} \partial_t{{\tilde{\boldsymbol{\mathsf{u}}}^{\prime,{\dagger}}}}+{{\tilde{\boldsymbol{\mathsf{U}}}}} \boldsymbol{\cdot}\boldsymbol{\nabla}\tilde{\boldsymbol{\mathsf{u}}}^{\prime,{\dagger}}=0. \end{equation}

Equation (3.18) is incorporated into $\boldsymbol{\mathsf{A}}^{\dagger}$ and $\boldsymbol{\mathsf{B}}^{\dagger}$ at the inflow, centreline, and outflow locations, where $\tilde {\boldsymbol{\mathsf{U}}}\neq \boldsymbol{\mathsf{0}}$.

Figure 10 shows the adjoint global modes corresponding to the direct global modes A and B. The adjoint modes can be used to quantify the sensitivity of the global direct modes to external forcing (Giannetti & Luchini Reference Giannetti and Luchini2007). We use both direct and adjoint modes for constructing the structural sensitivity tensor which can be used to guide the stabilisation of these modes.

Figure 10.

Adjoint global modes (a) A and (b) B, visualised by the streamwise velocity component. Dashed lines mark the extent of the local absolute instability. Approximate separation streamlines for the mean flow are drawn as solid black lines. Dark red and dark blue colours indicate the colour map extrema. Larger magnitudes correspond to higher sensitivities of the global modes to the underlying flow.

3.2.2. Structural sensitivity to a local feedback source

We follow the localised momentum feedback approach pioneered by Hill (Reference Hill1992) and Giannetti & Luchini (Reference Giannetti and Luchini2007). We perturb the momentum equation of the LNSE, after Fourier decomposition in time, by the linear operator $\delta \mathscr {M}$

(3.19)\begin{equation} \left.\begin{gathered} \tilde{\beta}^\prime \hat{{\boldsymbol{\mathsf{u}}}}^{\prime} + \mathscr{L}^+\{\tilde{\boldsymbol{\mathsf{U}}},Re_*\}\hat{{\boldsymbol{\mathsf{u}}}}^{\prime}+ \boldsymbol{\nabla}\hat{p}^{\prime}=\boldsymbol{\delta}\mathscr{M}(\hat{{\boldsymbol{\mathsf{u}}}}^{\prime},\hat{p}^{\prime}) \\ \boldsymbol{\nabla}\boldsymbol{\cdot}\hat{{\boldsymbol{\mathsf{u}}}}^{\prime}=0, \end{gathered}\right\}, \end{equation}

where the $\mathscr {L}^+$ operator is given as

(3.20)\begin{equation} \mathscr{L}^+\{{\tilde{\boldsymbol{\mathsf{U}}}}, Re_*\}\hat{{\boldsymbol{\mathsf{u}}}}^{\prime}= \hat{{\boldsymbol{\mathsf{u}}}}^{\prime}\boldsymbol{\cdot}\boldsymbol{\nabla} \tilde{\boldsymbol{\mathsf{U}}}+\tilde{\boldsymbol{\mathsf{U}}} \boldsymbol{\cdot}\boldsymbol{\nabla}\hat{{\boldsymbol{\mathsf{u}}}}^{\prime} -\frac{1}{Re_*}\nabla^2 \hat{{\boldsymbol{\mathsf{u}}}}^{\prime}. \end{equation}

Given the eigenvalue drift $\delta \tilde {{\beta }}$ and the perturbation eigenmode $(\delta \hat{{\boldsymbol{\mathsf{u}}}},\delta \hat {p})$, where $\hat{ {\boldsymbol{\mathsf{u}}}}^{\prime }=\hat{ {\boldsymbol{\mathsf{u}}}}+\delta \hat{{\boldsymbol{\mathsf{u}}}}$, $\hat {{p}}^{\prime }=\hat {{p}}+\delta \hat {{p}}$ and ${\tilde {\beta }}^\prime ={\tilde {\beta }}+\delta {\tilde {\beta }}$, we arrive at

(3.21)\begin{equation} \left.\begin{gathered} \tilde{\beta} \delta\hat{{\boldsymbol{\mathsf{u}}}} + \mathscr{L}^+\{\tilde{\boldsymbol{\mathsf{U}}},Re_*\}\delta\hat{{\boldsymbol{\mathsf{u}}}}+ \boldsymbol{\nabla}\delta\hat{p}={-}\delta{\tilde{\beta}} \hat{{\boldsymbol{\mathsf{u}}}} +{{{\delta}}}\mathscr{M}(\hat{{\boldsymbol{\mathsf{u}}}},\hat{p})\\ \boldsymbol{\nabla}\boldsymbol{\cdot}\delta\hat{{\boldsymbol{\mathsf{u}}}}=0, \end{gathered}\right\} \end{equation}

using the Lagrange identity (similar to (3.10)) for $(\delta \hat{{\boldsymbol{\mathsf{u}}}},\delta \hat {p})$ and $(\hat{ {\boldsymbol{\mathsf{u}}}}^{{\dagger} },\hat {p}^{{\dagger} })$. Taking $\delta \mathscr {M} = \boldsymbol{\mathsf{K}}\hat{{\boldsymbol{\mathsf{u}}}} \delta _{\boldsymbol{\mathsf{D}}}(\xi -\xi _+,\eta -\eta _+)$ and integrating over the entire computational domain, one obtains

(3.22)\begin{equation} \delta\tilde{\beta}=\frac{\langle \hat{{\boldsymbol{\mathsf{u}}}}^{{{\dagger}}} \boldsymbol{\cdot} \delta\mathscr{M}\rangle}{\langle\hat{{\boldsymbol{\mathsf{u}}}}\boldsymbol{\cdot} \hat{{\boldsymbol{\mathsf{u}}}}^{{{\dagger}}} \rangle}= \frac{\langle\hat{{\boldsymbol{\mathsf{u}}}}^{{{\dagger}}}\boldsymbol{\cdot}\boldsymbol{\mathsf{K}} \hat{{\boldsymbol{\mathsf{u}}}}\delta_{\boldsymbol{\mathsf{D}}}(\xi-\xi_+,\eta-\eta_+)\rangle} {\langle\hat{{\boldsymbol{\mathsf{u}}}}\boldsymbol{\cdot}\hat{{\boldsymbol{\mathsf{u}}}}^{{{\dagger}}} \rangle} = \frac{\kappa_0\hat{{\boldsymbol{\mathsf{u}}}}(\xi_+,\eta_+)\hat{{\boldsymbol{\mathsf{u}}}}^{{{\dagger}}} (\xi_+,\eta_+)}{\langle\hat{{\boldsymbol{\mathsf{u}}}}\boldsymbol{\cdot}\hat{{\boldsymbol{\mathsf{u}}}}^{{{\dagger}}} \rangle}. \end{equation}

In the above formulation, $\delta _{\boldsymbol{\mathsf{D}}}$ is the Dirac delta function, $\langle \cdot \rangle$ is the integral over the computational domain and $\boldsymbol{\mathsf{K}}$ and $\kappa _0=\boldsymbol{\mathsf{K}}(\xi _+,\eta _+)$ indicate the strength of the momentum feedback on the entire domain and at point $(\xi _+, \eta _+)$, respectively. It follows that the response of each component of momentum feedback on the growth rate and frequency of the global modes can be realised by assessing the real and imaginary parts of the sensitivity tensor $\boldsymbol{\mathsf{S}}=\hat{{\boldsymbol{\mathsf{u}}}} \otimes \hat{{\boldsymbol{\mathsf{u}}}}^{\dagger}$, which consists of the dyadic product of the direct and adjoint velocity modes. The imaginary part of the components of $\boldsymbol{\mathsf{S}}$ measure the sensitivity of the global mode with respect to the growth rate, while the real part determines the sensitivity with respect to the frequency. Here, we focus on the imaginary part, as we are interested in stabilisation of the modes, hoping to eliminate the vortex shedding between and downstream of the leaflets of the BMHV.

The components of the sensitivity tensor associated with modes A and B are given in figures 11 and 12. For mode A, the highest sensitivities are generally concentrated in the wave-maker zone. This is not surprising as a local absolute instability is defined by an instability with infinite impulse response at a fixed location. A feedback, if sufficiently small to preserve the instability characteristics of the base flow, could be seen in extreme cases as an input, which will be maximally amplified in the wave-maker zone. In terms of different components of sensitivity, $S_{12}=\hat {u}\hat {v}^{\dagger}$ appears to have a slightly higher magnitude than other components. This implies that, within a linear regime, a forcing in the cross-stream direction will trigger the largest response in the streamwise component of the direct mode. This could be explained by further investigating the characteristics of the base flow, such as the properties of the rate of strain tensor, but such an analysis will be postponed to a future effort. For mode B, the maximum structural sensitivity occurs upstream and almost completely outside the wave-maker zone. This observation signals a large influence of the base flow on the dynamics of this shear-layer mode, possibly through a feedback mechanism supported by recirculation.

Figure 11.

Components of the sensitivity tensor $\boldsymbol{\mathsf{S}}$ for mode A. Yellow areas denote the maximum, blue areas denote the minimum sensitivity. The domain between the dashed white lines is locally absolutely unstable.

Figure 12.

Components of the sensitivity tensor $\boldsymbol{\mathsf{S}}$ for mode B. Yellow areas denote the maximum, blue areas denote the minimum sensitivity. The domain between the dashed white lines is locally absolutely unstable (wave-maker region).

3.2.3. Passive control based on localised feedback

Motivated by the original experimental work of Strykowski & Sreenivasan (Reference Strykowski and Sreenivasan1990) and the follow-up theoretical study of Hill (Reference Hill1992) who tried to suppress the first instability of cylinder wake flow at $Re\approx 50$ by introducing localised roughness elements into the system, we use this technique to try to suppress the ILEV instabilities in the BMHVs. Several differences exist between the two cases, the most importantly, the two orders of magnitude larger Reynolds number of the BMHV flow compared with the cylinder wake flow. It is also not clear whether a localised feedback will strengthen, or weaken the instability. Besides, given that more than one global mode is involved, with high structural sensitivities in distinct locations, the likelihood that only one feedback source may eliminate the instability entirely is likely small. It is also important to note that, because the Reynolds number is relatively high in our case, a localised feedback, positioned relatively far from the wall may cause more instability, e.g. via introducing von-Kármán-type vortex streets. Despite these apprehensions, we nonetheless put our findings from the linear structural sensitivity analysis to the test by performing 2-D DNSs. We introduce small feedback sources in the form of small cylinders (far from the wall, for mode A), or in the form of a semicircular bump (close to the wall, for mode B), and perform 2-D DNS using the parameters from Zolfaghari & Obrist (Reference Zolfaghari and Obrist2019). For mode A, we observed that a localised feedback did not seem to improve the flow scenario (not shown). The inserted cylinder near the upstream end of the wave-maker zone created instabilities which evolved into travelling vortices downstream of the feedback source. Given that mode B was not triggered, it induced waviness in the ILEV which interacted with the instabilities created by the cylinder. For mode B, we anticipated a more positive outcome, as this mode has a more localised sensitivity compared with mode A. This localised area is closer to the wall which, due to lower velocities and a smaller surface, generates less drag force on the original base flow. Thus, the base flow characteristics are likely to change less than for mode A. Significant changes in the base flow $\tilde {\boldsymbol{\mathsf{U}}}$ may render the structural sensitivity analysis invalid, because the eigenvalue shift $\delta \beta$ based on (3.22) assumes small and thus negligible changes in the base flow. A better outcome is also expected because the area to be triggered lies outside the absolute instability zone. For this reason, the chances of triggering more instabilities by the large impulse response of the wave maker are rather low.

Figure 13 shows the outcome of passive control via local feedback associated with mode B. For the local feedback, two small and equal circular bumps with a radius of $R_c=0.126$ mm were added on both leaflet surfaces at $\xi _c=4.95$ and $\eta _c=\pm 1.68$ mm. The positions of these bumps were chosen to be approximately in the centre of the maximum sensitivity area for mode B. It can be observed that with passive control the vorticity production between the leaflets was diminished, but not eliminated. In more detail, as shown using instantaneous vorticity fields, the location where first vortices, forming due to ILEV instabilities, detach into the bulk flow could be moved downstream by nearly one third of a chord length. Thanks to this improvement, the mutual interaction between the vortices, generated by ILEVs and forming on the top and bottom leaflets, was nearly eliminated. In general, even though the ILEV instability could not be fully controlled, this passive control device resulted in a noticeably less chaotic wake flow past the valve. This is further quantified in figure 14 using area-normalised enstrophy defined as

(3.23)\begin{equation} \boldsymbol{E}_S(\pmb{\omega})= \frac{1}{S}\int_S\pmb{\omega}\boldsymbol{\cdot}\pmb{\omega}\,{\rm d} S, \end{equation}

where $\pmb {\omega }=\boldsymbol {\nabla }\times \boldsymbol{\mathsf{u}}$ is the vorticity and $S$ indicates the circular area where the integration is performed. Significantly lower enstrophy is found for the control case both within the central orifice and in the wake, which demonstrates the efficacy of the control.

Figure 13.

Instantaneous vorticity for (a) baseline ILEV flow and (b) ILEV flow with local feedback control corresponding to mode B.

Figure 14.

Effect of passive control based on momentum feedback on normalised enstrophy $\boldsymbol {E}_S(\pmb {\omega })$ at ((a), in central orifice) 5.76 mm and ((b), in the wake) 12.96 mm downstream of the leading edge. The integration area $S$ was a circle of radius 3.6 mm and centred on $r=0$ mm.

4.1.1. Governing equations

We again consider the incompressible Navier–Stokes equations for a disturbance $(\tilde {\boldsymbol{\mathsf{u}}}^\prime,{\tilde {p}}^\prime )$ around a global mean flow $(\tilde {\boldsymbol{\mathsf{U}}},\tilde {\boldsymbol{\mathsf{P}}})$,

(4.1a,b)\begin{equation} \partial_t{{\tilde{\boldsymbol{\mathsf{u}}}}^\prime}+ {{\tilde{\boldsymbol{\mathsf{u}}}}^\prime}\boldsymbol{\cdot} \boldsymbol{\nabla}{{\tilde{\boldsymbol{\mathsf{U}}}}}+{{\tilde{\boldsymbol{\mathsf{U}}}}} \boldsymbol{\cdot}\boldsymbol{\nabla}{\tilde{\boldsymbol{\mathsf{u}}}^\prime}={-}\boldsymbol{\nabla}\tilde{p}^\prime+ \frac{1}{Re}\nabla^2{{\tilde{\boldsymbol{\mathsf{u}}}}^\prime};\quad \boldsymbol{\nabla}\boldsymbol{\cdot}{{\tilde{\boldsymbol{\mathsf{u}}}}^\prime}=0, \end{equation}

and seek to recover the ILEV instabilities as disturbances superimposed on a mean flow $\tilde {\boldsymbol{\mathsf{U}}}$ at each time $t$. Note that we use a different mean flow for the nonlinear non-modal analysis than for the earlier linear modal approach give in (3.1a,b). Here, $\tilde {\boldsymbol{\mathsf{U}}}$ is obtained by temporally averaging of the velocity field $\tilde {\boldsymbol{\mathsf{u}}};$ no scaling is involved, and the averaging is performed over the entire flow domain. We have

(4.2)\begin{equation} {{{{\tilde{\boldsymbol{\mathsf{U}}}}}}}(x,r) = \frac{1}{\tilde{\boldsymbol{{T}}}}\int^{t_0+\tilde{\boldsymbol{T}}}_{t_0} {\tilde{\boldsymbol{\mathsf{u}}}(x,r,t)}\,{\rm d} t, \end{equation}

which also implies that the statistical averaging is performed over the global $(x,r)$ instead of the local $(\xi,\eta )$ coordinates. The time-averaged flow field around the ILEV is shown in figure 15.

Figure 15.

Time-averaged (a) streamwise and (b) cross-stream velocity components around the ILEV in the BMHV.

4.3.1. Maximum energy growth between the leading and trailing edges

We start by analysing the energy growth between the leading and trailing edges (referred to as mid-leaflet in what follows). This area $\varOmega _v$ is located downstream of the ILEV zone and overlaps with the wave maker. It is expected that the optimal initial condition will be mainly focused around the leading edge, particularly in the ILEV zone, due to its close proximity and small value of $T$.

A circular energy probe of radius $R=3.6$ mm is located at $(\xi,r)=(5.76,0\ \textrm {mm})$. DAL simulations are then performed for different times $T=0.01$, $T=0.05$ and $T=0.1$. In contrast to the DAL simulations for channel flow, (see Appendix C) the current simulations require larger numbers of iterations to converge (see figure 16). This may be a consequence of the presence of multiple globally unstable modes in the system, e.g. stemming from the ILEV, cavities and the wake.

Figure 16.

Growth (left axis) and residual (right axis) for the mid-leaflet mask and for $T=0.01$.

Figure 17(a,b) shows the initial ($t=0$) and terminal ($t=T$) states corresponding to maximum energy growth at the mid-leaflet mask and $T=0.01$ (smallest computed time). Due to the short terminal time, the optimal initial conditions are mainly focused within the mask. This suggests that instability growth for this mask remains almost entirely locally for small time scales. An exploration of instabilities outside the mask requires a larger time horizon $T.$ Panels (c,d) and (e, f) of figure 17 show results for $T=0.05$ and $T=0.1,$ respectively. As expected, the optimal initial condition moves towards the leading edge, including a part of the ILEV zone, as $T$ is increased. These results furthermore reveal an additional ‘actor’ besides the ILEV: the flow impingement zone located on the leaflet's thickness (area upstream of the green dashed line in the bottom row of figure 17). This area displays high potential for creating disturbances that are subsequently amplified by the ILEV. Figure 17(c–f) also shows that the initial and terminal states are slightly asymmetric despite the symmetry of the leaflets. This is due to minor differences in the depth of the cavities in the 2-D submodel (Zolfaghari & Obrist Reference Zolfaghari and Obrist2019), which results in slightly asymmetric ILEV profiles on the upper and lower leaflets.

Figure 17.

Streamwise component of (a,c,e) optimal initial and (b,d,f) final states are demonstrated for a mid-leaflet mask. From top to bottom, cases with $T=0.01$, $T=0.05$ and $T=0.1$ are shown, respectively. The domain within the circle identifies the mask $\varOmega _v$. Dark blue and light yellow areas signify minimum negative and maximum positive velocity disturbances, respectively.

4.3.2. Maximum energy growth at the trailing edge

Next, we displace the mask $\mathscr {G}$ farther downstream, such that it targets the area at the trailing edge of the leaflets. This area was identified by Zolfaghari & Obrist (Reference Zolfaghari and Obrist2019) as the interaction zone where the vorticity waves produced by the ILEV interacted with the wake structures, to ultimately force their breakdown. It is important to explore the possibility of this breakdown, which visually seemed to be caused by ILEV, being additionally driven by other instabilities originating in the cavities or the wake itself.

DAL simulations were performed using a circular energy probe of radius $R=3.6$ mm located at $(\xi,r)=(14.76,0\ \textrm {mm})$. Due to an increased distance from the leading edge, larger times may be needed for revealing the leading-edge signature. For this reason, we performed looping simulations for the mid-leaflet case for a fourth terminal time of $T=0.2$, in addition to the previous values of $T=0.01$, $T=0.05$ and $T=0.1$.

Figure 18(a,b) shows the initial and terminal conditions corresponding to a maximum gain $G$ at trailing edge for $T=0.01$. Similar to the mid-leaflet mask, for the smallest time, the optimal initial condition is concentrated on the mask itself. However, by increasing the time horizon $T$, the instability is progressively driven by disturbances created farther upstream (see, from top to bottom the second, third and fourth rows of figure 18). Figure 18(e, f) shows that the optimal initial condition mainly consists of Orr structures linked to the shear profile that is formed downstream of the ILEV zone. These structures become skewed within the ILEV zone, due to the specific shape of the velocity profiles in this area (Zolfaghari & Obrist Reference Zolfaghari and Obrist2019). Figure 18(g,h) is particularly relevant, as it shows that, for sufficiently large $T$, the initial conditions for maximum energy growth near the trailing edge are mostly influenced by the leading-edge structures, including the ILEV instability.

Figure 18.

Streamwise component of the (a,c,e,g) optimal initial and (b,d,f,h) final states are demonstrated for a trailing-edge mask. From top to bottom, cases with $T=0.01$, $T=0.05$, $T=0.1$ and $T=0.2$ are shown, respectively. The domain within the circle identifies the mask $\varOmega _v$. Dark blue and light yellow areas signify minimum negative and maximum positive velocity disturbances, respectively.

4.3.3. Maximum energy growth in the wake region

Finally, we study the optimal energy growth in the wake area of the 2-D BMHV submodel. To this end, we place $\varOmega _v$ downstream of the trailing edge where growth can be influenced by upstream wake structures as well as other instability mechanisms such as those arising from the cavities.

The DAL simulations were performed using a circular energy probe of radius $R=3.6\ \textrm {mm}$ located at $(\xi,r)=(23.7,0\ \textrm {mm})$. Figure 19 shows the initial condition for maximum energy growth at this location. Given a sufficiently long time $T=0.4$, the maximum growth for this area is driven by an optimal initial condition originating at the leading edge. This is an important observation, as it verifies the crucial role the leading-edge design of the leaflet plays on the disturbance energy growth, triggering or promoting turbulent flow features in the wake of the valve. This phenomenon has been observed in Zolfaghari & Obrist (Reference Zolfaghari and Obrist2019), but is verified quantitatively here. A 3-D transient growth analysis, regardless of prohibitive computational costs, is unlikely to change the key outcome here, that is, the energy growth in the wake is promoted mainly by ILEV mechanism. This is supported by DNS evidence (see figure 7) where elimination of ILEV significantly reduced the turbulence in the wake of the valve, while other instability mechanisms were present.

Figure 19.

Initial (a) and final (b) states for a wake mask and a terminal time $T=0.4$. The domain within the circle identifies the mask $\varOmega _v$. Dark blue and light yellow areas indicate the minimum negative and maximum positive velocity disturbances, respectively.