Slip-free and pressure corner boundary condition
OpenLB – Open Source Lattice Boltzmann Code › Forums › on OpenLB › General Topics › Slip-free and pressure corner boundary condition
- This topic has 17 replies, 4 voices, and was last updated 6 years, 10 months ago by jb.
-
AuthorPosts
-
January 10, 2018 at 3:16 pm #1958jbMember
Dear OpenLB team,
Thank you for all the work you have put in making this code available.
I have a question about boundary conditions, and I would be happy to hear your opinion.Let me first give a description of the setup.
I use an STL file for the geometry and define the boundary conditions as follows: https://ibb.co/kOaB2RI also made a figure of the material numbers designated to each voxel, as a verification that they are implemented correctly:https://ibb.co/ckvThR
However, I am unclear how to treat the corners and edges. For instance, right now I get a number of unwanted effects. This image shows the streamwise velocity in a plane through the middle: https://ibb.co/h55B2R
My main concern is about the top right corner, where the velocity seems to slow down. Is there any way to fix this? Am I using the right boundary conditions? Am I forgetting a special treatment of this edge/corner?
Another concern is the bottom left corner, where I suppose that the velocity accelerates locally because of the sudden step towards the bounceback boundary condition in the corner? Is there a better way to do this? Right now this seems to result in some instabilities at the inlet region.
I appreciate your help very much,
Thank you in advance.January 10, 2018 at 5:08 pm #2766albert.minkModeratorI like the way you prepared the data.
In https://ibb.co/kOaB2R (Front View) I would like to see the slip boundary, green bar, for the whole top wall. Right now there is one voxel at the ‘corners/edges’ that is a velocity or pressure boundary.
Change the code where function rename() is called and keep me informed about this artefacts on the boundary.Best wishes
Albert MinkJanuary 10, 2018 at 5:49 pm #2767jbMemberDear Albert,
Thank you very much for the reply.Good idea!
I just tried it, and this is the result: https://ibb.co/fYwtz6
(as a reference the material numbers: https://ibb.co/kkA8XR )It seems like this did not solve the problem. Actually it now also has a problem in the top left corner.
Do you have any other suggestions?Thanks,
JuliaanJanuary 10, 2018 at 6:35 pm #2768albert.minkModeratorHi Juliaan,
the slip boundary in the most recent OpenLB version is found to be incorrect. (The release candidate is already fixed.) So I suspect the slip boundary to cause the error, at least at the top. What happens for a simple bounce-back top wall?
Regards
Albert MinkJanuary 11, 2018 at 1:27 pm #2769jbMemberHi Albert,
Ok, that is good to know. Do you have an idea when the fixed slip boundary will become available? Is this fix something I can already implement?
(I should mention that when I used the slip boundary I manually set the normals.)The results with a bounce-back look correct to me: https://ibb.co/gfwgMm
I think that the effects in the left corners are a logical result of using two different prescribed velocities in those edges (inlet and no-slip). I guess I should either use a pressure driven inlet/outlet (With a simple try this seemed to work at first, although eventually the simulation crashed, I will look more into it.) or prescribe a Poisseuile velocity profile like in the provided example ‘Cylinder3D’.Thanks for your help,
JuliaanJanuary 11, 2018 at 4:20 pm #2770albert.minkModeratorHi Juliaan,
happy to see your results with bounce back. I do agree with you explanation of two different velocities.
Concerning the slip boundary. With a git diff localized the fix in src/boundary/boundaryPostProcessors3D.hh Let’s have a look. The developer version reads
January 11, 2018 at 4:21 pm #2773albert.minkModeratorYep. Strange code interface in forum.
Code:template<typename T, template<typename U> class Lattice>
SlipBoundaryProcessor3D<T,Lattice>::
SlipBoundaryProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int discreteNormalX, int discreteNormalY, int discreteNormalZ)
: x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_)
{
OLB_PRECONDITION(x0==x1 || y0==y1 || z0==z1);
reflectionPop[0] = 0;
for (int iPop = 1; iPop < Lattice<T>::q; iPop++) {
reflectionPop[iPop] = 0;
// iPop are the directions which pointing into the fluid, discreteNormal is pointing outwarts
if (Lattice<T>::c[iPop][0]*discreteNormalX + Lattice<T>::c[iPop][1]*discreteNormalY + Lattice<T>::c[iPop][2]*discreteNormalZ < 0) {
// std::cout << “—–” <<s td::endl;
int mirrorDirection0;
int mirrorDirection1;
int mirrorDirection2;
int mult = 2 / (discreteNormalX*discreteNormalX + discreteNormalY*discreteNormalY + discreteNormalZ*discreteNormalZ);mirrorDirection0 = (Lattice<T>::c[iPop][0] – mult*(Lattice<T>::c[iPop][0]*discreteNormalX + Lattice<T>::c[iPop][1]*discreteNormalY + Lattice<T>::c[iPop][2]*discreteNormalZ)*discreteNormalX);
mirrorDirection1 = (Lattice<T>::c[iPop][1] – mult*(Lattice<T>::c[iPop][0]*discreteNormalX + Lattice<T>::c[iPop][1]*discreteNormalY + Lattice<T>::c[iPop][2]*discreteNormalZ)*discreteNormalY);
mirrorDirection2 = (Lattice<T>::c[iPop][2] – mult*(Lattice<T>::c[iPop][0]*discreteNormalX + Lattice<T>::c[iPop][1]*discreteNormalY + Lattice<T>::c[iPop][2]*discreteNormalZ)*discreteNormalZ);// bounce back for the case discreteNormalX = discreteNormalY = discreteNormalZ = 1, that is mult=0
if (mult == 0) {
mirrorDirection0 = -Lattice<T>::c[iPop][0];
mirrorDirection1 = -Lattice<T>::c[iPop][1];
mirrorDirection2 = -Lattice<T>::c[iPop][2];
}// computes mirror jPop
for (reflectionPop[iPop] = 1; reflectionPop[iPop] < Lattice<T>::q ; reflectionPop[iPop]++) {
if (Lattice<T>::c[reflectionPop[iPop]][0]==mirrorDirection0 && Lattice<T>::c[reflectionPop[iPop]][1]==mirrorDirection1 && Lattice<T>::c[reflectionPop[iPop]][2]==mirrorDirection2) {
break;
}
}
//std::cout <<iPop << ” to “<< jPop <<” for discreteNormal= “<< discreteNormalX << “/”<<discreteNormalY <<std::endl;
}
}
}January 24, 2018 at 1:14 pm #2771jbMemberHi Albert,
Thanks again.
Except for a small difference on the line “int mult = 2 / (d…” I didn’t find anything wrong. I actually think the problem is not with the slip-free bc, but potentially with the pressure outlet.Because I want to run high Reynolds number cases I have opted for a pressure driven + periodic domain approach and I have implemented a fringe region to set the inflow.
However, I have two questions:
1. when i try to run my simulations with MPI they fail over the following line of code (that I use to initialize the fringe region). For a simulation on 2 processors, the code fails when iCloc = 1.
Code:for (int iCloc = 0; iCloc < noOfCuboids; iCloc++) {
BlockGeometryStructure2D<T>& tmp = superGeometry.getBlockGeometry(iCloc);
dom_origin = tmp.getOrigin();
…
}Something goes wrong when I ask for the origin. Any suggestion what could be wrong here?
2.
I have tried to make the code as fast as possible by calculating as much as possible for the fringe region only once in the initialization phase. However, depending on the number of lattices included in the fringe region, the code has a slow down of 40-50%. I use the following function in the file superLattice2D.hh. Any suggestion on how I can reduce the computational cost for the fringe region?Code:template< typename T, template<typename U> class Lattice>
void SuperLattice2D<T,Lattice>::defineFringe( int* fringe_iCloc, int* fringe_iX, int* fringe_iY, T* fringe_weight, int fieldBeginsAt, int sizeOfField, T velocity, int fringe_N)
{
T fringeF[1];
T output[2];
for (int iR = 0; iR < fringe_N ; ++iR){
_extendedBlockLattices[ fringe_iCloc[iR] ].get( fringe_iX[iR] , fringe_iY[iR]).computeU(output);
fringeF[0] = fringe_weight[iR]*(output[0]-velocity);
_extendedBlockLattices[ fringe_iCloc[iR] ].get( fringe_iX[iR] , fringe_iY[iR]).defineExternalField (
fieldBeginsAt, sizeOfField, fringeF);
}Thanks a lot, your help is really appreciated,
JuliaanJanuary 24, 2018 at 1:56 pm #2772jbMembersmall correction/addition:
– The code works fine without MPI, also for multiple Cuboids. The problem only arises with MPI.
– The fringe region sets the inflow, I don’t add any other external force to the fluid voxels.
– I am working this out first in 2D. Later I will go to 3D.January 24, 2018 at 2:47 pm #2774albert.minkModeratorI know that OpenLB has a fringe zone implementation. I wonder whether you have implemented by your own, since I can not find the code you showed me in any file.
January 24, 2018 at 3:00 pm #2775jbMemberOh ok, stupid of me! I didn’t know!
Ok I found it on the developers guide.
Thanks!January 24, 2018 at 3:47 pm #2776Markus MohrhardParticipantHello Juliaan,
Quote:Quote from juliaan on January 24, 2018, 13:14
Hi Albert,Thanks again.
Except for a small difference on the line “int mult = 2 / (d…” I didn’t find anything wrong. I actually think the problem is not with the slip-free bc, but potentially with the pressure outlet.Because I want to run high Reynolds number cases I have opted for a pressure driven + periodic domain approach and I have implemented a fringe region to set the inflow.
However, I have two questions:
1. when i try to run my simulations with MPI they fail over the following line of code (that I use to initialize the fringe region). For a simulation on 2 processors, the code fails when iCloc = 1.
Code:for (int iCloc = 0; iCloc < noOfCuboids; iCloc++) {
BlockGeometryStructure2D<T>& tmp = superGeometry.getBlockGeometry(iCloc);
dom_origin = tmp.getOrigin();
…
}Something goes wrong when I ask for the origin. Any suggestion what could be wrong here?
You are not allowed to just access all the data when you use MPI. In MPI mode your data is distributed across multiple processes so you need to use a concept like:
Code:for (int iC=0; iC<this->_loadBalancer.size(); ++iC) {
BlockGeometryStructure2D<T>& tmp = superGeometry.getBlockGeometry(iC);
…
}If you execute this code now in each MPI process you will process each BlockLattice by the correct MPI process.
January 24, 2018 at 5:50 pm #2777jbMemberThanks for the help! My implementation works now! (although it definitely slows down the simulation)
Albert,
If I understand correctly, the Fringe2D function in OpenLB returns the coefficients to be used with a ‘linear velocity force model’. Is it then correct that the coefficients need to be used with SmagorinskyLinearVelocityForcedBGKdynamics ?Thanks for your help,
JuliaanJanuary 25, 2018 at 9:56 am #2778mathiasKeymasterYes, that is correct! Have a look at our upcomming spring school – it could be interesting. I can also make contact to our turbulence expert developer. In any of the two cases, please contact me directly.
Best
Mathias J. (Krause)January 25, 2018 at 5:10 pm #2779jbMemberCool, i got the fringe region to work. Using the Fringe2D function with SmagorinskyLinearVelocityForcedBGKdynamics is definitely the best option, no speed reduction!
-
AuthorPosts
- You must be logged in to reply to this topic.