
[Sponsors] 
June 20, 2014, 05:16 
PimpleFoam: High Courant and Instability Problem

#1 
New Member
AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 4 
Good morning,
I'm trying to study the porosity problem with OpenFoam because I have to simulate the porosity in a real geometry of thoracic aorta. So, I've started with a circular cylinder made of two parts. The first part of fluid, and the second part solid, which is the porous part of the cylinder. I've set the simulation with a Re=50. A constant velocity profile in zdirection (the axis direction of the cylinder) and it becomes unstable after a few steps. At the beginning the Courant is very low and the deltaT is that one I set, but after a few steps the deltaT becomes infinitesimal (e5,e6 and so on) and the Co reaches values of 256!!! How can I try to solve this problem and reach a convergence? If you need (and you probably do) I can post my setup, the geometry, mesh and so on. Thanks Alessandro 

June 20, 2014, 08:46 

#2 
New Member
Hans Barósz
Join Date: May 2014
Posts: 22
Rep Power: 4 
My first thought: what happens when you disable adjustTimeStep in your contolDict?


June 20, 2014, 09:14 

#3 
Senior Member
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,564
Rep Power: 26 
Hi,
Can you post your... 1. checkMesh output 2. fvSchemes 3. fvSolution (surely it'll be easier if you just post case files) 

June 20, 2014, 09:26 

#4 
New Member
AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 4 
1.checkMesh output
Code:
Create time Create polyMesh for time = 0 Time = 0 Mesh stats points: 999 faces: 6208 internal faces: 5224 cells: 2696 faces per cell: 4.24036 boundary patches: 4 point zones: 0 face zones: 1 cell zones: 2 Overall number of cells of each type: hexahedra: 0 prisms: 648 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 2048 polyhedra: 0 Checking topology... Boundary definition OK. Cell to face addressing OK. Point usage OK. Upper triangular ordering OK. Face vertices OK. *Number of regions: 2 The mesh has multiple regions which are not connected by any face. <<Writing region information to "0/cellToRegion" Checking patch topology for multiply connected surfaces... Patch Faces Points Surface topology wall 72 26 multiply connected (shared edge) Wall1 840 560 ok (nonclosed singly connected) Outlet1 36 26 ok (nonclosed singly connected) Inlet1 36 26 ok (nonclosed singly connected) <<Writing 26 conflicting points to set nonManifoldPoints Checking geometry... Overall domain bounding box (1 0.999994 0) (1 0.999998 16) Mesh (nonempty, nonwedge) directions (1 1 1) Mesh (nonempty) directions (1 1 1) Boundary openness (4.12825e17 1.3048e17 9.79374e18) OK. Max cell openness = 2.77155e16 OK. Max aspect ratio = 4.92402 OK. Minimum face area = 0.0431419. Maximum face area = 0.254703. Face area magnitudes OK. Min volume = 0.00415344. Max volume = 0.0510398. Total volume = 48.7967. Cell volumes OK. Mesh nonorthogonality Max: 53.7663 average: 16.5653 Nonorthogonality check OK. Face pyramids OK. Max skewness = 0.505704 OK. Coupled point location match (average 0) OK. Mesh OK. End Code:
ddtSchemes { default Euler; } gradSchemes { default Gauss linear; grad(p) Gauss linear; } divSchemes { default none; div(phi,U) Gauss linear; div((nuEff*dev(T(grad(U))))) Gauss linear; } laplacianSchemes { default none; laplacian(nu,U) Gauss linear corrected; laplacian((1A(U)),p) Gauss linear corrected; laplacian(nuEff,U) Gauss linear corrected; laplacian(rAU,p) Gauss linear corrected; } interpolationSchemes { default linear; interpolate(HbyA) linear; } snGradSchemes { default linear; } fluxRequired { default no; p ; } Code:
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 1e06; relTol 0; } UFinal { solver PBiCG; preconditioner DILU; tolerance 1e06; relTol 0; } } PIMPLE { nOuterCorrectors 1; nCorrectors 2; nNonOrthogonalCorrectors 2; pRefCell 0; pRefValue 0; } @HanSolo123: I'll try and I'll tell you the result. 

June 20, 2014, 09:35 

