**1. Introduction**

An Unmanned Aerial Vehicle (UAV) is an unmanned aircraft operated by radio remote control equipment and its program control device [1]. It can be divided into fixed-wing UAV [2], rotor UAV [3], parawing UAV [4], and flapping-wing UAV [5] according to the flight platform configuration. The fixed-wing UAV has the advantages of fixed-wing gliding, stronger "robustness" to control technology errors, enormous carrying capacity, more extended flight, etc. [6]. It is widely used in various fields, which increases the requirements for survival and recovery. The takeoff technology of UAVs is relatively mature, whereas the landing technology is still facing many challenges, especially the more difficult emergency landings. According to statistics, 60% of UAV accidents are caused by human operation errors, among which takeoff and landing phases account for 50%, and landing accidents account for the majority [7]. Therefore, one of the critical problems currently for UAVs is how they can be recovered safely and accurately [8].

At present, the typical recovery methods of UAVs are net collision recovery [9], cable hook recovery [10], parachute recovery [11], and vertical takeoff and landing UAV [12]. Parachute recovery has been widely adopted due to its simple overall layout, low level of technical difficulty, and low dependence on other systems. In the parachute recovery of UAVs, the deployment model can be divided into deployment by a crown chute or a moving mechanism [13]. The crown chute will take a long time to open the parachute,

**Citation:** Shi, W.; Yue, S.; Li, Z.; Xu, H.; Du, Z.; Gao, G.; Zheng, G.; Zhao, B. Dynamic Simulation and Parameter Analysis of Weaved Composite Material for Unmanned Aerial Vehicle Parachute Recovery in Deployment Phase. *Crystals* **2022**, *12*, 758. https://doi.org/10.3390/ cryst12060758

Academic Editors: Yong He, Wenhui Tang, Shuhai Zhang, Yuanfeng Zheng, Chuanting Wang and Per-Lennart Larsson

Received: 28 April 2022 Accepted: 24 May 2022 Published: 25 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

which is not conducive to emergency release. The moving mechanism can shoot the traction device quickly to ge<sup>t</sup> a fast opening speed. Whether it is a normal or emergency situation, the parachute can be opened reliably and maintain the safe recovery of the UAV. In respect of parachutes, the deployment phase and the inflation phase are always research hotspots in the aviation field. Compared to the research on the parachute inflation phase, the research on the deployment phase is more complicated, so there is little research literature. However, the essence of the parachute deployment process is rope dynamics, where the research has been a hot topic in recent years [14].

Composite materials are widely used in the aviation field because of their excellent strength-weight ratio [15,16], and the weaved composite materials are widely used in flexible mechanisms such as parachutes. The parachute deployment phase is short and complicated, involving several movements. Studies on this process first appeared in the 1970s. Wolf proposed the continuous deployment model [17], which assumed that the suspension line was always straight during the deployment phase. The calculated snatch force was equal to the actual situation. Shen Guanghui et al. [18] established a 6-DOF mathematical model of the Mars deceleration system. Ortega Enrique et al. [19] accurately calculated the structural load and stress during parachute inflation based on the filling-time expansion model and Ludtke's area law. Yu Li et al. [20] analyzed the deployment phase of the parachute and established a multi-particle dynamics model, focusing on the rope's flexible and significant deformation characteristics. Zhang Qingbin, Wang Haitao, et al. [21,22] systematically studied the spring damping model of rope and the influence of wake on the trajectory of the parachute; they applied it to the calculation and analysis of the parachute deployment phase, obtaining good results.

The researchers studied theoretical models of the deployment and inflation phases, but few studies considered the initial folding of the parachute. In this paper, a method of deployment of parachutes by tractor rocket is proposed to improve the survivability of UAVs in emergencies. Based on the explicit dynamics method, the dynamics model of the parachute deployment phase is established. Experiments were designed and carried out to verify the accuracy of the dynamic model by comparing the parachute deployment time and height. The influence of parachute weight and the launch temperature on the deployment performance was further studied. The above research provides a new idea and method for the design and simulation of the parachute deployment stage in the parachute recovery process of unmanned aerial vehicles.

