Transport of nanoparticles and kinetics in packed beds: A numerical approach with lattice Boltzmann simulations and particle tracking

International Journal of Heat and Mass Transfer 72 (2014) 319–328

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Transport of nanoparticles and kinetics in packed beds: A numerical approach with lattice Boltzmann simulations and particle tracking Ngoc H. Pham a, Daniel P. Swatske a, Jeffrey H. Harwell a, Bor-Jier Shiau b, Dimitrios V. Papavassiliou a,⇑ a b

School of Chemical, Biological, and Materials Engineering, University of Oklahoma, Norman, OK 73019-1004, United States Mewbourne School of Petroleum and Geological Engineering, University of Oklahoma, Norman, OK 73019-1004, United States

a r t i c l e

i n f o

Article history: Received 24 May 2013 Received in revised form 10 December 2013 Accepted 28 December 2013 Available online 1 February 2014 Keywords: Porous media Nanoparticle transport Particle tracking

a b s t r a c t The lattice Boltzmann method is employed in conjunction with a Lagrangian particle tracking algorithm to investigate the fate and transport of nanoparticles as they propagate in porous columns that are packed with spherical particles. In this approach, physical phenomena that result in particle retention and remobilization are represented by a probability for adsorption and desorption, respectively. The method is validated with experiments where polymer-stabilized purified multi-walled carbon nanotubes propagate in a column packed with inert glass beads. Comparison of simulation results to the conventional filtration equation leads to the correlation of the simulation input parameters to macroscopically observed parameters, such as adsorption and desorption rate constants. It is found that adsorption and desorption nominal rates do not affect the diffusivity of the nanoparticles significantly (the change in particle diffusivity is less than 4%) and that particles with smaller size, where the Brownian motion is dominant, are retained more than larger particles. The difference in particle recovery for small versus large particles is of up to 36% when the probability for adsorption is as high as 0.01. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Engineered nano-sized particles, or nanoparticles (NPs), have recently attracted much attention due to their potential for applications in various sectors of engineering. This class of materials possesses unusual and unique physical, chemical, and mechanical properties that make them ideal for many significant applications [1–6]. Metal-based NPs, for example, have found applications in improving biosensors, in cancer therapy, cell labeling, and targeted drug delivery [7,8]. Likewise, carbon-based NPs have been intensively utilized in high-flux membranes, composite fillers, pollution prevention, and energy storage [9–13]. At the same time, concerns have been raised that the widespread use of NPs would lead to releases of significant amounts of them into the environment. Due to their tiny size, NPs released uncontrollably can migrate through soils and sentiments and can penetrate deeply in the subsurface, reaching aquifers and groundwater resources [14–16]. In other situations, however, the dispersion of nanoparticles in the pore space of a porous medium might be desirable, when one wants to track the path of injected hydraulic fracturing fluids or to potentially affect enhanced oil recovery in petroleum reservoirs [17]. To protect the aquifers and mitigate environmental risks, better prediction of NP transport in soil, or porous media in general, is ⇑ Corresponding author. Tel.: +1 (405) 325 5811; fax: +1 (405) 325 5813. E-mail address: [email protected] (D.V. Papavassiliou). 0017-9310/$ - see front matter Ó 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.12.075

highly desirable. Flow-through experiments with granular material packed columns are usually conducted to investigate the effect of various factors on the transport of NPs, as well as their interaction with the pore surfaces. It has been shown that propagation of NPs through the pores can be hindered by size exclusion (or straining), but most importantly by pore surface adsorption effects [18–21]. The former effect is simply the physical retention of particles by the porous medium due to the larger size of dispersed particles relative to pore throats [22]. In the case of nanoparticles, straining might occur due to agglomeration of the nanoparticles [18]. It has also been found that breakthrough or retention of a particle in the column depends on several parameters, such as the pH of the NP suspension, the ionic strength of the aqueous phase, the surface charge of the medium, the surface roughness, the pore fluid velocity, and the size of the particle [19,23–33]. In addition to experiments conducted in physical models, particle breakthrough behavior is predicted by utilizing the one-dimensional convectiondiffusion equation [34,35]. The terms of this mathematical equation indicate that the rate of change of particle concentration in the effluent is controlled by molecular diffusion (due to Brownian motion), convective diffusion (due to an external driving force), and the saturation level of adsorbed particles onto the pore surface area. Among these, the accurate modeling of the physics related to the pore surface saturation has attracted significant attention. As a result, the filtration equation has evolved to account for surface saturation by introducing a nonlinear term representing the

320

N.H. Pham et al. / International Journal of Heat and Mass Transfer 72 (2014) 319–328

Nomenclature c C C0 Cn dc dp D Da D0 ~ e f feq ff k ka kB kd k0 K L n n0 nx ny nz N N0 N(0, r) pa pd

lattice speed particle concentration in each pore volume interval initial particle concentration Cunningham correction factor collector diameter particle diameter effective diffusivity Damköhler number nominal molecular diffusivity microscopic velocity particle distribution function particle equilibrium distribution function forcing factor permeability of the porous medium adsorption rate constant Boltzmann constant desorption rate constant nominal adsorption rate constant equilibrium coefficient particle travel distance number of allowable directions ratio of the total amount of particles injected to the volumetric flow rate number of nodes in x direction number of nodes in y direction number of nodes in z direction number of particles captured number of particles injected normal distribution with zero mean and standard deviation r adsorption probability desorption probability

so-called blocking effect [36–38]. In addition, particles lost due to size exclusion can also be accounted for by adding a straining term into the convective diffusion equation [39]. In this work, we present a numerical approach to the transport and kinetics of NPs migrating through the pore space. We adopt the lattice Boltzmann method (LBM) to simulate the flow of an incompressible Newtonian fluid creeping through an infinite array of spheres packed in a random manner. Sequentially, a set of passive particles, representing the nanoparticles, are injected into the flow field and the trajectory of each individual particle is monitored in space and time using a particle tracking algorithm. The instantaneous position and velocity of each individual particle is calculated and recorded in each time step. To simulate the adsorption and desorption kinetics, particles are assigned a probability value related to their retention rate when they collide with the wall, and when attached a different probability related to their remobilization and detachment rate from the wall. To validate our simulation approach, laboratory experiments were conducted in a column packed with inert glass beads, and the simulation results were further compared with theoretical predictions. The nanoparticles were purified multi-walled carbon nanotubes (PMWCNTs) that were stabilized in suspension with the use of polymers. 2. Numerical and experimental methodology 2.1. Lattice Boltzmann method (LBM) We have employed LBM to simulate the flow through porous media composed of packed spheres. A custom-written code has been developed, details and validation of which can be found

Pe r S Sc t T u U US ~t U w ~ x x

