Electrostatic Effects on Convective Flux of Nanoparticles: A Mathematica/Matlab Implementation of W.M. Deen’s Model

ByKarl Smith

UPDATE: I have since ported the program into Matlab. If you open the following zip document(Archive) you can run Find_Sieving_Coefficient. You can select a results file and after ~5 minutes the expected sieving coefficient for the pore distribution is spat out.

Using Wolfram Alpha’s Mathematica programming language, I’ve implemented an amalgamation of different mathematical models to model the transport of charged gold nanoparticles moving through a charged pore. Infuriatingly, WordPress doesn’t allow Mathematica notebooks to be uploaded, so here’s a link to it hosted on my google drive (also a text file uploaded to the blog: Sieving model dynamic 3.14.14) EDIT:  code is now uploaded to the cloud @ http://drive.bbe.23a.myftpupload.com/index.php/apps/files?dir=/Shared/NRG/Code/Sieving%20Models . If you’re trying to understand the theory behind the separation, I produced a cleaner copy that sacrifices the dynamic updating and the ability to plot for ease of understanding (text file: Sieving Model static). If you do not have a copy of Wolfram, it’s available for free for everyone affiliated with HSEAS (and apparently @urmc.rochester.edu email addresses work too for the verification). They’ll make you jump through a few hoops, but it took me no more than twenty minutes to get it running on my computer. Wolfram allows for dynamic modulation of output, so I’ve crafted the following user interface:

Screen Shot 2014-03-14 at 4.04.20 PM

The sliders are useful for getting an intuition about what’s happening as various parameters are tweaked.

I chose mathematica because of a promotional video that told me that it was possible to upload mathematica code to the cloud, and anyone with the link (regardless if they had mathematica) could interact with a user interface similar to the one I put together. Unfortunately that functionality hasn’t been rolled out yet. So since not everyone will have access to the dynamic program, I’ve put together a few representative charts.

For the following, we assume a ‘standard separation’ with the following inputs:

particle diameter: 20 nm

particle concentration: .25 * 5.66E-2 (4 to 1 dilution of BBI values from their website)

average pore diameter: 50 nm

membrane thickness: 50 nm

transmembrane pressure: 1 psi

molar salt concentration: 10 mM

particle zetapotential: -0.015 V

membrane zetapotential: -0.020 V

Each parameter is then tweaked and plotted against the resultant sieving coefficient.

Average pore diameter v. sieving coefficient:

particle diameter

 

Transmembrane pressure v. sieving coefficient:

pressure

 

Membrane thickness v. sieving coefficient

membrane thickness

 

Molar salt concentration v. sieving coefficient:

saltconcentration

 

Membrane surface charge v. sieving coefficient:

zetaMembrane

 

The following comes (with some edits) from my qualifying exam proposal:

Filtrate concentration of gold can be understood as the ratio of solvent (water) flux to solute (gold) flux through the membrane. Solvent flux can be found from the applied pressure using the Dagan Equation:

Latex Q = \frac{\Delta P r^3}{\mu\left[ 3 + \frac{8}{\pi}\left(\frac{l}{r}\right)\right]}

