00001 { 00002 surfaceScalarField alphaf = fvc::interpolate(alpha); 00003 surfaceScalarField betaf = scalar(1) - alphaf; 00004 00005 volScalarField rUaA = 1.0/UaEqn.A(); 00006 volScalarField rUbA = 1.0/UbEqn.A(); 00007 00008 surfaceScalarField rUaAf = fvc::interpolate(rUaA); 00009 surfaceScalarField rUbAf = fvc::interpolate(rUbA); 00010 00011 Ua = rUaA*UaEqn.H(); 00012 Ub = rUbA*UbEqn.H(); 00013 00014 surfaceScalarField phiDraga = 00015 fvc::interpolate(beta/rhoa*dragCoef*rUaA)*phib + rUaAf*(g & mesh.Sf()); 00016 surfaceScalarField phiDragb = 00017 fvc::interpolate(alpha/rhob*dragCoef*rUbA)*phia + rUbAf*(g & mesh.Sf()); 00018 00019 forAll(p.boundaryField(), patchi) 00020 { 00021 if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi])) 00022 { 00023 phiDraga.boundaryField()[patchi] = 0.0; 00024 phiDragb.boundaryField()[patchi] = 0.0; 00025 } 00026 } 00027 00028 phia = (fvc::interpolate(Ua) & mesh.Sf()) + fvc::ddtPhiCorr(rUaA, Ua, phia) 00029 + phiDraga; 00030 phib = (fvc::interpolate(Ub) & mesh.Sf()) + fvc::ddtPhiCorr(rUbA, Ub, phib) 00031 + phiDragb; 00032 00033 phi = alphaf*phia + betaf*phib; 00034 00035 surfaceScalarField Dp("(rho*(1|A(U)))", alphaf*rUaAf/rhoa + betaf*rUbAf/rhob); 00036 00037 for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) 00038 { 00039 fvScalarMatrix pEqn 00040 ( 00041 fvm::laplacian(Dp, p) == fvc::div(phi) 00042 ); 00043 00044 pEqn.setReference(pRefCell, pRefValue); 00045 pEqn.solve(); 00046 00047 if (nonOrth == nNonOrthCorr) 00048 { 00049 surfaceScalarField SfGradp = pEqn.flux()/Dp; 00050 00051 phia -= rUaAf*SfGradp/rhoa; 00052 phib -= rUbAf*SfGradp/rhob; 00053 phi = alphaf*phia + betaf*phib; 00054 00055 p.relax(); 00056 SfGradp = pEqn.flux()/Dp; 00057 00058 Ua += (fvc::reconstruct(phiDraga - rUaAf*SfGradp/rhoa)); 00059 //Ua += rUaA*(fvc::reconstruct(phiDraga/rUaAf - SfGradp/rhoa)); 00060 Ua.correctBoundaryConditions(); 00061 00062 Ub += (fvc::reconstruct(phiDragb - rUbAf*SfGradp/rhob)); 00063 //Ub += rUbA*(fvc::reconstruct(phiDragb/rUbAf - SfGradp/rhob)); 00064 Ub.correctBoundaryConditions(); 00065 00066 U = alpha*Ua + beta*Ub; 00067 } 00068 } 00069 } 00070 00071 #include <finiteVolume/continuityErrs.H>