[lavdwall]

Hi everybody!

I am looking for a way to introduce particles via the fvOptions semiImplicitSource option. However, this doesn't work directly because fvOptions is not available in the alphaEqn (found in the file twoPhaseSystem.C).
There is however the following piece of code in twoPhaseEulerFoam.C:

            volScalarField contErr1
            (
                fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1)
              - (fvOptions(alpha1, rho1)&rho1)
            );

            volScalarField contErr2
            (
                fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2)
               - (fvOptions(alpha2, rho2)&rho2)
            );

How can I use the fvOptions that are given in these equations?

Kind regards,
Laurien

[GerhardHolzinger]

Originally Posted by lavdwall

Hi everybody!
I am looking for a way to introduce particles via the fvOptions semiImplicitSource option. However, this doesn't work directly because fvOptions is not available in the alphaEqn (found in the file twoPhaseSystem.C).

The solver twoPhaseEulerFoam is, as its name suggests, a two-phase Eulerian solver. The continuous (suspending) phase as well as the dispersed phase (bubbles, droplets or particles) are computationally both treated as continuous phases. Thus, there are no particles.

What fvOptions offers is a mass source to introduce the dispersed phase into the domain. Have a look on the injection tutorial case of the twoPhaseEulerFoam tutorials.

[lavdwall]

That was what I meant, I did know that there are no actual particles in the Eulerian models but I was looking for a way to introduce the dispersed phase, as in the tutorial Thanks for the tip!

In the tutorial the mass source is added in the rhoEqn, at least that is what I understand from the following piece of code in the fvOptions file.

    massSource1
    {
        type            scalarSemiImplicitSource;
        $injector1;

        scalarSemiImplicitSourceCoeffs
        {
            volumeMode      absolute;
            injectionRateSuSp
            {
                thermo:rho.air     (1e-3 0); // kg/s
            }
        }
    }

I don't want to introduce the dispersed phase via an injection of rho, but via an injection of alpha (the volume fraction of dispersed phase). But this doesn't seem to work. Any suggestions?

[GerhardHolzinger]

I get your point. However, in twoPhaseEulerFoam, there is no transport equation solved for rho. There is the alpha equation to compute the volume fraction, the momentum equation to compute the velocities, the pressure equation to resolve the pressure-velocity coupling and the energy equation to solve for the internal energy/enthalpy which in turn determines the temperature.
The density rho is computed from your selected equation of state, e.g. rho = rho(p,T).

Somehow, specifying a source for thermo:rho.air ends up being a source for the phase air. I frequenty run simulations which are based on the injection tutorial of twoPhaseEulerFoam.

[lavdwall]

Thanks for your explanation, I understand now

However, I still have some troubles with the simulations. When I introduce my dispersed phase using the mass source (injection of thermo::rho.particles) in fvOptions, the solver calculates negative temperatures already from the first time step. I want to add my dispersed phase with zero momentum so I don't have to add a momentum source in fvOptions. But how do I tell my solver that the dispersed phase is introduced at a certain temperature? I tried with an injection of thermo::T.particles, or thermo::h.particles, but this doesn't seem to work... Any help?

[GerhardHolzinger]

The injection tutorial of twoPhaseEulerFoam has a source for the internal energy. Depending whether you solve for the internal energy or the enthalpy, you need to specify an appropriate source.

[lavdwall]

Hi!

Thanks for the answer. I still have some problems. In order to speed up the calculation, I am not solving the energy equation because temperatures are not really important for my simulation. The case is running but the values for alpha are not bounded between 0 and 1, meaning that unphysical results are produced.

Below is a small part of my log-file.

Create time

Create mesh for time = 0.1


Reading g
Creating twoPhaseSystem

Selecting thermodynamics package 
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleInternalEnergy;
}

Reading face flux field phi.air
Selecting diameterModel for phase air: constant
Selecting turbulence model type RAS
Selecting RAS turbulence model kEpsilon
kEpsilonCoeffs
{
    Cmu             0.09;
    C1              1.44;
    C2              1.92;
    C3              0;
    sigmak          1;
    sigmaEps        1.3;
}

Selecting thermodynamics package 
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState rhoConst;
    specie          specie;
    energy          sensibleInternalEnergy;
}

