**Prediction of Solubility of Active Pharmaceutical Ingredients in Single Solvents and Their Mixtures — Solvent Screening**

Ehsan Sheikholeslamzadeh and Sohrab Rohani

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/60982

#### **Abstract**

In this chapter, the applicability of two predictive activity coefficient-based models will be examined. The experimental data from five different types of VLE (vaporliquid equilibrium) and VLLE (vapor-liquid-liquid equilibrium) systems that are common in industry are used for the evaluation. The nonrandom two-liquid segment activity coefficient (NRTL-SAC) and universal functional activity coefficient (UNI‐ FAC) were selected to model the systems. The various thermodynamic relations existing in the open literature will be discussed and used to predict the solubility of active pharmaceutical ingredients and other small organic molecules in a single or a mixture of solvents. Equations of states, the activity coefficient, and predictive models will be discussed and used for this purpose. We shall also present some of our results on solvent screening using a single and a mixture of solvents.

**Keywords:** Solubility prediction, Pharmaceuticals, , NRTL-SAC Thermodynamic model, Activity coefficient, Solvent screening, Single solvent, Solvent mixture

#### **1. Introduction**

The study of solutions and their properties is one of the important branches of thermodynam‐ ics. It is important to know the behavior of a mixture of components for a change in temper‐ ature, pressure, and composition. Knowledge in this area of thermodynamics helps engineers

© 2015 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

and scientists to better design, optimize, and operate the process units in the oil, gas, petro‐ chemicals, pharmaceuticals, agricultural, and other chemical-related industries. Knowledge of the thermodynamics is essential to improve process performance and product quality. For example, the use of phase behavior calculations to understand and estimate the production rate of solids from a crystallization process in a pharmaceutical industry is of paramount importance in modeling, optimization, and control of product quality.

A solution is called ideal if its mixture property is a linear combination of the properties of each of its constituents at the given temperature and pressure [1]:

$$m\_{\text{solution}} = \sum\_{i=1}^{N} \mathbf{x}\_i m\_i \tag{1}$$

where *msolution* and *mi* represent the molar property of the mixture and pure component *i*, respectively, and *xi* denotes the mole fraction of each constituent *i* in the solution. Not all solutions can be considered ideal. The interactions between the solute and solvent molecules in the solution renders a solution nonideal. The fundamental relationship which relates the thermodynamic properties is the Gibbs free energy [1]:

$$\mathbf{d}\,d\mathbf{G} = \mathbf{V}d\mathbf{P} - \mathbf{S}dT + \sum\_{i=1}^{N} \left.\frac{\partial \mathbf{G}}{\partial \mathbf{n}\_{i}}\right|\_{\mathbf{P}\_{i},\mathbf{T}\_{i},\mathbf{n}\_{j}} d\mathbf{x}\_{i} \tag{2}$$

where G, V, P, S, and T are, respectively, Gibbs free energy, volume, pressure, entropy, and temperature of the solution. The <sup>∂</sup> *<sup>G</sup>* ∂ *ni* <sup>|</sup> *<sup>P</sup>*,*<sup>T</sup>* ,*nj* represents the change in the Gibbs energy as a result of changes in the concentration of species *i* while the pressure, temperature, and the molar content of other species are kept constant. This is known as the chemical potential and in some textbooks is shown by *μi* . For a real solution, the chemical potential for each species is expressed by

$$
\mu\_i = \mu\_i^\circ + \mathbb{RT}lma\_i \tag{3}
$$

In which *μi <sup>o</sup>* and *ai* are the chemical potential of species *i* in its standard conditions and activity of the species *i*. The activity is defined as

$$\mathbf{a}\_{i} = \mathbf{y}\_{i}\mathbf{x}\_{i} \tag{4}$$

where *γ<sup>i</sup>* is the activity coefficient which is 1 for ideal solutions. Inserting Equation (4) into Equation (3), results in

Prediction of Solubility of Active Pharmaceutical Ingredients in Single Solvents and Their Mixtures — Solvent Screening http://dx.doi.org/10.5772/60982 3

$$
\mu\_i = \mu\_i^o + RTl\nu\nu\_i + RTl\nu x\_i \tag{5}
$$

The term *RTlnγ<sup>i</sup>* accounts for the nonideality of the solution (it is also referred to as the partial molar energy). As it can be seen from Equation (5), the terms on the right side, except *Tlnγ<sup>i</sup>* , are calculated from the pure properties. However, in order to have an accurate prediction of the chemical potential of any species in the solution, *RTlnγ<sup>i</sup>* should be known as well. In general, the activity coefficient is a function of temperature and composition and to a much less extent, the pressure. Because the activity coefficient is defined for a liquid solution, the pressure has very little effect on it. However, temperature and mole fractions of the species have significant effects on the activity coefficient of each species in a solution.

This chapter provides the reader with different and the most up-to-date thermodynamic models to estimate the activity coefficients of active pharmaceutical ingredients (API) in a solution.

#### **1.1. Thermodynamics of the solutions containing dissolved solids**

and scientists to better design, optimize, and operate the process units in the oil, gas, petro‐ chemicals, pharmaceuticals, agricultural, and other chemical-related industries. Knowledge of the thermodynamics is essential to improve process performance and product quality. For example, the use of phase behavior calculations to understand and estimate the production rate of solids from a crystallization process in a pharmaceutical industry is of paramount

A solution is called ideal if its mixture property is a linear combination of the properties of

1

solutions can be considered ideal. The interactions between the solute and solvent molecules in the solution renders a solution nonideal. The fundamental relationship which relates the

*N*

*<sup>G</sup> dG VdP SdT dx*

*o*

*i ii a x* = g

m m <sup>1</sup> , , *<sup>j</sup>*

*i i PTn*

<sup>=</sup> *n* ¶

where G, V, P, S, and T are, respectively, Gibbs free energy, volume, pressure, entropy, and

of changes in the concentration of species *i* while the pressure, temperature, and the molar content of other species are kept constant. This is known as the chemical potential and in some

<sup>=</sup> å (1)

represent the molar property of the mixture and pure component *i*,

*i*

. For a real solution, the chemical potential for each species is expressed

are the chemical potential of species *i* in its standard conditions and activity

is the activity coefficient which is 1 for ideal solutions. Inserting Equation (4) into

= -+å¶ (2)

represents the change in the Gibbs energy as a result

*ii i* = + *RTlna* (3)

(4)

denotes the mole fraction of each constituent *i* in the solution. Not all

*N solution i i i m xm* =

importance in modeling, optimization, and control of product quality.

each of its constituents at the given temperature and pressure [1]:

thermodynamic properties is the Gibbs free energy [1]:

∂ *ni* <sup>|</sup> *<sup>P</sup>*,*<sup>T</sup>* ,*nj*

where *msolution* and *mi*

2 Recent Advances in Thermo and Fluid Dynamics

temperature of the solution. The <sup>∂</sup> *<sup>G</sup>*

textbooks is shown by *μi*

*<sup>o</sup>* and *ai*

Equation (3), results in

of the species *i*. The activity is defined as

by

In which *μi*

where *γ<sup>i</sup>*

respectively, and *xi*

For a solid in equilibrium with itself in a solution, there is a thermodynamic relation [2]:

$$
\hat{f}\_{l}^{\ast} \left( T, P \right) = \hat{f}\_{l}^{\perp} \left( T, P\_{\prime} x\_{l} \right) \tag{6}
$$

where *f* ^ *i <sup>s</sup>*(*T* , *P*) denotes fugacity of component *i* in solid phase and *f* ^ *i <sup>L</sup>* (*T* , *P*, *xi* ) represents the fugacity in the liquid phase. Equation (7) below correlates the fugacity of a pure component *i* in solid and liquid state ( *f <sup>i</sup> <sup>s</sup>* and *f <sup>i</sup> <sup>L</sup>* , respectively).

$$f\_i^\* = \mathbf{x}\_i \mathbf{y}\_i f\_i^L \tag{7}$$

In order to find the ratio of the two fugacities, we need to consider all the thermodynamic processes that are involved from a solid to the liquid state. The three processes are shown in Figure 1. Path A in this figure shows the transformation from process conditions to the state where the solid starts melting. Path B shows the melting process at constant temperature and pressure. Path C indicates the change from melting to the process state. The sum of the three paths will give us the whole change from solid to liquid state of a pure component.

From the fundamental rules in thermodynamics, we have:

$$
\Delta G\_{s \to L} = RT \ln \frac{f\_i^L}{f\_i^S} \tag{8}
$$

And the Gibbs energy change is related to change of enthalpy and entropy:

$$\begin{array}{|l|l|}\hline \textbf{Solid State at} & \textbf{B} & \textbf{Liquid State at} \\ \hline \text{T}=\texttt{I\_{t\_{\text{blood}}}}\text{and P=} & \textbf{B} & \textbf{T}=\texttt{I\_{t\_{\text{blood}}}}\text{and P=} \\ \hline \text{P} & & & \textbf{P}\_{\text{boson}} \\ \hline \text{A} & & & \textbf{C} \\ \hline \text{A} & & & \textbf{C} \\ \hline \text{T}=\texttt{I\_{t\_{\text{Procast}}}}\text{and} & & & \textbf{P}\_{\text{Proages}} \\ \hline \text{P=P\_{\text{Process}}} & & & & \textbf{P}\_{\text{Proages}} \\ \hline \end{array}$$

*G H TS sL sL sL* ®® ® D =D - D (9)

**Figure 1.** Schematic diagram for finding the fugacity change from solid to liquid state of a pure substance.

From Figure 1, the whole change in enthalpy is:

$$\begin{aligned} \mathbf{C}\_{p,\text{solid}}dT + \Delta H\_{f\text{solid}} + \int\_{T\_f}^{T} \mathbf{C}\_{p,\text{liquid}}dT &= \int\_{T\_f}^{T} \Delta \mathbf{C}\_p dT + \Delta H\_{f\text{solid}} \\ \Delta H\_{s\to L} &= \Delta H\_A + \Delta H\_B + \Delta H\_C = \int\_{T} \end{aligned} \tag{10}$$

where Δ*Cp* =*Cp*,*liquid* −*Cp*,*solid* . In the same manner, the entropy change from solid to liquid state can be found from

$$
\Delta \mathbf{S}\_{s \to L} = \int\_{T\_f}^{T} \frac{\Delta \mathbf{C}\_p}{T} dT + \Delta \mathbf{S}\_{fusion} \tag{11}
$$

Also, Δ*<sup>S</sup> fusion* <sup>=</sup> <sup>Δ</sup>*<sup>H</sup> fusion <sup>T</sup> fusion* . If we substitute Equations (10) and (11) into Equation (9), then:

$$RT\ln\frac{f\_i^L}{f\_i^S} = \prod\_{T\_f}^T \Delta \mathcal{C}\_p dT + \Delta H\_{f\_{\text{fusion}}} + \int\_{T\_f}^T \frac{\Delta \mathcal{C}\_p}{T} dT + \frac{\Delta H\_{f\_{\text{fusion}}}}{T\_{f\_{\text{fusion}}}} \tag{12}$$

If we neglect the terms including change in the heat capacity (because of the large value of heat of fusion compared to the heat capacities), then we will get the following important equation:

$$
\ln \propto\_{s} = \frac{-\Delta H\_{fus}}{R} \left( \frac{1}{T\_w} - \frac{1}{T} \right) - \ln \gamma\_s \tag{13}
$$

where *xs* is the equilibrium mole fraction of the solute (dissolved solid) in a solution at temperature *T* . Equation (13) is the starting point in every solid-liquid equilibrium calculation procedure that relates the three important variables *xs*, *T* , and *γs*. However, we already know that *γs* is a function of solution temperature and the mole fraction of the species. Therefore, one can fix all the thermodynamic properties of a solution if the mole fraction of the species and the solution temperature are known. We recall that although the pressure has to be fixed as well, but due to its minimal effect on the activity coefficient and other properties of the system, we don't need to have it for calculations.

#### **1.2. Classification of the thermodynamic models for solutions**

*G H TS sL sL sL* ®® ® D =D - D (9)

**Figure 1.** Schematic diagram for finding the fugacity change from solid to liquid state of a pure substance.

*C dT H C dT C dT H*

*f f*

where Δ*Cp* =*Cp*,*liquid* −*Cp*,*solid* . In the same manner, the entropy change from solid to liquid state

*<sup>T</sup> fusion* . If we substitute Equations (10) and (11) into Equation (9), then:

D D

= D +D + + ò ò (12)

*i p fusion*

*f C H*

*f T T*

