
[Sponsors] 
Instalibity with pisoFoam in laminar shear flow 

LinkBack  Thread Tools  Search this Thread  Display Modes 
December 18, 2019, 11:47 
Instalibity with pisoFoam in laminar shear flow

#1 
New Member
Marcel Schrader
Join Date: Jan 2017
Posts: 3
Rep Power: 9 
Hello,
I would like to simulate the laminar flow between two planes with pisoFoam in OpenFOAM 4.x with a structured blockMesh. I am using pisoFoam instead of simpleFoam since I want to couple my CFD simulation with a DEM simulation later. I would like to simulate a cubic domain with a length of l=1 µm (80x80x80 cells) and shear rates in the range of 1 1/s up to 10.000 1/s with cyclic boundary conditions at four side walls of the cubic. I choosed the micro unit system to increase the absolute values in the simulation, which should avoid numerical errors due to small numbers. Nevertheless, my simulations are able to generate a linear velocity profile between the planes during the first timesteps with a constant courant number near to 0.05. If the simulation is at or near the steady state (zero iterations of the velocity in xdirection), the courant number starts to increase continuously and strange vortexes are developing. The effect occurs especially for low shear rates which means that the gradient of the velocity is low. I a already changed the unit system to nano, varied the schemes or the channel length without any positive effect. I noticed for higher shear rates (e.g. 1e6 1/s) that the simulation is stable over the time. Could small velocity changes between two cells be the reason for numerical errors which lead to instabilies? Here are my files for a shear rate of 1 1/s and a timestep of 1e3 s. I also attached a logfile. Code:
/** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 1.6   \\ / A nd  Web: http://www.OpenFOAM.org   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class volVectorField; object U; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 1 1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { backCyclic { type cyclic; } frontCyclic { type cyclic; } leftCyclic { type cyclic; } rightCyclic { type cyclic; } top { type fixedValue; value uniform (0.000001 0 0); } bottom { type fixedValue; value uniform (0 0 0); } } // ************************************************************************* // Code:
/** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 1.6   \\ / A nd  Web: http://www.OpenFOAM.org   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class volScalarField; object p; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 2 2 0 0 0 0]; internalField uniform 0; boundaryField { backCyclic { type cyclic; } frontCyclic { type cyclic; } leftCyclic { type cyclic; } rightCyclic { type cyclic; } top { type zeroGradient; } bottom { type zeroGradient; } } // ************************************************************************* // Code:
/** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 5.x   \\ / A nd  Web: www.OpenFOAM.org   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class polyBoundaryMesh; location "constant/polyMesh"; object boundary; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 6 ( top { type patch; nFaces 15000; startFace 4460000; } bottom { type patch; nFaces 15000; startFace 4475000; } leftCyclic { type cyclic; inGroups 1(cyclic); nFaces 10000; startFace 4490000; matchTolerance 0.0001; transform noOrdering; neighbourPatch rightCyclic; } rightCyclic { type cyclic; inGroups 1(cyclic); nFaces 10000; startFace 4500000; matchTolerance 0.0001; transform noOrdering; neighbourPatch leftCyclic; } frontCyclic { type cyclic; inGroups 1(cyclic); nFaces 15000; startFace 4510000; matchTolerance 0.0001; transform noOrdering; neighbourPatch backCyclic; } backCyclic { type cyclic; inGroups 1(cyclic); nFaces 15000; startFace 4525000; matchTolerance 0.0001; transform noOrdering; neighbourPatch frontCyclic; } ) // ************************************************************************* // Code:
/** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 1.6   \\ / A nd  Web: www.OpenFOAM.org   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSolution; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers { p { solver PCG; preconditioner DIC; tolerance 1e06; relTol 0; } pFinal { solver PCG; preconditioner DIC; tolerance 1e06; relTol 0; } U { solver PBiCG; preconditioner DILU; tolerance 1e05; relTol 0; } } PISO { nCorrectors 4; nNonOrthogonalCorrectors 0; pRefCell 0; pRefValue 100000; } // ************************************************************************* // Code:
/** C++ **\  =========    \\ / F ield  OpenFOAM: The Open Source CFD Toolbox   \\ / O peration  Version: 1.6   \\ / A nd  Web: www.OpenFOAM.org   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSchemes; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // ddtSchemes { default Euler; } gradSchemes { default Gauss linear; grad(p) Gauss linear; grad(U) Gauss linear; } divSchemes { default Gauss linear; div(phi,U) Gauss limitedLinearV 1; div((nuEff*dev(grad(U).T()))) Gauss linear; div(U) Gauss linear; } laplacianSchemes { default Gauss linear corrected; laplacian(nuEff,U) Gauss linear corrected; laplacian((1A(U)),p) Gauss linear corrected; laplacian(U) Gauss linear corrected; } interpolationSchemes { default linear; interpolate(U) linear; } snGradSchemes { default corrected; } fluxRequired { default no; p ; } // ************************************************************************* // 

December 19, 2019, 05:15 

#2 
Senior Member
Ruiyan Chen
Join Date: Jul 2016
Location: Hangzhou, China
Posts: 162
Rep Power: 9 
The case seems to be diverging, because the CFL number keeps increasing. I would think this has something to do with the boundary conditions. Why cyclic for the front and back, what about empty? In other words, because the flow is laminar, why 3D instead of 2D? I think a 2D simulation can quickly tell you if your boundary conditions are properly set up or not.


December 20, 2019, 04:12 

#3  
New Member
Marcel Schrader
Join Date: Jan 2017
Posts: 3
Rep Power: 9 
Quote:
Thank you for your reply. I have already tested 2D simulations of the flow, but similiar instabilities were observed. In long term I can not use 2D simulations since I would like to couple the CFD simulation with three dimensional particle structures in DEM. 

December 20, 2019, 08:57 

#4 
Senior Member
Troy Snyder
Join Date: Jul 2009
Location: Akron, OH
Posts: 219
Rep Power: 18 
A simple suggestion I have is to try initializing the entire velocity field to the same or 1/2 the value of the sliding wall velocity, rather than zero.


December 20, 2019, 09:49 

#5 
Senior Member
Peter Hess
Join Date: Apr 2011
Location: Austria
Posts: 250
Rep Power: 17 
I suppose here that the fluid is entering the domain from top
If sliding, then the following is not correct.  The pressure is no where defined!!! How the solver will know, which pressure to be used? The top pressure (at least) need to be defined in my opinion. zeroGraduent in top and bottom is not correct!  The velocity is defined at the top is given. Well, at the bottom is zero. The continuty equation is not satisfied! You need to change the velocity at the bottom to inletOutlet instead of fixed, to calculate the outlet velocity... If the velocity at the top refers to sliding velocity, then it is right defined. 

December 20, 2019, 10:09 

#6  
Senior Member
Troy Snyder
Join Date: Jul 2009
Location: Akron, OH
Posts: 219
Rep Power: 18 
Quote:
The incompressible solver will pick up reference value to set the pressure from pRefCell and pRefValue in the fvSolution file. It does not have to be defined on one of the boundaries. U and P at the top and bottom (solid walls) appear correct to me. 

December 20, 2019, 10:39 

#7 
Senior Member
Peter Hess
Join Date: Apr 2011
Location: Austria
Posts: 250
Rep Power: 17 
Yes! you are right.
Thanks for correcting me. 

December 22, 2019, 06:32 

#8 
Senior Member
Ruiyan Chen
Join Date: Jul 2016
Location: Hangzhou, China
Posts: 162
Rep Power: 9 
Ok, 2D simulations still have the instability problems, and I think your boundary conditions looks fine. What about the schemes, what if you change the div(phi, U) entry to simply Gauss linear? Also, you can check the conservation of mass flow rate between the inlet and the outlet, see if they are indeed conserved.
By the way, since you mentioned that the instability occours when you are near steadystate, what if you initialize the simulation from a steadystate solution(solved by simpleFoam for example), that could help you decide what the problem is. 

Tags 
channel flow, instable, pisofoam, shear 
Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Will the results of steady state solver and transient solver be same?  carye  OpenFOAM Running, Solving & CFD  9  December 28, 2019 05:21 
Calculate shear stress in a laminar fluid flow using simpleFoam  dieg01mp  OpenFOAM PostProcessing  4  October 14, 2019 08:00 
Pipe flow simulation problem with pisoFoam: not getting parabolic profile  sahmed  OpenFOAM Running, Solving & CFD  1  September 7, 2016 14:00 
Validation for laminar Disperse phase flow  shashwat  Main CFD Forum  0  April 4, 2008 02:20 
laminar and turbulent flow in one simulation  msna  FLUENT  0  January 27, 2007 17:35 