A simulation method for aircraft turbulence environment simulation

Abstract

The present invention discloses a simulation method for aircraft turbulence environment simulation, and the method includes the following steps: entering the turbulence module, and judging whether the aircraft speed reset flag is true; judging whether the aircraft flight freeze flag is true; judging whether the true airspeed V of the aircraft TR is greater than 15.0; judging whether the aircraft reset directory is 0; judging whether the aircraft wind profile index is not 0; judging whether the aircraft microburst activation flag is true; setting the required rough atmosphere flag of the aircraft to be true; judging whether the aircraft rough atmosphere flag is true; judging whether the aircraft turbulence type is 0; judging whether the aircraft turbulence intensity is 0; judging whether the aircraft repositioning directory is not 0. Beneficial effects: Turbulence models such as rough atmosphere, pebble turbulence, and discrete vertical gusts can be provided, and the turbulence intensity is level 10, and the gust intensity can be selected within the range of 0 to 25 m / s, so as to simulate the atmospheric disturbance environment by calculating the turbulence disturbance speed.

Description

Technical Field

[0001] The present invention relates to the technical field of simulation methods, and more specifically, to a simulation method for simulating the turbulent environment of an aircraft. Background Art

[0002] Turbulence, also known as turbulent flow, is a flow state of a fluid. When the flow velocity is very small, the fluid flows in layers without mixing, which is called laminar flow, or sheet flow; gradually increasing the flow velocity, the streamlines of the fluid begin to show wavy oscillations, and the frequency and amplitude of the oscillations increase with the increase of the flow velocity. This flow condition is called transitional flow; when the flow velocity increases to a very large value, the streamlines are no longer clearly distinguishable, and there are many small vortices in the flow field, which is called turbulent flow, also known as turbulent, disturbed or chaotic flow.

[0003] With the continuous improvement of living standards, traveling by plane is becoming more and more popular among people. When an aircraft flies in the atmosphere, it is often easily affected by atmospheric turbulence. The interference of atmospheric turbulence will cause the flight performance of the aircraft to decline, make it difficult to control, and even cause the pilot to induce oscillations, endangering flight safety. Therefore, the present invention proposes a simulation method for simulating the turbulent environment of an aircraft. Summary of the Invention

[0004] In view of the problems in the related art, the present invention proposes a simulation method for simulating the turbulent environment of an aircraft to overcome the above-mentioned technical problems existing in the existing related technologies.

[0005] For this purpose, the specific technical solution adopted by the present invention is as follows:

[0006] A simulation method for simulating the turbulent environment of an aircraft, the method comprising the following steps:

[0007] S1. Enter the turbulence module, and determine whether the aircraft speed reset flag is true. If it is, execute S19; if not, execute S2;

[0008] S2. Determine whether the aircraft flight freeze flag is true. If it is, execute S19; if not, execute S3;

[0009] S3. Determine whether the true airspeed V TR of the aircraft is greater than 15.0. If it is, execute S4; if not, execute S19;

[0010] S4. Determine whether the aircraft reset directory is 0. If it is, execute S5; if not, execute S19;

[0011] S5. Determine whether the aircraft wind profile index is not 0. If it is, execute S7; if not, execute S6;

[0012] S6. Determine whether the aircraft microburst activation flag is true. If it is, execute S7; if not, set the aircraft's required rough atmosphere flag to false and execute S8;

[0013] S7. Set the aircraft's required rough atmosphere flag to true and execute S8;

[0014] S8. Determine whether the aircraft rough atmosphere flag is true. If it is, execute S9; if not, execute S10;

[0015] S9. Determine whether the aircraft turbulence type is 0. If it is, execute S12; if not, execute S10;

[0016] S10. Determine whether the aircraft turbulence intensity is 0. If it is, execute S12; if not, execute S11;

[0017] S11. Determine whether the aircraft repositioning directory is not 0. If it is, execute S12; if not, execute S13;

[0018] S12. Determine whether the aircraft turbulence intensity factor η TU is 0. If it is, execute S14; if not, execute S13;

[0019] S13. Determine whether the aircraft turbulence type is not 3. If it is, call the clear vertical gust sub - process and execute S15; if not, call the calculate vertical gust sub - process and execute S15;

[0020] S14. First call the clear turbulence sub - process, then call the clear vertical gust sub - process, and exit the turbulence module;

[0021] S15. Determine whether the aircraft gust timing T JETUP is less than 30.0 seconds. If it is, execute S17; if not, execute S16;

[0022] S16. Determine whether the aircraft turbulence type is 3. If it is, execute S17; if not, execute S18;

[0023] S17. Call the clear turbulence sub - process and exit the turbulence module;

[0024] S18. Call the calculate other turbulence sub - process, call the generate random number sub - process, call the filter sub - process, call the calculate turbulence rate sub - process in sequence, and exit the turbulence module;

[0025] S19. Determine whether the aircraft re - positioning timing T REP is not 0. If it is, execute S21; if not, execute S20;

[0026] S20. Determine whether the aircraft true airspeed V TRIs it less than 15.0? If so, execute S21; if not, exit the turbulent flow module.

[0027] S21. Sequentially obtain the clear turbulent flow velocity U T = V T = W T = 0, η TU = η TUL = η TUD = 0, exit the turbulent flow module, where U T is the X-direction velocity output by the DRYDEN model, V T is the Y-direction velocity output by the DRYDEN model, W T is the Z-direction velocity output by the DRYDEN model, η TUL is the turbulent flow intensity factor of the previous cycle.

[0028] Furthermore, the vertical gust clearing sub-process in S13 includes the following steps:

[0029] S1301. Enter the vertical gust clearing sub-process;

[0030] S1302. Set both the gust distance factor and the gust timing to 0;

[0031] S1303. Set both the gust intensity factor and the vertical gust intensity to 0;

[0032] S1304. Set both the gust-required pebble turbulent flow intensity and the vertical velocity of the gust to 0;

[0033] S1305. Exit the vertical gust clearing sub-process.

[0034] Furthermore, the vertical gust calculation sub-process in S13 includes the following steps:

[0035] S1311. Enter the vertical gust calculation sub-process;

[0036] S1312. Determine whether the turbulent flow intensity is not 0. If so, obtain the vertical gust turbulent flow intensity factor where IAXTURIN represents the turbulent flow intensity. If not, obtain the vertical gust turbulent flow intensity factor η JETLOC = 0;

[0037] S1313. Determine whether the aircraft gust timing T JETUP is less than 30.0. If so, obtain T JETUP = T JETUP + Δt and execute S1314, where Δt represents the timing time interval. If not, obtain the distance factor = 0.008Mn + the distance factor, η JETLOC = 0 and execute S1315, where Mn represents the flight Mach number;

[0038] S1314. Determine T JETUP Whether it is greater than 40.0. If so, calculate the cobblestone turbulence intensity required for gust turbulence and execute S1317. If not, execute S1317;

[0039] S1315. Determine whether the distance factor is greater than or equal to 1.0. If so, execute S1316. If not, execute S1317;

[0040] S1316. Sequentially set the turbulence type number to 0, the turbulence intensity number to 0, the required gust flag to false, the light turbulence flag to false, the moderate turbulence flag to false, and the severe turbulence flag to false, and execute S1317;

[0041] S1317. Interpolate and calculate the vertical gust height factor η JUM = f(h1), where h1 represents the height;

[0042] S1318. Obtain the vertical gust velocity W JETUP = f(h2)·η JETLOC ·η JUM , where h2 represents the distance factor;

[0043] S1319. Exit the sub - process of calculating the vertical gust.

[0044] Furthermore, the sub - process of clearing turbulence in S14 includes the following steps:

[0045] S141. Enter the sub - process of clearing turbulence;

[0046] S142. Call the first fading sub - process (U RAL , U RA ) and the second fading sub - process (V RAL , V RA ), where U RAL represents the longitudinal component factor of the turbulent wind speed in the aircraft body system, U RA represents the longitudinal component of the turbulent wind speed in the aircraft body system, V RAL represents the lateral component factor of the turbulent wind speed in the aircraft body system, V RA represents the lateral component of the turbulent wind speed in the aircraft body system;

[0047] S143. Call the third fading sub - process (W RAL , W RA ) and the fourth fading sub - process (P RAL , P RA ), where W RAL represents the vertical component factor of the turbulent wind speed in the aircraft body system, W RA represents the vertical component of the turbulent wind speed in the aircraft body system, PRAL Indicates the pitch angular velocity factor generated by the turbulent wind speed, P RA Indicates the pitch angular velocity generated by the turbulent wind speed;

