2.1. Problem description for a spherical capsule in a channel flow

An initially spherical capsule of radius $a$ flows within a long microfluidic channel having a constant square section of side $2\ell$ (figure 1). The suspending fluid and capsule liquid core are incompressible Newtonian fluids with the same kinematic viscosity $\eta$.

Figure 1.

Sketch of the model geometry showing an initially spherical capsule of radius $a$ placed in a channel with a constant square section of side $2\ell$.

The capsule liquid core is enclosed by a hyperelastic isotropic membrane. Its thickness is assumed to be negligible compared to the capsule dimension. The membrane is thus modelled as a surface devoid of bending stiffness with surface shear modulus $G_S$. The two non-dimensional governing parameters of the problem are the size ratio $a/\ell$ and the capillary number

(2.1)\begin{equation} Ca= \eta V/G_S, \end{equation}

where $V$ is the mean axial velocity of the undisturbed external Poiseuille flow.

The flow Reynolds number is assumed to be very small. We solve the Stokes equations in the external ($\beta = 1$) and internal ($\beta = 2$) fluids, together with the membrane equilibrium equation to determine the dynamics of the deformable capsule within the microchannel.

For the fluid problem, we denote by $\boldsymbol {v}^{(\beta )}$, $\boldsymbol {\sigma }^{(\beta )}$ and $p^{(\beta )}$ the velocity, stress and pressure fields in the two fluids. These parameters are non-dimensionalized using $\ell$ as characteristic length, $\ell /V$ as characteristic time, and $G_S \ell$ as characteristic force. The non-dimensional Stokes equations

(2.2a,b)\begin{equation} \boldsymbol{\nabla} p^{(\beta)} = Ca\, \nabla^2\boldsymbol{v}^{(\beta)},\quad\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}^{(\beta)}=0,\quad\beta=1,2 \end{equation}

are solved in the domain bounded by the cross-sections $S_{in}$ at the tube entrance and $S_{out}$ at the exit. These cross-sections are assumed to be both located far from the capsule. The reference frame $(O,\boldsymbol {x}, \boldsymbol {y}, \boldsymbol {z})$ is centred at each time step on the capsule centre-of-mass $O$ in the HF code, but the displacement of the capsule centre-of-mass along the tube axis $\boldsymbol {Oz}$ is computed.

The boundary conditions of the problem are as follows.

1. (i) The velocity field is assumed to be the unperturbed flow field on $S_{in}$ and $S_{out}$, i.e. the flow disturbance vanishes far from the capsule.

2. (ii) The pressure is uniform on $S_{in}$ and $S_{out}$.

3. (iii) A no-slip boundary condition is assumed at the channel wall $W$ and on the capsule membrane $M$:

   (2.3a,b)\begin{equation} \forall \boldsymbol{x}\in W, \boldsymbol{v}(\boldsymbol{x})=\boldsymbol{0};\quad \forall \boldsymbol{x} \in M,\ \boldsymbol{v}(\boldsymbol{x})=\frac{\partial \boldsymbol{u}}{\partial t}. \end{equation}

4. (iv) The normal load $\boldsymbol {n}$ on the capsule membrane $M$ is continuous, i.e. the non-dimensionalized external load per unit area $\boldsymbol {q}$ exerted by both fluids is due to the viscous traction jump

   (2.4)\begin{equation} (\boldsymbol{\sigma}^{(1)}-\boldsymbol{\sigma}^{(2)})\boldsymbol{\cdot}\boldsymbol{n} = \boldsymbol{q}, \end{equation}
   where $\boldsymbol {n}$ is the unit normal vector pointing towards the suspending fluid.

To close the problem, the external load $\boldsymbol {q}$ on the membrane is deduced from the local equilibrium equation, which, in the absence of inertia, can be written as

(2.5)\begin{equation} \boldsymbol{\nabla}_s \boldsymbol{\cdot} \boldsymbol \tau + \boldsymbol{q} = \boldsymbol 0, \end{equation}

where $\boldsymbol {\tau }$ is the non-dimensionalized Cauchy tension tensor (forces per unit arc length in the deformed plane of the membrane), and $\boldsymbol {\nabla }_s \boldsymbol {\cdot }$ is the surface divergence operator. We assume that the membrane deformation is governed by the strain-softening neo-Hookean law. The principal Cauchy tensions can then be expressed as

(2.6)\begin{equation} \tau_1 = \frac{G_S}{\lambda_1 \lambda_2}\left[ \lambda_1^2 - \frac{ 1}{(\lambda_1 \lambda_2)^2}\right]\quad (\text{likewise for } \tau_2), \end{equation}

where $\lambda _1$ and $\lambda _2$ are the principal extension ratios measuring the in-plane deformation.

2.2. Numerical procedure

The FSI problem is solved by coupling a finite element method that determines the capsule membrane mechanics with a boundary integral method that solves for the fluid flows (Walter et al. Reference Walter, Salsac, Barthès-Biesel and Le Tallec2010; Hu et al. Reference Hu, Salsac and Barthès-Biesel2012). Thanks to the latter, only the boundaries of the flow domain, i.e. the channel entrance $S_{in}$ and exit $S_{out}$, the channel wall and the capsule membrane, have to be discretized to solve the problem. The mesh of the initially spherical capsule is generated by subdividing the faces of the icosahedron (regular polyhedron with 20 triangular faces) inscribed in the sphere until reaching the desired number of triangular elements. At the last step, nodes are added at the middle of all the element edges to obtain a capsule mesh with 1280 $P_2$ triangular elements and 2562 nodes, which correspond to a characteristic mesh size $\Delta h_C= 0.075\,a$. The channel mesh of the entrance surface $S_{in}$ and exit surface $S_{out}$, and of the channel wall, is generated using Modulef (INRIA, France). The central portion of the channel, where the capsule is located, is refined. The channel mesh comprises 3768 $P_1$ triangular elements and 1905 nodes.

At time $t=0$, a spherical capsule is positioned with its centre-of-mass $O$ on the channel axis. At each time step, the in-plane stretch ratios $\lambda _1$ and $\lambda _2$ are computed from the node deformation. The elastic tension tensor $\boldsymbol \tau$ is then deduced from the values of $\lambda _1$ and $\lambda _2$. The finite element method is used to solve the weak form of the membrane equilibrium equation (2.5) and determine the external load $\boldsymbol {q}$.

Applying the boundary integral method, the velocity of the nodes on the capsule membrane reads (Pozrikidis Reference Pozrikidis1992)

(2.7)\begin{align} \boldsymbol{v}(\boldsymbol{x}) \,{=}\, \boldsymbol{v}^\infty(\boldsymbol{x}) - \frac{1}{8{\rm \pi}\mu_F} \left[ \int_M \boldsymbol{J}(\boldsymbol{r})\boldsymbol{\cdot}\boldsymbol{q} \,{\rm d}S(\boldsymbol{y}) + \int_W \boldsymbol{J}(\boldsymbol{r})\boldsymbol{\cdot} \boldsymbol{f}\, {\rm d}S(\boldsymbol{y})\,{-}\,\Delta P \int_{S_{out}} \boldsymbol{J}(\boldsymbol{r})\boldsymbol{\cdot} \boldsymbol{n}\, {\rm d}S(\boldsymbol{y})\right] \end{align}

