# PimpleFoam: High Courant and Instability Problem

 Register Blogs Members List Search Today's Posts Mark Forums Read

 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 z-direction (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 (e-5,e-6 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,423 Rep Power: 25 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. <

June 20, 2014, 09:35
#5
New Member

AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 4
Quote:
 Originally Posted by HanSolo123 My first thought: what happens when you disable adjustTimeStep in your contolDict?
If I try the result is values of Courant number like e93!!!

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_64-linux-gnu/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_64-linux-gnu/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,423 Rep Power: 25 Try changing: 1. Code: `div(phi,U) Gauss linear;` to Code: `div(phi,U) Gauss upwind;` 2. Code: ``` p { solver PCG; preconditioner DIC; tolerance 1e-06; relTol 0; }``` to Code: ``` p { solver PCG; preconditioner DIC; tolerance 1e-06; relTol 0.1; }``` (the same for U) 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; }``` to Code: ```gradSchemes { default cellMDLimited Gauss linear 1.0; }``` Maybe you also need to check your BCs. Can you post them? (and if you disable time step adjustment, sure your time step won't be reducing, but it won't lead to convergence). Caio Martins likes this.

 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 post-processing 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,423 Rep Power: 25 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 { "(p|U)" { tolerance 1e-4; relTol 0; } }``` so the solver will be in pimple loop either 250 iterations (and at this point you know something very wrong), or while residuals for pressure and velocity become lower than 1e-4. 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,423 Rep Power: 25 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).
Attached Files
 TuboPoroso.zip (87.4 KB, 6 views)

 June 20, 2014, 10:16 #12 Senior Member   Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,423 Rep Power: 25 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 { "(p|U)" { tolerance 1e-4; 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,423 Rep Power: 25 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:
 Originally Posted by alexeym 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.
I read in an article (actually, the only article which talks about using porous media for blood simulation) that the solid part of the body (the porous part) must be meshed with a structured mesh, while the first part can be unstructured.

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,423
Rep Power: 25
Quote:
 Originally Posted by agiubilei And I can't use the TopoSet utility at all, so the simpliest way was that.
Using single region mesh and topoSet utility can help locate a problem. Why can't you use topoSet?

 Tags courant number, courant number increasing, instability, pimplefoam

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

All times are GMT -4. The time now is 07:45.