NOMENCLATURE

A, B, C

rigid bodies

b

span

b _A

span of centre wing

b _out

span of section outwards from hinge

c

chord length

C _Dparachute

drag coefficient of parachute

C_lp

rolling moment coefficient derivative by roll rate

C_lr

rolling moment coefficient derivative by yaw rate

C_L,W,max

maximum lift coefficient of main wing

C_mq

pitching moment coefficient derivative by pitch rate

C_np

yawing moment coefficient derivative by roll rate

C_nr

yawing moment coefficient derivative by yaw rate

C_Yr

side force coefficient due to yaw rate

DNS

direct numerical simulation

e

evaluation function

E _contributing

evaluation functions contributing to total minimum sigma level

E* _contributing

evaluation functions contributing to minimum sigma level

F

deployment torque scale; ratio of torques between deployment mechanism andaerodynamic force

F_{Z Aw}

z-direction force on wing of rigid body A

g

acceleration due to gravity

h

height

ID _Disp

identification number of the target dispersive input variable

ID _Eval

identification number of the target evaluation function

I_XX , I_YY , I_ZZ

moment of inertia

I_XZ

product of inertia

l

maximum lift force per unit wing span

LES

large eddy simulation

LSL

lower specification limit

m _total

total mass

M_b/2

bending moment on the centre of the main wing

M _{Reac Req}

requirement of the hinge reaction moment

n

sigma level

n _limit

limit load factor

n _min

minimum sigma level

n _{min Total}

total minimum sigma level

P, Q, R

angular rate

q

dynamic pressure

RANS

Reynolds-averaged Navier-Stokes equations

Re

Reynolds number

S _out

wing area of section outwards from hinge

S _parachute

reference area of parachute

t

time

T _aero

torque on hinge due to aerodynamic forces

T _deploy

torque due to deployment mechanism

T _spring|_Dep

torque of spring while wings are in deployed state

T_SL

left wing deployment delay time

T_SR

right wing deployment delay time

U, V, W

translational velocity

USL

upper specification limit

V _contributing

dispersive input variables contributing to total minimum sigma level

V _terminal

terminal velocity

X, Y, Z

coordinates of centre body

X ₀, Y ₀, Z ₀

inertial frame of reference

x

dispersive input variables vector

∆

difference for sensitivity calculation

θ _dep

deployment angle

μ_e

average of the evaluation function e

μ_x

average vector of dispersive input variables

μ_xi

average of i-th dispersive input variable

ρ

atmospheric density

σ_e

standard deviation of evaluation function e

σ_xi

standard deviation of i-th dispersive input variable

φ, θ, ψ

Euler angle

2.0 MISSION SCENARIO AND SIMULATION TARGET

Figure 1 shows a mission scenario and the simulation target. The Mars aircraft will be transported to Mars packed in an entry capsule. At Mars, the entry capsule will enter and descend into the Martian atmosphere with a parachute. When the entry capsule arrives at a predefined position, the Mars aircraft will be released from the entry capsule and will perform aerial deployment. The tail will be deployed first, and then the main wings will be deployed. The fully deployed Mars aircraft will control its attitude and start a horizontal flight. The Mars aircraft is able to conduct observations such as Martian magnetic field observations, thus capturing vast areas of the surface.

Figure 1. Mission scenario and simulation target.

Manoeuvres after the deployment, i.e. pull-up and observation flight, can be analysed using a usual flight simulation method that treats the aircraft as a single rigid body. However, the deployment process cannot be analysed using this method. The analysis of the deployment process requires a multibody dynamics simulation method. A motion analysis of the deployment process is important because the final condition of the deployment process coincides with the initial condition of the pull-up and observation flight. In addition, this study is focused on the design of the deployment mechanism of the main wing. Therefore, the simulation was performed from the end of the tail deployment to the end of the main wing deployment.

In the past, design of the deployment mechanism for a folding wing was carried out through a trial-and-error process. Therefore, the effect of the design on the deployment motion has not been clarified quantitatively. The simulation in this paper reveals how the design variables physically and quantitatively affect the aerial deployment behaviour.

3.0 SIMULATION METHOD

