Robust Design Optimization of Supersonic Biplane Airfoil Using Efficient Uncertainty Analysis Method for Discontinuous Problem

Abstract

Busemann’s supersonic biplane airfoil can reduce wave drag through shock interactions at its designed freestream Mach number. However, a choking phenomenon occurs with a decrease in the freestream Mach number, and the drag coefficient increases significantly, resulting in an aerodynamic problem with a discontinuous change in the performance function. In this study, an uncertainty analysis method, the divided inexpensive Monte Carlo simulation (IMCS), is proposed to solve discontinuous problems efficiently and is applied to Busemann’s biplane airfoil. In the divided IMCS, the discontinuity point is determined using a simple sampling method. The uncertainty input space is divided at the detected discontinuity point, and a surrogate model is constructed for each space. Uncertainty analysis was performed using the constructed surrogate models, and the results of the divided IMCS showed qualitative agreement with those of the conventional Monte Carlo simulation, which is the most straightforward uncertainty analysis method. Moreover, the divided IMCS significantly reduced the computational cost of the uncertainty analysis. A robust design optimization of the supersonic biplane airfoil was performed using the divided IMCS, yielding more robust designs than Busemann’s biplane airfoil. The usefulness of the divided IMCS for uncertainty analysis of discontinuous problems was confirmed.

1. Introduction

Several institutions and companies worldwide are engaged in the development of the next-generation supersonic transport (SST). Lockheed Martin X-59 Quiet Supersonic Technology is an SST model that aims to achieve low-boom flight. Boom Technology’s Boom Overture is also being developed to achieve low-drag and low-boom SST. Busemann’s supersonic biplane airfoil configuration is well known [1,2,3,4,5,6,7,8]. It comprises two triangle airfoils and can reduce wave drag by shock interactions at the design’s freestream Mach number (in this study, M_∞ = 1.7). However, a choking phenomenon occurs with a decrease in the freestream Mach number (for example, M_∞ = 1.6), and the drag coefficient (C_d) increases significantly, which results in an aerodynamic problem with a discontinuous change in the performance function (referred to as a discontinuous problem in this paper). Figure 1 shows pressure visualizations at M_∞ = 1.6 and 1.7. This discontinuity makes it difficult to accurately estimate the uncertainties in the aerodynamic performance with respect to uncertain inputs.

Considerable attention has been paid to uncertainty analysis [4,5,7,8,9,10,11,12,13,14,15,16], particularly in the field of computational fluid dynamics (CFD). Conventional numerical analyses deterministically evaluate performance with respect to defined input conditions. However, several sources include uncertainties in the input conditions, such as manufacturing tolerances and fluctuations in the operating conditions. Therefore, it is important to perform an uncertainty analysis to evaluate the mean, standard deviation, and probability density function (PDF) of the output functions with respect to uncertain inputs. Robust design optimization (RDO) can be performed by applying uncertainty analysis to design optimization [4,7,10,11,17,18,19,20,21,22,23]. Several studies have been conducted on RDO. Pisaroni et al. presented a novel approach for the robust optimization of aerodynamic shapes based on a combination of single- and multi-objective evolutionary algorithms and a continuation multilevel Monte Carlo methodology [20]. Zhao et al. provided a comprehensive review of robust aerodynamic design optimization methodologies, including uncertainty modeling, uncertainty quantification, and robust optimization [22]. Schaefer et al. presented a gradient-enhanced robust design strategy using spatially accurate polynomial chaos that was applied to the NASA X-59 QueSST aircraft to increase the robustness of its sonic boom performance [23].

