Changeset 449

Show
Ignore:
Timestamp:
03/18/2009 09:10:18 PM (3 years ago)
Author:
powell
Message:

Add unlimited number of Gaussian terms to free energy.

Location:
trunk/matml/src/ternary
Files:
4 modified

Legend:

Unmodified
Added
Removed

  • trunk/matml/src/ternary/freenergy.c

    r447 r449  
    4848 double T, double P, energy_params *eparams) 
    4949{ 
    50   /*+ This uses the Flory-Huggins free energy formula for a polymer solution, 
    51     which divides each species' entropy term by its molar volume and adds 
    52     binary interaction parameters. 
     50  /*+ This uses a free energy formula with a coefficient for each species' 
     51    entropy term, which allows one to use the Flory-Huggins model by setting 
     52    that coefficient equal to the reciprocal of the molar volume.  To this it 
     53    adds binary interaction parameters. 
    5354    +latex+The resulting regular solution expression for free energy is: 
    5455    +latex+\begin{eqnarray} 
    5556    +latex+  \label{eq:flory-huggins} 
    5657    +latex+  G &=& G_1 + (G_2-G_1) C_2 + (G_3-G_1) C_3 \\\ \nonumber 
    57     +latex+  &+& RT \left(\frac{1}{m_1}C_1\log C_1 
    58     +latex+    + \frac{1}{m_2}C_2\log C_2 
    59     +latex+    + \frac{1}{m_3}C_3\log C_3\right) \\\ \nonumber 
     58    +latex+  &+& RT \left(S_1C_1\log C_1 
     59    +latex+    + S_2C_2\log C_2 
     60    +latex+    + S_3C_3\log C_3\right) \\\ \nonumber 
    6061    +latex+  &+& \Omega_{12}C_1C_2 + \Omega_{13}C_1C_3 
    6162    +latex+  + \Omega_{23}C_2C_3 +\Omega_{123}C_1C_2C_3. 
    6263    +latex+\end{eqnarray} 
    6364    +*/ 
    64   return G1 + (G2-G1)*C2 + (G3-G1)*C3 
     65  int i; 
     66  double ret = G1 + (G2-G1)*C2 + (G3-G1)*C3 
    6567    + eparams->R*T*(  ((C1<=0.) ? 0. : eparams->S1*C1*log(C1)) 
    6668                    + ((C2<=0.) ? 0. : eparams->S2*C2*log(C2)) 
    … …  
    7072    + eparams->Omega12*C1*C2 + eparams->Omega13*C1*C3 + eparams->Omega23*C2*C3 
    7173    + eparams->Omega123*C1*C2*C3; 
     74 
     75  for (i=0; i<eparams->n_gauss; i++) 
     76    ret += (eparams->gauss) [i].G * 
     77      exp (-((C2-(eparams->gauss) [i].C2)*(C2-(eparams->gauss) [i].C2) + 
     78             (C3-(eparams->gauss) [i].C3)*(C3-(eparams->gauss) [i].C3)) / 
     79           (eparams->gauss) [i].w / (eparams->gauss) [i].w); 
     80 
     81  return ret; 
    7282} 
    7383 
    … …  
    115125    +latex+\end{equation} 
    116126    +*/ 
    117   return G2-G1 
     127  int i; 
     128  double ret = G2-G1 
    118129    + eparams->R*T*(((eparams->S1==0.) ? 0. : 
    119130                     ((C1<=0.) ? LOCAL_INFINITY : eparams->S1*(-log(C1)-1))) 
    … …  
    124135    + eparams->Omega12*(C1-C2) - eparams->Omega13*C3 + eparams->Omega23*C3 
    125136    + eparams->Omega123*(C1-C2)*C3; 
     137 
     138  for (i=0; i<eparams->n_gauss; i++) 
     139    ret -= 2. * (eparams->gauss) [i].G * (C2-(eparams->gauss) [i].C2) 
     140      / (eparams->gauss) [i].w / (eparams->gauss) [i].w 
     141      * exp (-((C2-(eparams->gauss) [i].C2)*(C2-(eparams->gauss) [i].C2) + 
     142               (C3-(eparams->gauss) [i].C3)*(C3-(eparams->gauss) [i].C3)) / 
     143             (eparams->gauss) [i].w / (eparams->gauss) [i].w); 
     144 
     145  return ret; 
    126146} 
    127147 
    … …  
    162182 double T, double P, energy_params *eparams, free_energy_flags flags) 
    163183{ 
    164   return G3-G1 
     184  int i; 
     185  double ret = G3-G1 
    165186    + eparams->R*T*(((eparams->S1==0.) ? 0. : 
    166187                     ((C2+C3>=1) ? LOCAL_INFINITY : eparams->S1*(-log(C1)-1))) 
    … …  
    171192    - eparams->Omega12*C2 + eparams->Omega13*(C1-C3) + eparams->Omega23*C2 
    172193    + eparams->Omega123*(C1-C3)*C2; 
     194 
     195  for (i=0; i<eparams->n_gauss; i++) 
     196    ret -= 2. * (eparams->gauss) [i].G * (C3-(eparams->gauss) [i].C3) 
     197      / (eparams->gauss) [i].w / (eparams->gauss) [i].w 
     198      * exp (-((C2-(eparams->gauss) [i].C2)*(C2-(eparams->gauss) [i].C2) + 
     199               (C3-(eparams->gauss) [i].C3)*(C3-(eparams->gauss) [i].C3)) / 
     200             (eparams->gauss) [i].w / (eparams->gauss) [i].w); 
     201 
     202  return ret; 
    173203} 
    174204 
    … …  
    211241 double T, double P, energy_params *eparams, free_energy_flags flags) 
    212242{ 
    213   return 
     243  int i; 
     244  double ret = 
    214245    eparams->R*T*(((eparams->S1==0.) ? 0. : 
    215246                   ((C2+C3>=1) ? LOCAL_INFINITY : eparams->S1/C1)) 
    … …  
    217248                  + ((C2>=1.) ? LOCAL_INFINITY : eparams->S4/(1.-C2))) 
    218249    - eparams->Omega12*2 - eparams->Omega123*2*C3; 
     250 
     251  for (i=0; i<eparams->n_gauss; i++) 
     252    ret += 2. * (eparams->gauss) [i].G 
     253      / (eparams->gauss) [i].w / (eparams->gauss) [i].w 
     254      * (2. * (C2-(eparams->gauss) [i].C2) * (C2-(eparams->gauss) [i].C2) 
     255         / (eparams->gauss) [i].w / (eparams->gauss) [i].w - 1.) 
     256      * exp (-((C2-(eparams->gauss) [i].C2)*(C2-(eparams->gauss) [i].C2) + 
     257               (C3-(eparams->gauss) [i].C3)*(C3-(eparams->gauss) [i].C3)) / 
     258             (eparams->gauss) [i].w / (eparams->gauss) [i].w); 
     259 
     260  return ret; 
    219261} 
    220262 
    … …  
    257299 double T, double P, energy_params *eparams, free_energy_flags flags) 
    258300{ 
    259   return 
     301  int i; 
     302  double ret = 
    260303    eparams->R*T*(((eparams->S1==0.) ? 0. : 
    261304                   ((C2+C3>=1) ? LOCAL_INFINITY : eparams->S1/C1)) 
    … …  
    263306                  + ((C3>=1.) ? LOCAL_INFINITY : eparams->S5/(1.-C3))) 
    264307    - eparams->Omega13*2 - eparams->Omega123*2*C2; 
     308 
     309  for (i=0; i<eparams->n_gauss; i++) 
     310    ret += 2. * (eparams->gauss) [i].G 
     311      / (eparams->gauss) [i].w / (eparams->gauss) [i].w 
     312      * (2. * (C3-(eparams->gauss) [i].C3) * (C3-(eparams->gauss) [i].C3) 
     313         / (eparams->gauss) [i].w / (eparams->gauss) [i].w - 1.) 
     314      * exp (-((C2-(eparams->gauss) [i].C2)*(C2-(eparams->gauss) [i].C2) + 
     315               (C3-(eparams->gauss) [i].C3)*(C3-(eparams->gauss) [i].C3)) / 
     316             (eparams->gauss) [i].w / (eparams->gauss) [i].w); 
     317 
     318  return ret; 
    265319} 
    266320 
    … …  
    303357 double T, double P, energy_params *eparams, free_energy_flags flags) 
    304358{ 
    305   return 
     359  int i; 
     360  double ret = 
    306361    eparams->R*T*((eparams->S1==0.) ? 0. : 
    307362                  ((C2+C3>=1) ? LOCAL_INFINITY : eparams->S1/C1)) 
    308363    - eparams->Omega12 - eparams->Omega13 + eparams->Omega23 
    309364    + eparams->Omega123*(2*C1-1); 
     365 
     366  for (i=0; i<eparams->n_gauss; i++) 
     367    ret += 4. * (eparams->gauss) [i].G 
     368      * (C2-(eparams->gauss) [i].C2) * (C3-(eparams->gauss) [i].C3) 
     369      / (eparams->gauss) [i].w / (eparams->gauss) [i].w 
     370      / (eparams->gauss) [i].w / (eparams->gauss) [i].w 
     371      * exp (-((C2-(eparams->gauss) [i].C2)*(C2-(eparams->gauss) [i].C2) + 
     372               (C3-(eparams->gauss) [i].C3)*(C3-(eparams->gauss) [i].C3)) / 
     373             (eparams->gauss) [i].w / (eparams->gauss) [i].w); 
     374 
     375  return ret; 
    310376} 
    311377 

  • trunk/matml/src/ternary/square.c

    r448 r449  
    3737  /* eparams: R,T0, G1@T0,G2,G3, C1,C2,C3, M1,M2,M3, O12,O13,O23,O123 */ 
    3838  energy_params /* Solid, liquid */ 
    39     eparams0 = {"Solid", 1.,1., 0.,-.1,-.2, -1.,-1.1,-1.2, 0.,1.,1.,1.,1., -.5,2.,.2, 0.}, 
    40     eparams1 = {"Liquid", 1.,1., .3,.1,-.1, -2.2,-2.,-1.8, 0.,1.,.67,1.,2., -.3,-.4,-.1, -.5}, 
     39    eparams0 = {"Solid", 1.,1., 0.,-.1,-.2, -1.,-1.1,-1.2, 0.,1.,1.,1.,1., -.5,2.,.2, 0., NULL,0}, 
     40    eparams1 = {"Liquid", 1.,1., .3,.1,-.1, -2.2,-2.,-1.8, 0.,1.,.67,1.,2., -.3,-.4,-.1, -.5, NULL,0}, 
    4141      allparams[2] = {eparams0, eparams1}; 
    4242 

  • trunk/matml/src/ternary/ternary.c

    r440 r449  
    3737  /* eparams: R,T0, G1@T0,G2,G3, C1,C2,C3, M1,M2,M3, O12,O13,O23,O123 */ 
    3838  energy_params /* Solid, liquid */ 
    39     eparams0 = {"Solid", 1.,1., 0.,-.1,-.2, -1.,-1.1,-1.2, 1.,.5,.2,0.,0., -.5,2.,.2, 0.}, 
    40     eparams1 = {"Liquid", 1.,1., .3,.1,-.1, -2.2,-2.,-1.8, 1.,1.,.67,0.,0., -.3,-.4,-.1, -.5}, 
     39    eparams0 = {"Solid", 1.,1., 0.,-.1,-.2, -1.,-1.1,-1.2, 1.,.5,.2,0.,0., -.5,2.,.2, 0., NULL,0}, 
     40    eparams1 = {"Liquid", 1.,1., .3,.1,-.1, -2.2,-2.,-1.8, 1.,1.,.67,0.,0., -.3,-.4,-.1, -.5, NULL,0}, 
    4141      allparams[2] = {eparams0, eparams1}; 
    4242 

  • trunk/matml/src/ternary/ternary.h

    r444 r449  
    7676                   +*/ 
    7777} ternary_point; /*+ Ternary "point" structure. +*/ 
     78 
     79typedef struct { 
     80  double C2; /*+ Composition at center of this Gaussian term +*/ 
     81  double C3; /*+ Composition at center of this Gaussian term +*/ 
     82  double G;  /*+ Free energy addition at the Gaussian center +*/ 
     83  double w;  /*+ Width of the Gaussian distribution +*/ 
     84  double n;  /*+ Exponent of the distribution, currently fixed at 2 +*/ 
     85} energy_gaussian; 
    7886 
    7987typedef struct { 
    … …  
    96104  double Omega23;  /*+ Species 2-3 regular solution interaction parameter +*/ 
    97105  double Omega123; /*+ Species 1-2-3 regular solution interaction parameter +*/ 
     106  energy_gaussian *gauss; /*+ Array of Gaussian energy parameter structs +*/ 
     107  int n_gauss;     /*+ Number of Gaussian energy terms +*/ 
    98108} energy_params;   /*+ Free energy parameters structure +*/ 
    99109