Peclet number particle radius adsorbed particle concentration Schmidt number time absolute temperature pore velocity macroscopic velocity in feq superficial velocity particle velocity at time t lattice specific weighing factor position location in the flow direction

Greek symbols Dt time interval of each time step Dx lattice constant x D~ travel distance by diffusion g single collector efficiency l fluid dynamic viscosity m fluid kinematic viscosity q fluid density qb bulk density of the column r standard deviation s relaxation time / porosity of the column X collision operator Subscript i lattice direction index

elsewhere [40,41]. In LBM, the geometry of the porous medium is discretized into lattice points (solid walls are described by logical ‘TRUE’ value and mesh nodes within the void pores are described as logical ‘FALSE’). The fluid flow is simulated by calculating the collisions and interactions between fluid particles that move on a rectangular lattice in phase space by solving the discrete Boltzmann equation for a single particle probability distribution function as a function of space and time as follows:

fi ð~ x þ~ ei Dt; t þ DtÞ ¼ fi ð~ x; tÞ þ X ð~ x; tÞ  ffi |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflffliffl{zfflfflfflffl} |{z} STREAMING

COLLISION

ð1Þ

FORCING

x is the position, t is the where f is the particle distribution function, ~ time, Dt is the time step, ~ e is the microscopic velocity, X is the collision operator, ff is the forcing factor (pressure drop over length) and the subscript ‘i’ is a lattice direction index (there are a total of 15 such directions in our model, which is commonly called D3Q15 in the LBM literature [42]). The terms on the right hand side of Eq. (1) constitute the three steps of the LBM algorithm, namely the streaming, collision and forcing steps. We use the single-relaxation time approximation of the collision term given by Bhatnagar, Gross and Krook (BGK) [43]. In the BGK model, the collision operator is approximated as

x; tÞ ¼  Xi ð~

 1 fi  fieq

ð2Þ

s

where the particle equilibrium distribution function, fieq , is given by

"

fieq ð~ xÞ

~ ~2 3U ~2 ~ ð~ ei  U ei  UÞ  ¼ wi q 1 þ 3 2 þ 9 c 2 c2 c4

# ð3Þ

321

N.H. Pham et al. / International Journal of Heat and Mass Transfer 72 (2014) 319–328

where c = Dx/Dt is the lattice speed, Dx is the lattice constant, q is the density, w is a lattice specific weighing factor and U is the macroscopic velocity. The time s appearing in Eq. (2) is the time scale with which the local particle distribution function relaxes to equilibrium, and is often referred to as the ‘‘relaxation time’’. It is related to the kinematic viscosity, m of the fluid as follows:

m¼

1 3



s

1 2

 ð4Þ

The final step in the LBM algorithm is to calculate the macroscopic properties of the fluid such as density, q, and velocity, U, at any instant, from the conservation equations of mass and momentum given by

q¼

n X

2.2. Lagrangian particle tracking (LPT)

fi

ð5Þ

i¼0

~¼ qU

Simulation boxes were all meshed such that nx  ny  nz = 101  101  101. The choice of this mesh size is a balance between grid independence [46] and reasonably fast convergence. Table 1 is a summary of simulation configurations and conditions implemented in the simulations presented in this work. Note that the reported permeabilities were directly calculated from simulations by employing Darcy’s law for flow in low Reynolds number regime. Briefly, a series of simulations with various pressure drops were performed for each sphere packing type, yielding different values of superficial velocity. A plot of superficial velocity versus pressure drop was then constructed and the permeability was obtained from the slope of the line.

n X

fi~ ei

ð6Þ

i¼0

where n is the number of directions that the fluid particles are allowed to move, in addition to the zero position, which is the rest position of a fluid particle. The simulation mesh consists of nx, ny, and nz nodes in the x, y, and z directions, respectively. Among these, fluid nodes are those within the flow field (i.e., within the empty pore space, given the logical value ‘‘FALSE’’) and wall nodes are those that make up the rigid wall (those given the logical value ‘‘TRUE’’). The velocity field generated by solving the above equations (Eqs. (5) and (6)) is equivalent to solving the Navier–Stokes equations for single-phase or multi-phase flows through the pore spaces with 2nd order accuracy [44]. Periodic boundary conditions were applied in all three directions (x, y, z). The no-slip boundary condition was applied at the wall faces using the bounce-back technique. In order to take advantage of the LBM parallelizability, the domain was decomposed using message passing interface. Geometries for the simulations in this study consist of typical ideal sphere packing arrays [body centered cubic (BCC), face centered cubic (FCC), and simple cubic (SC)] and arrays of spheres packed in a random manner. Fig. 1 is a visualization of these packing styles in the simulation box. Among these, arrays of randomly packed spheres were created by using event-driven molecular dynamics and a modified Lubachevsky–Stillinger algorithm [45]. The spheres were considered to be rigid and impermeable.

This method involves following the trajectory of point particles in the Lagrangian framework, as they travel in a certain flow field. The particles are assumed to be passive, i.e., they do not affect the flow field, and massless, because the dispersed NP suspensions are quite dilute in this study (on the order of 100 ppm), so the presence of the NPs is not expected to affect the flow. In the current LPT simulations, the particle-particle interactions are not taken into account. This means that the dispersion of a single particle is not affected by other particles, and that particles do not have a chance to agglomerate and get bigger in size. Given that no particle-particle interaction is considered, the only factor hindering the particle dispersion is adsorption. The NP motion, therefore, is a combination of convective transport (contributed by the velocity field that is obtained from the LBM flow simulation) and diffusive transport (contributed by Brownian motion). As a consequence, the equation of motion of particles is mathematically described as follows [41,47,48]

~ ~t þ D~ X tþDt ¼ ~ X t þ Dt U X

ð7Þ

~t is the particle velocity at time t, D~ where U X is the travel distance by diffusive transport, and Dt is the time step. As particles propagate in the flow field, Dt does not change, therefore we employ a static time-stepping approach. The chosen Dt is small enough in order to keep particles away from penetrating into the pore wall x within one time step. Typically, Dt must satisfy Dt < UDmax , where Umax is the maximum fluid velocity in the pore space and Dx is the resolution of the simulation grid. Eq. (7) provides the new particle position at time t + Dt from the old position at time t. It is very

Fig. 1. Types of sphere packing employed in this study. From left to right: FCC, BCC, SC, and randomly packed spheres.

Table 1 Simulation configurations and conditions. Simulation 1 2 3 4 5 6 7

Box size [l m3] 3

375 17983 17983 14143 11433 10003 3753

dP/L [Pa/m]

Packing manner

Dsphere [mm]

Porosity

Presented in

687 100 100–500 100 100 100 100

Random Random Random FCC BCC SC Random

0.15 1 1 1 1 1 0.15

0.35 0.55 0.55 0.267 0.325 0.476 0.35

