Numerical Investigation of Winglet Aerodynamics and Dimple Effect of NACA 0017 Airfoil for a Freight Aircraft

Abstract

Drag reduction is an ever-present challenge within the aeronautical engineering industry. This paper presents two substantial wing modifications: the addition of a winglet of a freighter aircraft and a dimpled wing on the NACA 0017 aerofoils. Studies on nine (9) different geometries of dimpled aerofoils were performed against a control model of an aerofoil without any dimple. Computational fluid dynamics (CFD) analysis was performed using two (2) commercial CFD platforms. This paper also explored two novel solutions of aircraft optimisation to mitigate the effects of drag and leading-edge pressure, while increasing the effect of lift. The optimised performance model of a freighter aircraft increased its aerodynamic efficiency. The study found that at take-off velocity of 82 m/s, winglets decreased pressure on the wing by 16.31%, through flow redirection and better flow integration into aerofoils wake. The study also analysed the separation layer and its effect through the appropriate use of the dimple effect. Increased lift effects were observed on a NACA 0017 aerofoil. Despite the low increase in drag of 6% from the modifications, the resultant L/D ratio was highly increased. This study also faced some challenges with validating the model. Hence some validation approaches were taken, and some other approaches suggested for future studies.

1. Introduction

Current trends in the research on dimples show aerodynamic effects on a variety of applications, such as golf balls, wind turbine blades, and aircraft wings (for military, cargo, passenger, or freighter aircrafts). Freighter aeroplanes serve public, government, and private sectors, delivering different goods and services across the continents. Air transportation is the quickest means of transport globally today. Aeroplanes, which are energy-powered flying vehicles whose wings are fixed, and the mass is greater than that of the displaced air [1] are not the only means of air transport. Others include helicopters, aircrafts, drones, stealth crafts, fighter jets, etc. These serve other purposes for logistics, military supplies, or distributing small, electronic components. Amazon, a global logistics company, has achieved current drone delivery applications for supplying goods [2]. Also, drones have been deployed for organ transplant in the University of Maryland Medical Center, USA [3,4] There are other applications of aviation crafts, such as military fighter jets, and stealth jets. Fast jets are typically required to respond to a combat situation as quickly and effectively as possible. Due to the application, they must be smaller with delicate parts to increase aerodynamic efficiency. They usually house a delta-shaped wing to increase agility in air-to-air combat [5], using an aerodynamically efficient aeroplane [6].

By size, the specification of the freighter aircraft results in a large cargo capacity, which is the most critical factor when considering their design. The Antonov An-225 Mriya’s wingspan of 88.4 metres [7] is nine times (9×) larger than the wingspan of the F-16 Fighting Falcon [8]. The Antonov’s wings are situated at the top of the fuselage to maximise the fuselage’s storage capacity. The wings slope to optimise the aircraft’s ability to glide; this consequently reduces fuel consumption. A smaller, fast jet has less exposed surface area and will experience less drag than a freighter aircraft when travelling at the same velocity and altitude [9]. One drag reduction method is to increase the aspect ratio of a wing. However, there are two immediate limitations on the length of an aircraft’s wingspan. One restriction is imposed by airports globally, for they do not have the room nor the runway width to adhere to the needs of an extended wing length. The second is imposed by the structural capability of the wing’s materials. Early aircrafts, such as the Wright Brothers Flyer, used wood as their primary construction material. Although this initially seems to be a primitive construct, an omen to the early days of material science, the use of wood to fabricate major aircraft components is present throughout the 1940s. This design feature is prominent in the De Havilland Mosquito, first conceived in November of 1940, flying throughout World War Two (WW2) [10]. For any full-size modern aircraft, aluminium is wood’s successor. While wood has a high strength to weight ratio, aluminium’s is considerably higher, without the long production process of drying and lamination. Aluminium alloys are also more resistant to erosion and more accessible to fabricate when compared to wood. Aluminium’s increased strength means it is structurally more capable than wood and can bare a larger wing load. In principle, aircraft wings can be described as cantilever beam; only supported at one end, the bending moment at the wing tip will be larger than the bending moment of the wing at the fuselage. The longer the wing, the larger the resulting bending moment and the more support the wing will need. Aluminium can provide the desired level of support, above most airport width restrictions. It, therefore, follows that adhering to airport demands is the primary restriction on aircrafts.

In regard to the effect of dimple on the aerofoil, different studies that include the dimple geometry, the angle of attack (AOA) and the inlet velocity have been considered by various researchers. The implementation of the dimple effect can be used to reduce the size of the vortex in the wake of the aerofoil can be reduced. The dimple effect is a drag reduction method, most notable for its appearance on golf balls. When the flow follows the impressions on the surface of the golf ball, the elevation changes induced by the dimple, thereby creating an extension to the boundary layer. As the flow progresses along with the golf ball, the low pressure in the boundary layer seeks equilibrium, effectively sucking in the surrounding high-pressure flow. This attempt to balance pressure allows the airflow surrounding the ball to travel closer to its surface and follow its shape for longer. The tighter airflow produces a reduction in the size of the vortex found in the wake of the ball [11]. However, the vortex phenomenon occurs in different fluid medium, both in air and water. An earlier study conducted in 1945 discovered that with a standard dimple depth of 0.254 mm, the ball experienced a 25% reduction in drag when compared to a smooth surfaced ball [12]. The author concluded that dimple depth also affected performance, with balls that house shallower impressions experiencing marginal or no aerodynamic improvement. Golf ball manufacturers have started to vary dimple geometry and patterns to optimise ball performance. Experiments have been conducted that compare the standard design to elliptical and hexagonal impressions.

One objective of this study is to determine the dimple effect phenomenon and to implement it by optimising aircraft wings. Despite little research on this, evidence from existing experiments on dimple spread, pattern, and the number of dimples present on the surface can vary drag results by up to 40% [13]. Arunkumar et al. [14] investigated dimples with a diameter of 0.2 c on the top of an aerofoil operating at 10° AOA. The study found that a semi-circular dimple geometry was optimum, increasing the lift efficiency coefficient by 0.73, and 30.8%. In another research by Soh and Al-Obai [15], simulations were conducted on a NACA 0012 aerofoil in which dimples of varying geometries were placed at 0.2 c, both on the top and the bottom of the aerofoil. When the results were compared against a smooth aerofoil with an inlet speed of 10 m/s, it was found that an elliptical dimple located on the bottom of the wing, having a depth of 0.03125 C increased aerodynamic efficiency by 6.42% at 5° AOA, and 3% at 10° AOA. This study concludes that dimples on the underside of the aerofoil are more effective than those placed on top. While both modifications resulted in aerodynamic improvement, the variation in results can be explained using Bernoulli’s principle, which states that ‘velocity is inversely proportional to pressure’. Thus, pressure underneath the aerofoil is more significant than that on top. It follows that a dimple located underneath the aerofoil will experience more pressure because of the surrounding flow. Consequently, the dimple has more potential to increase lift and decrease drag [15].

