2. Problem description

We consider a simplified, 2-D representation of inflow and outflow to/from a rectangular cavity, though we also briefly explore channel optimisation for different cavity geometries in § 4. This set-up is motivated by an idealisation of a ureteroscope tip within the renal pelvis, but more generally provides a toy geometry for considering shape optimisation on flows characterised by large vortical structures. We employ a Cartesian coordinate system $(x^{\star }, y^{\star })$ (stars indicate dimensional variables throughout) with corresponding coordinate directions $\boldsymbol {i}$ and $\boldsymbol {j}$.

The flow domain is represented by the light pink and light orange regions in figure 1, i.e. $\varOmega _{{t}}^{\star }\cup \varOmega _{{c}}^{\star }$. The $\varOmega _{{t}}^{\star }$ domain encompasses both inflow and outflow channels. The inflow boundary to the working channel is $\varGamma ^{\star }_{{in}}$, the upper and lower outflow boundaries are $\varGamma ^{\star }_{{out,u}}$ and $\varGamma ^{\star }_{{out,l}}$, respectively, $\varGamma ^{\star }_{{f}}$ is the solid boundary we will allow to vary in the shape optimisation and $\varGamma ^{\star }_{{wall}}$ denotes all other (fixed) solid boundaries. Following on from our previous experimental and numerical work in a similar geometry (Williams et al. Reference Williams, Castrejon-Pita, Turney, Farrell, Tavener, Moulton and Waters2020), we use the parameter values $l_1=1.2743\ \text {cm}$ and $l_2 = 1\ \text {cm}$ for the channel and cavity lengths, respectively, and set $a = 0.06\ \text {cm}$, $b = 0.25701\ \text {cm}$ and $h=0.11598\ \text {cm}$ in the initial geometry.

Figure 1.

Problem set-up. Fixed boundaries are $\varGamma ^{\star }_{{out,u}}$, $\varGamma ^{\star }_{{out,l}}$, $\varGamma ^{\star }_{{in}}$ and $\varGamma ^{\star }_{{wall}}$; the red boundary $\varGamma ^{\star }_{{f}}$ is the control boundary that we allow to vary in the shape optimisation.

The dimensional velocity field is given by $\boldsymbol {u}^{\star } = u^{\star }\boldsymbol {i} + v^{\star }\boldsymbol {j}$ and pressure $p^{\star }$. The flow is governed by the steady, incompressible Navier–Stokes and continuity equations, subject to the following boundary conditions. The flow is driven by an imposed fully developed parabolic profile at the inlet boundary with maximum velocity $U$. As outflow conditions, we impose zero normal stress and parallel velocity. On all other boundaries (solid black and red lines in figure 1) we assume no-slip and no-normal-velocity conditions.

A passive tracer of concentration $c^{\star }(\boldsymbol {x}^{\star }, t^{\star })$ over time $t^{\star }$ is advected by the flow and also diffuses with constant diffusion coefficient $\mathcal {D}$. We impose no total flux (which has both advective and diffusive components) of the tracer through the inflow or solid walls ($\varGamma ^{\star }_{{in}},\varGamma ^{\star }_{{wall}}, \varGamma ^{\star }_{{f}}$) and no diffusive flux of tracer through the outflows ($\varGamma ^{\star }_{{out,u}}$, $\varGamma ^{\star }_{{out,l}}$). We assume that initially the passive tracer has constant concentration $C$ that is uniform within the cavity and zero elsewhere.

2.1. Dimensionless governing equations

We non-dimensionalise as follows:

(2.1a–e)\begin{equation} \boldsymbol{x} =\frac{\boldsymbol{x}^{\star}}{a}, \quad t = \frac{Ut^{\star}}{a}, \quad \boldsymbol{u} = \frac{\boldsymbol{u}^{\star}}{U}, \quad p = \frac{a(p^{\star}-p_{{ref}})}{U\mu},\quad c = \frac{c^{\star}}{C}, \end{equation}

