Implementing a new interpolation scheme

chriss85 · May 23, 2016, 08:31

I want to solve a Laplace equation fvc::laplacian(sigma, phi)=0 which contains a nonsteady jump in sigma at a coupled boundary (baffle would also be fine, I just need the surface to mark the jump). If I use a standard interpolation scheme then sigma*fvc::grad(phi) won't be constant on both sides of the boundary because the interpolation scheme doesn't consider that the physically correct solution has a jump in the gradient at the boundary.

In the end I think this comes down to implementing a new surfaceInterpolation scheme that extrapolates the field value at the boundary face by considering only the cell centre value and the face values from other faces of the cell to the face at the boundary. I have some ideas how to calculate the face value, given the cell centre and cell face values but I'm not sure how to implement this. However, using extrapolation from both sides will yield a different solution depending on the side from which the face value is extrapolated, unless the physically correct solution for phi is obtained. This means that the framework which is commonly used in the surfaceInterpolation classes is not going to work for me as far as I can tell. Would it be better to implement a gradient operator instead?

I'm not really sure what I need to do to make OpenFOAM find the physically correct solution. Theoretically, one could use very small cells and smooth out the jump for a solution with minimal error but this is not possible because it will negatively affect the CFL number.

I would appreciate some input on this, be it ideas, literature references,...

chriss85 · May 25, 2016, 11:07

I made some progress. Most of my previous post was based on a misunderstanding of how the finite volume method works in regards to the conserved quantity. When calculating fvm::laplacian(sigma,phi)=0, the conserved quantity is not actually sigma*fvc::grad(phi) but rather the face flux returned by phiEqn.flux(). Looking in the gaussLaplacianScheme.C file one can find out that in the case of orthogonal meshes it is something like fvc::interpolate(sigma)*fvc::snGrad(phi). In the nonorthogonal case there are some explicit correctors employed which make things more difficult.

phiEqn.flux() returns a surfaceScalarField (as expected), so the quantity is defined on the surface instead of being located in the cell. One can use fvc::reconstruct(phiEqn.flux()) to get a volVectorField out of it. In the case of an orthogonal mesh this appears to be correct and I get a solution very close to the analytical one in a simple test case. However, doing the same on a nonorthogonal mesh with a corrected laplacian scheme results in an incorrect volVectorField containing atleast some cells that can contain some very wrong cells. Other parts on the other hand look ok, so I think this is very much mesh related and cells which are close to being orthogonal work fine.

Any idea about this? How to get it to work on nonorthogonal meshes?

chriss85 · May 30, 2016, 07:19

So I have finally solved this problem now.

The problem was caused by calculating the electrical current as j = - sigma * fvc::grad(phi). In the Finite Volume Method the conserved quantity is not the cell centered gradient but rather the flux across cell boundaries. This means that the current density should be calculated as j = - fvc::reconstruct(phiEqn.flux()). Same goes for the electric field, which is E = - fvc::reconstruct(phiEqn.flux() / fvc::interpolate(sigma)). This method yielded the best results for me so far and doesn't have issues at larger steps in the conductivity.

When you calculate the current density this way, you can also directly use the surfaceScalarField from the flux to evaluate the current through a patch.

Yujuan · July 23, 2020, 04:48

Hi, I think I meet similar problem as you. I don't totally understand the solution that you mentioned. I want to know how you exactly do when solving equations fvc::laplacian(gamma,phi). Can you explain in more detail? Thanks a lot!

Quote:

Originally Posted by chriss85

So I have finally solved this problem now.

The problem was caused by calculating the electrical current as j = - sigma * fvc::grad(phi). In the Finite Volume Method the conserved quantity is not the cell centered gradient but rather the flux across cell boundaries. This means that the current density should be calculated as j = - fvc::reconstruct(phiEqn.flux()). Same goes for the electric field, which is E = - fvc::reconstruct(phiEqn.flux() / fvc::interpolate(sigma)). This method yielded the best results for me so far and doesn't have issues at larger steps in the conductivity.

When you calculate the current density this way, you can also directly use the surfaceScalarField from the flux to evaluate the current through a patch.

chriss85 · July 23, 2020, 04:54

The equation is solved with fvm::laplacian(gamma, phi) = 0 (or some source term), it's just about how you calculate the grad(phi) or gamma * grad(phi) so that the conserved quantity, the flux across cell boundaries, is preserved.
For this you use eqn.flux() to get the surfaceScalarField with the face fluxes. You can then use the reconstruct() function to interpolate this to a volScalarField. In my tests I got the best results with this method compared to simply doing fvc::grad().

You should probably also solve the equation iteratively a few times when you have non orthogonal cells in your mesh (non-orthogonal corrector loops).

Yujuan · July 24, 2020, 13:29

Thanks a lot for your reply. So your procedure is that you solve equation fvm::laplacian(gamma,phi)=0, and then use eqn.flux() to get surfaceScalarField of phi, and finally use reconstruct to get volScalarField of phi? I'm still not clear how you change the way to calculate grad(phi). In my opinion, since the equation is constructed using fvm::laplacian(gamma,phi) = 0, then the calculation of grad(phi) and gamma*grad(phi) are already be done in fvm::laplacian class, any only when you use fvm::div(gamma(grad(phi)), you can change how you calculate grad(phi).

Quote:

Originally Posted by chriss85

