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/)
-   -   Instalibity with pisoFoam in laminar shear flow (https://www.cfd-online.com/Forums/openfoam-solving/223036-instalibity-pisofoam-laminar-shear-flow.html)

Marcel_BS December 18, 2019 11:47
Instalibity with pisoFoam in laminar shear flow
 
1 Attachment(s)
Hello,
I would like to simulate the laminar flow between two planes with pisoFoam in OpenFOAM 4.x with a structured blockMesh. I am using pisoFoam instead of simpleFoam since I want to couple my CFD simulation with a DEM simulation later. I would like to simulate a cubic domain with a length of l=1 µm (80x80x80 cells) and shear rates in the range of 1 1/s up to 10.000 1/s with cyclic boundary conditions at four side walls of the cubic. I choosed the micro unit system to increase the absolute values in the simulation, which should avoid numerical errors due to small numbers. Nevertheless, my simulations are able to generate a linear velocity profile between the planes during the first timesteps with a constant courant number near to 0.05. If the simulation is at or near the steady state (zero iterations of the velocity in x-direction), the courant number starts to increase continuously and strange vortexes are developing. The effect occurs especially for low shear rates which means that the gradient of the velocity is low.

I a already changed the unit system to nano, varied the schemes or the channel length without any positive effect. I noticed for higher shear rates (e.g. 1e6 1/s) that the simulation is stable over the time. Could small velocity changes between two cells be the reason for numerical errors which lead to instabilies?

Here are my files for a shear rate of 1 1/s and a timestep of 1e-3 s. I also attached a log-file.


Code:
/*--------------------------------*- C++ -*----------------------------------*\
| =========                |                                                |
| \\      /  F ield        | OpenFOAM: The Open Source CFD Toolbox          |
|  \\    /  O peration    | Version:  1.6                                  |
|  \\  /    A nd          | Web:      http://www.OpenFOAM.org              |
|    \\/    M anipulation  |                                                |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version    2.0;
    format      ascii;
    class      volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


dimensions      [0 1 -1 0 0 0 0];

internalField  uniform (0 0 0);

boundaryField
{
   
    backCyclic
    {
        type    cyclic;
    }
   
    frontCyclic
    {
        type    cyclic;
    }
   
    leftCyclic
    {
        type    cyclic;
    }

    rightCyclic
    {
        type    cyclic;
    }

    top
    {
        type    fixedValue;
        value    uniform (0.000001 0 0);
    }

    bottom
    {
        type    fixedValue;
        value    uniform (0 0 0);
    }

}

 // ************************************************************************* //
Code:
/*--------------------------------*- C++ -*----------------------------------*\
| =========                |                                                |
| \\      /  F ield        | OpenFOAM: The Open Source CFD Toolbox          |
|  \\    /  O peration    | Version:  1.6                                  |
|  \\  /    A nd          | Web:      http://www.OpenFOAM.org              |
|    \\/    M anipulation  |                                                |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version    2.0;
    format      ascii;
    class      volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField  uniform 0;

boundaryField
{
   
    backCyclic
    {
        type    cyclic;
    }
   
    frontCyclic
    {
        type    cyclic;
    }
   
    leftCyclic
    {
        type    cyclic;
    }

    rightCyclic
    {
        type    cyclic;
    }

    top
    {
        type    zeroGradient;
    }

    bottom
    {
        type    zeroGradient;
    }

}

// ************************************************************************* //
Code:
/*--------------------------------*- C++ -*----------------------------------*\
| =========                |                                                |
| \\      /  F ield        | OpenFOAM: The Open Source CFD Toolbox          |
|  \\    /  O peration    | Version:  5.x                                  |
|  \\  /    A nd          | Web:      www.OpenFOAM.org                      |
|    \\/    M anipulation  |                                                |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version    2.0;
    format      ascii;
    class      polyBoundaryMesh;
    location    "constant/polyMesh";
    object      boundary;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

6
(
    top
    {
        type            patch;
        nFaces          15000;
        startFace      4460000;
    }
    bottom
    {
        type            patch;
        nFaces          15000;
        startFace      4475000;
    }
    leftCyclic
    {
        type            cyclic;
        inGroups        1(cyclic);
        nFaces          10000;
        startFace      4490000;
        matchTolerance  0.0001;
        transform      noOrdering;
        neighbourPatch  rightCyclic;
    }
    rightCyclic
    {
        type            cyclic;
        inGroups        1(cyclic);
        nFaces          10000;
        startFace      4500000;
        matchTolerance  0.0001;
        transform      noOrdering;
        neighbourPatch  leftCyclic;
    }
    frontCyclic
    {
        type            cyclic;
        inGroups        1(cyclic);
        nFaces          15000;
        startFace      4510000;
        matchTolerance  0.0001;
        transform      noOrdering;
        neighbourPatch  backCyclic;
    }
    backCyclic
    {
        type            cyclic;
        inGroups        1(cyclic);
        nFaces          15000;
        startFace      4525000;
        matchTolerance  0.0001;
        transform      noOrdering;
        neighbourPatch  frontCyclic;
    }
)

  // ************************************************************************* //
Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                |                                                |
| \\      /  F ield        | OpenFOAM: The Open Source CFD Toolbox          |
|  \\    /  O peration    | Version:  1.6                                  |
|  \\  /    A nd          | Web:      www.OpenFOAM.org                      |
|    \\/    M anipulation  |                                                |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version    2.0;
    format      ascii;
    class      dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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-05;
        relTol          0;
    }

  }

PISO
{
    nCorrectors    4;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue      100000;
}


// ************************************************************************* //
Code:
/*--------------------------------*- C++ -*----------------------------------*\
| =========                |                                                |
| \\      /  F ield        | OpenFOAM: The Open Source CFD Toolbox          |
|  \\    /  O peration    | Version:  1.6                                  |
|  \\  /    A nd          | Web:      www.OpenFOAM.org                      |
|    \\/    M anipulation  |                                                |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version    2.0;
    format      ascii;
    class      dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default        Euler;
}

gradSchemes
{
    default        Gauss linear;
    grad(p)        Gauss linear;
    grad(U)        Gauss linear;
}

divSchemes
{
    default        Gauss linear;
    div(phi,U)      Gauss limitedLinearV 1;
    div((nuEff*dev(grad(U).T()))) Gauss linear;
    div(U)          Gauss linear;
}

laplacianSchemes
{
    default        Gauss linear corrected;
    laplacian(nuEff,U) Gauss linear corrected;
    laplacian((1|A(U)),p) Gauss linear corrected;
    laplacian(U) Gauss linear corrected;
}

interpolationSchemes
{
    default        linear;
    interpolate(U)  linear;
}

snGradSchemes
{
    default        corrected;
}

fluxRequired
{
    default        no;
    p              ;
}


// ************************************************************************* //

cryabroad December 19, 2019 05:15
The case seems to be diverging, because the CFL number keeps increasing. I would think this has something to do with the boundary conditions. Why cyclic for the front and back, what about empty? In other words, because the flow is laminar, why 3D instead of 2D? I think a 2D simulation can quickly tell you if your boundary conditions are properly set up or not.

Marcel_BS December 20, 2019 04:12
Quote:
Originally Posted by cryabroad (Post 752879)
The case seems to be diverging, because the CFL number keeps increasing. I would think this has something to do with the boundary conditions. Why cyclic for the front and back, what about empty? In other words, because the flow is laminar, why 3D instead of 2D? I think a 2D simulation can quickly tell you if your boundary conditions are properly set up or not.

Thank you for your reply. I have already tested 2D simulations of the flow, but similiar instabilities were observed. In long term I can not use 2D simulations since I would like to couple the CFD simulation with three dimensional particle structures in DEM.

tas38 December 20, 2019 08:57
A simple suggestion I have is to try initializing the entire velocity field to the same or 1/2 the value of the sliding wall velocity, rather than zero.

peterhess December 20, 2019 09:49
I suppose here that the fluid is entering the domain from top

If sliding, then the following is not correct.

--------------------------------------

The pressure is no where defined!!!

How the solver will know, which pressure to be used?

The top pressure (at least) need to be defined in my opinion.

zeroGraduent in top and bottom is not correct!

--------------------------------

The velocity is defined at the top is given.

Well, at the bottom is zero.

The continuty equation is not satisfied!

You need to change the velocity at the bottom to inletOutlet instead of fixed, to calculate the outlet velocity...

If the velocity at the top refers to sliding velocity, then it is right defined.

tas38 December 20, 2019 10:09
Quote:
Originally Posted by peterhess (Post 752970)
The pressure is no where defined!!!

How the solver will know, which pressure to be used?

The top pressure (at least) need to be defined in my opinion.

zeroGraduent in top and bottom is not correct!

--------------------------------

I suppose here that the fluid is entering the domain.

The velocity is defined at the top is given.

Well, at the bottom is zero.

The continuty equation is not satisfied!

You need to change the velocity at the bottom to inletOutlet instead of fixed, to calculate the outlet velocity...

If the velocity at the to refers to sliding velocity, then it is right defined.

The incompressible solver will pick up reference value to set the pressure
from pRefCell and pRefValue in the fvSolution file. It does not have to be defined on one of the boundaries.


U and P at the top and bottom (solid walls) appear correct to me.

peterhess December 20, 2019 10:39
Yes! you are right.

Thanks for correcting me.

cryabroad December 22, 2019 06:32
Ok, 2D simulations still have the instability problems, and I think your boundary conditions looks fine. What about the schemes, what if you change the div(phi, U) entry to simply Gauss linear? Also, you can check the conservation of mass flow rate between the inlet and the outlet, see if they are indeed conserved.

By the way, since you mentioned that the instability occours when you are near steady-state, what if you initialize the simulation from a steady-state solution(solved by simpleFoam for example), that could help you decide what the problem is.


All times are GMT -4. The time now is 06:43.