Changeset 363

Timestamp:

04/17/08 18:37:03 (7 months ago)

Author:

powell

Message:

Implementations (untested) of the free energy derivatives.

Files:

• trunk/matml/src/ternary/freenergy.c (modified) (16 diffs)

trunk/matml/src/ternary/freenergy.c

@@ -12 +12 @@
 /*++++++++++++++++++++++++++++++++++++++
   This local function calculates the free energy as a function of composition,
-  temperature, pressure, and the specified parameters.
+  temperature, pressure, and the specified parameters.  Note: if you change
+  this function, be sure you update all of the derivative functions as well!
 
   inline double _G It returns the free energy.
 
-  double C2 Concentration of species 2.
-
-  double C3 Concentration of species 3.
+  double C1 Concentration of species 1.
+
+  double C2 Concentration of species 2.
+
+  double C3 Concentration of species 3.
+
+  double G1 Free energy of pure species 1.
+
+  double G2 Free energy of pure species 2.
+
+  double G3 Free energy of pure species 3.
 
   double T Temperature.
@@ -28 +37 @@
 
 static inline double _G
-(double C2, double C3, double T, double P, energy_params *eparams)
-{
-  double C1 = (1.-C2-C3 < 0.) ? 0. : 1.-C2-C3,
-    G1 = eparams->G1_T0 + eparams->G1_C*(T-eparams->T0),
-    G2 = eparams->G2_T0 + eparams->G2_C*(T-eparams->T0),
-    G3 = eparams->G3_T0 + eparams->G3_C*(T-eparams->T0);
-
+(double C1, double C2, double C3, double G1, double G2, double G3,
+ double T, double P, energy_params *eparams)
+{
   /*+ This uses the Flory-Huggins free energy formula for a polymer solution,
     which divides each species' entropy term by its molar volume and adds
@@ -49 +54 @@
     +latex+\end{eqnarray}
     +*/
-  return G1 + (G2-G1)*C2 + (G3-G2)*C3
-    + eparams->R*T*(((C1<=0.)?0.:C1/eparams->M1*log(C1))
-                    + ((C2<=0.)?0.:C2/eparams->M2*log(C2))
-                    + ((C3<=0.)?0.:C3/eparams->M3*log(C3)))
-    + eparams->Omega12*C2*C1 + eparams->Omega13*C3*C1 + eparams->Omega23*C2*C3
+  return G1 + (G2-G1)*C2 + (G3-G1)*C3
+    + eparams->R*T*(((C1<=0.) ? 0. : C1/eparams->M1*log(C1))
+                    + ((C2<=0.) ? 0. : C2/eparams->M2*log(C2))
+                    + ((C3<=0.) ? 0. : C3/eparams->M3*log(C3)))
+    + eparams->Omega12*C1*C2 + eparams->Omega13*C1*C3 + eparams->Omega23*C2*C3
     + eparams->Omega123*C1*C2*C3;
 }
@@ -67 +72 @@
   inline double _G2 It returns the free energy derivative.
 
-  double C2 Concentration of species 2.
-
-  double C3 Concentration of species 3.
+  double C1 Concentration of species 1.
+
+  double C2 Concentration of species 2.
+
+  double C3 Concentration of species 3.
+
+  double G1 Free energy of pure species 1.
+
+  double G2 Free energy of pure species 2.
+
+  double G3 Free energy of pure species 3.
 
   double T Temperature.
@@ -79 +92 @@
 
 static inline double _G2
