TurbSim output for Unstable LOG wind profile

I am using the LOG wind profile option in TurbSim for different stability conditions depending on the RICH_NO.I am using TurbSim alpha version. I noticed on the program code that for Stable conditions: psiM = -5.0ZL and the wind speed at Ht can be obtained using U_{ref}( LOG( Ht / Z0 ) - psiM )/( LOG( z_Ref / Z0 ) - psiM ). I have used a calculator and checked it with the output from TurbSim and I got the same results. I have also got similar results for neutral atmospheric conditions.

I was trying to do the same to check the outputs for Unstable condition (I am using RICH_NO = -3.8) but the values I got using calculator were different from TurbSim output .

My understanding is that for unstable conditions, TurbSim uses the Businger-Dyer formula:
Ψ_m (z/L)=ln⁡[((1+x^2)/2) ((1+x)/2)^2 ]-2 tan^(-1)⁡x+π/2 Where x=〖(1-16 z/L)〗^(1/4)

this is written on the program code as:

tmp = (1.0 - 15.0*ZL)**0.25 (****15 is used here instead of 16 in the formula above. but this is not what I am worried about!!)

        !psiM = -2.0*LOG( (1.0 + tmp)/2.0 ) - LOG( (1.0 + tmp*tmp)/2.0 ) + 2.0*ATAN( tmp ) - 0.5 * PI
        psiM = -LOG( 0.125 * ( (1.0 + tmp)**2 * (1.0 + tmp*tmp) ) ) + 2.0*ATAN( tmp ) - 0.5 * PI

why is psiM multiplied by the negative sign on the code above?
most importantly,it would be great if you could point out the reason for not getting similar results using calculator. Let us take these values as an example: Rich_No=-3.8, z/L=-3.8, Z0=0.0001m,z_Ref=21m, Ht=90m, U_{ref}=3m/s and calculate Mean wind speed at hub height. Running TurbSim the value I obtained is 3.310 m/s. using a calculator I obtained 3.0294m/s. My working is shown below.

tmp = (1.0 - 15.0*(-3.8))**0.25 = 2.759669
psiM = -LOG( 0.125 * ( (1.0 + tmp)**2 * (1.0 + tmptmp) ) ) + 2.0ATAN( tmp ) - 0.5 * PI = 135.8691749 (but I used -135.8691749 according to the Businger-Dyer formula)
U_{ref}*( LOG( Ht / Z0 ) - psiM )/( LOG( z_Ref / Z0 ) - psiM ) = 3.0294 m/s

Please help.
With kind regards