[0048] S144. Call the fifth fading sub - process (Q RAL , Q RA ) and the sixth fading sub - process (R RAL , R RA ), where Q RAL Indicates the yaw angular velocity factor generated by the turbulent wind speed, Q RA Indicates the yaw angular velocity generated by the turbulent wind speed, R RAL Indicates the roll angular velocity factor generated by the turbulent wind speed, R RA Indicates the roll angular velocity generated by the turbulent wind speed;

[0049] S145. Obtain η TU = η TUL = η TUD = 0;

[0050] S146. Clear all filter outputs F UI = 0, F VI = 0, F WI = 0, I = 1, 2, 3, where F UI Indicates the longitudinal filter output, F VI Indicates the lateral filter output, F WI Indicates the vertical filter output;

[0051] S147. Exit the turbulent wind clearing sub - process.

[0052] Furthermore, the fading sub - process includes a first fading sub - process (U RAL , U RA ), a second fading sub - process (V RAL , V RA ), a third fading sub - process (W RAL , W RA ), a fourth fading sub - process (P RAL , P RA ), a fifth fading sub - process (Q RAL , Q RA ) and a sixth fading sub - process (R RAL , R RA ).

[0053] Furthermore, the fading sub - process includes the following steps:

[0054] Enter the fading sub - process and calculate using the formula RL = RAKFADER * RC, where RL represents U RAL , V RAL , WRAL , P RAL , Q RAL , R RAL , RAKFADER = 0.95, representing the fixed value of the factor for calculation speed, RC represents U RA , V RA , W RA , P RA , Q RA , R RA ;

[0055] Judge whether |RL| is less than 0.01. If so, then obtain RC = 0 and exit the fade - out sub - process. If not, obtain RC = RL and exit the fade - out sub - process.

[0056] Furthermore, the calculation of other turbulent sub - processes in S18 includes the following steps:

[0057] S1801. Enter the calculation of other turbulent sub - processes;

[0058] S1802. Judge whether the rough atmosphere requirement flag is true. If so, calculate the turbulence intensity factor and execute S1806. If not, execute S1803;

[0059] S1803. Judge whether the turbulence type is rough atmosphere. If so, obtain the turbulence intensity factor required for rough atmosphere where IAXTURIN represents the turbulence intensity and execute S1806. If not, execute S1804;

[0060] S1804. Judge whether the turbulence type is cobblestone turbulence. If so, obtain the turbulence intensity factor required for cobblestone atmosphere and execute S1805. If not, obtain η TUD = η JETLOC , and execute S1806;

[0061] S1805. Obtain the square - wave intensity and execute S1806;

[0062] S1806. Judge whether the turbulence intensity number is 0. If so, execute S1807. If not, execute 1808;

[0063] S1807. Set the rough atmosphere requirement flag to true. If it has been set, then obtain η TUD = 0 and execute S1808. If it has not been set, execute S1808;

[0064] S1808. Obtain the increment of the turbulence intensity factor Δη TU = η TUD - η TUL , and judge ΔηTU Is it greater than 0.01? If so, obtain Δη TU = 0.01, and execute S1810. If not, execute S1809; where η TUL is the turbulent intensity factor of the previous cycle;

[0065] S1809. Judge Δη TU Is it less than -0.05? If so, obtain Δη TU = -0.05. If not, execute S1810;

[0066] S1810. Obtain η TU = η TUL + Δη TU and η TUL = η TU , where η TUL represents the turbulent intensity factor of the previous cycle, and Δη TU represents the increment of the turbulent intensity factor;

[0067] S1811. Calculate the turbulent intensity height factor and the turbulent scale height factor by interpolation, and judge whether the turbulent type is rough atmosphere. If so, execute S1813. If not, execute S1812;

[0068] S1812. Judge whether the rough atmosphere requirement flag is true. If so, execute S1813. If not, execute S1814;

[0069] S1813. Calculate the rough atmosphere turbulent intensity and the turbulent scale in sequence, and obtain the clear cobblestone square wave intensity W COB = 0, A COB = 0, and execute S1821;

[0070] S1814. Judge whether the turbulent type is cobblestone turbulence. If so, execute S1816. If not, execute S1815;

[0071] S1815. Judge whether the turbulent type is vertical gust. If so, execute S1816. If not, execute S1821;

[0072] S1816. Obtain the cobblestone timing T CS = T CS + Δt, and judge whether the aircraft true airspeed V TR is greater than 700.0. If so, obtain Δη TU = 0.01, and execute S1818. If not, execute S1817;

[0073] S1817. Judge whether the true airspeed V TR is less than 200.0. If so, obtain the rough atmosphere timing TP = 0.333333, and proceed to the next step. If not, obtain T P = 1 / [0.006(V TR - 200.0) + 3.0], and proceed to the next step;

[0074] S1818. Judge T CS whether it is greater than T P . If so, obtain T CS = 0, and proceed to the next step. If not, proceed to the next step;

[0075] S1819. Judge T CS whether it is less than T P / 2.0. If so, obtain W COB = A COB . If not, obtain W COB = -A COB ;

[0076] S1820. Calculate the turbulence intensity and turbulence scale of the pebbles;

[0077] S1821. Exit the calculation of other turbulence subroutines.

[0078] Furthermore, the random number generation subroutine in S18 includes the following steps:

[0079] S1831. Enter the random number generation subroutine;

[0080] S1832. Generate 9 uniformly distributed random numbers between 0 and 1;

[0081] S1833. Generate 9 uniformly distributed random numbers between -0.5 and +0.5;

[0082] S1834. Save the current 9 normally distributed random numbers as the 9 normally distributed random numbers of the previous cycle;

[0083] S1835. Call the external retrieval subroutine and interpolate to obtain the normally distributed random numbers R U1 ~R U3 、R V1 ~R V3 、R W1 ~R W3 , where R U1 ~R U3 represents the U-direction normally distributed random number, R V1 ~R V3 represents the V-direction normally distributed random number, R W1 ~R W3 represents the W-direction normally distributed random number;

[0084] S1836. Exit the sub - process of generating random numbers.

[0085] Furthermore, the filter sub - process in S18 includes the following steps:

[0086] S1841. Enter the filter sub - process;

[0087] S1842. Calculate according to the time constant T1 = 2.0L U / V TR and the gain K = 2.8264η U L U / V TR where η U represents the turbulence intensity of the turbulent flow in the u - direction, and L U represents the scale of the turbulent flow in the u - direction;

[0088] S1843. Calculate according to EXP = e -Δt / T1 , the filter initial value CF0 = F U1 and the random input variable R f = R U1 where F UI represents the output of the longitudinal filter, I = 1, 2, 3;

[0089] S1844. Call the EQN200 sub - process to obtain F U1L = F U1 , F U1 = N F0, the gain K = 1.0, where F UIL represents the output of the longitudinal filter in the previous period, I = 1, 2, 3, and NF0 represents the filter output;

[0090] S1845. Calculate according to EXP = e -Δt / T1 , the filter initial value CF0 = F U2 and the random input variable R f = R U2 ;

[0091] S1846. Call the EQN200 sub - process to obtain F U2L = F U2 , F U2 = NF0, the time constant T1 = L U / V TR , the gain

[0092] S1847. Calculate according to EXP = e -Δt / T1 , CF0 = F U3 and R f = R U3 ;

[0093] S1848. Call the EQN200 sub - process to obtain F U3L = F U3 、F U3 = NF0, time constant T1 = 2.0L V / V TR 、Gain K = η U ·2.5133(L V / V TR ) 2 , where L V represents the scale of the turbulent flow in the v - direction;

[0094] S1849. According to EXP = e -Δt / T1 、Filter output value CF0 = F V1 and random input variable R f = R V1 for calculation, where F VI represents the lateral filter output, I = 1, 2, 3;

[0095] S1850. Call the EQN200 sub - process to obtain F V10L = F V1 、F V1 = N F0, time constant T1 = L V / V TR 、Gain K = 0.15915, where F VIL represents the lateral filter output of the previous cycle, I = 1, 2, 3;

[0096] S1851. According to EXP = e -Δt / T1 、Filter output value CF0 = F V2 、Filter output value DF0 of the previous cycle = F V2L 、Random input R f = R V2 and random input R of the previous cycle fL = R V2L for calculation;

