The numerical method employed in this study entails solving the general governing equations of fluid flow, the Navier-Stokes equations (NSE). By solving the equations numerically for a specific set of boundary conditions, the flow properties in a given geometry are predicted [2, 15]. In this study, the flowfield is apprximated as a steady, turbulent and viscous incompressible flow. Hence, the relevant form of the NSE is the Reynolds-average Navier-Stokes (RANS) equations, which are split into mean and fluctuating components [16]. The total velocity, \({u}_{i}\) is decomposed as a function of the mean velocity \(\stackrel{-}{{u}_{i}}\) and the fluctuating component \(\stackrel{´}{{u}_{i}}\) as shown in the Eq 1. In the same manner, averaging can be applied to all other flow variables as shown in Eq 2 for pressure.
$${u}_{i}=\stackrel{-}{{u}_{i}}+\stackrel{´}{{u}_{i}}$$
1
$$p=\stackrel{-}{p}+\stackrel{´}{p}$$
2
The averaged from of the continuity and momentum equations are given in Eq 3 as:
$$\frac{\partial }{\partial {X}_{i}}\left({\stackrel{-}{u}}_{i}\right)=0$$
3
$$\frac{D{\stackrel{-}{u}}_{i}}{Dt}=-\frac{1}{\rho }\frac{\partial P}{\partial {x}_{i}}+\frac{\partial }{\partial {x}_{j}}\left[\mu \left(\frac{\partial {\stackrel{-}{u}}_{i}}{\partial {x}_{j}}+\frac{\partial {\stackrel{-}{u}}_{i}}{\partial {x}_{i}}-\frac{2}{3}{\delta }_{ij}\frac{\partial {\stackrel{-}{u}}_{k}}{\partial {x}_{k}}\right)\right]+\frac{\partial }{\partial {x}_{j}}\left(-\stackrel{-}{\stackrel{´}{{u}_{i}}\stackrel{´}{{u}_{j}}}\right)$$
4
In The above equations, \(u\), \(\rho\) and \(\mu\) represents the velocity, density and dynamic viscosity of the fluid. The term \(-\stackrel{-}{\stackrel{´}{{u}_{i}}\stackrel{´}{{u}_{j}}}\) in Equation 4 is known as the Reynolds stress tensor (RST) and is defined in Equation 5. The RST accounts for the turbulent fluctuations and is required to be modelled to closed the RANS equation. Modeling the RST to obtain a solution requires some approximations in the governing differential equations at various levels of sophistication, leading to different turbulence models employed in predicting engineering flows. However, each turbulence model has some limitations, depending on the nature of the flowfield [11]. Today, numerous commercial and proprietary software have been developed for CFD analysis of various industrial and non-industrial flows. In this study, ANSYS-Fluent 14.0 CFD software is employed in modeling the flowfield of IU-18 UAV with the aim of predicting complex interaction of lift and drag forces involved.
The coefficients of drag and lift force \({C}_{d}\) and \({C}_{l}\) are defined in Equations 6 and 7, respectively.
$${C}_{d}=\frac{{F}_{d}}{0.5\rho {V}^{2}A}$$
6
$${C}_{l}=\frac{{F}_{l}}{0.5\rho {V}^{2}A}$$
7
where \({F}_{d}\) and \({F}_{l}\) are the drag and lift forces, respectively. \(\rho\) is the density of the air, \(A\) is the frontal area of the vehicle, and \(V\) is the relative speed.
The true aerodynamic efficiency of the wing is defined by the ratio of \({C}_{l}\) and \({C}_{d}\) [4] as.
$$\text{A}\text{e}\text{r}\text{o}\text{d}\text{y}\text{n}\text{a}\text{m}\text{i}\text{c} \text{e}\text{f}\text{f}\text{i}\text{c}\text{i}\text{e}\text{n}\text{c}\text{y}=\frac{{F}_{l} }{{F}_{d} }=\frac{{C}_{l}}{{C}_{d}}$$
8
Determining the Reynolds number of the flowfield, \(Re\) is vital in selecting the turbulence models in numerical analysis. \(Re\) is a non-dimensional parameter that represents the ratio of inertia and viscous forces [17]. It provides a quantitative representation of the degree of turbulence in a flow and is represented in Equation 9 [18, 19].
$$Re=\frac{Inertia force}{Viscous force}=\frac{\rho VD}{\mu }$$
9
where \(\rho , V, D \text{a}\text{n}\text{d} \mu\) is the air density, free stream velocity, characteristic length and dynamic viscosity, respectively. The D for an aerofoil is equal to the aerofoil length. For a wing, D is equal to the aerofoil chord length. While for a non-uniform wing, D is equal to the Mean Aerodynamic Chord, MAC length [20]. For fully developed turbulent flow, \({R}_{e}> \approx {10}^{6}\) [17].
For the IU-18 UAV at cruising altitude, the Reynolds number is obtained based on Equation 9 as follows:
$$Re=\frac{0.771\times 46\times 0.9}{1.66\times {10}^{-5}}=1.923\times {10}^{6}$$
Thus, the flow can be approximated as fully turbulent at the cruising altitude of the UAV. However, this does not eliminate the presence of low Reynolds number regions in the flowfield, particularly, at high AoA, where separation is highly likely.
Based on the Reynolds number of the IU-18 UAV computed above, the standard and realizable \(k-\epsilon\) turbulence models are selected for this study. The standard \(k-\epsilon\) model proposed by Launder and Spalding [21] is the most widely used turbulence model for industrial applications. The model has been acknowledged to be robust, economical and reasonably accurate for a wide range of industrial applications of fully developed turbulent flows. However, it performs poorly for flows with strong separation, large streamline curvature and large pressure gradient. Hence, in this study, the Realizable \(k-\epsilon\) turbulence model will be employed for the purpose of comparison. This model offers the good performance among the versions of \(k-\epsilon\) model. Hence, it is considered appropriate for modelling flows with separation and complex secondary flow feature [22, 23]. The numerical computation of the flowfield investigated in this study was carried out using Ansys-Fluent 14.0 CFD code. The details of the computer-aided design (CAD) of the wing geometry, numerical setup, boundary conditions and turbulence models employed are presented in the following subsections.
3.1 CAD Model of Wing of IU-18 UAV
The wing of the Ichoku-18 UAV has the following dimensions: a span of 11 m, a root chord of 1.20 m, a tip chord of 0.50 m, an aspect ratio of 11.68 and a gross area of 10.36 m2, Table 1 presents the remaining specifications of IU-18 UAV wing.
Table 1
Specifications of Wing of Ichoku-18 UAV.
Constraint
|
Dimensions
|
Trailing edge sweep angle
|
4.45o
|
Leading-edge sweep back angle
|
4.39o
|
MAC
|
0.90 m
|
Fitness ratio
|
13.0
|
Anhendral
|
3o
|
Twist
|
-2o
|
Taper ratio
|
0.42
|
Leading-edge sweep
|
0o
|
3.2 Numerical Setup
a) Geometry design modeller
The CAD model of the wing developed in CATIA V5 was situated in a computational domain (CD) with dimensions as shown in Figure 2. Subsequently, the CD was imported into ANSYS-Fluent 14.0 software for meshing and analysis.This study employed the C-type CD due to its simplicity and time savings [24].
Due to symmetric nature of the wing, only half of its CAD model was developed. To model variation in flight attitude, the study opted for changing the flow field velocity vectors and its orientation relative to the wing. This has been shown to present the same effect as changing the flow AoA as described in [24] and and employed by Roy et al. in [11].
b) Boundary conditions
The boundary conditions prescribed include ‘Inlet Velocity’ on the first face of the CD in the Z-direction, "Pressure Outlet" on the face opposite the inlet, and "wall” on the wing model under investigation. The rectangular walls of the CD were defined as ‘Symmetry’. Figure 3 shows the boundaries conditions prescribed on the CD.
c) Meshing
In this study, an unstructured mesh was generated on both the wing model and CD. To improve the grid quality, advance size function, relevance centre and smoothing in Ansys Fluent 14.0 were changed to ‘Proximity and Curvature’, ‘Fine’ and ‘High’, respectively. The Proximity and Curvature option was employed to refine the mesh over much of the model and control the mesh density in regions of the model where features are located more closely together without using numerous local controls. The fine smoothing option allows controlling the variation in the size of the mesh elements thereby improving the accuracy of the numerical analysis. The relevance centre provide finer mesh in the swept regions closer to the geometry [24]. Thereafter, the meshing is generated and done on both the wing model and CD.
To ensure the solution obtained is independent of the grid size, grid independence check was done before adopting the finer grid size. This was achieved through performing the simulation with finer grid size. The results of the check show that the analysis grid size has some effect on the CFD results. However, as show on Table 2, the effect of using the finer the grid on the values of Cd and Cl obtained was less than an order of magnitude, which is not significant. Figure 4 shows the mesh generated on the wing model and the CD. Subsequently, the finer mesh consisting of 86,546 nodes and 482,214 elements for the wing and the CD, respectively was adopted for the remaining calculations.
Table 2
Grid independence check results.
Grid size
|
AoA
|
Nodes
|
Elements
|
\(\mathbf{C}\mathbf{d}\)
|
\(\mathbf{C}\mathbf{l}\)
|
Coarse
|
0o
|
67,654
|
376,935
|
0.044
|
0.321
|
Finer
|
0o
|
86,546
|
482,214
|
0.026
|
0.246
|
d) Setup
The model was set up using a three-dimensional, steady-state, pressure-based coupled algorithm. The turbulence model implemented in this study is the standard and realizable \(k-\epsilon\) turbulence model. The input parameters are the cruising speed and altitude of the UAV, given as 46 m/s and 4,572 m, respectively. The remaining input parameters relate to properties of the surrounding air at altitude of 4,572 m, namely, pressure, viscosity, density and temperature, which are given as 57.182 kPa, \({1.66\times 10}^{-5} kg/m.s\), 0.771 kg/m3 and 258.43 K, respectively.
2.3 Solution initialization
To initialize the solution from the boundary conditions specified, a "Hybrid Initialization" was employed. Hybrid initialization is a collection of recipes and boundary interpolation methods. It solves the Laplace equation to produce a velocity field that conforms to complex domain geometries, and a pressure field that smoothly connects high and low-pressure values in the CD. All other variables (velocity, pressure, turbulence and VOF) were patched based on domain averaged values or a predetermined recipe [25].After initialization, seven (7) different cases namely; \({\text{A}\text{o}\text{A}=0}^{^\circ }, 4^\circ , {8}^{^\circ },12^\circ ,{16}^{^\circ }, 20^\circ \text{a}\text{n}\text{d} {24}^{^\circ }\) were simulated for the two turbulence models. Lift and drag coefficients values were monitored and saved for each iteration. The convergence history was also monitored. The residual Root-Mean-Square (RMS) Error default value was set to \({10}^{-6}\) for better convergence. In some simulation where convergence was not possible, calculations were run for a stabilized solution determined by low fluctuations in the residuals.