2.1 Dimensional formulation

2.1.1 Flow equations

We refer to the fluid/flow exterior to the tissue construct as the bulk fluid/flow, and the fluid/flow within the tissue construct as the interior fluid/flow or porous fluid/flow. In the rotating frame, the bulk flow satisfies the incompressible Navier–Stokes equations as follows:

(2.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}(\dot{\boldsymbol{u}}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}+2\dot{\unicode[STIX]{x1D711}}\boldsymbol{e}_{z}\times \boldsymbol{u}+\ddot{\unicode[STIX]{x1D711}}\boldsymbol{e}_{z}\times \boldsymbol{X}-\dot{\unicode[STIX]{x1D711}}^{2}\boldsymbol{X}^{\bot })=-\unicode[STIX]{x1D735}\tilde{p}+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\unicode[STIX]{x1D70C}\boldsymbol{g}, & \displaystyle\end{eqnarray}$$

(2.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}=u\boldsymbol{e}_{X}+v\boldsymbol{e}_{Y}+w\boldsymbol{e}_{z}$ is the bulk fluid velocity in the rotating frame, $\tilde{p}$ is the bulk fluid pressure, $\boldsymbol{X}^{\bot }=X\boldsymbol{e}_{X}+Y\boldsymbol{e}_{Y}$ is the projection of $\boldsymbol{X}$ onto the $(X,Y)$ plane, and an overdot represents differentiation with respect to time. The interior flow satisfies the incompressible Darcy equations

(2.3a ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{Q}-\dot{R}\boldsymbol{e}_{X}+(\dot{\unicode[STIX]{x1D711}}-\unicode[STIX]{x1D714})\boldsymbol{e}_{z}\times \boldsymbol{r}=-{\displaystyle \frac{k}{\unicode[STIX]{x1D707}}}(\unicode[STIX]{x1D735}\tilde{P}-\unicode[STIX]{x1D70C}\boldsymbol{g}), & \displaystyle\end{eqnarray}$$

(2.3b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{Q}=0, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{Q}=U\boldsymbol{e}_{X}+V\boldsymbol{e}_{Y}+W\boldsymbol{e}_{z}$ is the Darcy flow velocity in the rotating frame, $\tilde{P}$ is the interior fluid pressure, and $k$ is the permeability of the tissue construct. The transformation of the bulk flow and Darcy velocity from the rotating to the laboratory frame is given by

(2.4a,b ) $$\begin{eqnarray}\boldsymbol{u}_{LF}=\boldsymbol{u}+\dot{\unicode[STIX]{x1D711}}\boldsymbol{e}_{z}\times \boldsymbol{X},\quad \boldsymbol{Q}_{LF}=\boldsymbol{Q}+\dot{\unicode[STIX]{x1D711}}\boldsymbol{e}_{z}\times \boldsymbol{X},\end{eqnarray}$$

where $\boldsymbol{u}_{LF}$ and $\boldsymbol{Q}_{LF}$ are the bulk flow and Darcy velocity in the laboratory frame, respectively. Thus, the left-hand side of (2.3a ) is equivalent to $\boldsymbol{Q}_{LF}-\boldsymbol{V}_{TC}$ , the difference between the interior flow velocity in the laboratory frame and the velocity of the construct.

We denote the entire bioreactor surface as $S$ , and the interface between the tissue construct and the bulk fluid as $T$ . On the bioreactor surface, we impose a no-slip condition which, in the rotating frame, is given by

(2.5) $$\begin{eqnarray}\boldsymbol{u}=(\unicode[STIX]{x1D6FA}-\dot{\unicode[STIX]{x1D711}})\boldsymbol{e}_{z}\times \boldsymbol{X}\quad \text{on }S.\end{eqnarray}$$

On the tissue construct interface, we impose continuity of normal flux, continuity of pressure, and a no-slip condition (Levy & Sanchez-Palencia Reference Levy and Sanchez-Palencia1975; Carraro et al. Reference Carraro, Goll, Marciniak-Czochra and Mikelić2015) as follows

(2.6a-c ) $$\begin{eqnarray}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}=\boldsymbol{Q}\boldsymbol{\cdot }\boldsymbol{n},\quad \tilde{p}=\tilde{P},\quad \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{s}_{i}=(\dot{R}\boldsymbol{e}_{X}+(\unicode[STIX]{x1D714}-\dot{\unicode[STIX]{x1D711}})\boldsymbol{e}_{z}\times \boldsymbol{r})\boldsymbol{\cdot }\boldsymbol{s}_{i}\quad \text{on }T,\end{eqnarray}$$

where $\boldsymbol{n}$ is the unit normal pointing out of the tissue construct and $\boldsymbol{s}_{i}$ is any unit tangent vector to the surface.

2.1.2 Construct motion equations

In Dalwadi (Reference Dalwadi2014), the change in the linear and angular momentum of the construct due to the motion of fluid within and across the interface of the construct is found to be negligible when Darcy’s law is applicable, so we disregard these effects in this paper. Then, in the inertial frame, Newton’s second law for linear momentum yields

(2.7) $$\begin{eqnarray}\iint _{T}\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S-M_{c}g(\sin \unicode[STIX]{x1D711}(t)\boldsymbol{e}_{X}+\cos \unicode[STIX]{x1D711}(t)\boldsymbol{e}_{Y})=M_{c}\frac{\text{d}^{2}}{\text{d}t^{2}}(R(t)\boldsymbol{e}_{X}),\end{eqnarray}$$

where $\unicode[STIX]{x1D748}$ is the Newtonian stress tensor of the bulk fluid. Here, $M_{c}$ is the mass of the wetted construct, equal to the product of the effective construct density $\unicode[STIX]{x1D70C}_{c}$ and the construct volume

(2.8) $$\begin{eqnarray}M_{c}=\unicode[STIX]{x1D70C}_{c}\unicode[STIX]{x03C0}l^{2}d.\end{eqnarray}$$

In terms of experimentally measurable parameters, $\unicode[STIX]{x1D70C}_{c}$ is defined to be

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{c}=\unicode[STIX]{x1D70C}\unicode[STIX]{x1D719}+\unicode[STIX]{x1D70C}_{s}(1-\unicode[STIX]{x1D719}),\end{eqnarray}$$

where $\unicode[STIX]{x1D719}\in (0,1]$ is the porosity of the construct and $\unicode[STIX]{x1D70C}_{s}$ is the density of the solid matrix within the construct. We note that $\unicode[STIX]{x1D719}=0$ for a solid construct and $\unicode[STIX]{x1D719}=1$ for a fully fluid construct, but that the coupling conditions we derive in appendix A at first order are not valid for the lower limit. The first term on the left-hand side of (2.7) is the hydrodynamical force on the construct surface (and also contains the buoyancy forces), the second term on the left-hand side is the weight of the wet mass of the construct, and the right-hand side is the rate of change of linear momentum.

Similarly, Newton’s second law for angular momentum implies that

(2.10) $$\begin{eqnarray}\iint _{T}\boldsymbol{r}\times (\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n})\,\text{d}S={\displaystyle \frac{M_{c}l^{2}}{2}}\frac{\text{d}\unicode[STIX]{x1D714}}{\text{d}t}\boldsymbol{e}_{z},\end{eqnarray}$$

where the left-hand side is the hydrodynamical torque on the construct surface, and the right-hand side is the rate of change of angular momentum.

2.2 Dimensionless equations

We non-dimensionalize by setting

(2.11) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}(X,Y,z)=a(X^{\ast },Y^{\ast },\unicode[STIX]{x1D716}z^{\ast }),\quad (l,R)=a(b^{\ast },R^{\ast }),\\ t=\unicode[STIX]{x1D6FA}^{-1}t^{\ast },\quad \unicode[STIX]{x1D714}=\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D714}^{\ast },\\ (u,v,w)=a\unicode[STIX]{x1D6FA}(u^{\ast },v^{\ast },\unicode[STIX]{x1D716}w^{\ast }),\quad (U,V,W)=a\unicode[STIX]{x1D6FA}(U^{\ast },V^{\ast },\unicode[STIX]{x1D716}W^{\ast }),\\ (\tilde{p},\tilde{P})=\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D716}^{-2}(\tilde{p}^{\ast },\tilde{P}^{\ast }),\quad \unicode[STIX]{x1D748}=\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D716}^{-2}\unicode[STIX]{x1D748}^{\ast },\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x1D716}=h/a$ . In (2.11), asterisks denote dimensionless variables, and we henceforth drop asterisks for convenience. In dimensionless variables, we continue to use the vector notations $\boldsymbol{X}=X\boldsymbol{e}_{X}+Y\boldsymbol{e}_{Y}+z\boldsymbol{e}_{z}$ , $\boldsymbol{r}=\boldsymbol{X}-R\boldsymbol{e}_{X}$ , $\boldsymbol{u}=u\boldsymbol{e}_{X}+v\boldsymbol{e}_{Y}+w\boldsymbol{e}_{z}$ and $\boldsymbol{Q}=U\boldsymbol{e}_{X}+V\boldsymbol{e}_{Y}+W\boldsymbol{e}_{z}$ . The bulk flow equations (2.2) become

(2.12a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}Re(\dot{u}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})u-2\dot{\unicode[STIX]{x1D711}}v-\ddot{\unicode[STIX]{x1D711}}Y-\dot{\unicode[STIX]{x1D711}}^{2}X)=-p_{X}+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6FB}_{\bot }^{2}u+u_{zz}, & \displaystyle\end{eqnarray}$$

(2.12b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}Re(\dot{v}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})v+2\dot{\unicode[STIX]{x1D711}}u+\ddot{\unicode[STIX]{x1D711}}X-\dot{\unicode[STIX]{x1D711}}^{2}Y)=-p_{Y}+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6FB}_{\bot }^{2}v+v_{zz}, & \displaystyle\end{eqnarray}$$

(2.12c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}^{3}Re({\dot{w}}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})w)=-p_{z}+\unicode[STIX]{x1D716}^{4}\unicode[STIX]{x1D6FB}_{\bot }^{2}w+\unicode[STIX]{x1D716}^{2}w_{zz}, & \displaystyle\end{eqnarray}$$

(2.12d ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D6FB}_{\bot }^{2}=\unicode[STIX]{x2202}_{XX}+\unicode[STIX]{x2202}_{YY}$

$Re=\unicode[STIX]{x1D70C}h(a\unicode[STIX]{x1D6FA})/\unicode[STIX]{x1D707}$

$p=\tilde{p}+B(X\sin \unicode[STIX]{x1D711}+Y\cos \unicode[STIX]{x1D711})$

$B=\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D70C}ag/(\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FA})$

The interior flow equations (2.3) become

(2.13a ) $$\begin{eqnarray}\displaystyle & \displaystyle U-\dot{R}-(\dot{\unicode[STIX]{x1D711}}-\unicode[STIX]{x1D714})Y=-KP_{X}, & \displaystyle\end{eqnarray}$$

(2.13b ) $$\begin{eqnarray}\displaystyle & \displaystyle V+(\dot{\unicode[STIX]{x1D711}}-\unicode[STIX]{x1D714})(X-R)=-KP_{Y}, & \displaystyle\end{eqnarray}$$

(2.13c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}^{2}W=-KP_{z}, & \displaystyle\end{eqnarray}$$

(2.13d ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{Q}=0, & \displaystyle\end{eqnarray}$$

$K=k/(\unicode[STIX]{x1D716}a)^{2}$

$P=\tilde{P}+B(X\sin \unicode[STIX]{x1D711}+Y\cos \unicode[STIX]{x1D711})$

On the bioreactor boundary $S$ , the no-slip boundary condition (2.5) becomes

(2.14) $$\begin{eqnarray}\boldsymbol{u}=(1-\dot{\unicode[STIX]{x1D711}})\boldsymbol{e}_{z}\times \boldsymbol{X}\quad \text{on }S.\end{eqnarray}$$

On the tissue construct interface, the boundary conditions (2.6) become

(2.15a-c ) $$\begin{eqnarray}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}=\boldsymbol{Q}\boldsymbol{\cdot }\boldsymbol{n},\quad p=P,\quad \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{s}_{i}=(\dot{R}\boldsymbol{e}_{X}+(\unicode[STIX]{x1D714}-\dot{\unicode[STIX]{x1D711}})\boldsymbol{e}_{z}\times \boldsymbol{r})\boldsymbol{\cdot }\boldsymbol{s}_{i}\quad \text{on }T.\end{eqnarray}$$

Finally, we give the dimensionless form of the construct motion (2.7) and (2.10), which become

