UDF for Atmospheric Boundary Layer (ABL)

  • 671 Views
  • Last Post 04 March 2019
ahmadzaki posted this 26 September 2018

Hello,

I have been trying to simulate a homogenous ABL profile based on Richards and Hoxey's 1993 approach using FLUENT software, which requires a UDF at the inlet. I succeeded in  creating the velocity profile using the log law

Uref= u (log(href/z0)/log(h/z0)) [shown below in code]

Uref: reference velocity at refernce height (href), u: velocity at h (height), and z0: aerodynamic roughness height.

The result was successful as shown:

However, a homogenous turbulence kinetic energy (TKE) and dissipation rate (TDR) were not generated despite the inserted constant values (calculated from the mentioned apprach relative to friction velocity, u*) in the inlet display window as shown (the values were changed from 1 and 1 to differnt numbers of course).  The results changed due to the ground roughness (may be) that might be another possible influence (UDF?).

So I tried making a UDF code for both TKE and TDR as well to keep them constant based on the codes below, but it didn't work as well, the TKE gave high values at the top and minimum (zero) at the bottom:


#include "udf.h"

 #define z0 0.0000824          /* constants */
 #define z1 0.4
 #define uref 14.7
 #define CMU 0.09
 #define VKC 0.41


 /* profile for x-velocity */

 DEFINE_PROFILE(log_velocity,thread,index)
 {
   real y[ND_ND];
   real z2;
   face_t f;

   begin_f_loop(f,thread)
   {
        F_CENTROID(y,f,thread);
    z2 = y[1];
    F_PROFILE(f,thread,index) = uref*log(z2/z0)/log(z1/z0);
    }
 end_f_loop(f,thread)
   }



 /* profile for kinetic energy */
 
 DEFINE_PROFILE(k_profile,thread,index)
 {
    real y[ND_ND], lgr, ustr;
    real k;
    face_t f;

    lgr = log((z1+z0)/z0);
    ustr = (uref*VKC)/lgr;
      

    begin_f_loop(f,thread)
      {
         F_CENTROID(y,f,thread);

         k = pow(ustr,2.)/sqrt(CMU);
         F_PROFILE(f,thread,index) = k;
          
        }
      end_f_loop(f,thread)
   }



/* profile for dissipation rate */
 
 DEFINE_PROFILE(dissip_profile,thread,index)
 {
    real y[ND_ND], lgr, ustr;
    real eps;
    face_t f;

    lgr = log((z1+z0)/z0);
    ustr = (uref*VKC)/lgr;
   

    begin_f_loop(f,thread)
      {
         F_CENTROID(y,f,thread);

         eps = pow(ustr,3.)/(VKC*(z1+z0));
         F_PROFILE(f,thread,index) = eps;

        }
      end_f_loop(f,thread)
 }

 

what should I do?

Kindly notice that, I'm not good coder, in fact I don't know how to write codes, my efforts were done to try understand them, so this code is  a beginner trial. My question is How?, Am I missing something?

and yes, I read the manual.

 

Thank you for support

Ahmad

  • Liked by
  • EDAG
Order By: Standard | Newest | Votes
Kremella posted this 26 September 2018

Hello,

There are two things I like to do when I write my UDFs.

  • Fluent customization manual has multiple UDFs  examples, especially ones which use DEFINE_PROFILE macro. I like to first copy one of these UDFs and start there. I change only the parts of code when I start there and I am fairly confident that my UDF does not have any other syntax errors.
  • In addition to this, I add fprintf comments to my UDF. This would print data on my console window and helps me identify and debug a UDF.

Please let us know what you find. I hope this helps.

Thank you.

Best Regards,

Karthik

rwoolhou posted this 26 September 2018

Review the relevant sections of the UDF: you've not set   DEFINE_PROFILE(k_profile,thread,index) or DEFINE_PROFILE(dissip_profile,thread,index) as a constant. 

