CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Turbulence inflow generation - recycling method

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

Like Tree1Likes
  • 1 Post By eugene

Reply
 
LinkBack Thread Tools Display Modes
Old   June 12, 2010, 05:11
Default Turbulence inflow generation - recycling method
  #1
Senior Member
 
Jiang
Join Date: Oct 2009
Location: Japan
Posts: 186
Rep Power: 6
panda60 is on a distinguished road
Dear all,
I create this thread to discuss about recycling method to generate inflow turbulence. Now OpenFOAM only has "directMapped" method. That means the internal plane value is direct put back to inlet. This method may be fit for ground vehicles, but not appropriate in wind engineerings. In most cases of wind engineering, we not only need inflow turbulence, the inflow turbulence should also main a predetermined profile.
According to Lund's method 1998, Japanese researcher Kataoka found a method, making assumption that: in the driver region, the shear velocity and boundary layer thickness is not change, and we can get the wanted mean profile in recycling position. The following is Kataoka's method.

Because I am not familiar with foam. So if anyone can give me some suggestions ? And I hope others can discuss here, too.
Attached Images
File Type: jpg mesh.jpg (65.4 KB, 339 views)
File Type: jpg kataoka formula.jpg (40.2 KB, 236 views)
panda60 is offline   Reply With Quote

Old   June 12, 2010, 05:24
Default
  #2
Senior Member
 
Jiang
Join Date: Oct 2009
Location: Japan
Posts: 186
Rep Power: 6
panda60 is on a distinguished road
If I use Kataoka's formulation, if the "directMappedFixedValue" is the best utility we can modify ?

in directMappedFixedValueFvPatchField.C, we can find:

case directMappedPatchBase::NEARESTCELL:
{
if (mpp.sameRegion())
{
newValues = this->internalField();
}
else
{
newValues = nbrMesh.lookupObject<fieldType>
(
fldName
).internalField();
}
mapDistribute::distribute
(
Pstream::defaultCommsType,
distMap.schedule(),
distMap.constructSize(),
distMap.subMap(),
distMap.constructMap(),
newValues
);
break;
}
case directMappedPatchBase::NEARESTFACE:
{
Field<Type> allValues(nbrMesh.nFaces(), pTraits<Type>::zero);
const fieldType& nbrField = nbrMesh.lookupObject<fieldType>
(
fldName
);
forAll(nbrField.boundaryField(), patchI)
{
const fvPatchField<Type>& pf =
nbrField.boundaryField()[patchI];
label faceStart = pf.patch().patch().start();
forAll(pf, faceI)
{
allValues[faceStart++] = pf[faceI];
}
}
mapDistribute::distribute
(
Pstream::defaultCommsType,
distMap.schedule(),
distMap.constructSize(),
distMap.subMap(),
distMap.constructMap(),
allValues
);
newValues = this->patch().patchSlice(allValues);
break;
}
default:
{
FatalErrorIn
(
"directMappedFixedValueFvPatchField<Type>::updateC oeffs()"
)<< "Unknown sampling mode: " << mpp.mode()
<< nl << abort(FatalError);
}
}

in kataoka's formulaiton, u, v and w velocity need different recycling method, and we also need <u> ,<v> and <w> in recycling position, how can I call mean velocity ? how can I manipulate velocity component here ?
panda60 is offline   Reply With Quote

Old   June 16, 2010, 07:33
Default
  #3
Senior Member
 
Eugene de Villiers
Join Date: Mar 2009
Posts: 724
Rep Power: 11
eugene is on a distinguished road
Hi,

I would start from directMappedVelocityFluxFixedValueFvPatchField.

I suggest you create your own average of the mapped velocity inside the modified boundary condition. Depending on an external Umean in this instance would be unreliable.

You can then use the "newValues" variable to update the average and to derive the actual velocity on the patch.
eugene is offline   Reply With Quote

Old   June 17, 2010, 02:19
Default
  #4
Senior Member
 
Jiang
Join Date: Oct 2009
Location: Japan
Posts: 186
Rep Power: 6
panda60 is on a distinguished road
Dear eugene,
I am very grateful you come here.
Yes ,I think directMappedVelocityFluxFixedValueFvPatchField will be the best one to modify.
I think you are right, I should create my own averaged velocity, and then can be easily used.

