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<fv::CrankNicolsonDdtScheme<scalar> >(ddtAlpha()).ocCoeff();
36  }
37  else
38  {
40  << "Only Euler and CrankNicolson ddt schemes are supported"
41  << exit(FatalError);
42  }
43 
44  scalar cnCoeff = 1.0/(1.0 + ocCoeff);
45 
46  // Standard face-flux compression coefficient
47  surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()));
48 
49  // Add the optional isotropic compression contribution
50  if (icAlpha > 0)
51  {
52  phic *= (1.0 - icAlpha);
53  phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
54  }
55 
56  // Do not compress interface at non-coupled boundary faces
57  // (inlets, outlets etc.)
58  forAll(phic.boundaryField(), patchi)
59  {
60  fvsPatchScalarField& phicp = phic.boundaryField()[patchi];
61 
62  if (!phicp.coupled())
63  {
64  phicp == 0;
65  }
66  }
67 
68  tmp<surfaceScalarField> phiCN(phi);
69 
70  // Calculate the Crank-Nicolson off-centred volumetric flux
71  if (ocCoeff > 0)
72  {
73  phiCN = cnCoeff*phi + (1.0 - cnCoeff)*phi.oldTime();
74  }
75 
76  if (MULESCorr)
77  {
78  fvScalarMatrix alpha1Eqn
79  (
80  (
81  LTS
82  ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
83  : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
84  )
85  + fv::gaussConvectionScheme<scalar>
86  (
87  mesh,
88  phiCN,
89  upwind<scalar>(mesh, phiCN)
90  ).fvmDiv(phiCN, 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> talphaPhiUD(alpha1Eqn.flux());
102  alphaPhi = talphaPhiUD();
103 
104  if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
105  {
106  Info<< "Applying the previous iteration compression flux" << endl;
107  MULES::correct(alpha1, alphaPhi, talphaPhiCorr0(), 1, 0);
108 
109  alphaPhi += talphaPhiCorr0();
110  }
111 
112  // Cache the upwind-flux
113  talphaPhiCorr0 = talphaPhiUD;
114 
115  alpha2 = 1.0 - alpha1;
116 
117  mixture.correct();
118  }
119 
120 
121  for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
122  {
124 
125  tmp<surfaceScalarField> talphaPhiUn
126  (
127  fvc::flux
128  (
129  phi,
130  alpha1,
131  alphaScheme
132  )
133  + fvc::flux
134  (
136  alpha1,
138  )
139  );
140 
141  // Calculate the Crank-Nicolson off-centred alpha flux
142  if (ocCoeff > 0)
143  {
144  talphaPhiUn =
145  cnCoeff*talphaPhiUn + (1.0 - cnCoeff)*alphaPhi.oldTime();
146  }
147 
148  if (MULESCorr)
149  {
150  tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
151  volScalarField alpha10("alpha10", alpha1);
152 
153  MULES::correct(alpha1, talphaPhiUn(), talphaPhiCorr(), 1, 0);
154 
155  // Under-relax the correction for all but the 1st corrector
156  if (aCorr == 0)
157  {
158  alphaPhi += talphaPhiCorr();
159  }
160  else
161  {
162  alpha1 = 0.5*alpha1 + 0.5*alpha10;
163  alphaPhi += 0.5*talphaPhiCorr();
164  }
165  }
166  else
167  {
168  alphaPhi = talphaPhiUn;
169 
170  MULES::explicitSolve(alpha1, phiCN, alphaPhi, 1, 0);
171  }
172 
173  alpha2 = 1.0 - alpha1;
174 
175  mixture.correct();
176  }
177 
179  {
181  }
182 
183  if
184  (
185  word(mesh.ddtScheme("ddt(rho,U)"))
186  == fv::EulerDdtScheme<vector>::typeName
187  )
188  {
189  rhoPhi = alphaPhi*(rho1 - rho2) + phiCN*rho2;
190  }
191  else
192  {
193  if (ocCoeff > 0)
194  {
195  // Calculate the end-of-time-step alpha flux
196  alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
197  }
198 
199  // Calculate the end-of-time-step mass flux
200  rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
201  }
202 
203  Info<< "Phase-1 volume fraction = "
204  << alpha1.weightedAverage(mesh.Vsc()).value()
205  << " Min(" << alpha1.name() << ") = " << min(alpha1).value()
206  << " Max(" << alpha1.name() << ") = " << max(alpha1).value()
207  << endl;
208 }
tmp< surfaceScalarField > talphaPhiCorr0
Definition: createFields.H:138
tmp< surfaceScalarField > phiCN(phi)
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > flux(const surfaceScalarField &phi, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcFlux.C:45
dimensioned< scalar > mag(const dimensioned< Type > &)
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/dm2/henry/OpenFOAM/OpenFOAM-3.0.x/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:74
bool LTS
Definition: createRDeltaT.H:1
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
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
scalar cnCoeff
Definition: alphaEqn.H:44
surfaceScalarField alphaPhi(phi.name()+alpha1.name(), fvc::flux( phi, alpha1, alphaScheme ))
messageStream Info
tmp< surfaceScalarField > interpolate(const RhoType &rho)
dynamicFvMesh & mesh
scalar ocCoeff
Definition: alphaEqn.H:15
autoPtr< BasicCompressibleTurbulenceModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const typename BasicCompressibleTurbulenceModel::transportModel &transport, const word &propertiesName)
fvsPatchField< scalar > fvsPatchScalarField
surfaceScalarField phir("phir", phic *interface.nHatf())
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()))
label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")))
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
scalar icAlpha(alphaControls.lookupOrDefault< scalar >("icAlpha", 0))
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
rho1
Definition: pEqn.H:124
label patchi
bool alphaApplyPrevCorr(alphaControls.lookupOrDefault< Switch >("alphaApplyPrevCorr", false))
rho2
Definition: pEqn.H:125
volScalarField & alpha1
Definition: createFields.H:15
bool isType(const Type &t)
Check the typeid.
Definition: typeInfo.H:126
surfaceScalarField & phi
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
const word & executable() const
Name of executable without the path.
Definition: argListI.H:30
alpha2
Definition: alphaEqn.H:112
Foam::argList args(argc, argv)
label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")))
error FatalError
bool MULESCorr(alphaControls.lookupOrDefault< Switch >("MULESCorr", false))
volScalarField alpha10("alpha10", alpha1)
rhoPhi
Definition: alphaEqn.H:116
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
word alpharScheme("div(phirb,alpha)")
tmp< fv::ddtScheme< scalar > > ddtAlpha(fv::ddtScheme< scalar >::New( mesh, mesh.ddtScheme("ddt(alpha)") ))
U
Definition: pEqn.H:82
forAll(phic.boundaryField(), patchi)
Definition: alphaEqn.H:58