InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Computer and Information Science » Numerical Analysis and Scientific Computing » "Advances in Modeling of Fluid Dynamics", book edited by Chaoqun Liu, ISBN 978-953-51-0834-4, Published: November 7, 2012 under CC BY 3.0 license. © The Author(s).

Chapter 11

pMDI Sprays: Theory, Experiment and Numerical Simulation

By Ricardo F. Oliveira, Senhorinha Teixeira, José C. Teixeira, Luís F. Silva and Henedina Antunes
DOI: 10.5772/46099

Article top


Three commercial pMDI actuators with its drug canisters attached
Figure 1. Three commercial pMDI actuators with its drug canisters attached
Chemical structures of the two most common types of propellants (CFC-12 and HFA-134a) used in pMDI
Figure 2. Chemical structures of the two most common types of propellants (CFC-12 and HFA-134a) used in pMDI
Droplet velocity vectors represented in half of a mid-section plane of pMDI plume using a formulation only with HFA-134a, measured at the last moment of the actuation. Adapted from [13]
Figure 3. Droplet velocity vectors represented in half of a mid-section plane of pMDI plume using a formulation only with HFA-134a, measured at the last moment of the actuation. Adapted from [13]
Droplet mean velocity along the axial distance from the nozzle of a pMDI using a formulation only with HFA-134a. Adapted from [13]
Figure 4. Droplet mean velocity along the axial distance from the nozzle of a pMDI using a formulation only with HFA-134a. Adapted from [13]
Graphical representation of the pMDI HFA-134a salbutamol experimental data and its fitting for the Rosin-Rammler, Log-Normal and Nukiyama-Tanasawa distributions. Measurements obtained at 100 mm from the laser beam
Figure 5. Graphical representation of the pMDI HFA-134a salbutamol experimental data and its fitting for the Rosin-Rammler, Log-Normal and Nukiyama-Tanasawa distributions. Measurements obtained at 100 mm from the laser beam
Mass flux as a function of ALR
Figure 6. Mass flux as a function of ALR
Integration of the mass fluxes at constant liquid flow rate
Figure 7. Integration of the mass fluxes at constant liquid flow rate
Fresnel (a) and Fraunhofer (b) diffraction of light
Figure 8. Fresnel (a) and Fraunhofer (b) diffraction of light
View of the Malvern 2600 HSD
Figure 9. View of the Malvern 2600 HSD
Influence of light obscuration
Figure 10. Influence of light obscuration
Phase difference in scattered light from a spherical particle
Figure 11. Phase difference in scattered light from a spherical particle
Light scattering modes from a spherical particle
Figure 12. Light scattering modes from a spherical particle
Removing the 2π ambiguity
Figure 13. Removing the 2π ambiguity
a) Effect of collection angle on droplet size; (b) Effect of collection angle on mass flux
Figure 14. a) Effect of collection angle on droplet size; (b) Effect of collection angle on mass flux
Effect of collection angle upon the droplet probability density function
Figure 15. Effect of collection angle upon the droplet probability density function
High speed images of a puff taken from a Salbutamol HFA-134a pMDI (Ventolin®). These images were treated with greyscale and inverted colours after application of a threshold filter for easier visualisation of the plume
Figure 16. High speed images of a puff taken from a Salbutamol HFA-134a pMDI (Ventolin®). These images were treated with greyscale and inverted colours after application of a threshold filter for easier visualisation of the plume
Basic geometry of an impactor
Figure 17. Basic geometry of an impactor
A "testbox" representation constituted, the red plane (A) is the boundary condition ‘Velocity Inlet’ and the green plane (B) is the boundary ‘Outflow’. The other four outer walls of the domain have ‘Symmetry’ properties
Figure 18. A "testbox" representation constituted, the red plane (A) is the boundary condition ‘Velocity Inlet’ and the green plane (B) is the boundary ‘Outflow’. The other four outer walls of the domain have ‘Symmetry’ properties
Velocity magnitude contours of the air in the domain, along the longitudinal midsection plane (t=0.11 s)
Figure 19. Velocity magnitude contours of the air in the domain, along the longitudinal midsection plane (t=0.11 s)
Representation of the particle streams at the end of the injection (t=0.11 s), image shows the particles coloured by residence time inside the domain. The particle streams are draw as spheres with proportional size scaled 50 times more than the real diameter
Figure 20. Representation of the particle streams at the end of the injection (t=0.11 s), image shows the particles coloured by residence time inside the domain. The particle streams are draw as spheres with proportional size scaled 50 times more than the real diameter
Representation of the particle streams at the end of the injection (t=0.11 s), image shows the particles coloured by its velocity magnitude. The particle streams are draw as spheres with proportional size scaled 50 times more than the real diameter
Figure 21. Representation of the particle streams at the end of the injection (t=0.11 s), image shows the particles coloured by its velocity magnitude. The particle streams are draw as spheres with proportional size scaled 50 times more than the real diameter

pMDI Sprays: Theory, Experiment and Numerical Simulation

Ricardo F. Oliveira1, Teixeira Senhorinha2, José C. Teixeira1, Luís F. Silva1 and Henedina Antunes3, 4, 5

1. Introduction

The inhalation therapy is a cornerstone in asthma treatment. A disease characterized by episodes of wheezing, breathing difficulties, chest tightness and coughing. Essentially it is a chronic inflammatory disorder associated with airway hyper responsiveness [1]. This disease already affects over 300 million people worldwide, growing at a rate of 50% per decade, and causes the death of 220 thousand per year [1].

Therefore it is of utmost importance to improve the delivery effectiveness of the devices used in the inhalation therapy, which is the most used way of treatment, thanks to its greater efficiency, faster action, lower toxicity and minor drug waste [2].

Anti-inflammatory and bronchodilator drugs are used with the objective of reducing the inflammation of the pulmonary tissue, which cause the diameter reduction of the bronchus [1,3].

In the inhalation therapy procedure, the following devices are mostly used: nebulizers, Dry Powder Inhalers (DPI), and pressurized Metered-Dose Inhalers (pMDI), which can include a spacer attached (add-on device).

The nebulizer presents a simpler usage procedure but most of models available need to be attached to the power plug. They also present a greater waste of drug and a slower duration for each session of treatment and a lower efficiency when compared with other devices. The DPI is the most efficient device but also the most expensive, requiring a minimum inspiratory flux by the patient in order to work properly which is not expected in young children and elderly. Regarding the pMDI, it has a low cost and it is the most used by the medical community and its efficiency is acceptable. Its major drawbacks are the high spray velocities, which creates the so called “cold-Freon” sensation on the back of the throat, and the necessity of inspiratory coordination with the priming action, resulting that the younger and/or older patients may not be able to use the device properly. This problem is highlighted by the fact that the implementation of classes to teach the correct usage technique provide increased efficiency [46].

With the objective of solving some of the pMDI problems, the spacers were created and they are presented into three categories: simple tube, valved holding chamber and reverse flow. The most common type is the valved holding chamber, inside which the drug spray is released. This device leads to the reduction of the high velocity impact of the spray in the throat and the need for inspiratory coordination, allowing the patient to breath normally from the other side of the chamber. The one-way valve makes the normal respiratory cycle possible by not allowing the expiration flux to go back inside the chamber [58].

1.1. History of inhalation therapy

According to Crompton, inhalation therapy can be traced as far back as 4000 years ago in India, and many ‘hundreds of ingenious devices and hopeful medications’ (according to Sanders) have been discovered so far. Inhalation therapy is so important nowadays that is hard to imagine the treatment of asthma without a proper inhalation therapy procedure, inhalation device, and a necessary drug [9,10].

Briefly reviewing the antecedents of contemporary inhalation therapy and the therapies and medicines used for the management of asthma over the last 50 years, Crompton presents the 19th century hand-held glass bulb nebuliser up to the introduction of the first pMDI, which came into common clinical practice in 1956 [9]. Sanders, detailing the inhalation therapies since ancient times, refers several inhaler devices, some of which are the C. Bennet’s inhaler (dated from 1654), the J. Mudge’s inhaler (from 1778, and adapted from a pewter tankard), the first pressurised inhaler (presented in 1858, in Paris, France, by Sales-Girons), the Improved Nelson inhaler (proposed in 1865 by S. Maw & Sons, London) which is still being manufactured up to this day, and the Siegle’s steam spray inhaler (a German invention from early 1860’s that marked the beginning of the nebuliser therapy) [10].

Many contributions have been given so far, and many medications are available nowadays as inhaled treatments, as well as many inhaler devices (Rotahaler® and Accuhaler® are just two examples).

1.2. pMDI devices

The pMDI, shown in Figure 1, is the most commonly used aerosol delivery device in the world for the treatment of asthma. In 2000, an estimated production of 800 million units have been reported [5,11]. It has been the backbone of inhalation therapy for asthma for approximately 50 years [6]. The pMDI contains 100 - 400 doses in a small and (very) portable device that can be easily concealed in a pocket. The best feature of a pMDI device is that is always ready to be used, making it a magnificent example of an engineering solution.


