New Approaches to the Modeling and Calculation of Supercritical

Oct 19, 2015 - A 2D model of expanded supercritical fluid jet and unique calculation algorithms were developed and applied for the calculations...
0 downloads 0 Views 941KB Size
Download PDF
Article pubs.acs.org/IECR

New Approaches to the Modeling and Calculation of Supercritical Expansion and Micro- and Nanoparticle Formation V. I. Anikeev* and D. A. Stepanov Boreskov Institute of Catalysis, Siberian Branch of the Russian Academy of Sciences, Novosibirsk 630090, Russia ABSTRACT: Mathematical modeling of the rapid expansion of supercritical solution (nitrogen, nitrogen−water, propane− methanol−phenanthrene) and formation of particles was performed. A 2D model of expanded supercritical fluid jet and unique calculation algorithms were developed and applied for the calculations. The calculations allowed determination of jet boundaries, shape, and position of normal and barrel shock waves, condensation front, and particles size distribution. Parameter sensitivity analysis of the model helped to identify which parameters are of key importance for the production of particles with the required size and properties.

1. INTRODUCTION The present work is an indirect continuation of previous studies of rapid expansion of supercritical solutions (RESS), condensation of a solute, and formation of solid particles.1 In the latter work, the jet expansion was modeled in 1D approximation with preset shape of boundaries and shock waves. In the present work, the modeling is focused on the realistic jet expansion process. First, a 2D model for super- and subsonic jet expansion regions is considered. Second, the jet boundaries and shape and position of direct and hovering shock waves are calculated by solving the model equations. Third, the particle size distribution is considered, and a new model for particle growth during jet expansion is suggested. New approaches to the modeling of supercritical fluid expansion, particle nucleation and growth allowed more reliable optimization of the process parameters to provide the needed properties and size of the formed solid phase particles. In ref 1, we analyzed in detail the potentials, advantages, and shortcomings of existing methods to synthesize nanosized particles using supercritical fluids and explained the choice of the rapid expansion method. The problems related to particle formation by rapid expansion of supercritical solutions are thoroughly discussed in, for example, refs 2−5.

Figure 1. Principal scheme of submerged outflow of a fluid. I, reactor; II, capillary; III, supersonic expansion zone; IV, Mach disc; V, subsonic expansion zone.

sufficient velocity, including that in the radial direction which induces the formation of complex flux profile with expansion and contraction regions and complicated shock waves.6 The radial component of the flux velocity can vary considerably in the vicinity of the jet boundary and change its direction several times, until it becomes negligible because of dissipative effects. For this reason, the jet boundary downstream the capillary exit can form one or more barrel-shaped structures−hovering shock wave. At high outflow abnormality,6,7 one barrel shock is usually formed. The barrel-shaped nonisobaric jet boundaries enable the formation of similarly configured shockwaves, which, being reflected from the flow axis, might induce the formation of Mach configurations with straight trailing shock waves−normal shock waves. For convenience, the flow is conditionally divided into three zones: initial, III (supersonic expansion region); transient, IV (the region in vicinity of Mach disc); and the main flow, V (subsonic expansion region; Figure 1). The initial zone (supersonic expansion) is a barrel adjacent to the capillary cut. In this zone, the most nonuniform distribution of gas-

2. THEORETICAL ANALYSIS 2.1. Mathematical Model and Solution Methods. The general task is the modeling of the submerged efflux of supercritical fluid from the reactor through a capillary. The task is sequentially solved in an ideal gas approximation and then in real gas approximation. Figure 1 illustrates the rapid expansion of supercritical fluid through a capillary. Here, the supercritical state of the solution in the reactor is reached at temperature T0 and pressure p0. Then, the supercritical solution is throttled through the capillary to an external space with temperature T∞ and pressure p∞. When the fluid flows through a capillary and expands in the outer volume, its thermodynamic and hydrodynamic characteristics vary. Submerged efflux of the fluid is described by a nonisobaric underexpanded jet model. Just after capillary cut (exit), the nonisobaric flow acquires © 2015 American Chemical Society

Received: Revised: Accepted: Published: 11438

June 17, 2015 September 29, 2015 October 19, 2015 October 19, 2015 DOI: 10.1021/acs.iecr.5b02180 Ind. Eng. Chem. Res. 2015, 54, 11438−11447

Article

Industrial & Engineering Chemistry Research

The sets of eqs 1 and 2 in general form describe also the flow with condensation of saturated or supercooled vapor in the expansion wave sheaf. The main task to be solved in the case of a multicomponent solution expansion accompanied by spontaneous condensation of a component is to determine the position and shape of the front edge of condensation shock, its intensity, shock losses, and the volume of condensed phase downstream of the condensation shock wave. To solve this task, let us divide the whole flow field into several fluid tubes so that the flow in one tube can be assumed to be one-dimensional. The angle of the tube flow deflection and flow parameters upstream of expansion shock wave are known. According to ref 10, in general, any thermodynamic parameter Ω (p, ρ, or T) obeys the following relation:

dynamic parameters both in longitudinal and transverse directions is observed. In the main zone (subsonic expansion), the wave processes become weaker, the pressure equalizes, and the flow acquires isobaric character. Between the initial and main zones, a so-called transient region is located in slightly nonisobaric flow in the vicinity of Mach disc. Parametric instability here is much lower than in the initial region; the shockwaves are less intensive. In general, the state of the moving fluid with the known thermodynamic characteristics is determined by the fluid velocity, density, and pressure as functions of coordinates and time. In the present work, we compute nonisobaric nonisothermal jet using a method suggested in ref 7, which uses explicit difference scheme of the second order of accuracy and smoothing procedure. Consider a supersonic jet in the region of G(x) ≤ y ≤ F(x). In cross section x = xa (capillary cut), the gas-dynamic parameters are assumed to be known. Some boundary conditions are set at lines y1 = G(x) and y2 = F(x). Free surfaces, shockwaves, solid walls, symmetry axis, and so on may act as the y1 = G(x) and y2 = F(x) boundaries. Transform the region of x ≥ x0, G(x) ≤ y ≤ F(x) into half-strip 0 ≤ ξ ≤ 1 by replacing y with ξ using expression ξ=

