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

twoPhaseEulerSedFoam scour case diverges

Register Blogs Community New Posts Updated Threads Search

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 1, 2017, 21:25
Default twoPhaseEulerSedFoam scour case diverges
  #1
New Member
 
Lee
Join Date: May 2015
Location: AU
Posts: 2
Rep Power: 0
jylee4 is on a distinguished road
Hello,

Please pardon my ignorance, I am new to OpenFOAM. Could you please point me in the right direction?
I'm trying to model 2D scour beneath a pipeline, under a steady current condition, using twoPhaseEulerSedFoam (link below) on OpenFOAM-2.1.1. Although it was developed using OpenFOAM-2.1.0, the example cases ran fine on OpenFOAM-2.1.1.

https://github.com/csdms-contrib/twophaseeulersedfoam

I hope that the attached "alpha0.jpg" and "Ub0.jpg" can give you a good overview of the initial conditions. The problem is that there are very large time step continuity errors (as shown below)

Code:
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.1.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 2.1.1-221db2718bbb
Exec   : twoPhaseEulerSedFoam
Date   : May 02 2017
Time   : 09:59:46
Host   : "jylee4-VirtualBox"
PID    : 19239
Case   : /home/jylee4/OpenFOAM/jylee4-2.1.1/run/case2h
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Disallowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0


Reading g
Reading transportProperties

Calculating face flux field phia
Calculating face flux field phib
Reading field alpha

Reading field p

kEpsilonCoeffs
{
    Cmu             0.09;
    C1              1.44;
    C2              1.92;
    alphak          1;
    alphaEps        0.76923;
}

wallFunctionCoeffs
{
    kappa           0.41;
    E               9.8;
}

Reading field k

Reading field epsilon

Calculating field nutb

Calculating field nuEffa

Calculating field nuEffb

Calculating field DDtUa and DDtUb

Calculating field g.h

Selecting dragModel for phase a: GidaspowErgunWenYu
Selecting dragModel for phase b: GidaspowErgunWenYu
dragPhase is a
Selecting viscosityModel Gidaspow
Selecting conductivityModel Gidaspow
Selecting radialModel CarnahanStarling
Selecting granularPressureModel Lun
Selecting frictionalStressModel SrivastavaSundaresan
Courant Number mean: 0.00137713 max: 0.0190682

Reading transportProperties for gradP

Initializing with specified pressure gradient:(46.9899 0 0)


PIMPLE: max iterations = 50
    field p    : relTol 0, tolerance 0.0001
    field U    : relTol 0, tolerance 0.0001
    field alpha    : relTol 0, tolerance 0.0001


Starting time loop

Turbulence suspension term is included
Courant Number mean: 0.00137713 max: 0.0190682
Max Ur Courant Number = 0.0190682
Interface Courant Number mean: 0.00137713 max: 0.0190682
deltaT = 1.19904e-05 < 0.0005
Time = 1.19904e-05

