Loading...
Searching...
No Matches
janafThermoI.H
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9 Copyright (C) 2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "janafThermo.H"
30#include "specie.H"
31
32// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
33
34template<class EquationOfState>
36(
37 const EquationOfState& st,
38 const scalar Tlow,
39 const scalar Thigh,
40 const scalar Tcommon,
41 const typename janafThermo<EquationOfState>::coeffArray& highCpCoeffs,
42 const typename janafThermo<EquationOfState>::coeffArray& lowCpCoeffs,
43 const bool convertCoeffs
44)
45:
46 EquationOfState(st),
47 Tlow_(Tlow),
48 Thigh_(Thigh),
49 Tcommon_(Tcommon)
50{
51 if (convertCoeffs)
52 {
53 for (label coefLabel=0; coefLabel<nCoeffs_; coefLabel++)
54 {
55 highCpCoeffs_[coefLabel] = highCpCoeffs[coefLabel]*this->R();
56 lowCpCoeffs_[coefLabel] = lowCpCoeffs[coefLabel]*this->R();
57 }
58 }
59 else
60 {
61 for (label coefLabel=0; coefLabel<nCoeffs_; coefLabel++)
62 {
63 highCpCoeffs_[coefLabel] = highCpCoeffs[coefLabel];
64 lowCpCoeffs_[coefLabel] = lowCpCoeffs[coefLabel];
65 }
66 }
67}
68
69
70template<class EquationOfState>
72Foam::janafThermo<EquationOfState>::coeffs
73(
74 const scalar T
75) const
76{
77 if (T < Tcommon_)
78 {
79 return lowCpCoeffs_;
80 }
81 else
82 {
83 return highCpCoeffs_;
84 }
85}
86
87
88// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
89
90template<class EquationOfState>
92(
93 const word& name,
94 const janafThermo& jt
95)
96:
97 EquationOfState(name, jt),
98 Tlow_(jt.Tlow_),
99 Thigh_(jt.Thigh_),
100 Tcommon_(jt.Tcommon_)
101{
102 for (label coefLabel=0; coefLabel<nCoeffs_; coefLabel++)
103 {
104 highCpCoeffs_[coefLabel] = jt.highCpCoeffs_[coefLabel];
105 lowCpCoeffs_[coefLabel] = jt.lowCpCoeffs_[coefLabel];
107}
108
109
110// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
111
112template<class EquationOfState>
114(
115 const scalar T
116) const
117{
118 if (T < Tlow_ || T > Thigh_)
119 {
121 << "attempt to use janafThermo<EquationOfState>"
122 " out of temperature range "
123 << Tlow_ << " -> " << Thigh_ << "; T = " << T
124 << endl;
125
126 return clamp(T, Tlow_, Thigh_);
128
129 return T;
130}
131
132
133template<class EquationOfState>
135{
136 return Tlow_;
137}
138
139
140template<class EquationOfState>
142{
143 return Thigh_;
144}
145
146
147template<class EquationOfState>
148inline Foam::scalar Foam::janafThermo<EquationOfState>::Tcommon() const
150 return Tcommon_;
151}
152
153
154template<class EquationOfState>
158 return highCpCoeffs_;
159}
160
161
162template<class EquationOfState>
165{
166 return lowCpCoeffs_;
167}
168
169
170template<class EquationOfState>
172(
173 const scalar p,
174 const scalar T
175) const
176{
177 const coeffArray& a = coeffs(T);
178 return
179 ((((a[4]*T + a[3])*T + a[2])*T + a[1])*T + a[0])
180 + EquationOfState::Cp(p, T);
181}
182
183
184template<class EquationOfState>
186(
187 const scalar p,
188 const scalar T
189) const
190{
191 const coeffArray& a = coeffs(T);
192 return
193 (
194 ((((a[4]/5.0*T + a[3]/4.0)*T + a[2]/3.0)*T + a[1]/2.0)*T + a[0])*T
195 + a[5]
196 ) + EquationOfState::H(p, T);
197}
198
199
200template<class EquationOfState>
202(
203 const scalar p,
204 const scalar T
205) const
206{
207 return Ha(p, T) - Hc();
208}
209
210
211template<class EquationOfState>
212inline Foam::scalar Foam::janafThermo<EquationOfState>::Hc() const
213{
214 const coeffArray& a = lowCpCoeffs_;
215 return
216 (
217 (
218 (((a[4]/5.0*Tstd + a[3]/4.0)*Tstd + a[2]/3.0)*Tstd + a[1]/2.0)*Tstd
219 + a[0]
220 )*Tstd + a[5]
221 );
222}
223
224
225template<class EquationOfState>
227(
228 const scalar p,
229 const scalar T
230) const
231{
232 const coeffArray& a = coeffs(T);
233 return
234 (
235 (((a[4]/4.0*T + a[3]/3.0)*T + a[2]/2.0)*T + a[1])*T + a[0]*log(T)
236 + a[6]
237 ) + EquationOfState::S(p, T);
238}
239
240
241template<class EquationOfState>
243(
244 const scalar T
245) const
246{
247 const coeffArray& a = coeffs(T);
248 return
249 (
250 (
251 a[0]*(1 - log(T))
252 - (((a[4]/20.0*T + a[3]/12.0)*T + a[2]/6.0)*T + a[1]/2.0)*T
253 - a[6]
255 + a[5]
256 );
257}
258
259
260template<class EquationOfState>
262(
263 const scalar p,
264 const scalar T
265) const
266{
267 const coeffArray& a = coeffs(T);
268 return
269 (((4*a[4]*T + 3*a[3])*T + 2*a[2])*T + a[1]);
270}
271
272// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
273
274template<class EquationOfState>
276(
278)
279{
280 scalar Y1 = this->Y();
281
282 EquationOfState::operator+=(jt);
283
284 if (mag(this->Y()) > SMALL)
285 {
286 Y1 /= this->Y();
287 const scalar Y2 = jt.Y()/this->Y();
288
289 Tlow_ = max(Tlow_, jt.Tlow_);
290 Thigh_ = min(Thigh_, jt.Thigh_);
291
292 if
293 (
295 && !equal(Tcommon_, jt.Tcommon_)
296 )
297 {
299 << "Tcommon " << Tcommon_ << " for "
300 << (this->name().size() ? this->name() : "others")
301 << " != " << jt.Tcommon_ << " for "
302 << (jt.name().size() ? jt.name() : "others")
303 << exit(FatalError);
304 }
305
306 for
307 (
308 label coefLabel=0;
309 coefLabel<janafThermo<EquationOfState>::nCoeffs_;
310 coefLabel++
311 )
312 {
313 highCpCoeffs_[coefLabel] =
314 Y1*highCpCoeffs_[coefLabel]
315 + Y2*jt.highCpCoeffs_[coefLabel];
316
317 lowCpCoeffs_[coefLabel] =
318 Y1*lowCpCoeffs_[coefLabel]
319 + Y2*jt.lowCpCoeffs_[coefLabel];
320 }
321 }
322}
323
324
325// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
326
327template<class EquationOfState>
328inline Foam::janafThermo<EquationOfState> Foam::operator+
329(
332)
333{
334 EquationOfState eofs = jt1;
335 eofs += jt2;
336
337 if (mag(eofs.Y()) < SMALL)
338 {
340 (
341 eofs,
342 jt1.Tlow_,
343 jt1.Thigh_,
344 jt1.Tcommon_,
345 jt1.highCpCoeffs_,
346 jt1.lowCpCoeffs_
347 );
348 }
349 else
350 {
351 const scalar Y1 = jt1.Y()/eofs.Y();
352 const scalar Y2 = jt2.Y()/eofs.Y();
353
354 typename janafThermo<EquationOfState>::coeffArray highCpCoeffs;
355 typename janafThermo<EquationOfState>::coeffArray lowCpCoeffs;
356
357 for
358 (
359 label coefLabel=0;
360 coefLabel<janafThermo<EquationOfState>::nCoeffs_;
361 coefLabel++
362 )
363 {
364 highCpCoeffs[coefLabel] =
365 Y1*jt1.highCpCoeffs_[coefLabel]
366 + Y2*jt2.highCpCoeffs_[coefLabel];
367
368 lowCpCoeffs[coefLabel] =
369 Y1*jt1.lowCpCoeffs_[coefLabel]
370 + Y2*jt2.lowCpCoeffs_[coefLabel];
371 }
372
373 if
374 (
376 && !equal(jt1.Tcommon_, jt2.Tcommon_)
377 )
378 {
380 << "Tcommon " << jt1.Tcommon_ << " for "
381 << (jt1.name().size() ? jt1.name() : "others")
382 << " != " << jt2.Tcommon_ << " for "
383 << (jt2.name().size() ? jt2.name() : "others")
384 << exit(FatalError);
385 }
386
388 (
389 eofs,
390 max(jt1.Tlow_, jt2.Tlow_),
391 min(jt1.Thigh_, jt2.Thigh_),
392 jt1.Tcommon_,
393 highCpCoeffs,
394 lowCpCoeffs
395 );
396 }
397}
398
399
400template<class EquationOfState>
401inline Foam::janafThermo<EquationOfState> Foam::operator*
402(
403 const scalar s,
405)
406{
408 (
409 s*static_cast<const EquationOfState&>(jt),
410 jt.Tlow_,
411 jt.Thigh_,
412 jt.Tcommon_,
413 jt.highCpCoeffs_,
414 jt.lowCpCoeffs_
415 );
416}
417
418
419template<class EquationOfState>
420inline Foam::janafThermo<EquationOfState> Foam::operator==
421(
424)
425{
426 EquationOfState eofs
427 (
428 static_cast<const EquationOfState&>(jt1)
429 == static_cast<const EquationOfState&>(jt2)
430 );
431
432 const scalar Y1 = jt2.Y()/eofs.Y();
433 const scalar Y2 = jt1.Y()/eofs.Y();
434
435 typename janafThermo<EquationOfState>::coeffArray highCpCoeffs;
436 typename janafThermo<EquationOfState>::coeffArray lowCpCoeffs;
437
438 for
439 (
440 label coefLabel=0;
441 coefLabel<janafThermo<EquationOfState>::nCoeffs_;
442 coefLabel++
443 )
444 {
445 highCpCoeffs[coefLabel] =
446 Y1*jt2.highCpCoeffs_[coefLabel]
447 - Y2*jt1.highCpCoeffs_[coefLabel];
448
449 lowCpCoeffs[coefLabel] =
450 Y1*jt2.lowCpCoeffs_[coefLabel]
451 - Y2*jt1.lowCpCoeffs_[coefLabel];
452 }
453
454 if
455 (
457 && !equal(jt2.Tcommon_, jt1.Tcommon_)
458 )
459 {
461 << "Tcommon " << jt2.Tcommon_ << " for "
462 << (jt2.name().size() ? jt2.name() : "others")
463 << " != " << jt1.Tcommon_ << " for "
464 << (jt1.name().size() ? jt1.name() : "others")
465 << exit(FatalError);
466 }
467
469 (
470 eofs,
471 max(jt2.Tlow_, jt1.Tlow_),
472 min(jt2.Thigh_, jt1.Thigh_),
473 jt2.Tcommon_,
474 highCpCoeffs,
475 lowCpCoeffs
476 );
477}
478
479
480// ************************************************************************* //
scalar Ha(const scalar p, const scalar T) const
Definition EtoHthermo.H:32
#define R(A, B, C, D, E, F, K, M)
JANAF tables based thermodynamics package templated into the equation of state.
Definition janafThermo.H:90
static constexpr int nCoeffs_
Definition janafThermo.H:95
janafThermo(const EquationOfState &st, const scalar Tlow, const scalar Thigh, const scalar Tcommon, const coeffArray &highCpCoeffs, const coeffArray &lowCpCoeffs, const bool convertCoeffs=false)
Construct from components.
scalar limit(const scalar T) const
Limit the temperature to be in the range Tlow_ to Thigh_.
scalar Hs(const scalar p, const scalar T) const
Sensible enthalpy [J/kg].
scalar Tlow() const
Return const access to the low temperature limit.
scalar dCpdT(const scalar p, const scalar T) const
Temperature derivative of heat capacity at constant pressure.
scalar Cp(const scalar p, const scalar T) const
Heat capacity at constant pressure [J/(kg K)].
FixedList< scalar, nCoeffs_ > coeffArray
Definition janafThermo.H:96
scalar Hc() const
Chemical enthalpy [J/kg].
scalar Tcommon() const
Return const access to the common temperature.
scalar S(const scalar p, const scalar T) const
Entropy [J/(kg K)].
scalar Ha(const scalar p, const scalar T) const
Absolute Enthalpy [J/kg].
scalar Thigh() const
Return const access to the high temperature limit.
const coeffArray & highCpCoeffs() const
Return const access to the high temperature poly coefficients.
scalar Gstd(const scalar T) const
Gibbs free energy of the mixture in the standard state [J/kg].
const coeffArray & lowCpCoeffs() const
Return const access to the low temperature poly coefficients.
A class for handling words, derived from Foam::string.
Definition word.H:66
volScalarField & p
PtrList< volScalarField > & Y
const volScalarField & T
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition error.H:600
auto & name
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
#define WarningInFunction
Report a warning using Foam::Warning.
const scalar Tstd
Standard temperature: default in [K].
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition hashSets.C:40
bool equal(const T &a, const T &b)
Compare two values for equality.
Definition label.H:180
dimensionSet clamp(const dimensionSet &a, const dimensionSet &range)
dimensionedScalar log(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition Ostream.H:519
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition hashSets.C:26
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
Definition exprTraits.C:127
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition errorManip.H:125
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)