abenhadj posted this 26 September 2018

Correct

Best regards,

Amine

ahmadzaki posted this 26 September 2018

Hello Karthik,

Yes, Thats what I've been doing! (as I mentioned). However, my case deviates from the mentioned example of FLUENT UDF user manual.

Thank you

Ahmad

ahmadzaki posted this 26 September 2018

Hi Rwoolhou,

Yes I understand there is a problem, but for a guy who don't know how to code or just trip through try and error, can you be more specific?!

Thank you for support

Ahmad

ahmadzaki posted this 26 September 2018

hi Abenhadj,

No, it is incorrect because the results were incorrect!

 

Thank you

Ahmad

kkourbat posted this 27 September 2018

could you clarify where you see the deviation from prescribed boundary profile distributions? Is it right at the boundary or inside the domain?

ahmadzaki posted this 27 September 2018

Hi, kkourbat,

I have problems with the k and e. they constantly start at the inlet but eventually the change along the domain (inside the domain).. the following figures show the TKE contours (left) where it starts constant(red region) then starts to reduce along the domain (gradient blue regions).

I took three vertical line samples at the centre of the domain one at the inlet (white), middle (green), and outlet (red), then I plotted the TKE values at each line. You can see that the white line (inlet) is constant but as we go inside the domain, middle and output, the TKE reduces.

My objective is to keep both TKE and dissipation rates constant.

Thank you 

  • Liked by
  • EDAG
kkourbat posted this 27 September 2018

I don't understand why you think this is wrong. Fluent solves transport equations for TKE and dissipation rate, and they naturally evolve as the boundary layer develops downstream. They cannot be constant, as otherwise they wouldn't satisfy their governing equations and boundary conditions. You can turn off equations for TKE and dissipation so they are not solved to force constant TKE and dissipation throughout the domain. I'll let you decide how physical this would be for your particular need, as you would effectively be providing constant eddy viscosity which is never a constant in a turbulent boundary layer.  

ahmadzaki posted this 27 September 2018

I understand, however, in wind engineering we try to simulate a constant atmospheric boundary layer to match the actual boundary layer above cities and buildings. We do that, so when a building is placed in the domain we could be able to see its effect regardless of the domain effect since we are simulating large domains at which velocity, TKE, and TDR deteriorate which deviates from reality and wind environments around full-scale buildings. This is unlike domains made for aerofoils, flat plates, and ducts.

Thanks anyway 

abenhadj posted this 27 September 2018

I wrote correct to the reply provided by rwoolhou. He wrote that you need to review your UDF as you are using wrong profiles for dissipation and TKE.

But in your UDF you are then assigning TKE and Epsilon the value of Y-coordinate:

 

 begin_f_loop(f,thread)
      {
         F_CENTROID(y,f,thread);
         k=y[1];
     F_PROFILE(f,thread,index) = k;
           
        }
      end_f_loop(f,thread)

Generally for inlet you will need something like

 

teWind = (uStar^2)/sqrt(cmu=0.09)

 

edWind = (uStar^3)/(0.41*heightprof)

Best regards,

Amine

  • Liked by
  • ahmadzaki
abenhadj posted this 27 September 2018

teWind is the TKE Profile edWind is the dissipation profile. The profile height is the maximum f the of ground roughness and the actual height:

heightprofile=max(ground_roughness, height-height of the ground)

Best regards,

Amine

  • Liked by
  • ahmadzaki
ahmadzaki posted this 27 September 2018

@abenhadj

Thank you for your reply and concern.

I did write those equations but I used different symbols:

 

"k = pow(ustr,2.)/sqrt(CMU);" for TKE

and

"eps = pow(ustr,3.)/(VKC*(z1+z0));", where VKC=0.41

the code gave me the same results as when I define them as constant values in Fluent inlet user window (as shown in the question fig)

 

abenhadj posted this 28 September 2018

