
[Sponsors] 
December 12, 2016, 10:25 
Polyhedral mesh cyclic boundary problem

#1 
Member
Charlie Lloyd
Join Date: Feb 2016
Posts: 55
Rep Power: 5 
Hi all,
I have recently been using LES to simulate a turbulent flow at Re=180 on a hex mesh. The results have compared really well to DNS results for rms and mean profiles but I'm having issues when extending the problem to polyhedral meshes. The figures below show the profile (left is on the polymesh and right is the hexmesh). The results look reasonably similar in the ZY plane but there are large spikes in velocity in the XZ plane on the forced cyclic boundary. These spikes don't exist anywhere else so I think it is due to the fvOptions function that I am using. The mesh was generated using starccm+ and converted to OF using ccm26ToFoam and createPatch to match the two cyclic boundaries in Z and X. Has anyone else experienced something similar/ has any suggestions as to why this is occurring? My case file is set up almost exactly as the channel395 tutorial but I have changed the viscosity to 1/180 and the Ubar to 15.55556 to make Re_b = 2800, Re_tau = 180(ish) and the pressure gradient roughly equal to 1. I am also using the Lagrangian Dynamic LES model. Thanks in advance, Charlie FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object fvOptions; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // momentumSource { type meanVelocityForce; active yes; meanVelocityForceCoeffs { selectionMode all; fieldNames (U); Ubar (0 0 15.5555556); } } 

December 12, 2016, 18:35 

#2 
Senior Member
Santiago Lopez Castano
Join Date: Nov 2012
Posts: 254
Rep Power: 10 
Several could be the reasons why the 《red spots》 in the flow normal planes. The first is whether the two faces are topologically identical. Other if the volumes that share these faces are hexas or otherwise; if the facenormals of the internal faces are very slanted wrt the mean flow direction, there's a tendency to increase mass flux through the boundary face. Finally try to use the cubeVolRoot filter as the second filter, you need to be sure that the grid is not producing excessive oscillations over the bc's. What are you using as discretization for the convective term?
Sent from my GTI8190L using CFD Online Forum mobile app 

December 12, 2016, 19:05 

#3 
Member
Charlie Lloyd
Join Date: Feb 2016
Posts: 55
Rep Power: 5 
Thanks for the informative reply. I will provide more details when I am in the office tomorrow but I can answer some of those questions now:
the mesh faces are topologically identical  the mesh nodes/faces match on both boundaries to a higher tolerance. (Default is 10e3 I think?).  the volumes are likely to be mixed hex and poly but 95% of the mesh is poly. All the cells have a reasonable quality/skewness.  the face angles on the internal faces of the boundary is not something that I have looked at but are almost certainly not normal to the flow direction. I think this would be tough to specify in the meshing process? I will certainly investigate this tomorrow.  I will detail the fvSchemes tomorrow but I'm fairly sure the convective term is Gaussian linear. I'm not sure where the second filter is specified? Thanks for the reply Sent from my SMG920F using CFD Online Forum mobile app 

December 13, 2016, 05:26 

#4 
Senior Member
Santiago Lopez Castano
Join Date: Nov 2012
Posts: 254
Rep Power: 10 
The filter that you specify in LESProperties is, in fact, the second filter if dynamic procedures are used. Have in mind that FOAM is a library thought to be used for MILES, thus the primary filter is the grid itself
Sent from my GTI8190L using CFD Online Forum mobile app 

December 13, 2016, 05:26 

#5  
Member
Charlie Lloyd
Join Date: Feb 2016
Posts: 55
Rep Power: 5 
Quote:


December 13, 2016, 05:28 

#6 
Member
Charlie Lloyd
Join Date: Feb 2016
Posts: 55
Rep Power: 5 
Hi, I just got that message after I posted that reply. Should I swap vanDriest to cubeRootVol? Thanks for the help


December 13, 2016, 05:33 

#7 
Senior Member
Santiago Lopez Castano
Join Date: Nov 2012
Posts: 254
Rep Power: 10 
You need to have TDV/NVD algorithms for the convective term if you want to avoid oscillations due to the grid. Yes, delta defines the filter and 'delta' for the convolution
Sent from my GTI8190L using CFD Online Forum mobile app 

December 13, 2016, 08:46 

#8 
Member
Charlie Lloyd
Join Date: Feb 2016
Posts: 55
Rep Power: 5 
Ah okay, so should I specify Gauss filteredLinear as the div(phi,U) term?
I have now checked the mesh face orthogonality and it appears that most of the large spikes do occur at on cells that have faces not normal to the mean flow direction (see the attached photo). If this is the main reason why I am getting issues on the PBC I'm not sure I can fix it using meshing software! Eventually I will be simulating more complex geometries so polyhedral cells would be most useful for this. Is there a numerical trick I can use to correct these jumps? Perhaps if I specified a constant pressure gradient instead of the iterative one? Also, would you be able to provide a reference for the issues with nonorthogonality in PBC so I can read about it further? Or is it more of a speculative issue? Thanks in advance, Charlie 

