alphaEqn.H
Go to the documentation of this file.
1 {
2  word alphaScheme("div(phi,alpha)");
3  word alpharScheme("div(phirb,alpha)");
4 
5  tmp<fv::ddtScheme<scalar>> ddtAlpha
6  (
8  (
9  mesh,
10  mesh.ddtScheme("ddt(alpha)")
11  )
12  );
13 
14  // Set the off-centering coefficient according to ddt scheme
15  scalar ocCoeff = 0;
16  if
17  (
18  isType<fv::EulerDdtScheme<scalar>>(ddtAlpha())
19  || isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha())
20  )
21  {
22  ocCoeff = 0;
23  }
24  else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha()))
25  {
26  if (nAlphaSubCycles > 1)
27  {
29  << "Sub-cycling is not supported "
30  "with the CrankNicolson ddt scheme"
31  << exit(FatalError);
32  }
33 
34  ocCoeff =
35  refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha())
36  .ocCoeff();
37  }
38  else
39  {
41  << "Only Euler and CrankNicolson ddt schemes are supported"
42  << exit(FatalError);
43  }
44 
45  scalar cnCoeff = 1.0/(1.0 + ocCoeff);
46 
47  // Standard face-flux compression coefficient
48  surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()));
49 
50  // Add the optional isotropic compression contribution
51  if (icAlpha > 0)
52  {
53  phic *= (1.0 - icAlpha);
54  phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
55  }
56 
57  surfaceScalarField::Boundary& phicBf =
58  phic.boundaryFieldRef();
59 
60  // Do not compress interface at non-coupled boundary faces
61  // (inlets, outlets etc.)
62  forAll(phic.boundaryField(), patchi)
63  {
64  fvsPatchScalarField& phicp = phicBf[patchi];
65 
66  if (!phicp.coupled())
67  {
68  phicp == 0;
69  }
70  }
71 
72  tmp<surfaceScalarField> phiCN(phi);
73 
74  // Calculate the Crank-Nicolson off-centred volumetric flux
75  if (ocCoeff > 0)
76  {
77  phiCN = cnCoeff*phi + (1.0 - cnCoeff)*phi.oldTime();
78  }
79 
80  if (MULESCorr)
81  {
82  fvScalarMatrix alpha1Eqn
83  (
84  (
85  LTS
86  ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
87  : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
88  )
89  + fv::gaussConvectionScheme<scalar>
90  (
91  mesh,
92  phiCN,
93  upwind<scalar>(mesh, phiCN)
94  ).fvmDiv(phiCN, alpha1)
95  );
96 
97  alpha1Eqn.solve();
98 
99  Info<< "Phase-1 volume fraction = "
100  << alpha1.weightedAverage(mesh.Vsc()).value()
101  << " Min(" << alpha1.name() << ") = " << min(alpha1).value()
102  << " Max(" << alpha1.name() << ") = " << max(alpha1).value()
103  << endl;
104 
105  tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
106  alphaPhi = talphaPhiUD();
107 
108  if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
109  {
110  Info<< "Applying the previous iteration compression flux" << endl;
111  MULES::correct(alpha1, alphaPhi, talphaPhiCorr0.ref(), 1, 0);
112 
113  alphaPhi += talphaPhiCorr0();
114  }
115 
116  // Cache the upwind-flux
117  talphaPhiCorr0 = talphaPhiUD;
118 
119  alpha2 = 1.0 - alpha1;
120 
121  mixture.correct();
122  }
123 
124 
125  for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
126  {
128 
129  tmp<surfaceScalarField> talphaPhiUn
130  (
131  fvc::flux
132  (
133  phi,
134  alpha1,
135  alphaScheme
136  )
137  + fvc::flux
138  (
140  alpha1,
142  )
143  );
144 
145  // Calculate the Crank-Nicolson off-centred alpha flux
146  if (ocCoeff > 0)
147  {
148  talphaPhiUn =
149  cnCoeff*talphaPhiUn + (1.0 - cnCoeff)*alphaPhi.oldTime();
150  }
151 
152  if (MULESCorr)
153  {
154  tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
155  volScalarField alpha10("alpha10", alpha1);
156 
157  MULES::correct(alpha1, talphaPhiUn(), talphaPhiCorr.ref(), 1, 0);
158 
159  // Under-relax the correction for all but the 1st corrector
160  if (aCorr == 0)
161  {
162  alphaPhi += talphaPhiCorr();
163  }
164  else
165  {
166  alpha1 = 0.5*alpha1 + 0.5*alpha10;
167  alphaPhi += 0.5*talphaPhiCorr();
168  }
169  }
170  else
171  {
172  alphaPhi = talphaPhiUn;
173 
174  MULES::explicitSolve(alpha1, phiCN, alphaPhi, 1, 0);
175  }
176 
177  alpha2 = 1.0 - alpha1;
178 
179  mixture.correct();
180  }
181 
183  {
185  }
186 
187  if
188  (
189  word(mesh.ddtScheme("ddt(rho,U)"))
190  == fv::EulerDdtScheme<vector>::typeName
191  )
192  {
193  rhoPhi = alphaPhi*(rho1 - rho2) + phiCN*rho2;
194  }
195  else
196  {
197  if (ocCoeff > 0)
198  {
199  // Calculate the end-of-time-step alpha flux
200  alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
201  }
202 
203  // Calculate the end-of-time-step mass flux
204  rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
205  }
206 
207  Info<< "Phase-1 volume fraction = "
208  << alpha1.weightedAverage(mesh.Vsc()).value()
209  << " Min(" << alpha1.name() << ") = " << min(alpha1).value()
210  << " Max(" << alpha1.name() << ") = " << max(alpha1).value()
211  << endl;
212 }
surfaceScalarField & phi
fvsPatchField< scalar > fvsPatchScalarField
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
U
Definition: pEqn.H:83
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
word alpharScheme("div(phirb,alpha)")
tmp< surfaceScalarField > talphaPhiCorr0
Definition: createFields.H:132
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
tmp< surfaceScalarField > interpolate(const RhoType &rho)
forAll(phic.boundaryField(), patchi)
Definition: alphaEqn.H:62
label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")))
scalar ocCoeff
Definition: alphaEqn.H:15
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
rho1
Definition: pEqn.H:114
scalar icAlpha(alphaControls.lookupOrDefault< scalar >("icAlpha", 0))
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
dynamicFvMesh & mesh
rhoPhi
Definition: alphaEqn.H:116
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
rho2
Definition: pEqn.H:115
volScalarField alpha10("alpha10", alpha1)
bool isType(const Type &t)
Check the typeid.
Definition: typeInfo.H:126
surfaceScalarField phir("phir", phic *interface.nHatf())
tmp< surfaceScalarField > phiCN(phi)
alpha2
Definition: alphaEqn.H:112
bool LTS
Definition: createRDeltaT.H:1
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(psi *p+alphal *rhol0+((alphav *psiv+alphal *psil)-psi)*pSat, rhoMin);#1"/home/ubuntu/OpenFOAM-4.1/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:72
Info<< "Reading field p_rgh\n"<< endl;volScalarField p_rgh(IOobject("p_rgh", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);Info<< "Reading field U\n"<< endl;volVectorField U(IOobject("U", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);Info<< "Creating phaseChangeTwoPhaseMixture\n"<< endl;autoPtr< phaseChangeTwoPhaseMixture > mixture
Definition: createFields.H:33
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
volScalarField & alpha1
label patchi
bool alphaApplyPrevCorr(alphaControls.lookupOrDefault< Switch >("alphaApplyPrevCorr", false))
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()))
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
Definition: fvcFlux.C:32
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
surfaceScalarField alphaPhi(phi.name()+alpha1.name(), fvc::flux(phi, alpha1, alphaScheme))
tmp< fv::ddtScheme< scalar > > ddtAlpha(fv::ddtScheme< scalar >::New(mesh, mesh.ddtScheme("ddt(alpha)")))
scalar cnCoeff
Definition: alphaEqn.H:45
bool MULESCorr(alphaControls.lookupOrDefault< Switch >("MULESCorr", false))
label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")))