*Article* **A Data-Driven Space-Time-Parameter Reduced-Order Model with Manifold Learning for Coupled Problems: Application to Deformable Capsules Flowing in Microchannels**

**Toufik Boubehziz 1, Carlos Quesada-Granja 1, Claire Dupont 1, Pierre Villon 2, Florian De Vuyst 3,\* and Anne-Virginie Salsac <sup>1</sup>**


**Abstract:** An innovative data-driven model-order reduction technique is proposed to model dilute micrometric or nanometric suspensions of microcapsules, i.e., microdrops protected in a thin hyperelastic membrane, which are used in Healthcare as innovative drug vehicles. We consider a microcapsule flowing in a similar-size microfluidic channel and vary systematically the governing parameter, namely the capillary number, ratio of the viscous to elastic forces, and the confinement ratio, ratio of the capsule to tube size. The resulting space-time-parameter problem is solved using two global POD reduced bases, determined in the offline stage for the space and parameter variables, respectively. A suitable low-order spatial reduced basis is then computed in the online stage for any new parameter instance. The time evolution of the capsule dynamics is achieved by identifying the nonlinear low-order manifold of the reduced variables; for that, a point cloud of reduced data is computed and a diffuse approximation method is used. Numerical comparisons between the full-order fluid-structure interaction model and the reduced-order one confirm both accuracy and stability of the reduction technique over the whole admissible parameter domain. We believe that such an approach can be applied to a broad range of coupled problems especially involving quasistatic models of structural mechanics.

**Keywords:** data-driven model; model order reduction; proper orthogonal decomposition; manifold learning; diffuse approximation; microcapsule suspension; Hausdorff distance

#### **1. Introduction**

Numerical modeling and simulation today appear to be an indispensable science to analyze physics-coupled problems (e.g., micrometric and nanometric suspensions), but also for innovative design and optimization of complex three-dimensional systems in engineering and industry (health, automotive, aircraft, etc.). Although one can nowadays find robust and accurate open-source or commercial codes for the simulation of multiphysics systems, it is still hard to use them in the context of robust or optimal design because of the prohibitive computational time that does not match with engineering production horizons. In order to accelerate computations, one can make use of parallel HPC (High Performance Computing) facilities, but this can become financially costly. For most applications, even with HPC facilities, the evaluation of the solutions takes days or weeks for three-dimensional multi-coupled problems.

Alternatively, a current tendency is to use machine learning or artificial intelligence tools to capitalize knowledge stored into data and use it for future case studies. It leads to

#### **Citation:** Boubehziz, T.;

Quesada-Granja, C.; Dupont, C.; Villon, P.; De Vuyst, F.; Salsac, A.-V. A Data-Driven Space-Time-Parameter Reduced-Order Model with Manifold Learning for Coupled Problems: Application to Deformable Capsules Flowing in Microchannels. *Entropy* **2021**, *23*, 1193. https://doi.org/ 10.3390/e23091193

Academic Editors: Amine Ammar, Francisco Chinesta and Rudy Valette

Received: 29 July 2021 Accepted: 7 September 2021 Published: 9 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

less redundant or useless computations, while the database of results continues to grow with more and more relevant information contents. With machine learning, one can expect to explore the design spaces in an easier, faster and more efficient way. However, one usually needs expertise to design and train artificial neural networks (ANN) correctly. For high-dimensional data, the training stage may require large computational resources and issues of quality of data may also be raised. Machine learning can be extended to time dependent problems and dynamical system using e.g., recurrent neural networks [1].

However, for particular use cases like physical systems, machine learning algorithms are designed to return three-dimensional spatial fields. This means that the outputs of the networks are high-dimensional vectors, which may induce training convergence and accuracy issues. Over the last two years, one could observe the rise of so-called Physicsinformed neural networks (PINN) where the artificial neural networks are trained from a loss function that includes physical information (like partial differential equations), see e.g., [2].

Another class of 'machine learning' methods are the data-driven model order reduction (MOR) techniques, that use data generated from a (time-consuming) high-fidelity solver, also called the full-order model (FOM) [3,4]. Data-driven and non-intrusive reducedorder models (ROM) can be seen as supervised ANN [5,6]. For parametrized partial differential problems, ROMs usually perform a dimensionality reduction by means of suitable reduced bases. This can be achieved via different approaches such as the Proper Orthogonal Decomposition (POD) [7–9], piecewise tangential interpolation [10], Proper Generalized Decompositions (PGD) [11–13], Empirical Interpolation Methods (EIM) [14–16] or via different greedy procedures [17]. Then one has to find the manifold that maps the parameters to the coefficients of linear combination of the reduced basis functions. This can also be achieved in a supervised way, by means of universal approximation techniques like diffuse approximation [18] for example. Using ROMs may lead to substantial speedups as compared to FOM, from say 10 up to 10,000. One can even imagine real-time computations in some cases [19].

The 'ultimate' case is that of space-time-parameter problems involving spatial fields, timeline and design variables. This is of course of industrial importance, but still an issue and a current active field of research (see [20] for example). For such problems, the data are generally organized in data cubes (Figure 1). In this paper, a data-driven reducedorder modeling approach is proposed for space-time-parameter mechanical problems involving an equation of kinematics and a quasi-static law of equilibrium. As particular application, the physical problem that is addressed is the dynamics of dilute suspensions of micrometric capsules in microfluidic channels. Microcapsules can be used in Healthcare as innovative drug transportation vehicles into blood vessels and are expected to deliver drugs at identified targets [21,22]. They are composed of an elastic membrane protecting a liquid inner core and are used in suspension in another liquid. Testing them in microfluidic environments offers great potential to determine the capsule behavior and characterize the mechanical properties of the membrane [23–29], but also for sorting or enrichment of capsule suspensions [30–34]. We presently focus on the flow of a dilute suspension of initially spherical micrometric capsules in a microfluidic channel, which is a complex three-dimensional inertialess fluid-structure interaction problem that interestingly depends on only two independent design variables: the capillary number of the capsule *Ca*, which is a non-dimensional number that estimates the order of magnitude of the viscous forces acting on the capsule with respect to the elastic forces that build up in the membrane, and the confinement ratio *a*/ that provides a comparison between the initial capsule diameter and the channel width.

**Figure 1.** Space-time-parameter data cube.

Proper Orthogonal Decomposition has been shown to be particularly suitable to build reduced order models (ROM) of microcapsules [35], but so far no model capable of predicting capsule dynamics currently exists. The originality of the paper is to propose a ROM of the capsule-fluids interactions which provides the time-evolution of the capsule shape for any parameter values. From the capsule shape, it is indeed possible to deduce all the quantities of interest (viscous load, internal tensions within the membrane, membrane energy, etc.) in post-treatment. The ROM is inspired from the physical problem, in which the boundary condition stipulates that the fluid velocity equals the capsule membrane velocity. The correction of the capsule node position field can thus be obtained by integrating the velocity field over time. The challenge remains to correlate the position and velocity fields, which we propose to do with diffuse approximation and manifold learning [18,36,37] using the principal modes of both fields obtained by POD decomposition.

Numerical experiments will demonstrate the accuracy and efficiency of the approach, comparing reduced-order solutions to the full-order ones. We believe that the methodology proposed in this paper can be applied to a broad range of multiphysics problems such as fluid-structure interactions, structural dynamics using quasi-static structural mechanics models and related problems.

This paper is organized as follows. In Section 2, we describe the physical problem and its full order model solution. In Section 3.2, we construct parametric and spatial reduced-order modes using a set of pre-computed simulations. This allows to introduce a reduced-order model that expresses the displacement and velocity of a capsule at a selection of snapshots. In Section 3.4, we build a reduced model that corresponds to any parameter vector of *Ca* and *a*/ values by estimating its corresponding principal components. Then, with the use of a Diffuse Approximation (DA) method, we adopt a data-driven manifold learning to predict the deformation of the capsule in the flow for a chosen time discretization. Finally, in Section 4, we will validate the whole computational ROM approach with a comparison to the FOM solutions.

#### **2. Material and Methods**

#### *2.1. Problem Statement*

Let us consider a spherical micrometric capsule of radius *a* freely placed in a threedimensional microfluidic channel with square cross-section of length 2-(see Figure 2). The

capsule and the channel are filled with an incompressible Newtonian fluid of the same constant density *ρ* and dynamic viscosity *μ*. The capsule is enclosed by a thin hyperelastic isotropic membrane (surface shear modulus *Gs*, area expansion modulus *Ks* = 3*Gs*). It is subjected to a Poiseuille flow of mean velocity *V*.

**Figure 2.** Initial configuration considered in the FOM model: an initially spherical capsule is placed at the center of a square-section channel. The time-evolution of its dynamics is computed using a reference frame centred onto the capsule centre of mass.

The problem is governed by two dimensionless numbers:


The Reynolds number of the external flow is assumed to be very small (typically of order 10−<sup>2</sup> or less), inertia being negligible because of the spatial scales involved in the problem. As far as the internal flow is concerned, its velocity is induced by the motion of the capsule membrane, which is itself entrained by the external flow: it is thus of smaller amplitude than that of the external flow. Hence, the flow in the internal (*β* = *in*) and external (*β* = *ex*) fluids are described by the Stokes equations:

$$
\nabla \cdot \mathbf{v}^{\beta} = 0,\\
\nabla \cdot \sigma^{\beta} = 0,\\
\oint \mathbf{j} = \text{in}, \text{ex.} \tag{1}
$$

where *σ<sup>β</sup>* is stress tensor in the fluids. The Stokes equations are defined in the domains bounded by the capsule membrane for *β* = *in* and by the capsule membrane and the channel wall for *β* = *ex*. The inlet Γ*in* and outlet Γ*out* cross-sections of the channel are assumed to be far from the capsule (10 in the FOM model). The reference frame (*O*, *x*, *y*, *z*) is fixed on the capsule center of mass *O* at each time step. For the velocity vector field *v<sup>β</sup>* and the pressure field *pβ*, we consider the following boundary conditions (note that the boundary conditions include wall confinement effects, see [38] for more details):

• The flow perturbation induced by the capsule vanishes at Γ*in* and Γ*out*:

$$
\sigma^{ex}(\mathbf{x}, t) \to \mathfrak{v}^{\infty}(\mathbf{x}), \text{ when } \mathbf{x} \in \Gamma\_{in} \cup \Gamma\_{out} \tag{2}
$$