*i T T fusion*

*sL A B C*

D =D +D +D =

*H HHH* ®

*f*

*f f*

*L T T*

*RTln C dT H dT*

*S p fusion*

*T*

*T*

*p s L fusion*

*C S dT S T* ® D

*T T p solid fusion p liquid p fusion T T*

+D + = D +D

ò ò

*f*

(10)

*T*

D = +D ò (11)

ò

*T*

From Figure 1, the whole change in enthalpy is:

4 Recent Advances in Thermo and Fluid Dynamics

can be found from

Also, Δ*<sup>S</sup> fusion* <sup>=</sup> <sup>Δ</sup>*<sup>H</sup> fusion*

, ,

In order to predict the solubility of solids in pure and mixed solvents, the use of noncorrelative (predictive) thermodynamic models is of high importance. Since several years ago, there have been many thermodynamic models introduced to help scientists and engineers in the predic‐ tion of phase behaviors of various liquid-liquid and vapor-liquid systems, such as the Margules equation [3], the Wilson equation [4], the Van Laar equation [5], the nonrandom two-liquid (NRTL) equation [6], and the UNIQUAC equation [7]. As these models are applicable in the estimation of activity coefficients, they can also be utilized to predict solid-liquid equilibrium systems. In general, the models which govern the activity coefficients of the species in the solutions are grouped into two categories:


As it is apparent from the above two categories, the first one is not useful for the prediction of the phase behavior of mixtures, specifically for mixtures of more than two species. Because of the time-consumingandexpensivenatureofbinaryinteractionparameter evaluationofvarious chemical components, the first group of thermodynamic models (above) is not practical. As an example, the activity coefficient from Wilson's equation of state is found from [4]:

$$\ln \mathbf{y}\_i = -\ln \left[ \sum\_{j=1}^{n} \mathbf{x}\_j \mathbf{A}\_{ij} \right] + \mathbf{1} - \sum\_{z=1}^{n} \frac{\mathbf{x}\_z \mathbf{A}\_{kl}}{\sum\_{k=1}^{n} \mathbf{x}\_k \mathbf{A}\_{zk}} \tag{14}$$

The binary interaction parameter is:

$$\Lambda\_{\boldsymbol{\psi}} = \frac{V\_{\boldsymbol{\phi}}^{\perp}}{V\_{\boldsymbol{\iota}}^{\perp}} \exp\left[\frac{-\mathcal{\lambda}\_{\boldsymbol{\imath}\boldsymbol{\jmath}} - \mathcal{\lambda}\_{\boldsymbol{\imath}\boldsymbol{\varepsilon}}}{RT}\right] \tag{15}$$