The structure of this paper is as follows: Section 1 summarizes the relevant research and the main contributions of this paper. In Section 2, the dynamics model of the parachute deployment phase is established. In Section 3, the accuracy of the simulation model is verified by comparing it with the experiment. Parametric analysis of parachute weight and rocket launch temperature are included in Section 4. Section 5 concludes the whole research.

#### **2. Dynamic Model of the Parachute Deployment Phase**

#### *2.1. Overview of Parachute Ejecting Tractor Rocket*

Parachute recovery has a wide range of applications [23]. The mass of what is recovered ranges from a few kilograms to dozens of tons; the opening speed ranges from low to supersonic; and the height spans from thousands of meters to dozens of kilometers, which can safely and effectively slow the recovered UAVs down to the ground or water surface. In this paper, a typical fixed-wing UAV is selected, which adopts a conventional H-type double-tail-support layout, with a width of 7020 mm, a length of 4180 mm, a height of 1470 mm, and a weight of about 250 kg. It needs to be recycled within the height range of 150–5000 m. The temperature range is −40 ◦C–+45 ◦C. The carrying platform of the UAV recovery system cannot be subjected to excessive recoil force. The tractor rocket is used as the traction device to keep the overload within a reasonable range.

The working process of the tractor rocket is as follows: When the igniter receives instructions to light a solid propulsive agen<sup>t</sup> to produce a large amount of gas, the gas blow-out from the rocket nozzle jets makes the rocket accelerate. The rocket's move drives the extraction line, canopy, suspension lines, and sling to keep a straight direction. As the UAV descends, the parachute inflates until it lands smoothly.

#### *2.2. Model Description*

The parachute is the core part of the whole UAV parachute recovery system. The working process of the parachute includes two stages: parachute deployment, and inflation [24]. The UAV recovery system has strict weight and small swing angle requirements. It is necessary to select an appropriate parachute with a light weight to meet the specified swing angle due to its significant drag coefficient, small nominal area, and low volume. So the extended skirt parachute was chosen as the main parachute, which is shown in Figure 1.

**Figure 1.** Parachute structure.

There is also an extraction line and sling in the parachute system. The extraction line connects the tractor rocket to the parachute canopy, and the sling connects the suspension lines to the UAV. The canopy material is 1056 Anti-burning brocade silk material, the material of the extraction line and sling are 25–1000 Arlene belt, and the tensile force measured by the test is 19,860 N. The material of suspension lines is 2–200 high-strength polyethylene rope. The tensile strength is 2000 N. The main parameters of the parachute are shown in Table 1.

**Table 1.** Main parameters of parachute.


As the first action in the parachute recovery process, the deployment phase requires the parachute to complete a series of pre-set steps in a short time. As is shown in Figure 2. There are mainly two parachute deployment procedures: canopy-first deployment and lines-first deployment [25]. Under the same conditions, the snatch force of the canopy-first deployment method is larger than the lines-first deployment method, so the lines-first deployment method is adopted in this paper.

**Figure 2.** Deployment procedures. (**a**) The final state of a parachute deployment phase (**b**) Canopyfirst deployment (**c**) Lines-first deployment.

A finite element simulation model was needed to simulate the parachute deployment phase. The parachute uses a "Z" fold in the parachute container. Firstly, a single gore section model was established according to its size. Then, the individual gore sections were stacked and placed according to the constraints of each gore section position, with a certain angle and distance between each gore section. The partially enlarged view of the gore model is shown in Figure 3.

**Figure 3.** Gore model.

Then the gore model was divided, and the split length was consistent with the folding size of the parachute in the parachute container. The geometric parachute folded model was established by placing the parachute from the bottom to the top. The folded model was taken as the initial state of the parachute in the deployment phase. On this basis, an extractor rocket, suspension lines, and other components were added. The parachute in the final folded state is shown in Figure 4. The main parameters of the tractor rocket are shown in Table 2.

**Figure 4.** Simulation model.

**Table 2.** Main parameters of parachute.


#### *2.3. Establishment of Suspension Line*