The uncertainty analysis of Busemann’s biplane airfoil was also investigated. Tabata et al. optimized the 2D cross-sectional shape of Busemann’s biplane airfoil by applying the uncertainty analysis method of a divided difference filter (DDF) [7]. Optimal designs were obtained with a higher mean and smaller standard deviation of the lift-to-drag ratio (L/D) than those of Busemann’s biplane airfoil. However, the DDF is not a highly accurate method for discontinuous problems because the discontinuous change in the performance function of Busemann’s biplane airfoil cannot be captured. In [8], Kawai et al. proposed edge detection-based multielement polynomial chaos expansion (EDME-PCE), in which the uncertain space is adaptively decomposed by edge detection [24,25]. EDME-PCE can capture the discontinuity of the choking phenomenon; therefore, it is effective for the discontinuous problem of Busemann’s biplane airfoil. However, its computational cost is approximately 14-times higher than that of deterministic CFD analysis, and it is expensive to apply to RDO. In [26], Liem et al. proposed several methods to improve the accuracy of surrogate models for aerodynamic performance prediction (but not for Busemann’s biplane airfoil). The mixture-of-experts approach divides the input space into subspaces to build local surrogate models. In [27], Bettebghor et al. developed a surrogate modeling approximation approach using a mixture of experts for discontinuous problems. The surrogate models constructed using this approach were able to express discontinuities well; however, this approach required several high-fidelity computations to construct the surrogate models, which resulted in high computational costs.

This research aims to develop a highly accurate uncertainty analysis method with a low computational cost for discontinuous problems. The choking phenomenon of Busemann’s supersonic biplane airfoil was investigated as a discontinuous problem. For a highly accurate uncertainty analysis of discontinuous problems, detecting the discontinuity position with high accuracy is necessary. In this study, an uncertainty analysis method, divided inexpensive Monte Carlo simulation (IMCS), is proposed for efficiently solving the discontinuous problem. In the divided IMCS, the discontinuity position is detected using a simple sampling method, and the uncertainty input space is divided at the detected discontinuity position. The proposed sampling method for detecting the discontinuity position has the advantage of being applicable not only to Busemann’s biplane airfoil problem but also to other discontinuous problems, such as other SSTs or other engineering fields. Another advantage of the proposed sampling method is that uncertainty analysis of discontinuous problems can be performed at low computational cost. To confirm the usefulness of the proposed method, a one-dimensional uncertainty analysis problem of Busemann’s biplane airfoil was considered. The input uncertainty is included in M_∞ as a normal distribution. The results of the divided IMCS were compared with those of MCS, DDF, and IMCS. RDO of the 2D cross-sectional shape of the supersonic biplane airfoil was performed using the divided IMCS to obtain more robust airfoil shapes.

The remainder of this paper is organized as follows. In Section 2, the uncertainty analysis methods, including the proposed method, are described. The numerical schemes are described in Section 3. In Section 4, the proposed method is compared with other uncertainty analysis methods, and the characteristics and advantages of the proposed method are described. Section 5 presents the results of applying the proposed method to RDO. Finally, Section 6 presents the conclusions of this study. To demonstrate its extensibility to practical problems with multiple input uncertainties, Appendix A presents an example of applying the proposed method to a two-dimensional uncertainty analysis problem.

2. Uncertainty Analysis Methods

2.1. MCS (Monte Carlo Simulation)

The most straightforward approach for evaluating uncertainty is MCS; however, the number of samples required for statistical calculations is large, which results in a high computational cost. In the MCS, the mean value μ of a function f with respect to uncertain input variables $\overrightarrow{r}$ is calculated by the following equation:

$$\mu \left(f_{\overrightarrow{r}}\right)=\int \cdots \int \left(f_{\left(\overrightarrow{r}\right)}{\phi }_{\left(\overrightarrow{r}\right)}\right)d{r}_1\cdots d{r}_n$$

(1)

where n is the number of uncertain inputs, and all uncertain inputs are assumed to be independent. Equation (1) shows the general equation of the MCS; however, in this study, only one uncertain input was considered. φ is the PDF of $\overrightarrow{r}$. When samples ${\overrightarrow{r}}_i$ are generated based on the distribution of φ, the mean μ and standard deviation σ of f are predicted using the following equations:

$$\mu \left(f_{\overrightarrow{r}}\right)\cong \frac{1}{N}\sum _{i=1}^N{f}_{\left({\overrightarrow{r}}_i\right)}$$

(2)

$${\sigma }^2\left(f_{\overrightarrow{r}}\right)\cong \frac{1}{N}\sum _{i=1}^N{\left[f_{\left({\overrightarrow{r}}_i\right)}-\mu \left(f_{\overrightarrow{r}}\right)\right]}^2$$

(3)

