pEqn.H
Go to the documentation of this file.
1{
2 // rho1 = rho10 + psi1*p_rgh;
3 // rho2 = rho20 + psi2*p_rgh;
4
5 // tmp<fvScalarMatrix> pEqnComp1;
6 // tmp<fvScalarMatrix> pEqnComp2;
7
8 // //if (transonic)
9 // //{
10 // //}
11 // //else
12 // {
13 // surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi1);
14 // surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi2);
15
16 // pEqnComp1 =
17 // fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
18 // + fvc::div(phid1, p_rgh)
19 // - fvc::Sp(fvc::div(phid1), p_rgh);
20
21 // pEqnComp2 =
22 // fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
23 // + fvc::div(phid2, p_rgh)
24 // - fvc::Sp(fvc::div(phid2), p_rgh);
25 // }
26
27 PtrList<surfaceScalarField> alphafs(fluid.phases().size());
28 PtrList<volVectorField> HbyAs(fluid.phases().size());
29 PtrList<surfaceScalarField> phiHbyAs(fluid.phases().size());
30 PtrList<volScalarField> rAUs(fluid.phases().size());
31 PtrList<surfaceScalarField> rAlphaAUfs(fluid.phases().size());
32
33 phasei = 0;
34 for (phaseModel& phase : fluid.phases())
35 {
36 MRF.makeAbsolute(phase.phi().oldTime());
37 MRF.makeAbsolute(phase.phi());
38
39 HbyAs.set(phasei, new volVectorField(phase.U()));
40 phiHbyAs.set(phasei, new surfaceScalarField(1.0*phase.phi()));
41
42 ++phasei;
43 }
44
45 surfaceScalarField phiHbyA
46 (
47 IOobject
48 (
49 "phiHbyA",
50 runTime.timeName(),
51 mesh,
52 IOobject::NO_READ,
53 IOobject::AUTO_WRITE
54 ),
55 mesh,
56 dimensionedScalar(dimArea*dimVelocity, Zero)
57 );
58
59 volScalarField rho("rho", fluid.rho());
60 surfaceScalarField ghSnGradRho(ghf*fvc::snGrad(rho)*mesh.magSf());
61
62 phasei = 0;
63 for (phaseModel& phase : fluid.phases())
64 {
65 const volScalarField& alpha = phase;
66
67 alphafs.set(phasei, fvc::interpolate(alpha).ptr());
68 alphafs[phasei].rename("hmm" + alpha.name());
69
70 volScalarField dragCoeffi
71 (
72 IOobject
73 (
74 "dragCoeffi",
75 runTime.timeName(),
76 mesh
77 ),
78 fluid.dragCoeff(phase, dragCoeffs())/phase.rho(),
80 );
81 dragCoeffi.correctBoundaryConditions();
82
83 rAUs.set(phasei, (1.0/(UEqns[phasei].A() + dragCoeffi)).ptr());
84 rAlphaAUfs.set
85 (
86 phasei,
87 (
89 /fvc::interpolate(UEqns[phasei].A() + dragCoeffi)
90 ).ptr()
91 );
92
94
96 (
97 fvc::flux(HbyAs[phasei])
98 + MRF.zeroFilter
99 (
100 rAlphaAUfs[phasei]*fvc::ddtCorr(phase.U(), phase.phi())
101 )
102 );
103 MRF.makeRelative(phiHbyAs[phasei]);
104 MRF.makeRelative(phase.phi().oldTime());
105 MRF.makeRelative(phase.phi());
106
107 phiHbyAs[phasei] +=
109 *(
110 fluid.surfaceTension(phase)*mesh.magSf()
111 + (phase.rho() - fvc::interpolate(rho))*(g & mesh.Sf())
113 )/phase.rho();
114
115 auto dmIter = fluid.dragModels().cbegin();
116 auto dcIter = dragCoeffs().cbegin();
117
118 for
119 (
120 ;
121 dmIter.good() && dcIter.good();
122 ++dmIter, ++dcIter
123 )
124 {
125 const phaseModel *phase2Ptr = nullptr;
126
127 if (&phase == &dmIter()->phase1())
128 {
129 phase2Ptr = &dmIter()->phase2();
130 }
131 else if (&phase == &dmIter()->phase2())
132 {
133 phase2Ptr = &dmIter()->phase1();
134 }
135 else
136 {
137 continue;
138 }
139
140 phiHbyAs[phasei] +=
141 fvc::interpolate((*dcIter())/phase.rho())
142 /fvc::interpolate(UEqns[phasei].A() + dragCoeffi)
143 *phase2Ptr->phi();
144
145 HbyAs[phasei] +=
146 (1.0/phase.rho())*rAUs[phasei]*(*dcIter())
147 *phase2Ptr->U();
148
149 // Alternative flux-pressure consistent drag
150 // but not momentum conservative
151 //
152 // HbyAs[phasei] += fvc::reconstruct
153 // (
154 // fvc::interpolate
155 // (
156 // (1.0/phase.rho())*rAUs[phasei]*(*dcIter())
157 // )*phase2Ptr->phi()
158 // );
159 }
160
162
163 ++phasei;
164 }
165
166 surfaceScalarField rAUf
167 (
168 IOobject
169 (
170 "rAUf",
171 runTime.timeName(),
172 mesh
173 ),
174 mesh,
175 dimensionedScalar(dimensionSet(-1, 3, 1, 0, 0), Zero)
176 );
177
178 phasei = 0;
179 for (const phaseModel& phase : fluid.phases())
180 {
181 rAUf += mag(alphafs[phasei]*rAlphaAUfs[phasei])/phase.rho();
182
183 ++phasei;
184 }
185
186 // Update the fixedFluxPressure BCs to ensure flux consistency
187 {
188 surfaceScalarField::Boundary phib(phi.boundaryField());
189 phib = 0;
190 phasei = 0;
191 for (const phaseModel& phase : fluid.phases())
192 {
193 phib +=
194 alphafs[phasei].boundaryField()
195 *(mesh.Sf().boundaryField() & phase.U().boundaryField());
196
197 ++phasei;
198 }
199
201 (
202 p_rgh.boundaryFieldRef(),
203 (
204 phiHbyA.boundaryField() - MRF.relative(phib)
205 )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
206 );
207 }
208
209 while (pimple.correctNonOrthogonal())
210 {
211 fvScalarMatrix pEqnIncomp
212 (
213 fvc::div(phiHbyA)
214 - fvm::laplacian(rAUf, p_rgh)
215 );
216
217 pEqnIncomp.setReference(pRefCell, pRefValue);
218
219 solve
220 (
221 // (
222 // (alpha1/rho1)*pEqnComp1()
223 // + (alpha2/rho2)*pEqnComp2()
224 // ) +
225 pEqnIncomp,
226 mesh.solver(p_rgh.select(pimple.finalInnerIter()))
227 );
228
229 if (pimple.finalNonOrthogonalIter())
230 {
231 surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
232
233 phasei = 0;
234 phi = dimensionedScalar("phi", phi.dimensions(), Zero);
235
236 for (phaseModel& phase : fluid.phases())
237 {
238 phase.phi() =
240 + rAlphaAUfs[phasei]*mSfGradp/phase.rho();
241
242 phi +=
245 *mSfGradp/phase.rho();
246
247 ++phasei;
248 }
249
250 // dgdt =
251
252 // (
253 // pos0(alpha2)*(pEqnComp2 & p)/rho2
254 // - pos0(alpha1)*(pEqnComp1 & p)/rho1
255 // );
256
257 p_rgh.relax();
258
259 p = p_rgh + rho*gh;
260
261 mSfGradp = pEqnIncomp.flux()/rAUf;
262
263 U = dimensionedVector("U", dimVelocity, Zero);
264
265 phasei = 0;
266 for (phaseModel& phase : fluid.phases())
267 {
268 const volScalarField& alpha = phase;
269
270 phase.U() =
272 + fvc::reconstruct
273 (
275 *(
276 (phase.rho() - fvc::interpolate(rho))
277 *(g & mesh.Sf())
279 + mSfGradp
280 )
281 )/phase.rho();
282
283 //phase.U() = fvc::reconstruct(phase.phi());
284 phase.U().correctBoundaryConditions();
285
286 U += alpha*phase.U();
287
288 ++phasei;
289 }
290 }
291 }
292
293 //p = max(p, pMin);
294
295 #include "continuityErrs.H"
296
297 // rho1 = rho10 + psi1*p_rgh;
298 // rho2 = rho20 + psi2*p_rgh;
299
300 // Dp1Dt = fvc::DDt(phi1, p_rgh);
301 // Dp2Dt = fvc::DDt(phi2, p_rgh);
302}
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
const uniformDimensionedVectorField & g
volScalarField & p_rgh
surfaceScalarField & phi
const scalar pRefValue
const surfaceScalarField & ghf
const label pRefCell
IOMRFZoneList & MRF
const volScalarField & gh
pimpleControl & pimple
twoPhaseSystem & fluid
phaseModel & phase1
phaseModel & phase2
U
Definition: pEqn.H:72
volScalarField & p
phiHbyA
Definition: pcEqn.H:73
dynamicFvMesh & mesh
engineTime & runTime
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU))
autoPtr< multiphaseSystem::dragCoeffFields > dragCoeffs(fluid.dragCoeffs())
PtrList< fvVectorMatrix > UEqns(fluid.phases().size())
label phasei
Definition: pEqn.H:27
phib
Definition: pEqn.H:189
surfaceScalarField ghSnGradRho(ghf *fvc::snGrad(rho) *mesh.magSf())
PtrList< surfaceScalarField > phiHbyAs(fluid.phases().size())
PtrList< surfaceScalarField > rAlphaAUfs(fluid.phases().size())
PtrList< volVectorField > HbyAs(fluid.phases().size())
setSnGrad< fixedFluxPressureFvPatchScalarField >(p_rgh.boundaryFieldRef(),(phiHbyA.boundaryField() - MRF.relative(phib))/(mesh.magSf().boundaryField() *rAUf.boundaryField()))
PtrList< surfaceScalarField > alphafs(phases.size())
PtrList< volScalarField > rAUs
Definition: pEqn.H:4
volScalarField & alpha
CEqn solve()
static const char *const typeName
The type name used in ensight case files.