I also agree with you "Depending on an external Umean in this instance would be unreliable". But I can only do like this, no other method in our wind engineering field.
The following is a sketch map of wind tunnel, as you can see, in the inlet of wind tunnel is a big fan, can generate uniform velocity, then followed by a lot of roughness element put in ground , to generate the boundary layer profile which is nearly same as atmospheric boundary layer.
In our experiment, we measure the velocity profile in black triangle position, we put this profile as the inflow conditon for CFD simulation, then we can compare the CFD result with experiments.

So in most of our wind engineering field, this mean velocity profile is know before CFD simulation, and our generated inflow turbulence must maintain this mean profile, this is the most important thing for us. So kataoka's formulation is the best formulation in wind engineering field now, and is used widely and widely.

My professor said that not so many people use OpenFOAM in wind engineering now, so I think that is why not so many people are interesting in Kataoka's inflow formulation here. But as I know, more and more people begin to use this now.

I have two questions to ask:
first, becaue I want to use a pre-determined mean velocity profile, I should know the face centre of every face, How can I get face centre coordinate z here ?

vectorField allUValues(nbrMesh.nFaces(), vector::zero);
scalarField allPhiValues(nbrMesh.nFaces(), 0.0);
forAll(UField.boundaryField(), patchI)
{
const fvPatchVectorField& Upf = UField.boundaryField()[patchI];
const scalarField& phipf = phiField.boundaryField()[patchI];
label faceStart = Upf.patch().patch().start();
forAll(Upf, faceI)
{
allUValues[faceStart++] = Upf[faceI];
allPhiValues[faceStart] = phipf[faceI];
}
}

Second, how can I compile this new boundary condition ? compile to utility or lib , what should I do ?

Thank you very much.
Attached Images
File Type: jpg wind tunnel.jpg (62.0 KB, 188 views)
panda60 is offline   Reply With Quote

Old   June 18, 2010, 12:52
Default
  #5
Senior Member
 
n/a
Join Date: Sep 2009
Posts: 198
Rep Power: 6
deji is on a distinguished road
Hey there Jiang Guoyi. I am working on something similar as in I have mean experimental data for velocity and temperature that I feed at the inlet of a turbulent natural convection boundary layer. My issue is that the turbulent boundary layer is spatially developing, the displacment grows along the vertical plate. Is your flow fully developed turbulent flow, as in does the streamwise velocity stay unchanged or is yours a developing turbulent boundary layer. I mention this because I have been looking into using the directMapped tool as well.

So how does your flow evolve spatially? or it does not?
deji is offline   Reply With Quote

Old   June 19, 2010, 02:06
Default
  #6
Senior Member
 
Eugene de Villiers
Join Date: Mar 2009
Posts: 724
Rep Power: 11
eugene is on a distinguished road
I think you misunderstood me. There are 2 averaged velocities in the equation you provided. The first is the experimentally measure one <U>(z_inlet), the second is the average you should calculate internally inside the boundary condition <U>(yz_recyc).

I suggest you sepcify the experimental <U>(z_inlet) using the interpolation table format found in boundaries like timeVaryingUniformFixedValue. <U>(yz_recyc) on the other hand has to be found from an average of the mapped value. So what you need to do, is add another vector field to the boundary condition, (lets call it vectorField UMean_mapped). You should calculate UMean_patch as a running average of the mapped value - UMean_mapped = (N-1)/N * UMean_mapped + 1/N*U_mapped - where N is the number of iterations. You will have to store both N and UMean_patch as private member data of the boundary condition. .

The face centre inside a boundary condition is found from:
const vectorField& patchFaceCentres = patch().Cf();

I suggest you compile the boundary condition directly into your application. Just put the *.C and *.H files into the same directory as your application and add the *.C file to the top of Make/files.

I also suggest that you add a keyword like "up (0 0 1);" to specify which direction is "up". Having z hard-coded as the up direction is not good practice.

If possible, I would like a copy of this boundary once you are done.

Good luck,

Eugene
mgg likes this.
eugene is offline   Reply With Quote

Old   June 23, 2010, 10:09
Default
  #7
Senior Member
 
Jiang
Join Date: Oct 2009
Location: Japan
Posts: 186
Rep Power: 6
panda60 is on a distinguished road
Dear eugene,
Thanks, your suggestions are very helpfull to me . It is sure, If you want a copy of this, and if I can finish it. Now I have already try part of these. Because I don't have so much knowledge about Foam, so maybe I can only write some not code which is not so smart.

Could you have a look at if my code has problem or not. Because I still don't know how to store mean velocity, so I just give a value (Up 0 0 ) now.

