Flight simulation from takeoff to yawing of eVTOL airplane with coaxial propellers by fluid-rigid body interaction

In this study, we adopt a coupled fluid-rigid body simulation using the moving computational domain method and multi-axis sliding mesh method for the takeoff, hovering, and yawing flight of an electric vertical takeoff and landing aircraft (eVTOL). The aircraft model has four pairs of coaxial propellers, and the computational domain is divided into three domains to move the aircraft and eight propeller domains to rotate the propellers. As a result, we clarify the behavior and aerodynamic force of the aircraft when the input values are determined by the automatic control. The results in the flow field also show that the downwash spreads in a crisscross pattern on the ground, the wind reaches different ranges on the ground depending on the flight altitude, and that the coaxial propeller causes an asymmetry in the velocity field during yawing. Consequently, we conclude that this method is effective for the flight simulation of an eVTOL.


Introduction
Recently, the urban population has been increasing, and traffic congestion caused by private cars has become an important issue. To deal with the problem, the development of urban air mobility (UAM) has attracted attention around the world, and air taxi services using existing aircraft such as helicopters and light jets have become more popular [1]. In addition, the next generation of UAM, the electric vertical takeoff and landing vehicle (eVTOL), has gained attention as a means of transportation that is environmentally friendly and suitable for personal air travel. However, a large number of tests are needed to confirm the safety of flight, and thus the time and cost of the development increases. Aerodynamic performance is an important elemental technology of eVTOL, and various development methods have been studied. Experimental studies of eVTOL have been conducted to measure the lift and drag forces through wind tunnel tests [2] and to visualize the flow field around a rotating rotor using particle image velocimetry [3]. Numerical simulations have been adopted to study the aerodynamic interaction of coaxial propellers [4]. The simulations for an entire aircraft are conducted by modeling the propeller lift [5] or by computing the actual propeller rotation [6]. However, these are numerical wind tunnel contexts where the body is fixed, and they only focus on the steady state. This means that most of the previous studies have focused on only a few aspects of eVTOL aerodynamic performance. To reduce time and cost in the development process, reproducing flight tests by computer simulation and reducing the number of tests with an actual aircraft are required. For this purpose, development of a simulation method that can achieve flight-test scenario by coupling multiple factors is necessary. In particular, the unsteady flow generated by the flight and the aerodynamic effects exerted on the aircraft are currently evaluated by actual flight tests, and these are just a number of the issues to be addressed through eVTOL virtual flight tests. A simulation method has been developed that combines fluid and rigid-body dynamics using a moving grid [7], enabling various kinds of aircraft behavior to be simulated, such as takeoff, landing, cruising and crashing, as well as the transition between these flight phases. Thus, the number of prototypes and experiments can be reduced, which leads to various benefits including the reduction of cost and time in the development process. The objective of this study is to computationally examine the aerodynamic effect and its behavior in response to the effects of rotating the propellers and moving the body via the moving grid method. Since such examinations for eVTOL have not been carried out, basic flight behavior must be simulated first. Incidentally, the impact of downwash on objects on the ground during takeoff is regarded as a problem in the context of many kinds of aircraft [8,9]. In addition, the yawing method, which generates counter-torque by the speed difference of each propeller, is unique to multi-rotor aircraft, and simulating this behavior is effective in understanding the flow field around such aircraft. Furthermore, yawing flight is a suitable topic to confirm whether it is possible to simulate not only the translational but also the rotational motion of the aircraft. Therefore, in this study, we simulate takeoff and horizontal rotation at a constant altitude. An aircraft model with coaxial propellers [10] is used for the simulation. There are various types of eVTOLs, such as the tilt-rotor type [11] and one with propellers placed on the top of the main body [12]. The type used in this study is one of the options for practical use because of the compactness of the aircraft and the safety to keep the aircraft maneuverable in the event of propeller trouble, which is attributed to propeller redundancy [13]. To achieve free movement of the grid in three-dimensional (3D) space, we adopt the moving computational domain (MCD) method [14], which is based on the moving-grid finite volume (MGFV) method [15]. Since the MGFV method discretizes the governing equations in a space-time unified four-dimensional (4D) control volume, it automatically satisfies the geometric conservation law. We also introduce the multi-axis sliding mesh method [7], which divides computational domains using arbitrary boundary planes and moves a part of them. By combining the two methods, it is possible to compute the rotation of the propeller and the movement of the aircraft directly. Furthermore, by incorporating proportional-differential (PD) control into the computational program, the propeller speed can be automatically adjusted by providing only the target aircraft status (such as altitude, speed, and attitude). This computationally reproduces the flight as if the pilot controls the aircraft. This paper computationally demonstrates that flight simulation techniques considering fluid-rigid body interactions are capable of simulating an eVTOL's flight.

