# Modelling and numerics for respiratory aerosols

Laurent Boudin, Céline Grandmont, Alexander Lorz, Ayman Moussa

► To cite this version:

Laurent Boudin, Céline Grandmont, Alexander Lorz, Ayman Moussa. Modelling and numerics for respiratory aerosols. Communications in Computational Physics, 2015, 18 (3), pp.723-756. 10.4208/cicp.180714.200415a. hal-01044590

HAL Id: hal-01044590

https://hal.science/hal-01044590

Submitted on 29 Aug 2014

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
MODELLING AND NUMERICS FOR RESPIRATORY AEROSOLS

LAURENT BOUDIN, CÉLINE GRANDMONT, ALEXANDER LORZ, AND AYMAN MOUSSA

ABSTRACT. In this work, we first deal with the modelling and the discretization of an aerosol evolving in the air, in the respiration framework, within a domain which can be fixed or moving. We also investigate basic numerical properties of the numerical code which was developped, and also focus on the influence of the aerosol on the airflow.

# 1. INTRODUCTION

The evolution of droplets or particles in a surrounding fluid is a phenomenon encountered in several areas, ranging from medicine (aerosol therapy) to motor industry (transport and combustion of petrol). In the particular aerosol therapy case, in vivo observations of drug delivery in the airways induce several difficulties. For instance, aerosol deposition maps require heavy experimental protocols, which cannot be easily repeated, and the obtained measurements may not be accurate enough.

Consequently, the choice of physically relevant models, and thereafter the design of stable, efficient numerical methods allow *in silico* experiments which can provide a wide range of results for various physical situations and parameters (type of aerosols, surrounding fluids, pathological state...).

Several kinds of modelling are available to describe the aerosol movement in a fluid. Two of them are discussed quite in detail in [31]: one can consider individual particles on the one hand, or a collection of particles on the other hand.

Two-phase models consider the collection of droplets or particles as a fluid and study the evolution, for instance, of aerosol concentration in the ambient fluid. Those models are most certainly adapted in the case when the volume fraction occupied by the dispersed phase is not negligible with respect to the volume fraction occupied by the surrounding fluid [34, 7, 36, 17]. Unfortunately, such models do not allow an accurate description of particle deposition. This is why we shall only focus on spray models in which the discrete aspect of the dispersed phase is kept.

Following particles as individuals is the other classical strategy, see [9, 15, 16, 39, 46] for instance. Nevertheless, describing the behaviour of such a (very) large number of particles may lead to both technical and numerical difficulties if one tries to keep the track of each individual trajectory. For instance, the *Atomiser pocket aeroneb GO* from DTF<fnref n="1" /> Corporation has the following characteristics: airflow rate of 0.3 mL/min, average (in mass) diameter equal to 3.6 µm. Hence, this nebulizer allows the injection of 10¹⁰ particles in one minute.

In this context of very numerous particles, and since the volume occupied by the aerosol remains negligible in the human airways, the formalism of statistical physics and kinetic theory is especially well-fitted. This type of coupling was first introduced by O'Rourke [40] or Williams [43] and is now quite often used to model aerosol transport in the lung, see [10, 27, 13]. As for the interaction between the aerosol and the surrounding fluid, following a nomenclature introduced by O'Rourke (see also [22]), we assume the spray to be *thin*. This means that

* the aerosol volume fraction in the mixture remains negligible;

<footnote>Date: July 17, 2014.</footnote>

<footnote>This work was partially funded by the ANR-08-JCJC-013-01 project headed by C. Grandmont and the ANR-10-BLAN-1119 project headed by M. Filoche.</footnote>

<footnote n="1">¹See http://www.dtf.fr</footnote>
* there are no interactions between the aerosol particles;

* the aerosol can have an effect on the fluid, as a response to the drag force exerted by the fluid on the particles.

Note that, for two-phase models as well as ODE ones, the aerosol retroaction on the fluid is seldom taken into account, up to our knowledge. In the same context of aerosol kinetic modelling we consider hereafter, [27] presents a study of aerosol transport in the trachea where the retroaction is taken into account.

The aerosol is then described by a distribution function which satisfies a Vlasov-type equation. The fluid is assumed to be homogeneous, Newtonian, incompressible, and can be described using the Navier-Stokes equations, see [28] for instance. The physical domain can be either fixed or time dependent.

The aerosol and fluid are coupled through two terms: the particle acceleration, depending on the relative velocity of the particle in the fluid, and the retroaction force applied by the aerosol on the fluid. Consequently, we have to deal with a strong coupling of two types of equations: one at the macroscopic scale (Navier-Stokes), one at the mesoscopic scale (Vlasov). The existence of solutions of the obtained system in a fixed domain has been investigated by various authors [30, 6, 11, 45]. The main issue in these mathematical studies lies in the fact that the system is nonlinear and *strongly* coupled.

The numerical strategy faces the same difficulty together with the different level of descriptions for both phases (macro/mesoscopic). We choose an explicit time-advancing scheme, which allows to solve the fluid and aerosol parts in a staggered way. The system is thus uncoupled in the approximation procedure, reminiscent from the existence proof in [11]. For the space discretization, a finite element procedure and a particle-in-cell method are respectively used to approximate the fluid velocity and pressure, and the distribution function. The moving domain case is handled thanks to the arbitrary Lagrangian-Eulerian (ALE) method, see [38] in a finite element context.

This whole approximation scheme was implemented in the C++ library LifeV.

In this article, we aim to investigate the range of parameters for which our scheme is numerically stable and accurate. In particular, we study the influence of the retroaction force. Its explicit treatment may induce unphysical instabilities due to large particle velocities. In the human lung, if therapeutic aerosols seem not to require the retroaction term in most standard situations, it is not the case for polluting particles, whose average volume is larger. We also consider the moving domain case, which is usually not taken into account, and is a first step towards the bronchial wall motion.

Once the full aerosol-fluid model has been described, we present the considered numerical method and then focus on three cases: a fixed domain without or with retroaction and a moving domain without retroaction. In each case, we study the numerical sensitivity with respect to various parameters (time step, mesh size, initial datum, particle representativity...).

## 2. MODEL

In the upper airways, we can safely assume that the air is Newtonian and incompressible, thus governed by the incompressible Navier-Stokes equations. During the respiration process, some airway walls may be time-dependent. Thus our equations will be considered in a moving domain. When we focus on aerosols in human airways, the number of particles can be significant, whereas the volume occupied by the aerosol remains small. A classical strategy in statistical mechanics then consists in describing the spray behavior thanks to one single kinetic equation.

To model our problem, we are led to couple both types of equations to obtain a fluid/kinetic system we present in the next subsections.

2.1. **Geometries.** In our study, a typical fluid domain, denoted by $\Omega_t$, is a cylinder or a branch (see Figure 1), which can depend on time $t$. Its boundary $\Gamma_t = \partial\Omega_t$ is divided
into three (not necessarily connected) subsets: the inlet $\Gamma_t^{in}$, the outlet $\Gamma_t^{out}$ and the wall $\Gamma_t^{wall}$. In most situations, $\Gamma_t^{in}$ and $\Gamma_t^{out}$ do not depend on time, because they are chosen as artificial boundaries. On the contrary, the motion of $\Gamma_t^{wall}$ is driven by physiological phenomena.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Branch geometry with a moving wall.</figcaption></figure>

In the sequel $n_t$ will denote the unit vector going out and normal to $\partial\Omega_t$.

2.2. **Fluid equations.** Since, in our framework, the fluid is Newtonian, incompressible and homogeneous, the fluid mass density $\rho_{flu}$ remains constant. We denote $\nu$ the fluid kinematic viscosity and $\eta = \rho_{flu} \nu$ its dynamic one. The flow is classically described by its velocity field $u(t, x) \in \mathbb{R}^3$ and the pressure $p(t, x) \in \mathbb{R}$, where $t \ge 0$ is the time and $x = (x_1, x_2, x_3) \in \mathbb{R}^3$ is the position. It is governed by the following equations:

$$ \rho_{flu}[\partial_t u + (u \cdot \nabla_x)u] = -\nabla_x p + \eta \Delta_x u + F, \quad t \in \mathbb{R}_+, \quad x \in \Omega_t, \quad (2.1) $$

$$ \nabla_x \cdot u = 0, \quad t \in \mathbb{R}_+, \quad x \in \Omega_t, \quad (2.2) $$

where $F$ is a vector field representing the forces acting on the fluid (gravity, aerosol retroaction). Moreover, to take into account the fact that the domain itself can move, we consider given time-indexed open sets $\Omega_t$ of $\mathbb{R}^3$. If we denote, for any $t \ge 0$, the displacement $A_t: \Omega_0 \to \Omega_t$ from the initial (reference) configuration, we shall assume that $(t, x) \mapsto A_t(x)$ is as smooth as needed, see [38, 19].

2.3. **Aerosol equation.** The distribution function $f: \mathbb{R}_+ \times \mathbb{R}^3 \times \mathbb{R}^3 \to \mathbb{R}_+$ depends on time $t$ and position $x$ of the particles, but also on their velocity $v$. In fact, $f$ can also depend on particle radius, temperature or other relevant quantities, as seen in [40, 43]. Let us emphasize that we do not take into account any phenomenon modifying the aerosol distribution regarding radius (no collision, no abrasion, etc.) or other physical quantities (no significant temperature variation during the breathing process, for instance). This ensures that the initial radius distribution is conserved with respect to time. Therefore, for the sake of simplicity, the aerosol is chosen monodispersed in size, meaning that $r$ is a parameter, see Remark 2 below. Each particle is assumed to remain spherical and hence, its mass $m$ is constant and satisfies $m = 4\pi r^3 \rho_{aer}/3$, where $\rho_{aer}$ is the constant volume mass of each particle.

