Fluid-kinetic modelling for respiratory aerosols with variable size and temperature

. In this paper, we propose a coupled ﬂuid-kinetic model taking into account the radius growth of aerosol particles due to humidity in the respiratory system. We aim to numerically investigate the impact of hygroscopic eﬀects on the particle behaviour. The air ﬂow is described by the incompressible Navier-Stokes equations, and the aerosol by a Vlasov-type equation involving the air humidity and temperature, both quantities satisfying a convection-diﬀusion equation with a source term. Conservations properties are checked and an explicit time-marching scheme is proposed. Two-dimensional numerical simulations in a branched structure show the inﬂuence of the particle size variations on the aerosol dynamics. Résumé


Introduction
Aerosol therapy is one of the major curative way to treat chronic obstructive pulmonary diseases (COPD). Since in vivo observations of drug delivery in the human airways induce many difficulties, it appears crucial to model and to be able to numerically simulate the aerosol flow in the lung with a good enough accuracy, the main question being linked to the deposition phenomenon, and particularly its location.
The aerosol constitutes a set of very numerous particles. Those particles happen to have hygroscopic properties, i.e. abilities to exchange water with the bronchial air, which is full of humidity. Consequently, the aerosol particle sizes may vary in time, as it was assessed in [9][10][11]. Assuming that the droplets remain spherical, we intended to study the influence of radius growth on deposition (number of deposited particles, characteristic times of propagation/deposition inside a given realistic geometry, deposition areas...). It is worth noticing that those hygroscopic properties strongly rely on thermal effects, and the model we investigate involves both air and aerosol temperatures.
There are several kinds of models allowing to describe the aerosol motion in the air, two-phase models, agent-based ones and kinetic ones. In the first situation, the aerosol droplets are considered as a fluid which constitues a mixture with the ambient air in the lung. Then one focuses on the aerosol concentration in the air, see [1,5] for instance. Nevertheless, with those models, it may be difficult to determine the aerosol deposition areas. The second situation, which is used in [13,16] for instance, raises many difficulties to track the trajectory of numerous particles, in particular when there is a strong coupling between the particles and the air flow.
We here choose to use a kinetic approach, as in [2,7], which is relevant from the modelling point of view: there are lots of particles in the aerosol, but their volume is negligible compared to the airways volume. In this framework, the aerosol behaviour is described through a distribution function which depends on macroscopic variables (time, space position) as well as on microscopic ones (velocity, for instance). The Vlasov-type equation satisfied by the distribution function is usually coupled with the incompressible Navier-Stokes equations to describe the airflow [4]. The model we present here is an extension of the latter one, where the air temperature, the mass fraction of the water vapour in the air and the dependence of the distribution function on both the size and temperature of the particles are taken into account. Note that it is also investigated from the mathematical point of view in [12].
This article is divided into four main parts. The first one aims to present the model(s) we study. In the brief second one, we formally check some relevant balance properties of our model. Then we discuss the numerical method we use to discretize our model in the third section. Eventually, the last section is dedicated to numerical experiments which point out the aerosol size growth importance in the deposition phenomenon.

