(*Experimental Variables*)
averagePoreDiameter =1000.*^-9
particleDiameter = 20.*^-9
particleConcentration = .25*5.66*^-2 (*kg/m^3. 0.25 comes from a 4:1 dilution, other value from BBI website*)
membraneThickness = 50.*^-9
transmembranePressurePsi = 10.
molarSaltConcentration = .01 
zetaParticle = 0.015 (*volts*)
zetaMembrane = 0.020 (*volts*)

(*Converting our good practices into SI or radius*)
transmembranePressure = transmembranePressurePsi * 6894.75729 (* Pa, 34473 Pa = 5 psi*)
particleRadius = particleDiameter/2.
poreRadius = averagePoreDiameter/2.

(*Physical Constants*)
viscosity = 1.*^-3 (*N s/m^2*)
faradayConstant = 9.648*^4 (*C/mole*)
gasConstant = 8.31446 (*J/(mol*K)*)
vacuumPermittivity = 8.854*^-12 (*Farads/meter*)
densityOfWater = 1000 (*kg/m^3*)
boltsmannConstant = 1.38*^-23
dielectricConstWater = 80.
temp = 300.(*K*)

debyeLength = Sqrt[(vacuumPermittivity*dielectricConstWater*gasConstant*temp)/(2*faradayConstant^2*molarSaltConcentration)] 
(*for a symmetric monovalent electrolyte such as KCl*)

(*Stokes-Einstein Diffusion coefficient:*)
diffusionCoefficient = (boltsmannConstant * temp)/(6 Pi viscosity (particleRadius))

(*Dagan Equation. Volumetric Flow rate, m^3/s :*)
solventVolumeFlow = (transmembranePressure*poreRadius^3)/(viscosity(3+(8/Pi)(membraneThickness/poreRadius))) 

solventSpeed = solventVolumeFlow/(Pi poreRadius^2) (*m/s*)

(*Dimensionless variables*)
alpha := particleRadius/poreRadius (*alpha is the ratio of particle radius to pore radius*)
(*beta = r/poreRadius dimensionless radial position*)
tau := poreRadius/debyeLength 

sigmaC =  (zetaMembrane(1+tau alpha))/alpha (*surface charge density. I found this in Jess's thesis but I couldn't find where she found it*)
sigmaS = (zetaParticle tau BesselI[1,tau])/BesselI[0,tau] (*particle charge density*)

(*Paine and Scheer values interpolated to get g and k^-1. x is the diffusive drag, y is the convective drag:*)
o= Interpolation[{1.00000,1.04393,1.09178,1.14397,1.20096,1.26330,1.33159,1.40654,1.48892,1.57966,1.67980,1.79054,1.91328,2.04963,2.20150,2.37109,2.56100,2.77430,3.01464,3.28635,3.59464,3.94578,4.34741,4.808805,5.34141,5.95938,6.68043,7.52686,8.52705,9.71752,11.14580,12.87453,14.98751,17.59862,20.86550,25.01092,30.35733,37.38467,46.83132,59.87967,78.51925,106.31670,150.22684,225.51931,372.41035,737.25652}]
x[m_] := o[(m*45./.9)+1.]
l= Interpolation[{1.00000,1.04365,1.09062,1.14122,1.19584,1.25488,1.31881,1.38816,1.46351,1.54554,1.63500,1.73276,1.83981,1.95729,2.08651,2.22898,2.38646,2.56102,2.75506,2.97143,3.21351,3.48533,3.79174,4.13856,4.53292,4.98360 ,5.50116,6.09927,6.79476,7.60918,8.57025,9.71417,11.08878,12.75849,14.81147,17.37105 ,20.61390,24.80033,30.32686,37.82235,48.33522 ,63.72853,87.60633,127.82576,204.96024,393.56585}]
y[n_] := l[(n*45./.9)+1.]