The distribution function solves the Vlasov equation, i.e.

$$ \partial_t f + v \cdot \nabla_x f + \nabla_v \cdot (a f) = 0, \quad (2.3) $$

where $a(t, x, v)$ is the acceleration field undergone by the aerosol.

**Remark 1.** One can understand $f$ in two ways. It can be seen as

* a number density: $f(t, x, v) \, dx \, dv$ is the number of droplets located in the elementary volume of the phase space, centred at $(x, v)$ at time $t$;

* a probability density: if the zero-th moment of $f$ equals 1, $f(t, \cdot, \cdot)$ is the density function, at time $t$, of the probability of presence of the particles in the phase space.
**Remark 2.** We can take into account the case where the particles have different radii. Indeed, if the fluid velocity is given, all the terms involving the particles linearly depend on $f$. Hence, several distribution functions of different radii may be superponed to model a size-polydispersed distribution. This has implications in the numerical computations, see Remark 6 in Section 3.

2.4. **Interaction between the fluid and the aerosol.** The terms $F$ and $a$ must still be defined. In what follows, the gravitational effects will be neglected because we only focus on fluid-aerosol interaction. Both terms $F$ and $a$ model a drag force (or acceleration) between the fluid and the aerosol. We refer to Appendix A for a discussion about the drag force expression. We here use the Stokes law, which allows to write, for any $t$ and $x$,

$$a(t, x, v) = \frac{6\pi\eta r}{m} (u(t, x) - v), \qquad (2.4)$$

$$F(t, x) = -m \int_{\mathbb{R}^3} f(t, x, v) a(t, x, v) \, dv. \qquad (2.5)$$

Whether we take $F$ into account or not in this model is a real issue: it is called the aerosol retroaction on the fluid.

2.5. **Initial and boundary conditions.** Equations (2.1)–(2.5) must be supplemented by initial and boundary conditions. Recall that for any $t \ge 0$, the motion of the boundary is given by the displacement $A_t: \Omega_0 \to \Omega_t$ of the whole domain, from the initial (reference) configuration, $(t, x) \mapsto A_t(x)$ being as smooth as needed. The velocity of the domain at the point $x \in \partial\Omega_t$ is hence given by

$$w(t, x) := \dot{A}_t(A_t^{-1}(x)).$$

We assume that the fluid velocity fit $w$ on $\Gamma_t^{\text{wall}}$

$$u = w \text{ on } \Gamma_t^{\text{wall}}. \qquad (2.6)$$

Of course, when the domain is fixed, $w \equiv 0$ and we have a homogeneous Dirichlet boundary condition for $u$ on $\Gamma_0^{\text{wall}}$.

On both the inlet(s) and the outlet(s), we can assign Dirichlet or Neumann boundary conditions, for instance, if $u_t^{\text{in}}: \Gamma_t^{\text{in}} \to \mathbb{R}^3$, we can choose

$$u = u_t^{\text{in}} \text{ on } \Gamma_t^{\text{in}}, \quad (\nabla_x u + (\nabla_x u)^T) \cdot n_t - pn_t = 0 \text{ on } \Gamma_t^{\text{out}}. \qquad (2.7)$$

For a more realistic modelling in the airflow context, one can propose, as in [28, 8] the following strategy. In the proximal areas, the airflow is computed thanks to the Navier-Stokes equations whereas, in the distal part, it is described by a well chosen 0D boundary conditions, taking the diaphragm motion into account.

We choose an absorption boundary condition for the aerosol on the wall. In the kinetic formalism, it writes

$$(v - w) \cdot n_t < 0 \Rightarrow f = 0, \quad \text{on } \Gamma_t^{\text{wall}} \times \mathbb{R}^3. \qquad (2.8)$$

**Remark 3.** The boundary conditions on $\Gamma_t^{\text{wall}}$ are of course consistent with the respiration framework: indeed, the wall is coated with mucus and the aerosol particles deposit on the wall when they hit it.

Consider $u_0: \Omega_0 \to \mathbb{R}^3$ and $f_0: \Omega_0 \times \mathbb{R}^3 \to \mathbb{R}_+$ as initial data, i.e.

$$u(0, x) = u_0(x), \quad x \in \Omega_0, \quad f(0, x, v) = f_0(x, v), \quad x \in \Omega_0, \quad v \in \mathbb{R}^3. \qquad (2.9)$$

In a fixed domain, with homogeneous Dirichlet boundary conditions for the fluid and absorption for the spray, on $\partial\Omega_0$, we have the following proposition.

**Proposition 4.** Assume that $\Omega_t = \Omega_0$ for any $t \ge 0$, $u = 0$ on $\partial\Omega_0$ and $f = 0$ on $\partial\Omega_0 \times \mathbb{R}^3$ if $v \cdot n_0 < 0$. Then the energy of the whole aerosol-air system decreases.
*Proof.* Multiply respectively (2.3) by $mv^2/2$ and (2.1) by $u$. Then integrate the outcome equalities respectively on $\Omega_0 \times \mathbb{R}^3$ and $\Omega_0$, to get

$$ \frac{d}{dt} \left( \frac{m}{2} \iint_{\Omega_0 \times \mathbb{R}^3} f(t, \mathbf{x}, \mathbf{v}) \mathbf{v}^2 d\mathbf{v} d\mathbf{x} + \frac{1}{2} \int_{\Omega_0} \rho_{\text{flu}} |\mathbf{u}(t, \mathbf{x})|^2 d\mathbf{x} \right) \\ = - \iint_{\partial \Omega_0 \times \mathbb{R}^3} f(t, \mathbf{x}, \mathbf{v}) \mathbf{v} \cdot \mathbf{n}_0 d\mathbf{v} d\mathbf{x} \le 0, $$

since $f$ is nonnegative. The previous inequality ensures that the total kinetic energy of the coupled system decreases. □

**Remark 5.** If one considers other kinds of fluid boundary conditions, such as (2.7), the previous result does not hold and no energy bound can be derived. This lack of energy estimate comes from the Neumann boundary conditions for the Navier-Stokes system, and may also lead to numerical instabilities, see [25] for a review on this topic.

The nonlinearity in the Navier-Stokes equations and the strong coupling between the Vlasov and Navier-Stokes equations are the two major difficulties about the full system (2.1)–(2.7) from both mathematical and numerical viewpoints. System (2.1)–(2.9) was mathematically investigated in a fixed domain (global in time existence of weak solutions) in [30, 6, 11, 45]. In the two-dimensional case, uniqueness was also investigated [44]. Note that the approximation strategy used in the existence result established in [11] and especially the uncoupling process between the kinetic and fluid equations inspired the numerical scheme presented in the next section.

### 3. NUMERICAL SCHEME

We here propose a time-advancing scheme to solve the strongly coupled problem (2.1)–(2.7). We first uncouple the fluid and aerosol problems and solve the fluid part with a retroaction source term coming from the previous time step. Then we solve the kinetic part, using the updated fluid velocity to compute the drag force. Note that we may need time subcycling of the kinetic part to get an accurate value of the retroaction term.

The aerosol is computed thanks to a particle-in-cell (PIC) method [18, 20, 21, 42]. For the fluid, we use a Lagrange finite element method associated to an arbitrary Lagrangian-Eulerian (ALE) approach [32, 23], to handle the moving aspect. Those methods are briefly discussed in the next subsections. The most intricate part of this scheme probably lies in the coupling between (2.1) and (2.3) through (2.5).

Since the computational domain can move, we need to define the domain mapping $A_t: \Omega_0 \to \mathbb{R}^3$ and the associated velocity $w(t, \cdot): \Omega_t \to \mathbb{R}^3$. For instance, at each time $t$, $A_t$ can be computed from the boundaries movement as a solution to a Poisson problem set on $\Omega_t$.

Besides, our work is embedded in the C++ finite element library `LifeV2`, which previously owned numerical tools to handle both fixed and moving meshes, most classical finite element methods, and offered solvers for biological flows (Navier-Stokes, Darcy, etc.). If more details are needed, the reader is invited to refer to [37], see also [38, 19] about the fluid solvers.

In the following, we shall denote $T$ the final time of computation, and consider a regular subdivision $(t_n)_{0 \le n \le N}$ of $[0, T]$ with a step $\Delta t > 0$.

### 3.1. The Navier-Stokes equations.
Most features of the fluid solver we present below are standard. Nevertheless, we briefly explain how the whole computation is handled.

We discretize the Navier-Stokes equations (2.1)–(2.2) written in the ALE conservative form [38]. With our boundary conditions (2.6)–(2.7), it is given, for any $t \in [0, T]$, by

