A fractional order model of the COVID-19 outbreak in Bangladesh

1.   Introduction

Euler-Lagrange fluid-particle simulation has been extensively used as a surrogate for experiments in respiratory aerosol dynamics and inhalation drug delivery. This approach's popularity is rooted in its high accuracy in resolving fine-scale flow details and tracking individual droplets/particles in a time-efficient manner and thus providing an estimation of the dosimetry of released pharmaceuticals ^[1]. However, an important issue often neglected in numerical simulations of pulmonary drug delivery was the momentum exchange from the discrete phase to the carrier fluid or the reciprocal interactions between aerosols and flow (i.e., two-way coupling, or 2WC). One-way coupling (1WC) was adopted in most previous inhalation dosimetry simulations based on the assumption of dilute aerosol concentrations, where the presence of the aerosols did not affect the carrier flow ^[2]. This assumption was valid only when aerosols had a sufficiently low concentration and low exiting velocities, such as those from a mesh nebulizer. One-way coupling was also valid for aerosols from a valved holding chamber (VHC) or spacer connected to a metered-dose inhaler (MDI) due to the mixing of the spray droplets within it, which greatly reduced droplets' speed and concentration ^[3].

However, the local concentration of MDI spray aerosols can be high upon discharge close to the actuator orifice, which has a small diameter of 0.2–0.5 mm ^[4,5,6]. The pressurized propellant generates a high-speed jet flow from the orifice and atomizes the drug solution into micro-scale droplets, which develop into a spray plume by advancing along the streamwise direction and dispersing in transverse directions. Because of the large density ratio between the droplet and fluid (i.e., three orders of magnitude), the inertia and response time of the droplets differ notably from those of the carrier flow, which will induce large droplet-flow slip velocities. As a result, the momentum exchange from the spray droplets to the co-flow can be significant downstream of the actuator orifice and modify the flow behaviors there, which in turn will modify aerosol motions. These reciprocal interactions kick in from the onset of the actuation and occur throughout the actuation process (usually 0.2 s). Spatially, the 2WC effect is more significant near the orifice and persists till either the local mass loading or aerosol velocities become sufficiently low. Due to their higher inertia, aerosols often exhibit an integrative feature, i.e., the current aerosol distribution (or droplet location) depends on not only the instantaneous velocity distribution but also the prior histories of the droplets. By contrast, airflow can quickly respond to external changes and is predominantly determined by the instantaneous pressure field. As a result, the motions of individual droplets will be constantly modified by the carrier flow, and these modifications will accumulate throughout the time course, yielding even more different spray plume morphology and velocity than in the 1WC case.

Many experimental and numerical studies have considered the 2WC effects on particle-laden flows in general. Early observations by Torobin and Gauvin demonstrated that particles modified the drag force and heat transfer rate on the pipe, which couldn't be explained unless the flow turbulence had been modulated by the entrained particles ^[7]. Questions also arose about whether, at what conditions, and how the turbulence would be reduced or enhanced by the dispersed particles. It was known that when the aerosol droplets are comparable in size to the smallest fluid length scale or larger, unsteady flows can be important in affecting the aerosol motions ^[8,9]. Gore and Crowe reviewed a wide spectrum of experimental studies on particle-modulated turbulence in jets and pipes at varying orientations, with a Reynolds number of 8000–10⁵, volume fraction of 10^-6–0.2, and density ratio of 10^-3–2500 ^[10]. They first observed that particles smaller than 200 µm suppressed turbulence while particles larger than 200 µm enhanced turbulence. A non-dimensional parameter d_p/l_t (particle diameter to turbulence length scale) was then proposed, with turbulence attenuation at d_p/l_t < 0.1 and enhancement at d_p/l_t > 0.1 ^[10,11]. The percentage change in turbulence intensity varied from -50% to 400%, and exhibited a nonlinear relation to d_p/l_t ^[10,11]. Later, Elghobashi and Truesdell updated the criterion to dp/η, with η being the smallest turbulence length scale or Kolmogorov length scale ^[12].

Note that the particle-induced turbulence modulation can reversely modify particle motions at varying degrees, depending on the volume fraction, mass fraction, and/or diameter of the particles ^[13]. Using direct numerical simulations (DNS), Monchaux & Dejoan showed that 2WC had a weak effect on the global flow statistics when the particle loading is low (φ = 1.5×10^-5). They also reported that 2WC enabled the particles to drag the fluid down with them, which further enhanced the particle settling ^[14]. In a classic review, Elghobashi categorized the particle-laden turbulent flows into three regimes based on the particle volume fraction φ: one-way coupling (1WC) regime when φ = 10^-8–10^-6, two-way coupling (2WC) regime when φ = 10^-6–10^-3, and particle-particle collision regime when φ > 10^{-3 [15]}. The 1WC and 2WC regimes had dilute suspensions while the particle-particle collision regime had dense suspension. In the 2WC regime (φ = 10^-6–10^-3), a sub-division was made based on the time scale ratio, τ_p/τ_η where τ_p was the particle residence time, and τ_η was the Kolmogorov time scale. For particles with τ_p/τ_η > 100, 2WC enhanced turbulence production, and for particles with τ_p/τ_η < 100, 2WC enhanced dissipation ^[15].

Till now, studies of the 2WC effects on MDI spray flow and aerosol dynamics have not been reported. Considering the high droplet concentration close to the actuator orifice and the high-speed jet flow that generates small turbulence length scales, 2WC effects are expected to be significant in the orifice vicinity and decay from the orifice. It is unknown yet to what extent the spray plume will be affected and at what level it is affected in different directions (i.e., spray centerline vs. transverse directions) and under various conditions (i.e., mass loading, jet flow speed, and aerosol diameter). This study seeks to answer the above questions related to three popular MDIs, i.e., Ventolin, ProAir, and Qvar, which have different droplet size distributions and co-flow velocities. Specific aims include:

