1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
//! # Gas Struct
//! The `Gas` struct is built from a .json file, use `air.json` as reference. 
//! In this file, there are two main structs: `"phase"` and `"species_data"`.
//! `"speciesArray"` determines which species can exist in the gas. 
//! `"state"` determines the initial state of the gas upon instantiation.
//! `"species_data"` is a vector containing information about all species in `"speciesArray"`.
//! Therefore, every specie declared in `"speciesArray"` must have your basic information declared in
//! `"species_data"`. The order in which the data is added is not importante. There can be more
//! `"specie_data"` than `"speciesArray"` but never the opposite.

#![allow(non_snake_case)]
use crate::base::constants::{R, _P_REF};
use crate::reaction::{
    json_data,
    thermo::{ThermoInterp, ThermoProp},
};
use ndarray::prelude::*;
use std::collections::HashMap;

/// Struct of an ideal gas
#[derive(Debug, Clone)]
pub struct Gas {
    name: String,
    species: Vec<String>,
    species_molar_weight: Array1<f64>,
    mole_frac: Array1<f64>,
    species_atoms: Vec<HashMap<String, f64>>,
    thermo_interp: Vec<ThermoInterp>,
    thermo_prop: ThermoProp,
    num_species: usize,
}

impl Gas {
    /// Creates a `Gas` object from a file .json file
    pub fn new(gas_file: &str) -> Gas {
        let json_output = json_data::read_and_treat_json(gas_file);
        let num_species = json_output.species.len();

        let mut gas = Gas {
            name: json_output.name,
            species: json_output.species,
            species_molar_weight: json_output.species_molar_weight,
            mole_frac: json_output.mol_frac,
            species_atoms: json_output.species_atoms,
            thermo_interp: json_output.thermo_interp,
            thermo_prop: ThermoProp::new(),
            num_species,
        };
        gas.TP(json_output.ini_temp, json_output.ini_press);
        gas
    }