Building the model
As stated above, this section aims to explain and present a model describing the behaviour of a therapeutic aerosol in the respiratory system, where the size of the aerosol particles can vary because of the ambient water vapour in the airflow.
We use a simplified fluid domain, denoted by Ω: in our setting, Ω can be chosen as a cylinder or a branch, and is assumed not to depend on time. Its boundary Γ = ∂Ω is divided into three subsets, the wall Γ wall , the inlet Γ in and the outlet Γ out , see Figure 1.
As we already stated, we use a kinetic viewpoint to model the aerosol, which constitutes a dispersed phase. Let us consider a distribution function f representing the density of particles per unit volume in the domain of interest Ω ⊂ R 3 . It depends on time t ≥ 0, position x ∈ Ω, velocity v ∈ R 3 , radius r > 0 and temperature T > 0. The usual dependence of f with respect to the variables t, x, and v is supplemented by the one with respect to the size and temperature of the droplets. Indeed, the condensation phenomenon which induces the size variation involves matter and thermal exchanges between the aerosol and the air. Note that we assume that the particles remain spherical and do not interact with each other. The latter hypothesis relies on the fact that there are not enough particles for their collisions to be significant. In the respiratory system, the aerosol moves in the air, which can be assumed to be Newtonian and incompressible. In particular, this means that the air mass density air is constant. The air flow is then classically described through the air velocity u and its pressure p.
In order to consider the growth of particles and describe the matter and thermal exchanges with the air, we also need to work with the water vapour mass fraction in the air Y v,air and the air temperature T air . All those quantities indexed by "air" depend on t and x (except air , of course). Note that the water vapour mass density in the air considered as a gaseous mixture can be computed as air Y v,air .
Let us now focus on the various equations satisfied by the aerosol and air-related quantities, starting with the dispersed phase. The aerosol density f satisfies a Vlasov-type equation, which can be written as: where g is the gravitational field, α(u − v) the drag acceleration undergone by the aerosol from the air, and a and b respectively represent the radius and temperature growth evolutions of the particles. The Vlasov equation is supplemented with the following boundary and initial conditions where f in : given. The boundary condition on Γ wall is standard, it ensures that any droplet landing on the wall remains deposited on the wall afterwards.
We still need to define the functions α, a and b in (1). In order to do so, it is important to first understand how a particle behaves and what kind of phenomenon we deal with.
The main new phenomena taken into account in the model are the water vapour condensation on the droplet surface and the water evaporation from the droplet surface. In each particle, one can find active products (drug), excipient and water. Let us introduce some equivalent radii to simplify the situation. Denote by r drug the radius such that 4 3 πr 3 drug drug is the drug mass inside the droplet and by r ex the one such that 4 3 π(r 3 ex − r 3 drug ) ex is the excipient mass inside the droplet (see Figure 2). Obviously, drug and ex respectively are the (constant) drug and excipient mass densities. The latter equivalent radius r ex will also be named the particle dry radius, in the sense that there is no water in a dry particle.
Denoting by w the water mass density, assumed to be constant in the situations we investigate, it is now possible to define the mass and mass density of the particle, which both depend on r, i.e. This allows to properly define α. Indeed, we have, as in [4], where η is the constant air dynamic viscosity, so that the drag force satisfies the Stokes law. We rely on [10] to build the functions a and b. Since we want to take into account condensation and evaporation, the key physical quantities are the heat fluxes between the air and the droplets, the convective one Q d and the evaporating one L v N d , where L v is the latent heat of water vaporisation. The convective flux depends on r, T , but also on t and x through T air . The evaporating heat flux involves the water mass flux N d at the droplet surface. This flux depends on r, T , and also on t and x again, but this time through Y v,air , and also governs the size evolution of the droplets. The expressions of the functions a and b can then be written in terms of N d and Q d as where c P d is the constant specific heat of the droplet. Let us now detail the expressions of Q d and N d . We can first write where Nu is the droplet Nusselt number corresponding to the ratio of convective to conductive heat transfer between the droplet and the air, κ air the thermal conductivity of the air as a gaseous mixture, and C T the Knudsen correlation which is small in our case allowing to neglect non-continuum effects of the fluid, all of them being considered as constants in our case. Second, the water mass flux N d involves Y v,air and also the mass fraction Y v,surf of water vapour at the droplet surface, which depends on r and T . That quantity Y v,surf is quite intricate to be defined from the modelling viewpoint. We do not provide any detailed explanation on the model itself, the reader is invited to refer to [10] for more details. We first need the water vapour saturation pressure, which depends on T , and is here expressed in the cgs unit system: Then the influence of the Kelvin effect on the droplet surface concentration of water vapor is expressed thanks to where R v is the gas constant of water vapour and the droplet surface tension σ is assumed to be constant too in our framework (in particular not depending on T ). Eventually, we need the water activity coefficient given by where M w , M drug and M ex respectively denote the molar masses of water, drug and excipient, and the constants i drug and i ex are the so-called van't Hoff factors of the drug and the excipient, allowing to take into account the molecular dissociation during dissolution. Note that the previous expression S is different from the one in [10]. Indeed, S equals 0 when the particle is dry, i.e. when r = r ex , which was not possible in [10], but otherwise, of course, we recover the same value for S. We can now write the expression of Y v,surf , that is The expression of the water mass flux subsequently comes where Sh is the Sherwood number describing the water transfer between the air and the droplet, C m is the mass Knudsen number correction, and the binary diffusion coefficient D v of water vapour in the air is given, in the cgs unit system, by The definitions (6)- (8) and (13) of a, b, Q d and N d allow to determine lower bounds for the support of the distribution function with respect to both r and T . We know that a governs the time evolution of this support. Assume that, at initial time, the droplets all have a radius greater than r ex . In that situation, if r somehow reaches the value r ex , as we already pointed out, we have S(r ex ) = 0 and subsequently Y v,surf = 0, which implies N d ≤ 0 and a ≥ 0. That ensures that r cannot go below r ex in the model. For the support of f with respect to T , from the mathematical viewpoint, even if the formulae leading to (13) do not hold in the considered temperature range, one can check that if T somehow reaches 46.13 K (the value in P v,sat ), N d is non positive, since Y v,surf = 0 again. And Q d is clearly non positive too. That ensures that T cannot go below 46.13 K in the model. Anyway, if all the functions involved are smooth enough, if f init is chosen with a compact support in all its variables, that compactness property should hold at least for small times. Note that, from the physical viewpoint, it seems relevant to assume that all the temperatures remain around 300 K, and that Y v,air remains positive.
Let us now tackle the air equations. To begin with, the fluid velocity u(t, x) ∈ R 3 and its pressure p(t, x) ∈ R classically (see [4] for instance) satisfy the incompressible Navier-Stokes equations on R + × Ω, completed with the following boundary and initial conditions: where σ(u, p) = ∇ x u+(∇ x u) −p Id is the stress tensor, n is the outgoing normal vector from Γ, u in : R + ×Γ in → R 3 is chosen to be a Poiseuille flow and u init : Ω → R 3 is the initial datum. Eventually, the aerosol retroaction F on the air is given by In our numerical simulations, that latter term will be neglected (F = 0) because of its smallness related to the investigated particle size range, see [4], but we keep it for the time being, in order to remain consistent with respect to the total momentum conservation of our system.
We still need to write equations to describe the evolution of the air temperature T air and the water vapour mass fraction Y v,air in the air. The water vapor mass fraction Y v,air satisfies an advection-diffusion equation including a source term S Y which accounts for the water mass exchanges between the bronchial air and the aerosol. Other effects such as turbulence effects could also be taken into account (see [9,11]). Thus, we write, on R + × Ω, completed with the following boundary and initial conditions where Y v,air,init , Y in v,air , Y v,wall > 0 are given. The boundary condition on Γ wall ensures that the wall continuously provides water vapour to the air mixture. The source side term S Y of (18) is defined, very similarly to the one from [10], as Finally, the air temperature also satisfies an advection-diffusion equation on R + × Ω, which writes completed with the following boundary and initial conditions where T air,init , T in air , T wall > 0 are given. The source side term S T of (21) represents the heat transfer between the air and the aerosol through the water vapour, and is given, again very similarly to the one from [10], by In what follows, equations (1)-(23), which include all the effects related to aerosol size and temperature variation, will be referred to as the (A) model. If we drop the temperature effects, the (B) model is given by (1)- (6) and (9)-(20), with b = 0. And for the (C) model, where no aerosol size or temperature variation is allowed, we consider equations (1)- (5) and (15)-(17) with a = 0 and b = 0. To summarise the models under study, we rewrite below the PDEs for each one of them, without the boundary and initial conditions: varying size and temperature, varying size, no temperature variation, no size nor temperature variation.

