Winglet Geometry Impact on DLR-F 4 Aerodynamics and an Analysis of a Hyperbolic Winglet Concept

In this article, the growth of aerodynamic efficiency and the growth of the wing structural stress is studied for DLR-F4 typical transport aircraft wing-body, after installing classical Whitcomb winglets of different configurations and a delta wingtip fence. A new-concept curved-span winglet was mathematically developed and approved through Computational Fluid Dynamics (CFD) and static structural experiments, revealing the interaction of suband transonic air flow dynamics with the wingtip device geometry. The design space of the winglet geometry was explored briefly, and an evaluation of the lift-to-drag ratio increment depending on various winglet input parameters was performed. In particular, the winglet cant angle effect on lift and drag was thoroughly analyzed at various flow regimes and angles of attack, revealing an ambiguity and a conflicting character of results between highly canted winglets and nearly vertical ones. As a result of cant angle impact analysis, a curved winglet concept is suggested and mathematically parametrized, that could provide an innovative solution, alternative to a morphing winglet, but much simpler with a fixed structure. In conclusion, a multidisciplinary winglet efficiency estimation criterion is suggested for comparing the aerodynamic efficiency of different wingtip devices with respect to their structural weight penalty in real flight conditions.


Introduction
The potential advantage of using non-planar lifting surfaces to reduce induced drag has long been known.These include closed wings, C-wings and wings with winglets.Despite recent technological advances and the use of high rigidity Carbon Fiber Reinforced Plastic wing panels that allow very large aspect ratios, which greatly cut the induced drag at cruise, most airframe manufacturers are still applying winglets, both for further pushing the boundaries of the wing aerodynamic efficiency and due to their aesthetic value.Since their introduction in 1975 by National Aeronautics and Space Agency (NASA) engineer Richard Whitcomb [1,2], winglets have undergone major geometrical changes, driven by the quest for a sustainable efficiency throughout the entire flight envelope of modern airliners.Due to the largely different flow field nature at the wing tip area at different flight conditions, and the challenge of achieving a winglet shape optimal throughout the whole flight envelope, recent winglet geometry optimization studies in the last decade were mostly focused on cruise as a 'design case' for this wingtip device, resulting in trade-off solutions less optimal for take-off and climb at high angles of attack [3][4][5][6][7][8][9][10]. Few innovative solutions have been suggested including morphing winglets which adapt their cant and/or twist depending on the flow regime [11][12][13], or using an integrated moving device such as a winglet-integrated rudder [14,15] or a gust alleviating conventional aileron [16], as well as active vortex wake control with an oscillating mechanism [14,17,18].Bio-inspired devices were studied with a triangular leading-edge extension [19], with multiple elements, where each element is set at a given angle and contributes at a certain flight regime, thus sustaining overall efficiency [20,21].To date, none of these novel solutions could find application on an operational airliner, mainly due to their high complexity which both compromises reliability, and involves a weight penalty, potentially cancelling their aerodynamic gains.
Most winglet geometry optimization attempts [3][4][5][6][7][8][9][10] are purely experimental, based on the results of either computational and/or wind tunnel experiments on a large population of samples with a subsequent optimization using a statistical approach such as pareto front.Unfortunately, the general lack of a winglet local-flow-field theory makes it challenging to explain experimental results, and to uncover proper improvement directions for future designs.Keizo et al. [3] performed high fidelity CFD and wind tunnel multidisciplinary design exploration of a commercial jet winglet, which provides a fairly deep insight into the winglet design space, allowing an optimal design decision based on pareto front and simple trade studies.Although this is a pragmatic approach usually used by airframers to retrofit operational airliners with winglets [6,7], it results in an improved initial geometry, rather than a novel design, due to the lack of a closed-form relationship between the winglet geometry and overall lift-to-drag ratio.Given the winglet is a finite span lifting surface, classical wing theory has been extended to predict the winglet performance.For instance, a CFD combined with the lifting line theory approach has been used by Jan Himisch [4] and Streit et al. [6] to obtain an optimal cruise L/D ratio for a generic transonic airliner.An attempt of finding a winglet optimal twist was performed by Neal [10], by relating the winglet root and tip incidence angles with the induced drag reduction, obtained as output from potential flow solutions in the Trefftz plane behind the wing.However, there are virtually no recent theoretical or experimental studies of the winglet geometry impact on its local flow field and how this correlates with overall efficiency.Therefore, the current work is aimed at filling this gap through a quantitative analysis.
In this paper, a brief design space exploration of a wingtip fence and a classical winglet is presented and a detailed study has been performed on the winglet local flow field, where DLR-F4 lift and drag are correlated with changes in the winglet local angle of attack spanwise, which in turn changes according to the aircraft angle of attack and the winglet's own geometrical parameters, in particular its cant angle.CFD and analytic search for an optimal cant angle has led to considering a spanwise-variable cant (curved) winglet concept, alternative to cant-morphing solutions presented in [11][12][13].The curved winglet span was parametrized using a second order function, which was in turn optimized for getting the most beneficial local angle of attack distribution.Recently, the concept studied herein, of a curved-span winglet, has become relatively well-established with the introduction of blended winglets by Aviation Partners, Inc. [9,22], and Airbus wide body A350 XWB and very lately, the new long range A330 Neo, both featuring a double curvature (spanwise and plan form) winglet [23,24].At the concluding part, a comparison of the static structural stress growth on the wing structure is given for each of the CFD-analyzed wingtip devices, allowing an overall evaluation of their weight penalty, and its juxtaposition to their aerodynamic benefits.
An important contribution of this paper is the front-view parametrization and optimization method used herein for a non-planar lifting surface, applying as an objective function the spanwise local angle of attack distribution.Whereas most parametrization studies for planar wings seek either an optimal plan form, optimal span-wise airfoil shape variation and/or chord-wise airfoil curvature variation (see [25] for instance), the front projection optimization issue is paramount for non-planar lifting surfaces as wings with winglets, c-wings and closed lifting surfaces, where both the geometric and effective local angles of attack do not change uniformly spanwise, when the aircraft angle of attack changes.Examples of other approaches for winglets parametrization can be found in [25,26].The evolutionary geometry parametrization approach presented by Zingg et al. in [25] can be applied to parametrize a span-curved winglet geometry by starting from a classical straight winglet then using B-spline approximation to get a generic curved winglet, and by inserting additional knots at the winglet root and/or tip sections, based on gradual inclusion of different design variables.A comprehensive study by Luciano et al. [27] gives a method for predicting the minimum induced drag of non-planar wings through a configuration-invariant analytic formulation of the unknown circulation distribution.The authors present a numerical tool implementable in MATLAB, which allows the definition of non-planar wings' optimal circulation distribution with minimum induced drag.

Methodology and Tools for CFD and Structural Stress Simulation
The task of numerical simulation of air flow over the wing cantilever, and its stress-strain state under aerodynamic loads, was solved in ANSYS Workbench environment, allowing a combined simulation of the model using several modules from different disciplines, including a full analysis of the interaction of gas dynamics with the wing structure, and aeroelasticity phenomena.The project scheme, presented in Figure 1 shows the steps of a one-way data interchange between different modules: first, the initial solid 3D model is transmitted from the geometry CAD module to the aerodynamic calculations module (Fluent), where the mesh settings of airspace domain around the model are configured, and CFD calculations are performed.The solution, including air pressure field on the wing surface, is then transferred to the static structural strength calculation module, where the grid of the internal structure of the model is built, and its static stress-strain calculation is carried out.If any changes are applied to the original input geometry (a modified winglet twist angle for instance), subsequent modules automatically update the grid and recalculate.
Aerospace 2017, 4, 60 3 of 24 of non-planar wings through a configuration-invariant analytic formulation of the unknown circulation distribution.The authors present a numerical tool implementable in MATLAB, which allows the definition of non-planar wings' optimal circulation distribution with minimum induced drag.

Methodology and Tools for CFD and Structural Stress Simulation
The task of numerical simulation of air flow over the wing cantilever, and its stress-strain state under aerodynamic loads, was solved in ANSYS Workbench environment, allowing a combined simulation of the model using several modules from different disciplines, including a full analysis of the interaction of gas dynamics with the wing structure, and aeroelasticity phenomena.The project scheme, presented in Figure 1 shows the steps of a one-way data interchange between different modules: first, the initial solid 3D model is transmitted from the geometry CAD module to the aerodynamic calculations module (Fluent), where the mesh settings of airspace domain around the model are configured, and CFD calculations are performed.The solution, including air pressure field on the wing surface, is then transferred to the static structural strength calculation module, where the grid of the internal structure of the model is built, and its static stress-strain calculation is carried out.If any changes are applied to the original input geometry (a modified winglet twist angle for instance), subsequent modules automatically update the grid and recalculate.The continuity equation in the general form given above is valid for both incompressible (M ≤ 0.3), and compressible flows.Flow over the three-dimensional DLR-F4 wing-body configuration at a high subsonic velocity M = 0.6~0.85involves compressibility, therefore, an additional energy conservation equation is solved: .( ( )) .( ) ( . )