Fig. 3 Figs. 4 and 5 Figs. 7 and 8 Fig. 8 Fig. 8 Fig. 8 Figs. 6, 9 and 10

322

N.H. Pham et al. / International Journal of Heat and Mass Transfer 72 (2014) 319–328

rare that the position of a particle would coincide with a computational mesh node, so that it could assume the fluid velocity at that node. In the rest of the cases, a trilinear interpolation scheme is used to calculate the particle velocity, Ut, at the particle location from the eight neighboring grid nodes. The Brownian motion of each of the particles is considered to be a sequence of random jumps that are distributed following a normal distribution with zero mean and a standard deviation r, denoted as N(0,r). The standard deviation is directly related to the molecular diffusivity with Einstein’s theory for Brownian motion that results in the following equation in each of the three space dimensions [41,49]

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ¼ 2D0 Dt ¼

rffiffiffiffiffiffiffiffiffiffiffi 2mDt Sc

medium with a velocity equal to the mean fluid velocity (or pore velocity) without adsorption. The concentration ratio (C/C0) is then interpolated from the number of particles injected and the number of particles found in each PV interval as follows

C N number of PV injected ¼ C 0 N0 PV interval

ð9Þ

where C is the particle concentration in each PV interval, C0 is the initial particle concentration in the solution, N is the number of particles captured in each PV interval in LPT simulation, and N0 is the number of particles, initially released in the LPT simulation. 2.3. Experimental set up

ð8Þ

where D0 is the nominal molecular diffusivity, and Sc is the Schmidt number (i.e., the ratio of kinematic viscosity divided by molecular diffusivity). Eq. (8) accounts for the dependence of the Brownian motion on D0, which depends on the physical properties of the particle and the fluid. The LPT algorithm has been validated against Taylor–Aris diffusion theoretical predictions [41]. In order to simulate wall-particle interactions in our LPT implementation, particles are assigned a probability to be adsorbed by the pore surface upon collision with the pore wall. Thus, when a particle hits the pore wall, a random number between 0 and 1 is generated based on a uniform distribution and compared to a predefined adsorption probability, pa (0 6 pa 6 1). In order to keep the particle propagating, the random number must be larger than the adsorption probability pa. Otherwise, it will be retained by the pore wall and will remain adsorbed onto the wall. During each simulation, pa is assumed to be constant, since the injected particle suspensions in the experiments are quite dilute (100 ppm). This implies that the amount of injected particles is not enough to cover the pore surface significantly, and the adsorbed particles do not interfere with the adsorption of new particles. In this case, the change in the collector capacity to capture particles is small and can be negligible. Similarly, the same approach is applied to examine whether an adsorbed particle can escape from the pore wall and remobilize. Thus, desorption will only occur if the predefined desorption probability, pd, is larger than or equal to a randomly generated number that follows a uniform distribution between 0 and 1. This second random number is different from that generated in the adsorption process. In reality, the adsorption and desorption processes can happen at any instant, therefore it is necessary to check the conditions for particle adsorption and desorption at every time step of the LPT simulation. Utilizing the predefined pa values, we can simulate a wide range of pore surface conditions, covering both favorable and unfavorable conditions for particle adsorption. In other words, the factors that affect the particle adsorption process (such as surface roughness, surface charge, pH, etc.) are reflected by a single parameter, pa. In analogy, the detachment of adsorbed particles is determined by the particle desorption probability, pd, that can be used to predict the tail of the particle breakthrough curve. At the beginning of the LPT implementation, 100,000 particles were randomly and uniformly released at the inlet plane, which was perpendicular to the flow direction. These particles (with initially assigned pa and pd) often collided with the pore walls and underwent the adsorption–desorption process, leading to different residence time for the particles. At the outlet plane, one simulation box length away from the inlet, exit particles were captured and counted within an interval of 0.05 dimensionless particle residence time, that is related to the pore volume (PV). This dimensionless time scale is the ratio of the particle residence time over the average residence time at which a particle travels through the porous

