[panda60]

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.

[panda60]

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 ?

[eugene]

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.

[panda60]

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.

[deji]

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?

[eugene]

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

[panda60]

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]

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 .

[eugene]

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.

[panda60]

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]

And if the nearestFace can't run at last, I think modifing "directMapped" using nearCell is the only way I can do.

[panda60]

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]

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 ?

[eugene]

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)

[PerryLJohnson]

Originally Posted by eugene

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

[Pj.]

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