
[Sponsors] 
October 14, 2014, 09:34 
pimpleFoam with cyclic BCs

#1 
Senior Member
ArathoN
Join Date: Jul 2011
Posts: 137
Rep Power: 16 
After some times I retried to simulate a channel with cyclic BCs with the solver pimpleFoam and using RAS models (i want only to compute the mean velocity profile).
The problem is i can't get the pressure variable to converge it will always stay at around 0.40.6 and I really tried everything. First some info on the geometry and physical properties: Code:
H=0.025m Lx=5*H (i tried longer Lx but nothing changed) nu=1.5e5 Re(H,Uc)=15000 Re(H,Ubar)=13500 I=7% (from experimental data that i have) with Uc=velocity at centerline, Ub=Ubar and H = width of the channel. Now the initialization of teh turbulent flow are Code:
K=0.48 Code:
epsilon(l)=22 with l=10%H epsilon(nut/nu)=139 with nut/nu=10 Code:
Overall domain bounding box (0 0 0) (0.125 0.025 0.025) Mesh (nonempty, nonwedge) directions (1 1 0) Mesh (nonempty) directions (1 1 0) All edges aligned with or perpendicular to nonempty directions. Boundary openness (1.15504e18 7.86934e17 2.92744e16) OK. Max cell openness = 2.51663e16 OK. Max aspect ratio = 22.206 OK. Minimum face area = 7.0364e08. Maximum face area = 3.3667e05. Face area magnitudes OK. Min volume = 1.7591e09. Max volume = 4.20837e08. Total volume = 7.8125e05. Cell volumes OK. Mesh nonorthogonality Max: 0 average: 0 Nonorthogonality check OK. Face pyramids OK. Max skewness = 1.42118e06 OK. Coupled point location match (average 8.25619e18) OK. schemes and fvsolution parameters initial conditions yplus and celltocell expansion (or R) lenght of the channel Bcs Wall function or not Unfortunately this didn't help. Here the archive with all the files, I atatched also one graph of teh residual profile (done with pyFoam). I hope someone could figure out the problem. As a side note even with the LES simulation i couldn't get the right redisuals, the p and uy residuals were the problems. I used the same flow properties with a finer mesh. EDIT: Maybe i've found the error, if i choose in controldict to run the simulation with fixed deltat the courant number will explode after some steps. I don't know how to solve it, I tried changing the timeStep with no success. What velocity U and spacing do you consider when defining courant number? I usually try to get the lowest timestep so i use the highest velocity and the lowest spacing grade (I EDIT2: I tried also a simulation with no turbulence but again the residual floats around the unit. OT: I have 3 questions: 1 what relation do you use to define epsilon? 2 what do you put in the label pRefCell? I've chosen first 1001 but it gave me an error then i switched to 101, I know that it is teh identification on one cell of the domain where you initialize the pressure (so you'll not have an illconditioned problem) but it the choice of the cell important? 3 Reading teh code I noticed that the utility yPlusRAS only calls the calculation of yPlus done in the WF. Moreover studying the WF there two different definitions of the yPlus(utau,y,nu) and yStar(k). However i didn't find nothing RELIABLE about the comparison of the two definitions and on how they interact with each other (i know how to define them from the wall equations) These question are only my left doubts after going deep into the code and i couldn't answer them. Last edited by ArathoN; October 14, 2014 at 12:46. 

October 15, 2014, 04:35 

#2 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 27 
Hi.
I am pretty sure we can get this running! To your second question: pRef is not used for initialization as I understand it, but for setting a reference point in the pressure field for the complete simulation. Thus, the cell you choose is set to a fixed pressure (as you declare it with "pRefValue = 0"). This is needed for incompressible simulations, if there is no absolute pressure information at the boundaries. Anyway, it never changes the results of the simulation, because it is just an offset in the pressure field. You can choose any cell and any value. Ok, now to solve the problems. 1) The mesh you sent in the zipfile: What y+ do you expect here? I see you use RASmodel kepsilon, so since this is a highRe model I guess you have y+ >> ? 2) How can I create the mesh? I only import files from ICEM, so I don't know how to use the blockMeshDict. 3) Why do you use a transient solver? Why don't you use simpleFoam? 4) Is this 2D?
__________________
The skeleton ran out of shampoo in the shower. 

October 15, 2014, 07:17 

#3  
Senior Member
ArathoN
Join Date: Jul 2011
Posts: 137
Rep Power: 16 
Thank you for your help.
Quote:
Quote:
Code:
blockMesh Quote:
I tried one simulation with simpleFoam but i had the same results. yeah, I only want the mean velocity profile and not the perturbations fields. I think i'm doing something really stupid but can't figure out where, because even with the les simulation i had the same problem, the residuals of (p,Uy) were always too high around 0.51. 