(2.16a ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{f}_{c}+\boldsymbol{f}_{f}+\boldsymbol{f}_{g}=\unicode[STIX]{x1D716}^{2}Re\bar{\unicode[STIX]{x1D70C}}{\mathcal{V}}\frac{\text{d}^{2}}{\text{d}t^{2}}(R(t)\boldsymbol{e}_{X}), & \displaystyle\end{eqnarray}$$

(2.16b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D749}_{c}+\unicode[STIX]{x1D749}_{f}=\unicode[STIX]{x1D716}^{2}Re\bar{\unicode[STIX]{x1D70C}}{\mathcal{V}}{\displaystyle \frac{b^{2}}{2}}\frac{\text{d}\unicode[STIX]{x1D714}}{\text{d}t}\boldsymbol{e}_{z}, & \displaystyle\end{eqnarray}$$

where $\bar{\unicode[STIX]{x1D70C}}=\unicode[STIX]{x1D70C}_{c}/\unicode[STIX]{x1D70C}$ is the ratio between the effective density of the construct and the density of the nutrient fluid, ${\mathcal{V}}=\unicode[STIX]{x03C0}b^{2}(1-2\unicode[STIX]{x1D716}_{1})$ is the dimensionless volume of the construct, where $\unicode[STIX]{x1D716}_{1}=(h-d)/(2h)$ is the ratio of the gap distance between the flat sides of the tissue construct and bioreactor, and the thickness of the bioreactor, and we decompose the force and torque as follows:

(2.16c ) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{f}_{c}=\unicode[STIX]{x1D716}\iint _{T_{1}}\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S,\quad \unicode[STIX]{x1D749}_{c}=\unicode[STIX]{x1D716}\iint _{T_{1}}\boldsymbol{r}\times (\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n})\,\text{d}S,\quad \boldsymbol{f}_{f}=\iint _{T_{2}\cup T_{3}}\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S,\\ \displaystyle \unicode[STIX]{x1D749}_{f}=\iint _{T_{2}\cup T_{3}}\boldsymbol{r}\times (\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n})\,\text{d}S,\quad \boldsymbol{f}_{g}=-\unicode[STIX]{x1D716}B\bar{\unicode[STIX]{x1D70C}}{\mathcal{V}}(\sin \unicode[STIX]{x1D711}(t)\boldsymbol{e}_{X}+\cos \unicode[STIX]{x1D711}(t)\boldsymbol{e}_{Y}).\end{array}\right\}\end{eqnarray}$$

Each integral is defined over part of the tissue construct interface, decomposed as follows $T=T_{1}\cup T_{2}\cup T_{3}$ , where the component parts are defined to be

(2.17a ) $$\begin{eqnarray}\displaystyle & \displaystyle T_{1}=\{(X,Y,z)\in \mathbb{R}^{3}:(X-R)^{2}+Y^{2}=b^{2},\unicode[STIX]{x1D716}_{1}\leqslant z\leqslant 1-\unicode[STIX]{x1D716}_{1}\}, & \displaystyle\end{eqnarray}$$

(2.17b ) $$\begin{eqnarray}\displaystyle & \displaystyle T_{2}=\{(X,Y,z)\in \mathbb{R}^{3}:(X-R)^{2}+Y^{2}<b^{2},z=\unicode[STIX]{x1D716}_{1}\}, & \displaystyle\end{eqnarray}$$

(2.17c ) $$\begin{eqnarray}\displaystyle & \displaystyle T_{3}=\{(X,Y,z)\in \mathbb{R}^{3}:(X-R)^{2}+Y^{2}<b^{2},z=1-\unicode[STIX]{x1D716}_{1}\}. & \displaystyle\end{eqnarray}$$

$T_{1}$

$T_{2}$

$T_{3}$

$\boldsymbol{f}_{c}$

$\unicode[STIX]{x1D749}_{c}$

$\boldsymbol{f}_{f}$

$\unicode[STIX]{x1D749}_{f}$

$\boldsymbol{f}_{g}$

We have several dimensionless parameters in our model and we note that these parameter values can vary greatly in magnitude (table 1). We now make some further assumptions on the size of these parameters relevant to the HARV bioreactor. We first assume that the bioreactor thickness is much smaller than its radius, $\unicode[STIX]{x1D716}\ll 1$ , which is an intrinsic characteristic of the HARV bioreactor. We further assume that the flow is dominated by pressure and viscosity rather than inertial effects, so that $\unicode[STIX]{x1D716}Re\ll 1$ . However, as we are interested in the effect of a weak inertia on this system, we additionally assume that $Re\gg 1$ . We also make some size assumptions on the ratio between the width of the gap between the flat sides of the tissue construct and bioreactor, and the bioreactor thickness, defined as $\unicode[STIX]{x1D716}_{1}=(h-d)/(2h)$ . We first assume that this gap is much smaller than the bioreactor thickness, so that $\unicode[STIX]{x1D716}_{1}\ll 1$ . To focus on the effect of weak inertia, we additionally assume that the effect of the thin gap between the tissue construct and bioreactor on the flow is much smaller than the effect of inertia, and we take $\unicode[STIX]{x1D716}_{1}\ll \unicode[STIX]{x1D716}Re$ . This assumption is strictly stronger than the previous one, and allows us to neglect the role of the thin gap in the reduced flow problem we define below. This assumption will become relevant in appendix A when we derive the effective boundary conditions for the reduced flow problem. Although the gap has a small effect on the flow in the bioreactor, the gap remains important for the trajectory problem because it causes a significant friction between the construct and the bioreactor. Finally, we further assume that the thin gap is smaller than the boundary layers determined in Dalwadi et al. (Reference Dalwadi, Chapman, Waters and Oliver2016) and appendix A, requiring $\unicode[STIX]{x1D716}_{1}Re\ll 1$ ; thus, the gap does not affect the boundary layer structure and subsequent coupling results.

Table 1. Typical dimensional (left) and dimensionless (right) parameter values. The absolute upper bound of $1/12\approx 8.3\times 10^{-2}$ for $K$ can be obtained by considering unobstructed Poiseuille flow through a channel. Data for the HARV operating conditions and construct dimensions are taken from Wolf, Sams & Schwarz (Reference Wolf, Sams and Schwarz1992), Freed & Vunjak-Novakovic (Reference Freed and Vunjak-Novakovic1997), Ingram et al. (Reference Ingram, Techy, Saroufeem, Yazan, Narayan, Goodwin and Spaulding1997), Schwarz & Anderson (Reference Schwarz and Anderson1998), Yu et al. (Reference Yu, Botchwey, Levine, Pollack and Laurencin2004), Cummings et al. (Reference Cummings, Sawyer, Morgan, Rose and Waters2009). Data for construct permeability are taken from Šimáček & Advani (Reference Šimáček and Advani1996), Sucosky et al. (Reference Sucosky, Osorio, Brown and Neitzel2004), Nabovati, Llewellin & Sousa (Reference Nabovati, Llewellin and Sousa2009). We use the standard acceleration due to gravity, and the density and viscosity values of water at room temperature for $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D707}$ , respectively. Although the dimensionless groupings can vary significantly in magnitude across the range of rotating bioreactors, we are interested in the HARV bioreactor where $\unicode[STIX]{x1D716}$ and $\unicode[STIX]{x1D716}_{1}$ are small (Cummings et al. Reference Cummings, Sawyer, Morgan, Rose and Waters2009).

Additionally, tissue engineering applications require as large a permeability as possible to enhance nutrient delivery to cells, while retaining the structural integrity of the construct. Thus, $K$ is on the larger side of its range while still ensuring that Darcy’s law is applicable within the construct. The flow within our construct is governed by Darcy’s law because the underlying solid matrix of our porous medium is connected (Auriault Reference Auriault2009); a tissue construct requires the structural integrity provided by a connected solid matrix. We follow Dalwadi et al. (Reference Dalwadi, Chapman, Waters and Oliver2016) and consider small values of $K$ , but do not exploit the asymptotic limit in our analysis. That is, we treat $K$ as being of order unity in the asymptotic limits mentioned above for algebraic convenience, rather than scaling $K$ with one of these small parameters. The latter treatment would be highly unwieldy, as discussed in Dalwadi et al. (Reference Dalwadi, Chapman, Waters and Oliver2016). Finally, we note that whilst the Bond number can be very large, its large size does not affect our analysis because it is absorbed into the reduced pressure. The Bond number only plays a role in the tissue construct force balance in § 4 where it is multiplied by $\unicode[STIX]{x1D716}_{1}$ .

To summarize, the physically relevant asymptotic limit we consider is: $\unicode[STIX]{x1D716}\ll 1\ll Re$ , $\unicode[STIX]{x1D716}_{1}\ll \unicode[STIX]{x1D716}Re\ll 1$ , and $\unicode[STIX]{x1D716}_{1}Re\ll 1$ . We now investigate the coupled system of equations using an asymptotic analysis, requiring matched asymptotic expansions for the flow problem, and a multiple-scales analysis for the trajectory problem. We begin with the flow equations.

3.1 Solution at leading order: $O(1)$

The leading-order system is defined by (3.2) and (3.4). We follow Cummings & Waters (Reference Cummings and Waters2007), who consider a similar system, and transform to bipolar coordinates $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ , defined by

(3.8a,b ) $$\begin{eqnarray}X=\cosh \unicode[STIX]{x1D6FC}_{1}-{\displaystyle \frac{\sinh \unicode[STIX]{x1D6FC}_{1}\sinh \unicode[STIX]{x1D6FC}}{\cosh \unicode[STIX]{x1D6FC}-\cos \unicode[STIX]{x1D6FD}}},\quad Y={\displaystyle \frac{\sinh \unicode[STIX]{x1D6FC}_{1}\sin \unicode[STIX]{x1D6FD}}{\cosh \unicode[STIX]{x1D6FC}-\cos \unicode[STIX]{x1D6FD}}},\end{eqnarray}$$

where

(3.9a,b ) $$\begin{eqnarray}b={\displaystyle \frac{\sinh \unicode[STIX]{x1D6FC}_{1}}{\sinh \unicode[STIX]{x1D6FC}_{2}}},\quad R=\cosh \unicode[STIX]{x1D6FC}_{1}-\sinh \unicode[STIX]{x1D6FC}_{1}\coth \unicode[STIX]{x1D6FC}_{2}.\end{eqnarray}$$

The bipolar transformation is useful here as it transforms the domain of two non-concentric circles, one contained within the other, to a rectangular domain, allowing us to obtain an analytic solution. We show a schematic of the transformed domains in figure 3(b). The bulk flow region now corresponds to the finite rectangle defined by $0<\unicode[STIX]{x1D6FC}_{1}<\unicode[STIX]{x1D6FC}<\unicode[STIX]{x1D6FC}_{2}$ and $0<\unicode[STIX]{x1D6FD}<2\unicode[STIX]{x03C0}$ , and the porous region corresponds to the semi-infinite rectangle defined by $\unicode[STIX]{x1D6FC}_{2}<\unicode[STIX]{x1D6FC}<\infty$ and $0<\unicode[STIX]{x1D6FD}<2\unicode[STIX]{x03C0}$ . Additionally, both regions now require periodic boundary conditions at $\unicode[STIX]{x1D6FD}=0$ and $2\unicode[STIX]{x03C0}$ . Hence, the bioreactor boundary $\unicode[STIX]{x2202}D_{1}$ corresponds to $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FC}_{1}$ , and the tissue construct boundary $\unicode[STIX]{x2202}D_{2}$ corresponds to $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FC}_{2}$ . On both $\unicode[STIX]{x2202}D_{1}$ and $\unicode[STIX]{x2202}D_{2}$ , we find that $\boldsymbol{n}=-\boldsymbol{e}_{\unicode[STIX]{x1D6FC}}$ and $\boldsymbol{s}=-\boldsymbol{e}_{\unicode[STIX]{x1D6FD}}$ .

Laplace’s equation in (3.2) is invariant under the conformal bipolar transformation, and these governing equations must be solved subject to the boundary condition (3.4a ) on $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FC}_{1}$ , given by

(3.10) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}p_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}=0,\end{eqnarray}$$

along with the coupling conditions (3.4b ) and (3.4c ) on $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FC}_{2}$ , which become

(3.11a ) $$\begin{eqnarray}\displaystyle & \displaystyle p_{0}=P_{0}, & \displaystyle\end{eqnarray}$$

(3.11b ) $$\begin{eqnarray}\displaystyle & \displaystyle -{\displaystyle \frac{1}{12}}\frac{\unicode[STIX]{x2202}p_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}+K\frac{\unicode[STIX]{x2202}P_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}={\displaystyle \frac{\sinh \unicode[STIX]{x1D6FC}_{1}(R(1-\dot{\unicode[STIX]{x1D711}})\sin \unicode[STIX]{x1D6FD}\sinh \unicode[STIX]{x1D6FC}_{2}+\dot{R}(\cos \unicode[STIX]{x1D6FD}\cosh \unicode[STIX]{x1D6FC}_{2}-1))}{(\cosh \unicode[STIX]{x1D6FC}_{2}-\cos \unicode[STIX]{x1D6FD})^{2}}}. & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The system (3.2), (3.10)–(3.11) is solved by

(3.12a ) $$\begin{eqnarray}\displaystyle & \displaystyle p_{0}=\mathop{\sum }_{n=1}^{\infty }[A_{n}\cos n\unicode[STIX]{x1D6FD}+B_{n}\sin n\unicode[STIX]{x1D6FD}]\cosh n(\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FC}_{1}), & \displaystyle\end{eqnarray}$$