[0097] S1852. Call the EQN300 sub - process to obtain F V2L = F V2 、F V2 = NF0, time constant 1T1 = L V / V TR 、Gain K = η V L V / 2V TR and time constant 2 where η V represents the turbulence intensity of the turbulent flow in the v - direction;

[0098] S1853. According to EXP = e-Δt / T1 、The filter output value CF0 = F V3 、The filter output value of the previous cycle DF0 = F V3L 、The random input R f = R V3 and the random input R of the previous cycle fL = R V3L are used for calculation;

[0099] S1854. Call the EQN400 sub - process to obtain F V3L = F V3 、F V3 = NF0, time constant T1 = 2.0L W / V TR and gain K = η W · 2.5133 (L W / V TR ) 2 where η W represents the turbulence intensity of the turbulence in the w - direction, and L W represents the scale of the turbulence in the w - direction;

[0100] S1855. According to EXP = e -Δt / T1 、the filter output CF0 = F W1 and the random number R f = R W1 are used for calculation, and F WI represents the vertical filter output, I = 1, 2, 3;

[0101] S1856. Call the EQN200 sub - process to obtain F W1L = F W1 、F W1 = NF0, time constant T1 = 2.0L W / V TR 、gain K = 0.15915, where F WIL represents the vertical filter output of the previous cycle, I = 1, 2, 3;

[0102] S1857. According to EXP = e -Δt / T1 、the filter output value CF0 = F W2 、the filter output value of the previous cycle DF0 = F W2L 、the random input R f = R W2 and the random input R of the previous cycle fL = R W2L are used for calculation;

[0103] S1858. Call the EQN300 sub - process to obtain F W2L = F W2 、FW2 = NF0, time constant 1T1 = L W / V TR , gain K = η W L W / 2V TR and time constant 2

[0104] S1859. According to EXP = e -Δt / T1 , filter output value CF0 = F W3 , filter output of the previous cycle DF0 = F W3L , random number R f = R W3 and random number R of the previous cycle fL = R W3L for calculation;

[0105] S1860. Call the EQN400 sub - process to obtain F W3L = F W3 and F W3 = N F0;

[0106] S1861. Exit the filter sub - process.

[0107] Furthermore, the sub - process for calculating the turbulence rate in S18 includes the following steps:

[0108] S1871. Enter the sub - process for calculating the turbulence rate;

[0109] S1872. Determine whether the turbulence type is vertical gust. If so, obtain the turbulence factor η THL = 1.0, and execute S1878. If not, execute S1873;

[0110] S1873. Determine whether the flight altitude h R is less than the turbulence low - altitude limit h L - 0.5. If so, execute S1875. If not, execute S1874;

[0111] S1874. Determine whether h R is greater than h U + 0.5. If so, execute S1875. If not, execute S1876;

[0112] S1875. Obtain the turbulence factor η THL = 0, and execute S1878;

[0113] S1876. Determine whether h R is less than or equal to h L . If so, obtain η THL = 1.0 - (h L-h R ) · 2.0, and execute S1878; if not, then execute S1877;

[0114] S1877. Judge h R whether it is greater than or equal to the turbulent height limit h U , if so, then obtain η THL = 1.0 + (h U - h R ) · 2.0, and execute S1878; if not, then obtain η THL = 1.0, and execute S1878;

[0115] S1878. Calculate U RAL = F U1 · F U2 + F U3 、U RAL = F V1 · F V2 + F V3 、w RAL =

[0116] F W1 · F W2 + F W3 · W JETUP + W COB 、P RAL = 0.025 · l · (W RAL - W JETUP ) / L W 、

[0117] Q RAL = 0.025 (W RAL - W JETUP ) · 81.5 / L W 、R RAL = 0.025 · 81.5 · V RAL / L W 、U T = U RAL · η THL 、

[0118] V T = V RAL · η THL 、W T = W RAL · η THL 、P T = P RAL · η THL 、Q T = Q RAL · η THL 、R T = R RAL · ηTHL , U TOL = U T , V TOL = V T and W TOL = W T , where W JETUP represents the vertical gust velocity, and W COB represents the cobblestone turbulence square wave intensity, η THL represents the turbulence factor, P T represents the pitch angular velocity output by the DRYDEN model, Q T represents the yaw angular velocity output by the DRYDEN model, R T represents the roll angular velocity output by the DRYDEN model, U TL represents the X-direction velocity output by the DRYDEN model in the previous period, V TL represents the Y-direction velocity output by the DRYDEN model in the previous period, W TL represents the Z-direction velocity output by the DRYDEN model in the previous period, represents; represents the Y-direction acceleration output by the DRYDEN model, represents the Z-direction acceleration output by the DRYDEN model;

[0119] S1879. Exit the subroutine for calculating the turbulence rate.

[0120] The beneficial effects of the present invention are as follows: By using the present invention, turbulence models such as rough atmosphere, cobblestone turbulence, and discrete vertical gusts can be provided, and the turbulence intensity is 10 levels, and the gust intensity can be selected within the range of 0 - 25 m / s, so as to simulate the atmospheric disturbance environment by calculating the turbulence disturbance velocity. In addition, the present invention can be classified into vertical gust turbulence, rough atmosphere, and cobblestone turbulence according to the turbulence type, and can be classified into mild turbulence, moderate turbulence, and severe turbulence according to the turbulence intensity for a total of nine kinds of turbulence simulations. BRIEF DESCRIPTION OF THE DRAWINGS

[0121] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0122] Figure 1 is a schematic diagram of the subroutine for clearing vertical gusts in a simulation method for simulating the turbulence environment of an aircraft according to an embodiment of the present invention;

[0123] Figure 2Schematic diagram of the calculation vertical gust sub - process in a simulation method for aircraft turbulence environment simulation according to an embodiment of the present invention;

[0124] Figure 3 Schematic diagram of the turbulence clearing sub - process in a simulation method for aircraft turbulence environment simulation according to an embodiment of the present invention;

[0125] Figure 4 Schematic diagram of the fade - out sub - process in a simulation method for aircraft turbulence environment simulation according to an embodiment of the present invention;

[0126] Figure 5 Schematic diagram of the random number generation sub - process in a simulation method for aircraft turbulence environment simulation according to an embodiment of the present invention;

[0127] Figure 6 Schematic diagram of the EQN200 sub - process in a simulation method for aircraft turbulence environment simulation according to an embodiment of the present invention;

[0128] Figure 7 Schematic diagram of the EQN300 sub - process in a simulation method for aircraft turbulence environment simulation according to an embodiment of the present invention;

[0129] Figure 8 Schematic diagram of the EQN400 sub - process in a simulation method for aircraft turbulence environment simulation according to an embodiment of the present invention. Detailed implementation manners

[0130] To further illustrate each embodiment, the present invention provides accompanying drawings. These drawings are part of the disclosure of the present invention, mainly used to illustrate the embodiments, and can be combined with the relevant descriptions in the specification to explain the operation principle of the embodiments. With reference to these contents, those of ordinary skill in the art should be able to understand other possible implementation manners and the advantages of the present invention. The components in the figures are not drawn to scale, and similar component symbols are usually used to represent similar components.

[0131] According to an embodiment of the present invention, a simulation method for aircraft turbulence environment simulation is provided.

[0132] Now, the present invention will be further described in combination with the accompanying drawings and specific implementation manners. As Figure 1-8 shown, the simulation method for aircraft turbulence environment simulation according to an embodiment of the present invention includes the following steps:

[0133] S1. Enter the turbulence module, and determine whether the aircraft speed reset flag is true. If it is, execute S19; if not, execute S2;

[0134] S2. Determine whether the aircraft flight freeze flag is true. If it is, execute S19; if not, execute S3;

[0135] S3. Determine whether the true airspeed V of the aircraft TR is greater than 15.0. If so, execute S4; if not, execute S19;

[0136] S4. Determine whether the reset directory of the aircraft is 0. If so, execute S5; if not, execute S19;

[0137] S5. Determine whether the wind profile index of the aircraft is not 0. If so, execute S7; if not, execute S6;

[0138] S6. Determine whether the microburst activation flag of the aircraft is true. If so, execute S7; if not, set the required rough atmosphere flag of the aircraft to false and execute S8;

[0139] S7. Set the required rough atmosphere flag of the aircraft to true and execute S8;

[0140] S8. Determine whether the rough atmosphere flag of the aircraft is true. If so, execute S9; if not, execute S10;

[0141] S9. Determine whether the turbulence type of the aircraft is 0. If so, execute S12; if not, execute S10;

[0142] S10. Determine whether the turbulence intensity of the aircraft is 0. If so, execute S12; if not, execute S11;