1) Develop an LES-Lagrangian computational model for MDI spray plume development in open space and validate it with complimentary high-speed imaging of Ventolin sprays.

2) Compare the spatial and temporal evolutions of the Ventolin spray plumes between 2WC and 1WC predictions.

3) Evaluate the 2WC effects with mass loadings, droplet sizes, and inhaler types.

4) Understand diameter-dependent spray plume behaviors by correlating to flow Kolmogorov length/time scales and respective droplet-flow length/time ratios.

The results of this study will shed light on whether 2WC effects should be considered in MDI spray simulations and provide quantitative guidance in estimating the differences between 2WC and 1WC cases.

2.   Materials and methods

2.1.   Study design

The 2WC effects on MDI airflow and aerosol dynamics were numerically evaluated in an MDI-open-space model, as shown in Figure 1(a). The MDI was actuated, and the generated spray droplets were released into the open space. A cone-shaped space was adopted to account for the expansion of the spray plume and, at the same time, minimize the computational mesh size. An airspace geometry had a length of 0.5 m based on the observations that the maximal spray plume lengths from all three inhalers considered hereof were equal to or smaller than 0.5 m.

Figure 1(b) shows the model for Ventolin, including both the geometrical model and the model for the initial airflow and spray aerosols at actuation. The inhaler geometry was developed in SolidWorks (Dassault Systèmes, Concord, MA) based on the Ventolin inhaler with the most essential details included. Examples included the exact size of the mouthpiece and the actuator orifice indented in the cylindrical reservoir with an orifice diameter of 0.3 mm (Figure 1(b)). The orifice diameter of 0.3 mm followed Gabrio et al. ^[16], who evaluated the plume impact forces from several prototype MDIs (including Ventolin) with two different orifice diameters of 0.32 and 0.48 mm. They also measured the throat deposition from model actuators with orifice diameters ranging from 0.29 to 0.56 mm and reported that the throat deposition was sensitive to orifice diameters up to approximately 0.4 mm ^[16]. In this study, only the prototype actuator diameter of 0.32 mm (rounded to 0.3 mm) was considered. The exiting velocity of the airflow and aerosol was 40 m/s, which was based on previous measurements and computations ^[17,18,19]. A polydisperse aerosol size distribution was adopted (Figure 1(b)) based on previous measurements, with the mean droplet size being 9.1 µm and the standard geometrical derivation (GSD) being 1.32 ^[17,19]. The inhaler actuation started from 0.1 s and lasted for 0.2 s. After that, no further airflow and aerosols were released from the orifice, while the aerosols that had been released between 0.1–0.3 s would continue to be simulated till their complete decay.

An MDI is supplied with a canister that contains 8.5 g of formulations intended for 200 puffs ^[20,21]. Each actuation delivers 90–120 μg medications (i.e., 0.2–0.3% w/w), such as albuterol from Ventolin, albuterol sulfate from ProAir, and beclomethasone dipropionate from Qvar ^[22]. In addition to the 0.1–1.0% w/w medication, a formulation also contains 3–13% w/w ethanol to increase the drug solubility and propellent to atomize the drug formulation into micron-scaled droplets ^[23,24]. In this study, it was assumed that the propellant evaporated immediately and completely at the orifice, leaving the liquid droplets that contained only the medication and ethanol. As a result, the aerosol droplet mass per actuation was approximately 3.0 mg (8500 mg, 200 puffs, 7% medication plus ethanol). The computational model was first applied to simulate the Ventolin spray plume development and was validated against complimentary high-speed images. This validated case with 2WC was then used as the baseline case for later comparisons. Various spray properties were considered regarding their relative contributions to the 2WC effects. These included the presence of the propellant jet flow, the spray mass loading, and individual droplet sizes. All cases were simulated with both 2WC and 1WC to characterize the 2WC effects under different scenarios.

To evaluate the influences of the propellant flow, the spray plume behaviors were simulated with a 40-m/s velocity for aerosols but a zero velocity for the airflow from the orifice during actuation (0.1–0.3 s) under both two-way and one-way coupling. To evaluate the mass loading effects, three other drug masses were considered, which were one-tenth, ten times, and one hundred times (i.e., 0.30, 30, and 300 mg) that of the baseline (3.0 mg). To study the droplet diameter effects, seven monodisperse aerosols with a dose of 3.0 mg were considered, with the aerosol sizes being 2, 3.5, 5, 7.5, 10, 15, and 20 µm. Two other inhalers, ProAir and Qvar, were considered in addition to Ventolin, with their aerosol size distribution shown in Figure 1(c). Based on the measurements of Liu et al., the exiting spray velocity was 26 m/s for ProAir and 23 m/s for Qvar ^[17]. A summary of the test cases is listed in Table 1, with the corresponding figure numbers of the results.

2.2.   High-speed imaging

The high-speed camera used was Phantom VEO, with an acquisition rate of up to 11,000 frames per second (fps). In this study, the temporal and spatial evolutions of the Ventolin spray plume after actuation were recorded at 4000 fps in order to achieve the optimal quality of acquired images. A laser sheet (OXlaser, 488 nm, 100 mW) was implemented to enhance the spray-environment contrast by aligning the sheet with the spray plume centerline. To understand the spray plume development at the very beginning of the actuation, the mouthpiece was partially removed to reveal the orifice. Images of the spray from the orifice were thus directly captured without the blockage by the mouthpiece.

2.3.   Numerical methods

The LES-WALE model was used to simulate the multi-regime flows because of its ability to accurately capture turbulent-laminar transitions ^[25]. Droplet dynamics were simulated with a Lagrangian approach. An in-house MATLAB code was used to generate MDI droplets at the orifice, which specified the droplet diameter, size distribution, exit velocities, and plume angle ^[26,27]. The 2WC between the airflow and aerosol droplets were considered through the force below,