Figure 1.

Three commercial pMDI actuators with its drug canisters attached

1.2.1. Components and their role

The pMDI device essentially incorporates a disposable canister, where the drug formulation is stored, which can be replaced for a new one at any time, mounted in an actuator with a mouthpiece zone and usually a dust cap is included [6,12]. The four basic components that can be found in all pMDIs are: the formulation canister; the metering valve; the actuator; and the container [13]. On the following list are described the components of a hydrofluoroalkane (HFA) pMDI device, published in [5]:

  • Canister with O-Ring

  • Formulation

  • Drug substance

  • Propellant(s)

  • Surfactant(s)

  • Lubricant(s)

  • Co-Solvent(s)

  • Ferrule Gasket

  • Complete Valve

  • Mounting cup

  • Valve core assembly

  • Diaphragm

  • Stem
  • Spring

  • Seal
  • Metering tank

  • Retaining cup

  • Adaptor with dust cap

The spray pattern and Mass Median Diameter (MMD) are influenced by the ambient temperature, the design of the actuator nozzle and valve stem, also by the vapour pressure of the propellant in the formulation [5].

The formulation is composed by the drug, the propellants gases and it often contains surfactant and other excipients. In the formulation, the propellant is the component which creates vapour pressure inside the canister, allowing the drug to exit the pMDI and form a spray plume every time it is pressed. When the solid drug particles exit the pMDI they are encapsulated in a droplet made of propellant, which evaporates as it travels through air. Some of the components in the formulation of a common used drug in asthma treatment, Ventolin®, are listed in Table 1.

Ventolin® HFAVentolin® CFC
Active ingredientAlbuterol sulphateAlbuterol sulfate
ExcipientNoneOleic acid
Formulation typeSuspensionSuspension

Table 1.

Ventolin®’s formulation components present is two propellant types, HFA and chlorofluorocarbon (CFC). Adapted from [12]

The canisters for a pMDI are typically made of aluminium and they are designed to be light, compact and strong to hold the high internal vapour pressure of 3 to 5 bar made by the propellant in the formulation [6,12]. Its regular capacity ranges 15-30 mL [13].

The metering valve of a pMDI is crimped into the container. There are several designs of valves but they all work under de same principle: while the canister is in the inverted position, by gravity action, the valve fills with the formulation through a channel archiving the desired measured dose; at the priming moment, this channel closes and another one opens in the opposite side of the metering valve, allowing the formulation to rapidly expand into the expansion chamber. The actuator pit and the valve stem constitute the expansion chamber where the propellant begins to boil [12,13]. This process is often referred as the primary atomisation, where a flash evaporation of the propellant takes place [13,14]. After it exits the expansion chamber through the actuator nozzle it rapidly expands in the shape of a solid-cone plume [12,13]. In this way, the metering valve of a pMDI is a critical component in the effectiveness of the delivery system, due to the fact that its main functions are [46,1113]:

  • deliver, in an accurately and reproducibly way, a measured volume (20-100 µL), containing between 20 - 5000 µg of the dispersed drug;

  • form a propellant-tight seal for high pressure in the canister.

The actuator is normally a single plastic piece produced by injection moulding that consists of a mouthpiece, body and nozzle [6,11–13,15]. When the pMDI reached the market in 1956, it was composed of an elongated mouthpiece (around 8 cm); nowadays, it was reduced (from 2 to 3 cm) to a more compact size to improve the portability [12].

  • the mouthpiece is the interface part to the patient mouth;

  • the body provides support for the canister;

  • the nozzle has a very important role in controlling the atomisation process, to guarantee a spray plume formation. The nozzle diameter interferes directly with the particles size distribution.

1.2.2. Formulation: CFC versus HFA

On the 16th September of 1987, the Montreal Protocol was signed, with the objective to drastically reduce the CFC emission to atmosphere which was contributing to the destruction of the ozone layer. Being the CFC propellant the main component in the pMDI formulation, it was clear that some big changes were required to the pharmaceutical industry. The chosen propellants for this transition were the HFA. Although these greenhouse gases contribute to the global warming, their effect at this scale is considered negligible [6]. During this transition not every single drug type available with CFC propellant was transferred to the HFA formulation, reducing the types of drugs available to be prescribed [3].


Figure 2.

Chemical structures of the two most common types of propellants (CFC-12 and HFA-134a) used in pMDI

The CFC gas meets the criteria to be used as pMDI propellant, such as, the constant vapour pressure during the product’s life ensuring consistent dose; needs to be nontoxic, non-flammable; correct boiling point and densities and also needs to be compatible with the drug solution. These reasons made CFC the propellant of choice, it is formulated using the CFC-12 (Figure 2) as the major component, added with CFC-11 and/or CFC-114 to regulate the vapour pressure of the formulation [12]. Its substitute, the HFA-134a (Figure 2) or HFA-227, resulted in a reformulation of all the pMDI components due to the serious changes in the spray plume shape, particle size distribution and amount of drug dose delivered [7,11,16,17]. These HFA propellants forced the redesign of the actuator and valve components, essentially the gaskets [7]. The HFA-134a has similar thermodynamic properties to the CFC-12, but there is not a substitute HFA component to the CFC-11 or CFC-114. So there is the need to introduce excipients with lower volatility in the mixture to modify the vapour pressure, contributing to the development of new methods for the pressure-filling process [12]. The transition progressed in two different ways. Some companies decided to develop bioequivalent products using the HFA propellant, while others decided to take the opportunity to make a new product which is clinically more efficient in the delivery to the patient lungs [6,12].

Table 2 describes some of the most important physicochemical properties of the two main pMDI propellants, CFC-12 and HFA-134a. The HFA propellant products require a higher vapour pressure to reach the liquid state in comparison to the CFC, resulting in a higher speed spray [7,17]. This forced the redesign of the metering valves so they could deliver the same amount of formulation volume per puff, leading to changes in the particle size distribution, because the HFA formulation produces smaller particles [5,7].

PropertyCFC 12HFA 134a
Boiling point (ºC)-30-26
Vapor pressure (kPa)566572
Density (g/cm3)1.331.23
Viscosity (mPa.s)0.200.21
Surface tension (mN/m2)98
Dielectric constant2.19.5
Water solubility (ppm)1202200

Table 2.

Some physicochemical properties for two main pMDI propellants (CFC-12 and HFA-134a). Adapted from [11]

1.2.3. Mechanism and technique

The pMDI device contains the drug formulation inside the canister at a high pressure, which at moment of priming flows rapidly into the metering chamber and to the expansion chamber of the actuator, and then throughout the nozzle. At this point, the drug particles became airborne forming the spray plume, travelling through the air. The sudden pressure change makes the liquid propellant of the formulation to rapidly evaporate, going from a high diameter droplet into a smaller drug particle and/or droplet.

The pressing, also called priming, of the top of the metal canister against the plastic actuator, makes the spray plume discharge, also known as puff, into the air. This priming needs to follow some procedure steps to result in an efficient delivery to the patient. These steps are normally the following [18]:

  1. Shake inhaler

  2. Exhale

  3. Breathe in deeply and slowly

  4. Press on inhaler (co-ordinating actuation with inhalation)

  5. Hold breath

  6. Breath out slowly

  7. One actuation per inhalation

  8. Wait before second actuation

By studying the technique of the pMDI used by the patients, it was found that the majority of the patients have difficulty to execute all the steps correctly. The study included children (n= 26), adolescent (n= 23), adult (n= 209) and elderly (n= 66) patients. Among them, the children are the ones with more difficulty in the step 4, and the elderly have more difficulties executing the steps 7 and 8. The adolescents and adults are the patients with better abilities to use the pMDI device correctly [18].

2. Spray characterisation

In this section, all the parameters and variables needed to describe the spray will be presented and discussed, which will deal with its characterization, spray plume formation and position, droplet size distributions and mass flux.

2.1. Spray characteristics

The production of aerosol in a pMDI is a complicated process, which is influenced by several factors. The vapour pressure is the most important factor in the pMDI mechanics, and it influences the distribution and shape of the plume. Higher vapour pressure inside the canister results in a faster velocity plume and with faster evaporation rate [11]. Clark developed a mathematical equation based on empirical results from CFC formulation, although it performs well for HFA propellants too. Equation 1 allows the calculation of the MMD produced by the pMDI [11,13,14].


where the MMD is represented asD0.5, the qec is the quality of the flow from the Rosin-Rammler distribution, pecis the pressure in the expansion chamber and the p represents the ambient pressure.

The interaction of the formulation with the other pMDI components defines the final spray plume form and characteristics. An relationship can be used to predict the percentage of fine-particle fraction in a HFA-134a spray, see Eq. 2, by using the dimensions of its components, as described in [12].


Where, Arepresents the actuator nozzle diameter (in mm), Vis metered volume size (in µL) and C134 is the HFA content (in percent) [12].