-(double C2, double C3, double T, double P, energy_params *eparams)
-{
-  return 0.;
+(double C1, double C2, double C3, double G1, double G2, double G3,
+ double T, double P, energy_params *eparams)
+{
+  /*+
+    +latex+The derivative of the free energy in equation \ref{eq:flory-huggins}
+    +latex+with respect to $C_2$ is:
+    +latex+\begin{equation}
+    +latex+  \label{eq:flory-huggins-C2}
+    +latex+\end{equation}
+    +*/
+  return G2-G1
+    + eparams->R*T*(((C2+C3>=1) ? INFINITY : (-log(C1)-1)/eparams->M1)
+                    + ((C2<=0.) ?-INFINITY : (log(C2)+1)/eparams->M2))
+    + eparams->Omega12*(C1-C2) - eparams->Omega13*C3 + eparams->Omega23*C3
+    + eparams->Omega123*(C1-C2)*C3;
 }
 
@@ -94 +119 @@
   inline double _G3 It returns the free energy derivative.
 
-  double C2 Concentration of species 2.
-
-  double C3 Concentration of species 3.
+  double C1 Concentration of species 1.
+
+  double C2 Concentration of species 2.
+
+  double C3 Concentration of species 3.
+
+  double G1 Free energy of pure species 1.
+
+  double G2 Free energy of pure species 2.
+
+  double G3 Free energy of pure species 3.
 
   double T Temperature.
@@ -106 +139 @@
 
 static inline double _G3
-(double C2, double C3, double T, double P, energy_params *eparams)
-{
-  return 0.;
+(double C1, double C2, double C3, double G1, double G2, double G3,
+ double T, double P, energy_params *eparams)
+{
+  return G3-G1
+    + eparams->R*T*(((C2+C3>=1) ? INFINITY : (-log(C1)-1)/eparams->M1)
+                    + ((C3<=0.) ?-INFINITY : (log(C3+1)/eparams->M3)))
+    - eparams->Omega12*C2 + eparams->Omega13*(C1-C3) + eparams->Omega23*C2
+    + eparams->Omega123*(C1-C3)*C2;
 }
 
@@ -123 +161 @@
   inline double _G22 It returns the free energy derivative.
 
-  double C2 Concentration of species 2.
-
-  double C3 Concentration of species 3.
+  double C1 Concentration of species 1.
+
+  double C2 Concentration of species 2.
+
+  double C3 Concentration of species 3.
+
+  double G1 Free energy of pure species 1.
+
+  double G2 Free energy of pure species 2.
+
+  double G3 Free energy of pure species 3.
 
   double T Temperature.
@@ -135 +181 @@
 
 static inline double _G22
-(double C2, double C3, double T, double P, energy_params *eparams)
-{
-  return 0.;
+(double C1, double C2, double C3, double G1, double G2, double G3,
+ double T, double P, energy_params *eparams)
+{
+  return
+    eparams->R*T*(((C2+C3>=1) ? INFINITY : 1./C1/eparams->M1)
+                  + ((C2<=0.) ? INFINITY : 1./C2/eparams->M2))
+    - eparams->Omega12*2 - eparams->Omega123*2*C3;
 }
 
@@ -152 +202 @@
   inline double _G33 It returns the free energy derivative.
 
-  double C2 Concentration of species 2.
-
-  double C3 Concentration of species 3.
+  double C1 Concentration of species 1.
+
+  double C2 Concentration of species 2.
+
+  double C3 Concentration of species 3.
+
+  double G1 Free energy of pure species 1.
+
+  double G2 Free energy of pure species 2.
+
+  double G3 Free energy of pure species 3.
 
   double T Temperature.
@@ -164 +222 @@
 
 static inline double _G33
-(double C2, double C3, double T, double P, energy_params *eparams)
-{
-  return 0.;
+(double C1, double C2, double C3, double G1, double G2, double G3,
+ double T, double P, energy_params *eparams)
+{
+  return
+    eparams->R*T*(((C2+C3>=1) ? INFINITY : 1./C1/eparams->M1)
+                  + ((C3<=0.) ? INFINITY : 1./C3/eparams->M3))
+    - eparams->Omega13*2 - eparams->Omega123*2*C2;
 }
 
@@ -181 +243 @@
   inline double _G23 It returns the free energy derivative.
 
-  double C2 Concentration of species 2.
-
-  double C3 Concentration of species 3.
+  double C1 Concentration of species 1.
+
+  double C2 Concentration of species 2.
+
+  double C3 Concentration of species 3.
+
+  double G1 Free energy of pure species 1.
+
+  double G2 Free energy of pure species 2.
+
+  double G3 Free energy of pure species 3.
 
   double T Temperature.
@@ -193 +263 @@
 
 static inline double _G23
-(double C2, double C3, double T, double P, energy_params *eparams)
-{
-  return 0.;
+(double C1, double C2, double C3, double G1, double G2, double G3,
+ double T, double P, energy_params *eparams)
+{
+  return
+    eparams->R*T*((C2+C3>=1) ? INFINITY : 1./C1/eparams->M1)
+    - eparams->Omega12 - eparams->Omega13 + eparams->Omega23
+    + eparams->Omega123*(2*C1-1);
 }
 
@@ -218 +292 @@
 (double C2, double C3, double T, double P, energy_params *eparams)
 {
-  return _G (C2, C3, T, P, eparams);
+  double C1 = (1-C2-C3 < 0.) ? 0. : 1-C2-C3,
+    G1 = eparams->G1_T0 + eparams->G1_C*(T-eparams->T0),
+    G2 = eparams->G2_T0 + eparams->G2_C*(T-eparams->T0),
+    G3 = eparams->G3_T0 + eparams->G3_C*(T-eparams->T0);
+
+  return _G (C1,C2,C3, G1,G2,G3, T, P, eparams);
 }
 
@@ -244 +323 @@
 {
   int i;
+  double C1,
+    G1 = eparams->G1_T0 + eparams->G1_C*(T-eparams->T0),
+    G2 = eparams->G2_T0 + eparams->G2_C*(T-eparams->T0),
+    G3 = eparams->G3_T0 + eparams->G3_C*(T-eparams->T0);
 
   for (i=0; i<n; i++)
-    points[i].G = _G (points[i].C2, points[i].C3, T, P, eparams);
+    {
+      C1 = 1.-points[i].C2-points[i].C3;
+      points[i].G = _G (C1,points[i].C2,points[i].C3, G1,G2,G3, T, P, eparams);
+    }
 
   return 0;
@@ -274 +360 @@
 {
   int i;
+  double C1,
+    G1 = eparams->G1_T0 + eparams->G1_C*(T-eparams->T0),
+    G2 = eparams->G2_T0 + eparams->G2_C*(T-eparams->T0),
+    G3 = eparams->G3_T0 + eparams->G3_C*(T-eparams->T0);
 
   for (i=0; i<n; i++)
     {
-      points[i].G   = _G   (points[i].C2, points[i].C3, T, P, eparams);
-      points[i].G2  = _G2  (points[i].C2, points[i].C3, T, P, eparams);
-      points[i].G3  = _G3  (points[i].C2, points[i].C3, T, P, eparams);
-      points[i].G22 = _G22 (points[i].C2, points[i].C3, T, P, eparams);
-      points[i].G33 = _G33 (points[i].C2, points[i].C3, T, P, eparams);
-      points[i].G23 = _G23 (points[i].C2, points[i].C3, T, P, eparams);
+      C1 = 1.-points[i].C2-points[i].C3;
+      points[i].G  =_G  (C1,points[i].C2,points[i].C3, G1,G2,G3, T,P, eparams);
+      points[i].G2 =_G2 (C1,points[i].C2,points[i].C3, G1,G2,G3, T,P, eparams);
+      points[i].G3 =_G3 (C1,points[i].C2,points[i].C3, G1,G2,G3, T,P, eparams);
+      points[i].G22=_G22(C1,points[i].C2,points[i].C3, G1,G2,G3, T,P, eparams);
+      points[i].G33=_G33(C1,points[i].C2,points[i].C3, G1,G2,G3, T,P, eparams);
+      points[i].G23=_G23(C1,points[i].C2,points[i].C3, G1,G2,G3, T,P, eparams);
     }