by Malin Schulz, Richard Brinkhuis, Carol Crean, Richard P. Sear and Joseph L. Keddie
This Jupyter notebook has the Python code for a model that predicts when stratification will and will not occur, in the parameter space of the volume fraction of small colloid, $\phi_S$ and Peclet number for the small colloid, Pe$_S$. It also briefly summarises the equations, see main text for more details.
It relies on reading in the experimental data from a csv file in the same directory.
The model is a simple extension of that in Sear, J. Chem. Phys 148, 134909 (2018). Please see there for details. NB There are minor differences in notation between that work in JCP and this. We will highlight the novel part in this work.
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from scipy.optimize import brentq
#
filename="./Master_plot_table.csv"
print('NB needs csv with file =',filename)
NB needs csv with file = ./Master_plot_table.csv
Treatment of jamming is just as in Sear, J. Chem. Phys (2018)
Volume fraction of small colloid at descending interface, before the onset of jamming is
$$ \mbox{volume fraction of small colloid at descending air/water interface}=\phi_S(1+\mbox{Pe}_St^*) $$for Peclet number Pe$_S=v_{ev}H/D_S$ and reduced time $t^*=v_{ev}t/H$. Dropping 1 inside () and requiring that jamming occur at the latest at $t^*=1$ we have
$$ \phi_S\mbox{Pe_S}=\phi_{jam} $$and
$$ \phi_S>\phi_{jam}/\mbox{Pe_S}~~~~\mbox{Equation (1)} $$for jamming to occur during drying. This is equation (11) of Sear (2018).
From Equation A1 in the Appendix of Sear (2018) we have that immediately below a descending jamming front (of small colloidal particles) the volume fraction of small particles is
$$ \phi_S(\Delta z)=\phi_S+(\phi_{jam}-\phi_S)\exp[-\Delta z/\lambda_G] $$for $\Delta z$ the distance below the (moving) jamming front, and $\lambda_G=D_s/v_{jam}$. So the derivative
$$ \frac{\partial \phi_S}{\partial z} = -\frac{\phi_{jam}-\phi_S}{\lambda_G}\exp[-\Delta z/\lambda_G] $$The diffusiophoretic speed $U$ of the large colloid is given by Equation (18) of Sear (2018)
$$ U=-(9/4)D_S(\partial \phi_S/\partial z) $$Using above equation for the gradient
$$ U=-\frac{9}{4}D_s\frac{\phi_{jam}-\phi_S}{\lambda_G}\exp[-\Delta z/\lambda_G]=(-9/4)v_{jam}(\phi_{jam}-\phi_S)\exp[-\Delta z/\lambda_G] $$as $\lambda_G=D_s/v_{jam}$.
In the 2018 we considered a large colloid very close to the jamming front, i.e., $\Delta z=0$, and then required that the large colloid's diffusiophoretic speed be at least equal to $v_{jam}$ - the speed at which the jamming front descends. Then
$$ U = v_{jam}\frac{9}{4}(\phi_{jam}-\phi_S) \ge v_{jam} $$or
$$ \frac{9}{4}(\phi_{jam}-\phi_S) \ge 1 $$$$ \phi_S \le \phi_{jam}-\frac{4}{9}\simeq 0.64-\frac{4}{9}\simeq 0.196~~~~\mbox{Equation (2)} $$In summary: Equation (1) is a lower limit to the volume fraction of small colloid at which stratification can occur--- below (in either $\phi_S$ or Pe) that limit there is not enough small colloid to generate a jammed layer early enough in drying. Equation (2) is an upper limit to the volume fraction of small colloid at which stratification can occur --- above that limit the gradients in concentration of the small colloid are too small to drive diffusiophoresis of the large colloid that is fast enough to outrun the descending jammed layer of small colloid. Thus stratification can only occur between the values set by equations (1) and (2).
Also, note that both equations rely on numerous approximations, some of which are discussed in Sear and Warren (2017). Equation (1) is just one formulation of the fact that you need enough small colloid, while equation (2) is one formulation of the fact that if the initial concentration of small colloids is high then diffusiophoresis is weak to due weak gradients (eg in limit that $\phi_S=\phi_{jam}$ then no gradients at all).
phijam=0.64
print('small colloid jams at volume fraction ',phijam)
phi0upper=phijam-4.0/9.0
print('upper limit to phi0 at which strat occurs ',round(phi0upper,3))
# array for phiS
n_pts=1500
phi0_array=np.linspace(3.0e-4,phi0upper-3.4e-4,n_pts)
# curve for onset of jamming before end of drying
peclet_lowjam=phijam/phi0_array
#
#plot using matplotlib
plt.figure(figsize=(8,5))
plt.xticks([0.0,0.05,0.1,0.15,0.2,0.25,0.3],fontsize='24')
plt.yticks([1.0,5.0,10.0,15.0,20.0,25.0,30.0],fontsize='24')
#
plt.plot(phi0_array,peclet_lowjam,color='red',lw=4)
plt.plot([phi0_array[-1],phi0_array[-1]],[peclet_lowjam[-1],1e4],color='red',lw=4)
plt.ylim([0.99,300])
plt.yscale('log')
plt.xlim([0,0.25])
plt.xlabel('$\phi_{\mathrm{S}}$',fontsize='26')
plt.ylabel('$\mathrm{Pe}_{\mathrm{S}}$',fontsize='26')
plt.title('strat. only inside red curves',fontsize=26)
plt.tight_layout()
plt.show()
small colloid jams at volume fraction 0.64 upper limit to phi0 at which strat occurs 0.196
In Sear 2018, the phoretic velocity $U$ at $\Delta z=0$ was used. This assumes that the large colloids are much smaller than the width of the accumulation zone, $R_L\ll \lambda_G$. Here we relax this constraint but still assume that we can use the phoretic speed $U$ at the centre of the large colloid. In reality the diffusiophoretic stresses will vary over the surface of the large colloid when $R_L$ and $\lambda_G$ are comparable but we assume that we can just approximate $U$ by that at the large colloid centre.
Then the maximum velocity occurs at the closest point to the jamming front, that a large colloid can approach, this is $\Delta z=R_L$. So the
$$ \mbox{maximum diffusiophoretic speed}=(9/4)v_{jam}(\phi_{jam}-\phi_S)\exp[-R_L/\lambda_G] $$and if our condition for stratification is that this maximum value exceed the downward speed of the jamming front, $v_{jam}$, we have as a condition for stratification
$$ (9/4)v_{jam}(\phi_{jam}-\phi_S)\exp[-R_L/\lambda_G]>v_{jam} $$or
$$ (9/4)(\phi_{jam}-\phi_S)>\exp[R_L/\lambda_G] $$This is new to this work.
We want to work with the variables $\phi_S$ and Pe$_S$ so we use that
$$ \lambda_G=\frac{D_S}{v_{jam}}=\frac{D_S(1-\phi_S/\phi_{jam})}{v_{ev}}=\frac{H(1-\phi_S/\phi_{jam})}{\mbox{Pe}_S} $$where we used Equation (12) of Sear (2018) for $v_{jam}=v_{ev}/(1-\phi_S/\phi_{jam})$.
So we can rewrite our condition for stratification as
$$ (9/4)(\phi_{jam}-\phi_S)>\exp\left[\frac{R_L\mbox{Pe}_S}{H(1-\phi_S/\phi_{jam})}\right] $$for our systems, we have that approximately
$$ \frac{H}{R_L}=3000 $$and then our inequality is a just a function of $\phi_S$ and Pe$_S$. We set it to an equality, and solve for Pe$_S$ as a function of $\phi_S$.
Note that this is perhaps just as much the upper limit of validity of the model, as much as it is a prediction of an upper limit on stratification.
# below is function needed for root finding
def func(Pe):
expfunc=np.exp( Pe/(H_over_Rbig*(1-phi0/phijam)) )
f=(9.0/4.0)*(phijam-phi0)-expfunc
# print(Pe,f)
return f
#######
RLarge=80.0e-9
print('large colloid radius ',RLarge*1e9,' nm')
H0=500e-6
print('consider typical height for initial film ',H0*1e6,' um')
H_over_Rbig=H0/RLarge
print('ratio of film thickness to radius of bigger colloid H/R_L ',H_over_Rbig)
minPe=0.1
maxPe=10000
Pe_upper=np.zeros(n_pts)
for i in range(0,n_pts):
phi0=phi0_array[i]
Pe_upper[i]=brentq(func,minPe,maxPe)
#
#plot using matplotlib
plt.figure(figsize=(8,5))
plt.xticks([0.0,0.05,0.1,0.15,0.2,0.25,0.3],fontsize='24')
plt.yticks([1.0,5.0,10.0,15.0,20.0,25.0,30.0],fontsize='24')
plt.plot(phi0_array,Pe_upper,color='blue',lw=4)
plt.ylim([0.99,3000])
plt.yscale('log')
plt.xlim([0,0.3])
plt.xlabel('$\phi_{\mathrm{S}}$',fontsize='26')
plt.ylabel('$\mathrm{Pe}_{\mathrm{S}}$',fontsize='26')
plt.title('upper limit to strat. due to width of accumulation zone narrowing')
plt.tight_layout()
large colloid radius 80.0 nm consider typical height for initial film 500.0 um ratio of film thickness to radius of bigger colloid H/R_L 6250.0
It should be borne in mind that the above model is very crude, in particular, it is based on the initial Pe$_S$, which uses the diffusion constant $D_S$ at the start of drying. As the viscosity increases dramatically during drying then although we assume $\lambda_G$ is a constant during, in fact it will decrease due $D_S$ decreasing as the viscosity increases. However, the model is simple and gives a physically reasonable upper bound on stratification, above which stratification cannot occur due to the width of the accumulation zone of the small colloids (below a jammed layer of small colloids) being comparable to or smaller than the size of the large colloids. Above this limit the model breaks down and we do not predict stratification.
We now want to compare model predictions with experimental data. First, we read in experimental data from .csv file
df=pd.read_csv(filename)
#print(df)
Note that we calculate viscosity at the concentration of the thickener (ASE - alkali-swellable emulsion - our thickener) in the continuous phase, i.e., $\phi_{CP}$. This is
$$ \phi_{CP}=\phi_T/(1-\phi_{TOT})=\phi_T/0.7 $$for $\phi_{TOT}=0.3$ the total volume fraction of colloid in the suspension (both small and large colloid). In other words, we measure viscosity in ASE solutions without colloids and then need to use these values in suspensions with 30% solids. So we scale up ASE concentrations by a factor of $1/0.7$.
def Viscosity_fit(T,ASE):
if(ASE==0):
# in Pa s at T=60 C for water
if(abs(T-40-273.15) < 1e-5):
Viscosity=6.531e-4
elif(abs(T-60-273.15) < 1e-5):
Viscosity=4.658e-4
else:
# with ASE so use fit
# effective ASE % increased due to colloid taking up 30% of the volume
effectiveASE=ASE/0.7
intercept, gradient= -1.97591338472721 , 4.606782537301664
Viscosity=10.0**(intercept+gradient*np.log10(effectiveASE))
# print(intercept,gradient,Viscosity)
return Viscosity
#
def StokesEinstein(T,ASE):
#
viscosity=Viscosity_fit(T,ASE)
D=1.381e-23*T/(6.0*np.pi*viscosity*RSmall)
return D
#
def Peclet_calc(T,ASE,E,H0):
D_SE=StokesEinstein(T,ASE)
PeS=E*H0/D_SE
return PeS
just some calculations to check
RSmall=14.0e-9
print('small collid radius ',RSmall*1.0e9,' nm')
print('check first row at 0% ASE')
Tcheck=df.iloc[0,3]+273.15
print('T ',Tcheck,' K')
H0check=df.iloc[0,4]*1e-6
print('H0 ',H0check,' m')
Echeck=df.iloc[0,5]*1e-9
print('E ',Echeck,' m/s')
time_evapcheck=df.iloc[0,6]
print('evaporation time defined by time to evaporate to 0.3H0 ',np.round(time_evapcheck)\
,' min',np.round(0.7*H0check/Echeck/60))
D_SE=StokesEinstein(Tcheck,0.0)
print('D of small colloid ',D_SE,' m^2/s at viscosity (pure water)')
PeSmall=Peclet_calc(Tcheck,0.0,Echeck,H0check)
print('Pe number for small particles ',np.round(PeSmall,1))
#
print('')
#
n_row=14
print('check another row ',n_row)
ASEcheck=df.iloc[n_row,0]
print('ASE % ',ASEcheck,' %')
Tcheck=df.iloc[n_row,3]+273.15
print('T ',Tcheck,' K')
H0check=df.iloc[n_row,4]*1e-6
print('H0 ',H0check,' m')
Echeck=df.iloc[n_row,5]*1e-9
print('E ',Echeck,' m/s')
time_evapcheck=df.iloc[n_row,6]
print('evaporation time defined by H0/E ',H0check/Echeck/60,' min')
print('evaporation time defined by time to evaporate to 0.3H0 ',np.round(time_evapcheck,1)\
,' min',np.round(0.7*H0check/Echeck/60))
#
print('NB to calculate Pe viscosity is that at effective % ASE = % ASE/0.7 due to 30% solid colloids')
viscosity=Viscosity_fit(Tcheck,ASEcheck)
print('Viscosity at ',ASEcheck/0.7,' % SE = ',viscosity,' Pa s')
D_SE=StokesEinstein(Tcheck,ASEcheck)
print('D of small colloid ',D_SE,' m^2/s at viscosity ',viscosity, ' Pa s')
PeSmall=Peclet_calc(Tcheck,ASEcheck,Echeck,H0check)
print('Pe number for small particles ',np.round(PeSmall,1))
small collid radius 14.0 nm check first row at 0% ASE T 333.15 K H0 0.001 m E 5.550000000000001e-07 m/s evaporation time defined by time to evaporate to 0.3H0 21.0 min 21.0 D of small colloid 3.742870948647523e-11 m^2/s at viscosity (pure water) Pe number for small particles 14.8 check another row 14 ASE % 2 % T 333.15 K H0 0.001 m E 5.550000000000001e-07 m/s evaporation time defined by H0/E 30.030030030030026 min evaporation time defined by time to evaporate to 0.3H0 21.0 min 21.0 NB to calculate Pe viscosity is that at effective % ASE = % ASE/0.7 due to 30% solid colloids Viscosity at 2.857142857142857 % SE = 1.3318827273246456 Pa s D of small colloid 1.30899609410961e-14 m^2/s at viscosity 1.3318827273246456 Pa s Pe number for small particles 42398.9
Make plot of experimental data plus model predictions
n_pts=22
print('we have 22 data points')
print('')
plt.rc('xtick', labelsize=18)
plt.rc('ytick', labelsize=20)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3,sharex=True,sharey=True,figsize=(8,4.5))
ax1.set_title('(a)',fontsize=20)
ax2.set_title('(b)',fontsize=20)
ax3.set_title('(c)',fontsize=20)
n_pts_noASE=4
for i in range(0,n_pts_noASE):
PhiS=df.iloc[i,1]
T=df.iloc[i,3]+273.15
E=df.iloc[i,5]*1.0e-9
H0=df.iloc[i,4]*1.0e-6
Pe=Peclet_calc(T,0.0,E,H0)
size=0.125*df.iloc[i,4]
indicator=df.iloc[i,7]
viscosity=Viscosity_fit(T,0)
print(i,'0 % ASE Phi_S ',PhiS,' Pe_S ',round(Pe,1),' eta ',"{:.2e}".format(viscosity),\
'Pa s T ',T,' K E',np.round(E*1e9),' nm/s ','H ',round(H0*1e6),' um',indicator)
if(indicator=='no'):
no=ax1.scatter(PhiS, Pe,c='red',s=size,marker='s',label='no')
elif(indicator=='yes'):
yes=ax1.scatter(PhiS, Pe,c='green',s=size,marker='s',label='yes')
elif(indicator=='intermediate'):
inter=ax1.scatter(PhiS, Pe,c='orange',s=size,marker='s',label='inter')
#
n_pts_ASE1=10
for i in range(0,n_pts_ASE1):
PhiS=df.iloc[n_pts_noASE+i,1]
T=df.iloc[n_pts_noASE+i,3]+273.15
E=df.iloc[n_pts_noASE+i,5]*1.0e-9
H0=df.iloc[n_pts_noASE+i,4]*1.0e-6
ASEtmp=df.iloc[n_pts_noASE+i,0]
Pe=Peclet_calc(T,ASEtmp,E,H0)
sizes=0.125*df.iloc[n_pts_noASE+i,4]
indicator=df.iloc[n_pts_noASE+i,7]
viscosity=Viscosity_fit(T,ASEtmp)
print(n_pts_noASE+i,ASEtmp,'% ASE Phi_S ',PhiS,' Pe_S ',round(Pe,1),' eta ',np.round(viscosity,3),\
'Pa s T ',T,' K E',np.round(E*1e9),' nm/s ','H ',round(H0*1e6),' um ',indicator)
if(indicator=='no'):
no=ax2.scatter(PhiS, Pe,c='red',s=sizes,marker='s',label='no')
elif(indicator=='yes'):
yes=ax2.scatter(PhiS, Pe,c='green',s=sizes,marker='s',label='yes')
elif(indicator=='intermediate'):
inter=ax2.scatter(PhiS, Pe,c='orange',s=sizes,marker='s',label='inter')
#
n_pts_ASE2=8
for i in range(0,n_pts_ASE2):
PhiS=df.iloc[n_pts_noASE+n_pts_ASE1+i,1]
T=df.iloc[n_pts_noASE+n_pts_ASE1+i,3]+273.15
E=df.iloc[n_pts_noASE+n_pts_ASE1+i,5]*1.0e-9
H0=df.iloc[n_pts_noASE+n_pts_ASE1+i,4]*1.0e-6
ASEtmp=df.iloc[n_pts_noASE+n_pts_ASE1+i,0]
Pe=Peclet_calc(T,ASEtmp,E,H0)
sizes=0.125*df.iloc[n_pts_noASE+n_pts_ASE1+i,4]
indicator=df.iloc[n_pts_noASE+n_pts_ASE1+i,7]
viscosity=Viscosity_fit(T,ASEtmp)
print(n_pts_noASE+n_pts_ASE1+i,ASEtmp,'% ASE Phi_S ',PhiS,' Pe_S ',round(Pe,1),' eta ',np.round(viscosity,3),\
'Pa s T ',T,' K E',np.round(E*1e9),' nm/s ','H ',round(H0*1e6),' um',indicator)
if(indicator=='no'):
no=ax3.scatter(PhiS, Pe,c='red',s=sizes,marker='s',label='no')
elif(indicator=='yes'):
yes=ax3.scatter(PhiS, Pe,c='green',s=sizes,marker='s',label='yes')
elif(indicator=='intermediate'):
inter=ax3.scatter(PhiS, Pe,c='orange',s=sizes,marker='s',label='inter')
#
#
ax1.legend(handles=[yes,inter,no],fontsize=16,loc=(0.030,-0.38),ncol=3)
ax1.set_xticks([0.0,0.1,0.2])
ax2.set_xticks([0.0,0.1,0.2])
plt.yscale('log')
plt.xlim([0,0.205])
plt.ylim([1,55000])
ax1.set_ylabel('Pe$_{\mathrm{s}}$',fontsize=20)
ax1.set_xlabel('$\Phi_S$',fontsize=20)
ax2.set_xlabel('$\Phi_S$',fontsize=20)
ax3.set_xlabel('$\Phi_S$',fontsize=20)
#ax1.legend()
#
ax1.plot(phi0_array,peclet_lowjam,color='blue',linestyle=':',lw=4,zorder=-3)
ax1.plot(phi0_array,Pe_upper,color='blue',lw=4,zorder=-3)
ax1.fill_between(phi0_array,peclet_lowjam,Pe_upper,color='blue',alpha=0.1)
ax2.plot(phi0_array,peclet_lowjam,color='blue',linestyle=':',lw=4,zorder=-3)
ax2.plot(phi0_array,Pe_upper,color='blue',lw=4,zorder=-3)
ax2.fill_between(phi0_array,peclet_lowjam,Pe_upper,color='blue',alpha=0.1)
ax3.fill_between(phi0_array,peclet_lowjam,Pe_upper,color='blue',alpha=0.1)
ax3.plot(phi0_array,peclet_lowjam,color='blue',linestyle=':',lw=4,zorder=-3)
ax3.plot(phi0_array,Pe_upper,color='blue',lw=4,zorder=-3)
#
#plt.tight_layout()
plt.savefig('figure2abc.png',dpi=1000,bbox_inches='tight')
plt.show()
we have 22 data points 0 0 % ASE Phi_S 0.1 Pe_S 14.8 eta 4.66e-04 Pa s T 333.15 K E 555.0 nm/s H 1000 um yes 1 0 % ASE Phi_S 0.1 Pe_S 7.6 eta 6.53e-04 Pa s T 313.15 K E 190.0 nm/s H 1000 um yes 2 0 % ASE Phi_S 0.1 Pe_S 5.7 eta 6.53e-04 Pa s T 313.15 K E 190.0 nm/s H 750 um intermediate 3 0 % ASE Phi_S 0.1 Pe_S 3.8 eta 6.53e-04 Pa s T 313.15 K E 190.0 nm/s H 500 um no 4 1 % ASE Phi_S 0.1 Pe_S 1740.1 eta 0.055 Pa s T 333.15 K E 555.0 nm/s H 1000 um yes 5 1 % ASE Phi_S 0.1 Pe_S 1740.1 eta 0.055 Pa s T 333.15 K E 555.0 nm/s H 1000 um intermediate 6 1 % ASE Phi_S 0.1 Pe_S 1740.1 eta 0.055 Pa s T 333.15 K E 555.0 nm/s H 1000 um no 7 1 % ASE Phi_S 0.1 Pe_S 633.8 eta 0.055 Pa s T 313.15 K E 190.0 nm/s H 1000 um intermediate 8 1 % ASE Phi_S 0.1 Pe_S 475.3 eta 0.055 Pa s T 313.15 K E 190.0 nm/s H 750 um no 9 1 % ASE Phi_S 0.1 Pe_S 316.9 eta 0.055 Pa s T 313.15 K E 190.0 nm/s H 500 um intermediate 10 1 % ASE Phi_S 0.1 Pe_S 190.1 eta 0.055 Pa s T 313.15 K E 190.0 nm/s H 300 um intermediate 11 1 % ASE Phi_S 0.1 Pe_S 246.4 eta 0.055 Pa s T 303.15 K E 143.0 nm/s H 500 um yes 12 1 % ASE Phi_S 0.1 Pe_S 147.8 eta 0.055 Pa s T 303.15 K E 143.0 nm/s H 300 um yes 13 1 % ASE Phi_S 0.1 Pe_S 19.2 eta 0.055 Pa s T 293.15 K E 18.0 nm/s H 300 um no 14 2 % ASE Phi_S 0.1 Pe_S 42398.9 eta 1.332 Pa s T 333.15 K E 555.0 nm/s H 1000 um intermediate 15 2 % ASE Phi_S 0.1 Pe_S 15442.0 eta 1.332 Pa s T 313.15 K E 190.0 nm/s H 1000 um intermediate 16 2 % ASE Phi_S 0.1 Pe_S 7721.0 eta 1.332 Pa s T 313.15 K E 190.0 nm/s H 500 um no 17 2 % ASE Phi_S 0.1 Pe_S 4632.6 eta 1.332 Pa s T 313.15 K E 190.0 nm/s H 300 um no 18 2 % ASE Phi_S 0.1 Pe_S 6002.7 eta 1.332 Pa s T 303.15 K E 143.0 nm/s H 500 um no 19 2 % ASE Phi_S 0.1 Pe_S 3601.6 eta 1.332 Pa s T 303.15 K E 143.0 nm/s H 300 um no 20 2 % ASE Phi_S 0.1 Pe_S 781.4 eta 1.332 Pa s T 293.15 K E 18.0 nm/s H 500 um no 21 2 % ASE Phi_S 0.1 Pe_S 468.8 eta 1.332 Pa s T 293.15 K E 18.0 nm/s H 300 um no
The shaded area above is where there are enough small colloids combined with a high enough Peclet number, to create a jammed layer of small particles, but where the accumulation zone in front of this jammed layer is wide relative to the size of the large colloidal particles.
NB Dashed curve position does depend a bit on $H$ and experimental data for range of $H$ values, but this affect is small on the log axis for Pe$_S$.