<footnote n="2">²Free software under LGPL license, jointly developed in four institutions: École Polytechnique Fédérale de Lausanne (Switzerland), Politecnico di Milano (Italy), Inria Paris-Rocquencourt (France) and Emory University (USA), see http://www.lifev.org/</footnote>
$$
\begin{align*}
& \varrho_{\text{flu}} \frac{d}{dt} \int_{\Omega_t} \mathbf{u} \cdot \boldsymbol{\beta}(t, \mathbf{x}) \, d\mathbf{x} + \varrho_{\text{flu}} \int_{\Omega_t} [ (\mathbf{u} - \mathbf{w}) \cdot \nabla_x \mathbf{u} ] \cdot \boldsymbol{\beta} \, d\mathbf{x} \\
& - \varrho_{\text{flu}} \int_{\Omega_t} (\text{div } \mathbf{w}) \mathbf{u} \cdot \boldsymbol{\beta} \, d\mathbf{x} + \eta \int_{\Omega_t} \nabla_x \mathbf{u} : \nabla_x \boldsymbol{\beta} \, d\mathbf{x} + \int_{\Omega_t} p(\text{div } \boldsymbol{\beta}) \, d\mathbf{x} = \int_{\Omega_t} \mathbf{F} \cdot \boldsymbol{\beta} \, d\mathbf{x}, \\
& \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \int_{\Omega_t} (\text{div } \mathbf{u}) \, \mu(t, \mathbf{x}) \, d\mathbf{x} = 0,
\end{align*}
$$

where $\boldsymbol{\beta}$ and $\mu$ are suitable test functions (transported by the ALE mapping $\mathcal{A}_t$ from reference test functions), satisfying $\boldsymbol{\beta} = 0$ on the wall and on the inlets.

We choose a backward Euler scheme and a semi-implicit treatment of the convective term. In the simpler case of a fixed domain, it reduces to a standard semi-implicit Euler scheme.

For the space discretization, we use a Lagrange finite element method. For a characteristic size $h > 0$, consider a tetrahedric mesh of $\Omega_0$, denoted by

$$
\mathcal{T}_0^h = \bigcup_{i=1}^{M_h} K_{0,i}^h,
$$

where each $K_{0,i}^h$ is a tetrahedron. For each $t$, the tetrahedric mesh $\mathcal{T}_t^h$ of $\Omega_t$ is the union of $M_h$ tetrahedra whose vertices are transported from the reference configuration $\mathcal{T}_0^h$ by the discrete ALE mapping $\mathcal{A}_t^h$. This mapping should preserve the tetrahedric structure of the mesh, see [38] for details.

Then we define the basis functions, at each time $t^n$, $(\varphi_k^n)_{1\le l\le P_h}$ for velocities, and $(\psi_l^n)_{1\le l\le P_h}$ for pressures, transported from $\mathbb{P}1 - \mathbb{P}1$ reference basis functions on $\mathcal{T}_0^h$. Once again, $\mathcal{A}_t^h$ should preserve the chosen finite element space. Due to the choice of a $\mathbb{P}1 - \mathbb{P}1$ setting, we use a stabilized formulation, that we do not detail here.

Then, at time $t^n$, we approximate the unknowns on $\mathcal{T}_t^n$ by

$$
\mathbf{u}^n(\mathbf{x}) = \sum_{k=1}^{N_h} u_k^n \varphi_k^n(\mathbf{x}), \quad p^n(\mathbf{x}) = \sum_{\ell=1}^{P_h} p_\ell^n \psi_\ell^n(\mathbf{x}),
$$

for any $\mathbf{x} \in \Omega_{t^n}$. The unknowns then become the following column vectors

$$
U^n = (u_1^n, \dots, u_{N_h^n})^\mathrm{T}, \quad \Pi^n = (p_1^n, \dots, p_{P_h^n})^\mathrm{T}.
$$

Consequently, from $t^n$ to $t^{n+1}$, we have to solve the following linear system

$$
\begin{bmatrix} D^{n+1} & (B^{n+1})^{\mathrm{T}} \\ B^{n+1} & 0 \end{bmatrix} \begin{bmatrix} U^{n+1} \\ \Pi^{n+1} \end{bmatrix} = \begin{bmatrix} F^n \\ 0 \end{bmatrix} + \frac{\varrho_{\mathrm{flu}}}{\Delta t} \begin{bmatrix} M^n U^n \\ 0 \end{bmatrix}, \quad (3.1)
$$

where

$$
B^n = \left( - \int_{\Omega_{t^n}} \psi_i^n \nabla_x \boldsymbol{\varphi}_j^n \, d\mathbf{x} \right)_{1 \le i \le P_h, 1 \le j \le N_h},
$$

$$
M^n = \left( \int_{\Omega_{t^n}} \varphi_i^n \cdot \varphi_j^n \, d\mathbf{x} \right)_{1 \le i,j \le N_h},
$$

$$
F^n = \left( \int_{\Omega_{t^n}} \mathbf{F}^n \cdot \boldsymbol{\varphi}_i^n \, d\mathbf{x} \right)_{1 \le i \le N_h}^{\mathrm{T}},
$$

$$
D^{n+1} = \frac{\varrho_{\text{flu}}}{\Delta t} M^{n+1} - \varrho_{\text{flu}} C^{n+1/2} + \eta A^{n+1},
$$

with

$$
A^n = \left( \int_{\Omega_{t^n}} \nabla_x \varphi_i^n : \nabla_x \varphi_j^n \, d\mathbf{x} \right)_{1 \le i,j \le N_h},
$$

$$
C^{n+1/2} = \left( \sum_{k=1}^{N_h} (u_k^n - w_k^{n+1}) \left[ \int_{\Omega_{t^{n+1}}} \left( \varphi_k^{n+1} \cdot \nabla_x \varphi_j^{n+1} \right) \cdot \varphi_i^{n+1} \, d\mathbf{x} \right] \right)_{1 \le i,j \le N_h}.
$$
The definition of $F^n$ on $\Omega_{t^n}$ is given in 3.3 below.

The semi-implicit treatment of the convective term appears in the expression of $C^{n+1/2}$, where $u_k^n$ is used with the basis functions at time $t^{n+1}$. In the same way, note that (3.1) involves nodal quantities at time $t^{n+1}$. In particular, vectors in the right-hand side of (3.1), are defined using the finite element coordinates of quantities on $\Omega_{t^n}$, and further transported by the ALE mapping on $\Omega_{t^{n+1}}$ (on the corresponding nodes).

## 3.2. The Vlasov equation.
Unlike the finite element method, the particle method does not provide an approximation of $f$ on the mesh nodes. More precisely, the distribution function is computed as a weighted sum of Dirac masses in the positions and velocities of the numerical particles, i.e., in a measure sense,

$$f(t, \mathbf{x}, \mathbf{v}) = \sum_{p=1}^{N_{\text{num}}} \omega_p \delta_{\mathbf{x}_p(t)} \otimes \delta_{\mathbf{v}_p(t)}(\mathbf{x}, \mathbf{v}),$$

where the number of numerical particles $N_{\text{num}}$ is initially chosen by the user, $\omega_p$ is called the representativity of numerical particle $p$, and $t \mapsto (\mathbf{x}_p(t), \mathbf{v}_p(t))$ is the trajectory, in the phase space, of $p$. Hence, we just have to compute the characteristics to get an approximation of $f$ from its initial datum. Note that, as soon as $\mathbf{x}_p(t)$ does not belong to $\Omega_t$ any longer, that means that particle $p$ is deposited on $\Gamma_t^{\text{wall}}$ or went out of the domain through $\Gamma_t^{\text{in}}$ or $\Gamma_t^{\text{out}}$. From a numerical point of view, $N_{\text{num}}$ is not modified but the outgoing particles are not treated any more.

For practical purposes, $N_{\text{num}}$ is very small with respect to the average number of physical aerosol particles $N_{\text{aero}}$, which makes the computations much less costly. In fact, we can write

$$N_{\text{aero}} = \sum_{p=1}^{N_{\text{num}}} \omega_p,$$

which gives an order of magnitude of $\omega_p$. For instance, if we choose $N_{\text{aero}} = 10^{10}$ (see 2.3) and $N_{\text{num}} = 10^3$, we could set $\omega_p = 10^7$ for all numerical particles. The relevance of such a choice is discussed later in Section 5.

As long as the particle $p$ remains in the domain, its coordinates in the space phase solve the Cauchy problem

$$\dot{\mathbf{x}}_p(t) = \mathbf{v}_p(t), \quad \dot{\mathbf{v}}_p(t) = \mathbf{a}(t, \mathbf{x}_p(t), \mathbf{v}_p(t)), \quad 0 \le t \le T,$$

with chosen initial data $\mathbf{x}_p(0)$ and $\mathbf{v}_p(0)$. This problem is solved with a first-order semi-implicit Euler scheme

$$\mathbf{v}_p^{n+1} = \mathbf{v}_p^n + \Delta t \tilde{\mathbf{a}}^{n+1/2}, \quad \mathbf{x}_p^{n+1} = \mathbf{x}_p^n + \Delta t \mathbf{v}_p^{n+1}, \quad (3.2)$$

where $\tilde{\mathbf{a}}^{n+1/2}$ is the drag acceleration depending on quantities at both times $t^n$ and $t^{n+1}$. This is explained in the next subsection, because it is related to the coupling between both phases.

Eventually, we must emphasize that this method needs averaging, since it relies on statistical physics. Indeed, if the particles are injected with a uniform distribution at $\Gamma_t^{\text{in}}$, we must proceed with several initial numerical distributions and take the computations average to locate the deposition areas, for instance.