where $p_{{ref}}$ is the pressure at the outflow boundaries $\varGamma ^{\star }_{{out,u}}$ and $\varGamma ^{\star }_{{out,l}}$ of the truncated domain, $\mu$ is the dynamic viscosity of the fluid and we note that time has been non-dimensionalised using the advective time scale. The steady, incompressible Navier–Stokes equations are

(2.2a,b)\begin{equation} Re(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}= -\boldsymbol{\nabla} p + \boldsymbol{\nabla}\boldsymbol{\cdot}\mathcal{E}(\boldsymbol{u}), \quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0 \quad\text{in }\varOmega_{{c}}\cup\varOmega_{{t}} \end{equation}

where the Reynolds number $Re = \rho U a/\mu$, with $\rho$ the fluid density, $\mathcal {E}(\boldsymbol {u}) = \boldsymbol {\nabla }\boldsymbol {u}+(\boldsymbol {\nabla }\boldsymbol {u})^\top$ and $\boldsymbol {\nabla } = (\partial /\partial x, \partial /\partial y)$. (The long formulation of (2.2a), rather than reducing $\boldsymbol {\nabla }\boldsymbol {\cdot }\mathcal {E}(\boldsymbol {u})$ to $\nabla ^2\boldsymbol {u}$ using (2.2b), ensures zero stress as the natural boundary condition in the weak formulation of the equations.) The flow boundary conditions are

(2.3a,b)\begin{gather} u = \left[1-y^2\right], \quad v = 0 \quad \text{on }\varGamma_{{in}}, \end{gather}

(2.3c,d)\begin{gather} 2\frac{\partial u}{\partial x} - p = 0 \quad v = 0 \quad \text{on }\ \varGamma_{{out,l}}\cup\varGamma_{{out,u}}, \end{gather}

(2.3e)\begin{gather} u = v = 0 \quad \text{on }\varGamma_{{wall}}\cup\varGamma_{{f}}. \end{gather}

The advection–diffusion equation for tracer concentration is

(2.4)\begin{equation} \frac{\partial c}{\partial t} = Re^{-1}Sc^{-1}\nabla^2c - \boldsymbol{u}\boldsymbol{\nabla} c, \end{equation}

where the Schmidt number $Sc = \mu /(\rho \mathcal {D})$ represents the ratio of momentum diffusivity to tracer diffusivity and we recall that $\mathcal{D}$ is the tracer diffusion coefficient. Equation (2.4) is solved subject to the boundary conditions

(2.5a,b)\begin{equation} \boldsymbol{J}\boldsymbol{\cdot}\boldsymbol{n} = 0 \quad\text{on } \varGamma_{{in}},\varGamma_{{wall}}, \varGamma_{{f}}, \quad \boldsymbol{\nabla} c\boldsymbol{\cdot}\boldsymbol{n} = 0 \quad\text{on } \varGamma_{out,u}, \varGamma_{out,l}, \end{equation}

where $\boldsymbol {J} = -Sc^{-1}Re^{-1}\boldsymbol {\nabla } c + \boldsymbol {u}c$. Additionally, (2.4) is solved subject to the initial condition

(2.6)\begin{equation} c(\boldsymbol{x},0) = \begin{cases} 1 & \text{in }\varOmega_{{c}}, \\ 0 & \text{in }\varOmega_{{t}}. \end{cases} \end{equation}

2.2. Washout time

Following Williams et al. (Reference Williams, Castrejon-Pita, Turney, Farrell, Tavener, Moulton and Waters2020), to quantify the rate at which the passive tracer exits the entire domain $\varOmega _{{c}}\cup \varOmega _{{t}}$, we define the percentage loss at time $t$ as

(2.7)\begin{equation} \% \text{ loss} = \frac{\iint_{{\varOmega}_{{c}}\cup\varOmega_{{t}}}\left[ c({\boldsymbol{x}},0)-c(\boldsymbol{x}, t)\right]\,\mathrm{d}{\varOmega}}{\iint_{{\varOmega}_{{c}}\cup\varOmega_{{t}}} c(\boldsymbol{x}, 0)\,\mathrm{d}{\varOmega}}, \end{equation}

