PimpleFoam: High Courant and Instability Problem

agiubilei · June 20, 2014, 06:16

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

HanSolo123 · June 20, 2014, 09:46

My first thought: what happens when you disable adjustTimeStep in your contolDict?

alexeym · June 20, 2014, 10:14

Hi,

Can you post your...

1. checkMesh output
2. fvSchemes
3. fvSolution

(surely it'll be easier if you just post case files)

agiubilei · June 20, 2014, 10:26

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 · June 20, 2014, 10:35

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)

alexeym · June 20, 2014, 10:39

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).

agiubilei · June 20, 2014, 10:48

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 · June 20, 2014, 10:57

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

alexeym · June 20, 2014, 10:58

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 · June 20, 2014, 11:01

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.

agiubilei · June 20, 2014, 11:07

In the .zip file there are all the files of my simulation ("0", "constant" and "system" folders).

alexeym · June 20, 2014, 11:16

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;
        }
    }
}

agiubilei · June 20, 2014, 11:21

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.

alexeym · June 20, 2014, 11:29

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.

agiubilei · June 20, 2014, 11:35

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.

alexeym · June 20, 2014, 11:42

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?

June 20, 2014, 11:16		#12
alexeym Senior Member Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,938 Rep Power: 39	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, 06:16	PimpleFoam: High Courant and Instability Problem	#1
agiubilei New Member AGiub Join Date: Apr 2014 Posts: 10 Rep Power: 12	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, 09:46		#2
HanSolo123 New Member Hans Barósz Join Date: May 2014 Posts: 22 Rep Power: 12	My first thought: what happens when you disable adjustTimeStep in your contolDict?

June 20, 2014, 10:14		#3
alexeym Senior Member Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,938 Rep Power: 39	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, 10:39		#6
alexeym Senior Member Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,938 Rep Power: 39	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, 10:48		#7
agiubilei New Member AGiub Join Date: Apr 2014 Posts: 10 Rep Power: 12	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, 10:57		#8
agiubilei New Member AGiub Join Date: Apr 2014 Posts: 10 Rep Power: 12	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, 10:58		#9
alexeym Senior Member Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,938 Rep Power: 39	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, 11:01		#10
alexeym Senior Member Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,938 Rep Power: 39	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, 11:21		#13
agiubilei New Member AGiub Join Date: Apr 2014 Posts: 10 Rep Power: 12	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, 11:29		#14
alexeym Senior Member Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,938 Rep Power: 39	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.