The T3d2 element (two-node linear three-dimensional truss element) in ABAQUS was selected to construct the suspension lines model [26,27]. Element nodes connected the truss elements in the finite element model. When the truss elements are under compression loads, they will rotate around the nodes freely, making them unable to bear compression loads. Therefore, multiple truss elements connected by element nodes can simulate the dynamics performance of the rope. Each truss element can bear internal and external forces during the movement. The dynamic model of the suspension lines consists of a set of dynamic equations for a single truss element [28]. Figure 5 shows the force transmission of the truss element in the deployment process. The truss element is a double node element (node *i* and node *j*), and the unit axial force *Tij* of the element can be denoted as [29]:

$$T\_{ij} = k\_{ij} \left( l\_{ij} - l\_{ij}^0 \right) \tag{1}$$

where *lij* is the deformed length of the truss element *Hij*, parameter *lij*0 is the initial length, *Kij* is the equivalent stiffness of the truss element, which depends upon the material of the rope. The equivalent stiffness *Kij* is written as:

$$k\_{ij} = \frac{EA\_{ij}}{l\_{ij}^0} \tag{2}$$

where *E* is the elastic modulus of the rope material, *Aij* is the cross-section area of the rope.

The test was carried out on the SANS testing machine [30]. The measurement range for the force was between 1 KN and 100 KN, and the precision of force indication can be achieved within ±1%. Before the test, the end of the suspension line was tied with the upper and lower gripper of the testing machine, and the suspension line was in a relaxed state. The suspension line was stretched during the test, with the displacement and force sensor recording the real-time distance and tensile force. The initial length of the suspension line was 100 mm, and the maximum strain in the test was 0.2. The tensile stiffness tests for the suspension line are shown in Figure 6.

**Figure 6.** Tensile stiffness test of suspension line.

According to Figure 7, the line stiffness increases with the growth of external load. However, when the external load is under around 10 N, the relationship between the rope stiffness and the elastic modulus of the material can be simplified as linear:

$$K = \frac{EA}{L} \tag{3}$$

where *K* is the stiffness of the rope, and *E* is the elastic modulus of the rope material. *A* is the equivalent area of the rope, which is 0.196 mm2. *L* is the equal length, which is the original length of the rope (100 mm) in the tensile test.

**Figure 7.** The flowchart to determine line stiffness.

The procedure to determine the suspension line stiffness is listed as follows. Firstly, the initial maximum tension of the line segmen<sup>t</sup> is given, and the constant tensile stiffness is determined, based on Figure 6. After that, the dynamic model of suspension line deployment is executed, and the new maximum tension during deployment is obtained. Based on this, the new tensile stiffness was reselected, based on Figure 6. The above steps were repeated until the maximum tension between two iterations was lower than 0.1 N. Then, the current tensile stiffness was determined as the ultimate rope stiffness for the following analysis. The iterative flowchart is shown in Figure 7.

The initial maximum tension of the rope was first given as 10 N, in which the stiffness of the suspension line is 2.72 N/mm. Based on Equation (3), the elastic modulus was calculated as 1385.3 MPa. Based on Figure 7, the iteration was performed to obtain the tensile stiffness of the suspension line. After eight iterations, the ending condition was satisfied, and the whole loop was finished. The final force curve of the suspension line during deployment is shown in Figure 8.

From Figure 8, it can be seen that the maximum tension suspension line was 4.5 N. Therefore, the stiffness of the suspension line was 1.47 N/mm. The chosen stiffness was shown in Figure 9, and the time for maximum tension to appear ranged around 1.57 s. The suspension line was in full deployment and started straight state. Equation (3) obtains the elastic modulus of the suspension line, and the value is 748.7 MPa. The elastic moduli computed are equivalent to the corresponding structure, which can be applied to the suspension line structure.

**Figure 8.** Line segmen<sup>t</sup> tension during deployment.

**Figure 9.** Result of stiffness selection.

#### *2.4. Establishment of Canopy*

The membrane element is a thin plate that can withstand the membrane force without bending or transverse shear stiffness, so the only non-zero stress components in the membrane were those parallel to the middle surface of the membrane. The membrane was in a state of plane stress. Discretization and meshing obtains a combination of quadrilateral elements with four nodes. In the local coordinate system, according to the theoretical hypothesis of the finite element method, the displacement vector of any node of the membrane element is *d* [31,32]:

$$\{d\} = \{u, v, w\}^T\tag{4}$$

