Cantera  2.0
vcs_species_thermo.cpp
Go to the documentation of this file.
1 /**
2  * @file vcs_species_thermo.cpp
3  * Implementation for the VCS_SPECIES_THERMO object.
4  */
5 /*
6  * Copyright (2005) Sandia Corporation. Under the terms of
7  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
8  * U.S. Government retains certain rights in this software.
9  */
10 
11 
12 
13 #include "cantera/equil/vcs_solve.h"
14 #include "vcs_species_thermo.h"
15 #include "cantera/equil/vcs_defs.h"
16 #include "cantera/equil/vcs_VolPhase.h"
17 
18 #include "vcs_Exception.h"
19 #include "cantera/equil/vcs_internal.h"
20 
21 #include <cstdio>
22 #include <cstdlib>
23 #include <cmath>
24 
25 using namespace std;
26 
27 namespace VCSnonideal
28 {
29 
30 
31 VCS_SPECIES_THERMO::VCS_SPECIES_THERMO(size_t indexPhase,
32  size_t indexSpeciesPhase) :
33 
34  IndexPhase(indexPhase),
35  IndexSpeciesPhase(indexSpeciesPhase),
36  OwningPhase(0),
37  SS0_Model(VCS_SS0_CONSTANT),
38  SS0_feSave(0.0),
39  SS0_TSave(-90.0),
40  SS0_T0(273.15),
41  SS0_H0(0.0),
42  SS0_S0(0.0),
43  SS0_Cp0(0.0),
44  SS0_Pref(1.01325E5),
45  SS0_Params(0),
46  SSStar_Model(VCS_SSSTAR_CONSTANT),
47  SSStar_Params(0),
48  Activity_Coeff_Model(VCS_AC_CONSTANT),
49  Activity_Coeff_Params(0),
50  SSStar_Vol_Model(VCS_SSVOL_IDEALGAS),
51  SSStar_Vol_Params(0),
52  SSStar_Vol0(-1.0),
53  UseCanteraCalls(false),
54  m_VCS_UnitsFormat(VCS_UNITS_UNITLESS)
55 {
56  SS0_Pref = 1.01325E5;
57 }
58 
59 
60 /******************************************************************************
61  *
62  * destructor
63  */
64 VCS_SPECIES_THERMO::~VCS_SPECIES_THERMO()
65 {
66 }
67 
68 /*****************************************************************************
69  *
70  * Copy Constructor VCS_SPECIES_THERMO
71  */
72 VCS_SPECIES_THERMO::VCS_SPECIES_THERMO(const VCS_SPECIES_THERMO& b) :
73  IndexPhase(b.IndexPhase),
74  IndexSpeciesPhase(b.IndexSpeciesPhase),
75  OwningPhase(b.OwningPhase),
76  SS0_Model(b.SS0_Model),
77  SS0_feSave(b.SS0_feSave),
78  SS0_TSave(b.SS0_TSave),
79  SS0_T0(b.SS0_T0),
80  SS0_H0(b.SS0_H0),
81  SS0_S0(b.SS0_S0),
82  SS0_Cp0(b.SS0_Cp0),
83  SS0_Pref(b.SS0_Pref),
84  SS0_Params(0),
85  SSStar_Model(b.SSStar_Model),
86  SSStar_Params(0),
87  Activity_Coeff_Model(b.Activity_Coeff_Model),
88  Activity_Coeff_Params(0),
89  SSStar_Vol_Model(b.SSStar_Vol_Model),
90  SSStar_Vol_Params(0),
91  SSStar_Vol0(b.SSStar_Vol0),
92  UseCanteraCalls(b.UseCanteraCalls),
93  m_VCS_UnitsFormat(b.m_VCS_UnitsFormat)
94 {
95 
96  SS0_Params = 0;
97 }
98 
99 /*****************************************************************************
100  *
101  * Assignment operator for VCS_SPECIES_THERMO
102  */
103 VCS_SPECIES_THERMO&
104 VCS_SPECIES_THERMO::operator=(const VCS_SPECIES_THERMO& b)
105 {
106  if (&b != this) {
107  IndexPhase = b.IndexPhase;
108  IndexSpeciesPhase = b.IndexSpeciesPhase;
109  OwningPhase = b.OwningPhase;
110  SS0_Model = b.SS0_Model;
111  SS0_feSave = b.SS0_feSave;
112  SS0_TSave = b.SS0_TSave;
113  SS0_T0 = b.SS0_T0;
114  SS0_H0 = b.SS0_H0;
115  SS0_S0 = b.SS0_S0;
116  SS0_Cp0 = b.SS0_Cp0;
117  SS0_Pref = b.SS0_Pref;
118  SSStar_Model = b.SSStar_Model;
119  /*
120  * shallow copy because function is undeveloped.
121  */
122  SSStar_Params = b.SSStar_Params;
123  Activity_Coeff_Model = b.Activity_Coeff_Model;
124  /*
125  * shallow copy because function is undeveloped.
126  */
127  Activity_Coeff_Params = b.Activity_Coeff_Params;
128  SSStar_Vol_Model = b.SSStar_Vol_Model;
129  /*
130  * shallow copy because function is undeveloped.
131  */
132  SSStar_Vol_Params = b.SSStar_Vol_Params;
133  SSStar_Vol0 = b.SSStar_Vol0;
134  UseCanteraCalls = b.UseCanteraCalls;
135  m_VCS_UnitsFormat = b.m_VCS_UnitsFormat;
136  }
137  return *this;
138 }
139 
140 /******************************************************************************
141  *
142  * duplMyselfAsVCS_SPECIES_THERMO(): (virtual)
143  *
144  * This routine can duplicate inherited objects given a base class
145  * pointer. It relies on valid copy constructors.
146  */
147 
148 VCS_SPECIES_THERMO* VCS_SPECIES_THERMO::duplMyselfAsVCS_SPECIES_THERMO()
149 {
150  VCS_SPECIES_THERMO* ptr = new VCS_SPECIES_THERMO(*this);
151  return ptr;
152 }
153 
154 
155 /**************************************************************************
156  *
157  * GStar_R_calc();
158  *
159  * This function calculates the standard state Gibbs free energy
160  * for species, kspec, at the solution temperature TKelvin and
161  * solution pressure, Pres.
162  *
163  *
164  * Input
165  * kglob = species global index.
166  * TKelvin = Temperature in Kelvin
167  * pres = pressure is given in units specified by if__ variable.
168  *
169  *
170  * Output
171  * return value = standard state free energy in units of Kelvin.
172  */
173 double VCS_SPECIES_THERMO::GStar_R_calc(size_t kglob, double TKelvin,
174  double pres)
175 {
176  char yo[] = "VCS_SPECIES_THERMO::GStar_R_calc ";
177  double fe, T;
178  fe = G0_R_calc(kglob, TKelvin);
179  T = TKelvin;
180  if (UseCanteraCalls) {
181  AssertThrowVCS(m_VCS_UnitsFormat == VCS_UNITS_MKS, "Possible inconsistency");
182  size_t kspec = IndexSpeciesPhase;
183  OwningPhase->setState_TP(TKelvin, pres);
184  fe = OwningPhase->GStar_calc_one(kspec);
185  double R = vcsUtil_gasConstant(m_VCS_UnitsFormat);
186  fe /= R;
187  } else {
188  double pref = SS0_Pref;
189  switch (SSStar_Model) {
190  case VCS_SSSTAR_CONSTANT:
191  break;
192  case VCS_SSSTAR_IDEAL_GAS:
193  fe += T * log(pres/ pref);
194  break;
195  default:
196  plogf("%sERROR: unknown SSStar model\n", yo);
197  exit(EXIT_FAILURE);
198  }
199  }
200  return fe;
201 }
202 
203 /**************************************************************************
204  *
205  * VolStar_calc:
206  *
207  * This function calculates the standard state molar volume
208  * for species, kspec, at the temperature TKelvin and pressure, Pres,
209  *
210  * Input
211  *
212  * Output
213  * return value = standard state volume in m**3 per kmol.
214  * (VCS_UNITS_MKS)
215  */
216 double VCS_SPECIES_THERMO::
217 VolStar_calc(size_t kglob, double TKelvin, double presPA)
218 {
219  char yo[] = "VCS_SPECIES_THERMO::VStar_calc ";
220  double vol, T;
221 
222  T = TKelvin;
223  if (UseCanteraCalls) {
224  AssertThrowVCS(m_VCS_UnitsFormat == VCS_UNITS_MKS, "Possible inconsistency");
225  size_t kspec = IndexSpeciesPhase;
226  OwningPhase->setState_TP(TKelvin, presPA);
227  vol = OwningPhase->VolStar_calc_one(kspec);
228  } else {
229  switch (SSStar_Vol_Model) {
230  case VCS_SSVOL_CONSTANT:
231  vol = SSStar_Vol0;
232  break;
233  case VCS_SSVOL_IDEALGAS:
234  vol= Cantera::GasConstant * T / presPA;
235  break;
236  default:
237  plogf("%sERROR: unknown SSVol model\n", yo);
238  exit(EXIT_FAILURE);
239  }
240  }
241  return vol;
242 }
243 
244 /**************************************************************************
245  *
246  * G0_R_calc:
247  *
248  * This function calculates the naught state Gibbs free energy
249  * for species, kspec, at the temperature TKelvin
250  *
251  * Input
252  * kglob = species global index.
253  * TKelvin = Temperature in Kelvin
254  *
255  * Output
256  * return value = naught state free energy in Kelvin.
257  */
258 double VCS_SPECIES_THERMO::G0_R_calc(size_t kglob, double TKelvin)
259 {
260 #ifdef DEBUG_MODE
261  char yo[] = "VS_SPECIES_THERMO::G0_R_calc ";
262 #endif
263  double fe, H, S;
264  if (SS0_Model == VCS_SS0_CONSTANT) {
265  fe = SS0_feSave;
266  return fe;
267  }
268  if (TKelvin == SS0_TSave) {
269  fe = SS0_feSave;
270  return fe;
271  }
272  if (UseCanteraCalls) {
273  AssertThrowVCS(m_VCS_UnitsFormat == VCS_UNITS_MKS, "Possible inconsistency");
274  size_t kspec = IndexSpeciesPhase;
275  OwningPhase->setState_T(TKelvin);
276  fe = OwningPhase->G0_calc_one(kspec);
277  double R = vcsUtil_gasConstant(m_VCS_UnitsFormat);
278  fe /= R;
279  } else {
280  switch (SS0_Model) {
281  case VCS_SS0_CONSTANT:
282  fe = SS0_feSave;
283  break;
284  case VCS_SS0_CONSTANT_CP:
285  H = SS0_H0 + (TKelvin - SS0_T0) * SS0_Cp0;
286  S = SS0_Cp0 + SS0_Cp0 * log((TKelvin / SS0_T0));
287  fe = H - TKelvin * S;
288  break;
289  default:
290 #ifdef DEBUG_MODE
291  plogf("%sERROR: unknown model\n", yo);
292 #endif
293  exit(EXIT_FAILURE);
294  }
295  }
296  SS0_feSave = fe;
297  SS0_TSave = TKelvin;
298  return fe;
299 }
300 
301 /**************************************************************************
302  *
303  * eval_ac:
304  *
305  * This function evaluates the activity coefficient
306  * for species, kspec
307  *
308  * Input
309  * kglob -> integer value of the species in the global
310  * species list within VCS_GLOB. Phase and local species id
311  * can be looked up within object.
312  *
313  * Note, T, P and mole fractions are obtained from the
314  * single private instance of VCS_GLOB
315  *
316  *
317  * Output
318  * return value = activity coefficient for species kspec
319  */
320 double VCS_SPECIES_THERMO::eval_ac(size_t kglob)
321 {
322 #ifdef DEBUG_MODE
323  char yo[] = "VCS_SPECIES_THERMO::eval_ac ";
324 #endif
325  double ac;
326  /*
327  * Activity coefficients are frequently evaluated on a per phase
328  * basis. If they are, then the currPhAC[] boolean may be used
329  * to reduce repeated work. Just set currPhAC[iph], when the
330  * activity coefficients for all species in the phase are reevaluated.
331  */
332  if (UseCanteraCalls) {
333  size_t kspec = IndexSpeciesPhase;
334  ac = OwningPhase->AC_calc_one(kspec);
335  } else {
336  switch (Activity_Coeff_Model) {
337  case VCS_AC_CONSTANT:
338  ac = 1.0;
339  break;
340  default:
341 #ifdef DEBUG_MODE
342  plogf("%sERROR: unknown model\n", yo);
343 #endif
344  exit(EXIT_FAILURE);
345  }
346  }
347  return ac;
348 }
349 
350 /*****************************************************************************/
351 }
Generated by   doxygen 1.8.2