case directMappedPolyPatch::NEARESTFACE:
{
vectorField allUValues(nbrMesh.nFaces(), vector::zero);

forAll(UField.boundaryField(), patchI)
{
const fvPatchVectorField& Upf = UField.boundaryField()[patchI];

label faceStart = Upf.patch().patch().start();
const vectorField& faceCentres = Upf.patch().Cf(); // add by myself

scalar delta = 0.6; // boundary layer thickness
forAll(Upf, faceI)
{
scalar z = faceCentres[faceI].z();
scalar thita = z/delta;
// damping function, which impedes the growth of boundary layer
scalar phithita = 0.5*(1-tanh(8.0*(thita-1.0)/(0.7-0.4*(thita-0.3)))/tanh(8.0));
// predetermined Ux mean velocity profile, mostly come from experimental measurement
scalar Up = -595.5*z*z*z*z*z*z*z*z+2473.0*z*z*z*z*z*z*z-4372.4*z*z*z*z*z*z
+4299.1*z*z*z*z*z-2559.5*z*z*z*z+936.9*z*z*z-210.2*z*z+32.4*z+1.6;

scalar Ux = Up+phithita*(Upf[faceI].x()-Up);
scalar Uy = phithita*(Upf[faceI].y()-0);
scalar Uz = phithita*(Upf[faceI].z()-0);

allUValues[faceStart++] = vector(Ux, Uy, Uz);
}
}

mapDistribute::distribute
(
Pstream::defaultCommsType,
distMap.schedule(),
distMap.constructSize(),
distMap.subMap(),
distMap.constructMap(),
allUValues
);
newUValues = patch().patchSlice(allUValues);

}

These code has pass the compile, but I don't know if it is right or not.
Maybe I can define a volVectorField Umean, and average in the whole field, maybe this is not good.
panda60 is offline   Reply With Quote

Old   June 23, 2010, 10:19
Default
  #8
Senior Member
 
Jiang
Join Date: Oct 2009
Location: Japan
Posts: 186
Rep Power: 6
panda60 is on a distinguished road
And I also have trouble to use "directMappedVelocityFlux", I want to know if it is the same with "directMapped" ?

I have tried the example in OpenFOAM/run/incompressible/pisoFoam/les/pitzDailyDirectMapped .

The method "nearestCell" works fine, but the method "nearestFace" doesn't work ,it has the following message:

Starting time loop

Reading/calculating field UMean

Reading/calculating field pMean

Reading/calculating field UPrime2Mean

Reading/calculating field pPrime2Mean

fieldAverage: starting averaging at time 0

Time = 1e-05

Courant Number mean: 0 max: 0.2


