CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Verification & Validation

Doubt of fvc::laplacian in laplacianFoam

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

Like Tree2Likes
  • 1 Post By trailer
  • 1 Post By Tobermory

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   July 18, 2021, 17:17
Default Doubt of fvc::laplacian in laplacianFoam
  #1
Member
 
Join Date: Mar 2021
Posts: 38
Rep Power: 3
trailer is on a distinguished road
Hello to all,

I am new to CFD and OpenFOAM programming, and as part of my initial work I am taking a look into the OpenFOAM explicit operators (finite volume calculus aka fvc).

Using the laplacianFoam solver I am initializing a field to be:

T(x,y) = x^2


I am not solving the T equation, just computing the laplacian of the field.



Assuming DT to be unity, the result should be 2.


Here is the code I have for the T field:

Code:
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 1 0 0 0];


internalField  #codeStream
{
    codeInclude
    #{
        #include "fvCFD.H"
    #};

    codeOptions
    #{
        -I$(LIB_SRC)/finiteVolume/lnInclude \
        -I$(LIB_SRC)/meshTools/lnInclude
    #};

    codeLibs
    #{
        -lmeshTools \
        -lfiniteVolume
    #};

    code
    #{
        // Access internal mesh information
        const IOdictionary& d = static_cast<const IOdictionary&>(dict);
        const fvMesh& mesh = refCast<const fvMesh>(d.db());

        scalarField myInternalField(mesh.nCells(), 0.0);

        const volVectorField& cellCenter = mesh.C();

        forAll(myInternalField, cellI)
        {
            const scalar x = cellCenter[cellI].component(0);

            myInternalField[cellI] = x*x;
        }

        myInternalField.writeEntry("", os);

    #};

};

boundaryField
{
    "(left|right|top|bottom)"
    {
        type        codedFixedValue;
        value       uniform 1;

        name        Tboundary;  

        code
        #{

            const fvPatch& boundaryPatch = patch(); 

            const vectorField& Cf = boundaryPatch.Cf();

            scalarField& field = *this;

            forAll(Cf, faceI)
            { 
                const scalar x = Cf[faceI].x();
             
                const scalar solutionT = x*x;

                field[faceI] = solutionT;
            }

        #};
    } 

    frontAndBack
    {
        type            empty;
    }
}