for any $\boldsymbol {x}$ in the spatial domain when the suspending and internal fluids have the same viscosity. The vector $\boldsymbol {f}$ is the disturbance wall friction due to the capsule, $\Delta P$ is the additional pressure drop, and $\boldsymbol {r}=\boldsymbol {y}-\boldsymbol {x}$.

To update the position of the membrane nodes, the nodal displacement $\boldsymbol {u}$ is computed by integrating equation (2.3a,b) in time. The procedure is repeated until the desired non-dimensional time $VT/\ell$.

For later development, it is more convenient to work on the condensed abstract form of the system. The full-order semi-discrete FSI system to solve consists of the kinematics and the membrane equilibrium algebraic equations:

(2.8)$$\begin{gather} \{\dot{\boldsymbol{u}}\} = \{\boldsymbol{v}\}, \quad t\in [0,T], \end{gather}$$

(2.9)$$\begin{gather}\{\boldsymbol{v}\} = \varphi(\{\boldsymbol{u}\}), \end{gather}$$

where $\varphi$ is a nonlinear mapping from $\mathbb {R}^{3d}$ to $\mathbb {R}^{3d}$, and $d$ is the number of nodes on the membrane. Regarding time discretization, a Runge–Kutta Ralston scheme is used:

(2.10)$$\begin{gather} \{\hat{\boldsymbol{u}}^{n+2/3}\} = \{\boldsymbol{u}^n\} + \tfrac{2}{3}\,\Delta t \, \{\boldsymbol{v}^{n}\}, \end{gather}$$

(2.11)$$\begin{gather}\{\hat{\boldsymbol{v}}^{n+2/3}\} = \{\varphi\}(\{\hat{\boldsymbol{u}}^{n+2/3}\}), \end{gather}$$

(2.12)$$\begin{gather}\{\boldsymbol{u}^{n+1}\} = \{\boldsymbol{u}^n\} + \Delta t \, \left(\tfrac{1}{4}\{\boldsymbol{v}^{n}\}+\tfrac{3}{4}\{\hat{\boldsymbol{v}}^{n+2/3}\}\right), \end{gather}$$

(2.13)$$\begin{gather}\{\boldsymbol{v}^{n+1}\} = \{\varphi\}(\{\boldsymbol{u}^{n+1}\}), \end{gather}$$

(2.14)$$\begin{gather}\{\boldsymbol{u}^0\} = \{\boldsymbol{0}\},\ \{\boldsymbol{v}^{0}\} = \{\varphi\}(\{\boldsymbol{0}\}), \end{gather}$$

where $\Delta t>0$ is a constant time step, and $\{\boldsymbol{u}\}^n$ and $\{\boldsymbol{v}\}^n$ respectively represent the discrete membrane displacement field and the discrete membrane velocity field at discrete time $t^n=n\,\Delta t$. The initial condition is simply $\{\boldsymbol{u}\}^0=\{\boldsymbol{0}\}$.

The whole numerical scheme is subject to some Courant–Friedrichs–Lewy type stability condition on the time step (Walter et al. Reference Walter, Salsac, Barthès-Biesel and Le Tallec2010) because of its explicit nature. The numerical method is conditionally stable if the time step $\Delta t$ satisfies

(2.15)\begin{equation} \frac{V}{\ell}\,\Delta t < O\left( \frac{\Delta h_C}{\ell}\,Ca\right). \end{equation}

From the computational point of view, the resolution of (2.9) at each time step requires (i) the computation of the disturbance wall friction $\boldsymbol {f}$ at all the wall nodes, (ii) the additional pressure drop $\Delta P$, (iii) the traction jump $\boldsymbol {q}$ at the membrane nodes, and (iv) the boundary integrals for each node. The resulting numerical FOM may thus be time-consuming, depending on the membrane discretization and the number of time steps. Figure 2 shows that the evolution of the computational cost when $a/\ell =0.7$, considering the mesh discretization described above and a workstation equipped with two Intel Xeon Gold 6130 CPU (2.1 GHz) processors. A week of computation is sometimes necessary to simulate the dynamics of an initially spherical capsule in a microchannel over the non-dimensional time $VT/\ell =10$.

Figure 2.

Simulation time of the dynamics of the capsule over a non-dimensional time $Vt/\ell =10$ ($a/\ell =0.7$) according to the time step $\Delta t$. Simulations were performed on a workstation equipped with two Intel Xeon Gold 6130 CPU (2.1 GHz) processors.

For that reason, an MOR strategy is studied in this paper, in order to reduce the computational time by several orders of magnitude. ROMs try to approximate solutions of the initial problem by strongly lowering the dimensionality of the numerical model, generally using a reduced basis of suitable functions, then derive a low-order system of equations.

In the case of differential algebraic equations (DAEs) such as (2.8)–(2.9), the reduced system of equations to find should also be of DAE nature. We remark that it is often possible to reformulate DAEs as a system of ordinary differential equations (ODEs) (Ascher & Petzold Reference Ascher and Petzold1998). In the next section, we give details on the chosen ROM methodology for the particular case and context of the FSI capsule problem.

3.1. Dimensionality reduction and reduced variables for displacements and velocities

Assume first that for any $t\in [0,T]$, the discrete velocity field can be approximated accurately according to the expansion

(3.1)\begin{equation} \{\boldsymbol{v}\}(t) \approx \sum_{k=1}^K \beta_k(t)\, \{ \boldsymbol{\phi}_k\} \end{equation}

for some orthonormal modes $\{ \boldsymbol{\phi}_k\} \in \mathbb {R}^d$ and real coefficients $\beta _k(t)$. The truncation rank $K\leq d$ is, of course, expected to be far less than $d$ as expected in a general ROM methodology. From the kinematics equations, we have

(3.2)\begin{align} \{\boldsymbol{u}\}(t) &= \int_0^t \{\boldsymbol{v}\}(s)\, {\rm d}s\nonumber\\ &\approx \int_0^t \beta_k(s)\, \{ \boldsymbol{\phi}_k\} \, {\rm d}s, \end{align}

so the displacement field can be represented accurately by

(3.3)\begin{equation} \{\boldsymbol{u}\}(t) \approx \sum_{k=1}^K\alpha_k(t)\, \{ \boldsymbol{\phi}_k\}, \end{equation}

where ${\alpha _k(t)=\int _0^t \beta _k(s)\, {\rm d}s}$. The coefficients $(\alpha _k(t))_k$ and $(\beta _k(t))_k$ are called the reduced variables. For the sake of readability and mental correspondence between full-order unknowns and reduced ones, we will use the convenient notations

(3.4a,b)\begin{equation} \boldsymbol{\alpha}(t)=(\alpha_1(t),\ldots,\alpha_K(t))^{\rm T},\quad \boldsymbol{\beta}(t)=(\beta_1(t),\ldots,\beta_K(t))^{\rm T}, \end{equation}