Checking the physical conservations
In this section, we discuss the physical conservations of various quantities related to the (A) system, and we focus on two quantities which involve water vapour. The first one focuses on the water vapour mass exchange. More precisely, the water vapour coming from the air is supposed to lead to a size variation of the aerosol droplets. That property is stated in the following proposition.
Proof. On the one hand, multiplying (1) by m(r), integrating with respect to all the variables except t, and eliminating the conservative terms through integrations by parts, we obtain d dt Ω x)) f (t, x, v, r, T ) dv dr dT dx.
The exchanges of thermal energy associated to water transfers between the air and the aerosol are more intricate to understand. The following proposition holds.
Proposition 2. Assume that u = 0 and ∇ x T air · n = 0 on ∂Ω, and that Proof. On the one hand, we integrate (21) over Ω to obtain Then we multiply (1) by m(r)c P d T and integrate it with respect to all the variables except t to get d dt Ω Then we sum both previous equalities using (3), (4), (6) and (7) to recover (24), noticing that the term involving Q d vanishes, keeping two terms involving N d : one with L v to take the change of physical state into account and one with the added thermal energy in the aerosol due to the mass exchange.
For the total momentum conservation, when the retroaction term is involved, in fact, the formal proof is exactly the same as in [4]. We do not present any more details about the conservation of this quantity and the total kinetic energy decreasing.

