Aerodynamic Analysis of Camber Morphing Airfoils in Transition via Computational Fluid Dynamics

In this paper, the authors analyze an important but overlooked area, the aerodynamics of the variable camber morphing wing in transition, where 6% camber changes from 2% to 8% using the two airfoil configurations: NACA2410 and NACA8410. Many morphing works focus on analyzing the aerodynamics of a particular airfoil geometry or already morphed case. The authors mainly address "transitional" or "in-between" aerodynamics to understand the semantics of morphing in-flight and explore the linearity in the relationship when the camber rate is gradually changed. In general, morphing technologies are considered a new paradigm for next-generation aircraft designs with highly agile flight and control and a multidisciplinary optimal design process that enables aircraft to perform substantially better than current ones. Morphing aircraft adjust wing shapes conformally, promoting an enlarged flight envelope, enhanced performance, and higher energy sustainability. Whereas the recent advancement in manufacturing and material processing, composite and Smart materials has enabled the implementation of morphing wings, designing a morphing wing aircraft is more challenging than modern aircraft in terms of reliable numerical modeling and aerodynamic analysis. Hence, it is interesting to investigate modeling the transitional aerodynamics of morphing airfoils using a numerical analysis such as computational fluid dynamics. The result shows that the SST k-ω model with transition/curvature correction computes a reasonably accurate value than an analytical solution. Additionally, the CL is less sensitive to transition near the leading edge in airfoils. Therefore, as the camber rate changes or gradually increases, the aerodynamic behavior correspondingly changes linearly.


Introduction
The term morphing in aerospace refers to technologies that possibly enhance an aircraft's performance via alternating wing geometry corresponding to optimized wing shapes [1]. Recent airplane designs adapt fixed wings that could be said to be optimized in a single flight mode at a given task but far from aerodynamically optimized in others [2]. Meanwhile, adaptive or morphing wings could be more energy efficient [1,[3][4][5]. Biomimetics is an interdisciplinary field where materials and mechanisms resemble or mimic biological systems processes [6]. From biomimetic aspects, avian morphology inspires us in morphing wing aircraft designs and efficient flight control accordingly [5]. Morphing wings enable aircraft to optimize flight mode by changing wing shapes through active adjustment of wing shapes corresponding to external or internal flight conditions while flight [1,3,[7][8][9][10]. Most works in morphing focus on the mechanisms concerning conventional mechanical design or Smart materials and their proof-of-concept stages and material/structure-based internal mechanisms, which are far from realization and test flight of morphing aircraft .
The study of fluid flow around a morphing wing is essential to understanding aircraft performance and energy efficiency [46]. The aerodynamic loadings on the wing rely on the flow pattern. Also, it is coupled with the current shape of the wing [47][48][49][50][51]. A lifting-line theory could be used for well-defined airfoil shapes to analytically estimate and determine the L and D of the entire wing along the span direction [52]. Other standard techniques to study aerodynamics around the airfoils are experimental and numerical [53]. The significant advantage of experiments or testing is its accuracy; however, expensive in time and cost. For the case of a morphing wing, it would require replicating and controlling the behavior of a morphing wing inside an experimental apparatus such as a wind tunnel.
On the other hand, the numerical analysis offers flexibility in investigating various morphing geometries, providing the software's accurate boundary conditions and settings [54]. This article's related research works have numerical and experimental investigations [13,17,28,39,44]. However, it is necessary to eventually develop highly accurate and computationally inexpensive modeling and analysis tools to address morphing in transition states during flight to overcome the shortcomings of typical and conventional approaches and procedures [55]. As the first but most crucial step upward, authors investigate comparatively the numerical analysis of variable camber morphing airfoils and their behaviors in a transitional mode. In this present study, the authors have employed a CFD technique to establish the settings and conditions of airfoil models and adopted the National Advisory Committee for Aeronautics (NACA) [56]. First, the NACA2410 and NACA8410 are selected and analyzed where the entire wing has a taper ratio of one or a rectangular shape. Then, the results are compared with analytical results to validate them. Then, the NACA2410 on one side and the NACA8410 on the other side are modeled and studied to estimate their aerodynamics' transitional and linear/nonlinear behaviors morphing airfoils.