where the superscript ${\rm T}$ denotes the transpose of the matrix. The condensed matrix forms of (3.3) and (3.1), respectively, are

(3.5a,b)\begin{equation} \{\boldsymbol{u}\}(t) \approx \boldsymbol{\mathsf{Q}}\, \boldsymbol{\alpha}(t),\quad \{\boldsymbol{v}\}(t) \approx \boldsymbol{\mathsf{Q}}\, \boldsymbol{\beta}(t), \end{equation}

where $\boldsymbol{\mathsf{Q}}=[\{ \boldsymbol{\phi}_1\} ,\ldots,\{ \boldsymbol{\phi}_K\} ]\in \mathscr {M}_{dK}$. Since the modes $\{ \boldsymbol{\phi}_k\}$ are assumed to be orthonormal (for the standard Euclidean inner product), the matrix $\boldsymbol{\mathsf{Q}}$ is a semi-orthogonal matrix, i.e. $\boldsymbol{\mathsf{Q}}^{\rm T}\boldsymbol{\mathsf{Q}}=\boldsymbol{\mathsf{I}}_K$. In particular, we have $\boldsymbol {\alpha }(t)\approx \boldsymbol{\mathsf{Q}}^{\rm T}\, \{\boldsymbol{u}\}(t)$ and $\boldsymbol {\beta }(t)=\boldsymbol{\mathsf{Q}}^{\rm T}\, \{\boldsymbol{v}\}(t)$.

Note that the modes $\{ \boldsymbol{\phi}_k\}$ and reduced variables $\boldsymbol {\alpha },\boldsymbol {\beta }$ are determined for each parameter set $(Ca, a/\ell )$, but a common value of the truncation rank $K$ is chosen for all the sets. Its practical computation will be detailed in a next subsection, as well as that of the modes $\{ \boldsymbol{\phi}_k\}$.

3.2. ROM prototype

The expressions $\{ \tilde {\boldsymbol{u}}\} (t)=\boldsymbol{\mathsf{Q}}\, \boldsymbol {\alpha }(t)$ and $\{ \tilde {\boldsymbol{v}}\} (t)=\boldsymbol{\mathsf{Q}}\, \boldsymbol {\beta }(t)$ provide low-order representations of displacement and velocity fields, respectively. We can now write equations for the reduced vectors $\boldsymbol {\alpha }(t)$ and $\boldsymbol {\beta }(t)$. In this subsection, let us consider a projection Galerkin-type approach. Let us denote by $\langle \cdot,\cdot \rangle$ the standard Euclidean scalar product in $\mathbb {R}^d$. Considering a test vector $\{ \boldsymbol{w}\}$ in $W={\rm span}(\{ \boldsymbol{\varphi}_1\} ,\ldots,\{ \boldsymbol{\varphi}_K\} )$, we look for an approximate displacement field $\{\tilde{\boldsymbol{u}}\}(t)$ solution of the projected kinematics equations

(3.6)\begin{equation} \left\langle \frac{{\rm d}}{{\rm d}t}\{ \tilde {\boldsymbol{u}}\} (t),\{ \boldsymbol{w}\} \right\rangle = \langle \{ \tilde {\boldsymbol{v}}\} (t),\{ \boldsymbol{w}\} \rangle, \quad \forall\, \{ \boldsymbol{w}\} \in W. \end{equation}

By considering each test vector $\{ \boldsymbol{w}\} =\{ \boldsymbol{\varphi}_k\}$, we get the consistent reduced kinematics equation

(3.7)\begin{equation} \dot{\boldsymbol{\alpha}} = \boldsymbol{\beta}. \end{equation}

Consider now the projected field $\{ \tilde {\boldsymbol{v}}\} (t)$ that is a solution of the system of algebraic equations (Galerkin approach)

(3.8)\begin{equation} \langle \{ \tilde {\boldsymbol{v}}\} (t),\{ \boldsymbol{w}\} \rangle = \langle \varphi(\{ \tilde {\boldsymbol{u}}\} (t)),\{ \boldsymbol{w}\} \rangle, \quad \forall\, \{ \boldsymbol{w}\} \in W. \end{equation}

Again by taking the test vector $\{ \boldsymbol{w}\}=\{\boldsymbol{\phi}_k\}$, we have

(3.9)\begin{equation} \{\boldsymbol{\phi}_k\}^{\rm T} \boldsymbol{\mathsf{Q}}\,\boldsymbol{\beta}(t) = \{\boldsymbol{\phi}_k\}^{\rm T} \varphi(\boldsymbol{\mathsf{Q}}\,\boldsymbol{\alpha}(t)). \end{equation}

Considering all $k$ in $\{1,\ldots,K\}$, since $\boldsymbol{\mathsf{Q}}=[\{ \boldsymbol{\phi}_1\} ,\ldots,\{ \boldsymbol{\phi}_K\} ]$ and $\boldsymbol{\mathsf{Q}}^{\rm T}\boldsymbol{\mathsf{Q}}=\boldsymbol{\mathsf{I}}_K$ we get

(3.10)\begin{equation} \boldsymbol{\mathsf{Q}}^{\rm T} \boldsymbol{\mathsf{Q}}\,\boldsymbol{\beta}(t) = \boldsymbol{\beta}(t) = \boldsymbol{\mathsf{Q}}^{\rm T}\,\varphi(\boldsymbol{\mathsf{Q}}\,\boldsymbol{\alpha}(t)). \end{equation}

This is in the form

(3.11)\begin{equation} \boldsymbol{\beta}(t) = \boldsymbol{\varphi}_{rz}(\boldsymbol{\alpha}(t)), \end{equation}

with the mapping $\boldsymbol{\varphi}_r:\mathbb {R}^K\rightarrow \mathbb {R}^K$ defined by $\boldsymbol{\varphi}_r(\boldsymbol {\alpha })=\boldsymbol{\mathsf{Q}}^{\rm T}\,\varphi (\boldsymbol{\mathsf{Q}}\boldsymbol {\alpha })$. We get a reduced-order algebraic equilibrium equation. Unfortunately, because of the nonlinearities in $\varphi$, the computation of $\boldsymbol{\varphi}_r(\boldsymbol {\alpha })$ requires high-dimensional $O(d)$ operations, making this approach irrelevant from the performance point of view. A possible solution to deal with the nonlinear terms would be to use, for example, empirical interpolation methods (Barrault et al. Reference Barrault, Maday, Nguyen and Patera2004), but from the algorithm and implementation point of view, this would lead to an intrusive approach with specific code developments. We here rather adopt a linearization strategy in the following sense: by differentiating (3.11) with respect to time, we get

(3.12)\begin{equation} \dot{\boldsymbol{\beta}}(t) = \frac{\partial\boldsymbol{\varphi}_r}{\partial \boldsymbol{\alpha}}\,(\boldsymbol{\alpha}(t)) \, \dot{\boldsymbol{\alpha}}(t). \end{equation}