where µ is the fluid viscosity, ρ_p is the particle density, d_p is the particle diameter, R_e is the relative Reynolds number (R_e = ρ(u_p – u)l/µ), C_D is the drag coefficient, u_p is the particle velocity, u is the fluid velocity, ṁ_p is the aerosol mass flow rate, Δt is the time step, and F_other is a collection of other interaction forces, which are zero in this study. The mutual momentum transfer between the fluid and aerosols was numerically considered in each control volume for a given instant.

ANSYS Fluent 19.1 (Canonsburg, PA) was utilized to solve flow-particle governing equations. ANSYS ICEM CFD was utilized to generate the computational mesh (Figure 1(d)). The characteristic length scales in the MDI-airspace system differed by three orders of magnitude, i.e., 0.3-mm diameter for the actuator orifice, 2.54-cm width for the mouthpiece opening, and a 0.5-m length for the open space. To adequately resolve the flow details in all regions, multidomain meshes were generated by dividing the system into many zones, with each zone being filled with varying-sized cells according to the anticipated flow complexities. For instance, ultrafine meshes were generated close to the orifice and within the shear layers around the spray centerline, with the mesh density gradually decreasing from the centerline, which eventually became coarse in the far open space (Figure 1(d)). Grid independent studies were conducted by testing six meshes from 3.2 million to 12.0 million (Figure 1(d)). The parameter of interest was the airflow 3 cm from the mouthpiece at 0.5 s. The grid-independent mass-based DF was reached at 10.1 million. A total 2×10⁵ seed droplets were released for each actuation during 0.1–0.3 s. To approximate the continuous releasing of sprays during actuation, a group of 5000 droplets was released every 5 ms for 40 times, as opposed to a single injection in most previous numerical studies ^[28,29,30]. For an applied mass of 3.0 mg, the mass flow rate for each droplet parcel in the injection file was calculated to be 3.0×10^-9 kg/s based on a droplet parcel number of 5000 and a time step of 0.005 s per release for 40 releases. The mass loading of 0.30, 30, and 300 mg was obtained by scaling the baseline particle mass flow rate per parcel by a factor of 1/10, 10, and 100, respectively. A time step size of 5 ms was used for a flow time of 1.5 s, or a complete spray decay, which came first. One test case took 60 hours or so in a Ryzen 9 3960x workstation with 3.79 GHz frequency and 256G RAM.

2.4.   Kolmogorov scales

A deeper understanding of the fluid-aerosol interactions involves those at small scales, especially when the length and time scales are equivalent between the fluid and individual droplets. The Kolmogorov scales are the smallest in turbulent flows and describe the microscopic airflow behaviors being dominated by viscosity, such as how much energy is contained in the smallest eddies, how energy is transferred between small and large eddies, and how energy is dissipated by eddies ^[31]. The Kolmogorov length scale, η and time scale, τ_η of the airflow are modeled as ^[32],

where ν is the molecular viscosity and ɛ is the dissipation rate per unit mass, as shown below.

The length scale of a particle or droplet is obviously its diameter, and the time scale is described by the residence time, i.e., the duration required to react to a local flow change,

The particle/droplet Stokes number was also presented above, which is the ratio of the particle/droplet residence time over the characteristic macroscopic time scale of the flow (l₀/u₀). Here l₀ and u₀ are the macroscopic characteristic length and velocity scale in a local region. Considering the multiscale dimensions (three orders of magnitude differences) of the system in this study, the local characteristic scales also varied drastically, which further led to drastic changes in the particle/droplet Stokes number.

The aerosol-fluid length scale ratio (d_p/η) and time scale ratio (τ_p/τ_η) have been suggested as indexes to evaluate the effects of aerosol-to-fluid forces on the flow turbulence ^[15]. It has been demonstrated that the presence of aerosols suppressed turbulence when d_p/η < 0.1 and enhanced turbulence when d_p/η > 0.1 ^[10,11,12]. Both ratios were used to evaluate the spray plume morphology variations for different aerosol sizes.

3.   Results

3.1.   Ventolin spray plume visualization and simulation (baseline)

3.1.1.   High-speed imaging of Ventolin actuation

Figure 2(a) shows the recorded images of the spray plume development from Ventolin at two instants using the Phantom VEO high-speed camera at an acquisition rate of 4000 fps. Intensive turbulence and eddies were observed at the front of the spray plume and, to a lesser degree, at the transverse interface with the ambient air. Considering the plume front, an aerosol stripe developed and moved downward (red arrow, second panel, Figure 2(a)), which subsequently separated from the main plume and started to settle downward from the gravity (yellow arrow, third panel, Figure 2(a)). Simultaneously, the eddies at the transverse interface continued to grow by mixing with the ambient air due to large shear forces. Because of the front and shear forces, as well as the entrained mass of the air, the overall velocity of the spray plume decreased quickly, from 40 m/s at the orifice to nearly zero at 0.5 m. The morphology of the spray plume exhibited a highly dynamic nature even after the spray plume reached its maximal length, driven both by the decaying vortices and gravity.

Considering that the orifice was blocked by the mouthpiece, parts of the mouthpiece were removed, and the emanation process of the spray from the orifice was captured using the high-speed camera at 4000 fps (Figure 2(b)). It can be clearly seen that the spray plume exited the orifice as a very thin jet with a cross-section diameter comparable to that of the orifice (first recorded image, Figure 2(b)). The plume grew quickly in the axial direction, and at a slower speed, in the transverse direction as well (second and third recorded images, Figure 2(b)). However, no apparent recirculation or vortices were noticed close to the orifice and within the mouthpiece. These observations were used to validate the complimentary computational results in the later section.

3.1.2.   LES-Lagrangian simulations of Ventolin spray evolution