Check why the values are high or low at the extreme heights. Perhaps the issue is related to the mean speed itself. As I do not the formula I cannot comment here. ABL conditions are working for me. For the mean speed I use

speed = uStar*log(height_prof/groundRougness)/van karaman constant

 

You have to do some debugging to figure out the reason for the behavior. But Also please reformulate your question: On the beginning you said the profile are not constant and then you used a constant profile which you edited later. In other post you wrote :"My objective is to keep both TKE and dissipation rates constant." If you want to have constant profiles at the inlet then just give constant values here based on averaged wind profile.

Best regards,

Amine

  • Liked by
  • ahmadzaki
EDAG posted this 25 February 2019

Hi ahmadzaki, like you, I am also trying to simulate the horizontal homogeneous ABL, but i'm having problems with the TKE contour plot too. Have you solved this problem yet? Hope you can help me.

 

Best regards

ahmadzaki posted this 26 February 2019

Hello EDAG,

I'm currently trying to use CFX instead, since it is much easier to just implement the ABL equations as expressions, however, I'm having errors and problems that I'm still not able to achieve any results.

I can send you a UDF code, that was based on Richards and Hoxey approach, however, the code doesn't work with me, since its giving me a floating point and exits the solution.

 

 

 

  • Liked by
  • EDAG
EDAG posted this 26 February 2019

That would be very helpful! I was thinking on use CFX too, but I still want to suffer a little bit in Fluent. Please let me know if I can help you with something later.

ahmadzaki posted this 26 February 2019

Here is the code, as I said its based on Richards and Hoxey approach. I tried to modify it and compiled it, however it failed to inialise. If you succeded in understanding the second part of the code, starting from "Define_Source" and knew what to modify so it works with all cases, then please let me know.

 

https://drive.google.com/file/d/1Z0hKICnOkjFsdClfDvLsHYILxTcbiJNX/view?usp=sharing

  • Liked by
  • EDAG
EDAG posted this 28 February 2019

Hi ahmadzaki, I only have two questions about your code. I don't understand what that "HUGE" value is for and where did you find that x-momentum source equation?

I think that the units of the k-ϵ sources that you are using are not consistent with the units that Fluent needs, e.g. for the 'k' source term you are specifying that the units are m2/s2 and according to the source terms units for "k" on Fluent, these units should be in kg/m-s3. I'm not sure if a source term for "k" is really needed.

Besides this, your code has been very useful. Thanks.

 

Best regards

ahmadzaki posted this 28 February 2019

Its not my code, I got it from Dr. Hargreaves, he was to busy to explain it and I tried to figure it out myself. I didn't know HUGE term (added to , thats why I sent you the code.

I think the x-mon source represents the ground shear, or the degredation of TKE due to wall while the flow moving through x-directin (may be I'm wrong).

The k-e are correct based on Richard and Hoxey's equations. I think the whole point of adding the source loops is to maintain them constant, along the domain. I wrote a code without them and the TKE profile looked like this:

 

 

which was satisfying for me, I don't think it will make much a difference in my case (the degradation happens at the end, my building is still in the front region).. I'm not studying ABL, I'm just trying to implement it in my research..

I tried to find the HUGE trem in Hargreaves and Wright (2007) paper, but I couldn't find anything. I guess one of my issues is that I'm trying to implement the ABL codes on a scaled one (wind tunnel) and that causes my soulution to fail.. although I modified the dimensions. Therefore, I think the HUGE term is related to the domain itself, but I have to understand first..

Please let me know of any findings or whether you were able to implement the code to your case or not..I will be grateful to learn.

 

Thank you

Ahmad

  • Liked by
  • EDAG
EDAG posted this 01 March 2019

I would like to communicate with you more frequently so that we can progress more quickly. I'm also trying to just implement ABL in my research, not to study it deeper. My buildings are near the inlet too.

ahmadzaki posted this 04 March 2019

This is my email: azak215@aucklanduni.ac.nz

 

good luck

Close