Aerodynamic analysis and modeling of coaxial ducted fan aircraft with the ceiling effect

In this study, a novel coaxial ducted aircraft structure is considered, based on which the flow field with the ceiling effect under different vertical velocities and heights is obtained by exploiting the CFD (Computational Fluid Dynamics) numerical calculation methods. Taking into account the ceiling effect, a total lift mechanism model of the coaxial ducted aircraft structure is established based on the momentum theory and the blade element theory. In the proposed mechanism model, the relationship between the vertical climbing rate and the induced speed of the rotors and the induced velocity interaction between the coaxial rotors and the effect of the wake of the upper rotor on the lower rotor are considered. To verify the effectiveness and feasibility of the mechanism model, the numerical calculation and the bench test are conducted. The results indicated that, under the hovering condition, the error between the mechanism model and the bench test results falls into 15%; under the vertical climbing condition, the error between the mechanism model and the CFD numerical calculation results falls into 10%. It is concluded that the proposed mechanism model exhibits high adaptability and accuracy.


Introduction
UAVs (Unmanned Aerial Vehicles) have been extensively applied as an effective platform for various operations because of their strong maneuverability and high flexibility. For instance, UAVs have been used in numerous fields (e.g. aerial works, the communication relay, and the disaster monitoring). With the gradual maturity of the aircraft hardware development and the gradual improvement of relevant control algorithms, the overall practicability of UAVs is being improved day by day, so people hold increasingly higher requirements for its performance and the more stringent requirements for its environment (Ai et al., 2020b). The use scene of UAV has been progressively updated from a simple open environment to a complex environment, thereby significantly increasing the diversity of the UAV operation (Fumagalli et al., 2014;Ruggiero et al., 2018). However, when UAVs are in a complex environment, as restricted by the complex and changeable boundaries, UAVs will also encounter some challenges that are not available in the open environment.
Besides considering the collision between UAVs and obstacles, the biggest potential danger in a complex environment refers to the proximity effect: the ground effect, CONTACT Wei Fan fanweixx@bit.edu.cn the wall effect ( Robinson et al., 2014;Lee et al., 2012;Tanabe et al., 2018;Tavora, 2017), and the ceiling effect. When a UAV is near the environment for relevant operations, the proximity effect will impact the aerodynamic characteristics of the UAV, and even cause the UAV system to be unstable if the interaction is strong, which seriously reduces the safety and operation performance of the UAV (He et al., 2019). The main research on the proximity effect has currently been the ground effect, which is broadly used in UAV take-off and landing from the ground. There are various researches on lift empirical model of the ground effect and the related stable take-off and landing control (Cheeseman and Bennett, 1955;Davis & Pounds, 2016;Nonaka & Sugizaki, 2011;Powers et al., 2013;Li et al., 2015;. Compared with the ground effect, the ceiling effect has been less studied. When the rotor rotates close to the constraint surface (e.g. the bottom of the bridge or the interior ceiling), the rotor will undergo an aerodynamic interaction with the ceiling, which is termed as the 'ceiling effect'. For some detection applications (e.g. bridge bottom or the interior ceiling detection), manual detection is often timeconsuming and labour-consuming, and its efficiency is not high. Thus, aircraft can be used instead of the manual to complete the task, so it is of a significance for practical researches. The ceiling effect can make the rotor produce more lift at a constant power (Rossow, 1985;Wang and Huang, 2018), so the energy of the aircraft can be saved, and the endurance of the whole aircraft can be improved to a certain extent. However, under the additional lift of the ceiling effect, the prototype will obtain an upward acceleration, which causes the collision between the prototype and the ceiling. As a result, the aircraft will be unstable, and the prototype will be even caused to crash . Compared with the ground effect, the ceiling effect is more dangerous to some extent, so it should be studied in-depth, and its advantages of energy-saving should be exploited to avoid the possibility of collision between the prototype and the ceiling (Kocer, Kumtepeli et al., 2019).
Over the past few years, experts and scholars have conducted related researches on the ceiling effect. Cheeseman Bennett proposed a method to normalize the ratio of the rotor to the ceiling and the rotor radius; the ratio was taken as the ceiling effect coefficient and introduced into the total lift of the system (Cheeseman and Bennett, 1955). However, a large deviation on the lift prediction of a small aircraft with the ceiling effect by the follow-up research is found. Based on the uniform inflow control theory,  further refined the ceiling effect coefficient to further reduce the model error. According to Robinson et al. (2016), the numerical study was conducted on the time-varying aerodynamic effect of the system flying close to the ceiling. Jimenez-Cano et al. (2019) studied the thrust model fitted by CFD calculation when the distance between the prototype and the ceiling is different. However, the lift model was only related to the distance between the rotor and the ceiling, while the effect of an unsteady state was not considered. According to the research of Conyers, the ceiling effect of the single wing and multi-rotor UAVs was assessed by performing comprehensive experiments. The ceiling lift data of the prototype hovering under different parameter combinations were given and then compared with the theoretical value (Conyers et al., 2018). Besides the research on the empirical modeling under the ceiling effect, the study on the control of UAVs with the ceiling effect is also gradually advancing. Kocer, Tiryaki et al. (2019) proposed a novel framework of the prototype in the ceiling effect, which included the force estimation based on the constrained optimization and the nonlinear controller based on optimization. Thus, the prototype can be stable close to the ceiling. In the standard control, the change of rotor rotation speed is assumed to be proportional to the generated thrust. To compensate for the ceiling effect, Powers et al. (2013) adopt the real-time data to solve the rotor speed thrust equation instead of the standard control allocation. Kocer et al. (2018) formulated a control strategy combining NMPC and double loop PID to improve the Stability performance of the prototype with the ceiling effect. Gao et al. (2019) exploiting Ground and Ceiling Effects on MAV to plan a trajectory while minimizing energy consumption.
At present, the main research object of the ceiling effect refers to the multi-rotor aircraft with the single blade configuration. To the best of the author's knowledge, no research has been conducted on the ceiling effect of coaxial ducted fan aircraft.
However, the UAV with an open rotor structure (e.g. a helicopter or multi-rotor, which is extensively used at present) can't contact the target in a close range for its large safe flight space, which significantly limits its application scenarios and operation scopes. Some UAVs with special structures are capable of near the ceiling, and they fail to carry relevant tools to conduct relevant operations in the complex environment due to their sizes. Compared with the helicopter and the multi-rotor, the ducted vehicle exhibits a compact structure that can directly contact the external environment shortly, and it has fewer requirements for the flight space. As impacted by the existence of ducts, the rotor can avoid the direct contact with the outside world, so it also exhibits a better safety performance. Besides, under the effect of the airflow at the lip of the duct, the ducted fan aircraft will generate a greater lift force than the open rotor aircraft under the identical rotor size (Pflimlin et al., 2010). Figure 1 presents a comparison diagram of the space required for the safe flight of coaxial ducted aircraft, quadrotors, and helicopters. According to Figure 1, due to the circular expansion effect of duct (Johnson & Turbe, 2006), the safe flight space required for ducted fan aircraft is significantly less than that of a quadrotor and a helicopter, and it can contact targets or obstacles in a close range. Moreover, it exhibits higher trafficability and greater maneuverability under the identical payload. Given the mentioned advantages, the ducted fan aircraft is more suitable for related operations in a complex environment than the open rotor aircraft.
Existing studies largely focused on the lift model of the aircraft under the ceiling effect at the constant power or speed. Most of the mentioned empirical models only fit the total lift of the system to the expression of the normalized ratio of the distance from the rotor to the ceiling and the radius of the rotor. The model fails to consider the vertical velocity of the prototype, nor does it reveal the change of the system lift under the ceiling effect from the perspective of the mechanism modeling. Moreover, the multi-rotor lift model is not suitable for the novel coaxial duct structure aircraft. The main contributions of this paper include two parts. One is to study the flow field characteristics of a novel ducted fan structure aircraft with different vertical climbing rates with the ceiling effect. Another, instead of using the existing empirical model that only considered that the system lift variation is related to the distance from the UAV to the ceiling, more other factors are considered. The relationship between the vertical velocity and the induced velocity of the rotor is deduced based on the momentum theory, the blade element theory, and consider the interaction between the upper and lower rotor and the wake area of the upper to the lower rotor. Lastly, the lift mechanism model of the prototype with the climbing rates under the ceiling effect is built. It is worth noting that no lift mechanism model has been available for coaxial ducted aircrafts with the ceiling effect.
The structure of the rest of this study is as follows. In Section 2, the aerodynamic characteristics and flow field of the prototype with the ceiling at different positions and different velocities are analyzed. In section 3, the lift of the prototype considering the ceiling effect is modeled. In section 4, the accuracy of the lift mechanism model is verified by using the bench test and the CFD numerical calculation. In section 5, the main conclusions and prospects for the future of the study are summarized.

Aerodynamic analysis of the prototype with the ceiling effect
According to the existing team research, a series of related developments were achieved with the ducted structure aircraft as the objects (Ai et al., 2020b;Han et al., 2019;Zhang et al., 2016). Figure 2 illustrates the structure of the latest prototype.
The prototype mainly comprises a lift system, a control surface system, a support frame, and an airborne control system. The lift system of the prototype uses the composite structure of duct and coaxial twin rotors to ensure high dynamic performance and applicability of the prototype. On the premise that the prototype has the payload carrying relevant operation tools and various sensors, the rotor diameter and the structural size of the fuselage are reduced to make the structure of the prototype more compact. Under different requirements, the corresponding operating tools can be installed under the control system as well.
Moreover, as impacted by the characteristics of the duct structure, the prototype is usually applied for the low-velocity flight and hovering operation when it is close to objects in a complex environment. When the prototype is close to the ceiling, the ceiling effect will be generated. The ceiling effect significantly impacts the hovering and low-velocity maneuvering of the prototype, which will hinder the stable hovering and the related operations when the prototype is near the ceiling for a certain distance. Accordingly, the flow field characteristics of the prototype under the typical working conditions should be analyzed when the ceiling effect occurs. Obviously, the flow field characteristics of the prototype with the ceiling effect can facilitate the subsequent lift mechanism modeling. In the present section, the flow field of the prototype is analyzed under two typical conditions of a small vertical climbing rate and hovering at different positions under the ceiling effect.
As impacted by the symmetrical structure of the prototype, the major effect of the ceiling effect on the system refers to its lift (Conyers et al., 2018;Gao et al., 2019;Kocer, Kumtepeli et al., 2019). As impacted by the large distance between the front and rear ducts, the mutual interference of the rotors in the front and rear ducts is ignored (Zhou et al., 2017). Moreover, in the hovering and vertical climbing processes, the rudder surface is kept vertical without any deflection, and the total lift of the prototype originates from the composite lift system. Therefore, in the followup research, a single lift system is adopted to replace the whole prototype for the aerodynamic analysis and modeling.

Analysis of hovering condition of the prototype with the ceiling effect
When the UAV's rotor rotates, it will drive the surrounding airflow to move together. Accordingly, when the UAV is near the ceiling, the inflow process of the rotor will drive the relative static airflow around the ceiling; as a result, the interaction is generated. When the prototype with the ceiling effect, it will affect the rotor free inflow process, which should be analyzed according to the flow field of the prototype at different positions. As a type of CFD software, Fluent can be adopted to analyze the flow field and the aerodynamic characteristics when the rotors rotate close to the constraint; it is extensively used in aerospace, heat transfer and many other fields (Chen et al., 2019;Duan et al., 2020;Ez Abadi et al., 2020;Farzaneh-Gord et al., 2019;Galván et al., 2011;Ghalandari et al., 2019;Halfi et al., 2020;Khalil et al., 2019;Kong et al., 2018;Liu et al., 2018;Niu et al., 2018;Paz et al., 2020;Zhou et al., 2020).
The geometric model of a single duct with a coaxial twin-rotor is processed. Figure 3 illustrates the calculation domain and boundary conditions of the single duct with a coaxial twin-rotor under the ceiling effect. The computational domain is a cube of 20R × 20R × 20R, where R denotes the rotor radius, d denotes the distance from the upper rotor to the ceiling, and x, y and z are cartesian coordinates in the computational domain. Overall, the ceiling, duct and rotors are all set as the smooth wall without sliding. In addition, it should be noted that in Figure 3, the domain around the upper and lower rotors is a rotational domain, and other domains are stationary domains. The stationary domain and the rotational domain need to be connected by the interface, and the system interpolation calculation between the two domains should be carried out. Therefore, it is necessary to ensure that the meshes on the interface are consistent with each other and reduce the interpolation error.
The grid setting of the wall boundary layer is critical to the flow field simulation. In order to simulate the change of fluid shear stress in the boundary layer region, it is necessary to set enough grids with uniform transition in the boundary layer, and the height of the first layer grid on the wall is important. The dimensionless parameter y + is used to indicate whether the height is appropriate.
In this study, the classical (SST) k-w turbulence model which is suitable for rotor close to the constrained rotation was selected for the simulation (Habibnia & Pascoa, 2019;Duan et al., 2020). The SST k-ω turbulence model uses the default wall processing method when close to the wall. The recommended y + value range should meet: ( 1 ) In this paper, the model with the distance of 1.25R from the upper rotor to the ceiling is taken as an example, and the subsequent aerodynamic analysis steps are similar under the different distances between the upper rotor and the ceiling. Unstructured grid is used in the flow field simulation of the coaxial twin-rotor ducted fan. Moreover, ICEM software is used to generate the mesh files. The quality of the generated grids exceeds 0.3, and the grid half profile is shown in Figure 4. The number of all grids in the two rotational domains is about 2.5 million, while the number of all grids in the stationary domain is about 4 million, and the total number of all grids is about 9 million.
The spacing of the grid has a great influence on the numerical results. When the increase of the number of grids has little influence on the system numerical results, it shows that the numerical results are reasonable, and the results are considered to be independent of the grid (Kim et al., 2021;Shaheen et al., 2015).
To examine the grid dependence in terms of results, another three levels of grids with different grid numbers are generated by adjusting the scale factor. The grid numbers are 13 million, 10 million and 6.2 million, respectively. Moreover, the results of fine grid simulation are taken as the standard. Table 1 shows the total numbers of cells, computational time, and percentage error of lift coefficient with other levels grid and fine level grid. The error in the lift coefficient of the prototype is found to be less than 3% for the results from the fine and  1.0 * 10 7 14h 2.25% Medium2 9.0 * 10 6 12h 2.36% Coarse 6.2 * 10 6 8h 6.84% medium grids. Therefore, to balance the computational accuracy and the computational time, the medium2 grid is selected. When the grid-dependency of the prototype is studied, the mesh can be imported into ANSYS-Fluent 19.0 software for simulation. Next, related operations are performed on cell zone conditions, boundary conditions, and mesh interface, and corresponding simulation parameters are set.
After the given simulation settings, the simulation iteration is initiated with the lift of the duct and the upper and lower rotor as the output objects. When the set residual is reached, the simulation ends. Furthermore, the.cas file and .data file can be imported into Tecplot for post-processing.

Analysis of flow field characteristics
When the prototype generates the ceiling effect, a mutual interference is developed among the duct, the rotors and the ceiling. In this section, the CFD simulation method is adopted to analyze the flow field mechanism.
As indicated from Figure 5, the initial inflow, which almost parallels the ceiling, follows the streamline direction. Then the airflow begins to deflect axially above the rotor and finally passes through the rotor. With the decrease in the ceiling height (from Figure 5 (a) to Figure 5 (d)), the radial velocity of the airflow close to the duct lip tends to increase, so a larger range of lowpressure area is formed. As a result, the total lift of the system increases. When the distance between the upper rotor and the ceiling decreases to a certain extent, the tip vortex of the rotors will be generated, which will affect the aerodynamic performance of the system.

Lift characteristic analysis
As impacted by the symmetrical structure of the prototype, when the prototype hovers and vertically climbs under the ceiling effect, according to the existing research, the most significant change of the system is the lift, and the change of the force and the torque in other directions is noticeably small, which can be ignored (Conyers et al., 2018;Gao et al., 2019;Kocer, Tiryaki et al., 2019). Thus, this section analyzes the lift characteristics of the prototype with the ceiling effect. Figure 6 shows the curve of lift characteristics of the prototype varying with the distance from the upper rotor to the ceiling at a constant velocity of 4000 rpm/min, where T ICE denotes the lift of the prototype with the ceiling effect, T OCE denotes the lift of the prototype out of the ceiling effect. The fitting model of the data in Figure 5 is 0.1898 * (z/R) −2.415 + 1.008, which satisfies the model form of a * (z/R) −b + c, indicating that the simulation data decay exponentially with distance is reasonable (Jardin et al., 2017).
According to Figure 6, as the distance from the upper rotor to the ceiling increases, the total lift of the system will decrease, which is the change can be approximately fitted as an exponential function. As compared with the total lift in an unconstrained environment, the total lift increases by nearly 135% under the distance between the upper rotor to the ceiling reaching d/R = 0.55. When the ceiling distance exceeds 3R, the ceiling effect does not impact the lift characteristics of the prototype.

Analysis of vertical climbing condition of the prototype with the ceiling effect
To compare with the analysis of the hovering condition, this study also analyzes the vertical uniform climbing condition of the prototype under the ceiling effect. Because of the existence of the ceiling, the vertical velocity of the prototype can't be excessively large.  Accordingly, under the distance between the prototype to the ceiling reaching 0.55R and 0.75R, the aerodynamic analysis is conducted at different velocities of 0.05 m/s, 0.15 m/s and 0.3 m/s. Moreover, the rotor speed is identical to that under the hovering conditions since the given condition is vertical climbing at a uniform rate. Consistent with the simulation setting of hovering conditions, the difference is that the inflow is velocity inlet, and the inflow velocity is set according to different velocity requirements (Mi, 2020).

Analysis of flow field characteristics
Complying with the flow field analysis of hovering at different distances from the prototype to the ceiling, the streamline diagrams of the flow field at different vertical velocities at the distances of 0.55R and 0.75R are set below.
As indicated from Figures 7 and 8, the initial airflow almost parallels the ceiling along the streamline direction until the airflow begins to deflect axially above the rotor and finally passes through the rotor, which is consistent with Figure 5. With the increase in the distance from the upper rotor to the ceiling, the deflected airflow is closer to the axial flow. Under the identical distance between the upper rotor and the ceiling, the range of the highvelocity zone of the radial velocity tends to decrease with the increase in the vertical climbing rate. Since the prototype has a vertical velocity, the velocity at the duct outlet increases as compared with that under the hovering condition. On the whole, the streamlined distribution of the whole system is consistent with hovering under the low velocity.
To more specifically analyze the effect of the vertical rising velocity on the prototype with the ceiling effect, Figures 9 and 10 also present the pressure cloud map of different vertical velocities under the distance of 0.55R and 0.75R between the upper rotor to the ceiling.
According to Figures 9 and 10, with the increase in the vertical climbing rate, the comprehensive low-pressure area of the upper and lower rotors and the duct lip tends to decrease. It is indicated that when the distance from the prototype to the ceiling is the same, the total lift of the system slightly decreases with the increase in the climbing rate.
As suggested from the comparison between Figure 7 and Figure 8 and between Figure 9 and Figure 10, when the vertical velocity of the prototype is same and the distance from the upper rotor to the ceiling is different, the range of the low-pressure region is altered significantly. Under the identical vertical climbing rate, the total lift of the system increases obviously with the decrease in the distance from the upper rotor to the ceiling, which complies with the change of the lift characteristics when the prototype hovers.
As concluded from the mentioned analysis, when the prototype climbs at a small rate, the biggest factor of the lift of the system continues to be the distance from the prototype to the ceiling, instead of the climbing rate. However, it is noteworthy that when the distance from the upper rotor to the ceiling is excessively low, it may cause the tip vortex and the flow separation, which should be considered in future studies. Figure 11 presents the lift characteristics of the prototype at different distances from the upper rotor to the ceiling at a constant speed of 4000 rpm/min. As indicated from Figure 11, under the identical distance from the upper rotor to the ceiling, with the increase in the climbing rate of the prototype, the total lift of the system decreases, which also complies with the trend of the low-pressure range in Figure 9 and 10. Moreover, with the increase in the velocity, the total lift of the system will decrease more at the identical distance interval. Accordingly, under the same rotor speed and the same distance between the rotor and the ceiling, the small vertical velocity slightly impacts

Lift characteristic analysis
the total lift of the system. However, this analysis is only applicable when the prototype climbs at a low rate. When the vertical climbing rate of the prototype is large, it may cause a stall phenomenon, which does not comply with this changing trend.

Derivation of lift mechanism model for vertical climbing rate of the prototype with the ceiling effect
Since the aerodynamic characteristics of the ducted rotor are sophisticated, to simplify the derivation of the lift mechanism model of the prototype with the ceiling effect, the following assumptions are made in this study: 1) The gas flowing into the duct ignores its viscosity; 2) The flow into the duct is steady, and the flow state does not change over time; 3) The gas flowing via the rotor is approximately axial flow; 4) The tip vortex, the root vortex and the flow separation loss are not considered; 5) The rotor is an infinitely thin circular plane, which produces a stable and uniform induced velocity.
For the convenience of modeling, the prototype falls into two parts, i.e. ducted and coaxial twin rotors. The axial flow of the coaxial twin-rotor prototype in this paper exhibits two characteristics: 1. The interaction of the induced velocity exists when the upper and lower rotors are rotating; 2. The lower rotor is partially in the wake flow of the upper rotor.