Figure 3 shows the temporal and spatial evolution of the spray plume at different instants under 2WC. There were three stages during the plume development. The first stage occurs during 0.1–0.3 s (actuation, upper row of Figure 3), with different dynamics between the carrier flow and spray droplets. Due to the large density and high inertia, droplets traveled a much longer distance (~ 13 mm) in the first 5 ms than the airflow, with a vortex ring forming close to the mouthpiece exit (red arrow in Figure 2(a)). However, due to the drag force from the ambient air, the spray plume slowed down quickly, with further advancement of 7.5 mm at the next 15 ms (from 0.105 to 0.120 s) and 3 mm at the next 80 ms (from 0.120 to 0.200 s), as shown in Figure 2(a)–(c). By contrast, the coherent flow structures developed steadily. The single vortex ring generated from the mouthpiece at 0.105 s (red arrow in Figure 2(a)) moved forward, while new vortex rings were constantly generated at the mouthpiece (black and green arrows in Figure 2(b)). The green ring caught up with the red ring at 0.200 s and merged together, while younger rings (pink and black arrows, Figure 2(b)) were continuously generated and moved forward. At 0.3 s, the vortex rings caught up with the droplet plume and disturbed the shear-shaped spray plume into a more irregular morphology (Figure 3(d)).

In the second stage (0.3–0.8 s), vortex rings formed and moved forward, carrying the spray droplets with them. These vertex rings interacted with each other, as well as with the ambient air, leading to increasingly irregular patterns. Recirculating flows were induced from the interactions between the vortex and ambient air, which slowed down the advancement of the vortices and distributed the entrained droplets in transverse directions. For the aerosol droplets that fell below the spray centerline, the gravitational force became more dominant and started to settle in the negative y-direction (filled hollow arrow).

The third stage occurred after 0.8 s when no apparent vortices were generated at the mouthpiece. The vortices progressively decayed while all aerosol droplets started settling down under gravity, with larger droplets falling faster, as illustrated by the hollow arrow in Figure 3(h), (i).

The airflow development under the 2WC effects is further visualized in Figure 4 in terms of the velocity contour lines and streamlines. The jet flow was apparent in Figure 4(a) at the beginning of the actuation, with recirculating flows forming at the mouthpiece. The recirculating flow increased its size from 1.05 to 1.20 s by mixing with the surrounding air (red arrow, Figure 4(b)). It was also noted that jet flow pushed the ambient air at the front (green arrow) and sucked in the surrounding air due to low pressures induced by the high-velocity jet (pink arrow, Figure 4(b)). The recirculating flows above and below the jet centerline were approximately symmetrical, indicating an insignificant effect of the current mass loading (3.0 mg) on the airflow dynamics along the jet centerline (Figure 4(c)–(e)). The recirculating flow zones advanced much slowly than the initial discharge velocity (i.e., ~0.4 vs. 40 m/s), indicating the quick decay of the spray plume energy due to the strong mixing with the ambient air. Also, note that the recirculating flow below the jet centerline slightly lagged behind the one above the centerline; this was because there were more large droplets in the lower region, which slowed down the recirculating flow (red arrow, Figure 4(f)). Moreover, the streamlines in the left lower region (filled arrow) bent downward because of the predominate gravitational settling in this region.

3.2.   Ventolin spray simulations

3.2.1.   Ventolin spray simulations: 2WC vs. 1WC

Remarkable differences in spray dynamics were observed with and without the consideration of momentum exchange from aerosols to airflow during the first 0.3 s after actuation, as shown in Figure 5. The MDI was actuated at 0.1 s, and the applied dose was 3.0 mg. One obvious difference was the backflow into the mouthpiece when neglecting the aerosol-to-flow effects (i.e., 1WC from flow to aerosol only). This was reasonable because the spray droplets moved faster than the airflow and thus imparted a forward momentum onto the airflow, which effectively prevented the occurrence of backflows. By contrast, the airflow experienced a large resistance from the upfront pressure, which would push a fraction of airflow backward when neglecting the aerosol-to-flow momentum transfer. These backflows would entrain some droplets into the mouthpiece, as demonstrated in Figure 5(b). The upfront pressure resistance also gave rise to the mushroom-shaped vortex rings at the mouthpiece exit (blue arrow in Figure 5(b)), a phenomenon often reported in jet flows ^[33,34,35].

The second major difference in spray plume evolution between two-way and one-way coupling is the vortex generation and transportation along the centerline. With 2WC, the airflow had a higher momentum and, therefore, exhibited a higher rigidity as a bulk. A more coherent vortex formed (hollow arrow, 0.150 s, Figure 5(a)) than its counterpart with 1WC (hollow arrow, 0.15 s, Figure 5(b)). The vortex with 2WC moved faster, which caught up with or even surpassed the large droplets (red spheres) at 0.3 s (or 0.2 s after actuation), as denoted by the filled arrow in Figure 5(a). By comparison, the vortex under 1WC was still 26 mm short of the large droplets (0.3 s, Figure 5(b)). As alluded to in Figure 3, the vortex would persist for more than one second, which moved forward at a much-reduced speed (~0.4 m/s) than the discharge speed (40 m/s); it transported the polydisperse aerosol droplets along the axial direction, mixed with the surrounding air in transverse directions, and led to complex spray plume morphologies, as illustrated at 0.4 s in Figure 5(a), (b).

3.2.2.   Effects of the spray mass loading