#5  
New Member
AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 4 
Quote:
and then: Code:
DILUPBiCG: Solving for Ux, Initial residual = 0.99612, Final residual = 0.284478, No Iterations 1001 DILUPBiCG: Solving for Uy, Initial residual = 0.996631, Final residual = 0.0366921, No Iterations 1001 #0 Foam::error::printStack(Foam::Ostream&) at ??:? #1 Foam::sigFpe::sigHandler(int) at ??:? #2 in "/lib/x86_64linuxgnu/libc.so.6" #3 double Foam::sumProd<double>(Foam::UList<double> const&, Foam::UList<double> const&) at ??:? #4 Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:? #5 at ??:? #6 at ??:? #7 at ??:? #8 at ??:? #9 __libc_start_main in "/lib/x86_64linuxgnu/libc.so.6" #10 at ??:? Eccezione in virgola mobile (core dump creato) 

June 20, 2014, 09:39 

#6 
Senior Member
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,564
Rep Power: 26 
Try changing:
1. Code:
div(phi,U) Gauss linear; Code:
div(phi,U) Gauss upwind; Code:
p { solver PCG; preconditioner DIC; tolerance 1e06; relTol 0; } Code:
p { solver PCG; preconditioner DIC; tolerance 1e06; relTol 0.1; } 3. Either increase number of outer iterations or switch to residualControl for exiting PIMPLE loop. 4. Further you can try switching from Code:
gradSchemes { default Gauss linear; } Code:
gradSchemes { default cellMDLimited Gauss linear 1.0; } (and if you disable time step adjustment, sure your time step won't be reducing, but it won't lead to convergence). 

June 20, 2014, 09:48 

#7 
New Member
AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 4 
1. Done.
2. Done. 3. Which value of outer iterations should I set? Which files, referring to BCs, do you need? The "0" folder? Excuse me but it's my first time with OpenFoam and I'm still learning, so, thank you for your patience. 

June 20, 2014, 09:57 

#8 
New Member
AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 4 
With the 1. and 2. corrections, the problem reaches the convergence almost instantly and I can postprocessing the results.
But I don't really know if they're correct or not 

June 20, 2014, 09:58 

#9 
Senior Member
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,564
Rep Power: 26 
Well, something greater than 1 Start with 5. Or just switch to residualControl, i.e. you change PIMPLE dictionary to something like:
Code:
PIMPLE { nOuterCorrectors 250; nCorrectors 2; nNonOrthogonalCorrectors 2; pRefCell 0; pRefValue 0; residualControl { "(pU)" { tolerance 1e4; relTol 0; } } Yes, initial and boundary conditions are described in the files in 0 folder. It'll be easier if you just compress it and attach to the message. 

June 20, 2014, 10:01 

#10 
Senior Member
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,564
Rep Power: 26 
Well,
Correction no. 1 is just a switch to a lower order discretisation scheme. To be sure about results you have to use residualControl to check convergence. 

June 20, 2014, 10:07 

#11 
New Member
AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 4 
In the .zip file there are all the files of my simulation ("0", "constant" and "system" folders).


June 20, 2014, 10:16 

#12 
Senior Member
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,564
Rep Power: 26 
ICs and BCs seem to be reasonable. I lost one bracket it my previous post, so the correct PIMPLE dictionary (which is in fvSolution file) is
Code:
PIMPLE { nOuterCorrectors 250; nCorrectors 2; nNonOrthogonalCorrectors 2; pRefCell 0; pRefValue 0; residualControl { "(pU)" { tolerance 1e4; relTol 0; } } } 

June 20, 2014, 10:21 

#13 
New Member
AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 4 
The solutions seems to be reasonable for the velocity. The flux becomes instable almost instantly. The problem is with pression, I think. Because the result shows values of 1e14 for the pressure inside the first fluid part of the body.
So, I think there's something wrong. My purpose  in doing this simple simulation with the cylinder  is to investigate and try to understand the potential of the porous media in order to increase the Re number in my main simulation with the aortic vessel. But now, with your suggestions, the simulation converges. 

June 20, 2014, 10:29 

#14 
Senior Member
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,564
Rep Power: 26 
Well, there could be a problem with connection between two regions. I still can't get why do you need mesh with two regions. You can create single region mesh (if it's just a cylinder, you can even have fully hexagonal mesh), then you create cellZone with topoSet utility and that's it.


June 20, 2014, 10:35 

#15  
New Member
AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 4 
Quote:
And in the main simulation I have the real geometry of thoracic aorta (that must be meshed with an unstructured mesh) and 3 porous cylindrical media that must be structured. So I've followed this way. This preliminary study is only made in order to know if I can set an higher Re number in the real simulation. And I can't use the TopoSet utility at all, so the simpliest way was that. 

June 20, 2014, 10:42 

#16 
Senior Member
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,564
Rep Power: 26 

Tags 
courant number, courant number increasing, instability, pimplefoam 
Thread Tools  
Display Modes  