Numerical method
Before studying some relevant physical situations from the computational viewpoint, let us describe, in this section, the numerical scheme we implemented to solve system (A).
We proceed in the same way as in [4], by using a time-marching scheme and uncoupling the fluid and aerosol equations. First, we solve the three fluid equations (15)-(16), (18) and (21) with the source terms computed at the previous time step for both equations for Y v,air and T air (recall that F is chosen equal to 0). Then we solve the Vlasov equation (1) using the updated fluid quantities. For the fluid equations, we use a finite-element method, and the aerosol is computed thanks to a particle method. All the computations are performed in a two-dimensional setting with FreeFem++ [8], on a time interval [0, τ ], where τ > 0 is given. For a time step ∆t > 0 such that τ /∆t ∈ N * , we denote t n = n∆t for any n, with 0 ≤ n ≤ N . The domain Ω is discretized as a tetrahedric mesh Ω h .
Let us provide more details about the numerical solving. As already stated, we start by solving the fluid equations, and use a finite element method to do so. If we write weak formulations of those equations, we introduce the following test functions: χ ∈ L 2 (Ω) for (16), and ν, ψ, φ ∈ H 1 (Ω), vanishing on Γ in and Γ wall , respectively for (15), (18) and (21). Then, for the discretization, we are led to use P 2 functions for the velocities u and ν and P 1 functions for p, Y v,air , T air , χ, ψ and φ.
To go from t n to t n+1 , assume that u n , Y n v,air , T n air , S n T and S n Y are known. In order to handle the convective term in (15), we introduce the approximated characteristic flow X n , which approximates the solution X to the Cauchy problem, set on [t n , t n+1 ],Ẋ (s) = u n (s, X(s)), for any x ∈ Ω h . The approximation X n is computed thanks to the FreeFem++ command convect. Hence, we define u n+1 as the solution to the following discrete weak formulation We proceed in the same way to define Y n+1 v,air and T n+1 air as the solutions to the discrete weak formulations Note that, in the considered air temperature range, the value of D v , which theoretically depends on T air , only has a 2 % variation, so we choose not to update D v with respect to T air in this equation. Next comes the aerosol numerical solving. Since we use a particle method, we discretize the distribution function f as a weighted sum of Dirac masses in the position, velocity, radius and temperature variables. More precisely, we consider N num ∈ N * numerical particles, each of them having the representativity ω ∈ N * , so that the total number of physical particles is N aero = ωN num . Note that, for obvious computational cost reasons, N num must be chosen to be small with respect to N aero , but large enough to faithfully represent the distribution of the aerosol particles.
The distribution function is then approximated by where x p (t), v p (t), r p (t), T p (t) are the coordinates, in the space phase of f , of the numerical particle p ∈ {1, . . . , N num } at time t.
The particle coordinates solve the following differential system supplemented with initial data chosen to fit the distribution of droplets given by f init . The initial data chosen in our case are detailed in the next section.
The ODE on r p is solved thanks to a RK4 scheme, in order to remain very accurate for the radii computations, involving the newly computed value Y n+1 v,air at the current position of the particle x n p . The velocity and temperature ODEs are solved with a semi-implicit Euler scheme, respectively involving u n+1 , and Y n+1 v,air , T n+1 air , all of them again being valued at x n p . Eventually, the position x n+1 p is updated using the newly computed v n+1 p . Once we know (x n+1 p , v n+1 p , r n+1 p , T n+1 p ) for all p = 1, . . . , N num , we define the source terms for the next time step t n+1 as In order to take into account the deposition or exiting conditions from (2), in the previous definition of the discrete source terms at time t n+1 , we drop the indices p such that x n+1 p is outside the domain Ω h in the sum p ∈ {1, . . . , N num }. It means that, at time t n+1 , the corresponding particle p either deposited on the wall or went out of the domain through the outlet (depending on its distance to Γ wall and Γ out at time t n ) . Note that a particle p is also considered to be deposited if the distance between x p and Γ wall is smaller that r p . The same goes for Γ out for exiting the domain. Once the particle is deposited or out, it is of course not treated numerically anymore.
It appears mandatory to perform a time-subcycling for the aerosol computations. Without that subcycling, the particle would be able to go across various cells during the same (fluid) time step. Moreover, it is particularly relevant when dealing with the particle temperatures because, if we compute the coefficients involved, we can check that the temperature ODE is very stiff.