Effects of the mass loading on the airflow and aerosol dynamics are shown in Figure 6 by varying the applied dose mass by three orders of magnitude (0.25–250 mg). Both the streamlines, vortex structures, and aerosol droplets are shown at 0.4 and 0.9 s (i.e., 0.3 and 0.8 s after the actuation), with the vortices being 60% transparent to avoid blocking the aerosol droplets. At a very low mass loading (i.e., 0.25 mg, or one-tenth of the typical MDI dose), the airflow pattern and aerosol distribution resembled those under 1WC. The vortices above and below the axial centerline kept relatively symmetric at both 0.4 and 0.9 s, reflecting the negligible effects from the discrete droplets (Figure 6(a)). By contrast, a high mass loading (i.e., 25 mg, or ten times of typical MDI dose) led to a considerable delay of the vortex below the centerline because of the larger mass that, in turn, resulted from increased gravitational settling at this higher mass loading. Particularly, due to the faster motion of the upper vortex, the core flows bent downward, further escalating the spray plume asymmetry (Figure 6(b)). In Figure 6(c), increasing the applied spray dose to 250 mg resulted in even a higher level of asymmetry, with the majority of airflow streamlines bending downward at 0.9 s. As a result, the mass loading could have a significant impact on the spray plume evolution for a spray dose typical of MDI inhalers or higher; 2WC should be implemented in simulations of MDI drug deliveries.

3.3.   Effects of the spray droplet size

3.3.1.   Airflow and aerosol dynamics vs. droplet size

We further studied the airflow dynamics for monodisperse aerosols with varying droplet diameters, as shown in Figure 7. Significant differences were observed in the airflow pattern among different droplets (Figure 7(a)–(d)). For a given dose of 3.0 mg, the airflow with the 2-µm aerosol advanced a much short distance in comparison to the 1WC case (Figure 7(a) vs. 7(e)). The advancement distance increased with the droplet diameter of the monodisperse aerosol from 2 to 20 µm (Figure 7(a)–(d)). Compared to 1WC, where the airflow was not affected by the presence of aerosols, the spray plume length was shorter for 5-µm aerosols but longer for 10-µm aerosols. Moreover, the spray plume persistently expanded in transverse directions with the increasing aerosol droplet size, with the 20-µm spray aerosols exhibiting a much wider plume angle than other aerosols (Figure 7(d) vs. 7(a)–(c)).

The corresponding monodisperse aerosol dynamics under 2WC is shown in Figure 8(a)–(d) at 0.4 s after actuation for droplet diameter of 2, 5, 10 and 20 µm, all with a dose of 3.0 mg. For comparison purposes, the polydisperse aerosol dynamics under 1WC were also plotted, as shown in Figure 8(e). Very different patterns of aerosol distributions were noted for different droplet sizes, which, however, were generally consistent with the velocity contours in Figure 7. As expected, the largest differences occurred for the largest droplets hereof (i.e., 20 µm, Figure 8(d)).

3.3.2.   Kolmogorov scales and ratios (mid-plane)

To understand the strikingly distinct behaviors of the airflow and aerosols for different droplet sizes in Figures 7 and 8, the theory of Kolmogorov scales and ratios were extracted from the spray simulations at 0.5 s under 1WC (Figure 9). Examining 1WC simulations had the advantage of evaluating the tendency of flow/aerosol behaviors without complications caused by the fluid-aerosol interweaving under 2WC. For the spray flow at 0.5 s, the smallest Kolmogorov length scale, η, was 2.72 µm, which led to the maximal magnitude of the aerosol-flow length ratio, d_p/η, being 0.74 for 2-µm droplets, 3.68 for 10-µm droplets, and 7.35 for 20-µm droplets (Figure 9(a)). Note that both scales, at their minimal magnitudes, occurred in a very small region immediately downstream of the actuator orifice (first panel, Figure 9(a)). Both scales increased quickly in magnitude at other locations away from the orifice.

The distribution of the length ratio d_p/η is shown in the second to fourth panels in Figure 9(a) based on a color code range of 0–0.1. For 2-µm droplets, the d_p/η magnitude was less than 0.1 for almost all regions except for a very slim region close to the orifice (second panel, Figure 9(a)). Accordingly, the turbulence development was noticeable suppressed, leading to weaker turbulence and vortex intensities, which in turn yielded a shorter penetration distance than when aerosol-to-flow forces were neglected, as previously described in Figure 8(a) vs. 8(e). By contrast, the 10-µm droplets had a much larger region with d_p/η > 0.1 (red zone, third panel, Figure 9(a)). In this region, the forces from the aerosols would enhance the turbulent flows and increase the flow strength, as aforementioned in Figure 8(c) vs. 8(e). The d_p/η distribution largely resembled the morphology of the orifice jet flow, with the flow-enhancing zones along the axial centerline, which decayed quickly in the transverse directions (third panel, Figure 9(a)). Considering 20-µm droplets, the red zones enlarged in both axial and transverse directions, reflecting a much larger impact from the spray droplets. This observation was consistent with the highly expanded spray cloud for 20-µm droplets, as exhibited in Figures 7(d) and 8(d).

Figure 9(b) shows the Kolmogorov time scale and the droplet-specific time-scale-ratio, τ_p/τ_η for 2, 10, and 20-µm aerosols. For the spray flow at 0.5 s, the smallest Kolmogorov time scale, τ_η, was 5.1e-7 s, which was smaller than the residence time, τ_p, for all droplet sizes considered (2–20 µm), as listed in Table 2. It is also noted that the smallest Kolmogorov time scale (τ_η = 5.1e-7) existed in a very small region only; there was an appreciable region along the spray centerline with τ_η ranging from 1.0e-5 to 1.0e-3, which was equivalent to the residence time τ_p for droplets ranging from 2 to 20 µm. As a result, the airflow dynamics would be noticeably modified by the entrained droplets downstream of the actuator orifice, which was also demonstrated to modify the droplet behaviors and fates, as discussed in detail in Figures 7 and 8.

3.3.3.   Quantification of Kolmogorov length scale and ratio (1-D plot)