and we define a metric for the washout time, $T_{90}$, corresponding to the time required for 90 % of the tracer to exit the domain, i.e.

(2.8)\begin{equation} \frac{\iint_{{\varOmega}_{{c}}\cup\varOmega_{{t}}}\left[ c(\boldsymbol{x},0)-c(\boldsymbol{x},T_{90})\right]\,\mathrm{d}{\varOmega}}{\iint_{{\varOmega}_{{c}}\cup\varOmega_{{t}}} c(\boldsymbol{x},0)\,\mathrm{d}{\varOmega}} = 0.9. \end{equation}

2.3. Optimisation objectives

As discussed in § 1, vortices within a flow field are associated with regions where $\det \boldsymbol {\nabla }\boldsymbol {u} > 0$ (Jeong & Hussain Reference Jeong and Hussain1995). Therefore,

(2.9)\begin{equation} \iint_{\varOmega_{{c}}} \max(0, \det\boldsymbol{\nabla}\boldsymbol{u}) \,\mathrm{d}{\varOmega} \end{equation}

is a logical objective to associate with recirculation zones. However, to ensure differentiability of the objective and enhance optimisation performance, we utilise the smoothing function

(2.10)\begin{equation} g(s; \epsilon) = \begin{cases} 0, & s \leq 0, \\ s^3/(s^2 + \epsilon), & s > 0, \end{cases} \end{equation}

and consider minimising

(2.11)\begin{equation} \iint_{\varOmega_{{c}}} g(\det\boldsymbol{\nabla}\boldsymbol{u}; \epsilon) \,\mathrm{d}{\varOmega}, \end{equation}

as introduced by Kasumba & Kunisch (Reference Kasumba and Kunisch2010, Reference Kasumba and Kunisch2012). We choose $\epsilon = 0.01$ – as opposed to $\epsilon = 1$ used by Kasumba & Kunisch (Reference Kasumba and Kunisch2010, Reference Kasumba and Kunisch2012) – to well capture vortices within the flow. Contours of $g(\det \boldsymbol {\nabla }\boldsymbol {u}; 0.01)$ for $Re = 10$ and $Re = 20$ are plotted in figures 2(b) and 2(d), respectively. Comparison with corresponding streamlines for $Re = 10$ and $Re = 20$ shown in figures 2(a) and 2(c), respectively, demonstrates strong qualitative correlation between $g(\det \boldsymbol {\nabla }\boldsymbol {u}; 0.01)$ and areas of closed streamlines.

Figure 2.

(a,c) Streamlines in the cavity for $Re = 10$, $20$, respectively, and (b,d) corresponding contours of $g(\det \boldsymbol {\nabla }\boldsymbol {u}; 0.01)$, where $g(s;\epsilon )$ is given by (2.10): $(a)$ $Re = 10$, streamlines; $(b)$ $Re = 10$, $g(\det \boldsymbol {\nabla }\boldsymbol {u}; 0.01)$; $(c)$ $Re = 20$, streamlines; $(d)$ $Re = 20$, $g(\det \boldsymbol {\nabla }\boldsymbol {u}; 0.01)$.

2.4. Numerical details

2.4.1. State equation (Navier–Stokes) solver

