CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   PimpleFoam: High Courant and Instability Problem (https://www.cfd-online.com/Forums/openfoam-solving/137651-pimplefoam-high-courant-instability-problem.html)

agiubilei June 20, 2014 06:16

PimpleFoam: High Courant and Instability Problem
 
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 (Post 497962)
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 :confused:

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

1 Attachment(s)
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 (Post 497988)
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 (Post 497990)
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?


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