[dvcauwe]

Hello Foamers, I have been breaking my head over this for the past 2 days, so I would really appreciate it if I could get some feedback on this.

I'm simulating flow through a pipe using a modified version of rhoPimpleFoam with cyclic boundaries at the in- and outlet of my domain while the flow is heated with a constant temperature gradient imposed on the tube wall. The problem however is that when I impose a gradient that is twice as big, the rate at which my gas is heating is only slightly higher and takes a considerable amount of time (~1s) to reach the double value of my initial deltaT/time with half the gradient imposed.

I understand that it takes a while to develop a temperature profile but still energy should be conserved, right? Performing additional subiterations offers a slight improvement but still doesn't balance out the energy equation. I thought a good way to monitor the heat going into the volume would be this function:

const fvPatchField<scalar>& heatFluxPatchField = T.boundaryField()[patchI];

// Print patch total heat transfer rate
const tmp<volScalarField> tkappa = thermo.kappa();
const scalarField& kappa = tkappa().boundaryField()[patchI];

scalar patchHeatRate = gSum
(
	kappa * heatFluxPatchField.snGrad() *
	heatFluxPatchField.patch().magSf()
);

, which shows me the desired total heat input but unfortunately this energy is not actually absorbed by the gas.

The energy equation is standard I'd say:

    volScalarField& he = thermo.he();

    fvScalarMatrix EEqn
    (
        fvm::ddt(rho, he) + fvm::div(phi, he)
      + fvc::ddt(rho, K) + fvc::div(phi, K)
      - fvm::laplacian(turbulence->alphaEff(), he)
     ==
        fvOptions(rho, he)
    );
    EEqn.relax();
    fvOptions.constrain(EEqn);
    EEqn.solve();
    fvOptions.correct(he);
    thermo.correct();

My only remaining guess at this moment is that setting a constant gradient for temperature simply does not guarantee that amount of energy going into the fluid during this time interval but that would imply setting an energy flux for...h? Then I'm not entirely sure how to use this field again in my thermo calculations...

Anyway, thanks to anyone that might be able to shed some light on this!

Best regards,
David

[dvcauwe]

Ok, I isolated the problem to the term "fvm::div(phi, he)", which according to my limited experience with OpenFOAM uses the face values of he. The lack of energy conservation then stems from my inlet patch automatically being converted to a fixedEnergy bc with values that are not being updated according to the changing temperature field.

By running the simulation 10 timesteps at a time and then restarting, the inlet patch is updated each time but obviously that's not an option so I need some way to change the enthalpy values on this patch. I tried all of the stuff below but nothing seems to work, please let me know if there's something else that might!

    label patchID = mesh.boundaryMesh().findPatchID("in");

    scalarField Tp = T.boundaryField()[patchID];
    const scalarField& pp = p.boundaryField()[patchID];
    scalarField heFix(thermo.he(pp, Tp, patchID));
    fixedEnergyFvPatchScalarField& hPatch = refCast<fixedEnergyFvPatchScalarField>(he.boundaryField()[patchID]);
    hPatch==heFix;
    forAll(hPatch, faceI) { hPatch[faceI] = heFix[faceI];
    hPatch.evaluate(Pstream::blocking);
    hPatch.updateCoeffs();

I don't feel like recompiling fixedEnergy since that would imply recompiling pretty much all thermo classes... For now I just manually (or using groovyBC) fix the temperature values at the inlet patch to those at the outlet, causing the enthalpy bc to become a gradientEnergy or mixedEnergy, which dó allow their values to be manipulated. Obviously this is a rather lousy approach so unless someone can think of a solution, I guess I'll file a bug report for the fixedEnergy bc.

Best regards,
David