[0143] S11. Determine whether the repositioning directory of the aircraft is not 0. If so, execute S12; if not, execute S13;

[0144] S12. Determine whether the turbulence intensity factor η of the aircraft TU is 0. If so, execute S14; if not, execute S13;

[0145] S13. Determine whether the turbulence type of the aircraft is not 3. If so, call the clear vertical gust subroutine and execute S15; if not, call the calculate vertical gust subroutine and execute S15;

[0146] Among them, the clear vertical gust subroutine includes the following steps:

[0147] S1301. Enter the clear vertical gust subroutine;

[0148] S1302. Set both the gust distance factor and the gust timing to 0;

[0149] S1303. Set both the gust intensity factor and the vertical gust intensity to 0;

[0150] S1304. Set both the gust required cobblestone turbulence intensity and the vertical speed of the gust to 0;

[0151] S1305. Exit the vertical gust clearing sub - process.

[0152] The vertical gust calculation sub - process includes the following steps:

[0153] S1311. Enter the vertical gust calculation sub - process;

[0154] S1312. Determine whether the turbulence intensity is not 0. If so, obtain the vertical gust turbulence intensity factor where IAXTURIN represents the turbulence intensity. If not, obtain the vertical gust turbulence intensity factor η JETLOC = 0;

[0155] S1313. Determine the aircraft gust timing T JETUP whether it is less than 30.0. If so, obtain T JETUP = T JETUP + Δt, and execute S1314, where Δt represents the timing time interval. If not, obtain the distance factor = 0.008Mn + the distance factor, η JETLOC = 0, and execute S1315, where Mn represents the flight Mach number;

[0156] S1314. Determine whether T JETUP is greater than 40.0. If so, calculate the cobblestone turbulence intensity required for gust turbulence and execute S1317. If not, execute S1317;

[0157] S1315. Determine whether the distance factor is greater than or equal to 1.0. If so, execute S1316. If not, execute S1317;

[0158] S1316. Sequentially set the turbulence type number to 0, the turbulence intensity number to 0, the required gust flag to false, the light turbulence flag to false, the moderate turbulence flag to false, and the severe turbulence flag to false, and execute S1317;

[0159] S1317. Interpolate to calculate the vertical gust height factor η JUM = f(h1), where h1 represents the height;

[0160] S1318. Obtain the vertical gust velocity W JETUP = f(h2)·η JETLOC ·η JUM where h2 represents the distance factor;

[0161] S1319. Exit the vertical gust calculation sub - process.

[0162] S14. First call the turbulence clearing sub - process, then call the vertical gust clearing sub - process, and exit the turbulence module;

[0163] Among them, the turbulent flow elimination sub-process includes the following steps:

[0164] S141. Enter the turbulent flow elimination sub-process;

[0165] S142. Call the first attenuation sub-process (U RAL , U RA ) and the second attenuation sub-process (V RAL , V RA ), where U RAL represents the longitudinal component factor of the turbulent flow velocity in the aircraft body system, U RA represents the longitudinal component of the turbulent flow velocity in the aircraft body system, V RAL represents the lateral component factor of the turbulent flow velocity in the aircraft body system, V RA represents the lateral component of the turbulent flow velocity in the aircraft body system;

[0166] S143. Call the third attenuation sub-process (W RAL , W RA ) and the fourth attenuation sub-process (P RAL , P RA ), where W RAL represents the vertical component factor of the turbulent flow velocity in the aircraft body system, W RA represents the vertical component of the turbulent flow velocity in the aircraft body system, P RAL represents the pitch angular velocity factor generated by the turbulent flow velocity, P RA represents the pitch angular velocity generated by the turbulent flow velocity;

[0167] S144. Call the fifth attenuation sub-process (Q RAL , Q RA ) and the sixth attenuation sub-process (R RAL , R RA ), where Q RAL represents the yaw angular velocity factor generated by the turbulent flow velocity, Q RA represents the yaw angular velocity generated by the turbulent flow velocity, R RAL represents the roll angular velocity factor generated by the turbulent flow velocity, R RA represents the roll angular velocity generated by the turbulent flow velocity;

[0168] S145. Obtain η TU = η TUL = η TUD = 0;

[0169] S146. Clear all filter outputs F UI = 0, F VI = 0, F WI = 0, I = 1, 2, 3, where F UIRepresents the vertical filter output, F VI Represents the lateral filter output, F WI Represents the vertical filter output;

[0170] S147. Exit the turbulence clearing subroutine.

[0171] Specifically, the fading subroutine includes a first fading subroutine (U RAL , U RA ), a second fading subroutine (V RAL , V RA ), a third fading subroutine (W RAL , W RA ), a fourth fading subroutine (P RAL , P RA ), a fifth fading subroutine (Q RAL , Q RA ) and a sixth fading subroutine (R RAL , R RA ).

[0172] The fading subroutine includes the following steps:

[0173] Enter the fading subroutine and calculate using the formula RL = RAKFADER * RC, where RL represents U RAL , V RAL , W RAL , P RAL , Q RAL , R RAL , RAKFADER = 0.95, which is a fixed value of the factor representing the calculation speed, and RC represents U RA , V RA , W RA , P RA , Q RA , R RA ;

[0174] Judge whether |RL| is less than 0.01. If so, obtain RC = 0 and exit the fading subroutine. If not, obtain RC = RL and exit the fading subroutine.

[0175] S15. Judge whether the aircraft gust timing T JETUP is less than 30.0 seconds. If so, execute S17. If not, execute S16;

[0176] S16. Judge whether the aircraft turbulence type is 3. If so, execute S17. If not, execute S18;

[0177] S17. Call the turbulence clearing subroutine and exit the turbulence module;

[0178] S18. Call the other turbulence subroutine, the random number generation subroutine, the filter subroutine, and the turbulence rate calculation subroutine in sequence, and exit the turbulence module;

[0179] Among them, the other turbulence calculation subroutine includes the following steps:

[0180] S1801. Enter the other turbulence calculation subroutine;

[0181] S1802. Determine whether the rough atmosphere requirement flag is true. If it is, calculate the turbulence intensity factor and execute S1806. If not, execute S1803;

[0182] S1803. Determine whether the turbulence type is rough atmosphere. If it is, obtain the turbulence intensity factor required for the rough atmosphere where IAXTURIN represents the turbulence intensity, and execute S1806. If not, execute S1804;

[0183] S1804. Determine whether the turbulence type is pebble turbulence. If it is, obtain the turbulence intensity factor required for the pebble atmosphere and execute S1805. If not, obtain η TUD = η JETLOC , and execute S1806;

[0184] S1805. Obtain the square wave intensity and execute S1806;

[0185] S1806. Determine whether the turbulence intensity number is 0. If it is, execute S1807. If not, execute 1808;

[0186] S1807. Set the rough atmosphere requirement flag to true. If it has been set, obtain η TUD = 0, and execute S1808. If it has not been set, execute S1808;

[0187] S1808. Obtain the turbulence intensity factor increment Δη TU = η TUD - η TUL , and determine whether Δη TU is greater than 0.01. If it is, obtain Δη TU = 0.01, and execute S1810. If not, execute S1809; where η TUL is the previous cycle's turbulence intensity factor;

[0188] S1809. Determine whether Δη TU is less than -0.05. If it is, obtain Δη TU = -0.05. If not, execute S1810;

[0189] S1810. Obtain η successively to get η TU = η TUL + Δη TU and η TUL = η TU , where η TUL represents the turbulent intensity factor of the previous cycle, and Δη TU represents the increment of the turbulent intensity factor;

[0190] S1811. Calculate the turbulent intensity height factor and the turbulent scale height factor by interpolation, and determine whether the turbulent type is rough atmosphere. If so, execute S1813; if not, execute S1812;

[0191] S1812. Determine whether the rough atmosphere requirement flag is true. If so, execute S1813; if not, execute S1814;

[0192] S1813. Calculate the rough atmosphere turbulent intensity and turbulent scale successively to obtain the clear cobblestone square wave intensity W COB = 0, A COB = 0, and execute S1821;

[0193] S1814. Determine whether the turbulent type is cobblestone turbulence. If so, execute S1816; if not, execute S1815;

[0194] S1815. Determine whether the turbulent type is vertical gust. If so, execute S1816; if not, execute S1821;

[0195] S1816. Obtain the cobblestone timing T CS = T CS + Δt, and determine whether the aircraft true airspeed V TR is greater than 700.0. If so, obtain Δη TU = 0.01, and execute S1818; if not, execute S1817;

