CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

PimpleFoam: High Courant and Instability Problem

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

Reply
 
LinkBack Thread Tools Display Modes
Old   June 20, 2014, 05:16
Default PimpleFoam: High Courant and Instability Problem
  #1
New Member
 
AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 3
agiubilei is on a distinguished road
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
agiubilei is offline   Reply With Quote

Old   June 20, 2014, 08:46
Default
  #2
New Member
 
Hans Barósz
Join Date: May 2014
Posts: 18
Rep Power: 3
HanSolo123 is on a distinguished road
My first thought: what happens when you disable adjustTimeStep in your contolDict?
HanSolo123 is offline   Reply With Quote

Old   June 20, 2014, 09:14
Default
  #3
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,132
Rep Power: 20
alexeym will become famous soon enoughalexeym will become famous soon enough
Hi,

Can you post your...

1. checkMesh output
2. fvSchemes
3. fvSolution

(surely it'll be easier if you just post case files)
alexeym is offline   Reply With Quote

Old   June 20, 2014, 09:26
Default
  #4
New Member
 
AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 3
agiubilei is on a distinguished road
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 (non-closed singly connected)  
    Outlet1             36       26       ok (non-closed singly connected)  
    Inlet1              36       26       ok (non-closed singly connected)  
  <<Writing 26 conflicting points to set nonManifoldPoints

Checking geometry...
    Overall domain bounding box (-1 -0.999994 0) (1 0.999998 16)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (4.12825e-17 1.3048e-17 9.79374e-18) OK.
    Max cell openness = 2.77155e-16 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 non-orthogonality Max: 53.7663 average: 16.5653
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.505704 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End
2.fvSchemes

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((1|A(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               ;
}
3.fvSolution

Code:
 
solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0;
    }
pFinal
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0;
    }
    U
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-06;
        relTol          0;
    }
UFinal
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-06;
        relTol          0;
    }
}

PIMPLE
{
    nOuterCorrectors 1;
    nCorrectors     2;
    nNonOrthogonalCorrectors 2;
    pRefCell        0;
    pRefValue       0;
}

@HanSolo123: I'll try and I'll tell you the result.
agiubilei is offline   Reply With Quote

Old   June 20, 2014, 09:35
Default
  #5
New Member
 
AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 3
agiubilei is on a distinguished road
Quote:
Originally Posted by HanSolo123 View Post
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)
agiubilei is offline   Reply With Quote

Old   June 20, 2014, 09:39
Default
  #6
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,132
Rep Power: 20
alexeym will become famous soon enoughalexeym will become famous soon enough
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).
alexeym is offline   Reply With Quote

Old   June 20, 2014, 09:48
Default
  #7
New Member
 
AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 3
agiubilei is on a distinguished road
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.
agiubilei is offline   Reply With Quote

Old   June 20, 2014, 09:57
Default
  #8
New Member
 
AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 3
agiubilei is on a distinguished road
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
agiubilei is offline   Reply With Quote

Old   June 20, 2014, 09:58
Default
  #9
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,132
Rep Power: 20
alexeym will become famous soon enoughalexeym will become famous soon enough
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.
alexeym is offline   Reply With Quote

Old   June 20, 2014, 10:01
Default
  #10
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,132
Rep Power: 20
alexeym will become famous soon enoughalexeym will become famous soon enough
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.
alexeym is offline   Reply With Quote

Old   June 20, 2014, 10:07
Default
  #11
New Member
 
AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 3
agiubilei is on a distinguished road
In the .zip file there are all the files of my simulation ("0", "constant" and "system" folders).
Attached Files
File Type: zip TuboPoroso.zip (87.4 KB, 2 views)
agiubilei is offline   Reply With Quote

Old   June 20, 2014, 10:16
Default
  #12
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,132
Rep Power: 20
alexeym will become famous soon enoughalexeym will become famous soon enough
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;
        }
    }
}
alexeym is offline   Reply With Quote

Old   June 20, 2014, 10:21
Default
  #13
New Member
 
AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 3
agiubilei is on a distinguished road
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.
agiubilei is offline   Reply With Quote

Old   June 20, 2014, 10:29
Default
  #14
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,132
Rep Power: 20
alexeym will become famous soon enoughalexeym will become famous soon enough
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.
alexeym is offline   Reply With Quote

Old   June 20, 2014, 10:35
Default
  #15
New Member
 
AGiub
Join Date: Apr 2014
Posts: 10
Rep Power: 3
agiubilei is on a distinguished road
Quote:
Originally Posted by alexeym View Post
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.
agiubilei is offline   Reply With Quote

Old   June 20, 2014, 10:42
Default
  #16
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,132
Rep Power: 20
alexeym will become famous soon enoughalexeym will become famous soon enough
Quote:
Originally Posted by agiubilei View Post
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?
alexeym is offline   Reply With Quote

Reply

Tags
courant number, courant number increasing, instability, pimplefoam

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On



All times are GMT -4. The time now is 16:38.