CFD Online Discussion Forums - Implementing variable Diffusion Coefficient

Dear foamers,

I am trying to expand the functionality of twoLiquidMixingFoam so that it can handle Diffusion coefficients which depend on the volume fraction of one of the liquids:

$D = D(\alpha_1) = C_1 e^{C_2 \alpha_1},$

in which

and

are positive constants. Given this dependency and assuming

(neglecting advection) the advection-diffusion-equation looks as follows:

$\frac{\partial \alpha_1}{\partial t} = \nabla \cdot \left( D(\alpha_1) \nabla{\alpha_1} \right) = D(\alpha_1) \Delta{\alpha_1} + \nabla{D(\alpha_1)} \cdot \nabla{\alpha_1}$

The first term is already implemented in twoLiquidMixingFoam/alphaDiffusionEqn.H, whereas the second one is missing:

Code:

    fvScalarMatrix alpha1Eqn

    (

        fvm::ddt(alpha1)

      - fvc::ddt(alpha1)

      - fvm::laplacian

        (

            volScalarField("Dab", Dab + alphatab*turbulence->nut()),

            alpha1

        )

    );

I tried to implement it in two ways and tried to validate the resulting solver in a simple 1-D configuration initially containing a jump from alpha1=0 to alpha1=1.

The first try was taking the values from the last timestep:

Code:

    Dab = C1*Foam::exp(C2*alpha1);



    volScalarField Sc = fvc::grad(Dab) & fvc::grad(alpha1);



    fvScalarMatrix alpha1Eqn

    (

        fvm::ddt(alpha1)

      - fvc::ddt(alpha1)

      - fvm::laplacian

        (

            volScalarField("Dab", Dab + alphatab*turbulence->nut()),

            alpha1

        )

      ==

        Sc

    );

The results showed a growing total volume fraction of alpha1, meaning a production of phase 1. Additionally, the alpha1 values exceeded 1.

Code:

Starting time loop



...



Courant Number mean: 0 max: 0

Interface Courant Number mean: 0 max: 0

Z number max: 0

deltaT = 0.0119976

Time = 0.0119976



MULES: Solving for alpha1

Phase-1 volume fraction = 0.027  Min(alpha1) = 0  Max(alpha1) = 1

...

ExecutionTime = 0.54 s  ClockTime = 0 s



Courant Number mean: 0 max: 0

Interface Courant Number mean: 0 max: 0

Z number max: 3382.27

deltaT = 7.0944e-06

Time = 0.0120047



MULES: Solving for alpha1

Phase-1 volume fraction = 0.0681488  Min(alpha1) = 0  Max(alpha1) = 15.1184

MULES: Solving for alpha1

...

Both of these effects are dependent on the grid resolution. Nevertheless, even with poor resolution, this non-physical, non-conservative behavour persists.

The second way I tried to implement the additional term was by using Picard's method (https://www.cfd-online.com/Wiki/Sour...inearization):

Code:

    Dab = C1*Foam::exp(C2*alpha1);



    volScalarField Sp = C1*C2*C2*Foam::exp(C2*alpha1) * fvc::grad(alpha1) & fvc::grad(alpha1);

    volScalarField Sc = (C1*C2*Foam::exp(C2*alpha1) * fvc::grad(alpha1) & fvc::grad(alpha1))*(1-C2*alpha1);



    fvScalarMatrix alpha1Eqn

    (

        fvm::ddt(alpha1)

      - fvc::ddt(alpha1)

      - fvm::laplacian

        (

            volScalarField("Dab", Dab + alphatab*turbulence->nut()),

            alpha1

        )

      ==

        fvm::Sp(Sp, alpha1) + Sc

    );

In that case I was not able to get a stable calculation/converging solution. The solver diverged after some timesteps.

Does anyone have an Idea how to improve the implementation? I would be very grateful if anyone could give me a hint how to solve the problem.