where *u*, *v*, and *w* are all functions of three coordinate directions in the global coordinate system. Therefore, the displacement of this membrane element in the local coordinate system can be expressed as:

$$\mathbf{x}\{\mathbf{x}\} = \left\{ \mathbf{x}\_i, y\_i, z\_i, \mathbf{x}\_{j\prime} y\_j, z\_{j\prime}, \mathbf{x}\_{m\prime} y\_{m\prime} z\_{m\prime}, \mathbf{x}\_{n\prime} y\_{n\prime} z\_n \right\} \tag{5}$$

For quadrilateral elements, the form function can be expressed in local coordinates as:

$$N = \begin{bmatrix} IN\_{i\prime} \, IN\_{j\prime} \, IN\_{\text{m}} \, \end{bmatrix} \tag{6}$$

where *I* is the element matrix of size 4 × 4.

According to the finite element theory, the relation between the displacement function and the node displacement is:

$$\{d\} = \mathbf{N} \cdot \{\mathbf{x}\}^T \tag{7}$$

Due to the nonlinear characteristics of the membrane material, the green strain tensor was adopted, and the strain matrix of the membrane element [*B*] could be expressed in two parts:

$$[B] = [B\_0] + [B\_N] \tag{8}$$

where [*B*0] is the linear part of the strain matrix; [*BN*] is the nonlinear part of the strain matrix, which is a function of the node displacement vector {*d*}. On this basis, the virtual displacement principle was used, and the finite element was written in the incremental form as follows:

$$\{\mathbf{K}\} \cdot d\{d\} = d\{P\} = \int\_{V} [\mathbf{B}]^T \cdot d\{\boldsymbol{\sigma}\} d\boldsymbol{v} + \int\_{V} d[\mathbf{B}]^T \cdot \{\boldsymbol{\sigma}\} d\boldsymbol{v} \tag{9}$$

According to mechanical theory, [*B*0] has nothing to do with node displacement, so the following formula can be obtained:

$$d\{\varepsilon\} = d(\left[B\right] \cdot \{d\}) = \left[B\right] \cdot d\{d\} + d\left[B\_N\right] \cdot \{d\} \tag{10}$$

The final finite element incremental equation can be expressed as:

$$d\left(\left[K\right] + \left[K\_{\mathcal{T}}\right]\right) \cdot d\left\{d\right\} = d\left\{P\right\} \tag{11}$$

where [*<sup>K</sup>σ*] is the geometric stiffness matrix.

The tensile test of the canopy was carried out to obtain the weft and warp tensile modulus of the canopy material. The thickness of the specimen was 0.1 mm, the width was 10 mm, and the length was 450 mm. As is shown in Figure 10, the weft modulus was calculated to be 7.29 MPa, and the warp modulus 10.75 MPa.

**Figure 10.** Canopy tensile test.

#### *2.5. The Tractor Rocket Thrust Test*

In this paper, a tractor rocket is used to deploy the parachute. The rocket thrust tests are carried out at different temperatures to obtain the thrust parameters. The test system consists of a testbed, rocket, sensor, data transmission lines, and data collector, as is shown in Figure 11. The experiment proceeded three times, one at standard temperature (+20 ◦C), one at low temperature (−40 ◦C), and one at high temperature (+45 ◦C). Before the test, the rocket was kept in an incubator for more than twelve hours. During the trial, the rocket was taken out and installed quickly. The ambient temperature was 20 ◦C, and the sampling frequency was 10 kHz.

**Figure 11.** The tractor rocket thrust test layout.

The test results are shown in Figure 12 and Table 3. It can be seen that with the temperature increase, the peak thrust value, average thrust value, and peak time of the rocket all increase, and the working time of the rocket decreases.

**Figure 12.** Curves of rocket thrust.

**Table 3.** Main parameters of parachute.


#### **3. Experimental Validation**

#### *3.1. Experiment Design*

The parachute deployment experiment was designed to verify the correctness of the simulation results. The launch temperature of the tractor rocket was 20 ◦C. The test was carried out on a clear and windless day to facilitate the test record and avoid the impact of environmental factors. The test layout is shown in Figure 13.