Governing Equations
For steady flow external aerodynamics, ANSYS Fluent solves conservation equations for massthe continuity equation: ∇.(ρ → v ) = 0, and for momentum [28]: The continuity equation in the general form given above is valid for both incompressible (M ≤ 0.3), and compressible flows.Flow over the three-dimensional DLR-F4 wing-body configuration at a high subsonic velocity M = 0.6~0.85involves compressibility, therefore, an additional energy conservation equation is solved: The complete system of governing equations and their underlying physical principles are given by Anderson [29].The compressible viscous flow density was treated by the ideal-gas law: ρ = P abs /(R/M w T), while the viscosity was calculated through Sutherland approximation: 2. General Description of DLR-F4 Simulated Model and Mesh Convergence Study DLR-F4 prototype, consisting of a fuselage in the form of a solid of revolution and a large aspect ratio swept back wing of a tapered planform, featuring a supercritical airfoil DFVLR R-4, allows the aerodynamic simulation of a typical medium-range subsonic passenger aircraft configuration.Using the DLR-F4 CAD model shown in Figure 2 as initial input, the influence of different winglet configurations on overall lift, drag and aerodynamic loads on the wing were studied.The impact of different cant angles for a classical Whitcomb winglet was identified in Mach number range of M = 0.6~0.87,with Reynolds' number being equal to Re = 3 × 10 6 .
Aerospace 2017, 4, 60 4 of 24 The complete system of governing equations and their underlying physical principles are given by Anderson [29].The compressible viscous flow density was treated by the ideal-gas law: /( / ) , while the viscosity was calculated through Sutherland approximation:  As a reference area for calculating the aerodynamic coefficients of lift and drag, the wing plan area was taken including its ventral (non-wetted) part S = 0.148 m 2 .The dimensions of the simulated domain and the mesh resolution were chosen to match the materials presented at AIAA CFD Drag Predictions (DP) seminars, for instance [30,31].The simulated domain dimensions were five body lengths along the symmetry axis of the model, 2.6 body lengths along the vertical axis and 2.5 body lengths along the z-axis.The mesh size varied from seven to 20 million cells.Standard k-ϖ turbulence model provided acceptable accuracy with reasonable computational cost.The standard k-ϖ turbulence model is founded on the developed by Wilox [32] modifications for low-Reynoldsnumber effects, compressibility, and shear flow spreading, where the following transport equations allow obtaining the turbulence kinetic energy, k: and the specific dissipation rate, ϖ: For CFD calculations on ANSYS FLUENT, an unstructured grid was used, consisting of tetrahedral cells condensing towards the surface of the wing-body, around which a structured submesh was inflated allowing for a suitable boundary layer resolution.Parameters of the sub-mesh were set as 15 layers with a growth factor of 1.5, so that for each 10% body length, the resolution of the boundary layer could be ensured by at least 20 grid cells.Typical fragments of the domain mesh used in gas dynamics calculations are shown in Figure 3 below.As a reference area for calculating the aerodynamic coefficients of lift and drag, the wing plan area was taken including its ventral (non-wetted) part S = 0.148 m 2 .The dimensions of the simulated domain and the mesh resolution were chosen to match the materials presented at AIAA CFD Drag Predictions (DP) seminars, for instance [30,31].The simulated domain dimensions were five body lengths along the symmetry axis of the model, 2.6 body lengths along the vertical axis and 2.5 body lengths along the z-axis.The mesh size varied from seven to 20 million cells.Standard kturbulence model provided acceptable accuracy with reasonable computational cost.The standard kturbulence model is founded on the developed by Wilox [32] modifications for low-Reynolds-number effects, compressibility, and shear flow spreading, where the following transport equations allow obtaining the turbulence kinetic energy, k: and the specific dissipation rate, : For CFD calculations on ANSYS FLUENT, an unstructured grid was used, consisting of tetrahedral cells condensing towards the surface of the wing-body, around which a structured sub-mesh was inflated allowing for a suitable boundary layer resolution.Parameters of the sub-mesh were set as 15 layers with a growth factor of 1.5, so that for each 10% body length, the resolution of the boundary layer could be ensured by at least 20 grid cells.Typical fragments of the domain mesh used in gas dynamics calculations are shown in Figure 3   The influence of the inflated boundary-layer-resolution structured submesh on lift and drag computational results is illustrated by the comparison with wind-tunnel experimental results in Figure 4, where its effect on drag prediction accuracy becomes evident at high angles of attack up to 4°.An example of the wall Y-plus distribution is given at Figure 5    The influence of the inflated boundary-layer-resolution structured submesh on lift and drag computational results is illustrated by the comparison with wind-tunnel experimental results in Figure 4, where its effect on drag prediction accuracy becomes evident at high angles of attack up to 4 • .An example of the wall Y-plus distribution is given at Figure 5 for the baseline model without winglets.The influence of the inflated boundary-layer-resolution structured submesh on lift and drag computational results is illustrated by the comparison with wind-tunnel experimental results in Figure 4, where its effect on drag prediction accuracy becomes evident at high angles of attack up to 4°.An example of the wall Y-plus distribution is given at Figure 5    The influence of the inflated boundary-layer-resolution structured submesh on lift and drag computational results is illustrated by the comparison with wind-tunnel experimental results in Figure 4, where its effect on drag prediction accuracy becomes evident at high angles of attack up to 4°.An example of the wall Y-plus distribution is given at Figure 5    The results (which are available at open source website of AGARD (Advisory Group for Aerospace Research and Development)) of wind tunnel experiments conducted during 1994 at DLR (Germany), ONERA (France), DRA (UK) and National Aerospace Laboratory of the Netherlands [30,31], allowed a verification of the robustness of our research technique, and validate the selected boundary conditions, mesh resolution, and turbulence models for the Mach and Reynold's numbers' ranges corresponding to typical climb and cruising flight mode conditions, a comparison of ANSYS Fluent computational results using different turbulence models with selected AGARD experimental results from [33,34] is given at Figure 6 below for the lift and drag coefficients at flight angles of attack.
Aerospace 2017, 4, 60 6 of 24 (Germany), ONERA (France), DRA (UK) and National Aerospace Laboratory of the Netherlands [30,31], allowed a verification of the robustness of our research technique, and validate the selected boundary conditions, mesh resolution, and turbulence models for the Mach and Reynold's numbers' ranges corresponding to typical climb and cruising flight mode conditions, a comparison of ANSYS Fluent computational results using different turbulence models with selected AGARD experimental results from [33,34] is given at Figure 6  To calculate the increments of aerodynamic loads on the wing in ANSYS Static Structural finite element solver, an unstructured mesh for the inner structure of the model was also built from tetrahedral cells, where a high resolution mesh was set only for the wing structure, being the main object for stress-strain analysis in this paper, with proximity and curvature refinement.While for the fuselage body, a coarse grid was constructed and the boundary condition "fixed" was set, allowing zero-displacements in space under external loads.This simulates the weight of an aircraft fuselage, and allowed us to obtain aerodynamic wing loading values close to those loads encountered during a steady level flight.For static stress analysis, the total number of cells of the inner mesh of the wing and body combined is ~1 million.An example of a typical grid used in stress calculations is shown in Figure 7.There are virtually no open source winglet prototypes available for a known transport airplane configuration with a known optimal geometry, or known empirical equations relating the winglet efficiency to its geometrical parameters.Thus, for an initial assessment of the influence on overall To calculate the increments of aerodynamic loads on the wing in ANSYS Static Structural finite element solver, an unstructured mesh for the inner structure of the model was also built from tetrahedral cells, where a high resolution mesh was set only for the wing structure, being the main object for stress-strain analysis in this paper, with proximity and curvature refinement.While for the fuselage body, a coarse grid was constructed and the boundary condition "fixed" was set, allowing zero-displacements in space under external loads.This simulates the weight of an aircraft fuselage, and allowed us to obtain aerodynamic wing loading values close to those loads encountered during a steady level flight.For static stress analysis, the total number of cells of the inner mesh of the wing and body combined is ~1 million.An example of a typical grid used in stress calculations is shown in Figure 7.
Aerospace 2017, 4, 60 6 of 24 (Germany), ONERA (France), DRA (UK) and National Aerospace Laboratory of the Netherlands [30,31], allowed a verification of the robustness of our research technique, and validate the selected boundary conditions, mesh resolution, and turbulence models for the Mach and Reynold's numbers' ranges corresponding to typical climb and cruising flight mode conditions, a comparison of ANSYS Fluent computational results using different turbulence models with selected AGARD experimental results from [33,34] is given at Figure 6  To calculate the increments of aerodynamic loads on the wing in ANSYS Static Structural finite element solver, an unstructured mesh for the inner structure of the model was also built from tetrahedral cells, where a high resolution mesh was set only for the wing structure, being the main object for stress-strain analysis in this paper, with proximity and curvature refinement.While for the fuselage body, a coarse grid was constructed and the boundary condition "fixed" was set, allowing zero-displacements in space under external loads.This simulates the weight of an aircraft fuselage, and allowed us to obtain aerodynamic wing loading values close to those loads encountered during a steady level flight.For static stress analysis, the total number of cells of the inner mesh of the wing and body combined is ~1 million.An example of a typical grid used in stress calculations is shown in Figure 7.There are virtually no open source winglet prototypes available for a known transport airplane configuration with a known optimal geometry, or known empirical equations relating the winglet efficiency to its geometrical parameters.Thus, for an initial assessment of the influence on overall There are virtually no open source winglet prototypes available for a known transport airplane configuration with a known optimal geometry, or known empirical equations relating the winglet efficiency to its geometrical parameters.Thus, for an initial assessment of the influence on overall DLR-F4 aerodynamics of various geometric dimensions and setting angles of the winglet, and for their optimal choice in further work, a brief exploratory study of the design space of both a triangular wingtip fence and a classical winglet prototype was carried out, and a mathematical model was developed.

Design Space Exploration for Different Wingtip Device Geometries
In design optimization theory, the notion of design space of a system denotes the totality of all the values of geometrical dimensions and setting angles that are subject to free adjustment.A degree of freedom for a wingtip device geometry corresponds to a geometrical parameter; changing within certain limits affects the overall aerodynamic efficiency of the wing + wingtip device system.The growth in the lift-to-drag ratio L/D at cruising flight conditions as compared to the initial model without winglets was chosen as the main criterion for winglet geometry analysis and optimization.

