Shape optimization of Stokesian peristaltic pumps using boundary integral methods

TL;DR

This paper introduces a boundary integral method for optimizing the shape of peristaltic pumps that transport viscous fluids, enabling efficient shape updates without volume remeshing and demonstrating improved performance through numerical examples.

Contribution

It develops a boundary integral approach for shape optimization of Stokesian peristaltic pumps, deriving shape derivatives and avoiding volume remeshing during optimization.

Findings

Significant cost savings in shape optimization process.

Effective boundary integral formulas for shape derivatives.

Successful numerical demonstrations of optimized pump shapes.

Abstract

This article presents a new boundary integral approach for finding optimal shapes of peristaltic pumps that transport a viscous fluid. Formulas for computing the shape derivatives of the standard cost functionals and constraints are derived. They involve evaluating physical variables (traction, pressure, etc.) on the boundary only. By emplyoing these formulas in conjuction with a boundary integral approach for solving forward and adjoint problems, we completely avoid the issue of volume remeshing when updating the pump shape as the optimization proceeds. This leads to significant cost savings and we demonstrate the performance on several numerical examples.

Equations

$${\mathbf{x}}_{,s}=\mathbf{\tau },{\mathbf{\tau }}_{,s}=\kappa \mathbf{n},{\mathbf{n}}_{,s}=-\kappa \mathbf{\tau }\text{on Γ+ and Γ−}.$$

$$\mathbf{U}(s)=U\mathbf{\tau }(s),$$

$$\mathbf{u}(\mathbf{x})={\mathbf{u}}^{\text{D}}(\mathbf{x}):=\frac{{\ell }^{\pm }}{L}c{\mathbf{\tau }}^{\pm }(\mathbf{x}),\mathbf{x}\in {\Gamma }^{\pm },$$

$$\text{mbox}div\mathbf{\sigma }[\mathbf{u},p]$$

$$\mathbf{\sigma }[\mathbf{u},p]$$

$$\text{mbox}div\mathbf{u}$$

$$\mathbf{u}$$