The Kolmogorov length scale, η, and ratio, d_p/η, were further quantified in the axial and transverse directions downstream of the actuator orifice at 0.5 s (Figure 10(a)). As alluded to in Figure 9, the η magnitude increased with the distance from the orifice because of the decaying eddies or vortices away from the orifice (Figure 10(a)). Figure 10(b) shows the length scale ratio d_p/η vs. the axial distance (i.e., Y-direction). As expected, in most regions, the d_p/η magnitude for 2-µm droplets was lower than 0.1 (pink dashed line). By contrast, for 10-µm and 20-µm aerosols, the d_p/η magnitude was larger than 0.1 for Y < 12 cm. A bump was observed at Y = 1–3 cm for 10-µm and 20-µm aerosols (Figure 10(b)) because the spray actuation had been completed at 0.5 s and the core jet flow started to decay thenceforth.

The distributions of the Kolmogorov length scale in the circumneutral directions (Z- and X- directions) are shown in Figure 10(c), (d) at different distances from the orifice. For all cross-sections considered, the Kolmogorov length scale η magnitude (the eddy scale) increased quickly from the axial centerline and became larger than 200 µm (10 times of the 20-µm droplets) when the radius from the centerline was larger than 0.5 cm. For a given radius, the η magnitude also decreased with the axial distance. At the cross-sections of Y = 3 and 6 cm, the η distributions were approximately symmetric upper-lower (Z) or left-right (X). This symmetry was also observed in the X-direction (left-right) at Y = 9 and 12 cm (two solid lines, Figure 10(d)). However, apparent asymmetry was observed in the gravitational Z-direction, with a much higher value above the axial centerline than below.

3.3.4.   Relative flow-droplet velocities

The particle Reynolds number has been suggested as an index to differentiate the turbulence suppression and enhancement. In Figure 11, the flow and aerosol velocities were examined for the 2-µm aerosols under 2WC at t = 0.5 s. Figure 11(a) shows the airflow velocities at the locations that are coincident with the instantaneous particles. A wide range of velocities was observed, with low-speed flows in most of the regions downstream of the mouthpiece in contrast to a high-speed jet flow along the centerline. To better understand the velocity distributions, a slice view was presented in the lower panel of Figure 11(a), which only showed the aerosols within a 2-mm-thickness slice in the midplane. Without obscuring from surrounding droplets, the jet flow and vortex (pink square) were clearly displayed. The vortex was further visualized in a zoomed view, where flow recirculation was apparent.

The instantaneous velocities of the spray droplets within the midplane slice (2 mm thickness) are presented in Figure 11(b), which were much higher than the airflow velocities. Unlike the various directions of the airflow, most droplets exhibited a forward direction. The differences in the velocity and direction between airflow and aerosol droplets are shown in the right lower panel of Figure 11(a), where both the instantaneous velocity vectors of the flow and droplets are presented. Note that the flow and droplet velocity vectors were not to the scale in order to plot them in the same frame. It was observed that both the magnitude and direction could differ significantly between the airflow and an individual droplet, and that these differences varied among droplets. This made it challenging to use the particle Reynolds number as a practical criterion for turbulence-depression-enhancement.

3.3.5.   Decay of flow and aerosol velocities

The effects of the monodisperse aerosol diameter on the spray plume evolution were evaluated by comparing the aerosol velocities at two points (i.e., 3 and 6 cm downstream of the mouthpiece). Both the maximal and average velocities were quantified and presented in Figure 12(a) as a function of the droplet diameter (2, 3.5, 5, 7.5, 10, 15, 20 µm). The aerosol velocity peaked around 5–7 µm for both the maximal and average values at both sampling points considered. The velocity decreased thereafter and then increased again from 15 µm (Figure 12(a)). The velocity decay from 3 to 6 cm had a similar magnitude for both the maximal and average velocities (Figure 12(a)).

Figure 12(b) shows the velocities of the airflow along the axial direction (1–8 cm from the mouthpiece) for different aerosol diameters. Overall, the aerosol speeds are much higher than those of the carrier flow (Figure 12(a) vs. 12(b)). The airflow slowed down quickly after exiting the mouthpiece. However, the rate of decrease differed among aerosols, with 2-µm aerosols decaying most drastically and 5-µm aerosols decaying the slowest, which was consistent with the trend observed in Figure 12(a).

3.4.   Two-way vs. one-way for ProAir and Qvar

The effects of the 2WC on different MDI types were evaluated by considering the other two inhalers: ProAir with a mean droplet size of 5.2 μm emitted at 26 m/s, and Qvar inhaler with a mean droplet size of 2.9 μm emitted at 23 m/s ^[17]. As expected, the predicted airflow dynamics without the aerosol impact were very similar regardless of the inhaler types (Figure 13(a) vs. 13(b)). The differences in the one-way-predicted spray morphologies at 0.5 s between ProAir and Qvar were noticeable but not significant. These slight differences were presumably attributed to the differences in the droplet sizes (5.2 vs. 2.9 µm) and initial velocities (26 vs. 23 m/s). Remarkable differences were predicted with 2WC in both the airflow and aerosol dynamics between Proair and Qvar (Figure 13(a) vs. 13(b)). Even though the airflow was suppressed by the entrained spray aerosols from both inhalers, that from Qvar was significantly more suppressed due to its smaller droplet sizes and lower d_p/η ratio. As a result, the Qvar spray plume was much less dispersed and advanced a shorter distance (Figure 13(b)).

The airflow and aerosol dynamics under 2WC after 0.5 s are presented in Figure 14(a), (b) for ProAir and Qvar, respectively. For ProAir (with a mean aerosol diameter of 5.2 µm), a smaller vortex formed below the axial centerline, as denoted by the hollow red arrow in Figure 14(a). The airflow posterior to the vortex bent downward and became increasingly more apparent with time (hollow blue arrow in Figure 14(a)), indicating an increasing gravitational effect from large spray droplets. The bending-downward flow and the small vortex were not observed in Qvar. This was reasonable considering the much smaller aerosol mean diameter (2.9 µm) than that of ProAir (5.2 µm). Moreover, the spray plumes from Qvar were less dispersed than those from ProAir and advanced a shorter distance than their counterparts in ProAir. The dipole vortices below and above the axial centerline were slightly off asymmetric at t = 1.0 s, but overall, the bending-downward motions of airflow or aerosols were not observed (third panel, Figure 14(b)). The differences shown in this figure, together with Figure 7, indicated that the 2WC effects are sensitive to the spray droplet size, and thus it will be essential in numerical studies of MDI pulmonary drug delivery to consider the 2WC effects for different types of inhalers, which often had very different spray droplets and exit velocities.