where N is the number of samples for the MCS. A Latin hypercube sampling (LHS) method is utilized for generating the samples in this study.

2.2. DDF (Divided Difference Filter)

In DDF, μ and σ of f are predicted as follows [7,9,10,11]

$$\mu \left(f_{\left(\overrightarrow{r}\right)}\right)=G_0{f}_{\left({\overrightarrow{r}}_0\right)}+\sum _{i=1}^n{G}_i\left(f_{\left({\overrightarrow{r}}_{i+}\right)}+f_{\left({\overrightarrow{r}}_{i-}\right)}\right)$$

(4)

$${\sigma }^2\left(f_{\left(\overrightarrow{r}\right)}\right)=\frac{1}{2}\sum _{i=1}^n\left[G_i{\left(f_{\left({\overrightarrow{r}}_{i+}\right)}-f_{\left({\overrightarrow{r}}_{i-}\right)}\right)}^2+\left(G_i-2G_i^2\right){\left(f_{\left({\overrightarrow{r}}_{i+}\right)}+f_{\left({\overrightarrow{r}}_{i-}\right)}-2f_{\left({\overrightarrow{r}}_0\right)}\right)}^2\right]$$

(5)

where n is the number of uncertain inputs, and all uncertain inputs are assumed to be independent. Equations (4) and (5) show the general equations of the DDF; however, in this study, only one uncertain input was considered. ${\overrightarrow{r}}_0=\left({\mu }_{r_1}, \cdots , {\mu }_{r_i}, \cdots ,{\mu }_{r_n}\right)$ and ${{\overrightarrow{r}}_i}_{\pm }=\left({\mu }_{r_1}, \cdots , {\mu }_{r_i}\pm ∆r_i, \cdots ,{\mu }_{r_n}\right)$ are referred to as sigma points. The weight coefficients $G_0$, $G_i$, and $∆r_i$ are defined by the following equations:

$$G_0=1-\sum _{i=1}^n\frac{1}{K_{r_i}}$$

(6)

$$G_i=\frac{1}{2K_{r_i}}$$

(7)

$$∆r_i={\sigma }_{r_i}\sqrt{K_{r_i}}$$

(8)

where $K_{r_i}$ and ${\sigma }_{r_i}$ are, respectively, the kurtosis and standard deviation of the i-th uncertain input variable. The PDF of the i-th uncertain input variable determines $K_{r_i}$. For normal distributions, $K_{r_i}$ is set to 3. The PDF of f cannot be evaluated. The computational cost of DDF is significantly lower than that of MCS, because the number of evaluated samples is only 2n + 1.

2.3. IMCS (Inexpensive Monte Carlo Simulation)

IMCS replaces a large number of exact evaluations of MCS with functional estimations on a surrogate model. Therefore, IMCS can reduce the computational cost of uncertainty analysis [12,13,14,15,16]. Because the accuracy of IMCS depends on the accuracy of the surrogate model, it tends to be less accurate for discontinuous problems. In this study, a surrogate model was constructed using ordinary Kriging (OK).

2.4. Divided IMCS

The divided IMCS was proposed to solve the uncertainty analysis of discontinuous problems efficiently. In this study, to confirm the usefulness of the proposed method for the discontinuous problem, a one-dimensional uncertainty analysis problem for Busemann’s biplane airfoil was investigated. The one-dimensional uncertainty analysis problem can be explained as follows. Figure 2 shows a flowchart of the divided IMCS. In the divided IMCS, the discontinuity point is detected using a simple sampling method, and the uncertainty input space is divided at the detected discontinuity point. The simple sampling method for detecting discontinuity points has the advantage of being easily applicable, not only to Busemann’s biplane airfoil problem but also to other discontinuous problems. In the sampling method, two initial samples selected at both ends of the uncertainty input space were first evaluated. In this study, the freestream Mach number M_∞ is the uncertain input. The difference in the performance function (C_d) of the two samples is calculated, and if the difference is higher than a threshold value ΔC_d (in this study, ΔC_d = 0.05), a discontinuity exists between the two samples. Figure 3 shows an example of adding samples to detect the discontinuity point. An additional sample was defined at the center of the two samples and evaluated. The discontinuity range is defined as the distance between the two samples in which a discontinuity exists. Additional samples are iteratively generated until the discontinuity range becomes smaller than a threshold value ασ, where α is a user-defined coefficient, and σ is the standard deviation of the uncertain input. When the iterative calculation is completed, the center of the discontinuity range is set as the final discontinuity point. The uncertainty input space was divided into two clusters at the discontinuity point. If a cluster does not have three samples, additional samples are selected at outward locations of the initial sample (at intervals of ασ). Additional samples are evaluated until each cluster contains three samples. A surrogate model of OK is constructed in each divided uncertainty input space. Then, each sample point for the IMCS is evaluated using the surrogate model constructed for the region to which the sample point for the IMCS belongs. In this study, OK was utilized as a standard surrogate model; however, arbitrary surrogate models can be adopted based on the needs or preferences of users. The samples evaluated for the detection of the discontinuity point were also used to construct surrogate models; therefore, the proposed method has the potential to reduce computational cost.

