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(alphaPhic/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 surfaceScalarField::Boundary& phicBf = phic.boundaryFieldRef();
76
77 // Do not compress interface at non-coupled boundary faces
78 // (inlets, outlets etc.)
79 forAll(phic.boundaryField(), patchi)
80 {
81 fvsPatchScalarField& phicp = phicBf[patchi];
82
83 if (!phicp.coupled())
84 {
85 phicp == 0;
86 }
87 }
88
89 tmp<surfaceScalarField> phiCN(alphaPhic);
90
91 // Calculate the Crank-Nicolson off-centred volumetric flux
92 if (ocCoeff > 0)
93 {
94 phiCN = cnCoeff*alphaPhic + (1.0 - cnCoeff)*alphaPhic.oldTime();
95 }
96
98 {
99 #include "alphaSuSp.H"
100
101 fvScalarMatrix alpha1Eqn
102 (
103 (
104 LTS
105 ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alphac, alpha1)
106 : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
107 )
108 + fv::gaussConvectionScheme<scalar>
109 (
110 mesh,
111 phiCN,
112 upwind<scalar>(mesh, phiCN)
113 ).fvmDiv(phiCN, alpha1)
114 - fvm::Sp(fvc::ddt(alphac) + fvc::div(phiCN), alpha1)
115 ==
116 Su + fvm::Sp(Sp + divU, alpha1)
117 );
118
119 alpha1Eqn.solve();
120
121 Info<< "Phase-1 volume fraction = "
122 << alpha1.weightedAverage(mesh.Vsc()).value()
123 << " Min(" << alpha1.name() << ") = " << min(alpha1).value()
124 << " Max(" << alpha1.name() << ") = " << max(alpha1).value()
125 << endl;
126
127 tmp<surfaceScalarField> talphaPhi1UD(alpha1Eqn.flux());
128 alphaPhi10 = talphaPhi1UD();
129
130 if (alphaApplyPrevCorr && talphaPhi1Corr0.valid())
131 {
132 Info<< "Applying the previous iteration compression flux" << endl;
133 MULES::correct
134 (
135 alphac,
136 alpha1,
138 talphaPhi1Corr0.ref(),
139 zeroField(), zeroField(),
140 oneField(),
141 zeroField()
142 );
143
145 }
146
147 // Cache the upwind-flux
148 talphaPhi1Corr0 = talphaPhi1UD;
149
150 alpha2 = 1.0 - alpha1;
151
152 mixture.correct();
153 }
154
155
156 for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
157 {
158 #include "alphaSuSp.H"
159
160 surfaceScalarField phir(phic*mixture.nHatf());
161
162 tmp<surfaceScalarField> talphaPhi1Un
163 (
164 fvc::flux
165 (
166 phiCN(),
167 cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(),
168 alphaScheme
169 )
170 + fvc::flux
171 (
172 -fvc::flux(-phir, alpha2, alpharScheme),
173 alpha1,
175 )
176 );
177
178 if (MULESCorr)
179 {
180 tmp<surfaceScalarField> talphaPhi1Corr(talphaPhi1Un() - alphaPhi10);
181 volScalarField alpha10("alpha10", alpha1);
182
183 MULES::correct
184 (
185 alphac,
186 alpha1,
187 talphaPhi1Un(),
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 {
208 alphaPhi10 = talphaPhi1Un;
209
210 MULES::explicitSolve
211 (
212 alphac,
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 if
239 (
240 word(mesh.ddtScheme("ddt(rho,U)"))
241 == fv::EulerDdtScheme<vector>::typeName
242 )
243 {
244 #include "rhofs.H"
246 }
247 else
248 {
249 if (ocCoeff > 0)
250 {
251 // Calculate the end-of-time-step alpha flux
252 alphaPhi10 =
253 (alphaPhi10 - (1.0 - cnCoeff)*alphaPhi10.oldTime())/cnCoeff;
254 }
255
256 // Calculate the end-of-time-step mass flux
257 #include "rhofs.H"
258 rhoPhi = alphaPhi10*(rho1f - rho2f) + alphaPhic*rho2f;
259 }
260
261 Info<< "Phase-1 volume fraction = "
262 << alpha1.weightedAverage(mesh.Vsc()).value()
263 << " Min(" << alpha1.name() << ") = " << min(alpha1).value()
264 << " Max(" << alpha1.name() << ") = " << max(alpha1).value()
265 << endl;
266}
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)
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))
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