Although the Mars aircraft deploys not only the main wing but also the tail, this study deals with the deployment of the main wing. The analysis model consisted of three rigid bodies: a centre body, a right wing and a left wing. These bodies will be referred to as rigid bodies A, B and C in this paper. The centre body consisted of a centre wing, a fuselage and a tail. The right and left wings were assumed to be flat plates. Each body was connected to each other with spring-loaded hinges.

Figure 2 shows a wing deployment process. At first, the right and left wings were folded under the fuselage. Then both wings rotated around the hinges. The rotation stopped when the wings aligned with the centre wing. In this study, the left wing always started its deployment at the beginning of the simulation. However, to avoid contact between the wings, the right wing did not start its deployment simultaneously; the right wing started its deployment after the deployment delay time T_SR . The hinge axis was set along the chord direction because this configuration allows for the wings to be placed in parallel when in the folded state.

Figure 2. Wing deployment process.

This study used the aerial deployment motion simulator in Reference Fujita, Motoda and NagaiRef. 7. As described above, the model was a multibody system consisting of the three rigid bodies connected with hinges. The equations of motion for this multibody system were reported in Reference Fujita, Motoda and NagaiRef. 8. The equation of motion of the folding wing aircraft was derived using a velocity transformation method^{(Reference Tajima11)}. External forces were the aerodynamic force, the hinge torque and gravity. Therefore, the motion was affected by the aerodynamic force. However, the aerodynamic force is a function of the velocity relative to the air. Hence, the aerodynamic force was also affected by the motion.

From an aerodynamics perspective, the aerial deployment process is a complex problem. The points to consider include the low Reynolds number flow (Re = O(10⁴)), the flow separation, the process time (several seconds), the three-dimensional body and flow, and the wing motion. A potential-flow-based method and a method based on the Reynolds-averaged Navier-Stokes equations (RANS) method usually cannot treat the separated flow. On the other hand, although the large eddy simulation (LES) and the direct numerical simulation (DNS) can treat the separation, they cannot calculate the flow around the three-dimensional aircraft with wings in motion for such a long period of time.

Therefore, this study applied the experimental results of the aerodynamic characteristics of a flat plate to the wings and the tails as basic characteristics. The total aerodynamic characteristics (i.e. data from all 360° of angles-of-attack) of the flat plate were obtained by Okamoto et al^{(Reference Okamoto and Azuma12)} and were used for the longitudinal characteristics. The effect of a sideslip was evaluated based on Reference Hoerner and BorstRef. 13. It should be noted that the interaction of the flows around the bodies was ignored. Therefore, the aerodynamic characteristics were independent of the deployment angle. Although the deployment angle was independent, the deployment speed affected the aerodynamic force through the change of the air speed, angle-of-attack and sideslip angle.

The dynamic derivatives were not introduced for each flat plate. However, some of the derivatives were approximately introduced as a nature of this multibody dynamics. The aerodynamic force of each flat plate was independently calculated using the local air velocity vector at a reference point of each plate. On the other hand, one of the typical estimation methods for the dynamic derivatives of the aircraft calculates the dynamic derivatives based on the local difference of the aerodynamic force due to local air velocity vector using the strip theory^{(Reference Kato, Oya and Karasawa14)}. Therefore, the effect of such dynamic derivatives was automatically and approximately reflected in the simulation. The derivatives produced by the difference of the forces between the right and the left wing due to its rotation (i.e. C_lp , C_lr , C_np and C_nr ), and the derivatives produced by the change of the forces acting on the tail due to its rotation (i.e. C_mq , C_lr , C_nr and C_Yr ), were approximately effective. It should be noted that the design dihedral in this simulation is zero. However, if the dihedral exists, additional derivatives such as C_{l β} , etc. are also approximately effective.

In summary, this model can simulate a post-stall condition fast enough using the experimental result of the aerodynamic characteristics. Some of the dynamic derivatives were approximately introduced as a nature of this multibody dynamics. This model can also simulate the three-dimensional behaviour of the aircraft with moving wings.