PIMPLE: iteration 1
DILUPBiCG:  Solving for alpha, Initial residual = 1.49347e-07, Final residual = 3.21951e-25, No Iterations 3
Dispersed phase volume fraction = 0.182364  Min(alpha) = 0  Max(alpha) = 0.6
DILUPBiCG:  Solving for alpha, Initial residual = 5.486e-10, Final residual = 6.01472e-21, No Iterations 2
Dispersed phase volume fraction = 0.182364  Min(alpha) = 0  Max(alpha) = 0.6
DILUPBiCG:  Solving for Theta, Initial residual = 1, Final residual = 7.43123e-11, No Iterations 8
kinTheory: max(Theta) = 0.000137375 min(Theta) = 0
kinTheory: min(nua) = 0, max(nua) = 2.27681e+32
kinTheory: min(pa) = 0, max(pa) = 5.37874
Max gradPs =119910,Min gradPs =0
Max gradPf =573016,Min gradPf =0
swak4Foam: Allocating new repository for sampledGlobalVariables
GAMG:  Solving for p, Initial residual = 1, Final residual = 6.44256e-10, No Iterations 161
GAMG:  Solving for p, Initial residual = 0.0567147, Final residual = 5.77517e-10, No Iterations 73
time step continuity errors : sum local = 4.97447e-13, global = 3.86504e-14, cumulative = 3.86504e-14
DILUPBiCG:  Solving for alpha, Initial residual = 1.64902e-06, Final residual = 8.65104e-25, No Iterations 4
Dispersed phase volume fraction = 0.182364  Min(alpha) = 0  Max(alpha) = 0.6
DILUPBiCG:  Solving for alpha, Initial residual = 1.25664e-08, Final residual = 6.23708e-24, No Iterations 3
Dispersed phase volume fraction = 0.182364  Min(alpha) = 0  Max(alpha) = 0.6
GAMG:  Solving for p, Initial residual = 0.0058589, Final residual = 8.81012e-10, No Iterations 119
GAMG:  Solving for p, Initial residual = 0.794251, Final residual = 9.101e-10, No Iterations 89
time step continuity errors : sum local = 5.50086e-14, global = 1.83734e-15, cumulative = 4.04877e-14
DILUPBiCG:  Solving for alpha, Initial residual = 1.61092e-06, Final residual = 7.69513e-22, No Iterations 13
Dispersed phase volume fraction = 0.182364  Min(alpha) = 0  Max(alpha) = 0.6
DILUPBiCG:  Solving for alpha, Initial residual = 1.59849e-08, Final residual = 3.91697e-22, No Iterations 8
Dispersed phase volume fraction = 0.182364  Min(alpha) = 0  Max(alpha) = 0.6
PIMPLE: iteration 2
DILUPBiCG:  Solving for alpha, Initial residual = 8.06833e-10, Final residual = 2.4574e-21, No Iterations 10
Dispersed phase volume fraction = 0.182364  Min(alpha) = 0  Max(alpha) = 0.6
DILUPBiCG:  Solving for alpha, Initial residual = 4.38846e-10, Final residual = 4.53047e-22, No Iterations 10
Dispersed phase volume fraction = 0.182364  Min(alpha) = 0  Max(alpha) = 0.6
DILUPBiCG:  Solving for Theta, Initial residual = 6.31348e-17, Final residual = 6.31348e-17, No Iterations 0
kinTheory: max(Theta) = 0.000137375 min(Theta) = 0
kinTheory: min(nua) = 0, max(nua) = 1.81897e+32
kinTheory: min(pa) = 0, max(pa) = 5.37874
Max gradPs =113411,Min gradPs =0
Max gradPf =541822,Min gradPf =0
GAMG:  Solving for p, Initial residual = 1, Final residual = 6.25606e-10, No Iterations 70
GAMG:  Solving for p, Initial residual = 0.000466378, Final residual = 8.81131e-10, No Iterations 52
time step continuity errors : sum local = 5.33473e+07, global = 1.91681e+06, cumulative = 1.91681e+06
DILUPBiCG:  Solving for alpha, Initial residual = 3.45233e-08, Final residual = 3.45233e-08, No Iterations 1001
Dispersed phase volume fraction = 0.182515  Min(alpha) = 0  Max(alpha) = 0.6
DILUPBiCG:  Solving for alpha, Initial residual = 3.45166e-08, Final residual = 1.00473e-07, No Iterations 1001
Dispersed phase volume fraction = 0.374457  Min(alpha) = 0  Max(alpha) = 0.634
GAMG:  Solving for p, Initial residual = 0.999775, Final residual = 7.53284e-10, No Iterations 101
GAMG:  Solving for p, Initial residual = 0.000436472, Final residual = 8.85226e-10, No Iterations 47
time step continuity errors : sum local = 2.54219e+11, global = -2.31365e+10, cumulative = -2.31345e+10
#0  Foam::error::printStack(Foam::Ostream&) in "/home/jylee4/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1  Foam::sigFpe::sigHandler(int) in "/home/jylee4/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2   in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::DILUPreconditioner::calcReciprocalD(Foam::Field<double>&, Foam::lduMatrix const&) in "/home/jylee4/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#4  Foam::DILUPreconditioner::DILUPreconditioner(Foam::lduMatrix::solver const&, Foam::dictionary const&) in "/home/jylee4/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#5  Foam::lduMatrix::preconditioner::addasymMatrixConstructorToTable<Foam::DILUPreconditioner>::New(Foam::lduMatrix::solver const&, Foam::dictionary const&) in "/home/jylee4/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#6  Foam::lduMatrix::preconditioner::New(Foam::lduMatrix::solver const&, Foam::dictionary const&) in "/home/jylee4/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#7  Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/home/jylee4/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#8  Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/home/jylee4/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/lib/libfiniteVolume.so"
#9  Foam::fvMatrix<double>::solve() in "/home/jylee4/OpenFOAM/jylee4-2.1.1/platforms/linux64GccDPOpt/bin/twoPhaseEulerSedFoam"
#10  
 in "/home/jylee4/OpenFOAM/jylee4-2.1.1/platforms/linux64GccDPOpt/bin/twoPhaseEulerSedFoam"