2.2. Velocity

Within the pMDI spray plume, the droplet velocity varies accordingly to the position in it. In the middle of the plume cone the velocity is maximum, while in the outer zones it gets its lower values (see Fig. 3).

The droplet velocity also decelerates along the axial distance from the nozzle, due to momentum exchange with the air. As reported by Dunbar, using phase-Doppler particle analysis measurements of a HFA-134a spray plume during an actuation, values were taken at different distances from the nozzle of the actuator (see Fig. 4) [13].


Figure 3.

Droplet velocity vectors represented in half of a mid-section plane of pMDI plume using a formulation only with HFA-134a, measured at the last moment of the actuation. Adapted from [13]


Figure 4.

Droplet mean velocity along the axial distance from the nozzle of a pMDI using a formulation only with HFA-134a. Adapted from [13]

According to the measurements made by Dunbar, the author concluded that a HFA propellant formulation produces a spray with higher velocities than a CFC formulation. This fact is attributed to the higher vapour pressure used in the HFA formulation. The plume behaves like a spray up to the distance of 75 mm from the nozzle and as an aerosol afterwards that distance, where the droplet motion is being influenced by the gas [13].

2.3. Droplet size distributions

As the spray is formed downstream of the nozzle of the plastic actuator, it has suffered the influence of the drug formulation, valve and actuator design dimensions. These factors create a resulting plume of spray with a characteristic particle size distribution. The spray formation process causes an interesting variation of the MMD of the dose, having a value of 35 µm at the nozzle which reduces to 14 µm at just 10 cm from it, due to evaporation process of the propellant [19].

Particle size distributions can be represented in two forms, as a Probability Density Function (PDF) or as a Cumulative Distribution Function (CDF). Normally, an independent variable (x) and two more adjustable parameters, representing the particle size and the distribution of particle sizes are used [20]. There are several types of mathematical distributions to describe spray particles/droplets, being the Log-Normal, Rosin-Rammler and Nukiyama-Tanasawa, the most common.

It is a common idea that pharmaceutical aerosols can be accurately described by the Log-Normal function fitting the measured data to the cumulative mass distribution (see Eq. 4). The Log-Normal PDF (see Eq. 3) was developed from the normal distribution [20].


where σg is the geometric standard deviation (which shall be ≠ 0), the μ represents the mean diameter and erfc is the complementary error function.

Using data from coal powder, an empirical function was developed to describe the distribution, this function is called Rosin-Rammler, also known as Weibull distribution (see Eq. 5). It has been widely applied in the atomisation field, using an independent variable (x) and other two to describe the distribution, λrepresents the mean diameter and k is the distribution spread. The CDF is simply given by Eq. 6, making it easy for graphical representation [20,21].


In 1939, Nukiyama and Tanasawa developed a function for characterisation of twin-fluid atomizers (see Eq. 7). It makes use of the gamma distribution function, Γ(x), and two other parameters: b that represents the size and δ that is the distribution parameter. In the CDF form it is described by Eq. 8 [20,21].


Dunbar and Hickey studied the best distribution and fitting method for empirical data obtained from a pMDI spray using an Andersen eight-stage cascade impactor. They concluded that the best fitting method is through the nonlinear least squares of the PDF and that the Log-Normal and Nukiyama-Tanasawa PDFs produced better fitting results to the data than the Rosin-Rammler function [20].

The authors measured a pMDI HFA-134a formulation of salbutamol using the laser diffraction analysis technique (throughout a particle sizer Malvern 2600), using an independent model for data processing. The data were fitted to the three distribution models described above (see Eqs. 4, 6 and 8) using the least squares method. A comparison is shown in Figure 5.


Figure 5.

Graphical representation of the pMDI HFA-134a salbutamol experimental data and its fitting for the Rosin-Rammler, Log-Normal and Nukiyama-Tanasawa distributions. Measurements obtained at 100 mm from the laser beam

2.4. Mass flux

Mass (or volume) fluxes are an important characteristic of any spray. In particular for drug delivery applications, such knowledge is of utmost importance to the understanding of droplet/particle transport through the respiratory system and any delivery devices.

Space averaged measuring systems (such as laser diffraction based) only provide a crude estimative of droplet/particle concentration. The user has to specify the sample length illuminated by the laser beam and the concentration (estimated for the duration of the measurement throughout the entire sample) is calculated as the ratio of the light intensity collected by the sample over the total light intensity existing in the absence of any sample. One must remember that the light array detector includes a central detector where the undiffracted light is focused.

Because the technique only measures particle size, the procedure does not yield data for mass flux. Furthermore, the information is a time-space averaged quantity which is of limited value, particularly for transient phenomena such as that occurring during the operation of inhalator devices.

Phase Doppler anemometry may overcome such limitations, because it is a point measuring procedure and, crucially, includes data for velocity. In addition, because the sampling occurs over a period of time, one may capture time variations in properties such as the droplet mass flux.

The mass flux results from the integration of the total droplet volume (mass) over a certain area for a period of time. Various sources of errors make this calculation difficult. Crucial to this procedure is the correct determination of the effective area. The probe volume results from the intersection of two Gaussian laser beams (of diameter de) which has an ellipsoidal shape. The actual probe volume is the intersection of this ellipsoid with the slit length (L) in the direction of the receiving optics which is oriented at an angle from the forward direction. A simple approximation for the cross area is the projected area (deL) of the volume described above. However one should take into account the cross section area of the volume perpendicular to the particle velocity. This only coincides with the former in a one dimensional flow if the main axis is oriented with this velocity. Therefore one must take into account the actual velocity vector of each particle which is problematic because most of the phase Doppler anemometry configurations are 2D. Zhang and Ziada refer a method for overcoming this limitation [22].

In addition, the droplets used for the calculation of the mass flux may not cross the probe volume at the same rate throughout its area. Therefore, the droplets used for this calculation are allocated to a certain position in the probe volume by taking into account the Doppler burst length. The calculation of the number of droplets in a certain time is based on the number of droplets during the elapsed time of the measurement. It is argued that this is not accurate because in that time more particles have been through the probe volume than those actually accounted for and this ratio depends upon the validation rate of the experiment. The attempted/validated samples ratio should correct this estimate. Furthermore, other sources of uncertainty are caused by the fact that not all the particles may trigger the system. This may result in a bias against the contribution of small particles.

Because the particle counting is based on average values recorded in size bins, other sources of uncertainty may come from high turbulence and insufficient number of particles per size bin. This may be a serious drawback because bins for large particles are most likely to have small number of droplets and, on the other hand these droplets are those which most contribute to the total mass/volume fluxes. Concerns regarding this estimative have been raised by other authors [23].

Figure 6 shows an example of the limitations and uncertainties in the calculation of a mass flux. The data reports the variation of the local liquid mass fluxes of a twin fluid atomizer for various gas flow rates (the liquid flow rate is constant at 1 L/min), which shows a considerable scatter for the various Air to Liquid Ratio (ALR).