where *v*<sup>∞</sup> is the Poiseuille flow velocity of the suspending fluid in the absence of capsule. For a square channel we have the expression in expansion form

$$v^{\infty}(x,y) = \frac{\sum\_{n=1,3,\dots}^{\infty} \frac{\pi V}{n^3} \left[1 - \frac{\cosh(n\pi x/\ell)}{\cosh(n\pi/2)}\right] \sin(n\pi(y/\ell + 1/2))}{2\left[\frac{\pi^4}{96} - \sum\_{n=1,3,\dots}^{\infty} \frac{\tanh(n\pi\pi/2)}{n^5\pi/2}\right]}.\tag{3}$$

• Uniform pressure at Γ*in* and Γ*out*:

$$p^{\varepsilon x}(\mathbf{x}, t) = 0, \ \mathbf{x} \in \Gamma\_{\mathrm{in}\nu} \tag{4}$$

$$p^{\varepsilon x}(\mathbf{x}, t) = \Delta p(t) + \Delta p^{\infty}, \; \mathbf{x} \in \Gamma\_{\text{out}\_{\prime}} \tag{5}$$

where Δ*p*<sup>∞</sup> is the undisturbed suspending pressure drop in the absence of capsule and Δ*p* is the additional pressure drop due to the capsule.

• No slip boundary conditions on the channel wall *W*:

$$\mathbf{w}^{\mathrm{ex}}(\mathbf{x}, t) = 0, \text{ for } \mathbf{x} \in \mathcal{W} \tag{6}$$

• No slip boundary conditions on the capsule membrane *C*:

$$\boldsymbol{\sigma}^{in}(\mathbf{x},t) = \boldsymbol{\sigma}^{ex}(\mathbf{x},t) = \frac{\partial}{\partial t} \mathbf{x}(\mathbf{X},t), \text{ for } \mathbf{x} \in \mathbb{C} \tag{7}$$

where *<sup>∂</sup> <sup>∂</sup><sup>t</sup> <sup>x</sup>*(*X*, *<sup>t</sup>*) is the membrane velocity at position *<sup>x</sup>* at time *<sup>t</sup>*, and *<sup>X</sup>* is the reference position vector of the capsule membrane.

• The normal loading continuity indicates that the load *q* on the membrane is due to the viscous traction jump

$$\left[\sigma^{ex}(\mathbf{x}) - \sigma^{in}(\mathbf{x})\right] \cdot \mathfrak{n} = \mathfrak{q}\_{\prime} \text{ for } \mathbf{x} \in \mathbb{C} \tag{8}$$

where *n* is the outward unit normal vector.

As the membrane thickness is negligibly small compared to the capsule dimensions, the membrane can be considered as a hyperelastic surface devoid of bending stiffness. The in-plane deformation is then measured by the principal extension ratios *λ*<sup>1</sup> and *λ*2, that measure the in-plane deformation. Owing to the combined effects of hydrodynamic forces, boundary confinement, and membrane deformability, the capsule can be highly deformed. Consequently, the choice of membrane constitutive law is important. We consider the Neo-Hookean (NH) constitutive law that models the membrane as an infinitely thin sheet of a three-dimensional isotropic and incompressible material. It was indeed shown to adequately model microcapsules with a cross-linked proteic membrane [23,24,39]. The principal Cauchy in-plane tensions *τ<sup>i</sup>* (*i* = 1, 2) (forces per unit arc length of deformed surface curves) can be expressed as a function of the principal extension ratios:

$$\pi\_1 = \frac{G\_s}{\lambda\_1 \lambda\_2} \left[ \lambda\_1^2 - \frac{1}{\left(\lambda\_1 \lambda\_2\right)^2} \right], \quad \pi\_2 = \frac{G\_s}{\lambda\_1 \lambda\_2} \left[ \lambda\_2^2 - \frac{1}{\left(\lambda\_1 \lambda\_2\right)^2} \right]. \tag{9}$$

#### *2.2. Discrete Full Order Model (FOM)*

The Fluid-Structure Interaction (FSI) problem is numerically modeled by coupling the Boundary Integral Method (BIM) that solves the fluid Equations (2)–(8) with the Finite Element Method (FEM) that solves the membrane mechanical problem [38,40] using the Caps3D in-house code. The unknowns are the discrete displacement field {*u*} and the discrete velocity field {*v*} at the nodes of the membrane mesh. The equation of kinematics states that *<sup>d</sup> dt* {*u*} <sup>=</sup> {*v*}. The forces exerted onto the membrane are computed by the FEM. The deformation of the membrane is computed from the velocity vector field obtained at the membrane nodes by solving the Stokes equations with the BIM, leading to a nonlinear relation written in abstract form {*v*} <sup>=</sup> {N }({*u*}). For a given parameter vector *θ* = (*θ*1, *θ*2) *<sup>T</sup>*, where *θ*<sup>1</sup> = *Ca* and *θ*<sup>2</sup> = *a*/-, the time-continuous semi-discrete FSI scheme reads in abstract form

$$\begin{aligned} \frac{d}{dt} \{\mathfrak{u}\}(t) &= \{\mathfrak{v}\}(t), \\ \{\mathfrak{v}\}(t) &= \{\mathcal{N}\}(\{\mathfrak{u}\}(t), \mathfrak{f}), \quad t \in (0, T\_f], \\ \{\mathfrak{u}\}(0) &= \{\mathbf{0}\}, \; \{\mathfrak{v}\}(t) = \{\mathcal{N}\}(\{\mathbf{0}\}, \mathfrak{f}) \end{aligned}$$

where {*u*}(*t*) and {*v*}(*t*) represent the discrete FE displacement and velocity fields at continuous time *t*, and *Tf* is the final time. For time discretization, either a forward Euler scheme or a second order Runge-Kutta scheme is used with a suitable constant time step *δt* > 0. The Euler scheme reads

$$\begin{aligned} \{\mathfrak{u}^{i+1}\} &= \{\mathfrak{u}^{i}\} + \delta t \{\mathfrak{v}^{i}\}, \\ \{\mathfrak{v}^{i+1}\} &= \{\mathcal{N}\} (\{\mathfrak{u}^{i+1}\}, \mathfrak{f}), \\ \{\mathfrak{u}^{0}\} &= \{\mathbf{0}\}, \{\mathfrak{v}^{0}\} = \{\mathcal{N}\} (\{\mathbf{0}\}, \mathfrak{f}). \end{aligned}$$

where {*u<sup>i</sup>* } and {*v<sup>i</sup>* } represent the discrete FE displacement and velocity fields at discrete time *t <sup>i</sup>*,*<sup>δ</sup>* <sup>=</sup> *<sup>i</sup> <sup>δ</sup><sup>t</sup>* <sup>≤</sup> *Tf* . For second-order accuracy in time, a Runge-Kutta Ralston scheme is used:

$$\begin{aligned} \{\hat{\mathfrak{u}}^{i+2/3}\} &= \{\mathfrak{u}^{i}\} + \frac{2}{3}\delta t \{\mathfrak{v}^{i}\}, \\ \{\hat{\mathfrak{v}}^{i+2/3}\} &= \{\mathcal{N}\} (\{\hat{\mathfrak{u}}^{i+2/3}\}, \mathfrak{g}), \\ \{\mathfrak{u}^{i+1}\} &= \{\mathfrak{u}^{i}\} + \frac{\delta t}{4} \left(\{\mathfrak{v}^{i}\} + 3\{\hat{\mathfrak{v}}^{i+2/3}\}\right), \\ \{\mathfrak{v}^{i+1}\} &= \{\mathcal{N}\} (\{\mathfrak{u}^{i+1}\}, \mathfrak{g}), \\ \{\mathfrak{u}^{0}\} &= \{\mathbf{0}\}, \{\mathfrak{v}^{0}\} = \{\mathcal{N}\} (\{\mathbf{0}\}, \mathfrak{g}). \end{aligned}$$

Because of the explicit nature of the numerical schemes for the equation of kinematics, the time step is subject to a Courant-Friedrichs-Lewy (CFL)-like stability condition

$$
\dot{\gamma}\,\delta t < \mathbb{C}\,\frac{\Delta h\_{\mathbb{C}}}{\ell}\,\mathbb{C}a,\tag{10}
$$

,

where *γ*˙ = *V*/-, *C* > 0 is a constant and Δ*hC* is the typical mesh size (see [40]). In practice, we first use small time steps, and tune them to be big enough but not too close to the stability boundary. This process does not take too much time.

#### *2.3. Design of Experiment, Database of FOM Results*