The equation is solved with fvm::laplacian(gamma, phi) = 0 (or some source term), it's just about how you calculate the grad(phi) or gamma * grad(phi) so that the conserved quantity, the flux across cell boundaries, is preserved.
For this you use eqn.flux() to get the surfaceScalarField with the face fluxes. You can then use the reconstruct() function to interpolate this to a volScalarField. In my tests I got the best results with this method compared to simply doing fvc::grad().

You should probably also solve the equation iteratively a few times when you have non orthogonal cells in your mesh (non-orthogonal corrector loops).

May 23, 2016, 08:31	Implementing a new interpolation scheme	#1
chriss85 Senior Member Join Date: Oct 2013 Posts: 397 Rep Power: 19	I want to solve a Laplace equation fvc::laplacian(sigma, phi)=0 which contains a nonsteady jump in sigma at a coupled boundary (baffle would also be fine, I just need the surface to mark the jump). If I use a standard interpolation scheme then sigma*fvc::grad(phi) won't be constant on both sides of the boundary because the interpolation scheme doesn't consider that the physically correct solution has a jump in the gradient at the boundary. In the end I think this comes down to implementing a new surfaceInterpolation scheme that extrapolates the field value at the boundary face by considering only the cell centre value and the face values from other faces of the cell to the face at the boundary. I have some ideas how to calculate the face value, given the cell centre and cell face values but I'm not sure how to implement this. However, using extrapolation from both sides will yield a different solution depending on the side from which the face value is extrapolated, unless the physically correct solution for phi is obtained. This means that the framework which is commonly used in the surfaceInterpolation classes is not going to work for me as far as I can tell. Would it be better to implement a gradient operator instead? I'm not really sure what I need to do to make OpenFOAM find the physically correct solution. Theoretically, one could use very small cells and smooth out the jump for a solution with minimal error but this is not possible because it will negatively affect the CFL number. I would appreciate some input on this, be it ideas, literature references,...

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Constructing Surface Interpolation Scheme from Divergence Scheme Information	ngj	OpenFOAM Programming & Development	2	June 17, 2022 10:12
Flux Limiter and Interpolation Scheme	pablitobass	OpenFOAM Programming & Development	10	December 3, 2016 14:15
Implementing other interpolation method for mesh to mesh interpolation	Jacks	OpenFOAM Programming & Development	2	March 21, 2014 07:13
Use of upwind scheme for interpolation of u/v	quarkz	Main CFD Forum	6	August 30, 2011 05:10
On the harmonic interpolation scheme for the diffusivity coefficient	vkrastev	OpenFOAM	5	April 11, 2011 10:13

May 25, 2016, 11:07		#2
chriss85 Senior Member Join Date: Oct 2013 Posts: 397 Rep Power: 19	I made some progress. Most of my previous post was based on a misunderstanding of how the finite volume method works in regards to the conserved quantity. When calculating fvm::laplacian(sigma,phi)=0, the conserved quantity is not actually sigmafvc::grad(phi) but rather the face flux returned by phiEqn.flux(). Looking in the gaussLaplacianScheme.C file one can find out that in the case of orthogonal meshes it is something like fvc::interpolate(sigma)fvc::snGrad(phi). In the nonorthogonal case there are some explicit correctors employed which make things more difficult. phiEqn.flux() returns a surfaceScalarField (as expected), so the quantity is defined on the surface instead of being located in the cell. One can use fvc::reconstruct(phiEqn.flux()) to get a volVectorField out of it. In the case of an orthogonal mesh this appears to be correct and I get a solution very close to the analytical one in a simple test case. However, doing the same on a nonorthogonal mesh with a corrected laplacian scheme results in an incorrect volVectorField containing atleast some cells that can contain some very wrong cells. Other parts on the other hand look ok, so I think this is very much mesh related and cells which are close to being orthogonal work fine. Any idea about this? How to get it to work on nonorthogonal meshes?

May 30, 2016, 07:19		#3
chriss85 Senior Member Join Date: Oct 2013 Posts: 397 Rep Power: 19	So I have finally solved this problem now. The problem was caused by calculating the electrical current as j = - sigma * fvc::grad(phi). In the Finite Volume Method the conserved quantity is not the cell centered gradient but rather the flux across cell boundaries. This means that the current density should be calculated as j = - fvc::reconstruct(phiEqn.flux()). Same goes for the electric field, which is E = - fvc::reconstruct(phiEqn.flux() / fvc::interpolate(sigma)). This method yielded the best results for me so far and doesn't have issues at larger steps in the conductivity. When you calculate the current density this way, you can also directly use the surfaceScalarField from the flux to evaluate the current through a patch.

July 23, 2020, 04:54		#5
chriss85 Senior Member Join Date: Oct 2013 Posts: 397 Rep Power: 19	The equation is solved with fvm::laplacian(gamma, phi) = 0 (or some source term), it's just about how you calculate the grad(phi) or gamma * grad(phi) so that the conserved quantity, the flux across cell boundaries, is preserved. For this you use eqn.flux() to get the surfaceScalarField with the face fluxes. You can then use the reconstruct() function to interpolate this to a volScalarField. In my tests I got the best results with this method compared to simply doing fvc::grad(). You should probably also solve the equation iteratively a few times when you have non orthogonal cells in your mesh (non-orthogonal corrector loops).