$$\text{find }\bm{u}\in\text{\boldmath$\mathcal{V}$},\ p\in\mathcal{P},\ \bm{f}\in\mbox{\boldmath$\mathcal{F}$},\qquad\left\{\begin{aligned} \text{(a) \ }&&a(\bm{u},\bm{v})-b(\bm{v},p)-\big{\langle}\hskip 1.00006pt\bm{f},\bm{v}\hskip 1.00006pt\big{\rangle}_{\Gamma}&=0&\qquad&\forall\bm{v}\in\text{\boldmath$\mathcal{V}$}\\ \text{(b) \ }&&b(\bm{u},q)&=0&&\forall q\in\mathcal{P}\\ \text{(c) \ }&&\big{\langle}\hskip 1.00006pt\bm{u}^{\text{\scriptsize D}},\bm{g}\hskip 1.00006pt\big{\rangle}_{\Gamma}-\big{\langle}\hskip 1.00006pt\bm{u},\bm{g}\hskip 1.00006pt\big{\rangle}_{\Gamma}&=0&&\forall\bm{g}\in\mbox{\boldmath$\mathcal{F}$}\end{aligned}\right.$$

$$a(\mathbf{u},\mathbf{v})={\int }_{\Omega }2\mu \mathbf{D}[\mathbf{u}]:\mathbf{D}[\mathbf{v}]\text{d}V,b(\mathbf{v},q)={\int }_{\Omega }q\text{mbox}div\mathbf{v}\text{d}V$$

$$J_{\text{\scriptsize PL}}(\Omega):=\big{\langle}\hskip 1.00006pt\bm{f}_{\Omega},(\bm{u}^{\text{\scriptsize D}}+c\bm{e}_{1})\hskip 1.00006pt\big{\rangle}_{\Gamma}.$$

$$C_{\text{V}}(\Omega ):=|\Omega |-V_0=0$$

$$C_{\text{Q}}(\Omega ):=Q(\Omega )-Q_0=0\text{with}Q(\Omega )=\frac{1}{L}{\int }_{\Omega }({\mathbf{u}}_{\Omega }+c{\mathbf{e}}_1)\cdot {\mathbf{e}}_1\text{d}V.$$

$$J_{\text{\scriptsize PL}}(\Omega)=\big{\langle}\hskip 1.00006pt\bm{f}_{\Omega},(\bm{u}^{\text{\scriptsize D}}+c\bm{e}_{1})\hskip 1.00006pt\big{\rangle}_{\Gamma}=a(\bm{u}_{\Omega},\bm{u}_{\Omega}).$$

$$\mathbf{x}\in {\Omega }_{\text{all}}\mapsto {\mathbf{x}}^{\eta }=\mathbf{x}+\eta \mathbf{\theta }(\mathbf{x}),$$

$$\mathcal{O}=\big{\{}\Omega\subset\Omega_{\text{all}},\text{ $\Omega$ is periodic and simply connected, }|\Omega|=V_{0},\ Q(\Omega)=Q_{0}\big{\}}.$$

$$\Theta=\big{\{}\bm{\theta}\in W^{1,\infty}(\Omega_{\text{all}})\ \text{such that \ (i) }\bm{\theta}|_{\Gamma_{0}}=\bm{\theta}|_{\Gamma_{L}},\ \text{(ii) }\bm{\theta}\!\cdot\!\bm{e}_{1}=0\text{ on $\Gamma_{0}$},\ \text{(iii) }\bm{\theta}(\bm{z}^{-})=\mathbf{0}\big{\}},$$

$$\makebox[5.28589pt]{$\buildrel\scriptscriptstyle\star\over{\bm{a}}$}(\bm{x})=\lim_{\eta\to 0}\dfrac{1}{\eta}\big{[}\hskip 1.00006pt\bm{a}(\bm{x}^{\eta},\eta)-\bm{a}(\bm{x},0)\hskip 1.00006pt\big{]}\qquad\bm{x}\in\Omega,$$

$$\text{(a) \ }(\bm{\nabla}\bm{a})^{\star}=\bm{\nabla}\makebox[5.28589pt]{$\buildrel\scriptscriptstyle\star\over{\bm{a}}$}-\bm{\nabla}\bm{a}\!\cdot\!\bm{\nabla}\bm{\theta},\qquad\text{(b) \ }(\mbox{div}\,\bm{a})^{\star}=\mbox{div}\,\makebox[5.28589pt]{$\buildrel\scriptscriptstyle\star\over{\bm{a}}$}-\bm{\nabla}\bm{a}\!:\!\bm{\nabla}\bm{\theta}$$

$$\big{\langle}\hskip 1.00006ptJ^{\prime}(\Omega),\bm{\theta}\hskip 1.00006pt\big{\rangle}=\lim_{\eta\to 0}\dfrac{1}{\eta}\big{(}\hskip 1.00006ptJ(\Omega_{\eta}(\bm{\theta}))-J(\Omega)\hskip 1.00006pt\big{)}.$$

$$\text{(a) }{I}_{\text{V}}(\eta )={\int }_{{\Omega }_{\eta }(\mathbf{\theta })}F(\cdot ,\eta )\text{d}V,\text{(b) }{I}_{\text{S}}(\eta )={\int }_{S_{\eta }(\mathbf{\theta })}F(\cdot ,\eta )\text{d}s,$$

$$\text{(a) \ }\dfrac{\text{d}I_{\text{\scriptsize V}}}{\text{d}\eta}\Big{|}_{\eta=0}=\int_{\Omega}\big{[}\hskip 1.00006pt\makebox[7.81946pt]{$\buildrel\scriptscriptstyle\star\over{F}$}+F(\cdot,0)\,\mbox{div}\,\bm{\theta}\hskip 1.00006pt\big{]}\;\text{d}V,\qquad\text{(b) \ }\dfrac{\text{d}I_{\text{\scriptsize S}}}{\text{d}\eta}\Big{|}_{\eta=0}=\int_{S}\big{[}\hskip 1.00006pt\makebox[7.81946pt]{$\buildrel\scriptscriptstyle\star\over{F}$}+F(\cdot,0)\,\mbox{div}_{\!S}\bm{\theta}\hskip 1.00006pt\big{]}\,\text{d}s,$$

$$\mbox{div}_{\!S}\bm{\theta}=\big{(}\hskip 1.00006pt\bm{I}-\bm{n}\otimes\bm{n}\hskip 1.00006pt\big{)}\!:\!\bm{\nabla}\bm{\theta}=\partial_{s}\theta_{s}-\kappa\theta_{n}$$

$$\mathbf{\theta }={\theta }_s\mathbf{\tau }+{\theta }_n\mathbf{n}.$$

$$\big{\langle}\hskip 1.00006ptJ^{\prime}(\Omega),\bm{\theta}\hskip 1.00006pt\big{\rangle}=\int_{\partial\Omega}g\,\theta_{n}\,\text{d}s$$

$$\big{\langle}\hskip 1.00006ptC_{\text{\scriptsize V}}^{\prime}(\Omega),\bm{\theta}\hskip 1.00006pt\big{\rangle}=\int_{\Omega}\mbox{div}\,\bm{\theta}\;\text{d}V=\int_{\partial\Omega}\theta_{n}\,\text{d}s=\int_{\Gamma}\theta_{n}\,\text{d}s,$$

$$\displaystyle a(\makebox[5.72458pt]{$\buildrel\scriptscriptstyle\star\over{\bm{u}}$},\bm{v})-b(\bm{v},\makebox[5.03125pt]{$\buildrel\scriptscriptstyle\star\over{p}$})-\big{\langle}\hskip 1.00006pt\makebox[5.97226pt]{$\buildrel\scriptscriptstyle\star\over{\bm{f}}$},\bm{v}\hskip 1.00006pt\big{\rangle}_{\Gamma}$$

$$\displaystyle b(\makebox[5.72458pt]{$\buildrel\scriptscriptstyle\star\over{\bm{u}}$},q)$$

$$\displaystyle\big{\langle}\hskip 1.00006pt\makebox[5.72458pt]{$\buildrel\scriptscriptstyle\star\over{\bm{u}}$}{}^{\text{\scriptsize D}},\bm{g}\hskip 1.00006pt\big{\rangle}_{\Gamma}-\big{\langle}\hskip 1.00006pt\makebox[5.72458pt]{$\buildrel\scriptscriptstyle\star\over{\bm{u}}$},\bm{g}\hskip 1.00006pt\big{\rangle}_{\Gamma}$$

$$a^1(\mathbf{u},\mathbf{v},\mathbf{\theta })$$

$$b^1(\mathbf{u},q,\mathbf{\theta })$$

$$\text{(a) \ }\makebox[5.72458pt]{$\buildrel\scriptscriptstyle\star\over{\bm{u}}$}{}^{\text{\scriptsize D}}=\frac{c\makebox[4.16669pt]{$\buildrel\scriptscriptstyle\star\over{\ell}$}}{L}\bm{\tau}+\frac{c\ell}{L}(\partial_{s}\theta_{n}+\kappa\theta_{s})\bm{n},\qquad\text{with \ (b) \ }\makebox[4.16669pt]{$\buildrel\scriptscriptstyle\star\over{\ell}$}=-\int_{0}^{\ell}\kappa\theta_{n}\,\text{d}s$$

$$a^{1}(\bm{u},\hat{\bm{u}}{},\bm{\theta})-b^{1}(\bm{u},\hat{p},\bm{\theta})-b^{1}(\hat{\bm{u}}{},p,\bm{\theta})\\ =\int_{\Gamma}\Big{\{}\hskip 1.00006pt(\bm{\sigma}[\bm{u},p]\!:\!\bm{D}[\hat{\bm{u}}{}])\theta_{n}-\bm{f}\!\cdot\!\bm{\nabla}\hat{\bm{u}}{}\!\cdot\!\bm{\theta}-\hat{\bm{f}}{}\!\cdot\!\bm{\nabla}\bm{u}\!\cdot\!\bm{\theta}\hskip 1.00006pt\Big{\}}\,\text{d}s+\int_{\Gamma_{L}}\Delta\hat{p}\,(\partial_{2}u_{1})\theta_{2}\,\text{d}s.$$

$$\big{\langle}\hskip 1.00006ptJ^{\prime}_{\text{\scriptsize PL}}(\Omega),\bm{\theta}\hskip 1.00006pt\big{\rangle}=\int_{\Gamma}\big{[}\hskip 1.00006pt\bm{f}\!\cdot\!\makebox[5.72458pt]{$\buildrel\scriptscriptstyle\star\over{\bm{u}}$}{}^{\text{\scriptsize D}}+(\bm{u}^{\text{\scriptsize D}}+\bm{e}_{1})\!\cdot\!\big{(}\hskip 1.00006pt\makebox[5.97226pt]{$\buildrel\scriptscriptstyle\star\over{\bm{f}}$}+\bm{f}\mbox{div}_{\!S}\bm{\theta}\hskip 1.00006pt\big{)}\hskip 1.00006pt\big{]}\,\text{d}s$$

Full text

1 Introduction

Many physiological flows are realized owing to peristalsis, a transport mechanism induced by periodic contraction waves in fluid-filled (or otherwise) tubes/vessels [12, 9, 3]. This mechanism is used in engineering applications like microfluidics, organ-on-a-chip devices and MEMS devices for transporting and mixing viscous fluids at small scale (e.g., see [27, 30, 25, 15, 32, 7]). Owing to its importance in science and technological applications, numerous analytical and numerical studies were carried out in the past decades to characterize the fluid dynamics of peristalsis in various physical scenarios; some recent works include [26, 28, 14] for non-Newtonian flows and [5, 1, 13] for particle transport.

Several applications—optimal transport of drug particles in blood flow [19], understanding sperm motility in the reproductive tract [24], propulsion of soft micro-swimmers [8, 23]—require scalable numerical methods that can handle arbitrary shape deformations. One of the challenges of existing mesh-based methods (e.g., finite element methods) is the high computational expense of re-meshing, needed to proceed between optimization updates or in transient solution of the forward/adjoint problems. Boundary integral equation (BIE) methods, on the other hand, avoid volume discretization altogether for linear partial differential equations (PDEs) and are highly scalable even for moving geometry problems [22]. While BIE methods have been used widely for shape optimization problems, including in linear elasticity, acoustics, electrostatics, electromagnetics and heat flow (e.g., [17, 33, 2, 31, 10]), we are not aware of their application to optimization of peristaltic pumps transporting simple (or complex) fluids.

The primary goal of this work is to derive shape sensitivity formulas that can be used in conjunction with an indirect BIE method to enable fast numerical optimization routines. Our work is inspired by that of Walker and Shelley [29] who considered the shape optimization of a peristaltic pump transporting a Newtonian fluid at low to moderate Reynolds numbers ($Re$) and applied a finite element discretization. Here, we restrict our attention to problems in the zero $Re$ limit only. Following [29], we consider shapes that minimize the input fluid power under constant volume and flow rate constraints. The forward problem requires solving the Stokes equations in a tube with periodic flow conditions and prescribed slip on the walls while the adjoint problem has a prescribed pressure drop condition across the tube. For both problems, we employ the recently developed periodic BIE solver of Marple et al. [18]; the primary advantage compared to classical BIE methods that rely on periodic Green’s functions is that the required slip or pressure-drop conditions can be applied directly.

The paper is organized as follows. In Section 2, we define the shape optimization problem and the PDE formulation of the forward problem in strong and weak forms. In Section 3, we derive the shape sensitivity formulas for a specific objective function and the functionals required to impose the given constraints on the pump shape. Using these formulas, we present a numerical optimization procedure in Section 4 based on a periodic boundary integral equation formulation for the PDE solves. We present validation and shape optimization results in Section 5 followed by conclusions in Section 6.

2 Problem formulation

2.1 Formulation of the wall motion

Pumping is achieved by the channel wall shape moving along the positive direction $\bm{e}_{1}$ at a constant velocity $c$, as a traveling wave of wavelength $L$ (the wave period therefore being $T=L/c$). The quantities $L,c$ (and hence $T$) are considered as fixed in the wall shape optimization process. This apparent shape motion is achieved by a suitable material motion of the wall, whose material is assumed to be flexible but inextensible. Like in [29], it is convenient to introduce a wave frame that moves along with the traveling wave, i.e. with velocity $c\bm{e}_{1}$ relative to the (fixed) lab frame.

Let then $\Omega$ denote, in the wave frame, the fluid region enclosed in one wavelength of the channel (see Fig. 1), whose boundary is $\partial\Omega=\Gamma\cup\Gamma_{p}$. The wall $\Gamma:=\Gamma^{+}\cup\Gamma^{-}$, which is fixed in this frame, has disconnected components $\Gamma^{\pm}$ which are not required to achieve symmetry with respect to the $x_{1}$ axis and have respective lengths $\ell^{\pm}$. The remaining contour $\Gamma_{p}:=\Gamma_{0}\cup\Gamma_{L}$ consists of the periodic planar end-sections $\Gamma_{0}$ and $\Gamma_{L}$, respectively situated at $x_{1}=0$ and $x_{1}=L$; the endpoints of $\Gamma_{L}$ are denoted by $\bm{z}^{\pm}$ (Fig. 1). Both channel walls are described as arcs $s\mapsto\bm{x}^{\pm}(s)$ with the arclength coordinate $s$ directed “leftwards” as depicted in Fig. 1, whereas the unit normal $\bm{n}$ to $\partial\Omega$ is everywhere taken as outwards to $\Omega$. The position vector $\bm{x}(s)$, unit tangent $\bm{\tau}(s)$, unit normal $\bm{n}(s)$ and curvature $\kappa(s)$ obey the Frenet formulas

[TABLE]

The opposite orientations of $(\bm{\tau},\bm{n})$ on $\Gamma^{+}$ and $\Gamma^{-}$ resulting from our choice of conventions imply opposite sign conventions on the curvature, which is everywhere on $\Gamma$ taken as $\kappa:=\bm{\tau}_{,s}\!\cdot\!\bm{n}$ for overall consistency.

In the wall frame, the wall material velocity must be tangent to $\Gamma$ (wall material points being constrained to remain on the surface $\Gamma$); moreover the wall material is assumed to be inextensible. In the wave frame, the wall material velocities $\bm{U}$ satisfying both requirements must have, on each wall, the form

[TABLE]

where $U$ is a constant. Moreover, in the wave frame, all wall material points travel over an entire spatial period during the time interval $T=L/c$, which implies $U=c\ell/L$. Finally, the viscous fluid must obey a no-slip condition on the wall, so that the velocity of fluid particles adjacent to $\bm{x}(s)$ is $\bm{U}(s)$. Concluding, the pumping motion of the wall constrains on each wall the fluid motion through

[TABLE]

Remark 1**.**

*The no-slip condition $\bm{u}=\bm{\tau}$ proposed in [29], which corresponds with the present notations to setting $U=c$, is inconsistent with the period $T=L/c$ of the traveling wave, unless the channel is straight ($\ell^{\pm}=L$). This discrepancy alters the boundary condition (3) and hence the shape optimization problem (the value of the power loss functional for given wall shape being affected in wall shape-dependent fashion). *

2.2 PDE formulation of forward problem

The stationary Stokes flow in the wave frame [29] solves the governing PDE formulation

[TABLE]

where $\bm{u}$, $p$ are the velocity and pressure, $\bm{D}[\bm{u}]:=(\bm{\nabla}\bm{u}+\bm{\nabla}^{\text{\scriptsize T}}\bm{u})/2$ is the strain rate tensor, and $\bm{\sigma}[\bm{u},p]$ is the stress tensor. Equations (4a,b,c) respectively express the local balance of momentum (absent any body forces), the constitutive relation ($\mu$ being the dynamic viscosity) and the incompressibility condition. The prescribed wall velocity $\bm{u}^{\text{\scriptsize D}}$ is given by (3). Problem (4) defines $p$ up to an arbitrary additive constant.

2.3 Weak formulation of forward problem

In this study, flow computations rely on a boundary integral equation formulation, see Section 4.3. It is however convenient, for later establishment of shape derivative identities, to recast the forward problem (4) as a mixed weak formulation (e.g. [4], Chap. 6):

[TABLE]

where \big{\langle}\hskip 1.00006pt\cdot,\cdot\hskip 1.00006pt\big{\rangle}_{\Gamma} stands for the $L^{2}(\Gamma)$ duality product, and the bilinear forms $a$ and $b$ are defined by

[TABLE]

The function spaces in problem (5) are as follows: $\mathcal{V}$ is the space of all periodic vector fields contained in $H^{1}(\Omega;\mathbb{R}^{2})$, $\mathcal{P}$ is the space of all $L^{2}(\Omega)$ functions with zero mean (i.e. obeying the constraint \big{\langle}\hskip 1.00006ptp,1\hskip 1.00006pt\big{\rangle}_{\Omega}=0) and \mbox{\boldmath\mathcal{F}}=H^{-1/2}(\Gamma;\mathbb{R}^{2}). The chosen definition of $\mathcal{P}$ caters for the fact that $p$ would otherwise be defined only up to arbitrary additive constants. The Dirichlet boundary condition (4d) is (weakly) enforced explicitly through (5c), rather than being embedded in the velocity solution space $\mathcal{V}$, as this will facilitate the derivation of shape derivative identities. The unknown $\bm{f}$, which acts as the Lagrange multiplier associated with the Dirichlet boundary condition (4d), is in fact the force density (i.e stress vector) $\bm{\sigma}[\bm{u},p]\!\cdot\!\bm{n}$.

2.4 Formulation of optimization problem, objective functional

We address the problem of finding the channel shape that makes peristaltic pumping most efficient, i.e., leads to minimum power loss for given mass flow rate. The power loss cost functional [29] is defined by

[TABLE]

It is a shape functional in that its value is completely determined by the domain $\Omega$ (in a partly implicit way through $\bm{f}_{\Omega}$). The minimization of $J_{\text{\scriptsize PL}}(\Omega)$, as in [29], subject to two constraints, involvesd two other shape functionals:

1. (a)

The channel volume per wavelength has a prescribed value $V_{0}$:

[TABLE] 2. (b)

The mass flow rate per wavelength $Q(\Omega)$ has a prescribed value $Q_{0}$. Expressing this constraint in the wave frame gives:

[TABLE]

Remark 2**.**

The power loss functional $J_{\text{\scriptsize PL}}(\Omega)$ is insensitive to modification of $p$ by an additive constant. This is seen, for instance, by combining (5a) with $\bm{v}=\bm{u}+c\bm{e}_{1}$ and (5b) with $q=p$, which yields

[TABLE]

3 Shape sensitivities

In this section, we review available shape derivative concepts that address our needs. Rigorous expositions of shape sensitivity theory can be found in [6, Chaps. 8,9] or [11, Chap. 5].

3.1 Shape sensitivity analysis: an overview

Let $\Omega_{\text{all}}\subset\mathbb{R}^{2}$ denote a fixed domain chosen so that $\Omega\Subset\Omega_{\text{all}}$ always holds for the shape optimization problem of interest. Shape changes are described with the help of transformation velocity fields, i.e. vector fields $\bm{\theta}:\Omega_{\text{all}}\to\mathbb{R}^{2}$ such that $\bm{\theta}=\mathbf{0}$ in a neighborhood $\partial\Omega_{\text{all}}$. Then, shape perturbations of domains $\Omega\Subset\Omega_{\text{all}}$ can be mathematically described using a pseudo-time $\eta$ and a geometrical transform of the form

[TABLE]

which defines a parametrized family of domains $\Omega_{\eta}(\bm{\theta}):=(\bm{I}+\eta\bm{\theta})(\Omega)$ for any given “initial” domain $\Omega\Subset\Omega_{\text{all}}$. The affine format (11) is sufficient for defining the first-order derivatives at $\eta=0$ used in this work.

Assumptions on shape transformations, admissible shapes

The set $\mathcal{O}$ of admissible shapes for the fluid region $\Omega$ in a channel period (in the wave frame) is defined as

[TABLE]

Accordingly, let the space $\Theta$ of admissible transformation velocities be defined as

[TABLE]

ensuring that the shape perturbations (i) are periodic, (ii) prevent any deformation of the end sections $\Gamma_{p}^{\pm}$ along the axial direction, and (iii) prevent vertical rigid translations of the channel domain. The provision $\bm{\theta}\in W_{0}^{1,\infty}(\Omega_{\text{all}})$ ensures that (a) there exists $\eta_{0}>0$ such that $\Omega_{\eta}(\bm{\theta})\Subset\Omega_{\text{all}}$ for any $\eta\in[0,\eta_{0}]$, and (b) the weak formulation for the shape derivative of the forward solution is well defined if the latter belongs to the standard solution space. Note that the volume and flow rate constraints (8), (9) are not embedded in $\Theta$; they will be accounted for via a Lagrangian.

Material derivatives

In what follows, all derivatives are implicitly taken at some given configuration $\Omega$, i.e. at initial ”time” $\eta=0$. The “initial” material derivative $\buildrel\scriptscriptstyle\star\over{\bm{a}}$ of some (scalar or tensor-valued) field variable $\bm{a}(\bm{x},\eta)$ is defined as

[TABLE]

and the material derivative of gradients and divergences of tensor fields are given by

[TABLE]

Likewise, the first-order “initial” derivative $J^{\prime}$ of a shape functional $J:\mathcal{O}\to\mathbb{R}$ is defined as

[TABLE]

Material differentiation of integrals

Consider, for given transformation velocity field $\bm{\theta}\in\Theta$, generic domain and contour integrals

[TABLE]

where $\Omega_{\eta}(\bm{\theta})=(\bm{I}+\eta\bm{\theta})(\Omega)$ is a variable domain and $S_{\eta}(\bm{\theta}):=(\bm{I}+\eta\bm{\theta})(S)$ a (possibly open) variable curve. The derivatives of $I_{\text{\scriptsize V}}(\eta)$ and $I_{\text{\scriptsize S}}(\eta)$ are given by the material differentiation identities

[TABLE]

which are well-known material differentiation formulas of continuum kinematics. In (18b), $\mbox{div}_{\!S}$ stands for the tangential divergence operator, which in the present 2D context is given, with $\kappa$ denoting the curvature (see Sec. A.1) by

[TABLE]

where we have set $\bm{\theta}$ on $\Gamma$, using the notations and conventions introduced in (1), in the form

[TABLE]

Shape functionals and structure theorem

The structure theorem for shape derivatives (see e.g. [6, Chap. 8, Sec. 3.3]) then states that the derivative of any shape functional $J$ is a linear functional in the normal transformation velocity $\theta_{n}=\bm{n}\!\cdot\!\bm{\theta}|_{\partial\Omega}$. For PDE-constrained shape optimization problems involving sufficiently smooth domains and data, the derivative \big{\langle}\hskip 1.00006ptJ^{\prime}(\Omega),\bm{\theta}\hskip 1.00006pt\big{\rangle} has the general form

[TABLE]

where $g$ is the shape gradient of $J$. The structure theorem intuitively means that (i) the shape of $\Omega_{\eta}$ is determined by that of $\partial\Omega_{\eta}$, and (ii) tangential components of $\bm{\theta}$ leave $\Omega_{\eta}$ unchanged at leading order $O(\eta)$.

Example: shape derivative of channel volume constraint.

Since $|\Omega|$ is given by (17a) with $F=1$, this purely geometric constraint is readily differentiated using (18a) and Green’s theorem, to obtain

[TABLE]

where the last equality results from provision (ii) in (13).

3.2 Shape derivative of forward solution

The power loss and the mass flow rate constraint are expressed in terms of functionals $J_{\text{\scriptsize PL}}$ and $C_{\text{\scriptsize Q}}$ that depend on $\Omega$ implicitly through the solution $(\bm{u},\bm{f},p)$ of the forward problem. Finding their shape derivatives will then involve the solution derivative $(\bm{u},\bm{f},p)^{\star}$. Setting up the governing problem for $(\bm{u},\bm{f},p)^{\star}$ is then a necessary prerequisite.

Proposition 1**.**

The governing weak formulation for the shape derivative (\makebox[5.72458pt]{\buildrel\scriptscriptstyle\star\over{\bm{u}}},\makebox[5.97226pt]{\buildrel\scriptscriptstyle\star\over{\bm{f}}},\makebox[5.03125pt]{\buildrel\scriptscriptstyle\star\over{p}}) of the solution $(\bm{u},\bm{f},p)$ of problem (5) is

[TABLE]

with the trilinear forms $a^{1}$ and $b^{1}$ given by

[TABLE]

Proof.

The proposition is obtained by applying the material differentiation identities (18) to the weak formulation (5), assuming that the test functions are such that \makebox[5.20601pt]{\buildrel\scriptscriptstyle\star\over{\bm{v}}}=\mathbf{0}, \makebox[5.12846pt]{\buildrel\scriptscriptstyle\star\over{\bm{g}}}=\mathbf{0} and \makebox[4.8229pt]{\buildrel\scriptscriptstyle\star\over{q}}=0, i.e. are convected under the shape perturbation. The latter provision is made possible by the absence of boundary constraints in the adopted definition of \text{\boldmath\mathcal{V}},\mbox{\boldmath\mathcal{F}},\mathcal{P} (Sec. 2.3). $\Box$

Moreover, the material derivative \makebox[5.72458pt]{\buildrel\scriptscriptstyle\star\over{\bm{u}}}{}^{\text{\scriptsize D}} of the Dirichlet data $\bm{u}^{\text{\scriptsize D}}$, involved in (23c), is given in the following lemma, proved in Sec. A.1:

Lemma 3.1**.**

Let $\bm{u}^{\text{\scriptsize D}}$ be the prescribed wall velocity (3). On either component of the wall $\Gamma$, we have

[TABLE]

Remark 3.2**.**

*The provision $\bm{\theta}\in W^{1,\infty}(\Omega_{\text{all}})$ in (13) serves to ensure well-definiteness of the trilinear forms $a^{1},b^{1}$ introduced in Proposition 1 for any (\bm{u},p)\in\text{\boldmath\mathcal{V}}\times\mathcal{P}. *

Remark 3.3**.**

*To maintain the zero-mean constraint on $p$ as $\Omega$ is perturbed, the mean of $\buildrel\scriptscriptstyle\star\over{p}$ must be fixed through \big{\langle}\hskip 1.00006pt\makebox[5.03125pt]{\buildrel\scriptscriptstyle\star\over{p}},1\hskip 1.00006pt\big{\rangle}_{\Omega}+\big{\langle}\hskip 1.00006ptp\mbox{div}\,\bm{\theta},1\hskip 1.00006pt\big{\rangle}_{\Omega}=0 (this provision plays no actual role in the sequel). *

We now state an identity involving the trilinear forms $a^{1},b^{1}$ that will play a crucial role in the derivation of convenient shape derivative formulas for the functionals involved in this work. The proof of this result is given in Sec. A.2

Lemma 3.4**.**

Let $(\bm{u},p)$ and $(\hat{\bm{u}}{},\hat{p})$ respectively satisfy $\mbox{div}\,\bm{u}=0$, $\mbox{div}\,(\bm{\sigma}[\bm{u},p])=\mathbf{0}$ and $\mbox{div}\,\hat{\bm{u}}{}=0$, $\mbox{div}\,(\bm{\sigma}[\hat{\bm{u}}{},\hat{p}])=\mathbf{0}$ in $\Omega$, with $\bm{\sigma}[\bm{u},p],\bm{\sigma}[\hat{\bm{u}}{},\hat{p}]$ given for both states by the constitutive relation (4b). Moreover, assume that $\bm{u}$, $\hat{\bm{u}}{}$ and $p$ are periodic, and set $\bm{f}:=\bm{\sigma}[\bm{u},p]\!\cdot\!\bm{n}$, $\hat{\bm{f}}{}:=\bm{\sigma}[\hat{\bm{u}}{},\hat{p}]\!\cdot\!\bm{n}$ and $\Delta\hat{p}(x_{2}):=\hat{p}(L,x_{2})-\hat{p}(0,x_{2})$ (i.e. periodicity is not assumed for $\hat{p}$). Then, the following identity holds:

[TABLE]

Proof 3.5**.**

*See Sec. A.2 *

3.3 Shape derivative of power loss functional

The derivative of the cost functional (7) is given, using the material differentiation identity (18b), by

[TABLE]

in terms of the derivatives \makebox[5.72458pt]{\buildrel\scriptscriptstyle\star\over{\bm{u}}}{}^{\text{\scriptsize D}} of the Dirichlet data and $\buildrel\scriptscriptstyle\star\over{\bm{f}}$ of the stress vector. The former is given by Lemma 3.1. The latter is part of the solution of the derivative problem (23), whose solution therefore seems necessary for evaluating \big{\langle}\hskip 1.00006ptJ^{\prime}_{\text{\scriptsize PL}}(\Omega),\bm{\theta}\hskip 1.00006pt\big{\rangle} in a given shape perturbation $\bm{\theta}$. However, a more effective approach allows to bypass the actual evaluation of $\buildrel\scriptscriptstyle\star\over{\bm{f}}$ by combining the forward and derivative problems with appropriate choices of test functions:

Lemma 3.6**.**

The following identity holds:

[TABLE]

Proof 3.7**.**

Consider the forward problem (5) with (\bm{v},\bm{g},q)=(\makebox[5.72458pt]{\buildrel\scriptscriptstyle\star\over{\bm{u}}},\makebox[5.97226pt]{\buildrel\scriptscriptstyle\star\over{\bm{f}}},\makebox[5.03125pt]{\buildrel\scriptscriptstyle\star\over{p}}) and the derivative problem (23) with $(\bm{v},\bm{g},q)=(\bm{u}+\bm{e}_{1},\bm{f},p)$. Making these substitutions and forming the combination $(\ref{der:weak}a)-(\ref{der:weak}b)-(\ref{der:weak}c)-(\ref{forward:weak}a)+(\ref{forward:weak}b)+(\ref{forward:weak}c)$, one obtains

[TABLE]

*The Lemma follows from using $\bm{u}=\bm{u}^{\text{\scriptsize D}}$ on $\Gamma$ and rearranging the above equality. *

Then, using Lemma 3.6 in (26) and evaluating $a^{1}(\bm{u},\bm{u},\bm{\theta})-2b^{1}(\bm{u},p,\bm{\theta})$ by means of Lemma 3.4 with $(\hat{\bm{u}}{},\hat{\bm{f}}{},\hat{p})=(\bm{u}+\bm{e}_{1},\bm{f},p)$ (implying that $\Delta p=0$), the derivative of $J_{\text{\scriptsize PL}}$ is recast in the following form, which no longer involves the solution of the derivative problem:

[TABLE]

This expression is still somewhat inconvenient as it involves (through $\bm{D}[\bm{u}]$) the complete velocity gradient on $\Gamma$. This can be alleviated by using the decomposition

[TABLE]

of the velocity gradient (where $\bm{\nabla}_{\!\!S}\bm{u}$ and $\partial_{n}\bm{u}$ respectively denote the tangential gradient and the normal derivative of $\bm{u}$) and expressing $\partial_{n}\bm{u}$ in terms of $\bm{f}$ by means of the constitutive relation (4b). In view of the specific form (3) of the Dirichlet data, the latter step is here conveniently carried out explicitly (in Sec. A.3), using curvilinear coordinates, and yields:

Lemma 3.8**.**

Let $(\bm{u},\bm{f},p)$ solve the forward problem (5). On the channel wall $\Gamma$, we have

[TABLE]

*where $f_{s}:=\bm{f}\!\cdot\!\bm{\tau}$ (hence the viscous part of $\bm{f}$ is tangential to $\Gamma$). *

Finally, we evaluate the density of the integral in (27) using the formulas of Lemmas 3.1 and 3.8. On performing traightforward algebra and rearranging terms, we obtain the following final result for \big{\langle}\hskip 1.00006ptJ^{\prime}_{\text{\scriptsize PL}}(\Omega),\bm{\theta}\hskip 1.00006pt\big{\rangle}, which is suitable for a direct implementation using the output of a boundary integral solver:

Proposition 2**.**

The shape derivative of the power loss cost functional $J_{\text{\scriptsize PL}}$ in a shape perturbation whose transformation velocity field $\bm{\theta}$ satisfies asumptions (13) is given (with $f_{s}:=\bm{f}\!\cdot\!\bm{\tau}$) by

[TABLE]

*where $\epsilon\in\{+,-\}$ refers to the upper wall $\Gamma^{+}$ or the lower wall $\Gamma^{-}$. *

Remark 3.9**.**

The shape derivative \big{\langle}\hskip 1.00006ptJ^{\prime}_{\text{\scriptsize PL}}(\Omega),\bm{\theta}\hskip 1.00006pt\big{\rangle} as given by Proposition 2 is a linear functional on $\theta_{n}|_{\Gamma}$ (since $\buildrel\scriptscriptstyle\star\over{\ell}$ is one, see Lemma 3.1), as predicted by the structure theorem. Moreover, it is insensitive to the pressure being possibly defined up to an additive constant: indeed, replacing $p$ with $p+\Delta p$ adds

[TABLE]

*to \big{\langle}\hskip 1.00006ptJ^{\prime}_{\text{\scriptsize PL}}(\Omega),\bm{\theta}\hskip 1.00006pt\big{\rangle} as given by Proposition 2, and this quantity vanishes due to the assumed periodicity of $\theta_{n}$. *

3.4 Shape derivative of mass flow rate functional

The derivative of the mass flow rate functional $C_{\text{\scriptsize Q}}$, defined by (9), is given (recalling (22) for the shape derivative of $|\Omega|$) by

[TABLE]

To evaluate \big{\langle}\hskip 1.00006ptI^{\prime}(\Omega),\bm{\theta}\hskip 1.00006pt\big{\rangle}, with $\bm{u}$ solving the forward problem (4), we recast $I(\Omega)$ as an integral over the end section $\Gamma_{0}$ by means of the divergence theorem, recalling that $\bm{u}_{\Omega}$ is divergence-free, tangential on $\Gamma$ and periodic:

[TABLE]

Consequently, using identity (18b) with provision (ii) of (13), we find:

[TABLE]

The forward solution derivative $\buildrel\scriptscriptstyle\star\over{\bm{u}}$ in the above expression will now be eliminated with the help of an adjoint problem. Let the adjoint state $(\hat{\bm{u}}{},\hat{\bm{f}}{},\hat{p})$ be defined as the solution of the mixed weak formulation

[TABLE]

The adjoint state $(\hat{\bm{u}}{},\hat{\bm{f}}{},\hat{p})$ thus results from applying a unit pressure difference $\Delta\hat{p}=1$ between the channel end sections while prescribing a no-slip condition on the channel walls, i.e. Problem (2.14–2.18) of [18] with homogeneous Dirichlet data on the walls (see Remark 3.11).

Then, combining the adjoint problem (34) with (\bm{v},\bm{g},q)=(\makebox[5.72458pt]{\buildrel\scriptscriptstyle\star\over{\bm{u}}},\makebox[5.97226pt]{\buildrel\scriptscriptstyle\star\over{\bm{f}}},\makebox[5.03125pt]{\buildrel\scriptscriptstyle\star\over{p}}) and the derivative problem (23) with $(\bm{v},\bm{g},q)=(\hat{\bm{u}}{},\hat{\bm{f}}{},\hat{p})$ and recalling the no-slip condition $\hat{\bm{u}}{}=\mathbf{0}$ on $\Gamma$ yields the identity

[TABLE]

We next use this identity in (33) and substitute the resulting expression of \big{\langle}\hskip 1.00006ptI^{\prime}(\Omega),\bm{\theta}\hskip 1.00006pt\big{\rangle} into (32), to obtain

[TABLE]

We first apply Lemma 3.4 with $\Delta\hat{p}=1$ to the above expression, to obtain

[TABLE]

whose form is similar to (27). The first integral in (37) can be reformulated with the help of Lemma 3.1, Lemma 3.8 and its following counterpart for the adjoint solution (whose proof is given in Sec. A.3):

Lemma 3.10**.**

Let $(\hat{\bm{u}}{},\hat{\bm{f}}{},\hat{p})$ solve the adjoint problem (34). On the channel wall $\Gamma$, we have

[TABLE]

*where $\hat{f}_{s}:=\hat{\bm{f}}{}\!\cdot\!\bm{\tau}$ (hence the viscous part of $\hat{\bm{f}}{}$ is tangential to $\Gamma$). *

Carrying out these derivations, and evaluating the second integral of (37), establishes the following final, implementation-ready, formula for \big{\langle}\hskip 1.00006ptC_{\text{\scriptsize Q}}^{\prime}(\Omega),\bm{\theta}\hskip 1.00006pt\big{\rangle}:

Proposition 3**.**

The shape derivative of the mass flow rate functional $C_{\text{\scriptsize Q}}$ in a domain shape perturbation such that the transformation velocity field $\bm{\theta}$ satisfies asumptions (13) is given by

[TABLE]

*where $f_{s}$ and $(\hat{p},\hat{f}_{s})$ are components of the solutions of the forward problem (5) and the adjoint problem (34), respectively, and $\epsilon\in\{+,-\}$ refers to the upper wall $\Gamma^{+}$ or the lower wall $\Gamma^{-}$. *

Remark 3.11**.**

Let the Stokes solution $(\hat{\bm{u}}{},\hat{p})$, periodic up to a constant (unit) pressure drop, be defined by

[TABLE]

*(i.e. Problem (2.14–2.18) of [18] with homogeneous Dirichlet data on the walls). Taking the dot product of the first equation by \bm{v}\in\text{\boldmath\mathcal{V}}, integrating by parts and using the periodicity of $\bm{v}$ and $2\mu\bm{D}[\bm{u}]$, we obtain equation (34a), while equations (34b,c) express the incompressibility and no-slip conditions. *

4 Numerical algorithm

We now formally define the shape optimization problem and describe a numerical algorithm to solve it based on the shape sensitivity formulas derived in the previous section.

4.1 Optimization method

Our goal is to find the shape of the peristaltic pump in the wave frame that minimizes the power loss functional subject to constant volume and flow rate constraints, that is,

[TABLE]

with $\mathcal{O}$ defined as in (12). While there are numerous approaches for solving constrained optimization problems [20], we use an Augmented Lagrangian (AL) approach and avoid second-order derivatives of the cost functional, whose evaluation is somewhat challenging in the case of our shape optimization problem. It proceeds by forming an augmented Lagrangian, defined by

[TABLE]

where $\sigma=(\sigma_{1},\sigma_{2})$ are penalty coefficients that are positive and $\lambda=(\lambda_{1},\lambda_{2})$ are Lagrange multipliers. Setting the initial values $\sigma^{0}$ and $\lambda^{0}$ using heuristics, the augmented Lagrangian method introduces a sequence $(m=1,2,\dots)$ of unconstrained minimization problems:

[TABLE]

with explicit Lagrange multiplier estimates $\lambda^{m}$ and increasing penalties $\sigma^{m}$. We use the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [20], a quasi-Newton method, for solving  (43). Equation (22), Propositions 2 and 3 are used in this context for all gradient evaluations.

The overall optimization procedure is summarized in Algorithm 1.

4.2 Finite-dimensional parametrization of shapes

In view of both the structure theorem, see (21), and the fact that we rely for the present study on a boundary integral method, we only need to model shape perturbations of the channel walls. Here we consider, for each wall $\Gamma^{\pm}$, parametrizations of the form $\bm{x}^{\pm}=\bm{x}^{\pm}(t,\bm{\xi})$ with

[TABLE]

where $\phi_{k}(t)$ are the trigonometric polynomials $\{\cos(t),\cos(2t),\dots,\cos(Nt),\sin(t),\sin(2t),\dots,\sin(Nt)\}$. The set of admissible shapes as in (12) indicates $x_{1}^{\pm}(0)=0$ and $x_{1}^{\pm}(2\pi)=L$, which are enforced in (44). The constraint (iii) of (13) is fulfilled by a fixed value $\xi_{2,0}^{-}=x_{2}^{-}(0)$ in (45), which will be pre-assigned and excluded from the shape parameters.

Therefore, we take as design shape parameters the set $\bm{\xi}=\{\xi^{\pm}_{1,1},\cdots,\xi^{\pm}_{1,2N},\xi^{+}_{2,0},\xi^{\pm}_{2,1},\cdots,\xi^{\pm}_{2,2N}\}$ of dimension $8N+1$. Since the parametrization (44), (45) is linear in $\bm{\xi}$, transformation velocities $\bm{\theta}$ associated to perturbed parameters $\bm{\xi}(\eta)=\bm{\xi}+\eta\bm{p}$, i.e., to the $\eta$-dependent parametrization have the form

[TABLE]

4.3 Boundary integral solver

Solving the unconstrained optimization problem (43) requires evaluating the shape sensitivites (Propositions 2 and 3), which in turn require solving the forward and adjoint problems to obtain the corresponding traction and pressure on the channel walls. In both problems, the fluid velocity and pressure satisfy the Stokes equations:

[TABLE]

While the forward problem requires applying prescribed slip on the walls and periodic boundary conditions, the adjoint problem requires a no-slip on the walls and unit pressure drop across the channel. The boundary conditions can be summarized in a slightly modified form as follows.

[TABLE]

Here, $T$ is the traction vector, whose components are given by $T_{i}(\bm{u},p)=\sigma_{ij}(\bm{u},p)\bm{n}_{j}$ and it represents the hydrodynamic force experienced by any interface in the fluid with normal $\bm{n}$. These systems of equations correspond to (4) and (34) respectively by unique continuation of Cauchy data [18]. The standard boundary integral approach for periodic flows is to use periodic Greens functions obtained by summing over all the periodic copies of the sources [21]. The disadvantage of this approach is the slow convergence rate, specially, for high-aspect ratio domains; moreover, the pressure-drop condition cannot be applied directly. Instead, we use the periodization scheme developed recently in [18] that uses the free-space kernels only and enforces the inlet and outlet flow conditions in (48, 49) algebraically at a set of collocation nodes.

The free-space Stokes single-layer kernel and the associated pressure kernel, given a source point ${\bf y}$ and a target point ${\bf x}$, are given by

[TABLE]

The approach of [18] represents the velocity field as a sum of free-space potentials defined on the unit cell $\Omega$, its nearest periodic copies and at a small number, $K$, of auxiliary sources located exterior to $\Omega$ that act as proxies for the infinite number of far-field periodic copies:

[TABLE]

where

[TABLE]

Here, ${\bf d}$ is the lattice vector i.e. ${\bf d}=L{\bf e}_{1}$ and the source locations $\{{\bf y}_{m}\}_{m=1}^{K}$ are chosen to be equispaced on a circle enclosing $\Omega$. The associated representation for pressure is given by

[TABLE]

where

[TABLE]

The pair $(\bm{u},p)$ as defined by this representation satisfy the Stokes equations since $S$ and $Q$, defined in (50), are the Green’s functions. The unknown density function $\bm{\tau}$ and the coefficients $\{{\bf c}_{m}\}$ are then determined by enforcing the boundary conditions. For the forward problem, applying the conditions in (48) produces a system of equations in the following form [18]:

[TABLE]

The first row applies the slip condition on $\Gamma$ by taking the limiting value of $\bm{u}(\bm{x})$, defined in (51), as ${\bf x}$ approaches $\Gamma$ from the interior. The second row applies the periodic boundary conditions on velocity and traction as defined in (48). The operators $A,B,C$ and $D$ are correspondingly defined based on the representation formulas (51) and (53). In the case of the adjoint problem, the operators remain the same but the right hand side of (55) is modified according to the boundary conditions (49).

Remark 4.1**.**

*Note that the representation (51) implies that $A$ is a first-kind boundary integral operator and in general is not advisable for large-scale problems due to ill-conditioning of resulting discrete linear systems. However, for the problems we consider here, the dimension of linear system is usually small $\sim$100-200 and we always solve it using direct solvers. For large-scale problems requiring iterative solution (e.g., in peristaltic pumps transporting rigid or deformable particles), well-conditioned systems can be produced from second-kind operators, which can readily be constructed using double-layer potentials as was done in [18]. *

Finally, we use $M$ quadrature nodes each on $\Gamma^{+}$ and $\Gamma^{-}$ to evaluate smooth integrals using the standard periodic trapezoidal rule and weakly singular integrals, such as (52) evaluated on $\Gamma$, using a spectrally-accurate Nyström method (with periodic Kress corrections for the log singularity, see Sec. 12.3 of [16]). Thereby, a discrete linear system equivalent to (55) is obtained, which is solved using a direct solver for the unknowns $\bm{\tau}$ and $\{{\bf c}_{m}\}$. The traction vector for a given pipe shape, required for evaluating the shape derivatives (Props. 2 and 3), can then be obtained by a similar representation as (51) but with the kernel replaced by the traction kernel, given by,

[TABLE]

5 Results

In this section, we present validation results for our numerical PDE solver and test the performance of our shape optimization algorithm.

First, we demonstrate the performance of our forward problem (4) solver on an arbitrary pump shape as shown in Figure 2(a). To illustrate the convergence of the numerical scheme, we consider two scalar quatities: the objective function $J_{PL}$ and the mass flow rate $Q$, both of which depend on the traction vector obtained by solving the forward problem. In 2(b), we show the accuracy of the solver in computing these quantities as well as the CPU time it takes to solve as the number of quadrature points on each of the walls, $M$, is increased. Since we are using a spectrally-accurate quadrature rule, notice that the error decays rapidly and a small number of points are sufficient to achieve six-digit accuracy, which is more than enough for our application. The number of proxy points and the number of collocation points on the side walls are chosen following the analysis of [18]—both typically are small again owing to the spectral convergence of the method w.r.t these parameters. Consequently, each forward solve takes less than a second to obtain the solution to six-digit accuracy as can be observed from Figure 2(b).

Next we consider two different initial pump shapes and perform the shape optimization by imposing the same constraints on the flow rate $Q^{*}$ and the volume $V^{*}$. The initial and the intermediary shapes produced by our numerical optimization procedure are shown in Figure 3. Based on our earlier analysis of the forward solver, we set $M=64$ and the number of modes $N$ (e.g., in (44) and (45)) is set to 5, thereby, the number of the shape design parameters is 41. We use the Augmented Lagrangian approach described in Algorithm 1 with the values of the penalty parameters as listed. Each AL iteration requires solving an unconstrained optimization problem, which entails taking several BFGS iterations (steps 5-9 of Algo 1). Both the number of AL iterations and the number of BFGS iterations are shown for each shape update in Figure 3. In total, it costs 143 solves of the forward and adjoint problems (48)–(49) for the case in Fig. 3(a) and 197 for the case in Fig. 3(b).

We show the evolution of the objective function and the constraints with the iteration index, corresponding to these two test cases, in Fig. 4. Due to the fact that we are using an AL approach, the constraints are enforced progressively by increasing the penalty coefficients $\sigma$ and as a result, the objective function approaches a local minimum in a non-monotonic fashion. In the second test case (Fig. 3b), for instance, $J_{PL}$ increases significantly since the initial shape doesn’t satisfy the constraints. On the other hand, in the first test case, the initial values for the $Q$ and $V$, obtained by a forward solve, match the target values $Q^{*}$ and $V^{*}$, thereby $J_{PL}$ is reduced as the optimization proceeds, to nearly half of its initial value. While not guaranteed in general (due to the possibility of getting stuck in a local minimum), the final shapes in both cases coincide, which is another validation of our analytic shape sensitivity calculations.

The equilibrium shapes and their interior fluid flow shown in Figure 3 reaffirm the classical observation [12] of trapping i.e., an enclosed bolus of fluid particles near the center line indicated by the closed streamlines in the waveframe. As is also well-known, trapping occurs beyond a certain pumping range only; in Fig. 5, we show the optimal shapes obtained by our algorithm at different flow rates but containing the same volume of fluid. Here, we fixed the bottom wall to be flat. Notice that the bolus appears to form for shapes beyond $Q=0.9$. Moreover, as expected, the optimal value of power loss is higher for higher flow rates with the extreme case of a flat pipe transporting zero net flow with no power loss.

6 Conclusions

We derived new analytic formulas for evaluating the shape derivatives of the power loss and the mass flow rate functionals that arise in the shape optimization of Stokesian peristaltic pumps. While we restricted our attention to two-dimensional shapes, extension of these formulas to rotationally-symmetric shapes is rather straightforward. We applied the recently developed periodic boundary integral solver of [18] to solve the forward/adjoint PDE problems efficiently.

One of the main directions of our future investigation is to consider optimal shapes for transporting passive (e.g., colloids, bubbles or vesicles) and/or active (e.g., spermatozoa, bacteria) particles in Stokes flow. This is of current interest in both science and technological applications. Clearly, they are more computationally intensive since transient PDEs need to be solved as opposed to the quasi-static ones considered here. However, particulate flow solvers based on BIEs are proven to be efficient and scalable to simulating millions of deformable particles [22]; our work lays the foundation for applying these methods to shape optimization of peristaltic pumps transporting such complex fluids.

7 Acknowledgements

RL and SV acknowledge support from NSF under grants DMS-1719834 and DMS-1454010. The work of SV was also supported by the Flatiron Institute (USA), a division of Simons Foundation, and by the Fondation Mathématique Jacques Hadamard (France).

Appendix A Proofs

A.1 Proof of Lemma 3.1

We use the Frenet formulas (1) and associated conventions. To evaluate \makebox[5.72458pt]{\buildrel\scriptscriptstyle\star\over{\bm{u}}}{}^{\text{\scriptsize D}}, we let $\Gamma$ depend on the fictitious time $\eta$, setting

[TABLE]

(where $\Gamma$ stands for $\Gamma^{+}$ or $\Gamma^{-}$, and likewise for $\ell$) and seek the relevant derivatives w.r.t. $\eta$ evaluated at $\eta=0$. Note that for $\eta\not=0$, $s$ is no longer the arclength coordinate along $\Gamma_{\eta}$, and $\partial_{s}\bm{x}_{\eta}$ is no longer of unit norm; moreover, the length of $\Gamma_{\eta}$ depends on $\eta$. The wall velocity $\bm{U}=(c\ell/L)\bm{\tau}$ for varying $\eta$ is then given by

[TABLE]

having set $g_{\eta}=|\partial_{s}\bm{x}_{\eta}|$ (note that $g_{0}=1$). Our task is to evaluate $\text{d}/\text{d}\eta\,\bm{U}_{\eta}(s)$ at $\eta=0$. We begin by observing that the derivative of $g$ is (since $\partial_{s}\bm{x}_{\eta}=\bm{\tau}$ and $g=1$ for $\eta=0$)

[TABLE]

and the length $\ell_{\eta}$ of $\Gamma_{\eta}$ and its derivative $\buildrel\scriptscriptstyle\star\over{\ell}$ are given (noting that $s$ spans the fixed interval $[0,\ell]$ for all curves $\Gamma_{\eta}$) by

[TABLE]

The last equality in (ii), which results from the assumed periodicity of $\bm{\theta}$, proves item (b) of the lemma.

Using these identities in (58), the sought derivative \makebox[5.72458pt]{\buildrel\scriptscriptstyle\star\over{\bm{u}}}{}^{\text{\scriptsize D}} of the wall velocity is found as

[TABLE]

thus establishing item (a) of the lemma. The proof of Lemma 3.1 is complete.

A.2 Proof of Lemma 3.4

The proof proceeds by verification, and rests on evaluating $\mbox{div}\,\bm{A}$, with the vector function $\bm{A}$ defined by

[TABLE]

First, it is easy to check, e.g. using component notation relative to a Cartesian frame, that

[TABLE]

Next, invoking the constitutive relation (2.2b) for both states, one has

[TABLE]

i.e. (using incompressibility)

[TABLE]

Finally, using (64) in (61a) and the balance equation in (61b), together with the corresponding identities obtained by switching $(\bm{u},p)$ and $(\hat{\bm{u}}{},\hat{p})$, one obtains

[TABLE]

(with the second equality stemming from $\bm{\sigma}[\bm{u},p]\!:\!\bm{D}[\hat{\bm{u}}{}]=2\mu\bm{D}[\bm{u}]\!:\!\bm{D}[\hat{\bm{u}}{}]=2\mu\bm{D}[\hat{\bm{u}}{}]\!:\!\bm{D}[\bm{u}]=\bm{\sigma}[\hat{\bm{u}}{},\hat{p}]\!:\!\bm{D}[\bm{u}]$). Using definitions (24a,b) of $a^{1}$ and $b^{1}$, we therefore observe that

[TABLE]

The last step consists of applying the first Green identity (divergence theorem) to the above integral. The Lemma follows, with the contribution of the end section $\Gamma_{L}$ therein stemming from condition (ii) in (13) and the periodicity conditions at the end sections. The latter hold by assumption for both $\bm{u}$ and $\hat{\bm{u}}{}$, and the interior regularity of solutions in the whole channel then implies the same periodicity for $\bm{\nabla}\bm{u}$ and $\bm{\nabla}\hat{\bm{u}}{}$; moreover, periodicity is also assumed for $p$ (but not necessarily for $\hat{p}$) as well as for $\bm{\theta}$.

A.3 Proof of Lemmas 3.8 and 3.10

Let points $\bm{x}$ in a tubular neighborhood $V$ of $\Gamma$ be given in terms of curvilinear coordinates $(s,z)$, so that

[TABLE]

and let $\bm{v}(\bm{x})=v_{s}(s,z)\bm{\tau}(s)+v_{n}(s,z)\bm{n}(s)$ denote a generic vector field in $V$. Then, at any point $\bm{x}=\bm{x}(s)$ of $\Gamma$, we have

[TABLE]

Assuming incompressibility, the condition $\mbox{div}\,\bm{v}=0$ can be used for eliminating $\partial_{n}v_{n}$ and we obtain

[TABLE]

Recalling now that the forward and adjoint solutions respectively satisfy $\bm{u}=(c\ell/L)\bm{\tau}$ and $\hat{\bm{u}}{}=\mathbf{0}$ on $\Gamma$, and that $2\bm{D}[\bm{v}]=\bm{\nabla}\bm{v}+\bm{\nabla}\bm{v}^{\text{\scriptsize T}}$, we have

[TABLE]

The corresponding stress vectors $\bm{f}=-p\bm{n}+2\mu\bm{D}[\bm{u}]\!\cdot\!\bm{n}$ and $\hat{\bm{f}}{}=-\hat{p}\bm{n}+2\mu\bm{D}[\hat{\bm{u}}{}]\!\cdot\!\bm{n}$ on $\Gamma$ are found as

[TABLE]

which in particular prove items (c) of Lemmas 3.8 and 3.10. Finally, using the above in (73) establishes the remaining items (a), (b) of both lemmas.