[chriss85]

I would like to add some source terms to the density, impulse and energy equations in rhoCentralFoam. Unfortunately I don't have a complete understanding of the algorithm yet, so it's somewhat difficult for me.

Why are two equations solved for impulse and energy respectively (viscid case), i.e. for rhoU and for U? Same goes for rhoE and e...

for U:

solve(fvm::ddt(rhoU) + fvc::div(phiUp));

        U.dimensionedInternalField() =
            rhoU.dimensionedInternalField()
           /rho.dimensionedInternalField();
        U.correctBoundaryConditions();
        rhoU.boundaryField() = rho.boundaryField()*U.boundaryField();

        volScalarField rhoBydt(rho/runTime.deltaT());

        if (!inviscid)
        {
            solve
            (
                fvm::ddt(rho, U) - fvc::ddt(rho, U)
              - fvm::laplacian(muEff, U)
              - fvc::div(tauMC)
            );
            rhoU = rho*U;
        }

and for e:

// --- Solve energy
        surfaceScalarField sigmaDotU
        (
            "sigmaDotU",
            (
                fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
              + (mesh.Sf() & fvc::interpolate(tauMC))
            )
            & (a_pos*U_pos + a_neg*U_neg)
        );

        solve
        (
            fvm::ddt(rhoE)
          + fvc::div(phiEp)
          - fvc::div(sigmaDotU)
        );

        e = rhoE/rho - 0.5*magSqr(U);
        e.correctBoundaryConditions();
        thermo.correct();
        rhoE.boundaryField() =
            rho.boundaryField()*
            (
                e.boundaryField() + 0.5*magSqr(U.boundaryField())
            );

        if (!inviscid)
        {
            solve
            (
                fvm::ddt(rho, e) - fvc::ddt(rho, e)
              - fvm::laplacian(turbulence->alphaEff(), e)
            );
            thermo.correct();
            rhoE = rho*(e + 0.5*magSqr(U));
        }

To which of these equations can I add my source terms?
I would think to the latter ones, as the turbulence is added there, however, some terms also appear in the first equations.
Do I need to discretize them using a central(upwind) scheme somehow?

By the way, for anyone interested, I think this might be a good read about the pressure and density-based solvers in OF:
http://www.cimec.org.ar/ojs/index.ph...File/4231/4157

[qutadah.r]

I think this solver uses flux terms to calculate the field using some kind of flux splitting scheme, and therefore it needs the equations in the conservative form, which have d(rho*U)/dt = RHS and d(rho*E)/dt = RHS for both momentum and energy equations, therefore you find the rho*U and rho*E recalculated again using the appropriate boundary conditions.