Different numerical methods have been utilized on aerofoils, such as turbulence SST k-omega [16,17]. Simulations have also been conducted on the effects of a dimple pattern on the fuselage [18]. In this study, dimples were positioned at the second half of the fuselage after the flow had become more laminar. Therefore, the turbulent flow that the dimples induced resulted in drag being increased by 30%. There is potential for this concept to be improved if the dimple pattern could be relocated to the plane’s nose. To become more effective, dimples can be placed on the separation point of the aerofoil [19]. A separation point occurs when the boundary layer develops further, and the eddies contained within start to experience flow reversal [20]. After the separation point, the size of the boundary layer begins to increase, eventually leading to the formation of a turbulent wake. Consequently, this directs the flow reversal point towards the trailing edge, decreases the wake’s size behind the aerofoil, and subsequently decreases drag. For most symmetrical aerofoils at 0° AOA (angle of attack), a separation point occurs towards the trailing edge of the aerofoil, as AOA increases the separation towards the leading edge [21,22]. Figure 1 shows an illustration depicting the constituent forces acting on dimpled aerofoils, including its Angle of attack (AOA).

This paper investigates the flow behaviour on a winglet and the dimple effect on aircraft wings to determine separation points that can be applied to improve aircraft performance. To understand the flow progression over an aerofoil, the formation of a boundary layer and the subsequent separation point over a NACA 00017 aerofoil at 0° AOA has been conducted with details on the interaction of airflow around aerofoil. Section 1 introduces the research, while Section 2 presents some fundamental theories on the investigation. The study results are presented in Section 3, and some discussions are made in Section 4. Concluding remarks are made in Section 5.

2. Numerical Model

2.1. The Fundamental Theory

2.1.1. Aerodynamic Forces

There are different forces that can influence the performance and aerodynamics of the aircrafts. These range from the viscous force, motor force, drag force, weight of the aircraft, wing geometry, turbulence, and other forces. However, the effect of drag force will be analysed in this study, as this has been found to affect the performance greatly. The drag force is as a result of the combined effects of pressure and wall shear forces. The drag force can be expressed mathematically by Equation (1), where, C_d is the drag coefficient, ρ is the density of air, V is the flow velocity relative to the aircraft, and A is the cross-sectional area of the wing. However, there is a need to increase the power of the commercial van to make it more efficient. To obtain the power (P) required to overcome this drag, we consider the expression in Equation (2) where the drag coefficient is proportional to the cube of the velocity (V).

$$FD=\frac{1}{2}{C}_d\rho V^{2}A$$

(1)

$$P=\frac{1}{2}\rho C_{d}A{V}^3$$

(2)

2.1.2. Air Flow Theories

The airflow around the freight aircrafts is high-speed aerodynamics; thus, the airflow field can be regarded as three-dimensional (3D) compressible flow. However, in the present study, the focus is on an aerofoil and a weightless winglet. The fundamental governing equations for the three-dimensional incompressible flow are presented in this section. The Navier–Stokes equation was applied with the K-$\omega$ or k-$\epsilon$ turbulence model to solve the compressible airflow. Navier–Strokes equations consist of governing equations, such as the Continuity equation (see Equation (3)) and Momentum equations (see Equations (4)–(6)), where u, v, and z are the components of velocity vector in x, y, and z-components; ρ is the air density, P₀ is static pressure, and S_m,x, S_m,y, and S_m,z are source terms for x, y, and z components. On the left hand side of Equations (3)–(6), the first term is a density gradient. Also, on the right hand side of Equations (3)–(6), the first term is a pressure gradient. It should be noted that the Navier–Stokes Equations (3)–(6) are not incompressible since it involves density variation. Hence, they are expressed herein since they are compressible but should have the required energy equation as in Equations (4)–(6). Thus, we have:

$$\frac{\partial \rho }{\partial t}+𝛻\left(pV\right)=0$$

(3)

$$\frac{\partial \rho \left(u{E}_T\right)}{\partial t}+𝛻\left(\rho uV\right)=-\frac{\partial p}{\partial x}+𝛻\left(\mu 𝛻u\right)+S_{m,x}$$

(4)

$$\frac{\partial \rho \left({\mathit{vE}}_T\right)}{\partial t}+𝛻\left(\rho vV\right)=-\frac{\partial p}{\partial y}+𝛻\left(\mu 𝛻v\right)+S_{m,y}$$

(5)

$$\frac{\partial \rho \left(w{E}_T\right)}{\partial t}+𝛻\left(\rho wV\right)=-\frac{\partial p}{\partial z}+𝛻\left(\mu 𝛻w\right)+S_{m,z}$$

(6)

2.1.3. Governing Equations for Aerofoils

Equations (7) and (8) show the normal and axial forces on top of the aerofoil, respectively, and Equations (9) and (10) display the respective forces underneath the aerofoil [22]. Figure 1 shows an illustration depicting the constituent forces acting on dimpled aerofoils.

$$dN{\prime }_u=\left(-p_{u}cos\theta -{\tau }_{u}sin\theta \right)d{s}_u$$

(7)

$$dA{\prime }_u=\left(-p_{u}sin\theta -{\tau }_{u}cos\theta \right)d{s}_u$$

(8)

$$dN{\prime }_l=\left(p_{l}cos\theta -{\tau }_{l}sin\theta \right)d{s}_l$$

(9)

$$dA{\prime }_l=\left(p_{l}sin\theta -{\tau }_{l}cos\theta \right)d{s}_l$$

(10)

2.2. Turbulence Modelling

