In this research, a quadrotor eVTOL aircraft concept proposed by Silva et al. (2018), as shown in Figure 1A, is used to study the effect of wind on a multirotor eVTOL aircraft trajectory. This eVTOL has six-seater (up to 545 kg payload) capacity with a long range cruise airspeed (LRC) of 50.4 m/s (98 kts).

Table 1 lists the aircraft performance data for the lateral trajectory optimization in cruise phase.

To study the effect of wind on the cruise phase of the multirotor eVTOL aircraft, a lateral flight dynamics model (two dimensional in space and one dimensional in time) is considered. The four lateral states of the model are: λ, τ, V, ψ; where λ is the latitude, τ is the longitude, V is the true airspeed (assumed to have only horizontal component during the cruise) and ψ is the heading angle w. r.t north (Tsuchiya et al., 2009; Yomchinda et al., 2011). The three control variables of the model are: the net thrust (T), the rotor tip-path-plane pitch angle (θ) and the rotor tip-path-plane roll (bank) angle (ϕ). Therefore, the quasi-steady cruise flight dynamics and kinematics of the multirotor eVTOL aircraft considering wind in a vehicle-carried frame of reference are as follows (Erzberger and Lee, 1980; Tsuchiya et al., 2009; Yomchinda et al., 2011; Sridhar et al., 2015; Weitz, 2015; Pradeep, 2019; Pradeep and Wei, 2019): d V d t = T ⁡ cos ⁡ ϕ ⁡ sin ⁡ θ − D m − d W N d t cos ⁡ ψ − d W E d t sin ⁡ ψ (1) V d ψ d t = T ⁡ sin ⁡ ϕ m + d W N d t sin ⁡ ψ − d W E d t cos ⁡ ψ (2) R Earth + h d λ d t = V ⁡ cos ⁡ ψ + W N = V G S ⁡ cos ⁡ χ (3) R Earth + h cos ⁡ λ d τ d t = V ⁡ sin ⁡ ψ + W E = V G S ⁡ sin ⁡ χ (4) where D is the parasite drag, V _GS is the ground speed, χ is the course, h is the altitude above mean sea level and R _Earth is the mean radius of the Earth, and W _N and W _E are the components of the wind in north and east directions, respectively. The time derivative of wind components is assumed to be zero given the short duration flights in the UAM environment (Thipphavong et al., 2018; Uber-Elevate, 2019).

For the quadrotor eVTOL aircraft, assuming that the rotors have negligible interference with each other, the net thrust (T) produced by the four rotors is given by: T = ∑ n = 1 4 T rotor n (5) where ( T rotor ) n is the thrust produced by the nth rotor. Also, assuming all the rotors produce the same amount of thrust (T _rotor) in cruise phase, the net thrust (T) produced by the rotors is given by: T = 4 T rotor (6)

The parasite drag (D) on the multirotor eVTOL is calculated as follows (Silva et al., 2018): D = 1.1984 ρ V 2 2 (7) where ρ is the density of air, which is a function of altitude and 1.1984 is a product of drag coefficient and reference area (Silva et al., 2018).

Using momentum theory (Heyson, 1975; Hoffmann et al., 2007; Johnson, 2012), the hover induced velocity (v _h) is given by: v h = T rotor 2 ρ A rotor (8) where A _rotor is the rotor diVsk area (πR ²) and R is the radius of the rotor.

Consider an isolated rotor in forward motion at true airspeed (V), with angle of attack (α) between the air-stream and the rotor disk (tip-path-plane). Using momentum theory in forward flight, the solution for induced velocity (v _i) is given by (Heyson, 1975; Hoffmann et al., 2007; Johnson, 2012): v i = v h 2 V ⁡ cos ⁡ α 2 + V ⁡ sin ⁡ α + v i 2 (9) where the hover induced velocity (v _h) on the right-hand side of Eq. 9 is computed using Equation 8.

Equation 9 is a quartic polynomial that can be analytically solved for induced velocity (v _i). Out of the four roots of the quartic Eq. 9, the one with real-positive and value lower than v _h is the correct solution for cruise phase. Equation 9 can also be numerically solved using an iterative technique with initial guess for v _i as v _h.

Once v _i is computed, the induced power loss of an isolated rotor (P _{induced rotor}) in forward flight is computed as follows (Johnson, 2012; Johnson, 2015): P induced rotor = κ T rotor v i (10) where the induced power factor (κ) is assumed to be 1.75 (Silva et al., 2018).

