[danielr]

I try to understand why the flux conservation in a simple diffusion problem seems to be inaccurate with my following implementation:

I have modified the scalarTransportFoam solver to accept a nonuniform diffusivity k instead of the constant diffusivity DT. The velocity field U is set zero. So I solve the diffusion problem
d/dt T=div k grad T
in a unit square (16x16 cells) with k=0.01 in the upper half and k=1 in the lower half.
Boundary conditions are T=0 at the top and fixedGradient at the bottom.

The attached figure shows a color plot of T (above) at steady state (integration time t=100) and a probe line (below) of k,T and the flux magnitude |F| along the vertical axis.

I calculate the diffusive flux (F=k grad T) in the modified scalarTransportFoam as F=k*fvc::grad(T).
The calculated F shows a wiggle, where k jumps from 1 to 0.01, which can't be correct in steady state.

What is the correct expression to evaluate the flux k*grad T?

Thanks
Daniel

[olwi]

Hi Daniel,

When you solve a laplacian with a finite volume method, the fluxes you use are evaluated on cell *faces*. The flux you computed from a gradient in cell centres (using an fvc operator) therefore does'n have anything to do with the fluxes you actually use in the solution. So, I'm not surprised it looks strange. The fluxes across cell faces are certainly conserved in the code.

Another thing: When you have such great variation in the diffusivity, you should be careful with the interpolation scheme you choose for k. It's physically more sensible to use "harmonic" rather than "linear" interpolation (especially if you have interfaces where k has step changes). Remember that a laplacian is evaluated as a sum of all face fluxes k_f*snGrad(T), where k_f is the interpolated value of k at each face.

When you select a laplacian scheme like "Gauss linear corrected" in fvSchemes, then "Gauss" stands for this method of surface integrals on faces (only option, I think), "linear" is the interpolation scheme for k in your case (here I would use "harmonic"), and "corrected" specifies that the surface-normal gradient of T is computed with full non-orthogonal corrections applied.

/Ola

[danielr]

Originally Posted by olwi

Hi Daniel,

When you solve a laplacian with a finite volume method, the fluxes you use are evaluated on cell *faces*. The flux you computed from a gradient in cell centres (using an fvc operator) therefore does'n have anything to do with the fluxes you actually use in the solution. So, I'm not surprised it looks strange. The fluxes across cell faces are certainly conserved in the code.

I was suspecting something like this. I didn't doubt that the scalarTransportFoam was correct, but was not sure how to get the correct flux in my case.

Remember that a laplacian is evaluated as a sum of all face fluxes k_f*snGrad(T), where k_f is the interpolated value of k at each face.

Evaluating k_f*snGrad(T) in my modified solver as you describe:
surfaceScalarField k_f=fvc::interpolate(k);
surfaceScalarField snGradT_f=fvc::snGrad(T);
F=fvc::interpolate(k)*fvc::snGrad(T);
returns the expected fluxes without wiggles. Thank you for this solution.

Can I use OpenFOAMs flux() functions to get the flux used in the laplacian(k,T) scheme immediately, like fvc::flux(phi,T) gets the convective flux?
This is obviously nonsense:
fvc::flux(fvc::laplacian(k, T));
but does a similar function exist, which always returns the used flux independent of the chosen discretisation scheme?

Thanks
Daniel

[olwi]

Hi Daniel,

I don't think there's a flux function for the diffusive fluxes, at least I didn't find one in the Doxygen djungle (I looked through the fvc namespace).

I note that the flux F you compute above (F=fvc::interpolate(k)*fvc::snGrad(T)) is a surfaceScalarField. Are you ok with that? If you really want something in cell centres to look at in the post-pro, then you can "reconstruct" a volVectorField (!) from your face field F like this:
Fc = fvc::reconstruct(F*mesh.magSf()),
where the surfaceScalarField mesh.magSf() is the face areas.

If you make sure that you use exactly the same fvSchemes for interpolation and snGrad in both the laplacian and the explicit calculation of F, then this Fc will be a very nice vectorfield, even on unstructured grids! Probably the best you can ever hope for. It's in some way "consistent" with your conservative face fluxes, although the notion of "conservation" only makes sense when you look at fluxes in and out of a control volume (like a finite volume cell!), and not vector values in points...

Happy hacking!

/Ola