A Selection of Physical Chemistry Problems Solved using

A Selection of Physical Chemistry Problems Solved using Mathematica Housam BINOUS National Institute of Applied Sciences and Technology binoushousam@yahoo. com 1 - Applied Thermodynamics 2 - Chemical Kinetics

Azeotropes for Ternary Systems Case of Acetone - Chloroform - Methanol System Gas Constant and Total Pressure : Vapor Pressure Using Antoine Equation : R=1. 987; P=760; A 1=7. 11714; B 1=1210. 595; C 1=229. 664; A 2=6. 95465; B 2=1170. 966; C 2=226. 232; A 3=8. 08097; B 3=1582. 271; C 3=239. 726; PS 1=10^(A 1 -B 1/(C 1+T)); PS 2=10^(A 2 -B 2/(C 2+T)); PS 3=10^(A 3 -B 3/(C 3+T));

Liquid phase activity coefficients from Wilson model : l 12=116. 1171; l 21=-506. 8519; l 13=-114. 4047; l 31=545. 2942; l 23=-361. 7944; l 32=1694. 0241; V 1=74. 05; V 2=80. 67; V 3=40. 73; X 3=1. -X 1 -X 2; A 12=V 2/V 1 Exp[-l 12/(R*(273. 15+T))]; A 13=V 3/V 1 Exp[-l 13/(R*(273. 15+T))]; A 32=V 2/V 3 Exp[-l 32/(R*(273. 15+T))]; A 21=V 1/V 2 Exp[-l 21/(R*(273. 15+T))]; A 31=V 1/V 3 Exp[-l 31/(R*(273. 15+T))]; A 23=V 3/V 2 Exp[-l 23/(R*(273. 15+T))]; GAM 1=Exp[-Log[X 1+X 2*A 12+X 3*A 13]+1 -(X 1/(X 1+X 2*A 12+X 3*A 13)+ X 2*A 21/(X 1*A 21+X 2+X 3*A 23)+X 3*A 31/(X 1*A 31+X 2*A 32+X 3))]; GAM 2=Exp[-Log[X 1*A 21+X 2+X 3*A 23]+1 -(X 1*A 12/(X 1+X 2*A 12+X 3*A 13)+ X 2/(X 1*A 21+X 2+X 3*A 23)+X 3*A 32/(X 1*A 31+X 2*A 32+X 3))]; GAM 3=Exp[-Log[X 1*A 31+X 2*A 32+X 3]+1 -(X 1*A 13/(X 1+X 2*A 12+X 3*A 13)+ X 2*A 23/(X 1*A 21+X 2+X 3*A 23)+X 3/(X 1*A 31+X 2*A 32+X 3))];

Modified Raoult’s law : Y 1=X 1*PS 1*GAM 1/P; Y 2=X 2*PS 2*GAM 2/P; Y 3=X 3*PS 3*GAM 3/P; Using the Mathematica’s function Find. Root to solve a system of nonlinear equations using different initial guesses : In[25]: =Find. Root[{Y 1==X 1, Y 2==X 2, P==X 1*PS 1*GAM 1+X 2*PS 2*GAM 2+X 3*PS 3*GAM 3}, {X 1, 0. 3}, {X 2, 0. 4}, {T, 40}] Out[25]={X 1 ->0. 329313, X 2 ->0. 230369, T->57. 3763} In[26]: =Find. Root[{Y 1==X 1, Y 2==X 2, P==X 1*PS 1*GAM 1+X 2*PS 2*GAM 2+X 3*PS 3*GAM 3}, {X 1, 0}, {X 2, 0. 4}, {T, 57}] Out[26]={X 1 -> 4. 068181432585476 10^-27, X 2 -> 0. 6547858067091471, T -> 53. 89598911261277} In[27]: =Find. Root[{Y 1==X 1, Y 2==X 2, P==X 1*PS 1*GAM 1+X 2*PS 2*GAM 2+X 3*PS 3*GAM 3}, {X 1, 0. 3}, {X 2, 0}, {T, 57}] Out[27]={X 1 -> 0. 7895499074011907, X 2 -> 7. 732502667736949 10^-31, T -> 55. 37684581029206} In[28]: =Find. Root[{Y 1==X 1, Y 2==X 2, P==X 1*PS 1*GAM 1+X 2*PS 2*GAM 2+X 3*PS 3*GAM 3}, {X 1, 0. 3}, {X 2, 0. 6}, {T, 63}] Out[28]={X 1 ->0. 337271, X 2 ->0. 662729, T->64. 5366}

