[amwitt]

Hi all,

I am looking to add a variable diffusion coefficient to an advection diffusion equation I added in interFoam to model mass transport in turbulent free surface flows. The coefficient should be weighted based on the value of alpha1 in each cell, i.e. D = Dl*alpha1+Dg*(1-alpha1) where Dl and Dg are scalars. It should also include the addition of turbulent diffusion if alpha1 does not equal 0 or 1, i.e. not at an interface.

I included the following in the createFields.H header:

dimensionedScalar Dl
    (
        twoPhaseProperties.lookup("Dl")
    );

    dimensionedScalar Dg
    (
        twoPhaseProperties.lookup("Dg")
    );

I then created a volScalarField for D in createFields.H as :

Info<< "Reading field D\n" <<endl;
    volScalarField D
    (
      IOobject
      (
         "D",
         runTime.timeName(),
         mesh,
         IOobject::NO_READ,
         IOobject::AUTO_WRITE
      ),
       (Dl*alpha1+Dg*(scalar(1)-alpha1))
      );

Finally, I added the following to the CEqn.H file, which contains the transport equation:

volScalarField Dt = D+(turbulence->nut())/(scalar(0.7));

forAll(D, cellI)
   {
   if
   (
     alpha1[cellI]==scalar(1)
     ||
     alpha1[cellI]==scalar(0)
   )
      {
       D[cellI] = Dt[cellI];
      }
   }


    fvScalarMatrix CEqn
    (
         fvm::ddt(C)
         +fvm::div(phi,C)
         -fvm::laplacian(fvc::interpolate(D),C)
    );

    CEqn.solve();

Everything compiles and seems to run well. However, when I check the values of D, they aren't matching what I expect. In areas where alpha1=1 or 0, I am getting values of D that include turbulent diffusion. It's not clear to me why this is happening. Additionally, the D values shown in paraview look like they are not progressing in time. It appears they are adding to the existing D field each time step. However, the C values are clearly impacted by diffusion. What I want to do is change D each time step based on alpha1 values in each cell and write the results to the time directory to ensure I'm calculating correctly. Any advice on where I'm going wrong is appreciated.

Thanks,
Adam

[akidess]

volScalarField Dt = D+(turbulence->nut())/(scalar(0.7));

forAll(D, cellI)
   {
   if
   (
     alpha1[cellI]==scalar(1)
     ||
     alpha1[cellI]==scalar(0)
   )
      {
       D[cellI] = Dt[cellI];
      }
   }

What does this piece of code do? For all cells where alpha1 is either 1 or 0, you set D=DT. Earlier you set DT=D+turbD, so your output is exactly what you should be expecting?

[amwitt]

That is how I expect the code to work. However, the D field is not evolving in that manner. I have attached four pictures that show the alpha1, nut and D fields. I started computing the transport of C and D at 10s, and at 10.1s (pic 1) the D field is clearly not calculating correctly. There are a few patches of blue indicating turbulent diffusion, but those areas are not consistent with the alpha1 and nut fields as they should be if the code was working correctly. Advancing in time, the D field has patches of blue added to it but it is not evolving like any of the other variable fields. Small areas of high diffusion are added to the D field from the previous time step, and everywhere else in the domain the D field is static in time.

I'm not sure if the D field is just incorrectly written at each time step, which affects the post-processing, of if there is a fundamental mistake in the underlying code for the evolution of D.

[bigphil]

Hi Adam,

Just a couple of quick points which might be affecting things:

• you should set the values on the boundaryField as well as the internalField;
• comparing floats exactly is not a great idea.

Maybe try something like this:

volScalarField Dt = D+(turbulence->nut())/(scalar(0.7));

// set internal field
forAll(D, cellI)
   {
   if ( mag(alpha1[cellI] - 1) < SMALL
     ||
     mag(alpha1[cellI]) < SMALL )
      {
       D[cellI] = Dt[cellI];
      }
   }

// set boundary field
forAll(D.boundaryField(), patchI)
   {
   forAll(D.boundaryField()[patchI], faceI)
      {
       if ( mag(alpha1.boundaryField()[patchI][faceI] - 1) < SMALL
        ||
        mag(alpha1.boundaryField()[patchI][faceI]) < SMALL )
         {
          D.boundaryField()[patchI][faceI] = Dt.boundaryField()[patchI][faceI];
         }
      }
   }
// update coupled boundaries
D.correctBoundaryConditions();

