alphaEqn.H
Go to the documentation of this file.
1{
2 word alphaScheme("div(phi,alpha)");
3 word alpharScheme("div(phirb,alpha)");
4
5 surfaceScalarField phir
6 (
7 IOobject
8 (
9 "phir",
10 runTime.timeName(),
11 mesh,
12 IOobject::NO_READ,
13 IOobject::NO_WRITE
14 ),
15 mixture.cAlpha()*mag(phi/mesh.magSf())*mixture.nHatf()
16 );
17
18 for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
19 {
20 // Create the limiter to be used for all phase-fractions
21 scalarField allLambda(mesh.nFaces(), 1.0);
22
23 // Split the limiter into a surfaceScalarField
24 slicedSurfaceScalarField lambda
25 (
26 IOobject
27 (
28 "lambda",
29 mesh.time().timeName(),
30 mesh,
31 IOobject::NO_READ,
32 IOobject::NO_WRITE,
33 false
34 ),
35 mesh,
36 dimless,
37 allLambda,
38 false // Use slices for the couples
39 );
40
41
42 // Create the complete convection flux for alpha1
43 surfaceScalarField alphaPhi1
44 (
45 fvc::flux
46 (
47 phi,
48 alpha1,
49 alphaScheme
50 )
51 + fvc::flux
52 (
53 -fvc::flux(-phir, alpha2, alpharScheme),
54 alpha1,
56 )
57 + fvc::flux
58 (
59 -fvc::flux(-phir, alpha3, alpharScheme),
60 alpha1,
62 )
63 );
64
65 // Create the bounded (upwind) flux for alpha1
66 surfaceScalarField alphaPhi1BD
67 (
68 upwind<scalar>(mesh, phi).flux(alpha1)
69 );
70
71 // Calculate the flux correction for alpha1
72 alphaPhi1 -= alphaPhi1BD;
73
74 // Calculate the limiter for alpha1
75 if (LTS)
76 {
77 const volScalarField& rDeltaT =
78 fv::localEulerDdt::localRDeltaT(mesh);
79
80 MULES::limiter
81 (
82 allLambda,
83 rDeltaT,
84 geometricOneField(),
85 alpha1,
86 alphaPhi1BD,
88 zeroField(),
89 zeroField(),
90 oneField(),
91 zeroField()
92 );
93 }
94 else
95 {
96 MULES::limiter
97 (
98 allLambda,
99 1.0/runTime.deltaT().value(),
100 geometricOneField(),
101 alpha1,
102 alphaPhi1BD,
103 alphaPhi1,
104 zeroField(),
105 zeroField(),
106 oneField(),
107 zeroField()
108 );
109 }
110
111 alphaPhi1 = alphaPhi1BD + lambda*alphaPhi1;
112
113 // Reset allLambda to 1.0
114 allLambda = 1.0;
115
116 // Create the complete flux for alpha2
117 surfaceScalarField alphaPhi2
118 (
119 fvc::flux
120 (
121 phi,
122 alpha2,
123 alphaScheme
124 )
125 + fvc::flux
126 (
127 -fvc::flux(phir, alpha1, alpharScheme),
128 alpha2,
130 )
131 );
132
133 // Create the bounded (upwind) flux for alpha2
134 surfaceScalarField alphaPhi2BD
135 (
136 upwind<scalar>(mesh, phi).flux(alpha2)
137 );
138
139 // Calculate the flux correction for alpha2
140 alphaPhi2 -= alphaPhi2BD;
141
142 // Further limit the limiter for alpha2
143 if (LTS)
144 {
145 const volScalarField& rDeltaT =
146 fv::localEulerDdt::localRDeltaT(mesh);
147
148 MULES::limiter
149 (
150 allLambda,
151 rDeltaT,
152 geometricOneField(),
153 alpha2,
154 alphaPhi2BD,
155 alphaPhi2,
156 zeroField(),
157 zeroField(),
158 oneField(),
159 zeroField()
160 );
161 }
162 else
163 {
164 MULES::limiter
165 (
166 allLambda,
167 1.0/runTime.deltaT().value(),
168 geometricOneField(),
169 alpha2,
170 alphaPhi2BD,
171 alphaPhi2,
172 zeroField(),
173 zeroField(),
174 oneField(),
175 zeroField()
176 );
177 }
178
179 // Construct the limited fluxes
180 alphaPhi2 = alphaPhi2BD + lambda*alphaPhi2;
181
182 // Solve for alpha1
183 solve(fvm::ddt(alpha1) + fvc::div(alphaPhi1));
184
185 // Create the diffusion coefficients for alpha2<->alpha3
186 volScalarField Dc23(D23*max(alpha3, scalar(0))*pos0(alpha2));
187 volScalarField Dc32(D23*max(alpha2, scalar(0))*pos0(alpha3));
188
189 // Add the diffusive flux for alpha3->alpha2
190 alphaPhi2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);
191
192 // Solve for alpha2
193 fvScalarMatrix alpha2Eqn
194 (
195 fvm::ddt(alpha2)
196 + fvc::div(alphaPhi2)
197 - fvm::laplacian(Dc23 + Dc32, alpha2)
198 );
199 alpha2Eqn.solve();
200
201 // Construct the complete mass flux
202 rhoPhi =
203 alphaPhi1*(rho1 - rho3)
204 + (alphaPhi2 + alpha2Eqn.flux())*(rho2 - rho3)
205 + phi*rho3;
206
207 alpha3 = 1.0 - alpha1 - alpha2;
208 }
209
210 Info<< "Air phase volume fraction = "
211 << alpha1.weightedAverage(mesh.V()).value()
212 << " Min(" << alpha1.name() << ") = " << min(alpha1).value()
213 << " Max(" << alpha1.name() << ") = " << max(alpha1).value()
214 << endl;
215
216 Info<< "Liquid phase volume fraction = "
217 << alpha2.weightedAverage(mesh.V()).value()
218 << " Min(" << alpha2.name() << ") = " << min(alpha2).value()
219 << " Max(" << alpha2.name() << ") = " << max(alpha2).value()
220 << endl;
221}
Y[inertIndex] max(0.0)
surfaceScalarField & phi
const volScalarField & alpha1
volScalarField & rho2
const surfaceScalarField & alphaPhi1
volScalarField & rho1
const surfaceScalarField & alphaPhi2
rhoPhi
Definition: alphaEqn.H:6
alpha2
Definition: alphaEqn.H:9
dynamicFvMesh & mesh
engineTime & runTime
bool LTS
Definition: createRDeltaT.H:1
surfaceScalarField phir(fvc::flux(UdmModel.Udm()))
word alpharScheme("div(phirb,alpha)")
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
CEqn solve()
dimensionedScalar lambda("lambda", dimTime/sqr(dimLength), laminarTransport)
Info<< "Creating temperaturePhaseChangeTwoPhaseMixture\n"<< endl;autoPtr< temperaturePhaseChangeTwoPhaseMixture > mixture
Definition: createFields.H:39
label nAlphaCorr(alphaControls.get< label >("nAlphaCorr"))