October 15, 2014, 08:54 

#4  
New Member
Hans Barósz
Join Date: May 2014
Posts: 22
Rep Power: 12 
Quote:
BTW: your k, which is needed to solve the NSE, is build from the 3 velocity components. But in your case, there are only 2. When you have a look at DNS data or in those LES threads, you can easily see that the 3rd component has a huge impact. When you leave it, you will make a mistake. So my advice: go for 3D. 

October 15, 2014, 10:30 

#5  
Senior Member
ArathoN
Join Date: Jul 2011
Posts: 137
Rep Power: 16 
Quote:
Tbh the RAS simulation, given the simplicity of the case, will always give with the right condition and mesh a mean velocity profile similar to the experimental data. However the perturbation will be far away from the real data especially in the regions near the mixing layer because they are totally modeled with and assumption that in these region isn't so valid (hyp. boussinesq). About the averaging, the RAS equation are based on the mean value of the variable and its perturbation for every time step, so in my opinion the velocity field represent the mean part of that variable. Let me explain properly in the turbulence equation (eg. pEqn or Ueqn) you are using the mean value of the velocity (see every textbook about RANS eg. turbulent flow of Pope) so the output will be the mean value. But this can be wrong. 

October 15, 2014, 11:29 

#6 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 27 
This case will converge to steady state, so if you use a RAS model, 2d is absolutely sufficient. It will give exactly the same result as 3d with symmetry b.c. in the third dimension. You are right, arathon. Also, since it will converge to steady state, there is no need to "average the fluctuations for a long time" since there are no fluctuations. All fluctuations are captured by the RAS model.
About the turbulence model: You use pretty wrong boundary conditions. First you must decide whether you want to have an y+ of 1 or let's say 100. For the first case you need a lowRe model with lowRe boundary conditions. Let's say you go for the LaunderSharmaKE. For the second case you need a highRe model with highRe b.c. As I read in this forum the standard kepsilon is a bad choice since it is pretty unstable. Better use the realizableKE. Then decide on the solver. I think you want simpleFoam. You need to change the fvSolution file so it knows what to do when SIMPLE is started. What your fvSolution is also missing is the underrelaxation. You won't get simple stable without underrelaxation. I didn't get PIMPLE stable without it, too. What do I need to change in the blockMeshDict to get a high and a lowRe mesh? If you tell me, I will get you both cases running...
__________________
The skeleton ran out of shampoo in the shower. 

October 15, 2014, 11:46 

