rust_steam/
lib.rs

1const _R: f64 = 0.461526;
2
3const REGION_1_COEFFS: [[f64; 3]; 34] = [
4    [0.0, -2.0, 0.14632971213167],
5    [0.0, -1.0, -0.84548187169114],
6    [0.0, 0.0, -0.37563603672040e1],
7    [0.0, 1.0, 0.33855169168385e1],
8    [0.0, 2.0, -0.95791963387872],
9    [0.0, 3.0, 0.15772038513228],
10    [0.0, 4.0, -0.16616417199501e-1],
11    [0.0, 5.0, 0.81214629983568e-3],
12    [1.0, -9.0, 0.28319080123804e-3],
13    [1.0, -7.0, -0.60706301565874e-3],
14    [1.0, -1.0, -0.18990068218419e-1],
15    [1.0, 0.0, -0.32529748770505e-1],
16    [1.0, 1.0, -0.21841717175414e-1],
17    [1.0, 3.0, -0.52838357969930e-4],
18    [2.0, -3.0, -0.47184321073267e-3],
19    [2.0, 0.0, -0.30001780793026e-3],
20    [2.0, 1.0, 0.47661393906987e-4],
21    [2.0, 3.0, -0.44141845330846e-5],
22    [2.0, 17.0, -0.72694996297594e-15],
23    [3.0, -4.0, -0.31679644845054e-4],
24    [3.0, 0.0, -0.28270797985312e-5],
25    [3.0, 6.0, -0.85205128120103e-9],
26    [4.0, -5.0, -0.22425281908000e-5],
27    [4.0, -2.0, -0.65171222895601e-6],
28    [4.0, 10.0, -0.14341729937924e-12],
29    [5.0, -8.0, -0.40516996860117e-6],
30    [8.0, -11.0, -0.12734301741641e-8],
31    [8.0, -6.0, -0.17424871230634e-9],
32    [21.0, -29.0, -0.68762131295531e-18],
33    [23.0, -31.0, 0.14478307828521e-19],
34    [29.0, -38.0, 0.26335781662795e-22],
35    [30.0, -39.0, -0.11947622640071e-22],
36    [31.0, -40.0, 0.18228094581404e-23],
37    [32.0, -41.0, -0.93537087292458e-25],
38];
39
40const REGION_1_BACK_COEFFS_PH: [[f64; 3]; 20] = [
41            [0.0,  0.0, -0.23872489924521e+3],
42            [0.0,  1.0,  0.40421188637945e+3],
43            [0.0,  2.0,  0.11349746881718e+3],
44            [0.0,  6.0, -0.58457616048039e+1],
45            [0.0, 22.0, -0.15285482413140e-3],
46            [0.0, 32.0, -0.10866707695377e-5],
47            [1.0,  0.0, -0.13391744872602e+2],
48            [1.0,  1.0,  0.43211039183559e+2],
49            [1.0,  2.0, -0.54010067170506e+2],
50            [1.0,  3.0,  0.30535892203916e+2],
51            [1.0,  4.0, -0.65964749423638e+1],
52            [1.0, 10.0,  0.93965400878363e-2],
53            [1.0, 32.0,  0.11573647505340e-6],
54            [2.0, 10.0, -0.25858641282073e-4],
55            [2.0, 32.0, -0.40644363084799e-8],
56            [3.0, 10.0,  0.66456186191635e-7],
57            [3.0, 32.0,  0.80670734103027e-10],
58            [4.0, 32.0, -0.93477771213947e-12],
59            [5.0, 32.0,  0.58265442020601e-14],
60            [6.0, 32.0, -0.15020185953503e-16]
61
62];
63const REGION_1_BACK_COEFFS_PS: [[f64; 3]; 20] = [
64            [0.0,  0.0,  0.17478268058307e+03],
65            [0.0,  1.0,  0.34806930892873e+02],
66            [0.0,  2.0,  0.65292584978455e+01],
67            [0.0,  3.0,  0.33039981775489],
68            [0.0, 11.0, -0.19281382923196e-06],
69            [0.0, 31.0, -0.24909197244573e-22],
70            [1.0,  0.0, -0.26107636489332],
71            [1.0,  1.0,  0.22592965981586],
72            [1.0,  2.0, -0.64256463395226e-01],
73            [1.0,  3.0,  0.78876289270526e-02],
74            [1.0, 12.0,  0.35672110607366e-09],
75            [1.0, 31.0,  0.17332496994895e-23],
76            [2.0,  0.0,  0.56608900654837e-03],
77            [2.0,  1.0, -0.32635483139717e-03],
78            [2.0,  2.0,  0.44778286690632e-04],
79            [2.0,  9.0, -0.51322156908507e-09],
80            [2.0, 31.0, -0.42522657042207e-25],
81            [3.0, 10.0,  0.26400441360689e-12],
82            [3.0, 32.0,  0.78124600459723e-28],
83            [4.0, 32.0, -0.30732199903668e-30]
84];
85
86
87// Region 2
88
89const REGION_2_COEFFS_RES: [[f64; 3]; 43] = [
90    [1.0, 0.0, -0.0017731742473213],
91    [1.0, 1.0, -0.017834862292358],
92    [1.0, 2.0, -0.045996013696365],
93    [1.0, 3.0, -0.057581259083432],
94    [1.0, 6.0, -0.05032527872793],
95    [2.0, 1.0, -0.000033032641670203],
96    [2.0, 2.0, -0.00018948987516315],
97    [2.0, 4.0, -0.0039392777243355],
98    [2.0, 7.0, -0.043797295650573],
99    [2.0, 36.0, -0.000026674547914087],
100    [3.0, 0.0, 2.0481737692309E-08],
101    [3.0, 1.0, 4.3870667284435E-07],
102    [3.0, 3.0, -0.00003227767723857],
103    [3.0, 6.0, -0.0015033924542148],
104    [3.0, 35.0, -0.040668253562649],
105    [4.0, 1.0, -7.8847309559367E-10],
106    [4.0, 2.0, 1.2790717852285E-08],
107    [4.0, 3.0, 4.8225372718507E-07],
108    [5.0, 7.0, 2.2922076337661E-06],
109    [6.0, 3.0, -1.6714766451061E-11],
110    [6.0, 16.0, -0.0021171472321355],
111    [6.0, 35.0, -23.895741934104],
112    [7.0, 0.0, -5.905956432427E-18],
113    [7.0, 11.0, -1.2621808899101E-06],
114    [7.0, 25.0, -0.038946842435739],
115    [8.0, 8.0, 1.1256211360459E-11],
116    [8.0, 36.0, -8.2311340897998],
117    [9.0, 13.0, 1.9809712802088E-08],
118    [10.0, 4.0, 1.0406965210174E-19],
119    [10.0, 10.0, -1.0234747095929E-13],
120    [10.0, 14.0, -1.0018179379511E-09],
121    [16.0, 29.0, -8.0882908646985E-11],
122    [16.0, 50.0, 0.10693031879409],
123    [18.0, 57.0, -0.33662250574171],
124    [20.0, 20.0, 8.9185845355421E-25],
125    [20.0, 35.0, 3.0629316876232E-13],
126    [20.0, 48.0, -4.2002467698208E-06],
127    [21.0, 21.0, -5.9056029685639E-26],
128    [22.0, 53.0, 3.7826947613457E-06],
129    [23.0, 39.0, -1.2768608934681E-15],
130    [24.0, 26.0, 7.3087610595061E-29],
131    [24.0, 40.0, 5.5414715350778E-17],
132    [24.0, 58.0, -9.436970724121E-07],
133];
134
135const REGION_2_COEFFS_IDEAL: [[f64; 2]; 9] = [
136    [0.0, -0.96927686500217e1],
137    [1.0, 0.10086655968018e2],
138    [-5.0, -0.56087911283020e-2],
139    [-4.0, 0.71452738081455e-1],
140    [-3.0, -0.40710498223928],
141    [-2.0, 0.14240819171444e1],
142    [-1.0, -0.43839511319450e1],
143    [2.0, -0.28408632460772],
144    [3.0, 0.21268463753307e-1],
145];
146
147const REGION_4_SATURATION_COEFFS: [f64; 10] = [
148    0.11670521452767e4,
149    -0.72421316703206e6,
150    -0.17073846940092e2,
151    0.12020824702470e5,
152    -0.32325550322333e7,
153    0.14915108613530e2,
154    -0.48232657361591e4,
155    0.40511340542057e6,
156    -0.23855557567849,
157    0.65017534844798e3,
158];
159
160// Region 3
161
162
163const REGION_3_COEFFS: [[f64; 3]; 40] = [
164    [0.0,  0.0,    0.10658070028513e1],
165    [0.0,  0.0,   -0.15732845290239e2],
166    [0.0,  1.0,    0.20944396974307e2],
167    [0.0,  2.0,   -0.76867707878716e1],
168    [0.0,  7.0,    0.26185947787954e1],
169    [0.0,  10.0,  -0.28080781148620e1],
170    [0.0,  12.0,   0.12053369696517e1],
171    [0.0,  23.0,  -0.84566812812502e-2],
172    [1.0,  2.0,   -0.12654315477714e1],
173    [1.0,  6.0,   -0.11524407806681e1],
174    [1.0,  15.0,   0.88521043984318],
175    [1.0,  17.0,  -0.64207765181607],
176    [2.0,  0.0,    0.38493460186671],
177    [2.0,  2.0,   -0.85214708824206],
178    [2.0,  6.0,    0.48972281541877e1],
179    [2.0,  7.0,   -0.30502617256965e1],
180    [2.0,  22.0,   0.39420536879154e-1],
181    [2.0,  26.0,   0.12558408424308],
182    [3.0,  0.0,   -0.27999329698710],
183    [3.0,  2.0,    0.13899799569460e1],
184    [3.0,  4.0,   -0.20189915023570e1],
185    [3.0,  16.0,  -0.82147637173963e-2],
186    [3.0,  26.0,  -0.47596035734923],
187    [4.0,  0.0,    0.43984074473500e-1],
188    [4.0,  2.0,   -0.44476435428739],
189    [4.0,  4.0,    0.90572070719733],
190    [4.0,  26.0,   0.70522450087967],
191    [5.0,  1.0,    0.10770512626332],
192    [5.0,  3.0,   -0.32913623258954],
193    [5.0,  26.0,  -0.50871062041158],
194    [6.0,  0.0,   -0.22175400873096e-1],
195    [6.0,  2.0,    0.94260751665092e-1],
196    [6.0,  26.0,   0.16436278447961],
197    [7.0,  2.0,   -0.13503372241348e-1],
198    [8.0,  26.0,  -0.14834345352472e-1],
199    [9.0,  2.0,    0.57922953628084e-3],
200    [9.0,  26.0,   0.32308904703711e-2],
201    [10.0, 0.0,   0.80964802996215e-4],
202    [10.0, 1.0,  -0.16557679795037e-3],
203    [11.0, 26.0, -0.44923899061815e-4],
204];
205
206pub enum Region {
207    Region1,
208    Region2,
209    Region3,
210    Region4,
211    Region5,
212}
213
214#[derive(Debug)]
215pub enum IAPWSError {
216    OutOfBounds(f64, f64),
217    NotImplemented()
218}
219
220// ===============     Main API ===================
221
222/// Determines which region of the pT chart
223/// a point belongs to.
224///
225/// Returns an error if the point is outside the
226/// bounds of the IAPWS-IF97 correlations.
227///
228/// Temperature is assumed to be in K
229/// Pressure is assumed to be in Pa
230///
231/// Example
232///
233/// ```rust
234/// use rust_steam::{region};
235/// let region = region(300.0, 101325.0).unwrap();
236/// ```
237pub fn region(t: f64, p: f64) -> Result<Region,IAPWSError> {
238    if t < 273.15 || t > 2273.15{
239        return Err(IAPWSError::OutOfBounds(t, p));
240    }
241    if p > 100e6 || p < 0.0{
242        return Err(IAPWSError::OutOfBounds(t, p));
243    }
244
245    if t <= 623.15 {
246        if p > p_sat(t){
247            return Ok(Region::Region1);
248        }else if p < p_sat(t) {
249            return Ok(Region::Region2);
250        }
251        return Ok(Region::Region4);
252    }
253
254    if t <= 863.15 {
255        if p > p_boundary_2_3(t){
256            return Ok(Region::Region3);
257        }else if p < p_boundary_2_3(t) {
258            return Ok(Region::Region2);
259        }
260        return Ok(Region::Region4);
261    }
262
263    if t < 1073.15 {
264        return Ok(Region::Region2);
265    }
266
267    if t >= 1073.15 {
268        if p > 50e6 {
269            return Err(IAPWSError::OutOfBounds(t, p));
270        }
271        return Ok(Region::Region5);
272    }
273    
274
275    Err(IAPWSError::OutOfBounds(t, p))
276}
277
278/// Calculates the water enthalpy in kJ/kg at a given
279/// temperature and pressure.
280///
281/// Temperature is assumed to be in K
282/// Pressure is assumed to be in Pa
283///
284/// Example
285///
286/// ```rust
287/// use rust_steam::{h_tp};
288/// let h = h_tp(300.0, 101325.0).unwrap();
289/// ```
290pub fn h_tp(t: f64, p: f64) -> Result<f64, IAPWSError> {
291    let region = region(t, p)?;
292    return match region {
293        Region::Region1 => Ok(h_tp_1(t, p)),
294        Region::Region2 => Ok(h_tp_2(t, p)),
295        _ => Err(IAPWSError::NotImplemented())
296    }
297}
298
299/// Calculates the water internal energy in kJ/kg at a given
300/// temperature and pressure.
301///
302/// Temperature is assumed to be in K
303/// Pressure is assumed to be in Pa
304///
305/// Example
306///
307/// ```rust
308/// use rust_steam::{u_tp};
309/// let u = u_tp(300.0, 101325.0).unwrap();
310/// ```
311pub fn u_tp(t: f64, p: f64) -> Result<f64, IAPWSError> {
312    let region = region(t, p)?;
313    return match region {
314        Region::Region1 => Ok(u_tp_1(t, p)),
315        Region::Region2 => Ok(u_tp_2(t, p)),
316        _ => Err(IAPWSError::NotImplemented())
317    }
318}
319
320// ================    Region 1 ===================
321
322/// Returns the region-1 gamma
323/// Temperature is assumed to be in K
324/// Pressure is assumed to be in Pa
325fn gamma_1(t: f64, p: f64) -> f64 {
326    let tau = tau_1(t);
327    let pi = pi_1(p);
328    let mut sum = 0.0;
329    for i in 0..REGION_1_COEFFS.len() {
330        let ii = REGION_1_COEFFS[i][0] as i32;
331        let ji = REGION_1_COEFFS[i][1] as i32;
332        let ni = REGION_1_COEFFS[i][2];
333        sum += ni * (7.1 - pi).powi(ii) * (tau - 1.222).powi(ji);
334    }
335    sum
336}
337
338/// Returns the region-1 gamma_pi
339/// Temperature is assumed to be in K
340/// Pressure is assumed to be in Pa
341fn gamma_pi_1(t: f64, p: f64) -> f64 {
342    let tau = tau_1(t);
343    let pi = pi_1(p);
344    let mut sum = 0.0;
345    for i in 0..REGION_1_COEFFS.len() {
346        let ii = REGION_1_COEFFS[i][0] as i32;
347        let ji = REGION_1_COEFFS[i][1] as i32;
348        let ni = REGION_1_COEFFS[i][2];
349        sum += -ni * f64::from(ii) * (7.1 - pi).powi(ii - 1) * (tau - 1.222).powi(ji);
350    }
351    sum
352}
353
354/// Returns the region-1 gamma_pi_pi
355/// Temperature is assumed to be in K
356/// Pressure is assumed to be in Pa
357fn gamma_pi_pi_1(t: f64, p: f64) -> f64 {
358    let tau = tau_1(t);
359    let pi = pi_1(p);
360    let mut sum = 0.0;
361    for i in 0..REGION_1_COEFFS.len() {
362        let ii = REGION_1_COEFFS[i][0] as i32;
363        let ji = REGION_1_COEFFS[i][1] as i32;
364        let ni = REGION_1_COEFFS[i][2];
365        sum += ni * f64::from(ii*(ii-1)) * (7.1 - pi).powi(ii - 2) * (tau - 1.222).powi(ji);
366    }
367    sum
368}
369
370/// Returns the region-1 gamma_tau
371/// Temperature is assumed to be in K
372/// Pressure is assumed to be in Pa
373fn gamma_tau_1(t: f64, p: f64) -> f64 {
374    let tau = tau_1(t);
375    let pi = pi_1(p);
376    let mut sum = 0.0;
377    for i in 0..REGION_1_COEFFS.len() {
378        let ii = REGION_1_COEFFS[i][0] as i32;
379        let ji = REGION_1_COEFFS[i][1] as i32;
380        let ni = REGION_1_COEFFS[i][2];
381        sum += ni * f64::from(ji) * (7.1 - pi).powi(ii) * (tau - 1.222).powi(ji - 1);
382    }
383    sum
384}
385
386/// Returns the region-1 gamma_tau_tau
387/// Temperature is assumed to be in K
388/// Pressure is assumed to be in Pa
389fn gamma_tau_tau_1(t: f64, p: f64) -> f64 {
390    let tau = tau_1(t);
391    let pi = pi_1(p);
392    let mut sum = 0.0;
393    for i in 0..REGION_1_COEFFS.len() {
394        let ii = REGION_1_COEFFS[i][0] as i32;
395        let ji = REGION_1_COEFFS[i][1] as i32;
396        let ni = REGION_1_COEFFS[i][2];
397        sum += ni * f64::from(ji * (ji-1)) * (7.1 - pi).powi(ii) * (tau - 1.222).powi(ji - 2);
398    }
399    sum
400}
401
402/// Returns the region-1 gamma_pi_tau
403/// Temperature is assumed to be in K
404/// Pressure is assumed to be in Pa
405fn gamma_pi_tau_1(t: f64, p: f64) -> f64 {
406    let tau = tau_1(t);
407    let pi = pi_1(p);
408    let mut sum = 0.0;
409    for i in 0..REGION_1_COEFFS.len() {
410        let ii = REGION_1_COEFFS[i][0] as i32;
411        let ji = REGION_1_COEFFS[i][1] as i32;
412        let ni = REGION_1_COEFFS[i][2];
413        sum += -ni * f64::from(ii * ji) * (7.1 - pi).powi(ii - 1) * (tau - 1.222).powi(ji - 1);
414    }
415    sum
416}
417
418/// Returns the region-1 tau
419/// Temperature is assumed to be in K
420/// Pressure is assumed to be in Pa
421fn tau_1(t: f64) -> f64 {
422    1386.0 / t
423}
424
425/// Returns the region-1 pi
426/// Temperature is assumed to be in K
427/// Pressure is assumed to be in Pa
428fn pi_1(p: f64) -> f64 {
429    p / (16.53e6)
430
431}
432
433/// Returns the region-1 specific enthalpy
434/// Temperature is assumed to be in K
435/// Pressure is assumed to be in Pa
436pub fn h_tp_1(t: f64, p: f64) -> f64 {
437    _R * t * tau_1(t) * gamma_tau_1(t, p)
438}
439
440/// Returns the region-1 specific volume
441/// Temperature is assumed to be in K
442/// Pressure is assumed to be in Pa
443pub fn v_tp_1(t: f64, p: f64) -> f64 {
444    // The multiplication by 1000 is necessary to convert R from kJ/kg.K to J/kg.K
445    ((_R * 1000.0) * t / p) * pi_1(p) * gamma_pi_1(t, p)
446}
447
448/// Returns the region-1 specific internal energy
449/// Temperature is assumed to be in K
450/// Pressure is assumed to be in Pa
451pub fn u_tp_1(t: f64, p: f64) -> f64 {
452    _R * t * (tau_1(t) * gamma_tau_1(t, p) - pi_1(p) * gamma_pi_1(t,p))
453}
454
455/// Returns the region-1 specific internal energy
456/// Temperature is assumed to be in K
457/// Pressure is assumed to be in Pa
458pub fn s_tp_1(t: f64, p: f64) -> f64 {
459    _R * (tau_1(t) * gamma_tau_1(t, p) - gamma_1(t,p))
460}
461
462/// Returns the region-1 specific isobaric heat capacity
463/// Temperature is assumed to be in K
464/// Pressure is assumed to be in Pa
465pub fn cp_tp_1(t: f64, p: f64) -> f64 {
466    let tau = tau_1(t);
467    _R * (-tau.powi(2) * gamma_tau_tau_1(t, p))
468}
469
470/// Returns the region-1 specific isochoric heat capacity
471/// Temperature is assumed to be in K
472/// Pressure is assumed to be in Pa
473pub fn cv_tp_1(t: f64, p: f64) -> f64 {
474    let tau = tau_1(t);
475    let corr = (gamma_pi_1(t, p) - tau * gamma_pi_tau_1(t, p)).powi(2)/gamma_pi_pi_1(t, p);
476    _R * (-tau.powi(2) * gamma_tau_tau_1(t, p) + corr)
477}
478
479/// Returns the region-1 speed of sound
480/// Temperature is assumed to be in K
481/// Pressure is assumed to be in Pa
482pub fn w_tp_1(t: f64, p: f64) -> f64 {
483    let tau = tau_1(t);
484    let gamma_pi = gamma_pi_1(t, p);
485    let gamma_pi_tau = gamma_pi_tau_1(t, p);
486    let gamma_pi_pi = gamma_pi_pi_1(t, p);
487    let gamma_tau_tau = gamma_tau_tau_1(t, p);
488    let term = (gamma_pi - tau * gamma_pi_tau).powi(2)/(tau.powi(2) * gamma_tau_tau);
489
490    // The multiplication by 1000 is necessary to convert R from kJ/kg.K to J/kg.K
491    let square = (_R * 1000.0) * t * (gamma_pi.powi(2)/(term - gamma_pi_pi));
492    square.sqrt()
493}
494
495/// Returns the region-1 theta for backwards calculations
496/// Temperature is assumed to be in K
497fn theta_1_back(t: f64) -> f64 {
498    t
499}
500
501/// Returns the region-1 eta for backwards calculations
502/// Enthalpy is assumed to be in kJ/kg
503fn eta_1_back(h: f64) -> f64 {
504    h / 2500.0 
505}
506
507/// Returns the region-1 pi for backwards calculations
508/// Pressure is assumed to be in Pa
509fn pi_1_back(p: f64) -> f64{
510    p / 1e6
511}
512
513/// Returns the region-1 sigma for backwards calculations
514/// Entropy is assumed to be in kJ/kg.K
515fn sigma_1_back(s: f64) -> f64{
516    s
517}
518
519/// Returns the region-1 backward correlation for T(p,s)
520/// Entropy is assumed to be in kJ/kg.K
521/// Pressure is assumed to be in Pa
522pub fn t_ps_1(p: f64, s: f64) -> f64 {
523    let sig = sigma_1_back(s);
524    let pi = pi_1_back(p);
525    let mut sum = 0.0;
526    for i in 0..REGION_1_BACK_COEFFS_PS.len() {
527        let ii = REGION_1_BACK_COEFFS_PS[i][0] as i32;
528        let ji = REGION_1_BACK_COEFFS_PS[i][1] as i32;
529        let ni = REGION_1_BACK_COEFFS_PS[i][2];
530        sum += ni * pi.powi(ii) * (sig + 2.0).powi(ji);
531    }
532    sum
533}
534
535/// Returns the region-1 backward correlation for T(p,h)
536/// Enthalpy is assumed to be in kJ/kg
537/// Pressure is assumed to be in Pa
538pub fn t_ph_1(p: f64, h: f64) -> f64 {
539    let eta = eta_1_back(h);
540    let pi = pi_1_back(p);
541    let mut sum = 0.0;
542    for i in 0..REGION_1_BACK_COEFFS_PH.len() {
543        let ii = REGION_1_BACK_COEFFS_PH[i][0] as i32;
544        let ji = REGION_1_BACK_COEFFS_PH[i][1] as i32;
545        let ni = REGION_1_BACK_COEFFS_PH[i][2];
546        sum += ni * pi.powi(ii) * (eta + 1.0).powi(ji);
547    }
548    sum
549}
550
551// ================    Region 2 ===================
552
553/// Auxiliary equation separating regions 2 and 3
554fn t_boundary_2_3(p: f64) -> f64 {
555    let pi = p / 1e6;
556    let theta_star = 1.0;
557    let n3 = 0.10192970039326e-2;
558    let n4 = 0.57254459862746e3;
559    let n5 = 0.13918839778870e2;
560    let theta = n4 + ((pi - n5)/n3).sqrt();
561    theta * theta_star
562}
563fn p_boundary_2_3(t: f64) -> f64 {
564    let p_star = 1e6;
565    let theta = t / 1.0;
566    let n1 = 0.34805185628969e3;
567    let n2 = -0.11671859879975e1;
568    let n3 = 0.10192970039326e-2;
569    p_star * (n1 + n2 * theta + n3 * theta*theta)
570
571    }
572
573/// Returns the region-2 tau
574/// Temperature is assumed to be in K
575/// Pressure is assumed to be in Pa
576fn tau_2(t: f64) -> f64 {
577    540.0 / t
578}
579
580/// Returns the region-2 pi
581/// Temperature is assumed to be in K
582/// Pressure is assumed to be in Pa
583fn pi_2(p: f64) -> f64 {
584    p / 1e6
585}
586
587/// Returns the region-2 ideal gamma
588/// Temperature is assumed to be in K
589/// Pressure is assumed to be in Pa
590fn gamma_2_ideal(t: f64, p: f64) -> f64 {
591    let pi = pi_2(p);
592    let mut sum: f64 = 0.0;
593    let tau = tau_2(t);
594    for i in 0..REGION_2_COEFFS_IDEAL.len() {
595        let ji = REGION_2_COEFFS_IDEAL[i][0];
596        let ni = REGION_2_COEFFS_IDEAL[i][1];
597        sum += ni * tau.powf(ji);
598    }
599    pi.ln() + sum
600}
601
602/// Returns the region-2 residual gamma
603/// Temperature is assumed to be in K
604/// Pressure is assumed to be in Pa
605fn gamma_2_res(t: f64, p: f64) -> f64 {
606    let pi = pi_2(p);
607    let mut sum: f64 = 0.0;
608    let tau = tau_2(t);
609    for i in 0..REGION_2_COEFFS_RES.len() {
610        let ii = REGION_2_COEFFS_RES[i][0] as i32;
611        let ji = REGION_2_COEFFS_RES[i][1] as i32;
612        let ni = REGION_2_COEFFS_RES[i][2];
613        sum += ni * pi.powi(ii) * (tau - 0.5).powi(ji);
614    }
615    sum
616}
617
618/// Returns the region-2 ideal gamma_tau
619/// Temperature is assumed to be in K
620/// Pressure is assumed to be in Pa
621fn gamma_tau_2_ideal(t: f64, _: f64) -> f64 {
622    let mut sum: f64 = 0.0;
623    let tau = tau_2(t);
624    for i in 0..REGION_2_COEFFS_IDEAL.len() {
625        let ji = REGION_2_COEFFS_IDEAL[i][0] as i32;
626        let ni = REGION_2_COEFFS_IDEAL[i][1];
627        sum += ni * f64::from(ji) * tau.powi(ji - 1);
628    }
629    sum
630}
631
632/// Returns the region-2 ideal gamma_tau_tau
633/// Temperature is assumed to be in K
634/// Pressure is assumed to be in Pa
635fn gamma_tau_tau_2_ideal(t: f64, _: f64) -> f64 {
636    let mut sum: f64 = 0.0;
637    let tau = tau_2(t);
638    for i in 0..REGION_2_COEFFS_IDEAL.len() {
639        let ji = REGION_2_COEFFS_IDEAL[i][0] as i32;
640        let ni = REGION_2_COEFFS_IDEAL[i][1];
641        sum += ni * f64::from(ji * (ji - 1)) * tau.powi(ji - 2);
642    }
643    sum
644}
645
646/// Returns the region-2 ideal gamma_pi
647/// Temperature is assumed to be in K
648/// Pressure is assumed to be in Pa
649fn gamma_pi_2_ideal(_: f64, p: f64) -> f64 {
650    1.0 / pi_2(p)
651}
652
653/// Returns the region-2 residual gamma_tau
654/// Temperature is assumed to be in K
655/// Pressure is assumed to be in Pa
656fn gamma_tau_2_res(t: f64, p: f64) -> f64 {
657    let mut sum: f64 = 0.0;
658    let tau = tau_2(t);
659    let pi = pi_2(p);
660    for i in 0..REGION_2_COEFFS_RES.len() {
661        let ii = REGION_2_COEFFS_RES[i][0] as i32;
662        let ji = REGION_2_COEFFS_RES[i][1] as i32;
663        let ni = REGION_2_COEFFS_RES[i][2];
664        sum += ni * pi.powi(ii) * f64::from(ji) * (tau - 0.5).powi(ji - 1);
665    }
666    sum
667}
668
669/// Returns the region-2 residual gamma_tau_tau
670/// Temperature is assumed to be in K
671/// Pressure is assumed to be in Pa
672fn gamma_tau_tau_2_res(t: f64, p: f64) -> f64 {
673    let mut sum: f64 = 0.0;
674    let tau = tau_2(t);
675    let pi = pi_2(p);
676    for i in 0..REGION_2_COEFFS_RES.len() {
677        let ii = REGION_2_COEFFS_RES[i][0] as i32;
678        let ji = REGION_2_COEFFS_RES[i][1] as i32;
679        let ni = REGION_2_COEFFS_RES[i][2];
680        sum += ni * pi.powi(ii) * f64::from(ji * (ji - 1)) * (tau - 0.5).powi(ji - 2);
681    }
682    sum
683}
684
685/// Returns the region-2 residual gamma_pi
686/// Temperature is assumed to be in K
687/// Pressure is assumed to be in Pa
688fn gamma_pi_2_res(t: f64, p: f64) -> f64 {
689    let mut sum: f64 = 0.0;
690    let tau = tau_2(t);
691    let pi = pi_2(p);
692    for i in 0..REGION_2_COEFFS_RES.len() {
693        let ii = REGION_2_COEFFS_RES[i][0] as i32;
694        let ji = REGION_2_COEFFS_RES[i][1] as i32;
695        let ni = REGION_2_COEFFS_RES[i][2];
696        sum += ni * pi.powi(ii - 1) * f64::from(ii) * (tau - 0.5).powi(ji);
697    }
698    sum
699}
700
701/// Returns the region-2 residual gamma_pi_pi
702/// Temperature is assumed to be in K
703/// Pressure is assumed to be in Pa
704fn gamma_pi_pi_2_res(t: f64, p: f64) -> f64 {
705    let mut sum: f64 = 0.0;
706    let tau = tau_2(t);
707    let pi = pi_2(p);
708    for i in 0..REGION_2_COEFFS_RES.len() {
709        let ii = REGION_2_COEFFS_RES[i][0] as i32;
710        let ji = REGION_2_COEFFS_RES[i][1] as i32;
711        let ni = REGION_2_COEFFS_RES[i][2];
712        sum += ni * pi.powi(ii - 2) * f64::from(ii * (ii - 1)) * (tau - 0.5).powi(ji);
713    }
714    sum
715}
716
717/// Returns the region-2 residual gamma_pi_tau
718/// Temperature is assumed to be in K
719/// Pressure is assumed to be in Pa
720fn gamma_pi_tau_2_res(t: f64, p: f64) -> f64 {
721    let mut sum: f64 = 0.0;
722    let tau = tau_2(t);
723    let pi = pi_2(p);
724    for i in 0..REGION_2_COEFFS_RES.len() {
725        let ii = REGION_2_COEFFS_RES[i][0] as i32;
726        let ji = REGION_2_COEFFS_RES[i][1] as i32;
727        let ni = REGION_2_COEFFS_RES[i][2];
728        sum += ni * pi.powi(ii - 1) * f64::from(ii * ji) * (tau - 0.5).powi(ji - 1);
729    }
730    sum
731}
732
733/// Returns the region-2 specific volume
734/// Temperature is assumed to be in K
735/// Pressure is assumed to be in Pa
736pub fn v_tp_2(t: f64, p: f64) -> f64 {
737    ((_R * 1000.0) * t / p) * pi_2(p) * (gamma_pi_2_ideal(t, p) + gamma_pi_2_res(t, p))
738}
739
740/// Returns the region-2 enthalpy
741/// Temperature is assumed to be in K
742/// Pressure is assumed to be in Pa
743pub fn h_tp_2(t: f64, p: f64) -> f64 {
744    _R * t * tau_2(t) * (gamma_tau_2_ideal(t, p) + gamma_tau_2_res(t, p))
745}
746
747/// Returns the region-2 internal energy
748/// Temperature is assumed to be in K
749/// Pressure is assumed to be in Pa
750pub fn u_tp_2(t: f64, p: f64) -> f64 {
751    let tau = tau_2(t);
752    let pi = pi_2(p);
753    let tau_term = gamma_tau_2_ideal(t, p) + gamma_tau_2_res(t, p);
754    let pi_term = gamma_pi_2_ideal(t, p) + gamma_pi_2_res(t, p);
755    _R * t * (tau * tau_term - pi * pi_term)
756}
757
758/// Returns the region-2 entropy
759/// Temperature is assumed to be in K
760/// Pressure is assumed to be in Pa
761pub fn s_tp_2(t: f64, p: f64) -> f64 {
762    let tau = tau_2(t);
763    let tau_term = gamma_tau_2_ideal(t, p) + gamma_tau_2_res(t, p);
764    let pi_term = gamma_2_ideal(t, p) + gamma_2_res(t, p);
765    _R * (tau * tau_term - pi_term)
766}
767
768/// Returns the region-2 isobaric specific heat
769/// Temperature is assumed to be in K
770/// Pressure is assumed to be in Pa
771pub fn cp_tp_2(t: f64, p: f64) -> f64 {
772    let tau = tau_2(t);
773    - _R  * tau.powi(2) * (gamma_tau_tau_2_ideal(t, p) + gamma_tau_tau_2_res(t, p))
774}
775
776/// Returns the region-2 isochoric specific heat
777/// Temperature is assumed to be in K
778/// Pressure is assumed to be in Pa
779pub fn cv_tp_2(t: f64, p: f64) -> f64 {
780    let tau = tau_2(t);
781    let pi = pi_2(p);
782    let num = (1.0 + pi*gamma_pi_2_res(t, p) - tau * pi * gamma_pi_tau_2_res(t,p)).powi(2);
783    let den = 1.0 - pi.powi(2) * gamma_pi_pi_2_res(t,p);
784    cp_tp_2(t, p) - _R * num/ den
785}
786
787/// Returns the region-2 sound velocity
788/// Temperature is assumed to be in K
789/// Pressure is assumed to be in Pa
790pub fn w_tp_2(t: f64, p: f64) -> f64 {
791    let tau = tau_2(t);
792    let pi = pi_2(p);
793    let num = 1.0 + 2.0*pi*gamma_pi_2_res(t, p) +  pi.powi(2) * gamma_pi_2_res(t,p).powi(2);
794    let subnum = (1.0 + pi * gamma_pi_2_res(t, p) - tau * pi * gamma_pi_tau_2_res(t, p)).powi(2);
795    let subden = tau.powi(2)*(gamma_tau_tau_2_ideal(t, p) + gamma_tau_tau_2_res(t, p));
796    let den = 1.0 - pi.powi(2) * gamma_pi_pi_2_res(t,p) + subnum / subden;
797    ((_R * 1000.0 * t) * num / den).sqrt()
798}
799
800// ================    Region 4 ===================
801
802/// Returns the saturation pressure in Pa
803/// Temperature is assumed to be in K
804pub fn p_sat(t: f64) -> f64 {
805    let n1 = REGION_4_SATURATION_COEFFS[0];
806    let n2 = REGION_4_SATURATION_COEFFS[1];
807    let n3 = REGION_4_SATURATION_COEFFS[2];
808    let n4 = REGION_4_SATURATION_COEFFS[3];
809    let n5 = REGION_4_SATURATION_COEFFS[4];
810    let n6 = REGION_4_SATURATION_COEFFS[5];
811    let n7 = REGION_4_SATURATION_COEFFS[6];
812    let n8 = REGION_4_SATURATION_COEFFS[7];
813    let n9 = REGION_4_SATURATION_COEFFS[8];
814    let n10 = REGION_4_SATURATION_COEFFS[9];
815
816    let theta = t + n9 / (t - n10);
817
818    let coef_a = theta * theta + n1 * theta + n2;
819    let coef_b = n3 * theta * theta + n4 * theta + n5;
820    let coef_c = n6 * theta * theta + n7 * theta + n8;
821    (2.0 * coef_c / (-coef_b + (coef_b * coef_b - 4.0 * coef_a * coef_c).sqrt())).powi(4) * 1e6
822}
823
824/// Returns the saturation temperature in K
825/// Pressure is assumed to be in Pa
826pub fn t_sat(p: f64) -> f64 {
827    let n1 = REGION_4_SATURATION_COEFFS[0];
828    let n2 = REGION_4_SATURATION_COEFFS[1];
829    let n3 = REGION_4_SATURATION_COEFFS[2];
830    let n4 = REGION_4_SATURATION_COEFFS[3];
831    let n5 = REGION_4_SATURATION_COEFFS[4];
832    let n6 = REGION_4_SATURATION_COEFFS[5];
833    let n7 = REGION_4_SATURATION_COEFFS[6];
834    let n8 = REGION_4_SATURATION_COEFFS[7];
835    let n9 = REGION_4_SATURATION_COEFFS[8];
836    let n10 = REGION_4_SATURATION_COEFFS[9];
837
838    let beta = (p / 1e6).powf(0.25);
839
840    let coef_e = beta * beta + n3 * beta + n6;
841    let coef_f = n1 * beta * beta + n4 * beta + n7;
842    let coef_g = n2 * beta * beta + n5 * beta + n8;
843
844    let coef_d = 2.0 * coef_g / (-coef_f - (coef_f * coef_f - 4.0 * coef_e * coef_g).sqrt());
845
846    (n10 + coef_d - ((n10 + coef_d).powi(2) - 4.0 * (n9 + n10 * coef_d)).sqrt()) / 2.0
847}
848
849#[cfg(test)]
850mod tests {
851    use super::*;
852
853    fn is_close(x: f64, y: f64, tol: f64) -> bool {
854        (x - y).abs() < tol
855    }
856
857    #[test]
858    fn canary() {}
859
860    #[test]
861    fn region_1_enthalpy() {
862        let h1 = h_tp_1(300.0, 3e6);
863
864        assert!(is_close(h1 / 1000.0, 0.115331273, 1e-9));
865
866        let h1 = h_tp_1(300.0, 80e6);
867        assert!(is_close(h1 / 1000.0, 0.184142828, 1e-9));
868
869        let h1 = h_tp_1(500.0, 3e6);
870        assert!(is_close(h1 / 1000.0, 0.975542239, 1e-9));
871    }
872
873    #[test]
874    fn region_1_internal_energy() {
875        let u1 = u_tp_1(300.0, 3e6);
876
877        assert!(is_close(u1 / 1000.0, 0.112324818, 1e-9));
878
879        let u1 = u_tp_1(300.0, 80e6);
880        assert!(is_close(u1 / 1000.0, 0.106448356, 1e-9));
881
882        let u1 = u_tp_1(500.0, 3e6);
883        assert!(is_close(u1 / 1000.0, 0.971934985, 1e-9));
884    }
885
886    #[test]
887    fn region_1_entropy() {
888        let s1 = s_tp_1(300.0, 3e6);
889        assert!(is_close(s1 , 0.392294792, 1e-9));
890
891        let s1 = s_tp_1(300.0, 80e6);
892        assert!(is_close(s1 , 0.368563852, 1e-9));
893
894        let s1 = s_tp_1(500.0, 3e6);
895        assert!(is_close(s1 , 0.258041912e1, 1e-9));
896    }
897
898    #[test]
899    fn region_1_cp() {
900        let cp1 = cp_tp_1(300.0, 3e6);
901        println!("{}", cp1);
902        assert!(is_close(cp1/10.0 , 0.417301218, 1e-9));
903
904        let cp1 = cp_tp_1(300.0, 80e6);
905        assert!(is_close(cp1/10.0 , 0.401008987, 1e-9));
906
907        let cp1 = cp_tp_1(500.0, 3e6);
908        assert!(is_close(cp1/10.0 , 0.465580682, 1e-9));
909    }
910
911    #[test]
912    fn region_1_sound_velocity() {
913        let w1 = w_tp_1(300.0, 3e6);
914        assert!(is_close(w1/10000.0 , 0.150773921 , 1e-9));
915
916        let w1 = w_tp_1(300.0, 80e6);
917        assert!(is_close(w1/10000.0 , 0.163469054, 1e-9));
918
919        let w1 = w_tp_1(500.0, 3e6);
920        assert!(is_close(w1/10000.0 , 0.124071337, 1e-9));
921    }
922
923    #[test]
924    fn region_1_specific_volume() {
925        let v1 = v_tp_1(300.0, 3e6);
926        assert!(is_close(v1 * 100.0, 0.100215168, 1e-9));
927
928        let v1 = v_tp_1(300.0, 80e6);
929        assert!(is_close(v1 * 1000.0 , 0.971180894, 1e-9));
930
931        let v1 = v_tp_1(500.0, 3e6);
932        assert!(is_close(v1 * 100.0, 0.120241800, 1e-9));
933    }
934
935    #[test]
936    fn region_1_backwards_t_ph() {
937        let t = t_ph_1( 3e6, 500.0);
938        assert!(is_close(t / 1000.0, 0.391798509, 1e-9));
939
940        let t = t_ph_1( 80e6, 500.0);
941        assert!(is_close(t / 1000.0, 0.378108626, 1e-9));
942
943        let t = t_ph_1( 80e6, 1500.0);
944        assert!(is_close(t / 1000.0, 0.611041229, 1e-9));
945
946    }
947
948    #[test]
949    fn region_1_backwards_t_ps() {
950        let t = t_ps_1( 3e6, 0.5);
951        assert!(is_close(t / 1000.0, 0.307842258, 1e-9));
952
953        let t = t_ps_1( 80e6, 0.5);
954        assert!(is_close(t / 1000.0, 0.309979785, 1e-9));
955
956        let t = t_ps_1( 80e6, 3.0);
957        assert!(is_close(t / 1000.0, 0.565899909, 1e-9));
958
959    }
960
961    #[test]
962    fn region_2_specific_volume() {
963        let v = v_tp_2(300.0, 0.0035e6);
964        assert!(is_close(v / 100.0, 0.394913866, 1e-9));
965
966        let v = v_tp_2(700.0, 0.0035e6);
967        assert!(is_close(v / 100.0, 0.923015898 , 1e-9));
968
969        let v = v_tp_2(700.0, 30e6);
970        assert!(is_close(v / 0.01, 0.542946619, 1e-9));
971    }
972
973    #[test]
974    fn region_2_enthalpy() {
975        let h = h_tp_2(300.0, 0.0035e6);
976        assert!(is_close(h / 10000.0, 0.254991145, 1e-9));
977
978        let h = h_tp_2(700.0, 0.0035e6);
979        assert!(is_close(h / 10000.0, 0.333568375, 1e-9));
980
981        let h = h_tp_2(700.0, 30e6);
982        assert!(is_close(h / 10000.0, 0.263149474, 1e-9));
983    }
984
985    #[test]
986    fn region_2_internal_energy() {
987        let u = u_tp_2(300.0, 0.0035e6);
988        assert!(is_close(u / 10000.0, 0.241169160, 1e-9));
989
990        let u = u_tp_2(700.0, 0.0035e6);
991        assert!(is_close(u / 10000.0, 0.301262819, 1e-9));
992
993        let u = u_tp_2(700.0, 30e6);
994        assert!(is_close(u / 10000.0, 0.246861076, 1e-9));
995    }
996
997    #[test]
998    fn region_2_entropy() {
999        let s = s_tp_2(300.0, 0.0035e6);
1000        println!("{}", s);
1001        assert!(is_close(s / 10.0, 0.852238967, 1e-9));
1002
1003        let s = s_tp_2(700.0, 0.0035e6);
1004        assert!(is_close(s / 100.0, 0.101749996, 1e-9));
1005
1006        let s = s_tp_2(700.0, 30e6);
1007        assert!(is_close(s / 10.0, 0.517540298, 1e-9));
1008    }
1009
1010    #[test]
1011    fn region_2_cp() {
1012        let cp = cp_tp_2(300.0, 0.0035e6);
1013        assert!(is_close(cp / 10.0, 0.191300162, 1e-9));
1014
1015        let cp = cp_tp_2(700.0, 0.0035e6);
1016        assert!(is_close(cp / 10.0, 0.208141274 , 1e-9));
1017
1018        let cp = cp_tp_2(700.0, 30e6);
1019        assert!(is_close(cp / 100.0, 0.103505092, 1e-9));
1020    }
1021
1022    #[test]
1023    fn region_2_sound_velocity() {
1024        let w = w_tp_2(300.0, 0.0035e6);
1025        assert!(is_close(w / 1000.0, 0.427920172, 1e-9));
1026
1027        let w = w_tp_2(700.0, 0.0035e6);
1028        assert!(is_close(w / 1000.0, 0.644289068, 1e-9));
1029
1030        let w = w_tp_2(700.0, 30e6);
1031        assert!(is_close(w / 1000.0, 0.480386523, 1e-9));
1032    }
1033
1034    #[test]
1035    fn region_2_3_auxiliary_boundary() {
1036        let p = p_boundary_2_3(623.15);
1037        assert!(is_close(p / 1e8, 0.165291643, 1e-9));
1038
1039        let t = t_boundary_2_3(0.165291643e8);
1040        assert!(is_close(t / 1e3, 0.623150000, 1e-9));
1041
1042    }
1043
1044
1045    // Region 4
1046
1047    #[test]
1048    fn region_4_saturation_pressure() {
1049        let ps = p_sat(300.0);
1050        assert!(is_close(ps / 10000.0, 0.353658941, 1e-9));
1051
1052        let ps = p_sat(500.0);
1053        assert!(is_close(ps / 1e7, 0.263889776, 1e-9));
1054
1055        let ps = p_sat(600.0);
1056        assert!(is_close(ps / 1e8, 0.123443146, 1e-9));
1057    }
1058
1059    #[test]
1060    fn region_4_saturation_temperature() {
1061        let ts = t_sat(0.1e6);
1062        assert!(is_close(ts / 1000.0, 0.372755919, 1e-9));
1063
1064        let ts = t_sat(1e6);
1065        assert!(is_close(ts / 1000.0, 0.453035632, 1e-9));
1066
1067        let ts = t_sat(10e6);
1068        assert!(is_close(ts / 1000.0, 0.584149488, 1e-9));
1069    }
1070}
rust_steam/lib.rs

rust_steam/
lib.rs
