Skip to content
This repository was archived by the owner on Sep 14, 2022. It is now read-only.

Commit 1c7bfd9

Browse files
g-bauerg-bauer
committed
refactored reference for thermal conductivity
1 parent b34e004 commit 1c7bfd9

File tree

1 file changed

+82
-63
lines changed

1 file changed

+82
-63
lines changed

‎src/eos/mod.rs‎

Lines changed: 82 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -161,6 +161,23 @@ fn omega22(t: f64) -> f64 {
161161
- 6.435e-4 * t.powf(0.14874) * (18.0323 * t.powf(-0.76830) - 7.27371).sin()
162162
}
163163

164+
#[inline]
165+
fn chapman_enskog_thermal_conductivity(
166+
temperature: SINumber,
167+
molarweight: SINumber,
168+
m: f64,
169+
sigma: f64,
170+
epsilon_k: f64,
171+
) -> SINumber {
172+
let t = temperature.to_reduced(KELVIN).unwrap();
173+
0.083235 * (t * m / molarweight.to_reduced(GRAM / MOL).unwrap()).sqrt()
174+
/ sigma.powi(2)
175+
/ omega22(t / epsilon_k)
176+
* WATT
177+
/ METER
178+
/ KELVIN
179+
}
180+
164181
impl EntropyScaling<SIUnit> for PcSaft {
165182
fn viscosity_reference(
166183
&self,
@@ -261,48 +278,7 @@ impl EntropyScaling<SIUnit> for PcSaft {
261278
Ok(a + b * s - c * (1.0 - s.exp()) * s.powi(2) - d * s.powi(4) - e * s.powi(8))
262279
}
263280

264-
// fn thermal_conductivity_reference(
265-
// &self,
266-
// state: &State<SIUnit, E>,
267-
// ) -> EosResult<SINumber> {
268-
// if self.components() != 1 {
269-
// return Err(EosError::IncompatibleComponents(self.components(), 1));
270-
// }
271-
// let p = &self.parameters;
272-
// let res: Array1<SINumber> = (0..self.components())
273-
// .map(|i| {
274-
// let tr = (temperature / p.epsilon_k[i] / KELVIN)
275-
// .into_value()
276-
// .unwrap();
277-
// let cp = State::critical_point_pure(&state.eos, Some(state.temperature)).unwrap();
278-
// let s_res_cp_reduced = cp
279-
// .entropy(Contributions::ResidualNvt)
280-
// .to_reduced(SIUnit::reference_entropy())
281-
// .unwrap();
282-
// let s_res_reduced = cp
283-
// .entropy(Contributions::ResidualNvt)
284-
// .to_reduced(SIUnit::reference_entropy())
285-
// .unwrap();
286-
// let ref_ce = 0.083235
287-
// * ((temperature / KELVIN).into_value().unwrap()
288-
// / (p.molarweight[0]
289-
// / p.m[0]))
290-
// .sqrt()
291-
// / p.sigma[0]
292-
// / p.sigma[0]
293-
// / omega22(tr)
294-
// * p.m[0];
295-
// let alpha_visc = (-s_res_reduced / s_res_cp_reduced).exp();
296-
// let ref_ts = (-0.0167141 * tr / p.m[0] + 0.0470581 * (tr / p.m[0]).powi(2))
297-
// * (p.m[0] * p.m[0] * p.sigma[i].powi(3) * p.epsilon_k[0])
298-
// / 100000.0;
299-
// (ref_ce + ref_ts * alpha_visc) * WATT / METER / KELVIN
300-
// })
301-
// .collect();
302-
// Ok(res[0])
303-
// }
304-
305-
// Equation 11 of DOI: 10.1021/acs.iecr.9b03998
281+
// Equation 4 of DOI: 10.1021/acs.iecr.9b04289
306282
fn thermal_conductivity_reference(
307283
&self,
308284
temperature: SINumber,
@@ -313,38 +289,81 @@ impl EntropyScaling<SIUnit> for PcSaft {
313289
return Err(EosError::IncompatibleComponents(self.components(), 1));
314290
}
315291
let p = &self.parameters;
292+
let mws = self.molar_weight();
316293
let state = State::new_nvt(&Rc::new(Self::new(p.clone())), temperature, volume, moles)?;
317294
let res: Array1<SINumber> = (0..self.components())
318295
.map(|i| {
319296
let tr = (temperature / p.epsilon_k[i] / KELVIN)
320297
.into_value()
321298
.unwrap();
322-
let ce = 83.235
323-
* f64::powf(10.0, -1.5)
324-
* ((temperature / KELVIN).into_value().unwrap() / p.molarweight[0] * p.m[0])
325-
.sqrt()
326-
/ (p.sigma[0] * p.sigma[0])
327-
/ omega22(tr);
328-
ce * WATT / METER / KELVIN
329-
+ state.density
330-
* self
331-
.diffusion_reference(temperature, volume, moles)
332-
.unwrap()
333-
* self
334-
.diffusion_correlation(
335-
state
336-
.molar_entropy(Contributions::ResidualNvt)
337-
.to_reduced(SIUnit::reference_molar_entropy())
338-
.unwrap(),
339-
&state.molefracs,
340-
)
341-
.unwrap()
342-
* (state.c_v(Contributions::Total) - 1.5 * RGAS)
299+
let s_res_reduced = state
300+
.molar_entropy(Contributions::ResidualNvt)
301+
.to_reduced(RGAS)
302+
.unwrap()
303+
/ p.m[i];
304+
let ref_ce = chapman_enskog_thermal_conductivity(
305+
temperature,
306+
mws.get(i),
307+
p.m[i],
308+
p.sigma[i],
309+
p.epsilon_k[i],
310+
);
311+
let alpha_visc = (-s_res_reduced / -0.5).exp();
312+
let ref_ts = (-0.0167141 * tr / p.m[i] + 0.0470581 * (tr / p.m[i]).powi(2))
313+
* (p.m[i] * p.m[i] * p.sigma[i].powi(3) * p.epsilon_k[i])
314+
* 1e-5
315+
* WATT
316+
/ METER
317+
/ KELVIN;
318+
ref_ce + ref_ts * alpha_visc
343319
})
344320
.collect();
345321
Ok(res[0])
346322
}
347323

324+
// // Equation 11 of DOI: 10.1021/acs.iecr.9b03998
325+
// fn thermal_conductivity_reference(
326+
// &self,
327+
// temperature: SINumber,
328+
// volume: SINumber,
329+
// moles: &SIArray1,
330+
// ) -> EosResult<SINumber> {
331+
// if self.components() != 1 {
332+
// return Err(EosError::IncompatibleComponents(self.components(), 1));
333+
// }
334+
// let p = &self.parameters;
335+
// let state = State::new_nvt(&Rc::new(Self::new(p.clone())), temperature, volume, moles)?;
336+
// let res: Array1<SINumber> = (0..self.components())
337+
// .map(|i| {
338+
// let tr = (temperature / p.epsilon_k[i] / KELVIN)
339+
// .into_value()
340+
// .unwrap();
341+
// let ce = 83.235
342+
// * f64::powf(10.0, -1.5)
343+
// * ((temperature / KELVIN).into_value().unwrap() / p.molarweight[0] * p.m[0])
344+
// .sqrt()
345+
// / (p.sigma[0] * p.sigma[0])
346+
// / omega22(tr);
347+
// ce * WATT / METER / KELVIN
348+
// + state.density
349+
// * self
350+
// .diffusion_reference(temperature, volume, moles)
351+
// .unwrap()
352+
// * self
353+
// .diffusion_correlation(
354+
// state
355+
// .molar_entropy(Contributions::ResidualNvt)
356+
// .to_reduced(SIUnit::reference_molar_entropy())
357+
// .unwrap(),
358+
// &state.molefracs,
359+
// )
360+
// .unwrap()
361+
// * (state.c_v(Contributions::Total) - 1.5 * RGAS)
362+
// })
363+
// .collect();
364+
// Ok(res[0])
365+
// }
366+
348367
fn thermal_conductivity_correlation(&self, s_res: f64, x: &Array1<f64>) -> EosResult<f64> {
349368
if self.components() != 1 {
350369
return Err(EosError::IncompatibleComponents(self.components(), 1));

0 commit comments

Comments
 (0)