*2.4. Modelling Geometric Imperfections in Our DNS*

Section 3.6 presents results from DNS in which we mimiced the geometric perturbation of the experiments of Xu et al. [17]. Inspired by the optimal baffle designed by Marensi et al. [24], we here model the effect of geometric perturbations with an additional volume force in Equation (2) of the form

$$F\_p(r, \theta, z, t) = -A\_p \cdot f\_p(r, \theta, z) \cdot \mathbf{u}(r, \theta, z, t). \tag{3}$$

The body force *<sup>F</sup><sup>p</sup>* acts against the velocity field *<sup>u</sup>* and is localised in the radial, azimuthal, and axial direction by

$$f\_p(r, \theta, z) = f(r) \cdot \lg(\theta, z) \cdot h(z) \text{with} \tag{4}$$

$$f(r) = \frac{1}{2} + \frac{1}{\pi} \arctan(M\_7(r - r\_0)),\tag{5}$$

$$\log(\theta, z) = \frac{1}{\pi} (\arctan(M\_0(\theta - \pi(\theta\_0(z) - L\_0))) - \arctan(M\_0(\theta - \pi(\theta\_0(z) + L\_0)))),\tag{6}$$

$$h(z) = \frac{1}{\pi} \Big( \arctan \left( M\_z \left( z - z\_0 + \frac{L\_z}{2} \right) \right) - \arctan \left( M\_z \left( z - z\_0 - \frac{L\_z}{2} \right) \right) \Big) \text{and} \tag{7}$$

*Entropy* **2021**, *23*, 46

$$
\partial \partial \mathfrak{u}(z) = 1 + 2\Lambda \theta \frac{(z - z\_0)}{L\_z}. \tag{8}
$$

These localisation functions satisfy the constraints max(*f*) = 1, min(*f*) = 0, max(*g*) = 1, min(*g*) = 0, max(*h*) = 1, and min(*h*) = 0; the perturbation amplitude is given by *Ap*.

Due to the big parametric space in hand, we designed three simple body force set-ups and left further optimisation of parameters as future work. The first set-up is an axisymmetric force that models the effect of a small circumferential contraction similar to weak stenosis in blood vessels [25] or imperfect pipe joints in laboratory experiments [17] (see Figure 4a). The second set-up is a highly localised force that approximates the effect of a single bump or an individual roughness element (see Figure 4b). The third set-up is also a highly localised force that approximates the effect of a single bump or an individual roughness element, but this time, it is tilted with respect to the axial direction (see Figure 4c). The parameters defining the perturbations are given in Table 1. We studied the effect of the axisymmetric force on steady laminar Hagen–Poiseuille flow at *Re* = 2400 to select a suitable value of the force amplitude *Ap*. Our criterion was that the force must be strong enough to sufficiently disturb the flow without creating too long of a re-circulation region. For *Ap* = 0.25, we found a fair compromise between these two constraints.

**Figure 4.** Geometric representation of the perturbation force (*Fp*) in terms of iso-surfaces (black) of the localisation function for *fp* = 0.5. (**a**) Axisymmetric contraction. (**b**) Localised bump. (**c**) Tilted bump. See Table 1 for details. The direction of the mean bulk flow (*us*) is always from left to right.

**Table 1.** Parameters to control the body force term in Equation (3) to model the effect of geometric perturbations: Magnitude (*Ap*) and slope (*M*), size (*L*), and location in the radial (*r*), azimuthal (θ), and axial (*z*) direction. Geometric representations of the perturbations are shown in Figure 4.


The goal of this model is to serve as a proof of concept. Our hypothesis is that geometric imperfections employed in the experiments locally modify the flow pattern causing the instability. The model satisfies this requirement, as it represents a small perturbation to the flow. It is meant for testing such a hypothesis, whereas the precise shape of its geometry plays an ancillary role. In order to faithfully reproduce the experiments of Xu et al. [17], one would need to have a boundary-fitted mesh or use immersed boundary methods. We are, however, confident that if the DNS was exactly reproducing the precise imperfections of the experiments, the exact same behaviour would be observed in the DNS.

#### **3. Results**

We first tested the effect of the pulsation on spatiotemporal intermittency by performing a DNS initialised with a snapshot of the statistically steady pipe flow (SSPF), as shown in Figure 2a. We refer to these simulations as IC SSPF. Next, we followed Xu et al. [20] and performed a linear non-modal stability analysis to identify the optimal perturbation for the parameter values of interest (*Re* = 2400, *Wo* = 8, and several *A*). This method produces the geometry (radial shape and axial/azimuthal wavenumbers) and the initial time (*t*0) of the perturbation achieving the maximum energy amplification. We used these optimal perturbations on top of the Sexl–Womersley velocity profile as initial conditions for a second set of DNS in order to test whether puffs or helical waves were developed. We refer to these simulations as IC SWOP. In a last step, we performed a third set of DNS with the body force term in Equation (3) to mimic the experimental setup of Xu et al. [17]. All parameter combinations for which we have performed DNS are summarised in Figure 5b.

**Figure 5.** Turbulent fraction (*Ft*) in the computational pipe domain based on the axial vorticity data shown in Figure 2 and Figure 7. The threshold to distinguish turbulent from laminar regions is set to ω2 *z <sup>r</sup>*,<sup>θ</sup> <sup>=</sup> <sup>4</sup> <sup>×</sup> <sup>10</sup><sup>−</sup>2. (**a**) Time series of the turbulent fraction for several amplitudes *<sup>A</sup>* (line styles) and different numerical set-ups (symbols and colours from those in (**b**)). (**b**) Time-averaged turbulent fraction *Ftt*><sup>2</sup> for four different set-ups: The statistically steady pipe flow (SSPF) serves as reference data and as initial condition (IC) for the first set-up. The IC for the second set-up are composed out of the analytical Sexl–Womersley (SW) velocity profile superimposed with an optimal perturbation (OP). The third set-up is initialised with an unperturbed SW flow and then permanently perturbed using a localised body force (see Section 3.6).