start 27238 out of range 0 ... 29#0 Foam::error:rintStack(Foam::Ostream&) at ~/OpenFOAM/OpenFOAM-1.6/src/OSspecific/POSIX/printStack.C:203
#1 Foam::error::abort() at ~/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/error.C:242
#2 Foam::Ostream& Foam:perator<< <Foam::error>(Foam::Ostream&, Foam::errorManip<Foam::error>) at ~/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/errorManip.H:87
#3 Foam::UList<Foam::Vector<double> >::checkStart(int) const at ~/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/UListI.H:82
#4 Foam::SubList<Foam::Vector<double> >::SubList(Foam::UList<Foam::Vector<double> > const&, int, int) at ~/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/SubListI.H:61
#5 Foam::List<Foam::Vector<double> >::subList const Foam::fvPatch:atchSlice<Foam::Vector<double> >(Foam::List<Foam::Vector<double> > const&) const at ~/OpenFOAM/OpenFOAM-1.6/src/finiteVolume/lnInclude/fvPatch.H:191
#6 Foam::directMappedFixedValueFvPatchField<Foam::Vec tor<double> >::updateCoeffs() at ~/OpenFOAM/OpenFOAM-1.6/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C:267
#7 Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>::GeometricBoundaryField::updateCoef fs() at ~/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/GeometricBoundaryField.C:271
#8 Foam::fvMatrix<Foam::Vector<double> >::fvMatrix(Foam::GeometricField<Foam::Vector<doub le>, Foam::fvPatchField, Foam::volMesh>&, Foam::dimensionSet const&) at ~/OpenFOAM/OpenFOAM-1.6/src/finiteVolume/lnInclude/fvMatrix.C:230
#9 Foam::fv::gaussLaplacianScheme<Foam::Vector<double >, double>::fvmLaplacianUncorrected(Foam::GeometricFi eld<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) in "/home/panda60/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPDebug/libfiniteVolume.so"
#10 Foam::fv::gaussLaplacianScheme<Foam::Vector<double >, double>::fvmLaplacian(Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) at ~/OpenFOAM/OpenFOAM-1.6/src/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianSchemes.C:114
#11 Foam::fv::laplacianScheme<Foam::Vector<double>, double>::fvmLaplacian(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) at ~/OpenFOAM/OpenFOAM-1.6/src/finiteVolume/lnInclude/laplacianScheme.C:108
#12 Foam::tmp<Foam::fvMatrix<Foam::Vector<double> > > Foam::fvm::laplacian<Foam::Vector<double>, double>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&, Foam::word const&) at ~/OpenFOAM/OpenFOAM-1.6/src/finiteVolume/lnInclude/fvmLaplacian.C:220
#13 Foam::tmp<Foam::fvMatrix<Foam::Vector<double> > > Foam::fvm::laplacian<Foam::Vector<double>, double>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) at ~/OpenFOAM/OpenFOAM-1.6/src/finiteVolume/lnInclude/fvmLaplacian.C:252
#14 Foam::tmp<Foam::fvMatrix<Foam::Vector<double> > > Foam::fvm::laplacian<Foam::Vector<double>, double>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) at ~/OpenFOAM/OpenFOAM-1.6/src/finiteVolume/lnInclude/fvmLaplacian.C:264
#15 Foam::incompressible::LESModels::GenEddyVisc::divD evBeff(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) const at ~/OpenFOAM/OpenFOAM-1.6/src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.C:95
#16 Foam::incompressible::LESModel::divDevReff(Foam::G eometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) const at ~/OpenFOAM/OpenFOAM-1.6/src/turbulenceModels/incompressible/LES/LESModel/LESModel.H:240
#17 main at ~/OpenFOAM/OpenFOAM-1.6/applications/solvers/incompressible/pisoFoam/pisoFoam.C:70
#18 __libc_start_main in "/lib/libc.so.6"
#19 _start at /usr/src/packages/BUILD/glibc-2.9/csu/../sysdeps/i386/elf/start.S:122


From function UList<T>::checkStart(const label)
in file /home/panda60/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude/UListI.H at line 78.

FOAM aborting

I don't know the reason , Could anyone help me .
panda60 is offline   Reply With Quote

Old   June 25, 2010, 05:12
Default
  #9
Senior Member
 
Eugene de Villiers
Join Date: Mar 2009
Posts: 724
Rep Power: 11
eugene is on a distinguished road
To average the mapped field on the inlet patch, first define 2 new variables in the boundary condition .H file:

label N_;
vecterField UpatchMean_;

Then in the evaluation section you would add:

N_++;

UpatchMean_ = (N_ - 1)/N_ * UpatchMean + 1/N_*Umapped

To get the deviatoric velocity:

Udev = Umapped - UpatchMean_;

To set the new inlet velocity:

this->operator==(Up_specifiedMeanVelocityAtInlet + scalingFactor*Udev);

Of course you still have to calculate the scaling factor and specify the mean inlet velocity. Also, remember to write out the mean velocity and number of averaging steps in the write function and read these values in the dictionary constructor if they are available.

I don't know why the nearest face mapping isn't working.
eugene is offline   Reply With Quote

Old   June 26, 2010, 09:29
Default
  #10
Senior Member
 
Jiang
Join Date: Oct 2009
Location: Japan
Posts: 186
Rep Power: 6
panda60 is on a distinguished road
Dear Eugene,
Thanks,
Could you try this example OpenFOAM/run/incompressible/pisoFoam/les/pitzDailyDirectMapped using "nearestFace" method in your computer ? I am afraid I have make mistake somwwhere.
In this tutorials, I just change the nearestCell to nearestFace, but can run. The following is my boundary file in polyMesh directory.

5
(
inlet
{
//type patch;
//nFaces 30;
//startFace 27238;
type directMappedPatch;
nFaces 30;
startFace 27238;
//sampleMode nearestCell;
sampleMode nearestFace;
sampleRegion region0;
samplePatch none;
offset ( 0.04925005 0 0 );
}
outlet
{
type patch;
nFaces 57;
startFace 27268;
}
upperWall
{
type wall;
nFaces 275;
startFace 27325;
}
lowerWall
{
type wall;
nFaces 302;
startFace 27600;
}
frontAndBack
{
type empty;
nFaces 27570;
startFace 27902;
}
)
panda60 is offline   Reply With Quote