⎛ ∂Ω ⎞ ⎛ ∂Ω ⎞ ⎛ ∂Ω ⎞ ⎜ ⎟ = − rρW ⎜ ⎟ ⎟ +⎜ ⎝ ∂Θ ⎠r ⎝ ∂Θ ⎠ψ ⎝ ∂ψ ⎠ Θ

where ψ is the stream line, ρ is the flow density, W is the flow velocity, r is the radius vector of flow line, and Θ is the angle between r and the y axis. Because dl = r dΘ, eq 3 can be written as ⎛ ∂Ω ⎞ ⎛ ∂Ω ⎞ ⎛ ∂Ω ⎞ ⎜ ⎟ = − rρW ⎜ ⎟ ⎟ + r⎜ ⎝ ∂Θ ⎠r ⎝ ∂l ⎠ψ ⎝ ∂ψ ⎠

y − G (x ) F (x ) − G (x )

Θ

i = 1, 2, 3, 4

⎛ ∂p ⎞ 1 ⎛ ∂p ⎞ ⎜ ⎟ = ⎜ ⎟ ⎝ ∂l ⎠ψ r ⎝ ∂Θ ⎠r

(5)

where Ω = p. Below are given some more equations characterizing the variation of fluid parameters in the expansion wave that are needed for the calculation of expansion rate along the flow line. If we assume that the fluid expands isentropically before reaching condensation shock then the variation of dimensionless expansion rate along the stream line will be described by the known equation10

(1)

where Φ1 = ρuyl1; Q 1 = (ρv + ρul1)y Φ2 = (p + ρu 2)yl1; Q 2 = [ρvu + (p + pu 2)l1)]y

λ 2 = λr2 + λ Θ2 = 1 +

Φ3 = ρuvyl1; Q 2 = (p + pu 2 + ρuvl 2)y

2 sin 2(z Θ̅ ) k−1

(6)

where

Φ4 = (h + p)uyl1; Q 4 = [(h + p)v + (h + p)ul 2] B1 = B2 = B4 = 0; B3 = pl1

(4)

In the case when the expanding fluid contains a condensing component but its condensation is (still or already) negligible, the equation for calculation of the static pressure along the flow tube is written as

While solving the task, it is necessary in some cases (for example, depending on the gradient of the required functions) to condense or rarefy the grid of the finite-difference scheme as dictated by particular character of the flow and introduce new variable η = f(ξ). After the change of variables, the system of gas-dynamics equations in divergence form becomes ∂Q i ∂Φi + = Bi , ∂x ∂η

(3)

v 1 = sin z Θ̅ a z * u λΘ = = cos z Θ̅ a *

λr =

(2)

The algorithm for solving the sets of eqs 1 and 2 is given in detail in the Appendix. 2.2. Modeling of Two-Phase Flow. Equations of the above mathematical model describe the expansion of homogeneous flow in each of the regions in the scheme presented in Figure 1. The model becomes much more complicated if the expansion of multicomponent fluid is accompanied by the formation of new phase. In our earlier studies,1 a number of simplifying assumptions were introduced, which generally complied with the physical principles of the process. The assumptions were as follows: (1) formation of new phase only in the subsonic expansion region, (2) 1D flow of a given expansion geometry (triangle expansion in the superand subsonic flow regions), and (3) uniform distribution of the flow velocity, density, and pressure throughout the Mach disc.

z=

⎛ k − 1 ⎞0.5 ⎜ ⎟ ⎝ k + 1⎠

Here, Θ̅ = Θ + Θ0 and a* is the critical velocity. The constant of integration Θ0 is calculated at λ1 > 1 by equation: Θ0 =

1 arcsin m(λ1 sin Θ1) − Θ1 z

Here, Θ1 is the angle between the normal and the first characteristic. 11439

DOI: 10.1021/acs.iecr.5b02180 Ind. Eng. Chem. Res. 2015, 54, 11438−11447

Article

Industrial & Engineering Chemistry Research

2.3. Particle Nucleation and Growth. As discussed above, the expansion of a supercooled solution proceeds through a metastable state followed by condensation of the solute species into a new phase. Formation of new (second) phase is the main phenomenon that dictates the necessity to model the particle nucleation in the flow and to calculate the particle size distribution. A brief review presented in our earlier work1 concerning nucleation models (see also refs 11−14) accounting for, besides coagulation, continuous particle growth by condensation, i.e., transition of a solute from the gas phase into liquid (solid) particle. The metastable state of supercooled phase is usually characterized by supersaturation degree S = p/pS, or supercooling degree ΔT = TS − T. In terms of molecular kinetics, stable supersaturated state of the fluid appears because of the fact the phase-transition time is a finite quantity. As the fluid becomes more and more supersaturated, the probability of random molecules association into a structure (a so-called nucleus) increases. The structure is considered to be stable and can act as a condensation center. Nucleation proceeds usually at the condensation shock. The dependence between the nucleus critical size and the fluid supersaturation degree is described by the Kelvin equation:15

To calculate the variation of the relative static pressure and temperature along the flow axis, the following expressions are suggested in ref 10: ε=

⎛ 1 + cos 2z Θ̅ ⎞k / k − 1 p =⎜ ⎟ p0 ⎝ k+1 ⎠

(7)

and ⎛ p ⎞k − 1/ k T 1 + cos 2z Θ̅ = ⎜⎜ ⎟⎟ = T0 p k+1 ⎝ 0⎠

(8)

Equations 6−8 are needed for the calculation of the variation of fluid expansion rate and condensation zone. The fluid expansion rate is described by differential equation for the streamline in polar coordinates:

dr r dΘ = λr λΘ

(9)

By integrating eq 9, we obtain r = r0(cos z Θ0)k + 1/ k − 1(cos z Θ̅ )−(k + 1)/ k − 1

(10)

where r0 is the radius-vector of flow line at Θ = 0. The value of longitudinal pressure gradient along the flow line strongly affects the character of supersaturated vapor expansion. The higher is the expansion rate pW, defined by equation

pW ψ = −

W dp p dl