[0196] S1817. Determine whether the true airspeed V TR is less than 200.0. If so, obtain the rough atmosphere timing T P = 0.333333, and execute the next step; if not, obtain T P = 1 / [0.006(V TR - 200.0)+ 3.0], and execute the next step;

[0197] S1818. Determine whether T CS is greater than T P . If so, obtain T CS = 0, and execute the next step; if not, execute the next step;

[0198] S1819, Determine T CS Whether it is less than T P / 2.0. If so, obtain W COB = A COB If not, obtain W COB = -A COB ;

[0199] S1820, Calculate the turbulence intensity and turbulence scale of pebbles;

[0200] S1821, Exit the calculation of other turbulence subroutines.

[0201] The random number generation subroutine includes the following steps:

[0202] S1831, Enter the random number generation subroutine;

[0203] S1832, Generate 9 uniformly distributed random numbers between 0 and 1;

[0204] S1833, Generate 9 uniformly distributed random numbers between -0.5 and +0.5;

[0205] S1834, Store the current 9 normally distributed random numbers as the 9 normally distributed random numbers of the previous cycle;

[0206] S1835, Call the external retrieval subroutine and interpolate to obtain the normally distributed random numbers R U1 ~R U3 、R V1 ~R V3 、R W1 ~R W3 ,where R U1 ~R U3 represents the U-direction normally distributed random number, R V1 ~R V3 represents the V-direction normally distributed random number, R W1 ~R W3 represents the W-direction normally distributed random number;

[0207] S1836, Exit the random number generation subroutine.

[0208] The filter subroutine includes the following steps:

[0209] S1841, Enter the filter subroutine;

[0210] S1842, According to the time constant T1 = 2.0L U / V TR and the gain K = 2.8264η U L U / V TRPerform calculations, where η U represents the turbulence intensity of the turbulent flow in the u direction, and L U represents the scale of the turbulent flow in the u direction;

[0211] S1843. According to EXP = e -Δt / T1 , filter initial value CF0 = F U1 and random input variable R f = R U1 Perform calculations;

[0212] S1844. Call the EQN200 subroutine to obtain F U1L = F U1 , F U1 = NF0, gain K = 1.0, where F UIL represents the output of the longitudinal filter in the previous period, I = 1, 2, 3, and NF0 represents the filter output;

[0213] S1845. According to EXP = e -Δt / T1 , filter initial value CF0 = F U2 and random input variable R f = R U2 Perform calculations;

[0214] S1846. Call the EQN200 subroutine to obtain F U2L = F U2 , F U2 = NF0, time constant T1 = L U / V TR , gain

[0215] S1847. According to EXP = e -Δt / T1 , CF0 = F U3 and R f = R U3 Perform calculations;

[0216] S1848. Call the EQN200 subroutine to obtain F U3L = F U3 , F U3 = NF0, time constant T1 = 2.0L V / V TR , gain K = η U ·2.5133(L V / V TR ) 2 , where L V represents the scale of the turbulent flow in the v direction;

[0217] S1849. According to EXP = e -Δt / T1, the filter output value CF0 = F V1 and the random input variable R f = R V1 are calculated, where F VI represents the lateral filter output, I = 1, 2, 3;

[0218] S1850. Call the EQN200 subroutine to obtain F V10L = F V1 , F V1 = NF0, the time constant T1 = L V / V TR , the gain K = 0.15915, where F VIL represents the previous cycle's lateral filter output, I = 1, 2, 3;

[0219] S1851. According to EXP = e -Δt / T1 , the filter output value CF0 = F V2 , the previous cycle's filter output value DF0 = F V2L , the random input R f = R V2 and the previous cycle's random input R fL = R V2L are calculated;

[0220] S1852. Call the EQN300 subroutine to obtain F V2L = F V2 , F V2 = NF0, the time constant 1T1 = L V / V TR , the gain K = η V L V / 2V TR and the time constant 2 η V represents the turbulence intensity of the turbulence in the v direction;

[0221] S1853. According to EXP = e -Δt / T1 , the filter output value CF0 = F V3 , the previous cycle's filter output value DF0 = F V3L , the random input R f = R V3 and the previous cycle's random input R fL = R V3L are calculated;

[0222] S1854. Call the EQN400 subroutine to obtain F V3L = F V3 , F V3 = NF0, the time constant T1 = 2.0L W / VTR and gain K = η W ·2.5133(L W / V TR ) 2 , where η W represents the turbulence intensity of the turbulent flow in the w direction, and L W represents the scale of the turbulent flow in the w direction;

[0223] S1855. According to EXP = e -Δt / T1 , filter output CF0 = F W1 and random number R f = R W1 to perform calculations, where F WI represents the vertical filter output, and I = 1, 2, 3;

[0224] S1856. Call the EQN200 subroutine to obtain F W1L = F W1 , F W1 = NF0, time constant T1 = 2.0L W / VTR, gain K = 0.15915, where F WIL represents the vertical filter output of the previous cycle, and I = 1, 2, 3;

[0225] S1857. According to EXP = e -Δt / T1 , filter output value CF0 = F W2 , filter output value DF0 of the previous cycle = F W2L , random input R f = R W2 and random input R of the previous cycle fL = R W2L to perform calculations;

[0226] S1858. Call the EQN300 subroutine to obtain F W2L = F W2 , F W2 = NF0, time constant 1T1 = L W / V TR , gain K = η W L W / 2V TR and time constant 2

[0227] S1859. According to EXP = e -Δt / T1 , filter output value CF0 = F W3 , filter output DF0 of the previous cycle = F W3L , random number R f = R W3 and random number R of the previous cycle fL=R W3L Perform calculations;

[0228] S1860. Call the EQN400 sub - process to obtain F W3L =F W3 and F W3 =NF0;

[0229] S1861. Exit the filter sub - process.

[0230] Specifically, the EQN200 sub - process includes the following steps:

[0231] Enter the EQN200 sub - process;

[0232] Calculate NF0 = CF0e -Δt / T1 +R f K(1 - e -Δt / T1 );

[0233] Judge whether |NF0| is less than 0.00001. If so, obtain NF0 = 0. If not, execute the next step;

[0234] Exit the EQN200 sub - process.

[0235] The EQN300 sub - process includes the following steps:

[0236] Enter the EQN300 sub - process;

