<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRMECA,Unicode,biblatex,published]{cedram}
\addbibresource{crmeca20250172.bib}
\newenvironment{noXML}{}{}%
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\sfrac#1#2{{#1}/{#2}}
\def\sofrac#1#2{({#1}/{#2})}
\def\stbfrac#1#2{({#1})/({#2})}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}\hline}
\def\botline{\\\hline}
\def\xsection{}
\allowdisplaybreaks
\def\ndash{\text{--}}
\def\rmpi{\uppi}
\usepackage{etoolbox}
\usepackage{upgreek}
\usepackage{upgreek}
\usepackage[T1]{fontenc}
\def\jobid{crmeca20250172}
\def\morerows#1#2{#2}
%\graphicspath{{/tmp/\jobid_figs/web/}} 
\graphicspath{{./figures/}} 
\def\0{\phantom{0}}
\def\mn{\phantom{$-$}}
\def\refinput#1{}
\newcounter{runlevel}
\let\MakeYrStrItalic\relax
\csdef{Seqnsplit}{\\}
\def\Tg#1#2#3#4{}
\def\refinput#1{}
\def\back#1{}
\DOI{10.5802/crmeca.304}
\datereceived{2025-02-24}
\daterevised{2025-05-20}
%\datererevised{2025-05-22}
\dateaccepted{2025-05-21}
\ItHasTeXPublished
}
{\documentclass[crmeca]{article}
\usepackage{upgreek}
\usepackage[T1]{fontenc} 
\usepackage{stmaryrd}
\def\CDRdoi{10.5802/crmeca.304}
\def\newline{\break}
\def\hyperlink#1#2{#2}
\def\hypertarget#1#2{#2}
\def\sfrac#1#2{{#1}/{#2}}
\def\sofrac#1#2{({#1}/{#2})}
\def\stbfrac#1#2{({#1})/({#2})}
\def\CDRsupplementaryTwotypes#1#2{}
}
\makeatother

\dateposted{2025-06-25}
\begin{document}

\begin{noXML}

\CDRsetmeta{articletype}{review}

\title{Multiphysics simulations in biomedical: fluid--structure
interaction coupled with mass transfer}

