alphaEqn.H
Go to the documentation of this file.
1 {
2  #include "alphaScheme.H"
3 
4  // Set the off-centering coefficient according to ddt scheme
5  scalar ocCoeff = 0;
6  {
7  tmp<fv::ddtScheme<scalar>> tddtAlpha
8  (
10  (
11  mesh,
12  mesh.ddtScheme("ddt(alpha)")
13  )
14  );
15  const fv::ddtScheme<scalar>& ddtAlpha = tddtAlpha();
16 
17  if
18  (
19  isType<fv::EulerDdtScheme<scalar>>(ddtAlpha)
20  || isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha)
21  )
22  {
23  ocCoeff = 0;
24  }
25  else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha))
26  {
27  if (nAlphaSubCycles > 1)
28  {
30  << "Sub-cycling is not supported "
31  "with the CrankNicolson ddt scheme"
32  << exit(FatalError);
33  }
34 
35  if
36  (
38  || mesh.time().timeIndex() > mesh.time().startTimeIndex() + 1
39  )
40  {
41  ocCoeff =
42  refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha)
43  .ocCoeff();
44  }
45  }
46  else
47  {
49  << "Only Euler and CrankNicolson ddt schemes are supported"
50  << exit(FatalError);
51  }
52  }
53 
54  // Set the time blending factor, 1 for Euler
55  scalar cnCoeff = 1.0/(1.0 + ocCoeff);
56 
57  tmp<surfaceScalarField> phiCN(phi);
58 
59  // Calculate the Crank-Nicolson off-centred volumetric flux
60  if (ocCoeff > 0)
61  {
63  (
64  "phiCN",
65  cnCoeff*phi + (1.0 - cnCoeff)*phi.oldTime()
66  );
67  }
68 
69  if (MULESCorr)
70  {
71  #include "alphaSuSp.H"
72 
73  fvScalarMatrix alpha1Eqn
74  (
75  (
76  LTS
77  ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
78  : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
79  )
80  + fv::gaussConvectionScheme<scalar>
81  (
82  mesh,
83  phiCN,
84  upwind<scalar>(mesh, phiCN)
85  ).fvmDiv(phiCN, alpha1)
86  // - fvm::Sp(fvc::ddt(dimensionedScalar(dimless, 1), mesh)
87  // + fvc::div(phiCN), alpha1)
88  ==
89  Su + fvm::Sp(Sp + divU, alpha1)
90  );
91 
92  alpha1Eqn.solve();
93 
94  Info<< "Phase-1 volume fraction = "
95  << alpha1.weightedAverage(mesh.Vsc()).value()
96  << " Min(" << alpha1.name() << ") = " << min(alpha1).value()
97  << " Max(" << alpha1.name() << ") = " << max(alpha1).value()
98  << endl;
99 
100  tmp<surfaceScalarField> talphaPhi1UD(alpha1Eqn.flux());
101  alphaPhi10 = talphaPhi1UD();
102 
103  if (alphaApplyPrevCorr && talphaPhi1Corr0.valid())
104  {
105  Info<< "Applying the previous iteration compression flux" << endl;
107  (
108  geometricOneField(),
109  alpha1,
110  alphaPhi10,
111  talphaPhi1Corr0.ref(),
112  oneField(),
113  zeroField()
114  );
115 
116  alphaPhi10 += talphaPhi1Corr0();
117  }
118 
119  // Cache the upwind-flux
120  talphaPhi1Corr0 = talphaPhi1UD;
121 
122  alpha2 = 1.0 - alpha1;
123 
124  mixture.correct();
125  }
126 
127  for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
128  {
129  #include "alphaSuSp.H"
130 
131  // Split operator
132  tmp<surfaceScalarField> talphaPhi1Un
133  (
134  fvc::flux
135  (
136  phiCN(),
137  (cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime())(),
138  compressionScheme.rewind()
139  )
140  );
141 
142  if (MULESCorr)
143  {
144  tmp<surfaceScalarField> talphaPhi1Corr(talphaPhi1Un() - alphaPhi10);
145  volScalarField alpha10("alpha10", alpha1);
146 
148  (
149  geometricOneField(),
150  alpha1,
151  talphaPhi1Un(),
152  talphaPhi1Corr.ref(),
153  Sp,
154  (-Sp*alpha1)(),
155  oneField(),
156  zeroField()
157  );
158 
159  // Under-relax the correction for all but the 1st corrector
160  if (aCorr == 0)
161  {
162  alphaPhi10 += talphaPhi1Corr();
163  }
164  else
165  {
166  alpha1 = 0.5*alpha1 + 0.5*alpha10;
167  alphaPhi10 += 0.5*talphaPhi1Corr();
168  }
169  }
170  else
171  {
172  alphaPhi10 = talphaPhi1Un;
173 
175  (
176  geometricOneField(),
177  alpha1,
178  phiCN,
179  alphaPhi10,
180  Sp,
181  (Su + divU*min(alpha1(), scalar(1)))(),
182  oneField(),
183  zeroField()
184  );
185  }
186 
187  alpha2 = 1.0 - alpha1;
188 
189  mixture.correct();
190  }
191 
193  {
195  talphaPhi1Corr0.ref().rename("alphaPhi1Corr0");
196  }
197  else
198  {
199  talphaPhi1Corr0.clear();
200  }
201 
202  #include "rhofs.H"
203 
204  if
205  (
206  word(mesh.ddtScheme("ddt(rho,U)"))
207  == fv::EulerDdtScheme<vector>::typeName
208  || word(mesh.ddtScheme("ddt(rho,U)"))
209  == fv::localEulerDdtScheme<vector>::typeName
210  )
211  {
213  }
214  else
215  {
216  if (ocCoeff > 0)
217  {
218  // Calculate the end-of-time-step alpha flux
219  alphaPhi10 =
220  (alphaPhi10 - (1.0 - cnCoeff)*alphaPhi10.oldTime())/cnCoeff;
221  }
222 
223  // Calculate the end-of-time-step mass flux
225  }
226 
227  Info<< "Phase-1 volume fraction = "
228  << alpha1.weightedAverage(mesh.Vsc()).value()
229  << " Min(" << alpha1.name() << ") = " << min(alpha1).value()
230  << " Max(" << alpha1.name() << ") = " << max(alpha1).value()
231  << endl;
232 }
ITstream compressionScheme(compressionSchemes.found(alphaScheme) ? mesh.divScheme(divAlphaName) :ITstream(divAlphaName, tokenList { word("Gauss"), word("interfaceCompression"), alphaScheme, alphaControls.lookup< scalar >("cAlpha") }))
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
volScalarField alpha10("alpha10", alpha1)
zeroField Su
Definition: alphaSuSp.H:1
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
surfaceScalarField rho1f(fvc::interpolate(rho1))
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const bool alphaApplyPrevCorr(alphaControls.lookupOrDefault< Switch >("alphaApplyPrevCorr", false))
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const bool alphaRestart
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
surfaceScalarField alphaPhi10(alphaPhi10Header, phi *fvc::interpolate(alpha1))
alpha2
Definition: alphaEqn.H:115
phi
Definition: pEqn.H:104
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
surfaceScalarField rho2f(fvc::interpolate(rho2))
dynamicFvMesh & mesh
rhoPhi
Definition: alphaEqn.H:106
autoPtr< BasicCompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const typename BasicCompressibleMomentumTransportModel::transportModel &transport)
const fv::ddtScheme< scalar > & ddtAlpha
Definition: alphaEqn.H:15
volScalarField & alpha1(mixture.alpha1())
bool isType(const Type &t)
Check the typeid.
Definition: typeInfo.H:126
tmp< surfaceScalarField > phiCN(phi)
bool LTS
Definition: createRDeltaT.H:1
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(rho0+psi *p, rhoMin);# 1 "/home/ubuntu/OpenFOAM-8/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
Definition: pEqn.H:68
scalar ocCoeff
Definition: alphaEqn.H:5
const label nAlphaSubCycles(alphaControls.lookup< label >("nAlphaSubCycles"))
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
messageStream Info
const label nAlphaCorr(alphaControls.lookup< label >("nAlphaCorr"))
scalar cnCoeff
Definition: alphaEqn.H:55
zeroField divU
Definition: alphaSuSp.H:3
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
Definition: fvcFlux.C:32
phaseChangeTwoPhaseMixture & mixture
Definition: createFields.H:38
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
const bool MULESCorr(alphaControls.lookupOrDefault< Switch >("MULESCorr", false))
zeroField Sp
Definition: alphaSuSp.H:2
tmp< surfaceScalarField > talphaPhi1Corr0