[Hgholami]

Dear Foamer
The pUCoupledFoam solver coupled velocity and pressure as 4 quantities and solve them as coupled manner.
Calculating phi in this solver is some different from the segregate fluid solver such as icoFoam. In pUCoupled solver

phi = (fvc::interpolate(U) & mesh.Sf()) + pEqn.flux() + presSource;

but in icoFoam is

phi = (fvc::interpolate(U) & mesh.Sf());
adjustPhi(phi, U, p);
if (piso.finalNonOrthogonalIter())
{
phi -= pEqn.flux();
}

The problem is come when I want to add temperature equation as energy in solver. the TEqn is

fvScalarMatrix TEqn
(
fvm::ddt(T)
+fvm::div(phi,T)
-fvm::laplacian(DT,T)
);
TEqn.solve()

by locating the TEqn in pUCoupledFoam and running some cases, I found some diverge in temperature. I think the problem comes from phi in div(phi,T), but by change the location of TEqn in code the problem still exist.
Anybody have ideas?
Thanks