Wingtip Fence Design Space
For a two-dimensional delta planform wingtip fence (Figure 8b), based on theoretical considerations and wind tunnel experiments performed at [35,36], the key to L/D ratio growth of a wing equipped with a thin-airfoil tip fence is a large relative area of the wingtip fence S fence /S wing .While the remaining dimensions, the height and aspect ratio, have relatively little effect on its efficiency.Thus, we have a one-dimensional design space, where L/D ratio growth is directly proportional to the wingtip fence relative area S fence /S wing .However, S fence upper values are limited by losses due to viscous friction of large fences (S fence /S wing ≥ 0.3), that can well exceed the gains in induced drag.Further increasing the fence area leads to a smaller growth, and eventually to a fall in the L/D ratio (Figure 8a, solid line).CFD calculations of the aerodynamic coefficients of DLR-F4 model equipped with few oval planform (for simplicity) endplates, aimed at finding near optimal values of a wingtip fence plan area, showed that the optimal endplate relative area lies in the range: [S fence /S wing ] optimal ~0.25 to 0.3.These results were taken into account for subsequent CFD and structural stress investigations of this near-optimal delta planform wingtip fence.
Aerospace 2017, 4, 60 7 of 24 DLR-F4 aerodynamics of various geometric dimensions and setting angles of the winglet, and for their optimal choice in further work, a brief exploratory study of the design space of both a triangular wingtip fence and a classical winglet prototype was carried out, and a mathematical model was developed.

Design Space Exploration for Different Wingtip Device Geometries
In design optimization theory, the notion of design space of a system denotes the totality of all the values of geometrical dimensions and setting angles that are subject to free adjustment.A degree of freedom for a wingtip device geometry corresponds to a geometrical parameter; changing within certain limits affects the overall aerodynamic efficiency of the wing + wingtip device system.The growth in the lift-to-drag ratio L/D at cruising flight conditions as compared to the initial model without winglets was chosen as the main criterion for winglet geometry analysis and optimization.

Wingtip Fence Design Space
For a two-dimensional delta planform wingtip fence (Figure 8b), based on theoretical considerations and wind tunnel experiments performed at [35,36], the key to L/D ratio growth of a wing equipped with a thin-airfoil tip fence is a large relative area of the wingtip fence Sfence/Swing.While the remaining dimensions, the height and aspect ratio, have relatively little effect on its efficiency.Thus, we have a one-dimensional design space, where L/D ratio growth is directly proportional to the wingtip fence relative area Sfence/Swing.However, Sfence upper values are limited by losses due to viscous friction of large fences (Sfence/Swing ≥ 0.3), that can well exceed the gains in induced drag.Further increasing the fence area leads to a smaller growth, and eventually to a fall in the L/D ratio (Figure 8a, solid line).CFD calculations of the aerodynamic coefficients of DLR-F4 model equipped with few oval planform (for simplicity) endplates, aimed at finding near optimal values of a wingtip fence plan area, showed that the optimal endplate relative area lies in the range: [Sfence/Swing] optimal ~0.25 to 0.3.These results were taken into account for subsequent CFD and structural stress investigations of this near-optimal delta planform wingtip fence.

Winglet Design Space
A classical winglet is unlike the 'flat' wingtip fence shown above.Whitcomb winglet is a complex three-dimensional lifting surface with at least six geometrical parameters: sweep λ, cant ѱ and twist ξ angles, the height h, fillet radius R and the taper ratio B = b/b0 (Figure 9).In the coordinate system

Winglet Design Space
A classical winglet is unlike the 'flat' wingtip fence shown above.Whitcomb winglet is a complex three-dimensional lifting surface with at least six geometrical parameters: sweep λ, cant ψ and twist ξ angles, the height h, fillet radius R and the taper ratio B = b/b 0 (Figure 9).In the coordinate system z-axis denotes the wing span direction perpendicular to the aircraft (XY) symmetry plane, x-axis: the direction of flight and y-axis normal to XZ plane: the normal force direction.With over six degrees of freedom of the winglet design space, it is not viable to explore the whole design space by manual variations of each dimensional combinations.In order to best explore the vast multitude of possible combinations of input parameters, automated optimization simulations were carried out using ANSYS Design Exploration tool on relatively coarse grids for computational cost economy.In Figure 10 below, we show an example of surfaces estimating the effect on aerodynamic efficiency of winglet-equipped DLR-F4, of combined variations of both the winglet height (m) with twist (rad) and fillet radius (m) with taper ratio (dimensionless).Despite relatively low accuracy coarse grids, response surfaces in Figure 10 allowed us to unveil the optimal envelope for values of the winglet geometrical parameters.This in turn allowed us to narrow the winglet six-dimensional design space to only one-dimension, where for further highaccuracy CFD investigations and loads or structural stress evaluation, it is sufficient to independently only vary parameters that mostly influence the L/D ratio and wing root bending moment, while keeping the values of other dimensions fixed.As for airfoiled winglets, it was shown that L/D ratio, as well as the wing root bending moment and as a result, span wise stress magnitude and distribution, are mostly sensitive to the winglet cant angle ѱ (see Sections 3.1.1and 3.2).It is worth noting that using response surfaces, it is possible to optimize the winglet geometry for both maximum L/D ratio, and minimum wing root bending moment at the same time varying all six variables at once, but the computational cost would be unreasonably expensive, compromising the fidelity of With over six degrees of freedom of the winglet design space, it is not viable to explore the whole design space by manual variations of each dimensional combinations.In order to best explore the vast multitude of possible combinations of input parameters, automated optimization simulations were carried out using ANSYS Design Exploration tool on relatively coarse grids for computational cost economy.In Figure 10 below, we show an example of surfaces estimating the effect on aerodynamic efficiency of winglet-equipped DLR-F4, of combined variations of both the winglet height (m) with twist (rad) and fillet radius (m) with taper ratio (dimensionless).With over six degrees of freedom of the winglet design space, it is not viable to explore the whole design space by manual variations of each dimensional combinations.In order to best explore the vast multitude of possible combinations of input parameters, automated optimization simulations were carried out using ANSYS Design Exploration tool on relatively coarse grids for computational cost economy.In Figure 10 below, we show an example of surfaces estimating the effect on aerodynamic efficiency of winglet-equipped DLR-F4, of combined variations of both the winglet height (m) with twist (rad) and fillet radius (m) with taper ratio (dimensionless).Despite relatively low accuracy coarse grids, response surfaces in Figure 10 allowed us to unveil the optimal envelope for values of the winglet geometrical parameters.This in turn allowed us to narrow the winglet six-dimensional design space to only one-dimension, where for further highaccuracy CFD investigations and loads or structural stress evaluation, it is sufficient to independently only vary parameters that mostly influence the L/D ratio and wing root bending moment, while keeping the values of other dimensions fixed.As for airfoiled winglets, it was shown that L/D ratio, as well as the wing root bending moment and as a result, span wise stress magnitude and distribution, are mostly sensitive to the winglet cant angle ѱ (see Sections 3.1.1and 3.2).It is worth noting that using response surfaces, it is possible to optimize the winglet geometry for both maximum L/D ratio, and minimum wing root bending moment at the same time varying all six variables at Despite relatively low accuracy coarse grids, response surfaces in Figure 10 allowed us to unveil the optimal envelope for values of the winglet geometrical parameters.This in turn allowed us to narrow the winglet six-dimensional design space to only one-dimension, where for further high-accuracy CFD investigations and loads or structural stress evaluation, it is sufficient to independently only vary parameters that mostly influence the L/D ratio and wing root bending moment, while keeping the values of other dimensions fixed.As for airfoiled winglets, it was shown that L/D ratio, as well as the wing root bending moment and as a result, span wise stress magnitude and distribution, are mostly sensitive to the winglet cant angle ψ (see Sections 3.1.1and 3.2).It is worth noting that using response surfaces, it is possible to optimize the winglet geometry for both maximum L/D ratio, and minimum wing root bending moment at the same time varying all six variables at once, but the computational cost would be unreasonably expensive, compromising the fidelity of results.Thus, a high fidelity cant angle ψ CFD and mathematical optimization was prioritized over other geometrical parameters to get a cost effective, 'multi-fidelity' approach.

Mathematical Modeling of the Winglet Local Flow Field
Given the winglet is an airfoiled finite-span lifting surface, finite wing theory can be applied to describe the three-dimensional flow field in its vicinity.In this paper, the winglet surface local flow field has been quantitatively analyzed through the term of the winglet local angle of attack α winglet , which is different in each section of the winglet (z) and generally assumed to be a function of both the winglet geometry, mainly the cant ψ and twist ξ winglet angles, and the general angle of attack α wing of DLR-F4 model, as well as three-dimensional flow field variables: Here, the first and the second terms are geometric, while the third term arises from the three-dimensional flow field over the winglet span.In general: • ξ winglet (z): the winglet twist, a function of z coordinate (changes spanwise): ξ winglet (0) = 0 means no twist at the winglet root.By convention, an outward twist at some z n coordinate, which leads to an increased local α winglet (z n ) ≥ α winglet (0), is considered positive: ξ winglet (z n ) ≥ 0. • K ψ : Cant coefficient accounting for the winglet angle of attack α winglet sensitivity to changing the general angle of attack of the model (or the wing) α = α wing , which is linked to the cant ψ as follows: omitting both ξ winglet (z) and α f low_ f ield (definition below), K ψ can be found from (1) and the assumption that increasing the angle of attack of a wing equipped with a vertical winglet (ψ = 0, see Figure 9) does not increase the angle of attack of the winglet itself α vertical winglet = const (as a vertical winglet would be rotating over z axis in its own plane, experiencing a greater slip instead).From (1): K ψ=0 = 0.While for a horizontal winglet (ψ = π/2), the winglet angle of attack is directly equal to the wing angle of attack ∆α horizontal winglet ≈ ∆α wing (a horizontal winglet is just an extension of the wing itself), from (1) this means K ψ=π = 1.Solving the system of equations: Only positively canted winglets were considered in this paper, where ψ is ranged from 0 ≤ ψ ≤ π 2 , which implies: 0 ≤ K ψ ≤ 1. (2) can be verified for a half-canted winglet (ψ = π/4): π/2 = 0.5.This means that for any increase of the general angle of attack, half of that increase will contribute to geometrically increasing the winglet angle of attack: ∆α ψ=π/4 winglet = 0.5∆α wing .

