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 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 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 pub fn molar_isochoric_heat_capacity(&self, contributions: Contributions) -> MolarEntropy<D> {
52 self.temperature * self.ds_dt(contributions) / self.total_moles
53 }
54
55 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 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 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 pub fn molar_entropy(&self, contributions: Contributions) -> MolarEntropy<D> {
96 self.entropy(contributions) / self.total_moles
97 }
98
99 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 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 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 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 pub fn molar_enthalpy(&self, contributions: Contributions) -> MolarEnergy<D> {
154 self.enthalpy(contributions) / self.total_moles
155 }
156
157 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 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 pub fn molar_helmholtz_energy(&self, contributions: Contributions) -> MolarEnergy<D> {
182 self.helmholtz_energy(contributions) / self.total_moles
183 }
184
185 pub fn internal_energy(&self, contributions: Contributions) -> Energy<D> {
187 self.temperature * self.entropy(contributions) + self.helmholtz_energy(contributions)
188 }
189
190 pub fn molar_internal_energy(&self, contributions: Contributions) -> MolarEnergy<D> {
192 self.internal_energy(contributions) / self.total_moles
193 }
194
195 pub fn gibbs_energy(&self, contributions: Contributions) -> Energy<D> {
197 self.pressure(contributions) * self.volume + self.helmholtz_energy(contributions)
198 }
199
200 pub fn molar_gibbs_energy(&self, contributions: Contributions) -> MolarEnergy<D> {
202 self.gibbs_energy(contributions) / self.total_moles
203 }
204
205 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 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 pub fn isenthalpic_compressibility(&self) -> InvP<D> {
221 self.isentropic_compressibility() * Dimensionless::new(self.grueneisen_parameter() + 1.0)
222 }
223
224 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 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 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 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 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 pub fn specific_entropy(&self, contributions: Contributions) -> SpecificEntropy<D> {
288 self.molar_entropy(contributions) / self.total_molar_weight()
289 }
290
291 pub fn specific_enthalpy(&self, contributions: Contributions) -> SpecificEnergy<D> {
293 self.molar_enthalpy(contributions) / self.total_molar_weight()
294 }
295
296 pub fn specific_helmholtz_energy(&self, contributions: Contributions) -> SpecificEnergy<D> {
298 self.molar_helmholtz_energy(contributions) / self.total_molar_weight()
299 }
300
301 pub fn specific_internal_energy(&self, contributions: Contributions) -> SpecificEnergy<D> {
303 self.molar_internal_energy(contributions) / self.total_molar_weight()
304 }
305
306 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 pub fn speed_of_sound(&self) -> Velocity {
318 (self.density * self.total_molar_weight() * self.isentropic_compressibility())
319 .inv()
320 .sqrt()
321 }
322}