The SST k-$\omega$ (shear stress transport) turbulence model [16,17] was used to model the turbulent flow with automatic wall function. There are several considerations to make when selecting an appropriate turbulence model, the primary affecting factor of which is flow type. The simulations throughout this study were solved using the SST k-$\omega$ model. This model type utilizes the functions found in two independent models, k-$\epsilon$ and k-$\omega$, and combines their calculations to obtain a more stable and therefore, reliable result, as presented in the Equations (11)–(15). As the intensity of the eddy currents at the inlet of two independent simulations can never be modelled identically, varying flows in the free stream change and effect overall results.

The SST k-$\omega$ model uses the accuracy of the k-$\epsilon$ to resolve values in the free stream where there is a constant flow type and use the k-$\omega$ model to solve near the wall of the object where Reynold number changes [16,17]. This is achieved by adding a blending function, $F_1$ [23]. A value of $F_1=1$ is given to the cells within the boundary layer. This will cancel $\left(1-F_1\right)$ and the equation will be solved using the k-$\omega$ model, producing accurate, near-wall results. If $F_1=0$, then $\left(1-F_1\right)=1$ indicating flow is in the free stream, and the governing equations for the k-$\epsilon$ equation will be utilized [23]. The SST k-$\omega$ model also includes a near all viscosity limiter that limits the freeflow behaviour of the fluid. This replicates effects found within the Boussinesq’s turbulence hypothesis, which states that the prominent eddies found within a highly turbulent fluid alter the way momentum is transferred, causing the fluid to behave as if it were vicious.

$$\frac{\partial \tilde{\epsilon }}{\partial t}=\frac{\partial }{\partial y}\left\{\left(v_{kv}+\frac{v_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \tilde{\epsilon }}{\partial y}\right\}+C_1{f}_1{v}_t\frac{\tilde{\epsilon }}{k}{\left(\frac{\partial u}{\partial y}\right)}^2-C_2{f}_2\frac{{\tilde{\epsilon }}^2}{k}+E$$

(11)

$$f_1=1$$

(12)

$$f_2=1-0.3exp\left(-R_{eT}{}^2\right)$$

(13)

$$\frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho u_{j}k\right)}{\partial x}=p-{\beta }^{*}\rho \omega k+\frac{\partial }{\partial x}\left[\left(\mu +{\sigma }_k\frac{\rho k}{\omega }\right)\frac{\partial k}{\partial x}\right]$$

(14)

$$\frac{\partial \omega }{\partial t}+u_j\frac{\partial \omega }{\partial x}=\alpha S^2-\beta {\omega }^2+\frac{\partial }{\partial x}\left[\left({\nu }_{kv}+{\sigma }_{\omega }{\nu }_t\right)\frac{\partial \omega }{\partial x}\right]+2\left(1-F_1\right){\sigma }_{\omega }{}^2\frac{1}{\omega }\frac{\partial k_t}{\partial x}\frac{\partial \omega }{\partial x}$$

(15)

2.3. Aerofoil Geometry

To accurately model this wing and its curvature, aerofoil cross-sections were taken from several points along the wing’s length, from the root aerofoil to the wingtip. Data points for each NACA 0017 aerofoil were downloaded from [24], where they were lofted together to create a CAD model of the wing shown in Figure 2.

The three distinct winglets in Figure 2a–c were then designed in comparison to the normal Boeing 787-9 Dreamliner wings in Figure 2d, with the hypothesis that the delta-shaped winglets yield the highest percentage in drag reduction. Details of the parameters are presented in

Table 1. Details of the Freighter aircraft.

Parameter	Value
Total wingspan	60 m
Root chord length	13.5 m
Speed range	1350 km
Take-off weight	254,011 kg
Take-off speed chosen	82 m/s

and

Table 2. Parameters for the enclosure for the aerofoil dimple designs in 2D.

Parameter	Value (mm)
Leading Edge to Front Boundary	15,000
Trailing Edge to Rear Boundary	30,000
Total Boundary Length	46,000
Leading Edge to Horizontal Boundary	15,000
Total Boundary Height	30,000
Aerofoil Length	1000
Maximum Aerofoil Height	170
Dimple Length for ELLIPSE Design	46.835
Max. Dimple Depression Depth ELLIPSE	11.94
Dimple Length for SQUARE Design	46.835
Max. Dimple Depression Depth SQUARE	13.962
Dimple Length for TRIANGLE Design	46.835
Max. Dimple Depression Depth TRIANGLE	9.879

. The wing geometry for this study was modelled on the Boeing 787-9 Dreamliner aircraft and designed in Autodesk Inventor Professional 2019 CAD software. The Boeing 787-9 Dreamliner is a long-haul aircraft [25], with a total wingspan of 60 m, a root chord length of 13.5 m, a range of 14,010 km [26], and a take-off weight of 254,011 kg [27], as presented in Table 1.

The geometrical designs for the dimples and their location along the aerofoil are found in Figure 3. In turn each dimple geometry (ellipse, triangle, square) was placed atop, underneath and both on top and underneath the aerofoil. The dimensions of the dimples are shown in Table 2.

Aerofoil geometry is primarily governed by the structural requirements of the entire wing [28]. Thicker aerofoils are typically found on larger aircraft as they can accommodate a larger spar. A larger spar increases resistance to bending stress, something that is more prominent on longer wings, as shown in Equation (16), as in the literature [29]. Inherently, structural support is more beneficial for aerodynamic performance despite the elevated level of drag larger aerofoils create [28]. As a secondary point, aircraft fuel tanks are housed in the wing, and a larger spar has a larger area, so it can accommodate more fuel as needed for lengthy flights. Adhering to these parameters, the thicker NACA 0017 aerofoil was selected for this part of the study. The dimples were selected to be placed at 0.92 c, given that this is approximately where the separation point over the NACA 0017 aerofoil is located. This location was experimentally obtained and validated with literature from Anyoji, et al. [21].

$$M_{maximum}=EI\frac{d^{2}x}{d{x}^2}$$

(16)

2.4. Methodology of CFD Analysis