where *i* and *j* refer to the compounds present in the solution. From Equation (14), it can be seen that for obtaining the activity coefficient of a component 1 in a pure solvent 2, we need four interaction parameters (Λ12, Λ21, Λ<sup>11</sup> ∧Λ22, which are temperature dependent. In addition, from Equation (15), it is evident that for calculating the value of the binary interaction parameters, additional experimental data, such as molar volume is needed.

From the predictive category (second group), the universal functional activity coefficient (UNIFAC) model is a well-known example. The main application of the UNIFAC model is in systems showing nonelectrolytic and nonideal behaviors. Fredenslund et al. [8] developed the UNIFAC model. The NRTL segment activity coefficient (NRTL-SAC) model was first intro‐ duced in 2004 by Chen et al. [9]. This model was proposed in order to compensate for the weakness of the UNIFAC in predicting the solubility of complex chemical molecules with functional groups that had not been studied for the UNIFAC parameters. Also, in some cases, the UNIFAC group addition rule becomes invalid [2]. One of the main advantages of the NRTL-SAC model in comparison to the other predictive methods is its ability to predict organic electrolyte systems. The UNIFAC method identifies the molecule in terms of its functional groups, while the NRTL-SAC model divides the whole surface of the molecule to four segments. These segments' values are found from their interactions with other molecules in a solution. Based on Chen et al., three purely hydrophilic, hydrophobic, and polar segments have been selected as water, hexane, and acetonitrile, respectively.

#### **2. Ideal solutions**

An ideal solution is simply defined as a mixture of chemical components with its thermody‐ namic property related to the linear sum of each pure species thermodynamic property (Equation (1)). A common example is a solution which obeys Raoult's law. This law states that the total pressure of a system is a linear combination of the component's vapors pressure at the system's temperature, provided that the total pressure is less than 5 atm. In order to derive Raoult's law, we start from Equation (5) and assume that the liquid solution is ideal:

$$
\mu\_i = \mu\_i^o + \text{RTllux}\_i \tag{16}
$$

If Equation (16) is written for component *i* in both the liquid and vapor phases, we have:

Prediction of Solubility of Active Pharmaceutical Ingredients in Single Solvents and Their Mixtures — Solvent Screening http://dx.doi.org/10.5772/60982 7

$$
\mu\_{\text{i}}^{\text{L}} = \mu\_{\text{i}}^{\circ, \text{L}} + RT\text{lux}\_{\text{i}} \tag{17}
$$

$$
\mu\_i^V = \mu\_i^{o,V} + \text{RTlm}\frac{\ddot{f}\_i}{P} \tag{18}
$$

where *f* ^ *i* and *P* are the fugacity of species *i* in the vapor mixture and the total pressure, respectively. From the thermodynamic equilibrium criteria, if the two phases of liquid and vapor are in a state of equilibrium, then:

$$
\mu\_l^L = \mu\_l^V \tag{19}
$$

After substitution of Equations (17) and (18) into Equation (19), we get:

$$
\mu\_i \mu\_i^{o,L} + \text{RTllux}\_i = \mu\_i^{o,V} + \text{RTltu} \frac{\hat{f}\_i}{P} \tag{20}
$$

If Equation (19) is written for the case of pure component *i* in liquid phase at equilibrium with its vapor phase, then:

$$
\mu\_i^{o,L} = \mu\_i^{o,V} + \text{RTlm}\frac{f\_i}{P} \tag{21}
$$

By comparing Equations (20) and (21), it yields:

$$\text{RTllux}\_{i} = \text{RTl} \ln \frac{\hat{f}\_{i}}{f\_{i}} \tag{22}$$

Equation (22) is simplified to:

The binary interaction parameter is:

6 Recent Advances in Thermo and Fluid Dynamics

*L*

*ij L i*

additional experimental data, such as molar volume is needed.

have been selected as water, hexane, and acetonitrile, respectively.

**2. Ideal solutions**

*V*

*j ij ii*

where *i* and *j* refer to the compounds present in the solution. From Equation (14), it can be seen that for obtaining the activity coefficient of a component 1 in a pure solvent 2, we need four interaction parameters (Λ12, Λ21, Λ<sup>11</sup> ∧Λ22, which are temperature dependent. In addition, from Equation (15), it is evident that for calculating the value of the binary interaction parameters,

From the predictive category (second group), the universal functional activity coefficient (UNIFAC) model is a well-known example. The main application of the UNIFAC model is in systems showing nonelectrolytic and nonideal behaviors. Fredenslund et al. [8] developed the UNIFAC model. The NRTL segment activity coefficient (NRTL-SAC) model was first intro‐ duced in 2004 by Chen et al. [9]. This model was proposed in order to compensate for the weakness of the UNIFAC in predicting the solubility of complex chemical molecules with functional groups that had not been studied for the UNIFAC parameters. Also, in some cases, the UNIFAC group addition rule becomes invalid [2]. One of the main advantages of the NRTL-SAC model in comparison to the other predictive methods is its ability to predict organic electrolyte systems. The UNIFAC method identifies the molecule in terms of its functional groups, while the NRTL-SAC model divides the whole surface of the molecule to four segments. These segments' values are found from their interactions with other molecules in a solution. Based on Chen et al., three purely hydrophilic, hydrophobic, and polar segments

An ideal solution is simply defined as a mixture of chemical components with its thermody‐ namic property related to the linear sum of each pure species thermodynamic property (Equation (1)). A common example is a solution which obeys Raoult's law. This law states that the total pressure of a system is a linear combination of the component's vapors pressure at the system's temperature, provided that the total pressure is less than 5 atm. In order to derive

*ii i* = + *RTlnx* (16)

Raoult's law, we start from Equation (5) and assume that the liquid solution is ideal:

*o*

If Equation (16) is written for component *i* in both the liquid and vapor phases, we have:

m m é ù - l l

ê ú ë û

(15)

*exp <sup>V</sup> RT*

L = ê ú

$$
\hat{f}\_i = f\_i \mathbf{x}\_i \tag{23}
$$

Now, if the liquid and vapor mixtures are assumed to be ideal, then the fugacity values can be substituted by their corresponding pressure and thus:

$$P\_i = P\_i^{sat} \propto\_i \tag{24}$$

*Pi* is called the partial pressure of species *i* in the vapor mixture and *Pi sat* denotes the pure vapor pressure of the component *i* at the solution temperature.

#### **2.1. Ideal solution mixtures: VLE phase behavior**

In order to calculate the properties of a system of components that obey Raoult's law, Equation (24) is written for each of the species present in the solution,

$$\sum\_{l=1}^{N} \mathbf{y}\_{l} P\_{\text{tot}} = \sum\_{l=1}^{N} \mathbf{P}\_{l}^{\text{sat}} \mathbf{x}\_{l} \tag{25}$$

and by simplification:

pressures, then:

**Bubble point**

**Dew point**

cases are possible:

$$P\_{\text{tot}} = \sum\_{i=1}^{N} P\_i^{sat} \mathbf{x}\_i \tag{26}$$

by which the solution equilibrium temperature can be found. As we know, the vapor pressure of each species is temperature-dependent and by having the mole fraction of the components in the liquid phase, one can find the temperature which corresponds to the total pressure of the system. This problem can be solved quickly using a spreadsheet and by changing the solution temperature until Equation (26) is satisfied. Sometimes, the temperature of the solution is known and the total pressure needs to be calculated, which is a straightforward calculation from Equation (26). This way, there is no need to use the nonlinear solvers to find the temperature.

There are two examples of the binary systems which exhibit negative deviation from Raoult's law.

**Figure 1.** T‐x‐y diagram for the binary systems of hexane‐benzene (left) and ethyl acetate‐benzene (right) at a pressure of 101.33 kPa. **Figure 2.** T-x-y diagram for the binary systems of hexane-benzene (left) and ethyl acetate-benzene (right) at a pressure of 101.33 kPa.

In order to better demonstrate whether a system follows Raoult's law, a diagram of the phase equilibrium called T‐x‐y should be plotted. This plot (Figure 1) shows the equilibrium temperatures at which either a liquid solution will start bubbling (bubble curve) or a vapor mixture starts condensing (dew curve). The two systems with their experimental data and the calculation curve of the ideal solution is shown in Figure 1. In Figure 1, the system of hexane‐benzene at the pressure of 101.33 kPa

If one is interested in finding the bubble and dew curves for an ideal solution at low to moderate

From Equation (26), by knowing the mole fraction of each component in the liquid phase, two

**1. Pressure known.** If the total pressure of the system is given, then Equation (26) should be

��� � ���� ����

**2. Temperature known.** For this case, the problem is straightforward. The vapor pressure of each species can be found readily from Antoine's equation and then inserted into Equation

�

���

��������

�27�

solved for the temperature (if the vapor pressure follows Antoine's law):

�

��

(27) to find the total pressure.

��� ���

[10] and the system of ethylacetate‐benzene [11] show negative deviations from Raoult's law.

In order to better demonstrate whether a system follows Raoult's law, a diagram of the phase equilibrium called T-x-y should be plotted. This plot (Figure 2) shows the equilibrium temperatures at which either a liquid solution will start bubbling (bubble curve) or a vapor mixture starts condensing (dew curve). The two systems with their experimental data and the calculation curve of the ideal solution is shown in Figure 2. In Figure 2, the system of hexanebenzene at the pressure of 101.33 kPa [10] and the system of ethylacetate-benzene [11] show negative deviations from Raoult's law.

If one is interested in finding the bubble and dew curves for an ideal solution at low to moderate pressures, then:

#### **• Bubble point**

*Pi*

and by simplification:

the temperature.

pressures, then:

of 101.33 kPa.

**Bubble point**

**Dew point**

cases are possible:

law.

is called the partial pressure of species *i* in the vapor mixture and *Pi*

In order to calculate the properties of a system of components that obey Raoult's law, Equation

*sat i tot i i*

å å<sup>=</sup> (25)

<sup>=</sup> å (26)

1 1

*yP P x* = =

1

by which the solution equilibrium temperature can be found. As we know, the vapor pressure of each species is temperature-dependent and by having the mole fraction of the components in the liquid phase, one can find the temperature which corresponds to the total pressure of the system. This problem can be solved quickly using a spreadsheet and by changing the solution temperature until Equation (26) is satisfied. Sometimes, the temperature of the solution is known and the total pressure needs to be calculated, which is a straightforward calculation from Equation (26). This way, there is no need to use the nonlinear solvers to find

There are two examples of the binary systems which exhibit negative deviation from Raoult's

**Figure 1.** T‐x‐y diagram for the binary systems of hexane‐benzene (left) and ethyl acetate‐benzene (right) at a pressure of 101.33 kPa. In order to better demonstrate whether a system follows Raoult's law, a diagram of the phase equilibrium called T‐x‐y should be plotted. This plot (Figure 1) shows the equilibrium temperatures at which either a liquid solution will start bubbling (bubble curve) or a vapor mixture starts condensing (dew curve). The two systems with their experimental data and the calculation curve of the ideal solution is shown in Figure 1. In Figure 1, the system of hexane‐benzene at the pressure of 101.33 kPa

**Figure 2.** T-x-y diagram for the binary systems of hexane-benzene (left) and ethyl acetate-benzene (right) at a pressure

[10] and the system of ethylacetate‐benzene [11] show negative deviations from Raoult's law.

If one is interested in finding the bubble and dew curves for an ideal solution at low to moderate

From Equation (26), by knowing the mole fraction of each component in the liquid phase, two

**1. Pressure known.** If the total pressure of the system is given, then Equation (26) should be

��� � ���� ����

**2. Temperature known.** For this case, the problem is straightforward. The vapor pressure of each species can be found readily from Antoine's equation and then inserted into Equation

�

���

��������

�27�

solved for the temperature (if the vapor pressure follows Antoine's law):

�

��

(27) to find the total pressure.

��� ���

*N sat tot i i i P Px* =

*i i*

*N N*

pressure of the component *i* at the solution temperature.

(24) is written for each of the species present in the solution,

**2.1. Ideal solution mixtures: VLE phase behavior**

8 Recent Advances in Thermo and Fluid Dynamics

*sat* denotes the pure vapor

From Equation (26), by knowing the mole fraction of each component in the liquid phase, two cases are possible:

**1. Pressure known.** If the total pressure of the system is given, then Equation (26) should be solved for the temperature (if the vapor pressure follows Antoine's law):

$$P\_i^{sat} = A - \frac{B}{T + C} \, ^\prime P\_{\text{tot}} = \sum\_{l=1}^N P\_l^{sat} \left( ^t \Gamma \right) x\_l \tag{27}$$


$$y\_i = \frac{P\_i^{sat} \mathbf{x}\_i}{P\_{tot}} \tag{28}$$

by taking the sum for all the components in the liquid phase and knowing that the sum of the mole fractions in a phase is 1, we get:

$$\mathbf{1} = \sum\_{i=1}^{N} \mathbf{x}\_i = \frac{y\_i P\_{\text{tot}}}{P\_i^{\text{sat}} \left(T\right)} \tag{29}$$

If the total pressure is known, then Equation (29) should be solved for the temperature which determines the vapor pressure of each of the species in the mixture.

**2. Temperature known.** In this case, Equation (29) is isolated for the *Ptot* to find the total pressure of the system:

$$i = \frac{i}{\sum N \frac{P\_i^{sat} \left(T\right)}{\mathcal{Y}\_i}} \tag{30}$$

which has a straightforward solution.

#### **2.2. Ideal solution mixtures: SLE phase behavior**

For a solid-liquid equilibrium behavior, Equation (13) for the ideal solution (in which *γ<sup>s</sup>* =1) is simplified to:

$$
\hbar \text{m} \mathbf{x}\_s = \frac{-\Delta H\_{fus}}{\mathcal{R}} \left( \frac{1}{T\_m} - \frac{1}{T} \right) \tag{31}
$$

From Equation (31), one can find the mole fraction (or what is called *solubility* in the case of SLE) by inserting the appropriate values for the enthalpy of fusion, Δ*H fus*, melting tempera‐ ture, *Tm*, and the solution temperature, *T* . It can be found from this equation that for any solvent, the solubility of a solid will be the same regardless of the nature of the solvent. For very few cases, this assumption might be correct; however, for most of the practical and common applications, the above formula does not work well. Therefore, a term accounting for the nonideality of the system should be added to Equation (31).

#### **3. Nonideal solutions**

Almost all real VLE or SLE systems show nonideality. Equations that predict the behavior of different phase equilibria are divided into two:


We focus on the thermodynamic models that deal with the liquid mixtures in this chapter. From the two categories of activity coefficient models, the correlative one is not very useful for solubility prediction and solvent screening purposes. The main reason for this is the lack of experimental data for the binary interaction parameters of the solute-solvent, soluteantisolvent, and solvent-antisolvent systems. As an example, the activity coefficient from Wilson's equation of state is found from Equations (14) and (15). It can be seen that for obtaining the activity coefficient of a component 1 in a pure solvent 2, we need four interaction param‐ eters (Λ12, Λ21, Λ<sup>11</sup> ∧Λ22, which are temperature dependent. It is evident that for calculating the value of the binary interaction parameters, additional experimental data, such as molar volume is needed. Other models which belong to the first category have the same limitations as Wilson's method. The Wilson model was used in the prediction of various hydrocarbons in water in pure form and mixed with other solvents by Matsuda et al. [11]. In order to estimate the pure properties of the species, the Tassios method [12] with DECHEMA VLE handbook [13] were used. Matsuda et al. also took some assumptions in the estimation of binary interactions (because of the lack of data) that resulted in some deviations from the experimental data.

**2. Temperature known.** In this case, Equation (29) is isolated for the *Ptot* to find the total

1

*i i*

For a solid-liquid equilibrium behavior, Equation (13) for the ideal solution (in which *γ<sup>s</sup>* =1) is

*fus* 1 1

*H*

*m*

è ø

*R TT* -D æ ö <sup>=</sup> ç ÷ -

From Equation (31), one can find the mole fraction (or what is called *solubility* in the case of SLE) by inserting the appropriate values for the enthalpy of fusion, Δ*H fus*, melting tempera‐ ture, *Tm*, and the solution temperature, *T* . It can be found from this equation that for any solvent, the solubility of a solid will be the same regardless of the nature of the solvent. For very few cases, this assumption might be correct; however, for most of the practical and common applications, the above formula does not work well. Therefore, a term accounting for

Almost all real VLE or SLE systems show nonideality. Equations that predict the behavior of

**•** Equations of state (EOS) that are useful to predict the nonideal behavior of the vapor phase

**•** Equations which are applied to the liquid phase to predict the nonideal behavior of the liquid

We focus on the thermodynamic models that deal with the liquid mixtures in this chapter. From the two categories of activity coefficient models, the correlative one is not very useful for solubility prediction and solvent screening purposes. The main reason for this is the lack of experimental data for the binary interaction parameters of the solute-solvent, soluteantisolvent, and solvent-antisolvent systems. As an example, the activity coefficient from

*P T N*

*y*

*tot sat*

å

*P*

*s*

*lnx*

the nonideality of the system should be added to Equation (31).

=

*i*

=

( )

(30)

(31)

pressure of the system:

10 Recent Advances in Thermo and Fluid Dynamics

which has a straightforward solution.

simplified to:

**3. Nonideal solutions**

solutions

different phase equilibria are divided into two:

**2.2. Ideal solution mixtures: SLE phase behavior**

From the predictive category, we bring some examples of the application of the UNIFAC model. In one study, this model has been used to predict the solubility of naphtalene, anthra‐ cene, and phenanthrene in various solvents and their mixtures [8]. They showed the applica‐ bility of the UNIFAC model in prediction of the phase behavior of solutes in solvents. There have been efforts to make the UNIFAC model more robust and powerful in the prediction of phase behaviors [14]. In one study, the solubility of buspirone-hydrochloride in isopropyl alcohol was measured and evaluated by the modified UNIFAC model [15]. It was concluded that for highly soluble pharmaceutics, the modified form of the UNIFAC model was not suitable. In another study, the solubility of some chemical species in water and some organic solvents was predicted by the UNIFAC model [16]. For some unknown functional groups, they used other known groups which had chemical structures that were similar to unknown ones. In conclusion, it was stated that the UNIFAC model is not a proper model for use in crystal‐ lization and related processes. The UNIFAC model also has been utilized to predict the solubility of some aromatic components as well as long-chain hydrocarbons [17]. The results showed that the predictions for the linear hydrocarbons are not as good as the ones for the aromatics.

Chen et al. developed the NRTL-SAC model in 2004 [9]. One of the main reasons for developing this model was to enhance the predicting capability of the solubility of complex chemical molecules with functional groups that were not included in the UNIFAC model. Also, in some cases, the UNIFAC group addition rule becomes invalid [2]. One of the main advantages of the NRTL-SAC model in comparison to the other predictive methods is its ability to predict organic electrolyte systems [18]. The UNIFAC method identifies the molecule in terms of its functional groups, while the NRTL-SAC model divides the whole surface of the molecule to four segments. The hydrophobic segment (X) denotes parts of the molecule surface which don't participate in forming a hydrogen bond, such as hexane. The polar segments (Y- and Y+) do not belong to either hydrophobic or hydrophilic segment. The polar attractive segment (Y-) shows attractive interaction with hydrophilic segment, while the polar repulsive segment (Y +) has repulsive characteristic with hydrophilic segment. Hydrophilic segment (Z) contributes to the part of the molecule which tends to form a hydrogen bond, such as water.

Based on the interaction forces between the molecule surfaces, the four segments of the NRLT-SAC model have been identified. The three reference solvents mentioned before were used to determine the rest of chemical components' segment numbers [9]. As an example, the VLE and LLE data containing the component of interest in binary mixtures with hexane, water, and acetonitrile were collected and then used in order to find the segment numbers of the compo‐ nent. Sheikholeslamzadeh et al. have investigated various industrial case studies to show the applicability of the NRTL-SAC in processes consisting of solvent mixtures [19]. The segment numbers for each of the nonreference solvents are found from the nonlinear optimization methods which minimize the deviation from the model output and the experimental data for the whole range of data (VLE and LLE). Most of the solvents have some segment numbers that are high in value compared to other segment numbers. The reason is that it is less probable that a chemical component collects all four segments at the same level of value. Currently, more than a hundred solvent segment numbers have been evaluated and optimized which can be found in some commercial simulation packages. The two important predictive equations of UNIFAC and NRTL-SAC, as representatives of the activity coefficient models, are presented here.

#### **3.1. UNIFAC model**

In general, the activity coefficient models of the predictive category split the activity coeffi‐ cients into two segments:


With the above definition, the total activity of a component in the solution is the sum of the two parts:

$$
\ln \mathbf{\gamma}\_i = \ln \mathbf{\gamma}\_i^{\mathbb{C}} + \ln \mathbf{\gamma}\_i^{\mathbb{R}} \tag{32}
$$

In which *γ<sup>i</sup>* is the activity coefficient of component *i* in the solution, *γ<sup>i</sup> <sup>C</sup>* is the combinatorial part and *γ<sup>i</sup> <sup>R</sup>* is the residual part. Up to this point, all of the group contribution and activity coefficient methods (i.e. NRTL-SAC) have been the same, but the methods in which the activities have been calculated are different. In the UNIFAC model, the combinatorial part for component *i* is found from the following equation [8]:

$$\ln \mathbf{y}\_i^c = \ln \left[ \frac{\boldsymbol{\phi}\_i}{\boldsymbol{\chi}\_i} \right] + \frac{\boldsymbol{z}}{2} q\_i \ln \left[ \frac{\boldsymbol{\theta}\_i}{\boldsymbol{\phi}\_i} \right] + L\_i - \frac{\boldsymbol{\phi}\_i}{\boldsymbol{\chi}\_i} \sum\_{j=1}^n \boldsymbol{x}\_j L\_j \tag{33}$$

where

$$L\_i = \frac{\mathbb{Z}}{\mathbf{2}} (r\_i - q\_i) - \left(r\_i - \mathbf{1}\right) \tag{34}$$

*z* is the coordination number and is taken to be 10. In Equation (33), *ϕ<sup>i</sup>* is the segment fraction and *θ<sup>i</sup>* is the area fraction of component *i* and is related to the mole fraction of species *i* in the mixture:

determine the rest of chemical components' segment numbers [9]. As an example, the VLE and LLE data containing the component of interest in binary mixtures with hexane, water, and acetonitrile were collected and then used in order to find the segment numbers of the compo‐ nent. Sheikholeslamzadeh et al. have investigated various industrial case studies to show the applicability of the NRTL-SAC in processes consisting of solvent mixtures [19]. The segment numbers for each of the nonreference solvents are found from the nonlinear optimization methods which minimize the deviation from the model output and the experimental data for the whole range of data (VLE and LLE). Most of the solvents have some segment numbers that are high in value compared to other segment numbers. The reason is that it is less probable that a chemical component collects all four segments at the same level of value. Currently, more than a hundred solvent segment numbers have been evaluated and optimized which can be found in some commercial simulation packages. The two important predictive equations of UNIFAC and NRTL-SAC, as representatives of the activity coefficient models, are presented

In general, the activity coefficient models of the predictive category split the activity coeffi‐

**•** A part that includes the contribution of the chemical structure and the size of a compound

**•** The second part that includes the contribution of the functional sizes and binary interaction

With the above definition, the total activity of a component in the solution is the sum of the

*C R*

coefficient methods (i.e. NRTL-SAC) have been the same, but the methods in which the activities have been calculated are different. In the UNIFAC model, the combinatorial part for

<sup>1</sup> 2

qf

=

 g

*<sup>R</sup>* is the residual part. Up to this point, all of the group contribution and activity

*n*

ëû ëû <sup>å</sup> (33)

*<sup>z</sup> L rq r* = - -- (34)

= + (32)

*<sup>C</sup>* is the combinatorial

*ii i ln ln ln* gg

is the activity coefficient of component *i* in the solution, *γ<sup>i</sup>*

*c i ii i i i j j i ii j <sup>z</sup> ln ln q ln L x L x x*

éù éù = + +- êú êú

( ) ( <sup>1</sup>) <sup>2</sup> *i ii i*

f

f

between pairs of the functional groups (residual part)

component *i* is found from the following equation [8]:

g

here.

two parts:

In which *γ<sup>i</sup>*

part and *γ<sup>i</sup>*

where

**3.1. UNIFAC model**

cients into two segments:

12 Recent Advances in Thermo and Fluid Dynamics

(combinatorial part)

$$\phi\_i = \frac{r\_i \mathbf{x}\_i}{\sum\_{j=1}^{n} r\_j \mathbf{x}\_j} \tag{35}$$

$$\theta\_i = \frac{q\_i \mathbf{x}\_i}{\sum\_{j=1}^n q\_j \mathbf{x}\_j} \tag{36}$$

*qi* and *ri* are the pure component surface area and volumes (van der Waals), respectively. These parameters are not temperature-dependent and are only functions of the chemical structure of a functional group. In the UNIFAC model, for every functional group, there is a unique value for surface area and volume that can be found in common texts and handbooks [1]. The first step in modeling the UNIFAC for a specific binary or ternary system is to break down the chemical structure of a molecule into the basic functional groups. As it is suggested in thermodynamic textbooks [19], the optimum way of breaking it down is the one which results in the minimum number of subgroups, and with each subgroup having the maximum replicates. The *qi* and *ri* can be found from Equations (37) and (38):

$$r\_i = \sum\_k v\_k \, ^i \mathbb{R}\_k \tag{37}$$

$$\mathbf{q}\_{i} = \sum\_{k} \mathbf{v}\_{k}^{\top} \mathbf{Q}\_{k} \tag{38}$$

*vk <sup>i</sup>* is the number of subgroup *k* in component *i*. The residual part of the UNIFAC is found from the equation below:

$$\ln \mathbf{\hat{v}}\_i^R = \sum\_{k=1}^{n\_k} \mathbf{v}\_k^i \left[ \ln \boldsymbol{\Gamma}\_k - \ln \boldsymbol{\Gamma}\_k^{-i} \right] \tag{39}$$

In Equation (39), Γ*<sup>k</sup>* is the residual activity coefficient of subgroup *k* in the mixture and Γ*<sup>k</sup> <sup>i</sup>* is that value in a pure solution of the component *i*. This term is added so when the mole fraction approaches unity, the term *lnγ<sup>i</sup> <sup>R</sup>* tends to zero (*γ<sup>i</sup> <sup>R</sup>* →1). The residual activity coefficient of subgroup *k* in a solution is given by

$$
\ln \Gamma\_k = \mathbb{Q}\_k \left[ 1 - \ln \sum\_m \theta\_m \mathbb{W}\_{mk} - \sum\_m \left[ \frac{\theta\_m \mathbb{W}\_{mk}}{\sum\_n \theta\_n \mathbb{W}\_{nm}} \right] \right] \tag{40}
$$

Equation (40) is also applicable to the case of Γ*<sup>k</sup> <sup>i</sup>*, in which the parameters of the right-hand side of the equation are written based on the pure component *i*. *θ<sup>m</sup>* is the area fraction of the functional group *m* in the mixture:

$$\theta\_m = \frac{\mathbb{Q}\_w X\_w}{\sum\_n \mathbb{Q}\_n X\_n} \tag{41}$$

*Xm* is the mole fraction of subgroup *m* in the mixture. *ψnm* is the group interaction parameter between groups *n* and *m* and is dependent on the temperature:

$$\left[\left.\nu\right\vert\_{nm} = \exp\left[\frac{-\mu\_{nm} - \mu\_{nm}}{RT}\right] = \exp\left[\frac{-a\_{nm}}{T}\right] \tag{42}$$

The group interaction parameter *anm* is found from the large sets of VLE and LLE data in the literature, which are tabulated for many subgroups. It is worth noting that *anm* ≠*amn*. There are some modifications to the original UNIFAC equation in order to make the model robust for some complex systems. In the UNIFAC-DM method, the modification is made on the combi‐ natorial part:

$$\ln \mathbf{y}\_i^c = \ln \left[ \frac{\boldsymbol{\phi}\_i^\cdot}{\boldsymbol{\chi}\_i} \right] + \frac{\boldsymbol{z}}{2} q\_i \ln \left[ \frac{\boldsymbol{\theta}\_i}{\boldsymbol{\phi}\_i} \right] + L\_i - \frac{\boldsymbol{\phi}\_i}{\boldsymbol{\chi}\_i} \sum\_{j=1}^n \boldsymbol{\chi}\_j L\_j \tag{43}$$

In which the term *<sup>ϕ</sup><sup>i</sup>* ' *xi* is defined as

$$\frac{\phi\_i^{\cdot}}{\alpha\_i} = \frac{r^{3/4}}{\sum\_j \mathbf{x}\_j r\_j^{3/4}} \tag{44}$$

The algorithm for finding the solute solubility in the mixture of ternary solution (solvent/ cosolvent/solute) is shown in Figure 3. The algorithm starts with known values, such as the physical properties of the solute. After making an initial guess for the solubility, the program obtains the activity coefficients and the new solubility is found and is compared with the old value and the calculations are repeated to converge to a unique value for solubility. This procedure is done for all the experimental data points.

1 *m mk*

é ù é ù ê ú ê ú G= - - ê ú

q y

*m m <sup>n</sup> n nm*

ë û ê ú ë û

side of the equation are written based on the pure component *i*. *θ<sup>m</sup>* is the area fraction of the

*m m*

*Q X Q X*

*n n n*

*Xm* is the mole fraction of subgroup *m* in the mixture. *ψnm* is the group interaction parameter

*nm mm nm*

<sup>1</sup> 2

3/4

The algorithm for finding the solute solubility in the mixture of ternary solution (solvent/ cosolvent/solute) is shown in Figure 3. The algorithm starts with known values, such as the physical properties of the solute. After making an initial guess for the solubility, the program obtains the activity coefficients and the new solubility is found and is compared with the old

qf

=

*n*

= + +- ê ú ê ú ê ú ë û ë û <sup>å</sup> (43)

å (44)

*uu a exp exp RT <sup>T</sup>*

The group interaction parameter *anm* is found from the large sets of VLE and LLE data in the literature, which are tabulated for many subgroups. It is worth noting that *anm* ≠*amn*. There are some modifications to the original UNIFAC equation in order to make the model robust for some complex systems. In the UNIFAC-DM method, the modification is made on the combi‐

é ù éù -- - <sup>=</sup> ê ú êú <sup>=</sup>

q y

> q y

å å <sup>å</sup> (40)

å (41)

ë û ëû (42)

*<sup>i</sup>*, in which the parameters of the right-hand

*k k m mk*

*m*

q=

between groups *n* and *m* and is dependent on the temperature:

'

f

*c i ii i i i j j i ii j <sup>z</sup> ln ln q ln L x L x x*

é ù é ù

'3/4

*<sup>i</sup> j j <sup>j</sup>*

*r x x r*

*i*

f= f

*nm*

y

g

*xi* is defined as

'

*ln Q ln*

Equation (40) is also applicable to the case of Γ*<sup>k</sup>*

functional group *m* in the mixture:

14 Recent Advances in Thermo and Fluid Dynamics

natorial part:

In which the term *<sup>ϕ</sup><sup>i</sup>*

**Figure 3.** The algorithm of converging to the solubility of a ternary system using the UNIFAC model.

#### **3.2. Nonrandom two-liquid segment activity coefficient (NRTL-SAC)**

According to Chen et al. [9], the NRTL-SAC model is based on the derivation of the original NRTL model for polymers. From Equation (32), the activity coefficient is made up of two terms, combinatorial and residual. Like the UNIFAC model, the activity coefficients must be gener‐ ated in order to obtain solubility. In the NRTL-SAC model, the combinatorial part is calculated by Equation (45):

$$\ln \mathbf{y}\_i^{\mathbb{C}} = \ln \frac{\bigotimes\_{\mathbf{x}\_i}}{\mathbf{x}\_i} + \mathbf{1} - r\_i \sum\_{j} \frac{\bigotimes\_{\mathbf{x}\_j}}{\mathbf{x}\_j} \tag{45}$$

With the definitions:

$$r\_i = \sum\_j r\_{j,i} \tag{46}$$

$$\mathcal{D}\_{i} = \frac{r\_{i}\mathbf{x}\_{i}}{\sum r\_{i}\mathbf{x}\_{i}} \tag{47}$$

where xi is the mole fraction of component *i*, rm,i is the number of segment *m*, ri is the total segment number in component *i*, and ∅*i* is the segment mole fraction in the mixture.

The residual term is defined as

$$\ln \text{r} \mathbf{y}\_i^R = \ln \mathbf{y}\_i^{lc} = \sum\_m \mathbf{r}\_{m,i} \left( \ln \Gamma\_m{}^{\;k\;l} - \ln \Gamma\_m{}^{\;k,l} \right) \tag{48}$$

In Equation (48), there are two terms, *ln*Γ*mlc* and *ln*Γ*mlc*,*<sup>i</sup>*, which are the activity coefficients of segment *m* in solution and component *i*, respectively.

The two mentioned terms are found using Equations (49) and (50):

$$\ln \boldsymbol{\Gamma}\_{\boldsymbol{m}}{}^{k} = \frac{\sum\_{j} \mathbf{x}\_{j} \mathbf{G}\_{j,\boldsymbol{m}} \boldsymbol{\pi}\_{j,\boldsymbol{m}}}{\sum\_{k} \mathbf{x}\_{k} \mathbf{G}\_{k,\boldsymbol{m}}} + \sum\_{\boldsymbol{m}'} \frac{\mathbf{x}\_{\boldsymbol{m}'} \mathbf{G}\_{\boldsymbol{m},\boldsymbol{m}'}}{\sum\_{k} \mathbf{x}\_{k} \mathbf{G}\_{k,\boldsymbol{m}'}} \left( \boldsymbol{\pi}\_{\boldsymbol{m},\boldsymbol{m}'} - \frac{\sum\_{j} \mathbf{x}\_{j} \mathbf{G}\_{j,\boldsymbol{m}'} \boldsymbol{\pi}\_{j,\boldsymbol{m}'}}{\sum\_{k} \mathbf{x}\_{k} \mathbf{G}\_{k,\boldsymbol{m}'}} \right) \tag{49}$$

$$\ln \Pi\_{\;m} \stackrel{\text{l.c.}}{=} \frac{\sum\_{j} \mathbf{x}\_{j,l} \mathbf{G}\_{j,m} \mathbf{r}\_{j,m}}{\sum\_{k} \mathbf{x}\_{k,l} \mathbf{G}\_{k,m}} + \sum\_{m'} \frac{\mathbf{x}\_{m',l} \mathbf{G}\_{m,m'}}{\sum\_{k} \mathbf{x}\_{k,l} \mathbf{G}\_{k,m'}} \left( \mathbf{r}\_{m,m'} - \frac{\sum\_{j} \mathbf{x}\_{j,l} \mathbf{G}\_{j,m'} \mathbf{r}\_{j,m'}}{\sum\_{k} \mathbf{x}\_{k,l} \mathbf{G}\_{k,m'}} \right) \tag{50}$$

In the two above equations, *l* is referred to as the component and *j, k, m,* and *m*′ are referred to the segments in each component. xj,l is the segment-based mole fraction of segment species *j* in component *l* only. The mole fractions of segments in the whole solution and in components are defined as below:

$$\mathbf{x}\_{\rangle} = \frac{\sum\_{l} \mathbf{x}\_{l} r\_{j,l}}{\sum\_{z} \sum\_{l} \mathbf{x}\_{z} r\_{j,z}} \tag{51}$$

Prediction of Solubility of Active Pharmaceutical Ingredients in Single Solvents and Their Mixtures — Solvent Screening http://dx.doi.org/10.5772/60982 17

$$\mathbf{x}\_{j,l} = \frac{r\_{j,l}}{\sum\_{i} r\_{j,l}} \tag{52}$$

*Gi*, *<sup>j</sup>* and *τi*, *<sup>j</sup>* are the local binary values which can be related to each other based on NRTL nonrandom parameter *αij* , and are shown by their values in Table 1. *Gi*, *<sup>j</sup>* and *τi*, *<sup>j</sup>* have the following relation:

1 *<sup>C</sup> <sup>j</sup> <sup>i</sup> i i*

> *i ji*, *j*

segment number in component *i*, and ∅*i* is the segment mole fraction in the mixture.

, *R lc lc lc i i i mi m m m ln ln r ln ln*

In Equation (48), there are two terms, *ln*Γ*mlc* and *ln*Γ*mlc*,*<sup>i</sup>*, which are the activity coefficients of

, , , , , , ,, ,

t

 ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢

,

 ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢

å å <sup>å</sup> åå å (50)

t

æ ö

è ø å å <sup>å</sup> åå å (49)

æ ö

è ø

*kk k k km <sup>m</sup> k km k km*

, , , , , ,

*kk k kl km <sup>m</sup> kl km kl km*

*xG xG x G*

In the two above equations, *l* is referred to as the component and *j, k, m,* and *m*′ are referred to the segments in each component. xj,l is the segment-based mole fraction of segment species *j* in component *l* only. The mole fractions of segments in the whole solution and in components

> , ,

*l jl l*

*x r*

*x r* <sup>=</sup> <sup>å</sup>

*z jz z i*

*x G x G x G*

*j j j jm jm j jm jm lc m mm*

,,, ,, , , , ,

G= + ç ÷ -

*j j jl jm jm jl jm jm lc l m l mm*

*x G x G x G*

G= + ç ÷ -

*x G x G x G*

*i*

Æ =

*i i*

*r x r x*

*j j j*

is the mole fraction of component *i*, rm,i is the number of segment *m*, ri is the total

( ) ,

= = G -G å (48)

 t

> t

å å (51)

*ln ln r*

g

With the definitions:

16 Recent Advances in Thermo and Fluid Dynamics

The residual term is defined as

*ln*

*ln*

are defined as below:

g

segment *m* in solution and component *i*, respectively.

 g

The two mentioned terms are found using Equations (49) and (50):

t

t

*m m m*

*m m m*

*j*

*x*

where xi

*i j j*

= +- å (45)

*r r* <sup>=</sup> å (46)

å (47)

*x x*

Æ Æ

$$\mathbf{G}\_{i,j} = \mathbf{e}^{-\alpha\_{l,j}\mathbf{r}\_{l,j}} \tag{53}$$

Therefore, from fixed values of *τi*, *<sup>j</sup>* and *αi*, *<sup>j</sup>* one can find *Gi*, *<sup>j</sup>* . The segment numbers for the common solvents can be found from the literature [20]. After putting the values of segments for solvents and initial guess values for the solute segments, the written code for NRTL-SAC starts solving for the mole fractions at saturation for all of the species in the solu‐ tion (see Figure 3).

**Figure 4.** Algorithm flowchart for parameter estimation using NRTL-SAC model.

It is worth noting that the main difference in Figures 3 and 4 is the use of parameter estimation method for the calculation of the NRTL-SAC parameters, while for the UNIFAC model, the calculation is straightforward. Once the parameters (here, the segment numbers) are found, then they could be used for validation against other experimental data.

### **4. Application of solution thermodynamics in industry**

One of the main applications of the thermodynamic models is in the chemical industries which use solvent (or their mixtures) [19-22]. Two cases of the vapor-liquid equilibrium of common industrial solvent systems are discussed here.

#### **4.1. VLE study of two binary azeotropic systems**

There are some solvents within the chemical and pharmaceutical industries which are of importance to study, such as ethanol, isopropyl alcohol, and water in addition to some aromatic components, such as benzene and toluene. These solvents are sometimes used as additives to other valuable chemicals to maintain some performances or enhance their physical or chemical properties. The systems that form azeotropes make the distillation process calculations complex. As a result, having an accurate knowledge about the phase behaviors of such systems are important (such as toluene-ethanol or toluene-isopropyl alcohol). The four case studies of the VLE data for the mentioned azeotropic binary systems were used in a study by Chen et al. [24]. The operating pressures were at four distinct levels for the calculations. They have used three equations of state (from the correlative group) to fit the experimental data with the model outputs and found the necessary binary interaction parameters. Based on their work, the estimated parameters were found; however, they were not used to further validate the models for other operating conditions. In a study by Sheikholeslamzadeh et al., they used the NRTL-SAC and UNIFAC models to predict these azeotropic systems [20]. Table 1 (from their work) shows average relative deviation for the compositions and temperature for the two systems. It is seen from the table that the deviations from the experimental data in the vapor phase mole fractions are almost one-third relative to the NRTL-SAC model. On the other hand, the relative deviation for the saturation temperature is higher using the NRTL-SAC model. Another fact from this finding is that the deviations for both binary systems when the NRTL-SAC model is used are nearly the same. This is not the case when utilizing the UNIFAC model to predict the systems' behaviors. The final conclusion from Table 1 would be the better predictive capability of the UNIFAC model for light alcohols than the heavier ones.


Prediction of Solubility of Active Pharmaceutical Ingredients in Single Solvents and Their Mixtures — Solvent Screening http://dx.doi.org/10.5772/60982 19


**Table 1.** The results of the prediction using the NRTL-SAC and UNIFAC models with the experimental values of Chen et al. [20,24] 201.3 0.27 2.03 3.47 0.36 3.43 5.40

**Figure 4.** The results for the systems (A) ethanol–toluene at 101.3 kPa, (B) ethanol–toluene at 201.3 kPa, (c) isopropyl alcohol–toluene at 101.3 kPa, and (D) isopropyl alcohol–toluene at 201.3 kPa [20]. **Figure 5.** The results for the systems (A) ethanol-toluene at 101.3 kPa, (B) ethanol-toluene at 201.3 kPa, (c) isopropyl alcohol-toluene at 101.3 kPa, and (D) isopropyl alcohol-toluene at 201.3 kPa [20].

#### **4.2. VLE study of a ternary system**

It is worth noting that the main difference in Figures 3 and 4 is the use of parameter estimation method for the calculation of the NRTL-SAC parameters, while for the UNIFAC model, the calculation is straightforward. Once the parameters (here, the segment numbers) are found,

One of the main applications of the thermodynamic models is in the chemical industries which use solvent (or their mixtures) [19-22]. Two cases of the vapor-liquid equilibrium of common

There are some solvents within the chemical and pharmaceutical industries which are of importance to study, such as ethanol, isopropyl alcohol, and water in addition to some aromatic components, such as benzene and toluene. These solvents are sometimes used as additives to other valuable chemicals to maintain some performances or enhance their physical or chemical properties. The systems that form azeotropes make the distillation process calculations complex. As a result, having an accurate knowledge about the phase behaviors of such systems are important (such as toluene-ethanol or toluene-isopropyl alcohol). The four case studies of the VLE data for the mentioned azeotropic binary systems were used in a study by Chen et al. [24]. The operating pressures were at four distinct levels for the calculations. They have used three equations of state (from the correlative group) to fit the experimental data with the model outputs and found the necessary binary interaction parameters. Based on their work, the estimated parameters were found; however, they were not used to further validate the models for other operating conditions. In a study by Sheikholeslamzadeh et al., they used the NRTL-SAC and UNIFAC models to predict these azeotropic systems [20]. Table 1 (from their work) shows average relative deviation for the compositions and temperature for the two systems. It is seen from the table that the deviations from the experimental data in the vapor phase mole fractions are almost one-third relative to the NRTL-SAC model. On the other hand, the relative deviation for the saturation temperature is higher using the NRTL-SAC model. Another fact from this finding is that the deviations for both binary systems when the NRTL-SAC model is used are nearly the same. This is not the case when utilizing the UNIFAC model to predict the systems' behaviors. The final conclusion from Table 1 would be the better predictive capability of the UNIFAC model for light alcohols than the heavier ones.

> **%ARD (equilibrium temperature and vapor-phase mole fractions) System 1 System 2**

> > **(K) 2-Propanol Toluene**

**Ethanol Toluene Temperature,**

NRTL-SAC 101.3 0.46 5.14 9.81 0.41 4.91 7.19

then they could be used for validation against other experimental data.

**4. Application of solution thermodynamics in industry**

industrial solvent systems are discussed here.

18 Recent Advances in Thermo and Fluid Dynamics

**Thermodynamic model**

**Pressure, (kPa)**

**Temperature, (K)**

**4.1. VLE study of two binary azeotropic systems**

With the increasing prices of fossil fuels, the demand to obtain alternatives has received much attention. One of the candidates for this purpose is ethanol, as it decreases the air pollution when blended to the conventional fossil fuels and thereby, increasing the performance of burning fuels within the vehicle engines. It would also be cost-effective when adding other low-price additives to the fuel. The higher the purity of the alcohol being used as additive, the better the performance of the fuel. It was found that the addition of glycols to the mixture containing alcohols and water can improve the separation processes and utilize less energy to perform the process [25,26]. The experimental data containing the ternary system of ethylene glycol-water-ethanol and the performance of separation by varying the glycol concentration were performed by Kamihama et al. [27].

**Figure 5.** VLE diagrams of the systems (A) ethanol–water, (B) ethanol–ethylene glycol, and (C) water– ethylene glycol at 101.3 kPa [20]. **Figure 6.** VLE diagrams of the systems (A) ethanol-water, (B) ethanol-ethylene glycol, and (C) water-ethylene glycol at 101.3 kPa [20].

They performed binary system experiments for each of the pairs in the ternary system. It was found that glycol can move the azeotrope point and therefore, enhance the separation process. They performed binary system experiments for each of the pairs in the ternary system. It was found that glycol can move the azeotrope point and therefore, enhance the separation process. Sheikholeslamzadeh et al. have used both the NRTL-SAC and UNIFAC models to perform phase calculations and assess the capacity of the mentioned models in the prediction of binary and ternary systems containing glycol and alcohols [20].

Sheikholeslamzadeh et al. have used both the NRTL‐SAC and UNIFAC models to perform phase calculations and assess the capacity of the mentioned models in the prediction of binary and ternary From Figure 6, the UNIFAC model has the capability of predicting the system of ethanol-water, perfectly. However, this is not the case for the systems using ethylene glycol-water and

From Figure 5, the UNIFAC model has the capability of predicting the system of ethanol–water, perfectly. However, this is not the case for the systems using ethylene glycol–water and ethylene glycol– ethanol. On the other side, the NRTL‐SAC model gave satisfactory results for all three binary system, specifically, the systems that contain ethylene glycol. Also from Figures 5B and 5C, it can be seen that for nonazeotropic systems, the NRTL‐SAC model could best show capability in prediction comparable to the

systems containing glycol and alcohols [20].

ethylene glycol-ethanol. On the other side, the NRTL-SAC model gave satisfactory results for all three binary system, specifically, the systems that contain ethylene glycol. Also from Figures 6B and 6C, it can be seen that for nonazeotropic systems, the NRTL-SAC model could best show capability in prediction comparable to the UNIFAC model. The UNIFAC model could best locate the azeotrope point and the VLE behavior of the ethanol-water system.

low-price additives to the fuel. The higher the purity of the alcohol being used as additive, the better the performance of the fuel. It was found that the addition of glycols to the mixture containing alcohols and water can improve the separation processes and utilize less energy to perform the process [25,26]. The experimental data containing the ternary system of ethylene glycol-water-ethanol and the performance of separation by varying the glycol concentration

**Figure 5.** VLE diagrams of the systems (A) ethanol–water, (B) ethanol–ethylene glycol, and (C) water– ethylene glycol at 101.3 kPa [20].

**Figure 6.** VLE diagrams of the systems (A) ethanol-water, (B) ethanol-ethylene glycol, and (C) water-ethylene glycol at

They performed binary system experiments for each of the pairs in the ternary system. It was found that glycol can move the azeotrope point and therefore, enhance the separation process. Sheikholeslamzadeh et al. have used both the NRTL-SAC and UNIFAC models to perform phase calculations and assess the capacity of the mentioned models in the prediction of binary

<sup>0</sup> 0.2 0.4 0.6 0.8 <sup>1</sup> <sup>360</sup>

x,y (Water)

360

380

T em perature [K]

C

400

420

440

460

<sup>0</sup> 0.2 0.4 0.6 0.8 <sup>1</sup> <sup>340</sup>

Bubble points (experimental) Dew points (experimental) Bubble curve (UNIFAC) Dew curve (UNIFAC) Bubble curve (NRTL-SAC) Dew curve (NRTL-SAC)

x,y (Ethanol)

B

Bubble points (Experimental) Dew points (Experimental) Bubble curve (UNIFAC) Dew curve (UNIFAC) Bubble curve (NRTL-SAC) Dew curve (NRTL-SAC)

They performed binary system experiments for each of the pairs in the ternary system. It was found that glycol can move the azeotrope point and therefore, enhance the separation process. Sheikholeslamzadeh et al. have used both the NRTL‐SAC and UNIFAC models to perform phase calculations and assess the capacity of the mentioned models in the prediction of binary and ternary

From Figure 6, the UNIFAC model has the capability of predicting the system of ethanol-water, perfectly. However, this is not the case for the systems using ethylene glycol-water and

From Figure 5, the UNIFAC model has the capability of predicting the system of ethanol–water, perfectly. However, this is not the case for the systems using ethylene glycol–water and ethylene glycol– ethanol. On the other side, the NRTL‐SAC model gave satisfactory results for all three binary system, specifically, the systems that contain ethylene glycol. Also from Figures 5B and 5C, it can be seen that for nonazeotropic systems, the NRTL‐SAC model could best show capability in prediction comparable to the

were performed by Kamihama et al. [27].

20 Recent Advances in Thermo and Fluid Dynamics

A

Bubble points (Experimental) Dew points (Experimental) Bubble point (UNIFAC) Dew point (UNIFAC) Bubble point (NRTL-SAC) Dew point (NRTL-SAC)

<sup>0</sup> 0.2 0.4 0.6 0.8 <sup>1</sup> <sup>350</sup>

x,y (Ethanol)

380

and ternary systems containing glycol and alcohols [20].

400

Tem

perature [K]

420

440

460

480

355

101.3 kPa [20].

T e m pera tu re [K ]

360

365

370

systems containing glycol and alcohols [20].

For the ternary system of solvents consisting of ethanol, water, and ethylene glycol, Kamihama et al. [27] conducted the vapor-liquid measurements at a pressure 101.3 kPa. In order to use the correlative models (such as Wilson) for the ternary system, the binary interaction param‐ eters should be known for each pair of components in the mixture at that temperature and pressure. They used this method to find the ternary behavior of the system. Sheikholeslam‐ zadeh et al. used the NRTL-SAC model with the four conceptual segments of each solvent, which were already accessible in the literature [9, 18]. The results showed the high performance of the NRTL-SAC model in the prediction of this ternary system. The bubble point temperature as well as the vapor and liquid compositions could be estimated fairly well with the NRTL-SAC model. The predicted results are shown in Figure 7, giving the vapor phase mole fractions of the three species as well as the mixture temperature. It is apparent that the deviation from the experimental data for the case of water and ethanol in this ternary mixture is almost zero. The UNIFAC model could not match the experimental data as well as the NRTL-SAC model. For the ethylene glycol, as the experimental concentrations of the vapor phase are very small, with the inclusion of errors of the experimentation, the NRTL-SAC model could also give satisfactory results. Finally, for the bubble point temperature, Figure 7 illustrates the perfect predictions of the NRTL-SAC model compared with the UNIFAC model. The average relative deviation for the whole set of experimental data on the saturated temperature from the NRTL-SAC model is one-third that from the UNIFAC model.

**Figure 7.** Vapor phase mole fractions of ethylene glycol (bottom-left), ethanol (top-left), water (top-right), and the bub‐ ble temperature (bottom-right) prediction from the NRTL-SAC model [20].

#### **5. Conclusion**

There are several equations of state that describe the phase behavior of chemical components of a system at various temperatures, pressures, and compositions. From these models, the first group which needs various experimental data to predict the system behavior at other condi‐ tions is not very attractive. Instead, the second group (predictive models) is based on the activity coefficients that are found from the molecular structures with a few experimental data. In this chapter, the capacity of handling two binary and ternary systems of solvents using those predictive models was assessed. The NRTL-SAC and UNIFAC models were chosen for the modeling of those systems. The NRTL-SAC model showed relative advantage over the UNIFAC model in almost all cases, except for the systems containing light alcohols with water. The preference of using NRTL-SAC is due to its simplicity compared to the UNIFAC. If the four segment parameters of a specific component are known, then they can be set as a unique value for that component irrespective of the mixture conditions. This way, various operating conditions can be defined and the phase behavior of the components can be predicted accurately and rapidly. One of the main advantages of the NRTL-SAC model is that it can be written in various computer programming languages to be used in process simulation analysis. One good example of such work can be found in Sheikholeslamzadeh et al. [19,21]. They used the NRTL-SAC model to find the solid-liquid equilibrium for three pharmaceuticals. Then, the parameters they optimized for the given pharmaceuticals were used further to perform a solvent screening method. This method is really costly and time-consuming in industry. The authors [19,20] have developed an algorithm to find the best combination of solvents, temperatures, and pressures for the best yield of pharmaceutical production using the NRTL-SAC model. In conclusion, the NRTL-SAC model and other similar ones open a new window for the engineers and scientists to have a wider, more accurate, and rapid predictions of the solubility of active pharmaceuticals in mixtures of solvents.

#### **Author details**

Ehsan Sheikholeslamzadeh and Sohrab Rohani\*

\*Address all correspondence to: Srohani@uwo.ca

Department of Chemical and Biochemical Engineering, Western University, Canada

#### **References**

[1] J. M. Smith, Hendrick C Van Ness, Michael Abbott. Introduction to Chemical Engi‐ neering Thermodynamics, McGraw-Hill Science, Sixth Edition, 2004.

[2] J. M. Prausnitz, R. N. Lichtenthaler, E. G. Azevedo. Molecular Thermodynamics of Fluid-Phase Equilibria, Prentice-Hall International Series, Third Edition, 1999.

**5. Conclusion**

22 Recent Advances in Thermo and Fluid Dynamics

**Author details**

**References**

There are several equations of state that describe the phase behavior of chemical components of a system at various temperatures, pressures, and compositions. From these models, the first group which needs various experimental data to predict the system behavior at other condi‐ tions is not very attractive. Instead, the second group (predictive models) is based on the activity coefficients that are found from the molecular structures with a few experimental data. In this chapter, the capacity of handling two binary and ternary systems of solvents using those predictive models was assessed. The NRTL-SAC and UNIFAC models were chosen for the modeling of those systems. The NRTL-SAC model showed relative advantage over the UNIFAC model in almost all cases, except for the systems containing light alcohols with water. The preference of using NRTL-SAC is due to its simplicity compared to the UNIFAC. If the four segment parameters of a specific component are known, then they can be set as a unique value for that component irrespective of the mixture conditions. This way, various operating conditions can be defined and the phase behavior of the components can be predicted accurately and rapidly. One of the main advantages of the NRTL-SAC model is that it can be written in various computer programming languages to be used in process simulation analysis. One good example of such work can be found in Sheikholeslamzadeh et al. [19,21]. They used the NRTL-SAC model to find the solid-liquid equilibrium for three pharmaceuticals. Then, the parameters they optimized for the given pharmaceuticals were used further to perform a solvent screening method. This method is really costly and time-consuming in industry. The authors [19,20] have developed an algorithm to find the best combination of solvents, temperatures, and pressures for the best yield of pharmaceutical production using the NRTL-SAC model. In conclusion, the NRTL-SAC model and other similar ones open a new window for the engineers and scientists to have a wider, more accurate, and rapid predictions

of the solubility of active pharmaceuticals in mixtures of solvents.

Department of Chemical and Biochemical Engineering, Western University, Canada

neering Thermodynamics, McGraw-Hill Science, Sixth Edition, 2004.

[1] J. M. Smith, Hendrick C Van Ness, Michael Abbott. Introduction to Chemical Engi‐

Ehsan Sheikholeslamzadeh and Sohrab Rohani\*

\*Address all correspondence to: Srohani@uwo.ca


#### **Chapter 2**

## **Dynamics of Droplets**

[17] A. T. Kan and M. B. Tomson. UNIFAC prediction of aqueous and nonaqueous solu‐ bilities of chemicals with environmental interest. Environ. Sci. Technol. 1996, 30,

[18] C. C. Chen and Y. Song. Extension of nonrandom two-liquid segment activity coeffi‐

[19] Ehsan Sheikholeslamzadeh and Sohrab Rohani. Solubility prediction of pharmaceuti‐ cal and chemical compounds in pure and mixed solvents using predictive models.

[20] Ehsan Sheikholeslamzadeh and Sohrab Rohani. Vapour-liquid and vapour-liquidliquid equilibrium modeling for binary, ternary, and quaternary systems of solvents.

[21] Ehsan Sheikholeslamzadeh, Chau-Chyun Chen, and Sohrab Rohani. Optimal solvent screening for the crystallization of pharmaceutical compounds from multisolvent

[22] Ehsan Sheikholeslamzadeh and Sohrab Rohani. Modeling and optimal control of sol‐ ution mediated polymorphic transformation of l-glutamic acid. Ind. Eng. Chem. Res.

[23] Stichlmair, J., Fair, J., Bravo, J. L. Separation of azeotropic mixtures via enhanced dis‐

[24] Chen, R., Zhong, L., Xu, C. Isobaric vapor-liquid equilibrium for binary systems of toluene + ethanol and toluene + isopropanol at (1013.3, 121.3, 161.3, and 201.3) kPa. J.

[25] Pequenín, A., Carlos Asensi, J., Gomis V. Isobaric vapor−liquid−liquid equilibrium and vapor−liquid equilibrium for the quaternary system water−ethanol−cyclohexane

[26] Pequenín, A., Carlos Asensi, J., Gomis V. Experimental determination of quaternary and ternary isobaric vapor-liquid-liquid equilibrium and vapor-liquid equilibrium for the systems water-ethanol-hexane-toluene and water-hexane-toluene at 101.3

[27] Kamihama, N., Matsuda, H., Kurihara, K., Tochigi, K., Oba, S. Isobaric vapor-liquid equilibria for ethanol + water + ethylene glycol and its constituent three binary sys‐

cient model for electrolytes. Ind. Eng. Chem. Res. 2005, 44, 8909-8921.

Ind. Eng. Chem. Res. (2011) 51, 464-473.

tillation. Chem. Eng. Prog. B 1989, 85, 63.

kPa. J. Chem. Eng. Data, 2010, 56, 3991.

Chem. Eng. Data. 2011, 57, 155.

systems. Ind. Eng. Chem. Res. (2012) 51, 13792-13802.

−isooctane at 101.3 kPa. J. Chem. Eng. Data, 2010, 55, 1227.

tems. J. Chem. Eng. Data. dx.doi.org/10.1021/je2008704.

Fluid Phase Equilibria, 333, 97-105.

(2013) 52, 2633-2641.

1369-1376.

24 Recent Advances in Thermo and Fluid Dynamics

Hossein Yahyazadeh and Mofid Gorji-Bandpy

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/61901

#### **Abstract**

Capturing non-Newtonian power-law drops by horizontal thin fibers with circular cross‐ section in a quiescent media can be studied in this chapter. The case is simulated using volume of fluid (VOF) method providing a notable reduction of a computational cost. Open source OpenFOAM software is applied to conduct the simulations. This model is an extension of the one developed earlier by Lorenceau, Clanet, and Quéré [1]. To vali‐ date the model, water drops affecting a fiber of radius 350μm were simulated and thresh‐ old drop radiuses were obtained regarding to the impact velocity. These results agreed well with the experimental data presented by Lorenceau et al. [1]. In the next step, non-Newtonian power-law drops landing on thin fiber of radius 350μm were simulated. The final goal of this study was to obtain the threshold velocity and radius of a drop that is completely captured by the fiber. Threshold radiuses for both shear-thinning and shearthickening drops were obtained and compared with corresponding Newtonian drops. Results show that the threshold radius of drop increases in a fixed velocity as *n*, powerlaw index, increases. Furthermore, shear-thinning nature of the drop leads to instabilities in high Reynolds numbers (Re) as it influences the fiber.

**Keywords:** Impact, Wetting, Filtration, Non-Newtonian Drop, Power-Law Model

#### **1. Introduction**

The problem of the removal of aerosol particles from gas streams has become of increasing importance from the standpoint of public health and the recovery of valuable products. Technology of controlling the aerosol particles or improving the liquid phase of aerosol is very important in many industrial processes such as oil and petroleum, electronic, mining, and food, as well as waste products like noxious emission of aerosol in chemical plants. There are several ways for this purpose among which fibrous filters are more popular so that it is obvious to try to improve their efficiency. The efficiency of collection and the pressure drop are the most important practical considerations in the design of these fibrous filters [2]. Various

© 2015 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

experimental studies on liquid aerosol filtration have shown a filter clogging in several stages leading to a drainage stage with a constant pressure drop [3, 4]. Contal et al. [3] introduced four stages describing the clogging of fibrous filters by liquid droplets: In the first stage, the droplets impact the wet fibers. In the second stage, the amount of those droplets captured by fiber increases so that neighbor droplets coalesce. In the third stage, liquid shells form in the net that leads to a large increase of the pressure drop and to a dramatic decrease in the efficiency of the filter. Finally, there will be a pseudo-stationary state between the droplet collection and gravitational drainage of the liquid phase [1]. The first stage of clogging is studied in many researches experimentally. Patel et al. [5] investigated the drop breakup in a flow through fiber filters and the breakup probability for a drop. Hung and Yao [6] evaluated the impact of water droplets on a cylindrical object experimentally. Both these studies show that depending on the speed of the impacting droplets, the drops can be captured and broken into several pieces. Lorenceau et al. [1] investigated the threshold velocity and radius of drops captured by thin fibers. They have shown that the threshold impact velocity of drops decreases as the drop radius increases and vice versa. Following this study, Lorenceau et al. [7] investigated offcenter impact of droplets on horizontal fibers.

Most of liquid droplets in aerosols produced in different industrial and natural processes show non-Newtonian behavior. Therefore, studying the phenomenon of impaction of a non-Newtonian droplet on thin fibers is of prime importance. All mentioned studies have been devoted to Newtonian fluids. Colliding of non-Newtonian droplets on thin fibers, to our knowledge, has not been investigated previously. However, colliding of non-Newtonian drops on flat plates is evaluated in some experimental and numerical studies. Haeri and Hashema‐ badi [8] studied spreading of power-law fluids on inclined plates both numerically and experimentally. Saïdi et al. [9] investigated experimentally the influence of yield stress on the fluid droplet impact control. Son and Kim [10] studied experimentally the spreading of inkjet droplet of non-Newtonian fluid on a solid surface with controlled contact angle at low Weber and Reynolds numbers. Kim and Baek [11] investigated the parameters that govern the impact dynamics of yield stress on the fluid droplets on solid surfaces.

In this chapter, we focus on the first stage of clogging and investigate the impaction of non-Newtonian power-law fluids on thin fibers. We aim to obtain the threshold radius of impacting droplets in different impact velocities. Effect of shear-thinning and shear-thickening behavior of droplets is evaluated and compared with corresponding Newtonian fluids. For this purpose, volume of fluid method is used and open source OpenFOAM software is applied for simula‐ tions.

The order of contents in this chapter is as follows: In Section 2 the governing equations and numerical methodology are given. Validation of the obtained results for Newtonian fluids in comparison with experimental observations presented by Lorenceau et al. [1] is provided in Section 3.1. A general description of the observations for non-Newtonian droplets captured by horizontal fiber is presented in Section 3.2, and the behavior of non-Newtonian drops affecting thin fibers is explained. Finally, a summary of the results and conclusions is provided.

#### **2. Methodology**

experimental studies on liquid aerosol filtration have shown a filter clogging in several stages leading to a drainage stage with a constant pressure drop [3, 4]. Contal et al. [3] introduced four stages describing the clogging of fibrous filters by liquid droplets: In the first stage, the droplets impact the wet fibers. In the second stage, the amount of those droplets captured by fiber increases so that neighbor droplets coalesce. In the third stage, liquid shells form in the net that leads to a large increase of the pressure drop and to a dramatic decrease in the efficiency of the filter. Finally, there will be a pseudo-stationary state between the droplet collection and gravitational drainage of the liquid phase [1]. The first stage of clogging is studied in many researches experimentally. Patel et al. [5] investigated the drop breakup in a flow through fiber filters and the breakup probability for a drop. Hung and Yao [6] evaluated the impact of water droplets on a cylindrical object experimentally. Both these studies show that depending on the speed of the impacting droplets, the drops can be captured and broken into several pieces. Lorenceau et al. [1] investigated the threshold velocity and radius of drops captured by thin fibers. They have shown that the threshold impact velocity of drops decreases as the drop radius increases and vice versa. Following this study, Lorenceau et al. [7] investigated off-

Most of liquid droplets in aerosols produced in different industrial and natural processes show non-Newtonian behavior. Therefore, studying the phenomenon of impaction of a non-Newtonian droplet on thin fibers is of prime importance. All mentioned studies have been devoted to Newtonian fluids. Colliding of non-Newtonian droplets on thin fibers, to our knowledge, has not been investigated previously. However, colliding of non-Newtonian drops on flat plates is evaluated in some experimental and numerical studies. Haeri and Hashema‐ badi [8] studied spreading of power-law fluids on inclined plates both numerically and experimentally. Saïdi et al. [9] investigated experimentally the influence of yield stress on the fluid droplet impact control. Son and Kim [10] studied experimentally the spreading of inkjet droplet of non-Newtonian fluid on a solid surface with controlled contact angle at low Weber and Reynolds numbers. Kim and Baek [11] investigated the parameters that govern the impact

In this chapter, we focus on the first stage of clogging and investigate the impaction of non-Newtonian power-law fluids on thin fibers. We aim to obtain the threshold radius of impacting droplets in different impact velocities. Effect of shear-thinning and shear-thickening behavior of droplets is evaluated and compared with corresponding Newtonian fluids. For this purpose, volume of fluid method is used and open source OpenFOAM software is applied for simula‐

The order of contents in this chapter is as follows: In Section 2 the governing equations and numerical methodology are given. Validation of the obtained results for Newtonian fluids in comparison with experimental observations presented by Lorenceau et al. [1] is provided in Section 3.1. A general description of the observations for non-Newtonian droplets captured by horizontal fiber is presented in Section 3.2, and the behavior of non-Newtonian drops affecting thin fibers is explained. Finally, a summary of the results and conclusions is provided.

center impact of droplets on horizontal fibers.

26 Recent Advances in Thermo and Fluid Dynamics

tions.

dynamics of yield stress on the fluid droplets on solid surfaces.

#### **2.1. Volume of fluid method**

Hirt and Nichols [12] demonstrated the volume of fluid (VOF) method and started a new trend in multiphase flow simulation. It relies on the definition of an indicator function *γ*. This function allows us to know whether one fluid or another occupies the cell, or a mix of both. In the conventional volume of fluid method [12], the transport equation for an indicator function *γ*, representing the volume fraction of one phase, is solved simultaneously with the continuity and momentum equations as follows:

$$\nabla \mathbf{J} \mathbf{U} = \mathbf{0} \tag{1}$$

$$
\frac{\partial \chi}{\partial t} + \nabla \cdot \left(\mathbf{U}\boldsymbol{\gamma}\right) = 0 \tag{2}
$$

$$\frac{\partial \left(\rho \mathbf{U}\right)}{\partial t} + \nabla \cdot \left(\rho \mathbf{U} \mathbf{U}\right) = -\nabla \cdot p + \nabla .\mathbf{T} + \rho \mathbf{f}\_{\mathbf{b}} \tag{3}$$

where **U** represents the velocity field shared by the two fluids throughout the flow domain, *γ* is the phase fraction, **T**=2*μ* **S** − 2*μ* (∇.**U**)**I**/ 3 is the deviatoric viscous stress tensor, with the mean rate of strain tensor **S**=0.5 <sup>∇</sup>**<sup>U</sup>** <sup>+</sup> (∇**U**)*<sup>T</sup>* and **I**=*δi*, *<sup>j</sup>* , *ρ* is density, *p* is pressure, and **f***b* are body forces per unit mass. In VOF simulations, the latter forces include gravity and surface tension effects at the interface. The phase fraction *γ* can take values within the range 0≤*γ* ≤1, with the values of zero and one corresponding to regions accommodating only one phase, for example, *γ* =0 for gas and *γ* =1 for liquid. Accordingly, gradients of the phase fraction are encountered only in the region of the interface. Two immiscible fluids are considered as one effective fluid throughout the domain, the physical properties of which are calculated as weighted averages based on the distribution of the liquid volume fraction, thus being equal to the properties of each fluid in their corresponding occupied regions and varying only across the interface:

$$
\rho = \rho\_1 \,\gamma + \rho\_{\overline{\mathbf{g}}} \left( 1 - \gamma \right) \tag{4}
$$

$$
\mu = \mu\_1 \gamma + \mu\_\mathbf{g} \left( 1 - \gamma \right) \tag{5}
$$

where *ρ*l and *ρg* are densities of liquid and gas, respectively.

In this study, a modified approach similar to the one proposed in [13] is used.

#### **2.2. Power-law model**

Generally, the form of the function relating *τxy* to shear rate *γ*˙ is quite complicated. It has been found, however, that the two-parameter equation of state

$$
\tau\_{\infty y} = K \dot{\gamma}^{\mathcal{N}} \tag{6}
$$

is adequate for many non-Newtonian fluids [16], where *K* is the fluid consistency coefficient and *n* is power-law index. This is a well-known power-law model, which is used in this study.

#### **2.3. Computational method**

All computations are performed using the code OpenFOAM [17], an open source computa‐ tional fluid dynamics (CFD) toolbox, utilizing a cell-center-based finite volume method on a fixed unstructured numerical grid and employing the solution procedure based on the pressure implicit with splitting of operators (PISO) algorithm for coupling between pressure and velocity in transient flows [18].

A grid spacing equal to 1/20 of the droplet radius was used to discretize the computational domain. The mesh size was determined based on a mesh refinement study in which the grid spacing was progressively decreased until further reductions made no significant change in the predicted droplet shape during the impact. Bussmann et al. [19] have provided a detailed description of such a mesh refinement study earlier. The mesh size was uniform throughout the computational domain, which included the droplet and fiber.

#### **3. Results and discussion**

#### **3.1. Validation of the Solution**

In this section, we present the results obtained for water droplets impacting a fiber of radius 350 μm. Physical properties of water are given in Table 1. Threshold radiuses of the droplets in different impact velocities are obtained and compared with those exhibited by Lorenceau et al. [1]. All obtained results are shown in Figure 1, which demonstrates the threshold radiuses in different impact velocities.


**Table 1.** Physical properties of fluids

**Figure 1.** Threshold radiuses of the water droplets in different impact velocities

In this study, a modified approach similar to the one proposed in [13] is used.

t

found, however, that the two-parameter equation of state

the computational domain, which included the droplet and fiber.

Generally, the form of the function relating *τxy* to shear rate *γ*˙ is quite complicated. It has been

 g = & *n*

is adequate for many non-Newtonian fluids [16], where *K* is the fluid consistency coefficient and *n* is power-law index. This is a well-known power-law model, which is used in this study.

All computations are performed using the code OpenFOAM [17], an open source computa‐ tional fluid dynamics (CFD) toolbox, utilizing a cell-center-based finite volume method on a fixed unstructured numerical grid and employing the solution procedure based on the pressure implicit with splitting of operators (PISO) algorithm for coupling between pressure

A grid spacing equal to 1/20 of the droplet radius was used to discretize the computational domain. The mesh size was determined based on a mesh refinement study in which the grid spacing was progressively decreased until further reductions made no significant change in the predicted droplet shape during the impact. Bussmann et al. [19] have provided a detailed description of such a mesh refinement study earlier. The mesh size was uniform throughout

In this section, we present the results obtained for water droplets impacting a fiber of radius 350 μm. Physical properties of water are given in Table 1. Threshold radiuses of the droplets in different impact velocities are obtained and compared with those exhibited by Lorenceau et al. [1]. All obtained results are shown in Figure 1, which demonstrates the threshold radiuses

*n ρ* (*kg* / *m*3) *k* (*mPa* ⋅ *s <sup>n</sup>*) *σ* (*mN* / *m*) fluid 1 1000 1 70 water

*xy K* (6)

**2.2. Power-law model**

28 Recent Advances in Thermo and Fluid Dynamics

**2.3. Computational method**

and velocity in transient flows [18].

**3. Results and discussion**

**3.1. Validation of the Solution**

in different impact velocities.

**Table 1.** Physical properties of fluids

Defining *R*M as characteristic radius of the drop at zero impact velocity as follows:

$$R\_{\mathbf{M}} = 3^{1/3} b^{1/3} k^{-2/3} \tag{7}$$

where *k* <sup>−</sup><sup>1</sup> = *σ* / *ρg* is the capillary length and *b* is the radius of the fiber, and also characteristic velocity of the drop as:

$$\mathbf{M\_M} = \sqrt{4gR\_{\mathbf{M}}} \tag{8}$$

Variation of the dimensionless threshold radius versus dimensionless impact velocity is plotted in Figure 2, which shows a good agreement with the experimental and theoretical data presented by Lorenceau et al. [1].

**Figure 2.** Comparison of the obtained results with experimental and theoretical data

Errors of the achieved results in comparison with experimental results of Lorenceau et al. (2004) are presented in Table 2.


**Table 2.** Errors of the achieved results in comparison with the experimental and theoretical data

#### **3.2. Results for non-newtonian droplets**

Defining *R*M as characteristic radius of the drop at zero impact velocity as follows:

velocity of the drop as:

30 Recent Advances in Thermo and Fluid Dynamics

presented by Lorenceau et al. [1].

are presented in Table 2.

where *k* <sup>−</sup><sup>1</sup> = *σ* / *ρg* is the capillary length and *b* is the radius of the fiber, and also characteristic

Variation of the dimensionless threshold radius versus dimensionless impact velocity is plotted in Figure 2, which shows a good agreement with the experimental and theoretical data

**Figure 2.** Comparison of the obtained results with experimental and theoretical data

Errors of the achieved results in comparison with experimental results of Lorenceau et al. (2004)


<sup>=</sup> <sup>4</sup> M M *U gR* (8)

Four forces, capillary, gravity, inertia, and viscosity, are important during the impaction of droplet to the fiber. Capillary and gravity forces are in contrast. While the weight of the drop tends to detach it from the fiber, surface tension is responsible for the sustentation. In order to analyze the dynamics of the capture of a drop by a fiber, viscous and inertia effects must be considered. Here again, those effects are antagonistic. While inertia tends to make the drop cross the fiber, viscous effects will induce a dissipation and can thus possibly stop the drop [1]. Effects of variation of viscosity were not investigated in the previous works. In order to investigate these effects, both shear-thinning and shear-thickening fluids will be considered.

At first, we present general observations of threshold radiuses at different velocities for shearthinning droplets, *n* <1. For this purpose, physical properties of the fluid given in Table 1 were used except the power-law index, which is equal to 0.5. Figure 3 represents the threshold radiuses in different velocities for this type of droplet.

Next, we present our observations for the two kinds of shear-thickening droplets, *n*=1.5, 2, impacting the fiber which are shown in Figures 4 and 5.

Results show that, on the one hand, in shear-thinning drops, the threshold radius decreased in a fixed velocity in comparison with the corresponding Newtonian drops. Moreover, instabilities were observed during the impaction of the drop. This is because viscosity in shearthinning fluids decreases as it has an inverse relation with shear rate. That is, as a shearthinning droplet impacts a thin fiber, due to high shear rates, viscosity decreases in comparison with the corresponding Newtonian fluid, so that the viscosity force reduces (i.e., a force that helps the drop to stick on the fiber is reduced). Thus, in order to balance the antagonistic forces, drop size has to be decreased, which means the threshold radius has to be decreased. On the other hand, in the case where the drop is a shear-thickening fluid, the threshold radius increased in a fixed velocity in comparison with the corresponding Newtonian fluid, since the viscosity of shear-thickening fluids has a direct relation to the shear rate. It means that, viscosity increases as the drop impacts the fiber compared to Newtonian drops. Thereby, in shear-

**Figure 3.** Threshold radiuses of the shear-thinning (*n*=0.5) droplets in different impact velocities

thickening drops, viscous forces increase and the drop's tendency to stick to the fiber increases. Therefore, the threshold drop size will also increase. This tendency is obvious in Figures 4 and 5. Threshold radius of Newtonian, shear-thinning and shear-thickening droplets versus impact velocity is plotted in Figure 6.

Figure 6 shows that the threshold radius will decrease generally with the increase of an impact velocity or vice versa. In shear-thinning droplets, viscosity has an inverse relation with shear rate, as mentioned before, so that it will decrease with the increase of impact velocity. It means inertia forces increase whereas viscosity forces decrease. That is why instabilities are observed in shear-thinning droplets impacting fiber at high speeds. The threshold radius of shearthickening droplets, also, decreases with the increase of impact velocity. In shear-thickening droplets, although there is a direct relation between the viscosity and the shear rate, velocity increases with second power compared to capillary forces. That is, increase of inertia forces is

**Figure 4.** Threshold radiuses of the shear-thickening (*n*=1.5) droplets in different impact velocities

thickening drops, viscous forces increase and the drop's tendency to stick to the fiber increases. Therefore, the threshold drop size will also increase. This tendency is obvious in Figures 4 and 5. Threshold radius of Newtonian, shear-thinning and shear-thickening droplets versus impact

**Figure 3.** Threshold radiuses of the shear-thinning (*n*=0.5) droplets in different impact velocities

Figure 6 shows that the threshold radius will decrease generally with the increase of an impact velocity or vice versa. In shear-thinning droplets, viscosity has an inverse relation with shear rate, as mentioned before, so that it will decrease with the increase of impact velocity. It means inertia forces increase whereas viscosity forces decrease. That is why instabilities are observed in shear-thinning droplets impacting fiber at high speeds. The threshold radius of shearthickening droplets, also, decreases with the increase of impact velocity. In shear-thickening droplets, although there is a direct relation between the viscosity and the shear rate, velocity increases with second power compared to capillary forces. That is, increase of inertia forces is

velocity is plotted in Figure 6.

32 Recent Advances in Thermo and Fluid Dynamics

higher than increase of viscosity forces so that the threshold radius will decrease to make the drop stick to the fiber. However, in higher velocities, this rate will decrease as it can be seen in Figure 6. In other words, rate of the decreasing of threshold radius with the increasing of impact velocity will decrease due to high growth of shear rate.


**Figure 5.** Threshold radiuses of the shear-thickening (*n*=2) droplets in different impact velocities

There are some general observations that are common to Newtonian and non-Newtonian fluids. For all kinds of fluids, the threshold radius of the droplet has a reverse relation with the impact velocity of the droplet. In all cases, at a fixed impact velocity, droplets with a radius greater than the threshold radius have passed the fiber without breakup, and drops with a radius lower than the threshold radius clung to the fiber entirely.

#### **4. Conclusion**

Impaction of non-Newtonian power-law droplet to the horizontal fiber of circular cross section is investigated in this study. Volume of fluid technique is employed, significantly reducing the computational cost. Outcomes are divided into three parts: First, it has been observed that the threshold radius of the droplets decreased with the increase of impact velocity for New‐

**Figure 6.** Threshold radiuses of droplets versus impact velocity for Newtonian, shear-thinning and shear-thickening fluids

tonian, shear-thinning and shear-thickening fluids, due to the growth of inertia forces versus viscosity and capillary forces. Second, instabilities have been observed in shear-thinning droplets at high impact velocities due to severe reduction of viscosity. Finally, in shearthickening droplets, at high impact speeds, rate of reduction of threshold radius has decreased due to increase of the viscosity forces.

#### **Author details**

There are some general observations that are common to Newtonian and non-Newtonian fluids. For all kinds of fluids, the threshold radius of the droplet has a reverse relation with the impact velocity of the droplet. In all cases, at a fixed impact velocity, droplets with a radius greater than the threshold radius have passed the fiber without breakup, and drops with a

Impaction of non-Newtonian power-law droplet to the horizontal fiber of circular cross section is investigated in this study. Volume of fluid technique is employed, significantly reducing the computational cost. Outcomes are divided into three parts: First, it has been observed that the threshold radius of the droplets decreased with the increase of impact velocity for New‐

radius lower than the threshold radius clung to the fiber entirely.

**Figure 5.** Threshold radiuses of the shear-thickening (*n*=2) droplets in different impact velocities

**4. Conclusion**

34 Recent Advances in Thermo and Fluid Dynamics

Hossein Yahyazadeh and Mofid Gorji-Bandpy\*

\*Address all correspondence to: gorji@nit.ac.ir

Department of Mechanical Engineering, Babol University of Technology, Babol, Iran

#### **References**

[1] É. Lorenceau, C. Clanet, D. Quéré, Capturing drops with a thin fiber, Journal of Col‐ loid and Interface Science, vol. 279, 192–197, 2004.


[15] J.U. Brackbill, D. B. Kothe, C. Zemach, A continuum method for modeling surface tension, Journal of Computational Physics, vol. 100, 335–354, 1992.

[2] C.Y. Chenz, Filtration Of Aerosols By Fibrous Media,DA18-108-CML-4789 with the

[3] P. Contal, J. Simao, D. Thomas, T. Frising, S. Call, J.C. Appert-Collin, D. Bemer, Clog‐ ging of fibre filters by submicron droplets. Phenomena and influence of operating

[4] D.C. Walsh, J.I.T. Stenhouse, K.L. Scurrah, A. Graef, The effect of solid and liquid aerosol particle loading on fibrous filter material performance, Journal of Aerosol

[5] P. Patel, E. Shaqfeh, J.E. Butler, V. Cristini, J. Blawzdziewicz, M. Loewenberg, Drop breakup in the flow through fixed fiber beds: An experimental and computational in‐

[6] L.S. Hung, S.C. Yao, Experimental investigation of the impaction of water droplets on cylindrical objects, International Journal of Multiphase Flow, vol. 25, 1545–1559,

[7] E. Lorenceau, C. Clanet, D. Quéré, M. Vignes-Adler, Off-centre impact on a horizon‐

[8] S. Haeri, S.H. Hashemabadi, Three dimensional CFD simulation and experimental study of power law fluid spreading on inclined plates, International Communica‐

[9] Alireza Saïdi, Céline Martin, Albert Magnin, Influence of yield stress on the fluid droplet impact control, Journal of Non-Newtonian Fluid Mechanics, vol. 165, Issues

[10] Yangsoo Son, Chongyoup Kim, Spreading of inkjet droplet of non-Newtonian fluid on solid surface with controlled contact angle at low Weber and Reynolds numbers,

[11] Eunjeong Kim, Jehyun Baek, Numerical study of the parameters governing the im‐ pact dynamics of yield-stress fluid droplets on a solid surface, Journal of Non-New‐

[12] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free

[13] H. Rusche, Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions, Imperial College of Science, Technology and Medicine, Ph.D. thesis,

[14] G. Černe, S. Petelin, I. Tiselj, Coupling of the interface tracking and the two-fluid models for the simulation of incompressible two-phase flow, Journal of Computa‐

boundaries, Journal of Computational Physics, vol. 39, 201–225, 1981.

Journal of Non-Newtonian Fluid Mechanics, vol. 162, Issues 1–3, 78–87, 2009.

tal fibre, European Physical Journal Special Topics, vol. 166, 3–6, 2009.

tions in Heat and Mass Transfer, vol. 35, 1041–1047, 2008.

tonian Fluid Mechanics, vol. 173–174, 62–71, 2012.

tional Physics, vol. 171, 776–804, 2001.

Chemical Corps, U.S. Army, Washington 25, D.C, March 6, 1955.

conditions, Aerosol Science, vol. 35, 263–278, 2004.

vestigation, Physics of Fluids, vol. 15, 1146, 2003.

Science, vol. 27, 617–618, 1996.

36 Recent Advances in Thermo and Fluid Dynamics

1999.

2002.

11–12, 596–606, 2010.

