CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Manipulating source fvc::grad(vsf1+vsf2)

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 1 Post By andyru
  • 1 Post By cliffoi

Reply
 
LinkBack Thread Tools Display Modes
Old   January 10, 2012, 11:31
Question Manipulating source fvc::grad(vsf1+vsf2)
  #1
New Member
 
Andreas Ruopp
Join Date: Aug 2009
Location: Stuttgart / Germany
Posts: 29
Rep Power: 7
andyru is on a distinguished road
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 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
mm.abdollahzadeh likes this.
andyru is offline   Reply With Quote

Old   January 26, 2012, 20:31
Default
  #2
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
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.
mm.abdollahzadeh likes this.
cliffoi is offline   Reply With Quote

Old   February 3, 2012, 05:43
Default what about the correctBoundaryConditions
  #3
New Member
 
Andreas Ruopp
Join Date: Aug 2009
Location: Stuttgart / Germany
Posts: 29
Rep Power: 7
andyru is on a distinguished road
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);


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
andyru is offline   Reply With Quote

Old   February 3, 2012, 12:46
Default
  #4
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
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.
cliffoi is offline   Reply With Quote

Old   February 7, 2012, 17:00
Thumbs up summarising your last post
  #5
New Member
 
Andreas Ruopp
Join Date: Aug 2009
Location: Stuttgart / Germany
Posts: 29
Rep Power: 7
andyru is on a distinguished road
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

Last edited by andyru; February 7, 2012 at 17:24. Reason: Better explanation
andyru is offline   Reply With Quote

Old   February 7, 2012, 17:30
Default
  #6
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
That's pretty much it. Let me know if it actually works...
cliffoi is offline   Reply With Quote

Old   February 9, 2012, 15:14
Default I made progess
  #7
New Member
 
Andreas Ruopp
Join Date: Aug 2009
Location: Stuttgart / Germany
Posts: 29
Rep Power: 7
andyru is on a distinguished road
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
andyru is offline   Reply With Quote

Reply

Tags
fvc, gauss, grad, manipulate, source

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
swak4foam building problem GGerber OpenFOAM Installation 54 April 24, 2015 16:02
pisoFoam compiling error with OF 1.7.1 on MAC OSX Greg Givogue OpenFOAM Programming & Development 3 March 4, 2011 18:18
OpenFOAM on MinGW crosscompiler hosted on Linux allenzhao OpenFOAM Installation 127 January 30, 2009 20:08
DxFoam reader update hjasak OpenFOAM Post-Processing 69 April 24, 2008 01:24
DecomposePar links against liblamso0 with OpenMPI jens_klostermann OpenFOAM Bugs 11 June 28, 2007 17:51


All times are GMT -4. The time now is 00:33.