#11  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#12  
 in "/home/jylee4/OpenFOAM/jylee4-2.1.1/platforms/linux64GccDPOpt/bin/twoPhaseEulerSedFoam"
Floating point exception (core dumped)
How did I end up in this mess:

1. Copied the "Sumer1996" steady flow example case for twoPhaseEulerSedFoam, and replaced the mesh via fluent3DMeshToFoam (checkMesh.log provided below)

checkMesh.log:
Code:
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.1.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 2.1.1-221db2718bbb
Exec   : checkMesh
Date   : May 02 2017
Time   : 09:14:18
Host   : "jylee4-VirtualBox"
PID    : 18757
Case   : /home/jylee4/OpenFOAM/jylee4-2.1.1/run/case2h
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Disallowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           224000
    internal points:  0
    faces:            445612
    internal faces:   221612
    cells:            111204
    boundary patches: 8
    point zones:      0
    face zones:       1
    cell zones:       1

Overall number of cells of each type:
    hexahedra:     111204
    prisms:        0
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    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: 1 (OK).

Checking patch topology for multiply connected surfaces ...
    Patch               Faces    Points   Surface topology                  
    bed                 448      898      ok (non-closed singly connected)  
    fluidin             149      300      ok (non-closed singly connected)  
    fluidout            149      300      ok (non-closed singly connected)  
    pipe                200      400      ok (non-closed singly connected)  
    sedin               99       200      ok (non-closed singly connected)  
    sedout              99       200      ok (non-closed singly connected)  
    symsides            222408   224000   ok (non-closed singly connected)  
    top                 448      898      ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-1 0 0) (1 0.33 0.005)
    Mesh (non-empty, non-wedge) directions (1 1 0)
    Mesh (non-empty) directions (1 1 0)
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (3.79223e-21 8.50472e-18 1.63851e-15) OK.
    Max cell openness = 2.89168e-16 OK.
    Max aspect ratio = 20.0908 OK.
    Minumum face area = 5.62712e-07. Maximum face area = 8.99651e-05.  Face area magnitudes OK.
    Min volume = 2.81356e-09. Max volume = 3.6539e-07.  Total volume = 0.00329018.  Cell volumes OK.
    Mesh non-orthogonality Max: 44.0733 average: 10.1293
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.456248 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End
2. then edited the boundary conditions in the "0" folder. I used groovyBC to specify a log velocity profile for the fluid phase at the inlet:

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

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

internalField       uniform (0 0 0);