Old   June 26, 2010, 09:32
Default
  #11
Senior Member
 
Jiang
Join Date: Oct 2009
Location: Japan
Posts: 186
Rep Power: 6
panda60 is on a distinguished road
And if the nearestFace can't run at last, I think modifing "directMapped" using nearCell is the only way I can do.
panda60 is offline   Reply With Quote

Old   June 27, 2010, 08:47
Default
  #12
Senior Member
 
Jiang
Join Date: Oct 2009
Location: Japan
Posts: 186
Rep Power: 6
panda60 is on a distinguished road
Great, Today I see OpenFOAM-1.7 has been released, I will install and try in new version. And my professor just buy 3 new very good computer for me, next weak I will get these computers ,and Installation in these computers will use a lot of time.
panda60 is offline   Reply With Quote

Old   July 8, 2010, 01:50
Default
  #13
Senior Member
 
Jiang
Join Date: Oct 2009
Location: Japan
Posts: 186
Rep Power: 6
panda60 is on a distinguished road
Dear all,
I am trying to get turbulent flow using directMapped, nearestCell.
But after 10000 steps , I nearly can't see the turbulence.
My domain is 2m long. my recycling position is in the middle place.
mesh size dx=15mm. top and two sides using symmetry plane.
And I am using standard Smagorinsky model.

The following is my x velocity contours, but I nearly can't see turbulence.
I think my domain is difficult to achieve turbulent.
Anyone who has used directMapped can give me some suggestions ?
Attached Images
File Type: jpg recy.jpg (47.9 KB, 151 views)
panda60 is offline   Reply With Quote

Old   July 8, 2010, 08:32
Default
  #14
Senior Member
 
Eugene de Villiers
Join Date: Mar 2009
Posts: 724
Rep Power: 11
eugene is on a distinguished road
If you run a periodic channel flow with exactly the same numerical settings, do you get the same problem? Did you initialise the turbulence in a sensible way? (search for perturbU tool for a way of starting turbulence)
eugene is offline   Reply With Quote

Old   May 8, 2011, 02:36
Default
  #15
New Member
 
Perry L. Johnson
Join Date: Feb 2011
Location: Orlando, FL, USA
Posts: 17
Rep Power: 5
PerryLJohnson is on a distinguished road
Quote:
Originally Posted by eugene View Post
To average the mapped field on the inlet patch, first define 2 new variables in the boundary condition .H file:

label N_;
vecterField UpatchMean_;

Then in the evaluation section you would add:

N_++;

UpatchMean_ = (N_ - 1)/N_ * UpatchMean + 1/N_*Umapped

To get the deviatoric velocity:

Udev = Umapped - UpatchMean_;

To set the new inlet velocity:

this->operator==(Up_specifiedMeanVelocityAtInlet + scalingFactor*Udev);

Of course you still have to calculate the scaling factor and specify the mean inlet velocity. Also, remember to write out the mean velocity and number of averaging steps in the write function and read these values in the dictionary constructor if they are available.

I don't know why the nearest face mapping isn't working.

Quick note: does not N_ in this case need to be a floating point type?
Example: N_ = 2; (N-1)/N = 1/N = 1/2 = 0 if integer type number
PerryLJohnson is offline   Reply With Quote

Old   April 25, 2013, 01:34
Default
  #16
Pj.
Member
 
Luca
Join Date: Mar 2013
Posts: 49
Rep Power: 3
Pj. is on a distinguished road
Hi, I'm trying to use the Kataoka's method too in OF and i found this topic.

Since i see it's a bit old, is there any news? Panda, have you succeeded in the code writing?

Thank you,
Luca
Pj. is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
cartesian grid generation method Abu Taleb Main CFD Forum 7 April 14, 2001 09:49
Cartesian grid generation method Abu Taleb Main CFD Forum 0 April 8, 2001 12:03
Latest News in Mesh Generation Robert Schneiders Main CFD Forum 1 February 18, 2000 00:48
Latest new in mesh generation Robert Schneiders Main CFD Forum 0 February 16, 2000 07:12
Latest news in mesh generation Robert Schneiders Main CFD Forum 0 March 2, 1999 04:07


All times are GMT -4. The time now is 08:42.