(3.12b ) $$\begin{eqnarray}\displaystyle & \displaystyle P_{0}=\mathop{\sum }_{n=1}^{\infty }[C_{n}\cos n\unicode[STIX]{x1D6FD}+D_{n}\sin n\unicode[STIX]{x1D6FD}]\text{e}^{-n(\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FC}_{2})}, & \displaystyle\end{eqnarray}$$

$A_{n}$

$B_{n}$

$C_{n}$

$D_{n}$

(3.13) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle A_{n}=L_{n}\dot{R}\operatorname{sech}n(\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{1}),\quad B_{n}=L_{n}R(1-\dot{\unicode[STIX]{x1D711}})\operatorname{sech}n(\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{1}),\\ \displaystyle C_{n}=L_{n}\dot{R},\quad D_{n}=L_{n}R(1-\dot{\unicode[STIX]{x1D711}}),\quad L_{n}=-{\displaystyle \frac{24\sinh \unicode[STIX]{x1D6FC}_{1}\text{e}^{-n\unicode[STIX]{x1D6FC}_{2}}}{\tanh n(\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{1})+12K}}.\end{array}\right\}\end{eqnarray}$$

In Cummings & Waters (Reference Cummings and Waters2007), a solid tissue construct is considered. In the limit of $K\rightarrow 0$ , corresponding to a solid tissue construct, the coefficients $A_{n}$ and $B_{n}$ in (3.13) tend to those obtained in Cummings & Waters (Reference Cummings and Waters2007).

3.2 Solution at first correction: $O(\unicode[STIX]{x1D716}Re)$

We now solve the $O(\unicode[STIX]{x1D716}Re)$ system, given by (3.6)–(3.7), again using the bipolar transformation (3.8). We introduce

(3.14) $$\begin{eqnarray}\unicode[STIX]{x1D6F9}=p_{1}-{\displaystyle \frac{1}{2}}(X^{2}+Y^{2})+{\displaystyle \frac{3}{560}}|\unicode[STIX]{x1D735}p_{0}|^{2},\end{eqnarray}$$

and note that Laplace’s equation in (3.6) is invariant under the conformal bipolar transformation, yielding the governing equations

(3.15a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D6F9}=0\quad \unicode[STIX]{x1D6FC}_{1}<\unicode[STIX]{x1D6FC}<\unicode[STIX]{x1D6FC}_{2},\quad 0<\unicode[STIX]{x1D6FD}<2\unicode[STIX]{x03C0}, & \displaystyle\end{eqnarray}$$

(3.15b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}_{\bot }^{2}P_{1}=0\quad \unicode[STIX]{x1D6FC}_{2}<\unicode[STIX]{x1D6FC}<\infty ,\quad 0<\unicode[STIX]{x1D6FD}<2\unicode[STIX]{x03C0}. & \displaystyle\end{eqnarray}$$

(3.16a ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}{\unicode[STIX]{x2202}n}={\displaystyle \frac{(1-\dot{\unicode[STIX]{x1D711}})}{10}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}n}\frac{\unicode[STIX]{x2202}p_{0}}{\unicode[STIX]{x2202}s}-{\displaystyle \frac{1}{5}}\frac{\unicode[STIX]{x2202}p_{0}}{\unicode[STIX]{x2202}s}\quad \text{on }\unicode[STIX]{x2202}D_{1},\end{eqnarray}$$

and the continuity of pressure (3.7b ) and continuity of flux (3.7c ) conditions at the construct interface become

(3.16b ) $$\begin{eqnarray}\unicode[STIX]{x1D6F9}-P_{1}=-{\displaystyle \frac{X^{2}+Y^{2}}{2}}+{\displaystyle \frac{3|\unicode[STIX]{x1D735}p_{0}|^{2}}{560}}+\unicode[STIX]{x1D6F1}(s,t)(KP_{0n})^{2}H(-P_{0n})\quad \text{on }\unicode[STIX]{x2202}D_{2},\end{eqnarray}$$

(3.16c ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}{\unicode[STIX]{x2202}n}-12K\frac{\unicode[STIX]{x2202}P_{1}}{\unicode[STIX]{x2202}n} & = & \displaystyle {\displaystyle \frac{\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{B}(s,t)}{10}}\nonumber\\ \displaystyle & & \displaystyle -\,{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}}(12KP_{0n}\{p_{0s}\unicode[STIX]{x1D6EC}_{a}+\boldsymbol{s}\boldsymbol{\cdot }\boldsymbol{u}^{slip}\unicode[STIX]{x1D6EC}_{b}\})H(-P_{0n})\quad \text{on }\unicode[STIX]{x2202}D_{2},\qquad\end{eqnarray}$$

where $\boldsymbol{B}$ is defined in (A 31b ), $H(x)$ is the Heaviside function, $\boldsymbol{u}^{slip}$ is the velocity of the bioreactor walls relative to the construct defined in (A 9a ).

As with the $O(1)$ equations, it is convenient to solve this problem in terms of Fourier series by transforming to the bipolar coordinate system described by (3.8)–(3.9). Thus, we rewrite our boundary conditions as

(3.17a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}{\unicode[STIX]{x2202}n}=-{\displaystyle \frac{\cosh \unicode[STIX]{x1D6FC}_{1}-\cos \unicode[STIX]{x1D6FD}}{\sinh \unicode[STIX]{x1D6FC}_{1}}}\left(a_{10}+\mathop{\sum }_{n=1}^{\infty }(a_{1n}\cos n\unicode[STIX]{x1D6FD}+b_{1n}\sin n\unicode[STIX]{x1D6FD})\right)\quad \text{on }\unicode[STIX]{x2202}D_{1},\qquad & \displaystyle\end{eqnarray}$$

(3.17b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F9}-P_{1}=a_{20}+\mathop{\sum }_{n=1}^{\infty }(a_{2n}\cos n\unicode[STIX]{x1D6FD}+b_{2n}\sin n\unicode[STIX]{x1D6FD})\quad \text{on }\unicode[STIX]{x2202}D_{2}, & \displaystyle\end{eqnarray}$$

(3.17c ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}{\unicode[STIX]{x2202}n}-12K\frac{\unicode[STIX]{x2202}P_{1}}{\unicode[STIX]{x2202}n} & = & \displaystyle -{\displaystyle \frac{\cosh \unicode[STIX]{x1D6FC}_{2}-\cos \unicode[STIX]{x1D6FD}}{\sinh \unicode[STIX]{x1D6FC}_{1}}}\nonumber\\ \displaystyle & & \displaystyle \times \,\left(a_{30}+\displaystyle \mathop{\sum }_{n=1}^{\infty }(a_{3n}\cos n\unicode[STIX]{x1D6FD}+b_{3n}\sin n\unicode[STIX]{x1D6FD})\right)\quad \text{on }\unicode[STIX]{x2202}D_{2},\end{eqnarray}$$

$a_{in}(t)$

$b_{in}(t)$

$i\in \{1,2,3\}$

$a_{in}$

$b_{in}$

$\unicode[STIX]{x1D6FC}$

$\unicode[STIX]{x1D6FD}$

$a_{10}=a_{30}$

$\int _{\unicode[STIX]{x2202}D_{2}}\unicode[STIX]{x2202}P_{1}/\unicode[STIX]{x2202}n\,\text{d}s=0$

$a_{1n}(t)$

$b_{1n}(t)$

(3.18a-c ) $$\begin{eqnarray}a_{10}=0,\quad a_{1n}(t)={\displaystyle \frac{\dot{\unicode[STIX]{x1D711}}-3}{10}}nB_{n},\quad b_{1n}(t)={\displaystyle \frac{3-\dot{\unicode[STIX]{x1D711}}}{10}}nA_{n}.\end{eqnarray}$$

However, the expressions for $a_{2n}$ , $a_{3n}$ , $b_{2n}$ , and $b_{3n}$ are not as easily expressible in terms of analytically known coefficients; hence we determine these numerically.

The solution to (3.15) and (3.17) is

(3.19a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F9}=\unicode[STIX]{x1D6F9}_{1}+\unicode[STIX]{x1D6F9}_{2}, & \displaystyle\end{eqnarray}$$

(3.19b ) $$\begin{eqnarray}\displaystyle & \displaystyle P_{1}=\mathop{\sum }_{n=1}^{\infty }[\overline{E}_{n}\cos n\unicode[STIX]{x1D6FD}+\overline{F}_{n}\sin n\unicode[STIX]{x1D6FD}]\text{e}^{-n(\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FC}_{2})}, & \displaystyle\end{eqnarray}$$

(3.19c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F9}_{1}=\mathop{\sum }_{n=1}^{\infty }[\overline{A}_{n}\cos n\unicode[STIX]{x1D6FD}+\overline{B}_{n}\sin n\unicode[STIX]{x1D6FD}]\cosh n(\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FC}_{2}), & \displaystyle\end{eqnarray}$$

(3.19d ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F9}_{2}=a_{10}(\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FC}_{2})+a_{20}+\mathop{\sum }_{n=1}^{\infty }[\overline{C}_{n}\cos n\unicode[STIX]{x1D6FD}+\overline{D}_{n}\sin n\unicode[STIX]{x1D6FD}]\sinh n(\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FC}_{2}), & \displaystyle\end{eqnarray}$$

$\overline{A}_{n}$

$\overline{B}_{n}$

$\overline{C}_{n}$

$\overline{D}_{n}$

$\overline{E}_{n}$

$\overline{F}_{n}$

(3.20a ) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{A}_{n}=\overline{L}_{n}(a_{3n}-a_{1n}\operatorname{sech}n(\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{1})+12Kna_{2n}), & \displaystyle\end{eqnarray}$$

(3.20b ) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{B}_{n}=\overline{L}_{n}(b_{3n}-b_{1n}\operatorname{sech}n(\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{1})+12Knb_{2n}), & \displaystyle\end{eqnarray}$$

(3.20c ) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{C}_{n}=\overline{L}_{n}((a_{3n}+12Kna_{2n})\tanh n(\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{1})+12Ka_{1n}\operatorname{sech}n(\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{1})), & \displaystyle\end{eqnarray}$$

(3.20d ) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{D}_{n}=\overline{L}_{n}((b_{3n}+12Knb_{2n})\tanh n(\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{1})+12Kb_{1n}\operatorname{sech}n(\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{1})), & \displaystyle\end{eqnarray}$$

(3.20e ) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{E}_{n}=\overline{L}_{n}(a_{3n}-a_{1n}\operatorname{sech}n(\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{1})-na_{2n}\tanh n(\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{1})), & \displaystyle\end{eqnarray}$$

(3.20f ) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{F}_{n}=\overline{L}_{n}(b_{3n}-b_{1n}\operatorname{sech}n(\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{1})-nb_{2n}\tanh n(\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{1})), & \displaystyle\end{eqnarray}$$

(3.20g ) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{L}_{n}=1/(n(\tanh n(\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{1})+12K)), & \displaystyle\end{eqnarray}$$

$O(\unicode[STIX]{x1D716}Re)$

4.1 Periodic solutions for the trajectory

We recall that the tissue construct trajectories in HARV bioreactors appear to be periodic over a few orbits, but drift from this periodic orbit over the time scale of hours. Additionally, there are three types of periodic motion that are observed: steady, settling, and orbital. In steady motion, the tissue construct centre is stationary. In settling motion, the trajectory of the construct centre does not encircle the bioreactor centre. In orbital motion, the trajectory of the construct does encircle the bioreactor centre. In this subsection, we show that the solutions to the governing equation for the trajectory of the tissue construct centre of mass (4.10) are periodic up to $\unicode[STIX]{x1D716}_{1}$ . In later sections, we investigate the long-time drift from these periodic orbits using an asymptotic analysis, exploiting the small parameters $\unicode[STIX]{x1D716}_{1}$ and $\unicode[STIX]{x1D716}Re$ .

Up to $O(\unicode[STIX]{x1D716}_{1})$ , the system (4.10) only yields periodic solutions. This can be seen by noting that the steady points for this system are $(R,\unicode[STIX]{x1D711})=(R_{\ast },0)$ when $\overline{R}>0$ and $(R,\unicode[STIX]{x1D711})=(R_{\ast },\unicode[STIX]{x03C0})$ when $\overline{R}<0$ , where $R_{\ast }(1+\unicode[STIX]{x1D716}_{1}\unicode[STIX]{x1D6FE}_{1}(R_{\ast }))=|\overline{R}|$ . We note that, up to $O(\unicode[STIX]{x1D716}_{1})$ , $R_{\ast }\sim |\overline{R}|(1-\unicode[STIX]{x1D716}_{1}\unicode[STIX]{x1D6FE}_{1}(|\overline{R}|))$ . These steady points occur within the bioreactor when $R_{\ast }\leqslant |1-b|$ , with equality when the tissue construct touches the bioreactor edge. Moreover, a linear stability analysis suggests that these steady points are centres.

We can confirm that the steady points are true centres of the nonlinear system (4.10) up to the $O(\unicode[STIX]{x1D716}_{1})$ terms, by noting that there exists a first integral

(4.14) $$\begin{eqnarray}H+\overline{R}R\cos \unicode[STIX]{x1D711}={\displaystyle \frac{R^{2}}{2}}+\unicode[STIX]{x1D716}_{1}\int _{0}^{R}\hat{R}\unicode[STIX]{x1D6FE}_{1}(\hat{R})\,\text{d}\hat{R},\end{eqnarray}$$