A drive power was affected on the hinges as a torque owing to the deployment mechanism. Owing to the characteristics of the spring, the strength of a spring was at its weakest when the wing was in the deployed state. The torque of the spring while the wings are in the deployed state T _spring|_Dep shall be used as the standard in this simulation. In order to design an appropriate strength of the spring under a variety of deployment conditions, it was considered appropriate to relate the strength of the spring with a torque owing to the aerodynamic forces. A deployment torque scale F was, therefore, defined in this study as a dimensionless quantity expressed by the ratio of the torque due to the deployment mechanism T _deploy to the torque due to the aerodynamic forces T _aero

(1) $$\begin{equation} F \equiv \frac{{{T_{{\rm{deploy}}}}}}{{{T_{{\rm{aero}}}}}} = \frac{{{{\left. {{T_{{\rm{spring}}}}} \right|}_{{\rm{Dep}}}}}}{{{T_{{\rm{aero}}}}}} = \frac{{2{{\left. {{T_{{\rm{spring}}}}} \right|}_{{\rm{Dep}}}}}}{{l{b_{{\rm{out}}}}^2}} = \frac{{2{{\left. {{T_{{\rm{spring}}}}} \right|}_{{\rm{Dep}}}}}}{{qc{C_{L,\,W,\,{\rm{max}}}}{b_{{\rm{out}}}}^2}} = \frac{{2{{\left. {{T_{{\rm{spring}}}}} \right|}_{{\rm{Dep}}}}}}{{q{S_{{\rm{out}}}}{b_{{\rm{out}}}}{C_{L,\,W,\,{\rm{max}}}}}}, \end{equation}$$

where S _out and b _out are the wing area and span of the outer wing respectively, l is a maximum lift force per unit wing span, c is a chord length and q is a dynamic pressure. The maximum lift coefficient of the main wing C_{L,W, max} was assumed to be 1 based on the aerodynamic characteristics of the flat plate^{(Reference Okamoto and Azuma12)}. The deployment mechanism torque was obtained as a function of the deployment torque scale F. In this study, the deployment mechanism torque was treated as a value that was independent of the deployment angle to reduce the parameter because a change in the deployment mechanism torque due to deployment angle was not essential.

Atmospheric density and gravitational acceleration can be expressed as a function of altitude. The atmospheric density on Mars was calculated based on Ref. 15. The gravitational acceleration was calculated using the law of gravity. The incident angle of the main wing was 0°. The incident angle of the horizontal tail was set to 5° to make the trim angle-of-attack positive. Here, the positive direction of the incident angle of the horizontal tail is the direction in which the leading edge of the horizontal tail aimed downward. The trim angle-of-attack was 8.8° and the true airspeed at the equilibrium gliding condition was 65 m/s. The quaternion was used for the attitude expression in this simulation and the Euler angle was used only for the output. Figure 3 shows a definition of the Euler angle. It should be noted that, in this study, the downward direction was set to the direction to which all Euler angles were zero. It means that the X ₀, Y ₀ and Z ₀ directed to the ground, the east and the south, respectively. The reason why the reference of the Euler angle changed is that the initial condition of this simulation was mainly the downward dropping condition as it will be later described; this attitude was a singular attitude for the usual Euler angle definition. An Adams-Bashforth-Moulton PECE solver was used as the numerical integration method^{(Reference Shampine and Gordon16)}.

Figure 3. Definition of the Euler angle.

4.0 ROBUSTNESS EVALUATION

A sigma level methodReference Shimoyama, Oyama and Fujii ^{(17 ,Reference Shimoyama, Jeong and Obayashi18)} was adopted for the robustness evaluation. The sigma level is a function of the average and the standard deviation of the evaluation function. The average and the standard deviation were calculated using the sensitivity method. The sensitivity was numerically obtained using some of the results from the simulations. Evaluation functions in this study were margins for safe deployment conditions, which were defined in this section.

A sigma level method is one of the suitable ways to evaluate the robustness for the problem, which has constraints for the evaluation function. The sigma level n indicates the probability of the evaluation function staying within the constraints. The sigma level n was defined as follows:

(2) $$\begin{equation} n \equiv {\rm{min}}\left( {\frac{{{\mu _e} - LSL}}{{{\sigma _e}}},\,\,\,\frac{{USL - {\mu _e}}}{{{\sigma _e}}}} \right), \end{equation}$$

