Instalibity with pisoFoam in laminar shear flow

Marcel_BS · December 18, 2019, 12:47

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, 06: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, 05:12

Quote:

Originally Posted by cryabroad

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, 09: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, 10: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, 11:09

Quote:

Originally Posted by peterhess

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, 11:39

Yes! you are right.

Thanks for correcting me.

cryabroad · December 22, 2019, 07: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.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Will the results of steady state solver and transient solver be same?	carye	OpenFOAM Running, Solving & CFD	9	December 28, 2019 06:21
Calculate shear stress in a laminar fluid flow using simpleFoam	dieg01mp	OpenFOAM Post-Processing	4	October 14, 2019 09:00
Pipe flow simulation problem with pisoFoam: not getting parabolic profile	sahmed	OpenFOAM Running, Solving & CFD	1	September 7, 2016 15:00
Validation for laminar Disperse phase flow	shashwat	Main CFD Forum	0	April 4, 2008 03:20
laminar and turbulent flow in one simulation	msna	FLUENT	0	January 27, 2007 18:35

December 19, 2019, 06:15		#2
cryabroad Senior Member Ruiyan Chen Join Date: Jul 2016 Location: Hangzhou, China Posts: 162 Rep Power: 10	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.

December 20, 2019, 09:57		#4
tas38 Senior Member Troy Snyder Join Date: Jul 2009 Location: Akron, OH Posts: 220 Rep Power: 19	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.

December 20, 2019, 10:49		#5
peterhess Senior Member Peter Hess Join Date: Apr 2011 Location: Austria Posts: 250 Rep Power: 17	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.

December 20, 2019, 11:39		#7
peterhess Senior Member Peter Hess Join Date: Apr 2011 Location: Austria Posts: 250 Rep Power: 17	Yes! you are right. Thanks for correcting me.

December 22, 2019, 07:32		#8
cryabroad Senior Member Ruiyan Chen Join Date: Jul 2016 Location: Hangzhou, China Posts: 162 Rep Power: 10	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.