Reading face flux field phi.particles
Selecting diameterModel for phase particles: constant
Selecting turbulence model type RAS
Selecting RAS turbulence model kineticTheory
Selecting viscosityModel Gidaspow
Selecting conductivityModel Gidaspow
Selecting radialModel SinclairJackson
Selecting granularPressureModel Lun
Selecting frictionalStressModel JohnsonJackson
kineticTheoryCoeffs
{
    equilibrium     off;
    e               0.9;
    alphaMax        0.63;
    alphaMinFriction 0.5;
    residualAlpha   0.0001;
    viscosityModel  Gidaspow;
    conductivityModel Gidaspow;
    granularPressureModel Lun;
    frictionalStressModel JohnsonJackson;
    radialModel     SinclairJackson;
    JohnsonJacksonCoeffs
    {
        Fr              0.05;
        eta             2;
        p               5;
        phi             30;
    }
}

Selecting default blending method: none
Selecting dragModel for (particles in air): GidaspowErgunWenYu
Selecting swarmCorrection for (particles in air): none
Selecting swarmCorrection for (particles in air): none
Selecting swarmCorrection for (particles in air): none
Selecting heatTransferModel for (particles in air): RanzMarshall
Calculating field DDtU1 and DDtU2

Creating field dpdt

Creating field kinetic energy K

No MRF models present

Creating finite volume options from "system/fvOptions"

Selecting finite volume options model type scalarSemiImplicitSource
    Source: massSource1
    - applying source at time 0.1 for duration 0.15
    - selecting cells using cellSet c0
    - selected 27825 cell(s) with volume 9.04545e-05

Selecting finite volume options model type scalarSemiImplicitSource
    Source: energySource1
    - applying source at time 0.1 for duration 0.15
    - selecting cells using cellSet c0
    - selected 27825 cell(s) with volume 9.04545e-05

Courant Number mean: 0.894556 max: 10.2645
Max Ur Courant Number = 10.2644

PIMPLE: no residual control data found. Calculations will employ 2 corrector loops


Starting time loop

Courant Number mean: 0.894556 max: 10.2645
Max Ur Courant Number = 10.2644
Time = 0.10002

PIMPLE: iteration 1
MULES: Solving for alpha.air
MULES: Solving for alpha.air
smoothSolver:  Solving for alpha.air, Initial residual = 0, Final residual = 0, No Iterations 1
alpha.air volume fraction = 1  Min(alpha.air) = 1  Max(alpha.air) = 1
smoothSolver:  Solving for e.air, Initial residual = 1, Final residual = 1, No Iterations 0
smoothSolver:  Solving for e.particles, Initial residual = 1, Final residual = 1, No Iterations 0
min T.air 300
min T.particles 300
GAMG:  Solving for p, Initial residual = 0.00212471, Final residual = 4.74626e-06, No Iterations 2
PIMPLE: iteration 2
MULES: Solving for alpha.air
MULES: Solving for alpha.air
smoothSolver:  Solving for alpha.air, Initial residual = 0.000914294, Final residual = 3.41217e-11, No Iterations 3
alpha.air volume fraction = 0.999982  Min(alpha.air) = 0.999166  Max(alpha.air) = 1
smoothSolver:  Solving for e.air, Initial residual = 1, Final residual = 1, No Iterations 0
smoothSolver:  Solving for e.particles, Initial residual = 1, Final residual = 1, No Iterations 0
min T.air 300
min T.particles 300
GAMG:  Solving for p, Initial residual = 0.00176224, Final residual = 7.87692e-09, No Iterations 8
smoothSolver:  Solving for epsilon.air, Initial residual = 0.0013402, Final residual = 3.68387e-06, No Iterations 4
bounding epsilon.air, min: -1.71626e+06 max: 5.58129e+07 average: 418414
smoothSolver:  Solving for k.air, Initial residual = 0.000176795, Final residual = 4.68494e-06, No Iterations 3
smoothSolver:  Solving for Theta.particles, Initial residual = 0.98061, Final residual = 0.915433, No Iterations 50
ExecutionTime = 9.1 s


Courant Number mean: 0.894554 max: 10.1871
Max Ur Courant Number = 10.187
Time = 0.10004