[0237] Calculate NF0=(2CF0 - DF0·EXP)·EXP+Δt·[K·EXP(R f -R fL ) / T1 2 ;

[0238] Judge whether |NF0| is less than 0.00001. If so, obtain NF0 = 0. If not, execute the next step;

[0239] Exit the EQN300 sub - process.

[0240] The EQN400 sub - process includes the following steps:

[0241] Enter the EQN400 sub - process;

[0242] Calculate F1=(2.0·CF0 - DF0·EXP)·EXP, F 21 =KT2 / T1 2 , F 22 =T1 2 / T2, F 23 =(ΔtT1 / T 2) -Δt, F 31 =R f ·F21 ·[F 22 +(-F 23 -F 22 )·EXP], F 32 =R fL ·F 21 ·[F 22 (EXP - 1.0)+F 23 ·EXP, NF0 = F1 + F 31 +F 32 ;

[0243] Judge whether |NF0| is less than 0.00001. If so, then get NF0 = 0. If not, then execute the next step;

[0244] Exit the EQN400 sub - process.

[0245] The sub - process of calculating the turbulence rate includes the following steps:

[0246] S1871. Enter the sub - process of calculating the turbulence rate;

[0247] S1872. Judge whether the turbulence type is vertical gust. If so, then get the turbulence factor η THL =1.0 and execute S1878. If not, then execute S1873;

[0248] S1873. Judge whether the flight altitude h R is less than the low - altitude limit of turbulence h L -0.5. If so, then execute S1875. If not, then execute S1874;

[0249] S1874. Judge whether h R is greater than the high - altitude limit of turbulence h U +0.5. If so, then execute S1875. If not, then execute S1876;

[0250] S1875. Get the turbulence factor η THL =0 and execute S1878;

[0251] S1876. Judge whether h R is less than or equal to h L , if so, then get η THL =1.0-(h L -h R )·2.0 and execute S1878. If not, then execute S1877;

[0252] S1877. Judge whether h R is greater than or equal to h U , if so, then get η THL =1.0+(hU -h R )·2.0, and execute S1878, if not, then get η THL =1.0, and execute S1878;

[0253] S1878, calculate U in sequence RAL =F U1 ·F U2 +F U3 、V RAL =F V1 ·F V2 +F V3 , W RAL =

[0254] F W1 ·F W2 +F W3 ·W JETUP +W COB , P RAL =0.025·l·(W RAL -W JETUP ) / L W ,

[0255] Q RAL =0.025(W RAL -W JETUP )·81.5 / L W , R RAL =0.025·81.5·V RAL / L W , U T =U RAL ·η THL ,

[0256] V T =V RAL ·η THL , W T =W RAL ·η THL , P T =P RAL ·η THL , Q T =Q RAL ·η THL , R T =R RAL ·η THL , U TOL =U T 、V TOL =V T and W TOL =W T , where W JETUP Indicates the vertical gust speed, WCOB Represents the intensity of the pebble turbulent square wave, P T Represents the pitch angular velocity output by the DRYDEN model, Q T Represents the yaw angular velocity output by the DRYDEN model, R T Represents the roll angular velocity output by the DRYDEN model, U TL Represents the X - direction velocity output by the DRYDEN model in the previous cycle, V TL Represents the Y - direction velocity output by the DRYDEN model in the previous cycle, W TL Represents the Z - direction velocity output by the DRYDEN model in the previous cycle; Represents the Y - direction acceleration output by the DRYDEN model, Represents the Z - direction acceleration output by the DRYDEN model;

[0257] S1879. Exit the sub - process of calculating the turbulent flow rate.

[0258] S19. Judge the aircraft re - positioning timing T REP Whether it is not zero. If so, execute S21; if not, execute S20;

[0259] S20. Judge the aircraft true airspeed V TR Whether it is less than 15.0. If so, execute S21; if not, exit the turbulent flow module;

[0260] S21. Obtain the clear turbulent flow velocity U T =V T =W T =0, η TU =η TUL =η TUD =0, exit the turbulent flow module;

[0261] Wherein, U T Is the X - direction velocity output by the DRYDEN model, V T Is the Y - direction velocity output by the DRYDEN model, W T Is the Z - direction velocity output by the DRYDEN model, η TUL Is the turbulent flow intensity factor in the previous cycle.

[0262] For the convenience of understanding the above - mentioned technical solution of the present invention, the following will detail the relevant calculation formulas in the turbulent flow module of the present invention.

[0263] The relevant calculation formulas of the turbulent flow module are as follows:

[0264] Vvi = interp(Ttn)

[0265] Vhf = interp(Hg)

[0266] Hfu = erp(Hg) * Vvi

[0267] Hfv = erp(Hg) * Vvi - Ran * (2.8284 * Hfu * Cfu / Vr) * (1 - e -Tt / 2.0*Cfu / Vr) )

[0268] Hfw = erp(Hg) * Vvi + Ran * (1 - e -Tt / (2.0*Cfu / Vr) )

[0269]

[0270] Cfv = erp(Hg) Ran * (2.5133 * Hfv * (Cfv / Vr) 2 ) * (1 - e -Tr / 2.0*Cfv / Vr) )

[0271] Ffv2 = (2 * LFf2 - Yffv2 * e -Tt / (2.0*Cfv / Vr) ) * e -Tt / (2.0*Cfv / Vr) +0.15915 / (2.0 * Cfv / Vr) 2 * Tt * e -Tt / (2.0*Cfv / r) *(Lffv2 - Yffv2)

[0272] Ffv3 = (4 * Lffv3 - Yffv3 * e -Tt / (2.0*Cfv / Vr) ) * e -Tt / (2.0*Cfv / Vr) +Yffv3 * Hfv * Cfv / Vr / 2+LRan * Hfv * Cfv / Vr / 2 * ((e -Tt(2.0*Cfv / Vr) -1)+Tt / 0.74205) * e -Tt / (2.0*Cfv / Vr)

[0273] Ffw1 = LFfw1 + e -Tt / (2.0*Cfw / Vr) +Ran * (2.5133 * Hfw * (Cfw / Vr) 2 ) * (1 - e -Tt / (2.0*Cfw / Vr) )

[0274] Ffw2 = (2 * LFfw2 - Yffw2 * e -Tt / (2.0*Cfw / Vr) ) * e -Tt / (2.0*Cfw / Vr) +0.15915 / (2.0 * Cfw / Vr) 2 * Tt * e -Tt / (2.0*Cfw / Vr) *(Lffw2 - Yffw2)

[0275] Ffw3 = (4 * Lffw3 - Yffw3 * e -Tt / (2.0*Cfw / r) ) * e -Tt / (2.0*Cfw / Vr) +Yffw3 * Hfw * Cfw / Vr / 2+LRan * Hfwv * Cfw / Vr / 2 * ((e -Tt / (2.0*Cfw / Vr)-1)+Tt / 0.74205)*e -Tt / (2.0*Cfw / Vr)

[0276] Ssu = Ffu1 * Ffu2 + Ffu3

[0277] Ssv = Ffv1 * Ffv2 + Ff3

[0278] Ssw = Ffw1 * Ffw2 + Ffw3 + Vvv + Csi

[0279] Shp = 0.025 * Ws * (Ssw - Vvv)

[0280] Shq = 0.025 * Wt * (Ssw - Vvv)

[0281] Shr = 0.025 * Wt * Ssv

[0282] Sau = (Ssu - Lssu) / Tt

[0283] Sav = (Ssv - Lssv) / Tt

[0284] Saw = (Ssw - Lssw) / Tt

[0285] Among them, the corresponding table of variables in the turbulent flow module formula is as follows:

[0286]

[0287]

[0288] In summary, by means of the above technical solutions of the present invention, through the use of the present invention, turbulent flow models such as rough atmosphere, pebble turbulent flow, and discrete vertical gusts can be provided, and the turbulent flow intensity is 10 levels, and the gust intensity can be selected within the range of 0 - 25 m / s. Thus, the atmospheric disturbance environment can be simulated by calculating the turbulent disturbance velocity. In addition, the present invention can be divided into vertical gust turbulent flow, rough atmosphere, and pebble turbulent flow according to the turbulent flow type, and can be divided into mild turbulent flow, moderate turbulent flow, and severe turbulent flow according to the turbulent flow intensity, with a total of nine kinds of turbulent flow simulations.

[0289] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present invention shall be included within the protection scope of the present invention.

Claims

1. A simulation method for aircraft turbulence environment simulation, characterized in that, The method comprises the following steps: S1. Enter the turbulence module, and determine whether the aircraft speed reset flag is true. If it is, execute S19; if not, execute S2; S2. Determine whether the aircraft flight freeze flag is true. If it is, execute S19; if not, execute S3; S3. Determine whether the true airspeed V of the aircraft TR is greater than 15.

0. If so, execute S4; if not, execute S19; S3. Determine whether the aircraft reset directory is 0. If it is, execute S4; if not, execute S19; S4. Determine whether the aircraft wind profile index is not 0. If it is, execute S7; if not, execute S5; S5. Determine whether the aircraft microburst activation flag is true. If it is, execute S7; if not, set the aircraft required rough atmosphere flag to false and execute S6; S6. Set the aircraft required rough atmosphere flag to true and execute S7; S7. Set the aircraft required rough atmosphere flag to true and execute S8; S8. Determine whether the aircraft rough atmosphere flag is true. If it is, execute S9; if not, execute S10; S9. Determine whether the aircraft turbulence type is 0. If it is, execute S12; if not, execute S10; S10. Determine whether the aircraft turbulence intensity is 0. If it is, execute S12; if not, execute S11; S12. Determine the aircraft turbulence intensity factor η TU Whether it is 0. If so, execute S14; if not, execute S13; S11. Determine whether the aircraft repositioning directory is not 0. If it is, execute S12; if not, execute S13; S13. Determine whether the aircraft turbulence type is not 3. If it is, call the clear vertical gust sub - process and execute S15; if not, call the calculate vertical gust sub - process and execute S15; S15. Determine the aircraft gust timing T JETUP Whether it is less than 30.0 seconds. If so, execute S17; if not, execute S16; S14. First call the clear turbulence sub - process, then call the clear vertical gust sub - process, and exit the turbulence module; S15. Determine whether the aircraft turbulence type is 3. If it is, execute S16; if not, execute S17; S16. Call the clear turbulence sub - process and exit the turbulence module; S17. Call the clear turbulence sub - process and exit the turbulence module; S18. Call the calculate other turbulence sub - process, call the generate random number sub - process, call the filter sub - process, call the calculate turbulence rate sub - process in sequence, and exit the turbulence module; The calculate other turbulence sub - process in S18 comprises the following steps: S1803. Determine whether the turbulence type is a rough atmosphere. If so, obtain the turbulence intensity factor required for a rough atmosphere. where IAXTURIN represents the turbulence intensity, and execute S1806. If not, execute S1804. S1804. Determine whether the turbulence type is pebble turbulence. If so, obtain the turbulence intensity factor required for the pebble atmosphere and execute S1805. If not, obtain η and execute S1806. Here, η TUD = η JETLOC , and execute S1806, where η TUD represents the turbulence intensity factor, and η JETLOC represents the vertical gust turbulence intensity factor; S1805. Obtain the square wave intensity And execute S1806; S1801. Enter the calculate other turbulence sub - process; S1807. Set the rough atmosphere requirement flag to true. If it has already been set, obtain η TUD = 0 and execute S1808. If it has not been set, execute S1808; S1808. Obtain the increment of the turbulence intensity factor Δη TU = η TUD - η TUL , and determine whether Δη TU is greater than 0.