Our implementation is based on the Firedrake finite element library (Rathgeber et al. Reference Rathgeber, Ham, Mitchell, Lange, Luporini, McRae, Bercea, Markall and Kelly2016; Homolya, Kirby & Ham Reference Homolya, Kirby and Ham2017; Luporini, Ham & Kelly Reference Luporini, Ham and Kelly2017; Homolya et al. Reference Homolya, Mitchell, Luporini and Ham2018), which uses the unified form language (UFL) (Alnæs et al. Reference Alnæs, Logg, Ølgaard, Rognes and Wells2014). UFL is a domain specific language embedded in Python that enables the expression of partial differential equations in their weak form in close to mathematical notation. We discretise the equations using continuous, piecewise-quadratic Lagrange finite elements for the velocity field and discontinuous, piecewise-affine Lagrange finite elements for the pressure. This element pair is also referred to as the quadratic Scott–Vogelius element, with the velocity space denoted by $V_h$ and the pressure space by $Q_h$. While the Scott–Vogelius element pair is not known to be inf-sup stable on arbitrary triangulations (see Guzmán & Scott (Reference Guzmán and Scott2018) for a recent discussion), it has been shown to be inf-sup stable on meshes with a particular structure. To obtain an apposite mesh structure, we perform a barycentric refinement of an unstructured mesh generated by the mesh generator software Gmsh (Geuzaine & Remacle Reference Geuzaine and Remacle2009); this guarantees inf-sup stability (Qin Reference Qin1994). We choose the Scott–Vogelius element as it enforces the incompressibility constraint exactly, a property that more recently has been understood to be of importance for the accuracy of numerical solutions to incompressible flow (John et al. Reference John, Linke, Merdon, Neilan and Rebholz2017).

Since the Navier–Stokes equations are nonlinear, we employ a Newton scheme. The Jacobian required for this is derived automatically and symbolically by the UFL library. To assemble the linear systems arising in each step of the Newton iteration, Firedrake generates highly optimised C code at runtime. The system is then solved using the MUMPS (Amestoy et al. Reference Amestoy, Duff, L'Excellent and Koster2000) sparse direct solver. For higher accuracy, we perform a few (usually two) iterations of PETSc's (Balay et al. Reference Balay, Gropp, McInnes and Smith1997, Reference Balay2016, Reference Balay2018; Dalcin et al. Reference Dalcin, Paz, Kler and Cosimo2011) GMRES implementation using the LU factorisation as a preconditioner, until the residual is reduced to $10^{-10}$.

Starting from a zero initial guess, the solver finds a symmetric solution. However, for large enough Reynolds number, this solution is not stable. For these cases we perform a first solve with a no-slip boundary condition applied to one of the outlets (we use $\varGamma_{out, u}$ although identical, albeit mirrored, results would be obtained if we used $\varGamma_{out, l}$), and use this solution as the initial guess, to guarantee that the stable, asymmetric solution branch is found and subsequently used for the initial step of the optimisation routine (Williams et al. Reference Williams, Castrejon-Pita, Turney, Farrell, Tavener, Moulton and Waters2020).

2.4.2. Optimisation

For the optimisation we rely on the Fireshape shape optimisation library (Paganini & Wechsung Reference Paganini and Wechsung2020). Fireshape is an open source library that connects the Firedrake finite element library with the Trilinos Rapid Optimisation Library (Ridzal & Kouri Reference Ridzal and Kouri2014). Fireshape employs the perturbation-of-identity approach (Micheletti Reference Micheletti1972; Murat & Simon Reference Murat and Simon1976; Simon Reference Simon1980; Delfour & Zolésio Reference Delfour and Zolésio2011), which is based on deforming an initial domain $\varOmega _{{t}}$ using deformations of the form $\boldsymbol {T}=\boldsymbol {Id} + \boldsymbol {X}$, where $\boldsymbol {Id}(\boldsymbol {x}) = \boldsymbol {x}$ is the identity mapping. The optimisation is then formulated as finding an optimal deformation $\boldsymbol {X}$, and at each iteration the shape is given by $\varOmega _{{t}}^{(k)} = (\boldsymbol {Id} + \boldsymbol {X}^{(k)})(\varOmega _{{t}})$. Denoting by $X_h$ the space of continuous, piecewise-affine vector fields that vanish on $\partial \varOmega _{{t}} \setminus \varGamma _{f}$, the vanilla shape optimisation problem can then be expressed as

(2.12)\begin{equation} \left. \begin{aligned} \underset{{(\boldsymbol{u}, p)\in V_h\times Q_h, \boldsymbol{X}\in X_h}}{\text{minimise}} & \quad J_e(\varOmega_{{c}}) = \iint_{\varOmega_{{c}}} g(\det \boldsymbol{\nabla} \boldsymbol{u}, \epsilon) \,\mathrm{d}{\varOmega}\\ \text{subject to} & \quad (\boldsymbol{u}, p) \text{ satisfy } (2.2) \text{ in } (\boldsymbol{Id} + \boldsymbol{X})(\varOmega_{t}) \cup \varOmega_{c}. \end{aligned} \right\} \end{equation}