3. Numerical Schemes

Figure 4 shows the Busemann’s biplane airfoil configuration. The thickness ratio of the upper and lower airfoils was 0.05, and the distance ratio between the upper and lower airfoils was 0.5. The angle of attack was set to 2.0 degrees. The 2D compressible Euler equations were solved using a gridless method [28,29]. Computational points were distributed referring to the Mach angle calculated from the condition of M_∞. Therefore, the distribution of the computational points was changed with the variation in M_∞, as shown in Figure 5. Approximately 16,000 computational points were used in the study. In [6], the design optimization of a supersonic biplane airfoil was discussed, and a convergence study for the computational point resolution was performed. The resolution of the computational points in this study was selected by considering the computational load. The numerical flux at the intermediate position of the two computational points was calculated using Roe’s approximate Riemann solver, which has second-order spatial accuracy. Temporal discretization was performed using a four-stage Runge–Kutta scheme.

4. Comparison of Uncertainty Analysis Methods

In this section, the proposed divided IMCS is compared with the MCS, DDF, and IMCS. The input uncertainty is provided on M_∞ as a normal distribution. The mean μ of M_∞ is set to 1.70, and the standard deviation σ is set to 0.05. Therefore, the sigma points for the DDF are defined at M_∞ of 1.613, 1.70, and 1.787 from Equation (8). In the divided IMCS, three cases are investigated with the user-defined coefficient α of 0.1, 0.5, and 1.0. The two initial samples for the divided IMCS are defined at M_∞ of 1.55 and 1.85, which means μ ± 3σ of M_∞. In the conventional IMCS, a surrogate model is constructed over the uncertainty input space using the sample points obtained in the divided IMCS (α = 1.0). IMCS was performed using 1000 samples for both the conventional IMCS and divided IMCS. For comparison, MCS was also performed using 1000 samples.

Figure 6 shows the surrogate models of C_d and L/D obtained by the conventional IMCS and divided IMCS (α = 0.1, 0.5, and 1.0). The exact response and PDF of the uncertain input M_∞ are also shown in Figure 6. The surrogate model of the conventional IMCS showed an oscillatory response and did not show qualitative agreement with the exact response. By contrast, the surrogate models of the divided IMCS show qualitative agreement with the exact response and can successfully capture the discontinuous behavior. It can be confirmed that the accuracy of the surrogate model is improved by appropriately dividing the uncertainty input space.

Table 1. Comparison of mean and standard deviation (std) of C_d, L/D, and computational cost.

Method	Cd				L/D				Computational Cost 1
	Mean		Std		Mean		Std
MCS	1.46 × 10−2	0%	2.45 × 10−2	0%	15.3	0%	4.52	0%	1000
DDF	2.10 × 10−2	44%	3.22 × 10−2	31%	14.4	−6%	5.89	30%	3
IMCS	3.09 × 10−2	112%	2.98 × 10−2	22%	12.4	−19%	5.18	15%	6
divided IMCS (α = 0.1)	1.45 × 10−2	−1%	2.44 × 10−2	0%	15.3	0%	4.48	−1%	8
divided IMCS (α = 0.5)	1.55 × 10−2	6%	2.50 × 10−2	2%	15.1	−1%	4.56	1%	7
divided IMCS (α = 1.0)	1.79 × 10−2	22%	2.89 × 10−2	18%	14.1	−8%	5.00	11%	6