Nevertheless, the integration of the liquid fluxes over the diameter should yield a constant value for all the cases. As shown in Fig. 7, the results (normalized) may differ by a factor of 5. It should be stated that some factors may contribute to such differences: very high turbulence levels, dense sprays (particularly at high ALR's).

Figure 6.

Mass flux as a function of ALR

Figure 7.

Integration of the mass fluxes at constant liquid flow rate

In conclusion, this information must always be treated with caution and may only provide an indication of the relative value of the droplet mass fluxes.

3. Experimental techniques

The drug delivery through the oral system may be put in a simple manner as the transport of a dispersed phase through the upper and lower respiratory tracts. In addition it may include additional systems such as a spacer. Therefore it is widely accepted that particle size plays an important role in determining where the aerosol particles are deposited. The respiratory system acts as a filter that progressively filters particles of smaller size as the inhaled air flow reaches the lower tracts of the lungs. Particles with an aerodynamic diameter above 6 μm will deposit in the oropharyngeal region while those below in the range 2-4 μm will reach the alveoli being the most effective as drug delivery agents. In short, particle size is of the greatest relevance for accessing the effectiveness of the drug delivery systems.

In a brief sentence, droplet sizing is one of the most sensitive parameter to measure in a wide range of dispersed multiphase flows. The problem is compounded if the flow is confined by walls.

Over the last decades, a variety of techniques have been developed in order to overcome some basic limitations: the required accuracy for measuring small objects; the statistical requirement of very large samples; the fact that droplets may be fast moving in the flow stream; the spray may be dense; possible sampling interference.

Over the years, various literature reviews have been published. As an example, one should refer the work of Black et al, Hong et al and Mitchell and Nagel [24–26].

3.1. Laser diffraction techniques

Coherent light (wavelength λ) diffracted by a particle can be used to measure its size. Two kinds of diffraction must be considered. In Figure 8 (a) light diffracted by an angle θ interferes with the undiffracted light, forming an interference pattern on the screen. This is near field or Fresnel diffraction. However, it can be observed that particles with the same size but positioned at different places along the light axis (diffracting by the same angle) will give a different radial position on the detector, which can cause confusion.

Figure 8.

Fresnel (a) and Fraunhofer (b) diffraction of light

Alternatively, if the detector plane is moved away from the particle (ideally at an infinite distance) this problem will not occur. This is known as the far field or Fraunhofer diffraction. In practice this can be achieved by placing a lens in the light path so that the undiffracted light will be focused onto a central spot F and light diffracted by an angle θ is focused at the point D. The distance DF is related to the diffraction angle θ by the expressionDFfθ, where f is the lens focal length – Figure 8 (b). In the Fraunhofer approximation that only describes diffraction of light at the contour of the particle, it is assumed that all particles are much larger than λ and that only near forward scattering is considered (i.e. θ is small).

In the method proposed by Shifrin the light intensity pattern on the focal plane is measured at different radial positions by means of a photo multiplier, moving this unit along the focal plane [27]. The minimum particle diameter measurable is proportional to the light wavelength by a factor of 10/π. This method is more suitable for measuring the diffracted light at small angles.

The method developed by Swithenbank et al relates the measured light energy at the lens focal plane to the particle size [28]. This technique is commercially available in an instrument manufactured and distributed by Malvern Instruments Ltd. (England). The technology employed is based on the Fraunhofer diffraction of a parallel beam of monochromatic light by a moving particle. The diffraction pattern produced consists of a series of alternate light and dark concentric rings, the spacing of which is dependent on the drop diameter. This results in a series of overlapping diffraction rings, each of which is associated with a characteristic particle size range, depending on the focal length of the Fourier transform lens. This lens focuses the diffraction pattern onto a photo detector that measures the light energy distribution. Typically, the photo detector comprises 31 semi-circular rings surrounding a central circle. Each ring is therefore most sensitive to a specific size of droplets. The output of the photo detector is multiplexed through an A/D converter, which is connected to a computer. The computer provides an instant display of the measured distribution of drop sizes. Thus, a major advantage of the instrument is the speed at which data can be both accumulated and analysed. Another important attribute is that the diffraction pattern is independent of the position of the drop in the light beam. This allows measurement of size distributions to be made with drops moving at any speed. Lefebvre presents a very complete description of the most important characteristics of the instrument. Figure 9 shows such instrumentation [29].


Figure 9.

View of the Malvern 2600 HSD

The relationship between the size distribution [D] and the light intensity distribution [I] in each ring is of the form:


in which [A] is a square matrix whose coefficients depend upon the optical configuration of the instrument and of the sensitivity of each one of the photodetectors. By selecting an appropriate focusing lens, the sizing range can be adjusted. Typically the dynamic range is of approximately 100:1. For instance, by selecting a 63 mm lens the measuring range will fit in between 0.5 and 100 μm. Because of the nature of the coefficients, equation 9 is solved by fitting a probability distribution (see section 2) to the data. In addition, the occurrence of multiple diffraction due to the spray denseness, may result in a bias towards smaller diameters (outer rings).

Figure 10 shows such influence upon the droplet distribution. It is accepted that only if the fraction of light collected by the sample exceeds 50% of the incident, there is cause for concern. This is rarely the case for the aerosols such as those produced by pMDIs.


Figure 10.

Influence of light obscuration

3.2. Phase doppler anemometry

Phase Doppler Anemometry (PDA) is an extension of the Laser Doppler Anemometry (LDA) principle, capable of measuring simultaneously the size and velocity of individual droplets. The velocity measurement principle is the same used for LDA, making use of one detector for each velocity component to be measured. However, to perform size measurements two detectors are needed, in principle, as described below.

Figure 11 represents light scattered from a spherical particle, collected by two detectors. Since the detectors are at different positions, the optical path lengths for reflection from the two incident beams are not equal. Consequently, when a particle crosses the probe volume both detectors receive a Doppler burst of the same frequency, but with a phase difference.

The phase difference between the two Doppler bursts depends on the size, or its surface curvature, of the particle. If Δt is the time lag separating the wave fronts reaching the two detectors, the corresponding phase difference isΦ12=2πfΔt. For a spherical particle, the phase difference increase with its size.

The relationship between size and phase of a particular optical arrangement is described by the Mie scattering theory. However, a geometrical optics approximation can provide a good basis for analysis.


Figure 11.

Phase difference in scattered light from a spherical particle

The phase of the Doppler burst received at detector i is expressed as:


where n1 is the refractive index of the scattering medium and βi is a geometrical factor, which depends on the scattering mode and the collection arrangement.

As shown in Figure 12, when an incident light strikes a particle, several modes of scattering may occur. The light may be reflected or suffer 1st, 2nd or nth order of diffraction. The phase difference on the detector will depend upon the dominant scattering mode as well as the collection arrangement.


Figure 12.

Light scattering modes from a spherical particle

The angle of intersection between the two beams, θ, is determined by the beam separation, St, and the lens focal length. This angle determines the fringe separation. The collection angle, φ, and the elevation angle, ψ, define the direction towards the photo-detectors from the measuring volume.

For the reflection mode the geometrical factor is given by:


This expression shows that, for a particular optical arrangement, βi is a constant factor, which gives a linear relationship between phase and size. Furthermore, the refractive index of the particle, n2, does not appear, making this scattering mode useful in situations where the value of the refractive index is not known.

On the other hand for 1st order refraction it takes the form:



n2 = particle refractive index


For the other modes of scattering similar relationships may be derived [30]. As expressed in the equations for βi above, changing any of the angles θ, φ or , will affect the geometrical factor and, consequently, the sensitivity and range of the PDA.

In practice, there are some restrictions in the selection of the geometrical optical parameters. For instance, the selection of scattering angle, φ, is quite restricted, either to ensure a specific scattering mode or a sufficient signal-to-noise ratio, or from practical considerations of the measurement situation (optical access and working distance). The required working distance also affects the possible range of θ and .

The scattering angle, φ, and the elevation angle, ψ, are those with the most important effect upon the parameter βi. Once the elevation angle is selected (usually 0°), the scattering angle is crucial in defining the optical set-up. As explained before, although light scattering is fully described by the Mie theory, the PDA instrumentation relies upon the simpler approach based on the geometric optics principles. The question here is to decide about the reliability of this approach. Regarding this, the scattering angle plays an important role.

In a measurement situation all scattering modes are present. Depending on the collection angle, some may be dominant relatively to the others or may be of equal importance. A phase Doppler system should be set up so that the relationship between particle diameter and phase difference is linear. In general, this is obtained by choosing an angle where one single mode dominates the scattered light received by the receiving optics, and where the signal–to-noise ratio is as high as possible. Tayali and Bates have shown that, depending on the refractive index of both the medium and the particles, the scattering angle should be such that only one scattering mode is clearly dominant [31]. In this way the geometric optics principles will approach the exact theory and processing errors will be smaller. This shows that one condition for a reliable PDA set-up is the knowledge of the particle and medium refractive index.

Finally, if two detectors are present, the phase difference is calculated by subtracting the phase between each of the detectors (see Eq. 10). The phase difference of the larger particles may fall beyond the 2π range. Thus, there is no way one can tell whether the diameter is D3 or D3’. This is called the 2π ambiguity.

Therefore, with two photo-detectors, there is a compromise between, on the one hand, high sensitivity and small working range and, on the other hand, a larger working range at the expense of the sensitivity.

This problem is overcome by using three detectors asymmetrically positioned, as shown in Figure 13, and measuring more than one phase difference. Therefore, the phase difference between the two closer will extend the working range and the phase difference from the detectors further apart will give the necessary resolution.


Figure 13.

Removing the 2π ambiguity

This solution has another useful feature. The phase difference corresponding to each pair of detectors gives information about the curvature of the particle surface. This can be used to deal with non-spherical particles, either by measuring their sphericity ratio or by establishing a criterion for sphericity acceptance.

This technique was applied, by Liu et al, specifically to pMDI for measuring the droplets instantaneous velocity and mean diameter of the spray distribution [32]. A total of nine commercial pMDI products were tested, showing different droplet velocities for each one. It is also noticed the decrease in the spray droplets along the axial distance from the pMDI nozzle, as previously reported by Dunbar [13].

3.3. Uncertainties in PDA measurements

In applying a technique such as the PDA, one must be aware of the problems and potential limitations; errors may result of inappropriate configuration of the instrument. The major question concerns the fact that any sample will observe a limited slice window of the population; the question balances the various compromises necessary to carry out any experimentation.

In this topic one looks into the influence of PDA configuration and set up on the measured data. Particular attention is given to parameters that do affect the measurements and discussion is provided into these factors based on experimental data. For that purpose, the tests were carried out at constant gas/liquid flow rates in one nozzle. In addition, all the tests were carried out at the centre of the spray. Various parameters were slightly changed ant their effect upon the mean diameters, data rates and droplet fluxes are discussed. Insight is given into the best approach to obtain consistently reliable results.

Figure 14 (a) and (b) show the influence of the influence of the collection angle in the mean diameter and the mass flux, respectively. It must be stressed that in the case of a more complex flow (with walls for example) the difficulties should be greater and the results would show higher deviations. As a general conclusion, it can be concluded that considerable spread can be obtained unintentionally: a factor of two in diameter and a factor of six in the mass flux is possible. The various possible causes for such observations are discussed below.


Figure 14.

a) Effect of collection angle on droplet size; (b) Effect of collection angle on mass flux