We remark that while the integration domain $\varOmega _{{c}}$ in the objective stays fixed, it is the velocity field in the integrand that changes as the channel domain $\varOmega _{{t}}$ is deformed.

Both analysing and solving shape optimisation problems are known to be difficult for several reasons. First, shape objectives are usually highly nonlinear and non-convex in the control $\boldsymbol {X}$, so a unique optimal shape often does not exist. Second, when deforming meshes one needs to ensure that the mesh does not degenerate or overlap. Fireshape has several capabilities that allow us to mitigate these problems.

To avoid the first issue, we add a Tikhonov-type term to the objective, that is, a term of the form $\alpha \|\boldsymbol {X}\|_{X_h}^2$ (we choose $\alpha =0.01$ in our simulations). Such a term is commonly used for ill-posed optimisation problems as it is strictly convex. The norm $\|\cdot \|_{X_h}^2$ on the space $X_h$ is chosen to penalise non-smoothness of the deformations,

(2.13)\begin{equation} \left. \begin{gathered} \|\boldsymbol{X}\|_{X_h}^2 := (\boldsymbol{X}, \boldsymbol{X})_{X_h},\\ (\boldsymbol{X}, \boldsymbol{Y})_{X_h} := \iint_{\varOmega_{{t}}} \mathcal{E} (\boldsymbol{X})\colon \mathcal{E} (\boldsymbol{Y}) + 0.001 \Delta \boldsymbol{X} \colon \Delta \boldsymbol{Y} \,\mathrm{d}{\varOmega}. \end{gathered} \right\} \end{equation}

Since standard finite elements only admit a first (weak) derivative, we employ a standard interior penalty method (Brenner & Sung Reference Brenner and Sung2005; Wells, Kuhl & Garikipati Reference Wells, Kuhl and Garikipati2006). To counteract the second issue, we add a term $\beta \| \max (0, \sigma _{{max}}(\boldsymbol {\nabla }\boldsymbol {X}) - 0.9)\|_{L^2}^2$ to the objective, which penalises large values of the largest singular value of the gradient of $\boldsymbol {X}$ (we choose $\beta =10^3$). This modification is motivated by a result that guarantees that for a convex domain, the deformed mesh does not overlap if the largest singular value of $\boldsymbol {\nabla }\boldsymbol {X}$ is smaller than 1 (Delfour & Zolésio (Reference Delfour and Zolésio2011), § 3 and Theorem 2.15). Finally, we note that any shape objective (such as $J_e$) is invariant with respect to deformations of the interior mesh nodes: deforming the interior of the mesh does not change the overall shape and hence does not have an impact on the objective. To remove this invariance, we require that the deformation satisfy

(2.14)\begin{equation} (\boldsymbol{X}, \boldsymbol{Y})_{X_h} = 0 \quad \text{for all } \boldsymbol{Y} \in X_{h, 0} = \{\boldsymbol{Y}\in X_h: \boldsymbol{Y}\vert_{\partial\varOmega_{{t}}}=0\}. \end{equation}

This is the weak formulation of requiring $\boldsymbol {\nabla }\boldsymbol {\cdot } \mathcal {E} (\boldsymbol {X}) + \gamma \Delta ^2 \boldsymbol {X}=0$ inside the domain that is free to move. We refer to Wechsung (Reference Wechsung2019, § 8.2) for more details on these modifications. To summarise, the minimisation problem that we solve is given by

