JSBSim Turbulence Simulation #953

Dobid · 2023-09-01T15:57:14Z

Dobid
Sep 1, 2023

Hello!

I'm currently working on RL based control for fixed wing UAVs under disturbances and I would like to know if someone could explain to me how JSBSim generates turbulences.

Additional details coming from my investigations:
When looking at the code, I see that the scale lengths and amplitudes are implemented as shown in the cited reference (MIL-F-8785C) and in the Continuous gusts wikipedia page.

jsbsim/src/models/atmosphere/FGWinds.cpp

Lines 258 to 280 in e4aa2a9

    
           // Turbulence model according to MIL-F-8785C (Flying Qualities of Piloted Aircraft) 
        
           double b_w = in.wingspan, L_u, L_w, sig_u, sig_w; 
        
             if (b_w == 0.) b_w = 30.; 
        
           // clip height functions at 10 ft 
        
           if (h <= 10.) h = 10; 
        
           // Scale lengths L and amplitudes sigma as function of height 
        
           if (h <= 1000) { 
        
             L_u = h/pow(0.177 + 0.000823*h, 1.2); // MIL-F-8785c, Fig. 10, p. 55 
        
             L_w = h; 
        
             sig_w = 0.1*windspeed_at_20ft; 
        
             sig_u = sig_w/pow(0.177 + 0.000823*h, 0.4); // MIL-F-8785c, Fig. 11, p. 56 
        
           } else if (h <= 2000) { 
        
             // linear interpolation between low altitude and high altitude models 
        
             L_u = L_w = 1000 + (h-1000.)/1000.*750.; 
        
             sig_u = sig_w = 0.1*windspeed_at_20ft 
        
                           + (h-1000.)/1000.*(POE_Table->GetValue(probability_of_exceedence_index, h) - 0.1*windspeed_at_20ft); 
        
           } else { 
        
             L_u = L_w = 1750.; //  MIL-F-8785c, Sec. 3.7.2.1, p. 48 
        
             sig_u = sig_w = POE_Table->GetValue(probability_of_exceedence_index, h); 
        
           }

The next lines are dedicated to declaring variables following the Yeager report's conventions:

jsbsim/src/models/atmosphere/FGWinds.cpp

Lines 292 to 307 in e4aa2a9

    
           double 
        
             T_V = in.totalDeltaT, // for compatibility of nomenclature 
        
             sig_p = 1.9/sqrt(L_w*b_w)*sig_w, // Yeager1998, eq. (8) 
        
             //sig_q = sqrt(M_PI/2/L_w/b_w), // eq. (14) 
        
             //sig_r = sqrt(2*M_PI/3/L_w/b_w), // eq. (17) 
        
             L_p = sqrt(L_w*b_w)/2.6, // eq. (10) 
        
             tau_u = L_u/in.V, // eq. (6) 
        
             tau_w = L_w/in.V, // eq. (3) 
        
             tau_p = L_p/in.V, // eq. (9) 
        
             tau_q = 4*b_w/M_PI/in.V, // eq. (13) 
        
             tau_r =3*b_w/M_PI/in.V, // eq. (17) 
        
             nu_u = generator->GetNormalRandomNumber(), 
        
             nu_v = generator->GetNormalRandomNumber(), 
        
             nu_w = generator->GetNormalRandomNumber(), 
        
             nu_p = generator->GetNormalRandomNumber(), 
        
             xi_u=0, xi_v=0, xi_w=0, xi_p=0, xi_q=0, xi_r=0;

Then follows the implementation of the Tustin model.
Q1: Since we're already in a switch turbType; case ttTustin why do we need a if(turbType == ttTustin) condition ?

jsbsim/src/models/atmosphere/FGWinds.cpp

