uavWindTurbulence

Generate wind turbulence for UAV platform with discretized Von Kármán velocity spectra

Since R2024a

Description

The uavWindTurbulence value object uses the discretized Von Kármán spectral representation to generate turbulence for a UAV platform by passing band-limited white noise through appropriate forming filters. You can choose the mathematical representation of the Von Kármán turbulence from US Military Specification MIL-F-8785C, MIL-HDBK-1797, and MIL-HDBK-1797B. For more information, see Algorithms and .

The object discretizes the Von Kármán spectral representation with a rate similar to the UAV scenario UpdateRate property. To simulate wind turbulence in a UAV scenario:

1. Create the uavWindTurbulence object and set its properties.

2. Attach the wind model to a UAV platform by using the addWind function.

3. Obtain the wind velocity from the UAV scenario by using the windVelocity function, and the wind angular rate using the windAngularRate while the UAV scenario is running.

Note

Simulating a UAV scenario that contains a uavWindTurbulence object requires Aerospace Toolbox.

Creation

Syntax

turbulenceWind = uavWindTurbulence

turbulenceWind = uavWindTurbulence(Name=Value)

Description

turbulenceWind = uavWindTurbulence creates a wind turbulence object with default property values.

turbulenceWind = uavWindTurbulence(Name=Value) specifies properties using one or more name-value arguments. For example, uavWindTurbulence(StartTime=10) sets the wind start time to 10 seconds.

Properties

You can modify these value object properties before you add the uavWindTurbulence value object to the UAV platform using addWind function.

StartTime — Wind start time

0 (default) | nonnegative scalar

Wind start time, specified as a nonnegative scalar. This value is relative to the simulation start time of the UAV scenario. Units are in seconds

Example: 2

Data Types: double

StopTime — Wind stop time

Inf (default) | nonnegative scalar

Wind stop time, specified as a nonnegative scalar. This value is relative to the simulation start time of the UAV scenario. If specified as inf, the wind gust continues until the UAV scenario simulation ends. Units are in seconds

Example: 5

Data Types: double

Specification — Military reference

"MIL-F-8785C" (default) | "MIL-HDBK-1797" | "MIL-HDBK-1797B"

Military reference, specified as"MIL-F-8785C", "MIL-HDBK-1797", or "MIL-HDBK-1797B". The military reference affects the application of turbulence scale lengths in the lateral and vertical directions.

ModelType — Angular rate signs

"Von Karman (+q +r)" (default) | "Von Karman (+q -r)" | "Von Karman (-q +r)"

Angular rate signs, specified as one of these options:

• "Von Karman (+q +r)" — Positive vertical and lateral angular rate spectra.

• "Von Karman (+q -r)" — Positive vertical and negative lateral angular rate spectra.

• "Von Karman (-q +r)" — Negative vertical and positive lateral angular rate spectra.

Data Types: string

WindSpeedAt6m — Wind speed at 6 meters above ground

15 (default) | nonnegative scalar

Wind speed at a height of 6 meters above ground, specified as a nonnegative scalar in m/s.

Data Types: double

WindDirectionAt6m — Wind direction at 6 meters above ground

0 (default) | nonnegative scalar

Wind direction at a height of 6 meters above ground, specified as a nonnegative in degrees clockwise from north.

Data Types: double

ProbExceedanceHighAltitudeIntensity — Probability of exceeding turbulence intensity

"10^-2 - Light" (default) | "2x10^-1" | "10^-1" | "10^-3 - Moderate" | "10^-4" | "10^-5 - Severe" | "10^-6"

Probability of the turbulence intensity being exceeded, specified as "2x10^-1", "10^-1", "10^-2 - Light", "10^-3 - Moderate", "10^-4", "10^-5 - Severe", or "10^-6"

Data Types: string

ScaleLengthAtMediumHighAltitude — Turbulence scale length

762 (default) | positive scalar

Turbulence scale length above 2000 feet, specified as a real scalar. This length is assumed constant.

Data Types: double

WingSpan — Wingspan

