Fluids

Fluids

UDF for Atmospheric Boundary Layer (ABL)

    • ahmadzaki
      Subscriber

      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

    • Karthik R
      Administrator

      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

    • Rob
      Ansys Employee

      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. 

    • DrAmine
      Ansys Employee

      Correct

    • ahmadzaki
      Subscriber

      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
      Subscriber

      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
      Subscriber

      hi Abenhadj,


      No, it is incorrect because the results were incorrect!


       


      Thank you


      Ahmad

    • Konstantine Kourbatski
      Ansys Employee

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

    • ahmadzaki
      Subscriber

      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 


    • Konstantine Kourbatski
      Ansys Employee

      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
      Subscriber

      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 

    • DrAmine
      Ansys Employee

      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)

    • DrAmine
      Ansys Employee

      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)

    • ahmadzaki
      Subscriber

      @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)


       

    • DrAmine
      Ansys Employee

      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.

    • EDAG
      Subscriber

      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
      Subscriber

      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.


       


       


       

    • EDAG
      Subscriber

      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
      Subscriber

      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

    • EDAG
      Subscriber

      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
      Subscriber

      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

    • EDAG
      Subscriber

      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
      Subscriber

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


       


      good luck

    • amarti
      Subscriber

      Dear Ahmad,


       
      I  wrote you an email to the above direction, but since it seems to be an institutional university email and this is an old post I don't know if you are still using it so i will paste the main boty of the mail here:
       
      ---
       
      I am also working on this issue and I was experiencing some problems mantaining a constant turbulent kinetic energy profile.
       
      Then I used the code file you put there available to download. For what I, first of all, would like to thank you because after compiling it, it worked perfectly and worked as expected.
       
      I first tried it just with the PROFILES but without SOURCES. Geting still bad results. But then I implemented the SOURCE terms with magnific results. And here come the questions I would like to ask you, if you may have the kindness to answer me.
       


      1. I am trying to understand the equations that lead to those source terms. Do you have any paper (of Dr. Hargreaves, Richards or whoever) that explain those source addition? I am trying to understand where that code lines com from and I still don't find it in Dr. Hargreaves work.

      2. In the forum you say that this code comes from Dr. Hargreaves. Did he himself send it to you? Do you have somehow any direct contact to him?




      If you could help me with any information in those two points I would really be very thankful. If you didn't get any further in the comprehension of this code I would really be open to collaborate and share any new finding or step I do.
       
      Kind regards.
       
      Albert.
Viewing 23 reply threads
  • You must be logged in to reply to this topic.