Difference between revisions of "BubbleFoam"
From OpenFOAMWiki
(→Implementation of the phase continuity equation) |
(→Implementation of the pressure equation) |
||
| Line 329: | Line 329: | ||
==== Implementation of the pressure equation ==== | ==== Implementation of the pressure equation ==== | ||
| + | |||
| + | The details of the implementation of the pressure equation are commented in the code snippet below, where the implementation of the pseudo-staggered algorithm developed by Weller [10, 8] to obtain a stable solution in completely separated flows is adopted. | ||
| + | |||
| + | <cpp> | ||
| + | { | ||
| + | // The phase fraction fields are interpolated at cell faces | ||
| + | surfaceScalarField alphaf = fvc::interpolate(alpha); | ||
| + | surfaceScalarField betaf = scalar(1) - alphaf; | ||
| + | |||
| + | // The central coeffcients are computed | ||
| + | volScalarField rUaA = 1.0/UaEqn.A(); | ||
| + | volScalarField rUbA = 1.0/UbEqn.A(); | ||
| + | |||
| + | // The values of the central coefficients are interpolated at cell faces | ||
| + | surfaceScalarField rUaAf = fvc::interpolate(rUaA); | ||
| + | surfaceScalarField rUbAf = fvc::interpolate(rUbA); | ||
| + | |||
| + | // The phase velocities are predicted | ||
| + | Ua = rUaA*UaEqn.H(); | ||
| + | Ub = rUbA*UbEqn.H(); | ||
| + | |||
| + | // Accounting for the explict drag term and buoyancy effect. Both | ||
| + | // these terms are managed at cell faces and their effect is included in the phase flux. | ||
| + | // Note that the whole drag term is interpolated and then multiplied by the flux, | ||
| + | // instead of reusing the previously computed value of the central coefficient at cell | ||
| + | // faces (rUaAf). | ||
| + | |||
| + | surfaceScalarField phiDraga = | ||
| + | fvc::interpolate(beta/rhoa*dragCoef*rUaA)*phib + | ||
| + | rUaAf*(g & mesh.Sf()); | ||
| + | surfaceScalarField phiDragb = | ||
| + | fvc::interpolate(alpha/rhob*dragCoef*rUbA)*phia + | ||
| + | rUbAf*(g & mesh.Sf()); | ||
| + | |||
| + | forAll(p.boundaryField(), patchi) | ||
| + | { | ||
| + | if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi])) | ||
| + | { | ||
| + | phiDraga.boundaryField()[patchi] = 0.0; | ||
| + | phiDragb.boundaryField()[patchi] = 0.0; | ||
| + | } | ||
| + | } | ||
| + | |||
| + | // The phase fluxes and the total flux are predicted. | ||
| + | phia = (fvc::interpolate(Ua) & mesh.Sf()) + | ||
| + | fvc::ddtPhiCorr(rUaA, Ua, phia) + phiDraga; | ||
| + | phib = (fvc::interpolate(Ub) & mesh.Sf()) | ||
| + | + fvc::ddtPhiCorr(rUbA, Ub, phib) + phiDragb; | ||
| + | |||
| + | phi = alphaf*phia + betaf*phib; | ||
| + | |||
| + | // The mixture central coefficient is computed | ||
| + | surfaceScalarField Dp("(rho*(1|A(U)))", alphaf*rUaAf/rhoa + | ||
| + | betaf*rUbAf/rhob); | ||
| + | |||
| + | // The non-ortogonal correction loop is started | ||
| + | for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) | ||
| + | { | ||
| + | // The pressure equation is set up | ||
| + | fvScalarMatrix pEqn | ||
| + | ( | ||
| + | fvm::laplacian(Dp, p) == fvc::div(phi) | ||
| + | ); | ||
| + | |||
| + | // References for the pressure are set | ||
| + | pEqn.setReference(pRefCell, pRefValue); | ||
| + | |||
| + | // The pressure equation is solved | ||
| + | pEqn.solve(); | ||
| + | |||
| + | if (nonOrth == nNonOrthCorr) | ||
| + | { | ||
| + | // The surface normal pressure gradient is computed | ||
| + | surfaceScalarField SfGradp = pEqn.flux()/Dp; | ||
| + | |||
| + | // The phase fluxes and the total flux are corrected to account for the updated | ||
| + | // pressure field. | ||
| + | phia -= rUaAf*SfGradp/rhoa; | ||
| + | phib -= rUbAf*SfGradp/rhob; | ||
| + | |||
| + | phi = alphaf*phia + betaf*phib; | ||
| + | |||
| + | // Pressure is explicitly relaxed before correcting the phase velocity fields | ||
| + | p.relax(); | ||
| + | |||
| + | // The pressure surface normal gradient is recomputed. | ||
| + | SfGradp = pEqn.flux()/Dp; | ||
| + | |||
| + | // Velocities are not updated using directly the pressure, as done conventionally, | ||
| + | // but reconstructing the velocity update from the fluxes, which are the primary | ||
| + | // variable [8] | ||
| + | |||
| + | Ua += (fvc::reconstruct(phiDraga - rUaAf*SfGradp/rhoa)); | ||
| + | //Ua += rUaA*(fvc::reconstruct(phiDraga/rUaAf - | ||
| + | // SfGradp/rhoa)); | ||
| + | Ua.correctBoundaryConditions(); | ||
| + | |||
| + | Ub += (fvc::reconstruct(phiDragb - rUbAf*SfGradp/rhob)); | ||
| + | //Ub += rUbA*(fvc::reconstruct(phiDragb/rUbAf - | ||
| + | // SfGradp/rhob)); | ||
| + | Ub.correctBoundaryConditions(); | ||
| + | |||
| + | U = alpha*Ua + beta*Ub; | ||
| + | } | ||
| + | } | ||
| + | } | ||
| + | |||
| + | #include "continuityErrs.H" | ||
| + | |||
| + | </cpp> | ||
== bubbleFoam setup and solution strategy == | == bubbleFoam setup and solution strategy == | ||
