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  #include "alphaSuSp.H"
70 
71  if (MULESCorr)
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  );
87 
88  if (divU.valid())
89  {
90  alpha1Eqn -= Su() + fvm::Sp(Sp() + divU(), alpha1);
91  }
92 
93  alpha1Eqn.solve();
94 
95  Info<< "Phase-1 volume fraction = "
96  << alpha1.weightedAverage(mesh.Vsc()).value()
97  << " Min(" << alpha1.name() << ") = " << min(alpha1).value()
98  << " Max(" << alpha1.name() << ") = " << max(alpha1).value()
99  << endl;
100 
101  tmp<surfaceScalarField> talphaPhi1UD(alpha1Eqn.flux());
102  alphaPhi10 = talphaPhi1UD();
103 
104  if (alphaApplyPrevCorr && talphaPhi1Corr0.valid())
105  {
106  Info<< "Applying the previous iteration compression flux" << endl;
108  (
109  geometricOneField(),
110  alpha1,
111  alphaPhi10,
112  talphaPhi1Corr0.ref(),
113  oneField(),
114  zeroField()
115  );
116 
117  alphaPhi10 += talphaPhi1Corr0();
118  }
119 
120  // Cache the upwind-flux
121  talphaPhi1Corr0 = talphaPhi1UD;
122 
123  alpha2 = 1.0 - alpha1;
124 
125  mixture.correct();
126  }
127 
128  for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
129  {
130  // Split operator
131  tmp<surfaceScalarField> talphaPhi1Un
132  (
133  fvc::flux
134  (
135  phiCN(),
136  (cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime())(),
137  compressionScheme.rewind()
138  )
139  );
140 
141  if (MULESCorr)
142  {
143  tmp<surfaceScalarField> talphaPhi1Corr(talphaPhi1Un() - alphaPhi10);
144  volScalarField alpha10("alpha10", alpha1);
145 
146  if (divU.valid())
147  {
149  (
150  geometricOneField(),
151  alpha1,
152  talphaPhi1Un(),
153  talphaPhi1Corr.ref(),
154  Sp(),
155  (-Sp()*alpha1)(),
156  oneField(),
157  zeroField()
158  );
159  }
160  else
161  {
163  (
164  geometricOneField(),
165  alpha1,
166  talphaPhi1Un(),
167  talphaPhi1Corr.ref(),
168  oneField(),
169  zeroField()
170  );
171  }
172 
173  // Under-relax the correction for all but the 1st corrector
174  if (aCorr == 0)
175  {
176  alphaPhi10 += talphaPhi1Corr();
177  }
178  else
179  {
180  alpha1 = 0.5*alpha1 + 0.5*alpha10;
181  alphaPhi10 += 0.5*talphaPhi1Corr();
182  }
183  }
184  else
185  {
186  alphaPhi10 = talphaPhi1Un;
187 
188  if (divU.valid())
189  {
191  (
192  geometricOneField(),
193  alpha1,
194  phiCN,
195  alphaPhi10,
196  Sp(),
197  (Su() + divU()*min(alpha1(), scalar(1)))(),
198  oneField(),
199  zeroField()
200  );
201  }
202  else
203  {
205  (
206  geometricOneField(),
207  alpha1,
208  phiCN,
209  alphaPhi10,
210  oneField(),
211  zeroField()
212  );
213  }
214  }
215 
216  alpha2 = 1.0 - alpha1;
217 
218  mixture.correct();
219  }
220 
222  {
224  talphaPhi1Corr0.ref().rename("alphaPhi1Corr0");
225  }
226  else
227  {
228  talphaPhi1Corr0.clear();
229  }
230 
231  #include "rhofs.H"
232 
233  if
234  (
235  word(mesh.ddtScheme("ddt(rho,U)"))
236  == fv::EulerDdtScheme<vector>::typeName
237  || word(mesh.ddtScheme("ddt(rho,U)"))
238  == fv::localEulerDdtScheme<vector>::typeName
239  )
240  {
242  }
243  else
244  {
245  if (ocCoeff > 0)
246  {
247  // Calculate the end-of-time-step alpha flux
248  alphaPhi10 =
249  (alphaPhi10 - (1.0 - cnCoeff)*alphaPhi10.oldTime())/cnCoeff;
250  }
251 
252  // Calculate the end-of-time-step mass flux
254  }
255 
256  Info<< "Phase-1 volume fraction = "
257  << alpha1.weightedAverage(mesh.Vsc()).value()
258  << " Min(" << alpha1.name() << ") = " << min(alpha1).value()
259  << " Max(" << alpha1.name() << ") = " << max(alpha1).value()
260  << endl;
261 }
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
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:323
const bool alphaRestart
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
volScalarField & alpha1(mixture.alpha1())
surfaceScalarField alphaPhi10(alphaPhi10Header, phi *fvc::interpolate(alpha1))
rhoPhi
Definition: alphaEqn.H:25
alpha2
Definition: alphaEqn.H:115
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
surfaceScalarField rho2f(fvc::interpolate(rho2))
dynamicFvMesh & mesh
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
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(rho0+psi *p, rhoMin);# 1 "/home/ubuntu/OpenFOAM-9/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
bool isType(const Type &t)
Check the typeid.
Definition: typeInfo.H:126
phi
Definition: correctPhi.H:3
tmp< surfaceScalarField > phiCN(phi)
bool LTS
Definition: createRDeltaT.H:1
scalar ocCoeff
Definition: alphaEqn.H:5
const label nAlphaSubCycles(alphaControls.lookup< label >("nAlphaSubCycles"))
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
tmp< volScalarField > divU
Definition: alphaSuSp.H:9
const tmp< volScalarField::Internal > & Sp
Definition: alphaSuSp.H:7
messageStream Info
const label nAlphaCorr(alphaControls.lookup< label >("nAlphaCorr"))
scalar cnCoeff
Definition: alphaEqn.H:55
const tmp< volScalarField::Internal > & Su
Definition: alphaSuSp.H:6
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
Definition: fvcFlux.C:32
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))
tmp< surfaceScalarField > talphaPhi1Corr0