Based on the quasi-steady flight assumption, the instantaneous power required in forward cruise flight at a constant altitude is equal to the sum of the induced power, parasite power, and profile power as follows (Heyson, 1975; Leishman, 2002; Johnson, 2012; Pradeep, 2019; Pradeep and Wei, 2019): P required = P induced + P parasite + P profile (11)

The total induced power loss of the eVTOL aircraft (P _induced) is equal to the summation of the induced power loss of each rotor (P _{induced rotor}). Therefore, the induced power loss of the aircraft is given by Johnson (2012); Silva et al. (2018): P induced = ∑ n = 1 4 P induced rotor n = κ ∑ n = 1 4 T rotor n v i n (12)

The power required to propel the aircraft forward (the parasite power loss) at a constant altitude is given by Johnson (2012): P parasite = T V ⁡ sin ⁡ α (13)

The profile power loss is calculated from a mean blade drag coefficient (C _{d mean}) as follows (Leishman, 2002; Johnson, 2012; Silva et al., 2018): P profile = ρ A rotor Ω R 3 σ C d m e a n F P 8 (14) where Ω is the rotational velocity of the rotor blades, σ is the thrust weighted solidity ratio, C _{d mean} is the mean blade drag coefficient and F _P is the function that accounts for the increase of the blade section velocity with rotor edgewise and axial speed (Johnson, 2012; Johnson et al., 2018; Silva et al., 2018). However, F _P is assumed to be a constant (see Table 1) in this study because of cruise at constant altitude and nominal cruise speed.

Therefore, using the equations from (11) to (14), the instantaneous power required (P _required) in forward flight is given by: P required = κ ∑ n = 1 4 T rotor n v i n + T V ⁡ sin ⁡ α + ρ A rotor Ω R 3 σ C d m e a n F P 8 (15)

The performance index for the wind-optimal lateral trajectory optimization problem is constructed as follows: J = ∫ t 0 t f P required d t = ∫ t 0 t f κ ∑ n = 1 4 T rotor n v i n + T V ⁡ sin ⁡ α + ρ A rotor Ω R 3 σ C d m e a n F P 8 d t (16) where t ₀ is the initial flight time at the top of climb (TOC) placed directly above the origin vertiport at the cruise altitude and t _f is the final flight time to reach the top of descent (TOD) placed directly above the destination vertiport at the cruise altitude.

For a level flight (cruise) in the absence of vertical component of wind, the net vertical force on the multirotor eVTOL aircraft is zero; therefore, the following path constraint is imposed on the problem: T ⁡ cos ⁡ ϕ ⁡ cos ⁡ θ = m g (17) where m is the mass of the aircraft and g is the acceleration due to gravity.

The instantaneous power required (P_required in kw) is bounded by the total deliverable power (P _max) of the quadrotor eVTOL aircraft (Johnson, 2015; Silva et al., 2018): P required ≤ 494.25 (18) where P _required is defined in Equation 15.

The course angle (χ) for the great-circle trajectory between the two waypoints is calculated as follows (Chatterji et al., 1996): tan ⁡ χ = sin τ 2 − τ 1 cos λ 2 sin λ 2 ⁡ cos λ 1 − sin λ 1 ⁡ cos λ 2 ⁡ cos τ 2 − τ 1 (19) where the origin latitude-longitude is (λ ₁, τ ₁) and the destination latitude-longitude is (λ ₂, τ ₂). However, to generate a great-circle trajectory between the two waypoints using the optimal control framework developed in this research, a path constraint as a function of the wind components (W _N and W _E ), latitude-longitude coordinates of the two waypoints and heading angle (ψ) is required. Therefore, the course angle (χ) needs to be eliminated from Equation 19.

Using Eqs 3, 4, the heading angle (ψ) required to fly the course angle (χ) in the presence of the wind is computed as follows: tan ⁡ χ = V s i n ψ + W E V c o s ψ + W N (20) where W _N and W _E are the components of the wind in north and east directions, respectively.

Equations 19, 20 result in the following relation: V s i n ψ + W E V c o s ψ + W N = sin τ 2 − τ 1 cos λ 2 sin λ 2 ⁡ cos λ 1 − sin λ 1 ⁡ cos λ 2 ⁡ cos τ 2 − τ 1 (21)