rcr =

The stronger is the fluid supercooling before condensation shock, which reaches the maximum at the shock front. Pressure gradient in eq 5 can be expressed through other parameters of the fluid. For this purpose, we differentiate eq 7 on Θ and obtain

S=

As a result, eq 11 describing the variation of expansion rate along the flow axis converts to W dp W dp =− p dl rp dΘ λa k = * (1 + cos 2z Θ̅ )1/(k − 1) sin 2z Θ̅ rε (k − 1)(k + 1)k /(k − 1) (13)

pW = −

Substituting dimensionless rate λ (eq 6), relative static pressure ε (eq 7) and radius-vector (eq 10) into eq 13, we obtain

Φ(T , p , y) y(T , p)

⎛ ⎞ 2 2 σ * = σ∞⎜⎜1 − + 2 2 ⎟⎟ brp b rp ⎠ ⎝

0.5

r0(cos z Θ0)(k + 1)/(k − 1)

Φ(T , p , ye ) ye (Te , pe ) (16)

where Φ(T, p, y) is the dimensionless fugacity coefficient calculated by nonideal gas equation of state.18 The maximum theoretical order of supersaturation factor ranges within 106− 109. The selection of appropriate σ* value is one of the most problematic issues at calculating the nucleus critical radius and growth rate, J. Indeed, a threefold variation in σ* results in the change of J by 12 orders of magnitude. There are a number of approaches aimed at accounting for the fact that surface tension of nucleuses, comprised of tens of molecules, differs significantly from that of macro drops. In the present work, σ* = f(rp) is determined using the following expression suggested in19

kp0 dp =− (1 + cos 2z Θ̅ )1/(k − 1) sin 2z Θ̅ dΘ (k − 1)(k + 1)k /(k − 1) (12)

pW =

(15)

where σ* is the fluid surface tension, ρl is the fluid molar density (mol/m3), and S is the supersaturation factor calculated by the equation16,17

(11)

( k 2+k 1 RT0)

2σ * 1 ρl R gT ln S

Φ(k , Θ̅ )

(17)

where (1/b) ≈ (0.2−0.5) × 10−9m. Although a lot of original works and reviews devoted to detailed studies of the surface tension as a function of thermodynamic parameters are available in the literature (for example, refs 20 and 21), most authors agree that the Tolman equation19 provides the best approximation. To calculate the nucleation rate J in the fluid mass unit, we use the Frenkel−Zeldovich equation, which can be written in a form:

(14)

where ⎛ k ⎞ ⎟ sin 2z Θ̅ (cos z Θ̅ )(k + 1)/(k − 1) Φ(k , Θ̅ ) = ⎜ ⎝ k − 1⎠ ⎛ ⎞0.5 2 ⎜1 + sin 2 z Θ̅ ⎟ (1 + cos 2z Θ̅ )−1 ⎝ ⎠ k−1 11440

DOI: 10.1021/acs.iecr.5b02180 Ind. Eng. Chem. Res. 2015, 54, 11438−11447

Article

Industrial & Engineering Chemistry Research −2 ⎞ ⎛ ⎛ ln S ⎞ ⎟ J = z 0 exp⎜⎜ −⎜ ⎟ ⎟ ⎝ ⎝ z1 ⎠ ⎠

where ν is the kinematic viscosity coefficient. The coefficient “1/2” on the right side of eq 23 is introduced for correcting the double counting of the collisions between similar droplets. As a result of the collision, the larger particles grow in size, the amount of small droplets is reduced, and the total number of particles decreases as well.

(18)

where z0 =

p ρl R gT

2σ * πm 3

3. RESULTS AND DISCUSSION The demonstration of the calculated results is started from sequential complication of the modeled systems, from expansion of pure gas to expansion of complex mixtures with the formation of new phase. The calculations were performed at fixed parameters of the fluid in the reactor and the expansion medium (capillary diameter 0.1 mm, T0 = 650 K, p0 = 300 atm, T∞ = 293 K, and p∞ = 1 atm) with the exception of the studies of the effects of the temperature in the reactor (Figure 11) and fluid expansion pressure (Figure 12) on the particle size and particle size distribution. To demonstrate the correctness of the model and algorithm, Figure 2 presents the calculated boundaries of expanded

(19)

3 π ⎛⎜ σ * ⎞⎟ 1 z1 = 4 3m ⎜⎝ R gT ⎟⎠ ρl

(20)

where m is the mass of single molecule of a solute. More than ten physically meaningful equations based on various physical models are available for the calculation of the rate of condensation-driven droplets growth. The selection of a particular equation depends on the solute properties, supercooling temperature and some other parameters. In the present work, to calculate the rate of droplet growth on the nuclei by condensation mechanism, we used equation suggested in10 u

drp dl

= (8π )−1/2

RTS k+1 (R gT )1/2 p ln S k−1 ρl L h 2

(21)

where k is the adiabatic index and Lh is the evaporation (condensation) heat. The most important tasks in determining the structure of two-phase flow are to find a function that describes the particle size distribution and to study how this function forms and changes owing to particles movement in the flow. In accordance with refs 16 and 22, we assume that the particle size distribution function g(rp) obeys the normalization ∞ condition ∫ g (rp) drp = 1. Then, dg = g(rp) drp 0 Here, dg is a fraction of particles with radius in the range of (rp, rp + drp), and rp is the particle radius. Most often, the Gauss function is used to describe the particle size distribution g (rp) =

⎡ (r − r )2 ⎤ 1 p e ⎥ exp⎢ − 2 ⎢⎣ ⎥⎦ β 2π 2β

Figure 2. Boundaries of expanded nitrogen jet: 1, barrel shock; 2, jet boundary; 3, central (direct) shock wave or Mach disc; 4, reflected shock wave; 5, boundary of the main expansion zone. L, jet length; R, jet radius.

nitrogen jet. Obviously, the expansion geometry, position of barrel shocks and Max disc agree well with the real (classic) picture of an expanded gas jet described in detail in the literature, for example, that in ref 10. Figure 3 illustrates the expansion of a complex mixture (propane/methanol/phenanthrene = 90:7:3 mol %) containing

