CFD Online Discussion Forums - Manipulating source fvc::grad(vsf1+vsf2)

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Manipulating source fvc::grad(vsf1+vsf2) (https://www.cfd-online.com/Forums/openfoam-programming-development/96017-manipulating-source-fvc-grad-vsf1-vsf2.html)

Manipulating source fvc::grad(vsf1+vsf2)

Hello,

I need to manipulate the source term (fvc::grad(vsf1+vsf2)) of the momentum equation.
The physical behaviour should be, that every identified internal face should be treated as a zero gradient condition:

Code:

int icounterinternalfaces = 0;

forAll(phi,phii)

{

    label P = mesh.faceOwner()[phii];

    label N = mesh.faceNeighbour()[phii];   

    if(vsf1[P]<(vsf2[N]-vsf2[P]))

    {

        

        if(vsf1[N]<vsf3[N]){//faces which should be temporarily treated with a zero gradient condition for this timestep

            ....

        }

    }



      

}

Here I thought about two different approaches:
I'm using the gauss gradient scheme
...==fvc::grad(vsf1+vsf2)
First way I:
Should I copy the gaussGrad.H/C and manipulate the code? But then I need to use something like ...==fvc::grad(vsf1,vsf2,vsf3) to be able to distinguish between normal internal faces and the ones which should have a zero gradient condition. So I need to create also a new gradScheme-class.

Second way:
I use the fixedInternalValueFvPatchField, create temporarily an internal patch consisting of the specified faces and treating them with a zero gradient condition. But I'm not sure, if this is a good way.

I would be happy about every comment on this. I tend to use the first way, or do you think, there is a more easy way like manipulating directly the outcome of the grad-calculation :confused: like:

Code:

volVectorField sourceTerm = fvc::grad(vsf1+vsf2);

forAll(phi,phii)

{

    label P = mesh.faceOwner()[phii];

    label N = mesh.faceNeighbour()[phii];   

    if(vsf1[P]<(vsf2[N]-vsf2[P]))

    {

        

        if(vsf1[N]<vsf3[N]){//manipulating source term before giving it to the right sight of equation...

            sourceTerm...... = .....

        }

    }



      

}

Many thanks in advance

Andy

Try something like this in your code.

Code:

volScalarField sumVf = vf1+vf2;

surfaceScalarField sumVff = fvc::interpolate(sumVf);

volVectorField gradVsf = fvc::grad(sumVff);



const unallocLabelList& owner = mesh.owner();

const unallocLabelList& neighbour = mesh.neighbour();

const vectorField& Sf = mesh.Sf();

const scalarField& V = mesh.V();



scalarField& sumVfi = sumVf.internalField();

scalarField& issf = sumVff.internalField();

vectorField& igGrad = gradVsf.internalField();



forAll(owner, facei)

{

    label P = owner[facei];

    label N = neighbour[facei];

    

    if (vf1[P] < (vf2[N]-vf2[P])) & (vf1[N] < vf3[N])

    {

        igGrad[P] += Sf[facei]*(sumVfi[P] - issf[facei]) / V[P];

        igGrad[N] += Sf[facei]*(sumVfi[N] - issf[facei]) / V[N];

    }

}

This is based on the Gauss grad scheme. For every selected face it removes the original face contribution and adds the zero-gradient contribution. I'm not guaranteeing it will work but it's probably the easiest approach.

what about the correctBoundaryConditions

hi cliffoi,

many thanks for your hint.
Yes, I use the gaussgrad scheme and in the meantime I had a deeper look in that source code.
As you propose, I should calculate the gradient over all faces and subtracting the values of the specific faces, where neighour and owner fullfill some condition. Of course dividing the value with the volume of that cell.

With
volVectorField gradVsf = fvc::grad(sumVff);
I use the first template-class in gaussgrad-scheme, putting in the surfaceScalarField of the sumVf into the fvc::grad() and I get back the gradient (volVectorField gradVsf).
But what about the
correctBoundaryConditions(vsf, gGrad);
:confused:

In the original form, the gaussGrad does first the loop with the surfaceScalarField, and then the correctBoundaryConditions(vsf,gGrad)
with looping over the patches (except the coupled ones).
Of course, I have in the first template-class the function:
gGrad.correctBoundaryConditions();
But do I need the
correctBoundaryConditions(vsf, gGrad);
correction too?
And if so, how can I perform this calculation/correction?

Many thanks in advance!

Best regards

Andy

correctBoundaryConditions(vsf, gGrad) modifies the value of the gradient on all normal (non-coupled) boundaries so that the surface-normal component of gradient matches the surface-normal of the boundary condition. This correction has no influence on the internal field values of the gradient so, unless the boundary values are important to you, you can safely ignore this. Conveniently this function is static so you can call it yourself to correct the surface-normal component afterwards.

Code:

volScalarField sumVf = vf1+vf2;

surfaceScalarField sumVff = fvc::interpolate(sumVf);

volVectorField gradVsf = fvc::grad(sumVff);



const unallocLabelList& owner = mesh.owner();

const unallocLabelList& neighbour = mesh.neighbour();

const vectorField& Sf = mesh.Sf();

const scalarField& V = mesh.V();



//- Apply correction to internalField

scalarField& sumVfi = sumVf.internalField();

scalarField& issf = sumVff.internalField();

vectorField& igGrad = gradVsf.internalField();



forAll(owner, facei)

{

    label P = owner[facei];

    label N = neighbour[facei];

    

    if (vf1[P] < (vf2[N]-vf2[P])) & (vf1[N] < vf3[N])

    {

        igGrad[P] += Sf[facei]*(sumVfi[P] - issf[facei]) / V[P];

        igGrad[N] += Sf[facei]*(sumVfi[N] - issf[facei]) / V[N];

    }

}



//- Now apply correction to all coupled boundaries

forAll(mesh.boundary(), patchi)

{

    const unallocLabelList& pFaceCells =

        mesh.boundary()[patchi].faceCells();



    const vectorField& pSf = mesh.Sf().boundaryField()[patchi];



    const fvsPatchField<Type>& pssf = sumVff.boundaryField()[patchi];



    if (sumVf.boundaryField()[patchi].coupled())

    {

        const scalarField vf1P = vf1.boundaryField()[patchi].patchInternalField();

        const scalarField vf1N = vf1.boundaryField()[patchi].patchNeighbourField();



        const scalarField vf2P = vf2.boundaryField()[patchi].patchInternalField();

        const scalarField vf2N = vf2.boundaryField()[patchi].patchNeighbourField();



//        const scalarField vf3P = vf3.boundaryField()[patchi].patchInternalField();

        const scalarField vf3N = vf3.boundaryField()[patchi].patchNeighbourField();



        forAll(mesh.boundary()[patchi], facei)

        {

            if (vf1P[facei] < (vf2N[facei]-vf2P[facei]))

                & (vf1N[facei] < vf3N[facei])

            {

                label celli = pFaceCells[facei];



                igGrad[celli]]

                    += pSf[facei]*(sumVfi[celli] - pssf[facei]) / V[celli];

            }

        }

    }

}