**Remark 6.** As we already stated, we only describe the scheme (and the model) for a monodispersed (in radius) aerosol. It is of course possible to consider numerical particles with various radii. Hence, it is possible to reproduce the radius distribution of the particles in the aerosol. It is very useful when one a priori knows how an aerosol nebulizer generates particles. Of course, if one takes into account the size changes of aerosol particles due to the air humidity [33], this strategy must be adapted.
3.3. **The coupling.** Let us now focus on the coupling between the Vlasov and Navier-Stokes equations. We here choose to use an explicit time marching scheme, so that the fluid and aerosol parts are basically solved once per time step. It enables to reduce the computational cost. Nevertheless, it may lead to stability issues we investigate afterwards. Note that, as already stated, the proof of existence of weak solutions in [11] is based on the same kind of decoupling strategy. The whole process is quite similar to the ones in [40] or [27].

We assume that, at time $t^n$, $U^n$, $\Pi^n$, $(\mathbf{x}_p^n)_{1\le p\le N_{\text{num}}}$, $(\mathbf{v}_p^n)_{1\le p\le N_{\text{num}}}$ and $F^n$ are known. Of course, when $n = 0$, $F^0$ must be first computed from the initial data before starting the time loop. Then we solve (3.1) to obtain $U^{n+1}$ and $\Pi^{n+1}$. Next, to advance the aerosol in time, we compute a discrete particle acceleration denoted by $\tilde{\mathbf{a}}^{n+1/2}$ and defined by

$$ \tilde{\mathbf{a}}^{n+1/2} = \frac{6\pi\eta r}{m}(\tilde{\mathbf{u}}^{n+1}(\mathbf{x}_p^n) - \mathbf{v}_p^{n+1}), $$

where the quantity $\tilde{\mathbf{u}}^{n+1}(\mathbf{x}_p^n)$ is the fluid velocity at the position of particle $p$, at the previous time step. Since the fluid velocity is only known at mesh points, we need to interpolate it in order to evaluate $\tilde{\mathbf{u}}^{n+1}(\mathbf{x}_p^n)$.

Note that $\mathbf{x}_p^n \in \Omega_{t^n}$ and that $U^{n+1}$ lies in $\Omega_{t^{n+1}}$. First, particle $p$ may not be in $\Omega_{t^{n+1}}$ any more. In this case, the particle is considered to be deposited, hence taking the boundary condition of $f$ on $\Gamma_t^{\text{wall}}$ into account. Second, we have to find out in which cell particle $p$ is, and its associated barycentric coordinates. This is possible thanks to a locating algorithm [26]. Note that, apart from the initial locating processes, this algorithm should converge in a very small number of iterations, so it is not so costly. Third, we perform a linear interpolation of $U^{n+1}$ on the mesh points at time $t^{n+1}$ to get $\tilde{\mathbf{u}}^{n+1}(\mathbf{x}_p^n)$.

Once we know $\tilde{\mathbf{a}}^{n+1/2}$, we compute $\mathbf{v}_p^{n+1}$ and $\mathbf{x}_p^{n+1}$ thanks to (3.2), and eventually $\mathbf{F}^{n+1}$, as

$$ \mathbf{F}^{n+1}(\mathbf{x}) = -m \sum_{p=1}^{N_{\text{num}}} \omega_p \frac{6\pi\eta r}{m} (\tilde{\mathbf{u}}^{n+1}(\mathbf{x}_p^{n+1}) - \mathbf{v}_p^{n+1}) \delta_{\mathbf{x}_p^{n+1}}(\mathbf{x}), \quad \mathbf{x} \in \Omega_{t^{n+1}}. $$

Note that, to compute $\mathbf{F}$ at time $t^{n+1}$, we do not use $\tilde{\mathbf{a}}^{n+1/2}$. Consequently, the numerical total momentum is not conserved any more, which may imply stability issues.

**Remark 7.** The fluid time step may not be suitable to precisely follow the discrete particle trajectories. As seen on Figure 2, without subcycling, the particle goes across several cells during a same time step, whereas the high fluid velocity on the black nodes spins the particle trajectory. In this case, one should introduce a time subcycling strategy.

<figure><img src="image_3.png" /><figcaption>FIGURE 2. A particle with (solid line) and without (dashed line) time subcycling.</figcaption></figure>
**Remark 8.** In order for the retroaction term to be more accurate, we may use the smallest subcycle time step for the fluid phase too. Indeed, each particle may be able to exchange its momentum with the fluid in each cell it goes across.

In what follows, we numerically study accuracy and stability properties of our scheme. Since we aim to model the aerosol deposition in the human lung, the air is chosen as the ambient gas at temperature 310 K, so that $\rho_{flu} = 1.204 \times 10^{-3}$ g.cm$^{-3}$ and $\eta = 1.85 \times 10^{-4}$ g cm$^{-1}$ s$^{-1}$. The previous values can be found, for instance, in [1].

We first investigate in Section 4 the no-retroaction case in a fixed domain. We are interested, for instance, in sensitivity with respect to the mesh size, time step, initial particle locations. Then, in Section 5, we perform the same kind of tests in the retroaction case, again in a fixed domain. We moreover focus on representativity issues appearing when dealing with the retroaction term. Finally, Section 6 is dedicated to the moving domain case without retroaction. In particular, we exhibit an exact solution to (2.1)–(2.2) and compare it with our numerical solution. In every case, one of the main issues is to accurately predict the deposition map of the aerosol. We conclude our study by using our numerical code on a typical lung geometry.

## 4. NUMERICAL TESTS WITHOUT RETROACTION

We first focus on the behaviour of our numerical code when the action of the particles on the fluid is neglected. We consider an orthonormal system ($O, e_1, e_2, e_3$). The computational domain is a cylinder of axis ($O, e_3$), length $L = 5$ cm and radius $R = 0.2$ cm, centred at $O$. The airflow follows the Poiseuille law, with velocity $u$ going along the cylinder axis, i.e.

$$u(t, x) = U \left(1 - \frac{x_1^2 + x_2^2}{R^2}\right) e_3, \quad U \ge 0.$$

In particular, $u$ does not depend on $t$ and $x_3$, only on the transverse coordinates $x_1$ and $x_2$.

Particles evolve in the airflow. They have the same radius $r$, mass density $\rho_p$, and are injected in the fluid at time $T_{inj} = 0.2$ s with the same initial velocity $V^0$. Fig. 3 sums up the situation.

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Initial velocity $V^0$ for the computations : solid line for 4.1, dashed line for 4.2 and following.</figcaption></figure>

### 4.1. Behavior of the particle velocities.
In order to validate the characteristics method implementation in our code, let us consider here the case when there is only one numerical particle in the fluid. The Stokes law allows to obtain an analytic expression of the particle velocity. Since there is no action of the particle on the fluid, the particle asymptotically follows the streamlines of the fluid. From the analytic form of the fluid velocity, we can
deduce one for the particle velocity. Let us denote by $\mathbf{X}(t)$ and $\mathbf{V}(t)$ the location and velocity of the particle at time $t$, and choose $\mathbf{V}^0 = V^0 \mathbf{e}_3$, $V^0 \in \mathbb{R}$. Consequently, at each time $t$, $\mathbf{V}(t)$ is parallel to the line $\mathbb{R}\mathbf{e}_3$. Its coordinate $V(t)$ along this line is then given by

$$V(t) = \left[ V^0 + \mathbf{u}(\mathbf{X}(t)) \cdot \mathbf{e}_3 \int_{T_{\text{inj}}}^{t} \frac{e^{s/\tau}}{\tau} \mathrm{d}s \right] e^{-t/\tau},$$

where

$$\tau := \frac{2r^2 \rho_p}{9\eta}.$$

Hence, we get

$$V(t) = \mathbf{u}(X(t)) \cdot \mathbf{e}_3 + [V^0 - \mathbf{u}(X(t)) \cdot \mathbf{e}_3]e^{-(t-T_{\text{inj}})/\tau}, \quad (4.1)$$

as long as the particle remains in the cylinder.

4.1.1. *Motionless fluid.* We first assume that the fluid does not move, i.e. $U = 0$. Consequently, (4.1) becomes

$$V(t) = V^0 e^{-(t-T_{\text{inj}})/\tau},$$

the fluid slows down the particle. For standard physical parameters (water droplet in the air), the relaxation time $\tau$ is very small. Since we are only interested in the code validation, let us use a non realistic mass density for the particle, i.e. 10 g/cm$^3$, to allow, in the same time, a significant relaxation time and a small particle Reynolds number.

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Analytic/numerical particle velocities: (a) plots, (b) error.</figcaption></figure>

Choosing $V^0 = 5$ cm/s and $T_{\text{inj}} = 0.2$ s, we get $V(t) = 5e^{-3.33(t-0.2)}$ for any $t \ge T_{\text{inj}}$. In Fig. 4a, the analytic and numerical particle velocities are plotted with respect to time, and Fig. 4b shows the absolute error between both velocities.

4.1.2. *Nonzero Poiseuille profile.* We now assume that the maximum fluid velocity $U$ is positive, here $U = 3$ cm/s. To compute $V(t)$ for each particle thanks to (4.1), we need to find the value of the fluid velocity at the particle location. Thanks to the Poiseuille law, it is enough to know the distance between the particle and the cylinder axis.