Numerical simulations
In this section, we numerically investigate various situations. Let us first provide the specific setting of our numerical experiments. Note that the values of the physical constants appearing in our models are gathered in Appendix A.

Experimental context, reference situation
The final time is set at τ = 1 s. Our domain is supposed to represent the trachea and the first bifurcation in the human airways, see ( Figure 1). Its characteristic sizes and shape are the ones described in [14,15], taking into account a 3D-2D correction coefficient for each branch length. More precisely, the diameter of the trachea is set at D 0 = 1.80 cm, and its length at 0 = 7.52 cm. The right-hand bronchus has an angle of 25°with the tracheal axis, it is quite short ( 10 = 3.75 cm), but its diameter quite large, D 10 = 1.50 cm. The left-hand bronchus has an angle of 45°with the tracheal axis, it is longer than the first one ( 01 = 6.75 cm), but its diameter is smaller, with D 01 = 1.00 cm. Let us emphasize that the right-hand bronchus is the left branch on Figure 1, and conversely: we have the outsider's view, not the patient's.
The middle of the boundary inlet Γ in is the origin of the space coordinate system, with horizontal and vertical axes.
Let us now provide the boundary and initial conditions for the fluid equations. The velocity fluid is initialized at u init = 0, and, at the inlet, u in follows a Poiseuille law, i.e. it is vertically oriented from up to bottom and its modulus is given, for any x ∈ Γ in , by where u 0 = 50.0 cm.s −1 . The air temperature uses the following values T air,init = 37°C = 310 K, T in air = 24°C = 297 K, T wall = 37°C = 310 K.
The initial and boundary values of Y v,air uses the ones of the relative humidities in the airways, here chosen as RH lung = 0.990 and RH wall = 1.00. Then we write Figure 3 shows the values of |u|, Y v,air and T air at final time, where some kind of stationary state is reached for the fluid. Let us next describe the computational and experimental parameters for the aerosol. We consider 5 injections in the domain of 100 numerical particles each with representativity ω = 10 4 , periodically released between initial time and t = 0.25 s. Hence, we deal with N num = 500 numerical particles and N aero = 5 10 6 physical particles. All the numerical particles initially have the same vertical velocity v p,2 (0) = −100 cm.s −1 , the same radius r p (0) = 2.25 10 −5 cm, and the same temperature T p (0), equal to the air temperature T in air at the inlet. They are released from random positions x p (0) ∈ Γ in with its first coordinate in [−D 0 /4, D 0 /4], following a uniform law. We choose this latter interval instead of [−D 0 /2, D 0 /2] so that it allows a larger deposition phenomenon. Since we use a particle method, it is mandatory, in order to obtain meaningful results, to perform averaging computations over several initial randomly chosen distributions of droplets, in our case, 10 different distributions. Note that the reader can refer to [6], for instance, for more information related to numerical methods solving the Vlasov equation, including the particle methods.  This parameter set defines what we name the reference situation, which was first used to validate the code. Indeed, we checked, at the computational level, the total mass conservation in the water vapour exchange between air and aerosol and the thermal energy balance implied by the thermodynamic state change of water vapour. The next subsection also lies in the framework of the reference situation, in order to give an almost exhaustive overview of all the phenomena.

