Open access peer-reviewed chapter

A General-Purpose Multiphase/Multispecies Model to Predict the Spread, Percutaneous Hazard, and Contact Dynamics for Nonporous and Porous Substrates and Membranes

Written By

Navaz Homayun, Zand Ali, Gat Amir and Atkinson Theresa

Submitted: October 22nd, 2014 Reviewed: May 12th, 2015 Published: December 16th, 2015

DOI: 10.5772/60807

Chapter metrics overview

1,822 Chapter Downloads

View Full Metrics


A computational model to solve the coupled transport equations with chemical reaction and phase change for a liquid sessile droplet or the contact and spread of a sessile droplet between two approaching porous or non-porous surfaces, is developed. The model is general therefore it can be applied to toxic chemicals (contact hazard), drug delivery through porous organs and membranes, combustion processes within porous material, and liquid movements in the ground. The equation of motion and the spread of the incompressible liquid available on the primary surface for transfer into the contacting surface while reacting with other chemicals (or water) and/or the solid substrate are solved in a finite difference domain with adaptive meshing. The comparison with experimental data demonstrated the model is robust and accurate. The impact of the initial velocity on the spread topology and mass transfer into the pores is also addressed.


  • Liquid spread in porous materials
  • gaseous spread in porous material
  • evaporation in porous material
  • phase change in porous materials
  • multi-phase process in porous materials
  • multi-species process in porous materials

1. Introduction

Many hazardous chemicals are dispersed in their liquid phase as droplets and reside on the surface of the liquid shortly after dispersion. In this state they are referred to as sessile droplets. Many coupled processes commence after this stage, resulting in different forms of threats, depending on the surface and physical properties of the liquid droplet. For nonporous surfaces, if the vapor pressure of the liquid is high, the threat of exposure is through the respiratory system, and if the vapor pressure is low, a percutaneous or direct contact hazard can be expected. Furthermore, should liquid droplets enter a chemical reaction with the surface, the products or the remaining reactants can reduce or augment the threat. In general, the following processes can be expected on a nonporous surface:

  1. Surface evaporation

  2. Chemical reaction or adsorption at the contact surface only

  3. Formation of liquid bridge when a secondary surface come into contact with the sessile droplet on the nonporous or primary surface

  4. Spread of the liquid bridge as the gap between the two surfaces closes.

For porous surfaces the following processes can be simultaneously present:

  1. Surface evaporation

  2. Capillary transport

  3. Evaporation through the pores referred to as secondary evaporation

  4. Chemical reaction and adsorption inside the pores

  5. Possible formation of liquid bridge when a secondary surface comes into contact with any possible sessile droplet left due to capillary transport on the primary substrate in the postcontact phase

  6. Spread of the liquid bridge if the conditions for its formation are met

  7. Transfer of chemicals from the primary surface to the contacting surface occurs when the droplet on the primary substrate is observable or nonobservable. Therefore, contact hazard can exist even though the liquid droplet cannot be seen.

Although nonporous and porous surfaces share some processes, there are some differences in their outcome. Furthermore, in case of temperature variation, phase change can occur, thereby necessitating the inclusion of the energy equation in exposure risk models.

It is evident that the processes such as capillary transport (for porous surfaces), surface evaporation, in-pore or secondary evaporation, adsorption, and chemical reaction are initiated simultaneously as soon as a liquid droplet contacts a surface. In the case of a contact, a liquid droplet can be transported to the secondary (or contacting surface). This can occur when the droplet is sandwiched between the two surfaces, with the droplet forming a liquid bridge, or long after the droplet has disappeared from the surface. Therefore, an additional coupled process, which is mainly physical (liquid bridge spread), can be added to the above processes. The end result of all these processes determines the available chemical for atmospheric transport in the form of vapor or percutaneous hazard on surfaces available for transfer to a secondary surface upon contact.

When a wetting fluid is introduced into a porous medium, in the absence of external forces, the distribution of the liquid in pores is influenced by pressure differential, capillary, gravitational, and viscous forces. The main driver for the spread of this wetting fluid is capillary pressure. As the porous medium imbibes the liquid, the liquid flow rate decreases due to the domination of viscous forces, causing a reduction in the flow capillary number. The spread dynamics are influenced by the porous medium thickness, where two- and three-dimensional flows occur. For thin porous media, the liquid flows in a radial direction and can be modeled by the method of “common-lines” [1, 2]. In these solutions, two parameters need to be experimentally determined. In thick porous media, the flow becomes three-dimensional, although a cylindrical symmetry has been assumed in the past [3, 4] The experimental results of Denesuk et al. [5] verified this three-dimensional nature of the flow. Holman et al. [6] studied the surface spread in a porous medium as a function of permeability and have demonstrated that the surface spread is the dominant mechanism as compared to droplet penetration dynamics. In the above studies, a well-defined interface between the wetted and nonwetted volumes existed. It is known that spread dynamics becomes more complex when the saturation is less than unity, as demonstrated via MRI experiments by Mantle et al. [7]

When a droplet spreads into a porous substrate, saturation gets distributed between 0 and 1, necessitating a multiphase flow approach. Therefore, the governing equations need to account for the presence of all the phases [8, 9]. In these equations, two additional transport parameters are present: a relative permeability and a capillary pressure function that includes the interaction parameters between phases as shown by Chen et al. [10]. The relative permeability can be modeled using a power function as indicated by Brooks and Corey, [11] or van Genuchten equations, [12] in which each equation has a large number of parameters expressing the relative permeability. The capillary pressure is often modeled using the Leverett [13]and Udell [14] J-function.

In solving the droplet spread on or into a porous substrate, the shape of the droplet on the porous substrate surface (droplet free volume) and the shape of the wetted porous substrate imprint must be determined. Starov et al. [2, 15] have solved the governing equations to obtain the shape of a droplet free volume that spreads into a thin porous medium where he found a large variation of the droplet base radius. Using full numerical calculations for a three-dimensional porous medium, Alleborn et al. [4] have come to a similar conclusion. In respect to the imprint shape, the spread of the droplet is usually modeled by clearly separating the fully saturated and non-saturated regions. However, the spread dynamic is more complex due to the existence of partially saturated regions and varying transport parameters within the pores as indicated by Lenormand et al. [16]

The change of the flow mechanism, and the formation of partially saturated regions, requires solving the droplet spread into porous medium as a multiphase problem for which the multiphase flow parameters have to be determined. The multiphase parameters depend on the type of process and the predominance of gravity, viscous, or capillary forces. [16] Valavanides and Payatakes [17] summarized these influences on the multiphase parameters, where, in addition to phase content or phase saturation and emerging forces, solid–fluid contact angles, ratios of phase viscosities and flow rates, and the flow history (e.g., drainage or imbibition) altered the behavior of multiphase parameters. Discrete pore network models provide an alternative approach to elucidate the transfer phenomena and to evaluate transport parameters. In network models, an actual porous medium is represented as a network of pores that are connected by throats. [18] The phase progression is set using a potential threshold and the flow pattern is determined thereafter, as outlined by Prat. [19] From the phase distribution, the multiphase parameters are calculated (see, e.g., Constantinides and Payatakes [20]).

Further studies using a discrete pore network split the flow into primary and secondary regimes. The primary flow occurs when the volume of the sessile drop is greater than zero and the secondary flow happens after the sessile droplet volume becomes zero [21, 22, 23, 24]. In the secondary flow regime, the flow front at the surface becomes an irregular flow front, forming a heterogeneous porous medium at the surface. [25]

In spite of all previous progress in the field of porous media flow, there was still a need for a production code (or model) that included all the physiochemical processes in order to calculate the amount of a hazardous liquid droplet that poses a threat. The first problem encountered, in a realistic scenario, is the evaporation of a sessile droplet that is deposited onto an impermeable surface and exposed to outside wind conditions. This is called “convective evaporation” and its driving force is the outdoor wind velocity. As mentioned previously, a convective evaporation (laminar or turbulent) model needs to be developed in the absence of all other processes, i.e., for non-permeable non-reacting surfaces. An evaporation model for a sessile droplet on a non-porous surface with stagnant air above it is discussed by Popov [26]. Hu et al. [27] used finite element method to mesh the domain and then solved the mass conservation equations for the vapor concentration to find the evaporation rate. The above models assume no motion of air above the sessile droplet eliminating the convective effects. Baines et al. [28] developed a very simple convective evaporation model that did not compare with experimental results very well and in some cases overpredicted the evaporation rate by 90%.

A comprehensive literature review of the existing models is presented by Winter et al. [29] All of the previous studies focus on capturing the physics of evaporation and do not address the role of turbulence on evaporation. The effect of turbulence on evaporation of droplets that are moving in the core flow has been addressed by Navaz et al. [30] and Ward et al. [31] for spray combustion problems where high density and temperature gradients will reduce the evaporation time to the scale of turbulence fluctuations. However, for smaller or more moderate temperature gradients, the turbulence time scale will be smaller than the evaporation time scale. Therefore, a new approach should be developed to examine the effects of turbulence on evaporation of sessile drops. There are numerous works that examine atmospheric turbulence. The effect of the free stream Reynolds number on shear stress and friction velocity at a wall under outdoor conditions is discussed by Metzger and Klewicki. [32] Their study was performed on the Great Salt Lake Desert and spans over Re=2000 – 5e6 or Re=2000 - 5× 106. There is also a dataset for friction velocity measured over a period of several months compiled from the outdoor wind data of Klinger et al. [33]. It has been mentioned by Weber [34] that the main scaling parameter in similarity theory of atmospheric boundary layer is friction velocity. He concludes that there could be a significant difference in friction velocity depending on the type of air flow. Several studies have focused on the varying turbulence intensity in wind-driven flows under outdoor conditions [35, 36]. Their approach relies on the validity of some wall function that relates the friction velocity to the mean free stream velocity and height. Upon knowing the friction velocity, the shear stress and turbulence intensity are correlated. For the evaporation problem, the previous studies have correlated the free stream velocity (or Reynolds number) to the friction velocity. However, they have not addressed conditions where turbulence intensity can change while the mean free stream velocity stays constant, or the turbulence intensity can be altered without changing the free stream velocity (or mean velocity in a channel flow).

Once a droplet spreads into a porous substrate, secondary evaporation (from the porous substrate) occurs. This mass transport and vapor diffusion is hindered by granules of the porous medium, and the vapor may face liquid parcels of the agent that act as additional obstructions. These obstructions greatly alter the transport properties of the vapor inside the pores and the transport rate toward the surface. Basically, the transport properties for the vapor phase become a function of porosity and saturation. These properties are yet to be determined.

As mentioned earlier, capillary transport of sessile droplet of a chemical agent on a porous substrate can be further complicated by the presence of chemical reaction. The chemical reaction can occur between the chemical agent and the constituents of a porous substrate or the chemical agent and pre-existing chemicals inside the pores, such as water or other chemicals. It is known that some chemical agents undergo hydrolysis reaction when they encounter moisture in the soil or sand, as it will be shown later. A simple set of one-dimensional mass and energy equations for a medium with no motion (conductive only) with chemical reaction was solved by Sadig [37, 38] who also demonstrated how the local properties are affected by chemical reaction, as the composition of constituents vary in time. Xu et al., [39] Lichtner, [40] and Steefel [41] have solved a very simple reactive flow with diffusion in geological settings dealing with ion transport. To have a general-purpose model, it is required that the conservation equations for all the phases (solid, liquid, and gaseous) and existing constituents, in forms of reactants and products with a chemical reaction, be solved simultaneously in time. The variation of species concentration in any of the existing phases will alter the local properties of liquid, gas, or solid mixture. Furthermore, the production or destruction of any solid phase will also change the local porosity of the medium.