(where Latex Q, Latex \mu, Latex r, and Latex l are the solvent flux, solvent viscosity, pore radius, and pore length, respectively) which is an extension of the Hagen-poiseille equation to pores that cannot be approximated as infinitely long. Solutes at least a few times bigger than the solvent molecules (gold can be as small as 5 nm diameter; a water molecule is Latex \sim0.25 nm long) can be treated as Brownian particles subject to hydrodynamic resistance. For simple diffusion (i.e. spherical particles in an unbounded fluid) the solute flux is given by Fick’s Law of Diffusion as Latex N_s = -D_\infty \frac{dC}{dz}, where Latex D_\infty is the diffusivity of a solute molecule as given by the Stokes-Einstein equation, Latex C is concentration, and Latex z is distance from the pore. Because of steric and electrostatic interactions, the free diffusion of particles is `hindered’ within the narrow confines of the pore, so that the diffusive term is actually Latex N_s = -K^{-1}D_\infty \frac{dC}{dz}, with a diffusive hindrance factor (also called the enhanced drag) Latex K^{-1}. Similarly, purely convective transport can characterized by Latex N_s = GVC, with Latex V representing the bulk flow and Latex G the convective hindrance factor (also called the lag coefficient). Solute flux is a sum of the diffusive and convective components of flux through the membrane:

Latex \text{Solute Flux }= N = -K^{-1}D_\infty \frac{dC}{dz}+GVC

Deen (Hindered transport of large molecules in liquid-filled pores – pdf hosted by blog) derived the following system of equations to describe solute flux Latex N:

Latex Pe = \frac{W\langle V \rangle L}{HD_\infty}

Latex W = 4\int_0^{1-\alpha} G (1 - \beta^2)e^{-E(\beta)/kT} \beta d \beta

Latex H = 2\int_0^{1-\alpha} K^{-1} e^{-E(\beta)/kT} \beta d \beta

Latex \langle N \rangle = W\langle V \rangle C_{\text{feed}} \frac{[1-(C_{\text{filtrate}}/C_{\text{feed}})e^{-Pe}]}{1-e^{-Pe}}

where Latex \beta = \frac{r}{\text{pore radius}} is the dimensionless radial position, Latex Pe is the Peclet number (the ratio of convective to diffusive transport in the system), Latex \alpha = \frac{\text{particle radius}}{\text{pore radius}}, and Latex E(\beta) is the boltzman energy. Latex H and Latex W are two hydrodynamic functions that modify the diffusive and convective flux (respectively) due to steric and electrostatic interactions. Note that the limiting cases for the last equation are Latex \langle N \rangle = W\langle V \rangle C_{\text{feed}} for Latex Pe >> 1 and Latex \langle N \rangle = \frac{HD_\infty}{L}(C_{\text{feed}} -C_{\text{filtrate}}) for Latex Pe << 1.

You’ll notice that the filtrate concentration Latex C_{\text{filtrate}} on the right hand side of the equation, which is where the model begins to make assumptions that are more difficult to justify. To calculate the filtrate concentration of gold (which we need to calculate the diffusive flux) , we assume that there is no concentration polarization (no ‘traffic jam’ behind the membrane due to gold moving slower than water) and that the system is in steady state. If our solution is dilute, this is a fairly reasonable assumption – even though some gold is ‘rejected’ by the pore, it diffuses rapidly back into bulk solution and the concentration at the entrance approaches the concentration in the bulk. This rejection of gold is defined by the rejection coefficient Latex \sigma _r = 1 - W which represents the fraction of gold that does not pass through the membrane, and multiplying this by our feed, or bulk, concentration of gold gives us the expected concentration of gold in the filtrate due to convection. This is then used to calculate diffusive flux across the membrane in the above equation. This is why the sieving coefficient is so high in the limit of low transmembrane pressure – the model calculates the diffusive flux across the membrane assuming that the concentration of the feed is fixed, and that the concentration of the filtrate is also fixed at 0. But in the low-pressure limit, we should be using Jess’s model anyways.

Additionally, the model currently uses the average pore diameter for all calculations. As we’ve seen generally in Dagan flow, the larger pores are disproportionately responsible for transmembrane flux, due to the no-slip condition of flow and the r^4 dependence on flow rate. As a result, the average pore diameter will likely significantly underestimate the sieving behavior of the system. This is easily fixed by running the program for each pore in the pore histogram output, but it will necessitate re-writing the program in MATLAB so it can interface with the pore image processor.

In Electrostatic effects on the partitioning of spherical colloids between dilute bulk solution and cylindrical pores (pdf, hosted on blog), Smith et. al. find

Latex E(\beta)=R_0\epsilon(RT/F)^2\Delta G

where Latex R_0, \epsilon, R, F, \Delta G are the Cylinder radius, the solvent dielectric permittivity, the gas constant, the faraday constant, and the Gibbs free energy, respectively. Latex \Delta G at constant surface charge represents the change in free energy between a charged particle in solution and the charged particle within a charged pore, and is given as the following equation:

Latex \Delta G = \frac{\frac{8\pi\tau\alpha^4e^{\tau\alpha}}{1+\tau\alpha^2}\Lambda\sigma_s^2+\frac{4\pi^2\alpha^2I_0(\tau\beta)}{(1+\tau\alpha)I_1(\tau)}\sigma_s\sigma_c+\left(\frac{\pi I_0(\tau\beta)}{\tau I_1(\tau)}\right)^2\frac{(e^{\tau\alpha}-e^{-\tau\alpha})\tau\alpha L(\tau\alpha))}{1+\tau\alpha}\sigma_c^2}{\pi\tau e^{-\tau\alpha}-\frac{2(e^{\tau\alpha}-e^{-\tau\alpha})\tau\alpha L(\tau\alpha)\Lambda))}{1+\tau\alpha}}

where Latex \tau is the ratio of pore radius to debye length, Latex L(\tau\alpha)=\coth (\tau\alpha) - \frac{1}{\tau\alpha}, Latex I is the modified bessel function of the first kind, Latex \sigma_s and Latex \sigma_c are the surface charge of the sphere and the cylinder, respectively, and when the dielectric constant of the metal is much higher than that of the solvent, Latex \Lambda can be approximated:

Latex \Lambda \approx \frac{\pi}{2} I_0(\tau\beta) \sum_{t=0}^{\infty}\frac{\beta^t(2t)!}{2^{3t}(t!)^2}I_t(\tau\beta)\times\left[ \tau K_{t+1}(2\tau) + \frac{3}{4}K_t(2\tau)\right]

where Latex K is the modified bessel function of the second kind. Note that Latex \sigma_s and Latex \sigma_c are both exclusively in the numerator of equation, meaning as either surface’s charge increase the associated change in free energy likewise increases. These numbers are found from the experimentally derived zetapotentials via the following formulas:

Latex \sigma_C = (\zeta_{membrane} (1 + \tau \alpha))/ \alpha

Latex \sigma_S = \zeta_{particle} \tau I_1(\tau) /I_0(\tau)

which I found in Jess’s thesis, but I haven’t yet been able to figure out where she found it. Further note that because the surface charges are additive, even if one of the surfaces charges were zero there would still be an electrostatic component to the difference in free energy. Only if both of the particles had zero surface charge would this expression go to zero. This expression is derived by describing a charged sphere and a charged cylinder using the linearized Poisson-Boltzman equation in both spherical and cylindrical coordinates (respectively), identifying the boundary conditions of the two systems of equations, and then performing a coordinate transformation to apply one set of boundary conditions to the next.

Latex K^{-1} and G are interpolated from Paine and Scherr’s (Drag coefficients for the movement of rigid spheres through liquid-filled cylindrical pores – pdf hosted on blog) numerical simulations using mathematica’s standard Interpolation function. It’s important to note that Paine and Scheer assume that there is no electrostatic repulsion (Latex E(\beta) = 0) within the pore for the calculation of these values. Dechadilok and Deen examine electrostatic effects on both Latex G and Latex K^{-1} ; they found that accounting for electrostatics did marginally change both values, but the greater impact of charge-based repulsion came from whether or not the particles made it into the pore in the first place. Since we are accounting for electrostatics elsewhere, Paine and Sheer is a reasonable approximation. Note also that extensive experimental evidence indicates that the viscosity of water in pores down to a size of Latex 30 \AA does not deviate significantly from its bulk value, and we can thus treat viscosity as uniform everywhere in our system.

Similar Posts

2 Comments

  1. Wow. This will take a while to digest, but I’m glad its here!

    Can you host on owncloud instead of gdrive?

    Also talk to david about the trouble uploading to WP.

Comments are closed.