Numerical tests
We first showed, for the (A) model, the behaviour of the air velocity u and temperature T air at different times, as well as the movement in the domain of the various aerosol releases, in Figures 4-5. We do not provide snapshots for the water vapor mass fraction Y v,air because they would look like Figure 3 Figure 6. Local effects of the aerosol on the air temperature, here at time 0.25 + ∆t (in seconds). Figure 6, where the air temperature is shown just after the last aerosol release at time 0.25 s, plotting or not the particle positions, allows to point out an effect which could not be seen in the snapshots from Figure 5, because of the plotting of the particles in the domain. Indeed, for all aerosol releases except for the first one, there is a local air temperature increasing at the location of the particles. This effect is clearly not explained by direct thermal phenomena, hence it comes from the water vapour mass exchange between the humidified air and the droplets.
On Figure 7, we represent the particle trajectories obtained with the (A) model. In this example, most of the droplets (348 over 500) go out of the domain through the left branch. A less significant number of particles (98) go out through the right branch, and a smaller part of them (47) deposit on the wall (more precisely, on the "internal" wall of the left branch). There are still 7 droplets in the domain at final time. The major exit at the left branch is of course no surprise: because of its diameter, the branch is the most natural way out for the aerosol.
Let us now focus on Figure 8. The plots show the time evolution of (a) the radius and (b) the temperature of all the numerical particles from one of the initial distributions we used, in the (A) model. It is clear from Figure 8(a) that, at first, the droplets from the first release do not behave in the same way as the ones from the other releases: the size growth happens more slowly. The difference is confirmed by Figure 8(b): the temperature plots for the first release again behave very differently from the others. In particular, even if the temperatures of the first injected particles are initially 297 K, they almost instantaneously become close to 310 K because they move in an air at 310 K, see Figure 5(a). On the contrary, the other releases live in a cooler air, see Figure 5(b)-(e), and hence are not submitted to a thermal shock. Then, we can check on Figure 8 that each plot seems to have a characteristic time of behaviour change. It is definitively clear on Figure 8(b). For all releases except for the first one, there is a temperature jump which is, for the second one, located in time at approximately 0.25 s. If we observe the temperature snapshot 5(d), we can see that the corresponding time, which is 0.24 s, is approximately the one when the second release goes into the branches, where the diameters are significantly smaller than in the trachea, and where the wall effect consequently appears stronger. This explains the fact that the particle temperatures increase more. Then, on Figures 9-11, we compare the radius and temperature evolution, with respect to time, of particular droplets of the second release evolving in models (A), (B) and (C). Figure 9 focuses on a droplet exiting the domain through the left branch with the (A) model, Figure 10 through the right branch, and Figure 11 on a droplet depositing. Note that those same particles share the same future in models (A) and (C), but are all deposited in the (B) model.  Figure 11. Radius and temperature evolution of a particular droplet which deposits, comparison in the three models.
Let us first comment Figures 9(b)-11(b), i.e. the temperature plots. Of course, in the (B) and (C) models, the temperature remains constant, whereas, in the (A) model, the particle temperature grows until it (approximately) reaches T wall . This may seem peculiar. Indeed, the aerosol enters the domain at 297 K, and, except for the first aerosol release, it evolves in a cooled air. Since the droplet temperature variations cannot be explained because of the ambient air temperature, it means that they are triggered by hygroscopic phenomena. This leads us to study more carefully Figures 9(a)-11(a). As one can see, in the three situations, model (B) induces a larger size growth than model (A): this may explain the fact that the particle deposits in model (B). The hygroscopic effects imply radius variations for both models, but a part of this variation existing in the (B) model also has a temperature effect in the (A) model, which justifies that the radius for (A) is smaller than the one for (B).
Let us finally provide some statistical behaviours of the aerosol for each model. In Table 1, we first give the following values, all of them being computed from the 10 initial aerosol distributions used to solve the Vlasov equation: • the droplets mean radius and temperature at "final" time (in the sense that, when a particle exits or deposits during the computation, its radius and temperature do not change anymore until the end of the simulation), • the mean percentage of deposited particles and of particles reaching the boundaries Γ out (right and left) among the 500 particles, • the corresponding mean event time, i.e. the time when the particles exit or deposit. As seen in Figures 9(a)-11(a), the particle mean radius in the (A) model is intermediate between the ones from (C) and (B), and, as already emphasized, when we neglect the temperature effects, the particle radius must bear the whole hygroscopic effect in the (B) model. Besides, the size growth between (C) and (A) is significant. It seems to be the main reason for aerosol deposition in that case. Otherwise, (A) and (C) models have closer mean behaviours, which may imply that the (B) model is not relevant.
Finally, Table 2 gives the droplets mean radius at final time depending on their evolution (deposition, left/right exiting) for (A) and (B) models, since the radius remains constant to 2.25 10 −5 cm with the (C) model. There exists a significant variation of the radius, which can be linked to the mean event times in Table 1. In the (A) model, the particles going out through the right branch remain longer in the domain, thus undergoing a larger radius growth. For the (B) model, deposition or exiting happen more or less at the same time, leading to very similar radii for the particles.  Table 2. Statistics for the particles depending on their future (depositing/exiting).

