-
Notifications
You must be signed in to change notification settings - Fork 604
/
readInitialConditions.H
103 lines (89 loc) · 2.51 KB
/
readInitialConditions.H
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
const word constProp(initialConditions.lookup("constantProperty"));
if
(
(constProp != "pressure")
&& (constProp != "volume")
&& (constProp != "temperature")
)
{
FatalError
<< "in initialConditions, unknown constantProperty type "
<< constProp << nl
<< " Valid types are: pressure volume temperature."
<< exit(FatalError);
}
const word fractionBasis(initialConditions.lookup("fractionBasis"));
if ((fractionBasis != "mass") && (fractionBasis != "mole"))
{
FatalError << "in initialConditions, unknown fractionBasis type " << nl
<< "Valid types are: mass or mole."
<< fractionBasis << exit(FatalError);
}
const label nSpecie = Y.size();
const scalarList W(::W(thermo));
scalarList Y0(nSpecie, 0.0);
scalarList X0(nSpecie, 0.0);
dictionary fractions(initialConditions.subDict("fractions"));
if (fractionBasis == "mole")
{
forAll(Y, i)
{
const word& name = Y[i].name();
if (fractions.found(name))
{
X0[i] = readScalar(fractions.lookup(name));
}
}
scalar mw = 0.0;
const scalar mTot = sum(X0);
forAll(Y, i)
{
X0[i] /= mTot;
mw += W[i]*X0[i];
}
forAll(Y, i)
{
Y0[i] = X0[i]*W[i]/mw;
}
}
else // mass fraction
{
forAll(Y, i)
{
const word& name = Y[i].name();
if (fractions.found(name))
{
Y0[i] = readScalar(fractions.lookup(name));
}
}
scalar invW = 0.0;
const scalar mTot = sum(Y0);
forAll(Y, i)
{
Y0[i] /= mTot;
invW += Y0[i]/W[i];
}
const scalar mw = 1.0/invW;
forAll(Y, i)
{
X0[i] = Y0[i]*mw/W[i];
}
}
const scalar h0 = ::h0(thermo, Y0, p[0], T0);
forAll(Y, i)
{
Y[i] = Y0[i];
}
thermo.he() = dimensionedScalar(dimEnergy/dimMass, h0);
thermo.correct();
rho = thermo.rho();
scalar rho0 = rho[0];
scalar u0 = h0 - p0/rho0;
scalar R0 = p0/(rho0*T0);
Rspecific[0] = R0;
scalar integratedHeat = 0.0;
Info << constProp << " will be held constant." << nl
<< " p = " << p[0] << " [Pa]" << nl
<< " T = " << thermo.T()[0] << " [K] " << nl
<< " rho = " << rho[0] << " [kg/m3]" << nl
<< endl;