It has already been mentioned that the collection angle is by far the single most important factor in deciding an optical set up. However, in certain conditions, other angles could be more attractive and have been extensively used in the past. For example, if windows are necessary (Internal Combustion engines as an example) one single window would be very helpful in the rig construction. Therefore an angle of 150º (close to back scatter) is attractive. Other angles have also been employed, being the 30º one of the most common.

Results from a 10 μm water droplet show that reflection goes to zero at 72º (Brewster angle) and that 1st order refraction is over three orders of magnitude greater than any other non-zero mode. Thus, at this angle the geometric optics algorithm should not be affected by the presence of other scattering modes as they are of negligible intensity.

However, the detailed analysis of the droplet size PDFs also shows other interesting features, see Fig. 15. Firstly, it is obvious that the PDF at 150º is totally absent of small droplets. The problem is further compounded by the fact that, at that angle, the intensity of the scattered light is much lower than in the forward direction and therefore small droplets may not trigger the processors.

Although the mean diameters at 30 and 72º are close to each other, a close inspection in Figure 15 shows another story. Firstly, by altering the collection angle, the sizing range has dropped from 195 μm down to 134 μm which means that the distribution measured at 30º is truncated. If one removes the individual droplets above 134 μm in the 72º experiments, the mean diameter (d32 in particular) drops by approximately 10%! In Fig. 15 (c), the width of each bin in the distribution is made equal to that of the 72º collection angle. Therefore the counting in each bin can be directly compared. Although the 30º measurements should see more small droplets (close to forward collection where the signal intensity is higher which is supported by the higher data rates - 140 kHz as opposed to 84.5 kHz at 72º) the number of droplets at the lower end of the range is less than at 72º. This is clearly because at that angle the approximate solution is less accurate than at 72º. As expected, the data rates at 150º are much lower (33 kHz). However the mass flux shows that at this angle of collection, higher fluxes are observed (see Fig. 14 (b)). This is certainly due to the fact that the distributions contains a higher proportion of larger droplets and the flux (either on a mass or volume basis) are weighted by those large droplets.


Figure 15.

Effect of collection angle upon the droplet probability density function

3.4. Photography and high-speed camera

Photography can be used for particle sizing through its application to small and fast moving particles requires extreme care in order to ensure an acceptable level of accuracy. Because of the particles motion, a very short exposure time is required to “freeze” the image. This requires both a high level of illumination intensity and a wide lens aperture. This means that the depth of field is small and only a fraction of the total sample is likely to be in focus at any moment.

For smaller particles, the turbulence induced motion in a direction perpendicular to the focusing plane compounds the problem because they are prone to be out of focus. One technique proposed to overcome this problem consisted in producing multiple focusing planes that overlap each other. Another major hurdle is due to the limited amount of data that can be experimentally retrieved. The measuring process is very time consuming rendering the effective amount of data of limited statistical value. Furthermore, although the sizing range is high, the lower end of the limit (of approximately 5 μm) makes the technique of limited interest in drug aerosols measurements.

This technique can also be extended to measure the particle velocity by increasing the exposure time in such a way that the path described by the particle during the lens aperture is correlated with the particle velocity. Additionally, a double flash (with a controlled time delay) can be deployed to superimpose two photographs in the frame. This technique has the advantage of a more accurate determination of particle size.

Nonetheless the development of other techniques that provide better resolution (such as the phase Doppler) has rendered this technique of limited use for the last decades.

However the advent of high speed digital video camera can successfully be used to gather meaningful insight into spray dynamics. The problems associated with sample illumination are ever present but this technique, which can record up to 10,000 frames per second, is very useful for understanding and capturing details of transient phenomena. In particular for aerosol delivery of drugs the very nature of the spray resulting from a fast operating canister suits the potential of such technique. When combined with a plane illumination from a light sheet the particles can be easily identified and tracked to measure their velocity. Also, the frame high resolution coupled with filtering algorithms for identifying the contours of particles, are useful for measuring their size. However the most important feature of this technique is its ability to capture the transient nature of the aerosol formation over the delivery time.

Using a high-speed camera (FASTCAM-APX RS 250KC), the authors recorded a puff event from a Ventolin® HFA-134a pMDI. Results are shown in Figure 16, the images presented have an interval of 0.02 seconds. The images were taken at a rate of 6000 frames per second, allowing us to confirm the duration of the spray (0.1 seconds) and to calculate the approximated angle (10 degrees).

3.5. Impacting techniques and collectors

Impacting collectors have been used to measure particle size for a wide range of sprays. In particular for inhaler derived aerosols, it has been the most widely used method. Mitchell and Nagel have extensively reviewed this type of particle size analysers, their strengths and limitations [33].


Figure 16.

High speed images of a puff taken from a Salbutamol HFA-134a pMDI (Ventolin®). These images were treated with greyscale and inverted colours after application of a threshold filter for easier visualisation of the plume

The theoretical principles of an impactor are based on the solution of the two dimensional Navier Stokes equations of a fluid around a target oriented perpendicularly to the impinging flow. This flow field is subsequently coupled with Newton’s equation of motion to model the trajectories of particles/droplets through various stages and geometries.

In its basic form, the impactor consists of a jet of a known diameter that exits from an orifice diameter (D) located at a certain distance (S) from the flat target that acts as a collector (see Fig. 17).


Figure 17.

Basic geometry of an impactor

As the flow exits the orifice, its streamlines are diverged in the vicinity of the collection surface. The high inertia particles move through the streamlines and impact on the collection surface. The Stokes number (St) defines the critical particle size (aerodynamic diameter da) that will impact the collection surface through


where Ca is the Cunningham slip, ρa is the particle density, U is the gas velocity at the orifice exit and μ is the fluid viscosity.

The particle collection efficiency of an ideal impactor will increase in a step from 0 to 100% at a critical St. For each impactor stage, the corresponding cut-off diameter (d50) can be calculated through equation 14.


where n accounts for the number of circular orifices in the nozzle plate, Q refers to the flow rate and subscript 50 identify conditions at 50% efficiency.

A cascade impactor consists of several stages with progressively decreasing cut sizes, assembled in series, so that an aerosol is separated into various fractions according to the number of stages assembled. The design of the cascade can be tailored to measure the size down to sub-micron particles. The two stage impactor is used to separate the large from the small particles. The former are deposited in the upper tract of the respiratory system while the later are most likely to reach the alveoli. For inhaler testing, it is desirable that at least 5 stages are assembled with cut off aerodynamic diameters within the range 0.5-5.0 μm. One of the most widely used is the Anderson type impactor with 8 stages. The testing of aerosols usually involves the operation at a constant flow rate (28 or 60 L/min) although this does not represent the conditions typical of a tidal respiratory cycle. The combination of a breathing machine (Harvard type) and a vacuum pump to draw air through the impactor at a constant rate enables the measurement of aerosol deposition during the breathing tidal cycle [34].

4. Numerical studies

Some preliminary studies have been implemented in a Computational Fluid Dynamics (CFD) software code creating a computational model for a pMDI spray. In this way, a simple geometry test case (named “testbox”) was created and tested. These assumptions and simplifications are important in order to analyse spray flow without interference of the room walls.

The results given by the CFD analysis are normally a realistic approximation of a real-life system. The user has a wide choice regarding the level of detail of the results, since it is possible to simulate any fluid that cannot be reproduced experimentally due to economical or physical reasons. These advantages makes the CFD software to be a very powerful tool in the engineering research field [35,36].

The conservation equations for mass and momentum are solved in Fluent. Fluent is a finite volume based code, which uses the integral form of the conservation equations as its starting point. For this, a suitable grid for this simple geometry was generated. The solution domain is divided into a finite number of contiguous control volumes and the conservation equations are applied to each one. The governing equations are discretized on a curvilinear grid and Fluent uses a nonstaggered grid storage scheme to store the discrete values of dependent variables (velocities, pressure and scalars). A choice of interpolation schemes is then available.

By default, Fluent uses a Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm type to solve in a sequential way the discretized equations for mass and momentum. The algebraic equations are then solved iteratively using a Line-Gauss-Seidel solver and Multigrid acceleration techniques.

In the present application (spray in air flow), a dispersed second phase is included. Fluent models the dispersed second phase using a Lagrangian approach. The dispersed second phase is assumed sufficiently dilute and so, particle-particle interactions and the effects of the particles volume fraction on the gas phase are negligible. The calculations for the dispersed phase include the Lagrangian trajectory solution with stochastic tracking to account for the effects of gas.

