alphaEqn.H
Go to the documentation of this file.
1 {
2  word alphaScheme("div(phi,alpha)");
3  word alpharScheme("div(phirb,alpha)");
4 
6  (
7  IOobject
8  (
9  "phir",
10  runTime.timeName(),
11  mesh,
12  IOobject::NO_READ,
13  IOobject::NO_WRITE
14  ),
15  mixture.cAlpha()*mag(phi/mesh.magSf())*mixture.nHatf()
16  );
17 
18  for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
19  {
20  // Create the limiter to be used for all phase-fractions
21  scalarField allLambda(mesh.nFaces(), 1.0);
22 
23  // Split the limiter into a surfaceScalarField
25  (
26  IOobject
27  (
28  "lambda",
29  mesh.time().timeName(),
30  mesh,
31  IOobject::NO_READ,
32  IOobject::NO_WRITE,
33  false
34  ),
35  mesh,
36  dimless,
37  allLambda,
38  false // Use slices for the couples
39  );
40 
41 
42  // Create the complete convection flux for alpha1
44  (
45  fvc::flux
46  (
47  phi,
48  alpha1,
49  alphaScheme
50  )
51  + fvc::flux
52  (
54  alpha1,
56  )
57  + fvc::flux
58  (
59  -fvc::flux(-phir, alpha3, alpharScheme),
60  alpha1,
62  )
63  );
64 
65  // Create the bounded (upwind) flux for alpha1
66  surfaceScalarField alphaPhi1BD
67  (
68  upwind<scalar>(mesh, phi).flux(alpha1)
69  );
70 
71  // Calculate the flux correction for alpha1
72  alphaPhi1 -= alphaPhi1BD;
73 
74  // Calculate the limiter for alpha1
75  if (LTS)
76  {
77  const volScalarField& rDeltaT =
78  fv::localEulerDdt::localRDeltaT(mesh);
79 
81  (
82  allLambda,
83  rDeltaT,
84  geometricOneField(),
85  alpha1,
86  alphaPhi1BD,
87  alphaPhi1,
88  zeroField(),
89  zeroField(),
90  1,
91  0
92  );
93  }
94  else
95  {
97  (
98  allLambda,
99  1.0/runTime.deltaT().value(),
100  geometricOneField(),
101  alpha1,
102  alphaPhi1BD,
103  alphaPhi1,
104  zeroField(),
105  zeroField(),
106  1,
107  0
108  );
109  }
110 
111  // Create the complete flux for alpha2
113  (
114  fvc::flux
115  (
116  phi,
117  alpha2,
118  alphaScheme
119  )
120  + fvc::flux
121  (
123  alpha2,
125  )
126  );
127 
128  // Create the bounded (upwind) flux for alpha2
129  surfaceScalarField alphaPhi2BD
130  (
131  upwind<scalar>(mesh, phi).flux(alpha2)
132  );
133 
134  // Calculate the flux correction for alpha2
135  alphaPhi2 -= alphaPhi2BD;
136 
137  // Further limit the limiter for alpha2
138  if (LTS)
139  {
140  const volScalarField& rDeltaT =
141  fv::localEulerDdt::localRDeltaT(mesh);
142 
144  (
145  allLambda,
146  rDeltaT,
147  geometricOneField(),
148  alpha2,
149  alphaPhi2BD,
150  alphaPhi2,
151  zeroField(),
152  zeroField(),
153  1,
154  0
155  );
156  }
157  else
158  {
160  (
161  allLambda,
162  1.0/runTime.deltaT().value(),
163  geometricOneField(),
164  alpha2,
165  alphaPhi2BD,
166  alphaPhi2,
167  zeroField(),
168  zeroField(),
169  1,
170  0
171  );
172  }
173 
174  // Construct the limited fluxes
175  alphaPhi1 = alphaPhi1BD + lambda*alphaPhi1;
176  alphaPhi2 = alphaPhi2BD + lambda*alphaPhi2;
177 
178  // Solve for alpha1
179  solve(fvm::ddt(alpha1) + fvc::div(alphaPhi1));
180 
181  // Create the diffusion coefficients for alpha2<->alpha3
182  volScalarField Dc23(D23*max(alpha3, scalar(0))*pos0(alpha2));
183  volScalarField Dc32(D23*max(alpha2, scalar(0))*pos0(alpha3));
184 
185  // Add the diffusive flux for alpha3->alpha2
186  alphaPhi2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);
187 
188  // Solve for alpha2
189  fvScalarMatrix alpha2Eqn
190  (
192  + fvc::div(alphaPhi2)
193  - fvm::laplacian(Dc23 + Dc32, alpha2)
194  );
195  alpha2Eqn.solve();
196 
197  // Construct the complete mass flux
198  rhoPhi =
199  alphaPhi1*(rho1 - rho3)
200  + (alphaPhi2 + alpha2Eqn.flux())*(rho2 - rho3)
201  + phi*rho3;
202 
203  alpha3 = 1.0 - alpha1 - alpha2;
204  }
205 
206  Info<< "Air phase volume fraction = "
207  << alpha1.weightedAverage(mesh.V()).value()
208  << " Min(" << alpha1.name() << ") = " << min(alpha1).value()
209  << " Max(" << alpha1.name() << ") = " << max(alpha1).value()
210  << endl;
211 
212  Info<< "Liquid phase volume fraction = "
213  << alpha2.weightedAverage(mesh.V()).value()
214  << " Min(" << alpha2.name() << ") = " << min(alpha2).value()
215  << " Max(" << alpha2.name() << ") = " << max(alpha2).value()
216  << endl;
217 }
surfaceScalarField & phi
void limiter(scalarField &allLambda, const RdeltaTType &rDeltaT, const RhoType &rho, const volScalarField &psi, const surfaceScalarField &phiBD, const surfaceScalarField &phiCorr, const SpType &Sp, const SuType &Su, const scalar psiMax, const scalar psiMin)
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
surfaceScalarField & alphaPhi2
word alpharScheme("div(phirb,alpha)")
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
const surfaceScalarField & alphaPhi1
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")))
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
dynamicFvMesh & mesh
rhoPhi
Definition: alphaEqn.H:116
rhoEqn solve()
alpha2
Definition: alphaEqn.H:112
volScalarField scalarField(fieldObject, mesh)
dimensionedScalar pos0(const dimensionedScalar &ds)
bool LTS
Definition: createRDeltaT.H:1
Info<< "Reading field p_rgh\"<< endl;volScalarField p_rgh(IOobject("p_rgh", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);Info<< "Reading field U\"<< endl;volVectorField U(IOobject("U", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);Info<< "Creating phaseChangeTwoPhaseMixture\"<< endl;autoPtr< phaseChangeTwoPhaseMixture > mixture
Definition: createFields.H:33
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
volScalarField & alpha1
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
surfaceScalarField phir(IOobject("phir", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE), mixture.cAlpha() *mag(phi/mesh.magSf()) *mixture.nHatf())
dimensionedScalar lambda(laminarTransport.lookup("lambda"))
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
const dimensionedScalar & rho2
Definition: createFields.H:40
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
const dimensionedScalar & rho1
Definition: createFields.H:39
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
SlicedGeometricField< scalar, fvsPatchField, slicedFvsPatchField, surfaceMesh > slicedSurfaceScalarField