
[Sponsors] 
Variable diffusion coefficient for scalar transport in interFoam 

LinkBack  Thread Tools  Display Modes 
September 5, 2013, 20:11 
Variable diffusion coefficient for scalar transport in interFoam

#1 
New Member
Adam
Join Date: Mar 2012
Location: St. Paul, MN
Posts: 12
Rep Power: 5 
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*(1alpha1) 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") ); 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)) ); 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(); Thanks, Adam Last edited by amwitt; September 6, 2013 at 10:58. 

September 6, 2013, 02:59 

#2 
Senior Member
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 921
Rep Power: 17 
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]; } }
__________________
*On twitter @akidTwit *Spend as much time formulating your questions as you expect people to spend on their answer. *Help define the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oamtechnology 

September 6, 2013, 09:10 

#3 
New Member
Adam
Join Date: Mar 2012
Location: St. Paul, MN
Posts: 12
Rep Power: 5 
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 postprocessing, of if there is a fundamental mistake in the underlying code for the evolution of D. 

September 6, 2013, 10:48 

#4 
Senior Member
Philip Cardiff
Join Date: Mar 2009
Location: Dublin,Ireland
Posts: 572
Rep Power: 19 
Hi Adam,
Just a couple of quick points which might be affecting things:
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(); Philip 

September 6, 2013, 15:55 

#5 
New Member
Adam
Join Date: Mar 2012
Location: St. Paul, MN
Posts: 12
Rep Power: 5 
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]; (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.1E9, 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 

September 10, 2013, 19:12 

#6 
New Member
Adam
Join Date: Mar 2012
Location: St. Paul, MN
Posts: 12
Rep Power: 5 
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)) ); Code:
D = Dl*alpha1+Dg*(1alpha1); Dt = D+(turbulence>nut())/(scalar(0.7)); // set internal D field forAll(D, cellI) { if ( scalar(1)alpha1[cellI] < scalar(1e5)  alpha1[cellI] < scalar(1e5) ) { D[cellI] = Dt[cellI]; } } // set boundary field forAll(D.boundaryField(), patchI) { forAll(D.boundaryField()[patchI], faceI) { if ( (scalar(1)alpha1.boundaryField()[patchI][faceI]) < scalar(1e5)  (alpha1.boundaryField()[patchI][faceI]) < scalar(1e5) ) { 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(); } Adam 

Tags 
interface, interfoam, mass transfer, scalar transport 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
diffusion Term for add. Variable  Zaktatir  CFX  9  July 16, 2011 09:51 
inlet diffusion at the species transport  ahchoo  FLUENT  0  April 29, 2011 03:15 
Latest git 1.6.x: Crash when using inletOutlet for variable alpha1 in interFoam  carsten  OpenFOAM Bugs  6  September 23, 2009 09: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 