The CFD model of different geometries, as discussed in Section 2.2, was created using CAD software. The model was modified to suit the geometry for the CFD model by importing the model to the Autodesk CFD 2019 software. Some mesh refinements in the CFD model were performed to define the geometry due to the dimple type for the NACA 0017 aerofoil. Figure 3 shows the arrangements along with each dimple geometry. The modifications of small sizes, like edges, gaps, and shape were performed, for the three design concepts-triangular, square, and elliptical (round) shapes. 50 time-consuming simulations were conducted to obtain a stable value for static pressure. Regarding this study, three recognisable dimple geometries were cut out of the aerofoil: an ellipse (round), a square and a triangle. A common depression depth of 0.021 c and a width of 0.032 c was selected. This value was determined by considering the studies previously mentioned and a need to preserve the structural integrity of the aerofoil. Each dimple shape is located underneath the aerofoil and then compared to the effects of a dimple located on the top, and a dimple located on both the top and underneath of the aerofoil. This design concept was conducted to evaluate the combined effects of the dimples. Figure 4 represents the setup for the aerodynamic properties in SimScale platform, showing the reference axis.

2.5. Boundary Domain

For this study, two different models were considered: a 3D model for the winglet and a 2D model for the dimple effect. Boundary domains near the front of the aerofoil will produce unreliable results due to unstructured or still-developing flow traversing the aerofoil. A boundary domain that is too close to the trailing edge will not accurately capture the dispersion of turbulent air. For the CFD model, the geometry must include an inlet radius that fits into the CFD software, such as ANSYS Fluent, OpenFOAM, or Autodesk CFD 2019. Thus, the C-shaped domain can be exchanged for the square domain; however, the square domain was sufficient on this study.

For this investigation, the CFD model settings, as applied in Autodesk CFD, are presented in

Table 3. Parameters for the Initial Boundary Conditions.

Parameters	Values
Inlet Velocity, $({\mathit{v}}_{\mathit{t}\mathit{a}\mathit{k}\mathit{e}\mathbf{-}\mathit{o}\mathit{f}\mathit{f}}$)	82 m/s
Air Temperature, T	19.85 °C
Density of Air, ${\mathit{\rho }}_{\mathit{a}\mathit{i}\mathit{r}}$	1.225 Kg/m3
Viscosity of Air, ${\mathit{\mu }}_{\mathit{a}\mathit{i}\mathit{r}}$	1.81 × 10−5 Kg/m
Gas Constant	287.05 m2/s2–k
Compressibility Constant	1.4
Inlet Pressure (Gauge)	0 Pa
Outlet Pressure (Gauge)	0 Pa
Mesh Size	1 mm sub-domain (and auto)
Mesh Type	Tetrahedral
Flow Type	Incompressible (steady state)
Fluid Type	Inviscid
Turbulent Model	K-omega SST
Side Walls	Slip/Symmetry
Reynolds Number (2D)	5.55 × 106
Reynolds Number (3D)	76.314 × 106

. Figure 5 shows the boundary domain used in the 3D model, while Figure 6 illustrates the boundary conditions for the 2D model. This shows the main settings for the CFD simulation that was run to ensure a suitable boundary domain size. First, there is no reversed flow along the normal of the upper or lower boundary domain. Second, no pressure gradient is displayed along any of the five boundaries, not including the sixth boundary where the wing is walled. These factors indicate that the free-stream flows in these areas are unobstructed. This proves that the length of the suggested boundary domain is large enough to allow for turbulent air to reach a state of homogeneity.

2.6. Mesh Independence

Research into mesh quality highlighted the importance of obtaining results independent of the mesh size. A mesh size of 90 mm was deemed appropriate for the simulation in the 3D model. Refinement regions were created in

Table 4. Mesh Refinement Regions for the 2D Model.

Region Type	Box Region
X Offset	0.5286 m
Y Offset	−0.0096 m
Z Offset	0 m
X Length	1.2854 m
Y Length	0.5963 m
Z Length	0 m

, while

Table 5. Results of the Mesh Independence Study in 3D Winglet Analysis.

Mesh Scheme	Element Size (mm)	Number of Elements	Time (Minutes)	Pressure (Pa)
Auto-sized	2520	19,953	7	118,766.00
Refined wing edges	500	116,294	40	136,106.00
Refined wing edges with regions	90	505,390	210	139,333.00
Wing edges and regions are both refined	60	5,761,765	1274	139,331.00

shows the effects of element size on pressure and simulation time. The selected mesh is 28 times (28×) smaller than the auto-sized mesh. The mesh study shows the geometry of the aerofoil correctly and so accurately captures flow phenomena. This results in a pressure increase between these two parameters of 20,567 Pa.

The mesh size of 90 mm was solved approximately 6 times (6×) faster than the simulation with a mesh size of 60 mm but observed a pressure difference of 2 Pa. While the result of the smaller mesh can be considered more accurate, a pressure difference of 2 Pa is only 0.07% of the total pressure value. When making the selection, such a negligible decrease in pressure does not require an extended simulation time. Two regions were used in this simulation, one over the whole wing, and one solely around the leading edge. From the different simulation cases, some effects are seen between simulations 2 and 3. The addition of the region’s decreased mesh size by 410 mm, resulting in a 3227 Pa pressure increase. Figure 7 shows that after regions refine aerofoil shape, the pressure result becomes constant. It was concluded that at this stage, mesh independence had been achieved.

The results of the mesh independence study for the 2D study are presented in

Table 6. Effect of Mesh Size on Pressure in the 2D Dimple Effect Analysis.

Mesh Scheme	Mesh Size (mm)	Elements	Time (Minutes)	Pressure (Pa)
Auto-sized	20	2373	3	4000.28
Refined wing edges	5	7354	12	3983.80
Refined wing edges	2	30,604	27	3951.44
Refined wing edges	1	60,520	44	4091.87
Refined wing edges with regions	1	946,101	1179	4267.88
Refined wing edges with regions	0.5	3,711,713	57,454	4275.42

. The 2D model is less demanding of computational resources. Consequently, a finer mesh can be applied to this study. A mesh size of 1 mm was sufficiently acceptable in this model. Halving the mesh size to 0.5 mm increased simulation time by 76 h and 15 min, or 388.04%. Within this additional time period, the pressure was only affected by 0.178%. For the purposes of this study, such a negligible difference can be ignored. However, it is noted that for any application within the aviation industry, a more accurate result should always be used. For larger scale simulations to be applied on CFD models, results typically require accuracy of $\times {10}^{-9}$. However, these may consume more computational resources. In this analysis, one region was placed over the entire aerofoil; its effects can be seen in Figure 8. Simulation cases 4 and 5 both had a mesh size of 1 mm, but simulation 5 was run with a region of finer mesh surrounding the aerofoil. This action captured flow dynamics more clearly, leading to a pressure increase of 176.01 Pa when compared to simulation case 4. The graph on Figure 8 also shows that refinement of the region past an element size of 1 mm will give the same pressure value consistently and, as such, will be a waste of computational resources. This shows that the 1 mm mesh size is sufficiently acceptable for this CFD model.