    /// Set temperature and pressure. Thermo properties are recalculated  
    /// # Examples
    /// ```
    /// let temp = 350;
    /// let press = 2e5;
    /// gas.TP(temp, press);
    /// assert_eq!(350, gas.T());
    /// assert_eq!(2e5, gas.P());
    /// ```
    pub fn TP<'a>(&'a mut self, temp: f64, press: f64) -> &'a mut Self {
        self.thermo_prop.T = temp;
        self.thermo_prop.P = press;
        self.update_prop();
        self
    }

    /// Set temperature, pressure, and mole fraction of species. Thermo properties are recalculated  
    /// # Examples
    /// ```
    /// let temp = 350;
    /// let press = 2e5;
    /// let mole_frac = "O2:0.21, N2:0.79";
    /// gas.TPX(temp, press, mol_frac);
    /// assert_eq!(350, gas.T());
    /// assert_eq!(2e5, gas.P());
    /// assert_eq!(0.21, gas.mol_frac_of("O2"));
    /// ```
    pub fn TPX<'a>(&'a mut self, temp: f64, press: f64, mol_frac: &str) -> &'a mut Self {
        self.thermo_prop.T = temp;
        self.thermo_prop.P = press;
        self.X(mol_frac);
        self
    }

    /// Set mole fraction of species from `&str`. Thermo properties are recalculated  
    /// # Examples
    /// ```
    /// let mole_frac = "O2:0.21, N2:0.79";
    /// gas.X(mol_frac);
    /// assert_eq!(0.21, gas.mol_frac_of("O2"));
    /// assert_eq!(0.79, gas.mol_frac_of("N2"));
    /// ```
    pub fn X<'a>(&'a mut self, mole_frac: &str) -> &'a mut Self {
        let X = self.break_str_into_X_array(mole_frac);
        self.mole_frac = X;
        self.update_prop();
        self
    }

    /// Set mole fraction of species from ndarray::Array1<f64>, must be the same size.
    /// Thermo properties are recalculated
    pub fn X_array<'a>(&'a mut self, mol_frac: &Array1<f64>) -> &'a mut Self {
        if (mol_frac.sum() - 1.0).abs() > 1e-8 {
            println!(
                "Error!\n mol_fraction must sum 1.0: mol_frac = {}",
                mol_frac.sum()
            );
            std::process::exit(1);
        }
        self.mole_frac.assign(&mol_frac);
        self.update_prop();
        self
    }

    pub fn TPX_array<'a>(
        &'a mut self,
        temp: f64,
        press: f64,
        mol_frac: &Array1<f64>,
    ) -> &'a mut Self {
        self.thermo_prop.T = temp;
        self.thermo_prop.P = press;
        self.X_array(mol_frac);
        self
    }

    /// Returns the equivalent mole fraction, if a `mass`, in kg, of `self`
    /// were mixed with other gases in `add_gas<(mass, composition)>`
    /// # Examples
    /// ```
    /// gas.X("O2:1.0");
    /// let additional_gas = vec![(0.329337487, "N2:1.0")];
    /// let mix_mole_frac: Array1<f64> = gas.if_mixed_with(0.100, additional_gas);
    /// //mix_mole_frac has: O2:0.21, N2:0.79
    /// ```
    pub fn if_mixed_with(&self, mass: f64, add_gas: Vec<(f64, &str)>) -> Array1<f64> {
        if add_gas.iter().all(|(m, _)| *m == 0.0 as f64) {
            return self.mole_frac().clone();
        }
        let mut added_moles = Array::from_elem(self.num_species, 0.0);
        for (m, composition) in add_gas.iter() {
            let mole_frac = self.break_str_into_X_array(composition);
            let molar_weight = mole_frac.dot(&self.species_molar_weight);
            added_moles = added_moles + mole_frac * (*m) / molar_weight;
        }

        let new_mole_frac: Array1<f64>;
        let current_moles = self.mole_frac() * mass / self.M();
        let total_moles = current_moles.sum() + added_moles.sum();
        new_mole_frac = &(&current_moles + &added_moles) / total_moles;
        new_mole_frac
    }

    fn update_prop(&mut self) {
        let mut cp_array = Array::from_elem(self.species_molar_weight.len(), 0.);
        let mut h_array = Array::from_elem(self.species_molar_weight.len(), 0.);
        let mut s_array = Array::from_elem(self.species_molar_weight.len(), 0.);
        let mut cp_h_s: (f64, f64, f64);

        for (i, thermo_interp) in self.thermo_interp.iter().enumerate() {
            if self.T() < thermo_interp.Tmid() {
                cp_h_s = ThermoInterp::calc_thermo_properties(thermo_interp.coeffs_low(), self.T());
            } else {
                cp_h_s =
                    ThermoInterp::calc_thermo_properties(thermo_interp.coeffs_high(), self.T());
            }
            cp_array[i] = cp_h_s.0;
            h_array[i] = cp_h_s.1;
            s_array[i] = cp_h_s.2;
        }

        cp_array = R * cp_array; // [J/kmol/K]
        h_array = R * self.T() * h_array; // [J/kmol]
                                          // tmp: Array1<f64> = ln(mole_frac*P/P_ref)
        let tmp = self.mole_frac().mapv(|x| -> f64 {
            if x == 0.0 {
                0.0
            } else {
                (x * self.P() / _P_REF).ln()
            }
        });
        s_array = R * (s_array - tmp); // [J/kmol/K]
        let cv_array = &cp_array - R;
        // All properties in mass basis
        self.thermo_prop.M = self.mole_frac.dot(&self.species_molar_weight);
        self.thermo_prop.cp = self.mole_frac.dot(&(cp_array / self.M()));
        self.thermo_prop.cv = self.mole_frac.dot(&(cv_array / self.M()));
        self.thermo_prop.h = self.mole_frac.dot(&(h_array / self.M()));
        self.thermo_prop.s = self.mole_frac.dot(&(s_array / self.M()));
        self.thermo_prop.R = R / self.M();
        self.thermo_prop.k = self.cp() / self.cv();
        self.thermo_prop.rho = self.P() / (self.R() * self.T());
        self.thermo_prop.e = self.h() - self.P() / self.rho();
        self.thermo_prop.a = (self.k() * self.R() * self.T()).sqrt();
        self.thermo_prop.mu =
            (1.458e-6) * (self.T() * self.T() * self.T() / (self.T() + 110.4)).sqrt();
    }
    /// gas name
    pub fn name(&self) -> String {
        self.name.clone()
    }
    /// Vec with all species
    pub fn species(&self) -> &Vec<String> {
        &self.species
    }
    /// array of species mole fraction
    pub fn mole_frac(&self) -> &Array1<f64> {
        &self.mole_frac
    }
    /// array of all species molar weight [kg/kmol]
    pub fn mole_weight(&self) -> &Array1<f64> {
        &self.species_molar_weight
    }
    /// temperature [K]
    pub fn T(&self) -> f64 {
        self.thermo_prop.T
    }
    /// pressure [Pa]
    pub fn P(&self) -> f64 {
        self.thermo_prop.P
    }
    /// density [kg/m³]
    pub fn rho(&self) -> f64 {
        self.thermo_prop.rho
    }
    /// specific heat at constant pressure [J/kg.K]
    pub fn cp(&self) -> f64 {
        self.thermo_prop.cp
    }
    /// specific heat at constant volume [J/kg.K]
    pub fn cv(&self) -> f64 {
        self.thermo_prop.cv
    }
    /// ideal gas constant [J/kg.K]
    pub fn R(&self) -> f64 {
        self.thermo_prop.R
    }
    /// specific heat ratio (cp/cv)
    pub fn k(&self) -> f64 {
        self.thermo_prop.k
    }
    /// molar mass [kg/kmol]
    pub fn M(&self) -> f64 {
        self.thermo_prop.M
    }
    /// internal energy [J]
    pub fn e(&self) -> f64 {
        self.thermo_prop.e
    }
    /// enthalpy [J/kg]
    pub fn h(&self) -> f64 {
        self.thermo_prop.h
    }
    /// entropy [J/kg.K]
    pub fn s(&self) -> f64 {
        self.thermo_prop.s
    }
    /// sound speed [m/s]
    pub fn a(&self) -> f64 {
        self.thermo_prop.a
    }
    /// viscosity [m²/s]
    pub fn mu(&self) -> f64 {
        self.thermo_prop.mu
    }
    /// return the number os species
    pub fn num_species(&self) -> usize {
        self.num_species
    }
    /// check if `specie` contains in `species` vector
    pub fn contains_specie(&self, specie: &str) -> bool {
        self.species.contains(&specie.to_string())
    }
    /// get the index of `specie`
    pub fn get_specie_index(&self, specie: &str) -> usize {
        match self.species.iter().position(|s| s == specie) {
            Some(i) => i,
            None => {
                println!(
                    "Error at `get_specie_index()`. Specie \"{}\" not found",
                    specie
                );
                std::process::exit(1);
            }
        }
    }
    /// return the mole fraction of `specie`
    pub fn mole_frac_of(&self, specie: &str) -> f64 {
        let i = match self.species.iter().position(|s| s == specie) {
            Some(i) => i,
            None => {
                println!("Error at `mol_frac_of`. Specie \"{}\" not found", specie);
                std::process::exit(1);
            }
        };
        self.mole_frac[i]
    }
    /// return the mole weight of `specie`
    pub fn mole_weight_of(&self, specie: &str) -> f64 {
        let i = match self.species.iter().position(|s| s == specie) {
            Some(i) => i,
            None => {
                println!("Error at `mol_frac_of`. Specie {} not found", specie);
                std::process::exit(1);
            }
        };
        self.species_molar_weight[i]
    }
    /// return a `&HashMap` containing the `atoms_name` and `num_of_atoms` of `specie`
    pub fn atoms_of(&self, specie: &str) -> &HashMap<String, f64> {
        let i = match self.species.iter().position(|s| s == specie) {
            Some(i) => i,
            None => {
                println!("Error at `mol_frac_of`. Specie {} not found", specie);
                std::process::exit(1);
            }
        };
        &self.species_atoms[i]
    }

    fn break_str_into_X_array(&self, mole_frac: &str) -> Array1<f64> {
        let strings: Vec<String> = mole_frac
            .replace(&[',', '\"'][..], "")
            .split_whitespace()
            .map(|s| s.to_string())
            .collect();
        let mut X = Array::from_elem(self.species.len(), 0.);

        for word in strings.iter() {
            let specie: Vec<&str> = word.split(":").collect(); // specie should be like ["O2", "0.21"]
            if specie.len() != 2 {
                println!("Error setting mole fraction! \"{}\" is not valid", word);
                println!("Valid example: \"O2:0.21, N2:0.79\"");
                std::process::exit(1);
            }
            if self.species.contains(&specie[0].to_string()) {
                let i = self.species.iter().position(|s| s == specie[0]).unwrap();
                X[i] = specie[1].parse().unwrap();
            } else {
                println!(
                    "Error!\n Specie `{}` was not found in `species_data` in file `{}`",
                    specie[0],
                    self.name()
                );
                std::process::exit(1);
            }
        }
        if ((X.sum() - 1.0) as f64).abs() > 1e-8 {
            println!("Error!\n mol_fraction must sum 1.0: mol_frac = {}", X.sum());
            std::process::exit(1);
        }
        X
    }
}