alphaEqn.H
Go to the documentation of this file.
1 {
2  word alphaScheme("div(phi,alpha)");
3  word alpharScheme("div(phirb,alpha)");
4 
5  // Standard face-flux compression coefficient
6  surfaceScalarField phic(interface.cAlpha()*mag(phi/mesh.magSf()));
7 
8  // Add the optional isotropic compression contribution
9  if (icAlpha > 0)
10  {
11  phic *= (1.0 - icAlpha);
12  phic += (interface.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
13  }
14 
15  // Do not compress interface at non-coupled boundary faces
16  // (inlets, outlets etc.)
17  forAll(phic.boundaryField(), patchi)
18  {
19  fvsPatchScalarField& phicp = phic.boundaryField()[patchi];
20 
21  if (!phicp.coupled())
22  {
23  phicp == 0;
24  }
25  }
26 
27  tmp<surfaceScalarField> tphiAlpha;
28 
29  if (MULESCorr)
30  {
31  fvScalarMatrix alpha1Eqn
32  (
33  #ifdef LTSSOLVE
34  fv::localEulerDdtScheme<scalar>(mesh, rDeltaT.name()).fvmDdt(alpha1)
35  #else
36  fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
37  #endif
38  + fv::gaussConvectionScheme<scalar>
39  (
40  mesh,
41  phi,
42  upwind<scalar>(mesh, phi)
43  ).fvmDiv(phi, alpha1)
44  );
45 
46  alpha1Eqn.solve();
47 
48  Info<< "Phase-1 volume fraction = "
49  << alpha1.weightedAverage(mesh.Vsc()).value()
50  << " Min(alpha1) = " << min(alpha1).value()
51  << " Max(alpha1) = " << max(alpha1).value()
52  << endl;
53 
54  tmp<surfaceScalarField> tphiAlphaUD(alpha1Eqn.flux());
55  tphiAlpha = tmp<surfaceScalarField>
56  (
57  new surfaceScalarField(tphiAlphaUD())
58  );
59 
60  if (alphaApplyPrevCorr && tphiAlphaCorr0.valid())
61  {
62  Info<< "Applying the previous iteration compression flux" << endl;
63  #ifdef LTSSOLVE
65  #else
67  #endif
68 
70  }
71 
72  // Cache the upwind-flux
73  tphiAlphaCorr0 = tphiAlphaUD;
74 
75  alpha2 = 1.0 - alpha1;
76 
77  interface.correct();
78  }
79 
80  for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
81  {
82  surfaceScalarField phir(phic*interface.nHatf());
83 
84  tmp<surfaceScalarField> tphiAlphaUn
85  (
86  fvc::flux
87  (
88  phi,
89  alpha1,
90  alphaScheme
91  )
92  + fvc::flux
93  (
95  alpha1,
97  )
98  );
99 
100  if (MULESCorr)
101  {
102  tmp<surfaceScalarField> tphiAlphaCorr(tphiAlphaUn() - tphiAlpha());
103  volScalarField alpha10(alpha1);
104 
105  #ifdef LTSSOLVE
106  MULES::LTScorrect(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
107  #else
108  MULES::correct(alpha1, tphiAlphaUn(), tphiAlphaCorr(), 1, 0);
109  #endif
110 
111  // Under-relax the correction for all but the 1st corrector
112  if (aCorr == 0)
113  {
114  tphiAlpha() += tphiAlphaCorr();
115  }
116  else
117  {
118  alpha1 = 0.5*alpha1 + 0.5*alpha10;
119  tphiAlpha() += 0.5*tphiAlphaCorr();
120  }
121  }
122  else
123  {
124  tphiAlpha = tphiAlphaUn;
125 
126  #ifdef LTSSOLVE
128  #else
130  #endif
131  }
132 
133  alpha2 = 1.0 - alpha1;
134 
135  interface.correct();
136  }
137 
139 
141  {
143  }
144 
145  Info<< "Phase-1 volume fraction = "
146  << alpha1.weightedAverage(mesh.Vsc()).value()
147  << " Min(alpha1) = " << min(alpha1).value()
148  << " Max(alpha1) = " << max(alpha1).value()
149  << endl;
150 }