where LSL and USL are the lower and upper specification limits respectively, and μ_e and σ_e are the average and the standard deviation of the evaluation function e, respectively. A concept of the sigma level is shown in Fig. 4. The sigma level n is a measure of the distance between the average and the specification limits based on the standard deviation. Therefore, the sigma level n can be graphically understood as an inverse of the slope of the line connecting the lower/upper specification limits with the result of the average and the standard deviation of the evaluation function. For example, the solution of the lowest standard deviation is Sol. 3. However, because the average of Sol. 3 is close to the specification limit, the possibility of failure is not the lowest, and therefore the sigma level is not the highest. Solution 1 is worse than Sol. 2 because the average is closer to the limit, and the standard deviation is larger. The solution of the highest sigma level, and therefore of the highest robustness, is Sol. 2, owing to the balance of the average and the standard deviation. As described above, usage of the sigma level is effective to evaluate the robustness of the evaluation function, which has constraints. Table 1 shows the relation between the sigma level n and its safety probability^{(Reference Shimoyama, Jeong and Obayashi18)}.

Table 1 Relation between the sigma level and its safety probability^{(Reference Shimoyama, Jeong and Obayashi18)}

Figure 4. Concept of the sigma level n.

The sigma level method uses the average and the standard deviation of the evaluation function caused by the variation of the dispersive input variables. There are two calculation methods of the average and the standard deviation: the sensitivity method and the Monte Carlo method^{(Reference Shimoyama19)}. The sensitivity method numerically calculates the sensitivity to the evaluation function, and obtains the average and the standard deviation of the evaluation function using the Taylor expansion. The Monte Carlo method collected the result of the simulation by randomly changing the dispersive input variables, and then the average and the standard deviation were calculated based on the obtained data set. The sensitivity method is better from the standpoint of the balance of accuracy and calculation cost when the number of the dispersive input variables is low. Since the number of dispersive input variables for this study was low, as will be described later, the sensitivity method was adopted.

The j-th evaluation function e_j (x) of the dispersive input variables vector x can be written using the Taylor expansion, as follows:

(3) $$\begin{equation} {e_j}({{\bf x}} + \Delta {{\bf x}}) = {e_j}({{\bf x}}) + {\left( {\frac{{\partial {e_j}}}{{\partial {{\bf x}}}}} \right)^T}\Delta {{\bf x}} + \frac{1}{2}\Delta {{{\bf x}}^T}\frac{{{\partial ^2}{e_j}}}{{\partial {{{\bf x}}^2}}}\Delta {{\bf x}} + O(\Delta {{{\bf x}}^3}). \end{equation}$$

The average and the standard deviation of the evaluation function μ_{e j} and σ_{e j, i} respectively, can be obtained using Equation (3). On collecting the terms up to the first order, we obtain:

(4) $$\begin{equation} {\mu _{{e_j}}} = {e_j}({{{\bf \mu }}_{{\bf x}}}), \end{equation}$$

(5) $$\begin{equation} {\sigma _{{e_{j,\,i}}}} = \sqrt {{{\left( {\frac{{\partial {e_j}}}{{\partial {x_i}}}} \right)}^2}{\sigma _{{x_i}}}^2} , \end{equation}$$

where μ_x is the average vector of the dispersive input variables and σ_xi is the standard deviation of the i-th dispersive input variable x_i . Both μ_x and σ_xi are defined by the user. The first-order approximation was used to reduce the calculation cost.

The sensitivity of the dispersive input variables to the evaluation function was numerically obtained. It is expressed using the difference method as follows:

(6) $$\begin{equation} \frac{{\partial {e_j}}}{{\partial {x_i}}} \approx \frac{{{{\left. {{e_j}({{\bf x}})} \right|}_{{x_i} = {\mu _{{x_i}}} + \Delta {x_i}}} - {{\left. {{e_j}({{\bf x}})} \right|}_{{x_i} = {\mu _{{x_i}}}}}}}{{\Delta {x_i}}}, \end{equation}$$

where μ_xi is the average of the i-th dispersive input variable x_i and ∆x_i is a difference of x_i for the sensitivity calculation.

In this study, the sigma level n was a function of the design variables, evaluation functions and dispersive input variables. Therefore, it was calculated for each design point, evaluation function and dispersive input variable.

