CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Variable diffusion coefficient for scalar transport in interFoam

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree5Likes
  • 1 Post By bigphil
  • 4 Post By amwitt

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   September 5, 2013, 21:11
Default Variable diffusion coefficient for scalar transport in interFoam
  #1
New Member
 
Adam
Join Date: Mar 2012
Location: St. Paul, MN
Posts: 12
Rep Power: 14
amwitt is on a distinguished road
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:

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

    dimensionedScalar Dg
    (
        twoPhaseProperties.lookup("Dg")
    );
I then created a volScalarField for D in createFields.H as :

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

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

Last edited by amwitt; September 6, 2013 at 11:58.
amwitt is offline   Reply With Quote

Old   September 6, 2013, 03:59
Default
  #2
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,377
Rep Power: 29
akidess will become famous soon enough
Code:
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?
__________________
*On twitter @akidTwit
*Spend as much time formulating your questions as you expect people to spend on their answer.
akidess is offline   Reply With Quote

Old   September 6, 2013, 10:10
Default
  #3
New Member
 
Adam
Join Date: Mar 2012
Location: St. Paul, MN
Posts: 12
Rep Power: 14
amwitt is on a distinguished road
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.
Attached Images
File Type: jpg time1.jpg (86.0 KB, 143 views)
File Type: jpg time2.jpg (97.1 KB, 74 views)
File Type: jpg time3.jpg (97.3 KB, 65 views)
File Type: jpg time4.jpg (94.0 KB, 52 views)
amwitt is offline   Reply With Quote

Old   September 6, 2013, 11:48
Default
  #4
Super Moderator
 
bigphil's Avatar
 
Philip Cardiff
Join Date: Mar 2009
Location: Dublin, Ireland
Posts: 1,086
Rep Power: 34
bigphil will become famous soon enoughbigphil will become famous soon enough
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:
Code:
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
Z.Q. Niu likes this.
bigphil is offline   Reply With Quote

Old   September 6, 2013, 16:55
Default
  #5
New Member
 
Adam
Join Date: Mar 2012
Location: St. Paul, MN
Posts: 12
Rep Power: 14
amwitt is on a distinguished road
Thank you for the response Philip. Two questions related to your post:

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

Code:
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
Attached Images
File Type: jpg picture1.jpg (35.8 KB, 35 views)
File Type: jpg picture2.jpg (47.8 KB, 41 views)
File Type: jpg picture3.jpg (47.1 KB, 38 views)
File Type: jpg picture4.jpg (33.3 KB, 31 views)
File Type: jpg picture5.jpg (35.7 KB, 29 views)
amwitt is offline   Reply With Quote

Old   September 10, 2013, 20:12
Default
  #6
New Member
 
Adam
Join Date: Mar 2012
Location: St. Paul, MN
Posts: 12
Rep Power: 14
amwitt is on a distinguished road
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:

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

Code:
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
amwitt is offline   Reply With Quote

Reply

Tags
interface, interfoam, mass transfer, scalar transport

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
inlet diffusion at the species transport ahchoo FLUENT 1 May 8, 2017 01:19
diffusion Term for add. Variable Zaktatir CFX 9 July 16, 2011 10:51
Latest git 1.6.x: Crash when using inletOutlet for variable alpha1 in interFoam carsten OpenFOAM Bugs 6 September 23, 2009 10:46
Env variable not set gruber2 OpenFOAM Installation 5 December 30, 2005 05:27
Source terms for additional variable transport eqn Nandini Rohilla CFX 0 February 6, 2004 14:38


All times are GMT -4. The time now is 04:23.