(2.15)\begin{equation} \left. \begin{aligned} \underset{{(\boldsymbol{u}, p)\in V_h\times Q_h, \boldsymbol{X}\in X_h}}{\text{minimise}} & \quad J_e(\varOmega_{c}) + \alpha \|\boldsymbol{X}\|_{X_h}^2 + \beta \| \max(0, \sigma_{{max}}\boldsymbol{\nabla}\boldsymbol{X} - 0.9)\|_{L^2}^2\\ \text{subject to} & \quad (\boldsymbol{u}, p) \text{ satisfy } (2.2)\text{ in } (\boldsymbol{Id} + \boldsymbol{X})(\varOmega_{t})\cup\varOmega_{c},\\ & \quad \boldsymbol{X} \text{ satisfies } (2.14). \end{aligned} \right\} \end{equation}

In order to solve the optimisation problem efficiently the shape derivative is required. Fireshape relies on UFL to automatically derive the required adjoint equations and shape derivatives (Alnæs et al. Reference Alnæs, Logg, Ølgaard, Rognes and Wells2014; Ham et al. Reference Ham, Mitchell, Paganini and Wechsung2019). Finally, the resulting optimisation problem is solved using the L-BFGS implementation of the Rapid Optimization Library (Nocedal & Wright Reference Nocedal and Wright2006; Ridzal & Kouri Reference Ridzal and Kouri2014).

2.4.3. Tracer washout equations

We approximated the concentration function with piecewise-linear elements on the same mesh as the velocity field. The velocity field was solved for and subsequently restored to drive the advective term in (2.4). Equation (2.4) was discretised in time using a Crank–Nicolson scheme, and we added a mesh-dependent streamline upwind Petrov–Galerkin (SUPG) stabilisation term, as described in Franca, Frey & Hughes (Reference Franca, Frey and Hughes1992), as standard Galerkin discretisations of advection–diffusion equations are oscillatory in the advection-dominated regime (Silvester, Elman & Wathen Reference Silvester, Elman and Wathen2014). We validated our advection–diffusion code using the method of manufactured solutions, with details presented in Williams et al. (Reference Williams, Castrejon-Pita, Turney, Farrell, Tavener, Moulton and Waters2020, § B.2).

4. Summary and discussion

Many biomedical procedures require efficient clearance of debris, whether to improve visualisation in endoscopic surgeries, or to reduce accumulation of undesirable biological agents in a stented artery or ureter. Debris clearance typically relies on fluid flow – an injected saline solution, contrast agent or blood circulation – and debris aggregation sites can arise as a result of recirculation zones in the fluid. Modelling debris dynamics as the advection and diffusion of a passive tracer, the time required for the tracer to escape a region of closed streamlines depends nonlinearly on the velocity, diffusive properties of the debris and vortex size (Rhines & Young Reference Rhines and Young1983). This motivated the hypothesis that underpins the work in this article: for fixed flow speed and debris characteristics, a reduction in vortex size will promote more efficient debris clearance.

In this paper we considered flow in a 2-D Cartesian geometry, consisting of a rectangular cavity and inlet and outlet channels, and used shape optimisation to determine channel geometries that reduce the presence of vortices. We characterised recirculation zones via $\det (\boldsymbol {\nabla }\boldsymbol {u})>0$ – i.e. locally circular or spiral streamlines – thereby choosing the objective for our optimisation to be a smoothed form of $\max (0, \det \boldsymbol {\nabla }\boldsymbol {u})$ integrated over the cavity domain (Kasumba & Kunisch Reference Kasumba and Kunisch2012). We determined the relationship between washout time and the objective by solving an advection–diffusion equation for the motion of a passive tracer in both the initial and the optimised flow fields and calculating the resulting washout time. We used $T_{90}$ as the metric characterising washout time, corresponding to the time taken for 90 % of the tracer to exit the domain. We note here that alternatives to $T_{90}$ for the washout metric can be found in the literature: for example, Min et al. (Reference Min, Fischer and Pearlstein2020) track the location and magnitude of the peak of the remaining concentration. A comparison of washout metrics in relation to our objective function would be an interesting future study, but here we chose to restrict attention to $T_{90}$ as this is a commonly used standard measure of washout efficiency.