Thanks to the reduced kinematics equation (3.7), we get

(3.13)\begin{equation} \dot{\boldsymbol{\beta}}(t) = \frac{\partial\boldsymbol{\varphi}_r}{\partial \boldsymbol{\alpha}}\,(\boldsymbol{\alpha}(t)) \, \boldsymbol{\beta}(t). \end{equation}

Since $\boldsymbol{\varphi}_r$ is hard to evaluate, it is even harder to evaluate its differential. But the differential $({\partial \boldsymbol{\varphi}_r}/{\partial \boldsymbol {\alpha }})(\boldsymbol {\alpha }(t))$ can be approximated itself, or replaced by some matrix $\boldsymbol{\mathsf{A}}(t)$. Then we get a ROM structure (ROM prototype) in the form

(3.14)$$\begin{gather} \dot{\boldsymbol{\alpha}} = \boldsymbol{\beta}(t), \end{gather}$$

(3.15)$$\begin{gather}\dot{\boldsymbol{\beta}}(t) = \boldsymbol{\mathsf{A}}(t) \, \boldsymbol{\beta}(t) . \end{gather}$$

The differential system (3.14)–(3.15) is linear with variable coefficient matrix $\boldsymbol{\mathsf{A}}(t)\in \mathscr {M}_K(\mathbb {R})$. It can be written in matrix form as

(3.16)\begin{equation} \frac{{\rm d}}{{\rm d}t}\begin{pmatrix} \boldsymbol{\alpha}(t)\\ \boldsymbol{\beta}(t) \end{pmatrix} = \underbrace{\begin{pmatrix} [0] & \boldsymbol{\mathsf{I}}_K \\ [0] & \boldsymbol{\mathsf{A}}(t) \end{pmatrix}}_{=\boldsymbol{\mathsf{A}}(t)} \begin{pmatrix} \boldsymbol{\alpha}(t)\\ \boldsymbol{\beta}(t) \end{pmatrix}. \end{equation}

The spectral properties of the differential system (3.16) are related to the spectral properties of matrix $\boldsymbol{\mathsf{A}}(t)$. In particular, if all the (complex) eigenvalues $\lambda _k(t)$ of $\boldsymbol{\mathsf{A}}(t)$ are such that $\mathrm {Re}(\lambda _k(t))<0$ for all $k$ (uniformly distributed in time), then the system is dissipative.

3.3. Non-intrusive approach, singular value decomposition and POD modes

One of the requirements of this work is to achieve a non-intrusive ROM, meaning that the effective implementation of the ROM does not involve tedious low-level development of the FOM code. For that, a data-based approach is adopted: from the FOM code, it is possible to compute FOM solutions $(\{\boldsymbol{u}\}^n,\{\boldsymbol{v}\}^n)$ at discrete times $t^n$, $n=0,\ldots,N$ ($t^N=N\,\Delta t=T$), then store some snapshot solutions (called snapshots) into a database for data analysis and later design of a ROM. Proper orthogonal decomposition (POD) (Berkooz, Holmes & Lumley Reference Berkooz, Holmes and Lumley1993) is today a well-known dimensionality reduction approach to determine the principal components from solutions of partial differential equations. The Sirovich snapshot variant approach (Sirovich Reference Sirovich1987) is based on snapshot solutions from a FOM to get a posteriori empirical POD modes $\{ \boldsymbol{\varphi}_k\}$. For the sake of simplicity, assume that the snapshot solutions are all the discrete FOM solution at simulation instants. Applying a singular value decomposition (SVD) to the displacement snapshot matrix

(3.17)\begin{equation} \mathbb{S}^u = [\boldsymbol{u}^1,\boldsymbol{u}^2,\ldots,\boldsymbol{u}^N] \end{equation}

of size $d\times N$, we get the SVD decomposition

(3.18)\begin{equation} \mathbb{S}^u = \boldsymbol{\mathsf{U}} \varSigma \boldsymbol{\mathsf{V}}^{\rm T}, \end{equation}

with orthogonal matrices $\boldsymbol{\mathsf{U}}\in \mathscr {M}_d(\mathbb {R})$, $\boldsymbol{\mathsf{V}}\in \mathscr {M}_N(\mathbb {R})$ and the singular value matrix $\varSigma ={\rm diag}(\sigma _k)\in \mathscr {M}_{d\times N}(\mathbb {R})$, with $\sigma _k\geq 0$ for all $k$ organized in decreasing order: $\sigma _1\geq \sigma _2\geq \cdots \geq \sigma _{\min (d,N)}\geq 0$. From SVD theory, for a given accuracy threshold $\varepsilon >0$, the truncation rank $K=K(\varepsilon )$ is computed as the smallest integer such that the inequality

(3.19)\begin{equation} \frac{\displaystyle{\sum_{k=K+1}^{\min(d,N)} \sigma_k^2}}{\displaystyle{\sum_{k=1}^{\min(d,N)}\sigma_k^2}}\leq \varepsilon \end{equation}

holds (Shawe-Taylor & Cristianini Reference Shawe-Taylor and Cristianini2004). Proceeding like that, it is shown that the relative orthogonal projection error of the snapshots $\{\boldsymbol{v}\}^n$ onto the linear subspace $W$ spanned by the $K$ first eigenvectors of $\boldsymbol{\mathsf{U}}$ is controlled by $\varepsilon$. Denoting by ${\rm \pi} ^W$ the linear orthogonal projection operator over $W$, we have

(3.20)\begin{equation} \sum_{n=1}^N \|\{\boldsymbol{v}\}^n-{\rm \pi}^W\{\boldsymbol{v}\}^n\|^2 \leq \varepsilon \sum_{n=1}^N \|\{\boldsymbol{v}\}^n\|^2. \end{equation}

The matrix $\boldsymbol{\mathsf{Q}}$ is obtained as the restriction of $\boldsymbol{\mathsf{U}}$ to its first $K$ columns.

3.4. Data-driven identification of coefficient matrix

The system (3.14)–(3.15) is still not closed since the coefficient matrices $\boldsymbol{\mathsf{A}}(t)$ are unknowns. From FOM data, one can try to identify the matrices by minimizing some residual function that measures the model discrepancy. The simplest linear model corresponds to the case where $\boldsymbol{\mathsf{A}}(t)$ is searched as a constant-time matrix $\boldsymbol{\mathsf{A}}$. In this case, (3.15) becomes $\dot {\boldsymbol {\beta }}(t) = \boldsymbol{\mathsf{A}} \, \boldsymbol {\beta }(t)$. This is the scope of this paper. From the continuous-time problem, one could determine the matrix $\boldsymbol{\mathsf{A}}$ by minimizing the least squares functional

(3.21)\begin{equation} \min_{\boldsymbol{\mathsf{A}}\in\mathscr{M}_K(\mathbb{R})} \frac{1}{2}\int_0^{\rm T} \|\dot{\boldsymbol{\beta}}(t)-\boldsymbol{\mathsf{A}}\,\boldsymbol{\beta}(t)\|^2\, {\rm d}t. \end{equation}