(22)

where re is the modal radius of particle and β is the radius dispersion from re. Introduction of a polydisperse system of particles enables accounting for coagulation processes in supersonic flow. The two main mechanisms of droplets coagulation are known: inertial mechanism and the mechanism of turbulent diffusion. The inertial mechanism assumes that the drops are not fully entrained by the turbulent eddies. As a result, the relative velocities acquired by drops are strongly dependent on the drop masses. The diffusion mechanism suggests complete entrainment of the drops by turbulent eddies of a scale that enables the encounters between particles. Because the particles drift chaotically in turbulent eddies, this phenomenon is similar to diffusion and can be characterized by turbulent diffusion factor. For turbulent coagulation by inertial mechanism, the equation for particle number (n) variation is written as follows:23 2

ρ U 9/4 dn 1 = − πre4 l5/4 3/4 n2 dt 2 ρg νg d

Figure 3. Boundaries of expanded fluid jet: 1, barrel shock; 2 and 5, jet boundaries; 3, direct shock wave or Mach disc; and 4, reflected shock wave. L, jet length; R, jet radius.

(23) 11441

DOI: 10.1021/acs.iecr.5b02180 Ind. Eng. Chem. Res. 2015, 54, 11438−11447

Article

Industrial & Engineering Chemistry Research a condensing substance (phenanthrene) that exists in the reactor in supercritical state. Clearly, this picture is quite different from that presented in Figure 2 and should be considered in more detail. Both barrel shock (1) and jet boundary (2) in the zone of supersonic jet expansion demonstrate inflection points where the condensation (formation of second phase) of phenanthrene begins. The condensate formation in the jet causes the decrease in the gas phase volume and, consequently, changes the jet density and speed. The boundary of reflected shock (4) increases considerably in the transient expansion zone. Moreover, the distance to the Mach disc decreases almost twice. Figure 4 demonstrates the temperature profile (1) and the fraction of the formed second phase (phenanthrene) (2) along

Figure 5. Profile of flow temperature along the jet length at flow axis: 1, system “nitrogen and water”; 2, nitrogen. T, temperature; L, length.

Figure 4. 1, Temperature along the flow axis; 2, profile of the volume fraction of the formed liquid phase along the jet length. T, temperature; L, length; Dph, phase fraction.

Figure 6. Variation of flow pressure in the vicinity of Mach disc along the jet length at flow axis: 1, system “nitrogen and water”; 2, nitrogen. p, pressure; L, length.

the jet axis through the expanded jet length. Because the assumed 2D model implies the phase distribution along the jet radius, its fraction Dph (Figure 4) is calculated as the total amount of the formed phase in the jet cross section at every point of its length. It is seen (Figure 4) that as the supersonic jet expands its temperature monotonically decreases, and in the moment of jumplike formation, the new phase remains almost unchanged. This fact is due to the assumption of the same temperature for both phases. Condensation of the liquid phase from the flow proceeds until the Mach disc where the jet temperature jumps up and the liquid phase partially evaporates; its fraction decreases. The decrease in the jet temperature in the subsonic expansion zone causes again weak condensation of phenanthrene from the gas phase to the liquid one. The behavior of the main parameters of the expended jet along its length should depend on its qualitative and quantitative initial composition, i.e., on the composition of supercritical fluid in the reactor (Figure 5). Indeed, the addition of 10 mol % of water into nitrogen leads to considerable shortening of the distance to the Mach disc and its size reduction. Variation of initial composition (Figure 6; 1, nitrogen−water = 90:10 (mol %); 2, pure nitrogen) changes the pressure drop before/after the Mach disc and the distance of the direct shock wave formation from the capillary cut but affects insignificantly the pressure profile in the flow. As mentioned above, the RESS process is intended for the production of nano- and microparticles of the material dissolved in supercritical fluid in the reactor. The final step of RESS simulation is a complex task that includes the modeling of the solute condensation in the expanded fluid, selection of models for nucleation and particle growth, and determination of the particle size distribution through the jet cross section.

A 2D model of the flow expansion process makes it possible to study the characteristics of the flow in the cross section of the jet. Figure 7 illustrates how the size of phenanthrene

Figure 7. Particle size distribution throughout the jet radius: 1, size of nuclei; 2, particle radius. r, particle radius; R, jet radius.

particle (2) varies along the radius of a jet cross section where the formation of new phase proceeds. Here, the size of nuclei of the new phase remains constant, although under some other conditions the nuclei size will vary. Line 1 in Figure 7 shows that the radius of a nucleus is about 1.4 nm. Although the particle size in the jet flow seems to increase insignificantly, the total particle number strongly decreases considerably (by ∼5 orders of magnitude) in the zone of phenanthrene condensation (Figure 8). Further decrease in the 11442

DOI: 10.1021/acs.iecr.5b02180 Ind. Eng. Chem. Res. 2015, 54, 11438−11447

Article

Industrial & Engineering Chemistry Research particle number proceeds by other mechanism (collisions and growth in size).

Figure 10. Dependence of average particle size on the distance from the capillary cut. r, particle radius; L, jet length.

Figure 8. Variation of total number of particles along the expanded jet length. N, particle number; L, jet length.

Figure 9 presents the particle size distributions in two cross sections of expanded jet at various distances from the capillary

Figure 11. Effect of temperature in the reactor on the particle size distribution in the jet at a 35 mm distance from the capillary cut. Composition of the fluid (%): propane/methanol/phenanthrene = 90:7:3. External pressure, 5 atm. Temperature in the reactor: 1, 500 K; 2, 550 K; 3, 600 K. r, particle radius; R, jet radius.

Figure 9. Particle size distribution in two cross sections of the expanded jet. Initial composition of the fluid (mol %): propane/ methanol/ phenanthrene = 90:7:3. r, particle radius; R, jet radius. 1, distance of 25 mm from the capillary cut; 2, distance of 15 mm from the capillary cut.