For your reference, here is part of the TurbSim input file I used.
--------Turbine/Model Specifications-----------------------
31 NumGrid_Z - Vertical grid-point matrix dimension
31 NumGrid_Y - Horizontal grid-point matrix dimension
0.05 TimeStep - Time step [seconds]
600 AnalysisTime - Length of analysis time series [seconds] (program will add time if necessary: AnalysisTime = MAX(AnalysisTime, UsableTime+GridWidth/MeanHHWS) )
600 UsableTime - Usable length of output time series [seconds] (program will add GridWidth/MeanHHWS seconds unless UsableTime is “ALL”)
90 HubHt - Hub height [m] (should be > 0.5GridHeight)
165 GridHeight - Grid height [m]
165 GridWidth - Grid width [m] (should be >= 2
0 VFlowAng - Vertical mean flow (uptilt) angle [degrees]
0 HFlowAng - Horizontal mean flow (skew) angle [degrees]

--------Meteorological Boundary Conditions-------------------
“NONE” TurbModel - Turbulence model (“IECKAI”,“IECVKM”,“GP_LLJ”,“NWTCUP”,“SMOOTH”,“WF_UPW”,“WF_07D”,“WF_14D”,“TIDAL”,“API”,“USRINP”,“TIMESR”, or “NONE”)
“UsrSpec.inp” UserFile - Name of the file that contains inputs for user-defined spectra or time series inputs (used only for “USRINP” and “TIMESR” models)
“1-ED3” IECstandard - Number of IEC 61400-x standard (x=1,2, or 3 with optional 61400-1 edition number (i.e. “1-Ed2”) )
“B” IECturbc - IEC turbulence characteristic (“A”, “B”, “C” or the turbulence intensity in percent) (“KHTEST” option with NWTCUP model, not used for other models)
“NTM” IEC_WindType - IEC turbulence type (“NTM”=normal, “xETM”=extreme turbulence, “xEWM1”=extreme 1-year wind, “xEWM50”=extreme 50-year wind, where x=wind turbine class 1, 2, or 3)
“default” ETMc - IEC Extreme Turbulence Model “c” parameter [m/s]
“LOG” WindProfileType - Velocity profile type (“LOG”;“PL”=power law;“JET”;“H2L”=Log law for TIDAL model;“API”;“USR”;“TS”;“IEC”=PL on rotor disk, LOG elsewhere; or “default”)
“unused” ProfileFile - Name of the file that contains inputs profiles for WindProfileType=“USR” and/or TurbModel=“USRVKM” [-]
21 RefHt - Height of the reference velocity (URef) [m]
3 URef - Mean (total) velocity at the reference height [m/s] (or “default” for JET velocity profile)
“default” ZJetMax - Jet height [m] (used only for JET velocity profile, valid 70-490 m)
0.14 PLExp - Power law exponent [-] (or “default”)
0.0001 Z0 - Surface roughness length [m] (or “default”)

--------Non-IEC Meteorological Boundary Conditions------------
“default” Latitude - Site latitude [degrees] (or “default”)
-3.800 RICH_NO - Gradient Richardson number [-]
“default” UStar - Friction or shear velocity [m/s] (or “default”)
“default” ZI - Mixing layer depth [m] (or “default”)
“none” PC_UW - Hub mean u’w’ Reynolds stress [m^2/s^2] (or “default” or “none”)
“none” PC_UV - Hub mean u’v’ Reynolds stress [m^2/s^2] (or “default” or “none”)
“none” PC_VW - Hub mean v’w’ Reynolds stress [m^2/s^2] (or “default” or “none”)

Sorry the outputs are correct!!! I made a mistake of not using radian for the Acrtan in : psiM = -LOG( 0.125 * ( (1.0 + tmp)**2 * (1.0 + tmptmp) ) ) + 2.0ATAN( tmp ) - 0.5 * PI


I still have a problem with the sign of psiM that TurbSim code uses for unstable condition. I have checked the formula in the following books:

STULL, R. B. (1988): Boundary Layer Meteorology. Kluwer Academic Publishers (page 385)
Panofsky, H.A.; Dutton, J.A. (1984). Atmospheric Turbulence: Models and Methods for Engineering Applications. New York: Wiley-Interscience; 397 pp.

TurbSim uses : psiM = -LOG( 0.125 * ( (1.0 + tmp)**2 * (1.0 + tmptmp) ) ) + 2.0ATAN( tmp ) - 0.5 * PI
and U_{ref}*( LOG( Ht / Z0 ) - psiM )/( LOG( z_Ref / Z0 ) - psiM )

but according to the text books, psiM = LOG( 0.125 * ( (1.0 + tmp)**2 * (1.0 + tmptmp) ) ) - 2.0ATAN( tmp ) + 0.5 * PI
because negative psiM is already used in the main formula
U_{ref}*( LOG( Ht / Z0 ) - psiM )/( LOG( z_Ref / Z0 ) - psiM )

please, could you clarify this for me.

Hi, Ashenafi.

After taking a look at the code and some references, I think you are correct. PsiM appears to be specified incorrectly when z/L < 0. This part of the code hasn’t changed since it was called SNLWind3D, so the code has implemented this incorrectly for at least 15 years. Thank you for bringing this issue to my attention. I will make sure it gets fixed and updated on the web in the near future.

For others reading this post,
This bug has no effect on results using IEC spectral models. It affects only non-IEC models where Rich_No < 0 (and uses either the diabatic ustar or log profile).

Thank you very much for clearing the confusion. For z/L < 0, I used the user-defined wind profile option in TurbSim v2.0 alpha version. just need to manually calculate PsiM and wind speeds at different height. Every thing went well after that.

Thank you!!!