2.7. The y₊ Factor

An appropriate y₊ value to calculate fluid flow down to the sub-viscous area of the boundary layer is calculated in this section. First, the local Reynolds number is calculated in Equations (17) and (22), the skin friction of the aerofoil in Equations (18) and (23), followed by the wall shear stress in Equations (19) and (24), then the near-wall boundary fluid speed as a frictional velocity in Equations (20) and (25), and finally the distance to the first mesh node, $y^{+}{}_{cell}$, in Equations (21) and (26) [30]. Respective air density and dynamic viscosity values are taken at 15 °C and at sea level. This considers the increase in aircraft efficiency at take-off where lift and drag forces are at a maximum. The average take-off velocity of larger aircraft were taken, thus the calculations presented in Equations (17)–(26) [30]. The 2D study was also set up at take-off velocity, as presented in Table 3. Given this, the set-up parameters remain constant as for the 3D model, only aerofoil length changes. In the same vein, it was found that a $y^{+}$ value equal to one was suitable to study sub-viscous flow effectively.

$$R{e}_{Object}=\frac{{\rho }_{air}{v}_{air}l}{{\mu }_{air}}=\frac{1.225\times 82\times 13.751}{1.81\times {10}^{-5}}=76.314\times {10}^6$$

(17)

$$C_f=0.058\left(R{e}_{Object}{}^{-0.2}\right)=0.058[\left(76.314\times {10}^6{)}^{-0.2}\right]=0.00154$$

(18)

$${\tau }_{air}=0.5C_f{\rho }_{air}{v}_{air}{}^2=0.5\times 0.00154\times 1.225\times {82}^2=6.342 \mathrm{Pa}$$

(19)

$$v^{*}=\sqrt{\frac{{\tau }_{air}}{{\rho }_{air}}}=\sqrt{\frac{6.342}{1.225}}=2.275{ \mathrm{ms}}^{-1}$$

(20)

$$y^{+}{}_{cell}=\frac{y^{+}{\mu }_{air}}{v^{*}{\rho }_{air}}=\frac{1\times 1.81\times {10}^{-5}}{2.275\times 1.225}=6.4947\times {10}^{-6}$$

(21)

$$R{e}_{Object}=\frac{{\rho }_{air}{v}_{air}l}{{\mu }_{air}}=\frac{1.225\times 82\times 1}{1.81\times {10}^{-5}}=5.55\times {10}^6$$

(22)

$$C_f=0.058\left(R{e}_{Object}{}^{-0.2}\right)=0.058(\left(5.55\times {10}^6{)}^{-0.2}\right)=0.0026$$

(23)

$${\tau }_{air}=0.5C_f{\rho }_{air}{v}_{air}{}^2=0.5\times 0.0026\times 1.225\times {82}^2=10.708 \mathrm{Pa}$$

(24)

$$v^{*}=\sqrt{\frac{{\tau }_{air}}{{\rho }_{air}}}=\sqrt{\frac{10.708}{1.225}}=2.957{ \mathrm{ms}}^{-1}$$

(25)

$$y^{+}{}_{cell}=\frac{y^{+}{\mu }_{air}}{v^{*}{\rho }_{air}}=\frac{1\times 1.81\times {10}^{-5}}{2.957\times 1.225}=4.997\times {10}^{-6}$$

(26)

3. Results

3.1. Effect of Pressure on the Winglet

Table 7. Pressure comparison showing percentage difference between each winglet design and a conventional wing.

ID	Final Pressure (Pa)	Percentage Difference	Force (N) × 106
No-winglet	139,333.00	-	48.39
Winglet-1	123,199.00	11.57	44.36
Winglet-2	119,803.00	14.01	42.90
Winglet-3	116,606.00	16.31	41.61

shows the effect of three delta winglet designs when attached to a conventional wing. According to the table, all winglet designs support the proposed literature, decreasing the influence of pressure over the wing. An average pressure decrease of 13.97% was observed, with the third design being most optimum. Winglet 3 offers a pressure decrease of 16.31% when compared to a conventional wing. It can be stated that for most practical flows, pressure and friction are the two most relevant forces [31]. The decrease in pressure suggests that a reduction in the magnitude of the opposing frictional forces acting on the wing has been observed.

The general Reynold number that describes flow over the whole boundary domain is shown in Equation (27). When compared to the local Reynold number, found in Equation (28), the Reynold number over the domain is 14.017 × ${10}^6$ smaller. This decrease is to be expected given the non-slip conditions of the wall previously described. Both local and general flows are categorically turbulent, with a Reynold value less than ×${10}^3$ needed to make the flow laminar [32]. It is important to note that despite being turbulent, the flow over the boundary domain is closer to being laminar. A decrease in the pressure gradient found on the leading edge, as found in Table 7 and Figure 9, implies that the wing is separating flow more effectively. Consequently, the refined wing is more adapted to induce less disturbance in this region of more laminar flow. A wing that can maintain a more structured, laminar flow for 50% of its chord length or more, will create substantially less drag [33,34].

The total force on the wing has been calculated using Equation (28) [28] and displayed in Table 7. Inevitably, the ultimate extent to which these factors have been positively affected is a combination of the lift and drag forces. Figure 9 shows the effects of pressure on the force exerted on the wing. An $R^2$ value of 0.9852 shows that results have a near-perfect positive correlation, and that as pressure increases, so does the force exerted on the wing. This value is almost significant to a value of 0.1, meaning as pressure over the object increases, force will increase at 99% of the rate.

$$Re=\frac{vl}{{\mu }_{air}}=\frac{82\times 13.751}{1.81\times {10}^{-5}}=62.297\times {10}^6$$

(27)

$$F=pA=139333\times 347.358=48.398\times {10}^6$$

(28)

3.2. Effect of Pressure on the Winglet

The effect of the wake on the winglet for turbulence flow model is presented in this section. Figure 10 displays a pressure gradient on the wing before modifications were made and for the improved design. Two distinct differences can be observed. First, there is a supplementary pressure gradient on the leading edge of the aerofoil where the winglet has been attached to the main body. This helps to direct high-pressure airflow away from the top of the wing, as hypothesised. These effects are then seen in the wing’s wake where there is inherently less disturbed flow.