boundaryField
{
    bed
    {
        type               fixedValue;
        value              uniform (0 0 0);
    }
    fluidin
    {
    type         groovyBC;
        variables
    (
        "uf=0.043914724;" // friction velocity, m/s
        "z0=3e-05;" // bed roughness length, m
        "kappaa=0.40;" // von Karman's constant (Soulsby, 1997)
        //"height=pos().y-0.1+z0;" // elevation - sediment box height; CAN'T BE 0
        "height=sqrt(sqr(pos().y-0.1+z0));" // elevation - sediment box height; CAN'T BE 0
        "velx=(uf/kappaa)*log(height/z0)*normal();"
    );
        valueExpression "-velx";
    value         uniform (0 0 0);
    }
    fluidout
    {
        //type                zeroGradient;
        type                inletOutlet;
        inletValue          uniform (0 0 0); //no backflow
        value               uniform ( 0 0 0 ); //zero gradient for outflow
    }
    pipe
    {
        type               fixedValue;
        value              uniform (0 0 0);
    }
    sedin
    {
        type                zeroGradient;
        //type        fixedValue;
    //value        uniform (0 0 0);
    }
    sedout
    {
        //type                zeroGradient;
        type                inletOutlet;
        inletValue          uniform (0 0 0); //no backflow
        value               uniform ( 0 0 0 ); //zero gradient for outflow
    }
    symsides
    {
        type            empty;
    }
    top
    {
        type               symmetryPlane;
        //type                slip;
        //type            pressureInletOutletVelocity;
        //value           uniform (0.982 0 0);
    }
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
3. used setFieldsDict and funkySetFieldsDict to have non-uniform internal fields for alpha (volume fraction) and Ub (fluid velocity).

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

defaultFieldValues /* outside sediment box */
(
    volScalarFieldValue alpha 0
);

regions
(
    boxToCell
    {
        box ( -1 0 0 ) ( 1 0.1 0.005 ) ; // inside sediment box

        fieldValues
        (
            volScalarFieldValue alpha 0.6 //6.107060e-01
        );
    }
);
system/funkySetFieldsDict:
Code:
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.1.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version         2.0;
    format          ascii;
    //root            "~/OpenFOAM/jylee4-2.1.1/run";
    //case            "case2h";
    //instance        "system";
    //local           "";
    class           dictionary;
    location        "system";
    object          funkySetFieldsDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

expressions
(
    initialVelProfile
    {
        field Ub; // field to initialise
        condition "pos().y>0.1 && pos().y<0.33";
        keepPatches true; //keep the boundary conditions that were set before
        //expression "vector(1,0,0)*(0.043914724/0.4)*log((pos().y)/3e-05)";
        expression "vector(1,0,0)*(0.043914724/0.4)*log((sqrt(sqr(pos().y-0.1+3e-05)))/3e-05)";
        dimension [0 1 -1 0 0];
    }
);
// ************************************************************************* //
4. ran:
Code:
setFields
funkySetFields -time 0
twoPhaseEulerSedFoam
5. solution diverges; so I tried copy and pasting the mesh and "0" folder into a twoPhaseEulerFoam tutorial case (e.g. bed2), and it ran fine.

6. back to twoPhaseEulerSedFoam, I reduced the nOuterCorrectors in the PIMPLE loop to 1 (in fvSolution), and reduced the writeInterval to a value which is equal to deltaT (in controlDict), to visualise the output.
With reference to the attached "p1.jpg" and "Ub1.jpg" files, it seems that there are very large pressure and velocity values at the outlet, or more specifically, at the interface between the fluid and sediment phase. However, I have no idea how to fix this.
I've tried using zeroGradient and fixedValue for the outlet boundary for p, zeroGradient and inletOutlet for the outlet for Ua and Ub, but that didn't work.

Please find the files here:
https://www.dropbox.com/s/ka99wfmfnx...2h.tar.gz?dl=0

Many thanks for your patience and kind assistance.
Attached Images
File Type: jpg p1.jpg (27.1 KB, 50 views)
File Type: jpg Ub1.jpg (26.9 KB, 42 views)
File Type: jpg mesh.jpg (180.5 KB, 48 views)
File Type: jpg alpha0.jpg (23.6 KB, 51 views)
File Type: jpg Ub0.jpg (30.2 KB, 58 views)
jylee4 is offline   Reply With Quote

Old   October 16, 2017, 06:25
Default
  #2
New Member
 
Join Date: Nov 2016
Posts: 4
Rep Power: 9
asin is on a distinguished road
hey,I'm doing some research on scour recently,do you solve your problem?
asin is offline   Reply With Quote

Old   November 10, 2017, 10:50
Default
  #3
New Member
 
Tina Ebrahimi
Join Date: Oct 2017
Posts: 4
Rep Power: 8
Tina Ebrahimi is on a distinguished road
Quote:
Originally Posted by asin View Post
hey,I'm doing some research on scour recently,do you solve your problem?
Hello Asin. I want to simulate scouring in OpenFOAM too. Can I ask you for some help?
Tina Ebrahimi is offline   Reply With Quote

Old   November 12, 2017, 12:32
Default
  #4
Member
 
Declan
Join Date: Oct 2016
Location: Ireland
Posts: 40
Rep Power: 9
decah is on a distinguished road
Try checking out scourFoam

Erosion/sediment transport with OpenFOAM: Which model to use?

Or a Larangian solver depending on the scale of your problem
decah is offline   Reply With Quote

Reply

Tags
twophaseeulersedfoam


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


Similar Threads
Thread Thread Starter Forum Replies Last Post
[DesignModeler] DesignModeler Scripting: How to get Full Command Access ANT ANSYS Meshing & Geometry 53 February 16, 2020 15:13
Is Playstation 3 cluster suitable for CFD work hsieh OpenFOAM 9 August 16, 2015 14:53
Superlinear speedup in OpenFOAM 13 msrinath80 OpenFOAM Running, Solving & CFD 18 March 3, 2015 05:36
Transient case running with a super computer microfin FLUENT 0 March 31, 2009 11:20
Turbulent Flat Plate Validation Case Jonas Larsson Main CFD Forum 0 April 2, 2004 10:25


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