In this test, the parachute weight was 10 kg. Before the test, the rocket was kept at +20 ◦C for more than twelve hours. Next, test parts were trialed together and placed on the test site, and the high-speed camera was set up in front of the parachute container to record time and parachute track. The aerial drone hovered above to record the test process. When the rest of the preparation work was completed, the rocket was removed from the oven, installed in the rocket cabin, and fired quickly to ensure that the rocket temperature was around +20 ◦C.

**Figure 13.** Experiment layout.

#### *3.2. Comparison between Simulation and Experiment*

The simulation and experiment deployment process configuration at different times is obtained, as shown in Figure 14.

**Figure 14.** Comparison of deployment configuration between simulation and experiment.

The displacement of the tractor rocket directly reflected the deployment degree of canopy and suspension lines. As is shown in Figure 15, the tractor rocket reached its maximum displacement of 29.91 m in 0.8 s. Observation of the test video shows that the time for the parachute to be deployed was about 0.93 s, with a time error of 13.97%. The nominal radius of the canopy was 7.78 m; when it was deployed straight, the measured length was about 11 m. The length of the suspension lines, extraction line, and sling was 14 m, 1 m, and 5 m, respectively. Finally, the height of the rocket was 0.36 m. Therefore, the displacement of the rocket under ideal conditions was about 31.36 m, and it can be seen from the test video that both the canopy and the suspension lines were straight, so the perfect length can be used as the test length comparison. The total displacement error was 4.62%. The main reason for the error was that the friction between the canopy, suspension lines, parachute container, and the air resistance was not considered in the simulation. Figure 16 shows the change in the rocket velocity. It can be seen that the rocket began to accelerate, and with the parachute weight increase, the rocket speed decreased. Because the canopy was folded three times, there were three distinct fluctuations in the rocket's velocity. With the canopy and suspension lines fully straight, the rocket finished its work, and the gravity slowed it down until it fell back to the ground.

**Figure 15.** Displacement of tractor rocket.

**Figure 16.** Velocity of tractor rocket.

In general, the simulation results matched the test results in trends and magnitudes, validating the accuracy and reliability of the dynamic model, which provided a solid foundation for the following parameter analysis of deploying performance.

#### **4. Parameter Analysis of Parachute Deployment**

#### *4.1. Influence of Parachute Weight*

This simulation researched the deployment phase of parachutes with different weights of 10 kg and 20 kg under the launch temperature of +20 ◦C. The simulation process is shown in Figure 17. It can be concluded that the rocket deployed parachutes successfully with a relatively smooth phase and a good result.

The force curves of the extraction line and the sling are shown in Figure 18a,b. The force of the extraction line and the sling varied in the same way in the two conditions. The maximum snatch force of the parachute with a weight of 10 kg was greater than 20 kg, which is because the more weighted parachute could offset more rocket thrust, resulting in a smaller force acting on the extraction line and the sling. The strength of the extraction line and of the sling met the design requirements.

**Figure 18.** Deployment performance under different weight. (**a**) Extraction line force (**b**) Sling force (**c**) Suspension lines force (**d**) Tractor rocket displacement.

The force curves of the suspension lines are shown in Figure 18c. The maximum snatch force of the parachute with a weight of 20 kg was greater than 10 kg, and the lower the weight, the less the cord fluctuated. Meanwhile, the strength of the suspension lines met the design requirements.

The displacement curves of the rocket are shown in Figure 18d. The maximum displacement of the 10 kg parachute was 29.91 m at 0.8 s, and that of the 20 kg parachute was 24.58 m at 1.15 s. The error of time and displacement was 23.66% and 21.6%, respectively. Thus, with the increase in the parachute weight, the degree of deployment reduced.

To sum up, as the parachute weight increased, the maximum snatch force on the extraction line and the sling decreased, but that on the suspension lines increased, and the deployment effect reduced. Therefore, the 10 kg parachute is the most effective.

#### *4.2. Influence of Temperature*

This simulation studied the deployment phase of a rocket under −40 ◦C, +20 ◦C, and +45 ◦C, with a parachute weight of 10 kg. The deployment simulation phase is shown in Figure 19. The canopy and the suspension lines were deployed gradually from the initial folding state, and the whole process was relatively smooth. After the rocket was fired, a continuous thrust pulled the parachute and the suspension lines out of the container. As the rocket rose, the parachute and the suspension lines straightened. When the propellant was burned out, the parachute and rocket gradually fell to the ground due to gravity. During the test, the deployment process of the parachute was smooth and orderly, with a good parachute shape and no strong disturbance. As the temperature rose, the parachute took less time to become straight.

