[bijan darbari]

Hello

I'm solving flow and heat transfer of water between two parallel plates with 100 micrometer solid particles (eulerian - lagrangian ).
the flow is laminar. the solver is "reactingParcelFoam".


when the "BrownianMotionForce" activated. many particles escaped from inlet and courantnumber get very high value!!!! so the solver diverged.

the case file attached.

solutionlog:

Create time

Create mesh for time = 5000


Reading g

PIMPLE: Operating solver in PISO mode

Creating combustion model

Selecting combustion model PaSR<rhoChemistryCombustion>
Selecting chemistry type 
{
    chemistrySolver noChemistrySolver;
    chemistryThermo rho;
}

Selecting thermodynamics package 
{
    type            heRhoThermo;
    mixture         reactingMixture;
    transport       polynomial;
    thermo          hPolynomial;
    energy          sensibleEnthalpy;
    equationOfState icoPolynomial;
    specie          specie;
}

Selecting chemistryReader foamChemistryReader
chemistryModel: Number of species = 1 and reactions = 0
    using integrated reaction rate
Creating component thermo properties:
    multi-component carrier - 1 species
    liquids - 1 components
    solids - 1 components

Reading field U

Reading/calculating face flux field phi

Creating turbulence model

Selecting turbulence model type laminar
Creating field dpdt

Creating field kinetic energy K

Creating multi-variate interpolation scheme

Selecting radiationModel none

Constructing reacting cloud
Constructing particle forces
    Selecting particle force sphereDrag
    Selecting particle force gravity
    Selecting particle force BrownianMotion
Constructing cloud functions
    Selecting cloud function patchPostProcessing1 of type patchPostProcessing
Constructing particle injection models
Creating injector: model1
Selecting injection model patchInjection
    Constructing 2-D injection
    Choosing nParticle to be a fixed value, massTotal variable now does not determine anything.
Selecting distribution model fixedValue
Selecting dispersion model none
Selecting patch interaction model localInteraction
    Interaction fields will not be written
Selecting stochastic collision model none
Selecting surface film model none
Selecting U integration scheme Euler
Selecting heat transfer model none
Selecting T integration scheme analytical
Selecting composition model singleMixtureFraction
Selecting phase change model none
Selecting devolatilisation model none
Selecting surface reaction model none
No MRF models present

No finite volume options present

Courant Number mean: 0.0512453171 max: 0.1219978729

Starting time loop

faceSource faceSource1:
    total faces  = 40
    total area   = 0.001
    weight field = phi


Courant Number mean: 0.0512453171 max: 0.1219978729
Time = 5001


Solving 2-D cloud reactingCloud1

Cloud: reactingCloud1 injector: model1
    Added 95 new parcels

Cloud: reactingCloud1
    Current number of parcels       = 40
    Current mass in system          = 0.0001675516082
    Linear momentum                 = (-0.0896473776 -0.005854519442 0)
   |Linear momentum|                = 0.08983834208
    Linear kinetic energy           = 1649.634654
    model1:
        number of parcels added     = 95
        mass introduced             = 0.0003979350695
    Parcel fate (number, mass)      : patch frontAndBack
      - escape                      = 0, 0
      - stick                       = 0, 0
    Parcel fate (number, mass)      : patch top
      - escape                      = 0, 0
      - stick                       = 20, 8.37758041e-05
    Parcel fate (number, mass)      : patch bottom
      - escape                      = 0, 0
      - stick                       = 10, 4.188790205e-05
    Parcel fate (number, mass)      : patch inlet
      - escape                      = 55, 0.0002303834613
      - stick                       = 0, 0
    Parcel fate (number, mass)      : patch outlet
      - escape                      = 0, 0
      - stick                       = 0, 0
    Temperature min/max             = 273, 273
    Mass transfer phase change      = 0
    Mass transfer devolatilisation  = 0
    Mass transfer surface reaction  = 0

diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
rho min/max = 999.9999951, 1000.000006
PIMPLE: iteration 1
smoothSolver:  Solving for Ux, Initial residual = 0.9999999582, Final residual = 1.735509726e-07, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 0.9999999903, Final residual = 1.878729947e-07, No Iterations 1
DILUPBiCG:  Solving for h, Initial residual = 1, Final residual = 3.671871817e-07, No Iterations 2
T gas min/max   = -5675.506981, 3396.376588
GAMG:  Solving for p, Initial residual = 0.9999999949, Final residual = 0.01993556705, No Iterations 1
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
rho min/max = -24208271.61, 12617994.23
time step continuity errors : sum local = 300.6960653, global = -3.598339041e-06, cumulative = -3.598339041e-06
p min/max = -94020684.85, 11979377.26
GAMG:  Solving for p, Initial residual = 0.0331384453, Final residual = 0.001899684118, No Iterations 1
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
rho min/max = -3642313.277, 7998449.618
time step continuity errors : sum local = 59.72799648, global = -8.171689074e-05, cumulative = -8.531522978e-05
p min/max = -93995055.47, 11934471.79
GAMG:  Solving for p, Initial residual = 0.002625535523, Final residual = 6.848064135e-07, No Iterations 7
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
rho min/max = 823.5085624, 1120.188828
time step continuity errors : sum local = 0.02162576412, global = 0.0003361922804, cumulative = 0.0002508770506
p min/max = -93994845.23, 11910067.61
ExecutionTime = 0.68 s  ClockTime = 0 s

Courant Number mean: 23243.26381 max: 5386047.497
Time = 5002


Solving 2-D cloud reactingCloud1

Cloud: reactingCloud1 injector: model1
    Added 95 new parcels

Cloud: reactingCloud1
    Current number of parcels       = 55
    Current mass in system          = 0.0002303834613
    Linear momentum                 = (0.01392251322 -0.01478103055 0)
   |Linear momentum|                = 0.02030554698
    Linear kinetic energy           = 1486.082312
    model1:
        number of parcels added     = 190
        mass introduced             = 0.0007958701389
    Parcel fate (number, mass)      : patch frontAndBack
      - escape                      = 0, 0
      - stick                       = 0, 0
    Parcel fate (number, mass)      : patch top
      - escape                      = 0, 0
      - stick                       = 23, 9.634217471e-05
    Parcel fate (number, mass)      : patch bottom
      - escape                      = 0, 0
      - stick                       = 20, 8.37758041e-05
    Parcel fate (number, mass)      : patch inlet
      - escape                      = 135, 0.0005654866776
      - stick                       = 0, 0
    Parcel fate (number, mass)      : patch outlet
      - escape                      = 0, 0
      - stick                       = 0, 0
    Temperature min/max             = 273, 273
    Mass transfer phase change      = 0
    Mass transfer devolatilisation  = 0
    Mass transfer surface reaction  = 0

diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
rho min/max = 823.5085624, 1120.188828
PIMPLE: iteration 1
smoothSolver:  Solving for Ux, Initial residual = 0.5225870011, Final residual = 2.91586547e-07, No Iterations 5
smoothSolver:  Solving for Uy, Initial residual = 0.5192037232, Final residual = 5.095801868e-07, No Iterations 5
DILUPBiCG:  Solving for h, Initial residual = 0.9275706068, Final residual = 6.899198912e-08, No Iterations 5


--> FOAM FATAL ERROR: 
Maximum number of iterations exceeded

    From function Foam::scalar Foam::species::thermo<Thermo, Type>::T(Foam::scalar, Foam::scalar, Foam::scalar, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar) const) const [with Thermo = Foam::hPolynomialThermo<Foam::icoPolynomial<Foam::specie, 8>, 8>; Type = Foam::sensibleEnthalpy; Foam::scalar = double; Foam::species::thermo<Thermo, Type> = Foam::species::thermo<Foam::hPolynomialThermo<Foam::icoPolynomial<Foam::specie, 8>, 8>, Foam::sensibleEnthalpy>]
    in file /home/openfoam/OpenFOAM/OpenFOAM-3.0.1/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 66.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::error::abort() at ??:?
