## 2 Methodology
### Flow and Particle Transport Simulation  
The airflow field was developed by a CFD code using the ANSYS Workbench finite element algorithm. In this study, we used the laminar, incompressible and isothermal Navier-Stokes equations with air specifications at 37 \({}^{\circ}\)C for modelling the airflow conditions in the alveolar region. Because of the low Reynolds number and laminar flow regime, the velocity profile at the inlet of the model was assumed to be parabolic, and a no-slip boundary condition was applied at the wall boundaries. The volumetric flow rate \(\dot{Q}_{G}\) for various generations was scaled from the tracheal flow rate \(\dot{Q}_{T}=500\) mL/s as follows:  
\[\dot{Q}_{G}=\frac{\dot{Q}_{T}}{2^{G}} \tag{10}\]  
where \(G\) is the generation number. The application of Lagrangian tracking in the modelling of particle transport involves the integration of individual particles that are tracked from their injection point until they escape the domain. In addition, just as the fluid affects particle transport through forces, there is a counteracting influence of the particle on the fluid, i.e., there is a coupling between the fluid and the particles. One-way coupling occurs if the fluid is permitted to affect particle dynamics but the particles do not influence the fluid pathways, while two-way coupling exists if the particles influence the fluid behaviour reciprocally. The particle-tracking model implemented in this study, which is referred to as the one- and two-way particle-fluid discrete element method (DEM), uses the Euler-Lagrange approach. Therefore, the mutual effect of airflow and particles is calculated along with the mutual interactions between spherical particles by means of binary collisions [13]. The Lagrangian approach is used for tracking the particles that are computed with the post-processor of the CFD package. In order to compute the particle trajectories, we use the Langevin equation, which balances the mass acceleration of the particle with the forces acting on it if the particulate density is large in comparison to the fluid density:  
\[m\frac{\mathrm{d}\mathbf{u}_{\mathrm{P}}}{\mathrm{d}t}=\mathbf{F}_{\mathrm{D} }+m\mathbf{g}+\mathbf{F}_{\mathrm{Br}}+\mathbf{F}_{\mathrm{VM}}+\mathbf{F}_{ \mathrm{C},ij} \tag{11}\]  
where \(\mathbf{u}_{\mathrm{P}}\), \(\mathbf{F}_{\mathrm{D}}\), \(\mathbf{F}_{\mathrm{Br}}\), \(m\mathbf{g}\) and \(\mathbf{F}_{\mathrm{VM}}\) denote velocity of the particles, drag force, force resulting from Brownian motion, gravitational force and virtual (or added) force, respectively. Moreover, the inter-particle contact forces resulting from the DEM approach are determined by the \(\mathbf{F}_{\mathrm{C},ij}\) force. The governing equations for forces are described elsewhere [14; 15; 16; 17]. \[\mathbf{F}_{\mathrm{D}}=\frac{1}{2}C_{\mathrm{D}}\rho A_{\mathrm{P}}|\mathbf{ u}-\mathbf{u}_{\mathrm{P}}|(\mathbf{u}-\mathbf{u}_{\mathrm{P}}) \tag{12}\]  
where \(\mathbf{u}\), \(\rho\) and \(A_{\mathrm{p}}\) are fluid velocity, fluid density and cross section area of the particles and for \(Re_{\mathrm{P}}\leq 1000\):  
\[C_{\mathrm{D}}=\frac{24}{Re_{\mathrm{P}}}\bigg{(}1+\frac{1}{6}Re_{\mathrm{P}} \tfrac{2}{3}\bigg{)} \tag{13}\]  
The \(C_{\mathrm{D}}\) is the drag coefficient. \(Re_{\mathrm{P}}\) is particle Reynolds number and is expressed as follows:  
\[Re_{\mathrm{P}}=\frac{\rho|\mathbf{u}-\mathbf{u}_{\mathrm{P}}|d_{\mathrm{P}}} {\mu} \tag{14}\]  
In this equation, \(d_{\mathrm{P}}\) is the particle diameter and \(\mu\) is air viscosity. The Brownian term \(\mathbf{F}_{i}^{\mathrm{B}}\) defined by Equation (15) is stochastic random vector and \(\Delta t\) is the solution process time step.