4.   Discussion

4.1.   Mass loading and 2WC

The magnitude of the 2WC force (Eq (1)) was proportional to both the slip velocity and local mass loading. It was demonstrated in this study (Figures 4 and 5) that a moderate mass loading typical of an MDI spray dose per puff (e.g., 3.0 mg) could alter the fluid flow and aerosol behaviors. For an applied dose ten times smaller (i.e., 0.30 mg), no significant differences were predicted between two-way and one-way coupling (Figure 7). However, for an applied dose ten times larger (i.e., 30 mg), both the flow and aerosols were remarkably modified under two-way-coupling. The jet flow became apparently asymmetric above and below the spray centerline while the front of the spray plume started to bend downward. In other words, the gravitational effect became more pronounced at higher mass loadings, which drove the droplet-laden flow toward the ground (Figure 7). These observations were in agreement with Tong and Wang, who reported that moderate loadings of particles could induce noticeable flow modulation ^[36]. Similarly, Monchaux and Dejoan reported that 2WC strongly enhanced particle settling than the 1WC case and that these particles dragged the fluid downward ^[14]. These also conformed to the three regimes of particle-laden flows proposed by Elghobashi et al. ^[15]. At very low volume fractions (φ < 10^-6), the discrete phase did not affect the carrier flow and was regarded as the 1WC regime. In the second regime, the discrete phase had a low volume fraction but a moderate mass loading due to the high droplet-flow density ratio. In this case, droplets could perceivably modulate the co-flow (i.e., 2WC regime), but droplet-droplet interactions could be neglected. Note that the MDIs fell within this regime. In the third regime (i.e., dense flow regime), both the volume fraction and mass concentration were high, and collisions among droplets had to be included.

4.2.   Droplet size and Kolmogorov scales on two-way and one-way coupling

To evaluate the droplet size effect on 2WC, monodisperse aerosols of varying sizes (2–20 μm) were simulated and compared for a given applied mass of 3.0 mg. Strikingly different spray plume morphologies were observed for different droplet sizes, with 2-μm aerosols significantly suppressed and 20-μm aerosols dispersed in comparison to their counterparts with 1WC simulations (Figures 8 and 9). To find out the underlying reasons for the sticky behavior of small particles and dispersing behavior of large parties, Kolmogorov scales for the airflow and scale ratios between droplets and airflow were quantified at t = 0.5 s. Turbulence modification in particle-laden flows has been acknowledged for several decades ^[7,12,37]. Early experimental observations by Torobin and Gauvin showed that the presence of particles changed the wall drag in pipes, as well as the heat transfer and chemical reaction rates, which could not be explained unless the fluid turbulence had been modified by particles ^[7]. In a comprehensive review, Gore and Crowe listed a spectrum of potential effects of the carrier flow depending on the particle size and volumetric loading. The ratio of the particle size to the turbulence integral length was proposed as an indicator of whether the turbulence was increased or decreased ^[10]. In this study, the ratio of the droplet diameter over the Kolmogorov length scale (d_p/η) was quantified (Figures 10 and 11), which was found to be closely associated with the spray plume morphologies (i.e., advancement in the axial direction and dispersion in transverse directions) (Figures 8 and 9). For instance, the seemingly sticky aerosols of 2 μm had a d_p/η < 0.1 in near all regions of the spray, while the highly dispersed and advanced aerosols of 20 μm had a d_p/η < 0.1 only in the regions further than 15 cm (Figure 11(b)). This finding was consistent with Elghobashi and Truesdell, who reported that turbulence was suppressed when d_p/η < 0.1 and enhanced when d_p/η > 0.1 ^[12].

It was also observed that the transverse dispersion of 20-μm aerosols in Figure 9 was much larger than the regions with d_p/η > 0.1 (red zone) in Figure 10, indicating that the Kolmogorov scale ratio theory might not fully explain the spray suppression/expansion in the transverse directions. One salient feature of MDI sprays was mixing layers at the interface between the jet and ambient air, where intensive vortex formations occurred. Many experimental and numerical investigations have demonstrated that aerosol dynamics in the mixing layers differed significantly from that of the passive tracers if the particle/droplet response time was comparable to the characteristic time scale of the flow. Aerosols would accumulate in the outskirts of vortices, especially in the proximity of braid stagnation points where two counter-rotating vortices met ^[38]. Therefore, the strong mixing and vortex formations at the jet interfaces could have assisted the transverse dispersion in addition to the Kolmogorov-associated turbulence enhancement from the core flow of 20-μm aerosols.

Turbulence modification by particles was also correlated to the particle Reynolds number R_e,_p based on the interphase slip velocity, with a low R_e,_p decreasing the fluid turbulence, and a high R_e,_p increasing the turbulence ^[39,40,41]. In this study, the R_e,_p was also sought by characterizing the slip velocity between an individual droplet and the local flow (Figure 12). However, both the airflow and droplets in the MDI sprays were highly transient in time and heterogeneous in space (Figure 12), making it impractical to quantify the R_e,_p for individual droplets based on their instantaneous slip velocities relative to the local flow. By contrast, the Kolmogorov length scale of the spray flow was more practical to quantify, and its ratio to droplet diameters has been demonstrated to be useful in predicting suppression or enhancement of the carrier flow and aerosol dispersion (Figures 8 and 9 vs. Figures 10 and 11).