•
α f low_ f ield accounts for three-dimensional flow field near the winglet, involving a strong interference with the main wingtip vortex, which tends to increase local α root winglet at the winglet root sections (this can be clearly seen in Figure 11, where a substantially lower pressure on top winglet surface near the leading edge close to the wing junction).On the other hand, the winglet's own 'downwash' and its smaller own tip vortex, lead to a slight decrease in local α tip winglet at the winglet tip sections.The interaction between the tip vortex and angle of attack has been studied in [37,38], and a method for predicting local root and tip variations, using local induced velocities is given in [39].For a fixed winglet geometry at a fixed α wing , α f low_ f ield is a function of the location on the winglet span (z coordinate): Substituting ( 2) and ( 3) in ( 1), a non-twisted (ξ winglet (z) = 0) winglet local angle of attack would be: Finding the exact α f low_ f ield (z) function is beyond the scope of this paper, but its general shape has been approximately estimated for the non-twisted winglet shown in Figure 11, featuring a cant of ψ = π 4 , CFD-simulated at Mach number M = 0.8 and α wing = 2 • ~0.0349 rad.Substituting the values of ψ and α wing in (4), the local α winglet for this winglet depends on the location z, in the same fashion as α f low_ f ield (z) does, with a spanwise constant geometrical half α 'cant correction' of 0.0174: From ( 5) we can assume that the α winglet (z) and α f low_ f ield (z) should have the same graph shape.Computationally, this shape can be qualitatively obtained using the winglet surface pressure field and velocity vectors as shown in Figure 11 below.An optical method for experimental wind tunnel measurement of the local angle of attack is described in [40].
Aerospace 2017, 4, 60 10 of 24 induced velocities is given in [39].For a fixed winglet geometry at a fixed wing α , shape.Computationally, this shape can be qualitatively obtained using the winglet surface pressure field and velocity vectors as shown in Figure 11 below.An optical method for experimental wind tunnel measurement of the local angle of attack is described in [40].The mathematical model and Equation (4) described in this item will be particularly useful in interpreting the results of CFD simulation of few non-twisted winglets with different cant angles, and subsequently in the task of optimizing the cant angle for sustaining L/D ratio growth of DLR-F4 + winglet model at high angles of attack.The mathematical model and Equation (4) described in this item will be particularly useful in interpreting the results of CFD simulation of few non-twisted winglets with different cant angles, and subsequently in the task of optimizing the cant angle for sustaining L/D ratio growth of DLR-F4 + winglet model at high angles of attack.without winglets, we can confirm the positive effect of this wingtip device on DLR-F4 L/D ratio at flight angles of attack (up to ~4% gain).Winglets with larger cant angles generally perform better, due to higher lifting characteristics excepted at high and close to critical angles of attack where nearly vertical winglets (ψ = 25 • ) slightly prevail, due to a better induced drag reduction (acting as a vertical fence better isolating upper and lower wing tips with largely different pressure values helps weaken wingtip vortex and hence induced drag).As shown in Figure 12 below, dependence on the angle of attack of overall growth in L/D ratio (%) of the model equipped with winglets of different cant angles, as compared with initial DLR-F4 without winglets, we can confirm the positive effect of this wingtip device on DLR-F4 L/D ratio at flight angles of attack (up to ~4% gain).Winglets with larger cant angles generally perform better, due to higher lifting characteristics excepted at high and close to critical angles of attack where nearly vertical winglets (ѱ = 25°) slightly prevail, due to a better induced drag reduction (acting as a vertical fence better isolating upper and lower wing tips with largely different pressure values helps weaken wingtip vortex and hence induced drag).

Impact of the Winglet Cant Angle on DLR-F4 Lift-to-Drag Ratio
As shown in Figure 12 below, dependence on the angle of attack of overall growth in L/D ratio (%) of the model equipped with winglets of different cant angles, as compared with initial DLR-F4 without winglets, we can confirm the positive effect of this wingtip device on DLR-F4 L/D ratio at flight angles of attack (up to ~4% gain).Winglets with larger cant angles generally perform better, due to higher lifting characteristics excepted at high and close to critical angles of attack where nearly vertical winglets (ѱ = 25°) slightly prevail, due to a better induced drag reduction (acting as a vertical fence better isolating upper and lower wing tips with largely different pressure values helps weaken wingtip vortex and hence induced drag).From Figures 12 and 13, a major problem of winglet design is the proper choice of a balanced cant angle that would prevent the flow separation and loss of efficiency in the winglet tip area at high α during take-off and climb, while providing acceptable lift at moderate α during cruise.Geometric twist is a well-known solution in correcting local α of lifting surfaces.A negative twist of the highly canted winglet's tip sections, for instance, might be applied to geometrically reduce the tip local angle

Winglet Cant Optimization for a Sustainable Efficiency along α Range Using a Curved Span
From Figures 12 and 13, a major problem of winglet design is the proper choice of a balanced cant angle that would prevent the flow separation and loss of efficiency in the winglet tip area at high α during take-off and climb, while providing acceptable lift at moderate α during cruise.Geometric twist is a well-known solution in correcting local α of lifting surfaces.A negative twist of the highly canted winglet's tip sections, for instance, might be applied to geometrically reduce the tip local angle of attack and obtain a smoother pressure distribution, thus delaying flow separation at high α.However, exploratory studies performed at Figure 10 on relatively coarse grids revealed that twist correction is advantageous only at substantially high angles of attack.At smaller cruise angles of attack, a negatively twisted tip local angle of attack becomes too small, or even negative, leading to a substantial loss in L/D ratio.Thus, a twist-morphing concept might be considered, dynamically applying twist only at high α.Another, more efficient morphing is dynamic cant correction: shifting to smaller cant angles at higher angles of attack, thus dynamically reducing the cant coefficient K ψ (see Section 2.3.2,Equations ( 1) and ( 2)), making the winglet local angle of attack less sensitive to increasing the general angle of attack, and using the winglet as a vertical fence at high α.This can be achieved through a cant-morphing winglet.Several studies regarding morphing concepts were performed, see [11][12][13]26] for instance.In this paper, a much simpler and structurally lighter alternative is suggested that does not involve morphing, and therefore, no complex articulation mechanisms are needed: A curved winglet, where local cant angle is different at each section (highly canted root versus nearly vertical tip, with gradual transition in between, see Figure 14 below).

Winglet Cant Optimization for a Sustainable Efficiency along α Range Using a Curved Span
From Figures 12 and 13, a major problem of winglet design is the proper choice of a balanced cant angle that would prevent the flow separation and loss of efficiency in the winglet tip area at high α during take-off and climb, while providing acceptable lift at moderate α during cruise.Geometric twist is a well-known solution in correcting local α of lifting surfaces.A negative twist of the highly canted winglet's tip sections, for instance, might be applied to geometrically reduce the tip local angle of attack and obtain a smoother pressure distribution, thus delaying flow separation at high α.However, exploratory studies performed at Figure 10 on relatively coarse grids revealed that twist correction is advantageous only at substantially high angles of attack.At smaller cruise angles of attack, a negatively twisted tip local angle of attack becomes too small, or even negative, leading to a substantial loss in L/D ratio.Thus, a twist-morphing concept might be considered, dynamically applying twist only at high α.Another, more efficient morphing is dynamic cant correction: shifting to smaller cant angles at higher angles of attack, thus dynamically reducing the cant coefficient ψ K (see Section 2.3.2,Equations ( 1) and ( 2)), making the winglet local angle of attack less sensitive to increasing the general angle of attack, and using the winglet as a vertical fence at high α.This can be achieved through a cant-morphing winglet.Several studies regarding morphing concepts were performed, see [11][12][13]26] for instance.In this paper, a much simpler and structurally lighter alternative is suggested that does not involve morphing, and therefore, no complex articulation mechanisms are needed: A curved winglet, where local cant angle is different at each section (highly canted root versus nearly vertical tip, with gradual transition in between, see Figure 14 below).The curved winglet concept can be mathematically analyzed as follows: The curved winglet centerline front (YZ) projection should be parametrized using a second order function Y(z), where the graph is tangent to the wing at the wing-winglet joint.For instance, for a parabolic winglet as Figure 14.The concept of curved winglet provides an alternative to a much more complex and structurally heavier variable-cant, morphing winglet.
The curved winglet concept can be mathematically analyzed as follows: The curved winglet centerline front (YZ) projection should be parametrized using a second order function Y(z), where the graph is tangent to the wing at the wing-winglet joint.For instance, for a parabolic winglet as shown on Figure 15, the parametrization function is Y(z) = z 2 .For a curved winglet, the constant cant of a straight winglet ψ should be replaced by a cant function ψ(z), which gives the local cant spanwise at each point (From Figure 15, this is π/2 minus the slope of the tangent line of the parametrization function graph).Therefore, the cant function ψ(z) can be calculated by deriving the chosen parametrization function as follows: For the parabolic winglet case, the cant function is: z], which can be verified at the winglet root: ψ(0) = π 2 − arctan(0) = π 2 = ψ max , which means that the local cant at the winglet root is the maximum possible.This complies with wing-winglet tangential connection requirement.
For the parabolic winglet case, the cant function is: which can be verified at the winglet root:  CFD simulation of the parabolic winglet confirmed its reduced sensitivity to high angles of attack, where the gain in L/D ratio remains the highest up until α = 5°.The solid graph at Figure 16 below, corresponding to the parabolic winglet shows that a curved winglet provides a compromise, performing nearly as good as highly canted winglets at small α, with a sustainable L/D growth at higher angles of attack.CFD simulation of the parabolic winglet confirmed its reduced sensitivity to high angles of attack, where the gain in L/D ratio remains the highest up until α = 5 • .The solid graph at Figure 16 below, corresponding to the parabolic winglet shows that a curved winglet provides a compromise, performing nearly as good as highly canted winglets at small α, with a sustainable L/D growth at higher angles of attack.
For the parabolic winglet case, the cant function is: which can be verified at the winglet root:  CFD simulation of the parabolic winglet confirmed its reduced sensitivity to high angles of attack, where the gain in L/D ratio remains the highest up until α = 5°.The solid graph at Figure 16 below, corresponding to the parabolic winglet shows that a curved winglet provides a compromise, performing nearly as good as highly canted winglets at small α, with a sustainable L/D growth at higher angles of attack.Flow field visualization in Figure 17 clearly reveals the positive effect of the winglet curved span shape on the three-dimensional flow field, which now involves a much weaker tip vortex and a delayed flow separation.Pressure field on the top surface of the parabolic winglet, near the leading edge, reveals how the tip sections, being nearly vertical allow a better local angle of attack distribution spanwise, thus resulting in a smoother pressure distribution and a sustainable efficiency of the winglet tip sections.