where $H$ is a constant. To obtain (4.14), we form $\text{d}\unicode[STIX]{x1D711}/\text{d}R$ up to and including the $O(\unicode[STIX]{x1D716}_{1})$ terms in (4.10), then integrate with respect to $R$ . The first integral (4.14) defines closed orbits in $(R,\unicode[STIX]{x1D711})$ space, implying that the solutions to the nonlinear system (4.10) up to the $O(\unicode[STIX]{x1D716}_{1})$ terms must be periodic.

To calculate the long-time drift away from the periodic orbits due to inertia, we use the method of multiple scales on a two-variable system of nonlinear equations (4.10). However, before we tackle the full problem, we first interrogate a toy problem that captures the physical essence of the full problem, while still being amenable to an analytic solution. Although we do not expect the toy solution to completely describe the construct motion for the full problem, introducing the toy problem serves two purposes. Firstly, it will allow us to gain physical insight into the full problem and, secondly, it will allow us to set up the mathematical machinery used to solve the full problem.

4.2 Toy problem

The analytically difficult terms in the governing equations (4.10) are the infinite sums $\unicode[STIX]{x1D6FE}_{1}$ , $\unicode[STIX]{x1D6FE}_{2}$ , $\unicode[STIX]{x1D6FE}_{3}$ , $\unicode[STIX]{x1D6FE}_{4}$ , and $\unicode[STIX]{x1D6FE}_{5}$ (defined in (4.6b ), (4.8b,c ), (4.9b ), and (4.9c )) due to the fluid pressure and shear stress acting on the curved construct interface. If we ignore these infinite sums but include the inertial terms, the remaining forces acting on the system are due to gravity and friction on the flat sides of the construct. As the centrifugal force of the rotation on the construct (the $R\dot{\unicode[STIX]{x1D711}}^{2}$ term on the right-hand side of (4.10)) causes the construct to drift out to the bioreactor edge, throwing out the infinite sums we mention above means that we lose the potential balancing effect of the centrifugal force of the rotating fluid on the construct pushing the construct inwards for a construct that is lighter than the surrounding fluid. This is because the infinite sums contain all the information about the force exerted on the curved interface of the construct by the surrounding fluid.

For our toy model, we introduce an ad hoc centrifugal force acting on the construct due to the fluid. We assume that this force has the same form as the centrifugal force that appears in the $O(\unicode[STIX]{x1D716}Re)$ Navier–Stokes equations (A 21). That is, we take the centrifugal force in our toy model to be $-\unicode[STIX]{x1D716}Re\dot{\unicode[STIX]{x1D711}}^{2}R\boldsymbol{e}_{X}$ , which pushes the construct towards the bioreactor centre.

In summary, we use (4.10) for the basis of our toy model, but remove the infinite sums and surface integral, and replace them with $-\unicode[STIX]{x1D716}\unicode[STIX]{x1D716}_{1}Re\dot{\unicode[STIX]{x1D711}}^{2}R/2\boldsymbol{e}_{X}$ . We then obtain the following coupled ordinary differential equations in time for $R(t)$ and $\unicode[STIX]{x1D711}(t)$

(4.15a ) $$\begin{eqnarray}\displaystyle & \displaystyle -\dot{R}-\overline{R}\sin \unicode[STIX]{x1D711}={\displaystyle \frac{\unicode[STIX]{x1D716}_{1}\unicode[STIX]{x1D716}Re}{2}}(\bar{\unicode[STIX]{x1D70C}}\ddot{R}-(\bar{\unicode[STIX]{x1D70C}}-1)\dot{\unicode[STIX]{x1D711}}^{2}R), & \displaystyle\end{eqnarray}$$

(4.15b ) $$\begin{eqnarray}\displaystyle & \displaystyle R(1-\dot{\unicode[STIX]{x1D711}})-\overline{R}\cos \unicode[STIX]{x1D711}={\displaystyle \frac{\unicode[STIX]{x1D716}_{1}\unicode[STIX]{x1D716}Re\bar{\unicode[STIX]{x1D70C}}}{2}}(2\dot{\unicode[STIX]{x1D711}}\dot{R}+\ddot{\unicode[STIX]{x1D711}}R). & \displaystyle\end{eqnarray}$$

$O(1)$

$t=O(1)$

To solve the system (4.15), we introduce a polar coordinate system $(r(t),\unicode[STIX]{x1D703}(t))$ centred around the leading-order steady point to define the position of the construct centre. That is, we define

(4.16a,b ) $$\begin{eqnarray}R\cos \unicode[STIX]{x1D711}=\overline{R}+r\cos \unicode[STIX]{x1D703},\quad R\sin \unicode[STIX]{x1D711}=r\sin \unicode[STIX]{x1D703}.\end{eqnarray}$$

Differentiating (4.16) with respect to time allows us to transform (4.15) into

(4.17a ) $$\begin{eqnarray}\displaystyle & & \displaystyle -r\dot{r}-\overline{R}(\dot{r}\cos \unicode[STIX]{x1D703}+r(1-\dot{\unicode[STIX]{x1D703}})\sin \unicode[STIX]{x1D703})\nonumber\\ \displaystyle & & \displaystyle \quad =(\unicode[STIX]{x1D716}_{1}\unicode[STIX]{x1D716}Re){\displaystyle \frac{\bar{\unicode[STIX]{x1D70C}}}{2}}(r(\ddot{r}-r\dot{\unicode[STIX]{x1D703}}^{2})+\overline{R}((\ddot{r}-r\dot{\unicode[STIX]{x1D703}}^{2})\cos \unicode[STIX]{x1D703}-(2\dot{r}\dot{\unicode[STIX]{x1D703}}+r\ddot{\unicode[STIX]{x1D703}})\sin \unicode[STIX]{x1D703}))\nonumber\\ \displaystyle & & \displaystyle \qquad +\,{\displaystyle \frac{(\unicode[STIX]{x1D716}_{1}\unicode[STIX]{x1D716}Re)}{2}}\left({\displaystyle \frac{(r\dot{\unicode[STIX]{x1D703}}(r+\overline{R}\cos \unicode[STIX]{x1D703})+\overline{R}\dot{r}\sin \unicode[STIX]{x1D703})^{2}}{r^{2}+2r\overline{R}\cos \unicode[STIX]{x1D703}+\overline{R}^{2}}}\right),\end{eqnarray}$$

(4.17b ) $$\begin{eqnarray}\displaystyle & & \displaystyle r^{2}(1-\dot{\unicode[STIX]{x1D703}})+\overline{R}(r(1-\dot{\unicode[STIX]{x1D703}})\cos \unicode[STIX]{x1D703}-\dot{r}\sin \unicode[STIX]{x1D703})\nonumber\\ \displaystyle & & \displaystyle \quad =(\unicode[STIX]{x1D716}_{1}\unicode[STIX]{x1D716}Re){\displaystyle \frac{\bar{\unicode[STIX]{x1D70C}}}{2}}(r(2\dot{r}\dot{\unicode[STIX]{x1D703}}+r\ddot{\unicode[STIX]{x1D703}})+\overline{R}((\ddot{r}-r\dot{\unicode[STIX]{x1D703}}^{2})\sin \unicode[STIX]{x1D703}+(2\dot{r}\dot{\unicode[STIX]{x1D703}}+r\ddot{\unicode[STIX]{x1D703}})\cos \unicode[STIX]{x1D703})).\qquad\end{eqnarray}$$

We proceed by exploiting the small parameter $\unicode[STIX]{x1D6FF}:=\unicode[STIX]{x1D716}_{1}\unicode[STIX]{x1D716}Re$ . We note that the second-order system of ordinary differential equations (4.17) is singular in the limit $\unicode[STIX]{x1D6FF}\rightarrow 0$ . A standard boundary layer analysis for (4.17) when $t=O(\unicode[STIX]{x1D6FF})$ shows that the initial position is the correct matching condition for the $t=O(1)$ problem at leading order. We omit the analysis for brevity here, but details can be found in Dalwadi (Reference Dalwadi2014). We would only need to take the initial velocity into account if we wanted to calculate the $O(\unicode[STIX]{x1D6FF})$ correction terms in the construct trajectory.

Formally, we implement the standard procedure for the method of multiples scales (see, for example, Kevorkian & Cole (Reference Kevorkian and Cole1996)). We introduce the slow time scale $T=\unicode[STIX]{x1D6FF}t=O(1)$ , and refer to $t$ as the fast time scale. We treat each dependent variable as a function of both the fast and slow time scales, then remove the extra freedom this introduces by imposing periodicity in the fast time scale. Although this system is nonlinear, we do not need to use the full machinery of the method of Kuzmak (Reference Kuzmak1959) as we will find that the period of the fast time scale solution is not dependent on the slow time scale.

Figure 5. The three different periodic motion regimes at leading order, shown in the bioreactor domain, when $\overline{R}=0.25$ . The black cross denotes the position of the construct centre in the steady regime, when $r(0)=0$ , the blue line denotes a trajectory of the construct centre in the settling regime, when $r(0)=0.1<|\overline{R}|$ , and the yellow line denotes a trajectory of the construct centre in the orbital regime, when $r(0)=0.45>|\overline{R}|$ . The red dotted line denotes the boundary between orbital and settling motion; the periodic trajectory of the construct centre encloses the bioreactor centre for the former, but not for the latter. The dashed black line denotes the positions of the construct centre for which the construct hits the bioreactor edge, using $b=0.2$ .

The time derivative becomes

(4.18) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D6FF}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}T},\end{eqnarray}$$

and we seek an asymptotic expansion of the form

(4.19a ) $$\begin{eqnarray}\displaystyle & \displaystyle r(t,T)\sim r_{0}(t,T)+\unicode[STIX]{x1D6FF}r_{2}(t,T), & \displaystyle\end{eqnarray}$$

(4.19b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}(t,T)\sim \unicode[STIX]{x1D703}_{0}(t,T)+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}_{2}(t,T), & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D6FF}\rightarrow 0$

$1$

$O(1)$

(4.20a ) $$\begin{eqnarray}\displaystyle & \displaystyle r_{0}r_{0t}+\overline{R}(r_{0t}\cos \unicode[STIX]{x1D703}_{0}+r_{0}(1-\unicode[STIX]{x1D703}_{0t})\sin \unicode[STIX]{x1D703}_{0})=0, & \displaystyle\end{eqnarray}$$

(4.20b ) $$\begin{eqnarray}\displaystyle & \displaystyle r_{0}^{2}(1-\unicode[STIX]{x1D703}_{0t})+\overline{R}(r_{0}(1-\unicode[STIX]{x1D703}_{0t})\cos \unicode[STIX]{x1D703}_{0}-r_{0t}\sin \unicode[STIX]{x1D703}_{0})=0, & \displaystyle\end{eqnarray}$$

(4.21a,b ) $$\begin{eqnarray}r_{0}=r_{0}(T),\quad \unicode[STIX]{x1D703}_{0}=t+\unicode[STIX]{x1D719}(T),\end{eqnarray}$$

and the solutions are $2\unicode[STIX]{x03C0}$ -periodic on the short time scale $t$ . Thus, the leading-order solutions to the toy and full problems are circular orbits around the steady point (figure 5). Moreover, the construct motion regime over this short time scale is: steady if $r(0)=0$ (black cross), settling if $r(0)<|\overline{R}|$ (blue), and orbital if $r(0)>|\overline{R}|$ (yellow). Additionally, we note that the construct will hit the bioreactor wall if $r(0)+|\overline{R}|\geqslant 1-b$ (black dashed). This inequality provides parameter constraints to avoid the construct hitting the bioreactor wall on the fast time scale. Thus, we require $r(0)+|\overline{R}|<1-b$ for the construct orbit to be entirely contained within the bioreactor on the fast time scale.

Recalling that $\overline{R}=\overline{B}(\bar{\unicode[STIX]{x1D70C}}-1)$ , we may deduce from the inequality above how the values that $r(0)$ can take depend on the values of $\overline{B}$ and $\bar{\unicode[STIX]{x1D70C}}$ , and additionally that the steady point is only within the allowable physical domain when $\overline{B}|\bar{\unicode[STIX]{x1D70C}}-1|<1-b$ (figure 6). We further show in figure 6(b) how the construct motion regime depends on $r_{0}$ , $\overline{B}$ , and $\bar{\unicode[STIX]{x1D70C}}$ . As the steady point $\overline{R}=\overline{B}(\bar{\unicode[STIX]{x1D70C}}-1)$ increases in magnitude, there are more construct trajectories that will hit the bioreactor edge.

Figure 6. The parameter constraints to keep the construct from hitting the bioreactor edge. The white regions denote valid values and the grey regions denote invalid values. (a) $(\bar{\unicode[STIX]{x1D70C}},\overline{B})$ space: the white region is defined by $\overline{B}|1-\bar{\unicode[STIX]{x1D70C}}|<1-b$ . (b) $(\bar{\unicode[STIX]{x1D70C}},r_{0})$ space: the white region is defined by $r_{0}+\overline{B}|1-\bar{\unicode[STIX]{x1D70C}}|<1-b$ . The dotted red lines at $r_{0}=\overline{B}|1-\bar{\unicode[STIX]{x1D70C}}|$ define the boundaries between orbital and settling motion regimes, and the dashed blue line at $r_{0}=0$ denotes the steady regime.