(7) $$\begin{equation} n = f({\rm{Design}}\;{\rm{variables}},\;I{D_{{\rm{Eval}}}},\,I{D_{{\rm{Disp}}}}), \end{equation}$$

where ID _Eval and ID _Disp are the identification number of the target evaluation function and dispersive input variable, respectively.

In this paper, a minimum sigma level n _min was defined as the minimum value of the sigma level of all evaluation functions for each design point and each dispersive input variable. Thus, the minimum sigma level n _min was a function of the design variables and the dispersive input variables. Here, evaluation functions contributing to the minimum sigma level E* _contributing can be obtained as a function of the design variables and the dispersive input variables.

(8) $$\begin{equation} {n_{{\rm{min}}}} = f({\rm{Design}}\;{\rm{variables}},\;\,I{D_{{\rm{Disp}}}}).\vspace*{-20pt} \end{equation}$$

(9) $$\begin{equation} E_{{\rm{contributing}}}^ * = f({\rm{Design}}\;{\rm{variables}},\;\,I{D_{{\rm{Disp}}}}). \end{equation}$$

Furthermore, a total minimum sigma level n _{min Total} was defined as a minimum value of the sigma level of all evaluation functions and all dispersive input variables for each design point. Thus, the total minimum sigma level n _{min Total} was a function of only the design variables. Here the evaluation functions and the dispersive input variables contributing to the total minimum sigma level, E _contributing and V _contributing, can be obtained as a function of only the design variables.

(10) $$\begin{equation} {n_{{\rm{min}}\;\,{\rm{Total}}}} = f({\rm{Design}}\;{\rm{variables}}).\vspace*{-20pt} \end{equation}$$

(11) $$\begin{equation} {E_{{\rm{contributing}}}} = f({\rm{Design}}\;{\rm{variables}}).\vspace*{-20pt} \end{equation}$$

(12) $$\begin{equation} {V_{{\rm{contributing}}}} = f({\rm{Design}}\;{\rm{variables}}). \end{equation}$$

5.0 SIMULATION CONDITION

5.1 Design variables

Since fundamental specifications (e.g. planform of the aircraft) were designed depending on the total performance of the aircraft, such variables were constant in the analysis of the robustness. The design variables for this robustness analysis were the variables that were directly related to the aerial deployment: the deployment torque scale F and the right-wing deployment delay time T_SR . The deployment torque scale F has been defined by Equation (1), and it represents the ratio of torques between the deployment mechanism and the aerodynamic force. The range of the deployment torque scale F was set from 0 to 1.5. On the other hand, the right-wing deployment delay time T_SR is the time difference between the starting time of the deployment of the right and the left wing. The range of the right-wing deployment delay time T_SR was set from 0 to 2 s. It should be noted that the left wing always started its deployment at the beginning of the simulation. There were 21 design points for each design variable at regular intervals. The total number of the design points was 21 × 21 = 441.

5.2 Evaluation functions

In this study, the evaluation functions were margins for the safe deployment conditions. The safe deployment conditions are shown in Table 2. Three requirements were set as the safe deployment conditions to determine either the success or the danger of the aerial deployment motion. The margin was defined as follows:

(13) $$\begin{equation} {\rm{Margin}}\;{\rm{ = }}\;{\rm{min}}\;\left( {{\rm{min}}\left( {{\rm{Value}}} \right) - {\rm{Lower}}\;{\rm{limit,}}\;{\rm{Higher}}\;{\rm{limit}}\;{\rm{-}}\;{\rm{max}}\left( {{\rm{Value}}} \right)} \right) \end{equation}$$

Table 2 Safe deployment condition for analysis of robustness

Value is a time-series data of the target parameter. A minimum or a maximum value in the simulation was used. The margin must be positive. The higher the margin is, the safer the deployment. Therefore, LSL and USL for the margin were 0 and nothing, respectively. Meanwhile, the failure of the aerial deployment motion was defined as the wings not having been fully deployed at the end of the simulation, or the right wing tip contacting the left wing.