¹ Cost of one deterministic CFD analysis is considered as the unit.

shows a comparison of the mean and standard deviation of C_d and L/D, as well as the computational cost. The computational cost was defined by considering the cost of one deterministic CFD analysis as a unit. Computational costs other than those for CFD analyses (for example, construction of surrogate models) were not considered because such computational costs are significantly lower than those of CFD analyses. The numbers of samples for the divided IMCS at α = 0.1, 0.5, and 1.0 are, respectively, eight, seven, and six. Although the statistics provided by the DDF include the influence of the choking phenomenon as the choking phenomenon occurs at the sigma point at M_∞ = 1.613, the accuracy of the statistics is lower than that of the divided IMCS. The statistics provided by the conventional IMCS did not show qualitative agreement with those of the MCS because of the oscillatory response of the surrogate models. The statistics provided by the divided IMCS show qualitative agreement with those of MCS, demonstrating the effectiveness of the divided IMCS detecting the discontinuity point. The divided IMCS algorithm estimates the discontinuity point as the center of two samples between which a discontinuity exists. Therefore, the estimated discontinuity point includes a maximum error of ασ/2, which is one of the limitations of the divided IMCS. Smaller values of α mean that more samples are required, and the computational cost increases. However, the discontinuity point is detected more accurately. There is a trade-off relationship between the computational cost and accuracy of uncertainty analysis, and users can define the coefficient α depending on whether low computational cost or high accuracy is more important. Cases with smaller values of α show better agreement with the MCS. Even in the case of the largest α (α = 1.0), significantly better results than those with the conventional IMCS are obtained. Therefore, the user-defined coefficient α = 1.0 is utilized for RDO in the next section. The computational cost of the MCS was 1000, whereas that of the divided IMCS was 6–8. Therefore, the divided IMCS was 125-times more efficient than the MCS. It was confirmed that the divided IMCS is useful for uncertainty analysis problems including discontinuous functional variations.

5. Robust Design Optimization (RDO) Problem

5.1. Definition of the Optimization Problem

The 2D cross-sectional shape of Busemann’s biplane airfoil was optimized using the uncertainty analysis method of the divided IMCS to confirm its usefulness in RDO. The input uncertainty is provided on M_∞ as a normal distribution. The mean and standard deviation of the uncertain input M_∞ are set to 1.7 and 0.05, respectively. From Section 4, it is confirmed that cases with smaller values of α show better agreement with the MCS. Even in the case of the largest α (α = 1.0), significantly better results than that of the conventional IMCS were obtained. Therefore, the user-defined coefficient α was set to 1.0 for the RDO. The other parameter settings were the same as those described in the previous section. Because the objective of optimization is to design high-performance airfoil shapes with high robustness, the optimization problem was set to maximize the mean of L/D and minimize the standard deviation of L/D, as shown in Equation (9), to discuss the trade-off relationship. The geometrical constraints were also defined in Equation (9) and given for the total cross-sectional area and cross-sectional area of the upper airfoil to prevent the design of very thin airfoils. In Equation (9), S_Tnew, S_Torig, S_Unew, and S_Uorig, respectively, indicate the total sectional area of a newly designed airfoil, total sectional area of Busemann’s biplane airfoil, sectional area of a newly designed upper airfoil, and sectional area of the upper airfoil of Busemann’s biplane. Since it is known from Licher’s biplane airfoil configuration that the thickness of the upper airfoil tends to be reduced to improve aerodynamic performance [30], the cross-sectional area of the upper airfoil was set to be larger than half of S_Uorig.