The leading-order results (4.21) simplify the $O(\unicode[STIX]{x1D6FF})$ terms in (4.17), which become

(4.22a ) $$\begin{eqnarray}\displaystyle & \displaystyle r_{0}r_{2t}+\overline{R}(r_{2t}\cos \unicode[STIX]{x1D703}_{0}-r_{0}\unicode[STIX]{x1D703}_{2t}\sin \unicode[STIX]{x1D703}_{0})=\unicode[STIX]{x1D707}(r_{0},\unicode[STIX]{x1D703}_{0}), & \displaystyle\end{eqnarray}$$

(4.22b ) $$\begin{eqnarray}\displaystyle & \displaystyle r_{0}^{2}\unicode[STIX]{x1D703}_{2t}+\overline{R}(r_{0}\unicode[STIX]{x1D703}_{2t}\cos \unicode[STIX]{x1D703}_{0}+r_{2t}\sin \unicode[STIX]{x1D703}_{0})=\unicode[STIX]{x1D708}(r_{0},\unicode[STIX]{x1D703}_{0}), & \displaystyle\end{eqnarray}$$

(4.23a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}(r_{0},\unicode[STIX]{x1D703}_{0}) & = & \displaystyle {\displaystyle \frac{\bar{\unicode[STIX]{x1D70C}}}{2}}r_{0}(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0})-r_{0T}(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0})+r_{0}\unicode[STIX]{x1D703}_{0T}\overline{R}\sin \unicode[STIX]{x1D703}_{0}\nonumber\\ \displaystyle & & \displaystyle \quad -\,{\displaystyle \frac{(r_{0}(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0}))^{2}}{2(r_{0}^{2}+2r_{0}\overline{R}\cos \unicode[STIX]{x1D703}_{0}+\overline{R}^{2})}},\end{eqnarray}$$

(4.23b ) $$\begin{eqnarray}\unicode[STIX]{x1D708}(r_{0},\unicode[STIX]{x1D703}_{0})={\displaystyle \frac{\bar{\unicode[STIX]{x1D70C}}}{2}}r_{0}\overline{R}\sin \unicode[STIX]{x1D703}_{0}-r_{0}\unicode[STIX]{x1D703}_{0T}(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0})-r_{0T}\overline{R}\sin \unicode[STIX]{x1D703}_{0}.\end{eqnarray}$$

To determine equations for $r_{0}(T)$ and $\unicode[STIX]{x1D719}(T)$ , we need to impose the secularity condition that arises from the assumption that the solution to the system (4.22) is $2\unicode[STIX]{x03C0}$ -periodic in $t$ . This is equivalent to imposing orthogonality of the solutions of the adjoint problem under periodic boundary conditions via the Fredholm Alternative Theorem. Since this procedure is less well known for systems, we first describe it in abstract terms.

Consider the linear system for a 2D vector function $(x_{1},x_{2})$ , as follows

(4.24) $$\begin{eqnarray}L(x)=\left(\begin{array}{@{}cc@{}}L_{1} & L_{2}\\ L_{3} & L_{4}\end{array}\right)\left(\begin{array}{@{}c@{}}x_{1}\\ x_{2}\end{array}\right)=\left(\begin{array}{@{}c@{}}f_{1}\\ f_{2}\end{array}\right),\end{eqnarray}$$

where $f_{1},f_{2},x_{1},x_{2}$ are functions of $t\in (0,2\unicode[STIX]{x03C0})$ , $x_{1}$ and $x_{2}$ have periodic boundary conditions on $t=0$ , $2\unicode[STIX]{x03C0}$ , and the $L_{i}$ are linear operators. If there is a non-trivial solution of the adjoint problem

(4.25) $$\begin{eqnarray}\left(\begin{array}{@{}cc@{}}L_{1}^{\ast } & L_{3}^{\ast }\\ L_{2}^{\ast } & L_{4}^{\ast }\end{array}\right)\left(\begin{array}{@{}c@{}}y_{1}\\ y_{2}\end{array}\right)=\left(\begin{array}{@{}c@{}}0\\ 0\end{array}\right),\end{eqnarray}$$

with periodic boundary conditions, the original problem has solvability condition

(4.26) $$\begin{eqnarray}\int _{0}^{2\unicode[STIX]{x03C0}}(y_{1}f_{1}+y_{2}f_{2})\,\text{d}t=0.\end{eqnarray}$$

For the $O(\unicode[STIX]{x1D6FF})$ system defined in (4.22), the solvability condition (4.26) is

(4.27a,b ) $$\begin{eqnarray}\int _{0}^{2\unicode[STIX]{x03C0}}{\displaystyle \frac{(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0})\unicode[STIX]{x1D707}+\overline{R}\sin \unicode[STIX]{x1D703}_{0}\unicode[STIX]{x1D708}}{\overline{R}^{2}+2r_{0}\overline{R}\cos \unicode[STIX]{x1D703}_{0}+r_{0}^{2}}}\,\text{d}t=0,\quad \int _{0}^{2\unicode[STIX]{x03C0}}{\displaystyle \frac{(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0})\unicode[STIX]{x1D708}-\overline{R}\sin \unicode[STIX]{x1D703}_{0}\unicode[STIX]{x1D707}}{r_{0}(\overline{R}^{2}+2r_{0}\overline{R}\cos \unicode[STIX]{x1D703}_{0}+r_{0}^{2})}}\,\text{d}t=0.\end{eqnarray}$$

Substituting (4.23) into (4.27) yields the solvability conditions

(4.28a ) $$\begin{eqnarray}\displaystyle & \displaystyle r_{0T}={\displaystyle \frac{\bar{\unicode[STIX]{x1D70C}}}{2}}r_{0}-{\displaystyle \frac{r_{0}^{2}}{4\unicode[STIX]{x03C0}}}\int _{0}^{2\unicode[STIX]{x03C0}}{\displaystyle \frac{(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0})^{3}}{(r_{0}^{2}+2r_{0}\overline{R}\cos \unicode[STIX]{x1D703}_{0}+\overline{R}^{2})^{2}}}\,\text{d}t, & \displaystyle\end{eqnarray}$$

(4.28b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}_{0T}={\displaystyle \frac{\overline{R}r_{0}}{4\unicode[STIX]{x03C0}}}\int _{0}^{2\unicode[STIX]{x03C0}}{\displaystyle \frac{\sin \unicode[STIX]{x1D703}_{0}(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0})^{2}}{(r_{0}^{2}+2r_{0}\overline{R}\cos \unicode[STIX]{x1D703}_{0}+\overline{R}^{2})^{2}}}\,\text{d}t. & \displaystyle\end{eqnarray}$$

$t=\unicode[STIX]{x03C0}-\unicode[STIX]{x1D703}(0)$

(4.29) $$\begin{eqnarray}\unicode[STIX]{x1D703}_{0T}=0,\end{eqnarray}$$

so that $\unicode[STIX]{x1D703}_{0}$ is independent of $T$ . We evaluate the integral in (4.28a ) using contour integration, obtaining