#2  Foam::species::thermo<Foam::hPolynomialThermo<Foam::icoPolynomial<Foam::specie, 8>, 8>, Foam::sensibleEnthalpy>::THs(double, double, double) const at ??:?
#3  Foam::heRhoThermo<Foam::rhoReactionThermo, Foam::SpecieMixture<Foam::reactingMixture<Foam::polynomialTransport<Foam::species::thermo<Foam::hPolynomialThermo<Foam::icoPolynomial<Foam::specie, 8>, 8>, Foam::sensibleEnthalpy>, 8> > > >::calculate() at ??:?
#4  Foam::heRhoThermo<Foam::rhoReactionThermo, Foam::SpecieMixture<Foam::reactingMixture<Foam::polynomialTransport<Foam::species::thermo<Foam::hPolynomialThermo<Foam::icoPolynomial<Foam::specie, 8>, 8>, Foam::sensibleEnthalpy>, 8> > > >::correct() at ??:?
#5  ? at ??:?
#6  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#7  ? at ??:?
Aborted (core dumped)
bijan@bijan-virtual-machine:~/Desktop/mycase

reactingcloud1properties:

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

solution
{
    active          true;
    coupled         true;
    transient       yes;
 
    cellValueSourceCorrection on;
    maxCo           0.1;

    sourceTerms
    {
        schemes
        {
            rho             semiImplicit 1;
            U               semiImplicit 1;
            Yi              semiImplicit 1;
            h               semiImplicit 1;
            radiation       explicit 1;
        } 
      
    }

    interpolationSchemes
    {
        rho             cell;
        U               cellPoint;
        thermo:mu       cell;
        T               cell;
        Cp              cell;
        kappa           cell;
        p               cell;
    }

    integrationSchemes
    {
        U               Euler;
        T               analytical;
    }
}


constantProperties
{
    rho0            8000;
    T0              273;
    Cp0             710;

    constantVolume  true;
}


subModels
{
    particleForces
    {
        sphereDrag;
        gravity;
        BrownianMotion
    {
    lambda     1E-10;
    turbulence false;
    }

    }

    injectionModels
    {
         model1
        {
            type            patchInjection;
            SOI             5000;
            massTotal       0;
            parcelBasisType fixed;
            patchName       inlet;
            duration        10000;
            parcelsPerSecond 95;
            nParticle       1000;
            U0              (0.001 0 0);
            flowRateProfile constant 1;
            
            sizeDistribution
            {
                type        fixedValue;
                fixedValueDistribution
                {
                    value   100e-6;
                }
            }
        }
    }

    dispersionModel none;

    patchInteractionModel localInteraction;

    heatTransferModel none;

    compositionModel singleMixtureFraction;

    radiation       off;

    localInteractionCoeffs
    {
        patches
        (
            frontAndBack
            {
                type rebound;
                 e    1;
                mu   0;
            }

            top
            {
                type stick;
    
            }

            bottom
            {
                type stick;
              
            }
            
            inlet
            {
                type escape;
           
            }

            outlet
            {
                type escape;
            }
        );
    }

    singleMixtureFractionCoeffs
    {
        phases
        (
            gas
            {
            }
            liquid
            {
            
            }
            solid
            {
            C      1;
            }
        );
        YGasTot0        0;
        YLiquidTot0     0;
        YSolidTot0      1;
    }

   
}


cloudFunctions
{
    patchPostProcessing1
    {
        type            patchPostProcessing;
        maxStoredParcels 1000000;
        patches         ( outlet );
    }
}


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

controldict:

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

application     reactingParcelFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         10000;

deltaT          1;

writeControl    adjustableRunTime;

writeInterval  10;

purgeWrite      0;

writeFormat     ascii;

writePrecision  10;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

adjustTimeStep  no;

maxCo           6;

maxDeltaT       1;

functions
{
    faceSource1
    {
        type            faceSource;
        functionObjectLibs ("libfieldFunctionObjects.so");
        enabled         yes;
        outputControl   outputTime;
        log             yes;
        valueOutput     no;
        source          patch;
        sourceName      outlet;
        operation       weightedAverage;
        weightField     phi;

        fields
        (
            H2O
            T
        );
    }
}

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

fvschemes:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.1                                 |
|   \\  /    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;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss upwind;
    div(phid,p)     Gauss upwind;
    div(phi,K)      Gauss upwind;
    div(phi,h)      Gauss upwind;
    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;
    div(phi,omega) Gauss upwind;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
    div(phi,Yi_h)   Gauss upwind;
}