3.3. Effect of Flow Phenomenon on the Winglet

The effect of the wake on the winglet for turbulence flow model. There is a region of turbulent air on the winglet itself. It can be seen in Figure 11a that this area of low pressure is forming on the open side of the winglet and spreading upwards. Traces, shown in Figure 11b, were placed over this part of the wing to capture flow phenomena in this region, but a corresponding level of flow delamination cannot be identified. To deduce the cause of the pressure increase, it is recommended that a wind tunnel analysis be conducted, where smoke is released from the inlet. This alternative form of analysis will show the development of any smaller eddy currents that cannot be displayed with current computational power. Figure 11c,d shows the formation of the low-pressure region over the other winglets. Figure 11c shows the largest low-pressure area induced on winglet 1, while Figure 11d shows the largest low-pressure area induced on winglet 2, and has an area in between both winglets.

3.4. Effect of Mesh Size on Pressure Stabilization

It is paramount to ensure that the selected boundary domain does not affect experimental results. The 2D simulation demanded much less computational power than the 3D simulation. As a result, the initial auto-sized mesh will represent aerofoil geometry with far more accuracy. It follows that a boundary independence evaluation could be conducted with a more representative auto-sized mesh for this study. The results will fluctuate with boundary domain changes but not due to auto-sized mesh inaccurately capturing flow phenomena between models. The reason is that the mesh size is dependent on the stabilization of the pressure.

Table 8. Pressure fluctuation as boundary domain changes size.

Boundary Domain in Terms of c	Elements	Pressure Value (Pa)
Y value
6	2642	4117.41
12	2690	4291.67
25	2455	4330.38
30	2780	4337.83
35	2714	4259.73
X value in front of aerofoil
3	2676	4139.25
6	2816	4376.91
15	2777	4412.88
17.5	2847	4150.67
X value behind aerofoil
10	2751	4777.21
15	2777	4412.88
20	2869	4267.74
25	2939	4008.67
30	2373	4000.28
35	2576	4236.08

shows the results of the study based on the domain. Dimensions for each axis were selected once a point of relative stability had been reached; these points are 30 c, 15 c, and 30 c, respectively. The subsequent boundary domain is 32.17% bigger than the domain size suggested by Hassan et al. [35]. Results in the mesh size effect may have been affected by a slightly more constrictive boundary domain. However, given that the boundary domain was kept constant across all studies, it is reasonable to suggest that the effects of the flow phenomenon would still have been captured at the same percentage rate. Figure 12 identifies the stability points located for this dimple study.

3.5. Effect of Dimple on Separation Point

As a flow moves along an object, friction between the fluid and the object develops a surrounding boundary layer of turbulent flow. As the Boussinesq’s turbulence hypothesis has already explained, the turbulent flow will function as if it were vicious and effectively stick to the wall of the aerofoil. As shown in Figure 13 and Figure 14, this behaviour is replicated in CFD through the use of a non-slip wall, where velocity is made equal to zero. Given the large pressure gradient developed near the wall, y₊ is a metric developed for CFD to calculate how fluid flow behaves in this sub-viscous area. Accurate representation of the fluid mechanics in this sub-viscous layer allows the model to develop the boundary layer accurately as the flow progresses along the non-slip wall of the object [36]. Monitoring the boundary layer development will become important when locating the separation point between aerofoil and fluid flow.

3.6. Effect of Dimple on the Lift and Drag Coefficients

Figure 15 establishes that the dimple effect has proven to be temperamental when changing several flow characteristics for the better. Firstly, it can be noted that for all three dimple geometries, a dimple underneath the aerofoil performs the best out of the three locations. Figure 16 visualises this effect for the ellipse-shaped dimple and compares it to an unmodified aerofoil. It is observed that the dimple located on the bottom of the aerofoil produces the highest amount of lift, while any alternative positioning decreases lift and increases drag. When dimple geometries were placed on both the top and bottom of the aerofoil, the lift effects cancel each other out and create 5.36% more drag. The increase in drag can be explained through the study of Figure 15, Figure 16, Figure 17 and Figure 18. It is noticed that airflow collides with the dimple’s back instead of flowing directly over it. A plausible explanation is that the dimple is too large in length. As the near vacuum created in the dimple sucks in the surrounding air, the flow is pulled into the dimple at an unexpected rate. Consequently, as momentum is conserved in the x-direction, flow hits the back wall of the dimple, thus slowing down the total velocity and creating a region of higher pressure. It follows that an arrangement of two dimples will double this induction of drag. Experimental results support this finding to an accuracy of $\pm 0.005$.

When this disturbance occurs in a square-shaped dimple, the near perpendicular wall causes an area of even higher pressure to develop; hence, the square-shaped dimple is the least efficient. This drag increase is displayed in Figure 17. A drag increase is to be expected when the dimple is placed underneath the wing, provided that lift results are taken into consideration. In the aviation industry, it has been reported that as lift increases, so does drag [37]. As observed in Figure 17, this positioning induces more lift, and as a result of that, more drag. When the dimple was located underneath the aerofoil, lift efficiency across the three geometries increased by an average of 0.623, or 218.60%. In contrast, when the dimple was placed on the top, efficiency decreased by 0.522, or 183.09%. The increased level of force on the underside of the aerofoil is natural and explained through Figure 18; (see the governing equations in Section 2). When air flows underneath the aerofoil, the pressure summation from both axial and normal forces increases the total force exposure, hence the slightly larger pressure gradient.

3.7. Effect of Dimple on Vortex Induced Vibration

The results presented in Figure 19 represent the vorticity plots for the different designs considered for the aerofoil dimples and the locations. The results show the effect of vortex-induced vibration (VIV) on each of the designs, which shows different flow patterns around the ends of the aerofoils where the dimples are located. The effect of the vortex around the end of the aerofoil has a noticeable impact on the drag.

3.8. Effect of Dimple on Flow Induced Vibration

The results from Equation (29) show that the flow over the aerofoil is unconditionally turbulent. Although the progression of flow after the dimple is kept further towards a laminar state, as shown in Figure 19, the turbulent state can be attributed to the high inlet velocity. The small amount of flow reversal highlighted in Figure 20 and Figure 21 might increase the Reynolds number for the modified aerofoil. Flow reversal should be prevented if this geometric modification is completely optimised. However, despite the presence of flow reversal, the pressure gradient displayed in the figures above show minimal detriment to the overall direction or structure of the flow. The whole boundary Reynolds number for this 2D simulation is calculated as in Equation (29) [34]. Further work on flow around the dimples are given in Figure 22 and Figure 23.

