Acta Polytechnica doi:10.14311/AP.2017.57.0272 Acta Polytechnica 57(4):272–281, 2017 © Czech Technical University in Prague, 2017 available online at http://ojs.cvut.cz/ojs/index.php/ap AERODYNAMIC LOAD OF AN AIRCRAFT WITH A HIGHLY ELASTIC WING Pavel Schoř Institute of aerospace Engineering, Faculty of Mechanical Engineering ,Brno University of Technology , Technicka 2896/2, 616 69 Brno Czech Republic correspondence: schor@fme.vutbr.cz Abstract. In this article, a method for calculation of air loads of an aircraft with an elastic wing is presented. The method can predict a redistribution of air loads when the elastic wing deforms. Unlike the traditional Euler or Navier-Stokes CFD to FEM coupling, the method uses 3D panel method as a source of aerodynamic data. This makes the calculation feasible on a typical recent workstation. Due to a short computational time and low hardware demands this method is suitable for both the preliminary design stage and the load evaluation stage. A case study is presented. The study compares a glider wing performing a pull maneuver at both rigid and and elastic state. The study indicates a significant redistribution of air load at the elastic case. Keywords: aerodynamics, panel method, finite element method, wing load, fluid structure interaction. 1. Introduction 1.1. Problem overview A flight of a modern glider or a High Altitude Long Endurance (HALE) aircraft is a representative exam- ple of the fluid-structure interaction problem. These flying vehicles are a result of a search for the most aerodynamically effective shape. This leads to a need for a wing with a large wingspan and an aspect ratio. Most recent gliders have a wingspan ranging from 18 to 30 meters with an aspect ratio approximately from 30 to 50. A NASA Helios HALE aircraft had a wing with a wingspan of 75m and aspect ratio of 30 [1]. During a flight of such aircraft, the wing undergoes large deformations, still in the elastic regime, as shown on Figure 1. The Helios HP-03 experimental HALE aircraft crashed on June 26, 2003 after encountering atmo- spheric turbulence. The aircraft developed a wing Figure 1. The Helios aircraft undergoing large defor- mation [1]. tip displacement of more than 12m and consecutive pitch oscilations, resulting in structural failure of the leading edge structure, but leaving the main spar un- damaged. After the crash of Helios HP03 prototype, see Figure 2, a demand rose for new methods, which are capable to compute loads of highly elastic, “mor- phing” aircrafts [1]. In this article, a method for a solution of quasi-static load cases of such aircraft is presented. Moreover, recent airworthiness regulations explicitly require to calculate the load redistribution in case “if deflections under load would significantly change the distribution of external or internal loads” [2, 3]. 1.2. Methods overview Various methods are traditionally used in order to capture the change of the aerodynamic load generated in rigid and elastic state. The simplest option is to combine a vortex lattice method [4] with a simple beam finite element model [5] . The most sophisti- cated approach consists of Navier-Stokes CFD solver Figure 2. The Helios aircraft after structural failure during flight [1] 272 http://dx.doi.org/10.14311/AP.2017.57.0272 http://ojs.cvut.cz/ojs/index.php/ap vol. 57 no. 4/2017 Aerodynamic Load Of An Aircraft With A Highly Elastic Wing for aerodynamic input combined with a detailed FEM structural model with solid elements [6]. The aerody- namic load is then directly interpolated to the FEM model. Both of these approaches need to iterate be- tween the aerodynamic and structural model to obtain the equilibrium. The simple model with the lifting line theory provides the fastest solution times, how- ever, errors may be introduced because of an improper pitching moment distribution due to the simplified VLM aerodynamic model. The detailed model with the Euler or Navier-Stokes CFD solver with a detailed FEM structural model can provide the most accurate prediction of the load redistribution. Typically, the CFD volume mesh has to be updated according to the results from the FEM solution. This may cause issues with the deformation of the CFD volume mesh when the deformation of the wing is large. Recently an Immersed Boundary Method is used to solve this kind of fluid-structure interaction [7] When assuming the costs of the computational time of the CFD solver and also assuming that up to 10 iterations are usually required to obtain the converged solution for one particular angle of attack, the use of the CFD solver is restricted to resend high per- formance computer clusters, therefore making it very expensive and not feasible for a typical glider manufac- turer neither at a design stage nor at a load evaluation stage. In this article, a 3D panel method is used, combined with a geometrically nonlinear finite element method with beam elements to model the structure of the wing. The savings of computational time when compared to the CFD N-S solver and the differences in computed air loads justify the use of this low fidelity aerodynamic tool, making it possible to compute large number of load cases on an average workstation. 1.3. Eligibility Based on assumptions made in Sections 2.1 and 2.2, the method presented in this article is suitable for analysis of quasi-static maneuvers of an aircraft with elastic wings, operating at a low Mach number (M < 0.3) and sufficiently high Reynolds number (Re > 1000000). This includes the HALE aircraft during its low altitude flight or a glider performing a pull maneuver. 2. Methods 2.1. Aerodynamic model A fluid motion is traditionally described by set of equa- tions called Navier-Stokes equations [8]. The number of the exact solution of Navier-Stokes equations is small. For practical problems, these equations are solved numerically [9]. Solving the Navier-Stokes equations provides a de- tailed description of the flow, however the computa- tional costs are enormous. Therefore some assump- tions are made to reduce the complexity of the prob- lem. If the flow is assumed to be steady, inviscid, Figure 3. Potential flow domain. incompressible and irrotational, is called a potentital flow and it can be described by Laplace’s equation [8] : ∇2Φ = 0 (1) Where Φ is a scalar function called velocity potential and velocity −→q at each point can be obtained by its derivation: −→q = ∇Φ (2) An approach from [10] is used to solve the Laplace’s equation. Assume a cross section of a wing as shown in Figure 3. The surface S∞ encloses the problem at infinity, surface S represents a body (wing) and surface W represents its wake. The surface S+W divides the domain into two regions: the external region with the flow field of interest and velocity potential Φ and the internal region with the fictitious flow and velocity potential Φi. The surfaces are modeled by a doublet and source singularities. Green’s theorem is applied to both the outer and inner region and velocity potential ΦP is obtained by combining both expressions. [10] The velocity poten- tial ΦP gives a velocity potential anywhere in the two regions. It is expressed in terms of surface integrals in terms of velocity potential and its normal derivative over the boundary surface as ΦP = 1 4π ∫∫ S+W +S∞ (Φ− Φi)−→n · ∇ 1 r dS − 1 4π ∫∫ S+W +S∞ 1 r −→n · (∇Φ−∇Φi) dS. (3) Where r is the distance from the point P to the element dS on the surface, and −→n is a unit normal vector of the element dS. [10] The first integral in (3) represents the disturbance potential from a surface 273 Pavel Schoř Acta Polytechnica distribution of doublets with a density of (Φ− Φ∞). The second integral represents the contribution of sources with a density of −−→n · (∇Φ−∇Φi).[10] Next, the internal Dirichlet boundary condition is introduced. This boundary conditions sets the Internal flow equal to the Onset flow by Φi = Φ∞, (4) ΦP = Φ∞. (5) Following the procedure from [8, 10, 11], the surface is discretized into n surface panels and nw panels in the wake region with a constant distribution of singu- larities. The (3) is also discretized and the Dirichlet boundary condition is evaluated at centroid of each surface panel. These points are called Collocation points. The result of this procedure is a linear equa- tion [8] n∑ k=1 Ckµk + nw∑ l=1 Clµl + n∑ k=1 Bkσk = 0, (6) where Ck, Cl are doublet influence coefficients of body and wake panels, Bk is the source influence coefficient, µk is the doublet strength of a body panel, µl is the doublet strength of a wake panel, σk is the source strength of a body panel. The σk is set as σk = −→nk · −−→ Q∞, (7) where −→nk is the panel normal unit vector and −−→Q∞ is the freestream velocity. The method described above is known as the “first order panel method” [8]. It is used as the only source of aerodynamic data for the fluid-structure interac- tion. The formulation of the panel method is iden- tical to well-known and validated codes such as the VSAERO [10] and PMARC [11]. Boundary layer. The potential flow may give use- ful results for some applications, however for applica- tions where viscous effect become more important, the results obtained under assumption of potential flow may become inaccurate in terms of over predicted lift curve slope. A typical example of a flow with a signif- icant viscous effect is a glider wing at a low Reynolds number or a flow over a multi-element wing. The viscous effects can be introduced into the poten- tial flow model by displacing the actual body surface by a displacement thickness of the boundary layer [12]. The actual displacing of the body surface is per- formed by a modification of panels source strength in (7) to σk = −→nk · −−→ Q∞ + d ds (ueδ ∗), (8) where s is a distance on the body, ue is a panel edge velocity and δ∗ is a displacement thickness computed by the boundary layer analysis [12]. 2.2. Structural model The structure of the wing is modeled using the FEM with beam elements only. Each node of the element has six degrees of freedom. The stiffness properties are assumed constant for each element. The stiffness matrix of such element can be found in many FEM related textbooks [13]. However such elements are usually derived for isotropic materials. In fact, wings of modern gliders are made of a carbon fiber rein- forced composite orthotropic material. A wing made of such material exhibits coupling between bending and torsion. This causes issues, when one tries to reduce stiffness properties of a composite wing into single stiffness matrix of a beam element. A method from [14] was adopted to solve this issue. As the wing undergoes large deformations it is nec- essary to capture the geometrical non linearity. To solve this issue a Newton/Rhapson iterative method is used. The geometrical non-linearity, however, causes difficulties when calculating elements stiffness matri- ces — one must determine the derivation of the global stiffness matrix (tangent stiffness matrix). A method from [15] was adopted, which provides a numerically generated tangent stiffness matrix for truss and beam elements. 2.3. Aero-structural interaction Work flow. A brief overview of the interaction be- tween the aerodynamic and internal forces in the wing structure is shown on Figure 4. The interaction is done by using a modified Newton-Rhapson method. This method and steps called the “Compute aerody- namic forces” and “Compute wing displacement” from Figure 4 are explained in detail in the paragraph Load control scheme. Computing the rigid lift curve as a first step is not necessary, the main reason for this is to find the zero lift angle of the attack α0, where the computation of elastic lift curve starts. Starting from α0 makes the wing almost unloaded and the Newton-Rhapson method is likely to converge. Also, the increment of the angle of attack ∆α should be small enough to assure the convergence, it is usually from 0.5° to 2°. Interacting element. The interacting element is explained in Figure 5. It consists of aerodynamic panels, but there is only one element in the span-wise direction and one FEM beam element. Forces and moments are always integrated from the wingtip to plane of symmetry, at a point which is also a node of the finite element model. The load applied at the FEM node is a sum of forces moments from all panels belonging to the interacting element as shown on Figure 5: −→ f = n∑ i=1 −→ fi , (9) −→m = n∑ i=1 −→ fi × ri, (10) 274 vol. 57 no. 4/2017 Aerodynamic Load Of An Aircraft With A Highly Elastic Wing Figure 4. Aero-structural interaction — work flow. where fi is the force on ith panel in the global coor- dinate system x, y, z, as shown in Figure 5; ri is the vector between the reference node and control point of ith panel. The result of the FEM analysis is a displacement vector u at each node. Coordinates of aerodynamic panels are updated at each load step using nodal translational and rotational displacements, according to (11). Also, the nodal coordinates are updated as explained in [15]. PTSnew = ( Rx(∆φx)Ry(∆φy)Rz(∆φz) ) PTSold + ∆U, (11) where ∆φx, ∆φy, ∆φz, are rotational increments in the global coordinate system calculated from the last step of the Newton-Rhapson iteration, see bellow; ∆U is the displacement increment; Rx(φx), Ry(φy), Rz(φz) are rotational matrices defined by Rx(∆φx) = 1 0 0 0 cosφx − sinφx 0 sinφx cosφx  , (12) Ry(∆φy) =  cosφy 0 sinφy 0 1 0 − sinφy 0 cosφy  , (13) Rz(∆φz) = cosφz − sinφz 0 sinφz cosφz 0 0 0 1  . (14) Figure 5. Panel method — FEM interaction. Load control scheme. In the Newton/Rhapson method, the total load is divided into m load steps and for each load step, a Newton iteration is done until an equilibrium between the internal forces and external load is found. After that, the external load is increased by the load step. However, in the case of an elastic wing, the total external load is unknown because the load (lift) also depends on how much the wing is deformed. Therefore, the standard load control in the Newton/Rhapson method was changed in the following way: The load is controlled by in- creasing the angle of attack as explained in Figure 4. Usually the lift generated by the wing (sailplane) is the dependent variable in the u–F curve. To compute the load generated by the elastic wing, a following scheme is used: (1.) The first iteration starts at an angle of attack of zero lift of a rigid wing. The equilibrium between the external load computed from the panel method and internal forces is computed using the Newton iteration. (2.) The angle of attack is increased by a small step, e.g., 1° and the equilibrium between the aerody- namic load and the internal forces is computed by the Newton iteration. (3.) The angle of attack is further increased until the total load (lift) is equal to the desired load (lift). The result is a nonlinear lift curve. For each load step, the vector of internal forces is also a nonlinear function of previous load steps. Therefore when com- puting a load of an elastic wing it is convenient to compute a wing lift curve as described above with very small increment of the angle of attack, e.g., 0.5° and to save the vector of deformation and internal forces for each converged load step. When comput- ing large number of load cases, significant savings in the computational time are achieved when the pre-computed deformation vector and internal forces vector are loaded for the closest lower angle of at- tack. One then needs to compute only the difference between the pre-computed angle of attack and the desired angle of attack. 275 Pavel Schoř Acta Polytechnica Figure 6. Aero-structural interaction — load control in Newton-Rhapson method. Figure 7. Coordinate systems. 3. Numerical study 3.1. Overview A numerical study was conducted to demonstrate the importance of the elasticity of the wing for its load. Assuming a flight envelope from [2, par. CS22.333], and considering the point A of the envelope to be crit- ical, the aerodynamic load at this point is evaluated using three different methods: (1.) panel method with elastic wing; (2.) panel method with rigid wing; (3.) lifting line theory. The lifting line theory is used as a validation tool for the panel method. 3.2. Coordinate systems A global right-handed coordinate system OXY Z de- fined in Figure 7 is used. Additional local coordinate systems oixiyizi are used for evaluation of the local load. As the wing deforms, these coordinate systems Property Symbol Value Area A 11.55m2 Span L 21.0m Root chord cr 0.7m Tip chord ct 0.4m Dihedral Γ 0 rad Sweep Θ 0 rad Table 1. Wing geometry. Property Symbol Unit Cross section area AS m2 Second moment of area Ix m4 Second moment of area Iz m4 Torsional constant Jy m4 Offset of neutral axis ∆Xn m Offset of neutral axis ∆Zn m Offset of elastic axis ∆Xt m Offset of elastic axis ∆Zt m Young modulus E N/m2 Shear modulus G N/m2 Table 2. Wing cross sectional properties. follow the deformed geometry. The free stream veloc- ity V with positive angle of attack is also shown in Figure 7. 3.3. Wing specimen A simple tapered wing is examined in this study. The wing is designed as a planar wing in undeformed state. The wing has one straight main spar located at 0.25c, where c is local chord. The main reason for this geometry is the evaluation of the load by the lifting line theory, which can not be used for a non-planar wing geometry and nor swept wings. 276 vol. 57 no. 4/2017 Aerodynamic Load Of An Aircraft With A Highly Elastic Wing Figure 8. Spanwise variation of Iz, Ix, Jy. Figure 9. Airfoil E603 inviscid aerodynamic data. Geometry. Geometry of the wing is described in Table 1. An Eppler E603 airfoil is used for all sections of the wing. Stiffness properties. When analyzing the defor- mation of the wing by the FEM analysis described above, an instrumental step is the determining of all cross-sectional characteristics defined in Table 2 for all cross sections of the wing. In this example, the Young’s modulus E and shear modulus G are assumed as constant for the whole wing with following values: E = 60GPa, G = 7GPa. The wing has only one spar at 0.25c and the nodes are coincident with the spar. Therefore the offset of the elastic axis is almost negligible and both values are assumed to be zero. The neutral axis is located at the local centroid with non-zero offsets, but these offsets are not shown in this paper. The spanwise variation of the Second moments of area and torsional constant is shown on Figure 8. Finite element model. The wing is divided into 162 spanwise sections. Each section has one two- noded beam element described in Section 2. Nodal coordinates always correspond to a location of a main spar. The node located in the plane of symmetry has all degrees of freedom prescribed to zero. All coupling terms in the element stiffness matrices related to the bending-torsion coupling are intentionally set to zero for clarity. Airfoil data. The aerodynamic characteristics of the Eppler E603 airfoil were obtained from the XFOIL inviscid analysis, as shown in Figure 9. This aerody- namic data are used only for the lifting line analysis. 3.4. Loadcase VA Assuming a maximum lift coefficient of the wing CLmax = 1.4, from Figure 10, the flight conditions at point A of the maneuvering envelope are shown in Ta- ble 3 The total mass of the glider is MTOW = 750 kg, but inertia forces acting on the wing are not assumed in the analysis for clarity. 277 Pavel Schoř Acta Polytechnica Figure 10. Comparison of rigid and elastic wing — lift. Property Symbol Value Airspeed V 64.0m/s Wing lift coefficient CL 1.4 Load factor n 5.3 Table 3. Load case VA. Deformation sequence. Deformation sequence for various wing lift coefficient is shown in Figure 11. It must be noted that this deformation sequence is valid only for airspeed V = 64.0m/s. As the deformation of the wing depends on forces and moments applied on the wing surface, for a different airspeed, the cor- responding deformed shape will be different. Integral aerodynamic characteristics. Integral aerodynamic characteristics of the rigid and elastic wing are compared on Figures 10 and 12. A close agreement for rigid cases can be observed. For the elastic case a significant decrease of the slope of the wing lift curve can be seen in Figure 10. Additionally, the increase of pitching moment of the elastic wing can be seen in Figure 12. The differ- ences between rigid and elastic cases are discussed in Section 4 Load comparison. Actual comparison of the load for the rigid and elastic wing performing a steady ma- neuver with the same load factor n = 5.3 is shown in Figures 13–15. Only bending and torsional moments are shown. The differences between the rigid and elastic case are discussed in Section 4 4. Discussion 4.1. Validation of results There are two main options to validate the result computed by the presented method: Figure 11. Deformation sequence at VA. Numerical validation. A numerical validation means to prerform the case study using a RANS CFD solver as a source of aerodynamic data coupled to a wing modeled by shell elements or by 3D solid ele- ments. Experimental validation. Experimental valida- tion means to actually fly a glider, perform the pull maneuver and measure the displacement of the wing. Conclusion on validation. Unfortunately both numerical and experimental validation would require much more effort and funds than what was spent for this whole work. A literature review also did not give any example, which could be used as a validation benchmark. Therefore, this work is presented only as a numerical study without any other validation. 4.2. Decrease of wing lift curve slope As the lift generated by the elastic wing increases, the wing bents more, as shown in Figure 11. As the wing bents more, the slope of the lifting curve drops. This is demonstrated in Figure 10, which shows a reduced slope of the lift curve of the elastic wing, compared to the rigid wing. This phenomena could possibly be explained by the reduction of the projected area of 278 vol. 57 no. 4/2017 Aerodynamic Load Of An Aircraft With A Highly Elastic Wing Figure 12. Comparison of rigid and elastic wing — pitching moment. Figure 13. Comparison of rigid and elastic wing — Normal bending moment Mx. the wing and the change of direction of the lift force. This has following consequences: (1.) When the prescribed total lift is reached, the wing operates at a higher angle of attack. In the particular case presented in Section 3 the difference of the critical angle of attack is 6°. (2.) The normal bending moment is slightly reduced, see Figure 13. The tangential bending moment is, however, increased due to the increased angle of at- tack, see Figure 15. In the particular case presented in Section 3 the tangential bending moment of the elastic wing is 1.34 times higher, than in the rigid case. The normal bending moment is, however, 0.87 times lower. 4.3. Increase of wing pitching moment The explanation of the increase of the wing pitching moment requires more labor than the explanation of the decrease of the wing lift curve slope. First, assume a wing cross section located at the wing tip. The wing is of an undeformed shape and operates near the critical angle of attack. The section generates a pitching moment around the main spar MY 1 = 1 2ρv 2S1c1cm1(α). The resultant of the aerodynamic force is assumed to act at the main spar (0.25c). The direction of the force is upward and forward (negative X, positive Z). When the wing has zero dihedral and zero sweep, the contribution of this force to the pitching moment is zero, since the vector of the section reference point related to the wing reference point ~r1 = [0,−10.37, 0]. The wing reference point is defined as a point (node) on the main spar in the XZ plane. Next, assume a highly deformed wing shown in Figure 16. The wings tip section is displaced by a vector ~d1 = [−0.12, 2.5, 5.15], the vector of the section reference point related to the wing reference point is now ~r1 = [−0.3,−8, 5.15]. The resultant of the aerodynamic force at the wing 279 Pavel Schoř Acta Polytechnica Figure 14. Comparison of rigid and elastic wing — Torsional moment My. Figure 15. Comparison of rigid and elastic wing — Tangential bending moment Mz. reference point is computed as ~f1 = [−6.07, 49.38, 9.31], then, the resulting moment at the wing reference point generated by the force ~f1 is ~m1 = ~r1 × ~f1 = [−329.05,−28.51,−63.32]. When the pitching moment of the section is evalu- ated, the contribution of ∆MY 1 = −28.51Nm causes to increase the pitching moment coefficient of the sec- tion by a value of −0.55 mostly due to the vertical displacement of the section and the force resultant directing forward. 4.4. Wing maximum lift In this article, the maximum lift coefficient of the rigid wing was computed using the nonlinear lifting line theory as CLmax = 1.4. The airfoil sectional data used in the analysis were computed by the XFOIL inviscid analysis, assuming to be close to the experimental data. Since the purpose of this article is to present a Figure 16. Pitching moment of highly deformed wing. potential problems of an aircraft with a highly elastic wing, rather than to make a prescription based on practical results, the maximum lift coefficient of the elastic wing is also assumed to be CLmax = 1.4. 280 vol. 57 no. 4/2017 Aerodynamic Load Of An Aircraft With A Highly Elastic Wing For a practical application of the method, it may be in some cases necessary to limit the maximum lift of the elastic wing directly in the panel method. This is achieved by employing the viscous boundary layer coupling in the panel method, as in [10, 12]. Also, when including the viscous drag, the tangential force acting at a high angle of attack at wing sections may be reduced, resulting in a lower pitching moment. 5. Conclusions and future work 5.1. Future work The 3D panel method usually provides a fast and reliable source of aerodynamic data. However, when the viscous effects are significant as in the case of the HALE aircraft or a glider with bumps on its lift curve, the computed load may not be adequate. To solve this issues, it is necessary to extend the panel method with the calculation of the viscous boundary layer [10]. With this extension, the method may be applicable to solve practical problems. 5.2. Conclusions A simple and robust method for a calculation of a load of an aircraft with highly elastic wings was de- veloped. Two eventual problems of an aircraft with highly elastic wings were discovered: (1.) increase of tangential bending moment; (2.) increase of pitching moment. The main advantages of the presented method are simple definition of the wing stiffness, absence of the volume mesh around the wing and problems related with remeshing at large deformations. The main ad- vantage is a very fast computational time, one aeroe- lastic iteration can be computed within minutes on an average workstation. When the problem is not re- duced and the aerodynamic inputs are solved with the RANS CFD solver, it takes significantly more time on much more powerful computers, making such method unsuitable for a calculation of large numbers of load cases. Following the NASA Helios mishap report [1], a relation between a highly bent wing and its pitching moment was found by the present method. However, due to the fact that the Helios HP03 operated at a very low Reynolds number, making the present method unsuitable and due to the lack of any data regarding the Helios aircraft, no conclusions could be made. The demand for new analysis methods from the report [1] can be possibly fulfilled. Due to a very fast computational time, it is possible to extend the method into the time domain, and therefore also analyze the unsteady cases [8]. References [1] T. Noll, J. Brown, M. Perez-davis, et al. Investigation of the Helios Prototype Aircraft Mishap - Volume I Mishap Report. National Aeronautics and Space Administration, 2012. [2] EASA. Certification Specifications for Sailplanes and Powered Sailplanes CS-22. European Aviation Safety Agency, 2008. [3] EASA. Certification Specifications for Normal, Utility, Aerobatic, and Commuter Category Aeroplanes. European Aviation Safety Agency, 2012. [4] N.-L. R. Center. Vortex-lattice utilization. NASA SP-405. NASA, 1976. [5] M. Drela. Integrated simulation model for preliminary aerodynamic, structural, and control-law design of aircraft. AIAA Paper 99:1394, 1999. [6] J. Haider, C. H. Lee, A. Gil, et al. A first order hyperbolic framework for large strain computational solid dynamics in openfoam. In World Congress on Computational Mechanics and The Asia-Pacific Congress on Computational Mechanics, pp. 688–688. 2016. [7] C. S. Peskin. The immersed boundary method. Acta Numerica 11:479âĂŞ517, 2002. doi:10.1017/S0962492902000077. [8] J. Katz, A. Plotkin. Low-speed aerodynamics, from wing theory to panel methodsmcgraw-hill. New York 1991. [9] A. J. Chorin. Numerical solution of the navier-stokes equations. Mathematics of computation 22(104):745–762, 1968. [10] B. Maskew. Program VSAERO Theory Document: A Computer Program for Calculating Nonlinear Aerodynamic Characteristics of Arbitrary Configurations. NASA contractor report. National Aeronautics and Space Administration, 1987. [11] D. Ashby. Potential Flow Theory and Operation Guide for the Panel Code PMARC. National Aeronautics and Space Administration, 1999. [12] T. Cebeci, J. Shao, F. Kafyeke, E. Laurendeau. Computational Fluid Dynamics for Engineers: From Panel to Navier-Stokes Methods with Computer Programs. Springer, 2005. doi:10.1007/3-540-27717-x. [13] A. Píštěk, O. Grégr, V. Kahánek, R. Böhm. Pevnost a životnost letadel i. Vysoké učení technické v Brně, Brno 1988. [14] R. d. F. Lopez. A 3D Finite Beam Element for the Modelling of Composite Wind Turbine Wings. Master’s thesis, Royal Institute of Technology,Stockholm, Sweden, 2013. [15] S. Lebofsky. Numerically Generated Tangent Stiffness Matrices for Geometrically Non-Linear Structures. Ph.D. thesis, University of Washington, 2013. 281 http://dx.doi.org/10.1017/S0962492902000077 http://dx.doi.org/10.1007/3-540-27717-x Acta Polytechnica 57(4):272–281, 2017 1 Introduction 1.1 Problem overview 1.2 Methods overview 1.3 Eligibility 2 Methods 2.1 Aerodynamic model 2.2 Structural model 2.3 Aero-structural interaction 3 Numerical study 3.1 Overview 3.2 Coordinate systems 3.3 Wing specimen 3.4 Loadcase VA 4 Discussion 4.1 Validation of results 4.2 Decrease of wing lift curve slope 4.3 Increase of wing pitching moment 4.4 Wing maximum lift 5 Conclusions and future work 5.1 Future work 5.2 Conclusions References