Consider, for instance, two particles. One is located on the axis and the other one 0.1 cm away (in the radial direction). Let us take $V^0 = -0.1$ cm/s, opposite to the fluid flow. Then Fig. 5 shows again the very good agreement between the analytic and numerical velocities of the particle.

In the remainder of this section, we investigate the numerical sensitivity of our scheme with respect to various parameters of the computations, such as the mesh, the time step, etc. We study an aerosol with 3000 numerical particles sharing the same initial transverse velocity $\mathbf{V}^0 = V^0 \mathbf{e}_1$ with $V^0 = 6.5$ cm/s. Fig. 3 may again allow to understand the situation. The chosen final time is 1.4 s. Twenty different initial space distributions of
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Analytic/numerical velocities for a particle on the axis and another one 0.1 cm away from the axis.</figcaption></figure>

3000 numerical particles are used for the particle method: they are uniformly drawn inside a disk of radius 0.16 cm at $x_3 = -L/2$. All the particles have the same radius, i.e. 50 µm.

4.2. **Mesh sensitivity.** We consider seven different meshes, from approximately 10,000 tetrahedra to 640,000. More precisely, we give, in Table 1, the number of tetrahedra and a corresponding average cell size value ($h$) obtained as the cube root of the ratio between the cylinder volume and the number of cells.

<table><caption>TABLE 1. Various meshes under study.</caption><tr><td>Number of cells</td><td>9,954</td><td>28,800</td><td>67,200</td><td>80,136</td><td>130,200</td><td>295,416</td><td>645,150</td></tr><tr><td>$h$ (in cm)</td><td>0.0398</td><td>0.0279</td><td>0.0211</td><td>0.0199</td><td>0.0169</td><td>0.0128</td><td>0.0099</td></tr></table>

Due to computational costs, the range of explored cell sizes is narrowed between 0.01 and 0.04. We perform the simulations for all seven different meshes with 20 different initial distributions of particles, with a time step $\Delta t = 0.0027$ s. The finest mesh computation is taken as the reference one.

First, we consider a given initial particle distribution. For each numerical particle, we recover the time evolution of the distance between its reference position and the ones obtained using coarser meshes. Then we compute the average of this distance with respect to all the numerical particles. We iterate this process for each of the 20 initial particle distributions, which allows to define the average distance as a statistical quantity. For instance, in Fig. 6, the average particle distance between computations on meshes 80 and 640 is shown, as well as the associated standard deviation.

After an initial phase, we observe a linear growth of the average distance error. Hence, the growth ratios are quantities of interest. They can be seen as line slopes. In Fig. 7, we plot the slope values with respect to the typical cell size $h$, and once again with the statistical error.

Our numerical scheme behaves with a remarkable statistical stability. On the deposition phenomenon, the variation of the average fraction of deposited particles (and the standard deviation) at final time does not seem significant in terms of meshes, see Fig. 8. The percentage of 54.7% is clearly a satisfying value, and we can note that using a coarser mesh seems to overestimate the deposition effect a little.

Eventually, when focusing on the coarsest and the finest meshes, let us emphasize that the maximal distance between the positions of the same particle at deposition time in the two different meshes is only 0.015 cm.
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Average distance between the coarsest and finest meshes.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7. Growth ratios on the average distance error.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 8. Fraction of deposited particles with respect to $h$.</figcaption></figure>

4.3. **Time step sensitivity.** In order to study the sensitivity of the computations with respect to the time step, we use the finest mesh from Table 1 with 645,150 tetrahedra. We study three situations: the reference one with $\Delta t/4 = 0.000675$ s, and two others with
$\Delta t/2$ and $\Delta t$. We then compute the distance between the locations of each numerical particle in the reference situation and one of the two others, at each multiple of $\Delta t$, up to final time 1.4s, and take the average over all the particles.

<figure><img src="image_10.png" /><figcaption>FIGURE 9. Comparisons on averaged particle locations on (a) $\Delta t$ and $\Delta t/4$, (b) $\Delta t/2$ and $\Delta t/4$.</figcaption></figure>

Then Fig. 9a–b respectively show the average (over the 20 randomly drawn space initial distributions) position errors for $\Delta t/4$ and $\Delta t$, and $\Delta t/4$ and $\Delta t/2$, respectively. As usual, the average and the standard deviation are displayed. Both results are satisfactory because the errors seem to have a linear behaviour with respect to physical time $t$.

Furthermore, we can check the number of deposited particles does not depend on the time step in the previous tests. More precisely, the maximal distance between the impact positions in computations with $\Delta t/4$ and $\Delta t/2$ is 0.004 cm and in computations with $\Delta t/4$ and $\Delta t$ is 0.008 cm.

Notice that this stability behavior for the deposition of particles is clearly related to the subcycling strategy mentioned in Remark 7, that we systematically used in these computations. Without the refinement of the time step for the particles, some of them could ignore artificially the boundary layer of the fluid, going through several cells in one time step and go out of the mesh. This would indeed be an artefact since, reducing the time step, these particles would start to perceive more accurately the velocity of the fluid near the boundary, whence reducing the deposition ratio.

## 4.4. Sensitivity with respect to the initial particle locations.
This test focuses on the sensitivity with respect to the initial space distribution of the numerical particles. We choose $\Delta t = 0.0027$ s and consider every mesh from Table 1. For a given initial configuration of particles, we perturb the particle locations in the following way. Each particle is randomly (with a uniform law) relocated in a small disk of radius $h/2$ centred at the particle initial position. The small radius value ensures that the particle remains in the cylinder after the perturbation. Let us emphasize that the perturbation depends on the considered mesh, through $h$.

We perform the computations for 20 different perturbations of the same kind, and eventually compute the maximal distance between the locations of the particles at final time with respect to the reference configuration, averaged on the 20 draws. Fig. 10 shows a remarkable stability with respect to the perturbation of the initial distribution.

# 5. NUMERICAL TESTS WITH RETROACTION

In this section, the particles exert a retroaction force on the fluid, and the coupling between the Navier-Stokes and Vlasov equations is strong. The retroaction effect was already discussed in a two-dimensional setting [13]: there were situations when $F$ could not be neglected at all (big and numerous particles, high initial relative particle velocities).
<figure><img src="image_11.png" /><figcaption>FIGURE 10. Effect of a mesh-dependent perturbation of the initial particle distribution on the average distance between particles at final time.</figcaption></figure>

The discussion about *F* allows to classify respiratory aerosols following [40]. They can be very thin sprays (interaction between particles, spray volume fraction and *F* are negligible) or thin sprays (interaction between particles and spray volume fraction are negligible, but *F* is not).

In this section, the computations for both phases are performed at each time step, whereas, in Section 4, we were able to first compute the full time evolution of the fluid flow, and afterwards the movement of each numerical particle as a post-treatment. The computational cost of course significantly increases, since we now need to know the aerosol movement at each time to compute the retroaction source term in the Navier-Stokes equations (2.1). Let us emphasize that we should consider relevant situations in which the aerosol has a significant effect on the fluid flow.

We performed the same kind of sensitivity studies as in the previous section. We do not provide results regarding sensitivity with respect to the mesh, because they are very similar to the ones shown in 4.2, but we choose to present a result regarding the time step sensitivity. We also tackle the key issue on particles representativity.

5.1. **Time step sensitivity.** Let us mimic the test from 4.3, with the same value of $\Delta t = 0.0027$ s. We here use the 80,136-tetrahedron mesh for the computations. Fig. 11a–b can easily be compared to Fig. 9a. The average distance between the locations of particles computed with $\Delta t$ and $\Delta t/4$ is a little bit bigger when retroaction is on, for both small and high representativities. Nevertheless, it really remains of the same order of magnitude.

<figure><img src="image_12.png" /><figcaption>FIGURE 11. Comparisons on averaged particle locations on $\Delta t$ and $\Delta t/4$ with (a) $\omega_p = 1$, (b) $\omega_p = 10^4$.</figcaption></figure>
## 5.2. Representativity of the numerical particles.
In this subsection, we discuss the numerical approximation of the retroaction term, and mainly the influence of the particle representativity. First, we deal with a given number $N_{\text{num}}$ of numerical particles. We investigate the behaviour of the total kinetic energy (of the coupled system) for various values of the representativity. This implies that we do not use the same number of physical particles in each numerical experiment. Next, we aim to exhibit an optimal value of the representativity for a given realistic total number particles.

### 5.2.1. Representativity issues for a given number of numerical particles.
We inject 3000 numerical particles with axial velocity 3 cm/s at $z = 0.2$ in a disk with radius 0.12 cm, and we consider an initially motionless fluid. The boundary conditions are chosen as in Proposition 4. Since retroaction is taken into account, the fluid velocity should be non zero near the particles shortly after the injection. According to Proposition 4, the total kinetic energy should decrease. The computations are first performed for a chosen time step equal to 0.0027 s, and various representativities. In Fig. 12, we observe that the total kinetic energy does not decrease when representativity, and consequently the number of physical particles, grow.

<figure><img src="image_13.png" /><figcaption>FIGURE 12. Total kinetic energy for various representativities, with $\Delta t = 0.0027$.</figcaption></figure>

Moreover, performing the same test for the larger time step 0.004 s, we observe in Fig. 13 that the total kinetic energy does not decrease anymore for $\omega_p = 5.10^6$, whereas it did for $\Delta t = 0.0027$.