4.3.   Effect of inhaler type between two-way and one-way coupling

The high sensitivity of the MDI spray plume to its droplet sizes rested on several MDI-specific factors, such as the high-speed jet flow from the orifice that led to Kolmogorov length scales comparable to the droplet sizes (2–20 μm), and the applied mass loading (3.0 mg in this study) that was large enough to elicit perceivable aerosol-to-flow momentum exchanges. In this study, we demonstrated that reducing the mass loading by ten times (0.30 mg) led to negligible two-coupling effects while increasing the mass loading by ten times (30 mg) led to large variations in spray plume morphologies (Figure 7). It was reminded that the droplet-size-dependent spray behaviors were shown for monodisperse aerosols at a 3.0-mg mass loading (Figures 8 and 9). For polydisperse aerosols like an MDI spray, the transient morphology of the aerosol plume was an instantaneous collective manifestation of all droplets with varying diameters and probability densities. For Ventolin, 2WC provided a better approximation to the experimental aerosol behaviors close to the orifice (Figure 5 vs. Figure 2); it also predicted a more detached aerosol bolus from the mouthpiece and a slightly faster advancement than 1WC (Figure 5). By contrast, for ProAir and Qvar, 2WC predicted a spray plume advancing at lower speeds than 1WC (Figures 14 and 15). However, the Qvar spray development appeared more retarded in both axial advancement and transverse dispersion due to its smaller mean diameter (2.9 μm) and smaller length ratio, dp/η, than those of ProAir (5.2 μm).

Even though only three types of inhalers were considered hereof, the spray behaviors of other MDIs can be roughly estimated based on the droplet-Kolmogorov length scale, which can be readily calculated from the aerosol size and discharge flow speed from the actuator orifice. Another implication was that previous one-way-coupled predictions of the delivery efficiencies from MDI delivery might have errors of varying degrees, especially for Qvar with small aerosol sizes (2.9 μm) and suppressed plume developments. Future studies are warranted to include the 2WC effects in predictive inhalation dosimetry in the respiratory tract for different inhalers.

4.4.   Limitations

MDI actuation and the subsequent spray plume development are highly dynamic and involve multiscale and multi-physics. Multiple assumptions were made in this study to fill the gap of missing data in current literature or mitigate the prohibitive numerical requirements. These included the same discharge velocity for aerosol droplets and co-flow, a constant aerosol output during the 0.2-s actuation, a constant temperature, no evaporation, no electrostatic charge, and no droplet interactions. The MDI formulation contains approximately 80–90% w/w propellant, which vaporizes quickly to atomize the solution into micro-scale droplets; however, measurements of propellent flow speed at the actuator orifice are not yet available ^[24]. Liu et al. measured the spray (aerosol and co-flow) velocities at 3 and 6 cm from the mouthpiece of different inhalers, including Ventolin, ProAir, and Qvar ^[17]. These measurements were used to reversely estimate the discharge velocity at the orifice, as explained in Talaat et al. ^[18,30]. The spray output during actuation can also be time-dependent during the 0.2-s actuation, adding more flow fluctuations. Drastic temperature variation can occur during MDI actuation, leading to thermophoretic effects on droplets ^[42,43,44]. The ethanol droplets can evaporate and decrease in size, further altering the thermal environment and droplet behavior ^[45,46,47]. Particle-particle interactions can also be important at high mass loading, e.g., 300 mg, but should be insignificant due to the low volume fraction ratio of the liquid droplets in the open space. Electrostatic charges on the droplets can be another factor that complicates the trajectories of individual droplets and interactions among neighboring droplets ^{[48,49,50,51,52]}. An orifice diameter of 0.3 mm was used all three inhalers, which was larger than Qvar's actual diameter (0.25 mm) and smaller than Ventolin's diameter (0.48 mm). Its effects on aerosol dynamics were expected to be insignificant considering that the experimentally measured aerosol sizes and speeds were implemented for Qvar. For Ventolin, we acknowledged that our LES simulations could not be regarded as replicating the experimental observations, given the differences between the model and prototype orifice diameters. Rather, the general features of 2WC simulations seemed to be more representative of the Ventolin plume than 1WC simulations. Furthermore, this study considered only spray evolutions in the open space, which can be different from those in a closed space like the respiratory tract ^{[53,54,55,56]}. However, neglecting these compounding factors allowed us to isolate the 2WC effect and examine it vs. the 1WC counterpart systemically under controlled, parametrized conditions, i.e., with different mass loadings, droplet sizes, and inhaler types.

5.   Conclusions

In summary, spray plume evolutions in open space were simulated using the LES-Lagrangian approach under both 2WC and 1WC conditions. Differences in the temporal and spatial morphologies of the spray plumes were characterized, and the reasons underlying the differences were explored using Kolmogorov scale ratios relative to the spray droplets. Specific findings are listed below:

1) Numerical simulations with 2WC seemed to provide better agreement with complimentary high-speed imaging of the spray from the orifice than 1WC.

2) Increasing the mass loading altered the spray morphology, with both the airflow and aerosol droplets bending downward at the spray plume front.

3) The droplet size had a significant impact on the spray evolution in both axial (advancement) and transverse (dispersion) directions, with suppressed sprays for 2-μm aerosols and more dispersed sprays for 20-μm aerosols.

4) The Kolmogorov length scale, η of the spray jet flow were comparable to the MDI droplet sizes. The length scale ratio d_p/η correlated well with spray suppression when less than 0.1 and with spray enhancement when larger than 0.1.

5) Different levels of impact from 2WC were predicted on the spray plume evolution of Ventolin, ProAir, and Qvar in comparison to 1WC.

Acknowledgments

Amr Seifelnasr was gratefully acknowledged for reviewing this manuscript.

Conflict of interest

The authors declare there is no conflict of interest.