[saimat]

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:



in which and are positive constants. Given this dependency and assuming (neglecting advection) the advection-diffusion-equation looks as follows:



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

    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:

    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.

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):

    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.

[saimat]

Originally Posted by saimat

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

I figured out, that the implementation in OpenFoam is already suitable for variable coefficients (see Programmers Guide table 2.2 / Chapter 2.4.1). I was just misrouted by the name "laplacian(D,alpha)". Consequently adding an additional term is not necessary.