Lines 311 to 344 in e4aa2a9

    
           if (turbType == ttTustin) { 
        
             // the following is the Tustin formulation of Yeager's report 
        
             double 
        
               omega_w = in.V/L_w, // hidden in nomenclature p. 3 
        
               omega_v = in.V/L_u, // this is defined nowhere 
        
               C_BL  = 1/tau_u/tan(T_V/2/tau_u), // eq. (19) 
        
               C_BLp = 1/tau_p/tan(T_V/2/tau_p), // eq. (22) 
        
               C_BLq = 1/tau_q/tan(T_V/2/tau_q), // eq. (24) 
        
               C_BLr = 1/tau_r/tan(T_V/2/tau_r); // eq. (26) 
        
             // all values calculated so far are strictly positive, except for 
        
             // the random numbers nu_*. This means that in the code below, all 
        
             // divisors are strictly positive, too, and no floating point 
        
             // exception should occur. 
        
             xi_u = -(1 - C_BL*tau_u)/(1 + C_BL*tau_u)*xi_u_km1 
        
                  + sig_u*sqrt(2*tau_u/T_V)/(1 + C_BL*tau_u)*(nu_u + nu_u_km1); // eq. (18) 
        
             xi_v = -2*(sqr(omega_v) - sqr(C_BL))/sqr(omega_v + C_BL)*xi_v_km1 
        
                  - sqr(omega_v - C_BL)/sqr(omega_v + C_BL) * xi_v_km2 
        
                  + sig_u*sqrt(3*omega_v/T_V)/sqr(omega_v + C_BL)*( 
        
                        (C_BL + omega_v/sqrt(3.))*nu_v 
        
                      + 2/sqrt(3.)*omega_v*nu_v_km1 
        
                      + (omega_v/sqrt(3.) - C_BL)*nu_v_km2); // eq. (20) for v 
        
             xi_w = -2*(sqr(omega_w) - sqr(C_BL))/sqr(omega_w + C_BL)*xi_w_km1 
        
                  - sqr(omega_w - C_BL)/sqr(omega_w + C_BL) * xi_w_km2 
        
                  + sig_w*sqrt(3*omega_w/T_V)/sqr(omega_w + C_BL)*( 
        
                        (C_BL + omega_w/sqrt(3.))*nu_w 
        
                      + 2/sqrt(3.)*omega_w*nu_w_km1 
        
                      + (omega_w/sqrt(3.) - C_BL)*nu_w_km2); // eq. (20) for w 
        
             xi_p = -(1 - C_BLp*tau_p)/(1 + C_BLp*tau_p)*xi_p_km1 
        
                  + sig_p*sqrt(2*tau_p/T_V)/(1 + C_BLp*tau_p) * (nu_p + nu_p_km1); // eq. (21) 
        
             xi_q = -(1 - 4*b_w*C_BLq/M_PI/in.V)/(1 + 4*b_w*C_BLq/M_PI/in.V) * xi_q_km1 
        
                  + C_BLq/in.V/(1 + 4*b_w*C_BLq/M_PI/in.V) * (xi_w - xi_w_km1); // eq. (23) 
        
             xi_r = - (1 - 3*b_w*C_BLr/M_PI/in.V)/(1 + 3*b_w*C_BLr/M_PI/in.V) * xi_r_km1 
        
                  + C_BLr/in.V/(1 + 3*b_w*C_BLr/M_PI/in.V) * (xi_v - xi_v_km1); // eq. (25)

Q2: The next condition: if(turbType == ttMilspec is never being run no ? since it's in the case ttTustin ?

jsbsim/src/models/atmosphere/FGWinds.cpp

Lines 346 to 355 in e4aa2a9

    
           } else if (turbType == ttMilspec) { 
        
             // the following is the MIL-STD-1797A formulation 
        
             // as cited in Yeager's report 
        
             xi_u = (1 - T_V/tau_u)  *xi_u_km1 + sig_u*sqrt(2*T_V/tau_u)*nu_u;  // eq. (30) 
        
             xi_v = (1 - 2*T_V/tau_u)*xi_v_km1 + sig_u*sqrt(4*T_V/tau_u)*nu_v;  // eq. (31) 
        
             xi_w = (1 - 2*T_V/tau_w)*xi_w_km1 + sig_w*sqrt(4*T_V/tau_w)*nu_w;  // eq. (32) 
        
             xi_p = (1 - T_V/tau_p)  *xi_p_km1 + sig_p*sqrt(2*T_V/tau_p)*nu_p;  // eq. (33) 
        
             xi_q = (1 - T_V/tau_q)  *xi_q_km1 + M_PI/4/b_w*(xi_w - xi_w_km1);  // eq. (34) 
        
             xi_r = (1 - T_V/tau_r)  *xi_r_km1 + M_PI/3/b_w*(xi_v - xi_v_km1);  // eq. (35) 
        
           }

And by the way I find it odd that the case ttMilspec is empty:

jsbsim/src/models/atmosphere/FGWinds.cpp

Lines 247 to 250 in e4aa2a9

    
           case ttMilspec: 
        
           case ttTustin: { 
        
             // an index of zero means turbulence is disabled

Q3 : Where is the link to the Dryden spectrum model ? I couldn't find any part in the FGWinds.cpp where the Dryden power spectral densities are implemented as presented in the Dryden Wind Turbulence Model wikipedia page or the Matworks Simulink Documentation

So the question sythesizing all the previous ones could be: what is the difference between Tustin and Milspec turbulences and their relationship to the Dryden spectrum model ?

Sorry for the very long post, and thanks in advance !

seanmcleod · 2023-09-01T16:16:56Z

seanmcleod
Sep 1, 2023
Maintainer

Just a quick initial answer for some of your questions.

Q1: Since we're already in a switch turbType; case ttTustin why do we need a if(turbType == ttTustin) condition ?
Q2: The next condition: if(turbType == ttMilspec is never being run no ? since it's in the case ttTustin ?
And by the way I find it odd that the case ttMilspec is empty:

What you're missing is that the ttMispec case isn't empty, there is no break statement after it. So in effect the code below is common to both ttMilspec and ttTustin and hence the reason for those conditional checks since there are parts of the code that are different.

1 reply

Dobid Sep 1, 2023
Author

Alright, I see haha. Dang it 🥲

Dobid · 2023-09-04T15:59:30Z

Dobid
Sep 4, 2023
Author

As an extra question, to give more detail about my "what's the difference between Milspec and Tustin turbulences" question; here I am playing with the Tustin turbulence model, and according to the plots the turbulence linear velocities reach abnormally high velocities ~300 m/s (at least to me) in comparison with Milspec turbulences with same settings (severe, windspeed_at_20ft = 75 fps, + 22 m/s north east constant wind).
Tustin turbulences inner loop control (roll, pitch, Va | ignore two first plots):

Milspec turbulences inner loop control (roll, pitch, Va | ignore two first plots):

0 replies

Dobid · 2023-09-11T10:55:05Z

Dobid
Sep 11, 2023
Author

Friendly Up !

2 replies

agodemar Sep 11, 2023
Maintainer

@Dobid Hello, can you give us more context on your work with JSBSim?

Dobid Sep 11, 2023
Author

Hello,
I'm working on RL based / PID controller comparison for flight control under disturbances, hence me trying to understand precisely how turbulences are generated and the difference between the turbulence models that JSBSim offers.

1st question :
I'm trying to understand the taxonomy of turbulence generation and more particularly how the Milspec and Tustin turbulence models are related to the Dryden spectrum model since I don't see the code for Dryden spectrum power spectral densities implemented in the JSBSim code as in Dryden Wind Turbulence Model wikipedia page or the Matworks Simulink Documentation. The only common code I found are the scale lengths computation but it seems to come from the Continuous gusts formalism.
So I don't really get why there is (Dryden spectrum) next to the Milspec and Tustin models in the FGWinds in the documentation. Is the continuous gust model the same as Dryden spectrum model ? The namings kinda confuse me.

2nd question :
These plots describe the telemetry of the flight of a fixed wing UAV flying under severe Tustin turbulences while being PID controlled for its pitch, roll and airspeed. Looking at the wind speeds (bottom right hand corner) plot, it seems that the wind reaches abnormally high speeds (up to 300m/s for the down wind) which doesn't seem normal to me, especially when compared to the Milspec turbulence model (with the same settings) giving much more coherent orders of magnitude for turbulent wind speeds.

That's it, feel free to ask more details if needed, thanks for your reply !

seanmcleod · 2023-09-11T12:00:34Z

seanmcleod
Sep 11, 2023
Maintainer

I'm trying to understand the taxonomy of turbulence generation...

I don't have a definite answer for this off the top of my head. But did notice this bit of documentation in FGWinds.h which may partially help, or be a starting point.

jsbsim/src/models/atmosphere/FGWinds.h

Lines 56 to 160 in c7cd8ef

    
           /** Models atmospheric disturbances: winds, gusts, turbulence, downbursts, etc. 
        
               <h2>Turbulence</h2> 
        
               Various turbulence models are available. They are specified 
        
               via the property <tt>atmosphere/turb-type</tt>. The following models are 
        
               available: 
        
               - 0: ttNone (turbulence disabled) 
        
               - 1: ttStandard 
        
               - 2: ttCulp 
        
               - 3: ttMilspec (Dryden spectrum) 
        
               - 4: ttTustin (Dryden spectrum) 
        
               The Milspec and Tustin models are described in the Yeager report cited 
        
               below.  They both use a Dryden spectrum model whose parameters (scale 
        
               lengths and intensities) are modelled according to MIL-F-8785C. Parameters 
        
               are modelled differently for altitudes below 1000ft and above 2000ft, for 
        
               altitudes in between they are interpolated linearly. 
        
               The two models differ in the implementation of the transfer functions 
        
               described in the milspec. 
        
               To use one of these two models, set <tt>atmosphere/turb-type</tt> to 4 
        
               resp. 5, and specify values for 
        
               <tt>atmosphere/turbulence/milspec/windspeed_at_20ft_AGL-fps</tt> and 
        
               <tt>atmosphere/turbulence/milspec/severity</tt> (the latter corresponds to 
        
               the probability of exceedence curves from Fig.&nbsp;7 of the milspec, 
        
               allowable range is 0 (disabled) to 7). <tt>atmosphere/psiw-rad</tt> is 
        
               respected as well; note that you have to specify a positive wind magnitude 
        
               to prevent psiw from being reset to zero. 
        
               Reference values (cf. figures 7 and 9 from the milspec): 
        
               <table> 
        
                 <tr><td><b>Intensity</b></td> 
        
                     <td><b><tt>windspeed_at_20ft_AGL-fps</tt></b></td> 
        
                     <td><b><tt>severity</tt></b></td></tr> 
        
                 <tr><td>light</td> 
        
                     <td>25 (15 knots)</td> 
        
                     <td>3</td></tr> 
        
                 <tr><td>moderate</td> 
        
                     <td>50 (30 knots)</td> 
        
                     <td>4</td></tr> 
        
                 <tr><td>severe</td> 
        
                     <td>75 (45 knots)</td> 
        
                     <td>6</td></tr> 
        
               </table> 
        
               <h2>Cosine Gust</h2> 
        
               A one minus cosine gust model is available. This permits a configurable, 
        
               predictable gust to be input to JSBSim for testing handling and 
        
               dynamics. Here is how a gust can be entered in a script: 
        
               ~~~{.xml} 
        
               <event name="Introduce gust"> 
        
                 <condition> simulation/sim-time-sec ge 10 </condition> 
        
                 <set name="atmosphere/cosine-gust/startup-duration-sec" value="5"/> 
        
                 <set name="atmosphere/cosine-gust/steady-duration-sec" value="1"/> 
        
                 <set name="atmosphere/cosine-gust/end-duration-sec" value="5"/> 
        
                 <set name="atmosphere/cosine-gust/magnitude-ft_sec" value="30"/> 
        
                 <set name="atmosphere/cosine-gust/frame" value="2"/> 
        
                 <set name="atmosphere/cosine-gust/X-velocity-ft_sec" value="-1"/> 
        
                 <set name="atmosphere/cosine-gust/Y-velocity-ft_sec" value="0"/> 
        
                 <set name="atmosphere/cosine-gust/Z-velocity-ft_sec" value="0"/> 
        
                 <set name="atmosphere/cosine-gust/start" value="1"/> 
        
                 <notify/> 
        
               </event> 
        
               ~~~ 
        
               The x, y, z velocity components are meant to define the direction vector. 
        
               The vector will be normalized by the routine, so it does not need to be a 
        
               unit vector. 
        
               The startup duration is the time it takes to build up to full strength 
        
               (magnitude-ft_sec) from zero. Steady duration is the time the gust stays at 
        
               the specified magnitude. End duration is the time it takes to dwindle to 
        
               zero from the specified magnitude. The start and end transients are in a 
        
               smooth cosine shape. 
        
               The frame is specified from the following enum: 
        
               enum eGustFrame {gfNone=0, gfBody, gfWind, gfLocal}; 
        
               That is, if you specify the X, Y, Z gust direction vector in the body frame, 
        
               frame will be "1". If the X, Y, and Z gust direction vector is in the Wind 
        
               frame, use frame = 2. If you specify the gust direction vector in the local 
        
               frame (N-E-D) use frame = 3. Note that an internal local frame direction 
        
               vector is created based on the X, Y, Z direction vector you specify and the 
        
               frame *at the time the gust is begun*. The direction vector is not updated 
        
               after the initial creation. This is to keep the gust at the same direction 
        
               independent of aircraft dynamics. 
        
               The gust is triggered when the property atmosphere/cosine-gust/start is set 
        
               to 1. It can be used repeatedly - the gust resets itself after it has 
        
               completed. 
        
               The cosine gust is global: it affects the whole world not just the vicinity 
        
               of the aircraft. 
        
               @see Yeager, Jessie C.: "Implementation and Testing of Turbulence Models for 
        
                    the F18-HARV" (<a 
        
                    href="http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19980028448_1998081596.pdf"> 
        
                    pdf</a>), NASA CR-1998-206937, 1998 
        
               @see MIL-F-8785C: Military Specification: Flying Qualities of Piloted Aircraft 
        
           */

1 reply

Dobid Sep 11, 2023
Author

Ok, that's the documentation I took a screenshot of, in my last reply. I may not keep on searching about it since it took me a lot of time already and I'm not sure if diving for so long in the details of old military specification reports is going to be worth the time for my work.

As for my second question about the Tustin turbulence simulation. Can you reproduce similar results for the wind speeds under the severe turbulence setting and confirm that those turbulences are giving abnormally high wind speeds ?

seanmcleod · 2023-09-11T12:20:17Z

seanmcleod
Sep 11, 2023
Maintainer

I'm not sure if diving for so long in the details of old military specification reports is going to be worth the time for my work.

Short of one of the original coders happening to see your discussion here and answering it I'm not sure you have much of a choice, unless someone else decides to read through the specification documents to come up with answers.

As for my second question about the Tustin turbulence simulation. Can you reproduce similar results for the wind speeds....

Haven't tried it yet.

1 reply

Dobid Sep 11, 2023
Author

Thanks for your replies!
As an extra information for the reproductibilty of my 2nd question, I launched JSBSim by setting the simulation/randomseed property to 1.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

JSBSim Turbulence Simulation #953

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 5 comments 5 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

Select a reply

JSBSim Turbulence Simulation #953

Dobid Sep 1, 2023

Replies: 5 comments · 5 replies

seanmcleod Sep 1, 2023 Maintainer

Dobid Sep 1, 2023 Author

Dobid Sep 4, 2023 Author

Dobid Sep 11, 2023 Author

agodemar Sep 11, 2023 Maintainer

Dobid Sep 11, 2023 Author

seanmcleod Sep 11, 2023 Maintainer

Dobid Sep 11, 2023 Author

seanmcleod Sep 11, 2023 Maintainer

Dobid Sep 11, 2023 Author

Dobid
Sep 1, 2023

Replies: 5 comments 5 replies

seanmcleod
Sep 1, 2023
Maintainer

Dobid Sep 1, 2023
Author

Dobid
Sep 4, 2023
Author

Dobid
Sep 11, 2023
Author

agodemar Sep 11, 2023
Maintainer

Dobid Sep 11, 2023
Author

seanmcleod
Sep 11, 2023
Maintainer

Dobid Sep 11, 2023
Author

seanmcleod
Sep 11, 2023
Maintainer

Dobid Sep 11, 2023
Author