The limit load factor requirements n _limit were set to 5 and –3. The upper requirement was designed higher than conventional ones^{(Reference Takeuchi, Yonemoto, Iwata, Narumi, Matsumoto, Fukuda, Uchida and Kato20)}. Since the angle-of-attack of the deploying aircraft took on not only positive but also negative values, the lower limit of the load factor was also important. Therefore, the lower limit of the load factor was set to –3 to withstand the aerodynamic force at a negative angle-of-attack.

A hinge reaction moment was regulated to prevent the structure near the hinge from breaking. The strength of the structure was expressed based on a limit load factor n _limit. The required hinge reaction moment M _{Reac Req} was written as the following equation:

(14) $$\begin{equation} {M_{{\rm{Reac}}\,{\rm{Req}}}} = \frac{{{n_{{\rm{limit}}}}{m_{{\rm{total}}}}gb_{{\rm{out}}}^2}}{{2b}}, \end{equation}$$

where m _total, g, and b represent the total mass, acceleration due to gravity and span, respectively. It was assumed that the span of each rigid body was the same; therefore b _out = b/3. The sign of the moment was positive for the right hinge and negative for the left hinge. Reference Fujita, Motoda and NagaiReference 7 revealed that the most sensitive condition was the hinge reaction moment. Therefore, the structure was reinforced and the limitation of the hinge reaction moment was changed to the value based on a load factor of 8. Here, the additional mass of the reinforcement was ignored.

The main wing may be damaged by the bending moment caused by the deployment of the right and left wings. Therefore, the bending moment at the centre of the main wing span M_b/2 due to the Z-direction force on the centre wing F_{Z Aw} and the hinge moment M _hinge was determined using Equation (15).

(15) $$\begin{equation} \begin{array}{@{}l@{}} {\rm{for}}\,{M_{{\rm{hinge}},\,B}} \cdot {M_{{\rm{hinge}},\,C}} < 0\\[3pt] \,\,\,\,\,\,{M_{b/2}} = {\textstyle{{} \over {}}}\frac{{{F_{ZAw}}{b_A}}}{8} + {\mathop{\rm sgn}} \left( {{M_{{\rm{hinge}},\,B}}} \right)\min \left( {\left| {{M_{{\rm{hinge}},\,B}}} \right|,\left| {{M_{{\rm{hinge}},\,C}}} \right|} \right),\\ {\rm{for}}\,{M_{{\rm{hinge}},\,B}} \cdot {M_{{\rm{hinge}},\,C}} \ge 0\\[3pt] \,\,\,\,\,\,{M_{b/2}} = {\textstyle{{} \over {}}}\frac{{{F_{ZAw}}{b_A}}}{8}, \end{array} \end{equation}$$

where b _A is the span of the centre wing. The first term represents the bending moment due to the aerodynamic force. Here, a uniformly distributed aerodynamic force was assumed. The second term represents the bending moment due to the hinge moments. The bending moment varies depending on the combination of the moments on the right and left hinges. When the directions of the hinge moments are the same, the bending moment due to hinge is zero, and the body rotates according to the sum of the moments. On the other hand, when the directions of the hinge moments are opposite, a bending moment due to the hinge occurs, and the body rotates according to the difference of the moments. The magnitude of this bending moment due to the hinge moments is equal to the lowest hinge moment. Here, the limitation was set as a bending moment at a limit load factor n _limit.

5.3 Nominal condition

The specifications of the rigid bodies are shown in Table 3. The specifications were defined based on those of the Mars aircraft^{(Reference Nagai and Oyama21)}. The masses of the right and left wings were both set to 0.16 kg. The mass of the centre was obtained as the difference between the mass of the Mars aircraft and the sum of the right and left wings. Inertia matrices of the right and left wings are calculated as a 0.050 m thick and uniform-density flat plate. The X-coordinate of the reference point of the inertia calculation was set to the centre of gravity of the whole aircraft. The Y and Z coordinates of the reference point were set to the centre of the flat plate. Since there are plane symmetries with respect to the XY plane and the XZ plane, the product of inertia of the right and left wings around the reference point was zero. An inertia matrix of the centre body was calculated from the values of the Mars aircraft and the right and left wings using the parallel axis theorem.

Table 3 Specifications of rigid bodies

#The origin of the coordinates for CG position expression was the point on the leading edge, the upper surface, the centre in spanwise of the main wing.