Methodology
ANSYS FLUENT [54] has been adopted to analyze the aerodynamics of morphing airfoils of interest. There are a lot of numerical methods that can be used to study the aerodynamics of an airfoil within ANSYS FLUENT. However, as the simulated flow is laminar to turbulent transitional flow, a RANS turbulence model is used. The objective is to compute the associated stresses, RSM, Nonlinear eddy viscosity models, and Linear eddy viscosity models [46,56]. RSM, also known as an RST model, is a higher-level and elaborate turbulence model [57]. The eddy viscosity approach [58] was neglected, and the Reynolds stresses were computed directly. The exact Reynolds stress transport equation was related to the directional effects of the Reynolds stress fields. The nonlinear eddy viscosity models are a class of turbulence models, and an eddy viscosity coefficient relates the mean turbulence field to the mean velocity field in their nonlinearity. Nevertheless, as obtained from the RANS equation, the Reynolds stresses in this study are assumed in a linear constitutive relation with the mean flow straining field, which is what linear eddy viscosity models were designed.
Among various subcategories in linear eddy-viscosity models, the SST k−ω turbulence model-one of the best airfoil aerodynamics studies for laminar to turbulent transitional flow-was adopted. This model also predicts the adverse pressure gradients and separates flow. Furthermore, the model includes two additional equations to represent the turbulent flow properties that could account for historical effects, such as the convection and diffusion in turbulent energy. The turbulent kinetic energy, k, that determines the energy in the turbulent flows is the first variable. The second one is ω the specific turbulent dissipation that determines the scale of the turbulence.
The second-order upwind scheme is performed to discretize all spatial terms in conservation equations. The solver handles the pressure-velocity coupled algorithm, and the least-squares cell-based setting calculates the gradient term. Authors first solve three momentum equations and pressure-correction continuity equations, then calculate scalar values such as temperature and turbulence quantities. The residual criteria below 10 −5 for . Far-field pressure: 70 kPa( )

Geometry Properties
Some of geometrical explanation for the airfoil is shown in Figure 1. 01 below. Wingspan: 9. 89 m Chord length: 0. 54 m Re: 1,000,000 where = , , , and represent the density, free stream velocity, chord length, and molecular viscosity, respectively. So, the flow speed from the previous equation will be 24.68 m/s.

Geometry
An asymmetric airfoil NACA2410 for the baseline and NACA8410 for entirely morphed camber were used. In addition, the authors used MATLAB to generate the geometric profiles of airfoils for analysis in ANSYS Design Modeler. Figure 2 shows the basic geometry (NACA2410) of a morphing wing structure without any morphing. However, Figure 3 shows an entirely morphed on one end, no morphing on the other end NACA8410 and NACA2410. ρ, V ∞ , c, and µ represent the density, free stream velocity, chord length, and molecular viscosity, respectively. So, the flow speed from the previous equation will be 24.68 m/s.

Geometry
An asymmetric airfoil NACA2410 for the baseline and NACA8410 for entirely morphed camber were used. In addition, the authors used MATLAB to generate the geometric profiles of airfoils for analysis in the ANSYS Design Modeler. Figure 2 shows the basic geometry (NACA2410) of a morphing wing structure without any morphing. However, Figure 3 shows an entirely morphed on one end, no morphing on the other end NACA8410 and NACA2410.

NACA Four-Digit Airfoil Specification
This NACA airfoil is characterized by four-digit NACA MPXX, which determines the camber, maximum camber position, and airfoil thickness, respectively. M is the number of the maximum camber divided by 100, and P is the number of the position of the maximum camber divided by 10. The remaining two digits, denoted XX, indicate the thickness divided by 100. For instance, NACA2412, where M = 2, where camber is 0.02 or 2% of the chord, P = 4, the maximum camber is positioned at 0.4 or 40% of the chord, and XX = 12, the thickness is 0.12 or 12% of the chord [56]. The equations for camber geometry generation are tabulated in Table 1.
The equation gives the thickness distribution:

NACA Four-Digit Airfoil Specification
This NACA airfoil is characterized by four-digit NACA MPXX, which determines the camber, maximum camber position, and airfoil thickness. M is the number of the maximum camber divided by 100, and P is the number of the position of the maximum camber divided by 10. The remaining two digits denoted XX, indicate the thickness divided by 100. For instance, NACA2412, where M = 2, where camber is 0.02 or 2% of the chord, P = 4, the maximum camber is positioned at 0.4 or 40% of the chord, and XX = 12, the thickness is 0.12 or 12% of the chord [56]. The equations for camber geometry generation are tabulated in Table 1. Table 1. Camber rate and gradient of geometry for airfoil description.
The equation gives the thickness distribution: where a 0 = 0.2969, a 1 = −0.126, a 2 = −0.3516, a 3 = 0.2843, a 4 = −0.1015, or −0.1035 for the closed trailing edge. The constants a 0 to a 4 are a 20% thick airfoil and T 0.2 adjusts the constants to the desired thickness. At the trailing edge where x = 1, a finite thickness of 0.0021 chord width for a 20% airfoil exists. When a closed trailing edge is required, the a 4 could be adjusted correspondingly. The y t is a half thickness and applies to both sides in the camber line. Given x, it is straight-forward to calculate y c and the gradient and the thickness. The upper and lower surface position is perpendicular to the camber line.
The most efficient way to plot the airfoil is to iterate through equally spaced values in x and the upper and lower surface coordinates. The points are more widely spaced near the leading edge. The flatter sections can be seen on the plots. A cosine is used with uniform increments of β to group the points near the airfoil's ends.
The thickness description curves for airfoil description are summarized in Table 2. Table 2. Thickness description curves for airfoil description.
Upper Surface Lower Surface

Meshing
The far-domain is eight chord lengths away from the airfoil boundary, and a Design Modeler [54] adopted a C-type mesh. Then, a refined mesh is generated and applied next to the wall to capture details of y + close to 1. Last, the distance between a wall in the airfoil and the first layer, y is calculated as below: wall shear stress τ w counted by: A literature survey refers to a formula for Cf, skin friction on a plate thus: y value can be calculated by taking Re, µ, v air and ρ values and substituting them in previous equations.
A Mesh Independence study was conducted as well. Richardson Extrapolation was adopted to calculate the grid convergence index, and it was confirmed that an index of less than 3% was obtained for the various grids tested. Computational time also plays a role in deciding the mesh element size. The grid comprising 7,060,448 elements (Element size: 0.01 m) was the most refined mesh tested, but as shown in Figure 4, varying the element size and thus the number of cells shows only minor changes in the values of C D we obtained from the simulations. Therefore, the element size was adjusted to 0.05 m for the rest of the study.

Boundary Conditions
There are four boundary conditions imposed: (1) The inlet equals a velocity inlet, (2) The outlet equals a pressure outlet,  Computational models applied boundary conditions are shown in Figure 5a-e, as below. = = 1.89 × 10 mm A Mesh Independence study was conducted as well. Richardson Extrapolation was adopted to calculate the grid convergence index, and it was confirmed that an index of less than 3% was obtained for the various grids tested. Computational time also plays a role in deciding the mesh element size. The grid comprising 7,060,448 elements (Element size: 0.01 m), was the most refined mesh tested, but as shown in Figure 4, varying the element size and thus the number of cells shows only minor changes in the values of obtained from the simulations. Therefore, the element size was adjusted to 0.05 m for the rest of the study.

