abnormal temperature near interface when adding energy equation to interFoam
 

To be brief, I added energy equation to interFoam according to the modify tutorial in this thread:

http://www.cfd-online.com/Forums/ope...interfoam.html

The author mentioned two possible types of TEqn.H which is(in openfoam language):
fvScalarMatrix TEqn
(
fvm::ddt(rhoCp, T)
+ fvm::div(rhoPhiCpf, T)
- fvm::laplacian(kappaf, T)
);

and

fvScalarMatrix TEqn
(
fvm::ddt(rho, T)
+ fvm::div(rhoPhi, T)
- fvm::laplacian((kappa/Cp), T)
);

The author pointed out that with the first type of TEqn, temperature instability occurs near the interface when thermal property difference of two phase is large, on the other hand, with the second type TEqn, no instability will occur even thermal properties differs significantly.

I tested the behavior of interFoam with these two different kind of TEqn, exactly the same phenomena as the author pointed out occurred in my case. However I want go further and ask WHY this happens? As far as I am concerned, since the rhoPhiCpf*T term is also calculated with a Gauss upwind scheme, even the convection term becomes large because of Cp is multiplied, it should be limited to a reasonable level that at least will not cause instability. Unfortunately this obviously doesn't match with the test result.

Right now, the only way I could figure out to somehow find the problem is to look at the A matrix in A*T=b equation. By this I mean, since eventually the TEqn will be solved by openFoam's magic code and the LAST step must be solving the linear equation A*T=b(T is temperature field), by investigating A's value I could find what's happening near the interface and then distinguish the problem. However...I just don't know how to find out A, it seems that it hide very deeply because of the encapsulation structure of openfoam. Does anyone has any idea on how-to doing this mission? Or there is better way to figure out the problem?

I am very interesting in investigating the reason of this phenomena, please help and comment if you have any idea!

Have a look at gdbOF (http://openfoamwiki.net/index.php/Contrib_gdbOF).

Thanks for reply!
Actually I already knew the existence of gdbOF. Unfortunately I found the gdbOF manual was too brief and obviously there is a knowledge gap between my understanding and the lowest requirement in order to use gdbOF to debug...

Maybe I didn't give enough information?
I review the thread again
http://www.cfd-online.com/Forums/ope...interfoam.html

and "eberberovic" mentioned that face fluxes in the temperature equation should be consistent with the fluxes calculated for the momentum equation. I can't fully understand why it's important.

Anyway now I could some how concluded that:
1. If face fluxes in the temperature equation is NOT consistent with the fluxes calculated for the momentum equation, then totally temperature divergence happens.
2. If face fluxes in the temperature equation is consistent with the fluxes calculated for the momentum equation, but still the convection term is calculated in a form "rhophi*Cpf", a converged result could be obtained BUT instability occurs near interface.
3. Not calculate face flux in temperature equation at all. A fully converge and interface-stable result could be obtained.

But why is all that?
Please help!

ّI think you've define rhoPhiCpf in alphaEqn.H file? if yes, the results may be wrong because alphaEqn.H is included in alphaEqnSubCycle.H . In this condition you can get reasonable result by changing the value of nAlphaSubCycles in fvSolution (i.e nAlphaSubCycles=1) . in general you can act like definition of rhoPhi in alphaEqnSubCycle.H file.