[doug]

Hi all,

I am trying to write my own turbulence model and am currently studying the kEpsilon model in OF 1.4 linked with simpleFoam. In the simpleFoam.C source file, it says:

phi = fvc::interpolate(U) & mesh.Sf();

I'm not exactly sure what's happening here. My best guess is that it's interpolating the vector U from the cell centers to the cell faces and storing it as the transport property phi. (What exactly is the "& mesh.Sf()?)

In the kEpsilon.C source file, the dissipation equation is solved using phi_ as the item of transport. The divergence of phi_ is calculated as part of that equation. However, I thought you could only calculate the divergence of a vector. Is phi_ somehow a vector now instead of a scalar? If so, how was it set? I'm assuming that the phi defined in simpleFoam.C is the same as the phi_ in kEpsilon.C.

Since velocity is the item of transport, could the phi_ in the turbulence model be replaced by U_? This would be a true vector.

Thanks for your help.

Doug

[coops]

Hi Doug,

As I understand it fvc::interpolate(U) & mesh.Sf() is, as you say, the interpolation of the velocity from the cell centre to the cell face, the '&' is the dot product symbol, see programmer's guide, and meshSf() is the face area magnitudes.

So this is saying, interpolate the velocity field and take the dot product with the face area vector.

Not too sure about the rest,

Hope this helps,

Shaun

[gschaider]

@"However, I thought you could only calculate the divergence of a vector. Is phi_ somehow a vector now instead of a scalar?": According to Gauss's theorem (http://en.wikipedia.org/wiki/Divergence_theorem) the divergence in a volume can be calculated by summing the flux on the surfaces. phi is just that (the flux on the surface). So I guess (havn't looked at the code for that) div(phi) does just that: sum up the phis on the surface and pretend it's div(U)

[lam]

Hi Bernhard,

If you think (as Shaun) that div(phi) is something like div(U), do you think that I have well implemented:

nu*laplacian(U)=grad(div(U))/epsi
(I obtained this expression by considering:
div(U)=epsi*p, where epsi is very small
(epsi [0 -2 1 0 0 0 0] 1e-6))
and
nu*Laplacian(U)=grad(p)

by writting in Foam:

tmp<fvvectormatrix> UEqn(laplacian(nu,U));
solve(UEqn()==fvc::grad(fvc::div(phi)/epsi));

?

Lam

[isabel]

Originally Posted by gschaider

@"However, I thought you could only calculate the divergence of a vector. Is phi_ somehow a vector now instead of a scalar?": According to Gauss's theorem (http://en.wikipedia.org/wiki/Divergence_theorem) the divergence in a volume can be calculated by summing the flux on the surfaces. phi is just that (the flux on the surface). So I guess (havn't looked at the code for that) div(phi) does just that: sum up the phis on the surface and pretend it's div(U)

In OpenFOAM you can compute div(U) which is correct because U is a vector.

However, phi is a "surfaceScalarField" (not a vector). What div(phi) is? And why you can compute a divergence of an scalar? By definition, divergence is applied to a vector, for a scalar it would be a gradient.