Table 4 shows the nominal conditions. g and T_SL are the acceleration due to gravity and left wing deployment delay time, respectively; h is the height; U _A, V _A and W _A are the translational velocities; φ _A, θ _A and ψ _A are the attitude angles; P _A, Q _A and R _A are the angular rates; and θ_dep and (dθ_dep/dt) are the deployment angle and its speed. This condition means that the wing-folded aircraft drops aiming its nose downwards with a velocity of 65 m/s. Therefore, the initial flight-path angle was –90°. The left wing started its deployment when the simulation began. The right wing started its deployment after the right-wing deployment delay time T_SR had passed. The simulation time was limited to the value of the right-wing deployment delay time T_SR plus 2 s. This limitation is useful for the reduction of the total calculation cost. The control surfaces were fixed to check the effect of the input conditions on the behaviour.

Table 4 Nominal conditions

5.4 Dispersive input variables

A scope of the variation was defined considering the aerial deployment condition on Mars. The important dispersive input variables are as follows:

• • Initial drop velocity

• • Surrounding gust velocity

• • Initial pitch angle

• • Initial height

First, let us focus on the initial drop velocity. The current mission has a deceleration phase using a parachute before the deployment. Therefore, the initial drop velocity can be estimated from the terminal velocity of the parachuting condition. The terminal velocity V _terminal can be obtained from the balance of the gravity and drag.

(16) $$\begin{equation} {V_{{\rm{terminal}}}} = \sqrt {\frac{{2{m_{{\rm{total}}}}g}}{{\rho {S_{{\rm{parachute}}}}{C_{D\,{\rm{parachute}}}}}}} , \end{equation}$$

where ρ, S _parachute, and C _Dparachute are the atmospheric density and the reference area and the drag coefficient of the parachute, respectively. Since the average of the terminal velocity can be designed through the parachute, the value was set to 65 m/s based on the cruising velocity. The standard deviation was considered next. In Equation (16), the density and the drag coefficient can vary from the nominal condition. In this study, the standard deviation of the initial drop velocity was assumed to be 10% of the average velocity.

The surrounding gust in the vertical direction was then ignored because it was almost the same as the variation of the drop velocity. Therefore, only the surrounding gust in the horizontal direction was considered. Based on Reference Rafkin and MichaelsRef. 6, the average and the standard deviation of the surrounding gust velocity were defined as 30 m/s and 20% of the average, respectively. The direction of the surrounding gust in relation to the aircraft was unpredictable. In this study, the direction of the surrounding gust was set to the direction from the right-wing side to the left-wing side in the initial condition because the direction showed worst result.

A third dispersive input variable is the initial pitch angle. An initial flight-path angle was fixed to –90°, i.e. dropping condition. Here, we considered the situation that the initial pitch angle varies while keeping the flight-path angle. Therefore, the initial angle-of-attack of the centre body was the same as the initial pitch angle. The average of the pitch angle at the beginning of the deployment varies depending on the trim angle in the folded state, the condition of the control surfaces, etc. The standard deviation of the pitch angle varies depending on the surrounding wind, the unsteady aerodynamic force for folded-state aircraft, the attitude control using control surfaces, the natural frequency of the aircraft pitching motion, etc. Therefore, the average and the standard deviation of the initial pitch angle were assumed to be 0° and 10°, respectively. It should be noted that an attitude with the aircraft nose aiming toward the ground and the upper surface aiming toward the north had been set to the origin of the Euler angle of attitude to avoid the singularity.

The final dispersive input variable is the initial height. The change of the height has an impact on the atmospheric density. The average initial height was set to 8 km using the ARES^{(Reference Guynn, Croom, Smith, Parks and Gelhausen2)} as reference. The standard deviation of the initial height was defined to be 20% of the average. The atmospheric density at the height of 8 km is approximately 0.008 kg/m³.

6.0 RESULTS AND DISCUSSION

The sigma levels of the margins for the safe deployment conditions for each variation of the aerial deployment condition are discussed first. After that, the optimal value, the successful region and the contributing factors of the total minimum sigma level are evaluated.