December 13, 2016, 08:59 
References, boh!

#9 
Senior Member
Santiago Lopez Castano
Join Date: Nov 2012
Posts: 254
Rep Power: 10 
Ummm... On top of my head I have no specific reference in which they could mention this issue, in general good unstructured meshing is more of an "art". Anyway, I could recommend you to read "FiniteVolume CFD procedure and Adaptive Error Control Strategy for grids of arbitrary topology" of Muzaferija & Gosman as a starting point, then you could delve a bit on the references therein mentioned. As for how to "smoothout" those, you could either do a filtering of the resulting field (lowpass filter on postprocessing) or when running the simulations using dissipative (upwind) schemes for the convective terms; in the end it all depend on the degree of accuracy you wish to get, both options will "regularize" the solution in a way that you'll end up seeing the solution of a slightly different NS equation.


December 13, 2016, 13:01 

#10 
Member
Charlie Lloyd
Join Date: Feb 2016
Posts: 55
Rep Power: 5 
Thanks for the reference, I will be sure to go through it. I think adding dissipation into the numerics would reduce the accurately too much for this work so I think I will have to investigate the meshing further.
Would another solution be defining a constant body force to the Ueq in order to fix the mean pressure gradient? This would mean there would be no iterative forcing on the boundary so would this remove the boundary errors? 

December 13, 2016, 13:41 

#11 
Senior Member
Santiago Lopez Castano
Join Date: Nov 2012
Posts: 254
Rep Power: 10 
Mmmm, I honestly don't think that's the problem. The calculation of the driving force given a bulk flux/velocity returns a scalar, not a field. If you refer to what adjustPhi does, well, mass conservation on the boundaries must be satisfied at all times at each iteration after finishing the pressure correction step of PISO; but this function is only called when pressure on a given boundary are defined as zeroGradient but in your case those BC's should be 'cyclic'.


December 14, 2016, 09:06 

#12 
Member
Charlie Lloyd
Join Date: Feb 2016
Posts: 55
Rep Power: 5 
I have decided to run the problem with a constant forcing term just for a sanity check but I expect the continuity problem will still persist for the reasons you mentioned. I think it is time to investigate snappyHexMesh as an alternative for complex geometries!
Thanks for the advice. 

December 14, 2016, 09:12 

#13  
Senior Member
Santiago Lopez Castano
Join Date: Nov 2012
Posts: 254
Rep Power: 10 
Quote:


December 21, 2016, 08:53 

