-
Notifications
You must be signed in to change notification settings - Fork 1
/
ion_ion_code_both.m
389 lines (313 loc) · 15.8 KB
/
ion_ion_code_both.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
clear
clc
close all
%% Physical constants
kb=1.380649*10^-23; % boltzmann constant
permi= 78.38*8.854188*(10^-12); % value at T=298.15 and P=0.1 MPa
e=1.602177*10^-19; % electronic charge
%% Process parameters
c0 = 5*1000*6.022e23; % bulk electrolyte concentration in #/m^3
psiD=5; % wall potential
electrolyte='NaCl';
T=298.15; % system temperature in K
rho=4/(sqrt(3)*(1.42*sqrt(3)*1e-10)^2); % particles per unit area of the wall
%% solution parameters
N=201; % number of grid points
dx=0.1; % this dx is non-dimensional
Linf=(N-1)*dx; % Distance between the plates
%% declaration of variables
filename=[electrolyte, '_c0_', num2str(c0/(1000*6.022e23)), '_N_', num2str(N), '_dx_', num2str(dx),'.mat'];
nplusf=zeros(1,N); % initialization of solution array
nminusf=zeros(1,N);
psifinal=zeros(1,N);
%% Electrolyte properties
properties=systemprops(electrolyte, T);
zplus=properties(1);
zminus=properties(2);
epsilonpp=properties(3); % LJ parameter for cation-cation interactions
epsilonmm=properties(4); % LJ parameter for anion-anion interactions
epsilonpm=properties(5); % LJ parameter for cation-anion interactions
epsilonpw=properties(6); % LJ parameter for cation-wall interactions
epsilonmw=properties(7); % LJ parameter for anion-wall interactions
sigmapp=properties(8); % LJ parameter sigma for cation-cation interactions
sigmamm=properties(9); % LJ parameter sigma for anion-anion interactions
sigmapm=properties(10); % LJ parameter sigma for cation-anion interactions
sigmapw=properties(11); % LJ parameter sigma for cation-wall interactions
sigmamw=properties(12); % LJ parameter sigma for anion-wall interactions
aplus=properties(13); % cation diameter
aminus=properties(14); % anion diameter
aplus3=properties(15);
aminus3=properties(16);
sigmamm6=sigmamm^6;
sigmapp6=sigmapp^6;
sigmapm6=sigmapm^6;
sigmamm12=sigmamm6^2;
sigmapp12=sigmapp6^2;
sigmapm12=sigmapm6^2;
aplus3c0 =aplus3*c0;
aminus3c0 =aminus3*c0;
lambda_D = sqrt(permi*kb*T/(e^2*(zplus^2*c0*zminus+zminus^2*c0*zplus))); % definition of debye length
a=(aplus+aminus)/2;
a3c0=a^3*c0;
four_pi_a3c0 = 4*pi*a3c0;
factorpsi=(e^2*a^2*c0/(kb*T*permi));
aratio=(aplus3/aminus3)^(1/3);
%% Solve
xmesh=linspace(0,Linf*a/lambda_D,N); % definition of the 1-D geometry
guess = @(x) [psiD*exp(-x/lambda_D), -(psiD/lambda_D)*exp(-x/lambda_D)]; % intial guess for the solver
solinit=bvpinit(xmesh,guess); % guess for bvp solver
options= bvpset('stats','off','RelTol',0.001);
sol=bvp5c(@(x,y)odefcn(x,y,zplus,zminus,aplus3c0,aminus3c0,aratio),@(psia,psib) bcfcn(psia,psib,psiD),solinit,options);
y=deval(sol,xmesh); %obtaining solution in arrays
sol.y=y(1,:);
sol.x=xmesh;
figure;
psi=sol.y(1,:);
% plot(sol.x*lambda_D/a, psi);
xlabel('X');
ylabel({'\bf \psi'});
fofpsi= (1+(zminus*aplus3c0*(exp(-zplus*psi)-1)./(1-zplus*aminus3c0))).^(aratio^3-1); % equations obtained from chemical potential expressions expression
gofpsi= fofpsi + (zplus*aminus3c0*(exp(zminus*psi)-fofpsi)) + (zminus*aplus3c0*fofpsi.*(exp(-zplus*psi)-1));
nplus=fofpsi.*exp(-zplus*psi)./gofpsi;
nminus=exp(zminus*psi)./gofpsi;
psiInit=sol.y(1,1:N);
xlistp=linspace(0,Linf,10*N);
xlistn=linspace(0,Linf,10*N);
xmeshnew=linspace(0,Linf,N);
nplusold=nplus;
nminusold=nminus;
psiold=psiInit;
steptol=1e-6; % non-linear least squares solver options
functol=1e-6;
optol=1e-6;
maxfcneval=500000;
maxiter=80;
lb=zeros(1,3*N); % concentrations and electric potential cannot be negative
ub=[(1/(aplus3c0*zminus))*ones(1,N), (1/(aminus3c0*zplus))*ones(1,N), inf*ones(1,N)]; % upper bound constraint on concentrations to prevent the sum of volume fractions from exceeding unity
options = optimoptions('lsqnonlin','Display','iter','StepTolerance',steptol,'FunctionTolerance',functol,'MaxIterations',maxiter,'MaxFunctionEvaluations',maxfcneval,'optimalityTolerance',optol);
xfinal=lsqnonlin(@(x) residual(x,aplus3c0,aminus3c0,aplus,aminus,a,four_pi_a3c0,factorpsi,aratio,...
N,Linf,dx,rho,psiD,zplus,zminus,epsilonpp,epsilonmm,epsilonpm,epsilonpw,epsilonmw,...
sigmapw,sigmamw,sigmamm6,sigmapp6,sigmapm6,sigmamm12,sigmapp12,sigmapm12,xlistp,xlistn)...
,x0,lb,ub,options);
nplusf=xfinal(1:N);
nminusf=xfinal((N+1):(2*N));
psifinal=xfinal((2*N)+1:(3*N));
solx=sol.x*lambda_D/a;
beta=psifinal(1,1:N);
gamma=nplusf(1,1:N);
delta=nminusf(1,1:N);
figure;
plot(solx, nplusold,"Color",'green'); % plotting cation concentration profile with and without LJ interaction
hold on
plot(xmeshnew,nplusf,"Color",'red');
hold off
xlim([0 (N-1)*dx]);
xlabel('X');
ylabel({'\bf n_+'});
legend('No LJ','Ion LJ');
figure;
plot(solx, nminusold,"Color",'green'); % plotting anion concentration profile with and without LJ interaction
hold on
plot(xmeshnew, nminusf,"Color",'red');
hold off
xlim([0 (N-1)*dx]);
xlabel('X');
ylabel({'\bf n_-'});
legend('No LJ','Ion LJ');
figure;
plot(solx, psiold,"Color",'green');
hold on
plot(xmeshnew, psifinal,"color","red"); % plotting electric potential profile with and without LJ interaction
hold off
xlim([0 (N-1)*dx]);
xlabel('X');
ylabel({'\bf \Psi'});
legend('No LJ','Ion LJ');
data=[xmeshnew' ,beta' , gamma' , delta', psiold', nplusold', nminusold'];
writematrix(data,'new_ion_ion_both_c0_5.csv');
%% defining properties for different electrolytes
function properties=systemprops(electrolyte,T)
factoreps=1000/(8.314*T);
if (strcmpi(electrolyte,'MgCl2'))
% MgCl2 (2:1)
zplus=2;
zminus=1;
aplus=0.13e-9;
aplus3 = aplus^3;
aminus=0.362e-9;
aminus3 = aminus^3;
a=(aplus+aminus)/2;
epsilonpp=3.651900*factoreps; % value is 1.47398
epsilonmm=0.076923*factoreps; % value is 0.03104
epsilonpm=3.000*factoreps;
sigmapp=(0.116290e-9)/a; % value is 0.3743
sigmamm=(0.469906e-9)/a; % value is 1.5123
sigmapm=0.300e-9/a;
epsilonww = 0.2330488*factoreps;
epsilonpw=sqrt(epsilonpp*epsilonww); % epsilon of carbon is 0.2330488 kJ/mol
epsilonmw=sqrt(epsilonmm*epsilonww);
sigmaww=0.34e-9/a;
sigmapw=(sigmapp+sigmaww)/2; % sigma of carbon is 3.400 angstroms
sigmamw=(sigmamm+sigmaww)/2;
elseif (strcmpi(electrolyte,'NaCl'))
% NaCl (1:1)
zplus=1;
zminus=1;
aplus=0.19e-9;
aplus3 = aplus^3;
aminus=0.362e-9;
aminus3 = aminus^3;
a=(aplus+aminus)/2;
epsilonpp=1.472356*factoreps;
epsilonmm=0.076923*factoreps;
epsilonpm=1.438894*factoreps;
sigmapp=(0.221737e-9)/a;
sigmamm=(0.469906e-9)/a;
sigmapm=(0.300512e-9)/a;
epsilonww = 0.2330488*factoreps;
epsilonpw=sqrt(epsilonpp*epsilonww); % epsilon of carbon is 0.2330488 kJ/mol
epsilonmw=sqrt(epsilonmm*epsilonww);
sigmaww=0.34e-9/a;
sigmapw=(sigmapp+sigmaww)/2; % sigma of carbon is 3.400 angstroms
sigmamw=(sigmamm+sigmaww)/2;
elseif (strcmpi(electrolyte,'KCl'))
% KCl (1:1)
zplus=1;
zminus=1;
aplus=0.266e-9;
aplus3 = aplus^3;
aminus=0.362e-9;
aminus3 = aminus^3;
a=(aplus+aminus)/2;
epsilonpp=1.985740*factoreps;
epsilonmm=0.076923*factoreps;
epsilonpm=1.400000*factoreps;
sigmapp=(0.230140e-9)/a;
sigmamm=(0.469906e-9)/a;
sigmapm=(0.339700e-9)/a;
epsilonww = 0.2330488*factoreps;
epsilonpw=sqrt(epsilonpp*epsilonww); % epsilon of carbon is 0.2330488 kJ/mol
epsilonmw=sqrt(epsilonmm*epsilonww);
sigmaww=0.34e-9/a;
sigmapw=(sigmapp+sigmaww)/2; % sigma of carbon is 3.400 angstroms
sigmamw=(sigmamm+sigmaww)/2;
end
properties=[zplus,zminus,epsilonpp,epsilonmm,epsilonpm,epsilonpw,epsilonmw,sigmapp,sigmamm,sigmapm,sigmapw,sigmamw,aplus,aminus,aplus3,aminus3];
end
% calculating the residual which is then minimized using a least-squares solver
function y=residual(x,aplus3c0,aminus3c0,aplus,aminus,a,four_pi_a3c0,factorpsi,aratio,...
N,Linf,dx,rho,psiD,zplus,zminus,epsilonpp,epsilonmm,epsilonpm,epsilonpw,epsilonmw,...
sigmapw,sigmamw,sigmamm6,sigmapp6,sigmapm6,sigmamm12,sigmapp12,sigmapm12,xlistp,xlistn)
nplus=x(1:N);
nminus=x(N+1:(N*2));
psi=x((2*N+1):(3*N));
y=zeros(1,3*N);
aplusn=aplus/a; % non-dimensionalized cation diameter
awall=1.42e-10/a; % non-dimensionalized carbon diameter in the wall
apluswn=(aplusn+awall)/2; % average cation-wall diameter
aminusn=aminus/a; % non-dimensionalized anion diameter
aminuswn=(aminusn+awall)/2; % average anion-wall diameter
aavgn=(aplusn+aminusn)/2;
Np=length(xlistp);
Nn=length(xlistn);
xlist=linspace(0,Linf,N);
sigmapp=sigmapp6^(1/6);
sigmamm=sigmamm6^(1/6);
sigmapm=sigmapm6^(1/6);
for i=1:N
xcurr=(i-1)*dx;
deltax_curr=xcurr-xlistp; % terms used in the denominator of the ion-ion LJ interaction equation
deltax_Linf=(Linf/2) -xlistp;
deltax_curr_4= deltax_curr.^4;
deltax_curr_10=deltax_curr.^10;
deltax_Linf_4= deltax_Linf.^4;
deltax_Linf_10=deltax_Linf.^10;
lim1=(xlistp>=0 & xlistp<=xcurr-aplusn) | (xlistp>=(xcurr + aplusn) & xlistp<=Linf); % cation-cation limits to prevent cations from overlapping
lim2=(xlistp>=0 & xlistp<=(Linf/2-aplusn)) | (xlistp>=(Linf/2+aplusn) & xlistp<=Linf);
lim3=(xlistp>=0 & xlistp<=xcurr-aavgn) | (xlistp>=(xcurr + aavgn) & xlistp<=Linf); % cation-anion limits to prevent cation-anion from overlap
lim4=(xlistp>=0 & xlistp<=(Linf/2-aavgn)) | (xlistp>=(Linf/2+aavgn) & xlistp<=Linf);
% integrands corresponding to the 4 integral limits
integrand1=zeros(Np,1);
integrand2=integrand1;
integrand3=integrand2;
integrand4=integrand3;
% using the trapz function to calculate the integrals
integrand1_val = interp1(xlist,nplus,xlistp).*((sigmapp6./(2*deltax_curr_4)) - (sigmapp12./(5*deltax_curr_10)));
integrand1(lim1)=integrand1_val(lim1);
integral1=-four_pi_a3c0*epsilonpp*zminus*trapz(xlistp,integrand1);
integrand2_val =-interp1(xlist,nplus,xlistp).*((sigmapp6./(2*deltax_Linf_4)) - (sigmapp12./(5*deltax_Linf_10)));
integrand2(lim2)=integrand2_val(lim2);
integral2=-four_pi_a3c0*epsilonpp*zminus*trapz(xlistp,integrand2);
integrand3_val = interp1(xlist,nminus,xlistp).*((sigmapm6./(2*deltax_curr_4)) - (sigmapm12./(5*deltax_curr_10)));
integrand3(lim3)=integrand3_val(lim3);
integral3=-0.5*four_pi_a3c0*epsilonpm*zplus*trapz(xlistp,integrand3);
integrand4_val = -interp1(xlist,nminus,xlistp).*((sigmapm6./(2*deltax_Linf_4)) - (sigmapm12./(5*deltax_Linf_10)));
integrand4(lim4)=integrand4_val(lim4);
integral4=-0.5*four_pi_a3c0*epsilonpm*zplus*trapz(xlistp,integrand4);
logtermplus=log((aplus3c0*nplus(i)*zminus)/(1-aplus3c0*nplus(i)*zminus-aminus3c0*nminus(i)*zplus));
y(i)= zplus*psi(i) + logtermplus...
- log(aplus3c0*zminus/(1-aplus3c0*zminus-aminus3c0*zplus))...
+ integral1 + integral2+ integral3 + integral4; % includes electrostatic interaction and LJ interaction
end
%function for mu-
for i=1:N
xcurr=(i-1)*dx;
deltax_curr=xcurr-xlistn; % terms used in the denominator of the ion-ion LJ interaction equation
deltax_Linf=(Linf/2) -xlistn;
deltax_curr_4= deltax_curr.^4;
deltax_curr_10=deltax_curr.^10;
deltax_Linf_4= deltax_Linf.^4;
deltax_Linf_10=deltax_Linf.^10;
lim1=(xlistn>=0 & xlistn<=xcurr-aminusn) | (xlistn>=(xcurr + aminusn) & xlistn<=Linf); % anion-anion limits to prevent anion-anion overlap
lim2=(xlistn>=0 & xlistn<=(Linf/2-aminusn)) | (xlistn>=(Linf/2+aminusn) & xlistn<=Linf);
lim3=(xlistn>=0 & xlistn<=xcurr-aavgn) | (xlistn>=(xcurr + aavgn) & xlistn<=Linf); % anion-cation limits to prevent anion-cation overlap
lim4=(xlistn>=0 & xlistn<=(Linf/2-aavgn)) | (xlistn>=(Linf/2+aavgn) & xlistn<=Linf);
% integrands corresponding to the 4 integral limits
integrand1=zeros(Nn,1);
integrand2=integrand1;
integrand3=integrand2;
integrand4=integrand3;
% using the trapz function to calculate the integrals
integrand1_val = interp1(xlist,nminus,xlistn).*((sigmamm6./(2*deltax_curr_4)) - (sigmamm12./(5*deltax_curr_10)));
integrand1(lim1)=integrand1_val(lim1);
integral1=-four_pi_a3c0*epsilonmm*zplus*trapz(xlistn,integrand1);
integrand2_val =-interp1(xlist,nminus,xlistn).*((sigmamm6./(2*deltax_Linf_4)) - (sigmamm12./(5*deltax_Linf_10)));
integrand2(lim2)=integrand2_val(lim2);
integral2=-four_pi_a3c0*epsilonmm*zplus*trapz(xlistn,integrand2);
integrand3_val = interp1(xlist,nplus,xlistn).*((sigmapm6./(2*deltax_curr_4)) - (sigmapm12./(5*deltax_curr_10)));
integrand3(lim3)=integrand3_val(lim3);
integral3=-0.5*four_pi_a3c0*epsilonpm*zminus*trapz(xlistn,integrand3);
integrand4_val = -interp1(xlist,nplus,xlistn).*((sigmapm6./(2*deltax_Linf_4)) - (sigmapm12./(5*deltax_Linf_10)));
integrand4(lim4)=integrand4_val(lim4);
integral4=-0.5*four_pi_a3c0*epsilonpm*zminus*trapz(xlistn,integrand4);
logtermminus=log(aminus3c0*nminus(i)*zplus/(1-aminus3c0*nminus(i)*zplus))-((aminus3c0/aplus3c0)*log((1-aminus3c0*nminus(i)*zplus-aplus3c0*nplus(i)*zminus)/(1-aminus3c0*nminus(i)*zplus)));
y(N+i)= -zminus*psi(i) + logtermminus...
- log(aminus3c0*zplus/(1-aminus3c0*zplus))-((aminus3c0/aplus3c0)*log((1-aminus3c0*zplus-aplus3c0*zminus)/(1-aminus3c0*zplus)))...
+ integral1 + integral2+ integral3 + integral4; %includes electrostatic interaction and LJ interaction
end
y((N*2)+1)= psi(1)-psiD; % potential is psiD at the left wall (boundary condition)
%solving the poisson equation
for i=2:(N-1)
y((N*2)+i)= ((psi(i+1)- 2*psi(i) + psi(i-1))/(dx^2)) + factorpsi*(zplus*zminus*nplus(i) - zminus*zplus*nminus(i));
end
y(3*N)=psi(N)-psiD; % potential is psiD at the right wall (boundary condition)
end
% function to solve the poisson equation without LJ interaction
function differ=odefcn(x,y,zplus,zminus,aplus3c0,aminus3c0,aratio)
psi=y(1);
dpsidx=y(2);
fofpsi= (1+(zminus*aplus3c0*(exp(-zplus*psi)-1)./(1-zplus*aminus3c0))).^(aratio^3-1);
gofpsi= fofpsi + (zplus*aminus3c0*(exp(zminus*psi)-fofpsi)) + (zminus*aplus3c0*fofpsi.*(exp(-zplus*psi)-1));
nplus=fofpsi.*exp(-zplus*psi)./gofpsi;
nminus=exp(zminus*psi)./gofpsi;
differ=zeros(2,1);
differ(1)=dpsidx;
differ(2)=(nminus-nplus)./(zplus+zminus); % poisson equation
end
% boundary condition function for function named odefcn
function res=bcfcn(psia,psib,psiD)
res=[ psia(1)-psiD
psib(1)-psiD
];
% psib(2)+psib(1)/Linf];
end