(4.30) $$\begin{eqnarray}\displaystyle \int _{0}^{2\unicode[STIX]{x03C0}}{\displaystyle \frac{(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0})^{3}\,\text{d}t}{(r_{0}^{2}+2r_{0}\overline{R}\cos \unicode[STIX]{x1D703}_{0}+\overline{R}^{2})^{2}}} & = & \displaystyle \int _{0}^{2\unicode[STIX]{x03C0}}{\displaystyle \frac{(r_{0}+\overline{R}\cos t)^{3}\,\text{d}t}{(r_{0}^{2}+2r_{0}\overline{R}\cos t+\overline{R}^{2})^{2}}}\nonumber\\ \displaystyle & = & \displaystyle \left\{\begin{array}{@{}ll@{}}\displaystyle 0\quad & \displaystyle \text{if }r_{0}<|\overline{R}|,\\ \displaystyle {\displaystyle \frac{\unicode[STIX]{x03C0}}{r_{0}^{3}}}(2r_{0}^{2}-\overline{R}^{2})\quad & \displaystyle \text{if }r_{0}>|\overline{R}|.\end{array}\right.\end{eqnarray}$$

We note that the result depends on whether $r_{0}<|\overline{R}|$ or $r_{0}>|\overline{R}|$ , i.e. whether the construct is in settling or orbital motion. This is because the integrand in (4.30) stems from the centrifugal force of the fluid acting on the construct (which we introduced as proportional to $\dot{\unicode[STIX]{x1D711}}^{2}$ ), and the discontinuity of the integral arises because the effect of centrifugal force averaged over one periodic orbit is discontinuous as we move from $r_{0}<|\overline{R}|$ (settling) to $r_{0}>|\overline{R}|$ (orbital). To get a sense of why this is the case, note that the average of $\dot{\unicode[STIX]{x1D711}}$ over one orbit is $0$ for settling motion and $1$ for orbital motion. Although the introduced centrifugal force has a factor of $\dot{\unicode[STIX]{x1D711}}^{2}$ and not $\dot{\unicode[STIX]{x1D711}}$ , the discontinuity across motion regimes will transfer from $\dot{\unicode[STIX]{x1D711}}$ to $\dot{\unicode[STIX]{x1D711}}^{2}$ .

Using the integral result (4.30) in the secularity condition (4.28a ), we obtain the governing equation

(4.31) $$\begin{eqnarray}\frac{\text{d}r_{0}}{\text{d}T}=\left\{\begin{array}{@{}ll@{}}\displaystyle {\displaystyle \frac{\bar{\unicode[STIX]{x1D70C}}}{2}}r_{0}\quad & \displaystyle \text{if }r_{0}<|\overline{R}|,\\ \displaystyle {\displaystyle \frac{\bar{\unicode[STIX]{x1D70C}}-1}{2}}r_{0}+{\displaystyle \frac{\overline{R}^{2}}{4r_{0}}}\quad & \displaystyle \text{if }r_{0}>|\overline{R}|.\end{array}\right.\end{eqnarray}$$

As (4.30) is discontinuous at $r_{0}=|\overline{R}|$ , (4.31) is also discontinuous at the same point. We note that $r_{0T}$ scales linearly with $r_{0}$ for $r_{0}<|\overline{R}|$ (settling regime). In the orbital regime, when $r_{0}>|\overline{R}|$ , we see that $r_{0T}$ tends to a linear function of $r_{0}$ with gradient $(\bar{\unicode[STIX]{x1D70C}}-1)/2$ as $r_{0}$ increases, decaying to this function at a decreasing rate.

Though we can solve (4.31) analytically as follows

(4.32) $$\begin{eqnarray}r_{0}=\left\{\begin{array}{@{}ll@{}}\displaystyle r_{0}(0)\exp \left({\displaystyle \frac{\bar{\unicode[STIX]{x1D70C}}T}{2}}\right)\quad & \displaystyle \text{if }r_{0}<|\overline{R}|,\\ \displaystyle \left(\left(r_{0}^{2}(0)+{\displaystyle \frac{\overline{R}^{2}}{2(\bar{\unicode[STIX]{x1D70C}}-1)}}\right)\exp ((\bar{\unicode[STIX]{x1D70C}}-1)T)-{\displaystyle \frac{\overline{R}^{2}}{2(\bar{\unicode[STIX]{x1D70C}}-1)}}\right)^{1/2}\quad & \displaystyle \text{if }r_{0}>|\overline{R}|,\end{array}\right.\end{eqnarray}$$

noting that the construct hits the bioreactor wall if $r_{0}+|\overline{R}|\geqslant 1-b$ , it is constructive to analyse the long-time behaviour of the system by investigating the sign of the right-hand side of (4.31). This is also the route we take in the full problem, when an analytic solution is not feasible. We consider the system behaviour under the independent variation of $\overline{B}$ and $\bar{\unicode[STIX]{x1D70C}}$ . These parameters serve as proxies for the independent variation of $\unicode[STIX]{x1D6FA}$ and $\unicode[STIX]{x1D70C}_{c}$ , the bioreactor angular velocity and the construct density, respectively, which are the main experimentally controllable parameters. We note that $\overline{B}$ is inversely proportional to $\unicode[STIX]{x1D6FA}$ , and $\bar{\unicode[STIX]{x1D70C}}$ is directly proportional to $\unicode[STIX]{x1D70C}_{c}$ . As $\overline{R}=\overline{B}(\bar{\unicode[STIX]{x1D70C}}-1)$ , we emphasize that a change in either of these experimental parameters will change the steady point of the system. We now discuss the possible long-time behaviours for the toy problem.

Figure 7. Stability and bifurcation diagrams showing the possible long-time behaviours for the toy problem. In the unshaded region, the construct is in one of two general motion behaviours separated by the dotted blue lines, as discussed in the main text. In the shaded regions, the construct hits the bioreactor edge during the first attempted orbit, and our model breaks down. (a) The long-time construct behaviour in $(\bar{\unicode[STIX]{x1D70C}},\overline{B})$ parameter space. The shaded regions have boundaries defined by $\overline{B}|1-\bar{\unicode[STIX]{x1D70C}}|=1-b$ . (b) A bifurcation diagram of $r_{0}$ against $\bar{\unicode[STIX]{x1D70C}}$ , for fixed $\overline{B}\in ((1-b)/2,1-b)$ . Over the long time, $r_{0}$ will vary in the direction of the arrows. The solid black line represents a stable limit cycle and the dashed black line an unstable spiral. The constraints and short-time motion regimes are the same as those in figure 6. That is, the orbital and settling motion regimes are above and below the dotted red line, respectively, which is obscured by the stable limit cycle for behaviour Ia, and the shaded regions have boundaries defined by $r_{0}+\overline{B}|1-\bar{\unicode[STIX]{x1D70C}}|=1-b$ .

We present the bifurcation results for the toy problem in figure 7 and a schematic of the long-time behaviours within the bioreactor in figure 8. To briefly summarize figure 7 before going into the details, the dotted blue vertical line at $\bar{\unicode[STIX]{x1D70C}}=1$ denotes a change in the stability of the steady point at $r_{0}=0$ and the death of a stable limit cycle, the dotted yellow vertical line at $\bar{\unicode[STIX]{x1D70C}}=0.5$ arises from the location of the stable limit cycle moving fully into the orbital regime from the boundary between settling and orbital motions, and the remaining blue dotted line arises from the stable limit cycle moving outside the bioreactor domain, leading to all trajectories eventually hitting the bioreactor edge. To evaluate the long-time behaviour, it is helpful to first understand the system without the constraint of the construct hitting the bioreactor wall when $r_{0}+|\overline{R}|\geqslant 1-b$ . From (4.31) we see that while in settling motion ( $r_{0}<|\overline{R}|$ ), the construct will always spiral out until it reaches the boundary between settling and orbital motion at $r_{0}=|\overline{R}|$ in finite time. This boundary is the trajectory that goes through the centre of the bioreactor. Thus, the point $r_{0}=0$ is always an unstable spiral. Depending on the parameter values, several things can occur when the construct reaches the boundary between settling and orbital motion, which we now examine.

Figure 8. A schematic of the three possible motion regimes for the toy problem. In Regime Ia, the construct always tends to a stable limit cycle on the boundary between settling and orbital motion. In Regime Ib, the construct always tends to a stable limit cycle in orbital motion. In Regime II, the construct always spirals out until it hits the bioreactor edge. The crosses denote an unstable steady point and the blue solid lines denote a stable limit cycle.

When in orbital motion ( $r_{0}>|\overline{R}|$ ), the long-time construct trajectory depends on the value of $\bar{\unicode[STIX]{x1D70C}}$ , and is separated into two main behaviours. If $\bar{\unicode[STIX]{x1D70C}}<1$ , the construct will tend to a stable limit cycle (behaviour I in figure 7). If $\bar{\unicode[STIX]{x1D70C}}>1$ , the right-hand side of (4.31) is always positive, and therefore the construct will always spiral outwards (behaviour II in figure 7). The location of the stable limit cycle that appears for $\bar{\unicode[STIX]{x1D70C}}<1$ depends on the value of $\bar{\unicode[STIX]{x1D70C}}$ , and exhibits two sub-behaviours. If $\bar{\unicode[STIX]{x1D70C}}<1/2$ , then we see from (4.31) that $r_{0T}<0$ for all $r_{0}>|\overline{R}|$ . Thus, whether in orbital or settling motion, the construct will always tend towards the orbit which lies on the boundary between these two motion regimes when $\bar{\unicode[STIX]{x1D70C}}<1/2$ . We refer to this as behaviour Ia in figure 7. For $\bar{\unicode[STIX]{x1D70C}}\in (1/2,1)$ , we see from (4.31) that there is a critical value of $r_{0}$ for which $r_{0T}=0$ , occurring when

(4.33) $$\begin{eqnarray}r_{0}^{c}={\displaystyle \frac{|\overline{R}|}{\sqrt{2(1-\bar{\unicode[STIX]{x1D70C}})}}}=\overline{B}\sqrt{{\displaystyle \frac{1-\bar{\unicode[STIX]{x1D70C}}}{2}}},\end{eqnarray}$$

where we use (4.11) for the second equality in (4.33). As $\unicode[STIX]{x2202}r_{0T}/\unicode[STIX]{x2202}r_{0}<0$ at $r_{0}=r_{0}^{c}$ (which is always within the orbital regime for $\bar{\unicode[STIX]{x1D70C}}\in (1/2,1)$ ), this critical value corresponds to a stable limit cycle for the trajectory in the orbital regime. This behaviour is marked as Ib in figure 7.

Thus, there is a bifurcation in the system at $\bar{\unicode[STIX]{x1D70C}}=1$ ; a stable limit cycle grows out of the bioreactor centre for $\bar{\unicode[STIX]{x1D70C}}<1$ , and is annihilated as $\bar{\unicode[STIX]{x1D70C}}\rightarrow 1^{-}$ . As the scaling $r_{0}\sim \overline{B}$ removes $\overline{B}$ from the system (4.31), figure 7(b) is valid for all $\overline{B}$ , and the geometrical effect of varying $\overline{B}$ in figure 7(b) is to vary the intersections of both the grey boundary and the dotted red line with $\bar{\unicode[STIX]{x1D70C}}=0$ .

In this toy problem, the construct can reach the stable limit cycle at the boundary between settling and orbital motion for $\bar{\unicode[STIX]{x1D70C}}<1/2$ in finite rather than infinite time. This is because the right-hand side of (4.31), the equation governing the long-time behaviour of the construct, is discontinuous at this boundary. Rather than $\text{d}r_{0}/\text{d}T=0$ at $r_{0}=|\overline{R}|$ , we instead have $\text{d}r_{0}/\text{d}T>0$ for $r_{0}\rightarrow |\overline{R}|^{-}$ and $\text{d}r_{0}/\text{d}T<0$ for $r_{0}\rightarrow |\overline{R}|^{+}$ when $\bar{\unicode[STIX]{x1D70C}}<1/2$ .

We now consider the constraint $r_{0}+|\overline{R}|<1-b$ , required to keep the construct from hitting the bioreactor wall. More trajectories hit the bioreactor wall when $|\overline{R}|=\overline{B}|1-\bar{\unicode[STIX]{x1D70C}}|$ is larger, and this would cause the area of the grey region in figure 7(b) to increase as $\overline{B}$ is increased. This corresponds to the steady point moving towards the bioreactor edge, and would eventually cause the stable limit cycle to hit the bioreactor edge. Keeping $\bar{\unicode[STIX]{x1D70C}}$ constant, and increasing $\overline{B}$ from zero in figure 7(b), we can see when this would occur and quantify the critical parameter relationship. When the boundary of the leftmost grey region (with relationship $r_{0}+\overline{B}(1-\bar{\unicode[STIX]{x1D70C}})=1-b$ ) intersects the stable limit cycle in Ia (the boundary between settling and orbital motion with relationship $r_{0}=\overline{B}(1-\bar{\unicode[STIX]{x1D70C}})$ ), the parameter relationship is

(4.34a ) $$\begin{eqnarray}\overline{B}={\displaystyle \frac{1-b}{2(1-\bar{\unicode[STIX]{x1D70C}})}}\quad \text{for }\bar{\unicode[STIX]{x1D70C}}\in (0,1/2),\end{eqnarray}$$

yielding the boundary between Ia and II in figure 7(a). Similarly, the boundary of the leftmost grey region intersects the stable limit cycle in Ib (with relationship (4.33)) when

(4.34b ) $$\begin{eqnarray}\overline{B}={\displaystyle \frac{1-b}{1-\bar{\unicode[STIX]{x1D70C}}+\sqrt{{\displaystyle \frac{1-\bar{\unicode[STIX]{x1D70C}}}{2}}}}}\quad \text{for }\bar{\unicode[STIX]{x1D70C}}\in (1/2,1),\end{eqnarray}$$

yielding the boundary between Ib and II in figure 7(a).

Our toy model suggests that the construct will never remain within the settling regime, because it will instead either tend to the orbital regime or the bioreactor edge in finite time. Further, if $\bar{\unicode[STIX]{x1D70C}}>1$ , the construct will hit the bioreactor edge in finite time no matter where it starts. If $\bar{\unicode[STIX]{x1D70C}}<1$ , the long-time behaviour will depend on the parameter values of $\bar{\unicode[STIX]{x1D70C}}$ , $\overline{B}$ , and $b$ , either tending to a stable limit cycle or hitting the bioreactor edge, as described in figure 7. Notably, our toy model suggests that drift will always occur in the HARV bioreactor due to weak inertia, and is thus not solely down to cell growth varying the construct density. Additionally, the discontinuity in the long-time ordinary differential equations we derive (4.31), at the boundary between the settling and orbital regimes, suggests that we may encounter a sharp change in construct behaviour between the settling and orbital regimes.

In the next subsection, we analyse the full problem as defined in (4.10) using the method of multiple scales. We then discuss the results and how they compare to the model described in this subsection.

4.3 Full problem

We now analyse (4.10) using the method of multiple scales. As the equilibrium point $R=\overline{R}$ remains the same as for the toy problem at leading order, we continue to use the polar coordinate system defined in (4.16) for the full problem. The presence of an $O(\unicode[STIX]{x1D716}_{1})$ term in the full equations motivates the introduction of an intermediate time scale $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D716}_{1}t$ , along with the long time scale $T=\unicode[STIX]{x1D6FF}t$ used in the previous subsection. The time derivative then becomes

(4.35) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D716}_{1}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}+\unicode[STIX]{x1D6FF}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}T},\end{eqnarray}$$

and we seek asymptotic expansions of the form

(4.36a ) $$\begin{eqnarray}\displaystyle & \displaystyle r(t,\unicode[STIX]{x1D70F},T)\sim r_{0}(t,\unicode[STIX]{x1D70F},T)+\unicode[STIX]{x1D716}_{1}r_{1}(t,\unicode[STIX]{x1D70F},T)+\unicode[STIX]{x1D6FF}r_{2}(t,\unicode[STIX]{x1D70F},T), & \displaystyle\end{eqnarray}$$

(4.36b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}(t,\unicode[STIX]{x1D70F},T)\sim \unicode[STIX]{x1D703}_{0}(t,\unicode[STIX]{x1D70F},T)+\unicode[STIX]{x1D716}_{1}\unicode[STIX]{x1D703}_{1}(t,\unicode[STIX]{x1D70F},T)+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}_{2}(t,\unicode[STIX]{x1D70F},T), & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D716}_{1},\unicode[STIX]{x1D6FF}\rightarrow 0$

$\unicode[STIX]{x1D716}_{1}\ll \unicode[STIX]{x1D716}Re\ll 1$

$\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D716}_{1}\unicode[STIX]{x1D716}Re$

$\unicode[STIX]{x1D716}_{1}^{2}\ll \unicode[STIX]{x1D6FF}\ll \unicode[STIX]{x1D716}_{1}\ll 1$

$O(\unicode[STIX]{x1D716}_{1}^{2})$

At $O(1)$ , the governing equations (4.10) yield the same equations as in § 4.2, given by (4.20). Thus, we find that

(4.37a,b ) $$\begin{eqnarray}r_{0}=r_{0}(\unicode[STIX]{x1D70F},T),\quad \unicode[STIX]{x1D703}_{0}=t+\unicode[STIX]{x1D719}_{0}(\unicode[STIX]{x1D70F},T),\end{eqnarray}$$

and the trajectory for $t=O(1)$ is thus circular motion around the steady point, as discussed in § 4.2.

At $O(\unicode[STIX]{x1D716}_{1})$ , we find that

(4.38a ) $$\begin{eqnarray}\displaystyle & \displaystyle r_{0}r_{1t}+\overline{R}(r_{1t}\cos \unicode[STIX]{x1D703}_{0}-r_{0}\unicode[STIX]{x1D703}_{1t}\sin \unicode[STIX]{x1D703}_{0})=\unicode[STIX]{x1D707}_{1}(r_{0},\unicode[STIX]{x1D703}_{0}), & \displaystyle\end{eqnarray}$$

(4.38b ) $$\begin{eqnarray}\displaystyle & \displaystyle r_{0}^{2}\unicode[STIX]{x1D703}_{1t}+\overline{R}(r_{0}\unicode[STIX]{x1D703}_{1t}\cos \unicode[STIX]{x1D703}_{0}+r_{1t}\sin \unicode[STIX]{x1D703}_{0})=\unicode[STIX]{x1D708}_{1}(r_{0},\unicode[STIX]{x1D703}_{0}), & \displaystyle\end{eqnarray}$$

(4.39a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}_{1}(r_{0},\unicode[STIX]{x1D703}_{0})=\unicode[STIX]{x1D6FE}_{1}(R_{0})r_{0}\overline{R}\sin \unicode[STIX]{x1D703}_{0}-r_{0\unicode[STIX]{x1D70F}}(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0})+\overline{R}r_{0}\unicode[STIX]{x1D703}_{0\unicode[STIX]{x1D70F}}\sin \unicode[STIX]{x1D703}_{0}, & \displaystyle\end{eqnarray}$$

