alphaEqn.H
Go to the documentation of this file.
1{
2 word alphaScheme("div(phi,alpha)");
3 word alpharScheme("div(phirb,alpha)");
4
5 // Set the off-centering coefficient according to ddt scheme
6 scalar ocCoeff = 0;
7 {
8 tmp<fv::ddtScheme<scalar>> tddtAlpha
9 (
10 fv::ddtScheme<scalar>::New
11 (
12 mesh,
13 mesh.ddtScheme("ddt(alpha)")
14 )
15 );
16 const fv::ddtScheme<scalar>& ddtAlpha = tddtAlpha();
17
18 if
19 (
20 isType<fv::EulerDdtScheme<scalar>>(ddtAlpha)
21 || isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha)
22 )
23 {
24 ocCoeff = 0;
25 }
26 else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha))
27 {
28 if (nAlphaSubCycles > 1)
29 {
31 << "Sub-cycling is not supported "
32 "with the CrankNicolson ddt scheme"
33 << exit(FatalError);
34 }
35
36 if
37 (
39 || mesh.time().timeIndex() > mesh.time().startTimeIndex() + 1
40 )
41 {
42 ocCoeff =
43 refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha)
44 .ocCoeff();
45 }
46 }
47 else
48 {
50 << "Only Euler and CrankNicolson ddt schemes are supported"
51 << exit(FatalError);
52 }
53 }
54
55 // Set the time blending factor, 1 for Euler
56 scalar cnCoeff = 1.0/(1.0 + ocCoeff);
57
58 // Standard face-flux compression coefficient
59 surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()));
60
61 // Add the optional isotropic compression contribution
62 if (icAlpha > 0)
63 {
64 phic *= (1.0 - icAlpha);
65 phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
66 }
67
68 // Add the optional shear compression contribution
69 if (scAlpha > 0)
70 {
71 phic +=
72 scAlpha*mag(mesh.delta() & fvc::interpolate(symm(fvc::grad(U))));
73 }
74
75
76 surfaceScalarField::Boundary& phicBf =
77 phic.boundaryFieldRef();
78
79 // Do not compress interface at non-coupled boundary faces
80 // (inlets, outlets etc.)
81 forAll(phic.boundaryField(), patchi)
82 {
83 fvsPatchScalarField& phicp = phicBf[patchi];
84
85 if (!phicp.coupled())
86 {
87 phicp == 0;
88 }
89 }
90
91 tmp<surfaceScalarField> phiCN(phi);
92
93 // Calculate the Crank-Nicolson off-centred volumetric flux
94 if (ocCoeff > 0)
95 {
96 phiCN = cnCoeff*phi + (1.0 - cnCoeff)*phi.oldTime();
97 }
98
100 {
101 #include "alphaSuSp.H"
102
103 fvScalarMatrix alpha1Eqn
104 (
105 (
106 LTS
107 ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
108 : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
109 )
110 + fv::gaussConvectionScheme<scalar>
111 (
112 mesh,
113 phiCN,
114 upwind<scalar>(mesh, phiCN)
115 ).fvmDiv(phiCN, alpha1)
116 // - fvm::Sp(fvc::ddt(dimensionedScalar("1", dimless, 1), mesh)
117 // + fvc::div(phiCN), alpha1)
118 ==
119 Su + fvm::Sp(Sp + divU, alpha1)
120 );
121
122 alpha1Eqn.solve();
123
124 Info<< "Phase-1 volume fraction = "
125 << alpha1.weightedAverage(mesh.Vsc()).value()
126 << " Min(" << alpha1.name() << ") = " << min(alpha1).value()
127 << " Max(" << alpha1.name() << ") = " << max(alpha1).value()
128 << endl;
129
130 tmp<surfaceScalarField> talphaPhi1UD(alpha1Eqn.flux());
131 alphaPhi10 = talphaPhi1UD();
132
133 if (alphaApplyPrevCorr && talphaPhi1Corr0.valid())
134 {
135 Info<< "Applying the previous iteration compression flux" << endl;
136 MULES::correct
137 (
138 geometricOneField(),
139 alpha1,
141 talphaPhi1Corr0.ref(),
142 oneField(),
143 zeroField()
144 );
145
147 }
148
149 // Cache the upwind-flux
150 talphaPhi1Corr0 = talphaPhi1UD;
151
152 alpha2 = 1.0 - alpha1;
153
154 mixture.correct();
155 }
156
157
158 for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
159 {
160 #include "alphaSuSp.H"
161
162 surfaceScalarField phir(phic*mixture.nHatf());
163
164 alphaPhiUn =
165 fvc::flux
166 (
167 phi,
168 alpha1,
169 alphaScheme
170 )
171 + fvc::flux
172 (
173 -fvc::flux(-phir, alpha2, alpharScheme),
174 alpha1,
176 );
177
178 if (MULESCorr)
179 {
180 tmp<surfaceScalarField> talphaPhi1Corr(alphaPhiUn - alphaPhi10);
181 volScalarField alpha10("alpha10", alpha1);
182
183 MULES::correct
184 (
185 geometricOneField(),
186 alpha1,
188 talphaPhi1Corr.ref(),
189 Sp,
190 (-Sp*alpha1)(),
191 oneField(),
192 zeroField()
193 );
194
195 // Under-relax the correction for all but the 1st corrector
196 if (aCorr == 0)
197 {
198 alphaPhi10 += talphaPhi1Corr();
199 }
200 else
201 {
202 alpha1 = 0.5*alpha1 + 0.5*alpha10;
203 alphaPhi10 += 0.5*talphaPhi1Corr();
204 }
205 }
206 else
207 {
209
210 MULES::explicitSolve
211 (
212 geometricOneField(),
213 alpha1,
214 phiCN,
216 Sp,
217 (Su + divU*min(alpha1(), scalar(1)))(),
218 oneField(),
219 zeroField()
220 );
221 }
222
223 alpha2 = 1.0 - alpha1;
224
225 mixture.correct();
226 }
227
229 {
231 talphaPhi1Corr0.ref().rename("alphaPhi1Corr0");
232 }
233 else
234 {
235 talphaPhi1Corr0.clear();
236 }
237
238 #include "rhofs.H"
239
240 if
241 (
242 word(mesh.ddtScheme("ddt(rho,U)"))
243 == fv::EulerDdtScheme<vector>::typeName
244 || word(mesh.ddtScheme("ddt(rho,U)"))
245 == fv::localEulerDdtScheme<vector>::typeName
246 )
247 {
249 }
250 else
251 {
252 if (ocCoeff > 0)
253 {
254 // Calculate the end-of-time-step alpha flux
255 alphaPhi10 =
256 (alphaPhi10 - (1.0 - cnCoeff)*alphaPhi10.oldTime())/cnCoeff;
257 }
258
259 // Calculate the end-of-time-step mass flux
261 }
262
263 Info<< "Phase-1 volume fraction = "
264 << alpha1.weightedAverage(mesh.Vsc()).value()
265 << " Min(" << alpha1.name() << ") = " << min(alpha1).value()
266 << " Max(" << alpha1.name() << ") = " << max(alpha1).value()
267 << endl;
268}
scalar cnCoeff
Definition: alphaEqn.H:56
tmp< surfaceScalarField > phiCN(alphaPhic)
surfaceScalarField::Boundary & phicBf
Definition: alphaEqn.H:75
surfaceScalarField phic(mixture.cAlpha() *mag(alphaPhic/mesh.magSf()))
scalar ocCoeff
Definition: alphaEqn.H:6
const fv::ddtScheme< scalar > & ddtAlpha
Definition: alphaEqn.H:16
Y[inertIndex] max(0.0)
surfaceScalarField & phi
const volScalarField & alpha1
U
Definition: pEqn.H:72
alphaPhi10
Definition: alphaEqn.H:7
rhoPhi
Definition: alphaEqn.H:6
alpha2
Definition: alphaEqn.H:9
surfaceScalarField rho2f(fvc::interpolate(rho2))
surfaceScalarField rho1f(fvc::interpolate(rho1))
surfaceScalarField alphaPhiUn(IOobject("alphaPhiUn", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE), mesh, dimensionedScalar(phi.dimensions(), Zero))
const bool alphaRestart
tmp< surfaceScalarField > talphaPhi1Corr0
dynamicFvMesh & mesh
bool LTS
Definition: createRDeltaT.H:1
surfaceScalarField phir(fvc::flux(UdmModel.Udm()))
word alpharScheme("div(phirb,alpha)")
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
zeroField divU
Definition: alphaSuSp.H:3
zeroField Su
Definition: alphaSuSp.H:1
zeroField Sp
Definition: alphaSuSp.H:2
volScalarField alpha10("alpha10", alpha1)
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
Info<< "Creating temperaturePhaseChangeTwoPhaseMixture\n"<< endl;autoPtr< temperaturePhaseChangeTwoPhaseMixture > mixture
Definition: createFields.H:39
bool MULESCorr(alphaControls.getOrDefault("MULESCorr", false))
label nAlphaSubCycles(alphaControls.get< label >("nAlphaSubCycles"))
scalar scAlpha(alphaControls.getOrDefault< scalar >("scAlpha", 0))
label nAlphaCorr(alphaControls.get< label >("nAlphaCorr"))
scalar icAlpha(alphaControls.getOrDefault< scalar >("icAlpha", 0))
bool alphaApplyPrevCorr(alphaControls.getOrDefault("alphaApplyPrevCorr", false))
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333