
[Sponsors] 
July 18, 2021, 17:17 
Doubt of fvc::laplacian in laplacianFoam

#1 
Member
Join Date: Mar 2021
Posts: 38
Rep Power: 3 
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: 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 { "(leftrighttopbottom)" { 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; } } // ************************************************************************* // Code:
volScalarField fvcLaplacian ( IOobject ( "fvcLaplacian", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), fvc::laplacian(DT, T) ); Would anyone be so kind as to explain why? 

July 19, 2021, 09:11 

#2  
Senior Member
Join Date: Apr 2020
Location: UK
Posts: 462
Rep Power: 11 
Quote:
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. 

July 19, 2021, 10:55 

#3  
Member
Join Date: Mar 2021
Posts: 38
Rep Power: 3 
Quote:
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)); [https://i.ibb.co/HzBhJZZ/gradImage.png] with the utility: Code:
postProcess func "grad(T)" [https://i.ibb.co/Mhnq3Bs/gradImage2.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? 

July 19, 2021, 11:10 

#4  
Senior Member
Join Date: Apr 2020
Location: UK
Posts: 462
Rep Power: 11 
Quote:
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. 

July 19, 2021, 16:45 

#5 
Member
Join Date: Mar 2021
Posts: 38
Rep Power: 3 
So,
In an orthogonal mesh (blockMesh) the results with two Neumann boundary conditions are Ok for the fvc::laplacian operator. I created a simple nonstructured 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 (nonempty/wedge) directions (1 1 0) Mesh has 2 solution (nonempty) directions (1 1 0) All edges aligned with or perpendicular to nonempty directions. Boundary openness (7.228014483e18 1.445602897e18 4.902400823e16) OK. Max cell openness = 1.664421853e16 OK. Max aspect ratio = 2.326112477 OK. Minimum face area = 0.0003310493266. Maximum face area = 0.005626613724. Face area magnitudes OK. Min volume = 3.310493266e05. Max volume = 0.0001065683383. Total volume = 0.1. Cell volumes OK. Mesh nonorthogonality Max: 29.33089879 average: 7.205332936 Nonorthogonality check OK. Face pyramids OK. Max skewness = 0.500547513 OK. Coupled point location match (average 0) OK. Mesh OK. 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/laplacianInterior.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) 

July 20, 2021, 12:53 

#6 
Senior Member
Join Date: Apr 2020
Location: UK
Posts: 462
Rep Power: 11 
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. 

July 21, 2021, 11:15 

#7 
Senior Member
Join Date: Apr 2020
Location: UK
Posts: 462
Rep Power: 11 
This post might be interesting for you to read:
"Gauss linear" gradient makes OpenFOAM zerothorder 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 

Thread Tools  Search this Thread 
Display Modes  


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 
Multiregion problem using laplacianFoam  zfaraday  OpenFOAM PreProcessing  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 