Optimal Cant Function
For an arbitrarily shaped curved-span winglet, by replacing the constant cant ѱ with a cant function Substituting ( 6) in (7), we obtain the winglet local angle of attack dependency on z for a nontwisted, arbitrarily shaped curved winglet, parametrized through an arbitrary second order function Y(z): As was mentioned earlier, the winglet shape should be optimized for getting the highest possible, at all flight angles of attack, L/D ratio of DLR-F4 + winglet model.From (8), seeking the optimal shape of a curved winglet is seeking the optimal parametrization function Y(z) optim , the derivative of which Y'(z) optim , being substituted in (8) results in an optimal spanwise distribution of local

Optimal Cant Function
For an arbitrarily shaped curved-span winglet, by replacing the constant cant ψ with a cant function ψ(z), (4) becomes: Substituting ( 6) in (7), we obtain the winglet local angle of attack dependency on z for a non-twisted, arbitrarily shaped curved winglet, parametrized through an arbitrary second order function Y(z): which, given arctan[Y (z)] in rad, reduces to: As was mentioned earlier, the winglet shape should be optimized for getting the highest possible, at all flight angles of attack, L/D ratio of DLR-F4 + winglet model.From (8), seeking the optimal shape of a curved winglet is seeking the optimal parametrization function Y(z) optim , the derivative of which Y'(z) optim , being substituted in (8) results in an optimal spanwise distribution of local α winglet (z).
In the light of CFD results shown at Figures 12 and 13, Figures 16 and 17 and the definition of the suggested spanwise-curved winglet solution (Section 3.1.2),this optimal α winglet spanwise distribution can be achieved by ensuring, at the same time, the following three conflicting local winglet flow field conditions:

•
Ensuring the highest possible, yet under critical, local angle of attack of the winglet: for a curved winglet design, this is interpreted as the highest cant possible at root, in order to best harness the main wingtip vortex energy (see Figure 11) and use it to increase local α root winglet .Also, in order to get the highest possible local increase, should the wing angle of attack increases ∆α root winglet ≈ ∆α wing (see Section 2.3.2,Equation ( 2)).This highest cant ψ root max = π/2 is achieved through a tangent connection to the wing at z root = 0: ψ opt (0) = π 2 .Substituting in (6): • Ensuring the smallest possible pressure difference between upper and lower surfaces at the winglet tip, thus keeping the winglet tip vortex and its own downwash at the lowest possible intensity.This can be achieved by making the tip vertical (with zero cant), in order to prevent its angle of attack from increasing (and the pressure difference from rising), should the general angle of attack increases: ∆α tip winglet ≈ 0 (see Section 2.3.2).Hence ψ opt (z tip ) = 0, Substituting in (6): (10) means an infinite slope (vertical) tangent line to Y(z) optim.graph at the at the winglet tip.

•
Ensuring the most advantageous transition from α root winglet to α tip winglet .From (8), this requires a well-known flow field dependency α f low_ f ield (z); defining it is beyond the scope of this paper, see [37,39].We can only notice that α f low_ f ield (z) itself depends on the winglet geometry, meaning that the optimization process is iterative, starting from an arbitrary second order parametrization function, defining α f low_ f ield (z) for it (either mathematically or experimentally), identifying an improvement direction, and then gradually adapting the parametrization function up until the point where L/D ratio does not grow anymore.It is well-known, however, that optimal spanwise circulation distribution of a lifting surface with the minimum induced drag is elliptical.Therefore, we may assume that a similar elliptic cant parametrization which easily complies with (9) and (10), would result in the highest L/D ratio of the winglet.
Given the boundary conditions from ( 9) and (10), as well as a mathematically optimal elliptic transition, an optimal curved winglet may be parametrized by a quarter ellipse function (Figure 18): where: A and B correspond to the ellipse semi-major and semi-minor axes respectively.For a known maximum wing span, which is usually the case for an operational aircraft, it is easy to define the semi-major axis A. For illustration, let's assume that an operational DLR-F4 aircraft maximum allowable wing semi-span is L DLR−F4 max = 78.56cm.Initial DLR-F4 wing semi-span without winglet is L DLR−F4 without ≈ 58.56 cm, then A = L DLR−F4 max − L DLR−F4 without = 20 cm.The semi-minor axis B can be obtained from a fixed winglet height h = B. From response surfaces in Figure 10 let's choose: h = B = 0.125 m = 12.5 cm.
Substituting A and B in (11): The graph of this particular-case elliptic winglet is given at Figure 18 below.It is worth noting that the concept presented herein of a spanwise-curved winglet is not completely new, and can be viewed as an extrapolation of the (already known to be efficient, see [9,22] for details) blended winglet concept, featuring an increased fillet radius and tangent wing joint.Airbus's latest wide body airliner A350 features a similar curved winglet, while the Boeing 787 wing has its tip section elastically deformed upward in flight, resulting in a gradual transition to a lower local cant angle of its raked wingtips the closer we move to the tips [41].
The novelty of this research lies in the methodology used for quantitative analysis of the effect of various geometric parameters of the winglet, the cant angle in particular, through the analysis of their effect on the local angle of attack, which both geometrically and physically drastically varies from one section of the winglet to the other.This is important for winglet and non-planar wings design, because in comparison with the 3D flow field around a high aspect ratio wing such as DLR-F4 wing (which, far from the tips, can be replaced by a simpler 2D flow with constant-spanwise local α), the flow field around the winglet is deeply three-dimensional, dominated by a high intensity tip vortex from the wing, and the winglet's own downwash.
local cant angle of its raked wingtips the closer we move to the tips [41].
The novelty of this research lies in the methodology used for quantitative analysis of the effect of various geometric parameters of the winglet, the cant angle in particular, through the analysis of their effect on the local angle of attack, which both geometrically and physically drastically varies from one section of the winglet to the other.This is important for winglet and non-planar wings design, because in comparison with the 3D flow field around a high aspect ratio wing such as DLR-F4 wing (which, far from the tips, can be replaced by a simpler 2D flow with constant-spanwise local α), the flow field around the winglet is deeply three-dimensional, dominated by a high intensity tip vortex from the wing, and the winglet's own downwash.Figure 18.An example of an optimal elliptic parametrization of a curved-span winglet with tangent junction to the wing at its root, a vertical tangent line at the tip and an elliptic transition in between.

Stress-Strain State Simulation Results of Winglet-Equipped DLR-F4
Besides L/D ratio, another important result of CFD calculation is the air pressure distribution on the wing surface, which serves as the initial external loads input information for stress calculations on the wing cantilever.Physical experiments made by R. Whitcomb back in the 1970s [1,2], included an estimation of the growth of wing loading by integrating air pressure obtained from a large number of sensors over the wing and winglet upper and lower surfaces.In the course of computational experiments in this paper, a similar approach to sensors is used, where air pressure is read by the software for each of the individual surface mesh cells (Figure 19), and then integrated to calculate the aerodynamic loads.Figure 18.An example of an optimal elliptic parametrization of a curved-span winglet with tangent junction to the wing at its root, a vertical tangent line at the tip and an elliptic transition in between.