Governing equations
The governing equations are composed of the 3D Euler equations written in the conservation form: where q denotes the vector of conserved variables, and E, F , G represent the inviscid flux vectors. ρ is the density, u, v, w are the x, y, z components of the velocity vector, and e is the total energy per unit volume. Applying the perfect gas law, the pressure p is defined by where the ratio of specific heats γ is 1.4 in this study.

Flow simulation around moving objects
The conventional approach of placing objects in a uniform flow makes it difficult to simulate the unsteady flow generated around moving objects. Therefore, we adopt the MGFV method, which is one of the moving grid methods. This method satisfies the geometric conservation law because it uses the 4D control volume in a space-time unified domain, and thereby creating a situation where the flow is excited by the object's motion. The conserved quantities are defined in the cell center, and the flux vectors are evaluated by Roe's flux difference splitting method [16] with a monotonic upstream-centered scheme for conservation laws (MUSCL) scheme. The variables are reconstructed using the Green-Gauss method and Hishida's van Leer-like limiter [17]. The dual time stepping is adopted to solve the governing equations, where the two-stage Runge-Kutta method is used for pseudotime marching in internal iterations. Furthermore, we introduce the MCD method and multi-axis sliding mesh method to eliminate the restriction of the size of the computational domain and to enable free movement of the object. In the MCD method, the computational domain moves in line with the objects, as shown in Fig. 1. The multi-axis sliding mesh method moves only a part of the grid rotationally by dividing the computational domains, and the physical quantities are transferred through the boundaries of the divided domain. By combining the MCD and multi-axis sliding mesh methods, it is possible to simulate the flow generated by the relative motion of objects such as the propeller rotation and takeoff.

Combined simulation of fluid and rigid body
In this paper, the aerodynamic forces on the aircraft determine its position and orientation. The fluid-rigid body interaction proceeds in the following procedure at each time step: where m represents the mass, I is the inertia tensor, and r , F , ω , and T indicate the position, external force, angular velocity, and torque vectors, respectively.

Computational model
The aircraft model is based on the SD-03, which was developed by SkyDrive in Japan [10]. The SD-03 specifications are shown in Table 1. The shape is shown in Fig. 2(a), and the 3D model is illustrated in Fig. 2(b). The inertia tensor to solve the motion equations is calculated by assuming the density of the model is constant. The computational grids are generated by an unstructured grid generator MEGG3D [18,19]. The number of cells is approximately 3.2 million, and the minimum cell width at the front edge of the propeller is represented as 0.001L with the characteristic length L, which is 4 m . The computational domain is divided into 11 domains: one outer domain including the external and ground boundaries, one middle domain to change the altitude, one body domain to express the attitude change, and eight propeller domains to rotate the propellers. Figure 3 displays the divided domains and how they  Table 2 lists the number of cells and the radius of each spherical or cylindrical domain. The computational grid for each domain is shown in Fig. 4.

Computational conditions
The initial conditions are shown in Table 3. These variables are nondimensionalized by the characteristic speed and density U = 340 m/s and ρ = 1.25 kg/m 3 . The boundary conditions are the slip wall conditions on the ground and at the aircraft surface and the Riemann invariant boundary conditions at the outer boundary.

Flight conditions
The simulation is carried out from takeoff to climbing to the target altitude, and then the aircraft performs a horizontal rotation at a constant altitude, as shown in Fig. 5. The operations applied are as follows:       Two cases are demonstrated in this paper: the target altitude is fixed to 12 m and is changed from 12 m to 18 m while the aircraft is in the process of climbing.

