Data-driven kinematics-consistent model order reduction of fluid-structure interaction problems: application to deformable microcapsules in a Stokes flow

In this paper, we present a generic approach of a dynamical data-driven model order reduction technique for three-dimensional fluid-structure interaction problems. A low-order continuous linear differential system is identified from snapshot solutions of a high-fidelity solver. The reduced order model (ROM) uses different ingredients like proper orthogonal decomposition (POD), dynamic mode decomposition (DMD) and Tikhonov-based robust identification techniques. An interpolation method is used to predict the capsule dynamics for any value of the governing non-dimensional parameters that are not in the training database. Then a dynamical system is built from the predicted solution. Numerical evidence shows the ability of the reduced model to predict the time-evolution of the capsule deformation from its initial state, whatever the parameter values. Accuracy and stability properties of the resulting low-order dynamical system are analyzed numerically. The numerical experiments show a very good agreement, measured in terms of modified Hausdorff distance between capsule solutions of the full-order and low-order models both in the case of confined and unconfined flows. This work is a first milestone to move towards real time simulation of fluid-structure problems, which can be extended to non-linear low-order systems to account for strong material and flow non-linearities. It is a valuable innovation tool for rapid design and for the development of innovative devices.


Introduction
Fluid-structure interaction (FSI) problems often occur in Engineering (aircraft and automotive industries, wind turbines) as well as in medical applications (cardiovascular systems, artificial organs, artificial valves, medical devices, etc.). Today the design of such systems usually requires advanced studies and high-fidelity (HF) numerical simulations become an essential tool of computed-aided analysis. However, computational FSI is known to be very time-consuming even on high-performance computing facilities. Usually, engineering problems are parameterized and the search of suitable designs require numerous computer experiments leading to prohibitive computational times. For particular applications such as the tracking of drug carrier capsules flowing in blood vessels, it would be ideal to have real-time simulations for a better understanding of the behaviour of the dynamics and for efficiency assessment. Unfortunately, today high-fidelity real-time FSI simulations are far from being reached with current High Performance Computing (HPC) facilities.
A current trend is to use machine learning (ML) or artificial intelligence (AI) tools such as artificial neural networks (ANN). Such tools learn numerical simulations from HF solvers and try to map entry parameters with output criteria in an efficient way, with response times far less than HF ones, say 3 or 4 orders of magnitude smaller. In some sense, heavy HF computations and training stage are done in an offline stage, and learned ANNs can be used online for real time evaluations and analysis. However, ML and ANN today are not fully satisfactory for dynamical problems, and/or the training stage itself may be time consuming, thus requiring more Central Processing Unit (CPU) time. Another option is the use of model order reduction (MOR). Reduced-order modeling (ROM) can be seen as a 'grey-box' supervized ML methodology, taking advantage of the expected low-order dimensionality of the FSI mechanical problem. By 'grey-box' we mean that the low-dimensional encoding of the ML process is based on mechanical principles and a man-made preliminary dimensionality reduction study. This allows one for a better control of the ROM accuracy and behaviour. There are two families of MOR: intrusive and non-intrusive approaches. The intrusive approaches use physical equations. The low-order model is derived by setting the physical problem on a suitable low-dimensional space. The accuracy can be very good, but the price to pay is the generation of a new code which can be a tedious and long task. The non-intrusive approach does not require heavy code development. It is based on HF simulation results used as entry data. Although it is not based on high-fidelity physical equations, a non-intrusive approach can include a priori physical informations, like e.g. meaningful physical features, prototype of system of equations, pre-computed principal components, consistency with physical principles, etc.
In the recent literature, efficient intrusive ROMs for FSI have been proposed e.g. in (Quarteroni et al. 2016). But to our knowledge there are far less contributions in non-intrusive ROMs dedicated to FSI.
In this paper, we propose a data-driven model order reduction approach for FSI problem which is consistent with the equations of kinematics and is designed from a low-order meaningful system of equations. As case of study, we focus on the motion of a microcapsule, a droplet surrounded by a membrane, subjected to a confined and unconfined Stokes flow.
Artificial microcapsules can be used in various industrial applications such as in cosmetics (Miyazawa et al. 2000;Casanova & Santos 2016), food industry (Yun et al. 2021) and biotechnology, where drug targeting is a high potential application (Ma & Su 2013;Abuhamdan et al. 2021;Ghiman et al. 2022). Once in suspension in an external fluid, capsules are subjected to hydrodynamics forces, which may lead to large membrane deformation, wrinkle formation or damage. The numerical model must be able to capture the time-evolution of the nonlinear 3D large deformations of the capsule membrane. Different numerical strategies are possible to solve the resulting large systems of equations (Lefebvre & Barthès-Biesel 2007;Hu et al. 2012;Ye et al. 2017;Tran et al. 2020). However, they all have long computational times.
Different approaches have been used over the past decade to accelerate the computations, such as HPC (e.g. Zhao et al. (2010)) and Graphics Processing units (e.g. Matsunaga et al. (2014)). More recently, reduced order models have been proposed to predict the motion of capsules suspended in an external fluid flow. In Quesada et al. (2021), the authors used the large amount of data generated by numerical simulations to show how relevant it is to recycle these data to produce lower-dimensional problem using physics-based reduced order models. However, their method can only predict the steady-state capsule deformed shape. Boubehziz et al. (2021) show for the first time the efficiency of data-driven model-order reduction technique to predict the dynamics of the capsule in a microchannel. However, the method is cumbersome as it requires two POD bases, one to predict the velocity field, the other to capture the shape evolution over time. And then they reconstruct the solution in the parameter space thanks to diffuse approximation (DA) strategy. The proposed method serves different objectives. We have designed the method to be non intrusive for practical uses of existing high-fidelity FSI solver (also referred to as the Full-Order Model, or FOM). That means that the ROM methodology should be data-driven. We also want the ROM to be consistent with the equations of kinematics. The model must thus return the displacement { } and velocity { } fields from a few snapshots provided by the FOM. It must otherwise be able to predict the solution for any parameter vector in predefined admissible domain. Finally, the kinematics-consistent data-driven reduced-order model of capsule dynamics must ideally open the way to real-time simulations. To do so, we use a coupling between Proper Orthogonal Decomposition (POD) and Dynamic Mode Decomposition (DMD), as well as a Tikhonov regularization for robustness purposes and interpolation to predict solutions for any parameter value.
As indicated above, we mainly consider the case of an initially spherical capsule flowing in a microfluidic channel with a square cross-section. The corresponding FOM was developed by Hu et al. (2012) and used to get a complete numerical database of the three-dimensional capsule dynamics as a function of the parameters of the problem: the capsule-to-tube confinement ratio, hereafter referred to as size ratio /ℓ and the capillary number , which measures the ratio between the viscous forces acting onto the capsule membrane and the membrane elastic forces. For clarity reasons, different ROMs are introduced with increasing levels of generality, as detailed in Table 1. First, we consider a fixed parameter vector, and get a space-time ROM in the form of a low-order dynamical system. Next, we generate such ROMs for the parameter samples that fill the admissible parameter domain, and then assess the uniform accuracy (space-time accuracy over the whole sample set). Finally, we propose a strategy to derive a general space-time-parameter ROM for any value of the parameter vector ( , /ℓ) in the admissible space.
The paper is organized as follows. First, we present the physics of the problem and the FOM in Section 2. The strategy used to develop a non-intrusive space-time ROM is detailed in Section 3. The results are presented for a given configuration in Section 4. The accuracy is then estimated in Section 5 on the entire database, formed by all the cases that have reached a stationary state. Finally in Section 6 we present the methodology of space-time-parameter ROM. The ROM accuracy is confirmed by some numerical experiments.