laplacianSchemes
{
    default         Gauss linear uncorrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         uncorrected;
}

wallDist
{
    method meshWave;
}


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

fvsolution:

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

solvers
{
    rho
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-05;
        relTol          0.1;
    }

    rhoFinal
    {
        $rho;
        tolerance       1e-05;
        relTol          0;
    }

    "(U|k|omega)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-06;
        relTol          0.1;
    }

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

    p
    {
        solver          GAMG;
        tolerance       1e-6;
        relTol          0.1;
        smoother        DICGaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    pFinal
    {
        $p;
        tolerance       1e-06;
        relTol          0;
    }

    Phi
    {
        $p;
    }

    "(Yi|O2|N2|H2O)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0;
    }
    h
    {
        $Yi;
        relTol          0.1;
    }
    hFinal
    {
        $Yi;
    }
}

potentialFlow
{
    // Used for potentialFoam initialisation
    nNonOrthogonalCorrectors 5;
}

PIMPLE
{
    transonic       no;
    nOuterCorrectors 1;
    nCorrectors     3;
    nNonOrthogonalCorrectors 0;
    momentumPredictor yes;
}

relaxationFactors
{
    equations
    {
        ".*"            0.6;
    }
}


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

Best wishes.

[bijan darbari]

hello

I understand that if uncoupled the particles from continuous phase ( set "coupled false;" in "reactingcloud1properties" ) the solution run normally without any divergence.

but It seems that particles aggregated on cell faces when the "brownianmotionforce" activated,

I attached two pictures. in the first picture brownian force activated and in the second, brownian force deactivated. in both cases drag force activated.







what is the problem???
Is it a bug??

[chpjz0391]

Originally Posted by bijan darbari

hello

I understand that if uncoupled the particles from continuous phase ( set "coupled false;" in "reactingcloud1properties" ) the solution run normally without any divergence.

but It seems that particles aggregated on cell faces when the "brownianmotionforce" activated,

I attached two pictures. in the first picture brownian force activated and in the second, brownian force deactivated. in both cases drag force activated.







what is the problem???
Is it a bug??

Hi Bijan,

Have you found the answer to your question?
Why that happened? why the particles aggregated on cell faces ?

Regards,

Ping

[bijan darbari]

Dear chpjz0391

yes. I solved the problem.

Reduce solution time step in control dict ( may be 1e-6 or 1e-8 ). If the problem don't solved, then type "maxCo 0.05;" at the beginning of "kinematiccloud1properties"( or "reactingcloud1properties" ).

[chpjz0391]

Originally Posted by bijan darbari

Dear chpjz0391

yes. I solved the problem.

Reduce solution time step in control dict ( may be 1e-6 or 1e-8 ). If the problem don't solved, then type "maxCo 0.05;" at the beginning of "kinematiccloud1properties"( or "reactingcloud1properties" ).

Hi bijian,

Thanks for your reply. It's quite useful. I have another question, do you need to plot the concentration of particles(mg/m3). As I know, this could be achieved by void fraction. I mean void fraction multiply density of particles. However, I can only get the total concentration of particles in this way. If I inject different sized particles in my simulation, how can I get the concentration for each sized particles?

Kind Regards,

Ping

[bijan darbari]

Dear chpjz0391

Unfortunately I Don't know how to export concentration for each size of particles.
sorry.

[chpjz0391]

Originally Posted by bijan darbari

Dear chpjz0391

Unfortunately I Don't know how to export concentration for each size of particles.
sorry.

Dear bijian,

It's ok. Actually I am using the simpleReactingParcel to simulate the lagrangian particles for steady-state. I met the same problems as yours. Particles aggregated on cell faces. Do you have any suggestion?

Kind Regards,

Ping

[chpjz0391]

Originally Posted by chpjz0391

Dear bijian,

It's ok. Last time you told me ,reactingParcelFoam can simulation collisions between particles. But when I try to add collisionModel in reactingCloud1Properties. It does not work.

Hi bijian,

I have found that function, stochasticCollisionModel is used to simulate collision between particles.

Cheers.

Ping