01. If so, obtain Δη TU = 0.01 and execute S1810. If not, execute S1809; where η TUL is the turbulence intensity factor in the previous cycle. S1809. Determine Δη TU Whether it is less than -0.

05. If so, obtain Δη TU = -0.

05. If not, execute S1810; S1810, successively obtain η TU = η TUL + Δη TU and η TUL = η TU ; where η TUL represents the turbulent intensity factor of the previous cycle, and Δη TU represents the increment of the turbulent intensity factor; S1802. Determine whether the rough atmosphere requirement flag is true. If it is, calculate the turbulence intensity factor and execute S1806; if not, execute S1803; S1803. Determine whether the turbulence intensity number is 0. If it is, execute S1804; if not, execute S1805; S1813. Calculate the rough atmospheric turbulence intensity and turbulence scale in sequence to obtain the clear cobblestone square wave intensity W COB = 0, A COB = 0, and execute S1821; S1804. Calculate the turbulence intensity height factor using interpolation; S1805. Calculate the turbulence scale height factor using interpolation; S1816. Obtain the cobblestone timing T CS = T CS + Δt, and determine whether the true airspeed V of the aircraft TR is greater than 700.

0. If so, obtain Δη TU = 0.01, and execute S1818. If not, execute S1817; S1817. Determine the true airspeed V TR Check if it is less than 200.

0. If so, obtain the rough air time T P = 0.333333, and proceed to the next step. If not, obtain T P = 1 / [0.006(V TR - 200.0) + 3.0], and proceed to the next step; S1818, Determine T CS Whether it is greater than T P , if so, obtain T CS = 0, and execute the next step, if not, execute the next step; S1819, Determine T CS Is it less than T P / 2.0, if so, then obtain W COB = A COB , if not, then obtain W COB = -A COB ; S1806. Determine whether the turbulence intensity number is 0. If it is, execute S1807; if not, execute S1808; S1807. Determine whether the turbulence type is rough atmosphere. If it is, execute S1809; if not, execute S1810; S1808. Determine whether the turbulence type is rough atmosphere. If it is, execute S1809; if not, execute S1810; S1809. Calculate the turbulence intensity and turbulence scale; S1810. Determine whether the turbulence type is pebble turbulence. If it is, execute S1811; if not, execute S1812; S1811. Calculate the turbulence intensity height factor and turbulence scale height factor using interpolation, and determine whether the turbulence type is rough atmosphere. If it is, execute S1813; if not, execute S1812; S1812. Determine whether the rough atmosphere requirement flag is true. If it is, execute S1813; if not, execute S1814; S1813. Determine whether the turbulence type is pebble turbulence. If it is, execute S1815; if not, execute S1816; S1814. Determine whether the turbulence type is pebble turbulence. If it is, execute S1815; if not, execute S1816; S1815. Determine whether the turbulence type is vertical gust. If it is, execute S1816; if not, execute S1820; S1816. Determine whether the turbulence type is vertical gust. If it is, execute S1817; if not, execute S1818; S1817. Calculate the vertical gust intensity; S1818. Determine whether the turbulence type is vertical gust. If it is, execute S1819; if not, execute S1820; S1819. Calculate the vertical gust scale; S1820. Calculate the pebble turbulence intensity and turbulence scale; S1821. Exit the calculate other turbulence sub - process; S19. Determine whether the aircraft repositioning timing T REP is not zero. If so, execute S21; if not, execute S20; S20. Determine whether the true airspeed V of the aircraft TR is less than 15.

0. If so, execute S21; if not, exit the turbulence module. S21. Sequentially obtain the clear turbulent flow velocities U T = V T = W T = 0, η TU = η TUL = η TUD = 0, and exit the turbulent flow module; Among them, U T is the X-direction velocity output by the DRYDEN model, V T is the Y-direction velocity output by the DRYDEN model, W T is the Z-direction velocity output by the DRYDEN model, η TUL is the turbulence intensity factor of the previous cycle.

2. A simulation method for aircraft turbulence environment simulation according to claim 1, characterized in that, The vertical gust clearing sub-process in S13 includes the following steps: S1301. Enter the vertical gust clearing sub-process; S1302. Set both the gust distance factor and the gust timing to 0; S1303. Set both the gust intensity factor and the vertical gust intensity to 0; S1304. Set both the gust-required pebble turbulence intensity and the vertical velocity of the gust to 0; S1305. Exit the vertical gust clearing sub-process.

3. A simulation method for aircraft turbulence environment simulation according to claim 1, characterized in that, The vertical gust calculation sub-process in S13 includes the following steps: S1311. Enter the vertical gust calculation sub-process; S1312. Determine whether the turbulence intensity is not zero. If so, obtain the vertical gust turbulence intensity factor where IAXTURIN represents the turbulence intensity. If not, obtain the vertical gust turbulence intensity factor η JETLOC = 0; S1313. Determine the aircraft gust timing T JETUP Whether it is less than 30.

0. If so, obtain T JETUP = T JETUP + Δt, and execute S1314, where Δt represents the timing time interval. If not, obtain the distance factor = 0.008Mn + the distance factor, η JETLOC = 0, and execute S1315, where Mn represents the flight Mach number; S1314. Determine T JETUP Check if it is greater than 40.

0. If so, calculate the cobblestone turbulence intensity required for gust turbulence and execute S1317. If not, execute S1317; S1315. Determine whether the distance factor is greater than or equal to 1.

0. If so, execute S1316; if not, execute S1317; S1316. Sequentially set the turbulence type number to 0, the turbulence intensity number to 0, the required gust flag to false, the mild turbulence flag to false, the moderate turbulence flag to false, and the severe turbulence flag to false, and then execute S1317; S1317. Interpolation calculation of the vertical gust height factor η JUM = f(h1), where h1 represents the height; S1318. Obtain the vertical gust velocity W JETUP = f(h2)·η JETLOC ·η JUM , where h2 represents the distance factor; S1319. Exit the vertical gust calculation sub-process.

4. A simulation method for aircraft turbulence environment simulation according to claim 1, characterized in that, The turbulence clearing sub-process in S14 includes the following steps: S141. Enter the turbulence clearing sub-process; S142. Call the first fading sub-process (U RAL , U RA ) and the second fading sub-process (V RAL , V RA ), where U RAL represents the longitudinal component factor of the turbulent wind speed in the aircraft body system, and U RA represents the longitudinal component of the turbulent wind speed in the aircraft body system. V RAL represents the lateral component factor of the turbulent wind speed in the aircraft body system, and V RA represents the lateral component of the turbulent wind speed in the aircraft body system; S143. Call the third attenuation sub-process (W RAL ,W RA ) and the fourth attenuation sub-process (P RAL ,P RA ), where W RAL represents the vertical component factor of the turbulent wind speed in the aircraft body system, and W RA represents the vertical component of the turbulent wind speed in the aircraft body system, and P RAL represents the pitch angular velocity factor generated by the turbulent wind speed, and P RA represents the pitch angular velocity generated by the turbulent wind speed; S144. Call the fifth fade-out sub-process (Q RAL , Q RA ) and the sixth fade-out sub-process (R RAL , R RA ), where Q RAL represents the yaw angular velocity factor generated by the turbulent wind speed, and Q RA represents the yaw angular velocity generated by the turbulent wind speed, and R RAL represents the roll angular velocity factor generated by the turbulent wind speed, and R RA represents the roll angular velocity generated by the turbulent wind speed; S145. Obtain η TU = η TUL = η TUD = 0; S146. Clear all filter outputs F UI = 0, F VI = 0, F WI = 0, I = 1, 2, 3, where F UI represents the longitudinal filter output, F VI represents the lateral filter output, F WI represents the vertical filter output; S147. Exit the turbulence clearing sub-process.

