feos_core/state/
properties.rs

1use super::{Contributions, Derivative::*, PartialDerivative, State};
2use crate::equation_of_state::{IdealGas, Molarweight, Residual};
3use crate::ReferenceSystem;
4use ndarray::Array1;
5use quantity::*;
6use std::ops::Div;
7use typenum::P2;
8
9impl<E: Residual + IdealGas> State<E> {
10    fn get_or_compute_derivative(
11        &self,
12        derivative: PartialDerivative,
13        contributions: Contributions,
14    ) -> f64 {
15        let residual = match contributions {
16            Contributions::IdealGas => None,
17            _ => Some(self.get_or_compute_derivative_residual(derivative)),
18        };
19
20        let ideal_gas = match contributions {
21            Contributions::Residual => None,
22            _ => Some(match derivative {
23                PartialDerivative::Zeroth => {
24                    let new_state = self.derive0();
25                    self.eos.ideal_gas_helmholtz_energy(&new_state) * new_state.temperature
26                }
27                PartialDerivative::First(v) => {
28                    let new_state = self.derive1(v);
29                    (self.eos.ideal_gas_helmholtz_energy(&new_state) * new_state.temperature).eps
30                }
31                PartialDerivative::Second(v) => {
32                    let new_state = self.derive2(v);
33                    (self.eos.ideal_gas_helmholtz_energy(&new_state) * new_state.temperature).v2
34                }
35                PartialDerivative::SecondMixed(v1, v2) => {
36                    let new_state = self.derive2_mixed(v1, v2);
37                    (self.eos.ideal_gas_helmholtz_energy(&new_state) * new_state.temperature)
38                        .eps1eps2
39                }
40                PartialDerivative::Third(v) => {
41                    let new_state = self.derive3(v);
42                    (self.eos.ideal_gas_helmholtz_energy(&new_state) * new_state.temperature).v3
43                }
44            }),
45        };
46
47        match (ideal_gas, residual) {
48            (Some(i), Some(r)) => i + r,
49            (Some(i), None) => i,
50            (None, Some(r)) => r,
51            (None, None) => unreachable!(),
52        }
53    }
54
55    /// Chemical potential: $\mu_i=\left(\frac{\partial A}{\partial N_i}\right)_{T,V,N_j}$
56    pub fn chemical_potential(&self, contributions: Contributions) -> MolarEnergy<Array1<f64>> {
57        Quantity::from_reduced(Array1::from_shape_fn(self.eos.components(), |i| {
58            self.get_or_compute_derivative(PartialDerivative::First(DN(i)), contributions)
59        }))
60    }
61
62    /// Partial derivative of chemical potential w.r.t. temperature: $\left(\frac{\partial\mu_i}{\partial T}\right)_{V,N_i}$
63    pub fn dmu_dt(
64        &self,
65        contributions: Contributions,
66    ) -> <MolarEnergy<Array1<f64>> as Div<Temperature>>::Output {
67        Quantity::from_reduced(Array1::from_shape_fn(self.eos.components(), |i| {
68            self.get_or_compute_derivative(PartialDerivative::SecondMixed(DT, DN(i)), contributions)
69        }))
70    }
71
72    /// Molar isochoric heat capacity: $c_v=\left(\frac{\partial u}{\partial T}\right)_{V,N_i}$
73    pub fn molar_isochoric_heat_capacity(&self, contributions: Contributions) -> MolarEntropy {
74        self.temperature * self.ds_dt(contributions) / self.total_moles
75    }
76
77    /// Partial derivative of the molar isochoric heat capacity w.r.t. temperature: $\left(\frac{\partial c_V}{\partial T}\right)_{V,N_i}$
78    pub fn dc_v_dt(
79        &self,
80        contributions: Contributions,
81    ) -> <MolarEntropy as Div<Temperature>>::Output {
82        (self.temperature * self.d2s_dt2(contributions) + self.ds_dt(contributions))
83            / self.total_moles
84    }
85
86    /// Molar isobaric heat capacity: $c_p=\left(\frac{\partial h}{\partial T}\right)_{p,N_i}$
87    pub fn molar_isobaric_heat_capacity(&self, contributions: Contributions) -> MolarEntropy {
88        match contributions {
89            Contributions::Residual => self.residual_molar_isobaric_heat_capacity(),
90            _ => {
91                self.temperature / self.total_moles
92                    * (self.ds_dt(contributions)
93                        - self.dp_dt(contributions).powi::<P2>() / self.dp_dv(contributions))
94            }
95        }
96    }
97
98    /// Entropy: $S=-\left(\frac{\partial A}{\partial T}\right)_{V,N_i}$
99    pub fn entropy(&self, contributions: Contributions) -> Entropy {
100        Entropy::from_reduced(
101            -self.get_or_compute_derivative(PartialDerivative::First(DT), contributions),
102        )
103    }
104
105    /// Molar entropy: $s=\frac{S}{N}$
106    pub fn molar_entropy(&self, contributions: Contributions) -> MolarEntropy {
107        self.entropy(contributions) / self.total_moles
108    }
109
110    /// Partial molar entropy: $s_i=\left(\frac{\partial S}{\partial N_i}\right)_{T,p,N_j}$
111    pub fn partial_molar_entropy(&self) -> MolarEntropy<Array1<f64>> {
112        let c = Contributions::Total;
113        -(self.dmu_dt(c) + self.dp_dni(c) * (self.dp_dt(c) / self.dp_dv(c)))
114    }
115
116    /// Partial derivative of the entropy w.r.t. temperature: $\left(\frac{\partial S}{\partial T}\right)_{V,N_i}$
117    pub fn ds_dt(&self, contributions: Contributions) -> <Entropy as Div<Temperature>>::Output {
118        Quantity::from_reduced(
119            -self.get_or_compute_derivative(PartialDerivative::Second(DT), contributions),
120        )
121    }
122
123    /// Second partial derivative of the entropy w.r.t. temperature: $\left(\frac{\partial^2 S}{\partial T^2}\right)_{V,N_i}$
124    pub fn d2s_dt2(
125        &self,
126        contributions: Contributions,
127    ) -> <<Entropy as Div<Temperature>>::Output as Div<Temperature>>::Output {
128        Quantity::from_reduced(
129            -self.get_or_compute_derivative(PartialDerivative::Third(DT), contributions),
130        )
131    }
132
133    /// Enthalpy: $H=A+TS+pV$
134    pub fn enthalpy(&self, contributions: Contributions) -> Energy {
135        self.temperature * self.entropy(contributions)
136            + self.helmholtz_energy(contributions)
137            + self.pressure(contributions) * self.volume
138    }
139
140    /// Molar enthalpy: $h=\frac{H}{N}$
141    pub fn molar_enthalpy(&self, contributions: Contributions) -> MolarEnergy {
142        self.enthalpy(contributions) / self.total_moles
143    }
144
145    /// Partial molar enthalpy: $h_i=\left(\frac{\partial H}{\partial N_i}\right)_{T,p,N_j}$
146    pub fn partial_molar_enthalpy(&self) -> MolarEnergy<Array1<f64>> {
147        let s = self.partial_molar_entropy();
148        let mu = self.chemical_potential(Contributions::Total);
149        s * self.temperature + mu
150    }
151
152    /// Helmholtz energy: $A$
153    pub fn helmholtz_energy(&self, contributions: Contributions) -> Energy {
154        Energy::from_reduced(
155            self.get_or_compute_derivative(PartialDerivative::Zeroth, contributions),
156        )
157    }
158
159    /// Molar Helmholtz energy: $a=\frac{A}{N}$
160    pub fn molar_helmholtz_energy(&self, contributions: Contributions) -> MolarEnergy {
161        self.helmholtz_energy(contributions) / self.total_moles
162    }
163
164    /// Internal energy: $U=A+TS$
165    pub fn internal_energy(&self, contributions: Contributions) -> Energy {
166        self.temperature * self.entropy(contributions) + self.helmholtz_energy(contributions)
167    }
168
169    /// Molar internal energy: $u=\frac{U}{N}$
170    pub fn molar_internal_energy(&self, contributions: Contributions) -> MolarEnergy {
171        self.internal_energy(contributions) / self.total_moles
172    }
173
174    /// Gibbs energy: $G=A+pV$
175    pub fn gibbs_energy(&self, contributions: Contributions) -> Energy {
176        self.pressure(contributions) * self.volume + self.helmholtz_energy(contributions)
177    }
178
179    /// Molar Gibbs energy: $g=\frac{G}{N}$
180    pub fn molar_gibbs_energy(&self, contributions: Contributions) -> MolarEnergy {
181        self.gibbs_energy(contributions) / self.total_moles
182    }
183
184    /// Joule Thomson coefficient: $\mu_{JT}=\left(\frac{\partial T}{\partial p}\right)_{H,N_i}$
185    pub fn joule_thomson(&self) -> <Temperature as Div<Pressure>>::Output {
186        let c = Contributions::Total;
187        -(self.volume + self.temperature * self.dp_dt(c) / self.dp_dv(c))
188            / (self.total_moles * self.molar_isobaric_heat_capacity(c))
189    }
190
191    /// Isentropic compressibility: $\kappa_s=-\frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{S,N_i}$
192    pub fn isentropic_compressibility(&self) -> <f64 as Div<Pressure>>::Output {
193        let c = Contributions::Total;
194        -self.molar_isochoric_heat_capacity(c)
195            / (self.molar_isobaric_heat_capacity(c) * self.dp_dv(c) * self.volume)
196    }
197
198    /// Isenthalpic compressibility: $\kappa_H=-\frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{H,N_i}$
199    pub fn isenthalpic_compressibility(&self) -> <f64 as Div<Pressure>>::Output {
200        self.isentropic_compressibility() * (1.0 + self.grueneisen_parameter())
201    }
202
203    /// Thermal expansivity: $\alpha_p=-\frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_{p,N_i}$
204    pub fn thermal_expansivity(&self) -> <f64 as Div<Temperature>>::Output {
205        let c = Contributions::Total;
206        -self.dp_dt(c) / self.dp_dv(c) / self.volume
207    }
208
209    /// Grueneisen parameter: $\phi=V\left(\frac{\partial p}{\partial U}\right)_{V,n_i}=\frac{v}{c_v}\left(\frac{\partial p}{\partial T}\right)_{v,n_i}=\frac{\rho}{T}\left(\frac{\partial T}{\partial \rho}\right)_{s, n_i}$
210    pub fn grueneisen_parameter(&self) -> f64 {
211        let c = Contributions::Total;
212        (self.volume / (self.total_moles * self.molar_isochoric_heat_capacity(c)) * self.dp_dt(c))
213            .into_value()
214    }
215
216    /// Chemical potential $\mu_i$ evaluated for each contribution of the equation of state.
217    pub fn chemical_potential_contributions(
218        &self,
219        component: usize,
220        contributions: Contributions,
221    ) -> Vec<(String, MolarEnergy)> {
222        let new_state = self.derive1(DN(component));
223        let mut res = Vec::new();
224        if let Contributions::IdealGas | Contributions::Total = contributions {
225            res.push((
226                self.eos.ideal_gas_model(),
227                MolarEnergy::from_reduced(
228                    (self.eos.ideal_gas_helmholtz_energy(&new_state) * new_state.temperature).eps,
229                ),
230            ));
231        }
232        if let Contributions::Residual | Contributions::Total = contributions {
233            let contributions = self.eos.residual_helmholtz_energy_contributions(&new_state);
234            for (s, v) in contributions {
235                res.push((
236                    s,
237                    MolarEnergy::from_reduced((v * new_state.temperature).eps),
238                ));
239            }
240        }
241        res
242    }
243}
244
245impl<E: Residual + Molarweight + IdealGas> State<E> {
246    /// Specific isochoric heat capacity: $c_v^{(m)}=\frac{C_v}{m}$
247    pub fn specific_isochoric_heat_capacity(
248        &self,
249        contributions: Contributions,
250    ) -> SpecificEntropy {
251        self.molar_isochoric_heat_capacity(contributions) / self.total_molar_weight()
252    }
253
254    /// Specific isobaric heat capacity: $c_p^{(m)}=\frac{C_p}{m}$
255    pub fn specific_isobaric_heat_capacity(&self, contributions: Contributions) -> SpecificEntropy {
256        self.molar_isobaric_heat_capacity(contributions) / self.total_molar_weight()
257    }
258
259    /// Specific entropy: $s^{(m)}=\frac{S}{m}$
260    pub fn specific_entropy(&self, contributions: Contributions) -> SpecificEntropy {
261        self.molar_entropy(contributions) / self.total_molar_weight()
262    }
263
264    /// Specific enthalpy: $h^{(m)}=\frac{H}{m}$
265    pub fn specific_enthalpy(&self, contributions: Contributions) -> SpecificEnergy {
266        self.molar_enthalpy(contributions) / self.total_molar_weight()
267    }
268
269    /// Specific Helmholtz energy: $a^{(m)}=\frac{A}{m}$
270    pub fn specific_helmholtz_energy(&self, contributions: Contributions) -> SpecificEnergy {
271        self.molar_helmholtz_energy(contributions) / self.total_molar_weight()
272    }
273
274    /// Specific internal energy: $u^{(m)}=\frac{U}{m}$
275    pub fn specific_internal_energy(&self, contributions: Contributions) -> SpecificEnergy {
276        self.molar_internal_energy(contributions) / self.total_molar_weight()
277    }
278
279    /// Specific Gibbs energy: $g^{(m)}=\frac{G}{m}$
280    pub fn specific_gibbs_energy(&self, contributions: Contributions) -> SpecificEnergy {
281        self.molar_gibbs_energy(contributions) / self.total_molar_weight()
282    }
283
284    /// Speed of sound: $c=\sqrt{\left(\frac{\partial p}{\partial\rho^{(m)}}\right)_{S,N_i}}$
285    pub fn speed_of_sound(&self) -> Velocity {
286        (1.0 / (self.density * self.total_molar_weight() * self.isentropic_compressibility()))
287            .sqrt()
288    }
289}
feos_core/state/properties.rs

feos_core/state/
properties.rs