Figure 5 shows the sigma levels for the variation of the initial drop velocity. The legend “Success” means both wings were fully deployed and all safe deployment conditions were satisfied. The legend “Danger” means both wings were fully deployed, but at least one of the safe deployment conditions was violated. The legend “Failure” means that the right or left wing was not fully deployed or the right wing tip hit the left wing. The sigma level of the load factor showed relatively low values from (F, T_SR ) = (0.5, 0.3) to (F, T_SR ) = (1.5, 0.7), and from (F, T_SR ) = (0.8, 1.5) to (F, T_SR ) = (1.5, 1.8). This might be caused by the balance between the period of the motion and the right wing deployment delay time T_SR . The minimum sigma level n_min of more than 5 was obtained at large region of the design space.

Figure 5. Sigma levels for variation of initial drop velocity.

Figure 6 shows the sigma levels for the variation of the surrounding gust. The sigma level of the hinge reaction moment was lower than that of the initial drop velocity. This suggests that the wing deploys too rapidly owing to the surrounding gust. A high sigma level can be achieved at the design space where the deployment torque scale F is below 0.8, and the right-wing deployment delay time T_SR is from 0.2 to 1.2 s. Figure 5 indicates that the contributing factor for the minimum sigma level n_min for the variation of the surrounding gust is the hinge reaction moment.

Figure 6. Sigma levels for variation of surrounding gust.

Figure 7 shows the sigma levels for the variation of the initial pitch angle. The sigma level of the limit load factor is lower than that of the variations of the initial drop velocity and the surrounding gust. This reduction is caused by the change of the extremal value of the angle-of-attack due to the initial pitch angle change. The important result is that the lower limit of the successful deployment torque scale F increased to 0.5. This increase is the effect of the folding-direction aerodynamic force due to a bad initial pitch angle.

Figure 7. Sigma levels for variation of initial pitch angle.

Figure 8 shows the sigma levels for the variation of the initial height. Overall, the effect of the initial height is lower than other variations, and each sigma level has a relatively high value.

Figure 8. Sigma levels for variation of initial height.

Figure 9 shows the total minimum sigma level n _{min Total}, which is the minimum value among the minimum sigma level n _min of each variation. In addition, V _contributing and E _contributing, the dispersive input variables and the safe deployment conditions contributing to the total minimum sigma level n_{min Total} , are shown in Figs. 10 and 11, respectively. Globally, there are high sigma level regions at the deployment torque scale below 0.8, and the right-wing deployment delay time T_SR from 0.4 to 1.3 s. The result indicates that a low deployment-torque-scale F yields a high sigma level. However, the minimum value of the successful deployment torque scale F was approximately 0.5. This limitation depends on the aerial deployment performance in the variation of the initial pitch angle, as shown in Fig. 6. The change of the initial pitch angle caused a change in the angle-of-attack, which caused difficulties in the deployment of the wing. Therefore, the wings could not be deployed in the low deployment-torque design. In Fig. 9, the dispersive input variables contributing to the total minimum sigma level V _contributing, indicate that the variation of the surrounding gust has a high impact in general. On the other hand, in Fig. 10, the safe deployment conditions contributing to the total minimum sigma level E _contributing, indicate that the load factor and the hinge reaction moment have a high impact. As seen in the comparison of Fig. 9 to Fig. 10, the contributing factors at the high-sigma-level region are the condition of the load factor due to the initial pitch angle, and the safe deployment condition of the hinge reaction moment due to the surrounding gust. The results described above quantitatively clarified the effect of the deployment torque scale F and the right-wing deployment delay time T_SR on the robustness of the safe aerial deployment under a variety of conditions.

Figure 9. Total minimum sigma level n _{min Total}.

Figure 10. Dispersive input variables contributing to the total minimum sigma level V _contributing.

Figure 11. Safe deployment conditions contributing to the total minimum sigma level E _contributing.

The design point of the deployment torque scale F is 0.7, and the right-wing deployment delay time T_SR of 1.0 s was selected as a low deployment-torque-scale F and high sigma level n design. The total minimum sigma level n _{min Total} of this design point is more than 5. The sigma level around the design point is also relatively high. In addition, the selected design point has a margin of the deployment torque scale F of 0.1 to the failure condition of 0.5. The selected design point is expected to deploy safely and surely while keeping a low deployment-torque-scale F.