Hello all,
I'm a new user, hopefully you can assist me as I have work on this problem for quite sometime, I keep getting error of Error, (in dsolve/numeric/bvp/convertsys) unable to convert to an explicit first-order system
I copy my work here. If anyone need the original worksheet I can send this immadiately. Please help and I would like to say Thank you.
R:= 8.31434; # R is the gas constant (J mol^-1 K^-1)
F:= 96485; # F is the Faraday constant (A s mol^-1)
alpha1:= 0.5; # alpha is the transfer coefficient (l)
alpha2:= 0.5;
T:= 298.15; # T is the absolute temperature (K)
j0H:= 0.87; # j0 is the exchange current density (A m^-2)
j0Cl2:= 0.00001; # j0 is the exchange current density (A m^-2)
8.31434
96485
0.5
0.5
298.15
0.87
0.00001
Variables
Appliedpotential:= -1.0:
d1:= 0.0085; # d1 is the thickness of cathode boundry layer (m)
d2:= 0.02; # d2 is the distance cathode boundry layer and membrane (m)
d3:= 0.0018; # d3 is the thickness of the membrane (m)
d4:= 0.02; # d4 is the distance between membrane and anode boundry layer(m)
d5:= 0.00125; # d5 is the thickness of anode boundry layer(m)
n:=1;n2:= 2;z1:=1;z2:=2; # n is the number of electrons
alpha:= 0.5;
nu:= 10^(-5); # nu is the kinematic viscosity (m/s)
d:= 0.000015; # d is the diameter of graphite pores (m)
v:= 0.002; # v is the volume (m^3)
Ac:=0.0063; # Ac is the area of cathode (m^2)
Am:=0.0064; # Am is area of the membrane (m^2)
Aa:= Ac;
0.0085
0.02
0.0018
0.02
0.00125
1
2
1
2
0.5
1
------
100000
0.000015
0.002
0.0063
0.0064
0.0063
Concentrations
Bulk concentration of chloride, zinc(II) and nickel(II) ions in solution in mol/m3
cCl:= 5000;co2:=1e-3:
cZn2:= 30;
cNi2:= 18;pcl2:=1:ph2:=1:po2:=1;
5000
30
18
1
Equilibrium Electrode Potentials
ee1 (SHE / V) is the reversible potential for the reduction of zinc(II) an nickel(I) ions
eNi2 := evalf(-0.245+R*T/(2*F)*ln(10)*log10(abs(y1(x))));
eZn2 := evalf(-0.76+R*T/(2*F)*ln(10)*log10(abs(y2(x))));
eH:= evalf(-0.887-R*T/(2*F)*ln(10)*log10(ph2)+R*T/F*ln(10)*log10(y3(x)));
eCl2:= evalf(1.396+(0.059/2)*log10(0.001*(abs(y4(x))))-0.059*log10(0.006));
eo2:=1.229+R*T/(2*F)*ln(10)*log10(po2)+R*T/F*ln(10)*log10(abs(y3(x)));
ehocl:=1.494+R*T/(2*F)*ln(10)*log10(abs(C3(x)))+R*T/(2*F)*ln(10)*log10(abs(C5(x)))-0.0296*log10(abs(y4(x)));
-0.245 + 0.01284614433 ln(|y1(x)|)
-0.76 + 0.01284614433 ln(|y2(x)|)
-0.887 + 0.02569228866 ln(y3(x))
1.527089076 + 0.01281168722 ln(0.001 |y4(x)|)
1.229 + 0.02569228866 ln(|y3(x)|)
1.494 + 0.01284614433 ln(|C3(x)|) + 0.01284614433 ln(|C5(x)|)
0.0296 ln(|y4(x)|)
- ------------------
ln(10)
Rate Coefficients
Diffusion Coefficient (Electrochemical Systems Third Edition John Newman, Karen E Thomas-Alyea)
dCl:=2.03e-9:dCl2:=1.83e-9:dhocl:=1.28e-9:do2:=2e-9:dNi2:=5e-10:dZn2:=0.71e-9:doh:=5.26e-9:nu:=1e-6:dh:=9.312e-9:dNa:=1.334e-9:dso4:=1.065e-9:pH:=3;
3
Rate Coefficients
j0Cl:=1e-5:bcCl:=38:baCl:=77:bahocl:=39:bchocl:=39:bco:=19:bao:=38:bch:=19:bah:=19:baNi:19:bcNi:19:baZn:=50:bcZn:50:f:=25:
j0o2:=1e-8:j0hocl:=1e-6:
j0H:=0.87:j0Cl2:=1e-5:
j0Ni:=1e-4:j0Zn:=3.2e-7:
j(x):=j2(x)+j3(x)+j4(x)+j5(x):
Ionic Mobility,Ui
Error, missing operator or `;`
> mu[Ni2] := dNi2/(R*T); mu[Zn2] := dZn2/(R*T); mu[Cl2] := dCl2/(R*T); mu[H] := dh/(R*T); mu[Na] := dNa/(R*T); mu[Cl] := dCl/(R*T); mu[OH] := doh/(R*T); mu[SO4] := dso4/(R*T);
-13
2.017007023 10
-13
2.864149973 10
-13
7.382245705 10
-12
3.756473879 10
-13
5.381374738 10
-13
8.189048516 10
-12
2.121891389 10
-13
4.296224961 10
Diffusion Limiting Current Desnity
jL is the diffusion limiting current (A/m2)
jLCl2:=1.554*2*F*(dCl2)^(2/3)*(nu)^(-1/6)*f^(1/2)*1000*y4(x);
jLNi:=1.554*2*F*(dNi2)^(2/3)*(nu)^(-1/6)*f^(1/2)*1000*y1(x);
jLZn:=1.554*2*F*(dZn2)^(2/3)*(nu)^(-1/6)*f^(1/2)*1000*y2(x);
jLH:= 0.0193:
(1/2)
4486.505500 25 y4(x)
(1/2)
1889.096518 25 y1(x)
(1/2)
2386.600407 25 y2(x)
Charge Balance(Electroneutrality)
[H+] + 2* [Ni 2+] + 2* [Zn 2+] + [Na+] = [Cl-] + [OH-]+2*[SO42-]
eq1:=y3(x)+2*y1(x)+2*y2(x)+y7(x)-y4(x)-y5(x)-2*y6(x)=0;
y3(x) + 2 y1(x) + 2 y2(x) + y7(x) - y4(x) - y5(x) - 2 y6(x) = 0
>
Distribution of electrode and electrolyte phase potential, depending on j(x)
Sigma[e] is the particle phase's conductivity and Sigma[s] is the electrolyte phase's conductivity (S/m)
Phi[e] is the electrode phase potential (V), Phi[s] is the electrolyte phase potential (V), x is the distance in the direction of current flow normalised with respect to L
kappa[s](x)=evalf(F^2*((abs(2^2)*mu[Ni2]*y1(x))+(abs(2^2)*mu[Zn2]*y2(x))+(abs(1^2)*mu[H]*y3(x))+(abs(1^2)*mu[Na]*y7(x))-(abs(-1^2)*mu[Cl]*y4(x))-(abs(-1^2)*mu[OH]*y5(x))-(abs(-2^2)*mu[SO4]*y6(x))));
eq1:=y3(x)+2*y1(x)+2*y2(x)+y7(x)-y4(x)-y5(x)-2*y6(x)=0;
kappa[e]:=7.407*10^5;
ODE1:=(1000*(1-lambda)+1)*diff(phi[e](x),x)=j(x)/(kappa[e]);
ODE2:=(1000*(1-lambda)+1)*diff(phi[s](x),x)=-j(x)/(F^2*((abs(2^2)*mu[Ni2]*y1(x))+(abs(2^2)*mu[Zn2]*y2(x))+(abs(1^2)*mu[H]*y3(x))+(abs(1^2)*mu[Na]*y7(x))-(abs(-1^2)*mu[Cl]*y4(x))-(abs(-1^2)*mu[OH]*y5(x))-(abs(-2^2)*mu[SO4]*y6(x))));
kappa[s](x) = 0.007510813947 y1(x) + 0.01066535580 y2(x) + 0.03497034973 y3(x)
+ 0.005009712903 y7(x) - 0.007623476159 y4(x) - 0.01975344069 y5(x)
- 0.01599803371 y6(x)
y3(x) + 2 y1(x) + 2 y2(x) + y7(x) - y4(x) - y5(x) - 2 y6(x) = 0
5
7.40700000 10
/ d \
(1001 - 1000 lambda) |--- phi[e](x)| = 0.000001350074254 j2(x)
\ dx /
+ 0.000001350074254 j3(x) + 0.000001350074254 j4(x)
+ 0.000001350074254 j5(x)
/ d \ (j2(x) + j3(x) + j4(x) + j5(x))/
(1001 - 1000 lambda) |--- phi[s](x)| = -
\ dx /
/ / -13 -12
\9309355225 \8.068028092 10 y1(x) + 1.145659989 10 y2(x)
-12 -13
+ 3.756473879 10 y3(x) + 5.381374738 10 y7(x)
-13 -12
- 8.189048516 10 y4(x) - 2.121891389 10 y5(x)
-12 \\
- 1.718489984 10 y6(x)//
Relationship between the concentration c and the local current density
Eq2 is the H material balance, Eq3 is the Ni(II) material balance, Eq4 is the Zn(II) material balance, Eq5 is the Cl2 material balance
>
ODE3:=(1000*(1-lambda)+1)*v*n*F*((dh*diff(y3(x),x,x))+(z1*mu[H]*F*y3(x)*F*diff(phi[s](x),x,x))+(z1*mu[H]*F*diff(y3(x),x)*diff(phi[s](x),x)))=Ac*j2(x);
/ -9 / d / d \\
192.970 (1001 - 1000 lambda) |9.312 10 |--- |--- y3(x)||
\ \ dx \ dx //
/ d / d \\
+ 0.03497034973 y3(x) |--- |--- phi[s](x)||
\ dx \ dx //
-7 / d \ / d \\
+ 3.624433822 10 |--- y3(x)| |--- phi[s](x)|| = 0.0063 j2(x)
\ dx / \ dx //
ODE4:=(1000*(1-lambda)+1)*v*n2*F*((dNi2*diff(y1(x),x,x))+(z2*mu[Ni2]*y1(x)*F*diff(phi[s](x),x,x))+(z2*mu[Ni2]*F*diff(y1(x),x)*diff(phi[s](x),x)))=Ac*j3(x);
/ -10 / d / d \\
385.940 (1001 - 1000 lambda) |5 10 |--- |--- y1(x)||
\ \ dx \ dx //
-8 / d / d \\
+ 3.892218452 10 y1(x) |--- |--- phi[s](x)||
\ dx \ dx //
-8 / d \ / d \\
+ 3.892218452 10 |--- y1(x)| |--- phi[s](x)|| = 0.0063 j3(x)
\ dx / \ dx //
ODE5:=(1000*(1-lambda)+1)*v*n2*F*((dZn2*diff(y2(x),x,x))+(z2*mu[Zn2]*y2(x)*F*diff(phi[s](x),x,x))+(z2*mu[Zn2]*F*diff(y2(x),x)*diff(phi[s](x),x)))=Ac*j4(x);
/ -10 / d / d \\
385.940 (1001 - 1000 lambda) |7.1 10 |--- |--- y2(x)||
\ \ dx \ dx //
-8 / d / d \\
+ 5.526950203 10 y2(x) |--- |--- phi[s](x)||
\ dx \ dx //
-8 / d \ / d \\
+ 5.526950203 10 |--- y2(x)| |--- phi[s](x)|| = 0.0063 j4(x)
\ dx / \ dx //
ODE6:=(1000*(1-lambda)+1)*v*n2*F*((dCl2*diff(y4(x),x,x))+(z2*mu[Cl]*y4(x)*F*diff(phi[s](x),x,x))+(z2*mu[Cl]*F*diff(y4(x),x)*diff(phi[s](x),x)))=Ac*j5(x);
/ -9 / d / d \\
385.940 (1001 - 1000 lambda) |1.83 10 |--- |--- y4(x)||
\ \ dx \ dx //
-7 / d / d \\
+ 1.580240692 10 y4(x) |--- |--- phi[s](x)||
\ dx \ dx //
-7 / d \ / d \\
+ 1.580240692 10 |--- y4(x)| |--- phi[s](x)|| = 0.0063 j5(x)
\ dx / \ dx //
>
>
> ODE7 := (1000*(1-lambda)+1)*v*n*F*(doh*(diff(y5(x), x, x))+z1*mu[OH]*y5(x)*F*(diff(phi[s](x), x, x))+z1*mu[OH]*F*(diff(y5(x), x))*(diff(phi[s](x), x))) = 0;
/ -9 / d / d \\
192.970 (1001 - 1000 lambda) |5.26 10 |--- |--- y5(x)||
\ \ dx \ dx //
-7 / d / d \\
+ 2.047306907 10 y5(x) |--- |--- phi[s](x)||
\ dx \ dx //
-7 / d \ / d \\
+ 2.047306907 10 |--- y5(x)| |--- phi[s](x)|| = 0
\ dx / \ dx //
> ODE8 := (1000*(1-lambda)+1)*v*n*F*(dso4*(diff(y6(x), x, x))+z2*mu[SO4]*y6(x)*F*(diff(phi[s](x), x, x))+z2*mu[SO4]*F*(diff(y6(x), x))*(diff(phi[s](x), x))) = 0;
/ -9 / d / d \\
192.970 (1001 - 1000 lambda) |1.065 10 |--- |--- y6(x)||
\ \ dx \ dx //
-8 / d / d \\
+ 8.290425307 10 y6(x) |--- |--- phi[s](x)||
\ dx \ dx //
-8 / d \ / d \\
+ 8.290425307 10 |--- y6(x)| |--- phi[s](x)|| = 0
\ dx / \ dx //
> ODE9 := (1000*(1-lambda)+1)*v*n*F*(dNa*(diff(y7(x), x, x))+z1*mu[Na]*y7(x)*F*(diff(phi[s](x), x, x))+z1*mu[Na]*F*(diff(y7(x), x))*(diff(phi[s](x), x))) = 0;
/ -9 / d / d \\
192.970 (1001 - 1000 lambda) |1.334 10 |--- |--- y7(x)||
\ \ dx \ dx //
-8 / d / d \\
+ 5.192219416 10 y7(x) |--- |--- phi[s](x)||
\ dx \ dx //
-8 / d \ / d \\
+ 5.192219416 10 |--- y7(x)| |--- phi[s](x)|| = 0
\ dx / \ dx //
Butler-Volmer Relationship between reaction current density and overpotential
j2(x) is the current due to hydrogen evolution, j3(x) is the current due to Ni(II) reduction, j4(x) is the current due to Zn(II) reduction, j5(x) is the current due to Cl2 reduction
j2(x):=-j0H*exp(-17.85*(phi[e](x)-phi[s](x)-eH));
j3(x):=(j0Ni*(-exp(-(phi[e](x)-phi[s](x)-eNi2)*alpha*F/(R*T))))/(1+j0Ni/jLNi*exp(-alpha*F*(phi[e](x)-phi[s](x)-eNi2)/(R*T)));j4(x):=(j0Zn*(-exp(-(phi[e](x)-phi[s](x)-eZn2)*alpha*F/(R*T))))/(1+j0Zn/jLZn*exp(-alpha*F*(phi[e](x)-phi[s](x)-eZn2)/(R*T)));
j5(x):=(j0Cl2*(-exp(-(phi[e](x)-phi[s](x)-eCl2)*alpha*F/(R*T))))/(1+j0Cl2/jLCl2*exp(-alpha*F*(phi[e](x)-phi[s](x)-eCl2)/(R*T)));
-0.87 exp(-17.85 phi[e](x) + 17.85 phi[s](x) - 15.83295
+ 0.4586073526 ln(y3(x)))
- (0.0001 exp(-19.46109226 phi[e](x) + 19.46109226 phi[s](x) - 4.767967604
// 1 / -9 (1/2)
+ 0.2500000000 ln(|y1(x)|))) |1 + ----- \2.117414310 10 25 exp(
\ y1(x)
-19.46109226 phi[e](x) + 19.46109226 phi[s](x) - 4.767967604
\\
+ 0.2500000000 ln(|y1(x)|))/|
/
/ -7
- \3.2 10 exp(-19.46109226 phi[e](x) + 19.46109226 phi[s](x) - 14.79043012
\// 1 / -12 (1/2)
+ 0.2500000000 ln(|y2(x)|))/ |1 + ----- \5.363277388 10 25 exp(
\ y2(x)
-19.46109226 phi[e](x) + 19.46109226 phi[s](x) - 14.79043012
\\
+ 0.2500000000 ln(|y2(x)|))/|
/
- (0.00001 exp(-19.46109226 phi[e](x) + 19.46109226 phi[s](x) + 29.71882140
// 1 / -11 (1/2)
+ 0.2493294270 ln(0.001 |y4(x)|))) |1 + ----- \8.915624868 10 25
\ y4(x)
exp(-19.46109226 phi[e](x) + 19.46109226 phi[s](x) + 29.71882140
\\
+ 0.2493294270 ln(0.001 |y4(x)|))/|
/
Numerical solution of the differential equation system using Maple's dsolve routine
> dsys := {eq1, ODE1, ODE2, ODE3, ODE4, ODE5, ODE6, ODE7, ODE8, ODE9, y1(d1) = cNi2, y2(d1) = cZn2, y3(d1) = 5000, y4(d1) = 5000, y5(d1) = 96, y6(d1) = 1000, y7(d1) = 2000, phi[e](0) = Appliedpotential, phi[s](0) = 0};
> sol10 := dsolve(dsys, numeric, continuation = lambda, abserr = 0.1e-4, maxmesh = 8192);