10 (default) | positive scalar

Wingspan, specified in meters.

Data Types: double

BandLimitedNoiseSampleTime — Noise sample time

0.1 (default) | positive scalar

Noise sample time at which the unit variance white noise signal is generated, specified in seconds.

Noise sample time must be an integer multiple of the scene sample time.

• In MATLAB^®, the scene sample time is determined by UpdateRate property of uavScenario.

• In Simulink^®, the scene sample time is determined by the Sample time property of UAV Scenario Configuration block

Data Types: double

RandomNoiseSeeds — Noise seeds

[23341 23342 23343 23344] (default) | vector of the form [u_g v_g w_g p_g]

Random noise seeds, specified as a vector of the form [u_g v_g w_g p_g] which generates the turbulence for each of the three velocity components and the roll rate.

Data Types: double

Examples

Generate Turbulence Wind in UAV Scenario

Create a UAV scenario with default parameters.

scene = uavScenario;

Create a UAV platform. Specify the trajectory to 80 meters north of scenario origin for 80 seconds.

distance = 80;
duration = 80;
platform = uavPlatform("UAV",scene,Trajectory=waypointTrajectory([0 0 -50;distance 0 -50],[0 duration]));

Create a turbulence wind object which starts at 20 seconds. Add the wind object to the UAV platform.

turbulenceWind = uavWindTurbulence(StartTime=20, BandLimitedNoiseSampleTime=1/scene.UpdateRate);
addWind(platform,turbulenceWind);

Set up the scenario.

setup(scene);

% Create arrays of zeros to store longitudinal, horizontal, and vertical wind velocity
uw = zeros(scene.UpdateRate*duration, 1);
vw = zeros(scene.UpdateRate*duration, 1);
ww = zeros(scene.UpdateRate*duration, 1);

% Create arrays of zeros to store x position of the UAV
xpos = zeros(scene.UpdateRate*duration, 1);

% Start a counter index
idx = 1;

Run the scenario. Obtain the vertical wind velocity at each time step.

while scene.CurrentTime <= duration
    % Obtain resultant wind velocity
    vel = windVelocity(platform);
    
    % Obtain UAV position
    [motion,lla]=read(platform);
    
    % Obtain UAV X position
    xpos(idx) = motion(1);

    % Obtain wind velocities
    uw(idx) = vel(1);
    vw(idx) = vel(2);
    ww(idx) = vel(3);

    % Increment counter index
    idx = idx + 1;

    % Advance the scenario time step
    advance(scene);
end

Plot the wind velocities against distance.

tiledlayout(3,1)

nexttile
plot(xpos,uw)
title("Longitudinal Wind Velocity")
xlabel("Distance (m)")
ylabel("Wind Velocity (m/s)")

nexttile
plot(xpos,vw)
title("Lateral Wind Velocity")
xlabel("Distance (m)")
ylabel("Wind Velocity (m/s)")

nexttile
plot(xpos,ww)
title("Vertical Wind Velocity")
xlabel("Distance (m)")
ylabel("Wind Velocity (m/s)")

Limitations

The frozen turbulence field assumption that this object uses is valid when the aircraft ground speed is large compared to the root-mean-square turbulence velocity or intensity.

For an example on how to model turbulence effect on UAV which flies at slow speed or hovers, visit the Simulate Wind Turbulence Using CETI Model example.

More About

Algorithms

According to the military references, turbulence is a stochastic process defined by velocity spectra. For an aircraft flying at a speed V through a frozen turbulence field with a spatial frequency of Ω radians per meter, the circular frequency ω is calculated by multiplying V by Ω. The following table displays the component spectra functions:

	MIL-F-8785C	MIL-HDBK-1797 and MIL-HDBK-1797B