Stress-Strain State Simulation Results of Winglet-Equipped DLR-F4
Besides L/D ratio, another important result of CFD calculation is the air pressure distribution on the wing surface, which serves as the initial external loads input information for stress calculations on the wing cantilever.Physical experiments made by R. Whitcomb back in the 1970s [1,2], included an estimation of the growth of wing loading by integrating air pressure obtained from a large number of sensors over the wing and winglet upper and lower surfaces.In the course of computational experiments in this paper, a similar approach to sensors is used, where air pressure is read by the software for each of the individual surface mesh cells (Figure 19), and then integrated to calculate the aerodynamic loads.In Figure 20, a comparison of the wing root bending moment (WRBM) values and maximum equivalent Von Mises stress for winglets of cant angles: =25°, 45° and 75° is shown for a parabolic winglet, a delta planform wingtip fence and the initial DLR-F4 wing without wingtip device.Simulation was performed at positive angles of attack and Mach number of M = 0.75.From Figure 20a, we can conclude that the increase of wing root bending moment for winglets is directly proportional to the value of their cant angle.This can be explained by a greater additional lifting force, the highly canted the winglet is.Thus, the smallest increase in wing loading corresponds to the 'non-lifting' wingtip fence.At moderate angles of attack α~1°, all winglets cause virtually the same WRBM increment, but increasing the angle of attack, the bending moment increases at a slower pace the less canted the winglet is.This is a direct implication of the cant coefficient K described in In Figure 20, a comparison of the wing root bending moment (WRBM) values and maximum equivalent Von Mises stress for winglets of cant angles: ψ =25 • , 45 • and 75 • is shown for a parabolic winglet, a delta planform wingtip fence and the initial DLR-F4 wing without wingtip device.Simulation was performed at positive angles of attack and Mach number of M = 0.75.From Figure 20a, we can conclude that the increase of wing root bending moment for winglets is directly proportional to the value of their cant angle.This can be explained by a greater additional lifting force, the highly canted the winglet is.Thus, the smallest increase in wing loading corresponds to the 'non-lifting' wingtip fence.At moderate angles of attack α~1 • , all winglets cause virtually the same WRBM increment, but increasing the angle of attack, the bending moment increases at a slower pace the less canted the winglet is.This is a direct implication of the cant coefficient K ψ described in Section 2.3.2,Equations ( 1) and ( 2), which is greater for highly canted winglets, leading to a faster-pace growth in their angle of attack and subsequently greater lifting force acting on them, which in turn leads to a faster pace growth of WRBM.In this regard, the parabolic winglet, although causing relatively high WRBM values, offers a compromise solution up until α~3 • .At α~4 • and higher, WRBM growth becomes higher than for conventional straight winglets.This can be explained by the sustainable lift generated by its tip sections (compare top surface pressures at both the winglets' tips at Figure 17).
In Figure 20, a comparison of the wing root bending moment (WRBM) values and maximum equivalent Von Mises stress for winglets of cant angles: =25°, 45° and 75° is shown for a parabolic winglet, a delta planform wingtip fence and the initial DLR-F4 wing without wingtip device.Simulation was performed at positive angles of attack and Mach number of M = 0.75.From Figure 20a, we can conclude that the increase of wing root bending moment for winglets is directly proportional to the value of their cant angle.This can be explained by a greater additional lifting force, the highly canted the winglet is.Thus, the smallest increase in wing loading corresponds to the 'non-lifting' wingtip fence.At moderate angles of attack α~1°, all winglets cause virtually the same WRBM increment, but increasing the angle of attack, the bending moment increases at a slower pace the less canted the winglet is.This is a direct implication of the cant coefficient ψ K described in Section 2.3.2,Equations ( 1) and ( 2), which is greater for highly canted winglets, leading to a fasterpace growth in their angle of attack and subsequently greater lifting force acting on them, which in turn leads to a faster pace growth of WRBM.In this regard, the parabolic winglet, although causing relatively high WRBM values, offers a compromise solution up until α~3°.At α~4° and higher, WRBM growth becomes higher than for conventional straight winglets.This can be explained by the sustainable lift generated by its tip sections (compare top surface pressures at both the winglets' tips at Figure 17).WRBM is a central factor of wingtip device weight penalty, but the combined effect of all winglet-induced aerodynamic loadings increases in bending, torque moment and the shear force, can be more fully studied by evaluating the changes of both the value and distribution of structural stress along the wing span.For this purpose, structural steel was set as the material of a one-piece wing cantilever structure.As the study presented here is mainly for a comparison purpose, structural material choice and internal structural layout is not important so far.On Figure 21 below, an example is shown of 34% growth in maximum Von Mises stress due to the installation of a highly canted winglet with ψ = 75 • , at M = 0.8 and α = 1 • .According to experience of wing structural weight estimation for transport aircraft, a minimum resulted in wing box structural weight growth of 10% for an Al-alloy box structure, and at least 15% for a spared structure [42].
From Figure 20, we can assume that the optimal for structural weight penalty wingtip device is the wingtip fence.Thus, it can be retrofitted into most operational aircraft almost without the need for a wing structure rework.It should be noted, however, that ensuring the aerodynamically optimal fence relative area S fence ~0.27S wing (based on the dependence shown earlier in Figure 8a), results in considerable local stress and deformations of the wingtip fence itself (Figure 22).This is due to the large plan area combined with an extremely small thickness (structural depth).Counter measures to locally strengthen the wingtip fence and its junction area with appropriate material choice are required, and they may lead to significant weight penalties.cantilever structure.As the study presented here is mainly for a comparison purpose, structural material choice and internal structural layout is not important so far.On Figure 21 below, an example is shown of 34% growth in maximum Von Mises stress due to the installation of a highly canted winglet with ѱ = 75°, at М = 0.8 and α = 1°.According to experience of wing structural weight estimation for transport aircraft, a minimum resulted in wing box structural weight growth of 10% for an Al-alloy box structure, and at least 15% for a spared structure [42].From Figure 20, we can assume that the optimal for structural weight penalty wingtip device is the wingtip fence.Thus, it can be retrofitted into most operational aircraft almost without the need for a wing structure rework.It should be noted, however, that ensuring the aerodynamically optimal fence relative area Sfence~0.27Swing(based on the dependence shown earlier in Figure 8a), results in considerable local stress and deformations of the wingtip fence itself (Figure 22).This is due to the large plan area combined with an extremely small thickness (structural depth).Counter measures to locally strengthen the wingtip fence and its junction area with appropriate material choice are required, and they may lead to significant weight penalties.

Developing a Multidisciplinary Criterion for Wingtip Device Aerodyamic Efficiency Assessment Regarding Structural Weight Penalty
Taking into account the conflicting nature of the influence of different wingtip devices on aircraft fuel efficiency, where we have a perceptible improvement in Lift-to-Drag ratio (up to 4-5%) at flight angles of attack and cruising Mach numbers on one hand, and a considerable increase in wing loading, involving a structural weight penalty on the other hand (Figures 20-22).Also, given the vast design space and the multitude of choices available for a device geometry, decisions about wingtip device configuration made during the early stages of wing design should be guided by some practical  From Figure 20, we can assume that the optimal for structural weight penalty wingtip device is the wingtip fence.Thus, it can be retrofitted into most operational aircraft almost without the need for a wing structure rework.It should be noted, however, that ensuring the aerodynamically optimal fence relative area Sfence~0.27Swing(based on the dependence shown earlier in Figure 8a), results in considerable local stress and deformations of the wingtip fence itself (Figure 22).This is due to the large plan area combined with an extremely small thickness (structural depth).Counter measures to locally strengthen the wingtip fence and its junction area with appropriate material choice are required, and they may lead to significant weight penalties.

Developing a Multidisciplinary Criterion for Wingtip Device Aerodyamic Efficiency Assessment Regarding Structural Weight Penalty
Taking into account the conflicting nature of the influence of different wingtip devices on aircraft fuel efficiency, where we have a perceptible improvement in Lift-to-Drag ratio (up to 4-5%) at flight angles of attack and cruising Mach numbers on one hand, and a considerable increase in wing loading, involving a structural weight penalty on the other hand (Figures 20-22).Also, given the vast design space and the multitude of choices available for a device geometry, decisions about wingtip device configuration made during the early stages of wing design should be guided by some practical

Developing a Multidisciplinary Criterion for Wingtip Device Aerodyamic Efficiency Assessment Regarding Structural Weight Penalty
Taking into account the conflicting nature of the influence of different wingtip devices on aircraft fuel efficiency, where we have a perceptible improvement in Lift-to-Drag ratio (up to 4-5%) at flight angles of attack and cruising Mach numbers on one hand, and a considerable increase in wing loading, involving a structural weight penalty on the other hand (Figures 20-22).Also, given the vast design space and the multitude of choices available for a device geometry, decisions about wingtip device configuration made during the early stages of wing design should be guided by some practical design philosophy allowing to reach an acceptable level of wingtip device aerodynamic efficiency, keeping the weight penalty in mind.In this paper, a quantitative evaluation criterion is suggested, that measures the aerodynamic benefits of wingtip devices at cruising flight mode, responsible for most of the timespan of aircraft operation and where, as a result, even small improvements in aerodynamic efficiency yield significant long term savings for airlines.As a test condition for wingtip device weight penalty though, it is proposed to measure loads for mostly loaded flight modes, including off-design conditions.Since the wing structural strength is usually designed to withstand maximum, despite only occasionally occurring (or may never occur) aerodynamic loads, besides typical cruising loads.
Thus, for quantitative comparison of the aerodynamic benefits of different wingtip devices, regarding their structural cost, it is proposed to calculate and compare a dimensionless ratio denoted W Tip_device of the increase in L/D ratio during cruising flight ∆(L/D) Cruise.(%) to the maximum possible (chosen from an off-design condition) growth of structural stress ∆δ max (%): Cruising flight mode is selected within the flight envelope to ensure the most economical flight conditions possible.The required thrust to ensure a steady level flight, P req_Level is equal to the aircraft weight divided by the L/D ratio [43]: Consequently, given a relatively constant aircraft weight G Cruise Aircra f t. = const, from (13) the most economical flight mode with the minimum required thrust P Cruise req_Level = min, is achieved at the maximum lift-to-drag ratio (L/D) cruise = max.For illustration on the studied DLR-F4 model, this is achieved at M = 0.6 and α = 3 • , where (L/D) max_DLR−F4 ≈ 20.Therefore, evaluating the aerodynamic efficiency of each winglet configuration through cruising lift-to-drag ratio growth ∆(L/D) Cruise_DLR−F4 (%) should be performed for this, presumably cruising, air flow regime.However, the most stressed (off-design for an operational aircraft) case was achieved at the highest value of Mach number M max = 0.8 and angle of attack α max = 5 • experimented in this research.Hence, we used this flow regime to seek the maximum possible wingtip-device-induced growth of wing structural stress (%) and we obtained the criterion W Tip_device for each case.In Figure 23    Based on the bar chart presented above, we can see the relative efficiency of the suggested curved winglet concept which provides a good compromise between classical Whitcomb winglets with small cant angles or wingtip fences, having minimum weight penalties and moderate aerodynamic efficiency, and 'lifting' winglets of larger cant angles giving maximum aerodynamic L/D ratio growth at cruise, but are much less efficient at high α during take-off and climb, and resulting in a significant structural stress growth and weight increments.

Discussion of Results and Conclusions
In this study, we presented a solution to a multidisciplinary optimization problem, concerning the interaction of sub-and near transonic aerodynamics with the cantilever of DLR-F4 wing-body prototype, equipped with wingtip devices of different geometries.The solution algorithm begins with low-cost CFD design space exploration on relatively coarse grids, which allows us to evaluate the general aerodynamic efficiency and viability of using different wingtip devices, as well as their most critical geometrical parameters.In the light of CFD results and using simple geometrical analysis, a mathematical model of the local flow field around the winglet and immediately downstream was developed using the local angle of attack function along the winglet span.High fidelity CFD simulations and mathematical analysis are then launched to optimize the chosen critical geometrical parameter (cant angle for instance), while keeping other parameters fixed.Besides lift and drag, CFD simulations yield air pressure distribution over the wing surface, which is used by the stress calculation module as an external load input for numerical modeling of the stress-strain state of the wing cantilever structure.Based on the results of evaluating the growth of lift-to-drag ratio of DLR-F4 after equipping it with different wingtip devices, and its juxtaposing with the increase in structural stress, a tool has been developed for estimating and comparing the overall efficiency through calculating the L/D ratio increment at cruising flight regime, and dividing it by the maximum structural stress increment, encountered during most dangerous off-design flight modes which the wing structure is, nonetheless, supposed to withstand.Using this algorithm, the wingtip device geometry can be optimized in a multidisciplinary fashion at the earliest stages of wing design.As for Whitcomb winglet configurations described in this paper, it has been revealed that the most balanced choice would be the suggested concept of a curved (nearly elliptical) winglet, which provides a compromise solution for cant angle choice.
Future research may enrich the used herein mathematical model of the winglet local flow field by analyzing the impact on its local flow field of other geometrical parameters such as the sweep angle, which may lead to a similar optimal elliptic-leading-edge compromise solution.Combined with a curved span, this curved leading edge would result in a double-curvature winglet, where a three-dimensional parametrization function, depicting a spatial double-curved winglet centerline, is required to be deduced and thoroughly optimized.
The winglet efficiency estimation factor should be extended to include the winglets' effect on directional stability, by quantifying possible reductions to the vertical tail plan area and its structural weight, as well as by a quantified analysis of winglets' impact on the wing aeroelasticity and the overall aero-acoustic and environmental footprint of airliners in the long run.In this regard, it is advisable to use, as an experimental platform, a more realistic transport aircraft prototype as DLR-F6 or -F11, that includes what is typical for passenger aircraft high lift devices and engine nacelles.
Acknowledgments: This work was supported by the Ministry of Education and Science of the Russian Federation, project No. 9.7170.2017/8.9.Special thanks to the reviewers for their valuable advices and comments.
Author Contributions: Djahid Gueraiche performed the experiments as part of a Ph.D. thesis, wrote the manuscript, developed the mathematical model and conducted the analysis of results.Sergey Popov supervised the work and helped to guide the choice of boundary conditions and mesh resolution, as well as to enhance the mathematical rigor of the mathematical model, and provided quality feedback and valuable advices.

Conflicts of Interest:
The authors declare no conflict of interest.The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.Aircraft take-off weight

Figure 1 .
Figure 1.Data exchange scheme about DLR-F4 simulated model between the geometry, gas dynamics and structural stress modules inside ANSYS Workbench environment.

Figure 1 .
Figure 1.Data exchange scheme about DLR-F4 simulated model between the geometry, gas dynamics and structural stress modules inside ANSYS Workbench environment.
Description of DLR-F4 Simulated Model and Mesh Convergence Study DLR-F4 prototype, consisting of a fuselage in the form of a solid of revolution and a large aspect ratio swept back wing of a tapered planform, featuring a supercritical airfoil DFVLR R-4, allows the aerodynamic simulation of a typical medium-range subsonic passenger aircraft configuration.Using the DLR-F4 CAD model shown in Figure 2 as initial input, the influence of different winglet configurations on overall lift, drag and aerodynamic loads on the wing were studied.The impact of different cant angles for a classical Whitcomb winglet was identified in Mach number range of M = 0.6~0.87,with Reynolds' number being equal to

Figure 3 .
Figure 3. Fragments of the fluid domain unstructured mesh used for aerodynamic calculations with an inflated on the surface structured submesh allowing for a fine boundary layer resolution.
for the baseline model without winglets.

Figure 4 .Figure 5 .
Figure 4. Comparison of computed DLR-F4 lift (a) and drag (b) coefficients with AGARD wind tunnel experiment results at M = 0.6 and Re = 3 × 10 6 using a mesh with an inflated boundary layer resolution submesh and without it.(AGARD: Advisory Group for Aerospace Research and Development)

Figure 3 .
Figure 3. Fragments of the fluid domain unstructured mesh used for aerodynamic calculations with an inflated on the surface structured submesh allowing for a fine boundary layer resolution.

Figure 3 .
Figure 3. Fragments of the fluid domain unstructured mesh used for aerodynamic calculations with an inflated on the surface structured submesh allowing for a fine boundary layer resolution.
for the baseline model without winglets.

Figure 4 .Figure 5 .
Figure 4. Comparison of computed DLR-F4 lift (a) and drag (b) coefficients with AGARD wind tunnel experiment results at M = 0.6 and Re = 3 × 10 6 using a mesh with an inflated boundary layer resolution submesh and without it.(AGARD: Advisory Group for Aerospace Research and Development)

Figure 4 . 24 Figure 3 .
Figure 4. Comparison of computed DLR-F4 lift (a) and drag (b) coefficients with AGARD wind tunnel experiment results at M = 0.6 and Re = 3 × 10 6 using a mesh with an inflated boundary layer resolution submesh and without it.(AGARD: Advisory Group for Aerospace Research and Development) for the baseline model without winglets.

Figure 4 .Figure 5 .
Figure 4. Comparison of computed DLR-F4 lift (a) and drag (b) coefficients with AGARD wind tunnel experiment results at M = 0.6 and Re = 3 × 10 6 using a mesh with an inflated boundary layer resolution submesh and without it.(AGARD: Advisory Group for Aerospace Research and Development)

Figure 5 .
Figure 5.An example of Wall Y-plus parameter distribution over the surface of DLR-F4 wing and body.

Figure 6 .
Figure 6.(a) Comparison of computed DLR-F4 lift coefficient for different turbulence models with wind tunnel experiment results at M = 0.75 and Re = 3 × 10 6 ; (b) Comparison of computed DLR-F4 drag coefficient with experimental results at M = 0.75 and Re = 3 × 10 6 .

Figure 7 .
Figure 7. Internal DLR-F4 (shown equipped with a curved winglet) unstructured mesh used to estimate the wing loads, stress and displacements.

Figure 6 .
Figure 6.(a) Comparison of computed DLR-F4 lift coefficient for different turbulence models with wind tunnel experiment results at M = 0.75 and Re = 3 × 10 6 ; (b) Comparison of computed DLR-F4 drag coefficient with experimental results at M = 0.75 and Re = 3 × 10 6 .

Figure 6 .
Figure 6.(a) Comparison of computed DLR-F4 lift coefficient for different turbulence models with wind tunnel experiment results at M = 0.75 and Re = 3 × 10 6 ; (b) Comparison of computed DLR-F4 drag coefficient with experimental results at M = 0.75 and Re = 3 × 10 6 .

Figure 7 .
Figure 7. Internal DLR-F4 (shown equipped with a curved winglet) unstructured mesh used to estimate the wing loads, stress and displacements.

Figure 7 .
Figure 7. Internal DLR-F4 (shown equipped with a curved winglet) unstructured mesh used to estimate the wing loads, stress and displacements.

Figure 8 .
Figure 8.(a) Influence of oval planform end plates relative area on L/D ratio growth (%) of DLR-F4 model at M = 0.75 and α = 3°.Solid line graph reveals the negative effect of viscous friction for fences of large areas; (b) Surface pressure field and streamlines over a delta planform wingtip fence at M = 0.8 and α = 1°.

Figure 8 .
Figure 8.(a) Influence of oval planform end plates relative area on L/D ratio growth (%) of DLR-F4 model at M = 0.75 and α = 3 • .Solid line graph reveals the negative effect of viscous friction for fences of large areas; (b) Surface pressure field and streamlines over a delta planform wingtip fence at M = 0.8 and α = 1 • .
Aerospace 2017, 4, 60 8 of 24 z-axis denotes the wing span direction perpendicular to the aircraft (XY) symmetry plane, x-axis: the direction of flight and y-axis normal to XZ plane: the normal force direction.
Aerospace 2017, 4, 60 8 of 24 z-axis denotes the wing span direction perpendicular to the aircraft (XY) symmetry plane, x-axis: the direction of flight and y-axis normal to XZ plane: the normal force direction.
location on the winglet span (z coordinate): the scope of this paper, but its general shape has been approximately estimated for the non-twisted winglet shown in Figure11, featuring a cant of 4 π ψ = , CFD-simulated at Mach number M = 0.8 and °= 2 wing α ~0.0349 rad.Substituting the values of ψ and wing α in (4), the local winglet α for this winglet depends on the location z, a spanwise constant geometrical half α 'cant correction' of 0.0174:

Figure 11 .
Figure 11.Three-dimensional flow field over a non-twisted, half canted winglet with ѱ = π/4 and its local angle of attack approximate dependency on the winglet span location.CFD results shown at M = 0.8, α = 2°.

Figure 11 .
Figure 11.Three-dimensional flow field over a non-twisted, half canted winglet with ψ = π/4 and its local angle of attack approximate dependency on the winglet span location.CFD results shown at M = 0.8, α = 2 • .
Winglet-Equipped DLR-F4 Model 3.1.1.Impact of the Winglet Cant Angle on DLR-F4 Lift-to-Drag Ratio As shown in Figure 12 below, dependence on the angle of attack of overall growth in L/D ratio (%) of the model equipped with winglets of different cant angles, as compared with initial DLR-F4 Winglet-Equipped DLR-F4 Model 3.1.1.Impact of the Winglet Cant Angle on DLR-F4 Lift-to-Drag Ratio

Figure 12 .
Figure 12.Dependence of overall improvement (%) in L/D ratio of DLR-F4 model on DLR-F4 geometric angle of attack due to the use of winglets of different cant angles.M = 0.75.

12 .
Dependence of overall improvement (%) in L/D ratio of DLR-F4 model on DLR-F4 geometric angle of attack due to the use of winglets of different cant angles.M = 0.75.From Figure 12, the influence of the winglet cant on DLR-F4 L/D ratio changes depending on the angle of attack of the model, α wing = α.Highly canted winglets are more sensitive to increasing the angle of attack of the wing than nearly vertical winglets.This sensitivity is illustrated in Figure 13 below, where winglets with different cant angles are exposed to a high-subsonic air flow regime M = 0.8 at a high angle of attack.Streamlines reveal a flow separation, an intense tip vortex and a loss of efficiency for the winglet with ψ = 75 • and, to a lesser extent the one with ψ = 45 • .While the flow field behind the least canted winglet with only ψ = 25 • , remains relatively smooth despite high angle of attack of the model.Aerospace 2017, 4, 60 11 of 24

Figure 12 . 24 Figure 13 .
Figure 12.Dependence of overall improvement (%) in L/D ratio of DLR-F4 model on DLR-F4 geometric angle of attack due to the use of winglets of different cant angles.M = 0.75.From Figure 12, the influence of the winglet cant on DLR-F4 L/D ratio changes depending on the angle of attack of the model, α α = wing

Figure 13 .
Figure 13.Pressure field and streamlines over winglets with different cant at a high angle of attack α = 4 • and M = 0.8, Velocity vectors in section planes reveal attached flow at the winglet root versus a flow separation and a visible tip vortex occurring at tip sections of highly canted winglets.

Figure 13 .
Figure 13.Pressure field and streamlines over winglets with different cant at a high angle of attack α = 4° and M = 0.8, Velocity vectors in section planes reveal attached flow at the winglet root versus a flow separation and a visible tip vortex occurring at tip sections of highly canted winglets.

Figure 14 .
Figure 14.The concept of curved winglet provides an alternative to a much more complex and structurally heavier variable-cant, morphing winglet.
that the local cant at the winglet root is the maximum possible.This complies with wing-winglet tangential connection requirement.

Figure 15 .
Figure 15.An example of a parabolic parametrization of a curved-span winglet.The tangent line slope at each point allows us to define the local winglet cant ) ( n z ψ

Figure 16 .
Figure 16.Comparison of overall improvement (%) in L/D ratio of DLR-F4 model equipped with a curved parabolic winglet with conventional winglets at M = 0.75 and different geometric angles of attack of the model.

Figure 15 .
Figure 15.An example of a parabolic parametrization of a curved-span winglet.The tangent line slope at each point allows us to define the local winglet cant ψ(z n ) which in turn affects the local angle of attack.
that the local cant at the winglet root is the maximum possible.This complies with wing-winglet tangential connection requirement.

Figure 15 .
Figure 15.An example of a parabolic parametrization of a curved-span winglet.The tangent line slope at each point allows us to define the local winglet cant ) ( n z ψ

Figure 16 .
Figure 16.Comparison of overall improvement (%) in L/D ratio of DLR-F4 model equipped with a curved parabolic winglet with conventional winglets at M = 0.75 and different geometric angles of attack of the model.

Figure 16 .
Figure 16.Comparison of overall improvement (%) in L/D ratio of DLR-F4 model equipped with a curved parabolic winglet with conventional winglets at M = 0.75 and different geometric angles of attack of the model.

.
In the light of CFD results shown at Figures 12, 13, 16 and 17 and the definition of the suggested spanwise-curved winglet solution (Section 3.1.2),this optimal winglet α spanwise distribution can be achieved by ensuring, at the same time, the following three conflicting local winglet flow field conditions:

Figure 19 .
Figure 19.An example of spanwise air pressure distribution (Pascal) over DLR-F4 wing surface.Surface pressure field serves as the basic information for aerodynamic loads estimation.

Figure 19 .
Figure 19.An example of spanwise air pressure distribution (Pascal) over DLR-F4 wing surface.Surface pressure field serves as the basic information for aerodynamic loads estimation.

Figure 20 .
Figure 20.(a) DLR-F4 wing root bending moment (WRBM) values for different-configuration winglets and wingtip fence; (b) Values of the maximum equivalent Von Mises stress depending on the geometric angle of attack of DLR-F4 model equipped with different wingtip devices.

Figure 20 .
Figure 20.(a) DLR-F4 wing root bending moment (WRBM) values for different-configuration winglets and wingtip fence; (b) Values of the maximum equivalent Von Mises stress depending on the geometric angle of attack of DLR-F4 model equipped with different wingtip devices.

Figure 22 .
Figure 22.Distribution of wing spanwise stress after the installation of a 'flat' delta planform wingtip fence: despite low-stressed wingspan, significant local stress and deformations can be observed on the wingtip fence itself and its junction area.Air flow regime-high subsonic: М = 0.8, α = 1°.

Figure 22 .
Figure 22.Distribution of wing spanwise stress after the installation of a 'flat' delta planform wingtip fence: despite low-stressed wingspan, significant local stress and deformations can be observed on the wingtip fence itself and its junction area.Air flow regime-high subsonic: М = 0.8, α = 1°.

Figure 22 .
Figure 22.Distribution of wing spanwise stress after the installation of a 'flat' delta planform wingtip fence: despite low-stressed wingspan, significant local stress and deformations can be observed on the wingtip fence itself and its junction area.Air flow regime-high subsonic: M = 0.8, α = 1 • .

.
below, shown is a bar chart comparing the values of W Tip_device for the studied wingtip devices. is selected within the flight envelope to ensure the most economical flight conditions possible.The required thrust to ensure a steady level flight, Level req Р _ is equal to the aircraft weight divided by the L/D ratio [43]: For illustration on the studied DLR-F4 model, this is achieved at М = 0.6 and α = 3°, where 20 evaluating the aerodynamic efficiency of each winglet configuration through cruising lift-to-drag ratio growth should be performed for this, presumably cruising, air flow regime.However, the most stressed (off-design for an operational aircraft) case was achieved at the highest value of Mach number Мmax = 0.8 and angle of attack αmax = 5° experimented in this research.Hence, we used this flow regime to seek the maximum possible wingtip-device-induced growth of wing structural stress (%) and we obtained the criterion device Tip W _ for each case.In Figure 23 below, shown is a bar chart comparing the values of device Tip W _ for the studied wingtip devices.

Figure 23 .
Figure 23.Comparison of values of the dimensionless criterion WTip_device for different wingtip devices installed on DLR-F4, revealing the relative aerodynamics-to-weight efficiency for each of the devices.

Figure 23 .
Figure 23.Comparison of values of the dimensionless criterion W Tip_device for different wingtip devices installed on DLR-F4, revealing the relative aerodynamics-to-weight efficiency for each of the devices.
Effective conductivity (k + k t where k t is the turbulent thermal conductivity, defined according to the turbulence model being used) Drag ratio, C l /C d S fence Wingtip fence plan area S wing DLR-F4 wing plan area, including its ventral non-wetted part [S fence /S wing ] optim Optimal wingtip fence relative area, yielding maximum L/D ratio of (LR-F4 + fence) system Winglet effective AOA, including its geometric angle of attack and α f low_ f ield α f low_ f ield Winglet AOA increment due to three-dimensional flow field over the winglet span Winglet cant function, depicting local cant angle at each z coordinate spanwise ψ opt (z) Winglet optimal cant function, depicting local ψ in the case of an optimal winglet shape, parametrized through Y(z) optim.A Semi-major axis of an elliptic winglet span curve B Semi-minor axis of an elliptic winglet span curve W Tip_device Dimensionless wingtip device overall efficiency criterion ∆(L/D) Cruise.DLR-F4 L/D ratio growth at cruise conditions, due to the wingtip device ∆δ max Maximum structural stress growth on DLR-F4 wing due to wingtip device R req_Level Required thrust for level flight G Aircra f t.