A low-volatile chemical in its liquid form poses a threat upon contact with skin or other materials (vehicles, etc.). Mass transfer to the secondary surface occurs during this process. However, the contact may occur before a sessile droplet completely sinks in the primary substrate, if porous, or after a sessile droplet is not observable on the surface. In the first case, a liquid bridge is formed between the two approaching surfaces and may spread depending on the surface properties. Furthermore, if the primary surface or substrate is non-porous, a liquid bridge is definitely formed and spreads as the distance between the approaching surfaces decreases. The forces applied by liquid bridges connected to static supporting surfaces were studied extensively by Butt et al. [42] However, in many applications, it is common that at least one of the bodies is moving, and thus, the motion of the solid body may be influenced by the forces associated with the liquid bridge [43, 44, 45]. Pitois et al. [46] studied the forces applied by a liquid bridge connecting two spheres moving at a constant speed relative to each other. Similarly, Meurisse and Querry [47] studied the effects of liquid bridges connecting two parallel flat plates, moving perpendicularly to the plane of the plates, at a constant speed, or at a constant force applied on the liquid. Both Pitois et al. [46] and Meurisse and Querry [47] observed a rapid change from an attractive force due to capillary effects to a repulsive force due to viscous effects for the case of surfaces approaching each other at a constant speed. De Souza et al. [48] studied the effect of contact angle hysteresis on the capillary forces in the absence of significant viscous forces and obtained good agreement between the experimental data and the existing models. The capillary force measurement of liquid bridges was examined and shown to be a reliable method for estimating advancing and receding wetting angles.

The absorption of chemical into the skin is of primary interest due to threat from contact with hazardous materials. Experimental studies of Ngo, O’Malley, and Maibach [49] have attempted to quantify mass absorption into the skin. This is a complex problem because skin chemistry and perspiration rate may vary between individuals affecting the absorption rate. Although some simplified models based on linear regression exist to estimate the absorption rate, a more comprehensive and scientifically sound approach provides the solution of coupled equations for multi-species systems, as discussed hereinafter.

However, there was still a need for a model that can predict the available threat from a liquid chemical toxicant release. This issue has been addressed by Kilpatrick et al. [50], Savage [51, 52], and Munro et al. [53], and prompted the Department of Defense (DoD), Edgewood Chemical and Biological Center (ECBC), and Defense Threat Reduction Agency (DTRA) to initiate a comprehensive research project called Agent Fate Program that was geared toward chemical agents of warfare. This project led to the development of the first-generation production code COMCAD (COmputational Modeling of Chemical Agent Dispersion) by Navaz et al. [54, 55], and the next generation for contact hazard MOCHA (Modeling Of Contact HAzard) codes. [56] Since the model is general, it can not only be applied to chemical agents of the warfare, but also any other liquid that can be disseminated into the environmental substrates in a sessile droplet form (e.g., pesticides, toxic industrial chemicals (TICs), and toxic industrial materials (TIMs)). This chapter summarizes the result of this research and its application to a variety of problems.

1.1. Mathematical model

The species mass and momentum, and energy equations for a multicomponent system in a porous medium with all phases being active is rather complex and are given in Navaz et al. [57], Navaz et al., [55] Vafai, [9] Navaz et al., [58] and Kuo. [59] The continuity equation for a multicomponent system for all phases includes any loss of mass due to surface and secondary evaporation, and cross-over mass due to chemical reaction, phase change, and adsorption. The momentum equations in porous media are considered for liquid and gaseous phases. Local thermal equilibrium is considered for the energy equation. That is to say that the local temperatures of all phases are equal (Ts=T=Tg). Note that the solid phase is stagnant, but the porosity will be changing in time if any solid constituent is formed due to phase change or chemical reaction. The continuity equation for the porosity is eliminated, but the porosity is updated in each time step (lagged by one time step). Although for an open system under the usual atmospheric conditions the gaseous phase pressure normally stays at the environmental conditions, the pressure term is included in the equations to extend the applicability of the model for the future applications, in which the pressure may change (closed systems or existence of explosions). All the conservation equations are solved explicitly on a finite difference mesh as will be described later. The sub-models describing the source terms (e.g., evaporation), chemical kinetics, and adsorption models are described in separate sections. We start with explaining the evaporation rate term.

Mass Conservation for Multiple Solid System (Substrate – Solid "k")


Mass Conservation for Multiple Solid System (Liquid "i")

(φρisi)t+(φρisiVi)=(ρ˙iSecondary  Evaporationρ˙iSurface  Evaporationω˙iReactionω˙iAdsorption)E2
ρmixture=i=1N(Liquids)ρiCi       Ci=Mass  Fraction=ρisii=1Nρisi,  s=i=1N(Liquid)siE3
Mass Conservation for Multiple Gaseous System (Gas"j")(φρgjsg)t+(φρgjsgVgj)=(ρgDjmixCj)+ρ˙,i=jSecondary  Evaporationω˙gjReaction-ω˙gjAdsorptionE5
kgj=Relative  Permeability to jth  constituent of the gas phase
μg(gas   mixture  viscosity) =j=1M(Gases)μgjCgj,E7
ρgj=Density of jth  constituent of the gas phaseVgj=Velocity  component of jth  gas  phase  constituent


sg=1iN(Liquids)si and ρg=j=1M(Gases)ρgjE8
Djmix=Effective  Diffusion Coefficient of gas species "j" into the mixture
Cj=Mass  Fraction=ρgjρgE9
Pi=PPci (Capillary Pressure)+ρigh*E10
h*=Local height of the droplet as a function of time (hydrostatic pressure)ω˙=Inter-Phase  Chemical  Species production  or  destructionK=Saturation  permeabilityk =k(s)   Relative permeability is a function of the local saturation=s2φ=φ(ρsk,x,y,z)  Local porosity is a function of local density of the solid phasebb1

Energy Equation



φs=(1φ), φ=sφ, φg=(1s)φ, V2=u2+υ2+w2E13
q˙=Heat   flux=keffT, Δh=Latent  heat  of  phase  change     ρ˙=rate  of  phase  Change H= Total enthalpy = cpT+V22   cp=Specific  heat τ=Shear  stress tensor given by:
τxx=μ[2ux23(ux+υy+wz)],      τxy=τyx=μ[uy+υx],    τxz=τzx=μ[wx+uz] τyy=μ[2υy23(ux+υy+wz)],     τyz=τzy=μ[υz+wy], τzz=μ[2wz23(ux+υy+wz)]E14

1.2. Evaporation models

When surface evaporation of a sessile droplet on a non-porous substrate occurs, one or both of the following mechanisms determines the topology of the sessile droplet. Either surface evaporation takes place with the droplet’s base area (radius) remaining constant while the height of the droplet decreases, or the surface evaporation causes the shrinking of a droplet while the ratio h/r remains constant. It has been observed that for the evaporation of nerve agents (HD, VX, and GD) the process starts with the first mechanism and when the contact angle reaches about 12°, the second mechanism is initiated.

Irrespective of the evaporation mechanism (r=const. or λr= h/r =const.), the volume and surface area of a sessile droplet in a form of a spherical cap is given by:


The dynamics and topology of this model is shown in Figure 1. where the evaporation is initiated with the first mechanism and then is switched to the h/r=const. case.

Having the mass rate of evaporation (m˙) equal to the negative change of the mass (m) left on the surface at time (t) and it can be written as:

m˙=ddtm=ddt(ρ V)=ρπ2(r2+h2)dhdtE16

where ( ρ ) is the liquid density. From Equation (16) the height of the droplet is calculated as:


Figure 1.

Schematic for the Topology and Dynamics of HD Evaporation on a Non-Permeable Surface.

The instantaneous droplet height can be obtained by integrating Equation (17), provided that the forcing function m˙ is known. There are numerous expressions in the literature for the evaporation rate m˙ given by researchers in spray combustion field [58, 59, 60, 61] for spherical droplets following a gas trajectory in a rocket engine. The current model follows the same process with a few exceptions. Toxic chemicals are generally tested in small wind tunnels mainly for better control of contamination. The flow in these wind tunnels is basically a channel flow and a boundary layer plate cannot be installed (due to size constraints) to create a scalable model to the open air scenarios. To overcome this issue, the Reynolds number in Equation (18) was based on the friction velocity instead of the free stream velocity. Embedded in the friction velocity is the impact of turbulence on evaporation. This is discussed by Navaz et al. [57] and is the focus of Section 1.3. Further modification was necessary to define a new transfer number, B. The transfer number in combustion is defined as: cp(TgTboiling point)hfgwhere hfg is the latent heat of evaporation. Although this definition will work in a combustion chamber, it does not apply to the evaporation of toxic chemicals under environmental conditions that mainly involve the vapor pressure of the sessile droplet and convective mass transfer. It is assumed that the evaporation process takes place isothermally and any minimal change of temperature at the interface and its conduction through the sessile droplet can be ignored. Two different length scales are defined depending on the droplet topology. [57] The following surface evaporation model, Equation (18), provides the source term (ρ˙surface evaporation) in Equations (2) and (5) for sessile droplets, which is derived, calibrated, and used for this study.