Longitudinal
$$H_u\left(s\right)$$	$$\frac{2{\sigma }_u^2{L}_u}{\pi V}\cdot \frac{1}{{\left[1+{\left(1.339L_u\frac{\omega }{V}\right)}^2\right]}^{\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$6$}\right.}}$$	$$\frac{2{\sigma }_u^2{L}_u}{\pi V}\cdot \frac{1}{{\left[1+{\left(1.339L_u\frac{\omega }{V}\right)}^2\right]}^{\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$6$}\right.}}$$
$$H_p\left(s\right)$$	$$\frac{{\sigma }_w^2}{V{L}_w}\cdot \frac{0.8{\left(\frac{\pi L_w}{4b}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}{1+{\left(\frac{4b\omega }{\pi V}\right)}^2}$$	$$\frac{{\sigma }_w^2}{2V{L}_w}\cdot \frac{0.8{\left(\frac{2\pi L_w}{4b}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}{1+{\left(\frac{4b\omega }{\pi V}\right)}^2}$$
Lateral
$${\Phi }_v\left(\omega \right)$$	$$\frac{{\sigma }_v^2{L}_v}{\pi V}\cdot \frac{1+\frac{8}{3}{\left(1.339L_v\frac{\omega }{V}\right)}^2}{{\left[1+{\left(1.339L_v\frac{\omega }{V}\right)}^2\right]}^{\raisebox{1ex}{$11$}\!\left/ \!\raisebox{-1ex}{$6$}\right.}}$$	$$\frac{2{\sigma }_v^2{L}_v}{\pi V}\cdot \frac{1+\frac{8}{3}{\left(2.678L_v\frac{\omega }{V}\right)}^2}{{\left[1+{\left(2.678L_v\frac{\omega }{V}\right)}^2\right]}^{\raisebox{1ex}{$11$}\!\left/ \!\raisebox{-1ex}{$6$}\right.}}$$
$${\Phi }_r\left(\omega \right)$$	$$\frac{\mp {\left(\frac{\omega }{V}\right)}^2}{1+{\left(\frac{3b\omega }{\pi V}\right)}^2}\cdot {\Phi }_v\left(\omega \right)$$	$$\frac{\mp {\left(\frac{\omega }{V}\right)}^2}{1+{\left(\frac{3b\omega }{\pi V}\right)}^2}\cdot {\Phi }_v\left(\omega \right)$$
Vertical
$${\Phi }_w\left(\omega \right)$$	$$\frac{{\sigma }_w^2{L}_w}{\pi V}\cdot \frac{1+\frac{8}{3}{\left(1.339L_w\frac{\omega }{V}\right)}^2}{{\left[1+{\left(1.339L_w\frac{\omega }{V}\right)}^2\right]}^{\raisebox{1ex}{$11$}\!\left/ \!\raisebox{-1ex}{$6$}\right.}}$$	$$\frac{2{\sigma }_w^2{L}_w}{\pi V}\cdot \frac{1+\frac{8}{3}{\left(2.678L_w\frac{\omega }{V}\right)}^2}{{\left[1+{\left(2.678L_w\frac{\omega }{V}\right)}^2\right]}^{\raisebox{1ex}{$11$}\!\left/ \!\raisebox{-1ex}{$6$}\right.}}$$
$${\Phi }_q\left(\omega \right)$$	$$\frac{\pm {\left(\frac{\omega }{V}\right)}^2}{1+{\left(\frac{4b\omega }{\pi V}\right)}^2}\cdot {\Phi }_w\left(\omega \right)$$	$$\frac{\pm {\left(\frac{\omega }{V}\right)}^2}{1+{\left(\frac{4b\omega }{\pi V}\right)}^2}\cdot {\Phi }_w\left(\omega \right)$$

The variable b represents the aircraft wingspan. The variables L_u, L_v, L_w represent the turbulence scale lengths. The variables σ_u, σ_v, σ_w represent the turbulence intensities. For more information, see Turbulence Scale Length and Intensities

Keep in mind that the longitudinal turbulence angular rate spectrum, Ф_p(ω), is a rational function. The rational function is derived from curve-fitting a complex algebraic function, not the vertical turbulence velocity spectrum, Ф_w(ω), multiplied by a scale factor.

The spectral density definitions of turbulence angular rates are defined in the references as three variations, which are displayed in the following table:

$$p_g=\frac{\partial w_g}{\partial y}$$	$$q_g=\frac{\partial w_g}{\partial x}$$	$$r_g=-\frac{\partial v_g}{\partial x}$$
$$p_g=\frac{\partial w_g}{\partial y}$$	$$q_g=\frac{\partial w_g}{\partial x}$$	$$r_g=\frac{\partial v_g}{\partial x}$$
$$p_g=-\frac{\partial w_g}{\partial y}$$	$$q_g=-\frac{\partial w_g}{\partial x}$$	$$r_g=\frac{\partial v_g}{\partial x}$$

The variations affect only the vertical (q_g) and lateral (r_g) turbulence angular rates.

The turbulence angular rate spectra has varying definitions because it contribute less to the aircraft gust response than the turbulence velocity spectra. The following table summarizes the combinations of vertical and lateral turbulence angular rate spectra.

Vertical	Lateral
Фq(ω) Фq(ω) −Фq(ω)	−Фr(ω) Фr(ω) Фr(ω)

To generate a signal with the correct characteristics, a unit variance, band-limited white noise signal is passed through forming filters. The forming filters are approximations of the Von Kármán velocity. The following tables displays the filters.

	MIL-F-8785C	MIL-HDBK-1797 and MIL-HDBK-1797B
Longitudinal
$$H_u\left(s\right)$$	$$\frac{{\sigma }_u\sqrt{\frac{2}{\pi }\cdot \frac{L_u}{V}}\left(1+0.25\frac{L_u}{V}s\right)}{1+1.357\frac{L_u}{V}s+0.1987{\left(\frac{L_u}{V}\right)}^2{s}^2}$$	$$\frac{{\sigma }_u\sqrt{\frac{2}{\pi }\cdot \frac{L_u}{V}}\left(1+0.25\frac{L_u}{V}s\right)}{1+1.357\frac{L_u}{V}s+0.1987{\left(\frac{L_u}{V}\right)}^2{s}^2}$$
$$H_p\left(s\right)$$	$${\sigma }_w\sqrt{\frac{0.8}{V}}\cdot \frac{{\left(\frac{\pi }{4b}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$6$}\right.}}{L_w{}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\left(1+\left(\frac{4b}{\pi V}\right)s\right)}$$	$${\sigma }_w\sqrt{\frac{0.8}{V}}\cdot \frac{{\left(\frac{\pi }{4b}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$6$}\right.}}{{\left(2L_w\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\left(1+\left(\frac{4b}{\pi V}\right)s\right)}$$
Lateral
$$H_v\left(s\right)$$	$$\frac{{\sigma }_v\sqrt{\frac{1}{\pi }\cdot \frac{L_v}{V\text{}}}\left(1+2.7478\frac{L_v}{V}s+0.3398{\left(\frac{L_v}{V}\right)}^2{s}^2\right)}{1+2.9958\frac{L_v}{V}s+1.9754{\left(\frac{L_v}{V}\right)}^2{s}^2+0.1539{\left(\frac{L_v}{V}\right)}^3{s}^3}$$	$$\frac{{\sigma }_v\sqrt{\frac{1}{\pi }\cdot \frac{2L_v}{V}}\left(1+2.7478\frac{2L_v}{V}s+0.3398{\left(\frac{2L_v}{V}\right)}^2{s}^2\right)}{1+2.9958\frac{2L_v}{V}s+1.9754{\left(\frac{2L_v}{V}\right)}^2{s}^2+0.1539{\left(\frac{2L_v}{V}\right)}^3{s}^3}$$
$$H_r\left(s\right)$$	$$\frac{\mp \frac{s}{V}}{\left(1+\left(\frac{3b}{\pi V}\right)s\right)}\cdot H_v\left(s\right)$$	$$\frac{\mp \frac{s}{V}}{\left(1+\left(\frac{3b}{\pi V}\right)s\right)}\cdot H_v\left(s\right)$$
Vertical
$$H_w\left(s\right)$$	$$\frac{{\sigma }_w\sqrt{\frac{1}{\pi }\cdot \frac{L_w}{V}}\left(1+2.7478\frac{L_w}{V}s+0.3398{\left(\frac{L_w}{V}\right)}^2{s}^2\right)}{1+2.9958\frac{L_w}{V}s+1.9754{\left(\frac{L_w}{V}\right)}^2{s}^2+0.1539{\left(\frac{L_w}{V}\right)}^3{s}^3}$$	$$\frac{{\sigma }_w\sqrt{\frac{1}{\pi }\cdot \frac{2L_w}{V}}\left(1+2.7478\frac{2L_w}{V}s+0.3398{\left(\frac{2L_w}{V}\right)}^2{s}^2\right)}{1+2.9958\frac{2L_w}{V}s+1.9754{\left(\frac{2L_w}{V}\right)}^2{s}^2+0.1539{\left(\frac{2L_w}{V}\right)}^3{s}^3}$$
$$H_q\left(s\right)$$	$$\frac{\pm \frac{s}{V}}{\left(1+\left(\frac{4b}{\pi V}\right)s\right)}\cdot H_w\left(s\right)$$	$$\frac{\pm \frac{s}{V}}{\left(1+\left(\frac{4b}{\pi V}\right)s\right)}\cdot H_w\left(s\right)$$

In UAV scenario, these filters are discretized with the sampling rate similar to the UAV scenario update rate. To minimize the discretization error, you must choose a scenario update rate in accordance with the UAV velocity and altitude.

Turbulence Scale Length and Intensities

Low-Altitude Model (Altitude <1000 feet)

According to the military references, the turbulence scale lengths at low altitudes, where h is the altitude in feet, are represented in the following table:

MIL-F-8785C	MIL-HDBK-1797 and MIL-HDBK-1797B
$$\begin{array}{l}{L}_w=h\\ L_u=L_v=\frac{h}{{\left(0.177+0.000823h\right)}^{1.2}}\end{array}$$	$$\begin{array}{l}2L_w=h\\ L_u=2L_v=\frac{h}{{\left(0.177+0.000823h\right)}^{1.2}}\end{array}$$

The turbulence intensities are given below, where W₂₀ is the wind speed at 20 feet (6 m). Typically for light turbulence, the wind speed at 20 feet is 15 knots; for moderate turbulence, the wind speed is 30 knots; and for severe turbulence, the wind speed is 45 knots.

$$\begin{array}{l}{\sigma }_w=0.1W_{20}\\ \frac{{\sigma }_u}{{\sigma }_w}=\frac{{\sigma }_v}{{\sigma }_w}=\frac{1}{{\left(0.177+0.000823h\right)}^{0.4}}\end{array}$$

The turbulence axes orientation in this region is defined as follows:

• Longitudinal turbulence velocity, u_g, aligned along the horizontal relative mean wind vector.

• Vertical turbulence velocity, w_g, aligned with vertical relative mean wind vector.

Medium/High Altitudes (Altitude > 2000 feet)

For medium to high altitudes the turbulence scale lengths and intensities are based on the assumption that the turbulence is isotropic. In the military references, the scale lengths are represented by the following equations:

MIL-F-8785C	MIL-HDBK-1797 and MIL-HDBK-1797B
L u = L v = L w = 2500 ft	L u = 2 L v = 2 L w = 2500 ft

The turbulence intensities are determined from a lookup table that provides the turbulence intensity as a function of altitude and the probability of the turbulence intensity being exceeded. The relationship of the turbulence intensities is represented in the following equation: σ_u=σ_v=σ_w.

Between Low and Medium/High Altitudes (1000 feet < Altitude < 2000 feet)

At altitudes between 1000 feet and 2000 feet, the turbulence velocities and turbulence angular rates are determined by linearly interpolating between the value from the low altitude model at 1000 feet transformed from mean horizontal wind coordinates to body coordinates and the value from the high altitude model at 2000 feet in body coordinates.

Version History

Introduced in R2024a