10 Psi2 = rho_foam/p_rgh;
12 volScalarField rAU(
"rAU", 1.0/UEqn.A());
13 surfaceScalarField rAUf(
"rAUf", fvc::interpolate(rAU));
15 volVectorField HbyA(
"HbyA", U);
18 surfaceScalarField phiHbyA
21 (fvc::interpolate(HbyA) & mesh.Sf())
22 + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
25 surfaceScalarField phig
28 interface.surfaceTensionForce()
30 - ghf*fvc::snGrad(rho)
36 constrainPressure(p_rgh, U, phiHbyA, rAUf);
38 tmp<fvScalarMatrix> p_rghEqnComp1;
39 tmp<fvScalarMatrix> p_rghEqnComp2;
41 if (pimple.transonic())
43 surfaceScalarField phid1(
"phid1", fvc::interpolate(Psi1)*phi);
44 surfaceScalarField phid2(
"phid2", fvc::interpolate(Psi2)*phi);
47 fvc::ddt(rho_gas) + fvc::div(phi, rho_gas)
48 - fvc::Sp(fvc::div(phi), rho_gas)
52 + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
54 deleteDemandDrivenData(p_rghEqnComp1.ref().faceFluxCorrectionPtr());
55 p_rghEqnComp1.ref().relax();
58 fvc::ddt(rho_foam) + fvc::div(phi, rho_foam)
59 - fvc::Sp(fvc::div(phi), rho_foam)
63 + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
65 deleteDemandDrivenData(p_rghEqnComp2.ref().faceFluxCorrectionPtr());
66 p_rghEqnComp2.ref().relax();
73 fvc::ddt(alpha1, rho_gas)
74 + fvc::div(fvc::interpolate(alpha1)*phi, rho_gas)
75 - fvc::Sp(fvc::ddt(alpha1)
76 + fvc::div(fvc::interpolate(alpha1)*phi), rho_gas)
78 + (alpha1*Psi1/rho_gas)*correction(fvm::ddt(p_rgh));
82 fvc::ddt(alpha2, rho_foam)
83 + fvc::div(fvc::interpolate(alpha2)*phi, rho_foam)
84 - fvc::Sp(fvc::ddt(alpha2)
85 + fvc::div(fvc::interpolate(alpha2)*phi), rho_foam)
87 + (alpha2*Psi2/rho_foam)*correction(fvm::ddt(p_rgh));
91 volScalarField p_rgh_0(p_rgh);
93 while (pimple.correctNonOrthogonal())
95 fvScalarMatrix p_rghEqnIncomp
98 - fvm::laplacian(rAUf, p_rgh)
103 p_rghEqnComp1() + p_rghEqnComp2()
106 mesh.solver(p_rgh.select(pimple.finalInnerIter()))
108 if (pimple.finalNonOrthogonalIter())
110 p = max(p_rgh + (alpha1*rho_gas + (1.0 - alpha1)*rho_foam)*gh, pMin);
111 p_rgh = p - (alpha1*rho_gas + (1.0 - alpha1)*rho_foam)*gh;
112 dimensionedScalar dumdgdt (
"dumdgdt",dimensionSet(1,-3,0,0,0,0,0), 1.0);
116 pos((1.0 - alpha1))*(p_rghEqnComp2 & p_rgh)/rho_foam
117 - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho_gas
120 phi = phiHbyA + p_rghEqnIncomp.flux();
122 + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
123 U.correctBoundaryConditions();
126 p = max(p_rgh + (alpha1*rho_gas + (1.0 - alpha1)*rho_foam)*gh, pMin);
129 rho_gas += Psi1*(p_rgh - p_rgh_0);
130 rho_foam += Psi2*(p_rgh - p_rgh_0);
131 rho = alpha1*rho_gas + (1 - alpha1)*rho_foam;
133 Info<<
"max(U) " << max(mag(U)).value() << endl;
134 Info<<
"min(p_rgh) " << min(p_rgh).value() << endl;