cut. The average particle size is ∼100 nm at a 15 mm distance (curve 2) and ∼350 nm at a distance of 25 mm (curve 1). The particle size increases mainly because of coagulation. Figure 10 illustrates how the average particle size varies along the jet length. Points 1 and 2 designate particle size distributions corresponding to the curves in Figure 9. Points 1 and 2 correspond, respectively, to 15 and 25 mm distances from the capillary cut. The calculations showed that the temperature and composition of the fluid in the reactor, external pressure, and some other parameters affect the size and size distribution of the formed phenanthrene particles (Figures 11−13). It is seen (Figure 11) that as the fluid temperature in the reactor increased from 500 to 600 K, the jet expansion radius at the same distance from the capillary cut (35 mm) increased insignificantly whereas the maximum particle size decreased by half. External pressure produces an even stronger effect on the size and particle size distribution (Figure 12). As the external pressure was elevated from 2 to 10 atm, the maximum particle size increased more than sixfold. That is, the size of the formed

Figure 12. Effect of external pressure on the particle size distribution in the jet at a 35 mm distance from the capillary cut. Composition of the fluid (%): propane/methanol/phenanthrene = 90:7:3. Temperature in the reactor, 500 K. External pressure: 1, 10 atm; 2, 5 atm; and 3, 2 atm. r, particle radius; R, jet radius.

particles can be easily controlled by varying the external pressure. 11443

DOI: 10.1021/acs.iecr.5b02180 Ind. Eng. Chem. Res. 2015, 54, 11438−11447

Article

Industrial & Engineering Chemistry Research

u=

v=

k k−1

k Φ2 + ⎡⎣ k − 1 Φ22 +

k+1 (Φ32 k−1

2 − 2Φ1 Φ4 )⎤⎦

k+1 Φ k−1 1

Φ3 Φ − Φ1y Φ ,p= 2 Φ2 , ρ = 1 y uy Φ1

(1.1)

Expressions for l1 and l2 depend on the type of grid refinement η = f(ξ). For grid refinement at upper boundary G(x)7 (a − 1) F(x) − G(x) ; l2 ln a (a − 1)η + 1 ln[(a − 1)η + 1] (G′(x) − F ′(x)) − G′(x) = ln a

l1 =

Figure 13. Effect of fluid composition in the reactor on the particle size distribution in the jet at a 35 mm distance from the capillary cut. Composition of the fluid (%): propane/methanol/phenanthrene: 1, 88:9:3; 2, 90:7:3; and 3, 92:5:3. Temperature in the reactor, 500 K. r, particle radius; R, jet radius.

where η=

The dash indicates a derived function. For grid refinement at lower boundary F(x)

As the fraction of phenanthrene in the fluid was decreased by ∼4%, the maximum particle size reduced more than 1.5 times at the same jet dimensions (Figure 13). It should be noted at the end of this section that we are faced with serious difficulties with regard to the experimental verification of the calculated results, but we hold out hope to overcome them.

η=

