[trailer]

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:




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:

/*--------------------------------*- 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

            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?

[Tobermory]

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.

[trailer]

Originally Posted by Tobermory

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:

 volVectorField gradT(fvc::grad(T));

Has the following distribution (mathematically should be ):
[https://i.ibb.co/HzBhJZZ/grad-Image.png]

with the utility:

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?

[Tobermory]

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]

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:

    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)

[Tobermory]

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.

[Tobermory]

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