Simulations of the FOM problem have been run varying the governing parameters in the range [0; 0.2] for the capillary number *Ca* and [0.75; 1.2] for the confinement ratio *a*/-. For *a*/- ≥ 0.95, the capsule is initially pre-deformed into an ellipsoid of semi-minor axis equal to 0.9. This pre-deformation does not have any impact on the steady-state capsule shape and is enough to avoid contacts between the capsule membrane and the channel wall. The resulting numerical database, composed of *Nc* = 118 (*Ca*, *a*/*l*) samples (Figure 3), contains the time-evolution of the three-dimensional position (or displacements) vectors of the capsule membrane nodes. Only the configurations for which a steady-state shape has been reached are retained. No steady state is found above the dotted red line of Figure 3, the microcapsules exhibiting continuous elongation owing to the strain-softening behavior of the membrane law [41].

**Figure 3.** (**a**) Values of *Ca* and *a*/*l* included in the FOM database in the case of an initially spherical capsule with a Neo-Hookean membrane flowing in a square-section microfluidic channel. The parameter domain where a steady capsule deformation can be reached is delimited by the red dotted line. (**b**) Time evolution of capsule deformation along the microfluidic channel shown as illustration for the 6 cases indicated with numbers in figure (**a**). The capsule is pre-deformed into an ellipsoid when *a*/-≥ 0.95.

In the ROM model, we consider the capsule positions in the laboratory reference frame (and not the reference frame centred on the capsule centre of mass as in the FOM model). The capsule thus moves along the microchannel. For data generation, we pick up time snapshot solutions at coarser discrete times *t <sup>i</sup>* = *i* Δ*t*, where Δ*t* = *m δt* for some integer *<sup>m</sup>* <sup>≥</sup> 1. The total number of coarse discrete times is denoted *Nt*. Let {*X*}[*n*] and {*x*}[*n*](*<sup>t</sup> i* ) <sup>∈</sup> <sup>R</sup>3, *<sup>n</sup>* <sup>=</sup> 1, ... , *Nx*, be the coordinates of node number *<sup>n</sup>* of the capsule mesh in the reference configuration (at time *t* = 0) and at discrete time *t i* , *i* = 1, ... , *Nt*. The coordinates in the current configuration {*x*}(*<sup>t</sup> i* ) are function of the (*Ca*, *a*/-) parameter value denoted *θ <sup>j</sup>*, *j* = 1, ... , *Nc*. The database is thus stored as a datacube of 3D-space, 1D-time and 2D-parameter data. The displacement vector is then {*u*} {*X*}[*n*], *<sup>t</sup> i* , *θ j* = {*x*}[*n*](*<sup>t</sup> i* ) − {*X*}[*n*]. The velocity vector {*v*} is calculated by finite differences from the position vector. Typically, for a standard capsule FOM simulation, *Nx* is of order 103 and *Nc* of order 102. The time step Δ*t* is chosen such that *Nt* is of order 102.

#### **3. Reduced Order Model (ROM)**

Reduced order modeling aims at deriving a lightweight model of low-order dimension from solutions obtained by the FOM, while trying to keep the same order of accuracy. There are many reasons for doing that. In particular, parameter exploration and sensitivity analysis are made easier because of large speedups using the ROM compared to the prohibitive FOM computational time. One can also imagine real-time parameter exploration and visualization of capsule evolution.

#### *3.1. Overview*

We first give a general overview of the proposed data-driven model order reduction methodology. The approach is classically made of an offline stage for the search of the principal components and POD coefficient matrices of the FOM solutions, followed by an online stage where a parameter is chosen and a low-order dynamical system is run to get the solutions.

1. **Offline stage**. We reduce the data dimensionality by means of a double POD basis for space and parameter variables. The displacement field is represented as

$$\{\mathfrak{u}\}(\{X\},t,\mathfrak{G}) = \sum\_{k=1}^{K\_u^c} \sum\_{\ell=1}^{K\_u^c} A\_{k\ell}(t) \left\{ \Phi\_u^{\ell} \right\}\_k (\psi\_u(\mathfrak{G}))\_{\ell\prime} \tag{11}$$

where {**Φ***<sup>r</sup> <sup>u</sup>*}*<sup>k</sup>* <sup>∈</sup> <sup>R</sup>3*Nx* are the spatial POD modes, *<sup>ψ</sup>u*(*θ*) <sup>∈</sup> <sup>R</sup>*K<sup>c</sup> <sup>u</sup>* the parameter modes and *Ak*-(*t*) scalar coefficients depending on time *t*. The truncation ranks are *K<sup>x</sup> <sup>u</sup>* and *Kc <sup>u</sup>*, respectively (the '*x*' superscript stands for 'space' and the '*c*' superscript for 'configuration'). We use a similar representation for the velocity field:

$$\{\sigma\}(\{X\}, t, \boldsymbol{\mathfrak{G}}) = \sum\_{k=1}^{k\_v^x} \sum\_{\ell=1}^{k\_v^x} B\_{k\ell}(t) \left\{ \boldsymbol{\Phi}\_v^r \right\}\_k (\boldsymbol{\upmu}\_v(\boldsymbol{\upTheta}))\_\ell. \tag{12}$$

The determination of the double POD basis is achieved by singular value decomposition (SVD) from the datacube with different rearrangements of the data in stacked matrix form. The truncation ranks *K<sup>x</sup> <sup>u</sup>*, *K<sup>c</sup> <sup>u</sup>*, *K<sup>x</sup> <sup>v</sup>* , *K<sup>c</sup> <sup>v</sup>* are expected to be rather small while ensuring accuracy of the representations.

	- (a) Estimate the displacement field {*u*}({*X*}, *<sup>t</sup>*, *<sup>θ</sup>q*) from expression (11). This requires an interpolation process at *θ* = *θ <sup>q</sup>*. For that, we decide to use a diffuse approximation technique [18] that can be used for any parameter space dimension;
	- (b) From the estimated displacement field {*u*}({*X*}, *<sup>t</sup> i* , *θq*) computed at different instants *t <sup>i</sup>* <sup>∈</sup> [0, *Tf* ], compute a low-order reduced basis {*ϕk*}(*<sup>θ</sup> <sup>q</sup>*), *<sup>k</sup>* <sup>=</sup> 1, ... , *mu* by singular value decomposition. We then get the low-order representations of both displacements and velocities:

$$\left(\{\mu\}\left(\{X\},t,\theta\_q\right)\right) = \sum\_{k=1}^{m\_u} a\_k(t) \left\{\mathfrak{q}^k\right\}(\theta\_q),\tag{13}$$

$$\left(\{\boldsymbol{v}\}\left(\{\boldsymbol{X}\},\boldsymbol{t},\boldsymbol{\theta}\_{q}\right)\right) = \sum\_{k=1}^{m\_{v}} \zeta\_{k}(\boldsymbol{t}) \left\{\boldsymbol{\gamma}^{k}\right\}(\boldsymbol{\theta}\_{q}) ,\tag{14}$$

(c) Manifold learning online stage: using diffuse approximation, we determine the low-order manifold M that links displacements and velocities in the (reduced-order) state space:

$$\xi = \mathcal{M}(\mathfrak{a}, \mathfrak{b}\_{\mathbb{Q}});$$

(d) Derivation of a low-order dynamical system: we then derive a lightweight differential-algebraic dynamical system, easy to solve numerically: for *θ* = *θ <sup>q</sup>*, solve

$$\begin{aligned} \frac{d\mathfrak{a}}{dt} &= Q\,\mathfrak{F}(t),\\ \mathfrak{F}(t) &= \mathcal{M}(\mathfrak{a}(t), \mathfrak{G}\_{\mathfrak{q}}). \end{aligned}$$

The high-dimensional displacement and velocity fields can then be reconstructed according to (13) and (14).

In the next section, we give all the details of the ROM methodology.

#### *3.2. Offline Stage*

#### 3.2.1. Global Parametric Reduced Basis (GPRB)

This first step consists in computing a parametric reduced basis in the whole parameter domain from the database of FOM results (see Section 2.3). For simplification reasons, we use the subscript that can be either *u* or *v* to express displacements and velocity respectively in the formulas.

Let *S<sup>i</sup> <sup>u</sup>* ∈ M3*Nx*,*Nc* (R) be the matrix of capsule displacement fields {*u*} and *Si <sup>v</sup>* ∈ M3*Nx*,*Nc* (R) the matrix of the velocity fields {*v*} at time *<sup>t</sup> i* , *i* = 1, ... , *Nt* (Figure 4a), considering all the configurations *θ <sup>j</sup>* for *j* = 1, . . . , *Nc* of the database, i.e.,

$$S\_{\mathfrak{u}}^{\bar{i}} = \left[ \{ \mathfrak{u} \} \left( \{ X \} , t^{\bar{i}} , \mathfrak{G}\_1 \right) , \dots , \; \{ \mathfrak{u} \} \left( \{ X \} , t^{\bar{i}} , \mathfrak{G}\_{\mathbb{N}\_c} \right) \right] \_{\prime \prime}$$

and

$$\mathcal{S}\_v^i = \left[ \{ v \} \left( \{ X \} , t^i , \theta\_1 \right) , \dots , \; \{ v \} \left( \{ X \} , t^i , \theta\_{N\_c} \right) \right] .$$

Then we stack all the matrices *S<sup>i</sup>* for *<sup>i</sup>* <sup>=</sup> 1, . . . , *Nt* into a big matrix <sup>S</sup> ∈ M3*Nx*×*Nt*,*Nc* (R):

$$\mathcal{S}\_{\boldsymbol{\varrho}} = \begin{bmatrix} \mathbf{S}\_{\boldsymbol{\varrho}}^1 \\ \mathbf{S}\_{\boldsymbol{\varrho}}^2 \\ \vdots \\ \mathbf{S}\_{\boldsymbol{\varrho}}^N \end{bmatrix} \\ \text{for } \boldsymbol{\varrho} = \boldsymbol{\mu}, \boldsymbol{v}.$$

We then apply SVD [42] and get:

$$\mathcal{S}\_{\emptyset} = \mathcal{U}\_{\emptyset} \Sigma\_{\mathbb{S}\_{\emptyset}} \Psi\_{\emptyset}^{T}, \text{ for } \emptyset = u, v,\tag{15}$$

where *<sup>U</sup>* ∈ M3*NxNt*,*Nc* (R), <sup>Ψ</sup> ∈ M*Nc* (R) are semi-orthogonal and orthogonal matrices, respectively, and <sup>Σ</sup>*S* ∈ M*Nc* (R) is the diagonal singular value matrix. The matrix <sup>Ψ</sup> of discrete parameter modes can be truncated according to *K<sup>c</sup>* parameters, so we note:

$$\begin{cases} \boldsymbol{\Psi}\_{\boldsymbol{\varrho}}^{r} = \left[ (\boldsymbol{\Psi}\_{\boldsymbol{\varrho}})\_{1}, \dots, (\boldsymbol{\Psi}\_{\boldsymbol{\varrho}})\_{K\_{\boldsymbol{\varrho}}^{c}} \right] \in \mathcal{M}\_{N\_{\boldsymbol{\varrho}}, K\_{\boldsymbol{\varrho}}^{c}}(\mathbb{R}),\\ \text{with } (\boldsymbol{\Psi}\_{\boldsymbol{\varrho}})\_{k} \in \mathcal{M}\_{N\_{\boldsymbol{\varrho}}, 1}(\mathbb{R}) \text{ for } k = 1, \dots, K\_{\boldsymbol{\varrho}'}^{c} \text{ and } \boldsymbol{\varrho} = \boldsymbol{u}\_{r} v. \end{cases} \tag{16}$$

The orthogonality property ensures that Ψ*r <sup>T</sup>* Ψ*r* = *IKc* .

**Figure 4.** FOM data rearrangements for (**a**) parametric data set selection and (**b**) spatial data set selection.

3.2.2. Global Spatial Reduced Basis (GSRB)

*Tj <sup>u</sup>* = {*u*} {*X*}, *<sup>t</sup>* 1, *θ <sup>j</sup>* 

Similarly, we build a global spatial reduced basis that captures the spatial data of capsule displacements. Let *T<sup>j</sup> <sup>u</sup>* ∈ M3*Nx*,*Nt*(R) be the displacement matrix and *<sup>T</sup><sup>j</sup> <sup>v</sup>* the velocity matrix for the *j*-th configuration *θ <sup>j</sup>*, for *j* = 1, ... , *Nc* at all time instants *t i* , *i* = 1, . . . , *Nt* (Figure 4b):

> , ..., {*u*}

and

$$T^{j}\_{\mathbb{P}} = \left[ \{ \mathfrak{v} \} \left( \{ \mathbf{X} \} , t^{1} , \mathfrak{\theta}\_{j} \right) , \dots , \; \{ \mathfrak{v} \} \left( \{ \mathbf{X} \} , t^{N\_{l}} , \mathfrak{\theta}\_{j} \right) \right] .$$

Then we define the global matrix <sup>T</sup> ∈ M3*Nx*,*Nt*×*Nc* (R) that horizontally gathers all the matrices *T<sup>j</sup>* for *j* = 1, . . . , *Nc* and = *u*, *v*, respectively:

$$\mathcal{T}\_{\emptyset} = \left[ T\_{\emptyset}^1, T\_{\emptyset}^2, \dots, T\_{\emptyset}^{N\_c} \right], \text{ for } \emptyset = \mathfrak{u}, \mathfrak{v}.$$

The SVD decomposition is applied on T to get

$$\mathcal{T}\_{\emptyset} = \Phi\_{\emptyset} \Sigma\_{T\_{\emptyset}} V\_{\emptyset}^T,\text{ for } \emptyset = u, v,\tag{17}$$

{*X*}, *<sup>t</sup>Nt* , *<sup>θ</sup> <sup>j</sup>*

 ,

where <sup>Φ</sup> ∈ M3*Nx* (R), *<sup>V</sup>* ∈ M*NcNt*,3*Nx* (R) are orthogonal and semi-orthogonal matrices, respectively, and <sup>Σ</sup>*S* ∈ M3*Nx* (R) is the diagonal singular value matrix with singular values organized in decreasing order. We can also apply a spatial basis truncation at a range of *K<sup>x</sup>* for a specified accuracy threshold. The reduced spatial POD basis is stored in the matrix:

$$\Phi\_{\boldsymbol{\upphi}}^{\boldsymbol{r}} = \left[ \{ \boldsymbol{\upPhi}\_{\boldsymbol{\upphi}} \}\_{1 \times} , \dots , \{ \boldsymbol{\upPhi}\_{\boldsymbol{\upphi}} \}\_{K\_{\boldsymbol{\upphi}}^{\boldsymbol{x}}} \right] \in \mathcal{M}\_{\mathbf{3N}\_{\boldsymbol{x}}, \mathbb{K}\_{\boldsymbol{\upphi}}^{\boldsymbol{x}}}(\mathbb{R}) \tag{18}$$

with the orthogonality property Φ*r <sup>T</sup>* Φ*r* = *IKx* , = *u*, *v*.

#### *3.3. Data Dimensionality Reduction*

Once the POD modes of S and T for the displacement fields ( = *u*) and the velocity fields ( = *v*) are computed, one can summarize (approximate) capsule displacement and velocity fields of the database at any discrete time *t <sup>i</sup>* (*i* = 1..., *Nt*) as

$$\left(\{\boldsymbol{u}\}\left(\{\boldsymbol{X}\},\,\boldsymbol{t}^{i},\,[\boldsymbol{\theta}\_{1},\ldots,\boldsymbol{\theta}\_{N\_{\mathrm{c}}}]\right)\approx\boldsymbol{\Phi}\_{\boldsymbol{u}}^{r}\boldsymbol{A}\left(\boldsymbol{t}^{i}\right)\left(\boldsymbol{\Psi}\_{\boldsymbol{u}}^{r}\right)^{T}\in\mathcal{M}\_{\mathrm{3N}\_{\mathrm{r}},\mathrm{N}\_{\mathrm{c}}}\left(\mathbb{R}\right),\tag{19}$$

$$\left(\{\boldsymbol{v}\}\left(\{\boldsymbol{X}\},\boldsymbol{t}^{i},\left[\boldsymbol{\theta}\_{1},\ldots,\boldsymbol{\theta}\_{N\_{\mathcal{E}}}\right]\right)\approx\boldsymbol{\Phi}\_{\boldsymbol{v}}^{r}\boldsymbol{B}\left(\boldsymbol{t}^{i}\right)\left(\boldsymbol{\Psi}\_{\boldsymbol{v}}^{r}\right)^{T}\in\mathcal{M}\_{\operatorname{3N}\_{\mathcal{E}}N\_{\mathcal{E}}}\left(\mathbb{R}\right),\tag{20}$$

where *A*(*t i* ) ∈ M*K<sup>x</sup> <sup>u</sup>*,*K<sup>c</sup> <sup>u</sup>* (R) and *<sup>B</sup> t i* ∈ M*K<sup>x</sup> <sup>v</sup>* ,*K<sup>c</sup> <sup>v</sup>* (R) are some coefficient matrices depending on time *t i* . If the approximation is chosen as the orthogonal projection over the vector spaces spanned by the POD modes, the coefficient matrices are computed as follows for *i* = 1..., *Nt*:

$$A(t^{i}) = \underbrace{\left(\Phi\_{\boldsymbol{u}}^{r}\right)^{T}}\_{K\_{\boldsymbol{u}}^{x}\times\left(3N\_{\boldsymbol{x}}\right)}\underbrace{\{\boldsymbol{u}\}\left(\{\mathbf{X}\},t^{i},[\boldsymbol{\Phi}\_{1},\dots,\boldsymbol{\Phi}\_{N\_{\boldsymbol{c}}}]\right)}\_{\left(3N\_{\boldsymbol{x}}\right)\times N\_{\boldsymbol{c}}}\underbrace{\Psi\_{\boldsymbol{u}}^{r}}\_{N\_{\boldsymbol{c}}\times K\_{\boldsymbol{u}}^{c}}\tag{21}$$

$$B(t^{i}) = \underbrace{\left(\Phi\_{\boldsymbol{v}}^{r}\right)^{T}}\_{\mathbf{K}\_{\boldsymbol{v}}^{x} \times \left(3\mathbb{N}\_{\mathrm{x}}\right)} \underbrace{\left\{\boldsymbol{v}\right\} \left(\{\boldsymbol{X}\}, t^{i}, [\boldsymbol{\Phi}\_{\boldsymbol{1}}, \dots, \boldsymbol{\Phi}\_{\boldsymbol{N}\_{\boldsymbol{c}}}]\right)}\_{\left(3\mathbb{N}\_{\boldsymbol{x}}\right) \times \mathbb{N}\_{\boldsymbol{c}}} \underbrace{\Psi\_{\boldsymbol{v}}^{r}}\_{\mathbf{N}\_{\boldsymbol{c}} \times \boldsymbol{K}\_{\boldsymbol{v}}^{c}}.\tag{22}$$

The outputs of the offline stage are respectively the POD matrices Φ*<sup>r</sup> <sup>u</sup>*, Φ*<sup>r</sup> <sup>v</sup>*, Ψ*<sup>r</sup> <sup>u</sup>*, Ψ*<sup>r</sup> <sup>v</sup>* and the small matrices *A*(*t i* ), *B*(*t i* ), *i* = 1, ... , *Nt*. The next online stage will operate on the summarized data (19), (20) with coefficients matrices (21), (22). The algorithm of the offline phase is summarized in Algorithm 1.

#### **Algorithm 1** Offline phase

**Require:** database of *θ<sup>k</sup>* for *k* = 1, . . . , *Nc*, truncations *K<sup>c</sup>* , number of snapshots *Nt*. **for** *i* ← 1, . . . , *Nt* **do if** ( = *u*) **then** *Si <sup>u</sup>* <sup>←</sup> {*u*} {*X*}, *<sup>t</sup> i* , *θ*<sup>1</sup> , ..., {*u*} {*X*}, *<sup>t</sup> i* , *θ Nc* ; <sup>S</sup>*<sup>u</sup>* <sup>←</sup> [S*u*; *<sup>S</sup><sup>i</sup> u*]; **else** *Si <sup>v</sup>* <sup>←</sup> {*v*} {*X*}, *<sup>t</sup> i* , *θ*<sup>1</sup> , ..., {*v*} {*X*}, *<sup>t</sup> i* , *θ Nc* ; <sup>S</sup>*<sup>v</sup>* <sup>←</sup> [S*v*, *<sup>S</sup><sup>i</sup> v*]; **end if end for for** *j* ← 1, . . . , *Nc* **do if** ( = *u*) **then** *Tj <sup>u</sup>* <sup>←</sup> {*u*} {*X*}, *<sup>t</sup>* 1, *θ <sup>j</sup>* , ..., {*u*} {*X*}, *<sup>t</sup>Nt* , *<sup>θ</sup> <sup>j</sup>* ; <sup>T</sup>*<sup>u</sup>* <sup>←</sup> [T*u*, *<sup>T</sup><sup>j</sup> u*]; **else** *Tj <sup>v</sup>* <sup>←</sup> {*v*} {*X*}, *<sup>t</sup>* 1, *θ <sup>j</sup>* , ..., {*v*} {*X*}, *<sup>t</sup>Nt* , *<sup>θ</sup> <sup>j</sup>* ; <sup>T</sup>*<sup>v</sup>* <sup>←</sup> [T*v*, *<sup>T</sup><sup>j</sup> v*]; **end if end for** Φ ← SVD(S), Ψ ← SVD(T), for ← *u*, *v*; **for** *i* = 1, . . . , *Nt* **do** *A*(*t i* ) <sup>←</sup> (Φ*<sup>r</sup> u*) *<sup>T</sup>* {*u*} {*X*}, *<sup>t</sup> i* , [*θ*1,..., *θ Nc* ] Ψ*r u*; *B*(*t i* ) <sup>←</sup> (Φ*<sup>r</sup> v*) *<sup>T</sup>* {*v*} {*X*}, *<sup>t</sup> i* , [*θ*1,..., *θ Nc* ] Ψ*r v*; **end for**

*3.4. Online Stage: Search for an Approximate Solution at a Query Configuration θ <sup>q</sup>*

In the online stage, a user will ask for an approximate solution at a new (query) configuration *θ* = *θ <sup>q</sup>* that has not been already computed by the FOM solver or is not stored in the database. Ingredients of the online stage will be: (i) the data summarization of the previous offline stage; (ii) a first estimation of the spatio-temporal solution at *θ* = *θ <sup>q</sup>*; (iii) the computation of a low-dimensional spatial reduced basis suitable for *θ* = *θ <sup>q</sup>*; (iv) the construction of a manifold M that links variables of displacements and velocities in the low-order state space to solve the equation of membrane mechanics; (v) finally, the building of a low-order differential-algebraic (DAE) system of equations that defines the reducedorder model. Substeps (ii) and (iv) will make use of diffuse approximation (DA) as a universal approximator for multivariate functions.

#### 3.4.1. First Estimation of the Solutions at *θ* = *θ <sup>q</sup>*

As an introduction, let us assume that, from the parameter sampling {*θ*1, ... , *θ Nc*}, we consider a polynomial Lagrange interpolation with Lagrange polynomials denoted by *Li*(*θ*) such that the Lagrange property

$$L\_i(\mathfrak{G}\_j) = \delta\_{i\bar{j}\prime} \quad 1 \le i, j \le N\_{\mathfrak{C}\_j}$$

is fulfilled (*δij* is the standard Kronecker symbol). Let us denote by *<sup>L</sup>*(*θ*)=(*Lj*(*θ*))*j*=1,...,*Nc* <sup>∈</sup> <sup>R</sup>*Nc* the vector that stores the Lagrange polynomials. Then

$$\mathcal{X}\{\boldsymbol{u}\}\Big(\{\mathcal{X}\},\,\boldsymbol{t}^{i},\boldsymbol{\theta}\_{\boldsymbol{q}}\Big):=\{\boldsymbol{u}\}\Big(\{\mathcal{X}\},\,\boldsymbol{t}^{i},\,[\boldsymbol{\theta}\_{1},\ldots,\boldsymbol{\theta}\_{N\_{c}}]\Big)L(\boldsymbol{\theta}\_{\boldsymbol{q}})\in\mathbb{R}^{3N\_{x}}$$

is an interpolated displacement field at parameter *θ* = *θ <sup>q</sup>* and discrete time *t* = *t i* . One can of course do the same for the velocity field.

Unfortunately, Lagrange polynomial interpolation is not suitable for parameter spaces of arbitrary dimension because of the curse of dimensionality and because it may suffer from instability issues (Runge phenomenon). Rather than using polynomial interpolation, we propose to use a Diffuse Approximation (DA) technique [18,29] which is an approximation method based on local low-order polynomial reconstruction (of order one or two) using a compactly-supported kernel function and weighted least squares. The DA method

is known to be a robust and reliable approach which is less sensitive to the location of the sampling points. Moreover, it can be applied to multivariate functions of arbitrary dimensions, which is interesting for larger or more general parameter spaces. It is particularly suited for the current problem, for which the sampling is performed on a Cartesian grid. It may fail in the occurrence of local point alignment within the cloud points, which does not occur in the present study. The accuracy of the DA method may slightly decrease close to the boundary of the domain, the number of neighboring points being reduced.

To estimate the displacement field for *<sup>θ</sup>* <sup>=</sup> *<sup>θ</sup> <sup>q</sup>*, we look for a vector *<sup>ψ</sup>u*(*<sup>θ</sup> <sup>q</sup>*) <sup>∈</sup> <sup>R</sup>*K<sup>c</sup> u* such that

$$\left(\{\boldsymbol{u}\}\left(\{\mathbf{x}\},\boldsymbol{t}^{\bar{i}},\boldsymbol{\theta}\_{q}\right) = \boldsymbol{\Phi}\_{\boldsymbol{u}}^{\tau}\boldsymbol{A}\left(\boldsymbol{t}^{\bar{i}}\right)\boldsymbol{\Psi}\_{\boldsymbol{u}}(\boldsymbol{\theta}\_{q})\tag{23}$$

returns an approximation of the displacement field at *θ* = *θ <sup>q</sup>*. Similarly for the velocity field, we search for a vector *<sup>ψ</sup>v*(*<sup>θ</sup> <sup>q</sup>*) <sup>∈</sup> <sup>R</sup>*K<sup>c</sup> <sup>v</sup>* that gives

$$\{\mathfrak{v}\}(\{\mathbf{x}\},t^{\dot{i}},\mathfrak{\theta}\_{\dot{q}}) = \mathfrak{\Phi}\_{\upsilon}^{r}B(t^{\dot{i}})\,\mathfrak{\Psi}\_{\upsilon}(\mathfrak{\theta}\_{\dot{q}}).\tag{24}$$

Each vector *<sup>ψ</sup>*(*<sup>θ</sup> <sup>q</sup>*) <sup>∈</sup> <sup>R</sup>*K<sup>c</sup>* can be locally approximated by

$$\Psi\_{\emptyset}(\theta\_{\emptyset}) = \mathcal{A}\_{\emptyset} \ p(\theta\_{\emptyset}), \text{ for } \emptyset = \boldsymbol{\mu}, \boldsymbol{v}, \tag{25}$$

where the matrix A ∈ M*K<sup>c</sup>* ,*m*(R) (to be determined) is the approximation coefficient matrix and *p θ q* <sup>∈</sup> <sup>R</sup>*<sup>m</sup>* is a vector of independent polynomial functions, where

$$\begin{cases} p(\boldsymbol{\theta}) = \begin{bmatrix} 1 & \mathbb{C}a & a/\ell \end{bmatrix}^{\top}, m = 3 & \text{for first order DA}, \\ p(\boldsymbol{\theta}) = \begin{bmatrix} 1 & \mathbb{C}a & a/\ell & \mathbb{C}a \end{bmatrix} \begin{bmatrix} \mathbb{C}a \end{bmatrix}^{2} \quad \begin{pmatrix} a/\ell \end{pmatrix}^{2}, m = 6 & \text{for second order DA}. \end{cases} \tag{26}$$

To approximate *ψ*(*θ <sup>q</sup>*), let us consider a neighborhood S (*θ <sup>q</sup>*) centered on *θ <sup>q</sup>* containing *M* neighboring points (Figure 5a). It is an ellipse of equation

$$\left(\theta\_1 - (\theta\_q)\_1\right)^2 + \vec{r}^2 \left(\theta\_2 - (\theta\_q)\_2\right)^2 = R^2$$

where *r*˜ is fixed (equal to 1.9 in Figure 5a) and *R* is chosen such that the ellipse contains *M* points (*M* being chosen by the operator). In other words, the distance between *θ* = (*θ*1, *θ*2)*<sup>T</sup>* and *θ <sup>q</sup>* is

$$d = \left( \left( \theta\_1 - (\theta\_q)\_1 \right)^2 + \vec{r}^2 \left( \theta\_2 - (\theta\_q)\_2 \right)^2 \right)^{\frac{1}{2}} / R. \tag{27}$$

The compactly supported Wendland weight function shown in Figure 5b is classically used. It has appropriate high-order approximation properties ([43]):

$$\begin{cases} w(d) = 2 \, d^3 - 3 \, d^2 + 1, \quad d \le 1, \\ 0, & \text{otherwise.} \end{cases} \tag{28}$$

Diffuse approximation consists in minimizing the weighted least square problem

$$\min\_{\mathcal{A}\_{\boldsymbol{\theta}} \in \mathcal{M}\_{\mathcal{K}\_{\boldsymbol{\theta}^{\rm c}}}(\mathbb{R})} \quad \operatorname{I}\_{\boldsymbol{\theta}\_{\boldsymbol{\theta}}}(\mathcal{A}\_{\boldsymbol{\theta}}) := \frac{1}{2} \sum\_{\boldsymbol{\theta} \in \mathcal{J}'(\mathfrak{B}\_{\boldsymbol{\theta}})} w(d(\boldsymbol{\theta})) \left\| \mathcal{A}\_{\boldsymbol{\theta}} \, p(\boldsymbol{\theta}) - [\Psi\_{\boldsymbol{\theta}}^{\boldsymbol{\tau}}(\boldsymbol{\theta})]^{T} \right\|\_{\mathbb{R}^{K\_{\boldsymbol{\theta}}^{\rm c}}}^{2} \tag{29}$$

where [Ψ*<sup>r</sup>* (*θ*)] *<sup>T</sup>* is the truncated matrix of modes that correspond to couples *θk*, *k* = 1, . . . , *Nc*. The solution A ( = *u*, *v*) of the weighted least square problem (29) is then

$$\mathcal{A}\_{\emptyset} = (\Psi\_{\emptyset}^{r})^{T} \mathcal{W} \mathcal{P} \left(\mathcal{P}^{T} \mathcal{W} \mathcal{P}\right)^{-1} \in \mathcal{M}\_{\mathbb{K}\_{\emptyset}^{c}, \mathfrak{w}}(\mathbb{R}) \tag{30}$$

where the matrix P∈M*Nc*,*m*(R) and the diagonal matrix of weights W∈M*Nc* (R) are defined as ⎡ ⎤

$$\mathcal{P} = \begin{bmatrix} p(\mathbf{\theta\_1})^T \\ \vdots \\ \vdots \\ p(\mathbf{\theta\_{N\_c}})^T \end{bmatrix} \text{and } \mathcal{W} = \begin{bmatrix} w\_1 & 0 & \cdots & 0 \\ 0 & w\_2 & & \vdots \\ \vdots & & \ddots & & \\ & \vdots & & \ddots & & \\ 0 & \cdots & w\_{N\_c} & & & \\ \vdots & & & \ddots & & \\ 0 & \cdots & & & w\_N \end{bmatrix}. \tag{31}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$
 
$$\text{1.05}$$

**Figure 5.** (**a**) DA elliptical region of interest (dashed line) defined around the point *θ <sup>q</sup>* = (*Ca* = 0.055, *a*/- = 0.95) in the parametric space with *M* = 10 neighbors; (**b**) Weight function *w*(*d*).

3.4.2. Construction of a Low-Order Reduced Basis Suitable for *θ* = *θ <sup>q</sup>*, Data Generation

From (23) and (24), one can easily generate some pseudo-snapshot matrices U(*θ <sup>q</sup>*) and V(*θ <sup>q</sup>*) that gather the estimated displacements and velocities at *Nt* discrete times, respectively:

$$\begin{cases} \mathcal{U}(\boldsymbol{\theta}\_{\boldsymbol{q}}) = \left[ \{ \boldsymbol{u} \} \left( \{ \boldsymbol{X} \} , t^{1} , \boldsymbol{\theta}\_{\boldsymbol{q}} \right), \dots, \{ \boldsymbol{u} \} \left( \{ \boldsymbol{X} \} , t^{N\_{l}} , \boldsymbol{\theta}\_{\boldsymbol{q}} \right) \right], \\ \mathcal{V}(\boldsymbol{\theta}\_{\boldsymbol{q}}) = \left[ \{ \boldsymbol{v} \} \left( \{ \boldsymbol{X} \} , t^{1} , \boldsymbol{\theta}\_{\boldsymbol{q}} \right), \dots, \{ \boldsymbol{v} \} \left( \{ \boldsymbol{X} \} , t^{N\_{l}} , \boldsymbol{\theta}\_{\boldsymbol{q}} \right) \right]. \end{cases} \tag{32}$$

One can then apply a new SVD decomposition of matrices U(*θ <sup>q</sup>*) and V(*θ <sup>q</sup>*) respectively to get spatial POD modes {*ϕk*}(*<sup>θ</sup> <sup>q</sup>*) <sup>∈</sup> <sup>R</sup>3*Nx* , *<sup>k</sup>* <sup>=</sup> 1, ... , *mu* for {*u*} and velocity POD modes {*γk*}(*<sup>θ</sup> <sup>q</sup>*) <sup>∈</sup> <sup>R</sup>3*Nx* , *<sup>k</sup>* <sup>=</sup> 1, . . . , *mv* for {*v*}.

$$\text{POD}(\mathcal{U}(\mathfrak{G}\_{\emptyset})) \quad \rightarrow \quad \{\mathfrak{q}^{1}\}(\mathfrak{G}\_{\emptyset}), \ldots, \{\mathfrak{q}^{m\_{u}}\}(\mathfrak{G}\_{\emptyset}) \tag{33}$$

$$\text{POD}(\mathcal{V}(\mathfrak{g}\_q)) \quad \rightarrow \quad \{\gamma^1\}(\mathfrak{G}\_q), \dots, \{\gamma^{m\_v}\}(\mathfrak{G}\_q) \tag{34}$$

where *mu* and *mv* are the truncation ranks of displacement and velocity modes determined in the next section on numerical experiments. One can then search the displacement and velocity fields at *θ* = *θ <sup>q</sup>* as

$$\{\mu\}\left(\{X\}, t, \mathfrak{G}\_q\right) = \sum\_{k=1}^{m\_u} a\_k(t) \left\{\mathfrak{g}^k\right\} \left(\mathfrak{G}\_q\right),\tag{35}$$

$$\{\upsilon\}\left(\{X\}, t, \theta\_q\right) = \sum\_{k=1}^{m\_v} \zeta\_k(t) \left\{\gamma^k\right\}(\theta\_q). \tag{36}$$

By denoting

$$\Phi(\boldsymbol{\theta}\_{\boldsymbol{\theta}}) = \left[ \{ \boldsymbol{\varrho}^{1} \} \{ \boldsymbol{\theta}\_{\boldsymbol{\theta}} \} , \dots , \{ \boldsymbol{\varrho}^{m\_{u}} \} \{ \boldsymbol{\theta}\_{\boldsymbol{\theta}} \} \right] \in \mathcal{M}\_{3\text{N}\boldsymbol{x},\boldsymbol{m}\_{u}}(\mathbb{R}) , \tag{37}$$

$$\Gamma(\boldsymbol{\theta}\_{q}) = \left[ \{ \gamma^{1} \} (\boldsymbol{\theta}\_{q}), \dots, \{ \gamma^{m\_{v}} \} (\boldsymbol{\theta}\_{q}) \right] \in \mathcal{M}\_{\text{3Nx}, m\_{v}}(\mathbb{R}) \tag{38}$$

and *<sup>α</sup>*(*t*)=[*α*1(*t*), ... , *<sup>α</sup>mu* (*t*)]*<sup>T</sup>* <sup>∈</sup> <sup>R</sup>*mu* , *<sup>ξ</sup>*(*t*)=[*ξ*1(*t*), ... , *<sup>ξ</sup>mv* (*t*)]*<sup>T</sup>* <sup>∈</sup> <sup>R</sup>*mv* , we have the vector formulas

$$\left(\{\mu\}\left(\{X\},t,\theta\_{\emptyset}\right)=\Phi(\theta\_{\emptyset})\,\mathfrak{a}\left(t\right), \quad \left\{\mathfrak{v}\right\}\left(\{X\},t,\theta\_{\emptyset}\right)=\Gamma(\theta\_{\emptyset})\,\mathfrak{F}(t).\tag{39}$$

The mode matrices Φ(*θ <sup>q</sup>*) and Γ(*θ <sup>q</sup>*) are assumed to be orthonormal (w.r.t the natural Euclidean inner product), so we have [Φ(*θ <sup>q</sup>*)]*<sup>T</sup>* Φ(*θ <sup>q</sup>*) = *Imu* and [Γ(*θ <sup>q</sup>*)]*<sup>T</sup>* Γ(*θ <sup>q</sup>*) = *Imv* .

3.4.3. Toward a Physically Consistent Dynamical Reduced-Order Model

Consider now the forward Euler scheme on the FSI system with a ROM time step *δtROM* > 0: at time *t <sup>i</sup>*+1,*ROM* = *t <sup>i</sup>*,*ROM* + *δtROM*, the numerical scheme is

$$\{\mathfrak{u}^{i+1}\} = \{\mathfrak{u}^{i}\} + \delta t^{ROM} \{\mathfrak{v}^{i}\},\tag{40}$$

$$\{\boldsymbol{\sigma}^{i+1}\} = \{\mathcal{N}\} (\{\boldsymbol{u}^{i+1}\}, \boldsymbol{\theta}\_{\emptyset}).\tag{41}$$

Let us emphasize that the equation of local mechanical equilibrium depends on the parameter *θ <sup>q</sup>*. For the reduced-order model, we would like to have a similar algebraic structure to (40), (41) but formulated as a low-dimensional system. If {*u<sup>i</sup>* } and {*v<sup>i</sup>* } are searched in the form {*u<sup>i</sup>* } <sup>=</sup> <sup>Φ</sup>(*<sup>θ</sup> <sup>q</sup>*) *<sup>α</sup><sup>i</sup>* and {*v<sup>i</sup>* } <sup>=</sup> <sup>Γ</sup>(*<sup>θ</sup> <sup>q</sup>*) *<sup>ξ</sup><sup>i</sup>* , respectively, Equation (40) becomes

$$\Phi(\boldsymbol{\theta}\_{\boldsymbol{q}}) \,\boldsymbol{a}^{i+1} = \Phi(\boldsymbol{\theta}\_{\boldsymbol{q}}) \,\boldsymbol{a}^{i} + \delta t^{\boldsymbol{R}\boldsymbol{M}\boldsymbol{M}} \,\Gamma(\boldsymbol{\theta}\_{\boldsymbol{q}}) \,\boldsymbol{\xi}^{i}.$$

By multiplying by [Φ(*θ <sup>q</sup>*)]*<sup>T</sup>* on the left, we get the system of *mu* equations

$$\mathfrak{a}^{i+1} = \mathfrak{a}^i + \delta t^{ROM} \, \mathcal{Q}(\mathfrak{e}\_q) \, \mathfrak{z}^i,\tag{42}$$

where *Q*(*θ <sup>q</sup>*)=[Φ(*θ <sup>q</sup>*)]*T*Γ(*θ <sup>q</sup>*). Equation (41) is replaced by

$$
\Gamma(\mathfrak{G}\_q) \mathfrak{F}^{i+1} = \{ \mathcal{N} \} (\Phi(\mathfrak{G}\_q) \, \mathfrak{a}^{i+1}, \mathfrak{G}\_q) .
$$

By multiplying by [Γ(*θ <sup>q</sup>*)]*<sup>T</sup>* on the left, we get

$$\mathfrak{F}^{i+1} = \mathcal{M}(\mathfrak{a}^{i+1}, \mathfrak{G}\_{\mathfrak{q}})$$

where

$$\mathcal{M}(\mathfrak{a}^{i+1}, \mathfrak{G}\_q) = [\Gamma(\mathfrak{G}\_q)]^T \{\mathcal{N}\} (\Phi(\mathfrak{G}\_q) \,\mathfrak{a}^{i+1}, \mathfrak{G}\_q) \in \mathbb{R}^{m\_{\mathcal{C}}}.\tag{43}$$

#### 3.4.4. Manifold Learning

Because of nonlinear terms, the direct computation of <sup>M</sup>(*αi*<sup>+</sup>1, *<sup>θ</sup> <sup>q</sup>*) in (43) requires high-dimensional computations, which makes the ROM irrelevant from a performance point of view. To "identify" a low-order manifold M, we rather adopt a data-driven approach based once again of diffuse approximation. We link the entry data *α<sup>D</sup> <sup>k</sup>* (*t i* ), *k* = 1, ... , *mu*, *i* = 1, ... , *Nt* to the output data *ξ<sup>D</sup> <sup>k</sup>* (*t i* ), *k* = 1, ... , *mv*, *i* = 1, ... , *Nt* (*'D'* stands for *'data'*). For that, one can compute the orthogonal projections of the pseudosnapshots over the POD bases, leading to the formulas

> *αD <sup>k</sup>* (*t i* ) = {*u*}({*X*}, *<sup>t</sup> i* , *<sup>θ</sup> <sup>q</sup>*), {*ϕk*}(*<sup>θ</sup> <sup>q</sup>*) *ξD*(*t i* ) = {*v*}({*X*}, *<sup>t</sup> i* , *<sup>θ</sup> <sup>q</sup>*), {*γk*}(*<sup>θ</sup> <sup>q</sup>*)

and

at instants *t <sup>i</sup>* = *i*Δ*t*. Manifold learning consists in achieving a (nonlinear) regression method that links entry and output data. We are looking for a manifold representation *ξ* = M(*α*, *θ <sup>q</sup>*) in the form

$$\boldsymbol{\xi}\_{k} = \boldsymbol{p}(\boldsymbol{a})^{T}\boldsymbol{a}^{k}, \quad k = 1, \ldots, m\_{\mathbb{D}} \tag{44}$$

where *<sup>p</sup>*(*α*) is the vector made of monomials in *<sup>α</sup>* of order zero and one, and *<sup>a</sup><sup>k</sup>* <sup>∈</sup> <sup>R</sup>*mu*+<sup>1</sup> is a vector of coefficients to be determined from the data. This corresponds to a local linear embedding process. For each *<sup>k</sup>* <sup>=</sup> 1, ... , *mv*, one looks for a coefficient vector *<sup>a</sup>k*(*t*) <sup>∈</sup> <sup>R</sup>*mu*+<sup>1</sup> solution of the weighted least square problem

$$\mathbf{a}^{k}(t) = \arg\min\_{\mathbf{a}\in\mathbb{R}^{m\_{H}+1}} \quad \frac{1}{2} \sum\_{i=1}^{N\_{l}} w\left(\frac{|t - t^{i}|}{R}\right) \left(p(\mathbf{a}^{D}(t^{i}))^{T}\mathbf{a} - \xi\_{k}^{D}(t^{i})\right)^{2} \tag{45}$$

where *<sup>t</sup>* <sup>∈</sup> [0, *Tf* ], *<sup>w</sup>* <sup>=</sup> *<sup>w</sup>*(*d*) is the weight function defined in Figure 5b and *<sup>d</sup>* <sup>=</sup> <sup>|</sup>*t*−*<sup>t</sup> i* | *<sup>R</sup>* . This returns a regression function

$$\mathfrak{g}\_k = \mathfrak{g}\_k(t, \mathfrak{a}(t)) = \mathfrak{p}(\mathfrak{a}(t))^T \mathfrak{a}^k(t). \tag{46}$$

#### 3.4.5. Low-Order Dynamical Reduced Order Model

The resulting time-discrete reduced-order model is then

$$t^{i+1,ROM} = t^{i,ROM} + \delta t^{ROM},\tag{47}$$

$$\mathfrak{a}^{i+1} = \mathfrak{a}^i + \delta t^{\text{ROM}} \, \mathcal{Q}(\mathfrak{e}\_q) \, \mathfrak{F}',\tag{48}$$

$$\xi\_k^{i+1} = p(\mathfrak{a}^{i+1})^T \mathfrak{a}^k (t^{i+1,ROM}) \quad \forall k \in \{1, \ldots, m\_v\}. \tag{49}$$

High-dimensional displacement and velocity fields can be reconstructed as follows:

$$\sigma\left(\mathfrak{u}\right)\left(\{X\},t^{i+1,\text{ROM}},\mathfrak{\theta}\_q\right) = \Phi(\mathfrak{\theta}\_q)\,\mathfrak{a}^{i+1}, \quad \{\upsilon\}\left(\{X\},t^{i+1,\text{ROM}},\mathfrak{\theta}\_q\right) = \Gamma(\mathfrak{\theta}\_q)\,\mathfrak{z}^{i+1}.$$

The online stage of the reduced-order model is summarized in Algorithm 2.

#### **Algorithm 2** Online phase

**Require:** choose a query parameter *θ <sup>q</sup>*, choose a time step *δtROM* > 0. Initialization: *t* = *t* 0,*ROM* = 0, *α*<sup>0</sup> = 0, *ξ*<sup>0</sup> = *ξD*(0); Compute Ψ*u*(*θ <sup>q</sup>*) and Ψ*v*(*θ <sup>q</sup>*) from the diffuse approximation approach; **for** *i* = 1..., *Nt* **do** {*u*}({*x*}, *<sup>t</sup> i* , *<sup>θ</sup> <sup>q</sup>*) <sup>←</sup> <sup>Φ</sup>*<sup>r</sup> <sup>u</sup> A*(*t i* ) Ψ*u*(*θ <sup>q</sup>*); {*v*}({*x*}, *<sup>t</sup> i* , *<sup>θ</sup> <sup>q</sup>*) <sup>←</sup> <sup>Φ</sup>*<sup>r</sup> <sup>v</sup> B*(*t i* ) Ψ*v*(*θ <sup>q</sup>*); **end for** <sup>U</sup>(*<sup>θ</sup> <sup>q</sup>*) <sup>←</sup> [{*u*}({*x*}, *<sup>t</sup>* 1, *<sup>θ</sup> <sup>q</sup>*),..., {*u*}({*x*}, *<sup>t</sup>Nt* , *<sup>θ</sup> <sup>q</sup>*)]; <sup>V</sup>(*<sup>θ</sup> <sup>q</sup>*) <sup>←</sup> [{*v*}({*x*}, *<sup>t</sup>* 1, *<sup>θ</sup> <sup>q</sup>*),..., {*v*}({*x*}, *<sup>t</sup>Nt* , *<sup>θ</sup> <sup>q</sup>*)]; Compute Φ(*θ <sup>q</sup>*), Γ(*θ <sup>q</sup>*), *Q*(*θ <sup>q</sup>*), *αD*(*t i* ) and *ξD*(*θi*), *i* = 1, . . . , *Nt*; **while** *t* < *Tf* **do** *<sup>t</sup>* <sup>←</sup> *<sup>t</sup>* <sup>+</sup> *<sup>δ</sup>tROM*; *<sup>t</sup> <sup>i</sup>*+1,*ROM* = *t <sup>i</sup>*,*ROM* + *δtROM*; *αi*+<sup>1</sup> = *α<sup>i</sup>* + *δtROM Q*(*θ <sup>q</sup>*)*ξ<sup>i</sup>* ; Compute *ak*(*t <sup>i</sup>*+1,*ROM*), *k* = 1, . . . , *mv* from the diffuse approximation approach; *ξi*+<sup>1</sup> *<sup>k</sup>* <sup>=</sup> *<sup>p</sup>*(*αi*+1)*<sup>T</sup> <sup>a</sup>k*(*<sup>t</sup> <sup>i</sup>*+1,*ROM*); If needed, reconstruct the high-dimensional displacements/velocity fields: {*u*} {*X*}, *<sup>t</sup> <sup>i</sup>*+1,*ROM*, *θ <sup>q</sup>* = Φ(*θ <sup>q</sup>*) *αi*<sup>+</sup>1; {*v*} {*X*}, *<sup>t</sup> <sup>i</sup>*+1,*ROM*, *θ <sup>q</sup>* = Γ(*θ <sup>q</sup>*) *ξi*<sup>+</sup>1; **end while**

#### **4. Numerical Experiments**

#### *4.1. Study Case*

We consider a capsule flowing in a square-base microchannel of base edges of length 2-. We want to capture the capsule dynamics for capillary numbers *Ca* belonging to the interval [0.005, 0.2] and aspect ratios *a*/ in the interval [0.75, 1.2] for which a steady state shape is reached. The Caps3D code [38,40] is then used as FOM solver. The comparison of the FOM results with experimental ones using a square-base cylinder have been thoroughly described in other previous studies (see for example [23,24,27,39] from A.V. Salsac's research team). The total non-dimensional time for simulation is *T* = 20. For any capillary number and aspect ratio, the capsule is discretized with the same mesh resolution and connectivity, consisting of *Nx* = 2562 nodes (corresponding to 1280 triangular elements), with a capsule mesh size Δ*hC* = 0.075 *a* (see Figure 6). A second-order RK2 Ralston scheme is used for time integration. The dimensionless time step is *γδ*˙ *<sup>t</sup>* <sup>=</sup> <sup>5</sup> · <sup>10</sup>−<sup>4</sup> for *Ca* <sup>&</sup>gt; 0.02 and *γδ*˙ *<sup>t</sup>* <sup>=</sup> <sup>10</sup>−<sup>4</sup> for *Ca* <sup>≤</sup> 0.02 .

**Figure 6.** Three-dimensional representation of a capsule flowing in a square microchannel at *T* = 0.

#### *4.2. FOM Result Database Generation*

A database of FOM results is generated from a sampling of the parameter domain (see Figure 7). It is observed that configurations for which a shape steady state is reached before the non-dimensional final time of 20 correspond to couples (*Ca*, *a*/-) in the parameter plane below the dashed red line of Figure 7. Using a Cartesian parameter sampling with step sizes of 0.01 in *Ca* and 0.05 in *a*/-, plus few additional points at *Ca* = 0.005, we get a database made of *Nc* = 118 configurations. From Caps3D FOM solutions, we pick up timesnapshot solutions every time step Δ*t* = 0.2 in non-dimensional time scale, corresponding to *Nt* <sup>=</sup> 100. This makes a datacube made of 2 <sup>×</sup> <sup>3</sup>*NxNcNt* <sup>≈</sup> 1.81 · 108 double precision float numbers taking about 1.45 GB of memory.

#### Clustering Strategy

For the sake of memory storage complexity, we adopt a strategy of data clustering with two weakly-overlapping clusters chosen manually, represented in Figure 7. For each cluster, a data dimensionality reduction is done following the offline-stage algorithm presented in Section 3. That means that two families of reduced-order models are actually computed. In the online stage, for a new query parameter vector *θ <sup>q</sup>*, one has to determine the cluster of belonging.

**Figure 7.** Design of computer experiment with sampling in the admissible parameter domain. The parameter domain is splitted up into two overlapping clusters: cluster 1 (squares), cluster 2 (crosses) and overlapping region (mixed squares and crosses).

#### *4.3. Elements of Analysis—Accuracy Criteria*

In order to measure the approximation error generated by the data dimensionality process, we introduce the classical Relative Information Content (RIC) (see for example [9]), which is computed as:

$$\text{RIC}(K) = \frac{\sum\_{k=K+1}^{r} \sigma\_k^2}{\sum\_{k=1}^{r} \sigma\_k^2}, \quad K = 1, \dots, r,\tag{50}$$

where *σ*˜ *<sup>k</sup>* is the *k*-th singular value from the SVD decomposition, *r* is the rank of the matrix of study (S or T) and *K* is the truncation rank. A supplementary indicator is the ratio

$$K \mapsto \frac{\partial\_K}{\partial\_1} \tag{51}$$

that gives an idea of the decay rate of the singular values.

The second criterion directly measures the error between the shape predicted by the ROM and the shape computed by the FOM. This is achieved by using the so-called Modified Hausdorff distance *dMH* [44] that we normalize by the capsule radius *a*. The modified Hausdorff distance computes the distance between two finite sets F and G of a normed space of norm ., and is defined as

$$d\_{MH}(\mathcal{F}, \mathcal{G}) = \max(d\_h(\mathcal{F}, \mathcal{G}), d\_h(\mathcal{F}, \mathcal{G})),\tag{52}$$

with

$$d\_h(\mathcal{F}, \mathcal{G}) = \frac{1}{N\_{\mathcal{F}}} \sum\_{p\_{\mathcal{F}} \in \mathcal{F}} d\_s(p\_{\mathcal{F}}, \mathcal{G}) \tag{53}$$

where *N*<sup>F</sup> is the number of points in the set F and *ds*(*p*<sup>F</sup> , G) is the distance between *p*<sup>F</sup> and the set G, which is defined as

$$d\_{\delta}(p\_{\mathcal{F}}, \mathcal{G}) = \min\_{p\_{\mathcal{G}} \in \mathcal{G}} \|p\_{\mathcal{F}} - p\_{\mathcal{G}}\|. \tag{54}$$

#### *4.4. Dimensionality Reduction Analysis*

A singular value decomposition analysis is first performed on the matrices S*<sup>u</sup>* and S*v*, and then on T*<sup>u</sup>* and T*v*. In Figure 8a, we plot the indicator (1 − RIC) (see (50)), as a function of the truncation rank *K*, for S*<sup>u</sup>* and S*v*. What can be seen is that (1 − RIC) rapidly converges towards the value 0 in all cases. An expected (<sup>1</sup> <sup>−</sup> RIC) of 10−<sup>7</sup> is reached for a truncation rank *K<sup>c</sup> <sup>u</sup>* (resp. *K<sup>c</sup> <sup>v</sup>*) of 7 for the displacement (resp. 23 for the velocity). Similarly in Figure 8b, we plot the indicator (<sup>1</sup> <sup>−</sup> RIC) for <sup>T</sup>*<sup>u</sup>* and <sup>T</sup>*v*. The number of modes *<sup>K</sup><sup>x</sup> u* (resp. *K<sup>x</sup> <sup>v</sup>* ) needed to reach the threshold of 10−<sup>7</sup> is 7 for the displacement (resp. 56 for the velocity).

**Figure 8.** Behaviour of the relative information content of the matrices S*<sup>u</sup>* and S*<sup>v</sup>* (**a**) and T*<sup>u</sup>* and T*<sup>v</sup>* (**b**) shown in the form (<sup>1</sup> <sup>−</sup> RIC) as a function of the truncation rank *<sup>K</sup>*. The horizontal red line corresponds to (<sup>1</sup> <sup>−</sup> RIC) = <sup>10</sup>−7.

As supplementary indicators, the singular values *σ*˜*<sup>K</sup>* normalized by *σ*˜1 are plotted in Figure 9a (resp. Figure 9b) for both matrices S*<sup>u</sup>* and S*<sup>v</sup>* (resp. T*<sup>u</sup>* and T*v*) in log10 scale. One can first observe a lower decay rate for the velocity fields compared to the displacements, meaning a greater information complexity for the velocity. Secondly, the decay rate is lower for the global spatial mode than for the parametric modes, indicating a larger entropy of information on the whole parameter domain. That justifies the derivation of suitable lower order spatial basis at a query parameter *θ <sup>q</sup>* in the online stage.

At the beginning of the online stage, for a query parameter *θ <sup>q</sup>*, an interpolated approximate solution is computed thanks to a diffuse approximation reconstruction. This allows us to get pseudo-snapshots in time for both displacements and velocities, stored in matrices U(*θ <sup>q</sup>*) and V(*θ <sup>q</sup>*), respectively. We assess the RIC for the two matrices, from an experimental parameter vector *θ <sup>q</sup>* = (0.10, 0.90). The comparison of the time evolution of POD coefficients between FOM and ROM models shows a high accuracy (see Figure 10). and Figure 11 shows that the RIC rapidly converges to 1. An expected RIC greater than <sup>1</sup> <sup>−</sup> <sup>10</sup>−<sup>7</sup> returns a truncation rank *mu* (resp. *mv*) of value 3 (resp. 8).

**Figure 9.** (**a**) Parametric normalized singular values *σ*˜*K*/*σ*˜1 for S*<sup>u</sup>* and S*v*; (**b**) Spatial normalized singular values *σ*˜*K*/*σ*˜1 for T*<sup>u</sup>* and T*v*.

**Figure 10.** FOM versus ROM comparison of the time evolution of the first three displacement (**a**) and velocity (**b**) POD coefficients for the query parameter *θ <sup>q</sup>* = (0.10, 0.90).

**Figure 11.** Online stage: behaviour of the relative information content of the matrices U(*θ <sup>q</sup>*) and V(*θ <sup>q</sup>*) shown in the form (1 − RIC) for query parameter *θ <sup>q</sup>* = (0.10, 0.90). The red line corresponds to (<sup>1</sup> <sup>−</sup> RIC) = <sup>10</sup>−7.

#### *4.5. ROM Accuracy Analysis*

The reduced-order model algorithm is applied with the following parameters and options:


The resulting time-evolution of the three-dimensional capsule shape, that is reconstructed with the ROM model, is illustrated in Figure 12 for the query couple *θ <sup>q</sup>* = (0.10, 0.90). The steady-state is reached before *γ*˙*t* = 3, which explains that the capsule shape is the same for *γ*˙*t* = 3, 6, 9.

We now focus on the accuracy analysis of the proposed reduced-order model. The methodology for error measurement is based on a '*Leave-one-out*' cross-validation procedure, where each sample FOM solution is taken out from the database and then evaluated by the ROM model and compared to the original FOM one. The error is measured using the modified Hausdorff distance calculated on the capsule shapes at different instants.

Figure 13 shows the heat maps of the FOM-vs-ROM error computed over the parameter space at the time instants *γ*˙*t* = 1, 2, 4 and 8. Figure 13 shows that the predicted ROM solutions are very accurate with a mean relative error below 0.2%. The maximum relative errors are below 3.5%: they occur along the boundary of the parameter domain, which is the only location where the predictions slightly lose in accuracy. This is probably due to a lack of well-distributed neighbors close to the boundaries, which affects the accuracy of the DA reconstruction (off-centre approximation). One can also notice that the accuracy of predictions decreases in time.

**Figure 13.** Heat maps of the normalized Hausdorff Distance *dMH*/*a* of configuration prediction shapes over the parametric space at different transient states: (**a**) *γ*˙*t* = 1; (**b**) *γ*˙*t* = 2; (**c**) *γ*˙*t* = 4; and (**d**) *γ*˙*t* = 8. The maximum error is 3.26% in (**d**).

The capsule cross-section profiles predicted by the ROM (red dots) are compared to the reference FOM solutions (solid black line) in Figure 14 at different time instants (*γ*˙*t* = 0, 1, 2 and 8) for the 6 configurations, selected as illustration in Figure 3. We observe that the reduced-order model returns very accurate solutions in terms of capsule shape as well-as axial position in the channel.

From the computing performance point of view, ROM-vs-FOM speedups are observed to be of order 10,000 with almost the same accuracy, making interactive exploration and real-time visual rendering possible.

#### *4.6. CapsuleExplorer: Capsule Visualization/Exploration Software*

We have developed an in-house software tool CapsuleExplorer based on the proposed ROM to provide the three-dimensional microcapsule deformation/evolution at any time *γ*˙*t* and for any *θ <sup>q</sup>* in the admissible parameter domain. CapsuleExplorer allows one to select a particular couple (*Ca*, *a*/-) in the admissible parameter domain, then to visualize the capsule dynamics between initial and final times, either in three dimensions or two dimensions with longitudinal or transversal cross-sectional view. The ROM high performance feature allows real-time exploration/visualization. CapsuleExplorer has been developed as a web application. Figures 15 and 16 show some screenshots of the graphics user interface, which will be useful for applications such as identifying the capsule wall mechanical properties through comparison with experimental results.

**Figure 14.** Comparison between the ROM (red dots) and FOM solutions (black line) of the capsule cross-section shapes in the plane *y* = 0 at the times *γ*˙*t* = 0, 1, 2 and 8 respectively, for the 6 parameter couples selected in Figure 3. The horizontal lines correspond to the channel walls.

**Figure 15.** CapsuleExplorer: (**a**) parameter domain exploration; (**b**) dynamic 3D capsule view.

**Figure 16.** CapsuleExplorer: (**a**) dynamic 2D cross-section longitudinal view; (**b**) dynamic 2D cross-section transversal view.

#### **5. Concluding Remarks**

In this paper, we have presented an innovative data-driven reduced-order model that enables the dynamics of a deformable membrane flowing in a microchannel, from its initial state to the steady shape state. The ROM is built to be valid in a large domain of interest in the parameter plane (*Ca*, *a*/-). Our FSI-ROM model first starts with an offline procedure to build two global orthonormal bases (space+parameter) that return good approximations of the FOM solutions over the whole parameter domain. The rather small truncation ranks already lead to an appreciable data dimensionality reduction, which is important for complexity and memory storage purposes.

The online stage consists in predicting the space-time solution for any query couple *θ <sup>q</sup>* = (*Ca* , *a*/-) in the parameter domain. In a first step, we determine a low-order basis for both the displacement and velocity vector variables. This is achieved by the use of diffuse approximation that returns an interpolated space-time solution at the query vector *θ <sup>q</sup>*. Then an SVD analysis provides a suitable low-order spatial basis for final construction of the ROM in the second step. The physically-based ROM is made of the kinematics equation and the law of membrane quasi-static equilibrium in their reduced formulation. The unknown variables become the POD coefficient vectors of displacement and velocity fields. The reduced quasi-static equilibrium law is determined once again by the use of a diffuse approximation. The manifold learning is achieved by the use of time-snapshot data of the interpolated solution at *θ* = *θ <sup>q</sup>*.

Numerical experiments confirm the efficiency of the method. ROM-vs-FOM speedups are observed to be of order 10,000 with almost the same accuracy (with less than a 0.3% error measured in terms of Hausdorff distance inside the parameter domain). Larger errors are encountered at the boundary of the parameter domain, but they still remain reasonable (up to 3.3% in Hausdorff distance). This work tends to show that model-order reduction techniques are complementary and valuable tools for the rapid design and optimization of capsules in healthcare engineering such as drug delivery through blood vessels.

The case of more complex FSI configurations such as the deformation of capsules flowing through a bifurcated microchannel will be investigated in a future work.

**Author Contributions:** Conceptualization, T.B., P.V., F.D.V. and A.-V.S.; Formal analysis, F.D.V., P.V.; Funding acquisition, A.-V.S.; Investigation, T.B., F.D.V. and A.-V.S.; Methodology, T.B., P.V. and F.D.V.; Project administration, A.-V.S.; Resources, C.Q.-G., C.D. and A.-V.S.; Software, T.B., C.Q.-G., C.D. and A.-V.S.; Supervision, F.D.V. and A.-V.S.; Validation, T.B.; Visualization, T.B.; Writing—original draft, T.B.; Writing—review & editing, C.Q.-G., C.D., P.V., F.D.V. and A.-V.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. ERC-2017-COG—MultiphysMicroCaps).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