This numerical behaviour may be explained by the explicit treatment of the retroaction as a source term in the Navier-Stokes equations. This explicit treatment may induce an unphysical energy production. More precisely, when computing velocity and pressure at $t^{n+1}$, the source term writes

$$F^n(\mathbf{x}) = -m \sum_{p=1}^{N_{\text{num}}} \omega_p \frac{6\pi\eta r}{m} (\tilde{u}^n(\mathbf{x}_p^n) - \mathbf{v}_p^n) \delta_{x_p^n}(\mathbf{x}), \quad \mathbf{x} \in \Omega_{t^n},$$

where $\omega_p$ clearly appears as a key parameter. Furthermore, the comparison between both cases with two different time steps suggests the existence of a CFL-like condition involving the product $\omega_p \Delta t$ among others.
<figure><img src="image_14.png" /><figcaption>FIGURE 13. Total kinetic energy for two “critical” representativities, with $\Delta t = 0.004$.</figcaption></figure>

5.2.2. *Representativity issues for a given number of physical particles.* In this subsection, the number of physical particles is constant: $10^{10}$, which is realistic for a standard nebulizer. We aim to determine how many numerical particles should be involved to represent those physical particles. By choosing very few numerical particles with a huge representativity in our computations, the space distribution of the aerosol particles may not well taken into account. Consequently, the aerosol retroaction, which only has a local effect, could be underestimated, leading to a non physical behaviour. Hence, the number of numerical particles should be a compromise between both the numerical cost and the space distribution of the particles in the fluid.

Here we consider three cases:

• 100 numerical particles, each representing $10^8$ physical particles,

• 1000 numerical particles, each representing $10^7$ physical particles,

• 10,000 numerical particles, each representing $10^6$ physical particles.

The 1000 and 10,000-particle test cases are obtained thanks to perturbation of the 100-particle one. More precisely, we build a random distribution of 100 particles in a disk of radius 0.06 cm at $z = 0.2$. Then we create a new distribution by replacing each of the 100 particles by 10 particles randomly chosen in a disk of radius 0.06 cm around them. This process is repeated to get the 10,000-particle situation from the 1000 one.

We check the fluid velocity at two different mesh nodes close to the z-axis, near the initial location of the particles at $z = 0.238095$, and away from this same location at $z = 0.595238$. More precisely, in Fig. 14–15, we show the numerical error on the fluid velocity norm at both nodes, comparing the 1000-particle case with the other ones. In both situations, the results are quite similar when dealing with 1000 and 10,000 particles, whereas they are significantly different when dealing with 1000 and 100 particles. Consequently, it is first clear that the choice of 100 numerical particles is not relevant to represent $10^{10}$ physical particles. Second, since the error between the 1000 and 10,000-particle situations remains small (of order of magnitude $10^{-2}$ cm/s), the choice of 1000 numerical particles appears as a good compromise in terms of numerical cost and accuracy. Note that the tests have been performed for several particle draws, and the standard deviation in any case is upper-bounded by $10^{-6}$.

## 6. NUMERICAL TESTS IN A MOVING DOMAIN

Let us now investigate the moving domain situation. At initial time, the domain is the cylinder from Figure 3. We study a case where the domain motion and the fluid behaviour can be analytically computed, see Proposition 9 and its proof within Appendix B. The
<figure><img src="image_15.png" /><figcaption>FIGURE 14. Fluid velocity error near the particles w.r.t. time, comparison with 1000 numerical particles.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 15. Fluid velocity error away from the particles w.r.t. time, comparison with 1000 numerical particles.</figcaption></figure>

considered time dependent domain is given, for any $t$, by

$$ \Omega_t = \{ (\kappa(t)r \cos \theta, \kappa(t)r \sin \theta, \zeta_z(t)) \mid 0 \le r \le R, 0 \le \theta \le 2\pi, -L \le 2z \le L \}, $$

where we set

$$ \kappa(t) = \exp(\lambda t), \quad \text{with} \quad \lambda = \log(0.9)/1.4 \approx -0.075. $$

A solution to the Navier-Stokes equations writes, using cylindrical coordinates,

$$ \mathbf{u}(t, r, z) = u_r \mathbf{e}_r + u_z \mathbf{e}_z = \frac{\dot{\kappa}(t)}{\kappa(t)} r \mathbf{e}_r - 2 \frac{\dot{\kappa}(t)}{\kappa(t)} z \mathbf{e}_z, $$

$$ p(t, r, z) = -\frac{r^2 \ddot{\kappa}(t)}{2 \kappa(t)} - \frac{z^2}{\kappa(t)^2} [3\dot{\kappa}(t)^2 - \ddot{\kappa}(t)\kappa(t)]. $$

We checked that the Navier-Stokes ALE solver accurately reproduced the previous solution. Indeed, using the coarsest mesh and time step, the error on the fluid velocity is at most of order $7.10^{-3}$ cm/s.

We can also obtain the trajectory, in the phase space, of each numerical particle in a semi-analytical way, thanks, for instance, to the python function expm. Let us use again the notations $\mathbf{X}(t)$ and $\mathbf{V}(t)$ introduced in 4.1. Note that $\mathbf{V}$ does not remain parallel to $\mathbf{e}_z$ anymore. We can write

$$ \dot{\mathbf{X}}(t) = \mathbf{V}(t), \quad \dot{\mathbf{V}}(t) = -\frac{1}{\tau}(\mathbf{V}(t) - \mathbf{u}(t, \mathbf{X}(t))), $$

where we arbitrarily set $\tau = 2$: it may not be a physical value, but allows to study the particle trajectory in a reasonable time scale. The previous differential system can be
rewritten with matrix notations and solved by computing a matrix exponential, i.e.

$$ \begin{pmatrix} \mathbf{X}(t) \\ \mathbf{V}(t) \end{pmatrix} = \exp(At) \begin{pmatrix} \mathbf{X}(0) \\ \mathbf{V}(0) \end{pmatrix}, $$

where

$$ A = \begin{pmatrix} 0 & I_3 \\ A_1 & -I_3/2 \end{pmatrix}, \quad A_1 = \frac{\lambda}{2} \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & -2 \end{pmatrix}, $$

and $I_3$ denotes the identity matrix of $\mathbb{R}^3$.

In what follows, we perform various numerical experiments without retroaction to investigate the solver efficiency when taking the mesh motion into account. Our numerical scheme uses a decoupling method, where the fluid-aerosol system and the domain are sequently solved once per time step. At each time step, we proceed in the following way:

1. move the domain,

2. solve the fluid equation,

3. locate and move the particles thanks to (3.2).

Other implicit or explicit strategies may have been considered. In particular, the order of the different steps is a crucial issue with respect to the particle deposition. Hence, we check that the chosen strategy is efficient by comparing our numerical solution with the semi-explicit one, see Fig. 16.

<figure><img src="image_17.png" /><figcaption>FIGURE 16. Comparisons on averaged particle locations with the semi-analytical reference solution.</figcaption></figure>

We choose the same initial conditions and number of numerical particles as in the lower part of page 10 describing the test case in a fixed domain.

6.1. **Mesh sensitivity.** We choose 0.0027 s as the time step for this test case. We observe, on Fig. 17, an almost constant deposition ratio with respect to the mesh size.

Next we compare, on Fig. 18, the deposited particle locations at both final time and deposition time (for each particle). Note that this deposition time is the one obtained in the reference situation. This fact may explain the slight difference in the distances observed for each mesh, comparing cases a and b.

6.2. **Time step sensitivity.** Fig. 19 shows a very good agreement between the three plots, for the considered time steps. The test case is performed in the 80, 136-tetrahedron mesh.
<figure><img src="image_18.png" /><figcaption>FIGURE 17. Fraction of deposited particles for various mesh size.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>FIGURE 18. Comparisons on averaged particle locations with the semi-analytical reference solution at (a) final time, and (b) deposition time in the reference situation.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>FIGURE 19. Fraction of deposited particles for various time steps: $\Delta t/4$, $\Delta t/2$ and $\Delta t$.</figcaption></figure>

## 7. CONCLUSIONS AND PROSPECTS

In this paper, we proposed a model to describe the interaction of a thin spray in an incompressible viscous flow in order to accurately predict the aerosol deposition and possibly take into account the aerosol retroaction on the fluid and/or some wall motion. A non
linear coupled Vlasov-Navier-Stokes system was thus considered and explicitly discretized in time. For the space approximation, a PIC method was used for the kinetic equation, whereas an ALE finite element method was considered for the fluid part. We numerically demonstrated the stability of our scheme in various situations: fixed or moving domains, with or without retroaction. In particular, we highlighted the existence of a CFL-like condition in the retroaction case. Note that our numerical simulations showed very good convergence properties for the deposition rate, with respect to the mesh size for instance. One crucial point here lies in the sub-cycling strategy used to advance the particles in time, enabling an accurate prediction of the deposition phenomenon, which is a crucial issue of our study. We also want to underline that considering a kinetic modelling of the aerosol is appropriate in this context: a phase field description, even if numerically less costly, does not allow an easy description of deposition, whereas following each individual particles may be too expansive, in particular when retroaction is not negligible. Consequently, our model and our explicit numerical strategy represent a good compromise and are a good candidate to perform efficient in silico experiments of aerosol deposition in complex situations such as lung modelling.

