[Gerben]

I am trying to solve a fluid in a solid container. The fluid is heated by a volumetric source term (an electromagnetic field, assumed to be constant in time but not homogeneous in space), causing the fluid to flow. The fluid should lose part of its heat to the surrounding solid. The temperature differences are fairly small, so I use the Boussinesq approximation. However, it turns out energy is not conserved.

My initial attempt was to use the following temperature equation for the fluid region:

fvScalarMatrix TEqn
(
    fvm::ddt(rho, T)
  + fvm::div(phi, T)
  - fvm::laplacian(alpha, T)
  ==
    fvc::Sp(source)
);
TEqn.solve();

and boundary condition compressible::turbulentTemperatureCoupledBaffleMix ed

Now, if I use the zeroGradient boundary condition between fluid and solid, so there is no heat exchange, I find that after 90 s, the thermal energy within the fluid region is 188 J, while integrating the provided power yields only 141 J.

If I remove the term fvm::div(phi, T) from the equation, the energy balance is correct, but evidently the temperature distribution does not take fluid into account anymore. The energy balance is also correct, if I split the equation like this:

fvScalarMatrix TEqn
(
    fvm::ddt(rho, T)
  + fvm::div(phi, T)
  - fvm::laplacian(alpha, T)
);
TEqn.solve();
fvScalarMatrix TsEqn
(
    fvm::ddt(rho, T)
  ==
    fvc::Sp(source)
);
TsEqn.solve();

However, if I understand it correctly, in this case T is updated only at the end of the time step, so after the second equation, meaning that the field T used in the second equation uses the same values as before the first equation; that is to say, it is the as if the TEqn equation is not there at all. If that is indeed correct, is there a way to update T in the mean time?

Moreover, when I then use the original boundary condition (compressible::turbulentTemperatureCoupledBaffleMi xed), there is heat flowing from the fluid region to the solid region, but the temperature in the fluid region is not affected. When I use the first TEqn (including the source term), the heat flow from the fluid to the solid region does affect the temperature in the fluid region, but there is still a nett energy production within the fluid.

Now, part of the problem seem to be in the implementation of the boundary condition. When I disable the source term and start with a homogeneously heated fluid (delta T is 5 K, compared to the solid), then after 10 s the heat lost by the fluid region is -529 J, whereas the heat gained by the solid is only 283 J. The amount gained by the solid appears to be correct. When I comment out the debug message in the turbulentTemperatureCoupledBaffleMixed source file, the heat transfer rate written to the log file is correct (and consistent with the energy gain in the solid region).

I am using OpenFOAM-2.2.x.

If anyone could give me a pointer how to solve this, that would be much appreciated.