Numerical analysis of an elasto-flexible wing configuration using computational fluid structural simulations

356 VII International Conference on Textile Composites and Inflatable Structures STRUCTURAL MEMBRANES 2015 E. Oñate, K.-U.Bletzinger and B. Kröplin (Eds) Numerical analysis of an elasto-flexible wing configuration using computational fluid structural simulations VI International Conference on Textile Composites and Inflatable Structures STRUCTURAL MEMBRANES 2015 E. Oñate, K.-U.Bletzinger, and B. Kröplin (Eds) NUMERICAL ANALYSIS OF AN ELASTO-FLEXIBLE WING CONFIGURATION USING COMPUTATIONAL FLUID STRUCTURAL SIMULATIONS J. PIQUEE*AND C. BREITSAMTER† * Institute of Aerodynamics and Fluid Mechanics, Technische Universität München, 85748 Garching, Germany e-mail: [email protected], web page: https://www.aer.mw.tum.de † Institute of Aerodynamics and Fluid Mechanics, Technische Universität München, 85748 Garching, Germany e-mail: [email protected], web page: https://www.aer.mw.tum.de Key words: Elasto-flexible wing, Coupled RANS and FEM Simulations, Coupled Panel method and analytical membrane code, SST model, γ-Reθ transition model Summary: In the present work, different numerical approaches dealing with fluid and structural interaction models are introduced to simulate the aerodynamics and structural behavior of 2D elasto-flexible wings. In order to verify the models, the results are compared with the experimental data based on wind tunnel investigations conducted at the Institute of Aerodynamics and Fluid Mechanics of the Technische Universität München. The results show a good agreement while the findings are to some extent similar to sailwing configurations: the lift is higher than on the rigid counterpart geometry whereas the drag is nearly equivalent and the stall region is smoother and appears delayed. 1 INTRODUCTION With an average growth rate of nearly 30% each year over the past five years1, energy produced by wind has been the fastest growing source of energy in the world today. With global warming and rising fuel prices, the wind energy industry will still significantly play an important role in the future. However, it still remains important to improve the technology in order to keep economically wind energy competitive. Increasing the rotor diameter and therefore the turbine size could offer a more efficient way to produce energy with wind turbines: a larger turbine can indeed capture more energy throughout its lifetime. However larger turbines means higher structural and fatigue loads2, which can cause more possible damages and then more unscheduled maintenance. The control of loads of a turbine has then led to intensive research in order to keep the wind energy domain competitive. Among others Basualdo3 numerically investigated a rear-mounted flap concerning 357 J. Piquee and C. Breitsamter. 2 its impact on loads. The study was conducted using a potential theory method for a variablegeometry airfoil by analyzing the changes in the airfoil displacements. The displacements affect directly the fatigue of the material through the bending moment. It was found that the use of a flapped airfoil can lead to load alleviation by reducing the airfoil displacement variations. Therefore, a variable geometry airfoil appears as a possible solution for the wind energy domain to increase performance. Such a design also suggests an alleviation of the effects of gust winds, which can result in a reduction of loads and fatigues. For Micro Air Vehicles, an increasing amount of resources has already been spent to investigate vehicles able to alter their shape like wings featuring a flexible membrane, with the main goal of improving their efficiency and reducing loads4,5,6. The aerodynamic force measurements show that aerodynamics characteristics of elasto-flexible membrane airfoils/wings are comparable or even better than their rigid counterparts7,8. Due to the flexibility of the airfoil/wing, its shape can adapt itself to the free stream flow by changing its geometry, which will induce new aerodynamics and structural properties. This adaptation to the flow results also in a delayed stall region: it was experimentally shown that flexible membrane airfoils/wings stall at higher angle of attack. For example, typical rigid wing for Micro Air Vehicles have stall angles between 12° and 15°, whereas flexible wings, which reduce their effective angle of attack to the surface deformation, may have stall angles up to 30° and 45°5. This paper focuses on a 2D elasto-flexible airfoil and has the purpose to present a numerical investigation of such a configuration. The numerical work was made by developing a coupling between fluid and structure solvers in order to simulate the strong relation between the mechanical and aerodynamic properties. The paper will present two different models to simulate the behavior of the 2D elasto-flexible airfoil: on one hand a Panel-method code coupled with an analytical membrane description and on the other hand a CFD-U/RANS simulation coupled with a FEM program at low Reynolds number (0.23x106). 2 MODEL GEOMETRY AND NUMERICAL METHODS The geometry used in this work, which is illustrated in Fig. 1, has been developed within an experimental study9: the corresponding investigation considers an elasto-flexible membrane wrapped around rigid leading- and trailing-edge spars. The leading-edge spar is made of a double ellipse to reduce the gradient of pressure which facilitates flow separation9,10, whereas the trailing-edge spar was designed as a sharp element for 2D experiments. The membrane material used for the wing surface was chosen for its compliance to adapt itself to the flow free stream and for its resistance as it has to be strong enough to bear the aerodynamic loads. In the experiments9, the membrane is a highly extensible, anisotropic elastic fabric coat covered with an impermeable rubber layer. As the membrane is an anisotropic material, its mechanical properties are dependent on all the directions and are linked with each other: the nonlinear deformation of the membrane depends on the ratio of the applied loads in the warp and in the weft directions11. However, for some reasons of simplicity, it was decided to consider the material as an isotropic material for both numerical models with a Young modulus equal as the Young modulus in the weft direction with E=2.1MPa. The pre-stress of the membrane was also 358 J. Piquee and C. Breitsamter. 3 taken into account as it controls significantly the deformation of the membrane: the initial elongation of the membrane is equal to 2% in the chord direction9. Figure 1 – Geometry of the elasto-flexible wing configuration 2.1 Panel method coupled with an analytical membrane model The first method to simulate such a configuration was realized with the software Matlab by developing a Panel method code coupled with an analytical membrane model. Based on a discretization of the body, this method consists on finding the pressure distribution around the geometry divided into panels and then determining the deformation of the elasto-flexible membrane defined within the analytical membrane behavior. As the Panel method is a well-known method, a detailed description of the model will not be given here, but for more information the reference12 can be used. Nevertheless, a brief description of the membrane code will be given; for more details, the references10,13 can be used. The membrane code is based on the theoretical membrane deflection equation which is for a 2D case described by Eq. (1) 𝑡𝑡𝑥𝑥 𝜕𝜕²𝑧𝑧(𝑥𝑥) 𝜕𝜕𝑥𝑥²= −(𝑃𝑃(𝑥𝑥)− 𝑃𝑃 ∞), (1) where the parameter z(x) represents the deflection of the membrane at the position (x), 𝑃𝑃(𝑥𝑥)− 𝑃𝑃 ∞ the pressure difference at the position (x) and tx the tension component defined as: 𝑡𝑡𝑥𝑥= 𝑡𝑡𝑥𝑥𝑥𝑥 + 𝐾𝐾𝑥𝑥∆𝑙𝑙𝑥𝑥. (2) The parameter Kx is defined as material property and ∆lx is the elongation of the membrane in the x-direction. Eq. (1) represents the relationship between the curvature of the membrane, the tension and the difference of pressure which will induce the deformation. The membrane is considered to be fixed at the leading-edge (𝑥𝑥 = 0) and at the trailing-edge (𝑥𝑥 = 𝑐𝑐). The introduction of the point xs, defined as the point where the membrane will separate from the contour of the leading-edge spar, is then important to find the new shape of the airfoil. Therefore, between 𝑥𝑥 = 0 and 𝑥𝑥 ≤ xs, the airfoil geometry is given by the leading-edge geometry, whereas for x > xs, the shape is determined by the deflection of the membrane which is given by solving Eq. (1). The boundary conditions of Eq. (1) are given in the following: 𝑧𝑧(𝑥𝑥𝑠𝑠)= 𝑧𝑧𝐿𝐿(𝑥𝑥𝑠𝑠) (3) -0.1 0 0.1 x/c y/c Fig.1 - Initial Geometry of the Sailwing 00.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 359 J. Piquee and C. Breitsamter. 4 𝑧𝑧′(𝑥𝑥𝑠𝑠)= 𝑧𝑧𝐿𝐿′(𝑥𝑥𝑠𝑠) (4) 𝑧𝑧(𝑐𝑐)= 0. (5) Eq. (3) assures that the membrane contour is the same as the leading-edge contour up to station xs, then Eq. (4) assures that the separation between the membrane and the leading-edge will appear smoothly and finally Eq. (5) illustrates the fixing condition of the membrane at the trailing-edge. The solution of Eq. (1) is obtained by integrating Eq.(1) two times as it is shown by Eq. (6), obtaining an equation for the parameter xs: 𝑧𝑧(𝑥𝑥)=1 𝑡𝑡𝑥𝑥∫ ∫ ∆𝑃𝑃(𝑥𝑥).𝑑𝑑𝑥𝑥.𝑑𝑑𝑥𝑥 𝑥𝑥 𝑥𝑥𝑠𝑠 𝑥𝑥 𝑥𝑥𝑠𝑠 + 𝑑𝑑𝑧𝑧𝐿𝐿(𝑥𝑥) 𝑑𝑑𝑥𝑥 𝑥𝑥𝑠𝑠 . 𝑥𝑥 + (𝑧𝑧𝐿𝐿(𝑥𝑥𝑠𝑠)−𝑑𝑑𝑧𝑧𝐿𝐿(𝑥𝑥) 𝑑𝑑𝑥𝑥 𝑥𝑥𝑠𝑠 .𝑥𝑥𝑠𝑠), (6) with 𝑡𝑡𝑥𝑥=∫ ∫ ∆𝑃𝑃(𝑥𝑥).𝑑𝑑𝑥𝑥.𝑑𝑑𝑥𝑥 𝑐𝑐 𝑥𝑥𝑠𝑠 𝑐𝑐 𝑥𝑥𝑠𝑠 (𝑐𝑐−𝑥𝑥𝑠𝑠).𝑑𝑑𝑧𝑧𝐿𝐿(𝑥𝑥) 𝑑𝑑𝑥𝑥 𝑥𝑥𝑠𝑠+𝑧𝑧𝐿𝐿(𝑥𝑥𝑠𝑠) . (7) To entirely solve the system, xs needs to be found: the method is to express the elongation of the membrane in two different ways. The point xs where those expressions are equal (Eq. (8)) will give the final solution. 𝑙𝑙(𝑥𝑥𝑠𝑠)−𝑙𝑙𝑜𝑜 𝑙𝑙𝑜𝑜=𝑡𝑡𝑥𝑥 𝐸𝐸𝑚𝑚− 𝜀𝜀𝑜𝑜, (8) with lo as initial length of the membrane, Em the Young Modulus of the membrane and εo the initial pre-elongation. Figure 2 – S812 Airfoil Geometry In order to validate the Panel code developed with Matlab, it was decided to compare the pressure distribution along the airfoil S812 (Fig. 2) and the lift coefficient for different angles of attack with those obtained with the software Xfoil and with experimental data15. At an angle of attack of α=2°, different numbers of panels were tested (Tab. 1) to check the influence of the discretization on the results. Fig. 3 shows the evolution of the lift coefficient obtained at α=2° for different numbers of panels which vary from 200 to 5000 around the S812 profile. As it is expected, the lift coefficient value increases with the number of panels up to a certain value. In this case, the final value of the lift coefficient obtained at α=2° with 5000 panels is equal to 0.6028. The software Xfoil provides 0.6025 with 400 panels and within the experiment a lift 0.1 -0.1 x/c y/c Fig.2 - S812 Airfoil Geometry 00.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 360 J. Piquee and C. Breitsamter. 5 Figure 5 – Pressure coefficient distribution at 2° for the S812 profile Table 1 – Lift coefficient comparison between the Panel Method code and Xfoil for different angles of attack -0.2 0 0.2 0.4 0.6 0.8 11.2 -1 -0.5 0 0.5 1 1.5 x/c -Cp 200 Panels 500 Panels 2000 Panels Xfoil Experimental Data 0.95 0.96 0.97 0.98 0.99 11.01 1.02 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 x/c -Cp 200 Panels 500 Panels 2000 Panels Xfoil Experimental Data Figure 4 – Pressure coefficient distribution at 2° for the S812 profile coefficient of 0.49 at a Reynolds number of 106 is obtained. Considering Xfoil as a reference, the error of the Panel code is then 0.05%. In Figs. 4 and 5 the pressure distribution is plotted for the S812 profile at 2° obtained with the Panel method for 200, 500 and 2000 panels as well as with Xfoil and by the experiments. The pressure distribution along the profile is wellapproximated with the Panel method when it is compared with the distribution obtained with Xfoil. Nevertheless, a difference is noticeable at the trailing edge, where a reverse pressure and a peak can be observed with the Panel method code. However the reverse pressure disappears whereas the peak augments when the number of panels increases. These unphysical phenomena is caused by the interpolation of the profile which is based on a cubic spline interpolation between the given grid points of the profile. Considering the experimental data, the Panel code overestimates the pressure along the profile which explains the difference between the lift coefficient values at α=2°, presented in Tab. 1. For the following, it was decided to work with 2000 panels as the error (compared with Xfoil) obtained at α=2° for the S812 profile is already below 0.5%. The lift coefficient from α=-2° to 8° for the S812 profile obtained with Xfoil and the Panel code is listed in Tab. 1. The relative error remains below 3%. α/Lift Coefficient Panel Method α=2° Xfoil Error (%) Experimental Data 0° 0.350 0.353 0.91 0.26 2° 0.6 0.603 0.45 0.49 4° 0.849 0.851 0.26 0.71 6° 1.097 1.098 0.13 0.94 8° 1.343 1.344 0.04 1.03 200500 1000 1500 2000 2500 3000 4000 5000 0.56 0.57 0.58 0.59 0.6 Panels number Cl for the airfoil S812 at 2° Fig.3 - Evolution of the Lift coefficient with the number of panels Figure 3 – Evolution of the lift coefficient with the number of panels 0.95 0.96 0.97 0.98 0.99 11.01 1.02 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 x/c -Cp 200 Panels 500 Panels 2000 Panels Xfoil Experimental Data 361 J. Piquee and C. Breitsamter. 6 2.2 Fluid structure simulations performed with CFD ans FEM The second method concentrates on Fluid-Structure-Interaction simulations based on a coupling between the fluid solver ANSYS CFX (CFD-U/RANS method) and the structural solver STATIC STRUCTURAL (FEM method), directly controlled by the Workbench of ANSYS. As the structural deformations of the airfoil are significant to change the fluid flow itself, the so called 2-ways interaction simulations were performed using a mesh motion inside the fluid solver. The displacement of the mesh is controlled by a mesh stiffness parameter, which controls which regions deform to absorb the motion and which remain relatively rigid. The mesh motion is applied near the small volumes as the deformation of the fluid grid appears inside the boundary layer of the airfoil, which is composed of the thinnest nodes.  Fluid model The fluid solver software used in this work is ANSYS CFX. As the considered Reynolds number is small, namely Re = 0.23x106, according to the Mayle correlation14, the percentage of laminar flow relative to the chord length is estimated to be more than 5%. As the laminarturbulent transition phenomenon affects the pressure distribution, the deformation and the behavior of the elasto-flexible wing, two models concerning the fluid simulations were investigated within ANSYS CFX: the k-ω based SST turbulence model and the k-ω based SST coupled with the transition γ-Reθ model. Furthermore, as the Mach number is low (𝑀𝑀𝑀𝑀 < 0.3) the fluid is considered as incompressible and it is also supposed isothermal. The residual criteria are set in order that the maximal value of all the residuals should be lower than 10-4. A second order upwind scheme for the convection terms is used and all the diffusion terms are discretized with the second order central difference scheme.  Mesh generation and Grid resolution study Considering the 2D airfoil geometry, a C-Grid block method is used to create a structured mesh for the fluid domain with the software Icemcfd. The C-Grid is created by using a 10/20 chord size, which means that upstream of the leading-edge, the fluid size is around 10 times the chord length and downstream of the trailing-edge, the fluid domain has the size of 20 times the chord length. An O-Grid is also created around the profile to refine the boundary layer. A value of 𝑦𝑦+= 1 is used which corresponds to a distance of the first nodes of 0.05 mm from the profile. To avoid mesh issues dependency, a grid sensitivity study was performed. Four meshes are tested to better understand the influence of the mesh on the results and to find a suitable grid for the simulations. The non-dimensional wall distance 𝑦𝑦+ remained constant for all grids but the number of layers inside the boundary layer and in the circumferencial direction of the body are refined with a ratio of 1.5 in both directions. The different characteristics of the grids are described in Table 2, and Figs. 6 and 7 represent the coarse and the extra-fine meshes. 362 J. Piquee and C. Breitsamter. 7 Figure 8 – Elasto-flexible 2D model mounted in wind tunnel test section Table 3 – Characteristics of the four different grids Parameters Coarse Medium Fine Extrafine Total Nodes Number 74750 141210 461030 690730 Normal layer 80 120 180 270 Circumferential layer 90 140 210 345 Minimal Angle 25.5 28.8 28.8 28.3 Aspect Ratio 1266 1865 1005 633 2.3 Experimental analysis The experiments focusing on force measurements have been conducted in a Göttingen type wind tunnel with an elasto-flexible 2D wing model. The 2D model consists of an extension of the geometry described in Fig. 1 in the third direction set up on a support as it can be seen in Fig. 8. The model has an overall span of 564 mm and a nominal chord of 220 mm, giving an aspect ratio of AR = ~ 2.5. In order to limit 3D flow effects, end-plates of circular planform are used to reproduce twodimensional flow condition. These endplates have a diameter of three times the airfoil chord and are made of Plexiglas to provide an optical access to the elastoflexible airfoil. The force measurements were conducted for an initial elongation of ∆Lo = 2% at a Reynolds number of Re = 278750. Figure 6 – Coarse Mesh Figure 7 – Extra-fine Mesh 363 J. Piquee and C. Breitsamter. 8 3 RESULTS & DISCUSSION 3.1 Grid dependency The results of the grid sensitivity study for α=0°, 4°, 8°, 10° and 14° are given in Fig. 9, where the lift coefficient is plotted for the four different meshes. The numbers 1, 2, 3 and 4 of the abscissae axes correspond to the different meshes from the coarse to the extra-fine. For both models, the coarse and medium meshes overestimate or underestimate the lift coefficient for all angles of attack by sometimes 1.16% compared to the results obtained with fine and extra-fine meshes. For the SST turbulence model, the lift coefficient converges for α=2°, 4°, 10° and 14° with the fine mesh: the fine mesh appears fine enough to achieve a grid independency in the results. For the γ-Reθ transition model, the lift coefficient converges for the angles of attack α=4°, 8° and 14°. For α=0° and 10° there is a difference in the results between the fine and extra-fine meshes, which does not exceed 0.57%. Therefore, it was decided to perform the FSI simulations with the extra-fine mesh, which appears fine enough that there is not noticeable influence in the results. 3.2 γ-Reθ Transition model The results with the γ-Reθ transition model are presented in Fig. 10, where the lift and drag coefficients are plotted as function of the angle of attack. In Fig. 10, the results for the Panel method, the experimental data9 and the rigid counterpart of the elasto-flexible concept are also plotted. 1 - Coarse Mesh 2 – Medium Mesh 3 – Fine Mesh 4 – Extra-Fine Mesh Figure 9 – Grid sensitivity study for the SST model the γ-Re θ model for the angles of attack of α=0°, 4°, 8°, 10° and 14° 1 2 3 4 0.17 0.175 0.18  = 0deg Cl 1 2 3 4 0.55 0.6 0.65  = 4deg Cl SST Turbulent  -  1 2 3 4 0.95 1 1.05  = 8deg Cl 1 2 3 4 1.1 1.15 1.2  = 10deg Cl 1 2 3 4 1 1.05 1.1  = 14deg Cl 364 J. Piquee and C. Breitsamter. 9 Figure 10 – Lift and Drag coefficients as function of the angle of attack for the γ-Reθ model The results show that the lift coefficient of an elasto-flexible wing is much higher than the lift obtained for its rigid counterpart whereas the drag coefficients are nearly equivalent. As the elasto-flexible wing can adapt itself to the flow environment, the pressure distribution along the wing causes a positive deformation of the membrane for positive angles of attack which results in an augmentation of the curvature of the profile and therefore an increase in the lift coefficient. At α=0° the lift coefficient of the elasto-flexible wing is about three times higher than the lift coefficient of the rigid wing (0.544/0.175) due to the curvature of the elasto-flexible wing. The lift coefficient then linearly increases until α=7° where the lift reaches a maximum of 1.33; from α=7° to α=16° the lift coefficient of the elasto-flexible wing stays constant which is a characteristic feature of the adaptivity of such a concept. The stall phenomenon appears over a larger angle of attack range compared to a rigid wing resulting in a smoother and delayed stall region. The elasto-flexible wing lift coefficient drops at α=16° whereas the lift for its rigid counterpart drops at α=12°. In this case, Figs. 11 and 12 illustrate the turbulence intensity of the flow at α=7° and α=9°. At α=7°, a laminar separation bubble can be seen on the upper side of the profile at around 2/3 of the chord, whereas at α=9°, the laminar separation bubble migrates to the leading-edge of the profile. Therefore, the detachment of the turbulent boundary layer appears sooner at α=9° compared to α=7°. The pressure distribution at α=9° for the region of the turbulent boundary layer separation point up to the end of the profile is then constant and exhibits a reduced suction level compared to the pressure distribution at α=7°; this explains the decrease of the deformation of the upper-sided membrane and therefore the small decrease of the lift coefficient. From α=9° to α=16°, the point where the boundary layer detaches, moves to the leading-edge resulting in a slight and slow decrease of the deformation of the membrane (due to the constant pressure): therefore the stall phenomenon appears smooth until the boundary layer detaches at the leading-edge causing the loss of lift beyond α=16°. As the Panel theory does not consider the viscosity phenomena, the pressure distribution 0 5 10 15 20 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 alpha ° Cl Cl CFD Elastoflexible Cl Panel Elastoflexible Cl CFD Rigid Cl Experimental 0 5 10 15 20 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 alpha ° Cd α α