But practically, we have velocity snapshot data only at discrete times, and we do not have access to the time derivatives of the velocity fields. So the following numerical procedure is adopted. From the velocity snapshot matrix $\mathbb {S}^v=[\{\boldsymbol{v}\}^1,\ldots,\{\boldsymbol{v}\}^N]$, we compute first the reduced snapshot variables

(3.22)\begin{equation} \boldsymbol{\beta}^n = \boldsymbol{\mathsf{Q}}\, \{\boldsymbol{v}\}^n,\quad n=1,\ldots,N. \end{equation}

Next, we determine a matrix $\boldsymbol{\mathsf{A}}$ that minimizes the least squares cost function

(3.23)\begin{equation} \min_{\boldsymbol{\mathsf{A}}\in\mathscr{M}_K(\mathbb{R})}\ \frac{1}{2} \sum_{n=1}^{N-1} \left\|\frac{\boldsymbol{\beta}^{n+1}-\boldsymbol{\beta}^n}{\Delta t}-\boldsymbol{\mathsf{A}}\boldsymbol{\beta}^n \right\|^2. \end{equation}

In (3.23), the finite difference $({\boldsymbol {\beta }^{n+1}-\boldsymbol {\beta }^n})/{\Delta t}$ is a first-order approximation (in $\Delta t$) of $\dot {\boldsymbol {\beta }}$ at time $t^n$. In Appendix A, we provide a mathematical analysis of the effect of time discretization in (3.23) about the impact on the stability of the resulting identified differential system compared to the initial one.

The minimization problem (3.23) can be written in condensed matrix form

(3.24)\begin{equation} \min_{\boldsymbol{\mathsf{A}}\in\mathscr{M}_K(\mathbb{R})}\, \frac{1}{2}\,\|\mathbb{Y}-\boldsymbol{\mathsf{A}}\mathbb{X}\|_F^2, \end{equation}

with the two data matrices

(3.25a,b)\begin{equation} \mathbb{X} = [\boldsymbol{\beta}^1,\boldsymbol{\beta}^2,\ldots,\boldsymbol{\beta}^{N-1}],\quad \mathbb{Y} = \left[\frac{\boldsymbol{\beta}^2-\boldsymbol{\beta}^1}{\Delta t},\ldots,\frac{\boldsymbol{\beta}^N-\boldsymbol{\beta}^{N-1}}{\Delta t} \right]. \end{equation}

Because $\mathbb {X}$ and $\mathbb {Y}$ store reduced variables (of size $K$), for a sufficient number of discrete snapshot times, these two matrices are horizontal ones. We will assume that the rank of $\mathbb {X}$ is always maximal, i.e. equal to $K$. The least squares solution $\boldsymbol{\mathsf{A}}$ of (3.24) is then given by

(3.26)\begin{equation} \boldsymbol{\mathsf{A}} = \mathbb{Y}\mathbb{X}^{\dagger}, \end{equation}

where $\mathbb {X}^{\dagger} =\mathbb {X}^{\rm T}(\mathbb {X}\mathbb {X}^{\rm T})^{-1}$ is the Moore–Penrose pseudo-inverse matrix of $\mathbb {X}$. This least squares approach has close connections with SVD-based DMD (Schmid Reference Schmid2010; Kutz et al. Reference Kutz, Brunton, Brunton and Proctor2016).

3.5. Tikhonov least squares regularized formulation

From standard spectral theory arguments, it is expected that the POD coefficients decay rapidly when $k$ increases as soon as both displacement and velocity fields are smooth enough. A possible side effect is the bad condition number of the matrix $\mathbb {X}$, since the last rows of $\mathbb {X}$ have small coefficients (thus leading to row vectors close to zero ‘at the scale’ of the first row of $\mathbb {X}$). Even if the solution $\boldsymbol{\mathsf{A}}$ in (3.26) always exists, the solution may be sensitive to perturbations, noise or round-off errors. In order to get a robust identification approach, one can regularize the least squares problem (3.24) by adding a Tikhonov regularization term (see e.g. Aster, Borchers & Thurber Reference Aster, Borchers and Thurber2019)

(3.27)\begin{equation} \min_{\boldsymbol{\mathsf{A}}\in\mathscr{M}_K(\mathbb{R})}\, \frac{1}{2}\,\|\mathbb{Y}-\boldsymbol{\mathsf{A}}\mathbb{X}\|_F^2 + \frac{\mu}{2}\, \|\mathbb{X}\|_F^2\, \|\boldsymbol{\mathsf{A}}\|_F^2, \end{equation}

where the scalar $\mu >0$ is the regularization coefficient. The factor $\|\mathbb {X}\|_F^2$ in the regularization term has been added for scaling purposes. The solution $\boldsymbol{\mathsf{A}}_\mu$ of (3.27) is given by

(3.28)\begin{equation} \boldsymbol{\mathsf{A}}_\mu = \mathbb{Y} \mathbb{X}^{\rm T}(\mathbb{X}\mathbb{X}^{\rm T}+\mu \|\mathbb{X}\|_F^2\, \boldsymbol{\mathsf{I}}_K)^{{-}1}. \end{equation}

3.6. Reduced-order continuous dynamical system

Once the matrix $\boldsymbol{\mathsf{A}}_\mu$ has been determined, we get the reduced-order continuous dynamical system

(3.29)$$\begin{gather} \dot{\boldsymbol{\alpha}} = \boldsymbol{\beta}, \end{gather}$$

(3.30)$$\begin{gather}\dot{\boldsymbol{\beta}} = \boldsymbol{\mathsf{A}}_\mu\, \boldsymbol{\beta} \end{gather}$$

with initial conditions $\boldsymbol {\alpha }(0)=\boldsymbol {0}$, $\boldsymbol {v}(0)=\boldsymbol{\mathsf{Q}}^{\rm T}\,\varphi (\{\boldsymbol{0}\} )$. At any time $t$, one can go back to the high-dimensional physical space using the POD modes: $\{\boldsymbol{u}\}(t)=\boldsymbol{\mathsf{Q}}\,\boldsymbol {\alpha }(t)$, $\{\boldsymbol{x}\}(t)=\{ \boldsymbol{X}\} +\{\boldsymbol{u}\}(t)$, $\{\boldsymbol{v}\}(t)=\boldsymbol{\mathsf{Q}}\,\boldsymbol {\beta }(t)$. As already mentioned, the system can be written in condensed matrix form

(3.31)\begin{equation} \dot{\boldsymbol{w}} = \mathbb{A}_\mu\, \boldsymbol{w}, \end{equation}

where

(3.32a,b)\begin{equation} \boldsymbol{w}(t)=(\boldsymbol{\alpha}(t),\boldsymbol{\beta}(t))^{\rm T}\quad {\rm and}\quad \mathbb{A}_\mu = \begin{pmatrix}[0]_K & \boldsymbol{\mathsf{I}}_K\\ [0]_K & \boldsymbol{\mathsf{A}}_\mu. \end{pmatrix}. \end{equation}

The exact analytical solution of (3.31) is