Hope it helps,
Philip

[amwitt]

Thank you for the response Philip. Two questions related to your post:

(1) Excuse my naievety in C++ but the line

D.boundaryField()[patchI][faceI] = Dt.boundaryField()[patchI][faceI];

leads me to believe I’m calling values of Dt from a specified boundary patch. Since I haven't stated any boundary conditions for Dt, I’m not sure what this means. Could you elaborate?

(2) What is the purpose of correctBoundaryConditions()? It seems from the documentation that this is only relevant for volVectorFields.


I tried implementing your suggestions and while there are some general improvements in D, I am still seeing the strange behavior. For example, the initial D field is shown in picture 1. The D field after the first recorded time step is shown in picture 2. It apears that the very first calculation updates the D field to match closely with the alpha1 field, as expected, but subsequent calculations do not update the D field everywhere. In fact, it appears in general to remain the same as it did after the first calculation, with small additions of high diffusion near the lower wall and inlet. Picture 3 shows additional diffusion near the walls 10 simulated seconds later, but identical values of D the red and green regions shown in Picture 2. Picture 4 shows a rescaled view of D, with D values that appear are at least an order of magnitude greater than any nut values found in the domain. Picture 5 shows an enhanced view of D near the inlet, where there is a buildup of diffusion. These lead me to believe the code either (a) isn’t looping correctly, (b) isn’t calculating updated D values correctly or (c) isn’t implementing boundary conditions appropriately. Something tells me it is a combination of (b) and (c).

For the D field boundary conditions I copied those used for nut, which gave me good results for nut before the scalar transport equation was added. For those, I used a nutWallFunction, uniform 0 for the lower wall and zeroGradient conditions at the top and right outlets. The only change I made was to fix the value at the inlet to Dl, or 2.1E-9, where the alpha1 value is fixed.

Something tells me the boundary conditions are not correct, and the values are somehow not getting updated throughout the domain after the initial calculation. I’m not quite sure how to go about testing these individually. Any advice is appreciated.

Thanks,
Adam

[amwitt]

So I think I figured out what I was doing wrong. There was no update or record of D in the CEqn.H file. This meant D was initialized first through createFields.H, and then each iteration only added to the D field when the alpha inequality was met in a particular cell. To remedy, I changed my D declaration in createFields.H to:

    Info<< "Reading field D\n" <<endl;
    volScalarField D
    (
      IOobject
      (
         "D",
         runTime.timeName(),
         mesh,
         IOobject::NO_READ,
         IOobject::AUTO_WRITE
      ),
       mesh,
       dimensionedScalar("D",dimensionSet(0,2,-1,0,0,0,0),0)
       //(Dl*alpha1+Dg*(scalar(1)-alpha1))
    );

    Info<< "Reading field Dt\n" <<endl;
    volScalarField Dt
    (
      IOobject
      (
         "Dt",
         runTime.timeName(),
         mesh,
         IOobject::NO_READ,
         IOobject::AUTO_WRITE
      ),
        mesh,
        dimensionedScalar("Dt",dimensionSet(0,2,-1,0,0,0,0),0)
        //D+(turbulence->nut())/(scalar(0.7)) 
    );

and changed CEqn.H to correctly compute at each time step:

D =  Dl*alpha1+Dg*(1-alpha1);
Dt = D+(turbulence->nut())/(scalar(0.7));

// set internal D field
forAll(D, cellI)
   {
   if
   (
     scalar(1)-alpha1[cellI] < scalar(1e-5)
     ||
     alpha1[cellI]           < scalar(1e-5)
   )
      {
       D[cellI] = Dt[cellI];
      }
   }

// set boundary field
forAll(D.boundaryField(), patchI)
      {
         forAll(D.boundaryField()[patchI], faceI)
         {
            if ( (scalar(1)-alpha1.boundaryField()[patchI][faceI]) < scalar(1e-5)
            ||
               (alpha1.boundaryField()[patchI][faceI]) < scalar(1e-5) )
            {
             D.boundaryField()[patchI][faceI] = Dt.boundaryField()[patchI][faceI];
            }
         }
      }
/// update coupled boundaries
D.correctBoundaryConditions();

//record updated D and Dt
if(runTime.outputTime())
{
D.write();
Dt.write();
}

This appears to do the trick and give me the proper evolution of D I was hoping for. Additionally I set all BCs to calculated since D is computed from internal and boundary values. Thanks for the advice in this thread.

Adam