5. A simulation method for aircraft turbulence environment simulation according to claim 4, characterized in that, The weakening sub-process includes the following steps: Enter the fade-out sub-process and calculate using the formula RL = RAKFADER * RC, where RL represents U RAL , V RAL , W RAL , P RAL , Q RAL , R RAL , RAKFADER = 0.95, which is the fixed value of the factor representing the calculation speed, and RC represents U RA , V RA , W RA , P RA , Q RA , R RA ; Determine whether |RL| is less than 0.

01. If so, obtain RC = 0 and exit the weakening sub-process; if not, obtain RC = RL and exit the weakening sub-process.

6. A simulation method for aircraft turbulent environment simulation according to claim 1, characterized in that, The random number generation sub-process in S18 includes the following steps: S1831. Enter the random number generation sub-process; S1832. Generate 9 uniformly distributed random numbers between 0 and 1; S1833. Generate 9 uniformly distributed random numbers between -0.5 and +0.5; S1834. Store the current 9 normally distributed random numbers as the 9 normally distributed random numbers of the previous cycle; S1835, call external search subroutine and interpolate to obtain normal distribution random number R U1 ~R U3 , R V1 ~R V3 , R W1 ~R W3 , where R U1 ~R U3 U represents a normally distributed random number, R V1 ~R V3 V represents a normally distributed random number, R W1 ~R W3 Represents W-direction normally distributed random numbers; S1836. Exit the random number generation sub-process.

7. A simulation method for aircraft turbulent environment simulation according to claim 6, characterized in that, The filter sub-process in S18 includes the following steps: S1841. Enter the filter sub-process; S1842. Calculate according to the time constant T1 = 2.0L U / V TR and the gain K = 2.8264η U L U / V TR where η U represents the turbulence intensity of the turbulent flow in the u direction, and L U represents the scale of the turbulent flow in the u direction; S1843. Calculate according to EXP = e -Δt / T1 , the initial value of the filter CF0 = F U1 , and the random input variable R f = R U1 , where F UI represents the output of the vertical filter, I = 1, 2, 3; S1844. Call the EQN200 sub - process to obtain F U1L = F U1 、F U1 = NF0, gain K = 1.0, where F UIL represents the output of the longitudinal filter in the previous cycle, I = 1, 2, 3, and NF0 represents the filter output; S1845. Calculate according to EXP = e -Δt / T1 , the initial value of the filter CF0 = F U2 , and the random input variable R f = R U2 ; S1846. Call the EQN200 sub - process to obtain F U2L = F U2 、F U2 = NF0, time constant T1 = L U / V TR 、Gain S1847. Calculate according to EXP = e -Δt / T1 , CF0 = F U3 and R f = R U3 ; S1848. Call the EQN200 sub - process to obtain F U3L = F U3 、F U3 = NF0, time constant T1 = 2.0L V / V TR 、Gain K = η U ·2.5133(L V / V TR ) 2 , where L V represents the scale of the turbulent flow in the v - direction; S1849. Calculate according to EXP = e -Δt / T1 , filter output value CF0 = F V1 and random input variable R f = R V1 , where F VI represents the lateral filter output, I = 1, 2, 3; S1850. Call the EQN200 sub - process to obtain F V1L = F V1 、F V1 = NF0, time constant T1 = L V / V TR 、Gain K = 0.15915, where F VIL represents the output of the lateral filter in the previous cycle, I = 1, 2, 3; S1851. Calculate based on EXP = e -Δt / T1 and the filter output value CF0 = F V2 and the filter output value DF0 = F of the previous cycle V2L and the random input R f = R V2 and the random input R of the previous cycle fL = R V2L ; S1852. Call the EQN300 sub-process to obtain F V2L = F V2 、F V2 = NF0, time constant 1T1 = L V / V TR 、Gain K = η V L V / 2V TR and time constant 2 where η V represents the turbulence intensity of the turbulent flow in the v direction; S1853. Calculate according to EXP = e -Δt / T1 The filter output value CF0 = F V3 The filter output value DF0 of the previous cycle = F V3L The random input R f = R V3 And the random input R of the previous cycle fL = R V3L Perform calculations; S1854, call sub - process EQN400 to obtain F V3L = F V3 , F V3 = NF0, time constant T1 = 2.0L W / V TR and gain K = η W ·2.5133(L W / V TR ) 2 , where η W represents the turbulence intensity of the turbulent flow in the w - direction, and L W represents the scale of the turbulent flow in the w - direction; S1855, according to EXP = e -Δt / T1 and the filter output CF0 = F W1 and the random number R f = R W1 perform calculations, F WI represents the vertical filter output, I = 1, 2, 3; S1856. Call the EQN200 sub - process to obtain F W1L = F W1 、F W1 = NF0, time constant T1 = 2.0L W / V TR 、Gain K = 0.15915, where F WIL represents the output of the vertical filter in the previous cycle, I = 1, 2, 3; S1857. Calculate based on EXP = e -Δt / T1 The filter output value CF0 = F W2 The filter output value DF0 of the previous cycle = F W2L The random input R f = R W2 And the random input R of the previous cycle fL = R W2L Perform calculations; S1858, call the EQN300 sub-process to obtain F W2L = F W2 、F W2 = NF0, time constant 1T1 = L W / V TR 、gain K = η W L W / 2V TR and time constant 2 S1859. Calculate according to EXP = e -Δt / T1 and the filter output value CF0 = F W3 and the filter output DF0 of the previous cycle = F W3L and the random number R f = R W3 and the random number R of the previous cycle fL = R W3L for calculation; S1860. Call the EQN400 sub - process to obtain F W3L = F W3 and F W3 = NF0; S1861. Exit the filter sub-process.

8. A simulation method for aircraft turbulent environment simulation according to claim 4, characterized in that, The turbulence rate calculation sub-process in S18 includes the following steps: S1871. Enter the turbulence rate calculation sub-process; S1872. Determine whether the type of turbulence is a vertical gust. If so, obtain the turbulence factor η THL = 1.0, and execute S1878. If not, execute S1873; S1873. Determine the flight altitude h R Whether it is less than the low-turbulence altitude limit h L -0.

5. If so, execute S1875; if not, execute S1874; S1874. Determine h R Is it greater than h U + 0.

5. If so, execute S1875; if not, execute S1876; S1875. Obtain the turbulence factor η THL = 0, and execute S1878; S1876. Determine h R Whether it is less than or equal to h L . If so, obtain η THL = 1.0 - (h L - h R )·2.0, and execute S1878. If not, execute S1877; S1877. Determine h R Is it greater than or equal to the turbulent height limit h U ? If so, obtain η THL = 1.0 + (h U - h R ) · 2.0, and execute S1878. If not, obtain η THL = 1.0, and execute S1878; S1878. Calculate U successively RAL = F U1 ·F U2 + F U3 、V RAL = F V1 ·F V2 + F V3 、W RAL = F W1 ·F W2 +F W3 ·W JETUP +W COB 、P RAL =0.025·l·(W RAL -W JETUP ) / L W 、 Q RAL = 0.025 (W RAL - W JETUP )·81.5 / L W 、R RAL = 0.025·81.5·V RAL / L W 、U T = U RAL ·η THL 、 V T =V RAL ·η THL , W T =W RAL ·η THL , P T =P RAL ·η THL , Q T =Q RAL ·η THL , R T =R RAL ·η THL , U TOL =U T 、V TOL =V T and W TOL =W T , where W JETUP Indicates the vertical gust speed, W COB represents the intensity of the pebble turbulence square wave, η THL represents the turbulence factor, P T Denotes the pitch velocity output by the DRYDEN model, Q T Indicates the DRYDEN model output yaw rate, R T Indicates the roll angular velocity output by the DRYDEN model, U TL Indicates the X-direction velocity output by the DRYDEN model in the previous cycle, V TL Indicates the Y-direction velocity output by the DRYDEN model in the previous cycle, W TL Indicates the Z-direction velocity output by the DRYDEN model in the previous cycle. Indicates the DRYDEN model outputs the X-axis acceleration; Indicates that the DRYDEN model outputs the Y-axis acceleration. Indicates the Z-direction acceleration output by the DRYDEN model; S1879. Exit the turbulence rate calculation sub-process.

Images