// ************************************************************************* //
In a modified version of laplacianFoam I am asking to compute the fvc::laplacian
Code:
            volScalarField fvcLaplacian
            (
                IOobject
                (
                    "fvcLaplacian",
                    runTime.timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
                fvc::laplacian(DT, T)
            );
In a 64x64 mesh the field has a uniform value value of 2 except near the left and right patches.





Would anyone be so kind as to explain why?
Shibi likes this.
trailer is offline   Reply With Quote

Old   July 19, 2021, 09:11
Default
  #2
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 462
Rep Power: 11
Tobermory will become famous soon enough
Quote:
Would anyone be so kind as to explain why?
Well, the easy but probably not very helpful answer is: because of your boundary condition.

Let's think a little deeper - to calculate the finite volume of the laplacian, the code uses the Gauss theroem to convert the divergence of the gradient into a surface flux of the gradient ... for which it needs to calculate the gradient of the field at the cell faces. Not a problem for the internal cells, but what happens at the cells bordering the boundaries? You have specified a fixed value for T, and whilst this is not a bad guess (it results in a laplacian value of around 1.5 it seems), it is not exact because the field is varying as x^2, while the fv method assumes linear variation across a cell. To get an exact solution, try setting the gradient of the field at the boundary and not the value. That should help.
Tobermory is offline   Reply With Quote

Old   July 19, 2021, 10:55
Default
  #3
Member
 
Join Date: Mar 2021
Posts: 38
Rep Power: 3
trailer is on a distinguished road
Quote:
Originally Posted by Tobermory View Post
Well, the easy but probably not very helpful answer is: because of your boundary condition.

Let's think a little deeper - to calculate the finite volume of the laplacian, the code uses the Gauss theroem to convert the divergence of the gradient into a surface flux of the gradient ... for which it needs to calculate the gradient of the field at the cell faces. Not a problem for the internal cells, but what happens at the cells bordering the boundaries? You have specified a fixed value for T, and whilst this is not a bad guess (it results in a laplacian value of around 1.5 it seems), it is not exact because the field is varying as x^2, while the fv method assumes linear variation across a cell. To get an exact solution, try setting the gradient of the field at the boundary and not the value. That should help.
Hello Sir,


Thank you very much for the input.

Putting one Neumann boundary condition on the right side (codedMixedBoundary):
[https://i.ibb.co/brYmjwF/image1.png]


Putting a Neumann boundary on the left and right side.
[https://i.ibb.co/DD0xVvN/image2.png]


Having two Dirichlet and Neumann boundary conditions appears to give the correct laplacian value.

However, something else weird occurred.

The gradient computed by the solver:
Code:
 volVectorField gradT(fvc::grad(T));
Has the following distribution (mathematically should be 2x):
[https://i.ibb.co/HzBhJZZ/grad-Image.png]

with the utility:
Code:
postProcess -func "grad(T)"
It shows the correct distribution:
[https://i.ibb.co/Mhnq3Bs/grad-Image2.png]

Back to the main topic, is placing 4 Dirichlet boundaries going to degrade the accuracy of the fvc::laplacian operator? It is supposed to be this way?
trailer is offline   Reply With Quote

Old   July 19, 2021, 11:10
Default
  #4
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 462
Rep Power: 11
Tobermory will become famous soon enough
Quote:
Back to the main topic, is placing 4 Dirichlet boundaries going to degrade the accuracy of the fvc::laplacian operator? It is supposed to be this way?
Glad that that helped.

Regarding the above question - your initial solution was not "wrong", it was just the solution to a slightly different problem - one where the boundary condition was not in equilibrium with the smooth solution that you were seeking. Both solutions were "correct".

Thinking of a real application, eg specifying velocity on an inlet plane, what is the impact? Well, you may want to apply a fixed value, in which case according to the above example there might be a minor error in the laplacian/diffusion term ... but this error is typically negligible, firstly because diffusion is typically much smaller than advection at an inlet, and secondly because the velocity should not be varying quadratically away from the inlet, so the laplacian gradient estimate at the inlet boundary is likely to be much better.

Bottom line - think carefully about the physicality of your boundary conditions - they are not trivial, and I have found that they have been the root of most of my problems in OF! Good luck.
trailer likes this.
Tobermory is offline   Reply With Quote

Old   July 19, 2021, 16:45
Default
  #5
Member
 
Join Date: Mar 2021
Posts: 38
Rep Power: 3
trailer is on a distinguished road
So,

In an orthogonal mesh (blockMesh) the results with two Neumann boundary conditions are Ok for the fvc::laplacian operator.

I created a simple non-structured mesh with gmsh:
[https://i.ibb.co/N1ny4w3/mesh.png]

Mesh statitics:
Code:
    Overall domain bounding box (0 0 0) (1 1 0.1)
    Mesh has 2 geometric (non-empty/wedge) directions (1 1 0)
    Mesh has 2 solution (non-empty) directions (1 1 0)
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (-7.228014483e-18 1.445602897e-18 -4.902400823e-16) OK.
    Max cell openness = 1.664421853e-16 OK.
    Max aspect ratio = 2.326112477 OK.
    Minimum face area = 0.0003310493266. Maximum face area = 0.005626613724.  Face area magnitudes OK.
    Min volume = 3.310493266e-05. Max volume = 0.0001065683383.  Total volume = 0.1.  Cell volumes OK.
    Mesh non-orthogonality Max: 29.33089879 average: 7.205332936
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.500547513 OK.
    Coupled point location match (average 0) OK.

Mesh OK.
Using the same boundaries as before and:

gradSchemes: Gauss linear
laplacinaSchemes: Gauss linear Corrected
divSchemes: none
interpolationSchemes: linear
snGradSchemes: corrected

The laplacian gets:
[https://i.ibb.co/PcCWQC4/laplacian1.png]

If I remove the boundary cells:
[https://i.ibb.co/NNM61Vd/laplacian-Interior.png]

Is this supposed to behave like this?



The gradient with Gauss linear gets:
from the solver (laplacianFoam):
[https://i.ibb.co/cLJkRHw/grad1.png]


From postProcess -func "grad(T)"
[https://i.ibb.co/0Y5tZ5C/grad2.png]


Changing the gradient scheme to least squares:
laplacian:
[https://i.ibb.co/Lrg9z9X/gradLSQ.png]


laplacian no cells near the boundary:
[https://i.ibb.co/ZgzKcsH/gradLSQ2.png]


gradient from laplacianFoam
[https://i.ibb.co/s1XZ1bK/gradLSQ3.png]


gradient from postProcess -func "grad(T)"
[https://i.ibb.co/L6gSSyH/gradLSQ4.png]




(Note: is it possible to resize an image inside the post? This way I can avoid putting links and instead place the images)
trailer is offline   Reply With Quote

Old   July 20, 2021, 12:53
Default
  #6
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 462
Rep Power: 11
Tobermory will become famous soon enough
Funnily enough, I am just struggling with a "simple" case at the moment that I built using gmsh, and which is based on triangles (extruded in one layer, as an asymmetric wedge model). The stability is awful, even with fully stable settings on tghe snGrad, laplacian and grad schemes.

I have just found that extruding some orthogonal hex based cells at the inlet boundary (see attached) seems to calm things down a little - otherwise I see junk like you were showing in your latest pictures. I am not sure what the cure is ... other than to avoid triangles/tets in OF, which is a bit lame. Please let me know if you make progress!

By the way - to upload figures in this forum, just click on the paperclip at the top of the window where you compose your message (to the right of the "fonts" selection) - you should be able to select your files.
Attached Images
File Type: jpg mesh.jpg (159.6 KB, 17 views)
Tobermory is offline   Reply With Quote

Old   July 21, 2021, 11:15
Default
  #7
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 462
Rep Power: 11
Tobermory will become famous soon enough
This post might be interesting for you to read:
"Gauss linear" gradient makes OpenFOAM zeroth-order accurate on unstructured meshes

and I am just looking also at this one, where the author is suggesting some skewCorrected settings that I had not seen before:
Unstructured mesh made out of OF
Tobermory is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Doubt regarding MUSCL scheme implementation and associated limiters SangeethCFD Main CFD Forum 0 October 7, 2017 10:00
Add solar radiation to laplacianFoam aar_out OpenFOAM Programming & Development 2 April 30, 2017 16:12
laplacianFoam: problem with fvc::laplacian leroyv OpenFOAM Programming & Development 2 September 2, 2016 08:34
Multi-region problem using laplacianFoam zfaraday OpenFOAM Pre-Processing 8 April 11, 2015 06:54
LaplacianFoam: fvc::laplacian explodes : bug in Openfoam? hrushi.397 OpenFOAM Running, Solving & CFD 0 August 12, 2013 02:41


All times are GMT -4. The time now is 21:52.