Full-order microcapsule model, parameters and quantities of interest
2.1. Problem description An initially spherical capsule of radius flows within a long microfluidic channel having a constant square section of side 2ℓ (Figure 1). The suspending fluid and capsule liquid core are incompressible Newtonian fluids with the same kinematic viscosity .
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 modeled as a surface devoid of bending stiffness with surface shear modulus . The two non-dimensional governing parameters of the problem are the size ratio /ℓ and the capillary number = / (2.1) where 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 ( = 1) and internal fluids ( = 2), together with the membrane equilibrium equation to determine the dynamics of the deformable capsule within the microchannel.
For the fluid problem, we denote ( ) , ( ) and ( ) the velocity, stress and pressure fields in the two fluids. These parameters are non-dimensionalized using ℓ as characteristic length, ℓ/ as characteristic time and ℓ as characteristic force. The non-dimensional Stokes equations ∇ ( ) = ∇ 2 ( ) , ∇ · ( ) = 0, = 1, 2 (2.2) are solved in the domain bounded by the cross sections at the tube entrance and at the exit. These cross sections are assumed to be both located far from the capsule. The Focus on Fluids articles must not exceed this page length reference frame ( , , , ) is centered at each time step on the capsule center of mass in the high-fidelity code, but the displacement of the capsule center of mass along the tube axis is computed. The boundary conditions of the problem are the following ones: • The velocity field is assumed to be the unperturbed flow field on and , i.e. the flow disturbance vanish far from the capsule.
• The pressure is uniform on and .
• A no-slip boundary condition is assumed at the channel wall and on the capsule membrane : (2.3) • The normal load on the capsule membrane is continuous, i.e. the non-dimensionalized external load per unit area exerted by both fluids is due to the viscous traction jump: where is the unit normal vector pointing towards the suspending fluid.
To close the problem, the external load on the membrane is deduced from the local equilibrium equation, which, in absence of internia, can be written as where is the non-dimensionalized Cauchy tension tensor (forces per unit arclength in the deformed plane of the membrane) and ∇ · 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) where 1 and 2 are the principal extension ratios measuring the in-plane deformation.

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. 2010;Hu et al. 2012). Thanks to the latter, only the boundaries of the flow domain, i.e the channel entrance and exit , 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 2 triangular elements and 2562 nodes, which correspond to a characteristic mesh size Δℎ = 0.075 . The channel mesh of the entrance surface and exit surface 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 1 triangular elements and 1905 nodes. At time = 0, a spherical capsule is positioned with its center of mass on the channel axis. At each time step, the in-plane stretch ratio 1 and 2 are computed from the nodes deformation. The elastic tension tensor is then deduced from the values of 1 and 2 . The finite element method is used to solve the weak form of the membrane equilibrium equation (2.5) and determine the external load .
Applying the boundary integral method, the velocity of the nodes on the capsule membrane reads (Pozrikidis 1992): (2.7) for any in the spatial domain when the suspending and internal fluids have the same viscosity. The vector is the disturbance wall friction due to the capsule, Δ is the additional pressure drop and = − .
To update the position of the membrane nodes, the nodal displacement is computed by integrating equation (2.3) in time. The procedure is repeated until the desired nondimensional time /ℓ. 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: where is a nonlinear mapping from R 3 to R 3 and is the number of nodes on the membrane. Regarding time discretization, a Runge-Kutta Ralston scheme is used: where Δ > 0 is a constant time step and { } and { } respectively represent the discrete membrane displacement field and the discrete membrane velocity field at discrete time = Δ . The initial condition is simply The whole numerical scheme is subject to some Courant-Friedrichs-Lewy (CFL) type stability condition on the time step (Walter et al. 2010) because of its explicit nature. The numerical method is conditionally stable if the time step Δ satisfies (2.10) From the computational point of view, the resolution of (2.9) at each time step requires i) the computation of the disturbance wall friction at all the wall nodes, ii) the additional pressure drop Δ , iii) the traction jump 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 /ℓ = 0.7, considering the mesh discretization described above and a workstation equipped with 2 processors Intel ® Xeon ® Gold 6130 CPU (2.1 GHz). A week of computation is sometimes necessary to simulate the dynamics of an initially spherical capsule in a microchannel over the non-dimensional time /ℓ = 10. For that reason, a model-order reduction (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 (RB) of suitable functions, then derive a low-order system of equations.
In the case of differential algebraic equations (DAE) like (2.8)-(2.9), the reduced system of equations to find should also be of DAE nature. Remark that it is often possible to reformulate DAEs as a system of ordinary differential equations (ODEs) (Ascher & Petzold 1998). In the next section, we give details on the chosen ROM methodology for the particular case and context of FSI capsule problem.

Non-intrusive space-time model-order reduction strategy
In this section, the parameter couple = ( , /ℓ) is fixed, thus we omit the dependency of the solutions with respect to for the sake of simplicity. For the derivation of the ROM model, we consider the semi-discrete time-continuous version of the FOM, i.e. (2.8)-(2.9).

Dimensionality reduction and reduced variables for displacements and velocities
Assume first that, for any ∈ [0, ], the discrete velocity field can be accurately approximated according to the expansion for some orthonormal modes { } ∈ R and real coefficients ( ). The truncation rank is of course expected to be far less than as expected in a general ROM methodology. From the kinematics equations we have so that the displacement field can be accurately represented by The coefficients ( ( )) and ( ( )) 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 where the exponent denotes the transpose of the matrix. The condensed matrix forms of (3.2) and (3.1) respectively are Since the modes { } are assumed to be orthonormal (for the standard Euclidean inner product), the matrix is a semi-orthogonal matrix, i.e. = . In particular, we have ( ) ≈ { }( ) and ( ) = { }( ). Note that the modes { } and reduced variables , are determined for each parameter set ( , /ℓ), but a common value of the truncation rank is chosen for all the sets. Its practical computation will be detailed in a next subsection, as well as that of the modes { }.

ROM prototype
The expressions {˜ }( ) = ( ) and {˜ }( ) = ( ) provide low-order representations of displacement and velocity fields respectively. We can now write equations for the reduced vectors ( ) and ( ) respectively. In this subsection, let us consider a projection Galerkintype approach. Let us denote ., . the standard Euclidean scalar product in R . Considering a test vector { } in = ({ 1 }, ..., { }), we look for an approximate displacement field{ }( ) solution of the projected kinematics equations By considering each test vector { } = { }, we get the consistent reduced kinematics equation = .
(3.4) Consider now the projected field {˜ }( ) which is solution of the system of algebraic equations (Galerkin approach): (3.5)

Again by taking the test vector
.
It is in the form ( ) = ( ( )) (3.6) with the mapping : R → R defined by ( ) = ( ). We get a reducedorder algebraic equilibrium equation. Unfortunately, because of the nonlinearities in , the computation of ( ) requires high-dimensional ( ) 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 (EIM) ( Barrault et al. 2004) 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 derivating (3.6) with respect to time, we get Thanks to the reduced kinematics equation (3.4), we get Since is hard to evaluate, it is even harder to evaluate its differential. But the differential ( ( )) can be approximated itself, or replaced by some matrix ( ). Then we get a ROM structure (ROM prototype) in the form (3.9) The differential system (3.8)-(3.9) is linear with variable coefficient matrix ( ) ∈ M (R).
It can be written in matrix form (3.10) The spectral properties of the differential system (3.10) are related to the spectral properties of matrix ( ). In particular, if all the (complex) eigenvalues ( ) of ( ) are such that Re( ( )) < 0 for all (uniformly distributed in time), then the system is dissipative.
3.3. Nonintrusive approach, SVD decomposition and POD modes One of the requirements of this work is to achieve a non-intrusive reduced-order model, meaning that the effective implementation of the ROM does not involve tedious low-level code development into the FOM code. For that, a data-based approach is adopted: from the FOM code, it is possible to compute FOM solutions ({ } , { } ) at discrete times , = 0, ..., ( = Δ = ), then store some snapshot solutions (called snapshots) into a database for data analysis and later design of a ROM. Proper Orthogonal Decomposition (POD) (Berkooz et al. 1993) is today a well-known dimensionality reduction approach to determine the principal components from solutions of partial differential equations. The Sirovich's snapshot variant approach (Sirovich 1987) is based on snapshot solutions from a FOM to get a posteriori empirical POD modes { }. 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 0. From SVD theory, for a given accuracy threshold > 0, the truncation rank = ( ) is computed as the smallest integer such that the inequality holds (Shawe-Taylor & Cristianini 2004). Proceeding like that, it is shown that the relative orthogonal projection error of the snapshots { } onto the linear subspace spanned by the first eigenvectors of is controlled by . Denoting the linear orthogonal projection operator over , we have: The matrix is obtained as the restriction of to its first columns.
3.4. Data-driven identification of coefficient matrix The system (3.8)-(3.9) is still not closed since the coefficient matrices ( ) 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 ( ) is searched as a time-constant matrix . In this case, equation (3.9) becomes ( ) = ( ). This is the scope of this article. From the time continuous problem, one could determine the matrix by minimizing the least square functional But practically, we only have velocity snapshot data 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 S = [{ } 1 , ..., { } ], we compute first the reduced snapshots variables: Next, we determine a matrix that minimizes the least square cost function: In (3.13), the finite difference +1 − Δ is a first-order approximation (in Δ ) of at time . In appendix A, we provide a mathematical analysis of the effect of time discretization in (3.13) about the impact on the stability of the resulting identified differential system compared to the initial one.
The minimization problem (3.13) can be written in condensed matrix form with the two data matrices Because X and Y store reduced variables (of size ), for a sufficient number of discrete snapshot times, these two matrices are horizontal ones. We will assume that the rank of X is always maximal, i.e. equal to . The least-square solution of (3.14) is then given by 3.5. Tikhonov least-square regularized formulation From standard spectral theory arguments, it is expected that the POD coefficients rapidly decay when increases as soon as both displacement and velocity fields are smooth enough. A possible side effect is the bad condition number of the matrix X, since the last rows of X have small coefficients (thus leading to row vectors close to zero 'at the scale' of the first row of X). Even if the solution in (3.16) 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-square problem (3.14) by adding a Tikhonov regularization term (see e.g. (Aster et al . 2019) where the scalar > 0 is the regularization coefficient. The factor X 2 in the regularization term has been added for scaling purposes. The solution of (3.17) is given by (3.18)

Choice of optimal regularization coefficient
Of course, the solution matrix depends on the regularization coefficient and one can ask what is the optimal choice for . There is a trade-off between the approximation quality measured by the residual Y − X and the norm solution . The minimization of should ensure that unneeded features will not appear in the regularized solution. When plotted on the log-log scale, the curve of optimal values ↦ → versus the residual ↦ → Y − X often takes on a characteristic L shape (Aster et al. 2019). A design of experiment with the test of different values of (starting say from 10 −12 to 10 −5 ) generally allow to find quasi-optimal values of located at the corner of the L-curve, thus providing a good trade-off between the two criteria.
3.6. Reduced-order continuous dynamical system Once the matrix has been determined, we get the reduced-order continuous dynamical system . .
The exact analytical solution of (3.21) is The stability of the differential system depends on the spectral structure of A , or equivalently on the spectrum of . Because of the stability of the fluid-capsule coupled system and from accurate solutions of the FOM solver, one can hope that the solution of the least-square 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 ( + )/2.

Model consistency with steady states
A steady state in our context is defined by a capsule that reaches a constant velocity { } ∞ , so that the motion becomes a translation flow in time in the direction { } ∞ . From (3.1), this shows that ( ) also reaches a constant vector ∞ , and = 0 at steady state. As a consequence, from (3.20), we get ∞ = 0, meaning that 0 is an eigenvalue of with ∞ as eigenvector. As a conclusion, the matrix must have zero in its spectrum in order to be consistent with the existence of steady states.

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 Δ gives (3.24) By multiplying (3.23) by we get the space-time approximate solution so the ROM model is completely consistent with the kinematics equation. Stability properties of the discrete system are linked to the spectral properties of the matrix For unconditional stability in time, it is required for the eigenvalues of + Δ 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.

Accuracy criteria and similarity distances between ROM and FOM solutions
In order to quantify the error induced by approximations, we introduce 3 accuracy criteria. The first accuracy criterion is the relative information content (RIC), defined by quantifies the relative amount of neglected information when truncating the number of modes at rank . The truncation rank is determined such that the RIC is inferior to the accuracy threshold . The accuracy threshold is fixed to 10 −6 .
The second accuracy criterion is the relative time residual R. It quantifies the relative error induced by the minimization of the least square cost function (3.13) using . It is given by where X represents the ℎ column of X and Y the ℎ column of Y. The index is thus linked to the snapshots ( ∈ {1, ..., }). 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 /ℓ hereafter. The third accuracy criteria Shape ( /ℓ) measures the difference between the 3D reference capsule shape given by the FOM (S FOM ) and the 3D shape predicted by the ROM (S ROM ). It is defined at a given non-dimensional time /ℓ as the ratio between the modified Hausdorff distance (MHD) computed between S and S and non-dimensionalized by ℓ Shape ( /ℓ) = MHD(S FOM ( /ℓ), S ROM ( /ℓ)) ℓ The modified Hausdorff distance is the maximum value of the mean distance between S FOM and S ROM and the mean distance between S ROM and S FOM (Dubuisson & Jain 1994).

Numerical experimentation on a given configuration
The method is first applied to a given configuration, in order to set the model parameters and to study its stability and precision. We consider the dynamics of an initially spherical capsule flowing in a microchannel when = 0.17 and /ℓ = 0.8. The time step between each snapshot Δ equals to 0.04. The dynamics predicted by the FOM is illustrated in Fig.  3 up to a nondimensional time /ℓ = 10. As the capsule flows, its membrane is gradually deformed by the hydrodynamic forces inside the channel during a temporary time until a steady state is reached. We assume that the capsule has reached its steady-state shape, when the surface area of the capsule varies by less than 5 × 10 −4 × (4 2 ) over a non-dimensional time /ℓ = 1. For ( = 0.17, /ℓ = 0.8), the steady state is reached at /ℓ = 2.6 and is characterized by a parachute capsule shape (Figure 3).

Proper orthogonal decomposition, truncation and modes
The singular value decomposition 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 about 1% only with one mode. The more modes  is kept, the less information is neglected. In the following, we fix the number of modes to 20. The accuracy threshold is thus equal to 10 −6 . 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 ( = 0.17, /ℓ = 0.8). The first six modes are mostly dedicated to change 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 .

Dynamic Mode Decomposition: empirical regularization
Before determining the matrix , we check the condition number of the matrices X and XX T . They are respectively equal to 7.9 × 10 4 and 6.2 × 10 9 . The condition numbers of these matrices are very high and the matrix , determined by solving (3.16), may be sensitive to perturbations or noise. To improve the robustness, a Tikhonov regularization is applied to solve the least-square problem (3.13) and the matrix is computed using (3.18), which depends on the regularization coefficient . To determine the optimal value of , the relative least square error X − Y / Y is represented according to the norm solution when 20 modes are considered and when is varied between 10 −12 and 10 −5 (Figure 6). The least square error X − Y and the norm solution are minimal when = 10 −9 . In the following, is thus fixed to = 10 −9 and the number of modes to 20.

Validity check of the ROM: spectral study of the resulting matrix
In order to detect anomalies, a spectral analysis of the reduced-order model learned by the DMD method is carried out. The spectrum of the matrix is represented in Figure 7. All the eigenvalues of the matrix have non-positive real parts. The resulting linear ROM is thus stable.
The temporal evolution of the residual R (Figure 8) shows that the error is less than 4%. The maximal value is reached at the beginning of the simulation ( /ℓ < 0.3) and R decreases afterwards. When /ℓ 2, 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 /ℓ > 6, is of order 0.01%. The  stationary state is thus well predicted by the model and the error during the transient stage is more than acceptable.

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 ( = 0.17, /ℓ = 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 perfectly overlap at /ℓ = 0, 0.4, 2, 4, 6. The temporal evolution of Shape is shown in Figure 11a. The maximum value of the error committed on the 3D shape Shape equals to 0.022%. The error on the capsule shape Shape is thus negligible. The deformation of the capsule from its initially spherical shape to its steady state over an non-dimensional time /ℓ = 10 can thus be estimated very precisely with the developed reduced-order model.  4.5. Sensitivity of the ROM on the learning time The DMD method predicts the capsule displacement at time +1 from that at time . The model has been constructed until now by considering the dynamics of the capsule over a non-dimensional time /ℓ of 10.
In order to study its influence, we increase the learning , i.e. the non-dimensional time over which the model is trained, from 2 to 8 and estimate the capsule dynamics using the ROM up to a non-dimensional time /ℓ of 10. The number of modes is always equal to 20 and = 10 −9 . The estimated shape is then compared to the one simulated with the FOM. The time evolution of Shape is shown in Figure 11a. The error between the 3D shape increases, as soon as we exceed the learning time. The longer the learning time, the smaller the error Shape at /ℓ = 10 (Figure 11b). The comparison of the capsule profile estimated by the ROM at /ℓ = 10 with the one simulated with the FOM (Figure 11c) confirms that considering only the dynamics of the capsule up to a non-dimensional time /ℓ of 2 is not sufficient. In fact, the capsule is still in the transient phase at /ℓ = 2. When = 4, a small difference persists in the parachute at /ℓ = 10 and Shape = 1.5%. However, when the learning time is increased above 6, Shape falls below 0.19% at /ℓ = 10. This learning time is sufficient to estimate the overall shape of the capsule.

Space-time ROM accuracy assessment over the full parameter sample set
The capillary number and the aspect ratio /ℓ are now considered as variable parameters. A database of 119 simulations of the deformation of an initially spherical capsule in a microchannel has been generated using the FOM with the same time step and mesh size as in section 4. Figure 12 shows have been computed to create the training database. When the capsule initial radius is close to or larger than the microchannel cross-dimension ( /ℓ 0.90), the capsule is predeformed into a prolate spheroid to fit in the channel. For a given /ℓ, a limit value of exists beyond which a capsule does not reach a steady-state ( Figure 12). This is due to the softening behavior of the neo-Hookean law. For all the couples ( , /ℓ) of the training database, the capsule shape is reconstructed from the ROM results at given non-dimensional times and compared to the shape predicted by the FOM at the same non-dimensional time. The evolution of the error committed on the capsule shape ℎ on the full database is illustrated in Figure 13 at /ℓ = 0, 0.4, 1, 2, 5, 10. ℎ is null at /ℓ = 0. The ROM is therefore able to predict the initial capsule shape correctly, whether it is spherical or slightly ellipsoidal. Until /ℓ 2, ℎ essentially remains zero on the majority of the database. Otherwise, it is equal to 0.15% at maximum. At /ℓ = 5 and 10, the error ℎ slightly increases for most of the couples ( , /ℓ) of the database. It remains fully acceptable since it is equal to 0.35% at maximum. When considering 20 modes and = 10 −9 , the developed ROM allows us to estimate with great precision the dynamics of an initially spherical capsule in a microchannel with a square cross-section.
To respect the stability condition (see Equation 2.10), the time step imposed to simulate the capsule dynamics with the FOM decreases, when the decreases. The lower the , the longer the simulation lasts (Figure 2). The time needed to calculate the capsule shape and write the results was estimated on the same workstation used to simulate and generate the result files with the FOM (2-CPU Intel ® Xeon ® Gold 6130, 2.1 GHz). The speedup is the ratio between the FOM runtime and the ROM runtime. Its evolution according to the FOM time step is illustrated in Figure 14. It was estimated from the ROM and FOM simulation time obtained when /ℓ = 0.7. The speedup varies between 52106 for a FOM time step of 10 −4 (i.e for the lowest value of tested) and 4200 for 5 × 10 −4 (i.e 0.05). It is thus possible to estimate the capsule dynamics very precisely with the developed ROM, while considerably reducing the computational time.
Another significant advantage is the gain in storage of the simulation results. By storing only the reduced variables , , the modes { } and the initial position of the nodes of each couple = ( , /ℓ), the training database is reduced from 1.9 GB, when computed with the FOM, to 0.15 GB with the ROM. It can therefore be more easily shared.

Full space-time-parameter ROM (for any admissible parameter value)
6.1. General methodology It is here again assumed that a training database of precomputed FOM results is available. Now we would like to derive a ROM for any parameter couple = ( , /ℓ) 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 { }( ; ) 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 Section 3. Below we give a detailed explanation of the two steps.
Step 1: predictor step. Considering a parameter couple , we first search the three nearest neighbor parameters in the sample set that form a nondegenerate triangle in the plane ( , /ℓ). Let us denote them by 1 , 2 and 3 . We will define a linear operator in the triangle ( 1 , 2 , 3 ). For that, let us introduce the barycentric coordinates ( 1 , 2 , 3 ), ∈ [0, 1], = 1, 2, 3 such that 1 + 2 + 3 = 1, (6.1) 1 1 + 2 2 + 3 3 = . (6.2) The 3 × 3 linear system (6.1),(6.2) is invertible as soon as the triangle ( 1 , 2 , 3 ) is nondegenerate. Notice that the ( = 1, 2, 3) are actually functions of . Let us now denote by { 1 }, { 2 } and { 3 } the displacement fields for the parameter vectors 1 , 2 and 3 respectively. Then we can consider the predicted velocity fieldˆ ( ; ) defined by (6.3) 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 Section 3 can be applied to these data to get a reduced dynamical system in the form We also have a matrix ( ) of orthogonal POD modes and we can go back to the highdimensional physical space by the standard operations  6.2. Numerical experiments, ROM accuracy assessment A testing database is created using the FOM as in Section 5 and considering ( , /ℓ)couples which are not in the training database. A set of 110 ( , /ℓ)-couples are included in this database ( Figure 15). For all the ( , /ℓ)-couples of the testing database, the capsule dynamics is interpolated from the dynamics of the 3 closest neighbors at a given nondimensional time. Capsule shapes obtained by the ROM are compared to the ones predicted by the FOM at the same nondimensional time. Figure 16 represents the evolution of the error committed on the capsule shape ℎ on the training database at /ℓ = 0, 0.4, 1, 2, 5, 10. At initial time, Shape is zero. The interpolation method is therefore able to capture the initial capsule shape. When the time increases, Shape increases and greater than if we apply directly the POD-DMD method on the FOM results and reconstruct the dynamics. However, Shape remains less than 0.3% on the majority of the testing database. It remains fully acceptable.
Shape is more important near the steady-state limit and when we approach the lowest values of because we are close to the limits of the training base.

Discussion and conclusion
As a summary, in this paper we have considered a -parametrized reduced-order model of microcapsule dynamics in the form The vector = ( , /ℓ) contains the governing parameters, the coefficients ( , ) and ( , ) are spectral coefficients of POD decomposition for the displacement and velocity fields respectively, and the matrix ( ) is identified from data using a dynamic mode decomposition least-square procedure. We have numerically proven for a broad range of capillary numbers and aspect ratios /ℓ that it is able to capture the dynamics up to the steady state of a capsule flowing in a channel and its large deformations. As a first approach, we have presently chosen to use a DMD method that is linear in time to build the ROM model. Still the ROM model captures spatial non-linearity by means of the POD modes. The   resulting reduced-order model is of great fidelity, weak discrepancies being only observed in the early transient stage. We have also shown that the learning time need to be larger than the transient stage duration and that we can go beyond the FOM time window used for the training of the ROM model.
For generalisation, we have computed the capsule dynamics for any parameter set. The generalization algorithm is based on interpolation: we first pre-calculate the ROM dynamic model at a finite number of points in the parameter space domain and determine the , and (and thus the capsule displacement) at these points. For any other value of the parameters, we first predict the time-evolution of the capsule node displacements using a linear interpolation procedure in the parameter space and then build a dynamical system based the DMD methodology. The error is mostly below 0.3% over the entire domain, which proves the precision and utility of the ROM approach.
Like any other data-driven model, the model requires a certain number of high-fidelity simulations to provide accurate predictions. By discretizing the parameter space in a regular and homogeneous way (Figure 12), we have not presently tried to optimize the number of FOM simulations. But sampling strategies like the Latin Hypercube Sampling (LHS) exist and result in a net reduction in FOM simulation number. The empirical law, conventional among the data-driven model community, is that one needs between 10 × and 50 × points, where is the dimension of the problem ( = 2 in our case). This law shows that the number of high-fidelity simulations does not explode with the problem dimension, owing to the linear dependence of the law.
To prove the generality of the proposed approach, we additionally apply the ROM to a capsule in simple shear flow. This classical case was extensively studied over the past years (Ramanujan & Pozrikidis 1998;Lac & Barthès-Biesel 2005;Li & Sarkar 2008;Walter et al. 2010;Foessel et al. 2011;Barthès-Biesel et al. 2010;Dupont et al. 2015). We build a ROM model to predict the evolution of an initially spherical capsule subjected to a shear rate until = 10 with 15 modes, a learning time of = 10 and = 10 −6 . The time step Δ between each snapshot is equal to 0.04. We retrieve that the initial capsule is elongated in the straining direction by the external flow and that the membrane rotates around the deformed shape due to the flow vorticity ( Figure 17). The ROM is thus able to recover the tank-treading motion. The capsule contours in both shear and perpendicular planes predicted by the ROM and simulated by the FOM are in very good agreement ( Figure 18). Figure 19 shows the evolution of the maximum error on the capsule shape for different values of . At = 0.1, folds appear periodically on the capsule, which prevents the ROM from reproducing precisely the wrinkling phenomenon. But from 0.3, the error is reduced by an order of magnitude and is below 0.2%. The linear differential model is stable as soon as the eigenvalues of have nonpositive real parts, and is consistent with steady states as soon as zero is an eigenvalue. Numerical experiments show that identified matrices from data have eigenvalues with negative real parts and one of the eigenvalues is very close to zero.
As it is often the case with spectral-like methods, there is a trade-off between accuracy and ill-conditioning effects: when a large number of POD modes are used ( > 20), the data matrix X of snapshot POD coefficients is ill-conditioned. For the determination of , we have used a Tikhonov regularization in the least square cost function (see (3.17)) in order to have a better conditioned problem and a -curve procedure to determine the best regularization coefficient . Unfortunately we observe some limitations in the accuracy. A perspective would be to use a proximal approach: within an iterative procedure, at iteration ( + 1), compute the matrix ( +1) solution of ( +1) = arg min using (0) = 0. At convergence, one can observe that the regularization term vanishes, so that one can expect better accuracy with this approach. This will be investigated in a future work.
We have proposed a successful and very efficient ROM for FSI problems. It is an alternative to the use of HPC. It must be seen as a complimentary (and non-competing) approach to full-order models, and has many advantages. Among them, one can mention the easiness in implementation. It leads to a very handy set of ODEs, that are easy to determine from an algorithmic point of view. Furthermore, the system can be run on any computer. The size of the matrices is, indeed, reduced from (3 × 2562nodes × 250snapshots) to about (3 × 2562nodes × ( + 1)), where the number of modes is = 20. The computation required time is a few milliseconds for one parameter set. The current speedups are between 5 000 and 52 000, which out-performs any full-order model approach. We believe that this work is an encouraging milestone to move toward real time simulation of general coupled problems and to deal with high-level parametric studies, sensitivity analysis, optimization and uncertainty quantification.
The next milestone following this work would be to go toward nonlinear differential dynamical systems as reduced-order models. There is three natural ways for that. The first one is to use Kernel Dynamic Model Decomposition (KDMD) rather than DMD. But we have recently shown in De Vuyst et al. (2022) that a non-linear low-order dynamical model does not provide significant improvement. The second one is to use Extended Dynamic Model Decomposition (EDMD) (Williams et al. 2015). The EDMD method adds some suitable nonlinear observables (or features) in the data, so that a linear 'augmented' dynamical system is searched for. A third option is would be to directly use artificial neural networks (ANN), in particular recurrent neural networks (RNN) (Trischler & D'Euleuterio 2016). The RNN training would replace the DMD procedure, and would be trained with the same POD coefficient matrices X and Y. As shown in the recent study by Lin et al. (2021), artificial intelligent may prove to be efficient and precise to predict capsule deformation. As soon as X is a full-rank matrix (meaning that and we reasonably have linearly independent measurements of ), the matrix XX is invertible and we get the estimate Let us denote by (resp. ) the (complex) eigenvalues of (resp. ). We have . Suppose now that we use a small time step Δ . From a Taylor expansion, we observe that = + Δ 2 ( ) 2 + (Δ ).
We would like to study what is the effect of the first-order error term Δ 2 ( ) 2 on the stability of the reconstructed dynamical system = . Suppose that the complex number has real and imaginary parts and respectively. Then = + Δ 2 ( 2 − 2 ) + (1 + Δ ) + (Δ ).
So there is again a time step Δ ★ > 0 for which, for any Δ < Δ ★ we have Re( ) 0.
As a conclusion, starting from a stable linear dynamical system (in the sense that Re( ) 0 for all ), using a small enough time step Δ and the forward Euler time discretization, the identification method leads to an estimated dynamical system which is also stable. Let us underline that this could not be the case using another time discretization as e.g. for the backward Euler time discretization and the associated least square problem We observe that for a center with a pure imaginary eigenvalue = , ≠ 0, one gets = Δ 2 2 + (Δ ) therefore > 0 for a small enough Δ . This is a counter-intuitive result: for numerical simulations, it is known that the backward Euler scheme provide more stability than the forward one. For system identification with time discretization of the residual term, it is safer to use the forward Euler scheme for stability of the estimated dynamical model.

Appendix B. Kinetic energy conservation
Another quantity of interest is the capsule kinetic energy { } 2 . Since the capsules are expected to reach a steady state after a transient stage in the Stokes pipe flow, the kinetic energy should also reach a constant value. From the differential equations, semi-orthogonality of and symmetry property of the scalar product, we successively have . Dissipation property is linked to the non-positiveness of the (real) eigenvalues of .

Appendix C. Practical computation of the pseudo-inverse matrix
The Moore-Penrose pseudo-inverse X † of a matrix X of size × , , with rank(X) = is defined by X † = X (XX ) −1 . (C 1) For an ill-conditioned matrix X, the direct computation of X † by formula (C 1) is unsuitable because the condition number of XX is the square of that of X. A more robust procedure can be derived by help of the QR factorization. There exists a semi-orthogonal matrixˆ of size × and an upper triangular square matrix of size × such that X =ˆ . Moreover, is invertible because X is assumed to be a full-rank matrix. Since X † is the solution of the matrix system X † (XX ) = X , we get X † ˆ ˆ =ˆ .