Calculation of Binary Interaction Parameters for Wilson Model Case of Methanol-Water binary system Gas Constant and Total Pressure : P=760; R=1. 987; Modified Raoult’s law with Wilson’s model: A 12[T_]: =18. 07/40. 73 Exp[-d 1/(R (T+273. 15))]; A 21[T_]: =40. 73/18. 07 Exp[-d 2/(R (T+273. 15))]; y[x_, T_]: =x 10^(Aa-Ba/(T+Ca)) Exp[-Log[x+A 12[T] (1 -x)] +(1 -x) (A 12[T]/(x+A 12[T] (1 -x))-A 21[T]/(A 21[T] x+1 -x))]/P Aa=8. 08097; Ba=1582. 271; Ca=239. 726;

Experimental data : {Methanol liquid mole fraction, Temperature, Methanol Vapor Mole fraction} from P. C. Wankat, Equilibruim Staged Separations, Prentice Hall 1988 mydata={{0, 100, 0}, {0. 02, 96. 4, 0. 134}, {0. 04, 93. 5, 0. 23}, {0. 06, 91. 2, 0. 304}, {0. 08, 89. 3, 0. 365}, {0. 1, 87. 7, 0. 418}, {0. 15, 84. 4, 0. 517}, {0. 2, 81. 7, 0. 579}, {0. 3, 78, 0. 665}, {0. 4, 75. 3, 0. 729}, {0. 5, 73. 1, 0. 779}, {0. 6, 71. 2, 0. 825}, {0. 7, 69. 3, 0. 87}, {0. 8, 67. 6, 0. 915}, {0. 9, 66, 0. 958}, {0. 95, 65, 0. 979}, {1, 64. 5, 1}}; Use Mathematica’s function Find. Minimum to determine the binary interaction coefficients : sum. Of. Squares[data_]: = Apply[Plus, Map[{(y[#[[1]], #[[2]]] - #[[3]])^2}&, data]]] param 1=Find. Minimum[sum. Of. Squares[mydata], {d 1, 0. 1, 90}, {d 2, . 01, 1000}, Max. Iterations->300] {0. 000322454, {d 1 ->127. 624, d 2 ->484. 178}}

Isobar Vapor-Liquid Equilibrium Calculations Case of Ethanol-Water System at 760 mm. Hg Vapor Pressure Using Antoine Equation : A 1=8. 07131; B 1=1730. 630; C 1=233. 426; A 2=8. 11220; B 2=1592. 864; C 2=226. 184; PS 2=10^(A 1 -B 1/(C 1+T)); PS 1=10^(A 2 -B 2/(C 2+T)); Activity coefficients using the Van Laar Model : G 1[i_]: =Exp[A 12 (A 21 (1 -x[i])/(A 12 x[i]+A 21 (1 -x[i])))^2] G 2[i_]: =Exp[A 21 (A 12 x[i]/(A 12 x[i]+A 21 (1 -x[i])))^2] A 12=1. 6798; A 21=0. 9227;

Compute T and y for given x using a While loop and Mathematica’s function Find. Root : i=0; P=760; T=. While[i<101, {x[i]=i 0. 01, T[i]=Find. Root[P== PS 1 G 1[i] x[i]+PS 2 G 2[i] (1 -x[i]), {T, 80}][[1, 2]], y[i]=PS 1 G 1[i] x[i]/P/. T-> T[i], i++}] Create tables using Mathematica’s command Table : tb=Table[{x[i], y[i]}, {i, 0, 100}]; tb 2=Table[{x[i], T[i]}, {i, 0, 100}]; tb 3=Table[{y[i], T[i]}, {i, 0, 100}];

Mathematica’s commands List. Plot and Show are used to plot Bubble point and dew point temperatures on the same figure : plt 1=List. Plot[tb 2, Plot. Style->RGBColor[1, 0, 0], Plot. Joined-> True, Plot. Range->All] plt 2=List. Plot[tb 3, Plot. Style->RGBColor[0, 0, 1], Plot. Joined-> True, Plot. Range->All] T x, y

Mathematica’s commands List. Plot, Line and Epilog are used to plot the VLE data and the y=x line : plt 1=List. Plot[tb, Plot. Style->RGBColor[1, 0, 0], Plot. Joined-> True, Epilog-> {RGBColor[0, 1, 0], Line[{{0, 0}, {1, 1}}]}] y x

Isotherm Vapor-Liquid Equilibrium Calculations Case of Ethanol-Ethyl acetate System at 70°C Vapor Pressure Using Antoine Equation : A 1=7. 10179; B 1=1244. 951; C 1=217. 881; A 2=8. 11220; B 2=1592. 864; C 2=226. 184; PS 1=10^(A 1 -B 1/(C 1+T)); PS 2=10^(A 2 -B 2/(C 2+T)); Partial pressure using Raoult’s law : T=70; P 1[x_]: =PS 1 x; P 2[x_]: =PS 2 (1 -x); plt 3=Plot[{P 1[x], P 2[x]}, {x, 0, 1}, Plot. Style->{RGBColor[1, 0, 1], RGBColor[0, 1, 0]}]

Liquid activity coefficients using the Margules Model : A 12=0. 8557; A 21=0. 7476; G 3=Exp[(A 12+2 (A 21 -A 12) x) (1 -x)^2] G 4=Exp[(A 21+2 (A 12 -A 21) (1 -x)) x^2] plt 8=Plot[{G 3, G 4}, {x, 0, 1}, Plot. Style->{RGBColor[0, 0, 1], RGBColor[0, 0, 1]}] gi Expect positive deviation from ideality because activity coefficients are greater than 1 : x

Partial pressure using Modified Raoult’s law : P 1[x_]: =G 3 PS 1 x; P 2[x_]: =G 4 PS 2 (1 -x); plt 5=Plot[{P 1[x], P 2[x]}, {x, 0, 1}, Plot. Style->{RGBColor[0, 0, 1], RGBColor[0, 0, 1]}] Show[plt 3, plt 5] Pi x

P[x_]: =G 3 x PS 1+G 4 (1 -x) PS 2 plt 10=Plot[P[x], {x, 0, 1}, Plot. Style->RGBColor[1, 0, 1]] Plotting P versus x and y to get the isotherm VLE diagram : tbl=Table[{x G 3 PS 1/(G 3 x PS 1+G 4 (1 -x) PS 2), G 3 x PS 1+G 4 (1 -x) PS 2}, {x, 0, 1, 0. 01}]; plt 11=List. Plot[tbl, Plot. Style->RGBColor[1, 0, 1], Plot. Joined->True] Show[plt 10, plt 11] P x, y

Wei-Prater mechanism 3 -reactant triangle network : A 1=A 2=A 3=A 1 with rate constants k 12, k 21, k 23, k 32, k 31, k 13. Rate constants are not independent : k 32=k 23 (k 12/k 21) (k 31/k 13) k 12=. 5; k 21=0. 7; k 13=. 1; k 31=. 2; k 23=. 9; k 32=k 23 (k 12/k 21) (k 31/k 13) A 1 o=. 7; A 2 o=0; A 3 o=. 3; tf=10; sequil=Solve[{0==-(k 12 +k 13) A 1+k 21 A 2+k 31 A 3 , 0==k 12 A 1 -(k 21 +k 23)A 2+k 32 A 3, A 3==A 1 o+A 2 o+A 3 o-A 1 -A 2}, {A 1, A 2, A 3}]//Simplify {A 1+A 2+A 3, A 2/A 1, k 12/k 21, A 3/A 1, k 13/k 31, A 3/A 2, k 23/k 32}/. sequil Steady state solution obtained using Mathematica’s function Solve : {{A 1 ->0. 451613, A 2 ->0. 322581, A 3 ->0. 225806}} {{1. , 0. 714286, 0. 5, 0. 7}}

transient solution obtained using Mathematica’s function NDSolve : sol. WP=NDSolve[{A 1'[t]==-(k 12 +k 13) A 1[t]+k 21 A 2[t]+ k 31 (A 1 o+A 2 o+A 3 o-A 1[t]-A 2[t]) , A 2'[t]==k 12 A 1[t] -(k 21 +k 23)A 2[t]+ k 32 (A 1 o+A 2 o+A 3 o-A 1[t]-A 2[t]), A 1[0]==A 1 o, A 2[0]==A 2 o}, {A 1[t], A 2[t]}, {t, 0, tf}]; ({A 1[t], A 2[t], A 1 o+A 2 o+A 3 o-A 1[t]-A 2[t]})/. sol. WP/. t->tf {{0. 45162, 0. 322577, 0. 225802}} Plotting the solution using Mathematica’s functions Plot and Parametric. Plot : Plot[Evaluate[Table[{A 1[t], A 2[t], A 1 o+A 2 o+A 3 o-A 1[t]-A 2[t]}/. sol. WP]], {t, 0, tf}, Frame->True, Default. Font->{"Symbol-Bold", 14}, Frame. Label->{"t", "A 1, A 2, A 3"}, Plot. Range->{{0, tf}, {0, 1}}, Plot. Style->{RGBColor[1, 0, 0], RGBColor[0, 1, 0], RGBColor[0, 0, 1]}]; Parametric. Plot[Evaluate[Table[{A 2[t], A 1[t]}/. sol. WP]], {t, 0, tf}, Frame->True, Default. Font->{"Symbol-Bold", 14}, Frame. Label->{"A 2", "A 1"}, Plot. Range->{{0, 1}, {0, 1}}];

Transient solution of Wei-Prater problem :

Consecutive reactions : A 1=A 2=A 3=A 4=A 5 with rate constants k 12, k 21, k 23, k 32, k 31, k 13. . . transient solution obtained using Mathematica’s function NDSolve : k 12=k 23=k 34=k 45=1; k 21=k 32=k 43=k 54=. 1; A 1 o=1; tf=10; sol 5=NDSolve[{A 1'[t]==-k 12 A 1[t]+k 21 A 2[t], A 2'[t]==k 12 A 1[t] -(k 21 +k 23)A 2[t]+k 32 A 3[t], A 3'[t]==k 23 A 2[t] -(k 32 +k 34)A 3[t]+k 43 A 4[t], A 4'[t]==k 34 A 3[t] -(k 43 +k 45)A 4[t]+k 54 (A 1 o-A 1[t]-A 2[t]-A 3[t]-A 4[t]), A 1[0]==A 1 o, A 2[0]==0, A 3[0]==0, A 4[0]==0}, {A 1[t], A 2[t], A 3[t], A 4[t]}, {t, 0, tf}]; Plotting the transient solution using Mathematica’s functions Plot : Plot[Evaluate[Table[{A 1[t], A 2[t], A 3[t], A 4[t], A 1 o-A 1[t]-A 2[t]-A 3[t]-A 4[t]} /. sol 5]], {t, 0, tf}, Frame->True, Default. Font->{"Symbol-Bold", 14}, Frame. Label->{"t", "A 1, A 2, A 3, A 4, A 5"}, Plot. Range->{{0, tf}, {0, 1}}, Plot. Style->{RGBColor[1, 0, 0], RGBColor[0, 1, 0], RGBColor[0, 0, 1], RGBColor[1, 1, 0], RGBColor[0, 1, 1]}];

Steady state solution obtained using Mathematica’s function Solve : soleq=Solve[{0==-k 12 A 1+k 21 A 2, 0==k 12 A 1 -(k 21 +k 23)A 2+k 32 A 3, 0==k 23 A 2 -(k 32 +k 34)A 3+k 43 A 4, 0==k 34 A 3 -(k 43 +k 45)A 4+k 54 (A 1 o-A 1 -A 2 -A 3 -A 4)}, {A 1, A 2, A 3, A 4}] A 5=(A 1 o-A 1 -A 2 -A 3 -A 4)/. soleq Plotting the steady state solution using Mathematica’s functions List. Plot : List. Plot[Flatten[{A 1, A 2, A 3, A 4, (A 1 o-A 1 -A 2 -A 3 -A 4)}/. soleq], Plot. Style->{Point. Size[0. 015], RGBColor[1, 0, 0]}]

transient solution : Steady state solution :

Lotka-Volterra Mechanism Foxes and rabbits interactions : Governing equations :

NDSolve finds the solutions to the ODEs and Plot gives the figure with typical oscillationsfor the case A=3. 7, k 1=1. 2, k 2=1. 5 and k 3=1. 2 : A =3. 7; k 1=1. 2; k 2=1. 5; k 3=1. 2; LV=NDSolve[{X'[t]==k 1 A X[t]-k 2 X[t] Y[t], Y'[t]==k 2 X[t] Y[t] k 3 Y[t], X[0]==. 85, Y[0]==3. 2}, {X[t], Y[t]}, {t, 0, 10}]; Plot[Evaluate[Table[{X[t], Y[t]}/. LV]], {t, 0, 10}, Plot. Style->{RGBColor[1, 0, 0], RGBColor[0, 0, 1]}] x, y t

Oregonator model of the BZ reaction Main chemical reactions taking places : Governing equations :

NDSolve finds the solutions to the ODEs for the case s=100, f=1. 1, q=10 -6 and w=3. 835 : s=100; f=1. 1; q=10^-6; w=3. 835; sol 1=NDSolve[{x'[t]==s (x[t]+y[t]-x[t] y[t]-q x[t]^2), y'[t]==1/s (-y[t]-x[t] y[t]+f z[t]), z'[t]==w (x[t]-z[t]), x[0]==1, y[0]==1, z[0]==1}, {x, y, z}, {t, 0, 5000}, Working. Precision->25, Accuracy. Goal->10, Precision. Goal->10, Max. Steps->Infinity] Plotting the solution using Mathematica’s functions Plot and Parametric. Plot : pl 1=Plot[Evaluate[x[t]/. sol 1], {t, 0, 1000}, Plot. Range->All, Plot. Style->RGBColor[0, 0, 1]] Parametric. Plot[Evaluate[{z[t], y[t]}/. sol 1], {t, 500, 1000}, Plot. Range->{1, 1. 20}, Plot. Style->RGBColor[0, 1, 0]]

The solution shows regular oscillations such as those observed in the Belousov-Zhabotinski experiments. Limit cycle are obtained and a lapse of time is necessary before oscillations are observed. u 2 u 1 u 3 t

Lindermann-Hinshelwood Mechanism Quasi steady state approximation : Solve[{RA==k 1 A^2 -k 2 A Ac, 0==k 1 A^2 -k 2 A Ac-k 3 Ac}, {RA, Ac}]//Simplify If A large, reaction rate law is first order If A small, reaction rate law is second order

Continuous-Stirred Tank Reactor Ain , Bin d V A, B, C

k 1=1; k 2=10^-2; k 0=d/V; V=10; d=0. 6; sol=NDSolve[{ c'[t]==k 1 a[t] b[t]-k 2 c[t]-k 0 c[t], a'[t]==-k 1 a[t] b[t]+k 2 c[t]+k 0 (0. 01 -a[t]), b'[t]==-k 1 a[t] b[t]+k 2 c[t]+k 0 (0. 02 -b[t]), a[0]==10^-2, b[0]==2 10^-2, c[0]==0}, {a, b, c}, {t, 0, 300}] Plot[Evaluate[1000 c[t]/. sol], {t, 0, 300}, Plot. Range->{0, 8}, Plot. Style->RGBColor[1, 0, 0]] 103 C NDSolve and Plot are used to get the concentration profile of the product : t

Conclusion Mathematica’s graphical algebraic, capabilities can numerical be put and into advantage to solve several Physical Chemistry problems such as applied thermodynamics and chemical kinetics.

1/ 2/