**Figure 19.** Different temperature parachute deployment process.

The force curves of the extraction line, the sling, and the suspension lines are shown in Figure 20a–c. The snatch forces varied in the same way at different temperatures. For the extraction line, at +45 ◦C, the maximum force was 3814.98 N at 0.28 s; at +20 ◦C, the maximum force was 713.55 N at 0.16 s; and at −40 ◦C, the maximum force was 2857.3 N at 0.18 s. The force at +45 ◦C was the largest and was little different at the other two temperatures. This is mainly because with the temperature increase, the peak thrust value of the rocket changed significantly. The peak thrust value was 3199.6 N at +45 ◦C, and only 1991 N and 1318.3 N at +20 ◦C and −40 ◦C, respectively, leading to a change in the extraction line.

In the sling and suspension lines, it was evident that as the temperature increased, the force increased, although the maximum force appeared earlier, with higher requirements for the strength of the material of the sling and suspension lines. The reason is mainly that the average thrust of the rocket increased as the temperature rose. The results accord with the variation trend of thrust peak value and mean value of rocket at different temperatures. At the same time, the other parameters in the parachute system remained unchanged.

**Figure 20.** Deployment performance under different temperatures. (**a**) Extraction line force (**b**) Sling force (**c**) Suspension lines force (**d**) Tractor rocket displacement.

The rocket displacement directly reflected the deployment degree of the parachute. In Figure 20d, at −40 ◦C, the maximum displacement of the rocket reached 28.35 m at 0.91 s, and the error of displacement was 9.6%. At +45 ◦C, the maximum displacement of the rocket was 29.91 m at 0.79 s, and the error of displacement was 4.6%. It can be concluded that as the temperature increased, the deployment time became shorter, and the deployment length did not change by much. The results accord with the rocket's working time and thrust variation trend at different temperatures. At these three temperatures, the difference between the maximum and minimum values of rocket displacement was only 1.559 m, accounting for 4.97% of the ideal full length. Therefore, it can be found that temperature has little influence on the degree of parachute straightening, which provides a theoretical basis for the normal operation of the UAV parachute recovery system at different temperatures.

It can be seen from the above analysis that with the temperature increase, the deployment time became shorter, and the deployment length of the canopy and suspension lines only changed slightly. The maximum snatch force on the extraction line, sling, and suspension lines increased.

## **5. Conclusions**

This paper proposed deploying the parachute by the tractor rocket to make the unmanned aerial vehicle land safely in an emergency. The deployment process of unmanned aerial vehicle parachute recovery was analyzed by combining experiment and finite element methods. Then, the dynamic characteristics of deployment process were studied and researched. Finally, the effects of parachute weight and rocket launch temperature on the

deployment characteristics were discussed, based on the dynamics model. The conclusions are as follows:


**Author Contributions:** Conceptualization, W.S. and S.Y.; Software, W.S.; Validation, Z.L., H.X., Z.D., G.G., G.Z. and B.Z.; Writing—Original Draft Preparation, Z.L., H.X. and W.S.; Writing—Review and Editing, W.S.; Supervision, S.Y. and G.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (52102436), the Fundamental Research Funds for the Central Universities (30920021109), the Natural Science Foundation of Jiangsu Province (BK20200496), the China Postdoctoral Science Foundation (2020M681615), and the project of Key Laboratory of Impact and Safety Engineering (Ningbo University), the Ministry of Education(CJ202107), and the State Key Laboratory of Mechanics and Control of Mechanical Structures (Nanjing University of Aeronautics and astronautics) (Grant No. MCMS-E-0221Y01).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are contained within the article. The data presented in this study can be seen in the content above.

**Acknowledgments:** We wish to express our gratitude to the members of our research team, Wenhui Shi, Shuai Yue, Zhiqian Li, Hao Xu, Zhonghua Du, Guangfa Gao, Guang Zheng and Beibei Zhao.

**Conflicts of Interest:** The authors declare no conflict of interest.