#14 
Member
Charlie Lloyd
Join Date: Feb 2016
Posts: 55
Rep Power: 5 
As you suggested I have extruded the faces on both boundaries in order to align them normal to the flow direction. However, after running the code again I am still getting spikes on the boundary (see the image below).
The only issue that I can see from my case files are that I have not used the 'preservePatches' function but I have not used it for any of the other cases which have all worked well. Is there anything obviously wrong with the files below? I have also attached a few timesteps from the pimpleFoam log  the only issue I can see is that the bounding flm and fmm values for the LES model are very large, corresponding to the spikes on the boundary. Code:
/** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 3.0.1   \\ / A nd  Web: www.OpenFOAM.org   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object controlDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // application pimpleFoam; startFrom latestTime; startTime 0; stopAt endTime; endTime 15000; deltaT 0.00001; maxCo 1; adjustTimeStep yes; writeControl timeStep; writeInterval 500; purgeWrite 0; writeFormat ascii; writePrecision 8; writeCompression on; timeFormat general; timePrecision 8; runTimeModifiable true; // ************************************************************************* // Code:
/** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 3.0.1   \\ / A nd  Web: www.OpenFOAM.org   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSchemes; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // ddtSchemes { default backward; } gradSchemes { default Gauss linear; } divSchemes { default none; div(phi,U) Gauss filteredLinear; div(phi,flm) Gauss filteredLinear; div(phi,fmm) Gauss filteredLinear; div(phi,nuTilda) Gauss filteredLinear; div((nuEff*dev2(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default corrected; } wallDist { method meshWave; } // ************************************************************************* // Code:
/** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 3.0.1   \\ / A nd  Web: www.OpenFOAM.org   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSolution; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers { p { solver GAMG; tolerance 0; relTol 0.05; smoother GaussSeidel; nPreSweeps 0; nPostSweeps 2; cacheAgglomeration true; nCellsInCoarsestLevel 10; agglomerator faceAreaPair; mergeLevels 1; } pFinal { $p; smoother DICGaussSeidel; tolerance 1e08; relTol 0; } "(UknuTildaflmfmm)" { solver smoothSolver; smoother symGaussSeidel; tolerance 1e08; } "(UknuTildaflmfmm)Final" { $U; } } PIMPLE { nOuterCorrectors 1; nCorrectors 2; nNonOrthogonalCorrectors 0; pRefPoint (0 0 0); pRefValue 0; } // ************************************************************************* // Code:
/** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 3.0.1   \\ / A nd  Web: www.OpenFOAM.org   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class dictionary; note "mesh decomposition control dictionary"; object decomposeParDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // numberOfSubdomains 64; // Keep owner and neighbour on same processor for faces in zones: // preserveFaceZones (heater solid1 solid3); // Keep owner and neighbour on same processor for faces in patches: // (makes sense only for cyclic patches) //preservePatches (cyclic_half0 cyclic_half1); // Keep all of faceSet on a single processor. This puts all cells // connected with a point, edge or face on the same processor. // (just having face connected cells might not guarantee a balanced // decomposition) // The processor can be 1 (the decompositionMethod chooses the processor // for a good load balance) or explicitly provided (upsets balance). //singleProcessorFaceSets ((f0 1)); // Use the volScalarField named here as a weight for each cell in the // decomposition. For example, use a particle population field to decompose // for a balanced number of particles in a lagrangian simulation. // weightField dsmcRhoNMean; // method scotch; // method hierarchical; method simple; // method metis; // method manual; // method multiLevel; // method structured; // does 2D decomposition of structured mesh multiLevelCoeffs { // Decomposition methods to apply in turn. This is like hierarchical but // fully general  every method can be used at every level. level0 { numberOfSubdomains 64; //method simple; //simpleCoeffs //{ // n (2 1 1); // delta 0.001; //} method scotch; } level1 { numberOfSubdomains 4; method scotch; } } // Desired output simpleCoeffs { n (4 4 4); delta 0.001; } hierarchicalCoeffs { n (1 2 1); delta 0.001; order xyz; } metisCoeffs { /* processorWeights ( 1 1 1 1 ); */ } scotchCoeffs { //processorWeights //( // 1 // 1 // 1 // 1 //); //writeGraph true; //strategy "b"; } manualCoeffs { dataFile "decompositionData"; } structuredCoeffs { // Patches to do 2D decomposition on. Structured mesh only; cells have // to be in 'columns' on top of patches. patches (bottomPatch); } //// Is the case distributed? Note: commandline argument roots takes //// precedence //distributed yes; //// Per slave (so nProcs1 entries) the directory above the case. //roots //( // "/tmp" // "/tmp" //); // ************************************************************************* // Code:
/** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 3.0.1   \\ / A nd  Web: www.OpenFOAM.org   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object turbulenceProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // simulationType LES; LES { LESModel dynamicLagrangian; turbulence on; printCoeffs on; delta vanDriest; dynamicLagrangianCoeffs { filter simple; ce 1.048; theta 1.5; } cubeRootVolCoeffs { deltaCoeff 1; } PrandtlCoeffs { delta cubeRootVol; cubeRootVolCoeffs { deltaCoeff 1; } smoothCoeffs { delta cubeRootVol; cubeRootVolCoeffs { deltaCoeff 1; } maxDeltaRatio 1.1; } Cdelta 0.158; } vanDriestCoeffs { delta cubeRootVol; cubeRootVolCoeffs { deltaCoeff 1; } smoothCoeffs { delta cubeRootVol; cubeRootVolCoeffs { deltaCoeff 1; } maxDeltaRatio 1.1; } Aplus 26; Cdelta 0.158; } smoothCoeffs { delta cubeRootVol; cubeRootVolCoeffs { deltaCoeff 1; } maxDeltaRatio 1.1; } } // ************************************************************************* // Code:
Courant Number mean: 0.045302183 max: 1.0263954 deltaT = 8.5202102e05 Time = 4.069965095738113 PIMPLE: iteration 1 smoothSolver: Solving for Ux, Initial residual = 0.0032301636, Final residual = 7.6381712e09, No Iterations 6 smoothSolver: Solving for Uy, Initial residual = 0.0034618555, Final residual = 9.157908e09, No Iterations 10 smoothSolver: Solving for Uz, Initial residual = 0.00032816168, Final residual = 7.7291454e09, No Iterations 4 Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.92385471 GAMG: Solving for p, Initial residual = 0.040015976, Final residual = 0.00082851291, No Iterations 2 time step continuity errors : sum local = 4.160954e09, global = 6.2520222e10, cumulative = 7.3492066e05 Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.92407577 GAMG: Solving for p, Initial residual = 0.0031098961, Final residual = 7.7192676e09, No Iterations 23 time step continuity errors : sum local = 6.2523534e10, global = 6.2520222e10, cumulative = 7.3492691e05 Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.92426602 smoothSolver: Solving for flm, Initial residual = 0.029661143, Final residual = 5.2469417e09, No Iterations 5 bounding flm, min: 667100.81 max: 2002938.7 average: 9.1001611 smoothSolver: Solving for fmm, Initial residual = 0.013253498, Final residual = 3.8549385e09, No Iterations 5 bounding fmm, min: 34953685 max: 2.3359435e+09 average: 13780.572 ExecutionTime = 172032.13 s ClockTime = 172379 s Courant Number mean: 0.044137162 max: 1.0137825 deltaT = 8.4043768e05 Time = 4.070049139505807 PIMPLE: iteration 1 smoothSolver: Solving for Ux, Initial residual = 0.0030180317, Final residual = 7.48811e09, No Iterations 3 smoothSolver: Solving for Uy, Initial residual = 0.0032379585, Final residual = 9.068158e09, No Iterations 3 smoothSolver: Solving for Uz, Initial residual = 0.00030659573, Final residual = 5.8584048e09, No Iterations 2 Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.9253623 GAMG: Solving for p, Initial residual = 0.036520884, Final residual = 0.00062184272, No Iterations 2 time step continuity errors : sum local = 2.953574e09, global = 5.9219333e10, cumulative = 7.3506524e05 Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.92528735 GAMG: Solving for p, Initial residual = 0.0026411813, Final residual = 8.5933446e09, No Iterations 25 time step continuity errors : sum local = 5.9222588e10, global = 5.9219333e10, cumulative = 7.3507117e05 Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.92539805 smoothSolver: Solving for flm, Initial residual = 0.041895575, Final residual = 9.7797156e10, No Iterations 6 bounding flm, min: 848496.88 max: 904626.11 average: 7.4115014 smoothSolver: Solving for fmm, Initial residual = 0.012833257, Final residual = 2.5986852e09, No Iterations 5 bounding fmm, min: 29135461 max: 1.2818239e+09 average: 12415.778 ExecutionTime = 172069.39 s ClockTime = 172416 s Courant Number mean: 0.041233053 max: 0.98170613 deltaT = 8.1079342e05 Time = 4.071023715504507 PIMPLE: iteration 1 smoothSolver: Solving for Ux, Initial residual = 0.0030767242, Final residual = 9.0609907e09, No Iterations 3 smoothSolver: Solving for Uy, Initial residual = 0.0032985252, Final residual = 5.9680932e10, No Iterations 4 smoothSolver: Solving for Uz, Initial residual = 0.00031225705, Final residual = 6.2052363e09, No Iterations 2 Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.92538027 GAMG: Solving for p, Initial residual = 0.040505436, Final residual = 0.00061931897, No Iterations 2 time step continuity errors : sum local = 3.0254994e09, global = 6.0425345e10, cumulative = 7.3507721e05 Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.92538532 GAMG: Solving for p, Initial residual = 0.0029228641, Final residual = 9.70328e09, No Iterations 23 time step continuity errors : sum local = 6.0429122e10, global = 6.0425345e10, cumulative = 7.3508325e05 Pressure gradient source: uncorrected Ubar = 15.555556, pressure gradient = 0.92541771 smoothSolver: Solving for flm, Initial residual = 0.041116896, Final residual = 1.3082178e09, No Iterations 6 bounding flm, min: 981020.44 max: 809933.1 average: 7.3560911 smoothSolver: Solving for fmm, Initial residual = 0.01373234, Final residual = 3.4180578e09, No Iterations 5 bounding fmm, min: 52698474 max: 1.3968272e+09 average: 12280.821 ExecutionTime = 172072.28 s ClockTime = 172419 s Courant Number mean: 0.042001431 max: 0.97433007 deltaT = 8.3215477e05 Time = 4.071106930981983 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
My radial inflow turbine  Abo Anas  CFX  27  May 11, 2018 02:44 
Wrong flow in ratating domain problem  Sanyo  CFX  17  August 15, 2015 07:20 
[Commercial meshers] icem fluent mesh with cyclic boundary condition  jiejie  OpenFOAM Meshing & Mesh Conversion  1  April 5, 2011 04:36 
[snappyHexMesh] snappyHexMesh won't work  zeros everywhere!  sc298  OpenFOAM Meshing & Mesh Conversion  2  March 27, 2011 22:11 
problem when I import mesh with cyclic graduated BC  Cyp  OpenFOAM  0  March 3, 2011 11:38 