Shrinkage boundary of wake
When the open coaxial twin rotors and the ducted coaxial twin rotors are rotating, the wake of the upper rotor will shrink (Figure 12). However, due to the existence of the duct, the shrinkage of the wake boundary of the two configurations will be different, and the duct will cause the wake expansion.
When the prototype is climbing under the ceiling, its vertical velocity is usually small for safety reasons, so the effect of vertical velocity is ignored. In addition, compared with the prototype that hovers in different positions under the ceiling effect, the shrinkage of the wake boundary is found to be basically unchanged. Thus, when the lower rotor is in the wake flow of the upper rotor, the wake boundary R' of the rotor is related to the distance h between the upper and lower rotor according to the vortex theory. It is assumed that h' = h / R, which can be fitted as a polynomial: As indicated from Figure 13, when the distance between the upper and lower rotors exceeds 2R, the shrinkage of the wake of the upper rotor will be stable and remain unchanged. Compared with the open rotor, the shrink radius of the wake will be stabilized at 0.707r, and the existence of the duct will alter the shrink boundary of the wake.

Interaction between upper and lower rotors
When the upper and lower rotors of the prototype are rotating at a constant speed, their induced speed will be affected by each other. Likewise, for safety reasons, the vertical velocity is generally given smaller, and the effect of vertical velocity on mutual interference is ignored. Besides, compared with the prototype hovering in different positions under the ceiling effect, the interaction between the upper and lower blades is found to be basically unchanged. Accordingly, when the distance between the upper and lower rotors is set, the mutual interference between the upper and lower rotors can be determined.
The effect of induced velocity on the lower rotor is : The effect of induced velocity on the upper rotor is : Where a and b respectively denote the influence coefficients of the lower rotor to the upper rotor and the upper rotor to the lower rotor, which are associated with the distance between the upper and lower rotor and the radius ratio h'. a = 0.9945 + 1.5132h − 0.8823h 2 − 0.2645h 3 + 0.6008h 4 + 0.
To simplify the modeling process, the prototype is decomposed into the duct and coaxial twin rotors. At present, no clear modeling method has been available for the lift of the duct which will be obtained by the identification method. The blade element momentum theory can be used to calculate the rotor aerodynamic characteristics. Based on the conventional blade element momentum theory, the lift model of coaxial twin rotors in this study is remodeled by considering several factors (e.g. the ceiling effect, the interaction between upper and lower rotors, and the shrinkage of upper blade wake). Overall, the effect of the existence of the duct on the modeling of coaxial twin rotors largely consists of three aspects: 1) the existence of duct alters the velocity of the part of the lower rotor that is not in the wake area of the upper rotor; 2) the duct regulates the shrinking area of the upper rotor wake; 3) the existence of the duct diffuser alters the outflow area. Figure 14 illustrates the inflow model of the prototype with the coaxial twin rotors at a constant climbing rate under the ceiling effect, where V 0 denotes the vertical climbing rate of the prototype, V 1 is the upper surface velocity of the upper rotor, V 2 expresses the lower surface velocity of the upper rotor, V 3 is the velocity of the upper surface of the lower rotor inside the wake boundary, V 4 represents the velocity of the lower surface of the lower rotor inside the wake boundary, and P 1 , P 2 , P 3 and P 4 denote the pressures corresponding to different positions, z refers to the distance from the upper rotor to the ceiling, V ∞ expresses the outflow velocity at infinity of the duct diffuser, V 5 represents the velocity of the upper surface of the lower rotor outside the wake boundary, V 6 is the velocity of the lower surface of the lower rotor outside the wake boundary, P 0 denotes the standard atmospheric pressure, P (r) expresses the air pressure between the upper rotor plane and the ceiling, and V (r) is the radial velocity of flow between the upper rotor plane and the ceiling. Moreover, R is the radius of the rotor, R' is the radius of the area where the lower rotor is subjected to the wake of the upper rotor, and r represents the radial position of the center of the rotor outward.
The lower surface velocity of the upper rotor consists of the vertical climbing rate of the prototype, the induced velocity of the upper rotor, as well as the interference velocity of the lower rotor to the upper rotor. Thus, it is written as: Besides, the lower surface velocity of the lower rotor includes the induced velocity of the lower rotor, the interference velocity of the upper rotor to the lower rotor, as well as the wake velocity of the upper rotor.
The airflow through the upper rotor is assumed to be the approximately axial flow. Thus, the wake velocity V' of the upper rotor can be written as follows: V 3 is the same as V'. Likewise, V 4 can be written as: In addition, as impacted by the existence of the duct and ceiling, the velocity of the upper surface of the lower rotor outside the wake boundary should be V 5 instead of V 0 . As impacted by the significantly small gap between the duct and the upper rotor, it is approximately considered that the velocity from the tip of the upper rotor to the inner wall of the duct is equated with the velocity V 1 . Accordingly, V 5 is expressed as: V 5 = V 1 * [(R duct 2 -R 2 )/ (R 2 -R' 2 )], and the range of R is [R,R duct ]. Where R duct is the radius of the inner wall of the duct.
The momentum theory can be adopted to calculate the lift which is generated by the mass flow and induced velocity of the upper rotor. The mass flow via the upper rotor is expressed as:ṁ 1 = ρA 1 V 2 (9) Where A 1 denotes the area of the upper rotor. By adopting the free inflow to the outlet of the upper rotor as the control body, the total lift of the upper rotor can be determined by complying with the momentum theory.
To solve the total lift of the lower rotor by employing the momentum theory, it needs to be divided into two parts, i.e. the lift of the lower rotor in the wake area of the upper rotor, and the lift of the lower rotor not in the wake area of the upper rotor.
Moreover, the momentum theory can be used to determine the lift which is generated by the mass flow and induced velocity when the lower rotor is in the wake area of the upper rotor. The mass flow rate from the wake of the upper rotor to the lower rotor is expressed as: Where A 2 represents the area of the lower rotor in the wake area of the upper rotor. With the rotor inflow to the outlet of the duct as the control body, the lift of the lower rotor in the wake area of the upper rotor can be determined in accordance with the momentum theory.
In addition, the lift force of the lower rotor not in the wake area of the upper rotor can also be calculated: Where A 3 represents the area of the lower rotor which is not in the wake area of the upper rotor. As compared with the lift of the prototype in an unconstrained environment, the pressure between the upper rotor and the ceiling will change under the existence of the ceiling. It can also be seen from Figure 14 that the speed of the rotor deflects to a large extent. When the prototype with the ceiling effect, the upper rotor is directly affected by the ceiling effect, which will change the induced velocity of the upper rotor. The lower rotor is in the wake area of the upper rotor, which is not directly affected by the ceiling effect, whereas it is affected by the outflow of the upper rotor due to the ceiling effect, thereby altering the induced velocity. The flow rate of the air deflected from the ceiling to the upper rotor is equated with the flow rate of the airflow into the upper rotor, and the following formula is yielded: By complying with the Bernoulli equation of kinetic energy conservation, the relationship about the pressure and the velocity between the ceiling and the upper rotor can be listed: Accordingly, the lift increased T ICE by the upper rotor due to the ceiling effect can be determined.
The relationship between the lower surface of the lower rotor and the duct diffusion port can be listed by the conservation of mass flow.
Where A 4 represents the area of the duct diffuser. Therefore, the lift model of the upper and lower rotors is derived based on the classical momentum theory and by considering the ceiling effect. There are five unknowns: the lift T 1 and T 2 of the upper and lower rotors, the induced velocities V i1 and V i2 of the upper and lower rotors, and the outflow velocity of duct V ∞ . To solve the mentioned unknowns, a new expression of the upper and lower rotor lift with respect to the induced velocity is derived by complying with the blade element theory.

Blade element theory
The existence of the ceiling will also affect the calculation of the blade element theory of the upper and lower rotors. The blade element theory of an unconstrained environment assumes that the inflow is an axial flow, and only the flow along with the rotation and axial directions of the rotor should be considered. However, as impacted by the existence of the ceiling, the airflow into the upper rotor deflects significantly, which will lead to the radial velocity, and the radial flow should be considered additionally.
Complying with the derivation of momentum theory, the radial velocity should be considered directly in the blade element theory calculation of the upper rotor, while the inflow of the lower rotor in the wake area of the upper rotor can be considered the axial flow. In Dreier (2007) and Cao (2015), the blade element lift of a rotating rotor at the position r from the rotation axis is defined as: Where U(r) represents the relative velocity of airflow relative to the rotor, ρ denotes the air density, C L (α) expresses the lift coefficient of the blade element, and c(r) is the chord length of the blade. When the prototype generates the ceiling effect, U(r), as the overall velocity vector of the airflow, covers the rotor rotation speed U = r, the vertical speed U z and the radial speed U r . If U acts as a vector of U(r), it can be written as U = [U r U U z ] T . Figure 15 presents a schematic diagram of the blade element theory of the upper rotor with the ceiling. In Figure 15,r,φ andẑ represent the direction of the hub coordinate system respectively.n(r) expresses a unit vector perpendicular to the rotor surface at the radial r from the centre of the rotor, and the angle of attack α is enclosed by U andn(r). Specific to a fixed pitch rotor,n(r) can be projected into the hub coordinate system by two small angles, i.e.n(r) = [sin θ φ sin θ r cos θ φ cos θ r cos θ φ ] T . To calculate the angle of attack α, the small angle approximation principle is used, we can write α as cos(90 • − α) = U · n(r)/(|U||n(r)|). Accordingly: Moreover, in the case of small angle approximation: sin θ φ (r) ≈ θ φ (r), sin θ r (r) ≈ θ r (r), cos θ φ (r) ≈ 1, cos θ r (r) ≈ 1.
Accordingly: In Dreier (2007) and Cao (2015), the lift coefficient C L is proportional to the angle of attack.
Accordingly, Equation (1) can be written as follows: Since U U z and U U r , U(r), U and U are all replaced by r, U z is replaced by V 2 , and U r (r) is replaced by v r (r) = (r/2z)V 1 .
By applying three lumped parameters (a,b,c), a,b,c express the dimensionless propeller's blade coefficients related to the airfoil structure. According to the streamline diagram analysis of Figure 5 and Figures 7 and 8, the whole lower rotor is in the axial flow approximately. Accordingly, when the blade element theory is followed to calculate the lift of the lower rotor, it is identical to the conventional calculation. However, due to the different inflow velocities, the angle of attack of the part in the wake area of the upper rotor is not identical to that of the part outside the wake area, so it should fall to two parts and be calculated separately. The schematic diagram of the conventional blade element theoretical model is presented in Figure 16, where the distance from the element to the hub centre is r, the blade pitch is expressed as θ r1 , the vertical velocity component is denoted as U z1 , and the tangential velocity component is U 1 .
The vertical velocity component and the tangential velocity component of the blade element in the wake area of the upper rotor can be written below: Assuming that the blade pitch is linear torsion, the blade element pitch can be expressed as a function of the According to the small angle hypothesis, the following conclusions are drawn: The inflow angle of the blade element can be approximately written as follows: Next, the angle of attack α r1 of the blade element is defined as: The aerodynamic lift dL of the blade element is: The direction of aerodynamic lift dL is perpendicular to the inflow direction.
According to the hypothesis of small angle, cosϕ ≈ 1, it yields: The number of blades of the rotor is set to Nr. Subsequently, the average total lift generated by the rotor in a blade rotation cycle is written as: Likewise, the lift T 2−b of the lower rotor outside the wake area of the upper rotor can be obtained from the conventional blade element theory. Notably, the vertical and tangential velocity components of the blade element outside the wake area of the upper rotor can be written as: Thus, the calculated lift T 2−b of the lower rotor outside the wake area of the upper rotor is as follows: The lift T 1−b determined by the blade element theory that considers the ceiling effect is equated with the sum of the lift T 1−m calculated according to the momentum theory of the upper rotor and the extra lift of the upper rotor due to the ceiling effect.
Moreover, for the part of the lower rotor in the wake area of the upper rotor, the lift T 2−b1 calculated by the blade element theory is equated with the lift T 2−m1 calculated by the momentum theory. At the same time, for the part of the lower rotor outside the wake area of the upper rotor, the lift T 2−b2 is calculated by the blade element theory is equated with the lift T 2−m2 calculated by the momentum theory.
The lift T 1 and T 2 of the upper and lower rotors can be solved by adopting the simultaneous Equation (19) (41) (42) (43) and the iterative solution when the prototype vertical climbs at a constant rate with the ceiling effect. The simultaneous equation is written as Equation (44). The relevant derivation diagram is shown in Figure 17.
However, the lift of the prototype with the ceiling effect includes the lift of the upper and lower rotors, and the lift of the duct. No accurate theoretical analysis has been conducted on the lift of ducts. On the whole, the ratio of the duct lift to the prototype lift is expressed as a coefficient k aug , which is identified with the relevant methods.

Identification of lift coefficient of duct
To estimate the lift ratio of the duct, a static bench test can be carried out. We take the rotational speed as the control variable to analyze the alteration of the ratio of the duct to total lift with the rotational speed. For the convenience of analysis, the normalized speed is defined below: Where max is the maximum speed, and here is 7200 rpm/min. Figure 18 shows that the lift ratio of the duct is basically unchanged with the change of rotor speed when the prototype is without the ceiling effect. Since this study considers the uniform vertical velocity of the prototype under the ceiling effect, the rotor speed always equals to 4000 rpm/min in the hovering process. When the prototype is unconstrained, the lift proportion of the duct is obtained at the hovering speed ( Figure 18). k aug ≈ 0.43 (46) Figure 18 shows how the lift ratio of the duct varies with the distance from the rotor speed when the prototype is hovering in an unconstrained environment. The change of the duct lift ratio should be analysed when the prototype is vertically climbing at a constant rate under the ceiling effect. Likewise, take the hovering speed of 4000 rpm/min, and compare the lift ratio of different distances from the prototype to the ceiling and different vertical climbing rates of the prototype. According to Figure 19, the duct lift ratio will change significantly with the distance change from the prototype to the ceiling. When the prototype is close to the ceiling, the duct lift ratio with the vertical velocity is larger than in the hovering process. When the prototype is far away from the ceiling, the lift ratio of the duct with vertical velocity exceeds that in the hovering process. Due to the small range of vertical velocity variation in this study, it is reported that under the identical distance from the prototype to the ceiling, the effect of velocity on the lift ratio of the duct is significantly small. Therefore, this paper primarily considers the duct lift ratio is affected by the distance from the prototype to the ceiling. The change of lift ratio with distance can be approximately fitted as follows:  Where T is the total lift of the prototype system.

Model validation
The Simulink module of MATLAB software can be adopted to model the lift mechanism model in the third section. For the lift mechanism model, its accurate performance should be verified. In this section, the lift mechanism model of the prototype is verified by conducting the CFD numerical calculation and the bench test. The mentioned results can verify the following items: the lift variation process of the distance from the prototype to the ceiling; the effect of vertical velocity on the lift variation of the prototype under the ceiling effect.

Simulation settings
The specific settings of the CFD simulation are illustrated in Section 2.1, and the prototype lift characteristics with different distances from the upper rotor to ceiling and different vertical velocities are taken as the output.

Bench test settings
To verify the accuracy of the lift mechanism model more effectively, the bench test of the physical prototype is also considered to conduct a better comparative analysis. The built duct bench is shown in Figure 20. The whole test bench comprises a test duct, twinrotors, brushless DC motor, connecting bracket, fixed base, a series of sensors and data collection systems. The test duct is made by 3D printing technology, and its structure is identical to that of the flight prototype. To study the aerodynamic force of the prototype with the proximity effect, the test bench is also equipped with a constrained board to simulate the constraint environment. In the case of the ceiling effect, the constrained board is set in front of the duct lip and parallel to the plane of the rotor. At the same time, the distance between the constrained board to the prototype can be regulated by test requirements.
Four types of sensors are installed on the test bench to collect the relevant data. The force sensor is installed on both sides of the duct body to directly measure the lift of the duct. Two-component balance is capable of measuring the force and moment of a single shaft. They are installed in the lower mount of the coaxial counterrotating blade drive motor to measure the lift and reaction torque of the rotor. The rotary encoder is installed in the drive motor mounting base, which is adopted to measure the real-time rotor speed. Hall current sensor is penetrated by the power line of the driving motor, which is used to measure the current of the driving motor in real-time.
The measuring range and accuracy of each sensor are given in Table 2.  There are two kinds of output signals of the sensors configured on the bench, one is the force sensor, the twocomponent balance, and the Hall current sensor, which outputs 0 ∼ 5 V analog voltage signal; the other is the pulse signal output by the rotary encoder. To process the output signal of the sensors, the test bench is equipped with a data collection system, as shown in Figure 21. The data collection system is composed of the test bench and the remote monitoring terminal. The test bench is used to process the sensor signal. At the same time, the remote monitoring terminal is used to monitor and store the data. The communication between the test bench and remote monitoring terminal is completed by using the wireless transmission module.

Validation of model results
In this section, the simulation results of the lift model are compared with the CFD numerical results and the bench test results. The hover and vertical climbing rate of the prototype are also selected to explain.

Hovering condition verification
Under the hovering conditions, the rotational speed of upper and lower rotors is set at 4000 rpm/min in CFD numerical simulation and bench test, and V0 is given as 0 m / s in the built lift mechanism model. It is worth noting that besides the real lift signal, various noise signals are also measured in the bench test (e.g. the highfrequency noise signal in the circuit, the mechanical noise signal of the rotor rotation frequency, and the lowfrequency vibration noise signal of the test bench). Circuit noise and mechanical noise can be removed by the filter, as shown in Figure 22. The fixed base of the test bench in Figure 20 is made of aluminum alloy with a total mass of 50 kg. As a result of the weight of the fixed base, it can ensure that the low-frequency vibration noise can be ignored when the rotor is rotating in the normal range.  This study selects the bench test data of the prototype hovering for 14s in an unconstrained environment, filters the measured lifts of the upper and lower paddles and ducts to calculate their average values, and finally sums them as the total lift of the system. Likewise, the total lift of the system at different distances from the rotor to the ceiling can be calculated when the prototype generates the ceiling effect.
The lift coefficient of the prototype under the hovering condition can be obtained by the following formula: Where T is the total lift of the system. ρ = 1.225 kg / m3 is the air density, A = π R 2 is the area of the rotor, and = ω π / 30 is the speed of the rotor. All units are the international system of units. In Figure 22, C T t = C T u + C T l + C T d is the total pulling force, and the subscripts u, l, d, and t represent the upper rotor, lower rotor, duct and overall system respectively.
The measurement error of the sensor is the main error in the bench test, which includes the measurement error of rotor speed and lift. According to Equation (43), using Kline McClintock equation (Kine & McClintock, 1953), the uncertainty of lift coefficient can be obtained: Equation (44) can be simplified as: Given the test data, the uncertainty of the lift coefficient is within 4%, which satisfies the requirements of test accuracy.
Compare the test data with CFD numerical calculation results and lift model simulation results, and the comparison chart is shown in Figure 23.
As shown in Figure 23, the changing trends of lift characteristic curves obtained using three different methods are nearly consistent. In this study, the simulation results of the CFD numerical simulation of the prototype hovering in an unconstrained environment were taken as T OCE . Compared with the lift mechanism model simulation results, the CFD numerical results are closer to the bench test data. This result is explained as some assumptions and approximations are made when the lift mechanism model is being built, thereby causes some differences between the calculation and the reality. However, under the different distances between the upper rotor and the ceiling, the error rate of the bench test and the model simulation remains basically at 15%, which has a high degree of confidence, indicating that the built mechanism model exhibits a high accuracy.

Vertical climbing condition verification
To more effectively compare the accuracy of the model when the vertical velocity V0 is insufficient. Besides, to be consistent with the hovering condition, the CFD numerical simulation calculation and the bench dynamic test of the prototype with different vertical velocities under the ceiling effect should be conducted. However, the bench test in this study can only be employed for the static test. When the distance between the constrained board and the duct is regulated manually and dynamically, the vertical velocity can't be ensured to be accurate and uniform constantly. Moreover, the external noise of the system will increase, and the response of the sensor will be delayed, causing significant errors in the measured lift data and a decrease in the degree of confidence. As indicated from Figure 23, the trends of the lift characteristic curves of different distances from the upper rotor to the ceiling obtained through the numerical calculation and the static test bench are the same, and the errors are within 10%. It is suggested that the CFD simulation results are basically consistent with the test results, which demonstrates that the numerical calculation results are reasonable. Thus, the CFD numerical calculation method is adopted to verify the accuracy of the lift model under the subsequent vertical climbing conditions. As indicated from Figure 23 and Figure 24, when the vertical climbing rate of the prototype and the distance from the upper rotor to the ceiling are different, the CFD numerical simulation results and the mechanism model simulation results display the same changing trends. At the same time, the error between the numerical simulation results and the lift model simulation results is maintained within 10%, indicating a high consistency.
Complying with the hovering condition, the total lift of the system decreases with the increase in the distance from the upper rotor to the ceiling at different vertical climbing rates, and the larger the distance is, the less the lift decreases.
Under the identical distance from the upper rotor to the ceiling, with the increase in the vertical climbing rate of the prototype, the total lift of the system will decrease. This is because, with the increase in the vertical velocity, the induced velocity will decrease, thereby leading to a decrease in the total lift of the system. At the constant vertical velocity, the variation of the total lift of the prototype with the distance from the upper rotor to the ceiling is consistent with that in the hovering process, and the total lift increases with the decrease in the distance.

Conclusion
Based on several reasonably simplified assumptions, the momentum theory and the blade element theory are used to build a lift mechanism model for the prototype with the ceiling effect is vertically climbing at a constant rate. Moreover, the interaction between the upper and lower rotors, the effect of the upper rotor wake on the lower rotor, and the additional lift attributed to the ceiling effect are considered. The accuracy of the proposed model is verified based on CFD numerical calculation and the bench test results under hovering and vertical climbing conditions. The lift mechanism model is capable of indicating the relationship between the lift of the prototype and the distance from the upper rotor to the ceiling and the vertical climbing rate. When the prototype with constant vertical climbing rate, with the decrease in the distance from the upper rotor to the ceiling, the total lift of the prototype increases exponentially. At the same time, if the distance between the upper rotor and ceiling is set as a constant, the overall lift of the prototype decreases as increasing vertical climbing rate. By comparing the influences of distance and vertical climb rate on the prototype lift, we find that the influence of distance on the lift is much greater than the vertical climb rate. As a whole, we can conclude that the distance between the upper rotor and ceiling plays a more critical role in the influence of the lift when the prototype with the ceiling effect.
In general, this study presents a lift mechanism model of coaxial twin-rotor ducted aircraft for the cases of small vertical climbing rates with the ceiling effect. The characteristics of the flow field and lift variation for the prototype that moves close to the obstacle are researched and disclosed thoroughly, which can also contribute to future studies of the prototype motion control.
In subsequent studies, the dynamic slide rail with variable velocity can carry different constrained boards by complying with the task requirements to realize the dynamic process bench analysis. Moreover, the prototype can carry IMU and laser rangefinder sensors to measure its vertical climbing rate and real-time distance to the ceiling in the subsequent actual flight experiments. Then, under the same controller with or without the proposed lift mechanism model to verify that the model can be used to improve the control performance, which is also one of the future studies.
Furthermore, the lift mechanism modeling of the prototype in other constrained environments (e.g. the ground and wall environment, and the uniform vertical descent under the ceiling effect) is a field worth investigating.