[Origami]

I am modelling complex terrain with varying elevations. I want to apply atmBoundaryLayerInletVelocity (see https://cpp.openfoam.org/v3/a00073.html) to the inlet.
Please see below. zGround and Zref are defined with reference to the origin of z axis, which is applicable for a flat terrain situation. does anyone know how to apply atmBoundaryLayerInletVelocity terrain with varying elevations?

inlet
{
type atmBoundaryLayerInletVelocity;
Uref 10.0;
Zref 500;
zDir (0 0 1);
flowDir (1 0 0);
z0 uniform 0.1;
zGround uniform 0.0;
}

[Tobi]

Yes,

out of the box it is not possible. I might investigate here in future and patch the stuff for individual terrains (e.g., inlet is not flat). We had this question 2 weeks ago and there is a workaround.

Specifying a prescribed velocity inlet

[Origami]

Hi Tobi,
Many thanks. Much appreciated for your reply. I will have a look at the link you are referring.
Cheers,
Hongbin

[Tobi]

Could you provide me your case?

[Origami]

Hi Tobi,
Sorry for responding to you this late as I have been very busy lately testing things.
I am using turbineSiting case coming with OpenFoam 8 for testing purpose.
The method you suggested as workaround worked. but the downside is that for running multiple processors (e.g. 8) I had to produce ununiform list for each processor for the inlet files. It is not practical for running large amount of modeling scenarios.

Therefore, I am thinking about using codedFixedValue in inlet files applying equations in the link. https://cpp.openfoam.org/v3/a00073.html

I managed to write a code below for U. the code works for flat terrain. the code basically, assign velocity (x direction) at cell centres at inlet patch. z direction is height. zg is assumed 0.

The next step is to add a code to calculate the height of the cell centres from terrain elevation at the boundary line intersecting with the inlet. I wrote a code using foamNewApp setting up folders. The logic in the code is that i loop though mesh nodes on terrain patch and put those also on the inlet to a List. so for a cell centre with y, I can find the corresponding zg by interpolating the points in the list. a very linear two points interpolation is used. my c++ skill is very rusty now. not sure if there is c++ library which can achieve this easily.

I attach the code. I need to migrate the code to the code below. the trouble is when i try to run simpleFoam. error pop up. mesh class is not recognised.

any idea how to migrate the attached code into the code below?

Many thanks,

Hongbin


inlet
{
type codedFixedValue;
value $internalField;
name codedInletBC;
code
#{
label id_ter = patch().boundaryMesh().findPatchID("terrain");
const polyPatch& ter = patch().boundaryMesh()[id].patch();
const fvPatch& patch = this->patch();
const vectorField& cf = patch.Cf();
scalar x=0, y=0, z=0, ustar=0.1, zg=0,k=0.41;
const scalar z0=0.2, uref = 3.4, zref = 10;
ustar = uref /log((zref+z0)/z0);
vectorField& field = *this;
forAll(cf, i){
x=cf[i].x();
y=cf[i].y();
z=cf[i].z();
field[i] =vector((ustar/k*log((z-zg+z0)/z0)),0,0);
}
#};
}

[Origami]

another thought to write a new boundary condition which inherit atmBoundaryLayer class. I need to rekindle my c++ knowledge. this is a pain for me right now.
https://cpp.openfoam.org/v3/a00073.html

i found a note how to write a new BC>
https://www.dam.brown.edu/people/ntrask/CFD/bctut.pdf

[Tobi]

Originally Posted by Origami

Hi Tobi,
The next step is to add a code to calculate the height of the cell centres from terrain elevation at the boundary line intersecting with the inlet.

Is it in the code snippet you already provided? If so, we could incorporate it into the real ATM boundary but there need to me more logic as the gravity direction can be x y z (means the ground can be in x-y plane, x-z plane or y-z plane). This needs to be incorporated.

unfortunately, I do have too much work at the moment to investigate into that particaluar problem. But it is on my todo list.

[Origami]

Thanks Tobi for your response. Yes. the calculation is in the code snippet I provided. I tested with turbineSiting case. it worked, although i have to admit my c++ skills is quite basic. The code assumes x is the wind direction and y is the cross-wind direction and z is the direction of height. yes, it would be ideal to consider a generic coordinate system in a similar way as atmBoundaryLayer describes . Many thanks for putting this on your to-do list. I think this inlet wind condition for non-flat terrain would be very useful to many wind engineering situation involving real terrain.

[Origami]

Hi Tobi,
Just to let you know I have add some more lines in the snippet to calculated the difference in height of the centres at the inlet patch to its intersecting line with the terrain. I put centres with adjusted z being the difference in a list called pcfAdj.

I have looked at atmBoundaryLayer class and atmBoundaryLayerInletVelocity class. the U profile is defined in atmBoundaryLayer functon atmBoundaryLayer::U (const vectorField& p)
p in the function is patch cell centres (i think)

U function is used in one of the constructor in
atmBoundaryLayerInletVelocityFvPatchVectorField

and refValue() = U(patch().Cf()); I guess it would might if patch().Cf() is replaced with the list of centres with adjusted z.

anyway, just let you know i have updated the snippet. just ignore my message if you are busy.

Cheers, Hongbin

[Origami]

Hi Tobi,
I've finally got the approach to work after so much try and error although the code is very messy needing a lot of tidying up. Thanks for your attention anyway.
Kind regards, Hongbin

[Tobi]

Can you share your code?

Does it work with all directions?

[Origami]

I am still wrapping up the code. Yes it works for all directions. as i mentioned before, regarding atmBoundaryLayerInletVelocityFvPatchVectorField, I only need to change refValue() = U(patch().Cf()) to something like refValue() = U(patchCentreRevised). no need to write code to handle directions. the same applies to epsilon and k.