(3.33)\begin{equation} \boldsymbol{w}(t) = \exp(\mathbb{A}_\mu t)\,\boldsymbol{w}(0). \end{equation}

The stability of the differential system depends on the spectral structure of $\mathbb {A}_\mu$, or equivalently on the spectrum of $\boldsymbol{\mathsf{A}}_\mu$. Because of the stability of the fluid–capsule coupled system, and from accurate solutions of the FOM solver, one can hope that the solution $\boldsymbol{\mathsf{A}}_\mu$ of the least squares identification problem has the expected spectral properties. This will be studied and discussed in the numerical experimentation section. From the kinetic energy point of view, it is shown in Appendix B that the stability of the kinetic energy is linked to the property of the (real) spectrum of the symmetric matrix $(\boldsymbol{\mathsf{A}}_\mu +\boldsymbol{\mathsf{A}}_\mu ^{\rm T})/2$.

3.7. Reduced-order discrete dynamical system

Of course, it is also possible to derive a discrete dynamical system from the continuous one by using a standard time advance scheme. For example, the explicit forward Euler scheme with a constant time step $\Delta t$ gives

(3.34)$$\begin{gather} \boldsymbol{\alpha}^{n+1} = \boldsymbol{\alpha}^n + \Delta t\, \boldsymbol{\beta}^n, \end{gather}$$

(3.35)$$\begin{gather}\boldsymbol{\beta}^{n+1} = \boldsymbol{\beta}^n + \Delta t\, \boldsymbol{\mathsf{A}}_\mu \boldsymbol{\beta}^n. \end{gather}$$

By multiplying (3.34) by $\boldsymbol{\mathsf{Q}}$, we get the space–time approximate solution

(3.36)\begin{equation} \{\boldsymbol{u}\}^{n+1} = \{\boldsymbol{u}\}^n + \Delta t\, \{\boldsymbol{v}\}^n, \end{equation}

so the ROM is completely consistent with the kinematics equation. Stability properties of the discrete system are linked to the spectral properties of the matrix

(3.37)\begin{equation} \boldsymbol{\mathsf{A}}_\mu^\Delta = \begin{pmatrix} \boldsymbol{\mathsf{I}}_K & \Delta t\, \boldsymbol{\mathsf{I}}_K \\ [0]_K & (\boldsymbol{\mathsf{I}}_K+\Delta t\, \boldsymbol{\mathsf{A}}_\mu) \end{pmatrix}. \end{equation}

For unconditional stability in time, it is required for the eigenvalues of $\boldsymbol{\mathsf{I}}_K+\Delta t\,\boldsymbol{\mathsf{A}}_\mu$ to stay in the unit disk of the complex plane.

More generally, it is possible to use any other time advance scheme, according to the expected order of accuracy or stability domain.

3.8. Accuracy criteria and similarity distances between ROM and FOM solutions

In order to quantify the error induced by approximations, we introduce three accuracy criteria. The first accuracy criterion is the relative information content (RIC), defined by

(3.38)\begin{equation} {RIC}(K) = \frac{\displaystyle{\sum_{k=K+1}^{\min(d,N)} \sigma_k^2}}{\displaystyle{\sum_{k=1}^{\min(d,N)}\sigma_k^2}}, \end{equation}

which quantifies the relative amount of neglected information when truncating the number of modes at rank $K$. The truncation rank is determined such that the RIC is inferior to the accuracy threshold $\varepsilon$. The accuracy threshold $\varepsilon$ is fixed at $10^{-6}$.

The second accuracy criterion is the relative time residual $\mathcal {R}$. It quantifies the relative error induced by the minimization of the least squares cost function (3.23) using $\boldsymbol{\mathsf{A}}_{\mu }$, and is given by

(3.39)\begin{equation} \mathcal{R}(j)=\frac{\Vert \boldsymbol{\mathsf{A}}_{\mu} \mathbb{X}_j-\mathbb{Y}_j \Vert^2}{\Vert \mathbb{Y}_j \Vert^2}, \end{equation}

where $\mathbb {X}_j$ represents the $j$th column of $\mathbb {X}$, and $\mathbb {Y}_j$ the $j$th column of $\mathbb {Y}$. The index $j$ is thus linked to the snapshots ($j\in \{1,\ldots, N\}$). To better draw a parallel between the evolution of this parameter and the capsule dynamics, this parameter will be represented as a function of the non-dimensional time $Vt/\ell$ hereafter.

The third accuracy criterion, $\varepsilon _{{Shape}}(Vt/\ell )$, measures the difference between the 3-D reference capsule shape given by the FOM ($\mathcal {S}_{{FOM}}$) and the 3-D shape predicted by the ROM ($\mathcal {S}_{{ROM}}$). It is defined at a given non-dimensional time $Vt/\ell$ as the ratio between the modified Hausdorff distance (MHD) computed between $\mathcal {S}_{FOM}$ and $\mathcal {S}_{ROM}$, and non-dimensionalized by $\ell$:

(3.40)\begin{equation} \varepsilon_{{Shape}}(Vt/\ell) =\frac{{MHD}(\mathcal{S}_{{FOM}}(Vt/\ell), \mathcal{S}_{{ROM}}(Vt/\ell))}{\ell}. \end{equation}

The MHD is the maximum value of the mean distance between $\mathcal {S}_{{FOM}}$ and $\mathcal {S}_{{ROM}}$, and the mean distance between $\mathcal {S}_{{ROM}}$ and $\mathcal {S}_{{FOM}}$ (Dubuisson & Jain Reference Dubuisson and Jain1994).

4.1. Proper orthogonal decomposition, truncation and modes

The SVD is first applied to the displacement snapshot matrix. To determine the truncation rank, the evolution of $1- RIC$ is illustrated in figure 4 as a function of number of modes considered. The RIC is approximately 1 % only with one mode. The more modes kept, the less information is neglected. In the following, we fix the number of modes to 20. The accuracy threshold $\varepsilon$ is thus equal to $10^{-6}$.

Figure 4.

Evolution of the relative amount of neglected information $1-RIC$, as a function of the number of modes ($Ca=0.13$, $a/\ell =0.8$).

The modes are determined from the displacement snapshot matrix. They are added to the sphere of radius 1 and amplified by a factor 2 to be visualized. The first six modes are represented in figure 5 for $Ca=0.13$, $a/\ell =0.8$.

Figure 5.

Representation of the first six modes of the capsule dynamics when $a/\ell =0.80$ and $Ca=0.13$. To be visualized, the modes of displacement were added to the sphere of radius 1 and amplified by a factor 2.

The first six modes are mostly dedicated to changing the shape of the rear of the capsule. The following modes appear to become noisy (not shown). However, these modes are not negligible if one wants to get an accuracy of $10^{-6}$.

4.2. Dynamic mode decomposition: empirical regularization

