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

Face pressure interpolation scheme based on linear interpolation

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

Like Tree1Likes
  • 1 Post By Shibi

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 17, 2021, 18:48
Default Face pressure interpolation scheme based on linear interpolation
  #1
Member
 
Join Date: Feb 2020
Posts: 90
Rep Power: 4
Shibi is on a distinguished road
Hello to all,


I would like to create a new interpolation scheme for the face pressure.

Traditionally this is estimated by a linear interpolation scheme, where:
p_f = (1-\lambda_f)p_P + \lambda_f p_N

where \lambda_f is an interpolation factor.

Can anyone point me to which class(classes) I should change in OpenFoam to implement something based on this linear interpolation scheme for the pressure?


Best Regards!
Shibi is offline   Reply With Quote

Old   June 18, 2021, 10:39
Default
  #2
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 469
Rep Power: 11
Tobermory will become famous soon enough
Typically pressure enters the transport equations through the laplacian term, and the discretisation of this is linear because that is consistent with the finite volume approach.

Using anything other than linear is pretty unorthodox, but I guess you could try setting the laplacianScheme for pressure to something other than linear ... i would strongly recommend rethinking this, though. Can I ask why you want to apply something other than linear?
Tobermory is offline   Reply With Quote

Old   June 18, 2021, 11:21
Default
  #3
Member
 
Join Date: Feb 2020
Posts: 90
Rep Power: 4
Shibi is on a distinguished road
Quote:
Originally Posted by Tobermory View Post
Typically pressure enters the transport equations through the laplacian term, and the discretisation of this is linear because that is consistent with the finite volume approach.

Using anything other than linear is pretty unorthodox, but I guess you could try setting the laplacianScheme for pressure to something other than linear ... i would strongly recommend rethinking this, though. Can I ask why you want to apply something other than linear?

Hi,



Sure, I want to define a density weighted interpolation scheme for the face pressure:


Where:



p_f = \frac{(1-\lambda_f) \rho_N p_P + \lambda_f \rho_P p_N}{ (1-\lambda_f)\rho_N + \lambda_f \rho_P}


This is for multi-phase flow where the linear interpolation should introduce considerable numerical errors if the pressure gradient at cell P and N differ greatly.
Shibi is offline   Reply With Quote

Old   June 19, 2021, 19:21
Default
  #4
Member
 
Join Date: Feb 2020
Posts: 90
Rep Power: 4
Shibi is on a distinguished road
So,


I guess a good starting point could be the linear interpolation scheme at:


Code:
/OpenFOAM/OpenFOAM-v2012/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linear

and change the weights() function.


Since this->mesh().surfaceInterpolation::weights() will get the linear interpolation weights I think I just need to access the density at cell P and N. But I am not sure on how to do this.



Any suggestions would be highly appreciated.
Shibi is offline   Reply With Quote

Old   June 20, 2021, 11:34
Default
  #5
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 469
Rep Power: 11
Tobermory will become famous soon enough
The cell values on either side of a given face are addressed using the lduAddressing functions lower() and upper(), with the face number being the index. Take a look at fvmDiv for an example of them in action, and good luck.
Tobermory is offline   Reply With Quote

Old   June 22, 2021, 14:57
Default
  #6
Member
 
Join Date: Feb 2020
Posts: 90
Rep Power: 4
Shibi is on a distinguished road
Hello to all,

Some new updates.

When you do fvc::grad(p), the solver will go though several files, namely: fvcGrad.C, gradScheme.C, linear.H , surfaceInterpolation.C, basicFvGeometryScheme.C (if you do not have the linear weights), surfaceInterpolationScheme.C.

Of special importance to this problem are the surfaceInterpolationScheme.C and linear.H files.

In linear.H the weights will be computed as:

Code:
    

forAll(owner, facei)    

 {        

 // Note: mag in the dot-product.         

// For all valid meshes, the non-orthogonality will be less than         

// 90 deg and the dot-product will be positive.  For invalid        

 // meshes (d & s <= 0), this will stabilise the calculation        

 // but the result will be poor.         

 
scalar SfdOwn = mag(Sf[facei] & (Cf[facei] - C[owner[facei]]));       

scalar SfdNei = mag(Sf[facei] & (C[neighbour[facei]] - Cf[facei])); 

w[facei] = SfdNei/(SfdOwn + SfdNei);     }
and in the surfaceInterpolationScheme.C the dotInterpolate function will be called for the interpolation procedure, where the values at the face will be calculated as:


Code:
    for (label fi=0; fi<P.size(); fi++)
     {        

          sfi[fi] = Sfi[fi] & (lambda[fi]*(vfi[P[fi]] - vfi[N[fi]]) + vfi[N[fi]]); 

    }
With:
Sfi[fi] being a vector of ones (for the time being I do not know what this operation exists),
vfi[P[fi]] the value of the variable at cell P (owner)
vfi[N[fi]] the value of the variable at cell N (neighbor)


The above expression, mathematically, reads:
p_f = \lambda \phi_P + (1-\lambda)\phi_N
\Leftrightarrow \lambda \phi_P +\phi_N - \lambda \phi_N
\Leftrightarrow \lambda (\phi_P - \phi_N) +\phi_N


For my case, to implement a density weighted interpolation scheme it seems that openFoam is "too" hard-coded into the linear interpolation (where a unique weighting parameter is defined).

I could of course rewrite the dotInterpolate function in the source code to accommodate my needs, but I would really prefer not to.

Would really appreciate some suggestions on this issue!


Edit: this is solved. Just had to redefine the correct template functions



Best Regards
Tibo99 likes this.

Last edited by Shibi; June 23, 2021 at 09:36.
Shibi is offline   Reply With Quote

Old   February 10, 2022, 13:53
Default
  #7
Senior Member
 
René Thibault
Join Date: Dec 2019
Location: Canada
Posts: 103
Rep Power: 5
Tibo99 is on a distinguished road
Hi Shibi!

Really glad that you were able to find a solution to your challenge.

I am trying to implement a new interpolation scheme also, but for the velocity U.

My friend and I tried a different approach then yours, but it doesn't work for some reason and I cant find the bug. It compile and everything without error, but the results are not affected at all.

Here is the link to the GitHub repository:
https://github.com/RTibo/OpenFOAM-Ne...ation-Scheme-U

So, I was thinking about changing the method by using a 'weighting' factor and then I found your post.

It basically the same thing that I'm trying to implement. The only thing that will differ probably is the 'weighting' factor equation (eq. 45, see image).

Since I needed the same equation in my previous code, the 'weight' equation factor is already coded, but I don't know where to implement this 'weight' factor code snippet in the file linear.H.

Since you said in your last post that you make it work by using the file linear.H as a starting point, which template functions you redefine exactly?


Thank you very much for your help and consideration.

Best Regards,

René

Link of the paper (Eq. 44 & 45): https://onlinelibrary.wiley.com/doi/10.1002/fld.2709
Attached Images
File Type: png IntU.png (42.1 KB, 18 views)
Tibo99 is offline   Reply With Quote

Old   February 13, 2022, 14:27
Default
  #8
Member
 
Join Date: Feb 2020
Posts: 90
Rep Power: 4
Shibi is on a distinguished road
If you only want to change the weights, write a new weights function.

See if this is useful: https://github.com/UnnamedMoose/Basi...discretisation
Shibi is offline   Reply With Quote

Old   February 14, 2022, 09:20
Default
  #9
Senior Member
 
René Thibault
Join Date: Dec 2019
Location: Canada
Posts: 103
Rep Power: 5
Tibo99 is on a distinguished road
Thanks for answering!

From what I look quickly, the scheme of the tutorial interpolate the actual weights from the center of cells to the face of the cells (horizontal and vertical).

The weights function that I need to code is depending of 'z' (see eq. 45 in my previous post), but it only apply at the horizontal face of the cells. That means, the weights of the vertical face of the cell is 1.

So, what I conclude is, I don't need to interpolate the weights from the center of the cells, but assign the weights of the horizontal face of the cell directly. Then, use this 'weights' as 'lambda' in the file 'surfaceInterpolationScheme.C'.

I don't know if I can make it work that way? Or is there a better approach?

What is your thoughts?

Thank you again for replying!

René

Last edited by Tibo99; February 14, 2022 at 16:39.
Tibo99 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
how to set periodic boundary conditions Ganesh FLUENT 15 November 18, 2020 07:09
pimpleFoam in OF1612 shows same time step twice in log file shang OpenFOAM Bugs 10 January 24, 2018 11:43
[blockMesh] Axisymmetrical mesh Rasmus Gjesing (Gjesing) OpenFOAM Meshing & Mesh Conversion 10 April 2, 2007 15:00
Explicit Pressure Based Scheme Apurva Main CFD Forum 4 June 3, 2004 21:01
Hydrostatic pressure in 2-phase flow modeling (long) DS & HB Main CFD Forum 0 January 8, 2000 16:00


All times are GMT -4. The time now is 14:00.