## APPENDIX A. ABOUT THE DRAG FORCE

In [40], O'Rourke gives a general expression of the drag force of the air on the aerosol. In terms of acceleration, it reads, for any $t$, $x$ and $v$,

$$ \boldsymbol{a}(t, \boldsymbol{x}, \boldsymbol{v}) = \frac{\pi r^2}{8m} \rho_{\text{flu}} C_{\text{D}} |\boldsymbol{u}(t, \boldsymbol{x}) - \boldsymbol{v}| (\boldsymbol{u}(t, \boldsymbol{x}) - \boldsymbol{v}), \quad (A.1) $$

where $C_D$ is a dimensionless parameter called the drag coefficient. This coefficient is semi-empirically computed using the particle Reynolds number, which allows to measure the ratio between inertia and viscosity forces on the particle:

$$ \mathrm{Re}_p := \frac{2}{\eta} \rho_{\text{flu}} r |\boldsymbol{u} - \boldsymbol{v}|. $$

Then $C_D$, which also depends on $\boldsymbol{u}$ and $\boldsymbol{v}$, is given by Schiller and Naumann's law [14], also used in [24, 27], in the KIVA codes [5, 4, 2, 3, 41] and for complex fluid-particle mixtures [12, 35]:

$$ C_D = \begin{cases} \frac{24}{\mathrm{Re}_p} \left( 1 + \frac{1}{6} \mathrm{Re}_p^{2/3} \right), & \text{if } \mathrm{Re}_p < 1000, \\ 0.424, & \text{if } \mathrm{Re}_p \ge 1000. \end{cases} $$

When $\mathrm{Re}_p$ is small, which is the case in the respiration framework, one can also use $C_D = 24/\mathrm{Re}_p$, see [29]. This formula immediately allows to recover the Stokes law from (A.1).

## APPENDIX B. ANALYTIC SOLUTIONS TO THE INCOMPRESSIBLE NAVIER-STOKES EQUATIONS

In this section, we derive analytic solutions to the incompressible Navier-Stokes equations without source terms in a prescribed moving domain, which is initially a cylinder.

Using the cylindrical coordinates $(r, \theta, z)$ and the associated basis $(e_r, e_\theta, e_z)$ (where $e_3 = e_z$), the cylinder at initial time is given by

$$ \mathcal{C} = \{(r \cos \theta, r \sin \theta, z) | 0 \le r \le R, 0 \le \theta \le 2\pi, -L \le 2z \le L\}, $$

and we can write $\boldsymbol{u}$ with the form $\boldsymbol{u}(t, r, \theta, z) = u_r e_r + u_\theta e_\theta + u_z e_z$. Assuming that $\boldsymbol{u}$ does not depend on $\theta$ anymore and $u_\theta = 0$, $\boldsymbol{u}$ satisfies

$$ \boldsymbol{u}(t, r, z) = u_r e_r + u_z e_z. \quad (B.1) $$

The incompressible Navier-Stokes equations then become

$$ \partial_t u_r + u_r \partial_r u_r + u_z \partial_z u_r = -\partial_r p + \left[ \frac{1}{r} \partial_r (r \partial_r u_r) + \partial_{zz}^2 u_r - \frac{u_r}{r^2} \right], \quad (B.2) $$
$$ \partial_t u_z + u_r \partial_r u_z + u_z \partial_z u_z = -\partial_z p + \left[ \frac{1}{r} \partial_r (r \partial_r u_z) + \partial_{zz}^2 u_z \right], \quad (B.3) $$

$$ \frac{1}{r} \partial_r (r u_r) + \partial_z u_z = 0. \quad (B.4) $$

Those equations hold in the time-depending domain $C(t)$ obtained from $C(0) := C$ thanks to a transformation to be determined. Let $\kappa \in C^2(\mathbb{R}_+; \mathbb{R}_+^*)$ such as $\kappa(0) = 1$. Consider a point located at cylindrical coordinates $(r, \theta, z)$ at initial time. We look for the displacement of this point taking the form $(\kappa(t)r, \theta, \zeta_z(t))$ (cylindrical coordinates), where $\zeta_z$ is smooth enough and $\zeta_z(0) = z$. The latter function $\zeta_z, -L/2 \le z \le L/2$, must still be defined. With that kind of displacements, $C(t)$ is no longer a cylinder but keeps the symmetry axis $(Oz)$ and remains centred at the origin.

In order to simplify again the computations, we assume that $u$ is such that $u_r$ does not depend on $z$ and remains zero on the axis. Let us now write the boundary condition on the wall. The fluid and the wall share the same velocity. This ensures that

$$ u_r(t, \kappa(t)R) = \dot{\kappa}(t)R, \quad (B.5) $$

$$ u_z(t, \kappa(t)R, \zeta_z(t)) = \dot{\zeta}_z(t). \quad (B.6) $$

The simplest possibility for (B.5) to hold is then

$$ u_r(t, r) = \frac{\dot{\kappa}(t)}{\kappa(t)}r. \quad (B.7) $$

We still have to specify $u_z$, $p$, and $\zeta_z$ for each $z \in [-L/2, L/2]$. Equation (B.4) implies

$$ \partial_z u_z = -\frac{1}{r} \partial_r (r u_r) = -2 \frac{\dot{\kappa}}{\kappa}. $$

A convenient $u_z$ is then not depending on $r$, and has the form

$$ u_z(t, r, z) = -2 \frac{\dot{\kappa}(t)}{\kappa(t)} z + v(t), \quad (B.8) $$

where $v$ is a $C^1$ arbitrary function of $t$ and has the dimension of a velocity. Note that $v$ can also depend on $r$, but that dependence is not useful here. We can now use (B.7)–(B.8) in (B.2) and obtain an equation on $p$

$$ -\partial_r p = r \frac{\ddot{\kappa}}{\kappa}. $$

Consequently, $p$ has the form

$$ p(t, r, z) = -\frac{r^2}{2} \frac{\ddot{\kappa}(t)}{\kappa(t)} + q(t, z), $$

where $q$ depends on $t$ and $z$, and has the dimension of a pressure, up to the air mass density. Eventually, (B.3) allows to identify a convenient function $q$. Indeed, we can write

$$ \partial_z q(t, z) + \frac{2z}{\kappa(t)^2} [3\dot{\kappa}(t)^2 - \ddot{\kappa}(t)\kappa(t)] = -\dot{v}(t) + 2\frac{\dot{\kappa}(t)}{\kappa(t)}v(t). $$

Hence, up to a constant, $p$ satisfies

$$ p(t, r, z) = -\frac{r^2 \ddot{\kappa}(t)}{2 \kappa(t)} - \frac{z^2}{\kappa(t)^2} [3\dot{\kappa}(t)^2 - \ddot{\kappa}(t)\kappa(t)] + \left[ 2\frac{\dot{\kappa}(t)}{\kappa(t)}v(t) - \dot{v}(t) \right] z. \quad (B.9) $$

Eventually, (B.6) gives the following Cauchy problem on $\zeta_z$, which has a unique solution on $\mathbb{R}_+$,

$$ \dot{y} = -2 \frac{\dot{\kappa}}{\kappa} y + v, \quad y(0) = z. \quad (B.10) $$

Let us now sum up the result we just proved.
**Proposition 9.** Let $v \in C^1(\mathbb{R}_+; \mathbb{R})$, $\kappa \in C^2(\mathbb{R}_+; \mathbb{R}_+^*)$ with $\kappa(0) = 1$. For any $z \in [-L/2, L/2]$, consider the unique solution $\zeta_z$ to Cauchy problem (B.10). Set, for any $t \ge 0$,

$$C(t) = \{(\kappa(t)r \cos \theta, \kappa(t)r \sin \theta, \zeta_z(t)) \mid 0 \le r \le R, 0 \le \theta \le 2\pi, -L \le 2z \le L\}.$$

Then $u$ and $p$, defined by (B.1) and (B.7)-(B.9), solve the incompressible Navier-Stokes equations (B.2)-(B.4), with boundary conditions (B.5)-(B.6) on the wall of $C(t)$.

**Acknowledgements.** Alexander Lorz wants to thank both ANR Samovar project and FSMP which provided his former post-doctoral positions. The authors are also grateful to Laurent Desvillettes, with whom the discussions were very fruitful, and to the Reo project-team members who also worked on the LifeV finite element library, and particularly Matteo Astorino, Miguel A. Fernández, Jean-Frédéric Gerbeau and Michaël Grasseau.

## REFERENCES

[1] J. Agnew, D. Pavia, and S. Clarke. Aerosol particle impaction in the conducting airways. *Phys. Med. Biol.*, 29:767-777, 1984.

[2] A. A. Amsden. KIVA-3: a KIVA program with block-structured mesh for complex geometries. Technical Report LA-12503-MS, Los Alamos National Laboratory, 1993.

[3] A. A. Amsden. KIVA-3V: a block-structured KIVA program for engines with vertical or canted valves. Technical Report LA-13313-MS, Los Alamos National Laboratory, 1997.

[4] A. A. Amsden, P. J. O'Rourke, and T. D. Butler. KIVA-II: A computer program for chemically reactive flows with sprays. Technical Report LA-11560-MS, Los Alamos National Laboratory, 1989.

