[jagan1mohan]

Hi FOAMers, I have solver for phi (a scalar) which is more-like scalar transport with extra terms.

Consider only fvm::ddt(phi)-fvm::laplacian(k,phi)=0 and in my tests, I write-out both transient and laplacian as fvc::div(k*fvc::grad(phi)) into all time folders. For square 2-D region, periodic in x with k = 0.001, after long time, transient term goes to 1e-38 but fvc::div(k*fvc::grad(phi)) in Paraview never goes below 1e-7. This did not change with mesh or time-step refinement.

Moving forward I initialize phiNOT in createFields as x-coordinate of cell and modify equation as fvm::ddt(phi)-fvm::laplacian(k,phi)-fvc::laplacian(k,phiNOT)=0. We know Nabla Square phiNOT=0 and including this should not affect solution.

But fvc::div(k*fvc::grad(phiNOT)) in Paraview is order of 0.01 and simulation does not reach steady-state due to this numerical error (I think). This error is especially high at boundaries.

How to address this situation as order of phi field is 0.0001? Any choice of fvSchemes or fvSolution can improve this? I can share more details to troubleshoot this behavior.

Thank you,

[Tobermory]

Okay, well the first thing is that your expression fvc::div(k*fvc::grad(phi)) is something rather different to fvm::laplacian(k,phi). To see this, you need to dive into the code. Deep breath.

fvm::laplacian(k,phi) calculates the face-centred gradients grad(phi), from the cell-centred phi values, and then calculates the divergence by multiplying these face gradients by the face-centre value of k (which is interpolated from the cell-centre values ... but since it's constant ignore this) and the corresponding face area and summing over the cell faces.

Now, what about fvc::div(k*fvc::grad(phi))? Well this starts with fvc::grad(phi). This returns the cell-centred gradient of phi based on the run-time selected gradient (see system/fvSchemes/gradSchemes) ... typically this is "Gauss linear" which means using linear interpolation to calculate the face-centred value of phi, and then using these face values to calculate the cell-centred gradient (this ignores corrections for non-orthogonality). Next, we calculate the divergence using this cell-centre gradient and k ... this means linearly interpolating k times the gradient to the face-centre, multiplying it by the face area and summing over the faces.

So, you see, your fvc expression introduces two sets of interpolation, which will give a different result on a stretched or non-orthogonal grid, and that's without considering the effects at the boundaries. If you were to use fvc::laplacian(k, phi), then this would use snGrad to calculate the face-centred gradient, and fv:dic to calculate the laplacian, and the result would be identical to the fvm expression.

Not sure whether that's the region for your fvc expression failing to drop to zero, but I suspect it is not helping.

Finally, just a stylistic thing - I would suggest avoiding calling your scalar phi, since that's openFoam terminology for the face mass flux ... your call.

[Tobermory]

In fact, even on a uniform mesh, the cell-centred gradient is different to the face-centred gradient. In 1D, for a cell N with neighbours W and E, the w-face and e-face values of phi are:



and so the gradient of phi for cell N is:

and for cell E is:


So the face gradient is:


whereas snGrad (and the calculation in fvmLaplacian) would simply have:

[jagan1mohan]

Hi, thank you for reply.

I'm solving scalar-transport-like equation which has terms containing field (known field) and (unknown field). Here, is scalar field (not to be confused with flux, ) and is solved for, is constant and . As we see here with , .
In my tests, results are dependent on inclusion and not-inclusion of fvc::laplacian(, ) into fvScalarMatrix.

Using blockMeshDict, for square 2-D region, periodic in , please find attached contours of in which same contour values are compared. I'm initializing in createFields.H and writing-out this Laplacian term for Paraview post-processing as fvc::laplacian(, ) which is not zero at four boundaries of this 2-D square.

What is going wrong in this implementation. Also, contours as shown in 'Without LaplacianPhiNot.png' match with analytical results.

Thank you in advance,