[floquation]

Hello all,

As a post-processing task, I'm using swak4foam during runtime to output the convective and diffusive flux over a faceSet. Initially, this resulted in strange non-conservative results. The reason for this error was eventually traced to the way in which I asked swak4foam to compute those fluxes: although my calculation was correct in an exact mathematical sense, it was wrong in a discretized numerical sense. Different fvSchemes were used, which resulted in the erroneous results. - But that has been solved now.

The reason for posting this question is to test my understanding. Is it correct what I am saying?

---

To obtain the diffusive flux, OpenFoam uses the "laplacian(alphaEff,T)" fvScheme, which requires two schemes: (1) interpolation scheme for the coefficient and (2) snGrad scheme for the field. OpenFoam integrates the Laplacian operator in an exact sense and applies Gauss' theorem to obtain the equation that will be discretized.
Therefore, to obtain the diffusive flux that is correct in a discretized sense, one should do exactly the same, so:

createField_gradT
{
    type expressionField;
    fieldName AgradT;
    verbose true;
    autowrite false;
    expression "snGrad(T)*area()"; // units: K*m
}

DiffFlux0_0
{
    type swakExpression;
    valueType faceSet;
    setName plane0_0;
    verbose true;
    expression "sum(-nu/Pr*AgradT)/sum(area())";
    accumulations (
        average
    );

    autoInterpolate true;
    warnAutoInterpolate false;
}

This appears to work (i.e. it gives conservative results if summed with the convective flux). However, I'm not sure if this still works with a non-constant thermal diffusivity coefficient? Perhaps the coefficient must then be moved to the "createField_gradT" expression to ensure a correct interpolation to the faceSet?

Also, is "snGrad(T)*area()" the correct flux? That is, simply the surface normal gradient multiplied by the area? Or is OpenFoam doing something more complicated than that?

---

The convective term was the root cause of my initial problems. OpenFoam uses the "div(phi,T)" fvScheme, which is again integrated in an exact sense to arrive at a flux representation which is then discretized. I eventually came up with the following:

createField_fluxT
{
    type expressionField;
    fieldName fluxT;
    verbose true;
    autowrite false;
    expression "flux(phi,T)"; // units: K*m^3/s
}

ConvFlux0_0
{
    type swakExpression;
    valueType faceSet;
    setName plane0_0;
    verbose true;
    expression "sum(fluxT)/sum(area())";
    accumulations (
        average
    );

    autoInterpolate true;
    warnAutoInterpolate false;
}

// Additional entry for fvSchemes.divSchemes:
flux(phi_copyBySwak4Foam,T_copyBySwak4Foam)      Gauss vanLeer; // should be same as div(phi,T)

Which, again, gives a conservative result if summed with the diffusive flux.
Now, the reason one cannot use e.g. "phi*interpolate(T)", is that it will then not select the appropriate scheme (here: vanLeer).

Is the flux()-function exclusively meant for the convective flux? That is, can't it be used somehow for the diffusive flux as well?