[5] A. A. Amsden, J. D. Ramshaw, P. J. O'Rourke, and J. K. Dukowicz. KIVA: A computer program for two- and three-dimensional fluid flows with chemical reactions and fuel spray. Technical Report LA-10245-MS, Los Alamos National Laboratory, 1985.

[6] O. Anoshchenko and A. Boutet de Monvel-Berthier. The existence of the global generalized solution of the system of equations describing suspension motion. *Math. Methods Appl. Sci.*, 20(6):495-519, 1997.

[7] B. Asgharian, W. Hofmann, and R. Bergmann. Particle deposition in a multiple-path model of the human lung. *Aerosol Sci. Tech.*, 34(4):332-339, 2001.

[8] L. Baffico, C. Grandmont, and B. Maury. Multiscale modeling of the respiratory tract. *Math. Models Methods Appl. Sci.*, 20(1):59-93, 2010.

[9] I. Balásházy and W. Hofmann. Deposition of aerosols in asymmetric airway bifurcations. *J. Aerosol Sci.*, 26(2):273-292, 1995.

[10] C. Baranger, L. Boudin, P.-E. Jabin, and S. Mancini. A modeling of biospray for the upper airways. *ESAIM Proc.*, 14:41-47, 2005.

[11] L. Boudin, L. Desvillettes, C. Grandmont, and A. Moussa. Global existence of solutions for the coupled Vlasov and Navier-Stokes equations. *Differential Integral Equations*, 22(11-12):1247-1271, 2009.

[12] L. Boudin, L. Desvillettes, and R. Motte. A modeling of compressible droplets in a fluid. *Commun. Math. Sci.*, 1(4):657-669, 2003.

[13] L. Boudin, C. Grandmont, B. Grec, and D. Yakoubi. Influence of the spray retroaction on the airflow. In *CEMRACS 2009 – Mathematical Modelling in Medicine*, volume 30 of *ESAIM Proc.*, pages 153-165. EDP Sci., Les Ulis, 2010.

[14] R. Clift, J. R. Graces, and M. E. Weber. *Bubbles, drops and particles*. Academic Press, 1978.

[15] J. K. Comer, C. Kleinstreuer, and Z. Zhang. Flow structures and particle deposition patterns in double bifurcation airway models. part 1. Air flow fields. *J. Fluid Mech.*, 435:25-54, 2001.

[16] J. K. Comer, C. Kleinstreuer, and Z. Zhang. Flow structures and particle deposition patterns in double bifurcation airway models. part 2. Aerosol transport and deposition. *J. Fluid Mech.*, 435:55-80, 2001.

[17] A. Comerford, G. Bauer, and W. A. Wall. Nanoparticle transport in a realistic model of the tracheobronchial region. *Int. J. Numer. Method Biomed. Eng.*, 26(7):904-914, 2010.

[18] G.-H. Cottet and P.-A. Raviart. On particle-in-cell methods for the Vlasov-Poisson equations. *Transport Theory Statist. Phys.*, 15(1-2):1-31, 1986.

[19] A. Decoène. *Hydrostatic model for three-dimensional free surface flows and numerical schemes*. PhD thesis, Université Pierre et Marie Curie, France, 2007.

[20] P. Degond and S. Mas-Gallic. The weighted particle method for convection-diffusion equations. I. The case of an isotropic viscosity. *Math. Comp.*, 53(188):485-507, 1989.

[21] P. Degond and S. Mas-Gallic. The weighted particle method for convection-diffusion equations. II. The anisotropic case. *Math. Comp.*, 53(188):509-525, 1989.

[22] L. Desvillettes. Some aspects of the modeling at different scales of multiphase flows. *Comput. Methods Appl. Mech. Engrg.*, 199(21-22):1265-1267, 2010.
[23] J. Donéa. Arbitrary Lagrangian Eulerian methods. In *Computational Methods for Transient Analysis*, volume 1 of *Computational Methods in Mechanics*, pages 473–516. North Holland, Elsevier, 1983.

[24] G. Dufour. Modélisation multi-fluide eulérienne pour les écoulements diphasiques à inclusions dispersées. PhD thesis, Université Paul-Sabatier Toulouse-III, France, 2005.

[25] J. Fouchet-Incaux. Artificial boundaries and formulations for the incompressible Navier-Stokes equations: applications to air and blood flows. *SéMA J.*, pages 1–40, 2014.

[26] P. J. Frey and P.-L. George. *Mesh generation*. Hermes Science Publishing, Oxford, 2000. Application to finite elements, Translated from the 1999 French original by the authors.

[27] T. Gemci, T. E. Corcoran, and N. Chigier. A numerical and experimental study of spray dynamics in a simple throat model. *Aerosol Sci. Technol.*, 36:18–38, 2002.

[28] C. Grandmont, Y. Maday, and B. Maury. A multiscale/multimodel approach of the respiration tree. In *New trends in continuum mechanics*, volume 3 of *Theta Ser. Adv. Math.*, pages 147–157. Theta, Bucharest, 2005.

[29] E. Guyon, J.-P. Hulin, and L. Weber. *Hydrodynamique physique*. CNRS Éditions, 2001.

[30] K. Hamdache. Global existence and large time behaviour of solutions for the Vlasov-Stokes equations. *Japan J. Indust. Appl. Math.*, 15(1):51–74, 1998.

[31] W. Hofmann. Modelling inhaled particle deposition in the human lung – a review. *J. Aerosol Sci.*, 42(10):693–724, 2011.

[32] T. J. R. Hughes, W. K. Liu, and T. K. Zimmermann. Lagrangian-Eulerian finite element formulation for incompressible viscous flows. *Comput. Methods Appl. Mech. Engrg.*, 29(3):329–349, 1981.

[33] P. W. Longest and M. Hindle. Numerical model to characterize the size increase of combination drug and hygroscopic excipient nanoparticle aerosols. *Aerosol Sci. Tech.*, 45(7):884–899, 2011.

[34] T. B. Martonen. Mathematical model for the selective deposition of inhaled pharmaceuticals. *J. Pharm. Sci.*, 82:1191–1199, 1993.

[35] J. Mathiaud. *Étude de systèmes de type gaz-particules*. PhD thesis, École normale supérieure de Cachan, France, 2006.

[36] C. Mitsakou, C. Helmis, and C. Housiadas. Eulerian modelling of lung deposition with sectional representation of aerosol dynamics. *J. Aerosol Sci.*, 36:74–95, 2005.

[37] A. Moussa. *Mathematical and numerical study of the transport of sprays in the human lung*. PhD thesis, École Normale Supérieure de Cachan, France, 2009.

[38] F. Nobile. *Numerical approximation of fluid-structure interaction problems with applications to haemodynamics*. PhD thesis, École Polytechnique Fédérale de Lausanne, Switzerland, 2001.

[39] J. M. Oakes, A. L. Marsden, C. Grandmont, S. C. Shadden, C. Darquenne, and I. E. Vignon-Clementel. Airflow and particle deposition simulations in health and emphysema: From in vivo to in silico animal experiments. *Annals of biomedical engineering*, 42(4):899–914, 2014.

[40] P. J. O'Rourke. *Collective drop effects on vaporizing liquid sprays*. PhD thesis, Los Alamos National Laboratory & Princeton University, USA, 1981.

[41] D. J. Torres and M. F. Trujillo. KIVA-4 An unstructured ALE code for compressible gas flow with sprays. *J. Comput. Phys.*, 219:943–975, December 2006.

[42] T. Westermann. Particle-in-cell simulations with moving boundaries — adaptive mesh generation. *J. Comput. Phys.*, 114(2):161–175, 1994.

[43] F. Williams. *Combustion Theory*. Benjamin Cummings, second edition, 1985.

[44] C. Yu. Uniqueness of global strong solutions to the 2D Navier-Stokes-Vlasov equations. *ArXiv e-prints*, July 2012.

[45] C. Yu. Global weak solutions to the incompressible Navier-Stokes-Vlasov equations. *J. Math. Pures Appl. (9)*, 100(2):275–293, 2013.

[46] Z. Zhang, C. Kleinstreuer, and C. S. Kim. Airflow and nanoparticle deposition in a 16-generation tracheobronchial airway model. *Ann. Biomed. Eng.*, 36(12):2095–2110, 2001.

L.B., A.M.: SORBONNE UNIVERSITÉS, UPMC UNIV PARIS 06 & CNRS, UMR 7598 LJLL, PARIS, F-75005, FRANCE & INRIA, ÉQUIPE-PROJET REO, BP 105, F-78153 LE CHESNAY CEDEX, FRANCE
*E-mail address:* laurent.boudin@upmc.fr, ayman.moussa@upmc.fr

C.G.: INRIA, ÉQUIPE-PROJET REO, F-78153 LE CHESNAY CEDEX, FRANCE & SORBONNE UNIVERSITÉS, UPMC UNIV PARIS 06 & CNRS, UMR 7598 LJLL, PARIS, F-75005, FRANCE
*E-mail address:* celine.grandmont@inria.fr

A.L.: SORBONNE UNIVERSITÉS, UPMC UNIV PARIS 06 & CNRS, UMR 7598 LJLL, PARIS, F-75005, FRANCE & INRIA, ÉQUIPE-PROJET MAMBA, BP 105, F-78153 LE CHESNAY CEDEX, FRANCE
*E-mail address*, Corresponding author: alexander.lorz@upmc.fr