[yeharav]

Hi all,

I am trying to calculate the individual terms of the convection terms in post-process.

The convection term is the vector

fvc::div(phi,U).

The convection term evaluated on the surface works fine

fvc::div(phi*fvc::interpolate(U))

.

However, when I use

fvc::grad(phi*fvc::interpolate(U))

and take components 0,1 and 2
the results are not related to anything. Infact even just component 0 is not related to anything.

Any ideas?

[HPE]

gradient of a vector is a tensor, not another vector.

[yeharav]

Hi,

I have solved this a while ago but I think its a good thing to post the solution here if somebody will ever need it.

So, as HPE mentioned the result is indeed a tensor. This is what we are trying to calculate.

That is we want a tensor whose [i.j] item /.


As I mentioned, just taking fvc::grad(phi*fvc::interpolate(U)) is erroneous.

There are several reasons for that:
1. phi is already multiplied by the surface area of the cell and therefore calling fvc::grad
multiplies it twice.
2. There are some correction for the boundary conditions.

Using the gauss term we try to calculate the surface integral of
/. So the question is now to evaluate these
fields on the cell face that is consistent with the way OF calculate the div(phi,U) term.

My solution is:
First calcualte the tensor terms on the face of each cell.

 surfaceTensorField interpolateUU    =   (phi*mesh.Sf()/mesh.magSf()))*fvc::interpolate(U);

Note that phi is scalar and

mesh.Sf()/mesh.magSf()

is the vector that is normal to the face direction (whose length is 1).

Now we need to correct the boundary fields.

forAll(mesh.boundary(), patchi)
{   
		auto& boundary_interpolateUU 	= interpolateU2.boundaryFieldRef()[patchi];

		const auto& boundary_phi	     = phi.boundaryField()[patchi];
		const auto& boundary_interpolateU    = interpolateU.boundaryField()[patchi];
		const auto& boundary_mesh 	     = mesh.boundary()[patchi];

		forAll(boundary_mesh, facei)
        	{
		    vector n_inside = -(boundary_mesh.Sf()[facei]/boundary_mesh.magSf()[facei]);
		    tensor tt = (boundary_phi[facei]*n_inside )*boundary_interpolateU[facei];
		    for (int i=0; i<9;i++) { 
				boundary_interpolateUU[facei].component(i) = tt.component(i);
		    }
	        }
}

Finally, we can calculate the value of the flux using th surface integral.

fvc::surfaceIntegrate(interpolateUU);