The coupling between the phases and its impact on both the dispersed phase trajectories and the continuous phase flow can be included.

4.1. Gas flow and spray simulation

The air flow inside the “testbox” is solved as transient, incompressible, Newtonian and viscous turbulent. The governing equations of conservation of mass and momentum are solved with appropriate modelling procedures to describe the effects of turbulence fluctuations.

The spray modelling includes the calculation of the dispersed phase trajectories with the turbulence effects. The dispersed phase is included in the model by defining the initial conditions which will be used to initiate the second phase calculations. Initial and boundary conditions also need some attention.

The main equations are now presented, as well as, a brief description of the boundary conditions and numerical solution.

4.1.1. Mathematical modelling

The fluid conservation equations are simplifications to the Navier-Stokes equations, which are not possible to solve by conventional numerical means, each of them implies several considerations according to the specifics of the problem.

Mass conservation implies that the mass entering a control volume equals the mass flowing out, creating a balance between input and the output flows for a certain volume. This concept is mathematically expressed by Eq. 15, assuming constant the fluid properties [35–37].


where ρ stands for density, t for time, xi (i = 1, 2, 3) or (x, y, z) are the three-dimensional Cartesian coordinates and ui or (ux, ui, uz) are the Cartesian components of the velocity vector u.

The conservation equations used in Fluent for turbulent flows are obtained from those for laminar flows using a time averaging procedure usually known as Reynolds averaging.

For the momentum equations, the velocity at a point is considered as a sum of the mean and the fluctuating components:


Dropping the overbar on the mean velocity, u¯, the ensemble-averaged momentum equations for predicting unsteady state turbulent flows are described by Eq. 17 [35-37]:


Whereμrepresents the viscosity, prepresents static pressure, grepresents gravitational acceleration and Fthe source term for the momentum equation. The effect of turbulence is incorporated through the "Reynolds stresses" terms,ρui'uj'¯. Fluent relates the Reynolds stresses to mean flow quantities via a turbulence model. The most commonly used model for turbulence calculations is the k-ε with applications in several fields of engineering. It is classified as a Reynolds-Averaged Navier–Stokes (RANS) based turbulence model, considering the eddy viscosity as linear, it is a two equation model. This model accounts for the generation of turbulent kinetic energy (k) and for the turbulent dissipation of energy (ε). The model formulation described is typically called the Standard k-ε, and it is calculated by Eq. 18 and 19 [37].


In these equations, Gkrepresents the generation of turbulent kinetic energy due to the mean velocity gradients,σk and σε are the turbulent Prandtl numbers for the k and ε, respectively.C1ε andC2εare model constants. Tubulent viscosity, μt, is modelled according to the Eq. 20, by combining the values of k and ε.


This model uses, by default, the following values for the empirical constants:


Fluent allows for inlet, outlet, symmetry, cyclic and wall boundary conditions for treatment of the continuous phase. Two different inlet boundaries can be considered: inlet (velocity) boundary or pressure boundary. An inlet boundary is a boundary at which flow enters (or exits) at a known velocity, and turbulence parameters. At the pressure boundary, the fluid pressure is defined instead of velocity. At wall boundaries, the normal velocity component is zero.

The CFD software favours the use the Euler-Lagrange approach to track the particle trajectory through the fluid domain. The dispersed phase can exchange momentum and mass with the fluid phase. The particle or droplet trajectories are computed individually at specified intervals during the fluid phase calculation by integrating its equation of motion, which one is written in a Lagrangian reference frame. A force balance on the droplet immersed in a turbulent gas flow equates the particle inertia with the forces acting on the droplet, formulated (for the x-direction in Cartesian coordinates) in Eq. 21 [37].


where the first term in the right side accounts for the drag force, the second one for the gravity effects and Fx includes additional forces which can be important only in certain circumstances. In the first term,


where Re is the relative Reynolds number defined by:


When the gas flow field is turbulent, Fluent can predict the trajectories of the second phase using only the mean fluid phase velocity in the trajectory equations or optionally it can include the instantaneous value of the fluctuating gas flow velocity in order to predict the dispersion of the droplets due to turbulence.

In the present application, the inclusion of the turbulence effect on the droplets has been considered important. By computing the droplets trajectories for a sufficient number of representative droplets, the stochastic nature of turbulence can be modelled. The stochastic tracking (random walk) model includes the effect of instantaneous turbulent velocity fluctuations on the particle trajectories.

The dispersed second phase is introduced in the computational domain at a certain position (with initial conditions) and its trajectory is then computed along the continuous gas flow field.

Fluent provides different ways for input of the initial conditions for the dispersed phase. These initial conditions provide the starting values for all the dependent dispersed phase variables above mentioned.

Each set of these initial conditions represents an 'injection' of 'droplet type' identified in Fluent by a label 'injection number'. Each injection or parcel means thousands of droplets assumed to have the same size, diameter, temperature, starting at the same spatial position with the same velocity components. The mass flow rate indicates the amount of liquid that will follow the same calculated droplet trajectory.

The user can define any number of different sets of initial conditions for dispersed phase droplets provided that sufficient memory has been allocated.

The treatment of the dispersed phase near a boundary of the computational domain is a very difficult task. In Fluent, different dispersed phase boundary conditions can be applied when a droplet reaches a physical boundary. These boundary conditions are defined by the user at each cell-type and they can be: 'reflect', 'trap', 'escape' and 'saltation'. The 'reflect' condition rebounds the particle/droplet off the boundary with a change in its momentum defined by a coefficient of restitution; the 'trap' condition terminates the trajectory calculation and if it is a case of an evaporating second phase, their entire mass is converted into vapour phase; with the 'escape' condition the trajectory is simply terminated; with 'saltation', the particle/droplet is placed in the gas field at a small distance from the wall and again a coefficient of restitution is used to modified the particle momentum.

4.1.2. Numerical solution

As mentioned before, Fluent uses a control volume based technique to solve the conservation equations for mass, momentum, and turbulence quantities. The domain is divided into discrete control volumes where the governing equations are integrated to obtain the algebraic equations for the unknowns (velocities, pressure and scalars). The integration of the differential equations in each control volume yields a finite-difference equation that conserves each quantity on a control-volume basis.

For the temporal discretization of the first term in the conservation equations, the implicit Euler scheme has been used.

Because Fluent defines the discrete control volumes using a non-staggered storage scheme (all variables are stored at the control volume cell center), interpolation schemes are needed to determine the face values of the unknowns from the stored values at the cell center.

The discretized equations are solved sequentially and the SIMPLE algorithm has always been used in the present application. This type of algorithm is based on using a relationship between velocity and pressure corrections in order to recast the continuity equation in terms of a pressure correction calculation. In this way, the calculated velocity and pressure fields satisfy the linearized momentum and continuity equations at any point.

Fluent does not solve each equation at all points simultaneously and so an iterative solution procedure is used with iterations continuing until the convergence criteria specified has been achieved.

The algebraic equation for each variable is solved using a Line Gauss-Seidel procedure (LGS) and the user can specify the direction in which the lines are solved (the direction of the flow or alternate directions) and the number of times the lines are solved in order to update a given variable within each global iteration cycle. To speed up the convergence achieved by the LGS procedure, Fluent uses a Multigrid acceleration technique by default to solve the pressure and enthalpy equations.

The new calculated values of a given variable obtained in each iteration by the approximate solution of the finite difference equations are then updated with the previews values of the variable using a under relaxation technique. The user can choose the best relaxation factors for each variable in order to achieve a better convergence.

From the above, it can be noted that the trajectory equation for each parcel droplet is solved by step-wise integration over discrete time steps. Small steps have to be used to integrate the equations of motion for the droplet and Fluent allows the user to control the integration time step size when the equations are analytically solved.

The second phase trajectory calculations are based on the local gas flow field conditions, but Fluent does not include the direct effect of the droplets on the generation or dissipation of turbulence in the continuous phase. It keeps track of mass and momentum gained or lost by each droplet injection, as its trajectory is calculated. These quantities can be incorporated in the subsequent continuous phase calculations and, in this way, to study the effect of the dispersed phase on the continuous gas flow field. This two-way coupling is essential to calculate an accurate solution for the two-phase flow.

The momentum transfer from the continuous phase to the dispersed phase is computed in Fluent by examining the change in momentum of a droplet as it passes through each control volume in the computational domain. This momentum change is computed by the Eq. 24 [37].


The mass flow rate of each injection, which has no impact on the droplet trajectory calculation, is now used in the calculation of the second phase effect on the gas phase. It must be noted that when different trajectories are calculated for the same injection, in order to simulate the turbulence effects of the gas phase, the mass flow rate of that injection is divided equally by the number of stochastic tracks computed for that injection and so the exchange terms of momentum, mass and heat are calculated using the divided mass flow rate. This momentum exchange acts as a momentum sink in the gas momentum conservation equation (see Eq. 16).

