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  forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
35  {
36  phaseModel& phase = iter();
37 
38  MRF.makeAbsolute(phase.phi().oldTime());
39  MRF.makeAbsolute(phase.phi());
40 
41  HbyAs.set(phasei, new volVectorField(phase.U()));
42  phiHbyAs.set(phasei, new surfaceScalarField(1.0*phase.phi()));
43 
44  phasei++;
45  }
46 
48  (
49  IOobject
50  (
51  "phiHbyA",
52  runTime.timeName(),
53  mesh,
54  IOobject::NO_READ,
55  IOobject::AUTO_WRITE
56  ),
57  mesh,
58  dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0)
59  );
60 
61  volScalarField rho("rho", fluid.rho());
63 
64  phasei = 0;
65  forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
66  {
67  phaseModel& phase = iter();
68  const volScalarField& alpha = phase;
69 
70  alphafs.set(phasei, fvc::interpolate(alpha).ptr());
71  alphafs[phasei].rename("hmm" + alpha.name());
72 
73  volScalarField dragCoeffi
74  (
75  IOobject
76  (
77  "dragCoeffi",
78  runTime.timeName(),
79  mesh
80  ),
81  fluid.dragCoeff(phase, dragCoeffs())/phase.rho(),
82  zeroGradientFvPatchScalarField::typeName
83  );
84  dragCoeffi.correctBoundaryConditions();
85 
86  rAUs.set(phasei, (1.0/(UEqns[phasei].A() + dragCoeffi)).ptr());
87  rAlphaAUfs.set
88  (
89  phasei,
90  (
91  alphafs[phasei]
92  /fvc::interpolate(UEqns[phasei].A() + dragCoeffi)
93  ).ptr()
94  );
95 
97 
98  phiHbyAs[phasei] =
99  (
100  (fvc::interpolate(HbyAs[phasei]) & mesh.Sf())
101  + rAlphaAUfs[phasei]*fvc::ddtCorr(phase.U(), phase.phi())
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())
112  - ghSnGradRho
113  )/phase.rho();
114 
115  multiphaseSystem::dragModelTable::const_iterator dmIter =
116  fluid.dragModels().begin();
117  multiphaseSystem::dragCoeffFields::const_iterator dcIter =
118  dragCoeffs().begin();
119  for
120  (
121  ;
122  dmIter != fluid.dragModels().end() && dcIter != dragCoeffs().end();
123  ++dmIter, ++dcIter
124  )
125  {
126  const phaseModel *phase2Ptr = NULL;
127 
128  if (&phase == &dmIter()->phase1())
129  {
130  phase2Ptr = &dmIter()->phase2();
131  }
132  else if (&phase == &dmIter()->phase2())
133  {
134  phase2Ptr = &dmIter()->phase1();
135  }
136  else
137  {
138  continue;
139  }
140 
141  phiHbyAs[phasei] +=
142  fvc::interpolate((*dcIter())/phase.rho())
143  /fvc::interpolate(UEqns[phasei].A() + dragCoeffi)
144  *phase2Ptr->phi();
145 
146  HbyAs[phasei] +=
147  (1.0/phase.rho())*rAUs[phasei]*(*dcIter())
148  *phase2Ptr->U();
149 
150  // Alternative flux-pressure consistent drag
151  // but not momentum conservative
152  //
153  // HbyAs[phasei] += fvc::reconstruct
154  // (
155  // fvc::interpolate
156  // (
157  // (1.0/phase.rho())*rAUs[phasei]*(*dcIter())
158  // )*phase2Ptr->phi()
159  // );
160  }
161 
163 
164  phasei++;
165  }
166 
168  (
169  IOobject
170  (
171  "rAUf",
172  runTime.timeName(),
173  mesh
174  ),
175  mesh,
176  dimensionedScalar("rAUf", dimensionSet(-1, 3, 1, 0, 0), 0)
177  );
178 
179  phasei = 0;
180  forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
181  {
182  phaseModel& phase = iter();
183  rAUf += mag(alphafs[phasei]*rAlphaAUfs[phasei])/phase.rho();
184 
185  phasei++;
186  }
187 
188  // Update the fixedFluxPressure BCs to ensure flux consistency
189  {
190  surfaceScalarField::GeometricBoundaryField phib(phi.boundaryField());
191  phib = 0;
192  phasei = 0;
193  forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
194  {
195  phaseModel& phase = iter();
196 
197  phib +=
198  alphafs[phasei].boundaryField()
199  *(mesh.Sf().boundaryField() & phase.U().boundaryField());
200 
201  phasei++;
202  }
203 
205  (
206  p_rgh.boundaryField(),
207  (
208  phiHbyA.boundaryField() - MRF.relative(phib)
209  )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
210  );
211  }
212 
213  while (pimple.correctNonOrthogonal())
214  {
215  fvScalarMatrix pEqnIncomp
216  (
218  - fvm::laplacian(rAUf, p_rgh)
219  );
220 
221  pEqnIncomp.setReference(pRefCell, pRefValue);
222 
223  solve
224  (
225  // (
226  // (alpha1/rho1)*pEqnComp1()
227  // + (alpha2/rho2)*pEqnComp2()
228  // ) +
229  pEqnIncomp,
230  mesh.solver(p_rgh.select(pimple.finalInnerIter()))
231  );
232 
233  if (pimple.finalNonOrthogonalIter())
234  {
235  surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
236 
237  phasei = 0;
238  phi = dimensionedScalar("phi", phi.dimensions(), 0);
239  forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
240  {
241  phaseModel& phase = iter();
242 
243  phase.phi() =
245  + rAlphaAUfs[phasei]*mSfGradp/phase.rho();
246  phi +=
248  + mag(alphafs[phasei]*rAlphaAUfs[phasei])
249  *mSfGradp/phase.rho();
250 
251  phasei++;
252  }
253 
254  // dgdt =
255 
256  // (
257  // pos(alpha2)*(pEqnComp2 & p)/rho2
258  // - pos(alpha1)*(pEqnComp1 & p)/rho1
259  // );
260 
261  p_rgh.relax();
262 
263  p = p_rgh + rho*gh;
264 
265  mSfGradp = pEqnIncomp.flux()/rAUf;
266 
267  U = dimensionedVector("U", dimVelocity, vector::zero);
268 
269  phasei = 0;
270  forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
271  {
272  phaseModel& phase = iter();
273  const volScalarField& alpha = phase;
274 
275  phase.U() =
276  HbyAs[phasei]
278  (
279  rAlphaAUfs[phasei]
280  *(
281  (phase.rho() - fvc::interpolate(rho))
282  *(g & mesh.Sf())
283  - ghSnGradRho
284  + mSfGradp
285  )
286  )/phase.rho();
287 
288  //phase.U() = fvc::reconstruct(phase.phi());
289  phase.U().correctBoundaryConditions();
290 
291  U += alpha*phase.U();
292 
293  phasei++;
294  }
295  }
296  }
297 
298  //p = max(p, pMin);
299 
300  #include "continuityErrs.H"
301 
302  // rho1 = rho10 + psi1*p_rgh;
303  // rho2 = rho20 + psi2*p_rgh;
304 
305  // Dp1Dt = fvc::DDt(phi1, p_rgh);
306  // Dp2Dt = fvc::DDt(phi2, p_rgh);
307 }
PtrList< fvVectorMatrix > UEqns(fluid.phases().size())
phaseModel & phase2
Definition: createFields.H:13
PtrList< surfaceScalarField > phiHbyAs(fluid.phases().size())
multiphaseSystem & fluid
Definition: createFields.H:10
phi
Definition: pEqn.H:20
dimensioned< scalar > mag(const dimensioned< Type > &)
label phasei
Definition: pEqn.H:37
setSnGrad< fixedFluxPressureFvPatchScalarField >(p_rgh.boundaryField(),( phiHbyA.boundaryField() -MRF.relative(mesh.Sf().boundaryField()&U.boundaryField()) *rho.boundaryField() )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField()))
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU))
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> reconstruct(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
autoPtr< multiphaseSystem::dragCoeffFields > dragCoeffs(fluid.dragCoeffs())
tmp< GeometricField< typename flux< Type >::type, fvsPatchField, surfaceMesh > > ddtCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
Definition: fvcDdt.C:155
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
volScalarField & p_rgh
tmp< surfaceScalarField > interpolate(const RhoType &rho)
dynamicFvMesh & mesh
const dictionary & pimple
p
Definition: pEqn.H:59
rhoEqn solve()
const scalar pRefValue
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
surfaceScalarField ghSnGradRho(ghf *fvc::snGrad(rho)*mesh.magSf())
IOMRFZoneList & MRF
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
PtrList< volVectorField > HbyAs(fluid.phases().size())
forAllIter(PtrDictionary< phaseModel >, fluid.phases(), iter)
Definition: pEqn.H:34
phiHbyA
Definition: pEqn.H:21
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
PtrList< surfaceScalarField > rAlphaAUfs(fluid.phases().size())
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
const dimensionedVector & g
void correctBoundaryConditions()
Correct boundary field.
const label pRefCell
PtrList< volScalarField > rAUs(fluid.phases().size())
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
rho
Definition: pEqn.H:1
PtrList< surfaceScalarField > alphafs(phases.size())
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
U
Definition: pEqn.H:82
const surfaceScalarField & ghf
const volScalarField & gh
phib
Definition: pEqn.H:191