CoolProp  6.6.0
An open-source fluid property and humid air property database
Ice.cpp
Go to the documentation of this file.
1 
2 #ifndef __powerpc__
3 # include <math.h>
4 # include <complex>
5 # include <iostream>
6 static std::complex<double> t1(0.368017112855051e-1, 0.510878114959572e-1);
7 static std::complex<double> r1(0.447050716285388e2, 0.656876847463481e2);
8 static std::complex<double> t2(0.337315741065416, 0.335449415919309);
9 static std::complex<double> r20(-0.725974574329220e2, -0.781008427112870e2);
10 static std::complex<double> r21(-0.557107698030123e-4, 0.464578634580806e-4);
11 static std::complex<double> r22(0.234801409215913e-10, -0.285651142904972e-10);
12 #endif
13 
14 #include "Ice.h"
15 
16 static double T_t = 273.16,
17  p_t = 611.657,
18  p_0 = 101325;
19 
20 // Complex Constants for EOS
21 static double g00 = -0.632020233449497e6;
22 static double g01 = 0.655022213658955;
23 static double g02 = -0.189369929326131e-7;
24 static double g03 = 0.339746123271053e-14;
25 static double g04 = -0.556464869058991e-21;
26 static double s0 = -0.332733756492168e4;
27 
28 double IsothermCompress_Ice(double T, double p) {
29 #ifndef __powerpc__
30  // Inputs in K, Pa, Output in 1/Pa
31  return -dg2_dp2_Ice(T, p) / dg_dp_Ice(T, p);
32 #else
33  return 1e99;
34 #endif
35 }
36 double psub_Ice(double T) {
37 #ifndef __powerpc__
38  double a[] = {0, -0.212144006e2, 0.273203819e2, -0.610598130e1};
39  double b[] = {0, 0.333333333e-2, 0.120666667e1, 0.170333333e1};
40  double summer = 0, theta;
41  theta = T / T_t;
42  for (int i = 1; i <= 3; i++) {
43  summer += a[i] * pow(theta, b[i]);
44  }
45  return p_t * exp(1 / theta * summer);
46 #else
47  return 1e99;
48 #endif
49 }
50 
51 double g_Ice(double T, double p) {
52 #ifndef __powerpc__
53  std::complex<double> r2, term1, term2;
54  double g0, theta, pi, pi_0;
55  theta = T / T_t;
56  pi = p / p_t;
57  pi_0 = p_0 / p_t;
58  g0 = g00 * pow(pi - pi_0, 0.0) + g01 * pow(pi - pi_0, 1.0) + g02 * pow(pi - pi_0, 2.0) + g03 * pow(pi - pi_0, 3.0) + g04 * pow(pi - pi_0, 4.0);
59  r2 = r20 * pow(pi - pi_0, 0.0) + r21 * pow(pi - pi_0, 1.0) + r22 * pow(pi - pi_0, 2.0);
60  // The two terms of the summation
61  term1 = r1 * ((t1 - theta) * log(t1 - theta) + (t1 + theta) * log(t1 + theta) - 2.0 * t1 * log(t1) - theta * theta / t1);
62  term2 = r2 * ((t2 - theta) * log(t2 - theta) + (t2 + theta) * log(t2 + theta) - 2.0 * t2 * log(t2) - theta * theta / t2);
63  return g0 - s0 * T_t * theta + T_t * real(term1 + term2);
64 #else
65  return 1e99;
66 #endif
67 }
68 
69 double dg_dp_Ice(double T, double p) {
70 #ifndef __powerpc__
71  std::complex<double> r2_p;
72  double g0_p, theta, pi, pi_0;
73  theta = T / T_t;
74  pi = p / p_t;
75  pi_0 = p_0 / p_t;
76  g0_p = g01 * 1.0 / p_t * pow(pi - pi_0, 1 - 1.0) + g02 * 2.0 / p_t * pow(pi - pi_0, 2 - 1.0) + g03 * 3.0 / p_t * pow(pi - pi_0, 3 - 1.0)
77  + g04 * 4.0 / p_t * pow(pi - pi_0, 4 - 1.0);
78  r2_p = r21 * 1.0 / p_t * pow(pi - pi_0, 1 - 1.0) + r22 * 2.0 / p_t * pow(pi - pi_0, 2 - 1.0);
79  return g0_p + T_t * real(r2_p * ((t2 - theta) * log(t2 - theta) + (t2 + theta) * log(t2 + theta) - 2.0 * t2 * log(t2) - theta * theta / t2));
80 #else
81  return 1e99;
82 #endif
83 }
84 
85 double dg2_dp2_Ice(double T, double p) {
86 #ifndef __powerpc__
87  std::complex<double> r2_pp;
88  double g0_pp, theta, pi, pi_0;
89  theta = T / T_t;
90  pi = p / p_t;
91  pi_0 = p_0 / p_t;
92  g0_pp = g02 * 2.0 * (2.0 - 1.0) / p_t / p_t * pow(pi - pi_0, 2.0 - 2.0) + g03 * 3.0 * (3.0 - 1.0) / p_t / p_t * pow(pi - pi_0, 3.0 - 2.0)
93  + g04 * 4.0 * (4.0 - 1.0) / p_t / p_t * pow(pi - pi_0, 4 - 2.0);
94  r2_pp = r22 * 2.0 / p_t / p_t;
95  return g0_pp + T_t * real(r2_pp * ((t2 - theta) * log(t2 - theta) + (t2 + theta) * log(t2 + theta) - 2.0 * t2 * log(t2) - theta * theta / t2));
96 #else
97  return 1e99;
98 #endif
99 }
100 
101 double dg_dT_Ice(double T, double p) {
102 #ifndef __powerpc__
103  std::complex<double> r2, term1, term2;
104  double theta, pi, pi_0;
105  theta = T / T_t;
106  pi = p / p_t;
107  pi_0 = p_0 / p_t;
108  r2 = r20 * pow(pi - pi_0, 0.0) + r21 * pow(pi - pi_0, 1.0) + r22 * pow(pi - pi_0, 2.0);
109  // The two terms of the summation
110  term1 = r1 * (-log(t1 - theta) + log(t1 + theta) - 2.0 * theta / t1);
111  term2 = r2 * (-log(t2 - theta) + log(t2 + theta) - 2.0 * theta / t2);
112  return -s0 + real(term1 + term2);
113 #else
114  return 1e99;
115 #endif
116 }
117 
118 double h_Ice(double T, double p) {
119 #ifndef __powerpc__
120  // Returned value is in units of J/kg
121  return g_Ice(T, p) - T * dg_dT_Ice(T, p);
122 #else
123  return 1e99;
124 #endif
125 }
126 
127 double rho_Ice(double T, double p) {
128 #ifndef __powerpc__
129  // Returned value is in units of kg/m3
130  return 1 / dg_dp_Ice(T, p);
131 #else
132  return 1e99;
133 #endif
134 }
135 
136 double s_Ice(double T, double p) {
137 #ifndef __powerpc__
138  // Returned value is in units of J/kg/K
139  return -dg_dT_Ice(T, p);
140 #else
141  return 1e99;
142 #endif
143 }