Proportional-differential control of aircraft
In this model, propellers 2, 3, 5, and 8 in Fig. 2(b) rotate clockwise, and propellers 1, 4, 6, and 7 rotate counterclockwise. The rotation of the propellers is determined by the superposition of four operation variables: the throttle to control the altitude, the aileron to control the roll, the rudder to control the yaw, and the elevator to control the pitch. The operations are shown in Table 4. These operating values are determined by applying the PD controller Eq. 7 and 8 to the climb speed v c [km/h] and the angular velocities around the x, y, and z-axes ω roll , ω yaw , and ω pitch [rad/s]: (6) v target = C a (y target − y current ). where h(t), r target (t) , and r current (t) represent the handling, target state, and the current state values, respectively. The proportional and differential gain, K p and K d , for each operation (Table 5) are determined by tuning the linear single-input single-output (SISO) transfer function model while checking their step responses. Here, the kinematic model of the aircraft is linearized by assuming that the changes in the propeller speed are small in the equations of motion with six degrees of freedom. By applying the mixing process [20] to this model, it can be considered as a SISO system without coupling for v c , ω roll , ω yaw , and ω pitch . The transfer functions in the Laplace domain from the input propeller rotation h to v c , ω roll , ω yaw , and ω pitch are regarded as follows: where D is the distance from the center of gravity of the aircraft to the propeller, and the inertia tensor I is calculated by the 3DCAD. In addition, K T , K A are constant values expressed as where g is the gravitational acceleration. A, B are constants found in Eqs. 15     where ω rot represents the propeller speed, F prop , T prop are the propeller thrust and torque, respectively. All physical quantities used in this determination are nondimensionalized by the characteristic length L, velocity U, and density ρ. Figure 6 shows the comparison of the responses between the transfer function control model and the flight simulation for the aircraft climbing and yaw angular velocities. Note that the yaw angular velocity has a limit on the input propeller speed, so the response when the limit is removed is also plotted. Furthermore, since the control model does not consider the propeller acceleration time nor the delay until the flow field is generated and aerodynamic forces are applied to the aircraft, the plots are overlaid to eliminate the time lag until the start of the motion. This figure shows that the aircraft reached the target values without overshooting in the flight simulation as designed by the control model. However, the time to reach the target value is longer in the flight simulation, and this difference is thought to be caused by not considering aerodynamic lag and drag in the control model. In addition, the motion model is linearized so that the change in propeller speed is assumed to be small, which could affect the curve of convergence near the target values. Elevator r 2 × 10 3 8 × 10 2 Fig. 6 Comparison of responses between flight simulation and control model   Figure 7 shows the aircraft altitude and yaw angular velocities. In both altitude targets, both altitude and angular velocities increase linearly, reach the target without exceeding it, and keep those states. Figure 8 shows the input operating values. Figure 9 shows the climbing force and horizontal rotational torque exerted on the aircraft. These indicate that the operating values fluctuate sharply when the target values change. While the throttle input during climbing is about 1.04 times that of hovering except for a momentarily large amount of input, the rudder input during the increase in angular velocity is about 3.73 times greater than that of constant rotation. Similar patterns can be seen in the force and torque graphs.  Figure 10 shows the velocity isosurface at 10 m/s where the altitude is 12 m . The flows generated by the propellers interact with each other and spread out in a crisscross pattern. When the airflow is spread the widest, the 10 m/s wind reaches the area about 16.9 m away from the aircraft. Figure 11 shows the velocity isosurface at 10 m/s during its yawing. These results indicate that during the yawing, an asymmetrical flow field is generated along the rotation. When the flow on the ground is nearly steady, the wind over 10 m/s spreads about 9.1 m and 8.0 m from the center of the aircraft in the cases when the target altitude is 12 and 18 m , respectively. It is clear that the impact of downwash on the ground becomes larger at lower altitudes. Figure 12 shows the velocity distribution in the plane across the propeller on the right and rear sides when the target altitude is 12 m . Here, the flight state in (a) is just after takeoff, (b) is during climbing, and (c) is during yawing. The flow velocity around the propellers and under the body is roughly related to propeller rotation, and a high-speed flow is observed only with large inputs. Figure 12(c) also shows that there is a difference in velocity distribution depending on whether the upper or lower propeller's rotation is increased.

Flow field around aircraft
In addition, Fig. 13 shows the flow vectors around the downwash during takeoff. In plane a, the flow below the propeller is rotating along with the rotation of the upper propellers. In plane b, the effect of the rotation is decreased and the flow is directed toward the center, but in plane c near the ground surface, the direction of the flow changes to spread outward. This phenomenon can be explained by the effects of the large vortex near the ground, as shown in Fig. 14.

Conclusion
This study presented a flight simulation of an eVTOL aircraft with coaxial propellers. The MGFV method was introduced to enable aircraft motion, and the MCD method enabled the aircraft to move freely in a 3D domain. The rotation of the propellers was simulated without any simplification by applying the multi-axis sliding mesh approach. The flight conditions of the aircraft were a takeoff toward a target altitude and yawing that keep the altitude. As a result, we can visualize the flow field generated in a series from the behavior. Although the classical PD controller with constant target values provides the appropriate operations to achieve the target behavior, the difference from the linear simplified dynamic model is apparent, which confirms the need to simulate the motion considering the effects of the fluid dynamics. These indicate that the combined flight simulation of fluid and rigid body with the MCD and multi-axis sliding mesh methods is an effective approach to simulate the flight-test scenario in a computer, to examine the aerodynamics, and to visualize flow fields of the eVTOL aircraft. Furthermore, the results show that the difference in altitude (12 and 18 m) did not affect the aircraft's behavior significantly, but made a difference in the wind strength on the ground although the intensity is not as great as the 1.5 times altitude difference. Future work will introduce atmospheric conditions and examine the behavior of eVTOLs at various flight attitudes.