ln[(a − 1)ξ + 1 aη − 1 ,ξ= ln a a−1

For a uniformly spaced grid: η = ξ , l1 = F(x) − G(x), l 2 = η(G′(x) − F ′(x)) − G′(x)

System (1) is approximated by difference equations on the grid defined by scheme: xn = nτ ; ηm = mζ

4. CONCLUSIONS

where n = 0, 1, ...; m = 0, 1, ..., Nχ; Nχ is the number of points on η, determined experimentally; and τ, ζ are the grid spacing. To determine the parameters of the internal nodes of the grid, the second-order-accurate explicit-difference scheme of predictor−corrector type is used. At first, using the known values on n layer in nodes (xn, ηm−1), (xn, ηm), and (xn, ηm+1) and the Lax scheme,7,8 the values on intermediate layer n + σ, 0 < σ ≤ 1, in nodes (xn+σ, ηm−1/2) and (xn+σ, ηm+1/2) are determined. (Here, σ is weight determined experimentally.) Then, using the “cross” scheme and the values in nodes (xn, ηm), (xn+σ, ηm−1/2), and (xn+σ, ηm+1/2), the values of the required functions in nodes (xn+1, ηm) on (n + 1) layer are calculated. Under this scheme, the difference equations of system (1) are as follows:

Mathematic modeling of the rapid expansion of supercritical solution and formation of particles by condensation was performed. Mathematical model and calculation algorithms were developed and applied for the calculations. A 2D model of expanded supercritical fluid jet was written for super- and subsonic zones; unique algorithms for solving the model equations were developed. As a result, the calculations allowed determination of jet boundaries, shape and position of normal and barrel shock waves, and condensation front. The calculation of initial nuclei size distribution in the jet cross section enabled determination of the particle size distribution at the end of the jet expansion. Parameter sensitivity analysis of the model helped to identify which parameters are of key importance for the production of particles with the required size and properties. Calculation of supercritical fluid expansion (nitrogen, nitrogen−water, propane−methanol−phenanthrene) and modeling of particle nucleation and growth provided efficient (more reliable) optimization of the process parameters.

■

ln[(a − 1)η + 1] aξ − 1 ;ξ= ln a a−1

for intermediate layer n + σ: (Φi)nm++σ1/2 =

1 [(Φi)nm + 1 + (Φi)mn ] 2 τσ − [(Q i)nm + 1 + (Q i)nm ] n τσ + [(Bi )nm + 1 + (Bi )nm ] 2

(2.1)

for layer (n + 1):

APPENDIX

We present a method for calculation of stationary flow of supersonic nonisobaric jet with the central shock wave. The equations of nonisobaric nonisothermal jet expansion (1) are solved by using finite-difference schemes. Given the known functions Φ1, Φ2, Φ3, Φ4 (2), variables u, v, p, and ρ are calculated by the following equations:7

(Φi)nm+ 1 = (Φi)nm − +

τ [(Q i)nm++σ1/2 + (Q i)nm+−σ1/2 ] ζ

ττ [(Bi )nm++σ1/2 + (Bi )nm+−σ1/2 ] 2n

(3.1)

To determine the parameters at the boundaries G(x) and F(x), differential equations of the initial system with respective 11444

DOI: 10.1021/acs.iecr.5b02180 Ind. Eng. Chem. Res. 2015, 54, 11438−11447

Article

Industrial & Engineering Chemistry Research

All values in node (xn+σ, η1/2) are calculated by eqs 2.1 at m = 0. By symmetry, we obtain

boundary conditions are used. In this case, the differential equations are approximated by a rectangle difference scheme with weights σ1 and 1 − σ1 for η derivatives on layers n + 1 and n. Consider in detail the algorithm for finding solutions in the boundary nodes for submerged jet, i.e., at the outer edge of the jet and the axis of symmetry. Let y = F(x) be a function of jet boundary. The pressure and entropy are kept constant along it. On the basis of these conditions and Bernoulli equations, four boundary expressions can be written: p = const, σ = const, W = const, and F′(x) = v/u. Here, W (= u + v) is the velocity vector. Then, we add the equation of system (1) at i = 2, which in differential form looks like (Φ2)nM+ 1

−

+σ n+σ σ n+σ ρ−n1/2 = ρ1/2 = −(l 2)1/2 ; (l 2)−n +1/2

−

(Φ2)nN+χ −1 1

+

2τ + σ1(Q 2)nNχ − 1 ζ

2τ n+σ (2vρ + uρl 2)1/2 ; ζ

(Φ̅ 2)n0+ 1 = (Φ̅ 2)n0 +

2τ n+σ (2vuρ + (p + u 2ρ)l 2)1/2 ; ζ

(Φ̅4 )n0+ 1 = (Φ̅4 )n0 +

ρW 2 ⎞ 2τ ⎛ k p+ ⎟ ⎜ 2 ⎠1/2 ζ ⎝k − 1

n+σ (2v + ul 2)1/2

(4.1)

k+1 k A1 = d1; A 2 = d 2 ; A 2 = d4 2(k − 1) (k − 1)

⎧ 2τ ⎨(p0 + ρ0 W02 cos Θ) F(x) l1 + σ1[ρ0 W02 cos Θ sin Θ ⎩ ζ ⎫n + 1 + (p0 + ρ0 W02 cos Θ)l 2] F(x)⎬ ⎭N

Then, the values of all boundary parameters on layer n + 1 are as follows: u 0n + 1 =

χ

= A2

d2 +

p0n + 1 =

The value of the boundary coordinate on layer n + 1 is

A 2 − A1u 0n + 1 (l1)n0+ 1

n

(x) = F (x) + [F ′(x)] τ

∂Q̅ i ∂Φ̅ i + = 0, i = 1, 2, 4 ∂x ∂η

; ρ0n + 1 =

A1 u(l1)n0+ 1

(10.1)

k=n

(Φ̅ i0)m =

∑ k =−n

where

ak (Φ̅ i)m − k

(11.1)

where ak is a member of descending geometric progression n satisfying the condition: ∑kk==−n (ak = 1), ak = a−k. The first member of the progression is determined by equation q0 = (q − 1)/(2qk+1 − q − 1). The geometric ratio q and the number of averaging points 2n + 1 are fixed. The smoothing procedure is performed for five nodes in internal points and three nodes in boundary points. No smoothing in the boundary points was applied. Various approximation models and methods suggested in literature for determining the position of direct shock wave9 produce significantly different results. Moreover, they are inappropriate for calculating the shape and size of direct shock wave; the parameters of the shock wave are independent of the character and structure of the upstream jet. The calculation of stationary flow of supersonic nonisobaric jet with central shock wave reduces to solving mixed elliptic− hyperbolic boundary value problem. Let us use the integral method7 in the approximation of axially symmetric flow of ideal gas in divergence form:

Φ̅ 1 = uρl1; Φ̅ 2 = (p + u 2ρ)l1; Φ̅4 ⎛ k ρW 2 ⎞ p+ =⎜ ⎟ul1; 2 ⎠ ⎝k − 1 Q̅ 1 = 2vρ + uρl 2 ; Q̅ 2 = 2vuρ + (p + u 2ρ)l 2 ; ⎛ k ⎛ k ρW 2 ⎞ ρW 2 ⎞ Q̅ 4 = 2⎜ p+ p+ ⎟ul 2 ⎟v + ⎜ 2 ⎠ 2 ⎠ ⎝k − 1 ⎝k − 1 (6.1)

All values (eqs 6.1) in node (xn+1, η0) are calculated by crossscheme using equation:

i = 1, 2, 4

; v0n + 1 = 0;

To extinguish oscillations at solution discontinuities, the smoothing procedure is applied. It is used for smoothing complexes Φi rather then the required functions. After smoothing, the final values of gas-dynamic parameters on layer n + 1 are determined using eq 1.1. The smoothed values of the complexes are calculated by equations:

Equation 5.1 allows a determination of the angle Θ of inclination of the velocity vector to the jet axis. Assume y = G(x) = 0 to be the symmetry axis. There is a single boundary condition v = 0 on this axis. After adding three equations of system (1), obtain

(Φ̅ i)n0+ 1 = (Φ̅ i)n0 −

d 22 − 4d1d4 2d1

(5.1)

n

(9.1)

Let us denote

On the basis of expression (2) and boundary conditions p = p0 = const, ρ = ρ0 = const, W = W0 = const, u = W0 cos Θ, and v = W0 sin Θ, we write eq 4.1 in expanded form as

F

(Φ̅ 1)n0+ 1 = (Φ̅ 0)n0 +

n+σ

(Φ2)nNχ − 1

2τ (1 − σ1)[(Q 2)nNχ + (Q 2)nNχ − 1] = A 2 ζ

n+1

(8.1)

In view of eqs 8.1, eqs 7.1 are written as

2τ + σ1(Q 2)nN+χ 1 ζ

(Φ2)nNχ

=

+σ n+σ +σ n+σ +σ n+σ = u1/2 = −v1/2 = p1/2 u−n 1/2 ; v−n 1/2 ; p−n1/2 ;

τ n+σ σ − (Q̅ i)−n +1/2 [(Q̅ i)1/2 ] = Ai , ζ (7.1) 11445

DOI: 10.1021/acs.iecr.5b02180 Ind. Eng. Chem. Res. 2015, 54, 11438−11447

Article

Industrial & Engineering Chemistry Research

index 1 are determined by calculating the triple point if its position on the barrel shock and the parameters of incoming flow are known. By symmetry condition, the central shock wave in the intersection with the stream axis is the direct one. If the parameters of the incoming flow are known, then the initial parameters with index 0 will be clear at the preset coordinate of direct shock wave position. Integrating eq 15.1, obtain

∂ρuy ∂ρvy + =0 ∂x ∂y ∂(yuvρ) ∂ ⎡ ⎛ 2 k − 1 ⎞⎤ =0 p⎟⎥ + ⎢y⎜⎝ρu + ⎠ ⎣ ⎦ ∂x ∂y 2k p W2 + = 1 ρ

(12.1)

xD − xB =

Here, u and v are the projection of velocity vector in cylindrical coordinates related to the maximum velocity; W is the velocity module; p and ρ are pressure and density related to the pressure and density, respectively, upstream of the shock. Divide the flow zone between the symmetry axis and tangential discontinuity line into Nn equidistant bands. Integrate the first and second equations of the system (eqs 12.1) with respect to y from the symmetry axis to each band boundary. As a result, we obtain an approximating system of ordinary differential equations. Write equation for Nn = 1. After selecting the desired functions as a pressure p and Q, the angle between the velocity vector and the positive flow axis, we obtain from the integrated equations

where ωB is the angle of direct shock inclination determined by the triple point. So far, the position of the triple point on the barrel shock yB and quantity (dΘ1/dx)x=xB in eqs 13.1 and 14.1 remain unknown. Equations 13.1 and 14.1 have the singular points. Expression in the second brace in the right member of eq 14.1 goes to zero in the singular point where W0 becomes equal to sonic speed. Expression in the second brace in the right member of eq 13.1 also becomes zero in the singular point where W1 attains the sonic speed. For continuous flow to exist in the vicinity of the sonic line (Mach disc), it is necessary that the expressions in the first braces in the right members of eqs 13.1 and 14.1 become zero as well. To meet these conditions, two free parameters are available at

⎧ ⎛ W ⎞ dΘ = ⎨2y1ρ1W12 sin Θ1 ⎜2 cos Θ1 − 0 ⎟ 1 dx W1 ⎠ dx ⎝ ⎩ ⎡k − 1 ⎛W ⎞⎤⎫ + 2tg Θ1⎢ (p1 − p0 ) + 2ρ1W12 cos Θ1 ⎜ 0 − cos Θ1⎟⎥⎬ ⎢⎣ 2k ⎝ W1 ⎠⎥⎦⎭

dp1

1 y ctgωB 2 B

⎪

⎪

⎪

x = x B/yB and

⎪

⎧ ⎡ ρ W2 ⎛ ⎪ W⎞ ⎨2y1⎢ 1 1 cos Θ1 ⎜cos Θ1 − 0 ⎟ ×⎪ ⎢ kp W1 ⎠ ⎝ ⎩ ⎣ 1 +

⎫−1 W ⎞⎤⎪ k − 1⎛ ⎬ ⎜1 − 2 cos2 Θ1 + cos Θ1 0 ⎟⎥⎪ 2k ⎝ W1 ⎠⎥⎦⎭

According to estimations, the variation of parameters in the initial cross section of subsonic jet downstream the direct shock is insignificant. Therefore, 1D approximation can be used to describe the jet expansion downstream the shock wave, which is assumed to be a zero-order approximation of integral method. 1D approximation of gas dynamics equations provides an equation for the pressure variation downstream the direct shock wave:

(13.1)

⎧ ⎛k − 1 W ρ W0W1 ⎞ dp1 ⎪ 0 ⎟ =⎨ − 1 2y cos Θ1 ⎜⎜ ⎪ 1 kp1 ⎟⎠ dx dx ⎝ 2k W1 ⎩

dp0

−1 dy ⎡ ⎛ 1 dp k − 1 1 ⎞⎤ ⎥ = −2 1 ⎢y1⎜ − ⎟ dx dx ⎢⎣ ⎝ kp 2k ρW 2 ⎠⎥⎦

⎫ ⎪ dΘ1 ⎬ + 2y1ρ1W1W0 sin Θ1 − 2tg Θ1(ρ0 W02 + 2ρ1W1W0 cos Θ1)⎪ dx ⎭ ⎧ ⎛ ρ W2 ⎫−1 k − 1 ⎞⎟⎪ 0 0 ⎜ ⎨y1⎜ ⎬ ×⎪ − 2k ⎟⎠⎪ ⎩ ⎝ kp0 ⎭

(14.1)

Indexes 0 and 1 are attributed to parameters on the zero (symmetry axis) and first (Nn = 1) lines. Add the equation for the flow line passing through triple point (Figure 1) and separating the flows downstream of reflected and direct shock waves. dx

= tg Θ1

(15.1)

■

One more equation is obtained from the boundary conditions for this flow line. Assume a completely supersonic jet at y > y1. Using the method of characteristics, write the consistency condition along the second family characteristic relating p and Θ in the form dΘ −

sin α cos α sin Θ sin α dx = 0 dp − kp cos(Θ − α) y

(17.1)

where y1 = y1(x) is the equation for the line of tangential discontinuity dividing the flow downstream the reflected and direct shock waves. Quantities p, ρ, and W in eq 17.1 are dimensionless (similar to those in eq 12.1). Note that the static pressure (p∞) is the same on the both sides of the tangential discontinuity line. The bracketed expression in the right member of eq 17.1 becomes zero in the singular point with M = 1. Consequently, in order for continuous flow to exist in the vicinity of the sonic line, it is necessary to satisfy the condition of M = 1. On the basis of these conditions, the position of the triple point on the direct shock wave can be determined by iterating.

⎪

dy1

dΘ1 dx

AUTHOR INFORMATION

Corresponding Author

*[email protected]. Notes

The authors declare no competing financial interest.

■

(16.1)

Solving the systems of eqs 13.1−16.1, determine the required functions p0(x), p1(x), Θ1(x) and y1(x). Initial parameters with 11446

NOMENCLATURE F, G = dimensionless functions l = coordinate, m DOI: 10.1021/acs.iecr.5b02180 Ind. Eng. Chem. Res. 2015, 54, 11438−11447

Article

Industrial & Engineering Chemistry Research

(2) Li, J.; Matos, H. A.; de Azevedo, E. G. Two-phase homogeneous model for particle formation from gas-saturated solution processes. J. Supercrit. Fluids 2004, 32, 275−286. (3) Turk, M.; Lietzow, R. Formation and stabilization of submicron particles via rapid expansion processes. J. Supercrit. Fluids 2008, 45, 346−355. (4) Tom, J. W.; Debenedetti, P. G. Particle formation with supercritical fluids  a review. J. Aerosol Sci. 1991, 22, 555−584. (5) Kwauk, X.; Debenedetti, P. G. Mathematical modelling of aerosol formation by rapid expansion of supercritical solutions in a converging nozzle. J. Aerosol Sci. 1993, 24, 445−469. (6) Dulov, V. G.; Luk’yanov, G. A. Gas-dynamics of outflow processes. Nauka: Novosibirsk, Russia, 1984 (in Russian). (7) Avduevskii, V. S.; Ashratov, E. A.; Ivanov A. V.; Pirumov, U. G. Supersonic non-isobaric gas jets. Mashinostroenie: Moscow, 1985 (in Russian). (8) Rouch, P. Computational fluid dynamics; Mir: Moscow, 1976 (in Russian). (9) Bobyshev, S. V.; Dobroserdov, I. L. The Principles of Construction of Algorithms for the Calculation of Non-Isobaric Turbulent Jets. Leningrad, Tutorial Leningradskiy Mechanical Institute: Leningrad, 1988 (in Russian). (10) Saltanov, G. A. In Supersonic Two-Phase Flows; Deich, M.E., Stepanchuk, V.F., Eds.; Vysshya Shkola: Minsk, 1972 (in Russian). (11) Erriguible, A.; Marias, F.; Cansell, F.; Aymonier, C. Monodisperse model to predict the growth of inorganic nanostructured particles in supercritical fluids through a coalescence and aggregation mechanism. J. Supercrit. Fluids 2009, 48, 79−84. (12) Mountain, R. D.; Mulholland, G. W.; Baum, H. Simulation of aerosol agglomeration in the free molecular and continuum flow regimes. J. Colloid Interface Sci. 1986, 114, 67. (13) Kaplan, C. R.; Gentry, J. W. Agglomeration of chain-like combustion aerosols due to brownain motion. Aerosol Sci. Technol. 1988, 8, 11−28. (14) Mansoori, G.A.. Principles of Nanotechnology; World Science: Hackensack, NJ, 2005. (15) Sutyagin, A. G. Spontaneous condensation of steam and the formation of condensation aerosols. Russ. Chem. Rev. 1969, 38, 79. (in Rusian). (16) Fuchs, N. A. The Mechanics of Aerosols; Academy of Sciences of the USSR: Moscow, 1964 (in Russian). (17) Weber, M.; Russell, L. M.; Debenedetti, P. G. Mathematical modeling of nucleation growth formed by the rapid expansion of a supercritical solution under subsonic conditions. J. Supercrit. Fluids 2002, 23, 65−80. (18) Anikeev, V. I.; Ermakova, A. The thermodynamic characteristics of formation of nanoparticles from supercritical solvents. Russ. J. Phys. Chem. 2007, 81, 2024−2029. (19) Tolman, R. C. The defect of droplet size on surfaces tension. J. Chem. Phys. 1949, 17, 333−337. (20) Zhukhovitskiy, D. I. The surface tension of the vapor-liquid boundary with a finite curvature. Kolloid. Zh. 2003, 65, 480−494 (in Russian). (21) Turk, M. Influence of thermodynamic behaviour and solute properties on homogeneous nucleation in supercritical solutions. J. Supercrit. Fluids 2000, 18, 169−184. (22) Fuchs, N. A. Evaporation and growth of droplets in a gaseous medium. [publisher-name]: Moscow, 1958 (in Russian). (23) Sinaiskiy, E. G.; Lapta, E.Ya.; Zaitsev, Yu.V. Separation of multiphase multi-component systems. Nedra-Business Center: Moscow, 2002 (in Russian).

l1, l2 = functions of special type x, y, ξ, η = dimensionless coordinates p = pressure, Pa ρ = density, kg/m3 u, v = longitudinal and transverse velocity, m/s U = dimensionless flow rate T = temperature, K h = enthalpy, J/kg d = diameter of particle, nm k = adiabatic exponent ak = member of descending geometric progression q = progression ratio Φi, Qi, Ii = functions a = grid condensation parameter W = velocity modulus, m/s Mχ = number of points on η N = number of points on x τ, ζ = grid spacing σ = weight factor σ1 = boundary weight factor ωB = angle of central shock inclination rs̅ = dimensionless jet radius S = oversaturation factor, dimensionless value ψ = flow line α = expansion angle, degrees Θ = angle between velocity vector and flow axis σ* = surface tension, J/m2 ε0 = specific dissipation energy ν = kinematic viscosity factor, m2 s M = Mach number, dimensionless Φ = fugitivity coefficient n = number of particles, 1/m3 rp = particle radius, nm re = modal radius of particles, nm β = radius dispersion g(r) = particle size distribution J = nucleation rate Mw = molecular weight, g/mol Rg = gas constant, J/(° mol) Lh = evaporation heat Indexes

* = parameter of a particle at critical radius up = upstream Mach disc after = downstream Mach disc 0 = conditions in the reactor (initial conditions) ∞ = external conditions e = oversaturation conditions fr = freezing (temperature) in = conditions at the capillary inlet out = conditions at the capillary outlet liq, l = liquid g = gas cr = critical values p = particle a = capillary end

■

REFERENCES

(1) Anikeev, V. I.; Stepanov, D. A.; Ermakova, A. Modeling and

calculation of the process of rapid expansion of supercritical fluid with the formation of nanoparticles. Theor. Found. Chem. Eng. 2011, 45, 141−155. 11447

DOI: 10.1021/acs.iecr.5b02180 Ind. Eng. Chem. Res. 2015, 54, 11438−11447