2.3.1. Nanoparticles The NPs used in this work are purified multi-walled carbon nanotubes (PMWCNTs), stabilized by a mixture of polyvinyl–pyrrolidone (PVP with molecular weight = 40,000, Sigma Aldrich Lot#080M0242 V) and a stabilizing polymer obtained from DOW Chemical Corporation (Lot#B-1222LABS) subject to a confidentiality clause that does not allow us to determine its composition. The PMWCNTs were obtained after purification of MWCNTs-Al2O3 particles with the process described in detail in Baez et al. [50]. The obtained PMWCNTs were stabilized by sequentially adding PVP and the DOW proprietary polymer in 10% standard API brine (i.e., 8 wt% NaCl and 2 wt% CaCl2). Briefly, the procedures involved first sonication (sonicator model FB505, Fisher Scientific) of a solution of 100 ppm PMWCNTs and 250 ppm PVP in API brine at 25% amplitude for 2 h. After 750 ppm of the DOW proprietary polymer was added to the solution, a second run of sonication took place for another 30 min. Finally, the dispersion was centrifuged for 1 h at 2000 rpm to remove any excessively large particle agglomerates that were not adequately exfoliated during multiple sonications [50]. 2.3.2. Column packing and propagation experiment Inert glass beads of 1 mm diameter were packed into a low pressure glass column (Kimble Chase, Kontes Chromaflex) of 2.54 cm inner diameter. Dry packing was employed until the column height reached two targeted heights: either 2.54 cm or 7.62 cm (one or three-inch height), yielding a porosity of 0.55 (see Fig. 2 for experimental set up). The packed column was then saturated by injecting 10 PV of the 10% brine to ensure equal ion concentration in the effluent and the injection. A Masterflex L/S peristaltic pump (Cole Parmer) was used to inject experimental fluid from a reservoir into the column. Sequentially, either a pulse injection (3in-column) or a 2 PV of continuous injection (1in-column) of the dispersed PMWCNT solution was injected from bottom to top at the volumetric flow rate of 0.3 mL/min, followed by a 5 PV post-flush with the particle-free API brine solution. During the injection with the particle suspension, the column remained saturated with prior API brine. Effluent samples were automatically collected at every 0.5 PV (for the 1in-column) or 0.25 PV (for the 3in-column) using a timed sample collector. Particle concentration in each PV interval was then measured by a UV–vis spectrophotometer (Thermoscientific, Genesys10s) using 3.5 mL quartz cuvette by comparing the readings with a calibration curve prepared daily. 3. Results and discussion As described earlier in the LPT implementation section, particle breakthrough behavior is subject to change upon different values of pa and pd. The numerically-obtained particle breakthrough curves are presented in Fig. 3 for different pa and pd at two different

323

N.H. Pham et al. / International Journal of Heat and Mass Transfer 72 (2014) 319–328

Pressure Gauges 30psi

Sample collector

Glass Column

Sand Packing

Pump

Fig. 2. Schematic diagram of column experiments with glass bead packing.

injection modes: a pulse of 1.76  105 PV and a 5 PV injection. The properties of water at room temperature were used to describe the fluid that was passed over the sphere packing by a pressure difference of 687 Pa/m. Released nanoparticles were assumed to be spherical with 10 nm in diameter, resulting in the Sc = 21,880 based upon the Stokes–Einstein relation for particle diffusion with the Cunningham correction [51]

D0 ¼

C n kB T 6plr

ð10Þ

where Cn is the Cunningham correction factor, kB is the Boltzmann constant, T is the absolute temperature, l is the fluid dynamic viscosity, and r is the particle radius. In order to simulate a continuous injection into the porous medium, data from the instantaneous injection were used instead of continually adding NPs in the simulation in each time step. For example, in order to determine how many particles migrate through the column from the initial injection at time t0 to time t0 + t, the number of particles that go through the column for pulse injections at time t0, t0 + Dt, t0 + 2Dt, t0 + 3Dt . . . t0 + t are calculated and summed up.

Pulse injection

5PV injection

2.5e-5

1.0

2.0e-5

0.8

C/C0

C/C0

3.0e-5

1.5e-5

0.6

1.0e-5

0.4

5.0e-6

0.2 0.0

0.0 0

1

2

3

4

5

6

0

2

4

6

PV

10

100

Cumulative recovery [%]

100

Cumulative recovery [%]

8

PV

80 60 40 20

80 60 40 20 0

0 0

1

2

3

PV

4

5

6

0

2

4

6

8

10

PV

Fig. 3. Breakthrough and cumulative recovery curves of particles passing through a randomly packed array of spheres with 0.35 porosity. . pa = 0, pd = 0; d pa = 0.01, pd = 0; h pa = 0.01, pd = 0.0001. The panels on the left are for the case of an instantaneous pulse injection, and the panels on the right are for a 5 pore volume (PV) injection.

324

N.H. Pham et al. / International Journal of Heat and Mass Transfer 72 (2014) 319–328

In either pulse or 5 PV injection scenario, the most favorable case in which particles were free from adsorption (pa = 0) is shown. It is obvious that in this case no particle retention is recorded (100% recovery of the particles is observed, as can be seen in the cumulative recovery curve – filled triangles) irrespective of what the injection manner was, and the delay of particle breakthrough is likely due to some particles that move through the column in regions with slow velocity. This case corresponds to the breakthrough of conservative tracer particles in common column experiments. In contrast, more than 74% of particles are lost if pa = 0.01 with no desorption (pd = 0, Fig. 3 – filled circles). This case is equivalent to conducting a column experiment without post-flushing the column with particle-free solution. The case in which particles are allowed to detach from the wall and remobilize is also presented, utilizing a specified desorption probability for attached particles (pa = 0.01 and pd = 0.0001). Similar to the case of tracer particles, 100% particles are recovered in this case, and the breakthrough curve is skewed to the right with a long tail. It can be anticipated that, for non-zero pa and pd, the larger the pa/pd ratio is, the longer the right tail of the breakthrough curve. This occurs when the desorption rate is significantly smaller than the adsorption rate, causing a delay in the elution of the particles. Ideally, if the particles are all in constant contact with the pore walls and the adsorption rate is of first order, pa in our model can be linked to the nominal adsorption rate constant k0 as expressed below [41]

  1 1 k0 ¼ ln 1  pa Dt

ð11Þ

3.1. Validation with experiments In Figs. 4 and 5, breakthrough curves of tracer particles generated by LBM/LPT simulations (pa = 0) are compared to those of PMWCNTs obtained by experimentally injecting the dispersed PMWCNT through the bead-packed columns. Fig. 4 is a presentation of results for a pulse injection of 0.087 PV of suspension through the 7.62 cm high column (the 3in column), whereas Fig. 5 is a presentation of results for 2 PV injection of the same suspension. It is noticed from Fig. 4 that PMWCNTs are retained by neither pore surface deposition nor size exclusion, and the cumulative recovery curves plateau at 100% in both cases. Early breakthrough in the experimental column is observed when the experimental breakthrough curve reaches a peak level (C/C0) at about 0.12 against 0.1 from the simulation prediction. In contrast to this difference, the position of peak concentration is almost identical. A deviation of the peak values is expected from the pulse injection in the bead-packed column, mainly due to instrumental constrains. In the experiment, samples were collected at every quarter of a PV and the particle concentration was measured with

UV-vis light absorbance. Due to the pulse injection, the highest particle concentration in the effluent sample is fifteen times smaller than the injected concentration and collecting samples at smaller PV intervals would lead to measurements around the detection limit of the equipment. This average of the results over a fourth of a PV in the experiment can lead to differences with the simulation. However, in general, the agreement between the simulation data and the experimental results is quite reasonable. Later, to overcome the measurement uncertainty, 2 PV of suspension was injected in a continuous manner into the column, instead of a pulse input. The comparison between these experimental results and simulations is depicted in Fig. 5. Note that the column length in this case (1-in column) was one-third of the above-described experiment (i.e., 3-in column). The results in Fig. 5 reveal that the agreement between experiments and simulations is improved, and both the peak position and the concentration value at the peak are identical. 3.2. Validation with theoretical predictions In addition to comparisons to experiments, validation of our simulation results can be done by using the conventional onedimensional convective-diffusion equation [34,35]. The governing equation is based on Fick’s second law of diffusion with the convective and particle filtration terms added

8 < @C ¼ D @ 2 C2  u @C  qb @t @x / @x : @S @t

ð12Þ

¼ q ka C  kd S b

where C is the concentration of particles, D is the effective diffusivity of the particles in the porous media, u is the pore velocity, qb is the bulk density of the column, / is the porosity of the column, S is the adsorbed particle concentration, ka is the adsorption rate constant, and kd is the desorption rate constant. Exact solution for this nonlinear unsteady state equation with a pulse input boundary condition at the inlet is known for specific cases. For a pulse injection of particles and conditions of irreversible adsorption, the distribution of the particle concentration at different times and column locations downstream is given by [52,53]

" # x ðx  utÞ2 Cðx; tÞ ¼ n0 pffiffiffiffiffiffiffiffiffiffiffi expðka tÞ exp  4Dt 2 pDt3

ð13Þ

where n0 is the ratio of the total volumetric amount of particle solution injected to the volumetric flow rate, and x is the location along the flow direction. This equation is useful in order to upscale simulation results by fitting them with this equation and obtaining, thus, the two macroscopic parameters, D and ka, accounting for the irreversible adsorption kinetics. For constant fluid velocity and particle diffusivity, the particle concentration in the fluid at a certain

100

Cumulative recovery [%]

0.14 0.12 Experiment Simulation

0.10

C/C 0

@S @t

/

0.08 0.06 0.04 0.02 0.00

80 Experiment Simulation

60 40 20 0

0

1

2

3

PV

4

5

0

1

2

3

4

5

PV

Fig. 4. Simulation of conservative tracer versus dispersion of stabilized PMWCNTs in the inert glass bead packed column. The experimental column is of 3in-height and 0.55 porosity. NPs were introduced into the pore matrix as a pulse input.

325

N.H. Pham et al. / International Journal of Heat and Mass Transfer 72 (2014) 319–328

100 Experiment Simulation

0.8

C/C0

Cumulative recovery [%]

1.0

0.6 0.4 0.2 0.0

80 Experiment Simulation

60 40 20 0

0

1

2

3

4

5

6

7

0

1

2

3

PV

4

5

6

7

PV

Fig. 5. Simulation of conservative tracer (i.e., pa = pd = 0) versus dispersion of stabilized PMWCNTs in the inert glass bead packed column. The experimental column is of 1inheight and 0.55 porosity. NPs were introduced into the pore matrix as a 2 PV input.

L ¼  ln

C C0

pa = 0.005 Eqn. (10)

1.5e-6 1.0e-6 5.0e-7 0.0

ð14Þ

0.5

1.0

1.5

2.0

2.5

PV

where g is the single collector efficiency, US is the superficial fluid velocity through the column, and dc is the collector diameter. Knowing the value of g is helpful in determining the travel distance, L, in which a certain fraction (1  C/C0) of particles are captured by the porous medium



pa = 0.002

2.0e-6

0.0

3 ð1  /Þ US g 2 dc /



pa = 0.001

2.5e-6



2dc 3ð1  /Þg

ð15Þ

At this point, it is of importance to emphasize that Eqs. (14) and (15) allow the upscaling of the mesoscopic simulation results to apparent parameters that are macroscopically observed in the column experiment. Fig. 6a is an illustration of how well our simulation results fit to Eq. (13). As seen in Fig. 6a, simulation predictions agree with Eq. (13), and out of each variation of pa, a pair of ka and D is obtained. The relationship between ka and pa is presented in Fig. 6b for two different Sc number cases: 21880 and 233545, corresponding to two different sizes of NPs with diameters 10 nm and 100 nm. The pulse experiments and the data-fit to Eq. (13) also determine how the effective diffusivity, D, behaves when the deposition rates change. It turns out that there is a rather small deviation from the mean value of D at various pa in the two cases of Sc simulated (Dave = 5.83  107 ± 1.06  108 and 6  107 ± 1.03  108 cm2/ s at Sc = 21880 and 233545, respectively). This finding suggests that the effective diffusivity of the NPs is affected mostly by their physical properties (e.g., particle size) and flow conditions rather than the deposition rate. Thus, one can conduct a tracer experiment (with no adsorption) for the particle of interest and determine D, which would still hold for other cases where particle adsorption is involved. In the conventional filtration theory, particle adsorption is formulated with a first order kinetics term, as seen in Eq. (12). In our simulation, this first order rate law is simulated by creating random numbers and comparing them to the predefined pa, and this approach is found to be consistent with the theory, as evidenced by a monotonic increase of ka versus pa in the small pa

(b) 0.018 0.016 0.014

Sc = 233545 Sc = 21880

0.012

ka [1/s]

ka ¼

(a) 3.0e-6

C/C0

column location is affected by ka, and the higher the adsorption rate the smaller the particle concentration. The adsorption rate, ka and the predefined adsorption probability, pa, in our simulation have the same physical meaning, and establishing their relation is now straightforward. Additionally, so long as ka is known, calculating the single collector efficiency, g, is also straightforward. The term single collector efficiency represents the ratio of the number of nanoparticles striking a single glass bead (the collector) over the number of particles approaching the collector in a unit of time [54]. This parameter can be directly correlated to ka as [55,56]

0.010 0.008 0.006 0.004 0.002 0.000 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012

pa Fig. 6. (a) Agreement between simulation data and the predictions from Eq. (10) for a pulse release. There is irreversible adsorption at Sc = 233,545. (b) Variation of ka at various pa at two different Sc.

region on Fig. 6b. It is also seen in Fig. 6b that the slope of the ka-pa relation for small pa decreases with increasing Sc, which implies that, at a particular pa, the adsorption rate is higher at smaller Sc. This finding is consistent with the role of the Brownian motion of particles at small Sc. The smaller the Sc, the farther the particles can jump off the fluid streamlines due to random Brownian movement and the more often they can collide with the pore walls, thereby, giving particle deposition a higher chance to occur. Note that the Brownian motion is explicitly accounted in the LPT algorithm through the standard deviation r of the probability distribution that the random jumps follow. If the mass retention is only via particle adsorption as assumed in LBM/LPT simulations, it can be concluded from Eq. (10) that large particles are unfavorable for adsorption [39]. A similar observation was reported in the work of Ko and Elimelech [25] when they experimentally investigated the surface coverage efficiency under various flow rates, ionic strengths, and particle sizes. In Fig. 6b, the decrease of ka with increasing Sc and vise versa implies that the particle deposition

N.H. Pham et al. / International Journal of Heat and Mass Transfer 72 (2014) 319–328

process is reaction limited at large Sc, and that is convection limited with small Sc. The role of these two factors can be revealed with the Damköhler number, the dimensionless quantity that accounts for the characteristic fluid time over the characteristic chemical reaction time

Da ¼

pffiffiffi k0 k u

ð16Þ

where Da is the Damköhler number and k is the medium permeabilpffiffiffi ity. In the expression of Da, k is used as the length scale instead of the column length, since the column length is unreasonable for considering an infinite medium. Another possibility would be to use the bead diameter, dc, as the length scale. However, for the case of a porous column with monodisperse beads the use of the medium permeability in Eq. (16) is equivalent to the use of the representapffiffiffi tive collector diameter, since k  dc /3 =ð1  /Þ2 , found by substituting Darcy’s law into the Blake–Kozeny equation. In addition, using the medium permeability to obtain a length scale helps to differentiate between different configurations of the porous medium, including polydisperse bead diameter columns and unconsolidated porous media, and provides generality to Eq. (16). It should be further noted that the Da - type (I) is used instead of the Da - type (II), since no significant interphase mass transfer limitation exists in our results. The calculated diffusion layer thickness corresponding to each case of simulations is less than 1.5 nm (data not shown). Fig. 7 is a log–log plot of ka at various Da, when the porous medium is represented by the randomly packed spheres array at three different pa. Notice that four different simulations at constant pressure drop (constant u) but different Sc (k0 changed) and another five simulations with changing u but constant Sc were carried out to yield a series of Da and ka. In our model, k0 is determined by the simulation conditions (see Eq. (11)), so different flow velocities also indicate changes in k0. It is apparent from Fig. 7 that ka is logarithmically proportional to Da and the slope is qualitatively unchanged despite different pa. The fitting equations are recorded to be ka = 0.0031 Da1.4, ka = 0.0017 Da1.39, and ka = 0.0013 Da1.36 from the smallest to the biggest pa, respectively. The same relation is logarithmically depicted in Fig. 8 with different sphere packing morphologies and constant pa of 0.001. Series of Da were obtained by keeping u unchanged and Sc varied in five distinct simulations for each type of sphere array. The data in the plot are seen to be correlated with a power correlation. Particularly the fitting equation to data points of random packing spheres turns out to be

ka ¼ 0:0031 Da1:4 :

ð17Þ

Substitution of Eq. (16) in Eq. (17) indicates that the non-linear correlation between Da and ka can be due to a non-linear relation between the adsorption rate constant ka and the nominal adsorption rate constant k0, and/or due to a non-linear relation of ka and the pore velocity u. A non-linear relationship between ka and u would be in agreement with the theory of Tufenkji and Elimelech [56], where it is argued that ka  Us g0 , with g0 being a non-linear function of Us and, thus, of u. The value of ka depends on the flow and on the pore structure. The effect of the medium configuration on particle dispersion and deposition is remarkable, and that might be due to either the tortuosity and/or the porosity of the medium. Similarly, the desorption probability in the LBM/LPT simulation can be correlated to the desorption rate constant, kd, by matching the solution of Eq. (12) to the simulation data when particle desorption is activated in the simulation. Unfortunately, solving this equation analytically with no sink term is not trivial. We instead used the method of lines, a general finite difference technique for the numerical solution of partial differential equations, in which the central difference representation was employed for the first and second derivatives in space [57]. The accuracy of the central differencing scheme is second order [58]. Fig. 9 is a depiction of the relation of pd versus kd at identical simulation conditions as in Fig. 6, except pa stays unchanged at 0.001 (ka = 0.00631 s1) and pd is varying from 3  106 to 1. It is evident from Fig. 9 that kd monotonically increases with pd at pd < 103.

0.1 FCC BCC

0.01

SC Random packing

ka [1/s]

326

0.001 ka = 0.0031Da1.4 R² = 0.9907

0.0001

0.00001 0.01

0.1

1

10

Da Fig. 8. Dependence of adsorption rate constant on Damköhler number under different geometric morphologies. Simulations were conducted under a constant pressure drop (dP/L = 100 Pa/m), a constant adsorption probability (pa = 0.001), and a range of Sc (Sc = 384, 1000, 5000, 10000, and 20000). The fitting equation presented is for randomly packed spheres.

0.01

100

ka [1/s]

0.001 pa = 0.001 pa = 0.003 pa = 0.005

kd

10-1

0.0001

Sc = 21880 Sc = 233545

10-2

0.00001 0.01

0.1

1

10

Da Fig. 7. Relation between the Damköhler number and the adsorption rate constant in a randomly packed array of spheres. Data obtained from simulations with constant pressure drop (dP/L = 100 Pa/m) and variable Sc (Sc = 384, 1000, 5000, 10000, and 20000) and simulations with constant Sc (Sc = 384) and different pressure drops (dP/L = 200, 300, 400, and 500 Pa/m).

10-3 10-6

10-5

10-4

10-3

10-2

10-1

100

pd Fig. 9. Variation of kd at various pd for two different Sc at reversible adsorption condition. The value of the pa was fixed at 0.001 and pd was varied from 3  106 to 1. Other than that, the simulation conditions were identical to those depicted in Fig. 6.

327

N.H. Pham et al. / International Journal of Heat and Mass Transfer 72 (2014) 319–328

20

20

Sc = 21880 Sc = 233545

Elapsed time [PV]

Elapsed time [PV]

25

15 10 5 0 10-6

10-5

10-4

10-3

10-2

10-1

100

Sc = 21880 Sc = 233545

15

10

5

0 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

K = ka/kd

pd

Fig. 10. Elapsed time for complete particle recovery versus pd (left plot) and equilibrium constant, K (right plot).

This behavior is expected, as first order kinetics are applied to simulate the desorption process. Notice that, in our LPT algorithm, the desorption process is allowed to take place after the first particle gets adsorbed on the sphere, and desorption of a particle is not hindered by others that are moving toward the same collector. There is no doubt that the delay in the particle breakthrough in the absence of size exclusion effects is controlled by the probability of desorption, pd. In other words, particle retention in this case is desorption limited. Fig. 10 is a presentation of the dimensionless elapsed time (presented in terms of PV) at which all the particles are recovered for various desorption rates (left panel) and equilibrium coefficients (right panel) at two different Sc numbers. The equilibrium coefficient is the ratio of adsorption rate constant over desorption rate constant, and is denoted by K. The plot on the left indicates that the delay in particle breakthrough is pronounced at pd < 103 for the two studied Sc. Recall that kd is almost constant at pd > 103 as seen in Fig. 9, and thereby the data presented in the left plot of Fig. 10 are not contradictory. It is also seen in this plot that particles with higher Sc cause an earlier breakthrough than those with smaller Sc. Longer retention of smaller Sc particles is expected due to higher collision frequency with the pore surfaces with the same probability of adsorption and desorption. Likewise, the plot on the right demonstrates the dependence of the elapsed time scale on different K. In contrast to the tendency recorded on the left one, particles with smaller Sc demonstrate an earlier breakthrough at a certain K, and the difference is more significant when K increases. 4. Conclusions A numerical approach using the Lattice Boltzmann method in association with a Lagrangian particle tracking algorithm is presented herein to simulate the transport and kinetics of nanoparticles inside the pore space of packed porous media. The packing is represented by packing solid, impermeable spheres in a periodic computational box. By applying the probability of adsorption and desorption concept, three different cases of particle-pore surface interaction, namely no adsorption, irreversible adsorption, and reversible adsorption have successfully been simulated. Simulation results have been found to be in agreement with those predicted by the conventional filtration equation in all cases of the interaction. This fact allows simulation results to be upscaled and to determine macroscopic parameters by fitting simulation data to equations. Additionally, particle breakthrough of PMWCNTs in inert glass bead packed column experiments agree with the breakthrough predicted by the simulation of a conservative tracer. Therefore, this approach can provide a rapid qualitative and quantitative investigation into microfluidics porous systems, and can be generalized to other more sophisticated porous media like sandstone.

Acknowledgments This work was supported by the Advanced Energy Consortium: http://www.beg.utexas.edu/aec/ Member companies include BP America Inc., BG Group, Petrobras, Schlumberger, Statoil, Shell, and Total. The use of computing facilities at the University of Oklahoma Supercomputing Center for Education and Research (OSCER) and at XSEDE (under allocation CTS-090025) is also gratefully acknowledged. We would also like to acknowledge useful discussions with Professor Daniel Resasco at the University of Oklahoma.

References [1] R.J. Aitken, M.Q. Chaudhry, A.B.A. Boxall, M. Hull, Manufacture and use of nanomaterials: current status in the UK and global trends, Occup. Med. 56 (5) (2006) 300–306. [2] N.M. Franklin, N.J. Roqers, S.C. Apte, G.E. Batley, G.E. Gadd, P.S. Casey, Comparative toxicity of nanoparticulate ZnO, bulk ZnO, and ZnCl2 to a freshwater microalga (Pseudokirchneriella subcapitata): the importance of particle solubility, Environ. Sci. Technol. 41 (24) (2007) 8484–8490. [3] A. Corma, P. Atienzar, H. García, J.Y. Chane-Ching, Hierarchically mesostructured doped CeO2 with potential for solar-cell use, Nat. Mater. 3 (6) (2004) 394–397. [4] Q. Fu, A. Weber, M. Flytzani-Stephanopoulos, Nanostructured Au–CeO2 catalysts for low-temperature water–gas shift, Catal. Lett. 77 (1) (2001) 87–95. [5] V.D. Kosynkin, A.A. Arzgatkina, E.N. Ivanov, M.G. Chtoutsa, A.I. Grabko, A.V. Kardapolov, N.A. Sysina, The study of process production of polishing powder based on cerium dioxide, J. Alloys Compd. 303–304 (2000) 421–425. [6] F.E. Livingston, H. Helvajian, Variable UV laser exposure processing of photosensitive glass-ceramics: maskless micro- to meso-scale structure fabrication, Appl. Phys. A: Mater. Sci. Process. 81 (8) (2005) 1569–1581. [7] A.M. Schrand, M.F. Rahman, S.M. Hussain, J.J. Schlager, D.A. Smith, A.F. Syed, Metal-based nanoparticles and their toxicity assessment, WIREs Nanomed. Nanobiotechnol. 2 (5) (2010) 544–568. [8] H.C. Huang, K. Rege, J.J. Heys, Spatiotemporal temperature distribution and cancer cell death in response to extracellular hyperthermia induced by gold nanorods, ACS Nano 4 (5) (2010) 2892–2900. [9] M.S. Mauter, M. Elimelech, Environmental applications of carbon-based nanomaterials, Environ. Sci. Technol. 42 (16) (2008) 5843–5859. [10] D.S. Sholl, J.K. Johnson, Making high-flux membranes with carbon nanotubes, Science 312 (5576) (2006) 1003–1004. [11] J. Theron, J.A. Walker, T.E. Cloete, Nanotechnology and water treatment: applications and emerging opportunities, Crit. Rev. Microbiol. 34 (1) (2008) 43–69. [12] W.M. Prickett, B.D. Van Rite, D.E. Resasco, R.G. Harrison, Vascular targeted single-walled carbon nanotubes for near-infrared light therapy of cancer, Nanotechnology 22 (45) (2011) 455101. [13] N. Behabtu, C.C. Young, D.E. Tsentalovich, O. Kleinerman, X. Wang, A.W.K. Ma, E.A. Bengio, R.F. ter Waarbeek, J.J. de Jong, R.E. Hoogerwerf, S.B. Fairchild, J.B. Ferguson, B. Maruyama, J. Kono, Y. Talmon, Y. Cohen, M.J. Otto, M. Pasquali, Strong, light, multifunctional fibers of carbon nanotubes with ultrahigh conductivity, Science 339 (182) (2013) 182–186. [14] B. Nowack, T.D. Becheli, Occurrence, behavior and effects of nanoparticles in the environment, Environ. Pollut. 150 (1) (2007) 5–22. [15] G. Bystrzejewska-Piotrowska, J. Golimowski, P.L. Urban, Nanoparticles: their potential toxicity, waste and environmental management, Waste Manage. 29 (9) (2009) 2587–2595. [16] R. Zhang, Y. Bai, B. Zhang, L. Chen, B. Yan, The potential health risk of titania nanoparticles, J. Hazard. Mater. 211–212 (2012) 404–413. [17] A. Yethiraj, A. Striolo, Fracking: what can physical chemistry offer?, J Phys. Chem. Lett. 4 (4) (2013) 687–690.

328

N.H. Pham et al. / International Journal of Heat and Mass Transfer 72 (2014) 319–328

[18] X. Jiang, M. Tong, R. Lu, H. Kim, Transport and deposition of ZnO nanoparticles in saturated porous media, Colloids Surf., A 401 (2012) 29–37. [19] Z. Li, E.S. Demessie, A.A. Hassan, G.A. Sorial, Transport and deposition of CeO2 nanoparticles in water-saturated porous media, Water Res. 45 (15) (2011) 4409–4418. [20] N. Solovitch, J. Labille, J. Rose, P. Chaurand, D. Borschneck, M.R. Wiesner, J. Bottero, Concurrent aggregation and deposition of TiO2 nanoparticles in a sandy porous media, Environ. Sci. Technol. 44 (13) (2010) 4897–4902. [21] X. Jiang, M. Tong, H. Kim, Influence of natural organic matter on the transport and deposition of zinc oxide nanoparticles in saturated porous media, J. Colloid Interface Sci. 386 (1) (2012) 34–43. [22] S. Sirivithayapakorn, A. Keller, Transport of colloids in saturated porous media: a pore-scale observation of the size exclusion effect and colloid acceleration, Water Resour. Res. 39 (4) (2003) 1109. [23] M. Elimelech, Effect of particle size on the kinetics of particle deposition under attractive double layer interactions, J. Colloid Interface Sci. 164 (1) (1994) 190–199. [24] S.R. Kanel, S.R. Al-Abed, Influence of pH on the transport of nanoscale zinc oxide in saturated porous media, J. Nanopart. Res. 13 (9) (2011) 4035–4047. [25] C.H. Ko, M. Elimelech, The ‘‘shadow effect’’ in colloid transport and deposition dynamics in granular porous media: measurements and mechanisms, Environ. Sci. Technol. 34 (17) (2000) 3681–3689. [26] N. Saleh, H.J. Kim, T. Phenrat, K. Matyjaszewski, R.D. Tilton, G.V. Lowry, Ionic strength and composition affect the mobility of surface-modified Fe0 nanoparticles in water-saturated sand columns, Environ. Sci. Technol. 42 (9) (2008) 3349–3355. [27] D. Liu, P.R. Johnson, M. Elimelech, Colloid deposition dynamics in flow through porous media: role of electrolyte concentration, Environ. Sci. Technol. 29 (12) (1995) 2963–2973. [28] C. Shen, B. Li, C. Wang, Y. Huang, Y. Jin, Surface roughness effect on deposition of nano- and micro-sized colloids in saturated columns at different solution Ionic strengths, Vadose Zone J. 10 (3) (2011) 1071–1081. [29] D. Grolimund, M. Elimelech, M. Borkovec, Aggregation and deposition kinetics of mobile colloidal particles in natural porous media, Colloids Surf., A 191 (1,2) (2001) 179–188. [30] X. Liu, D.M. O’Carroll, E.J. Petersen, Q. Huang, C.L. Anderson, Mobility of multiwalled carbon nanotubes in porous media, Environ. Sci. Technol. 43 (21) (2009) 8153–8158. [31] D.P. Jaisi, M. Elimelech, Single-walled carbon nanotubes exhibit limited transport in soil columns, Environ. Sci. Technol. 43 (24) (2009) 9161–9166. [32] Y. Li, Y. Wang, K.D. Pennell, L.M. Abriola, Investigation of the transport and deposition of fullerene (C60) nanoparticles in quartz sands under varying flow conditions, Environ. Sci. Technol. 42 (19) (2008) 7174–7180. [33] D.P. Jaisi, N.B. Saleh, R.E. Blake, M. Elimelech, Transport of single-walled carbon nanotubes in porous media: filtration mechanisms and reversibility, Environ. Sci. Technol. 42 (22) (2008) 8317–8323. [34] J. Smith, B. Gao, H. Funabashi, T.N. Tran, D. Luo, B.A. Ahner, T.S. Steenhuis, A.G. Hay, M.T. Walter, Pore-scale quantification of colloid transport in saturated porous media, Environ. Sci. Technol. 42 (2) (2008) 517–523. [35] Y. Tian, B. Gao, C.S. Batista, K.J. Ziegler, Transport of engineered nanoparticles in saturated porous media, J. Nanopart. Res. 12 (7) (2010) 2371–2380. [36] P.R. Johnson, M. Elimelech, Dynamics of colloid deposition in porous media: blocking based on random sequential adsorption, Langmuir 11 (3) (1995) 801– 812. [37] J.E. Saiers, J.M. Hornberger, First- and second-order kinetics approaches for modeling the transport of colloidal particles in porous media, Water Resour. Res. 30 (9) (1994) 2499–2506.

[38] Z. Adamczyk, B. Siwek, M. Zembala, P. Belouschek, Kinetics of localized adsorption of colloid particles, Adv. Colloid Interface Sci. 48 (1994) 151–280. [39] S.A. Bradford, J. Simunek, M. Bettahar, M.TH. VanGenuchten, S.R. Yates, Modeling colloid attachment, straining, and exclusion in saturated porous media, Environ. Sci. Technol. 37 (10) (2003) 2242–2250. [40] R.S. Voronov, S.B. VanGordon, V.I. Sikavitsas, D.V. Papavassiliou, Computational modeling of flow-induced shear stresses within 3D saltleached porouss caffolds imaged via micro-CT, J. Biomech. 43 (7) (2010) 1279–1286. [41] R.S. Voronov, S.B. VanGordon, V.I. Sikavitsas, D.V. Papavassiliou, Efficient Lagrangian scalar tracking method for reactive local mass transport simulation through porous media, Int. J. Numer. Methods Fluids 67 (4) (2010) 501–517. [42] Y.H. Qian, D. D’Humières, P. Lallemand, Lattice BGK models for Navier–Stokes equation, Europhys. Lett. 17 (6) (1992) 479–484. [43] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (3) (1954) 511–525. [44] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364. [45] M. Skoge, A. Donev, F.H. Stillinger, S. Torquato, Packing hyperspheres in highdimensional Euclidean spaces, Phys. Rev. E 74 (4) (2006) 041127. [46] R.S. Voronov, Local fluid stress and nutrient transport effects via simulations of perfusion through bone tissue engineering scaffolds, (Ph.D. thesis), University of Oklahoma, Norman, OK, 2010. [47] G.V. Kolmakov, R. Revanur, R. Tangirala, T. Emrick, T.P. Russell, A.J. Crosby, A.C. Balazs, Using nanoparticle-filled microcapsules for site-specific healing of damaged substrates: creating a ‘‘repair-and-go’’ system, ACS Nano 4 (2) (2010) 1115–1123. [48] R. Verberg, A. Alexeev, A.C. Balazs, Modeling the release of nanoparticles from mobile microcapsules, J. Chem. Phys. 125 (22) (2006) 224712. [49] D.V. Papavassiliou, Scalar dispersion from an instantaneous line source at the wall of a turbulent channel for medium and hight Prandtl number fluids, Int. J. Heat Fluid Flow 23 (2) (2002) 161–172. [50] J.L. Baez, M.P. Ruiz, J. Faria, J.H. Harwell, B. Shiau, D.E. Resasco, Stabilization of interfactically-active-nanohybrids/polymer suspensions and transport through porous media, in: 18th SPE Improved Oil Recovery Symposium, Society of Petroleum Engineers, Tulsa, OK, 2012, p. 154052. [51] M. Shapiro, H. Brenner, Dispersion/reaction model of aerosol filtration by porous filters, J. Aerosol Sci. 21 (1) (1990) 97–125. [52] D. Grolimund, M. Elimelech, M. Borkovec, K. Barmettler, R. Kretzschmar, H. Sticher, Transport of in situ mobilized colloidal particles in packed soil columns, Environ. Sci. Technol. 32 (22) (1998) 3562–3569. [53] D.P. Jaisi, M. Elimelech, Single-walled carbon nanotubes exhibit limited transport in soil columns, Environ. Sci. Technol. 43 (24) (2009) 9161–9166. [54] P.R. Johnson, N. Sun, M. Elimelech, Colloid transport in geochemically heterogeneous porous media: modeling and measurements, Environ. Sci. Technol. 30 (11) (1996) 3284–3293. [55] R. Kretzschmar, M. Borkovec, D. Grolimund, M. Elimelech, Mobile subsurface colloids and their role in contaminant transport, Adv. Agron. 66 (1999) 121– 193. [56] N. Tufenkji, M. Elimelech, Correlation equation for predicting single-collector efficiency in physicochemical filtration in saturated porous media, Environ. Sci. Technol. 38 (2) (2004) 529–536. [57] M.B. Cutlip, M. Shacham, The numerical method of lines for partial differential equations, CACHE News 47 (1998) 18–21. [58] D.A. Anderson, J.C. Tannehill, R.H. Pletcher, Computational Fluid Mechanics and Heat Transfer, first ed., McGraw-Hill, New York, 1984.