$$Re=\frac{vl}{{\mu }_{air}}=\frac{82\times 1}{1.81\times {10}^{-5}}=4.530\times {10}^6$$

(29)

As seen on Figure 20 and Figure 21, the value of the R² < 1, for both the velocity and pressure profiles, which means the values are consistent with very low variance and were increasing based on different geometric features, as shown on

Table 9. The pressure and velocity distribution for the 10 dimple designs.

Pressure Plot		Velocity Plot		Grouped Plot
Descriptors	Pressure (Pa)	Descriptors	Velocity (m/s)	Descriptors	Velocity (m/s)	Pressure (Pa)
Triangle Dimple Both	106,800	No Dimple	109.453	Round Dimple Both	157.581	106,880
No Dimple	106,828	Square Dimple Bottom	109.706	Round Dimple Bottom	122.988	106,847
Square Dimple Bottom	106,845	Round Dimple Top	109.724	Round Dimple Top	109.724	107,061
Round Dimple Bottom	106,847	Square Dimple Both	109.79	Triangle Dimple Both	136.903	106,800
Triangle Dimple Bottom	106,872	Triangle Dimple Bottom	118.837	Triangle Dimple Bottom	118.837	106,872
Round Dimple Both	106,880	Round Dimple Bottom	122.988	Triangle Dimple Top	133.365	106,886
Triangle Dimple Top	106,886	Triangle Dimple Top	133.365	Square Dimple Both	109.79	106,893
Square Dimple Both	106,893	Triangle Dimple Both	136.903	Square Dimple Bottom	109.706	106,845
Round Dimple Top	107,061	Round Dimple Both	157.581	Square Dimple Top	248.701	107149
Square Dimple Top	107,149	Square Dimple Top	248.701	No Dimple	109.453	106,828

. It was observed that the design with no dimple had the least velocity distribution. It was also observed that the dimple design with a triangled bottom had the average velocity and pressure distribution. This is because it does not hold onto any flow to recirculate or create any C-pillar, as noticed in Figure 16, Figure 17 and Figure 18. Details of the exploded view are shown in Figure 22 and Figure 23. On closer inspection, it can be observed that the dimples create some high-vortex effects.

3.9. Validation Results

Validation of computational models from commercial tools, such as Autodesk CFD, is quite a challenge due to limited trust and robustness in CFD modelling. Table 9 shows the result of the mesh analysis from convergence conducted in Autodesk CFD 2019. As presented in

Table 10. Results of a dimpled wing in comparison to a smooth one.

Dimple Type	Design Case Reference Point	Pressure (Pa)	Lift	Drag	Lift Coefficient	Drag Coefficient	Percentage Difference (P.D) of This Model Compared to a Conventional Aerofoil		Efficiency
							Lift P.D.	Drag P.D.
No Dimple	-	4267.88	19.486	68.355	0.154	0.5390	-	-	0.2850723
Ellipse	Top	4182.86	−49.581	69.810	−0.389	0.5480	−254.444	2.1290	−0.710229
	Bottom	4264.90	70.458	72.289	0.554	0.5679	261.580	5.7565	0.9746588
	Both	4192.92	0.363	71.994	0.003	0.5633	−0.732	5.3243	0.0050412
Square	Top	4188.69	−31.975	71.785	−0.250	0.5609	−164.090	5.0185	−0.44542
	Bottom	4265.64	60.069	71.926	0.469	0.5620	208.266	5.2253	0.8351424
	Both	4174.82	9.543	73.120	0.074	0.5662	−51.024	6.9720	0.1305168
Triangle	Top	4190.9	−29.110	70.967	−0.229	0.5586	−149.391	3.8214	−0.410197
	Bottom	4265.28	65.827	71.921	0.518	0.5661	237.815	5.2171	0.9152678
	Both	4275.29	26.732	73.003	0.210	0.5735	37.188	6.8009	0.3661807

, the drag forces remain relatively constant, with a maximum deviation of 6.97%, while the resultant lift force deviated from the normal by a maximum of positive 261.58% and negative 254.44%. Given the large fluctuations from the unmodified aerofoil, the study was replicated in SimScale platform to ensure validity. It is important to note that the results from the SimScale verification study are not directly comparable to the Autodesk CFD analysis, due to some refinements in the geometry and difference in enclosure by a few millimetres. However, the results both show good agreement.

Second, the boundary domain for the SimScale analysis had to be made substantially smaller for the simulation to create a well-structured mesh, as shown in Figure 24. Third, different software uses different scale factors, hence the larger drag coefficient. Effects of the scale factor can be nullified when looking at the percentage differences for lift and drag. As indicated in

Table 11. Results from the SimScale verification study.

Dimple Type	Design Case Reference Point	Lift Coefficient	Drag Coefficient	Percentage Difference (P.D.) of This Model Compared to a Conventional Aerofoil
				Lift P.D.	Drag P.D.
No Dimple	-	0.00103	13.3150	-	-
Ellipse	Top	−0.00071	13.3212	−168.547	0.1282
	Bottom	0.00409	13.3196	296.316	0.1399
	Both	0.00273	13.3455	164.954	0.0541
Square	Top	−0.00113	13.4501	−109.759	0.8384
	Bottom	0.00343	13.5110	232.368	1.4720
	Both	0.00183	13.3470	77.715	0.0658
Triangle	Top	−0.00131	13.3143	−126.939	−0.1794
	Bottom	0.00353	13.3167	242.058	0.1613
	Both	0.00017	13.2998	16.473	0.2887

, the ellipse-shaped dimple located at the bottom of the aerofoil performs the best, yielding a 296.316% lift increase, and a 0.0541% drag increase when compared to an unmodified aerofoil. This comparative study shows the same trend, with the ellipse-shaped dimple underneath the aerofoil yielding the best aerodynamic performance, while the square-shaped dimple on top of the aerofoil was the worst.

