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 == |