The coupling between the two phases can be automatically simulated in Fluent by solving alternately the dispersed and continuous phase equations, until the solutions in both phases have stopped changing. Usually, the second phase calculations are made after a certain number of gas phase iterations.

Also, the exchange of momentum and mass terms calculated during the second phase simulation is not introduced directly in the gas phase conservation equations because Fluent allows for the under relaxation of the interphase exchange terms during the subsequent calculations. A small under relaxation factor seems to improve the convergence of the solution

4.2. Spray plume simulation

The present simulation aims to highlight the configuration parameters that better fits the real-life case of a pMDI spray, by using commercially available CFD software, such as Fluent version 14.0 from ANSYS®.

In this research work, the Ventolin® was used because it is one of the most common drugs used in developed countries to treat asthma in children and adults. It mainly consists of salbutamol, which is the most frequently prescribed Short-Acting β-Agonist (SABA) [9,38,39]. The used pMDI actuator was the one that is typically sold with the Ventolin®, both produced by the GlaxoSmithKline® company.

4.2.1. “Testbox” geometry and grid

The “testbox” consists of a simple parallelepiped form with the dimensions of 0.2 x 0.2 x 0.3 (m) representing a sample of a room environment, with a pMDI actuator and canister in the middle of it. The spray injection point, the exit of the actuator’s nozzle, is located in the origin point (0,0,0), see Figure 18.

The geometry was drawn using an external design program and then loaded into the ANSYS® meshing software. The generated mesh for the numeric calculations was composed by tetrahedral and wedge elements, with sizes ranging from 0.1 mm to 20.0 mm, resulting in a 3D computational grid with 3060339 elements and 1022403 nodes. Several refinements near to wall zones of high proximity and curvature were included. The quality report showed that mesh had a good quality according to the skewness parameter, with an average values located at 0.21 (see Table 3).


Figure 18.

A "testbox" representation constituted, the red plane (A) is the boundary condition ‘Velocity Inlet’ and the green plane (B) is the boundary ‘Outflow’. The other four outer walls of the domain have ‘Symmetry’ properties

SkewnessElement Quality
Standard Deviation0.12920.2641

Table 3.

Mesh quality criteria parameters: skewness and element quality

The boundary conditions were defined in two opposite faces, one as a ‘Velocity Inlet’ (see Fig. 2 A), forcing air to move inside the “testbox” at 0.01 m/s and the other as an ‘Outflow’ (see Fig. 2 B), enabling the freely motion of the air, as well as particles. For the remaining four faces, a ‘Symmetry’ boundary condition was assumed. The pMDI actuator and canister boundaries were considered ‘Wall’, trapping all the particles that collide with them.

4.2.2. Drug properties and spray characterisation

As mentioned previously, Ventolin® was used since it is a common SABA drug applied in the treatment of asthma with a pMDI. The Ventolin® most important characteristics for the simulation are listed in Table 4.

PropellantHFA 134a[11,40]
Salbutamol density (kg/m3)1230[11]
Actuation dose (µg)100[7,40,41]
Actuation time (s)0.1[15]

Table 4.

Ventolin® drug properties

The particle size data were measured using the laser diffraction technique, using an independent model for data processing. Data were fitted to the cumulative function of Rosin-Rammler distribution, using the least squares method [20]. The CFD solver accepts this distribution type by inserting its corresponding parameters, see Table 5 [37].

Diameter distributionRosin-Rammler
Minimum diameter (µm)1.22
Maximum diameter (µm)49.50
Mean diameter (µm)16.54
Spread parameter1.86

Table 5.

Particles diameter distribution parameters

4.2.3. Numerical parameters in Fluent

The spray parameters used to configure the solver were obtained from various references, though some caution is required, such as, the angle of the spray was considered to be 10 degrees, estimated by high speed image analysis. Assumptions were made, such as, the spray particles to be solid instead of the well-known liquid droplets. The reason for this consideration is simple: when the drug exits the metering chamber of the pMDI, it undergoes a flash evaporation. Because this is an instantaneous process, it is assumed that no heat transfer between the gas and liquid phases [5,13,14]. Using the drug amount of a puff (100 µg) and divide it by the duration of a puff (0.1 s), it is possible to calculate the spray flow rate, see Table 1. The main spray parameters used to configure the solver are summarized in Table 6, such as the spray angle, maximum velocity of the spray plume, shape of the plume, as well as, the dimension of the actuator’s nozzle.

Spray typesolid-cone[12,37]
Angle (º)10
Velocity (m/s)100[12,42]
Radius (m)0.00025[12,15,42]
Flow rate (kg/s)1e-6

Table 6.

Solver spray configuration parameters

The solution of the differential equations for mass and momentum was done in a sequential manner, using the SIMPLE algorithm [3537,43]. The standard discretisation scheme was used for the pressure and the second order upwind scheme for the momentum, turbulent kinetic energy and turbulent dissipation rate equations. Convergence was reached in the simulation by using a criterion value of 1.0e-5 for the continuity (pressure), x, y and z velocity, and for k and ε turbulence parameters.

The time step used in the simulation was of 0.01 s, during 11 time steps with a maximum limit of 3000 iterations for each time step. Also the gravitational acceleration was assumed 9.81 m/s2 in the y axis direction, although for these particles diameter range, this effect might be negligible.

Some parameters used to configure the Discrete Phase Model (DPM) model in the solver are listed in Table 7. Although different configuration parameters have been tested, the ones presented here seemed to be the most appropriate to be used in the simulation. Due to the fact that the injection occurs in a limited period in time (0.1 s), the unsteady particle tracking was used. The interaction with the continuous phase was also included to approximate the simulation to the real event, because the velocity difference between the two phases is high.

Interaction with Continuous PhaseOn
Unsteady Particle TrackingOn
Inject Particles atParticle Time Step
Particle Time Step Size (s)0.001
Drag LawSpherical
Two-way coupling turbulenceOn

Table 7.

Parameters for the DPM model

The drag law used (spherical) is the simplest and most used law, once it is the one that better fits the presented model amongst the four different options available in the solver. This drag law considers the particle as a sphere, which is an acceptable simplification for the drug particles that exit the pMDI nozzle [14,37].

The total number of particle streams injected during the simulation was approximately 164000.

4.2.4. Results

The CFD results obtained are represented in the form of images, showing the contours of velocity along a longitudinal mid-section plane and as streams of tracked particles along the fluid domain at the last injection moment.

Figure 19 shows the velocity magnitude of air in the domain at the time of 0.11 seconds, which coincides with the end of the injection of the spray. It is possible to observe the high velocity that air shows at the nozzle zone, accelerated by the discharge of particles, reaching a velocity of 10 m/s and lowering its velocity until reaches an equilibrium, the “testbox” air velocity.

In Figures 20 and 21, the representation of the particles as spheres scaled by its diameter provides a direct reading of the particles diameter and positioning at the end of the injection. The bigger particles are at the end of the plume and the smaller ones are next to the nozzle. Although the particle size distribution configured provides the expected result, there are several factors that might influence numerical shape of the spray plume, such as the solver injects always the same quantity of particle streams at each instant of the injection, which is known not to be the spray real behaviour.


Figure 19.

Velocity magnitude contours of the air in the domain, along the longitudinal midsection plane (t=0.11 s)


Figure 20.

Representation of the particle streams at the end of the injection (t=0.11 s), image shows the particles coloured by residence time inside the domain. The particle streams are draw as spheres with proportional size scaled 50 times more than the real diameter


Figure 21.

Representation of the particle streams at the end of the injection (t=0.11 s), image shows the particles coloured by its velocity magnitude. The particle streams are draw as spheres with proportional size scaled 50 times more than the real diameter

In Figure 21, it can be observed that particles velocity drastically reduces after the injection into the still air, which was expected, but might not be fully accurate due to the simplification of the drag law used. In Figure 20 the residence time of each particle stream in the domain is coloured, showing that all diameters take the maximum time of 0.1 seconds, the duration of the injection, which means that particle streams are inside the domain since the beginning of the injection.

5. Conclusions

The authors believe this chapter will provide a meaningful insight into the pMDI spray thematic, with a significant description that covers several pMDI spray topics, such as, the components, mechanisms, distributions, different measurement techniques and numerical simulations, regarding the gas flow and spray simulation as well as the spray plume simulation. This review clearly highlights that a wide variety of techniques can be applied to the study of the drug delivery systems through the respiratory system. The experimental and numerical techniques complement each other to enhance the understanding of drug transportation into the respiratory tract.

The computational techniques can be successfully used to study complex flows in a cost effective manner. The availability of high performance computing tools enables the introduction of complex physical models into even more complex geometries that are representative of in vivo systems.

Nonetheless, experimentation is still a cornerstone approach to validate models, to turn them into more reliable and accurate simulation tools that should closely characterize the physical events they represent. The most recent techniques can be used to provide very detailed data of flow patterns, such as particle size and motion, amongst other parameters. Different techniques have been referred and summarized in the study herein reported, regarding laser diffraction techniques, phase Doppler anemometry and the use of photography and high-speed cameras to determine particle sizing and moving, as well as impacting techniques, to measure spray particle sizes.