#7  
Senior Member
ArathoN
Join Date: Jul 2011
Posts: 137
Rep Power: 16 
Quote:
I've confirmed that in OF 2.3 simplefoam supports cyclic BCs so i'll switch to it. I know about the underrelaxation, reading the OF code of solver and turbulence models was really a good journey of discovering how the concept learn on uni lectures are applied. Last about the yplus, i evaluate the y of the first cell center from wall using the desired yplus (if you want i can write all the process). Then i calculate the spacing deltay of the first cell and using this utility i define the number of grid points and the grading. In block/hex you have the definition of the block then the grid point in (x,y,z) afterward the grading (in OF it's the total expansion ratio deltalast/deltafirst). 

October 15, 2014, 11:51 

#8  
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 27 
Quote:
__________________
The skeleton ran out of shampoo in the shower. 

October 15, 2014, 12:41 

#9  
Senior Member
ArathoN
Join Date: Jul 2011
Posts: 137
Rep Power: 16 
Quote:
Code:
Ny Ryb Ryt 55 10.37 0.0964 40 16.25 0.0615 30 23.92 0.0418 25 30.28 0.0330 Code:
Y+=30 Ny Ryb Ryt 7 1.116 0.8961 7 1 1.0000 Code:
Y+=100 Ny Ryb Ryt 2 1.232 0.8116883117 Code:
Ny=Nyb=Nyt. Code:
blockMesh EDIT: I tried using simpleFoam with kOmegaSST but omega after some timesteps will explode and now the same happens with pimpleFoam. Last edited by ArathoN; October 16, 2014 at 07:03. 

October 16, 2014, 11:47 

#10 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 27 
These channels are crap... just give me some more time. I don't know why this "just runs" in Fluent and I always have stability problems in OpenFoam like hell...
__________________
The skeleton ran out of shampoo in the shower. 

October 16, 2014, 12:37 

#11  
Senior Member
ArathoN
Join Date: Jul 2011
Posts: 137
Rep Power: 16 
Quote:
If i choose simplefoam Ubar and gradP (and residual of turbulence variable) will rise until the solver will crash, the same if you use pimpleFoam with ddtschemes as steadyState. But if i choose the ddtscheme as Euler the simulation will not crash. I'm really curious on why this happens, with steady state the Ubar and gradp will go aroud e10 at the first step. I tried to change schemes and bcs with no success, I even changed deltaT. I'm looking at the code: Code:
Info<< "Pressure gradient source: uncorrected Ubar = " << magUbarAve << ", pressure gradient = " << gradP << endl; Code:
dGradP_ = (mag(Ubar_)  magUbarAve)/rAUave; scalar gradP = gradP0_ + dGradP_; Code:
magUbarAve += (flowDir_ & U[cellI])*volCell; (for all the cells) magUbarAve /= V_; Code:
Pressure gradient source: uncorrected Ubar = 2.27933, pressure gradient = 56090 Code:
magUbarAve = Ubar IMO the probem lies on the instabilities caused by the absence of the temporal derivative term that is one of the contributor to the convergence and stability of the solver (idea made after i've read "Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows " by Hrvoje Jasak). For this reason in simplefoam you need the underelaxaion factors to mitigate the absence of the temporal derivate contribution. Then this maybe affect the rAU matrix, so that the instabilities will propagate to the variable dGradP, infact there is no correction introduced for such a situation. However i can't explain how the Ubar value is like that inthe first timestep. Edit: maybe the root of the problem: Code:
magUbarAve =/= Ubar Edit 2: changing URF doesn't change the value of magUbarAve. Last edited by ArathoN; October 16, 2014 at 13:50. 

October 16, 2014, 12:41 

#12  
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 27 
Quote:
__________________
The skeleton ran out of shampoo in the shower. 

October 16, 2014, 13:17 

#13  
Senior Member
ArathoN
Join Date: Jul 2011
Posts: 137
Rep Power: 16 
Quote:
I even though maybe the pressure predictor is the cause and i've set it as 1 (so pisoFoam) but the simulation is running normally. I changed schemes, BCs, turbulence model and mesh properties but the only thing that make the simulation crash is the temporal derivative scheme. From what i gathered in the previous post with steady solver it computes wrongly the Ubar averaged over the domain and this leads to the crash. Unfortunately i can't understand why because the evaluation of volCell and V() is the same with steady or unsteady solvers as for U[CellI]. 

October 16, 2014, 15:40 

#14 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 27 
You did something wrong with simplefoam. I will send you the lowre case tomorrow. I don't know if the high re case works because it just such few cells in y direction.
__________________
The skeleton ran out of shampoo in the shower. 

October 17, 2014, 03:50 

#15  
Senior Member
ArathoN
Join Date: Jul 2011
Posts: 137
Rep Power: 16 
Quote:
By the way how do you choose the grading of the grid after fixing yplus? I know that it should stay under 1.2 (cell.tocell expansion), the aspect ratio should better be under 20 (read it in some book advised in openfoamwiki) and for wall resolved we need at least 10 pints into the boundary layer. Are there any other consideration you account off? EDIT: i plotted p with paraFoam and i get such results, aren't they quite strange? The pressure should always decay, why the is a rise in p? In this case Lx=20H shouls i lower it to filter this. Last edited by ArathoN; October 17, 2014 at 04:52. 

October 17, 2014, 09:41 

#16 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 27 
ArathoN,
Unfortunately, I didn't manage to get it converged. I get a picture quite similar to the residual plot of your first post. Surprisingly, the results look good and match pretty good with Fluent's results anyway. Now, I find yvelocity values (absolute) of about +/ 1e14 in my domain. I don't understand why these yvelocities of nearly zero can produce such a high residual. And why the pressure residual remains that high. Anyway, for the LaunderSharmaKE model, you can go with these b.c. (lowRe mesh!!!): k: Code:
internalField uniform 1; boundaryField { "(bottomWalltopWall)" { type fixedValue; value uniform 1.0e12; } "(inlet_topoutlet_topinlet_botoutlet_bot)" { type cyclic; } frontAndBack { type empty; } } Code:
internalField uniform 0.1; boundaryField { "(bottomWalltopWall)" { type fixedValue; value uniform 1e8; } "(inlet_topoutlet_topinlet_botoutlet_bot)" { type cyclic; } frontAndBack { type empty; } } Code:
internalField uniform 1e3; boundaryField { "(bottomWalltopWall)" { type nutLowReWallFunction; value uniform 0; } "(inlet_topoutlet_topinlet_botoutlet_bot)" { type cyclic; } frontAndBack { type empty; } } fvSolution: Code:
solvers { "(ppFinal)" { solver GAMG; tolerance 1e99; relTol 1e5; smoother DICGaussSeidel; nPreSweeps 0; nPostSweeps 1; nFinestSweeps 2; scaleCorrection true; directSolveCoarsestLevel false; cacheAgglomeration on; agglomerator faceAreaPair; nCellsInCoarsestLevel 500; mergeLevels 1; maxIter 15; } "(UkepsilonomegaUFinalkFinalepsilonFinalomegaFinal)" { solver PBiCG; preconditioner DILU; tolerance 1e99; relTol 1.0e6; maxIter 15; }; } SIMPLE { nNonOrthogonalCorrectors 0; residualControl { p 1e10; U 1e10; k 1e10; epsilon 1e10; omega 1e10; } pRefCell 1; pRefValue 0; } relaxationFactors { fields { p 0.3; } equations { "(U)" 0.7; "(k)" 0.5; "(epsilon)" 0.5; "(omega)" 0.8; } } Code:
ddtSchemes { default steadyState; } gradSchemes { default Gauss linear; grad(sqrt(k)) faceMDLimited leastSquares 1; } divSchemes { div((nuEff*dev(T(grad(U))))) Gauss linear; div((nuEff*T(grad(U)))) Gauss linear; div(phi,U) bounded Gauss linearUpwind grad(linUpw); div(phi,k) bounded Gauss linearUpwind grad(linUpw); div(phi,epsilon) bounded Gauss linearUpwind grad(linUpw); } laplacianSchemes { default Gauss linear uncorrected; } interpolationSchemes { default linear; } snGradSchemes { default uncorrected; } fluxRequired { default no; p ; }
__________________
The skeleton ran out of shampoo in the shower. 

October 23, 2014, 12:05 

#17 
Senior Member
ArathoN
Join Date: Jul 2011
Posts: 137
Rep Power: 16 
Finally I was able to obtain a converged solution with kOmagaSST and the uprofile is in accordance with the experimental data. I tried the laudershamaKE and while it didn't converge, res(p) ~0.1 & res(uy) ~0.01, the u profile was good and quite similar to komegasst.
By the way I reduced the tolerances on fvsolution yours were too restrictive and I noticed that after e9 the solution will not change. Finally when I see the flow in parafoam I noticed a weird behavior of p/uy at the inlet for komegasst and at the upper wall for laundersharma , there are pressure/uy spikes for no reason typical of a periodic condition on paper. See the images maybe this is the cause of the residual profile and maybe the pressure gradient option is still a bit bugged. I couldn't upload the images with the Android app, here a link to imgur http://imgur.com/a/t9lJy 

October 24, 2014, 05:18 

#18 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 27 
ArathoN, you can find a lowRe version of the kOmegaSST model here (post #5):
http://www.cfdonline.com/Forums/ope...tml#post491414 I implemented the same version that is implemented in Fluent, when you check the box "lowReblabla" in Fluent. The regular SST model in openFoam does not include the damping functions near the wall. Thus, even if you set lowRe b.c. in openFoam, you will not get matching profiles necessarily.
__________________
The skeleton ran out of shampoo in the shower. 

October 24, 2014, 16:30 

#19 
Senior Member
ArathoN
Join Date: Jul 2011
Posts: 137
Rep Power: 16 
I'll try it and report, i noticed that the laundersharmaKE gives better evaluation of wallshearstress hence better utau however the uprofile Is a bit worser.
OT: can you share with me the pdf of menter i couldn't find it. 

November 2, 2014, 17:38 

#20  
Senior Member
ArathoN
Join Date: Jul 2011
Posts: 137
Rep Power: 16 
Quote:
I tried the lowre version and it improved a little bit the data, however this would work only for the channel in a backward facing step it worsened the solution. I confirm that the pressure is the main problem when using cyclic condition with the simpleFoam solver, if you try to visualize the pressure field in parafoam after simulation you'll see that it is nonsense and unphysical (random high and low pressure spots). I tried to use a different Reynolds number and the same thing happened the pressure and Uy residual are a mess however the Uprofile is quite good. 

Tags 
channel, cyclic, openfoam, pimplefoam 
Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
How to setup cyclic BCs in simpleFOAM  hhuang84  OpenFOAM Running, Solving & CFD  14  December 8, 2016 12:39 
Possible createPatch/createBaffles bug?  simpomann  OpenFOAM Bugs  2  July 15, 2014 08:07 
renumberMesh and cyclic BC's  Jonathan  OpenFOAM PreProcessing  3  November 23, 2013 06:25 
Adjusting velocity using cyclic bc's  gregdB  Main CFD Forum  0  January 25, 2012 10:08 
Cyclic BCs using createPatch in OF 1.6.x  SunnyPP  OpenFOAM  2  August 6, 2010 11:21 