[jameswilson620]

All,

I am testing div() schemes for the transport of temperature. I am testing the solvers ability to maintain the resolution of a sharp interface (temperature front) advected by a divergence free velocity field in 1D. Upwind differencing breeds stability but introduces inconsistent diffusion terms due to the discretization; this drives my search for a more accurate scheme and also one that is robust.

I have tested variations of scalar transport equations without using interface compression such as mules in the VOF implementation in interFoam and have witnessed success using: div(phi,scalar) Gauss SuperBee.

So now I consider my new case while systematically testing each term of the energy equation in a modified two phase solver porousInterFoam. My current test is the convection term since it results in more difficulties than any other.

Assumptions:
1. 1D in x
2. Single phase (all phase 1)
3. convection dominant (k1=k2=1e-5[W/m/K])
4. constant inlet velocity

So here is the code:

cp 		= cp1*alpha1 + cp2*(1-alpha1);
k 		= k1*alpha1 +k2*(1-alpha1);
rho  	= rho1*alpha1 +rho2*(1-alpha1);

kEff = porosity*(alpha1*k1 + (1-alpha1)*k2) + (1-porosity)*kPorous;

rhoCpEff = porosity*(alpha1*rho1*cp1 + (1-alpha1)*rho2*cp2) + (1-porosity)*rhoPorous*cpPorous; 

rhoCpEff.oldTime() = porosity*(alpha1.oldTime()*rho1*cp1 + (1-alpha1.oldTime())*rho2*cp2) + (1-porosity)*rhoPorous*cpPorous; 
 
surfaceScalarField alphaPhi = (rhoPhi-phi*rho2)/(rho1-rho2+rhoEps);

surfaceScalarField rhoCpPhi = alphaPhi*(rho1*cp1-rho2*cp2)+phi*rho2*cp2;

fvScalarMatrix TEqn 
( 
  fvm::ddt(rhoCpEff,T) // Transient term Mean
  + fvm::div(rhoCpPhi,T) // Convective term fluid
  - fvm::laplacian(kEff, T) // Diffusion term Mean
//  - fvm::Sp(sensibleCorrectSource,T)
//  - fvm::Sp(hfgSource,T) //Split heat of vaporization to improve Diag.Dom.  
  ==
//  - hfgSource*Tsat //Split heat of vaporization to improve Diag.Dom.
);

the product rhoCp is constant and the inclusion of porosity (constant everywhere) is calculated and is considered in rhoCpEff and kEff. The convection of thermal energy in porous media only considers the interstitial fluid confined within the pores and therefore uses the product rhoCp of the fluid and NOT the effective rhoCpEff of the solid and fluid like the transient term and diffusion term.

Here is the fun part (all of the below utilize "div(corresponding argument) Gauss SuperBee" in system/fvSchemes):

fvm::div(rhoCpPhi,T) --> results in oscillations upstream of the sharp interface
fvm::rhoCp*div(phi,T) --> results in a bounded interface with expected numerical diffusion (this is perfectly acceptable but non-conservative)

rhoCp is a constant! div(rhoCpPhi,T) Gauss SuperBee vs div(phi,T) Gauss SuperBee differs only by a constant! Does this really matter in the SuperBee discretization scheme (or any other for that matter)?

I tested this also:

PHI = 2*phi;

rhoCp*fvm::div(PHI,T) --> using "div(PHI,T) Gauss SuperBee" results in oscillations upstream of the sharp interface
rhoCp*fvm::div(PHI,T) --> using "div(phi,T) Gauss SuperBee" results in oscillations upstream of the sharp interface

This simply confirms that the discretization of the div() operator broke because of the constant 2.

Then I tested an inlet velocity twice that of the previous cases and left phi unmodified but made a dummy variable PHI just to test passing the variable to another:

PHI = phi; // seems silly, i know

rhoCp*fvm::div(PHI,T) --> using "div(PHI,T) Gauss SuperBee" results in oscillations upstream of the sharp interface
rhoCp*fvm::div(PHI,T) --> using "div(phi,T) Gauss SuperBee" results in a bounded interface with expected numerical diffusion!!!!!

The system/fvschemes is very sensitive it appears, and in the previous case, sensitive to case! When manually modifying phi, i.e. PHI = 2*phi, this yields qualitatively accurate results but with oscillations upstream of the interface. When doing the equivalent (doubling the velocity) through the inlet BC, the results are smooth as is expected since the modification was very formal. My purpose in exploring this problem is I would like to use a porosity modified velocity, e.g. v = porosity*V where v is the darcian velocity and V is the interstitial pore velocity as this is the velocity that would advect a tracer particle through porous media or alternatively, an interface!

It seems like any modifications to phi, other than modifying the velocity properly through the BC's or something silly like PHI=phi, results in a broken system/fvschemes entry for div(something,T) Gauss SuperBee.

The figures attached are: Initial condition, Conservative divergence term (div(rhoCpPhi,T)) (oscillation), and non-conservative term (rhoCp*div(phi,T)) (smooth).