(4.39b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D708}_{1}(r_{0},\unicode[STIX]{x1D703}_{0})=\unicode[STIX]{x1D6FE}_{1}(R_{0})\overline{R}(\overline{R}+r_{0}\cos \unicode[STIX]{x1D703}_{0})-r_{0}\unicode[STIX]{x1D703}_{0\unicode[STIX]{x1D70F}}(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0})-r_{0\unicode[STIX]{x1D70F}}\overline{R}\sin \unicode[STIX]{x1D703}_{0},\qquad & \displaystyle\end{eqnarray}$$

(4.40) $$\begin{eqnarray}R_{0}(t,\unicode[STIX]{x1D70F},T)=\sqrt{\overline{R}^{2}+2r_{0}(\unicode[STIX]{x1D70F},T)\overline{R}\cos \unicode[STIX]{x1D703}_{0}(t,\unicode[STIX]{x1D70F},T)+r_{0}^{2}(\unicode[STIX]{x1D70F},T)}.\end{eqnarray}$$

We note that the linear operators acting on $r_{1}$ and $\unicode[STIX]{x1D703}_{1}$ in (4.38) are the same as for those acting on $r_{2}$ and $\unicode[STIX]{x1D703}_{2}$ in the toy problem in the previous subsection. Thus, the solvability conditions at this order can be obtained by substituting $\unicode[STIX]{x1D707}_{1}$ and $\unicode[STIX]{x1D708}_{1}$ from (4.39) into (4.27) in place of their unsubscripted counterparts. This yields

(4.41a ) $$\begin{eqnarray}\displaystyle & \displaystyle 0=\int _{0}^{2\unicode[STIX]{x03C0}}(-r_{0\unicode[STIX]{x1D70F}}+\overline{R}\unicode[STIX]{x1D6FE}_{1}(R_{0})\sin \unicode[STIX]{x1D703}_{0})\,\text{d}t, & \displaystyle\end{eqnarray}$$

(4.41b ) $$\begin{eqnarray}\displaystyle & \displaystyle 0=\int _{0}^{2\unicode[STIX]{x03C0}}(-r_{0}\unicode[STIX]{x1D703}_{0\unicode[STIX]{x1D70F}}+\overline{R}\unicode[STIX]{x1D6FE}_{1}(R_{0})\cos \unicode[STIX]{x1D703}_{0})\,\text{d}t. & \displaystyle\end{eqnarray}$$

$t=\unicode[STIX]{x03C0}+\unicode[STIX]{x1D719}_{0}(\unicode[STIX]{x1D70F},T)$

(4.42a,b ) $$\begin{eqnarray}r_{0}=r_{0}(T),\quad r_{0}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}={\displaystyle \frac{\overline{R}}{2\unicode[STIX]{x03C0}}}\int _{0}^{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D6FE}_{1}(R_{0})\cos \unicode[STIX]{x1D703}_{0}\,\text{d}t.\end{eqnarray}$$

Hence, as deduced in § 4.1, there is no drift at $t=O(1/\unicode[STIX]{x1D716}_{1})$ . There is, however, a change in the angular velocity of the construct centre around the steady point. Furthermore, using (4.40) we note that $R_{0}$ is $2\unicode[STIX]{x03C0}$ -periodic in $t$ with the same phase as $\cos \unicode[STIX]{x1D703}_{0}$ , and using (4.6b ) we note that $\unicode[STIX]{x1D6FE}_{1}(R_{0})>0$ is a monotonically increasing function of $R_{0}$ . Making the integral substitution $u=\cos t$ , we thus determine that the integral in (4.42b ) is always positive. Hence, we deduce that $\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{0}/\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}$ has the same sign as $\overline{R}$ , and vanishes when $\overline{R}=0$ .

At $O(\unicode[STIX]{x1D6FF})$ , the governing equations (4.10) give

(4.43a ) $$\begin{eqnarray}\displaystyle & \displaystyle r_{0}r_{2t}+\overline{R}(r_{2t}\cos \unicode[STIX]{x1D703}_{0}-r_{0}\unicode[STIX]{x1D703}_{2t}\sin \unicode[STIX]{x1D703}_{0})=\unicode[STIX]{x1D707}_{2}(r_{0},\unicode[STIX]{x1D703}_{0}), & \displaystyle\end{eqnarray}$$

(4.43b ) $$\begin{eqnarray}\displaystyle & \displaystyle r_{0}^{2}\unicode[STIX]{x1D703}_{2t}+\overline{R}(r_{0}\unicode[STIX]{x1D703}_{2t}\cos \unicode[STIX]{x1D703}_{0}+r_{2t}\sin \unicode[STIX]{x1D703}_{0})=\unicode[STIX]{x1D708}_{2}(r_{0},\unicode[STIX]{x1D703}_{0}), & \displaystyle\end{eqnarray}$$

(4.44a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}_{2}(r_{0},\unicode[STIX]{x1D703}_{0})=R_{0}(\unicode[STIX]{x1D6FE}_{2}-\unicode[STIX]{x1D6FE}_{4})+\left({\displaystyle \frac{\bar{\unicode[STIX]{x1D70C}}}{2}}r_{0}-r_{0T}\right)(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0})+r_{0}\unicode[STIX]{x1D703}_{0T}\overline{R}\sin \unicode[STIX]{x1D703}_{0}, & \displaystyle\end{eqnarray}$$

(4.44b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D708}_{2}(r_{0},\unicode[STIX]{x1D703}_{0})=R_{0}(\unicode[STIX]{x1D6FE}_{3}+\unicode[STIX]{x1D6FE}_{5})+\left({\displaystyle \frac{\bar{\unicode[STIX]{x1D70C}}}{2}}r_{0}-r_{0T}\right)\overline{R}\sin \unicode[STIX]{x1D703}_{0}-r_{0}\unicode[STIX]{x1D703}_{0T}\left(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0}\right).\qquad & \displaystyle\end{eqnarray}$$

As we expect, the linear operators acting on $r_{2}$ and $\unicode[STIX]{x1D703}_{2}$ in (4.43) are the same as those in the previous subsection. Thus, the solvability conditions can be obtained by substituting $\unicode[STIX]{x1D707}_{2}$ and $\unicode[STIX]{x1D708}_{2}$ into (4.27) in place of their unsubscripted counterparts, to yield

(4.45a ) $$\begin{eqnarray}\displaystyle & \displaystyle r_{0T}-{\displaystyle \frac{\bar{\unicode[STIX]{x1D70C}}}{2}}r_{0}={\displaystyle \frac{1}{2\unicode[STIX]{x03C0}}}\int _{0}^{2\unicode[STIX]{x03C0}}{\displaystyle \frac{(\unicode[STIX]{x1D6FE}_{2}-\unicode[STIX]{x1D6FE}_{4})(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0})+(\unicode[STIX]{x1D6FE}_{3}+\unicode[STIX]{x1D6FE}_{5})\overline{R}\sin \unicode[STIX]{x1D703}_{0}}{R_{0}}}\,\text{d}t, & \displaystyle\end{eqnarray}$$

(4.45b ) $$\begin{eqnarray}\displaystyle & \displaystyle r_{0}\unicode[STIX]{x1D703}_{0T}={\displaystyle \frac{1}{2\unicode[STIX]{x03C0}}}\int _{0}^{2\unicode[STIX]{x03C0}}{\displaystyle \frac{(\unicode[STIX]{x1D6FE}_{4}-\unicode[STIX]{x1D6FE}_{2})\overline{R}\sin \unicode[STIX]{x1D703}_{0}+(\unicode[STIX]{x1D6FE}_{5}+\unicode[STIX]{x1D6FE}_{3})(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0})}{R_{0}}}\,\text{d}t, & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D6FE}_{2}$

$\unicode[STIX]{x1D6FE}_{3}$

$\unicode[STIX]{x1D6FE}_{4}$

$\unicode[STIX]{x1D6FE}_{5}$

$R$

$\unicode[STIX]{x1D711}$

$r$

$\unicode[STIX]{x1D703}$

$r_{0}$

$\unicode[STIX]{x1D703}_{0}$

Of the two solvability conditions in (4.45), we are most interested in (4.45a ), as this condition determines whether the construct spirals inwards or outwards. To facilitate discussion of our results, we rewrite (4.45a ) as

(4.46a ) $$\begin{eqnarray}r_{0T}={\displaystyle \frac{\bar{\unicode[STIX]{x1D70C}}}{2}}r_{0}+{\mathcal{I}},\end{eqnarray}$$

where ${\mathcal{I}}={\mathcal{I}}_{p}^{x}+{\mathcal{I}}_{p}^{y}+{\mathcal{I}}_{s}^{x}+{\mathcal{I}}_{s}^{y}$ and

(4.46b ) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\mathcal{I}}_{p}^{x}={\displaystyle \frac{1}{2\unicode[STIX]{x03C0}}}\int _{0}^{2\unicode[STIX]{x03C0}}{\displaystyle \frac{\unicode[STIX]{x1D6FE}_{2}(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0})}{R_{0}}}\,\text{d}t,\quad {\mathcal{I}}_{p}^{y}={\displaystyle \frac{1}{2\unicode[STIX]{x03C0}}}\int _{0}^{2\unicode[STIX]{x03C0}}{\displaystyle \frac{\unicode[STIX]{x1D6FE}_{3}\overline{R}\sin \unicode[STIX]{x1D703}_{0}}{R_{0}}}\,\text{d}t,\\ \displaystyle {\mathcal{I}}_{s}^{x}=-{\displaystyle \frac{1}{2\unicode[STIX]{x03C0}}}\int _{0}^{2\unicode[STIX]{x03C0}}{\displaystyle \frac{\unicode[STIX]{x1D6FE}_{4}\overline{R}\sin \unicode[STIX]{x1D703}_{0}}{R_{0}}}\,\text{d}t,\quad {\mathcal{I}}_{s}^{y}={\displaystyle \frac{1}{2\unicode[STIX]{x03C0}}}\int _{0}^{2\unicode[STIX]{x03C0}}{\displaystyle \frac{\unicode[STIX]{x1D6FE}_{5}(r_{0}+\overline{R}\cos \unicode[STIX]{x1D703}_{0})}{R_{0}}}\,\text{d}t.\end{array}\right\}\end{eqnarray}$$

We note that ${\mathcal{I}}_{p}^{x}$ and ${\mathcal{I}}_{p}^{y}$ arise from the first correction to pressure acting on the construct interface in the $\boldsymbol{e}_{X}$ and $\boldsymbol{e}_{Y}$ directions, respectively, and ${\mathcal{I}}_{s}^{x}$ and ${\mathcal{I}}_{s}^{y}$ arise from the shear stress acting on the construct interface in the $\boldsymbol{e}_{X}$ and $\boldsymbol{e}_{Y}$ directions, respectively. We solve (4.46) numerically in the next subsection to efficiently determine the long-time behaviour of the construct.

4.3.1 Evaluating the long-time drift

We determine the long-time behaviour of the construct by evaluating the solvability condition (4.46) to determine $\text{d}r_{0}/\text{d}T$ . The right-hand sides of (4.46b ) are numerically evaluated using the trapezium rule. To do this, we must evaluate the infinite sums $\unicode[STIX]{x1D6FE}_{2}$ and $\unicode[STIX]{x1D6FE}_{3}$ , defined in (4.8b,c ), and the integrals $\unicode[STIX]{x1D6FE}_{4}$ and $\unicode[STIX]{x1D6FE}_{5}$ , defined in (4.9).

The functions $\unicode[STIX]{x1D6FE}_{2}$ and $\unicode[STIX]{x1D6FE}_{3}$ arise from the first correction to the fluid pressure acting on the construct interface. They are obtained by determining $\overline{E}_{n}$ and $\overline{F}_{n}$ , from (3.20e ) and (3.20f ), respectively. This requires the knowledge of $a_{in}$ and $b_{in}$ (where $i\in \{1,2,3\}$ ) obtained by equating (3.16) and (3.17). These coefficients can be determined numerically from the calculation of $p_{0}$ , given in (3.12a ), and the functions $\unicode[STIX]{x1D6F1}$ , $\unicode[STIX]{x1D6EC}_{a}$ , and $\unicode[STIX]{x1D6EC}_{b}$ we derive in appendix A to couple the outer regions. These last three functions are calculated beforehand for a range of arguments, and then we use numerical interpolation to evaluate each function at the parameter value we require. The functions $\unicode[STIX]{x1D6FE}_{4}$ and $\unicode[STIX]{x1D6FE}_{5}$ arise from the bulk fluid shear acting on the construct interface. They are obtained by substituting (4.4), the leading-order result for the bulk flow relative to the construct movement, into the shear results given in (4.9), using (3.12a ) to evaluate $p_{0}$ . To calculate these functions, we use 10 modes to approximate each infinite series, which suffices since each term is exponentially smaller than the last. We increase this to 40 modes when evaluating an orbit where the construct passes within 0.2 of the bioreactor edge, near which the convergence is slower. We use 100 points to discretize $\unicode[STIX]{x1D6FD}\in [0,2\unicode[STIX]{x03C0}]$ and 100 points to discretize $t\in [0,2\unicode[STIX]{x03C0}]$ and evaluate (4.46). This is enough because the system is periodic in both $\unicode[STIX]{x1D6FD}$ and $t$ , providing fast convergence. We checked that there were enough modes in the infinite sums and points in the grid by observing no notable change in the results when each was increased.

