[mayank.dce2k7]

Hi Foamers,

I am trying to calculate transient flow field in a single bifurcating tube(using OF-2.1.1 and swak4Foam) with inlet boundary condition varying with space and time:

I will start illustrating my case by giving 'checkMesh' report:

Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
points: 140311
faces: 1127971
internal faces: 1083988
cells: 521531
boundary patches: 4
point zones: 0
face zones: 1
cell zones: 1

Overall number of cells of each type:
hexahedra: 0
prisms: 125835
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 395696
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
INLET 910 612 ok (non-closed singly connected)
OUTLET1 562 390 ok (non-closed singly connected)
OUTLET2 566 394 ok (non-closed singly connected)
WALL 41945 21065 ok (non-closed singly connected)

Checking geometry...
Overall domain bounding box (-0.108847 -0.161672 -0.0190498) (0.108851 0.127 0.0190498)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (2.61084e-17 -1.48172e-16 3.15641e-17) OK.
Max cell openness = 2.91142e-16 OK.
Max aspect ratio = 7.00156 OK.
Minumum face area = 1.85341e-07. Maximum face area = 3.64501e-06. Face area magnitudes OK.
Min volume = 5.00262e-11. Max volume = 2.11586e-09. Total volume = 0.000347759. Cell volumes OK.
Mesh non-orthogonality Max: 63.6578 average: 14.781
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 2.94579 OK.
Coupled point location match (average 0) OK.

Mesh OK.

End

Boundary Conditions:

0/U

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

internalField  uniform (0 0 0);

boundaryField
{
    INLET
    {
                  
  
     type            groovyBC;

     variables "vavg=(-0.000163-0.00336*cos(time()*0.3254)-0.3416*sin(time()*0.3254));r2=(pow(pos().x,2)+pow(pos().z,2));R2=sum(area())/pi;para=-((R2-r2)/R2)*normal();";
     valueExpression "2*vavg*para";
     value        uniform (0 0 0);
 
    }

    OUTLET1
    {
       type  zeroGradient;
    }

    OUTLET2
    {
        type            zeroGradient;
     
}

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

   
}

0/p

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

internalField   uniform 0;

boundaryField
{
    INLET
    {
        type            zeroGradient;
    }

    OUTLET1
    {
      type            fixedValue;
      value           uniform 0;
    }
   
    OUTLET2
    {
      type            fixedValue;
      value           uniform 0;  
    }  

    WALL
    {
        type            zeroGradient;
    }

    
}

controlDict

libs ( "libOpenFOAM.so"
      "libsimpleSwakFunctionObjects.so"
      "libswakFunctionObjects.so"
      "libgroovyBC.so" ) ;
application     pimpleFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         20;

deltaT          0.02;

writeControl    adjustableRunTime;

writeInterval   0.4;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression uncompressed;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

adjustTimeStep  no;

maxCo           100;

functions {
    initSwak {
        // needed to allow DataEntry to work
        type initSwakFunctionObject;
        region region0;
    }
};

functions
{
    volFlowInlet
    {
        type                swakExpression;
        valueType           patch;
        patchName           INLET;
        
        verbose             true;
        expression          "phi";
        accumulations       ( sum );
    }
    
    volFlowOutlet1
    {
        type                swakExpression;
        valueType           patch;
        patchName           OUTLET1;
        
        verbose             true;
        expression          "phi";
        accumulations       ( sum );
    }

    volFlowOutlet2
    {
        type                swakExpression;
        valueType           patch;
        patchName           OUTLET2;
        
        verbose             true;
        expression          "phi";
        accumulations       ( sum );
    }
}

fvSchemes

ddtSchemes
{
    default         CrankNicholson 1;
}

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

divSchemes
{
    default         none;
    div(phi,U)      Gauss linearUpwind grad(U);
    div(phi,k)      Gauss limitedLinear 1;
    div(phi,epsilon) Gauss limitedLinear 1;
    div(phi,R)      Gauss limitedLinear 1;
    div(R)          Gauss linear;
    div(phi,nuTilda) Gauss limitedLinear 1;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         none;
    laplacian(nuEff,U) Gauss linear corrected;
    laplacian((1|A(U)),p) Gauss linear corrected;
    laplacian((1|A(U)),p_rgh) Gauss linear corrected;
    laplacian(rAU,p_rgh) Gauss linear corrected;
    laplacian(DkEff,k) Gauss linear corrected;
    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
    laplacian(DREff,R) Gauss linear corrected;
    laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
}

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

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p               ;
    p_rgh           ;
}

fvSolution

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0.01;
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    pFinal
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0.01;
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    "(U|kl|kt|omega)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-06;
        relTol          0;
    }

    "(U|kl|kt|omega)Final"
    {
        $U;
        tolerance       1e-06;
        relTol          0;
    }
}

PIMPLE
{
    nOuterCorrectors    50;
    nCorrectors         1;
    nNonOrthogonalCorrectors 2;
    pRefCell            0;
    pRefValue           0;
    residualControl
{
"(p|U)"
{
tolerance 1e-4;
relTol 0;
}
}
}

relaxationFactors
{
    fields
    {
    }
    equations
    {
        "U.*"           0.5;
        "kl.*"          1;
        "kt.*"          1;
        "omega.*"       1;
    }
}

The simulation blows blows at time=0.9s or little higher. The pulse lasts 20s and the reynolds number peaks(Re=2200) at 5s and 15s. If I choose adjustable time step then the it blows by showing time step in range of 10^-13 or so. If I choose a fixed time step of 0.01 or 0.02 or 0.005 it blows showing high courant number or unusually high flow rates. I have employed meshes of various kinds.....even pure tetrahedral mesh doesn't seem to work.

Here is the link to download the case directory along with the mesh:

https://www.dropbox.com/sh/u7sqlfvqtc6uynq/1MMKNryk7Z

Following are the error message snippets:

Expression volFlowInlet : sum=0.000155198
Expression volFlowOutlet1 : sum=-1.01936e+12
Expression volFlowOutlet2 : sum=1.01938e+12
Courant Number mean: 4.53531e+15 max: 9.33679e+17
Time = 1.23

Also the time step continuity errors goes exponentially high:

PIMPLE: iteration 50
DILUPBiCG: Solving for Ux, Initial residual = 0.0124649, Final residual = 3.77817e-07, No Iterations 5
DILUPBiCG: Solving for Uy, Initial residual = 0.00973289, Final residual = 3.76167e-07, No Iterations 5
DILUPBiCG: Solving for Uz, Initial residual = 0.0198487, Final residual = 4.315e-08, No Iterations 6
GAMG: Solving for p, Initial residual = 0.0537055, Final residual = 0.00027507, No Iterations 3
GAMG: Solving for p, Initial residual = 0.00744695, Final residual = 6.84879e-05, No Iterations 3
GAMG: Solving for p, Initial residual = 0.00143221, Final residual = 8.22524e-06, No Iterations 3
time step continuity errors : sum local = 3.47333e+09, global = -7.64618e+08, cumulative = -2.18106e+10
PIMPLE: not converged within 50 iterations
ExecutionTime = 55306.7 s ClockTime = 88055 s

Velocity fields have exploded and there is a flow reversal too. I tried avoiding flow reversal using outletinlet boudary condition since my outlets should behave as inlets in the beginning:

Code:
type  outletInlet;
outletValue uniform 0;
value uniform 0;
but it didn't help either.

Please suggest me how to make it work .