The CFD spray study proved to be useful in the study of add-on devices to the pMDI, such as holding chambers (or spacers) design. The velocity contours determined and the streams of tracked particles, obtained at a time corresponding to the end of the spray injection, suggests interesting aspects about the behaviour and performance of the spray and particles positioning and distribution, although other simulations need to be carried out in the future with the introduction of the spacer device in the “textbox” geometry. These aspects are currently under investigation by the authors and they represent a further contribution given to understand drug delivery through the oral system and to help the design of even more efficient delivery instruments.


1 - Global Initiative for Asthma, 2010 Global strategy for asthma management and prevention, GINA.
2 - J. C. Virchow, G. K. Crompton, Negro. R. Dal, S. Pedersen, A. Magnan, J. Seidenberg, P. J. Barnes, 2008 “Importance of inhaler devices in the management of airway disease,” Respiratory medicine, 102(1), 1019.
3 - M. B. Dolovich, 2004 “In my opinion- Interview with the Expert,” Pediatric asthma Allergy & Immunology, 17(4), 292300.
4 - D. P. Tashkin, 1998 “New devices for asthma,” Journal of Allergy and Clinical Immunology, 101(2), S409S416.
5 - M. B. Dolovich, J. B. Fink, 2001 “Aerosols and devices,” Respiratory care clinics of North America, 7(2), 131173, v.
6 - S. P. Newman, 2006 “Aerosols,” Encyclopedia of Respiratory Medicine, G.J. Laurent, and S.D. Shapiro, eds., Elsevier, 5864.
7 - C. Terzano, 2001 “Pressurized metered dose inhalers and add-on devices.,” Pulmonary pharmacology & therapeutics, 14(5), 351366.
8 - S. P. Newman, 2004 “Spacer devices for metered dose inhalers.,” Clinical pharmacokinetics, 43(6), 349360.
9 - G. Crompton, 2006 “A brief history of inhaled asthma therapy over the last fifty years.,” Primary Care Respiratory Journal, 15(6), 326331.
10 - M. Sanders, 2007 “Inhalation therapy: an historical review,” Primary care respiratory journal, 16(2), 7181.
11 - H. Smyth, 2003 “The influence of formulation variables and the performance of alternative propellant-driven metered dose inhalers,” Advanced Drug Delivery Reviews, 55, 807828.
12 - S. P. Newman, 2005 “Principles of metered-dose inhaler design,” Respiratory Care, 50(9), 11771190.
13 - C. A. Dunbar, 1997 “Atomization mechanisms of the pressurized metered dose inhaler,” Particulate Science and Technology, 15(3-4), 253271.
14 - W. H. Finlay, 2001 The mechanics of inhaled pharmaceutical aerosols: an introduction, Academic Press.
15 - H. Smyth, A. J. Hickey, G. Brace, T. Barbour, J. Gallion, J. Grove, 2006 “Spray pattern analysis for metered dose inhalers I: Orifice size, particle size, and droplet motion correlations,” Drug development and industrial pharmacy, 32(9), 10331041.
16 - K. J. Mc Donald, G. P. Martin, 2000 “Transition to CFC-free metered dose inhalers- into the new millennium,” International Journal of Pharmaceutics, 201(1), 89107.
17 - D. Tiwari, D. Goldman, W. a. Malick, P. L. Madan, 1998 “Formulation and evaluation of albuterol metered dose inhalers containing tetrafluoroethane (P134a), a non-CFC propellant.,” Pharmaceutical development and technology, 3(2), 163174.
18 - L. P. Pereira, Y. Clement, D. Simeon, 2001 “Educational intervention for correct pressurised metered dose inhaler technique in Trinidadian patients with asthma,” Patient Education and Counseling, 42, 9197.
19 - S. P. Newman, S. W. Clarke, 1983 “Therapeutic aerosols 1- physical and practical considerations,” Thorax, 38(12), 881886.
20 - C. A. Dunbar, A. J. Hickey, 2000 “Evaluation of probability density functions to aporoximate particle size distributions of representative pharmaceutical aerosols,” Journal of Aerosol Science, 31(7), 813831.
21 - R. A. Mugele, H. D. Evans, 1951 “Droplet Size Distribution in Sprays,” Industrial & Engineering Chemistry, 43(6), 13171324.
22 - Z. Zhang, S. Ziada, 2000 “PDA measurements of droplet size and mass flux in the three-dimensional atomisation region of water jet in air cross-flow,” Experiments in Fluids, 28(1), 2935.
23 - L. G. Dodge, J. A. Schwalb, 1989 “Fuel Spray Evolution: Comparison of Experiment and CFD Simulation of Nonevaporating Spray,” Journal of Engineering for Gas Turbines and Power, 111(1), 15
24 - Black. D. Lee, M. Q. Mc Quay, M. P. Bonin, 1996 “Laser-based techniques for particle-size measurement: A review of sizing methods and their industrial applications,” Progress in Energy and Combustion Science, 22(3), 267306.
25 - M. Hong, A. Cartellier, E. J. Hopfinger, 2004 “Characterization of phase detection optical probes for the measurement of the dispersed phase parameters in sprays,” International Journal of Multiphase Flow, 30(6), 615648.
26 - J. P. Mitchell, M. W. Nagel, 2004 “Particle Size Analysis of Aerosols from Medicinal Inhalers,” KONA, 22, 3265.
27 - L. P. Bayvel, 1980 “Application of the laser beam scattering technique method for multi-phase flow investigations,” Proceedings of the conference on multiphase transport, fundamentals and reactor safety applications, Y.N. Veziroglu, ed., New York.
28 - J. Swithenbank, J. M. Beer, D. S. Taylor, D. Abbot, G. C. Mc Creath, 1976 “A laser diagnostic technique for the measurement of droplet and particle size distribution,” American Institute of Aeronautics and Astronautics, 14th Aerospace Sciences Meeting, Washington, D.C.
29 - A. H. Lefebvre, 1989 Atomization and Sprays, Hemisphere Publishing Corp., New York and Washington, D.C.
30 - Dantec, 1995 PDA user’s guide, Dantec information.
31 - N. E. Tayali, C. J. Bates, 1989 Scattered light intensities and phase difference derived using geometrical optics.
32 - X. Liu, W. H. Doub, C. Guo, 2012 “Evaluation of metered dose inhaler spray velocities using Phase Doppler Anemometry (PDA),” International journal of pharmaceutics, 423(2), 235239.
33 - J. P. Mitchell, M. W. Nagel, 2003 “Cascade impactors for the size characterization of aerosols from medical inhalers: their uses and limitations.,” Journal of Aerosol Medicine, 16(4), 341377.
34 - S. A. Foss, J. W. Keppel, 1999 “In Vitro Testing of MDI Spacers : A Technique for Measuring Respirable Dose Output with Actuation In-Phase or Out-of-Phase with Inhalation,” Respiratory Care, 44(12), 14741485.
35 - J. H. Ferziger, M. Peric, 2003 Computational methods for fluid dynamics, Springer.
36 - H. K. Versteeg, W. Malalasekera, 1995 An introduction to computational fluid dynamics: the finite volume method, Longman, Harlow, England.
37 - ANSYS, 2009 ANSYS FLUENT Theory Guide, ANSYS Inc, Canonsburg, PA, USA.
38 - G. Jepson, T. Butler, D. Gregory, K. Jones, 2000 “Prescribing patterns for asthma by general practitioners in six European countries,” Respiratory medicine, 94(6), 578583.
39 - M. G. Zuidgeest, H. A. Smit, M. Bracke, A. H. Wijga, B. Brunekreef, M. O. Hoekstra, J. Gerritsen, M. Kerkhof, J. C. Jongste, H. G. Leufkens, 2008 “Persistence of asthma medication use in preschool children,” Respiratory Medicine, 102, 14461451.
40 - J. C. Dubus, C. Guillot, M. Badier, 2003 “Electrostatic charge on spacer devices and salbutamol response in young children,” International Journal of Pharmaceutics, 261(1-2), 159164.
41 - S. Verbanck, C. Vervaet, D. Schuermans, W. Vincken, 2004 “Aerosol Profile Extracted from Spacers as a Determinant of Actual Dose,” Pharmaceutical Research, 21(12), 22132218.
42 - A. R. Clark, 1996 “MDIs: physics of aerosol formation.,” Journal of aerosol medicine: the official journal of the International Society for Aerosols in Medicine, 9 Suppl 1(s1), S19S26.
43 - S. Abreu, L. F. Silva, H. Antunes, S. F. C. F. Teixeira, 2008 “Multiphase flow inside the Volumatic spacer: a CFD approach,” Proceedings of the 10th International Chemical and Biological Engineering Conference- CHEMPOR 2008, E.C. Ferreira, and M. Mota, eds., Braga, 6