Figure 24a was obtained from the SimScale simulation while Figure 24b was obtained from Autodesk CFD 2019. These mesh plots show the flow phenomenon and its effect over the aerofoil. Postprocessing of results visualises the cavities’ impact. Further post-processing was conducted using TecPlot 360 [38]. When looking at the surface where the dimple is located, there is considerably less low velocity after the dimple. Evidence suggests that flow remains in a state of equilibrium for longer, meaning the separation point has been successfully delayed. Due to the nature of the smaller boundary domain, results for this analysis might be slightly misrepresented, despite the small mesh shown in Figure 24a. Having proved the validity of the results, the presentation of these results will continue by only considering the most suitable CFD results.

4. Discussion

A CFD investigation was conducted on the aerofoil of an aircraft in this paper. Two different design modifications were performed on a conventional wing design to optimise the performance capabilities of freighter aircraft by increasing the lift-to-drag ratio. The findings from this study were for freighter aircrafts and those of similar aircraft classification. Thus, this study has a high positive impact on the aeronautical industry. Finite element analysis was performed by primarily presenting the mesh analysis and the convergence results obtained to compare results. Adding a winglet to a Boeing Dreamliner 787-9 wing reduced the pressure gradient found on the leading edge by a maximum of 16.31%. The design changes made in the recommendations section could help to decrease the pressure value by another 2–4%. To ensure the validity and stability of the results, the SimScale platform was also used, as the results should be stable to the power of ×10⁻⁹ compared to the CFD results.

The dimple effect study observed flow in the viscous boundary layer formed around the aerofoil. The separation point was delayed, leading to an increased aerodynamic performance. Results from the dimple study, having been evaluated, showed a high efficiency increase of 241%, and was found when an elliptical dimple was placed underneath the aerofoil. The triangle-shaped dimple model presented the second-best results at a 221% increase, and square-shaped ones were the worst at 193%. The position of the dimple on the wing also affected the results. Dimples placed underneath the wing are optimum, while dimples placed solely on top were found to decrease aerodynamic performance by up to 254%. Both Autodesk CFD and SimScale OpenCFD simulations deviated in a comparable way and to a similar magnitude. It cannot be ignored that some deviations in the efficiency results in this paper are considerably high. The results are to be taken conservatively, subject to further experimental validations, such as wind-tunnel analysis and additional CFD simulations, are recommended.

5. Concluding Remarks

The numerical investigation of a NACA 0017 with investigation on the effect of dimples in addition to the aerodynamic effect of a winglet of a freighter aircraft was performed in this paper. Overall, the paper is application-oriented and the specific outcome would be useful in the aeronautical society. This study presents some novelty in the application of aerofoils and winglets in the aviation industry. It also presents an understanding of the aerodynamic effect of a dimpled wing on the NACA aerofoil. In this study, much attention is spent to the design of surface dents close to the trailing edge. The manuscript describes a detailed study of aerofoil modifications that can decrease drag, thus increasing flow efficiency.

This research presents the following highlights from the models and methods: First, it assesses the effects of a dimple at a separation point where the aircraft is operating at 0° angle of attack (AOA). Second, the dimple effects were analysed when a dimple was both on top of and underneath the wing, in addition to one dimple either on top or underneath. Third, this study looked at the effects of the dimple on a larger aircraft wing with a more considerable spar length of aerofoils. Fourth, there is novelty in the design seen in the choice of dimple geometry with intricate vorticity plots for 10 model geometries of the NACA 0017 aerofoil. Finally, the methodology includes the novel application of an OpenFOAM platform called SimScale for the CFD analysis and verification of the model.

From this study, the following conclusions and recommendations were drawn:

(1)

It was observed that the dimple effect study presented a flow profile in the viscous boundary layer formed around the aerofoil. It was also noticed that each end of the dimple had a different characteristic effect in the flow distribution around the dimple. Also, there was some delay in the separation point that led to increased aerodynamic performance. Further studies on the different angle of attack (AOA) and the dimple effect at varying velocities, such as cruising and taxi, could provide insight as to whether each modification yields the same performance at lower velocities or Reynold numbers. In addition, the power efficiency and energy demand to study the fuel efficiency due to the winglet optimisation can be further investigated.

(2)

The investigation on the effect of wing extensions on aircraft, the dimple effect, and the precise effect of winglets on aircraft that already have larger aspect ratios showed good results for decreased drag effects on freighter aircraft, presented in Section 3 and Section 4. The addition of a winglet reduced the pressure gradient found on the leading edge by a maximum of 16.31%. More significantly, a maximum efficiency increase of 241% was found when an elliptical dimple was placed underneath the aerofoil. It cannot be ignored that some deviations in the efficiency results in this paper are considerably high. The results are to be taken conservatively and are subject to further experimental validations, such as wind-tunnel analysis and additional CFD simulations.

(3)

From this investigation, it was noticed that the location of the dimple on the wing also affected results. Dimples placed underneath the wing were optimum, while dimples placed solely on top were found to decrease aerodynamic performance by up to 254%. Both Autodesk CFD and SimScale’s OPENFOAM simulations deviated in a comparable way and to a similar magnitude. Note that most of the investigations are done with CFD, but the problem is very challenging for CFD, particularly computing the influence of dent details in such a flow. The dent is near the trailing edge of the foil, so the oncoming flow is highly turbulent. Fortunately, as seen here, the flow on the aerofoil until the dent is attached, and the modelling of attached flow, is not very critical. However, the situation changes completely when the flow starts to separate. There is no RANS model that can reliably model separated flow, despite the popularity of SST-k-omega, one of the comparatively better models. Hence, further research is recommended using the modern generation of LES models, which are not yet available in commercial codes.

(4)

The effect of the modification to a typical freighter aircraft geometry can restructure the airflow surrounding the wing to a more laminar state at take-off. In both cases, flow phenomenon surrounding the aerofoil can be observed in a positive light; hence optimisation has been achieved. However, the problem is that RANS models are too generous in adding turbulent eddy viscosity, such that effectively a flow at a much lower Reynolds number is solved. With more detailed CFD specification, the study can be enhanced in future research to include numerical robustness and the discretization of the convective terms, which does add viscosity effect.

(5)

In addition, an extensive grid refinement study should be performed, but at effectively a too low Reynolds number, where grid requirements are much lower. The problem is that it does not warn against errors caused by a too-diffusive turbulence model. CFD models required high validation levels to obtain high confidence index in the CFD results, such as (a) the use a central discretization (at least in the surroundings of the dent), and (b) changing the turbulence model (lower the amount of eddy viscosity) and see how much the results are changing. Also, detailed validation is highly desirable, especially inside the dent.