m˙surface  evaporation=2πLcμPr(Cf+C1RecmPrcn) n[1+B]Lc=Length scale {Mechanism 1: = 2πRsλrMechanism 2: =2πRs[1-sin(π2θa)]  
Re=Reynolds number based on the radius of curvature and friction velocityρg= Gas mixture densityB=Transfer Number = (y1y)ε,     y=PυjPE18
Pυj=Vapor pressure of species "j"where: m is the mass of liquid species "j" at each nodeCf,C1,cm,cn,ε=Model constants to be determined experimentally

The evaporation model for flow through the pores is also based on the methodology discussed by Navaz et al. [55, 58]. A numerical investigation of this process on hot surfaces has been published by Nikolopoulos et al. [62] However, in the context of finite difference or finite volume framework, a new model needs to be developed. It is assumed that the liquid phase at each node is spherical and the instantaneous length scale is found by knowing the mass and density of each liquid species at this node. The evaporation rates in the species continuity equation (Eq. 2) are obtained as indicated by Equation (19).

m˙Secondary evaporation=4πρgDjMixLc(C1+C2RemScn) n[1+(PυjP)ε1(PυjP)ε]ρ˙Secondary  Evaporation=m˙/cell
Sc=Local Schmidt  Number=νjDjMixρg= Gas mixture densityDjMix= Effective Diffusion Coefficient of species "j" into the mixturePυj=Vapor Pressure of species "j"P= Gas mixture pressureE19
Lc=Characteristic Length = (mjρj)13Re= Reynolds Number =Lcu2+υ2+w2νj   where: mj is the mass of liquid species "j" at each nodeC1,C2,C3,m,n,ε=Model constants to be determined experimentally

1.3. Scalability of surface evaporation model

As mentioned earlier, the Reynolds number in the evaporation model is based on friction velocity. This is due to the fact that a similar solution can be found for all cases involving the evaporation of a sessile droplet. Furthermore, the turbulence also speeds up the evaporation and it was postulated that the turbulence effects will be embedded in the friction velocity.

The model from Equation (18) suggests that the evaporation is influenced by geometry (because of the exposed surface to convective motion), convective transfer (turbulence and free stream velocity), and driving force (temperature and/or mole fraction). Based on analysis performed by similitudes, non-dimensional groups that lump all these effects and can influence the evaporation rate were derived. It was further postulated that the characteristic velocity representing the convective evaporation would be the friction velocity and not the traditional mean flow velocity. Based on these analyses, a geometrical factor of V/r2 defining the length scale was defined. Basically, a sessile droplet will “see” the friction velocity rather than the mean velocity of the flow. Furthermore, the turbulence effect that has an impact on evaporation rate is embedded in the friction velocity. Therefore, expressing the Reynolds number, which is the driver for convective evaporation in terms of the friction velocity, seems like a plausible proposition. Based on the similitude analysis that was explained earlier, a non-dimensional time scale is defined as:

tu*(V/r2)1n(1+B)=tu*(r2/V) n(1+B)E20

The model was used to calculate the evaporation rate and the normalized amount of mass left as a function of the non-dimensional time shown in Equation (20) for a variety of droplet initial volume, friction velocity, and initial contact angle (or base radius). The HD properties were used for this study. Figure 1.3.1 shows the collective effects of this study and it is observed that regardless of the droplet initial topology, air velocity, and turbulence intensity, a similar solution exists. This verifies the validity of our earlier postulate about the friction velocity and its ability to scale wind tunnel data to open air evaporation results. An identical experimental study was conducted with HD and the same similar solutions were calculated and plotted analogous to Figure 1.3.1 to create Figure 2.. It is observed that the experimental results verify the results obtained by the mathematical model.

Figure 2.

Similar solutions in the form of the percentage of mass left as a function of dimensionless time that is reduced with respect to geometrical, convective, and driving force parameters.

Figure 3.

Similar solutions for the available experimental data in the form of the percentage of mass left as a function of dimensionless time.

The significance of this result is in identifying the similar solution that lumps the effects of free stream momentum and turbulence (u*), the geometrical factor (r2/V), and heat or mass transfer (transfer number) into a single scalable variable. The variation in the initial contact angle makes a very small difference in the curve. Therefore, if experiments are conducted in different wind tunnels with different droplet sizes and different transfer numbers (liquids), the same evaporation rate curve should be obtained for all wind tunnels as long as the results are expressed in the non-dimensional time scales. In doing so, the friction velocity should be measured in each experiment. This is an important conclusion and should be directly tested to verify the results of Figures 2 and 3.

It is evident that the free stream velocity and turbulence intensity should be treated as independent variables, i.e., one can be altered while the other remains constant. A series of boundary layer measurement made at Caltech’s 6ʹ× 6 ʹ wind tunnel altered the turbulence intensity while the mean velocity was maintained constant. The turbulence intensity was altered by the number of bungee cords in the wind tunnel. A schematic of the wind tunnel and the cords for changing the turbulence level is shown in Figure 4. The wall shear stress was measured by the oil film technique in which a drop of olive oil is placed on a surface under a monochromatic light source. The combined reflection from oil and the glass surface create an interference image pattern as shown in Figure 5. The fringe spacing growth rate can be correlated to the wall shear stress as τw=2nμλwavedwdt where n is the index of refraction for oil, μ is the dynamic viscosity, and λwave is the light wavelength.

The turbulence intensity was measured by hot wire anemometry. Particle image velocimetry (PIV) can also be used to find the components of velocity fluctuations that will lead into the turbulence intensity calculations. Table 1 shows the wall shear stress as a function of the free stream velocity and turbulence intensity. The friction velocity (u*=τw/ρ) is calculated by finding the wall shear stress through interpolation from Table 1. The density can be obtained from the pressure and temperature data.

Figure 4.

Schematic of Caltech’s 6ʹ x 6 ʹ Wind Tunnel

Figure 5.

Fringe Spacing Growth Rate for Oil

In a group of wind tunnel experiments, HD droplet evaporation rates were measured for varying droplet volume, mean air speed, and temperature conditions. The turbulence intensity was estimated using the Edgewood Chemical and Biological Center (ECBC) wind tunnel data for the wall shear stress at the velocities that are mentioned hereinto. The wind tunnel tests were conducted for four different free stream velocities, u(m/s)={0.26, 1.77, 3.00, 3.66}, three initial droplet volumes, V(μL)={1, 6, 9}, and five different temperatures, T(°C)={15, 25, 35, 50, 55}. A numerical procedure that solves the governing Equation (17) along with Equation (18) was developed. The scheme is based on the Runge–Kutta fourth order integration algorithm. For the sake of brevity, the model predictions are compared with the experimental data for nine of approximately 60 cases in Figure 6. An excellent agreement is seen for all droplet sizes, air velocities, and temperatures. Detailed comparisons can be found in Navaz et al. [54, 55]

The analysis revealed that there is a difference in how each parameter influences the evaporation rate. Having the analytical formulation validated, the model can be used for parametric studies. To do this, the evaporation rate had to be calculated as the mean flow velocity and turbulence intensities were changed independently. A small (V0=1 μL) and a large (V0=9 μL) droplet were exposed to a constant mean boundary layer velocity as the turbulence intensity varied from 0 (laminar) to 6%. Then the mean velocity changed and the same calculations were repeated. Figure 7. (a and b) depict the results of this study and emphasizes the fact that that turbulence intensity is a driving force in convective cooling and the most viable method to account for its contribution is through friction velocity. A similar analysis was performed to examine the significance of the free stream velocity on the evaporation by considering the smallest and largest droplets as shown in Figure 7. (c and d). In this case, a free stream turbulence intensity of 2% with a constant prescribed free stream temperature of 35°C showed that if the free stream velocity is increased by a factor of about 14, the evaporation rate will decrease by a factor of about 4. By the same token, each free stream velocity can be correlated to the evaporation time.

In the last parametric study, the free stream velocity and turbulence intensity were maintained constant for the two droplet sizes and the temperature was varied. It appeared the air temperature (T) had the most significant effect on evaporation as seen from Figure 7. (e and f). For a 14% increase in absolute temperature, the evaporation rate is increased by a factor of 24 due to the vapor pressure being a strong function of temperature (the vapor pressure of HD at 30°C is about 0.25 mmHg).

There are two uncertainties that are encountered in the process of model validation with outdoor data: pertinent free stream velocity and turbulence intensity of the atmospheric flow. In the free atmospheric air flow, the wind velocity varies as a function of height (or vertical direction) and time, i.e., V(y, t). Furthermore, two velocities can be used as the ‘free stream velocity’ in the prediction of the droplet evaporation rate, one being the instantaneous wind velocity at the specific height V(yspec,t), and the other the height (y) averaged velocity V¯(t). In the outdoor experiments, the instantaneous velocity data were taken at seven or eight vertical locations up to approximately H=5 m high.

Figure 6.

Surface evaporation model prediction compared with wind tunnel experimental data for the following HD droplet sizes, wind speeds, and temperatures: (a) 1 μL, 1.77 m/s, 15°C, (b) 1 μL, 0.26 m/s, 35°C, (c) 1 μL, 1.77 m/s, 35°C, (d) 1 μL, 3.66 m/s, 50°C, (e) 6 μL, 1.77 m/s, 15°C, (f) 6 μL, 3.00 m/s, 35°C, (g) 9 μL, 3.00 m/s, 15°C, (h) 9 μL, 1.77 m/s, 35°C, and (i) 9 μL, 3.00 m/s, 55°C.

Free Stream Velocity (m/s) Turbulence Intensity
0% 0.3–0.4% 2.6% 4.1–5.4%
τw (Shear Stress at the wall, Pascal)
0 0.0 0.0 0.0 0.0
2 0.0037 0.0249 0.0251 0.0253
5 0.0147 0.0626 0.0663 0.0723
10 0.0466 0.2048 0.2102 0.2483
15 0.0765 0.4096 0.4450 0.5115

Table 1.

Wall Shear Stress Data as a Function of Turbulence Intensity and Free Stream Velocity

Figure 7.

Performing parametric study to explore the effect of turbulence intensity (a and b), free stream velocity (c and d), and air temperature (e and f) on evaporation.

The instantaneous velocity profile was used to find the transient average velocity, which is given by:


The same procedure was performed for the temperature profile and the graphical representations of these averaged values are embedded in all figures related to the outdoor data.

The second uncertainty encountered in model application is that the free stream turbulence intensity levels could be a function of time. In the atmospheric flow, one expects the existence of smaller turbulence intensity. In this work, the turbulence intensity for the outdoor experiments was assumed to be about 2% based on the “averaged” value of velocity measurements. Note that the turbulence intensity is an independent variable and needs to be extracted from the actual velocity measurement with high sampling rate.

Numerous model validations with open air test data are compiled by Navaz et al. [54, 55]. Two of these cases are presented here and an excellent agreement between outdoor data and model prediction are observed and summarized in Figure 8.(a and b).

Figure 8.

Model predictions for real-time wind velocity and temperature data as compared with the measurements for 1 μL average HD drop size. The embedded figures represent the instantaneous wind velocity and temperature.

1.4. Porous media models

The porous medium consists of solid inclusions with void spaces in between. The porosity is the ratio of the void volume to the total volume. It is an averaged property of a porous medium for “relatively” homogeneous substrates. For heterogeneous substrates, the porosity becomes a local value expressed as ϕ =ϕ(x,y,z). Furthermore, in the presence of chemical reaction, the local porosity changes in time as any of the reaction products may appear in solid phase making the porosity a function of time (or local solid phase density), i.e.. ϕ =ϕ(x,y,z,t) =ϕ(x,y,z,ρsk). For the time-dependent porosity, the gradient of the porosity should be included in the governing equations. This will make the equations more complex. However, since the reactions are generally slow, the porosity distribution and change in time is lagged by one time step (taken from the previous time step). This makes the algorithm less complicated and still very accurate.

Another property of the porous medium is the saturation permeability (K in Equation 10) that needs to be known. It is a momentum transport property of the porous medium. It is basically a measure of how “easily” fluid flows through fully saturated porous material. Therefore, this permeability is referred to as single-phase permeability, saturation permeability, or simply permeability. An experiment can be performed with any nonvolatile liquid including water to measure this property. [54, 63].

We have measured the saturation permeability for different porosities using glass beads. This was a necessary step for nonhomogeneous porous media. Saturation permeability is a function of porosity and changes locally as the porosity will change, especially due to those chemical reactions that produce an additional solid phase as a product. This relationship is:


where K is the saturation permeability in m2or ft2 and ϕ is the instantaneous local porosity (0.3, 0.5, etc.). This equation is general and valid for porosities up to 60% or 0.60. This procedure is detailed by Zand et al. [63]

Contrary to the single-phase (fully saturated) flow, multiphase flow can occur when a fluid phase prefers a specific path, thereby leaving some of the porous medium regions filled with the originally present phase, e.g., vapor–air mixture. The measure of how “easily” a phase flows through the porous medium in the multiphase flow is referred to as phase permeability. The relative permeability for each constituent in liquid and gaseous phase (ki,kgj) is defined as a ratio of the phase permeability to the saturation permeability and 0ki and kgj1 (zero for immobile phase and one for fully saturated flow). We eliminate the index “i” and “j” for simplicity, keeping in mind that these properties belong to each constituent in the liquid and/or gas phase. Hence, it is evident that the relative permeability is a function of the phase content or saturation, i.e.,k=k(s)   and  kg=kg(sg).

A hybrid approach (combined modeling and experimental studies) was used to find expressions for relative permeability. [54, 55] This hybrid approach needs to be performed separately for agents and the substrates of interest.

k=sk,     k=2kg=1+s2(2s3)E23

The capillary pressure function (see Equation 10) is not known for sessile droplets. Therefore, it is postulated that it has the following form:

pci=fi(*)σicosθiK/φJ(si),  where  J(si)=1.417(1si)2.120(1si)2+1.263(1si)3E24

where σi is the surface tension, θi is the contact angle inside the pores, K is the saturation permeability, siis the saturation for each constituent, and fi(*) is a function to be determined. This form is based on the original equation proposed by Leverett [13] and Udell [14]. In Equation (24), a geometrical scale is defined as K/φ, which is not necessarily the correct value as capillary pressure can be scaled with a different geometrical scale. Therefore, the capillary pressure is corrected for a potential error in geometrical scale using the functionfi(*).

Initially, the Buckingham theorem was used to identify the basic non-dimensional groups and the computational model developed was utilized to combine these basic Π groups into a more meaningful similar solution. A value of fi(*)=1 was assumed and numerical simulations were carried out over a wide range of initial droplet volume and wetted area (maintaining a constant initial contact angle), porous medium saturation permeability, liquid viscosity, and surface tension. For a constant initial contact angle and a specific substrate, a similar solution can be obtained with the following Π groups:

Π1=function(Π2)   for a constant contact angleΠ1=tσK0.5μ(Von  the  surfacer2)2                           Π2=Von  the  surfacer3E25

It should be noted again that the properties must be calculated for each constituent in the liquid phase marked by subscript “i” earlier. In Equation (25) r and V represent the wetted radius and the volume of the droplet left on the surface, respectively, and t is the time taken by a droplet to disappear from the surface. Figure 9.(a) shows a family of similar solution for these non-dimensional groups. Although this curve proves that the proposed non-dimensional groups are correct scalable quantities, it becomes a singular function when Von  the  surface0. Note that the properties relate to each liquid constituent on the surface only.

The above functions were used in an alternate method to eliminate this singularity. Again, the value f(*)=1 was kept, and the porosity and contact angle were changed to generate a family of curves with returning the instantaneous volume and wetted radius to their initial values, i.e., VInitial and rInitial. Figure 9.(b) is a collection of modified Π1 group as a function of contact angle for substrate porosity ranging from 0.1 to 0.9. The following procedure can be used to find the unknown f(*). The details of this method and all the experimental results and validations are compiled and discussed by Navaz et al. [64]. The process of finding f(*) is as follows:

  • Measure the initial contact angle, and based on the value of porosity, find the modified Π1 group representing tref (reference time).

  • Determine the function f(*) based on the following equation:

1fi(*)Unknown=nondimensional  measured time for the droplet to disppear from the surfacenondimensional  reference time (Π1  Read from Figure 9b)E26
  • Using the calculated value from step 2 should yield the correct f(*) value that can be used in Equation (24) for simulating any other scenario.

Figure 9.

Nondimensional groups and the similar solutions for the Π groups (a), and the nondimensional reference time as a function of the initial contact angle and porosity (b).

Table 2. shows the results of this study for five different liquids (including HD and VX, known as Mustard and nerve agents, respectively). More validations for this method is compiled by Navaz et al. [54]

Substance Experiment (time, s) – Extracted from Video Model (time, s)
1-2 Propandiol 1.45 1.50
Castor Oil 27.50 28.30
Glycerin 12.30 11.95
VX 1.50 1.45
HD 0.90 0.87

Table 2.

Measured and calculated time for the disappearance of a droplet from the surface of a porous material

Another transport property that needs to be known is the effective diffusion coefficient for the gas phase as indicated in Equations (5) and (6). The evaporation of a liquid inside a porous material generates a vapor phase concentration profile throughout the porous substrate, and the transport is described by a specie diffusion equation. This equation requires the diffusion coefficient as a transport parameter. In a porous medium, the molecular diffusion coefficient needs to be reduced as the vapor transport is affected by a presence of solid and liquid phases. This results in a definition of an effective diffusion coefficient, which is a function of porosity and saturation, as shown below:


The molecular diffusion coefficient as a function of temperature and molecular structure was taken from Treyball. [65] The function that describes how porosity and saturation attenuate this quantity was derived experimentally and is detailed by Navaz et al. [54, 55]

The effective diffusion coefficient, as given by Equation (27), is measured on a porous substrate bed where the boundary concentrations are given and transport takes place in the principal direction L. With reference to Figure 10., if the vapor is traveling from the lower surface to the upper surface (or vice versa), a one-dimensional diffusion problem governs this transport. The measurement technique is based on the gas chromatography and mass spectroscopy (GC/MS) method when the steady state of one-dimensional diffusion is achieved. At steady state, a constant evaporation rate is achieved, and therefore a constant slope mass loss or zero derivative of mass loss in time (m˙=dm/dt) should be observed. Then the equation in Figure (23) is used to find Deff.

Figure 10.

Experimental concept for the effective diffusivity measurements.

After the data reduction, the effective diffusion coefficient (Eq. 28) is obtained.


where DM is the molecular diffusion coefficient.

1.5. Contact dynamics and governing equations for the liquid bridge

Many parameters affect the amount or the concentration of transferred liquid into a secondary surface coming into contact with a sessile droplet on a substrate. The approach velocity between two surfaces affects the amount of mass transfer due to the resulting footprint after the contact. The mass transfer between the two surfaces, after the contact, occurs through the area of this footprint. There are basically two pathways by which a contact transfer can occur. These are schematically shown in Figures 29 and 30.

Figure 11.

Pathway 1: (a) Spread into the bottom surface, (b) contact dynamics, and (c) contact transfer

Figure 12.

Pathway 2: (a) Spread into the bottom surface, (b) contact dynamics, and (c) contact transfer

In this document, the bottom and top surfaces are referred to as the primary and contacting or secondary surfaces, respectively. The model is purely driven by the physiochemical properties of the medium, solids, liquids, and gases. When a droplet is placed on a primary surface (Figures 11a and 12a), capillary transport, surface evaporation (due to wind and temperature), and chemical reaction with the surface or other preexisting chemicals in the pores (such as moisture) are simultaneously initiated and the chemical available for transfer through contact is the end product or resultant of all these processes. We refer to this phenomenon as the precontact process. After this point, two possible scenarios could occur: creation of a liquid bridge between the two surfaces and its spread as the two surfaces get closer (Figure 11.b), or complete disappearance of the liquid droplet from the surface of the substrate (Figure 12.b). The last stage or post-contact period is when both surfaces come together and capillary transport occurs at the interface. The mass transfer to the contacting surface always occurs regardless of which pathway is taken. It should be noted when the upper surface comes into contact with or without the formation of the liquid bridge, chemical reaction, and evaporation can also be initiated in the upper surface.

The upper surface moves according to Newton’s second law assuming constant acceleration. The y-coordinate of the lower boundary of the upper surface (in motion) can be calculated by the equation of motion:

Equation  of  Motion:    y=F2mt2+Vot+yoy=Coordinate system (Vertical)        Vo=Initial Velocityyo=Initial space between the two surfacesF=Force exerted on the upper surface to move it down, m=mass of the upper surfaceE29

The liquid bridge has some free surface. It has been reported by Coblas et al. [66] that the shape of the free surface and its evolution varies as a function of inertial effects. For hydrophilic liquid bridges formed upon contact with small droplets (that is usually the case in chemical agent contamination and spray of pesticides on common surfaces), the liquid bridge is influenced by capillary flow and it generally assumes the shape of a hyperbloid when the velocity of approach is not significantly high (Longley et al.) [67]. The shape of the free surface of the liquid bridge (bulging in or out) has a very minimal effect on the amount of mass transfer into the contacting surfaces. Rather, the footprint area of the liquid bridge defines the surface available for the mass transfer. In the current work, this footprint is assumed to be in the shape of a circle, although the algorithm developed here may be modified to explore other shapes.

The contact angle of this hyperboloid and the surfaces can be found through experiments, although they will not affect the overall transfer of the mass. We have taken this angle to be π/8. The mass transfer into the porous media is calculated in a time step. Then the remaining mass of the liquid bridge is updated and its new volume is calculated. The separating distance is known from the equation of motion. By knowing the volume and the height of the liquid bridge, the contacting surface area (assumed to be circular) is calculated. The radius increase in time will determine the spread rate. Figure 13. shows that the arc FMG revolves around the C–Cʹ axis to produce a hyperboloid. We need to find the volume of this body of revolution analytically. If point A is the centroid of the surface engulfed by the arc FMG, the volume of the partial torus shaped resulting from this revolution according to Pappus centroid theorem is:

Volume of the partial torus = Vpt=2πAOA¯ where A is the area of the arc. Then the volume of the liquid between the two surfaces will be:


where t represents the time. Note that both the footprint radii and the liquid bridge height both are functions of time and are calculated in each time step. Ifθ^=GBF, then from basic trigonometry:

Figure 13.

Schematic of the hyperboloid geometry


The equation correlating instantaneous r(t) and h(t) with the liquid volume Liquid  Bridge(t) can be obtained as:

r2(t)h(t)G(β)4cos2βr(t)h2(t)G(β)tanβ8cos2β+h2(t)6Liquid  Bridgeπh(t)=0G(β)=π2βsin(π2β)E32

β is the contact angle between the surface and the liquid bridge. The droplet base radius (footprint), r(t), is calculated from the above equation. Both domains, or media, are remeshed and adapted to the new footprint of the liquid bridge. The remeshing is done after each time step.

As the gap between the two surfaces is closing, the height will reduce. At the same time, the liquid is absorbed into the porous media. Therefore, h(t) and Liquid  Bridge(t) are updated in each time step. Then the new footprint radius is obtained by Equation (32). It should be noted that when the height of the liquid bridge becomes less than a minimum threshold, the spread is stopped. This minimum threshold is a function of surface roughness and viscosity of the fluid. However, it is quite possible to assume a value of 10–100 μm for this minimum height. At this point the radius of the footprint remains constant and the absorption is driven solely by capillary pressure causing a reduction in h(t). This value of h(t) can be calculated by solving the above equation again with the “final” liquid footprint radius remaining a constant. This equation is given as:

(π6πG(β)tanβ8cos2β)h3(t)πr(t)G(β)4cos2βh2(t)+πr2(t)h(t)Liquid  Bridge=0E33

This will be a cubic equation in h(t). As the absorption into the porous media occurs, the volume of the very thin liquid bridge is updated and the new h(t) is calculated. In our model, this switching occurs when h(t) < 0.08 mm. This value is a function of surface roughness and can be specified by the user.

1.6. Chemistry model

The chemistry model is based on the Joint Army–Navy–NASA–Air Force (JANNAF) standard methodology. This methodology has three components: a reaction rate processor, species production destruction term calculator, and NASA thermodynamic file (provides specific heat, enthalpy, entropy, and Gibbs free energy) as a function of temperature. [68, 69, 70] This is shown in Equation (34). The outcome of this section will provide the term ω˙Reaction in the governing equations.


1.6.1. Bidirectional reactions

A general chemical reaction can be written in terms of its stoichiometric coefficients νij and νij' as:


where M¯i represents the i-th chemical species name (i=1, 2,...., NSP) and j represents the j-th reaction (j=1, 2,..., L). These reactions proceed (forward or reverse) according to the law of mass action, which states: “The rate at which an elementary reaction proceeds is proportional to the product of the molar concentrations of the reactant each raised to a power equal to its stoichiometric coefficient in the reaction equation.”

Let [M¯] denote the molar concentration of species i. The forward (left to right) reaction rate, Ratej (LR), for reaction j can be written as:


where kfj is the forward reaction rate constant. For species i, νijmoles on the left side of the reaction becomes νij' moles on the right side of the reaction. Consequently, the forward reaction for reaction j yields a time rate of change in the molar concentration of species i as follows:


Similarly, the reverse reaction for reaction j yields:


where kbj is the backward reaction rate constant. Thus, the net rate of change in the molar concentration of species i for reaction j (denoted by Xij) is as follows:


The “species production rate” is the time rate of change for the species density. For reaction j the net species production rate for species i is MwiXij where Mwi is the molecular weight of species i. Summing over all reactions gives the net species production rate ωi for the reaction set


The molar concentration of species i can be written as ρiMwi, which is the species i mass density divided by the species i molecular weight.


It follows that in terms of species i mass fraction Equation (39) becomes:

Xij=(νij'νij)[kfji=1NSP(ρCiMwi)νijkbji=1NSP(ρCiMwi)νij']         orXij=(νij'νij)[Kji=1NSP(ρCiMwi)νiji=1NSP(ρCiMwi)νij']kjKj=kfjkbj      and      kj=kbjE42

The reaction rate kj is from right to left (reverse) in the above equation and is often represented by the Arrhenius form:


where aj is the preexponential coefficient, nj is the temperature dependence of the preexponential factor, and bj is the activation energy. The ratio of forward to backward rate, Kj in Equation (42) is related to the equilibrium constant, Keql,by the following equation:




The quantity |λj|+1 is known as the order of the j-th reaction and the equilibrium constant is:

Keql=exp(ΔFR¯T),              ΔF=i=1NSPfiνiji=1NSPfiνij'fi=Gibbs Free Energy = Chemical Potential = h¯iTS¯ih¯i=Molar specific enthalpyS¯i=Molar specific entropyE46


Xij=(νij'νij)[Kji=1NSPC¯iνijρλji=1NSPC¯iνij']kjρi=1NSPνijC¯i=CiMwi( Molar mass fraction)E47

1.6.2. Third body dissociation recombination reactions

Reactions involving a third body have a distinct reaction rate for each particular third body. Benson and Fueno [71] have shown theoretically that the temperature dependence of recombination rates is approximately independent of the third body. Available experimental recombination rate data also indicates that the temperature dependence of recombination rates is independent of the third body within the experimental accuracy of the measurements. Assuming that the temperature dependence of recombination rates is independent of the third body, the recombination rate associated with the k-th species (third body) can be represented as:


where only the constant akj is different for different species (third bodies). Assuming that the reference species (third body) whose rate is specified in the program input has index k = kref, we may write:

kkj=akjakrefakrefTnjexp(bjR¯T)=mkjkref(j)     for kth speciesmkj=akjakrefE49

For the k-th species (third body) the stoichiometric coefficient is 1 on both sides of reaction j giving:

Xij(for 3rd body k)=(νij'νij)[KjρCkMwki=1NSP(ρCiMwi)νijρCkMwki=1NSP(ρCiMwi)νij'] =(νij'νij)[Kji=1NSP(ρCiMwi)νiji=1NSP(ρCiMwi)νij']mkjkref(j)E50

Summing over all third bodies k we get the total Xij:

Xij=k=1NSPXij(for 3rd body k)=(νij'νij)[Kji=1NSP(ρCiMwi)νiji=1NSP(ρCiMwi)νij']Mjkref(j)Mj=ρk=1NSPmkj(CkMwk)=ρk=1NSPmkjC¯k=ρk=1NSPmkj[M¯k]E51

This is the standard production rate modified by the factor Mj, which is an effective third body molar concentration. In terms of C¯i we may rewrite Eq. (51) as:

Xij=k=1NSPXij(for 3rd body k)=(νij'νij)[Kji=1NSPC¯iνijρλji=1NSPC¯iνij']Mjkjρi=1NSPνijMj={ρi=1NSPmijC¯k   for reactions involving 3rd bodies1                      for all other reactionsE52

1.6.3. Unidirectional reactions

For unidirectional reactions the reaction only proceeds from left to right and has the assumed form


where the chemical species B¯ does not occur on the right side of the equation. Let B¯ have species index K-UNI. A unidirectional reaction is obtained by setting the backward reaction rate constant to zero, so that for reaction j:

Xij=(νij'νij)(ρCKUNIMwKUNI)kfjSo that:d[B¯]dt=βρC¯KUNIkfj    and       d[A¯]dt=αiρC¯KUNIkfjE54

2. Numerical scheme and algorithm

The governing equations are transformed into the computational domain by rewriting them into a vector form with U being the conserved variable vector, E, F, and H are the inviscid flux vectors, and Ev, Fv, and Hv the viscous fluxes in each three directions, respectively. Vector G contains all source terms.


The following transformation is used to cast the conservation equations in the computational domain of ξ, η, and ζ that correspond to x, y, and z.


According to the chain rule of differentiation:



ξx=yηzζyζzη|n|,       ξy=xηzζxζzη|n|,     ξz=xηyζxζyη|n|ηx=yξzζyζzξ|n|,     ηy=xξzζxζzξ|n|,        ηz=xξyζxζyξ|n| ζx=yξzηyηzξ|n|,       ζy=xξzηxηzξ|n|,     ζz=xξyηxηyξ|n|  E58

With the norm of:


Equation (55) will now be transformed into:

U˜t˜+E˜ξ+F˜η+H˜ζ+E˜vξ+F˜vη+H˜vζ+G˜=0        with  U˜=1JU, G˜=1JGE˜=1J(ξxE+ξyF+ξzH),  F˜=1J(ηxE+ηyF+ηzH),  H˜=1J(ζxE+ζyF+ζzH)E˜v=1J(ξxEv+ξyFv+ξzHv),  F˜v=1J(ηxEv+ηyFv+ηzHv),  H˜v=1J(ζxEv+ζyFv+ζzHv)  E60

The Jacobian of the transformation is:


The explicit form of all of the flux vectors are described by Navaz et al. [55,56]. The above sets of equations are solved by finite difference method on a structured mesh. The explicit Runge–Kutta integration algorithm was used to find the distribution of all variables (gas, liquid, and solid concentrations, capillary pressure, velocities, etc.) in time. The program can accept externally generated mesh by commercial software such as GRIDGEN™ [72] or can generate a grid internally for one or two sessile droplet(s) residing on a porous on nonporous substrate. In the presence of a liquid bridge and its spread, a built-in and automated adaptive mesh generator is embedded in the computer model. The mesh or grid point distribution is based on Coon’s Patch [73]. First, the mesh for the substrate is generated and then the mesh defining the droplet(s) is overlaid on the substrate geometry. The mesh points defining the droplets are collapsed on the surface of the substrate except in the area(s) that droplet(s) reside on the surface. Figure 14 shows the top view of the adaptive mesh when liquid bridge expands during the contact process. Figures 11 and 12 shows snapshots of the absorption process during the contact. It displays the precontact and postcontact configurations.

Figure 14.

Top view of the adaptive mesh when liquid bridge still expanding, until the conditions for Equation (33) are reached.

The calculations start with the specified initial conditions and then the continuity and momentum equations are numerically integrated in time to find the distribution of liquid inside the pores. The saturation is equal to unity (s=1) at the interface between the droplet or liquid bridge with the porous media. The hydrostatic pressure is also added to the capillary pressure for more accurate calculations. That is, P=PPci+ρgh*where h* is the local height of the droplet above the surface. Mass is being transported into porous medium according to (ρυ˜φ)/J; where J is the Jacobian for the transformation and υ˜ is the contra-variant vertical velocity given by: υ˜=ηxu+ηyυ+ηzwwith u,υ,w being the three components of the velocity and ηx,ηy,ηz being the metrics of the transformation. The mass transfer is calculated in each time step and the instantaneous remaining mass yields the liquid bridge volume. The remaining mass divided by the liquid density will yield the instantaneous volume of the liquid bridge. Since the distance gap or the height of the liquid bridge is known from the equation of motion (Eq. 29), the temporal (or instantaneous) base or footprint radius can be found from Equation(32). When the distance between surfaces become less than a threshold, the spread is stopped, that is, the base radius stays constant at its last value and the capillary transport into both surfaces continues, causing a reduction in the height of the liquid bridge. At this point the instantaneous or temporal height of the liquid bridge can be calculated from Equation (33). This process continues until the amount of liquid between the two surfaces becomes zero.


3. Test cases

The computational model developed is called modeling of chemical hazard—dispersion and contact (MOCHA-DC). Test cases are designed to validate the physiochemical phenomena that are formulated in this chapter. More details are available in Navaz et al. [55, 56]. A test case describing the spread of VX nerve agent into sand is described. The capillary pressure, which is the main driver for the spread and transport of liquid into porous media, was obtained by the method outlined in Section 1.4.

3.1. Capillary transport and secondary evaporation

The topology of the penetration of a 6 μL VX droplet into the UK sand (ϕ = 0.44) after about 120 minutes is shown in Figure 15. VX has a very low vapor pressure at room temperature and so evaporation could be ignored. The shape of the liquid front is obtained experimentally at ECBC and is compared to the MOCHA-DC code prediction (saturation contours) in Figure 15. A comparison for the prediction of the depth of spread in time is observed from Figure 16. as compared with experimental data.

Figure 15.

Comparison of COMCAD prediction and experiment for VX on UK sand.

Figure 16.

Progression of depth in time for VX spread/transport through UK sand as predicted by the model and experiment

A similar validation was performed for the penetration of 1,2-propandiol, glycerin, and castor oil into a ceramic tile (ϕ = 0.26), coarse sand (ϕ = 0.31), and glass beads (ϕ = 0.28). The penetration was recorded after 2, 3, 5, 10, 20, 30, 60 minutes, for up to a week (for castor oil). The model prediction was excellent and they are detailed by Navaz et al. [55]. A sample case for glycerin spread into sand is shown in Figure 17.

Figure 17.

Penetration of glycerin into the coarse sand after 3 minutes (depth ~ 6 mm) and wetted pad from the experimental study

The secondary evaporation of several nerve agents are measured by Navaz et al.55. The validation of the model was performed using data from Reis et al. [74] and Mantle et al. [7] They studied the evaporation of water inside porous media using magnetic resonance imaging (MRI). Figure 18. shows an excellent comparison for the evaporation of water in sand and 120- and 400-μm diameter glass beads.

Data collected in collaboration with Kiple Acquisition Science Technology Logistics & Engineering Inc. (KASTLE), B.O.I.S.-Filtry, Brno, Czech Republic and VOP-026 Sternberk, Brno, Czech Republic, was used to validate the model. The data was taken in open air and the model was run under the same conditions (wind speed and temperature) for EDB (ethylene dibromide). An excellent agreement between the model prediction and data exists, as can be seen from Figure 19.

Figure 18.

Model validation for the secondary evaporation of water inside sand and glass beads compared with MRI data

Figure 19.

Model/experiment comparisons with Czech Data for 6 µL Ethylene Dibromide (EDB) droplet spread inside porous media (sand and glass beads) with secondary evaporation.

Figure 20.

Model/experiment comparisons with Czech Data for 6 µL Malathion (MAL) droplet spread inside porous media with secondary evaporation.

Wind tunnel test were performed for the evaporation of Malathion in dry concrete and glass beads. There was no (or minor) evaporation of Malathion inside dry concrete or glass beads due to the low vapor pressure of the species, as seen from Figure 20. However, it will be shown that in the presence of moisture or water this picture will change.

3.2. Chemical reaction in porous media

In this section series of validation, test cases are presented that involve chemical reaction(s) in different phases.

We chose to simulate two cases for a liquid interacting with a solid. The first case was a 50 µL drop of sulfuric acid placed on a brick of solid sodium sulfide. The second case was a 50 µL drop of sulfuric acid on a mixture of 75% sand and 25% sodium sulfide. In both cases the samples were open to the air to allow the product gas to escape. Mass loss was monitored and rates of product formation/reactant degradation inferred. Earlier experiments documented that during the time scale for these reactions the sulfuric acid did not react with the sand (100% sand case) and adsorption of atmospheric water was insignificant. Therefore, the mass change was due to the reaction alone. The reaction stoichiometric equation is given below. The model simulation of the liquid on the solid brick of sodium sulfide and the sand/sodium sulfide mixture are given in Figures 21a and b. Both results are validated with experimental data.

Na2S (s) + H2SO4 (aq) → Na2SO4 (s) + H2S (g)

(sodium sulfide + sulfuric acid → sodium sulfate + hydrogen sulfide)

Figure 21.

Model/ experiment comparison for (a) 100% sodium sulfate and 50 μL sulfuric acid reaction, and (b) 3:1 mixture of sand: sodium sulfate and 50 μL sulfuric acid droplet

In the next test case, glass beads were chosen as the porous substrate and droplets of cyclohexanol and phosphoric acid were deposited about 1 mm apart. The reaction started when the two droplets diffused together due to capillary forces, and the reaction proceeded thereafter. Figure 22. shows the contours of remaining reactants and the amount of product (cyclohexene) as predicted by the MOCHA-DC code from the initiation of the process (diffusion and reaction). The overall destruction of cyclohexanol is similar to experimental data and is described in Figure 22..

Figure 22.

Development of reaction zone between cyclohexanol and phosphoric acid Modeled by the MOCHA-DC Code (from top left to lower right) t=0, t=12 s, t=144 s, t=3600 s.

Figure 23.

Comparison of destruction of cyclohexanol to experimental data

The reaction of chemical warfare agents with water is of special interest in this section. This is called a hydrolysis reaction. VX and Malathion are among the chemicals that react with water. Water can be present as moisture in sand, concrete, and soil, or can be added to a preexisting chemical agent in a dry porous medium as rain. The model is capable of simulating all possible scenarios. The hydrolysis of VX proceeds as multiple chemical reaction processes that result in several breakdown products. The evolution of these products and the spread of a droplet therefore provide a reasonable validation of the interaction of the reaction and diffusion functions in the code. A simulation of a 6 μL droplet of O-ethyl-s-[2-N,N-(diisopropylamino)ethyl] methylphosphonothioate (VX) spread into sand with porosity of 35%. The hydrolysis reaction rates were defined using those for VX in moist sand given in Brevett et al. [75, 76] Brevett et al.’s [75, 76] work indicates that VX degrades into several products in the presence of water. All of these reactions were taken as first order in VX, due to the reported 25-32 fold molar excess of water present in the experiments. The droplet size was taken from the same study. The porous field in the model extended beyond the region of droplet spread. All model boundaries were set as impermeable, as the experimental study involved degradation under sealed conditions. The quantity of VX and the appearance of breakdown products indicated by the models throughout the simulation were consistent with data from the experimental study. VX spread to a radius of 7 mm in the model, similar to the spread reported by Brevett et al. [75, 76] and Wagner et al. [77] The match between the amount of VXH+ in the experiment and in the model during the degradation process results from the model’s ability to handle multiple reactions. The rate of production/destruction of each chemical species is a local event, calculated at each node. Species concentration and, therefore, the reaction rate vary across the saturation region within the sand, yielding a complex mixture of breakdown products. If the fluid motion predicted in the model was inappropriate it would lead to unrealistic local concentrations, which would cause the local rates of production/destruction (based on local concentrations in the reaction model) to be inappropriate. The agreement between the model and the experimental data, therefore, demonstrates robustness in both the fluid motion and reaction modeling methods. The result of this study is shown in Figure 24.

Figure 24.

A 6 μL VX droplet undergoes hydrolysis in damp sand at 50°C. Both agent and breakdown product amounts present in the sand, as measured in Brevett et al. [75, 76], are characterized by the model

3.3. Contact hazard

In this section, some test cases are presented to illustrate the model prediction for a sessile droplet spread dynamics between two moving nonporous or porous surfaces. The presented test cases will examine the gas and liquid mass transfer and the effect of chemical reactions during the droplet spread process.

In the first case, the spread of de-ionized water between a nonporous and porous glass plates as a function of velocity of approach was examined. A linear stage actuator (Thorlabs™ LNR50SEK1) controls the distance between the two parallel plates and a load cell (Interface™ MB-LBF5) is designated for measuring the force acting on the solid body. Two optical windows (Newport™ 20BW40-30) are used as the upper and lower plates as shown in Figure 25.. A feedback control loop was created, connecting the load-cell force input, the actuator position, and speed inputs and outputting the upper plate acceleration and speed. This feedback enabled to simulate a freely moving solid body with a given mass, m1, reacting to external forces and the force applied by the liquid bridge connected to a solid body with mass m2. The plates were cleaned before each experiment using acetone, isopropanol, and nitrogen gas. The advancing and receding wetting angles were estimated from quasi-static force measurements. The experiment utilized a 200 μL deionized (DI) water as a sessile droplet deposited onto a fused silica optical window (an impermeable surface), and a porous glass (VykorTM 7930) was used as the permeable contacting surface. The physical properties of water are μ=0.001Pa.s,ρ=1000kg/m3,σ=0.072N/m. The porous glass had a porosity of 28% with permeability of 2.08 × 10–19 m2/s (values provided by manufacturer). The surfaces were brought into contact with constant relative speeds of 0.5, 1.5, 2.0, and 2.5 μm/s, and the liquid bridge radius was measured in the experiment and compared to that calculated by the model (Figure 26). In all cases, there was a rapid initial change in geometry upon contact of the liquid with both plates, and a transition from sessile droplet geometry to the hyperboloid geometry of the liquid bridge. For smaller approach velocities, the absorption into the porous medium proceeded faster than the radial spreading throughout the liquid penetration process, causing a steady reduction in the footprint radius in time. However, for faster approach speeds (2 μm/s and above) the radius of the liquid bridge was only reducing in the initial stages due to rapid penetration into the porous medium. After some penetration of the liquid into the porous medium, and thus increasing the viscous resistance of the flow, the liquid began to spread radially and the radius of the liquid bridge increased with time. This finding demonstrated that the depth of liquid penetration and the imprint radius of droplet into the porous medium is determined by the speed of the approaching surfaces. We have also modeled two additional approach velocities of 1.0 and 3.0 μm/s to demonstrate the consistency of the results (Figure 26.). The model and experimental results indicated that the spreading increases as the approach velocity increases. This causes a larger footprint to be issued by the process.

Figure 25.

The experimental setup: two fused-silica plates connected by a DI water liquid bridge. (a) Plates in fixture and (a*) a closeup view of the liquid bridge (b) shows the porous glass used in the experiment

Figure 26.

200 μL of water on glass contacted by porous glass. The approach velocities of 0.5, 1.5, 2.0, and 2.5 μm/s are compared with the model. Two additional velocities of 1.0 and 3.0 μm/s are also modeled. For lower velocities the absorption into porous glass is at a faster pace than the spread speed, causing a continuous reduction in radius in time. As the approach velocity increases, the spread will proceed at a faster pace than the absorption, causing a continuous increase of radius in time. It appears that at velocity of 1.5 μm/s both processes take place at the same “speed”

An experiment conducted to measure the amount of mass transfer for two surfaces coming into contact with a liquid droplet in between. A 20 µL glycerin droplet was deposited on a porous media composed of play sand (porosity = 35%). Then a cloth at an initial distance of 11 cm with the speed of 35 cm/s was brought into contact with the glycerin. The amount of mass absorbed during the contact process into the cloth was measures and is compared to the model prediction. These tests were conducted with three additional chemicals and the results were compared to the model’s response. These conditions are shown in Figure 27.a and the agreement between the model and experimental data is seen to be excellent.

Another experiment was conducted with the same setup with the exception of varying the time of the contact. The purpose of this experiment is to demonstrate that even after the disappearance of the sessile droplet from the surface, mass transfer to a secondary porous material will take place. That is to say, even though a hazardous liquid substance may not be visible on a surface of a porous material, it can still pose a threat through contact. However, kitchen tile (porosity of 24%) with a lower porosity was selected for this experiment to have more control over the disappearance time of the sessile droplet from the surface. A wafer, fabricated in our laboratory from filter paper pulp, with a porosity of 60%, served as the secondary surface. Five separate experiments (three repetitions for each experiment) were conducted with a contact time to be, respectively, 1, 10, 20, 30, and 40 minutes after the droplet is placed on the kitchen tile. Every experiment was repeated three times at an initial gap of 2.5 cm between the two surfaces and a downward force of 1 N. The two surfaces were brought into contact. The amount of glycerin transferred into the secondary surface (wafer) was measured and compared favorably with the model predictions (Figure 27.b). These comparisons indicate that the model was fairly accurate in predicting the amount of mass being absorbed into a contacting porous surface.

Figure 27.

(a) Mass transferred to contacting surface as a function of the type of chemical and comparison of the model prediction with experimental data for each case, and (b) mass transfer of glycerin, that was initially deposited on sand, to a contacting cloth, where the amount transferred was a function of the time elapsed between deposition and the contact, the figure also shows the model prediction for each case.

In order to provide data for cases where contact between two surfaces results in a chemical reaction in one or both of the contacting surfaces, experiments were conducted measuring the rate of formation of product where one reactant resides on one surface and the other on the second, contacting, surface. Since there are many permutations for these types of experiments, we began our investigation for a case where the primary surface is nonporous, a glass slide, and the contacting surface is porous, an adsorbent pad. The reaction needed to be fast enough so that the product formation could be measured accurately and in a relatively short time span to avoid complications associated with evaporation of chemicals. The reaction chosen was addition of bromine to alkenes, since the progress of this reaction could be monitored both quantitatively and qualitatively. A 20 μL of styrene solution was placed on the surface of a glass slide and an adsorbent pad saturated with 1 mL of bromine in carbon tetrachloride solution was placed on top of the slide. The contact was stopped at various time intervals. The reaction was quenched by addition of excess cyclohexene to the pad, after removal from the slide, and adsorbent pad and the glass slide were each extracted with 2 mL of dichloromethane. The solutions were then analyzed by GC/MS and GC/FID. No traces of styrene or the addition product were observed on the glass slide even for the smallest time step, indicating complete adsorption into the pad. Figure 28. shows the model/experiment comparison of the product form according to:

The styrene–bromine reaction is as follows: C8H8(L)+Br2(L) = C8H8Br2(L)

Figure 28.

Product formation in time for the styrene and bromine reaction from the experiment and as predicted by MOCHA-DC and the model showing the styrene droplet just prior to contact with the bromine saturated upper layer

The second experiment is similar to the previous experiment but the glass substrate is replaced with permeable sand. A 60 μL styrene droplet was placed on sand and after completely disappearing from the surface was contacted by filtered paper saturated by bromine. The mass exchange occurs at the interface such that the reaction takes place in both media. The results are shown in Table 3. The comparison is quite satisfactory and the model accurately predicts the mass transfer between the two contacting surfaces and the amount of product in both zones.

Styrene Sand (mg) Styrene Pad (mg) Product Sand (mg) Product Pad (mg)
MOCHA-DC 44.50 7.12 4.90 1.75
Data 43.74 7.40 5.13 1.56

Table 3.

Data and model reactant and product mass comparison in two porous media

In another test case, a 6 μL VX droplet is situated on moist sand (5% saturation). The reaction occurs until a secondary moist porous medium (cloth) is brought into contact with the first surface after about 0.5 hour. The remaining VX at the surface of the sand starts to react with water content in or from the contacting surface. The reaction in both zones continues for about 3 hours. Figure 29. shows the amount of VX left and all products during this process in real time. We do not have any experimental data to compare with these results, however; although there are no data available for this case, it is a very challenging case to run and it demonstrates the robustness and generality of the code.

Figure 29.

VX on wet sand and contact with another wet porous surface

In the following test case, liquid and vapor transport through a dual layer fabric is examined when each of the layers has a different permeability. The layer with small permeability (pores of 0.1 μm) is defined hereafter as the outer later and the layer with the larger permeability (pores of 5 μm) is defined as the inner later. The first experiment started with a 40 μL de-ionized water droplet on the inner-layer. Figure 30.a shows the saturation distribution and its migration to the outer layer. In the second experiment the droplet was placed on the outer layer and it did not penetrate to the inner layer as seen in Figure 30.b.

Computer simulations were performed with the same conditions. The evaporation module was activated to monitor the vapor concentration as well as liquid saturation. The MOCHA-DC prediction is shown in Figure 31. for the case of a droplet deposited on the low-porosity outer layer. The model indicates the liquid spreads laterally, instead of continuing to spread vertically, when it reaches the interface with the high-porosity inner layer. Due to the smaller pores, the capillary pressure at the outer layer is greater than that of the inner layer. It is seen that this spread pattern virtually eliminates the possibility of liquid passing through the higher porosity inner layer. Since the double layer material is thin, the aspect ratio of vertical direction in these plots is taken to be eight times greater than other directions to make the visuals clearer.

On the other hand, when the droplet is placed on the inner layer as seen from Figures 30a and 32, its preferred direction for movement is vertical when it reaches the interface between the layers, rather than spreading horizontally. In the top image of Figure 32., the droplet is initially placed on the inner layer and gradually is transported vertically, due to capillary pressure, toward the outer layer. The liquid keeps moving upward toward the outer layer saturating the upper portion of the fabric. The fluid therefore travels through both the layers. This behavior would facilitate the transport of fluid, such as sweat, through a garment constructed of such a material. As inner layer fluids reach the outer layer, they might evaporate in the surrounding environment. Conversely, when fluids contact the outer layer they do not reach the inner layer surface, thus preventing contamination of the inner surface. This result suggests it may be possible to eliminate the likelihood of hazardous liquid contact with skin if such a fabric were used for protective clothing.

Figure 30.

Sessile droplet positioned on the inner layer (a) and the outer layer (b) on a double-layered fabric. The outer layer pore size is order of magnitude smaller compared with the inner layer. (Adapted from Gat et al. 2013.)

Figure 31.

Liquid transport and spread from the upper layer to the lower layer in the model. This time series shows the saturation contours in the center plane of the domain as the droplet is absorbed and moves and spreads laterally, but not into the lower layer of the dual layer fabric. The saturation within the lower layer remains zero.

Figure 32.

Liquid transport and spread from the lower layer to the upper layer as progressed in time. The figures show the saturation contours at the center plane of the domain. The top picture is the initial condition and as the droplet is absorbed it moves toward the upper layer rather immediately and then spread laterally.

The vapor phase also followed the same pattern. Although a small amount of vapor was predicted at the inner layer, it was significantly smaller than the amount of vapor present in the outer layer. Figure 33. shows the distribution of mole fraction of water vapor inside both layers. The vapor inside the bottom layer is one to two orders of magnitude less than the upper layer, implying that the capillary pressure difference between the two layers can be utilized as a tool to prevent volatile chemicals from reaching skin.

Figure 33.

Image of the vapor transport and spread in the upper and lower layers when a droplet is placed on the upper surface. The vapor phase demonstrates a behavior similar to that of the liquid phase. The droplet is 40 μL water. The evaporation occurs somewhat faster at the edges of the wetted region. The mole fraction of the vapor phase is shown.

Given that the pore size influences capillary pressure, a range of porous glasses were coupled to study the transport behavior. Figure 34. demonstrates the experimental and model prediction for all combinations of the two pore sizes for the evaporation of deionized water placed on the bottom surface of the dual layer.

Figure 34.

Evaporation vs. time for 200 μL ethanol droplet, positioned on the inner layer, into various combinations of nylon membranes (Scientific Tisch, TM) with 5 and 0.1 µm average pore size. Both model and experiment predict smaller evaporation time for the smaller inner pore size. For all cases the droplet was placed on the lower surface.

When the outer (upper) layer’s pore size is an order of magnitude smaller than the inner (lower) layer’s, the time for evaporation of a droplet on the lower surface was reduced by about 3040%. This suggests that such a combination can be used to increase the liquid evaporation rate from the inner side of a protective fabric so as to enhance the pass through of perspiration. For the smallest outer layer pore size, the model predicts somewhat slower evaporation. This could be due to the mesh size, which may need to be refined for faster processes. However, since the results are fairly close to the experimental values, we decided not to repeat this calculation and emphasize the fact that there are a number of values to the dual layer fabric design that make it suitable for a variety of applications.


4. Nomenclature

4.1. General

As Surface area

a Pre-exponential coefficient

b Activation energy

C Mass fraction

cp Specific heat at constant pressure

D Diffusion coefficient

E, F, H Flux vectors in x, y, and z directions

F Force

f Gibb’s Free Energy

f(*) Unknown quantity for capillary pressure

G Source term vector

g Gravitational acceleration

h Height of apex of a sessile droplet

h¯ Molar enthalpy

h* Local height of the droplet

H Enthalpy

hfg Latent heat of evaporation

H Maximum vertical distance for boundary layer measurements

J Jacobian for the transformation

K Saturated permeability for isotropic medium

Keql Equilibrium reaction rate

k Reaction rate

keff Effective conductivity

kg,k Relative permeability in gas or liquid phases

M Molecular weight

Molar concentration

m Mass

n Index of refraction

P Pressure

Pv Vapor pressure

q˙ Heat flux

r Wetted area radius

R or Rs Radius of curvature

Re Reynolds number, uR/ν

S Entropy

s Saturation, fraction of the pore volume occupied by a phase (or state)

Sc Schmidt number

t Time

T Temperature

U Conserved variables vector

u* Friction velocity

V Velocity

v and v’ Stoichiometric coefficients

υ˜ Contra-variant vertical velocity

u,υ,w Velocity components

x, y, z Cartesian coordinate system

X Net rate of change in molar concentration

Greek Symbols

α Thermal diffusivity

β Angle defined in Figure 12.

ϕ Porosity

λ or λr h/r

λj or λij Reaction order

μ Viscosity

ρ Density

ρ˙ Rate of phase change

σ Surface tension

τ Shear stress


ω Species production/destruction rate

ϴ Contact angle

ϴa Contact angle at which the evaporation topology changes

ϴi Contact angle inside the pores

ξ, η, ζ Computational domain coordinates


A Air

b Backward

c Capillary

eql Equilibrium

f Forward

fg Fluid to gas (evaporation)

g Related to the gaseous phase

i Liquid phase constituent

j Gas phase constituent

k Solid phase constituent

, liq Related to the liquid phase

m,n, Refer to xyz for shear stress tensor

M Related to mass transfer

mix Mixture of liquid and gas phase

n Temperature dependence

o Initial

pt Partial torus

r Relative

s Surface of the droplet, or related to the solid phase

s Solid to liquid phase change

g Liquid to gas phase change

sg Solid to gas phase change

T Related to heat transfer

vap Vapor phase


w At the wall

Far field – Free stream


~ Vector

UNI Unidirectional

per unit time



This project was collectively supported by The Air Force Research Laboratory, Human Effectiveness Directorate, Biosciences and Protection Division, Wright-Patterson AFB, US Army's Edgewood Chemical and Biological Center under contract, and The Defense Threat Reduction Agency (DTRA) under the contract HDTRA1-10-C-0064. The authors wish to thank all the above agencies and our CORs Drs. William Ginely, Sari Paikoff, and Brian Pate for their support and guiding this effort. Special thank also goes to Mr. Joseph B. Kiple of KASTLE Corporation for his involvement with the open air testing and sharing his experiences and vision with the entire team and Dr. D’Onofrio of ECBC for providing some of test data. We also appreciate Dr. Miroslav Skoumal and his team in Czech Republic for participating and sharing their state-of-the art open air testing facilities and their laboratories. A final word of thanks goes to the many Kettering University students who helped with the experimental and numerical studies.


  1. 1. Davis, S. H., Hocking, L. M., (1999), “Spreading and imbibition of viscous liquid on a porous base,” Phys. Fluids, 11, 48-57
  2. 2. tarov, V. M., (2004), “Surfactant solutions and porous substrates: spreading and imbibition,” Adv. Colloid Interface Sci., 111, 3–27
  3. 3. Alleborn, N., Rasziller, H., Anthonissen, K., and Lievens, O. (2003), “Spreading and sorption of a droplet on a porous substrate,” Fifth Europian Coating Symposium, Fribourg, Switzerland, pp. 246–251
  4. 4. Alleborn N. and Rasziller, H. (2004), “Spreading and sorption of a droplet on a porous substrate,” Chem. Eng. Sci., 59, 2071–2088
  5. 5. Denesuk, M., Smith, G. L., Zelinski, B. J. J., Kreidl, N. J., and Uhlmann, D. R., (1993), “Capillary penetration of liquid droplets into porous materials,” J. Colloid Interface Sci., 158, 114–120
  6. 6. Holman, R. K., Cima, M. J., Uhland, S. A., and Sachs, E., (2002), “Spreading and infiltration of inkjet-printed polymer solution droplets on a porous substrate”,J. Colloid Interface Sci., 249, 432–440
  7. 7. Mantle M. D., Reis N. C., Griffiths R. F., and Gladden, L. F. (2003), “MRI studies of the evaporation of a single liquid droplet from porous surfaces,” Magnetic Res. Imag., 21, 293–297.
  8. 8. Bear, J. (1988), Dynamics of Fluids in Porous Media, Dover
  9. 9. Vafai, K. (2000), Handbook of Porous Media, Marcel Dekker, Chapter 17
  10. 10. Chen, Q. and Balcom, B. J. (2005), “Measurement of rock core capillary pressure curves using a single-speed centrifuge and one dimensional magnetic resonance imaging,” Journal of Chemical Physics, Vol. 122, 214720
  11. 11. Brooks, R. H. and Corey, A. T. (1964), “Hydraulic Properties of Porous Media,” Colorado State University Hydrology Paper No. 3
  12. 12. Van Genuchten, M. T. (1980), “A closed-form equation for predicting the hydraulic conductivity of unsaturated soils,” Soil Sci. Soc. Am. J., 44, 892–898
  13. 13. Leverett, M. C. (1941), “Capillary behavior in porous solids,” AIME Trans., 142, 152–169
  14. 14. Udell, K. S. (1985), “Heat transfer in porous media considering phase change and capillarity—the heat pipe effect,” Int. J. Heat Mass Transfer, 28, 485–495
  15. 15. Starov, V. M., Kostvintsev, S. R., Sobolev, V. D., Velarde, M. G., and Zhdanov, S. A., (2002), “Spreading of liquid drops over dry porous layers: Complete wetting case,” J. Colloid Interface Sci., 252, 397–408
  16. 16. Lenormand, R., Touboul, E., and C. Zarcone, (1988), “Numerical models and experiments on immiscible displacement in porous media,” J. Fluid Mech., 189, 165–187
  17. 17. Valavanides, M. S., and Payatakes, A. C. (2001), “True-to-mechanism model of steady-state two-phase flow in porous media, using decomposition into prototype flows,” Adv. Water Resour., 24, 385–407
  18. 18. Fatt, I. (1956), “The network model of porous media III. Dynamic properties of networks with tube radius distribution,” Trans AIME, 207, 164–181
  19. 19. Prat, M. (1993), “Percolation model of drying under isothermal conditions in porous media,” Int. J. Multiphase Flow, 19,691–704
  20. 20. Constantinides, G. N. and Payatakes, A. C. (1996), “Network simulation of steady-state two-phase flow in consolidated media,” AIChE J., 42, 369–382
  21. 21. Markicevic, B., D’Onofrio, T. G., and Navaz, H. K. (2010), “On spread extent of sessile droplet into porous medium: Numerical solution and comparisons with experiments,”, Phys. Fluids, 22, Art. No. 012103
  22. 22. Markicevic, B. and Navaz, H. K. (2010), “Primary and secondary infiltration of wetting liquid sessile droplet into porous media,” Trans. Porous Media, 85, 953–974
  23. 23. Markicevic, B. and Navaz, H. K. (2010), “The influence of capillary flow on the fate of evaporating wetted imprint of the sessile droplet in porous medium,” Phys. Fluids, 22, Art. No. 122103
  24. 24. Lago, M., and Araujo, M. (2001), “Capillary rise in porous media,” J. Colloid Interface Sci., 234, 35–43
  25. 25. Markicevic, B. and Navaz, H. K. (2009), “Numerical solution of wetting fluid spread into porous media”, Int. J. Num. Meth. Heat Fluid Flow., 19, 521–534
  26. 26. Popov, Y. O. (2005), “Evaporation deposition pattern: Spatial dimensions of the deposit,” Phys. Rev., E 71, 036313-1–036313-17
  27. 27. Hu, H. and R. G. Larson, (2002), “Evaporation of a sessile droplet on a substrate,” J. Phys. Chem., B, 106, 1334–1344
  28. 28. Baines, W. D. and James, D. F. (1994), “Evaporation of a droplet on a surface,” Ind. Eng. Chem. Res., 33, 411–416
  29. 29. Winter, S., Karlsson, E., and Nyholm, S. (1990). Models for the Evaporation of Chemical Warfare Agents and Other Tracer Chemicals on the Ground. A Literature Review, Department of NBC Defense Report, S-901 82 Umea, Sweden, Art Hin Fokke, Oeseburg: TNO, Prins Maurits Laboratorium, The Netherland..
  30. 30. Navaz, H. K., Dang, A. L., and Rangel, R. H. (1992), “Numerical analysis of bipropellant combustion in liquid thrust chambers by an eulerian-eulerian approach,” AIAA Paper No. 92-3735.
  31. 31. Ward, P., Collings, and Hay, N. (1985), “The comparison of simple models of turbulent droplet diffusion suitability for use in computation of spray flames,” J. Engng Gas Turbines Power, Transactions of ASME, 107, 690–694
  32. 32. Metzger, M. M. and Klewicki, J. C. (2001),“A comparative study of near wall turbulence in high and low Reynolds number boundary layers,” Phys. Fluids, 13(3), 692–701
  33. 33. Klinger, B. A., Huang, B., Kirtman, B., Schop, P., and Wang, J. (2004), Monthly Climatology of Friction Velocity Cubed from SSM/I Wind Data, COLA Technical Report
  34. 34. Weber, R. O. (1999), “Remarks on the definition and estimation of friction velocity,” Boundary-Layer Meteorology 93,197–209
  35. 35. Frandsen, S. and Madsen, P. H. (2005), Spatially Average Of Turbulence Intensity Inside Large Wind Turbine Arrays, Riso National Laboratory, 4000 Roskilde, Denmark
  36. 36. D’Asaro, E. A. and Geoffrey, G. T. (1997), “Turbulence intensity measurements in a wind-driven mixed layer,” J. Phys. Oceanogr., 27,2009–2022
  37. 37. Zarrouk, S. J. (2004), Simulation of Complex Multi-Phase, Multi-Component Reacting Flows in Porous Media, Ph.D. thesis, The University of Auckland
  38. 38. Zarrouk, S. J. (2008), Reacting Flows in Porous Media: Complex Multi-Phase, Multi-Component Simulation, VDM Verlag
  39. 39. Xu, T., Sonnenthal, E. L., Spycher, N., and Pruess, K. (2004), TOUGHREACT User's Guide: A simulation program for non-isothermal multiphase reactive geochemical transport in variably saturated geologic media. Report LBNL-55460, Lawrence Berkeley National Laboratory, Berkeley, CA
  40. 40. Lichtner, P.C. (1985), “Continuum model for simultaneous chemical reactions and mass transport in hydrothermal systems,” Geochem. Cosmochem. Acta49, 779–800
  41. 41. Steefel, C. I. (2001), GIMRT, Version 1.2: Software for modeling multicomponent, multidimensional reactive transport, User's Guide. Report UCRL-MA-143182, Lawrence Livermore National Laboratory, Livermore, California
  42. 42. Butt, H.-J. and Kappl, M. (2008), “Normal capillary forces,” Adv. Colloid Interface Sci.146, 48
  43. 43. Carter, W. C. (1988), “The forces and behavior of fluids constrained by solids,” Acta Metall.36, 2283
  44. 44. Mazzon, D. N., Tardos, G. I., and Pfeffe, R. (1987), “The behavior of liquid bridges between two relatively moving particles,” Powder Technol.51, 71
  45. 45. Dodds, S., Carvalho, M. D. S., and Kumar S. (2009), “Stretching and slipping of liquid bridges near plates and cavities,” Phys. Fluids,21, 092103
  46. 46. O. Pitois, O., Moucheront, P., and Chateau, X. (2000), “Liquid bridge between twomoving spheres: An experimental study of viscosity effects,” J. Colloid Interface Sci. 231, 26
  47. 47. Meurisse, M.-H. and Querry, M. (2006), “Squeeze effects in a flat liquid bridgebetween parallel solid surfaces,” J. Tribol.128, 575
  48. 48. De Sauza, E. J., Gao, L., McCarthy, T. J., Arzt, E., and Crosby, A. J. (2008), “Effectof contact angle hysteresis on the measurement of capillary forces,” Langmuir24, 1391
  49. 49. Ngo, M.A., O’Malley, M., and Maibach, H. I. J. (2010), “Percutaneous absorption and exposure assessment of pesticides," Appl. Toxicol., 30(2),91–114
  50. 50. Kilpatrick, W. T., Ling, E. E., Hin, A. R. T., and Brevett, C. A. S., Fagan, M. W., and Murdock, W. P., Jr. (2004), “Testing requirements for predictive model development using environmental wind tunnels,” (HEDs),” AFRL-HE-WP-TR-2004-DRAFT, Air Force ResearchLaboratory Report
  51. 51. Savage, J. J. (2006). “Agent Fate Program Overview,” 2006 Scientific Conference on Chemical Biological Defense Research, Hunt Valley, MD
  52. 52. Savage, J. J. (2006), “Agent Fate Program Overview,” 74th Military Operations Research Society Symposium (MORSS), Colorado Springs, CO
  53. 53. Munro, N. B., Talmage, S. S., Griffin, G. D., Waters, L. C., Watson, A. P., King, J. F., and Hauschild, V. (1999),“The sources, fate, and toxicity of chemical warfare agent degradation products, ”Environ. Health Perspec., 107, 933
  54. 54. Navaz, H. K., Zand, A., Markicevic, B., Herman, M., Atkinson, T., Li, Nowakowski, A., Gheresus II, P., Rothstein, M., Kiple, J., and Middleton, V. (2011) Agent Fate Model Integration, Final Report, Prepared for Edgewood Chemical and Biological Center (ECBC)
  55. 55. Navaz, H. K., Zand, A., Atkinson, T., Nowakowski, A., Kiple, J., Kamensky, K., and Jovic, Z. (2013), Agent Fate Modeling-COMCAD Theory and Analysis, Code Structure, and Users Guide (Volumes I, II, and III), Final Report prepared for Edgewood Chemical and Biological Center (ECBC), January 7
  56. 56. Navaz, H. K., Zand, A., Dang, A. L., Gat, A., Gharib, M., Atkinson, T., Nowakowski, Kamensky, K., and Jovic, Z. (2014), Predictive Model for Assessment of Chemicals on and in Surfaces vs. Chemicals Available for Contact and Transport—Volume I: Theory and Analysis,Final Report prepared for Defense Threat Reduction Agency (DTRA)
  57. 57. Navaz, H. K., Chan, E., and Markicevic, B. (2008), “Convective evaporation of sessile droplets in a turbulent flow—comparison with wind tunnel data,” International Journal of Thermal Sciences 47, 963–971
  58. 58. Navaz, H. K., and Dang, A. D. (1994),“The Development of the Liquid Thrust Chamber Performance (LTCP) Code for Turbulent Two-Phase Flow Combustion of Dense Sprays,” Technical Report prepared for NASA/MSFC, Contract No. NAS8-38798
  59. 59. Kuo, K. K. (1986), Principles of Combustion, Wiley, Chapters 6 and 8
  60. 60. Chiu, H. H. and Zhou, X. Q., (1983) “GEMCHIP-CODE: Spray Combustion Process in Air Breathing Propulsion Application,” University of Illinois at Chicago
  61. 61. Chiu, H. H. (1989), “Theory of Droplet: I Renormalization Laws of Droplet Vaporization in Non-Dilute Sprays,” Constitutive Relationships and Models in Continuum Theories of Multiphase Flows, NASA Conference Publication 3047
  62. 62. Nikolopoulos, N., Theodorakakos, A., and Bergeles, G., (2007), “A numerical investigation of the evaporation process of a liquid droplet impinging onto a hot substrate”, Int. J. Heat Mass Transfer, 50, 303–319
  63. 63. Zand, A., Sikorski, Y., Sanders, M., and Navaz, H. K. (2007) “A simple laboratory experiment for the measurement of single phase permeability,” J. Physical Nat Sci, 1(2), pp 1-10
  64. 64. Navaz, H. K., Markicevic, B., Zand, A., Sikorski, Y., Chan, E., Sanders, M., and T. D’Onofrio (2008), “Sessile Droplet spread into porous substrates—Determination of capillary pressure using a continuum approach,” J. Colloid Interface Sci., 325, 440–446
  65. 65. Treyball, R. E. (1980), Mass-Transfer Operation, McGraw Hill.
  66. 66. Coblaş, D., Broboană, D., Bălan, C., and Hejjam, M. (2013) “ Numerical Simulation of Constant Velocity Squeeze Flow,” U. P. B. Bul., Series D, Vol. 75, Issue 2
  67. 67. Longley, L. E., Dooley, E., Givler, D. M., Napier, III, W. J., Chaudhury, M. K., and Daniel, S. (2012) “Drop motion induced by repeated stretching and relaxation on a gradient surface with hysteresis,” Langmuir 28 (39), pp 13912–13918
  68. 68. Gordon, S., and McBride, B. (Oct. 1994) Computer Program for Calculation of Complex Chemical Equilibrium Composition and Applications, NASA Reference Publication 1311
  69. 69. McBride, B. and Gordon, S. (Nov. 1992) Computer Program for Calculating and Fitting Thermodynamic Functions, NASA Reference Publication 1271
  70. 70. McBride, B., Gordon, S., and Reno, M. (Oct. 1993) Coefficients for Calculating Thermodynamic and Transport Properties of Individual Species, NASA Technical Memorandum 4513
  71. 71. Benson, S. W., and Fueno, T. (1962) “Mechanism of atom recombination by consecutive vibrational deactivation,” J. Chem. Phys., 36, 1597
  72. 72. GRIDGENTM, (2009), Reliable CFD Meshing, Pointwise Corporation, Fortworth, TX
  73. 73. Coons, S. A. (1966), “Computer graphic and innovative engineering design—super-Sculptor,” Datamation,12(5) 32–34
  74. 74. Reis, N., and Griffiths, R. (2003), "Investigation of the evaporation of the embedded liquid droplets from the porous surfaces using magnetic resonance imaging," Int J Heat Mass Transfer
  75. 75. Brevett, C. A. S., Kenneth B., Sumpter, J. P., and Robert G. (2009), “Evaporation and degradation of VX on silica sand,” J. Phys. Chem.,113(16), pp 6622-6633.
  76. 76. Brevett, C, Sumpter, K, Pence, J. et al. (2009) “Evaporation and degradation of VX on silica sand”, J. Phys. Chem., 113(16), 6622–6633
  77. 77. Wagner, G. W., O’Connor, R. J., Edwards, J. L., and Brevett, C. A. (2004) “Effect of drop size on the degradation of VX in concrete,” Langmuir, 20, 7146-7150

Written By

Navaz Homayun, Zand Ali, Gat Amir and Atkinson Theresa

Submitted: October 22nd, 2014 Reviewed: May 12th, 2015 Published: December 16th, 2015