We first discuss the behaviour of ${\mathcal{I}}$ , which is fundamental to understanding the long-time behaviour of the system. Then, we investigate how our system behaves in terms of experimentally controllable parameters. As predicted by the toy model in § 4.2, there is a sharp change in behaviour as a construct passes between settling and orbital motion (figure 9). The general trend is for ${\mathcal{I}}$ to increase at an approximately linear rate with $r_{0}$ in the settling regime, then to decay at a decreasing rate as $r_{0}$ increases in the orbital regime (figure 9 a). This behaviour is qualitatively the same as that in the toy problem, but without the discontinuities we saw previously. Additionally, we see from figure 9(a) that the main contribution to ${\mathcal{I}}$ is from the first-order correction to pressure ( ${\mathcal{I}}_{p}^{x}$ and ${\mathcal{I}}_{p}^{y}$ ) rather than the shear stress acting on the construct interface ( ${\mathcal{I}}_{s}^{x}$ and ${\mathcal{I}}_{s}^{y}$ ), even though both terms arise at the same asymptotic order. This does not mean that the shear stress acting on the construct at a given point in time is negligible at this order, rather that the average effect of shear stress over one orbit is negligible compared to the average effect of fluid pressure over one orbit.

Although ${\mathcal{I}}<0$ for some values of $r_{0}$ in figure 9(a), we can never have $r_{0T}<0$ for the parameter values used in figure 9. This is because choosing $\overline{R}=\overline{B}(\bar{\unicode[STIX]{x1D70C}}-1)>0$ implies that $\bar{\unicode[STIX]{x1D70C}}>1$ , and thus the $\bar{\unicode[STIX]{x1D70C}}r_{0}/2$ term in (4.46) must be larger than $r_{0}/2$ , which is larger than the magnitude of the negative values of ${\mathcal{I}}$ . We see that $r_{0T}$ is positive in figure 9(b) for a variety of values of $\bar{\unicode[STIX]{x1D70C}}$ . This is true in general, and thus our model suggests that a construct that is heavier than the surrounding fluid will always spiral out until it hits the bioreactor edge. However, we will show that if the solid part of the construct is made from material that is lighter than the surrounding fluid, stable limit cycles may exist within the bioreactor.

Figure 9. (a) The function ${\mathcal{I}}$ , defined in (4.46), and its component parts. We see that ${\mathcal{I}}$ , and hence the contribution to the long-time drift, has a marked change in behaviour at the boundary between these two regions. (b) The value of $\text{d}r_{0}/\text{d}T$ , defined in (4.46), using ${\mathcal{I}}$ from (a), and $\bar{\unicode[STIX]{x1D70C}}=1.01,1.1,1.2,1.3,1.4,1.5$ . Additionally, the boundary between the settling and orbital regimes is denoted by a black dotted line, with settling on the left and orbital motion on the right. For each subfigure, we use parameter values $\overline{R}=0.2$ , $b=0.2$ , and $K=0.05$ .

Henceforth, we consider the system behaviour under the independent variation of $\overline{B}$ and $\bar{\unicode[STIX]{x1D70C}}$ , as for the toy problem. We recall that these parameters serve as proxies for the independent variation of the main experimentally controllable parameters $\unicode[STIX]{x1D6FA}$ and $\unicode[STIX]{x1D70C}_{c}$ , the bioreactor angular velocity and the construct density, respectively. Again, $\overline{B}$ is inversely proportional to $\unicode[STIX]{x1D6FA}$ , $\bar{\unicode[STIX]{x1D70C}}$ is directly proportional to $\unicode[STIX]{x1D70C}_{c}$ , and a change in either of these experimental parameters will change the steady point of the system, since $\overline{R}=\overline{B}(\bar{\unicode[STIX]{x1D70C}}-1)$ .

Figure 10. Stability and bifurcation diagrams showing the three possible long-time behaviours (separated by the dotted blue lines) for the full problem, as discussed in the main text. (a) The long-time construct behaviour in $(\bar{\unicode[STIX]{x1D70C}},\overline{B})$ parameter space. (b) A bifurcation diagram of $r_{0}$ against $\bar{\unicode[STIX]{x1D70C}}$ , for $\overline{B}=0.1$ . Over the long time scale, $r_{0}$ will vary in the direction of the arrows. The solid black line represents a stable limit cycle, the dashed black line represents an unstable spiral, and the dash-dotted line represents an unstable limit cycle. The orbital and settling motion regimes are above and below the dotted red line, respectively. The shaded region has boundary defined by $r_{0}+\overline{B}|1-\bar{\unicode[STIX]{x1D70C}}|=1-b$ and represents the construct trajectory hitting the bioreactor edge on the first orbit. We use $b=0.2$ and $K=0.05$ .

Figure 11. A schematic of the three possible motion regimes for the full problem. In Regime Ib, the construct always tends to a stable limit cycle in orbital motion. In Regime II, the construct always spirals out until it hits the bioreactor edge. In Regime III, there exists both a stable and an unstable limit cycle, so the construct will either tend to the stable limit cycle or spiral out until it hits the bioreactor edge, depending on its initial position. The crosses denote an unstable steady point, the blue solid lines denote a stable limit cycle and the dash-dotted red lines denote an unstable limit cycle.

Calculating (4.46) as described above, we find that there are three main long-time behaviours in this system (figures 10 and 11). The first two behaviours are similar to those found in the toy problem, behaviours Ib and II. For $\overline{B}$ and $\bar{\unicode[STIX]{x1D70C}}$ small enough, there is a stable limit cycle in the system which attracts all initial positions (behaviour Ib in the toy problem). For larger $\overline{B}$ and $\bar{\unicode[STIX]{x1D70C}}$ , the construct will always spiral out towards the bioreactor edge, in this case the only steady point in the system is $r_{0}=0$ , which is always unstable (behaviour II in the toy problem). The final possible behaviour is new to the full system, and thus we define this as behaviour III. In this case, there are both unstable and stable limit cycles in the system (in addition to the ever-present unstable point at $r_{0}=0$ ), so the long-time behaviour of the system depends strongly on the initial conditions. The transition from behaviours II to III is a saddle–node (or fold) bifurcation of cycles (Strogatz Reference Strogatz2014). That is, the unstable and stable limit cycles appear in the system apropos of nothing, and the stability of the origin is unchanged. The transition from behaviours III to Ib is due to the unstable cycle falling off the physical domain.

We now consider the limit of a quickly rotating bioreactor in more detail. This will allow us to obtain some asymptotic results for the system, such as an analytic result for the parameter values at which the system stability changes. A quickly rotating bioreactor corresponds to $\overline{B}\rightarrow 0$ and thus $\overline{R}\rightarrow 0$ , from which we additionally attain $\dot{R}_{0}\rightarrow 0$ , and $\dot{\unicode[STIX]{x1D711}}_{0}\rightarrow 1$ . Thus the construct orbits the bioreactor centre with $R_{0}=r_{0}$ in this limit, and the system is in near rigid-body rotation with no flow across the construct interface relative to the construct motion. From these limiting behaviours, we find that all terms in the coupling conditions (3.16) vanish, apart from the $-(X^{2}+Y^{2})/2$ term in (3.16b ). It is a simple task to turn this remaining term into a Fourier series in $\unicode[STIX]{x1D6FD}$ on the construct interface, and subsequently determine

(4.47a-d ) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{2}=-2R_{0}\sinh ^{2}\unicode[STIX]{x1D6FC}_{2}\mathop{\sum }_{n=1}^{\infty }{\displaystyle \frac{n\tanh n(\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{1})\text{e}^{-2n\unicode[STIX]{x1D6FC}_{2}}}{\tanh n(\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{1})+12K}},\quad \unicode[STIX]{x1D6FE}_{3}=0,\quad \unicode[STIX]{x1D6FE}_{4}=0,\quad \unicode[STIX]{x1D6FE}_{5}=0.\end{eqnarray}$$

Figure 12. Comparison of numerical and asymptotic results for $b=0.2$ and $K=0.05$ in the asymptotic limit of $\overline{B}\rightarrow 0$ . (a) The solid grey line is ${\mathcal{I}}$ , obtained by resolving the solvability condition (4.46) using quadrature. The dashed black line in the asymptotic result (4.49) for ${\mathcal{I}}$ for small $R_{0}$ . (b) A stability diagram, showing the three possible behaviours of the construct for varying $\bar{\unicode[STIX]{x1D70C}}$ . The boundaries between behaviours are denoted with dotted blue lines, and the line $\bar{\unicode[STIX]{x1D70C}}=\bar{\unicode[STIX]{x1D70C}}^{\ast }$ ( ${\approx}0.606$ for the parameter values used in this figure), is obtained using the asymptotic result (4.50) for small $R_{0}$ . The solid black line represents a stable spiral, the dashed black line represents an unstable spiral, and the dash-dotted line represents an unstable limit cycle. Over the long time scale, $R_{0}$ will vary as discussed in the main text, and the arrows denote the direction of change in $R_{0}$ .

In the further limit of $R_{0}\rightarrow 0$ , we use the small $R_{0}$ results in appendix B to deduce that

(4.48) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{2}\sim -{\displaystyle \frac{R_{0}}{2\displaystyle \left(1+12K\left({\displaystyle \frac{1+b^{2}}{1-b^{2}}}\right)\right)}}+O(R_{0}^{3}).\end{eqnarray}$$

Hence, using $r_{0}=R_{0}$ , we obtain the asymptotic result

(4.49) $$\begin{eqnarray}{\mathcal{I}}\sim {\mathcal{I}}_{p}^{x}\sim -{\displaystyle \frac{R_{0}}{2\displaystyle \left(1+12K\left({\displaystyle \frac{1+b^{2}}{1-b^{2}}}\right)\right)}},\end{eqnarray}$$

which shows excellent agreement with the numerical resolution of ${\mathcal{I}}$ for $\overline{R}=0$ (figure 12 a). Substituting (4.49) into (4.46), we obtain

(4.50a ) $$\begin{eqnarray}{\displaystyle \frac{\text{d}R_{0}}{\text{d}T}}=(\bar{\unicode[STIX]{x1D70C}}-\bar{\unicode[STIX]{x1D70C}}^{\ast }){\displaystyle \frac{R_{0}}{2}}+O(R_{0}^{3})\quad \text{as }R_{0}\rightarrow 0,\end{eqnarray}$$

where

(4.50b ) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D70C}}^{\ast }={\displaystyle \frac{1}{1+12K\displaystyle \left({\displaystyle \frac{1+b^{2}}{1-b^{2}}}\right)}}.\end{eqnarray}$$

Thus, the bifurcation occurs when $\bar{\unicode[STIX]{x1D70C}}$ passes through $\bar{\unicode[STIX]{x1D70C}}^{\ast }$ , as we see in figure 12(b). As $\bar{\unicode[STIX]{x1D70C}}$ decreases through $\bar{\unicode[STIX]{x1D70C}}^{\ast }$ , $R_{0}=0$ changes from an unstable to a stable spiral, and an unstable limit cycle shoots out. Hence, there is a subcritical Hopf bifurcation at this point. We note that, in the limit $\overline{B}\rightarrow 0$ , the settling regime is annihilated, and the stable limit cycle in the orbital regime tends to and coalesces with the steady point at the origin. This explains why we have a Hopf bifurcation instead of a saddle–node bifurcation of cycles when $\overline{B}=0$ . As there is no settling regime in this limit, we do not delineate between behaviours Ia and Ib, and we refer to the behaviour where the construct tends to a steady point within the domain as behaviour I. Our asymptotic result (4.50) can also be used to determine where the boundary between behaviours II and III occurs for small $\overline{B}$ . As $\bar{\unicode[STIX]{x1D70C}}$ is decreased, the first limit cycles appear when $\overline{B}=0$ (see figure 10 a). Thus, (4.50b ) also provides a supremum for the values of $\bar{\unicode[STIX]{x1D70C}}$ for which stable limit cycles are possible.

Therefore, we conclude that the two possible long-time behaviours of the construct are to spiral out to the bioreactor edge, or towards a stable limit cycle within the bioreactor. In general, stable limit cycles can only exist when the bioreactor is rotating quickly and the tissue construct is lighter than the surrounding fluid. We are able to quantify when each of these behaviours is possible as a function of the system parameters and, moreover, we are able to quantify when the long-time behaviour of the construct depends on its initial position. Moreover, our work shows that spiralling can occur due to weak inertia and no cell growth over a time scale of $t=O((\unicode[STIX]{x1D716}_{1}\unicode[STIX]{x1D716}Re)^{-1})$ , and therefore construct spiralling may not be solely due to tissue growth as hypothesized by experimentalists. Using experimental parameters, the time scale of drift due to inertia is on the order of hours, a similar time scale as for tissue growth.