PIMPLE: iteration 1
MULES: Solving for alpha.air
MULES: Solving for alpha.air
smoothSolver:  Solving for alpha.air, Initial residual = 0.0517544, Final residual = 7.89374e-10, No Iterations 24
alpha.air volume fraction = 0.999964  Min(alpha.air) = 0.998058  Max(alpha.air) = 1
smoothSolver:  Solving for e.air, Initial residual = 1, Final residual = 1, No Iterations 0
smoothSolver:  Solving for e.particles, Initial residual = 0.999999, Final residual = 0.999999, No Iterations 0
min T.air 300
min T.particles 300
GAMG:  Solving for p, Initial residual = 0.00209044, Final residual = 3.44639e-06, No Iterations 2
PIMPLE: iteration 2
MULES: Solving for alpha.air
MULES: Solving for alpha.air
smoothSolver:  Solving for alpha.air, Initial residual = 0.00804932, Final residual = 2.84439e-10, No Iterations 6
alpha.air volume fraction = 0.999964  Min(alpha.air) = 0.998016  Max(alpha.air) = 1.00002
smoothSolver:  Solving for e.air, Initial residual = 1, Final residual = 1, No Iterations 0
smoothSolver:  Solving for e.particles, Initial residual = 1, Final residual = 1, No Iterations 0
min T.air 300
min T.particles 300
GAMG:  Solving for p, Initial residual = 0.00174703, Final residual = 9.91089e-09, No Iterations 7
smoothSolver:  Solving for epsilon.air, Initial residual = 0.00114242, Final residual = 3.54856e-06, No Iterations 4
smoothSolver:  Solving for k.air, Initial residual = 0.000172632, Final residual = 5.04051e-06, No Iterations 3
smoothSolver:  Solving for Theta.particles, Initial residual = 0.136171, Final residual = 0.0133347, No Iterations 50
ExecutionTime = 14.52 s


Courant Number mean: 0.894521 max: 10.1377
Max Ur Courant Number = 10.1377
Time = 0.10006

PIMPLE: iteration 1
MULES: Solving for alpha.air
MULES: Solving for alpha.air
smoothSolver:  Solving for alpha.air, Initial residual = 0.00836022, Final residual = 8.33284e-10, No Iterations 6
alpha.air volume fraction = 0.999946  Min(alpha.air) = 0.996837  Max(alpha.air) = 1
smoothSolver:  Solving for e.air, Initial residual = 1, Final residual = 1, No Iterations 0
smoothSolver:  Solving for e.particles, Initial residual = 0.999999, Final residual = 0.999999, No Iterations 0
min T.air 300
min T.particles 300
GAMG:  Solving for p, Initial residual = 0.00211449, Final residual = 3.87352e-06, No Iterations 2
PIMPLE: iteration 2
MULES: Solving for alpha.air
MULES: Solving for alpha.air
smoothSolver:  Solving for alpha.air, Initial residual = 0.0108336, Final residual = 4.31534e-10, No Iterations 6
alpha.air volume fraction = 0.999946  Min(alpha.air) = 0.997032  Max(alpha.air) = 1.00002
smoothSolver:  Solving for e.air, Initial residual = 1, Final residual = 1, No Iterations 0
smoothSolver:  Solving for e.particles, Initial residual = 0.999999, Final residual = 0.999999, No Iterations 0
min T.air 300
min T.particles 300
GAMG:  Solving for p, Initial residual = 0.00179842, Final residual = 6.10722e-09, No Iterations 8
smoothSolver:  Solving for epsilon.air, Initial residual = 0.00182793, Final residual = 9.80215e-06, No Iterations 6
bounding epsilon.air, min: -762705 max: 4.41792e+07 average: 418205
smoothSolver:  Solving for k.air, Initial residual = 0.000195834, Final residual = 8.3906e-06, No Iterations 3
smoothSolver:  Solving for Theta.particles, Initial residual = 0.0373695, Final residual = 0.006441, No Iterations 50
ExecutionTime = 19.92 s


Courant Number mean: 0.894453 max: 10.275
Max Ur Courant Number = 10.2749
Time = 0.10008

PIMPLE: iteration 1
MULES: Solving for alpha.air
MULES: Solving for alpha.air
smoothSolver:  Solving for alpha.air, Initial residual = 0.0107024, Final residual = 3.82048e-10, No Iterations 6
alpha.air volume fraction = 0.999929  Min(alpha.air) = 0.995972  Max(alpha.air) = 1
smoothSolver:  Solving for e.air, Initial residual = 1, Final residual = 1, No Iterations 0
smoothSolver:  Solving for e.particles, Initial residual = 0.999999, Final residual = 0.999999, No Iterations 0
min T.air 300
min T.particles 300
GAMG:  Solving for p, Initial residual = 0.00216317, Final residual = 3.32473e-06, No Iterations 2
PIMPLE: iteration 2
MULES: Solving for alpha.air
MULES: Solving for alpha.air
smoothSolver:  Solving for alpha.air, Initial residual = 0.0122464, Final residual = 3.83729e-10, No Iterations 6
alpha.air volume fraction = 0.999929  Min(alpha.air) = 0.99613  Max(alpha.air) = 1.00024
smoothSolver:  Solving for e.air, Initial residual = 1, Final residual = 1, No Iterations 0
smoothSolver:  Solving for e.particles, Initial residual = 0.999999, Final residual = 0.999999, No Iterations 0
min T.air 300
min T.particles 300
GAMG:  Solving for p, Initial residual = 0.00184023, Final residual = 5.20187e-09, No Iterations 8
smoothSolver:  Solving for epsilon.air, Initial residual = 0.0012377, Final residual = 4.01932e-06, No Iterations 4
bounding epsilon.air, min: -1.93643e+06 max: 7.05816e+07 average: 417772
smoothSolver:  Solving for k.air, Initial residual = 0.000181734, Final residual = 5.3518e-06, No Iterations 3
smoothSolver:  Solving for Theta.particles, Initial residual = 0.0203782, Final residual = 0.00309259, No Iterations 50
ExecutionTime = 25.41 s