Before determining the matrix $\boldsymbol{\mathsf{A}}$, we check the condition numbers of the matrices $\mathbb {X}$ and $\mathbb {XX}^{\rm T}$. They are respectively equal to $6.5\times 10^4$ and $4.3 \times 10^9$. The condition numbers of these matrices are very high, and the matrix $\boldsymbol{\mathsf{A}}$, determined by solving (3.26), may be sensitive to perturbations or noise. To improve the robustness, a Tikhonov regularization is applied to solve the least squares problem (3.23), and the matrix $\boldsymbol{\mathsf{A}}_{\mu }$ is computed using (3.28), which depends on the regularization coefficient $\mu$. To determine the optimal value of $\mu$, the relative least squares error $\Vert \boldsymbol{\mathsf{A}}_{\mu }\mathbb {X}-\mathbb {Y}\Vert _F/\Vert \mathbb {Y}\Vert _F$ is represented according to the norm solution $\Vert \boldsymbol{\mathsf{A}}_{\mu }\Vert _F$ when 20 modes are considered and when $\mu$ is varied between $10^{-12}$ and $10^{-3}$ (figure 6). The least squares error $\Vert \boldsymbol{\mathsf{A}}_{\mu }\mathbb {X}-\mathbb {Y}\Vert _F$ and the norm solution $\Vert \boldsymbol{\mathsf{A}}_{\mu }\Vert _F$ are minimal when $\mu =10^{-9}$. In the following, $\mu$ is thus fixed to $\mu =10^{-9}$ and the number of modes to 20.

Figure 6.

Evolution of the norm solution $\Vert \boldsymbol{\mathsf{A}}_{\mu }\Vert _F$ as a function of the least squares error $\Vert \boldsymbol{\mathsf{A}}_{\mu }\mathbb {X}-\mathbb {Y}\Vert _F / \Vert \mathbb {Y}\Vert _F$ when the number of modes is fixed to 20 and $(Ca=0.13, a/\ell =0.8)$.

4.3. Validity check of the ROM: spectral study of the resulting matrix

In order to detect anomalies, a spectral analysis of the ROM learned by the DMD method is carried out. The spectrum of the matrix $\boldsymbol{\mathsf{A}}_{\mu }$ is represented in figure 7. All the eigenvalues $\lambda _k$ of the matrix $\boldsymbol{\mathsf{A}}_{\mu }$ have non-positive real parts. The resulting linear ROM is thus stable.

Figure 7.

Eigenvalues $\lambda _k$ of $\boldsymbol{\mathsf{A}}_{\mu }$ ($Ca=0.13$, $a/\ell =0.8$, 20 modes, $\mu =10^{-9}$). Note that the maximum real part of the eigenvalues is exactly equal to zero.

The temporal evolution of the residual $\mathcal {R}$ (figure 8) shows that the error is less than 0.7 %. The maximal value is reached at the beginning of the simulation ($Vt/\ell < 0.3$), and $\mathcal {R}$ decreases afterwards. When $Vt/\ell \lesssim 6$, i.e. before the capsule has reached its steady state, high-frequency oscillations are observed. This probably means that a high-frequency mode is neglected, even if 20 modes are considered. For $Vt/\ell >6$, $\mathcal{R}$ is of order $10^{-7}$. The stationary state is thus well predicted by the model, and the error during the transient stage is more than acceptable.

Figure 8.

Temporal evolution of the normalized time residual with $Ca=0.13$, $a/\ell =0.8$, 20 modes and $\mu =10^{-9}$.

4.4. ROM online stage and accuracy assessment

The displacement of all the nodes of the capsule mesh estimated by the ROM is then added to the corresponding node of the sphere of radius 1 to visualize the temporal evolution of the capsule shape in three dimensions. Figure 9 shows the capsule dynamics for the reference case ($Ca=0.13$, $a/\ell =0.8$). The ROM allows us to reproduce the capsule deformation from the initial state up to the parachute-shaped steady state. For the FOM and the ROM, the capsule profile is then determined in the cutting plane passing through the middle of the microchannel. Figure 10 shows that the two capsule profiles overlap perfectly at $Vt/\ell =0, 0.4, 2, 4, 6$. The temporal evolution of $\varepsilon _{{Shape}}$ is shown in figure 11(a). The maximum value of the error committed on the 3-D shape $\varepsilon _{{Shape}}$ equals 0.004 %. The error on the capsule shape $\varepsilon _{{Shape}}$ is thus negligible. The deformation of the capsule from its initially spherical shape to its steady state over a non-dimensional time $Vt/\ell =10$ can thus be estimated very precisely with the developed ROM.

Figure 9.

Dynamics of a microcapsule flowing in microchannel with a square cross-section predicted by the ROM at the non-dimensional times $Vt/\ell = 0.4$, 2.8, 5.2, 7.6, 10, with $Ca=0.13$, $a/\ell =0.8$, 20 modes and $\mu =10^{-9}$. The initial spherical capsule is shown on the left by transparency.

Figure 10.

Comparison of the capsule contours given by the FOM (dotted line) and estimated by the ROM (orange line). The capsule is shown for $Ca=0.13$, $a/\ell =0.8$ at the non-dimensional times $Vt/\ell = 0$, 0.4, 2, 4, 6. The horizontal lines represent the channel borders. The number of modes is fixed at 20, and $\mu =10^{-9}$.

Figure 11.

(a) Comparison of the capsule contours given by the FOM (dotted line) and estimated by the ROM (orange line) for the different learning times $VT_{L}/\ell$. (b) Evolution of $\varepsilon _{{Shape}}$ measured at $Vt/\ell =10$ as a function of the learning time $VT_{L}/\ell$. (c) Influence of the learning time $VT_{L}/\ell$ on the temporal evolution of the error on the capsule shape $\varepsilon _{{Shape}}$. The error during the learning time is shown with a solid line. For this case, the parameters are 20 modes, $\mu =10^{-9}$, $Ca = 0.13$ and $a/\ell =0.8$.

The DMD method predicts the capsule displacement at time $t^{n+1}$ from that at time $t^n$. The model has been constructed until now by considering the dynamics of the capsule over a non-dimensional time $Vt/\ell =10$.

In order to study the sensitivity of the ROM to the learning time $VT_{L}/\ell$, i.e. the non-dimensional time over which the model is trained, we modify it with values between 2 and 8, knowing that the time to reach the steady state is in this case $VT_{SS}/\ell =6.2$. We estimate the capsule dynamics using the ROM up to a non-dimensional time $Vt/\ell =10$. The number of modes is always equal to 20, and $\mu =10^{-9}$. The comparison of the estimated shape at $Vt/\ell =10$ with the one simulated with the FOM (figure 11a) shows that $VT_{L}/\ell \geq 4$ is sufficient to predict very well the capsule dynamics. Figure 11(b) confirms that the error on the capsule shape is negligible when $VT_{L}/\ell \geq 4$. It is interesting that the ROM could predict the steady state even when $T_L < T_{SS}$. This could be due to the fact that the maximum real part of the eigenvalues is exactly equal to zero from $VT_{L}/\ell \geq 4$. The maximum real part of the eigenvalues is negative and close to zero for $VT_{L}/\ell =2$. Figure 11(c) shows that in the worst case ($VT_{L}/\ell =2$), the error on the capsule shape increases with time but reaches a plateau from time $Vt/\ell =8$. The system is thus stable without exponential drift, as proven by the negative value of the maximum real eigenvalue. In the inset, the error also increases for the other learning times but remains very small (below 0.2 %).

