feos_core/state/
properties.rs

1use super::{Contributions, State};
2use crate::equation_of_state::{Molarweight, Total};
3use crate::{ReferenceSystem, Residual};
4use nalgebra::allocator::Allocator;
5use nalgebra::{DefaultAllocator, OVector};
6use num_dual::{Dual, DualNum, Gradients, partial, partial2};
7use quantity::*;
8use std::ops::{Div, Neg};
9
10type InvP<T> = Quantity<T, <_Pressure as Neg>::Output>;
11type InvT<T> = Quantity<T, <_Temperature as Neg>::Output>;
12
13impl<E: Total<N, D>, N: Gradients, D: DualNum<f64> + Copy> State<E, N, D>
14where
15    DefaultAllocator: Allocator<N>,
16{
17    /// Chemical potential: $\mu_i=\left(\frac{\partial A}{\partial N_i}\right)_{T,V,N_j}$
18    pub fn chemical_potential(&self, contributions: Contributions) -> MolarEnergy<OVector<D, N>> {
19        let residual = || self.residual_chemical_potential();
20        let ideal_gas = || {
21            quantity::ad::gradient_copy(
22                partial2(
23                    |n, &t, &v| self.eos.ideal_gas_helmholtz_energy(t, v, &n),
24                    &self.temperature,
25                    &self.volume,
26                ),
27                &self.moles,
28            )
29            .1
30        };
31        Self::contributions(ideal_gas, residual, contributions)
32    }
33
34    /// Partial derivative of chemical potential w.r.t. temperature: $\left(\frac{\partial\mu_i}{\partial T}\right)_{V,N_i}$
35    pub fn dmu_dt(&self, contributions: Contributions) -> MolarEntropy<OVector<D, N>> {
36        let residual = || self.dmu_res_dt();
37        let ideal_gas = || {
38            quantity::ad::partial_hessian_copy(
39                partial(
40                    |(n, t), &v| self.eos.ideal_gas_helmholtz_energy(t, v, &n),
41                    &self.volume,
42                ),
43                (&self.moles, self.temperature),
44            )
45            .3
46        };
47        Self::contributions(ideal_gas, residual, contributions)
48    }
49
50    /// Molar isochoric heat capacity: $c_v=\left(\frac{\partial u}{\partial T}\right)_{V,N_i}$
51    pub fn molar_isochoric_heat_capacity(&self, contributions: Contributions) -> MolarEntropy<D> {
52        self.temperature * self.ds_dt(contributions) / self.total_moles
53    }
54
55    /// Partial derivative of the molar isochoric heat capacity w.r.t. temperature: $\left(\frac{\partial c_V}{\partial T}\right)_{V,N_i}$
56    pub fn dc_v_dt(
57        &self,
58        contributions: Contributions,
59    ) -> <MolarEntropy<D> as Div<Temperature<D>>>::Output {
60        (self.temperature * self.d2s_dt2(contributions) + self.ds_dt(contributions))
61            / self.total_moles
62    }
63
64    /// Molar isobaric heat capacity: $c_p=\left(\frac{\partial h}{\partial T}\right)_{p,N_i}$
65    pub fn molar_isobaric_heat_capacity(&self, contributions: Contributions) -> MolarEntropy<D> {
66        match contributions {
67            Contributions::Residual => self.residual_molar_isobaric_heat_capacity(),
68            _ => {
69                self.temperature / self.total_moles
70                    * (self.ds_dt(contributions)
71                        - (self.dp_dt(contributions) * self.dp_dt(contributions))
72                            / self.dp_dv(contributions))
73            }
74        }
75    }
76
77    /// Entropy: $S=-\left(\frac{\partial A}{\partial T}\right)_{V,N_i}$
78    pub fn entropy(&self, contributions: Contributions) -> Entropy<D> {
79        let residual = || self.residual_entropy();
80        let ideal_gas = || {
81            -quantity::ad::first_derivative(
82                partial2(
83                    |t, &v, n| self.eos.ideal_gas_helmholtz_energy(t, v, n),
84                    &self.volume,
85                    &self.moles,
86                ),
87                self.temperature,
88            )
89            .1
90        };
91        Self::contributions(ideal_gas, residual, contributions)
92    }
93
94    /// Molar entropy: $s=\frac{S}{N}$
95    pub fn molar_entropy(&self, contributions: Contributions) -> MolarEntropy<D> {
96        self.entropy(contributions) / self.total_moles
97    }
98
99    /// Partial molar entropy: $s_i=\left(\frac{\partial S}{\partial N_i}\right)_{T,p,N_j}$
100    pub fn partial_molar_entropy(&self) -> MolarEntropy<OVector<D, N>> {
101        let c = Contributions::Total;
102        -(self.dmu_dt(c) + self.dp_dni(c) * (self.dp_dt(c) / self.dp_dv(c)))
103    }
104
105    /// Partial derivative of the entropy w.r.t. temperature: $\left(\frac{\partial S}{\partial T}\right)_{V,N_i}$
106    pub fn ds_dt(
107        &self,
108        contributions: Contributions,
109    ) -> <Entropy<D> as Div<Temperature<D>>>::Output {
110        let residual = || self.ds_res_dt();
111        let ideal_gas = || {
112            -quantity::ad::second_derivative(
113                partial2(
114                    |t, &v, n| self.eos.ideal_gas_helmholtz_energy(t, v, n),
115                    &self.volume,
116                    &self.moles,
117                ),
118                self.temperature,
119            )
120            .2
121        };
122        Self::contributions(ideal_gas, residual, contributions)
123    }
124
125    /// Second partial derivative of the entropy w.r.t. temperature: $\left(\frac{\partial^2 S}{\partial T^2}\right)_{V,N_i}$
126    pub fn d2s_dt2(
127        &self,
128        contributions: Contributions,
129    ) -> <<Entropy<D> as Div<Temperature<D>>>::Output as Div<Temperature<D>>>::Output {
130        let residual = || self.d2s_res_dt2();
131        let ideal_gas = || {
132            -quantity::ad::third_derivative(
133                partial2(
134                    |t, &v, n| self.eos.ideal_gas_helmholtz_energy(t, v, n),
135                    &self.volume,
136                    &self.moles,
137                ),
138                self.temperature,
139            )
140            .3
141        };
142        Self::contributions(ideal_gas, residual, contributions)
143    }
144
145    /// Enthalpy: $H=A+TS+pV$
146    pub fn enthalpy(&self, contributions: Contributions) -> Energy<D> {
147        self.temperature * self.entropy(contributions)
148            + self.helmholtz_energy(contributions)
149            + self.pressure(contributions) * self.volume
150    }
151
152    /// Molar enthalpy: $h=\frac{H}{N}$
153    pub fn molar_enthalpy(&self, contributions: Contributions) -> MolarEnergy<D> {
154        self.enthalpy(contributions) / self.total_moles
155    }
156
157    /// Partial molar enthalpy: $h_i=\left(\frac{\partial H}{\partial N_i}\right)_{T,p,N_j}$
158    pub fn partial_molar_enthalpy(&self) -> MolarEnergy<OVector<D, N>> {
159        let s = self.partial_molar_entropy();
160        let mu = self.chemical_potential(Contributions::Total);
161        s * self.temperature + mu
162    }
163
164    /// Helmholtz energy: $A$
165    pub fn helmholtz_energy(&self, contributions: Contributions) -> Energy<D> {
166        let residual = || self.residual_helmholtz_energy();
167        let ideal_gas = || {
168            quantity::ad::zeroth_derivative(
169                partial2(
170                    |t, &v, n| self.eos.ideal_gas_helmholtz_energy(t, v, n),
171                    &self.volume,
172                    &self.moles,
173                ),
174                self.temperature,
175            )
176        };
177        Self::contributions(ideal_gas, residual, contributions)
178    }
179
180    /// Molar Helmholtz energy: $a=\frac{A}{N}$
181    pub fn molar_helmholtz_energy(&self, contributions: Contributions) -> MolarEnergy<D> {
182        self.helmholtz_energy(contributions) / self.total_moles
183    }
184
185    /// Internal energy: $U=A+TS$
186    pub fn internal_energy(&self, contributions: Contributions) -> Energy<D> {
187        self.temperature * self.entropy(contributions) + self.helmholtz_energy(contributions)
188    }
189
190    /// Molar internal energy: $u=\frac{U}{N}$
191    pub fn molar_internal_energy(&self, contributions: Contributions) -> MolarEnergy<D> {
192        self.internal_energy(contributions) / self.total_moles
193    }
194
195    /// Gibbs energy: $G=A+pV$
196    pub fn gibbs_energy(&self, contributions: Contributions) -> Energy<D> {
197        self.pressure(contributions) * self.volume + self.helmholtz_energy(contributions)
198    }
199
200    /// Molar Gibbs energy: $g=\frac{G}{N}$
201    pub fn molar_gibbs_energy(&self, contributions: Contributions) -> MolarEnergy<D> {
202        self.gibbs_energy(contributions) / self.total_moles
203    }
204
205    /// Joule Thomson coefficient: $\mu_{JT}=\left(\frac{\partial T}{\partial p}\right)_{H,N_i}$
206    pub fn joule_thomson(&self) -> <Temperature<D> as Div<Pressure<D>>>::Output {
207        let c = Contributions::Total;
208        -(self.volume + self.temperature * self.dp_dt(c) / self.dp_dv(c))
209            / (self.total_moles * self.molar_isobaric_heat_capacity(c))
210    }
211
212    /// Isentropic compressibility: $\kappa_s=-\frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{S,N_i}$
213    pub fn isentropic_compressibility(&self) -> InvP<D> {
214        let c = Contributions::Total;
215        -self.molar_isochoric_heat_capacity(c)
216            / (self.molar_isobaric_heat_capacity(c) * self.dp_dv(c) * self.volume)
217    }
218
219    /// Isenthalpic compressibility: $\kappa_H=-\frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{H,N_i}$
220    pub fn isenthalpic_compressibility(&self) -> InvP<D> {
221        self.isentropic_compressibility() * Dimensionless::new(self.grueneisen_parameter() + 1.0)
222    }
223
224    /// Thermal expansivity: $\alpha_p=-\frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_{p,N_i}$
225    pub fn thermal_expansivity(&self) -> InvT<D> {
226        let c = Contributions::Total;
227        -self.dp_dt(c) / self.dp_dv(c) / self.volume
228    }
229
230    /// 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}$
231    pub fn grueneisen_parameter(&self) -> D {
232        let c = Contributions::Total;
233        (self.dp_dt(c) / (self.molar_isochoric_heat_capacity(c) * self.density)).into_value()
234    }
235
236    /// Chemical potential $\mu_i$ evaluated for each contribution of the equation of state.
237    pub fn chemical_potential_contributions(
238        &self,
239        component: usize,
240        contributions: Contributions,
241    ) -> Vec<(&'static str, MolarEnergy<D>)> {
242        let t = Dual::from_re(self.temperature.into_reduced());
243        let v = Dual::from_re(self.density.into_reduced().recip());
244        let mut x = self.molefracs.map(Dual::from_re);
245        x[component].eps = D::one();
246        let mut res = Vec::new();
247        if let Contributions::IdealGas | Contributions::Total = contributions {
248            res.push((
249                self.eos.ideal_gas_model(),
250                self.eos.ideal_gas_molar_helmholtz_energy(t, v, &x),
251            ));
252        }
253        if let Contributions::Residual | Contributions::Total = contributions {
254            res.extend(
255                self.eos
256                    .lift()
257                    .molar_helmholtz_energy_contributions(t, v, &x),
258            );
259        }
260        res.into_iter()
261            .map(|(s, v)| (s, MolarEnergy::from_reduced(v.eps)))
262            .collect()
263    }
264}
265
266impl<E: Total<N, D> + Molarweight<N, D>, N: Gradients, D: DualNum<f64> + Copy> State<E, N, D>
267where
268    DefaultAllocator: Allocator<N>,
269{
270    /// Specific isochoric heat capacity: $c_v^{(m)}=\frac{C_v}{m}$
271    pub fn specific_isochoric_heat_capacity(
272        &self,
273        contributions: Contributions,
274    ) -> SpecificEntropy<D> {
275        self.molar_isochoric_heat_capacity(contributions) / self.total_molar_weight()
276    }
277
278    /// Specific isobaric heat capacity: $c_p^{(m)}=\frac{C_p}{m}$
279    pub fn specific_isobaric_heat_capacity(
280        &self,
281        contributions: Contributions,
282    ) -> SpecificEntropy<D> {
283        self.molar_isobaric_heat_capacity(contributions) / self.total_molar_weight()
284    }
285
286    /// Specific entropy: $s^{(m)}=\frac{S}{m}$
287    pub fn specific_entropy(&self, contributions: Contributions) -> SpecificEntropy<D> {
288        self.molar_entropy(contributions) / self.total_molar_weight()
289    }
290
291    /// Specific enthalpy: $h^{(m)}=\frac{H}{m}$
292    pub fn specific_enthalpy(&self, contributions: Contributions) -> SpecificEnergy<D> {
293        self.molar_enthalpy(contributions) / self.total_molar_weight()
294    }
295
296    /// Specific Helmholtz energy: $a^{(m)}=\frac{A}{m}$
297    pub fn specific_helmholtz_energy(&self, contributions: Contributions) -> SpecificEnergy<D> {
298        self.molar_helmholtz_energy(contributions) / self.total_molar_weight()
299    }
300
301    /// Specific internal energy: $u^{(m)}=\frac{U}{m}$
302    pub fn specific_internal_energy(&self, contributions: Contributions) -> SpecificEnergy<D> {
303        self.molar_internal_energy(contributions) / self.total_molar_weight()
304    }
305
306    /// Specific Gibbs energy: $g^{(m)}=\frac{G}{m}$
307    pub fn specific_gibbs_energy(&self, contributions: Contributions) -> SpecificEnergy<D> {
308        self.molar_gibbs_energy(contributions) / self.total_molar_weight()
309    }
310}
311
312impl<E: Total<N> + Molarweight<N>, N: Gradients> State<E, N>
313where
314    DefaultAllocator: Allocator<N>,
315{
316    /// Speed of sound: $c=\sqrt{\left(\frac{\partial p}{\partial\rho^{(m)}}\right)_{S,N_i}}$
317    pub fn speed_of_sound(&self) -> Velocity {
318        (self.density * self.total_molar_weight() * self.isentropic_compressibility())
319            .inv()
320            .sqrt()
321    }
322}