Conclusion
The previous section allowed to point out the relevance of model (A) compared to (B) and (C) to properly take into account the hygroscopic effects on aerosols in the airways.
Besides, in our opinion, there are plenty of situations remaining to investigate, which should be addressed in further works. Let us mention some of them below.
First, observe that the source term modelled by function b given by (7), which drives the behaviour of the aerosol particles, is naturally composed of two terms, which we may denote by b 1 and b 2 (the first one involves Q d and the second one N d ). In fact, b 1 and b 2 happen to have opposite signs and the same order of magnitude, around 2 10 5 K/s at initial time, whereas b approximately equals 400 K/s (see Figure 12). Even if the model seems to behave nicely with respect to temperatures (the temperature of the corresponding numerical particle is given on Figure 13), it is only because we used a very fine subcycling time step to guarantee numerical accuracy in the description of the thermal effects. It allows us to have numerical preservation of mass conservation and energy balance. However, from the computational viewpoint, this can probably be improved.  Figure 13. Temperature evolution of the chosen particle in Figure 12.
We did not provide any numerical tests with the presence of excipient in the aerosol, i.e. we chose r ex = r drug , though we wrote the (A) model involving the excipient. Nevertheless, since standard values of drug and ex are close to each other, it is unlikely that the excipient implies major behaviour changes on the aerosol.
Moreover, the air flow we here studied was only related to the inspiration part of the respiration. However, it seems difficult to tackle the expiration part, because boundary conditions on f at Γ out are then unclear.
The computational domain here represents the upper part of the airways, including the trachea. It would be relevant to study the model behaviour within other domains, not necessarily with a vertical main axis, to understand the effect of the geometrical variability.
Eventually, as it was emphasized in [4] where 3D computations are performed, compared to [3] in a 2D domain, two-dimensional simulations tend to increase the aerosol deposition phenomenon. In order to study the model very faithfully, it is probably mandatory to work in a 3D setting.  Table 3. Value of the physical constants.
1 Corporation designing, developing and manufacturing drug delivery devices