Figure 7 shows a flowchart of the RDO. For RDO, the Kriging response surface approach [31,32,33] and genetic algorithm (GA) were utilized. We used an in-house real-coded multi-objective GA developed by referring to [34,35]. First, uncertainty analyses for the 30 initial samples were performed using the divided IMCS. Among the 30 initial samples, 1 was the original Busemann’s biplane airfoil, and the other 29 were generated using the LHS method. In the design variables space, surrogate models for μ of L/D, σ of L/D, total cross-sectional area, and cross-sectional area of the upper airfoil were constructed by OK. The GA was utilized to search for global optimal solutions in the surrogate models. Additional samples were searched using expected improvement (EI) [36] and expected hypervolume improvement (EHVI) [37]. The surrogate models were updated using the results of the uncertainty analyses of the additional samples. This process was repeated until the total number of samples reached 200. The airfoil shapes were represented by Bezier curves using 13 control points and 15 design variables. Figure 8 shows an example of a deformed airfoil and the range of each design variable.

$$\mathrm{R}\mathrm{o}\mathrm{b}\mathrm{u}\mathrm{s}\mathrm{t} \mathrm{D}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n} \mathrm{O}\mathrm{p}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\hspace{1em}\hspace{1em}\left\{\begin{array}{c}\begin{array}{c}maximize {\mu }_{L/D}\\ minimize {\sigma }_{L/D}\end{array}\\ subject to \begin{array}{c}{S}_{T_{new}} \ge S_{T_{orig}}\\ S_{U_{new}} \ge \frac{S_{U_{orig}}}{2}\end{array}\end{array}\right.$$

(9)

$$\mathrm{D}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{c} \mathrm{D}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n} \mathrm{O}\mathrm{p}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n} \left\{\begin{array}{c}maximize L/D\\ subject to \begin{array}{c}{S}_{T_{new}} \ge S_{T_{orig}}\\ S_{U_{new}} \ge \frac{S_{U_{orig}}}{2}\end{array}\end{array}\right.$$

(10)

5.2. Results and Discussion

Figure 9 shows the performance of all samples obtained in the RDO. In Figure 9, two types of additional samples are shown: one for samples in which the choking phenomenon occurs and discontinuity exists (Figure 2, route A), and the other for samples in which the choking phenomenon does not occur and discontinuity does not exist (Figure 2, route B). The results for samples with the choking phenomenon showed that the standard deviation of L/D was larger than that for samples without the choking phenomenon. This was because the choking phenomenon caused a massive increase in C_d (that is, a massive decrease in L/D), which increased the standard deviation. There were several samples with better μ and σ of L/D than those of the original Busemann’s biplane airfoil. This indicates that the designs explored using RDO are more robust than the original Busemann’s biplane airfoil.

From the results of the RDO, five characteristic optimal designs were selected and referred to as designs A–E. Design A has the smallest σ of L/D, and design E has the largest μ of L/D. Designs B, C, and D exhibited intermediate performances between designs A and E. Figure 10 shows the surrogate models of L/D obtained by the divided IMCS. In the optimal designs A, B, and C, the choking phenomenon was not detected. In the divided IMCS, the choking phenomenon is detected from the two initial samples defined at M_∞ of 1.55 and 1.85. Therefore, the choking phenomenon is not detected when it occurs in the range of M_∞ < 1.55. The choking phenomenon occurred at M_∞ = 1.634 with Busemann’s biplane airfoil, and at M_∞ = 1.569 with the optimal designs D and E. Designs D and E have the same discontinuity point, because the divided IMCS algorithm estimates the discontinuity point as the center of the two samples between which a discontinuity exists. The limitation of the divided IMCS is that the estimated discontinuity point includes a maximum error of ασ/2. In designs D and E, the discontinuity point is located at a lower M_∞ where the PDF of M_∞ is sufficiently small so that the negative impact of the choking phenomenon on the robustness is also sufficiently small. Figure 11 shows the PDFs of L/D obtained using the divided IMCS. It can be confirmed that the PDF distributions shift to higher L/D from designs A to E. In the PDF of Busemann’s biplane airfoil, a peak attributed to the choking phenomenon can be observed in the area of L/D < 5. With respect to the representative designs A–E, the peak attributed to the choking phenomenon cannot be observed or is negligibly small. For comparison, deterministic design optimization was also performed, as defined in Equation (10). In this optimization, M_∞ was set to 1.7, and the fluctuation in M_∞ was not considered. The objective of deterministic design optimization was to maximize L/D. The geometrical constraints and design variables of the airfoil shapes were identical to those of the RDO. The optimal design obtained by deterministic design optimization is called the deterministic optimal. The exact off-design performance of the optimal designs is investigated by additional CFD evaluations with intervals of 0.01 in M_∞, and the results are shown in Figure 12. The choking phenomenon occurred with the designs A, B, and C in the region of M_∞ < 1.55. The discontinuity points of each optimal design moved towards larger M_∞ from designs A to E, and the deterministic optimal design had the right-most discontinuity point. The value of L/D at M_∞ of 1.7 also increased from designs A to E, and the deterministic optimal design had the largest L/D at M_∞ of 1.7. The deterministic optimal had the highest L/D at the mean condition of uncertain input (M_∞ = 1.7), while it had the poorest performance considering robustness. Additional IMCSs (referred to as IMCS’) were performed using the results of Figure 12, in which the surrogate models were defined by the linear interpolation of the CFD results shown in Figure 12.

Table 2. Comparison of mean and standard deviation (std) of L/D.

Airfoil Shape	Divided IMCS		MCS		IMCS’
	Mean of L/D	Std of L/D	Mean of L/D	Std of L/D	Mean of L/D	Std of L/D
Busemann’s biplane	14.1	5.00	15.3	4.52	15.2	4.59
Design A	15.0	0.0841	-	-	14.9	0.132
Design B	19.5	0.131	-	-	19.5	0.151
Design C	20.5	0.521	-	-	21.6	0.772
Design D	21.7	1.42	-	-	21.8	1.26
Design E	22.1	1.84	-	-	22.2	1.80
Deterministic Optimal	21.0	4.95	-	-	21.3	4.34

summarizes the mean and standard deviation values of L/D obtained by divided IMCS, MCS, and IMCS’. In Figure 9b, the results of the MCS are indicated for Busemann’s biplane airfoil, and the results of the IMCS’ are indicated for representative designs. Table 2 and Figure 9b indicate that the results of the optimal designs obtained by the divided IMCS agree well with the results of the IMCS’, regardless of whether the choking phenomenon is captured. The divided IMCS achieved reasonable uncertainty analysis, with a lower computational cost of 3–7 compared with the computational cost of 41–46 for IMCS’. The deterministic optimal had a high mean value of L/D and a large standard deviation value of L/D. Therefore, in deterministic design optimization, the highest performance can be realized under the design conditions; however, the robustness is low. By contrast, it can be confirmed that airfoils with high robust performance can be designed inexpensively by the RDO using the divided IMCS.

Figure 13 shows the 2D cross-sectional shapes of the optimal designs. The upper airfoils of the optimal designs were thinner than the lower airfoils, and the central vertices of the lower airfoil were rounded. The trend of thinner upper airfoils was the same as that of Licher’s biplane airfoil configuration [30], which was proposed as an efficient biplane airfoil configuration under lifted conditions. Figure 14, Figure 15, Figure 16, Figure 17 and Figure 18, respectively, show pressure distributions around Busemann’s biplane airfoil, deterministic optimal design, and optimal designs A, C, and E at M_∞ of 1.65, 1.70, and 1.75. In Busemann’s biplane airfoil at M_∞ = 1.7, the shock wave from the leading edge of the upper airfoil hits the upstream side of the center vertex of the lower airfoil, and the static pressure increases significantly, as shown in Figure 14b. At smaller M_∞ (Figure 14a), the shock wave hits further upstream of the center vertex of the lower airfoil. In the deterministic optimal design, in an off-design condition of M_∞ = 1.65, the shock wave also hits the upstream side of the center vertex of the lower airfoil, and the static pressure increases significantly, as shown in Figure 15a. Considering the optimal designs A, C, and E, the shock wave from the leading edge of the upper airfoil hits the rounded upper surface of the lower airfoil, and the distribution of the static pressure is not changed significantly, even with the variation in M_∞, as shown in Figure 16, Figure 17 and Figure 18. In the optimal designs A, C, and E, the upper airfoils are thinner, and the wedge angles of the leading edge are smaller. The angles of the shock wave from the leading edge of the upper airfoil decrease, and the shock wave hits the downstream side of the lower airfoil, decreasing the drag force. Figure 19 shows the relationship between the performance values obtained by IMCS’ and the minimum distance between the upper and lower airfoils. In Figure 19, the minimum distance is normalized to that of the original Busemann’s biplane airfoil, as shown in Figure 4. Additional uncertainty analyses with variations in the minimum distance of Busemann’s biplane airfoil were also performed for comparison. From Figure 19, it can be confirmed that the distances of the optimal designs are larger than those of the original Busemann’s biplane airfoil, and the distance increases in the following order for the deterministic optimal designs, E, D, C, B, and A. The mean values of Busemann’s biplane airfoil configurations were approximately constant with variations in the distance. Although the standard deviation values of Busemann’s biplane airfoil configurations decreased as the distance increased, the standard deviation values of the optimal designs obtained by RDO were significantly smaller. This means that the robustness of the optimal airfoils increased not only with the increase in the distance between the airfoils but also with the appropriate geometrical deformation of the airfoils, such as the rounded upper surface of the lower airfoil.

6. Conclusions

An efficient uncertainty analysis method, divided IMCS, was proposed for discontinuous problems. The choking phenomenon of the Busemann’s biplane airfoil was investigated as a discontinuous problem. To confirm the effectiveness of the proposed method, a one-dimensional uncertain input was considered. In the divided IMCS, the discontinuity point was detected using a sampling method, and the uncertainty input space was divided at the detected discontinuity point. A surrogate model of OK was constructed in each divided uncertainty input space, and each sample point for the IMCS was evaluated using the surrogate model generated for the region to which the sample point belonged. An uncertainty analysis of the supersonic flow around the 2D sectional shape of Busemann’s biplane airfoil was performed using the MCS, DDF, conventional IMCS, and divided IMCS. The input uncertainty was provided on M_∞ as a normal distribution. The user-defined coefficients for the divided IMCS α were set to 0.1, 0.5, and 1.0. The surrogate models of the divided IMCS showed qualitative agreement with the exact response and could successfully capture discontinuous behavior, whereas the surrogate model of the conventional IMCS did not show qualitative agreement with the exact response. In the divided IMCS, the estimated discontinuity point includes a maximum error of ασ/2, and this is a limitation of the divided IMCS. Smaller values of α required more samples, while the discontinuity point was more accurately detected. From the viewpoint of computational cost, the divided IMCS is significantly more efficient than MCS.

Subsequently, RDO for the 2D sectional shape of Busemann’s biplane airfoil was performed using the divided IMCS. The objectives of the RDO were maximizing μ of L/D and minimizing σ of L/D. More robust optimal designs than the original Busemann’s biplane airfoil were successfully explored using the RDO. A deterministic design optimization was also performed for comparison, in which L/D at M_∞ =1.7 was maximized. The exact off-design performance of the optimal designs was investigated using additional CFD evaluations, and additional IMCSs (referred to as IMCS’) were performed using the results of the exact off-design performance. The results obtained using the divided IMCS were in good agreement with those obtained using IMCS’. Leveraging deterministic design optimization, the highest performance could be realized under the design conditions, whereas the robustness was significantly lower. However, airfoils with highly robust performance can be designed inexpensively by RDO using the divided IMCS. In the robust optimal designs, the upper airfoil was thinner than the lower airfoil, and the center vertex of the lower airfoil was rounded to achieve highly robust performance.

In conclusion, the divided IMCS can perform highly accurate uncertainty analysis for discontinuous problems with a low computational cost and is useful for realizing practical methods for solving complex uncertainty analysis problems including discontinuities. Since only a one-dimensional uncertainty analysis problem was considered in this study, additional algorithm modifications are required for practical problems with multiple input uncertainties. Appendix A presents an example of applying the divided IMCS to a two-dimensional uncertainty analysis problem. In future studies, the divided IMCS will be extended to more practical problems with three- (or more) dimensional uncertainty input spaces. It can also be applied to more realistic RDO problems such as three-dimensional Busemann’s biplane wing configurations. Discontinuous problems have also been reported for structural optimization problems that consider buckling [27] and transonic airfoil shape optimization problems [38]. The divided IMCS method is also beneficial for such engineering problems.