Hence, the path constraint for the great-circle trajectory between the two waypoints is given by: V s i n ψ + W E sin λ 2 ⁡ cos λ 1 − sin λ 1 ⁡ cos λ 2 ⁡ cos τ 2 − τ 1 − V c o s ψ + W N sin τ 2 − τ 1 cos λ 2 = 0 (22)

The trajectory optimization problems can be numerically solved, either using the direct or indirect methods (Betts and Huffman, 1998). Direct methods typically discretize the trajectory optimization problem, and convert the original trajectory optimization problem into a nonlinear program. On the other hand, indirect methods are characterized by explicitly solving the optimality conditions stated in terms of the adjoint differential equations, the maximum principle, and associated boundary (transversality) conditions (Bryson, 1975). Using the calculus of variations, the optimal control necessary conditions can be derived by setting the first variation of the Hamiltonian function to zero. A common way to distinguish these two methods is that a direct method discretizes and then optimizes, while an indirect method optimizes and then discretizes (Betts and Huffman, 1998; Kelly, 2017).

In this research, PSOPT has been used to solve the optimal control problem to generate wind-optimal lateral trajectories for the multirotor eVTOL aircraft. PSOPT is an open-source optimal control software package written in C++ that uses direct collocation methods such as pseudospectral methods (Becerra, 2010). Pseudospectral methods directly discretize the original optimal control problem to formulate a nonlinear programming problem, which is then solved numerically using a sparse nonlinear programming solver to find approximate local optimal solutions. IPOPT is an open-source C++ package for large-scale nonlinear optimization, which uses an interior point method (Wächter and Biegler, 2006; Becerra, 2010). IPOPT is the default nonlinear programming algorithm used by PSOPT. Approximation theory and practice show that pseudospectral methods are well suited for approximating smooth functions, integration, and differentiation, i.e., all of which are relevant to optimal control problems (Becerra, 2010). Nearly all trajectory optimization techniques require a good initial guess to begin the optimization. In the best case, a good initialization ensures that the solver rapidly arrives at the optimal solution (Kelly, 2017).

The Rapid Refresh (RR) operational weather prediction system, available on an hourly basis at 13 km spatial resolution from the National Center for Environmental Prediction (Benjamin et al., 1998), has been used to extract the wind data for analysis. The wind data are extracted for the entire month of January 2019 (i.e., 31 days) during peak traffic hours in the morning and evening, i.e., i) 7 a.m.–11 a.m. (local time), and ii) 3 p.m.–7 p.m. (local time) at cruise altitude (1,600 ft above MSL) in 1-min intervals at grid points equispaced 10 km apart. The grid points cover 10⁴ km² area at the Dallas-Fort Worth and New York metropolitan areas of the United States.

NASA’s Future Air Traffic Management (ATM) Concepts Evaluation Tool (FACET), designed and developed over the last 20 years, provides a flexible simulation environment for the exploration, development and evaluation of advanced ATM concepts (Bilimoria et al., 2001). FACET models four-dimensional (4D) aircraft trajectories in the presence of winds using round-earth kinematic equations. Aircraft can be flown along flight plan routes or great-circle routes as they climb, cruise and descend per their aircraft-type performance models. Performance parameters of the aircraft are obtained from a data look-up table. Therefore, a look-up table is created for the multirotor eVTOL aircraft based on the performance data of the conceptual multirotor eVTOL proposed by Silva et al. (2018). The flight duration results of the optimal control framework using PSOPT solver are validated using the neighboring-optimal wind routing algorithm (Jardin and Bryson, 2001) implemented in FACET. FACET implementation performs a bi-linear interpolation in spatial and temporal dimensions on the available gridded wind data.

The following two metrics are used to compute operational benefits of flying wind-optimal lateral trajectory compared to great-circle trajectory for a given origin and destination.

Similar to the Dallas-Fort Worth case studies, two types of wind scenarios are considered: 1) strongest wind (as shown in Figure 6) and ii) highest spatial variability of the wind (as shown in Figure 8) in the New York metropolitan area. The vertiports and direct route lengths of 30 and 50 nautical miles are based on the Air Traffic Management - eXploration (ATM-X), Urban Air Mobility (UAM) Experiment two conducted at NASA Ames (Verma et al., 2020).