FreeFOAM The Cross-Platform CFD Toolkit
pEqn.H
Go to the documentation of this file.
1 {
2  volScalarField rUA = 1.0/UEqn.A();
4 
5  tmp<fvScalarMatrix> p_rghEqnComp;
6 
7  if (transonic)
8  {
9  p_rghEqnComp =
10  (
12  + fvm::div(phi, p_rgh)
14  );
15  }
16  else
17  {
18  p_rghEqnComp =
19  (
21  + fvc::div(phi, p_rgh)
23  );
24  }
25 
26 
27  U = rUA*UEqn.H();
28 
30  (
31  "phiU",
32  (fvc::interpolate(U) & mesh.Sf())
33  + fvc::ddtPhiCorr(rUA, rho, U, phi)
34  );
35 
36  phi = phiU +
37  (
38  fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
40  )*rUAf*mesh.magSf();
41 
42  for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
43  {
44  fvScalarMatrix p_rghEqnIncomp
45  (
46  fvc::div(phi)
47  - fvm::laplacian(rUAf, p_rgh)
48  );
49 
50  solve
51  (
52  (
53  max(alpha1, scalar(0))*(psi1/rho1)
54  + max(alpha2, scalar(0))*(psi2/rho2)
55  )
56  *p_rghEqnComp()
57  + p_rghEqnIncomp,
58  mesh.solver
59  (
60  p_rgh.select
61  (
62  oCorr == nOuterCorr-1
63  && corr == nCorr-1
64  && nonOrth == nNonOrthCorr
65  )
66  )
67  );
68 
69  if (nonOrth == nNonOrthCorr)
70  {
71  dgdt =
72  (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
73  *(p_rghEqnComp & p_rgh);
74  phi += p_rghEqnIncomp.flux();
75  }
76  }
77 
78  U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
79  U.correctBoundaryConditions();
80 
81  p = max
82  (
83  (p_rgh + gh*(alpha1*rho10 + alpha2*rho20))
84  /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
85  pMin
86  );
87 
88  rho1 = rho10 + psi1*p;
89  rho2 = rho20 + psi2*p;
90 
91  Info<< "max(U) " << max(mag(U)).value() << endl;
92  Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
93 
94  // Make the fluxes relative to the mesh motion
95  fvc::makeRelative(phi, U);
96 }