6.1. General methodology

It is here again assumed that a training database of $N$ pre-computed FOM results is available. Now we would like to derive a ROM for any parameter couple $\boldsymbol {\theta }=(Ca,a/\ell )$ in the admissible parameter domain. The proposed space–time-parameter ROM is made of two steps. The first step consists in predicting the space–time solution $\{\boldsymbol {u}\}(t;\boldsymbol {\theta })$ by means of a robust interpolation procedure. The second step consists in deriving a ROM in the form of a low-order dynamical system by using the predicted solutions of the first step as training data. Then we apply the former procedure detailed in § 3. Below, we give a detailed explanation of the two steps.

6.1.1. Step 1: predictor step

Considering a parameter couple $\boldsymbol {\theta }$, we first search the three nearest neighbour parameters in the sample set that form a non-degenerate triangle in the plane $(Ca,a/\ell )$. Let us denote them by $\boldsymbol {\theta }_1$, $\boldsymbol {\theta }_2$ and $\boldsymbol {\theta }_3$. We will define a linear operator in the triangle $(\boldsymbol {\theta }_1,\boldsymbol {\theta }_2,\boldsymbol {\theta }_3)$. For that, let us introduce the barycentric coordinates $(\lambda _1,\lambda _2,\lambda _3)$, $\lambda \in [0,1]$, $i=1,2,3$, such that

(6.1)$$\begin{gather} \lambda_1 + \lambda_2 + \lambda_3 = 1, \end{gather}$$

(6.2)$$\begin{gather}\boldsymbol{\theta}_1 \lambda_1 + \boldsymbol{\theta}_2 \lambda_2 + \boldsymbol{\theta}_3\lambda_3 = \boldsymbol{\theta}. \end{gather}$$

The $3\times 3$ linear system (6.1) and (6.2) is invertible as soon as the triangle $(\boldsymbol {\theta }_1,\boldsymbol {\theta }_2,\boldsymbol {\theta }_3)$ is non-degenerate. Notice that the $\lambda _i$ ($i=1,2,3$) are actually functions of $\boldsymbol {\theta }$. Let us now denote by $\{\boldsymbol{u}_1\}$, $\{\boldsymbol{u}_2\}$ and $\{\boldsymbol{u}_3\}$ the displacement fields for the parameter vectors $\boldsymbol {\theta }_1$, $\boldsymbol {\theta }_2$ and $\boldsymbol {\theta }_3$, respectively. Then we can consider the predicted velocity field $\hat {\boldsymbol {u}}(t;\boldsymbol {\theta })$ defined by

(6.3)\begin{equation} \{\hat{\boldsymbol{u}}\}(t,\boldsymbol{\theta}) = \lambda_1 \{\boldsymbol{u}_1\}(t) + \lambda_2 \{\boldsymbol{u}_2\}(t) + \lambda_3 \{\boldsymbol{u}_3\}(t). \end{equation}

6.1.2. Step 2: low-order dynamical system ROM

Expression (6.3) can be evaluated at some discrete instants in order to generate new training data. Then the SVD-DMD ROM methodology presented in § 3 can be applied to these data to get a reduced dynamical system in the form

(6.4)$$\begin{gather} \dot{\boldsymbol{\alpha}}(\boldsymbol{\theta}) = \boldsymbol{\beta}(\boldsymbol{\theta}), \end{gather}$$

(6.5)$$\begin{gather}\dot{\boldsymbol{\beta}}(\boldsymbol{\theta}) = \boldsymbol{\mathsf{A}}_\mu(\boldsymbol{\theta})\, \boldsymbol{\beta}(\boldsymbol{\theta}). \end{gather}$$

We also have a matrix $\boldsymbol{\mathsf{Q}}(\boldsymbol {\theta })$ of orthogonal POD modes, and we can go back to the high-dimensional physical space by the standard operations

(6.6a,b)\begin{equation} \{\hat{\boldsymbol{u}}\}(t,\boldsymbol{\theta}) \approx \boldsymbol{\mathsf{Q}}(\boldsymbol{\theta})\, \boldsymbol{\alpha}(t,\boldsymbol{\theta}),\quad \{\hat{\boldsymbol{v}}\}(t,\boldsymbol{\theta}) \approx \boldsymbol{\mathsf{Q}}(\boldsymbol{\theta})\, \boldsymbol{\beta}(t,\boldsymbol{\theta}). \end{equation}

Notice that the capsule position field $\{\boldsymbol {x}\}(t,\boldsymbol {\theta })$ is given by

(6.7)\begin{equation} \{ \boldsymbol{x}\}(t;\boldsymbol{\theta}) = \{\boldsymbol{X}\}(\boldsymbol{\theta}) + \{\hat {\boldsymbol{u}}\}(t,\boldsymbol{\theta}), \end{equation}

with an initial capsule position $\{\boldsymbol{X}\}(\boldsymbol {\theta })$ that may depend on $\boldsymbol {\theta }$ because of the pre-deformation pre-processing if $a/\ell \geq 0.90$.

A testing database is created using the FOM as in § 5 and considering $(Ca,a/\ell )$ couples that are not in the training database. A set of 110 $(Ca,a/\ell )$ couples are included in this database (figure 15). For all the $(Ca,a/\ell )$ couples of the testing database, the capsule dynamics is interpolated from the dynamics of the three closest neighbours at a given non-dimensional time. Capsule shapes obtained by the ROM are compared to the ones predicted by the FOM at the same non-dimensional time. Figure 16 represents the evolution of the error committed on the capsule shape $\varepsilon _{Shape}$ on the training database at $Vt/\ell =0, 0.4, 1, 2, 5, 10$. At the initial time, $\varepsilon _{{Shape}}$ is zero. The interpolation method is therefore able to capture the initial capsule shape. When the time increases, $\varepsilon _{{Shape}}$ increases and is greater than if we apply directly the POD-DMD method on the FOM results and reconstruct the dynamics. However, $\varepsilon _{{Shape}}$ remains less than 0.3 % on the majority of the testing database. It remains fully acceptable. Here, $\varepsilon _{{Shape}}$ is more important near the steady-state limit and when we approach the lowest values of $Ca$ because we are close to the limits of the training base.

Figure 15.

Values of $Ca$ and $a/\ell$ included in the testing database (open circles). The filled squares represent the cases in the training database. The dashed line delimits the domain where a steady-state capsule deformation exists for capsules following the neo-Hookean law.

Figure 16.

Heat maps of $\varepsilon _{{Shape}}$ on the testing database as functions of $Ca$ and $a/\ell$ at $Vt/\ell$ values (a) 0, (b) 0.4, (c) 1, (d) 2, (e) 5, ( f) 10. The dashed line delimits the domain for which a steady-state capsule deformation exists.