Courant Number mean: 0.894361 max: 10.2392
Max Ur Courant Number = 10.2392
Time = 0.1001

PIMPLE: iteration 1
MULES: Solving for alpha.air
MULES: Solving for alpha.air
smoothSolver:  Solving for alpha.air, Initial residual = 0.011217, Final residual = 2.80122e-10, No Iterations 6
alpha.air volume fraction = 0.999911  Min(alpha.air) = 0.995283  Max(alpha.air) = 1
smoothSolver:  Solving for e.air, Initial residual = 1, Final residual = 1, No Iterations 0
smoothSolver:  Solving for e.particles, Initial residual = 0.999999, Final residual = 0.999999, No Iterations 0
min T.air 300
min T.particles 300
GAMG:  Solving for p, Initial residual = 0.00221309, Final residual = 3.29251e-06, No Iterations 2
PIMPLE: iteration 2
MULES: Solving for alpha.air
MULES: Solving for alpha.air
smoothSolver:  Solving for alpha.air, Initial residual = 0.0130764, Final residual = 3.23792e-10, No Iterations 6
alpha.air volume fraction = 0.999911  Min(alpha.air) = 0.995308  Max(alpha.air) = 1.00058
smoothSolver:  Solving for e.air, Initial residual = 1, Final residual = 1, No Iterations 0
smoothSolver:  Solving for e.particles, Initial residual = 0.999999, Final residual = 0.999999, No Iterations 0
min T.air 300
min T.particles 300
GAMG:  Solving for p, Initial residual = 0.00189597, Final residual = 6.3676e-09, No Iterations 8
smoothSolver:  Solving for epsilon.air, Initial residual = 0.00165239, Final residual = 8.80676e-06, No Iterations 6
bounding epsilon.air, min: -557530 max: 4.30925e+07 average: 417572
smoothSolver:  Solving for k.air, Initial residual = 0.000191587, Final residual = 8.49417e-06, No Iterations 3
smoothSolver:  Solving for Theta.particles, Initial residual = 0.0137528, Final residual = 0.0016478, No Iterations 50
ExecutionTime = 31.1 s


Courant Number mean: 0.894234 max: 10.2039
Max Ur Courant Number = 10.2038
Time = 0.10012

PIMPLE: iteration 1
MULES: Solving for alpha.air
MULES: Solving for alpha.air
smoothSolver:  Solving for alpha.air, Initial residual = 0.0117368, Final residual = 2.3449e-10, No Iterations 6
alpha.air volume fraction = 0.999893  Min(alpha.air) = 0.994339  Max(alpha.air) = 1
smoothSolver:  Solving for e.air, Initial residual = 1, Final residual = 1, No Iterations 0
smoothSolver:  Solving for e.particles, Initial residual = 0.999999, Final residual = 0.999999, No Iterations 0
min T.air 300
min T.particles 300
GAMG:  Solving for p, Initial residual = 0.00227574, Final residual = 3.40033e-06, No Iterations 2
PIMPLE: iteration 2
MULES: Solving for alpha.air
MULES: Solving for alpha.air
smoothSolver:  Solving for alpha.air, Initial residual = 0.0137336, Final residual = 2.894e-10, No Iterations 6
alpha.air volume fraction = 0.999893  Min(alpha.air) = 0.994004  Max(alpha.air) = 1.00123

I already tried to solve the issue by using other iterpolation schemes, by using a smaller time step, by changing my mass source term, by changing the value of alphaSubCycles, ... but nothing seems to work and the maximum value of alpha.air keeps growing until the calculation eventually crashes.
Have you already encountered this issue?

[GerhardHolzinger]

Do you intend to have a Courant number of 10?

Try a smaller time step, to reach a Courant number smaller than one half, or one.

Also solver output like the following indicates that something unpleasant is going on:

bounding epsilon.air, min: -1.71626e+06 max: 5.58129e+07 average: 418414

[lavdwall]

I tried that already, but the same thing happens: alpha.air becomes larger than one already after a few time steps...