\alttitle{Simulations multiphysiques en biom\'edical : interaction
fluide-structure coupl\'ee au transfert de masse}

\author{\firstname{Badr} \lastname{Kaoui}\CDRorcid{0000-0003-2478-3858}}
\address{Equipe Interactions Fluides Structures Biologiques (IFSB),
Laboratoire Biom\'ecanique et Bioing\'enierie (BMBI),
UMR 7338 - Centre National de la Recherche Scientifique (CNRS) \&amp;
Universit\'{e} de Technologie de Compi\`{e}gne (UTC),
60200 Compi\`{e}gne, France}
\email{badr.kaoui@utc.fr}

\keywords{\kwd{CFD}
\kwd{FSI} 
\kwd{Multiphysics simulations}
\kwd{HPC}
\kwd{Biofluids}
\kwd{Biomedical engineering}}

\altkeywords{\kwd{M\'ecanique des fluides num\'erique}
\kwd{IFS}
\kwd{Simulations multiphysiques}
\kwd{Calcul haute performance}
\kwd{Biofluides}
\kwd{G\'enie biom\'edicale}}

\thanks{National Research Agency for funding the MechLymph
(ANR-20-CE45-0008-01),  UC-DC (ANR-23-CE51-0057), SOCA
(ANR-24-CE06-6245-02), GENCI-IDRIS (Grant 2024-[AD012A13687R2])}

\begin{abstract} 
We propose numerical schemes to study complex multiphysics problems
involving fluid--structure interaction (FSI) coupled with mass/heat
transfer (HMT). We primarily use the Lattice Boltzmann Method (LBM) to
compute fluid flow and the advection--diffusion--reaction of chemical
entities. The dynamics and deformation of structures are computed using
the Lattice Spring Method (LSM), and we implement the Immersed Boundary
Method (IBM) to achieve two-way coupling of fluid--structure
interaction. For biomedical applications, we study controlled drug
delivery from particles subjected to flow, passage of soft particles
through microfluidic constrictions, performance of artificial
pancreas-on-chip, as well as contraction--relaxation of lymphatic
vessels and the opening--closing of their valves.
\end{abstract}

\begin{altabstract} 
On propose des sch\'emas num\'eriques pour \'etudier des probl\`emes
multiphysiques complexes impliquant l'interaction fluide-structure
(IFS)  coupl\'ee au transfert de masse et de chaleur (TMC). On utilise
principalement la m\'ethode de Boltzmann sur r\'eseau (LBM) pour
simuler l'\'ecoulement des fluides ainsi que 
l'advection-diffusion-r\'eaction d'entit\'es chimiques.
La dynamique et la d\'eformation des structures sont calcul\'ees \`a
l'aide de la m\'ethode des ressorts sur r\'eseau (LSM), et on met en
\oe uvre la m\'ethode de la fronti\`ere immerg\'ee (IBM) pour assurer
le couplage bidirectionnel de l'interaction fluide-structure.
Dans le cadre des applications biom\'edicales, on \'etudie la
lib\'eration contr\^ol\'ee de m\'edicaments \`a partir de particules
soumises \`a un \'ecoulement, le passage de particules d\'eformables
\`a travers des constrictions microfluidiques, la performance de
pancr\'eas artificiels sur puce, ainsi que la contraction-relaxation
des vaisseaux lymphatiques et l'ouverture-fermeture de leurs valves.
\end{altabstract}

\maketitle

\end{noXML}

\section{Introduction}\label{sec1}
In our modern world, scientists and engineers are increasingly faced
with solving complex problems that involve multiple phenomena occurring
at multiple scales and simultaneously. In biomedical engineering, the
use of modeling and computer simulations has gained growing trust among
researchers. These simulations are increasingly employed to assist in
experiment planning and to elucidate the relationships and interactions
between different elements and factors in biological, biomedical and
biomimetic systems. While powerful commercial software solutions are
available on the market, they are primarily designed to meet the needs
of the aerospace and automotive industries. This highlights the
necessity of developing in-house numerical codes tailored to the
specific requirements of biomedical engineering research, and enabling
high-performance multiphysics simulations. The following sections
present the numerical techniques we use for biomedical applications,
which can also be adapted and applied to other scientific domains and
industrial sectors. Most biological and biomedical systems involve
transport phenomena, including fluid flow, biochemical species
transport, and interactions with deformable structures (see 
Figure~\ref{fig:schema}). Additionally, these systems operate under the
influence of external stimuli, such as applied flow or variations in
ambient temperature.

\begin{figure}[h!]
\includegraphics{fig01}
\caption{\label{fig:schema}Coupling of the numerical methods used for
solving each part and module of a multiphysics problem.}
\end{figure}

\section{Transport phenomena}\label{sec2}
The flow of an incompressible Newtonian Fluid is governed by the
Navier--Stokes equations,
{\begin{eqnarray}
&amp;\displaystyle
\frac{\partial \mathbf{u}}{\partial t} + 
\mathbf{u} \boldsymbol{\cdot} \nabla \mathbf{u} = 
-\frac{1}{\rho}\nabla p + \nu \nabla ^2 \mathbf{u} + \mathbf{f},
\label{eq:NS1}\Seqnsplit
&amp;\displaystyle
\nabla \boldsymbol{\cdot} \mathbf{u} = 0,
\label{eq:NS2}
\end{eqnarray}}\unskip
and the transport of chemical entities by the
advection--diffusion--reaction equation,
{\begin{equation}
\frac{\partial c}{\partial t} + 
\mathbf{u} \boldsymbol{\cdot} \nabla c = 
\nabla \cdot (D \nabla c) + S,
\label{eq:conv_diff}
\end{equation}}\unskip
where $\rho$ and $\nu$ are respectively the density and the kinematic
viscosity of the studied fluid (the solvent). $D$ is the diffusion
coefficient of the solute. The unkown solutions of these equations are
the velocity field $\mathbf{u}(\mathbf{x},t)$, the pressure field
$p(\mathbf{x},t)$ and the concentration field $c(\mathbf{x},t)$, where
$\mathbf{x}$ is the position and $t$ the time. $\mathbf{f}$ represents
a body force acting on the fluid, and $S$ is a source term.

Instead of solving these partial differential equations directly, for
example, using the Finite Element Method (FEM), we use an alternative
approach: the Lattice Boltzmann Method (LBM) to compute the 
solutions~\cite{Succi2001,Sukop2006,Aidun2010}. This method is well
suited for predicting the dynamics of complex fluids (e.g., particle
suspensions, multiphase fluids) and computing flows in complex
geometries (e.g., porous media, microfluidic systems).

Initially developed by physicists to overcome the drawbacks of the
Lattice Cellular Automata method~\cite{Sukop2006}, such as the lack of
Galilean invariance, the Lattice Boltzmann Method is now widely used by
scientists and engineers. Here, we will not present the theoretical
foundations of the method, but rather explain its principles and how it
works.

The distinctive feature of LBM lies in the use of a distribution
function that represents the probability of finding a fluid packet
(pseudo-particles) at a position $\mathbf{x}$, at time $t$ and with a
discrete velocity direction  $\mathbf{e}_i$. In LBM, both position and
velocity (momentum) spaces are discretized, with $\Delta x = \Delta y =
\Delta z = 1$ and $\Delta t = 1$. The evolution over time of the
distribution $f_i$ is given by the so-called Lattice Boltzmann
equation,
{\begin{equation}
f_i(\mathbf{x} + \mathbf{e}_i, t + 1) - 
f_i(\mathbf{x},t) = \Omega_i (\mathbf{x},t),
\label{eq:lbea}
\end{equation}}\unskip
with the left-hand side describing the streaming step of $f_i$, and the
right-hand side representing the collision step through the operator
$\Omega_i$. Two types of operators are widely used: BGK and MRT. The
BGK (Bhatnagar--Gross--Krook) operator is given by,
{\begin{equation}
\Omega_i (\mathbf{x},t)= 
-\frac{f_i - f_i^{\mathrm{eq}}}{\tau _{f}},
\end{equation}}\unskip
where $\tau_f$ is the relaxation time required for $f_i$ to tend toward its equilibrium function,
{\begin{equation}
f^{\mathrm{eq}} _{i}(\mathbf{x},t)=
\omega _i \rho [1 + 3(\mathbf{u}\boldsymbol{\cdot} \mathbf{e}_i) + 
\tfrac{9}{2}(\mathbf{u}\boldsymbol{\cdot}\mathbf{e}_i)^2 -
\tfrac{3}{2}(\mathbf{u})^2 ],
\end{equation}}\unskip
where $\omega_i$ are the weights associated to the type of the used
lattice, see  Figure~\ref{fig:lattices}. The microscopic physical
quantity $\tau_f$ is related to the macroscopic physical quantity, the
kinematic viscosity of the fluid $\nu$, through the formula,
{\begin{equation}
\nu = \tfrac{1}{3}
(\tau _{f} - \tfrac{1}{2}).
\end{equation}}\unskip
Thus, $\tau_f$ can only take values greater than $1/2$, and when
$\tau_f$ approaches $1/2$, numerical instabilities arise. The remedy
for this issue is the use of the Multiple-Relaxation-Time (MRT)
operator, which introduces different relaxation times for each velocity
direction $\mathbf{e}_i$.

\begin{figure}
\includegraphics{fig02}
\caption{\label{fig:lattices}D2Q9 lattice (two-dimensional space, 2D,
with 9 velocity directions, Q9) and D3Q19 lattice (three-dimensional
space, 3D, with 19 velocity directions, Q19).}
\end{figure}

The non-slip boundary condition on the solid walls of the domain,
$\mathbf{u} = 0$, is set using the bounce-back step in the LBM: the
distribution function is reflected in the opposite direction
($\mathbf{e}_{-i}$) to the incident velocity direction ($\mathbf{e}_i$)
if the position $\mathbf{x}$ is a solid node,
{\begin{equation}
f_{-i}(\mathbf{x} + \mathbf{e}_{-i}, t + 1) = 
f_i(\mathbf{x},t).
\label{eq:lbe}
\end{equation}}\unskip
There are also boundary conditions to be set at the inlet and the
outlet of the domain. These can be periodic or of the Zou and He
type~\cite{Zou1997}, which allow for setting either velocities or
pressures at the boundaries of the computational domain.

The local pressure and velocity of the fluid are evaluated respectively
by the equation of state of the fluid and the first moment of the
distribution $f_i$,
{\begin{eqnarray}
p(\mathbf{x},t) 
&amp;= &amp;
\frac{1}{3}\sum _{i} f_i (\mathbf{x},t) = 
\frac{1}{3}\rho(\mathbf{x},t),
\Seqnsplit
\mathbf{u}(\mathbf{x},t)
&amp;=&amp;
\frac{1}{\rho(\mathbf{x},t)}
\sum _{i} f_i (\mathbf{x},t)\mathbf{e}_i.
\end{eqnarray}}\unskip

The LBM is also used to compute numerically the solution of the
advection--diffusion--reaction 
equation~\cite{Kaoui2017,Kaoui2018,Kaoui2020}. To achieve this, an
additional Lattice Boltzmann equation is introduced to govern the
dynamics of another distribution function, $g_i$, associated with mass
transfer,
{\begin{equation}
g_i(\mathbf{x} + \mathbf{e}_i, t + 1) - 
g_i(\mathbf{x},t) =
-\frac{g_i - g_i^{\mathrm{eq}}}{\tau _{g}},
\label{eq:lbe2}
\end{equation}}\unskip
whose relaxation time $\tau _{g}$ is related to the solute diffusion
coefficient $D$ through the formula,
{\begin{equation}
D = \tfrac{1}{3}
(\tau_ {g} - \tfrac{1}{2}),
\end{equation}}\unskip
and the concentration is given by the zero-th moment of $g_i$,
{\begin{equation}
c(\mathbf{x},t) = 
\frac{1}{3}\sum _{i} g_i (\mathbf{x},t).
\end{equation}}\unskip

In 2D, the D2Q9 lattice is used with 9 velocity directions ($i=0\ndash8$),
and in 3D, the D3Q19 lattice is considered with 19 velocity directions
($i=0\ndash18$), see  Figure~\ref{fig:lattices}. The weights of these two
lattices are,
{\begin{eqnarray}
&amp;\displaystyle
\mathrm{2D}{:} \quad 
\omega _i=
\begin{cases}
4/9,  &amp; \mbox{for}\;i=0,\\
1/9,  &amp; \mbox{for}\;i=1\ndash4,\\
1/36, &amp; \mbox{for}\;i=5\ndash8.
\end{cases}
\Seqnsplit
&amp;\displaystyle
\mathrm{3D}{:} \quad
\omega _i=
\begin{cases}
1/3,  &amp; \mbox{for}\;i=0,\\
1/18, &amp; \mbox{for}\;i=1\ndash6,\\
1/36, &amp; \mbox{for}\;i=7\ndash18.
\end{cases}
\end{eqnarray}}\unskip
The numerical stability of the LBM solver for the mass transfer part is
further enhanced by using the same lattice as that of the LBM flow 
solver~\cite{Kaoui2017}. Moreover, this strategy enables the
development of modular and ergonomic computer codes.

The body force $\mathbf{f}$ in the Navier--Stokes equations, the source
term $S$ in the advection--diffusion--reaction equation, and the boundary
conditions vary depending on the studied problem. In LBM, $\mathbf{f}$
is accounted for by adding the term $3\omega_i (\mathbf{f} \boldsymbol{\cdot}
\mathbf{e}_i)$ to the right-hand side of 
Equation~(\ref{eq:lbe})~\cite{Kaoui2016}, while $S$ is incorporated by
adding the term $\omega_i S$ to the right-hand side  of
Equation~(\ref{eq:lbe2})~\cite{Kaoui2020}.

\section{Deformable structures}\label{sec3}
The mechanics of deformable structures are governed by different
elasticity laws. We are particularly interested in structures whose
thickness is negligible compared to their size. We model them as
two-dimensional structures of zero thickness. These structures evolve
in three-dimensional space and are subjected to three deformation modes
under the action of applied stresses: compression--dilation, shear and
bending. For each mode, the structure exerts in turn a restoring force
to recover its equilibrium shape at rest.

For small deformation regimes (e.g., Ref.~\cite{Bielinski2022}), we use
the Lattice Spring Method (LSM)~\cite{Ostoja-Starzewski2002}, which
models the structure as a network of points interconnected by springs,
see  Figure~\ref{fig:ressorts}. The mechanics of this network, at the
microscopic scale, is capable of reproducing faithfully the mechanical
behavior of the structure at the macroscopic scale. Each point $i$ in
the network experiences a force from its neighboring points $j$,
governed by Hooke's law,
{\begin{equation}
\mathbf{F}(i) = 
-\kappa \sum _{j \neq i}
[ d(i,j) - d_0(i,j) ] 
\frac{\mathbf{r}(i,j)}{\| \mathbf{r}(i,j) \|},
\label{eq:hook}
\end{equation}}\unskip
where $\kappa$ is the spring rigidity, $d$ and $d_0$ are the current
and initial distances between points $i$ and $j$, respectively, and
$\mathbf{r}$ the position vector connecting these points. Other types
of spring networks can be used, such as the one generated by triangular
meshing of the domain. The constant $\kappa$ can be related and
estimated from the macroscopic elastic moduli of a material, such as
Young's modulus, depending on the lattice  type~\cite{Omori2011}.

\begin{figure}
\includegraphics{fig03}
\caption{\label{fig:ressorts}Two-dimensional network of points
interconnected by Hookean springs of 
\mbox{rigidity} $\kappa$.}
\end{figure}

For problems involving large deformation regimes (e.g.,
Refs.~\cite{Kusters2014,Krueger2014,Bielinski2021}), models derived
from continuum mechanics are more suitable, as they allow for the
consideration of the nonlinear elasticity of structures. For instance,
the Neo-Hookean law, which describes the energetic cost of the shear
deformation mode,
{\begin{equation}
E_{\mathrm{s}} = 
\frac{\kappa _{\mathrm{s}}}{2} 
\int _{A} 
\left[I_1 - 1 + \frac{1}{I_2 + 1}\right] \mathrm{d}A,
\end{equation}}\unskip
where $\kappa_{\mathrm{s}}$ is the shear modulus, $I_1 = \lambda_1^2 +
\lambda_2^2 - 2$ and $I_2 = \lambda_1^2 \lambda_2^2 - 1$ are the two
strain invariants, $\lambda_1$ and $\lambda_2$ are the principal
stretch ratios. $A$ and $\mathrm{d}A$ are the surface area and surface
element of the structure, respectively. We evaluate this energy using
FEM.

\section{Fluid--structure interaction}\label{sec4}
Coupling between structure dynamics and fluid flow, known as
fluid--structure interaction (FSI), is achieved using the Immersed
Boundary Method  (IBM)~\cite{Peskin2002,Verzicco2023,Mittal2023}. In
this method, the structure is represented by a mobile Lagrangian mesh,
immersed in a stationary Cartesian Eulerian mesh, where the flow is
computed. See the illustration in Figure~\ref{fig:ibm} for the 2D case.

\begin{figure}
\includegraphics{fig04}
\caption{\label{fig:ibm}(A)~Deformable structure, here a vessel with a
bi-leaflet valve. (B)~Discretization of the structure into points
interconnected by springs (moving Lagrangian mesh), immersed in a fixed
grid on which the flow is computed (fixed Eulerian mesh). (C)~Nodes
used by the Immersed Boundary Method for interpolating velocities onto
the structure and extrapolating forces to the fluid nodes.}
\end{figure}

Physical quantities calculated at a point $\mathbf{r}_1 \equiv (x_1,
y_1, z_1)$ of one mesh are communicated to a point $\mathbf{r}_2 \equiv
(x_2, y_2, z_2)$ of the other mesh by interpolation or extrapolation
using a Dirac-type function $\Delta$, see  Figure~\ref{fig:dirac}. In
3D, we use,
{\begin{equation}
\Delta(\mathbf{r}_1,\mathbf{r}_2) = 
\delta(x_1,x_2)\delta(y_1,y_2)\delta(z_1,z_2),
\end{equation}}\unskip
with,
{\begin{equation}
\delta (a,b) =
\begin{cases} 
\displaystyle
\frac{1}{4}
\left(1+ \cos \frac{\rmpi (a - b)}{2}\right) &amp;
\mbox{if}\;|a - b | \leq 2, \\
0 &amp; \mbox{otherwise}.
\end{cases}
\end{equation}}\unskip
The velocity $\mathbf{v}$ of a point $\mathbf{r}$ on the structure is
thus evaluated by interpolating the velocities $\mathbf{u}$ computed at
the surrounding fluid points $\mathbf{x}$ (see 
Figure~\ref{fig:ibm}C---advection),
{\begin{equation}
\mathbf{v}(\mathbf{r},t)  = 
\sum_{\mathbf{x}} \Delta(\mathbf{r},\mathbf{x})\mathbf{u}(\mathbf{x},t).
\end{equation}}\unskip
This is the velocity used to move each point of the Lagrangian mesh,
for example, using the Euler scheme, under the assumption that inertial
effects are negligible,
{\begin{equation}
\mathbf{r}(t+1) = \mathbf{r}(t) + \mathbf{v}(\mathbf{r},t).
\end{equation}}\unskip
The displacements of all the points cause the structure to evolve as a
whole towards an out-of-equilibrium shape. As a result, the structure
exerts a reaction force on the fluid, thereby perturbing and altering
the flow in its vicinity. The force $\mathbf{F}$ exerted by a point
$\mathbf{r}$ of the structure is distributed by extrapolation to the
nearby fluid points $\mathbf{x}$ (see 
Figure~\ref{fig:ibm}C---reaction),
{\begin{equation}
\mathbf{f}(\mathbf{x},t) = 
\sum_{\mathbf{r}}
\Delta(\mathbf{r},\mathbf{x})\mathbf{F}(\mathbf{r},t).
\label{eq:force}
\end{equation}}\unskip

\begin{figure}
\includegraphics{fig05} 
\caption{\label{fig:dirac}The $\Delta$ function used by IBM in 2D.}
\end{figure}

\section{Multiphysics coupling}\label{sec5}
The three numerical methods, LBM, LSM (or FEM) and IBM, are coupled and
executed following an explicit time-stepping scheme. Each iteration in
the main code is structured into three distinct blocks: the LBM block
for the transport phenomena part (CFD and/or HMT), the LSM block for
the structure mechanics part (Structure), and the IBM block for the
two-way FSI coupling, organized into two levels: advection and
reaction.

The LBM and LSM blocks operate independently without direct
interaction.  It is the IBM block that ensures the coupling between the
LBM and LSM.  This two-way coupling enables the flow to advect and
deform the structure (fluid--structure interaction), while also
accounting for the flow disturbance caused by the presence of the
structure (structure--fluid interaction).  The various execution steps
of a typical computer code of multiphysics coupling are outlined in the
flowchart in  Figure~\ref{fig:algo}. The main challenge in this
multiphysics coupling explicit scheme lies in identifying suitable
numerical values for the key parameter set, particularly the coupling
coefficients, to ensure numerically stable simulations.  This involves
trial-and-error tuning of physical and numerical parameters, estimation
of key dimensionless numbers relevant to the problem, and analysis of
the relative importance of different physical effects.

\begin{figure}
\includegraphics{fig06}
\caption{\label{fig:algo}Flowchart summarizing the algorithm describing
the sequence of the execution steps for each block of the code: LBM,
LSM and IBM, as well as the coupling levels between them.}
\end{figure}

\section{Numerical library}\label{sec6}
We have developed a numerical library containing various modules that,
when coupled, enable multiphysics simulations (see 
Figures~\ref{fig:schema},~\ref{fig:algo} and 
Refs.~\cite{Kaoui2009,Kaoui2011,Kaoui2016,Kaoui2017,Kaoui2018,Kaoui2020,Kaoui2022,Bielinski2022,Bielinski2023}).
This library, unlike existing commercial and open-source software, is
better suited for biomedical applications involving transport phenomena
at the mesoscopic scale. It allows for coupling between the dynamics of
complex fluids and deformable structures, as well as mass or heat
transfer and the consideration of external stimuli.

For the simulation of fluid flow and mass or heat transfer, we have
developed the \textit{lbm2d} (2D simulations:
\url{https://gitlab.utc.fr/kaouibad/lbm2d}) and \textit{lbm3d} (3D
simulations: \url{https://gitlab.utc.fr/kaouibad/lbm3d}) modules based
on the LBM method. For the fluid--structure interaction, we have
developed the \textit{ibm2d}
(\url{https://gitlab.utc.fr/kaouibad/ibm2d}) and \textit{ibm3d}
(\url{https://gitlab.utc.fr/kaouibad/ibm3d}) modules based on the IBM
method. For computing the structure deformation, we have developed
modules based mainly on the LSM method. For the post-processing, we
have developed the \textit{io3d} module
(\url{https://gitlab.utc.fr/kaouibad/io3d}), which allows reading and
writing output files in VTK format (\url{https://vtk.org}). These files
are readable by the open-source software ParaView
(\url{https://www.paraview.org}) that allows 3D visualization of the
results (e.g., structure, streamlines, or concentration fields). The
functions in the library are programmed in modern Fortran language. The
2D codes are parallelized using CAF (Co-Array Fortran), which is better
suited for work stations with shared memory, while the 3D codes use MPI
(Message Passing Interface), a more appropriate parallelization method
for supercomputers with distributed memory. 

Figure~\ref{fig:speed} shows the performance of the \textit{lbm3d}
module on the Jean Zay supercomputer of the IDRIS institute of the
CNRS. The case study involves simple shear flow between two walls
(\url{https://gitlab.utc.fr/kaouibad/examples}),
computed in a 3D domain of size $400 \times 42 \times 42$ (a total of
705,600 points), with a total of 2000 time iterations. The performance
of the \textit{lbm3d} solver on the Jean Zay supercomputer demonstrates
strong scalability across a range of CPU counts.  The measured speed-up
closely follows the ideal linear trend, indicating efficient
parallelization and effective utilization of computational resources. 
While the solver maintains near-linear speed-up at lower to moderate
CPU counts, a slight deviation from ideal performance is observed as
the number of CPUs increases, which is typical due to factors such as
communication overhead.  Overall, the \textit{lbm3d} solver exhibits
high parallel efficiency, making it well-suited for large-scale
simulations on modern high-performance computing architectures.

\begin{figure}
\includegraphics{fig07}
\vspace*{-4pt}
\caption{\label{fig:speed}Performance of the \textit{lbm3d} solver on
Jean Zay supercomputer of IDRIS institute of CNRS. The straight line
indicates linear speed-up: perfect scaling.}
\vspace*{-4pt}
\end{figure}

\section{Applications}\label{sec7}
We limit ourselves to show only three applications from our ongoing
research projects:

\medskip\noindent
(i)~3D simulation of drug release from a capsule subjected to
external shear flow (cf.\  Figure~\ref{fig:caps}).  This simulation
allows for predicting the spatial distribution of the released
molecules, which is altered by the capsule dynamics and the resulting
flow.  This pharmaceutical application involves coupling between the
flow of the internal and the external fluids of the capsule, the
advection--diffusion of the solute, and the fluid--structure--solute 
coupling~\cite{Bielinski2022,Bielinski2023}.  The fluid transports the
solute by advection, while the solute itself has no direct influence on
the flow.

\begin{figure}
{\vspace*{-2pt}}
\includegraphics{fig08}
{\vspace*{-2pt}}
\caption{\label{fig:caps}Drug release from a capsule subjected to shear
flow. The simulation allows for analyzing the spatial distribution of
the released molecules and the flow developed around the capsule.}
{\vspace*{-2pt}}
\end{figure}

\medskip\noindent
(ii)~3D simulation of the performance of an organ (pancreas) on a
microfluidic chip (Organ-On-Chip), cf.\  Figure~\ref{fig:ooc}.  A
spheroid composed of pancreatic cell aggregates, modeled by a suspended
sphere, is trapped and maintained in one of the wells of a microfluidic
chip.  This spheroid secretes insulin, which is transported by the
ambient flow within the device. This bioengineering problem requires
computing fluid flow in a complex geometry: a microfluidic system with
wells designed to trap the  spheroids~\cite{Essaouiba2020}.  It relies
also on simulating fluid interaction with the sphere, and two
advection--diffusion--reaction equations: one for insulin dynamics and
the other for glucose.  These last two equations are coupled by source
terms based on a feedback mechanism: the glucose level carried by the
ambient fluid stimulates and controls the degree of insulin secretion
by the spheroid.

\begin{figure}
\includegraphics{fig09}
\caption{\label{fig:ooc}Suspended spheroid trapped by one of the wells
of a microfluidic chip operating as an artificial pancreas. The
simulation predicts the dynamics of the spheroid under flow within the
well, as well as the spatial distribution of the secreted insulin.} 
\end{figure}

\medskip\noindent
(iii)~3D simulation of the dynamics of a valve composed of two
crescent-shaped leaflets (cf.\  
\mbox{Figure~\ref{fig:valve}}).  The valves
described here are typical of the vessels in the lymphatic system. 
This biological system relies on strong coupling between several
phenomena: the deformation of structures (vessels and lymphatic
valves), fluid--structure interactions, and the transport of 
\mbox{biochemical}
entities such as calcium ions  (Ca$^{2+}$) and nitric oxide (NO). 
Unlike other fluid mechanics applications, the flow in this context is
not imposed.  It is instead generated by a complex mechanism involving
the synchronization between the opening and closing of the lymphatic
valves and the contraction and relaxation of the walls of the lymphatic
vessels, which are themselves regulated by biochemical signals.  This
mechanism ensures the drainage of lymph, the fluid circulating in the
lymphatic system~\cite{Padera2016}.  The simulation in
Figure~\ref{fig:valve} shows the opening and closing of a lymphatic
valve, as well as the prediction of the distribution of wall shear
stresses on the surface of one of the leaflets.

\begin{figure}
\includegraphics{fig10}
\caption{\label{fig:valve}Opening (on the left) and closing (on the
right) of a bicuspid lymphatic valve under the effect of oscillatory
flow. The simulation also allows for analyzing the spatial distribution
of wall shear stress on the surface of one of the valve's leaflets.}
\end{figure}

The source codes used to produce the results for the three applications
can be provided upon request by contacting the author via email.

{\vspace*{-2pt}}

\section{Conclusions}\label{sec8}

{\vspace*{-2pt}}

We have presented the development and application of modern numerical
methods to study complex fluid--structure interaction problems coupled
with transport phenomena, with a particular focus on biomedical
applications.  The combination of LBM, LSM and IBM methods allows for
an efficient approach to model multiphysics problems at the
mesoscopic-scale. The applications, such as drug release from
particles, organ-on-chip and lymphatic valves, demonstrate the
versatility and potential of these tools in addressing the challenges
posed by modern biomedical applications.  Moreover, the developed
numerical library offers a flexible, modular and high-performance
alternative to commercial software. These works open promising
perspectives for the study of complex biological, biomedical and
biomimetic systems, and could also find applications in other
scientific and industrial fields such as materials science, chemical
engineering, and environmental engineering.  In the future, the
integration of other numerical methods, along with the use of machine
learning approaches, could further complement and expand the range of
applications.

\section*{Declaration of interest}
The author does not work for, advise, own shares in, or receive funds
from any organization that could benefit from this article, and has
declared no affiliations other than his research institution.

\section*{Acknowledgments}
The author thanks his students, colleagues, collaborators, as well as
the National Research Agency for funding the MechLymph
(ANR-20-CE45-0008-01), UC-DC (ANR-23-CE51-0057) and SOCA
(ANR-24-CE06-6245-02) projects.  High-performance simulations were
carried out on the PILCAM2 supercomputers
(\url{https://pilcam2.utc.fr}) at UTC and Jean-Zay at CNRS, thanks to
the allocation of computing time provided by GENCI-IDRIS (Grant
2024-[AD012A13687R2]).

\CDRGrant[ANR]{ANR-20-CE45-0008-01}
\CDRGrant[ANR]{ANR-23-CE51-0057}
\CDRGrant[ANR]{ANR-24-CE06-6245-02}
\CDRGrant[GENCI-IDRIS]{2024-[AD012A13687R2]}

\back{}

\printbibliography
\refinput{crmeca20250172-reference.tex} 

\end{document}
</pre></body></html>