(*The following three equations come from Electrostatic effects on the partitioning of spherical colloids between dilute bulk solution and cylindrical pores, Smith & Deen.
lambda corresponds to equation, delta G to equation, and e to .*)
lambda[betaLambda_] := (Pi/2)BesselI[0,tau*betaLambda]Sum[((betaLambda^t((2. t)!))/(2.^(3t)*(t!)^2))*BesselI[t,tau betaLambda]*(tau BesselK[t+1, 2 tau] +(3/4)*BesselK[t,2 tau]),{t, 0, 20(*Should be Infinity, but on recommendation from Jess's thesis I do to 10*)}]
deltaG[betaDG_] :=(((*part1:*)((8. Pi tau alpha^4 Exp[tau alpha])/(1.+tau alpha^2))lambda[betaDG]sigmaS^2)(*part 2:*)+((4.Pi^2alpha^2 BesselI[0,tau betaDG])/((1.+tau alpha)BesselI[1,tau]))sigmaS sigmaC (*part3:*) + ((Pi BesselI[0,tau betaDG])/(tau*BesselI[1,tau]))^2 (((Exp[tau alpha]-Exp[-tau alpha])tau alpha (Coth[tau alpha]-(1./ (tau alpha))))/(1. + tau alpha))sigmaC^2)/(*part 4:*)(Pi tau Exp[-tau alpha]-(2.(Exp[tau alpha]-Exp[-tau alpha])tau alpha(Coth[tau alpha]-(1./ (tau alpha))) lambda[betaDG])/(1.+tau alpha))
e[betaE_] := poreRadius vacuumPermittivity dielectricConstWater ((gasConstant*temp)/faradayConstant)^2 * deltaG[betaE]
(*H and W come from Hindered Transport of Large Molecules in Liquid-Filled Pores by W.M. Deen. H is eq. 15, W eq. 16*)
(*Integrals solved by sampson's rule with n subdivisions as done in Jess's thesis. a can't be 0 because it will result in a 0^0 indeterminante form in the infinite sum in lambda*)
a =0.0000001
b = 1-alpha
n = 30
f [r_] :=x[alpha]^(-1)Exp[-e[r]/(boltsmannConstant*temp)]*r
g[u_] := y[alpha]^(-1)(1-u^2)Exp[-e[u]/(boltsmannConstant*temp)]*u
s = (b-a)/n
h = 2 N[(s/3)Sum[f[(a+(2j-2) s)]+4f[(a+(2j-1) s)]+f[(a+(2j) s)],{j, 0, n/2}]]
w = 4 N[(s/3)Sum[g[(a+(2p-2) s)]+4g[(a+(2p-1) s)]+g[(a+(2p) s)],{p, 0, n/2}]] 
peclet = (w solventSpeed membraneThickness)/(h diffusionCoefficient)
(*soluteSpeed = (((1-Exp[-peclet])/(w solventSpeed particleConcentration))+(1/(solventSpeed particleConcentration))Exp[-peclet])^(-1)*)
reflectionCoefficient = 1-w  (*should be 1-w but this seems to work. Hell if I know why*)
filtrateConcentrationConvection = particleConcentration*reflectionCoefficient (*filtrate concentration is first calculated assuming only convection is happening, and then that concentration gradient is used to calculate the diffusive flux and get the actual filtrate concentration*)
soluteFlux = w solventSpeed particleConcentration  ((1 - (filtrateConcentrationConvection/particleConcentration)Exp[-peclet])/(1-Exp[-peclet]))(*solute mass flux: kg/(m^2 s)*)
(*solventMassFlux = solventSpeed * densityOfWater (*1 m^3 of water is 1000kg*)*)

actualFiltrateConcentration = soluteFlux/solventSpeed
sieving = actualFiltrateConcentration/particleConcentration
soluteFlux
w
1.*10^-6
2.*10^-8
0.01415
5.*10^-8
10.
0.01
0.015
0.02
68947.6
1.*10^-8
5.*10^-7
0.001
96480.
8.31446
8.854*10^-12
1000
1.38*10^-23
80.
300.
9.74182*10^-8
2.19634*10^-11
2.64804*10^-12
3.37159
1.10265
0.0690093
InterpolatingFunction[{{1.,46.}},<>]
InterpolatingFunction[{{1.,46.}},<>]
1.*10^-7
0.98
30
0.0326667
0.91425
0.94752
7954.8
0.0524804
0.000742598
0.0452043
0.0134074
0.94752
0.0452043
0.94752