We first solved for washout time in the initial uniform channel geometry for Reynolds numbers $Re = 10$, $20$, $30$, $40$ and $50$ and confirmed that washout time depends non-trivially on Reynolds number, flow structure and the diffusion coefficient associated with the passive tracer. When the Schmidt number $Sc = 1$ – i.e. when tracer motion is governed by a balance of momentum diffusivity of the fluid and diffusivity of the concentration – washout time decreases monotonically with increasing Reynolds number. However, for decreased tracer diffusivity, we found that the case of fastest flow ($Re = 50$) in fact provided the slowest washout time. This seemingly counterintuitive result is easily explained by the fact that the fastest flow also creates the largest vortex, and in an advection-dominated regime vortex size has a dominant impact on washout time.

Running the optimisation algorithm for each of the five Reynolds numbers, we computed the optimal shape and found that washout time is indeed lower in the optimised geometries than in the initial geometries in each case. This partially confirmed the hypothesis that our objective serves as a proxy for washout time. To explore the correlation more thoroughly, we computed washout time for several geometries (and corresponding flow fields) along the optimisation pathway. With $Re = 20$ and $Re = 30$, washout time decreases monotonically with each shape iteration. For $Re = 10$, $40$ and $50$, washout time is not monotonic with iteration number in the optimisation algorithm, although after $200$ iterations the washout time is reduced to approximately $90\,\%$ of its initial value. This strong, if imperfect, correlation suggests that minimising washout may be well accomplished via a proxy objective measuring the size of vortical structures, a choice that creates significant computational advantage. This reflects the intuitive notion that if tracer transport is not diffusion dominated, tracer can remain ‘stuck’ in vortices for long periods of time, and thus reducing vortical structures can significantly reduce washout time. By construction, $Sc$ does not appear in the optimisation algorithm, but the connection between vortex size and washout suggests that the optimal shape for a given $Re$ will not vary significantly with varying $Sc$, at least in the advection-dominated regime. The connection between washout reduction and the size of vortical structures is less strong in the high-diffusion regime; here maximising flow in the non-vortical ‘racetrack’ regions may be more important than minimising the size of vortices, and a different proxy metric may be more suitable. However, such a regime is less interesting from an optimisation point of view, as it is characterised by very low washout times in any case.

Also of note on the computational side is the fact that, for two of the Reynolds numbers ($Re = 20$ and $Re = 30$), solution branch switching occurred during the optimisation routine. This feature was shown to occur as an artefact of the evolving broken pitchfork bifurcation structure for the flow problem during the evolution of the channel shape following the optimisation pathway. Such a feature may naturally arise in an optimisation scheme involving multi-solution flow fields, and highlights an important point, particularly when translating to a clinical or industrial setting: the flow field associated with a computed optimal shape may not be the only stable solution, so care may be needed in physically realising the desired optimal flow.

We also explored the washout time corresponding to the optimised channel shapes as a function of the diffusion of the passive tracer, captured by the Schmidt number $Sc$. As was observed in the base geometry, when plotting washout time against $Sc$ across the set of Reynolds numbers, numerous crossings of the lines are observed. This further underscores the complex relation between washout time, flow speed, vortices and diffusivity. Of course, from a device design point of view, the key question is what shape should the inflow and outflow channels be? Here, in lieu of a single answer, our analysis provides several points of insight: we have shown that the optimal shape is somewhat dependent on the Reynolds number, and the optimal Reynolds number is dependent on the tracer (debris) diffusivity. The diffusivity may not be well known, but we see the importance of a good estimate. In terms of Reynolds number, in practice this can be controlled by changing the upstream driving pressure or fluid properties, in which case figure 5 (or its equivalent in a different geometry) enables us to determine the best combination of shape and Reynolds number. A natural extension would be to optimise for both the Reynolds number and the shape simultaneously. However, in addition to non-trivial extensions of the Fireshape library being necessary for such a combined optimisation, we also expect more difficulties navigating the multiple solution branches if the Reynolds number changes during the optimisation. For these reasons, we do not explore this combined optimisation at this stage.

