FreeFOAM The Cross-Platform CFD Toolkit
pEqn.H
Go to the documentation of this file.
1 {
2  surfaceScalarField alphaf = fvc::interpolate(alpha);
3  surfaceScalarField betaf = scalar(1) - alphaf;
4 
5  volScalarField rUaA = 1.0/UaEqn.A();
6  volScalarField rUbA = 1.0/UbEqn.A();
7 
8  phia == (fvc::interpolate(Ua) & mesh.Sf());
9  phib == (fvc::interpolate(Ub) & mesh.Sf());
10 
11  rUaAf = fvc::interpolate(rUaA);
13 
14  Ua = rUaA*UaEqn.H();
15  Ub = rUbA*UbEqn.H();
16 
17  surfaceScalarField phiDraga =
18  fvc::interpolate(beta/rhoa*K*rUaA)*phib + rUaAf*(g & mesh.Sf());
19 
20  if (g0.value() > 0.0)
21  {
22  phiDraga -= ppMagf*fvc::snGrad(alpha)*mesh.magSf();
23  }
24 
25  if (kineticTheory.on())
26  {
27  phiDraga -= rUaAf*fvc::snGrad(kineticTheory.pa()/rhoa)*mesh.magSf();
28  }
29 
30  surfaceScalarField phiDragb =
31  fvc::interpolate(alpha/rhob*K*rUbA)*phia + rUbAf*(g & mesh.Sf());
32 
33  // Fix for gravity on outlet boundary.
34  forAll(p.boundaryField(), patchi)
35  {
36  if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi]))
37  {
38  phiDraga.boundaryField()[patchi] = 0.0;
39  phiDragb.boundaryField()[patchi] = 0.0;
40  }
41  }
42 
43  phia = (fvc::interpolate(Ua) & mesh.Sf()) + fvc::ddtPhiCorr(rUaA, Ua, phia)
44  + phiDraga;
45 
46  phib = (fvc::interpolate(Ub) & mesh.Sf()) + fvc::ddtPhiCorr(rUbA, Ub, phib)
47  + phiDragb;
48 
49  phi = alphaf*phia + betaf*phib;
50 
52  (
53  "(rho*(1|A(U)))",
54  alphaf*rUaAf/rhoa
55  + betaf*rUbAf/rhob
56  );
57 
58  for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
59  {
60  fvScalarMatrix pEqn
61  (
62  fvm::laplacian(Dp, p) == fvc::div(phi)
63  );
64 
65  pEqn.setReference(pRefCell, pRefValue);
66  pEqn.solve();
67 
68  if (nonOrth == nNonOrthCorr)
69  {
70  surfaceScalarField SfGradp = pEqn.flux()/Dp;
71 
72  phia -= rUaAf*SfGradp/rhoa;
73  phib -= rUbAf*SfGradp/rhob;
74  phi = alphaf*phia + betaf*phib;
75 
76  p.relax();
77  SfGradp = pEqn.flux()/Dp;
78 
79  Ua += fvc::reconstruct(phiDraga - rUaAf*SfGradp/rhoa);
80  Ua.correctBoundaryConditions();
81 
82  Ub += fvc::reconstruct(phiDragb - rUbAf*SfGradp/rhob);
83  Ub.correctBoundaryConditions();
84 
85  U = alpha*Ua + beta*Ub;
86  }
87  }
88 }
89