FreeFOAM The Cross-Platform CFD Toolkit
alphaEqns.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  interface.cAlpha()*mag(phi/mesh.magSf())*interface.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
43  surfaceScalarField phiAlpha1 =
44  fvc::flux
45  (
46  phi,
47  alpha1,
48  alphaScheme
49  )
50  + fvc::flux
51  (
52  -fvc::flux(-phir, alpha2, alpharScheme),
53  alpha1,
54  alpharScheme
55  )
56  + fvc::flux
57  (
58  -fvc::flux(-phir, alpha3, alpharScheme),
59  alpha1,
60  alpharScheme
61  );
62 
63  // Create the bounded (upwind) flux for alpha1
64  surfaceScalarField phiAlpha1BD =
65  upwind<scalar>(mesh, phi).flux(alpha1);
66 
67  // Calculate the flux correction for alpha1
68  phiAlpha1 -= phiAlpha1BD;
69 
70  // Calculate the limiter for alpha1
72  (
73  allLambda,
74  geometricOneField(),
75  alpha1,
76  phiAlpha1BD,
77  phiAlpha1,
78  zeroField(),
79  zeroField(),
80  1,
81  0,
82  3
83  );
84 
85  // Create the complete flux for alpha2
86  surfaceScalarField phiAlpha2 =
87  fvc::flux
88  (
89  phi,
90  alpha2,
91  alphaScheme
92  )
93  + fvc::flux
94  (
95  -fvc::flux(phir, alpha1, alpharScheme),
96  alpha2,
97  alpharScheme
98  );
99 
100  // Create the bounded (upwind) flux for alpha2
101  surfaceScalarField phiAlpha2BD =
102  upwind<scalar>(mesh, phi).flux(alpha2);
103 
104  // Calculate the flux correction for alpha2
105  phiAlpha2 -= phiAlpha2BD;
106 
107  // Further limit the limiter for alpha2
109  (
110  allLambda,
111  geometricOneField(),
112  alpha2,
113  phiAlpha2BD,
114  phiAlpha2,
115  zeroField(),
116  zeroField(),
117  1,
118  0,
119  3
120  );
121 
122  // Construct the limited fluxes
123  phiAlpha1 = phiAlpha1BD + lambda*phiAlpha1;
124  phiAlpha2 = phiAlpha2BD + lambda*phiAlpha2;
125 
126  // Solve for alpha1
127  solve(fvm::ddt(alpha1) + fvc::div(phiAlpha1));
128 
129  // Create the diffusion coefficients for alpha2<->alpha3
130  volScalarField Dc23 = D23*max(alpha3, scalar(0))*pos(alpha2);
131  volScalarField Dc32 = D23*max(alpha2, scalar(0))*pos(alpha3);
132 
133  // Add the diffusive flux for alpha3->alpha2
134  phiAlpha2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);
135 
136  // Solve for alpha2
137  fvScalarMatrix alpha2Eqn
138  (
139  fvm::ddt(alpha2)
140  + fvc::div(phiAlpha2)
141  - fvm::laplacian(Dc23 + Dc32, alpha2)
142  );
143  alpha2Eqn.solve();
144 
145  // Construct the complete mass flux
146  rhoPhi =
147  phiAlpha1*(rho1 - rho3)
148  + (phiAlpha2 + alpha2Eqn.flux())*(rho2 - rho3)
149  + phi*rho3;
150 
151  alpha3 = 1.0 - alpha1 - alpha2;
152  }
153 
154  Info<< "Air phase volume fraction = "
155  << alpha1.weightedAverage(mesh.V()).value()
156  << " Min(alpha1) = " << min(alpha1).value()
157  << " Max(alpha1) = " << max(alpha1).value()
158  << endl;
159 
160  Info<< "Liquid phase volume fraction = "
161  << alpha2.weightedAverage(mesh.V()).value()
162  << " Min(alpha2) = " << min(alpha2).value()
163  << " Max(alpha2) = " << max(alpha2).value()
164  << endl;
165 }