Difference between revisions of "BubbleFoam"
From OpenFOAMWiki
(→Implementation of the phase momentum equation) |
(→Implementation of the phase continuity equation) |
||
| Line 254: | Line 254: | ||
==== Implementation of the phase continuity equation ==== | ==== Implementation of the phase continuity equation ==== | ||
| + | |||
| + | <cpp> | ||
| + | { | ||
| + | // A word is defined to store the discretisation scheme label used to compute the | ||
| + | // divergence of the dispersed phase volume fraction. The text between “ “ is searched | ||
| + | // for in the fvSchemes dictionary | ||
| + | word scheme("div(phi,alpha)"); | ||
| + | |||
| + | // The relative flux is calculated | ||
| + | surfaceScalarField phir = phia - phib; | ||
| + | |||
| + | // The Courant number is computed on the base of the relative velocity | ||
| + | Info<< "Max Ur Courant Number = " | ||
| + | << ( | ||
| + | max | ||
| + | ( | ||
| + | mesh.surfaceInterpolation::deltaCoeffs()*mag(phir) | ||
| + | /mesh.magSf() | ||
| + | )*runTime.deltaT() | ||
| + | ).value() | ||
| + | << endl; | ||
| + | |||
| + | // The phase fraction correction loop is started | ||
| + | for (int acorr=0; acorr<nAlphaCorr; acorr++) | ||
| + | { | ||
| + | // The dipersed phase continuity equation object is created using | ||
| + | // implicit discretization for all its terms | ||
| + | fvScalarMatrix alphaEqn | ||
| + | ( | ||
| + | fvm::ddt(alpha) | ||
| + | + fvm::div(phi, alpha, scheme) | ||
| + | + fvm::div(-fvc::flux(-phir, beta, scheme), alpha, scheme) | ||
| + | ); | ||
| + | |||
| + | // Under-relaxation is applied to the dispersed phase continuity equation | ||
| + | alphaEqn.relax(); | ||
| + | |||
| + | // The dispersed phase continuity equation is solved | ||
| + | alphaEqn.solve(); | ||
| + | |||
| + | // The continuous phase continuity equation is defined as done for the dispersed | ||
| + | // phase. Notice that by default this equation is commented out because only the | ||
| + | // equation for the dispersed phase is solved. | ||
| + | |||
| + | /* | ||
| + | fvScalarMatrix betaEqn | ||
| + | ( | ||
| + | fvm::ddt(beta) | ||
| + | + fvm::div(phi, beta, scheme) | ||
| + | + fvm::div(-fvc::flux(phir, scalar(1) - beta, scheme), | ||
| + | beta, scheme) | ||
| + | ); | ||
| + | betaEqn.relax(); | ||
| + | betaEqn.solve(); | ||
| + | |||
| + | alpha = 0.5*(scalar(1) + sqr(scalar(1) - beta) | ||
| + | - sqr(scalar(1) - alpha)); | ||
| + | */ | ||
| + | |||
| + | // The continuous phase volume fraction is computed as the complement to one of | ||
| + | // the dispersed phase fraction. | ||
| + | beta = scalar(1) - alpha; | ||
| + | } | ||
| + | |||
| + | Info<< "Dispersed phase volume fraction = " | ||
| + | << alpha.weightedAverage(mesh.V()).value() | ||
| + | << " Min(alpha) = " << min(alpha).value() | ||
| + | << " Max(alpha) = " << max(alpha).value() | ||
| + | << endl; | ||
| + | } | ||
| + | |||
| + | rho = alpha*rhoa + beta*rhob; | ||
| + | </cpp> | ||
==== Implementation of the pressure equation ==== | ==== Implementation of the pressure equation ==== | ||