Even in the absence of a global optimal shape, figure 5 shows a clear qualitative trend that curving the jet of fluid via bending the inlet channel leads to a decrease in the size of vortical structures, and we anticipate that this qualitative finding could carry over to more complex (and possibly 3-D) domain geometries. To test this idea, in figure 7 we have run the optimisation algorithm on a set of three distinct geometries: in figure 7(a,d) we replace the square domain with a smooth curved boundary; in figure 7(b,e) we replace the back edge with an inwardly curved boundary; and in figure 7(c,e) we break the symmetry of the domain through the addition of a circular bulge in the corner. Here, figure 7(a–c) shows the flow pattern and computed $T_{90}$ for the straight inlet/outlet channels, while figure 7(d–f) shows the optimised shape and corresponding $T_{90}$. For each domain, the total cavity area was preserved, allowing direct comparison of the washout times. Unsurprisingly, the smooth domain in figure 7(a,d) exhibits a much faster washout than the domain in figure 7(c,f), in which debris remains stuck in the small vortex in the additional ‘mini-cavity’. Nevertheless, in all cases we observe a significant reduction in washout time. Perhaps more importantly, the optimised shape for each of the domains is qualitatively very similar, following the same trend of curving the tip of the inlet channel as was found in figure 5. While this is clearly not exhaustive in terms of possible geometric features, it does demonstrate a robustness to the framework and reinforces a design principle that curving the fluid can minimise the size of vortical structures and thus reduce the recirculation of debris.

Figure 7.

Optimisation on different geometries: (a–c) the flow pattern and computed $T_{90}$ before optimisation; (d–f) the same after optimisation. All computations are for $Re = 20$ and $Sc = 20$: $(a)$ $T_{90} = 1449.3$; $(b)$ $T_{90} = 1624.2$; $(c)$ $T_{90} = 2297.4$; $(d)$ $T_{90} = 896.0$; $(e)$ $T_{90} = 1028.9$; $(f)$ $T_{90} = 1815.2$.

An obvious extension to the study reported in this manuscript is to consider 3-D – rather than 2-D – geometries representative of real medical devices. However, a salient factor to consider is computational cost: solving shape optimisation problems involves many (${\sim }200$) state and adjoint equation solves. In two dimensions this is possible on a standard workstation, but the computational cost is several orders of magnitude larger in three dimensions. For this reason we focused on a 2-D case in this study, but are cognizant of the significance of the extension to three dimensions. We note that the choice of vortex metric is more open-ended in three dimensions, and a thorough comparison of results using different available metrics remains necessary future work (Jeong & Hussain Reference Jeong and Hussain1995; Haller Reference Haller2005).

Finally, we observe that in the optimisation scheme, the choice of which boundaries are fixed and which are free is important. In our set-up, we fixed the walls at the front of the scope, thereby preventing the inflow channel from flaring or tapering as it enters the cavity. Relaxing this restriction potentially introduces competing effects: a flared tip geometry may reduce the size of vortices, but it also reduces the flow speed of the open streamlines, or ‘racetrack’, as described in § 3.1. For high $Sc$ (tracer with low diffusivity), debris can get stuck in areas of slow velocity, inducing longer washout times, even though the presence of vortices has been reduced. In this instance, an alternative objective that promotes fast velocity in non-vortex regions – the racetrack – might be suitable, e.g. one could subtract a term from the objective of the form $\int _{\varOmega _{{c}}}|\boldsymbol {u}|\mathcal {H}(-\det \boldsymbol {\nabla }\boldsymbol {u})\,\textrm {d}\varOmega _{{c}}$. However, numerical complexities arise with determining a smoothed version of the indicator function $\mathcal {H}$, and further work is required to implement such a modified metric in practice.