Result
The k-ε (k-epsilon) model with the EWT (enhanced wall treatment), SST (shear stress transport) k-ω (k-omega), and SST k-ω with transition and curvature correction have been performed and analyzed. All models use a NACA2410 at 8° of AoA to verify the best model. The result shows the comparative study of and with analytical solutions as shown in Table 3. As a result, SST k-ω with IT/CC (intermittency transition/curvature

Result
The k-ε (k-epsilon) model with the EWT (enhanced wall treatment), SST (shear stress transport) k-ω (k-omega), and SST k-ω with transition and curvature correction have been performed and analyzed. All models use a NACA2410 at 8 • of AoA to verify the best model. The result shows the comparative study of C L and C D with analytical solutions as shown in Table 3. As a result, SST k-ω with IT/CC (intermittency transition/curvature correction) performs best in comparison. The result's deviation in the analytical solution comes because the model presumes turbulent flow around the airfoil starts from the front of the leading edge. However, the flow around the airfoil surface is laminar in practice and similar to flow over a flat plate. Then, the flow changes to turbulent downstream of the airfoil. Since conventional turbulent models induce the turbulent C D for the entire airfoil, the simulation is larger than the experiment or analytics. Table 3. Turbulence model effects on C L and C D . On the other hand, the SST k-ω with IT/CC outperforms in perspective. A detailed description is found in [59]. The TBL carries more energy and the C D is more significant than in the VBL. Meanwhile, the C L is less sensitive to the laminar and turbulence transition. Thus, the C L is higher in its value than experimental or analytical solutions [38].

NACA2410
The SST k-ω IT/CC is adopted to calculate the C L and C D where the NACA2410 3D airfoil with AoA changes from 0 • to 8 • , as shown in Table 4. Both C D and C L are compared with the analytical solutions, and the errors are below 13.8%. The C D and C L are graphed while varying AoA as shown in Figure 6. The C L obtained from ANSYS FLUENT match the analytical solution's most overlapping values. However, the C D from ANSYS is more significant than the one from the analytical solution. It is noted that the C D deviates from the analytical solution induced from the model setting and near-wall mesh resolution. Another reason for the deviation is that the analytical solution assumes the flow is laminar, but in reality, it is turbulent when Re = 1 million. So, the      Figure 9 shows the contour of the turbulent kinetic energy, k for AoA of attacks of 1 • and 8 • . It is noted that turbulent energy is generated near the trailing edge and downstream. Thus, the mesh is recommended to be sufficiently refined near the downstream region. Figure 9 shows the contour of the turbulent kinetic energy for AoA of attacks of 1° and 8°. It is noted that turbulent energy is generated near the trailing edge and downstream. Thus, the mesh is recommended to be sufficiently refined near the downstream region.  The pressure distribution is shown in Figure 10 for AoA 1 • , 5 • , and 8 • , respectively. As the AoA increases, it is noted that the pressure distribution at the top and bottom walls become separated from each. On average, both the C D and C L increase along with the AoA. The pressure distribution is shown in Figure 10 for AoA 1°, 5°, and 8°, respectively. As the AoA increases, it is noted that the pressure distribution at the top and bottom walls become separated from each. On average, both the and increase along with the AoA.

NACA8410
The simulations have been performed for the NACA8410 or 8% camber morphing case. The MATLAB generates a new mesh by applying the NACA8410 profiles at both ends. The same model setup and boundary conditions are applied for the NACA2410. The result of and for NACA8410 3D airfoil for the AoA varying from 0° to 8° is shown in Table 5. The and of NACA8410 and NACA2410 are graphed while varying AoA as shown in Figure 11. It is noted that both configurations have a similar trend.

NACA8410
The simulations have been performed for the NACA8410 or 8% camber morphing case. The MATLAB generates a new mesh by applying the NACA8410 profiles at both ends. The same model setup and boundary conditions are applied for the NACA2410. The result of C D and C L for NACA8410 3D airfoil for the AoA varying from 0 • to 8 • is shown in Table 5. The C D and C L of NACA8410 and NACA2410 are graphed while varying AoA as shown in Figure 11. It is noted that both configurations have a similar trend.  Figures 12 and 13 show the pressure and velocity distribution of the NACA8410. Compared with Figures 7 and 8 for the NACA2410, more significant pressure, and velocity variations are found, as shown in Figures 12 and 13. The distribution stretches up to farther distances, similar to the and plotted in Figure 11. Increasing the camber rate leads to more significant in and as well.  Figures 12 and 13 show the pressure and velocity distribution of the NACA8410. Compared with Figures 7 and 8 for the NACA2410, more significant pressure, and velocity variations are found, as shown in Figures 12 and 13. The distribution stretches up to farther distances, similar to the C D and C L plotted in Figure 11. Increasing the camber rate leads to more significant in L and D as well.
(a) (b) Figure 11. and of NACA2410 and NACA8410 vs. varying AoA. Figures 12 and 13 show the pressure and velocity distribution of the NACA8410. Compared with Figures 7 and 8 for the NACA2410, more significant pressure and velocity variations are found, as shown in Figures 12 and 13. The distribution stretches up to farther distances similar to the and plotted in Figure 11. Increasing the camber rate leads to more significant in and as well.  The effects of turbulent kinetic energy are also explored by generating turbulent kinetic energy contours for angles of attacks of 1° and 8°, as shown in Figure 14. It is also noted that turbulent energy is mainly dense near the trailing edge and downstream in low and high AoA. The effects of turbulent kinetic energy, k are also explored by generating turbulent kinetic energy, k contours for the AoA of 1 • and 8 • , as shown in Figure 14. It is also noted that turbulent energy is mainly dense near the trailing edge and downstream in low and high AoA. The effects of turbulent kinetic energy are also explored by generating turbulent kinetic energy contours for angles of attacks of 1° and 8°, as shown in Figure 14. It is also noted that turbulent energy is mainly dense near the trailing edge and downstream in low and high AoA. For a more direct and precise comparison, the pressure distribution plots around NACA2410 and NACA8410 airfoils are presented in Figure 15. The pressure difference between the top and the bottom surface of the NACA8410 shown in purple is more significant than the pressure differences of the NACA2410 in orange. That correlates to higher and for NACA8410 than for NACA2410. For a more direct and precise comparison, the pressure distribution plots around NACA2410 and NACA8410 airfoils are presented in Figure 15. For a more direct and precise comparison, the pressure distribution plots around NACA2410 and NACA8410 airfoils are presented in Figure 15. The pressure difference between the top and the bottom surface of the NACA8410 shown in purple is more significant than the pressure differences of the NACA2410 in orange. That correlates to higher and for NACA8410 than for NACA2410.

From NACA2410 to NACA8410 Transition
The NACA2410 on one side and the NACA8410 on the other for transitional morphing have been analyzed to simulate morphing wings in transition. As a result, the cross- The pressure difference between the top and the bottom surface of the NACA8410 shown in purple is more significant than the pressure differences of the NACA2410 in orange. That correlates to higher C D and C L for NACA8410 than for NACA2410.

From NACA2410 to NACA8410 Transition
The NACA2410 on one side and NACA8410 on the other have been analyzed to simulate morphing wings in transition. As a result, the cross-sectional geometry is linearly varied in between.
The computational methodology was compared against analytical data and validated. Similar to the previous step, the same boundary conditions and settings were also applied to this case. Simulations were run for NACA2410 to 8410 transition 3D morphing wing geometry over a range of AoA from 0 • to 8 • . The C D and C L results are summarized in Table 6. Table 6. NACA2410-8410 end simulated results.

Model
AoA  Figures 16-19 show the NACA2410-NACA8410 transitional morphing scenario in their multi-axial input/output characteristics on deformation and stress. sectional geometry is linearly varied in between. The computational methodology was compared against analytical data and validated. Similar to the previous step, the same boundary conditions and settings were also applied to this case. Simulations were run for NACA2410 to 8410 transition 3D morphing wing geometry over a range of AoA from 0° to 8°. The and results are summarized in Table 6. Figures 16-19 show the NACA2410-NACA8410 transitional morphing scenario in their multi-axial input/output characteristics on deformation and stress.   sectional geometry is linearly varied in between. The computational methodology was compared against analytical data and validated. Similar to the previous step, the same boundary conditions and settings were also applied to this case. Simulations were run for NACA2410 to 8410 transition 3D morphing wing geometry over a range of AoA from 0° to 8°. The and results are summarized in Table 6. Figures 16-19 show the NACA2410-NACA8410 transitional morphing scenario in their multi-axial input/output characteristics on deformation and stress.    The C D and C L vs. the AoA plots of NACA8410, NACA2410, and NACA2410 to 8410 transition cases are shown in Figure 16. As shown in Figure 20, the plot for the C D and C L of the transitional case (NACA2410 to NACA8410) lies in between the L plots for NACA2410 and NACA8410 airfoils. Similar trends are seen in the plots for the C D . Whereas the C L curve for the transitional case lies in the middle of the NACA2410 and NACA8410 C L curves, the C D curve for the transitional case lies closer to NACA8410 C D curve, which could be interpreted as the C L being linear overall in the wing spar direction, the C D is more sensitive to structural deformation from coupling. In addition, since the overall trend is maintained over the wing span direction, C L and C D can be interpolated in-between. The pressure and velocity distribution of the transitional case with NACA8410 configuration on one side and NACA2410 configuration on the other side are shown in Figures 21 and 22. A slightly more significant pressure difference is observed for the NACA8410 configuration end than for the NACA2410 configuration end, making sense as an 8% cambered airfoil configuration should generate a more significant pressure difference between the top and the bottom surfaces.       Figure 23 shows the contour of the turbulent kinetic energy at 5° AoA. It is noted that the turbulent energy is highly dense near the trailing edge and downstream.  Figure 23 shows the contour of the turbulent kinetic energy, k at 5 • AoA. Again, it is noted that the turbulent energy is highly dense near the trailing edge and downstream.

Summary
This paper presents an aerodynamic analysis of variable camber morphing wings via the CFD model with ANSYS FLUENT. The main focus lies on the transitional camber rate to resemble the camber morphing during flight. A refined mesh near the airfoil wall simulates the flow details. These suggested models perform based on the pre-defined initial and boundary conditions around a morphed wing. The results of the NACA2410, while varying AoA, are benchmarked and compared with analytical ones. The SST k-ω turbulence model with transition and curvature correction features has been adopted. It is noted that the C L is less sensitive to the flow transition near the airfoil's leading edge. The C D and C L of NACA8410 airfoil are also calculated using the benchmarked ANSYS FLUENT setup. Lastly, a variable camber wing with end configurations matching the geometric profile of the NACA2410 airfoil on one side and the NACA8410 airfoil on the other side were analyzed to simulate the behavior of a camber morphing wing in transition. The C D and C L obtained from simulations are studied and analyzed. Furthermore, some essential flow parameters are graphed and analyzed for all cases, including that the flow characteristics change far from the airfoil wall and downstream. Therefore, it is suggested that fine mesh be considered to capture all the details. The significant contribution of this study is to explore and suggest a methodological approach to understand an aerodynamic analysis of morphing wings in flight more accurately.

Conclusions
A study of a numerical model and analysis of the aerodynamics of a variable camber morphing wing is presented in this paper. Contrary to a conventional wing, where the shape of the airfoil or wing is fixed, and the aerodynamic parameters can be computed for the fixed geometry, a morphing wing changes shape during the flight, which makes the task of computing the aerodynamic parameters of a shape-changing wing in transition challenging. Morphing Wing design concepts proposed morphing the wing geometry to a specific configuration at one end and another at the other end. However, aerodynamic analysis computing the C L and C D for such a morphing wing transition is not readily available. Hence, the authors in this paper present methodologies for transitional aerodynamics in a morphing wing by adopting a geometry; NACA2410 on one side; the NACA8410 on the other. At the same time, the intermediate geometry linearly transforms in the middle of the airfoil. The SST k-ω turbulence model with intermittency transition/curvature correction was used to compute the C L and C D for the morphing condition. The computational model used is validated against the analytical model, and an assuring benchmark was obtained with an error of less than 13.8% for the C D and 0.5% for C L . Plots of C L and C D obtained for NACA2410, NACA8410, and NACA2410 to 8410 transition cases are compared. The results obtained for the morphing in transition case show that the variable camber morphing wing still performs better performance than the conventional airfoil. This study also indicates that morphing wings is a promising alternative to achieving a high-efficiency level on traditional aircraft. More importantly, it presents a method to analyze the aerodynamic performance of a scenario of a morphing wing in a transitional mode where a gradually changing morphing rate has been applied to a single body wing along the spanwise direction. Acknowledgments: Authors appreciate Muhammad Akbar and Anwar Ansulami at Tennessee State University, Nashville, TN, and Jacobs Somnic at SUNY Korea for their sincere efforts and technical advice.

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