fv::gaussGrad<scalar>::correctBoundaryConditions(sumVf, gradVsf);

Again, I haven't checked that this works but the concept is there.

summarising your last post

Hello cliffoi,

thank you for your very fast reply.

Please correct me, if I summaries something wrong:
1) First we go through the internal field, calculating the gradients with:
volVectorField gradVsf = fvc::grad(sumVff);

Then we manipulate the internal field with the first loop. Ok.

In file gaussGrad.C, after dividing the cell gradient with the volume, then comes:

Code:

gGrad.correctBoundaryConditions();

2)
This, we have to manipulate again with your suggestion:

Code:

//- Now apply correction to all coupled boundaries

forAll(mesh.boundary(), patchi)

{

    const unallocLabelList& pFaceCells =

        mesh.boundary()[patchi].faceCells();



    const vectorField& pSf = mesh.Sf().boundaryField()[patchi];



    const fvsPatchField<Type>& pssf = sumVff.boundaryField()[patchi];



    if (sumVf.boundaryField()[patchi].coupled())

    {

        const scalarField vf1P = vf1.boundaryField()[patchi].patchInternalField();

        const scalarField vf1N = vf1.boundaryField()[patchi].patchNeighbourField();



        const scalarField vf2P = vf2.boundaryField()[patchi].patchInternalField();

        const scalarField vf2N = vf2.boundaryField()[patchi].patchNeighbourField();



//        const scalarField vf3P = vf3.boundaryField()[patchi].patchInternalField();

        const scalarField vf3N = vf3.boundaryField()[patchi].patchNeighbourField();



        forAll(mesh.boundary()[patchi], facei)

        {

            if (vf1P[facei] < (vf2N[facei]-vf2P[facei]))

                & (vf1N[facei] < vf3N[facei])

            {

                label celli = pFaceCells[facei];



                igGrad[celli]]

                    += pSf[facei]*(sumVfi[celli] - pssf[facei]) / V[celli];

            }

        }

    }

}

3) Then we have to do finally: correctBoundaryConditions(vsf, gGrad);
Which is similar to your call:

Code:

fv::gaussGrad<scalar>::correctBoundaryConditions(sumVf, gradVsf);

Is it correct?

Thank you very much. You are a great help!

Best regards

Andy

That's pretty much it. Let me know if it actually works...

Hello Cliffoi,

the implementation of the source worked and I could manipulate the source term of the momentum equation.
Here I did one correction to your suggestion:

Code:

...

igGrad[P] += Sf[facei]*(sumVfi[P] - issf[facei]) / V[P];

igGrad[N] += Sf[facei]*(sumVfi[N] + issf[facei]) / V[N];

...

with using + instead of - for the neighbour cell gradient value, hence the surface normal vector in gaussGrad is pointing in opposite direction (from the neighbour cell point of view) and, I think, must use then opposite sign.

Then I got almost the expected result in my testcase.

However, I manipulated the source term of the momentum equation before the piso-loop. So, I have a flux-field, which is not divergence free, but approximately satisfies momentum (because this grad term has the same function as "fvc::grad(p)" in my incompressible flow).

Now I have to do the same corrections in the piso-loop [having a correction (?) of the grad(p)/grad(vsf1+vsf2) values as I see in Ferziger & Peric] because I have fluctuations in the flux-values after the time-step ...

So first, I must understand the piso-loop.

I think, some new interesting questions will arise from my side soon. Sorry about that in advance...

Best regards,

Andy