Abstract
This paper presents closed form solutions to the classical beam elasticity differential equation in order to effectively model the displacement of standard aerodynamic geometries used throughout a number of industries. The models assume that the components are constructed from in-plane generally anisotropic (though shown to be quasi-isotropic) composite materials. Exact solutions for the displacement and strains for elliptical and FX66-S-196 and NACA 63-621 aerofoil approximations thin wall composite material shell structures, with and without a stiffening rib (shear-web), are presented for the first time. Each of the models developed is rigorously validated via numerical (Runge-Kutta) solutions of an identical differential equation used to derive the analytical models presented. The resulting calculated displacement and material strain fields are shown to be in excellent agreement with simulations using the ANSYS and CATIA commercial finite element (FE) codes as well as experimental data evident in the literature. One major implication of the theoretical treatment is that these solutions can now be used in design codes to limit the required displacement and strains in similar components used in the aerospace and most notably renewable energy sectors.
1. Introduction
Currently there is great deal of interest in the design of structurally efficient aerodynamic structures constructed of advanced composite materials. In the aerospace sector these are defined as a number of chemically distinct materials being macroscopically bonded with the resultant amalgamation being superior than original chemically distinct constituents [1]. The emphasis in the renewable and aerospace communities is to embed a strong (though usually brittle) ceramic material such as glass or carbon fibre in a polymeric material such as epoxy or more recently the biodegradable polypropylene [2]. The resulting composite material is a heterogeneous anisotropic material with a reality high strength-to-weight ratio. The consequence is that the polymeric material is made a couple of orders of magnitude stronger in addition to significantly increasing the displacement at failure when compared to the original ceramics employed as the reinforcement.
Consequently, there has been much development of manufacturing processes in order to produce components such as wind turbine blades from these rapidly becoming inexpensive materials. Since the macroscopic bonding results in a heterogeneous more flexible material then the conventional assumptions of material failure tend to be inappropriate. It is known that most composite material wind-turbine blades employed in the industry fail via a combination of strain and displacement modes, where fibre-glass is concerned; it is now widely accepted [1] that such materials fail via maximum stain. In order to reduce both the cost of physical testing and time prohibitive computationally expensive simulations in the design phase, the need for accurate and effective modelling of the displacement and strain fields in materials which are used to construct wind turbine blades, under standard testing conditions, has very recently become of utmost importance.
It is generally accepted that the geometry of a blade is determined by aerodynamic properties, which are defined by aerofoil characteristics such as chord, twist, and thickness distribution. The structural design should incorporate material selection and spar/shear design with a purpose of minimizing the tip deflexion, thus limiting the chance of blade/tower collision [3]. Also of importance from a design viewpoint is to gain an understanding of the loads that a wind turbine blade undergoes during a design life of approximately 20 years [4].
Moreover, it is known that most structural failures of wind turbine blades are found in the blade root section, which has been confirmed via the application of a three-dimensional analytical model developed by El Chazly [5]. The results showed that maximum strains occurred at the root of the blade in the flapwise direction and that tapering blades diminishe strain values. In addition to saving material weight and cost, a successful blade design must satisfy a wide range of objectives:(1)maximize annual energy yields;(2)stall regulated machines;(3)resist extreme and fatigue loads;(4)restrict tip deflexion to avoid tower collisions;(5)avoid resonances.Objectives 1 and 2 are taken into consideration when a blade is designed to suit an aerodynamic design, while objectives 3 to 5 are taken into consideration with regard to structural design. When designing from an aerodynamic perspective the blade geometry is of vital importance. Though these two design processes are separated, it is often found that an iteration process will be needed between the two. This is to accommodate adequate blade thickness to house a spar/shear-web required to significantly reduce the displacement field and hence strain in the shell and consequently the layers of the composite material.
In short the work presented throughout this paper lays significant groundwork to the development of generic design-codes (expert system) for the development of structural aerodynamic components such as wind-turbine blades. The design solutions presented are in principle easily extended to allow all necessary geometrical and mechanical variations to be considered.
1.1. Governing Equations
From a structural (static) viewpoint any composite material wind-turbine blade can be modelled as cantilever with the root located at the hub. This explains the aforementioned observations [5] of structural failure at the root evident in the literature. Thus, it is obvious that the basic displacement equations of elasticity are applicable, albeit in anisotropic form for the heterogeneous materials considered [6].
Consider the following: where is the second moment of area about the -axis, is the force moment acting on the blade about the -axis, and is the elastic modulus of the composite material parallel to the -axis. Moreover, since the standard testing protocols [7] only consider loads in a single direction at right angles to the fibres then this simplifies the familiar expression for beam displacements; namely: Furthermore, if the fibres are assumed to align on application of load to the direction of maximum principal strain, the elastic modulus can be determined from the rule of mixtures, where is the total volume of the turbine blade, is the volume of the the reinforcement fibres, is the volume of the aforementioned polymeric matrix material, is the elastic modulus of the reinforcement phase, and is the elastic modulus of the matrix material. Incidentally for most practical design purposes so (3) reduces to the following: where is the volume faction of the the fibres. Moreover classical lamination theory (CLT) dictates [1] that some of the stiffness of the fibres is utilized for reinforcement in the traverse direction also, meaning that this expression has a more general form: where is a fibre orientation factor which is usually 0.5-0.6 for quasi-isotropic continuous layered composite materials and is also applicable for short-fibre composites, where it can be as low as 0.25 [1].
Additionally, since wall sections are generally small in comparison to the other cross-section dimensions then every infinitesimal second moment of area is proportional to the thickness of the section and an infinitesimal arc length; therefore any second moment area of the shell, symmetric about the abscissa (as simulated in this paper), can be determined from the following: where is the thickness of the material and which suitably describes the cross-section that is in the simulations that follow this, takes the form of a circle, an ellipse, or aerofoil as a function of the distance from the root.
Substitution of (5) and (7) into (2) renders the general governing elasticity shell differential equation as follows This equation is also applicable for quasi-isotropic and isotropic materials subjected to pure single direction bending.
With the displacement of the beam found, the global strain field in the shell from the simpler algebraic relationship: we note that the modelling presented in this paper investigates only cross-section variations with the material assumed to be constant; in principle variations in the material thickness and hence the effective modulus of the components can also be incorporated. Moreover, realistic empirical and/or computational fluid dynamics (CFD) load data can also be simulated. Thus, the models presented and their associated results are considered a representative benchmarking studies in order to demonstrate general purpose methods for preliminary design purposes.
2. Simulation Methods
Industrial standard design codes [8] dictate that for aerofoil sections, the thickness to chord ratio should be approximately 40% at the root and inboard region of the blade as this allows bending moments that are generated to be absorbed. This can be lowered to a value of 16% towards the tip of the blade. However, the blade thickness is in constant conflict between aerodynamic efficiency and stiffness and strength requirements. High performance aerofoils are desired for very specific aerodynamic reasons [9], that is, the thinnest blade possible, whereas structural requirements look for a thick enough cross-section to counter the load bearing elements. The maximum thickness is provided by the height of the spar-box or shear-web and has a direct influence on the section modulus on spar-box cross sections [10].
2.1. Benchmarking
In this section we present benchmarking studies required to validate the analytical and numerical modelling methods used prior to a more complete formulation of the tapered aerofoil section wind turbine blade. We will show that analytical solutions can be found for thin-wall tapered beams with elliptical sections. These solutions use extensions of the general double integration (antiderivative) method for the determination of beam deflexion. For example consider, the displacement field of circular shell section beam. This problem is being an extension of a typical advanced undergraduate mechanics of materials problem [11]. The key geometrical property, the second moment of area can be obtained via application of (7): where is the radius at the root of the beam, is the radius at the tip of the beam, and is a characteristic length with the length of the blade. Given that the aforementioned testing protocol [7] applies a point load at the tip of the blade this implies that the governing equation (8) becomes where is the force applied to the tip during testing. Initially this equation can be solved via direct double integration, and applying appropriate zero slope and deflexion boundary conditions at the root of the beam, to obtain the so-called analytical antiderivative solution: where is the effective elastic modulus of composite in the -direction and is the second moment of area at the root of the beam. This anti-derivative solution was initially validated using analog finite element and finite difference methods evident in the literature [12, 13]. We observed excellent agreement as indicated in Figures 1 and 2 with the finite difference (FD) and finite element (FE) numerical approaches. It is noted from the insert in Figure 1 that in general the shell approximation slightly overestimates the value of second moment of an area of the circular section investigated in this work [14]. This results in a minimal underestimation of the displacement field when compared to the solid anti-derivative formulation as depicted in the insert of Figure 2; also evident here is an underapproximation of the FE predictions from the commercial ANSYS software [15]. This initial benchmarking process is therefore offered here as a proof of methodologies employed throughout this paper.
We will show in the proceeding section that as the geometry of the beam cross-section becomes only slightly more complicated, the methods of integration required to obtain the necessary anti-derivatives become much more involved [12]. Moreover, for even the most idealized aerofoils (as those that follow) no explicit closed-form analytical solutions are possible, meaning that analogue numerical solutions must be obtained, which justifies the use of the numerical approaches outlined.
2.1.1. Elliptical Tapered Beam General Analytical Model
It will be shown in the proceeding section that from a structural design viewpoint most aerofoil sections have some elliptical geometrical nature (e.g., FX66-S-196 and NACA 63-621 [16, 17]). However, the inherent mathematical complexity with regard to the derivation of anti-derivatives of elliptical functions [18] has resulted in a lack of treatment of particularly elliptical shell beams and hence wind-turbine blades in the literature. For this reason an elliptical beam benchmarking regime was performed as a precursor to a complete and rigorous treatment of specific aerofoil sections: firstly as a verification of the second moment of area and then of the displacement and strain fields and finally and most importantly to provide the first general analytical treatment of the elliptical shell tapered beam problem absent in the literature.
The calculation of the second moment of area as a function of the the span was evaluated by two analytical and two independent numerical methods, namely, one-thousand point trapezoidal integration and adaptive Simpson quadrature. In the case of continuum solid formulation, the analytical model consisted of casting the exact formula as a function of the span variable; thus [11], where and are characteristic lengths based on the major and minor radii of the ellipse. Note here that this expression is now a function of only as the expression assumes a continuum solid cross-section. For reasons that became obvious later in this section (most notably the development of an explicit analytical model for elliptical shell sections) for elliptical shells (i.e., ), by expanding (13) and retaining only linear factors of we have Now we notice that the term in brackets is a constant as is the first term in the parenthesis and both are again effective length terms. This leads to the more lucid form of the second moment of area for a tapered elliptical shell: where
These formulae were validated via numerical integration of (7) using both the aforementioned integration methods and values extrapolated from suitable positions along a computer aided three-dimensional interactive application ( [19]) model. Figure 3 shows that for the elliptical shells investigated in this work (15) evaluates the required second moment of area well and hence was used in all subsequent formulations included as part of the blade model described in Section 2.2.
Hence substitution of (15) in the displacement differential field equation (8) gives where the partial fraction constants () can be evaluated from expansion and then comparison of the identity to render the fiollowing: We can find the slope equation by direct integration of (17) and imposing zero slope at the root This implies that the required displacement field can be found from the antiderivative of this expression and by imposing zero displacement at the root to render the following: with the final terms in the brackets being the constant of integration, meaning the analytical anti-derivative reduces to the following:
Therefore (16), (18), and (21) form the components of what will be referred to from here on in as the analytical shell model formulation; its solution evaluates the displacement of an elliptical shell tapering beam. As mentioned previously, to our knowledge this is the first time that such an explicit solution to this specific case has been presented in the literature.
The strain field can be found from the aforementioned algebraic definition (see (9)) and substitution of (13) and (15) to render continuum solid and shell models formulations; in the case of the latter this gives:
The validity of these models was verified using an explicit Shampine and Gordon [18] Runge-Kutta, resident in the Mathworks MATLAB [13] mathematical software. Further verification of these solutions was performed using two analogous finite element (FE) models, employing two separate but complementary commercial FE codes, namely, ANSYS and .
In the case of ANSYS the model was constructed from 1369 shell -type (approximate size range being ) elements of 5 mm thickness permitting both bending and membrane capabilities. Both in-plane and normal loads were permitted by the software, which implied that each of 1382 nodes produced had six degrees of freedom: these being translations in the principal nodal , , and directions and rotations about these axes. The mesh produced used advanced on-curvature sizing constraints with course relevance and span angle center with smooth relevance transition. In the solution set-up phase a fixed support (clamp) was applied to the the large end ellipse in order to simulate the root of a cantilever, while a force was applied to smaller end ellipse, thus simulating the blade tip load. A standard (default) linear static solution solution was then sought.
As a final comparative study a similar FE model was constructed using the advanced meshing tools resident in the code [19]. Here the Denizen frontal parabolic quadrangle shell -elements with an average size 66 mm and a thickness of 5 mm were employed. Identical mathematical boundary conditions were applied as in all the modelling methods described in this section; that is, a clamp was applied at the root of the blade and the generative structural analysis distributive force was applied at the tip.
2.2. Aerofoil Section Models
In this section we present analytical and numerical models for the deflexion of a tapering beam with specific aerofoil sections with and without the presence of a stiffening (shear-web) rib. We begin with predictions of the second moment of area for the aerofoil sections using the methods outlined in the previous section and then move on to the displacement and strain calculations which are compared with FE analogue models produced using the commercial ANSYS and codes.
2.2.1. Idealized Tapered Beam Model
It was necessary to fit the aerofoil sections to an elliptical and triangular section as demonstrated in Figure 4. Here, the end elevations of the idealized aerofoil sections without and with a reinforcement stiffening rib are shown in Figures 4(a) and 4(b), respectively, while an isometric MATLAB representation of the idealized tapered aerofoil shell beam used in the work is presented in Figure 4(c). The geometries modelled consisted of a root chord length of 600 mm and depth of 178 mm in line with specific aerofoils under investigation (FX66-S-196 and NACA 63–621 [16]). Here the geometry remains similar at the tip with a chord length of 110 mm being apparent and a resulting tip section depth of 33 mm.
The required second moment of area can be found exactly by noting that the section is half an ellipse plus the difference between two triangles; thus, application of (13) together with the standard result for that of the triangle renders the continuum solid section formula: Alternatively we could employ (15) and then the definition of the second moment of area to obtain a shell formula, this being These formulae were again validated using the previously described numerical integration method. What is noticeable here is that all the approximate methods tend to overestimate the second moment of area of the thin-wall aerofoil section with greater convergence being observed at the root. This said, all of the methods do predict the the second moment of area well.
2.2.2. Numerical Modelling
It is obvious from (23) and (24) (most notably the radical) that the extra terms imply that no closed form solution is possible for the displacement field; however, numerical finite-difference (Runge-Kutta [18]) solutions as discussed in the previous section are possible and hence were sought and verified using FE models created in the ANSYS and commercial codes. Examples of the FE model analogues to the analytical models are depicted in Figure 4 are shown in Figure 6 using (without the stiffening rib) and ANSYS (with a stiffening rib of 16 mm thickness); these models being typical. In the case of , the model shown consisted of 1471 quadrilateral 8-noded (Q8D) -type elements with parabolic interpolation functions. The average -value is being 50 mm, internal quadrilateral angle ranges ~, and an average aspect ratio is 1.168. The ANSYS model was analogous containing some 1175 -type isoparametric quadrilateral parabolic shell elements (average -value is 50 mm) with the stiffening rib as shown in Figure 6 or 1008 elements without the stiffening rib (not shown). Each ANSYS mesh utilized a course relevance center with a high degree of smoothing, fine span angel center, on curvature advanced size functionality with a default normal curvature angel of . Both models were used to simulate a standard design-code case as detailed in the literature [7, 8] (See Figure 5).
(a)
(b)
The finite-difference solutions were obtained via direct substitution of (23) and (24) into (8), while the solution of shell strain field was obtained from by substituting (23) and (24) into (9). Each of the solutions was obtained without and then with a stiffening rib. The stiffening (shear) web consisted of a further tapered solid rectangular beam located in the centre of the aerofoil section; see Figure 4(b) for details.
3. Results
The results of the FE displacement calculations from the ellipise, with an outside major axis of 300 mm, an outside minor axis of 89 mm, and shell of thickness 5 mm cantilever of length 4.2 m (These dimensions being analogous to the forthcoming analysis of the NACA 63-621 aerofoil section as detailed in [17]) benchmark case are shown in Figure 7 with an applied load of 8.5 kN. Here good agreement is demonstrated between the ANSYS and codes with Figure 7(a) showing a maximum deflexion of 359 mm calculated by the ANSYS code, while identical software calculations render a value of 357 mm as shown in Figure 7(b). The disparity of these values is being attributed, as one would expect, to the greater default mesh density imposed by the ANSYS software. The coarser code mesh render a stiffer model, indeed such overestimation of the stiffness is consequence of the FE numerical method. These maximum deflexion calculations are compared with the analytical model, of the previous section, and FD (Runge-Kutta) predictions in Table 1. In general excellent agreement between all the modelling methods is demonstrated. Best agreement is being observed between the Runge-Kutta and anti-derivative shell formulations which themselves are in close parity with the Runge-Kutta FD solution to the solid formulation with an average difference of just over 0.08% between the calculations. Less agreement is shown between FD and analytical (anti-derivative) solutions when compared with the FE calculations: discrepancies of 1.5% and 1.8%, respectively, for Finite element analysis (FEA) performed using the ANSYS and codes. This said, these differences are well within the default convergence settings of 5% (ANSYS) and 10% ().
(a)
(b)
Further demonstration of the suitability of the analytical elliptical shell beam model of the previous section is shown via the calculations of the displacement field in Figure 8. Here the displacement of the neutral axis of each the elliptical shell models is plotted as a function of the distance from the root (i.e., the span variable). Generally excellent agreement is shown between all the modelling procedures utilized in this work particularly with regard to the anti-derivative and the Runge-Kutta FD calculations; ANSYS-FE calculations though suffering from the aforementioned overestimation of stiffness, as more clearly indicated in the insert of Figure 8, these discrepancies are still well within the model convergence tolerance. However, it must therefore be inferred from the results presented in Figures 7 and 8 and Table 1 that the analytical anti-derivative model predicts the primary displacement field remarkably well; thus, implying that (16), (15), (18), and (21) are components of a mathematical analytical model of the structural response of a tapering elliptical beam shell which is offered for the first time in this work.
3.1. Elliptical Beam Failure Predictions
Of particular interest from a design viewpoint are the resulting strain-fields shown in Figure 9. The analytical (anti-derivative) method strain-field is being found from application of (22), while the Runge-Kutta values were found from double numerical differentiation of the numerical solutions shown in Figure 8. The ANSYS-FE strain-field data were obtained from a path on the top surface connecting an apex of the ellipse minor axis at the root to an apex on the minor axis at the tip. The models show more discrepancies than the predication of the displacement field (Figure 8). In particular the data from ANSYS begins to diverge from that obtained from the other modelling procedures at about 1.6 m from the root. Notwithstanding the initial inaccuracies of the FE data at the root due to Saint-Venant edge effects, the FE strain predictions remain continuously higher than all other similar model calculations. Both Runge-Kutta numerical finite difference calculations remain in agreement with analytical anti-derivative model over a greater distance from the root than the FE data, with divergence occurring at a value of 1.85 m. Saint-Venant effects are much more prevalent at the tip of the FE model where a difference of over 54% is evident when compared to the FD shell formulation solutions. Significantly none of the modelling procedures predict maxima strain values, where one may expect, at the root of the beam as shown in Table 2, which may indicate the position of delamination and eventual failure of the structure. The positions shown in Table 2 indicate potential failure of the beam just below 50% of the span. An apparent underapproximation from the Runge-Kutta methods being explained by the lack of solution points over the most significant changes in slope of the beam.
3.2. Idealized Blade Modelling
The normalized displacement field as a function of proportion of span of the idealized blade without and with a stiffening rib is shown, respectively, in Figures 11(a) and 11(b). In absence of the the stiffening component, Figure 10(a) is in good agreement up until about 40% of the span. Thereafter the agreement deteriorates the ANSYS-FE models predicting smallest magnitudes of the displacement field which throughout are in close synergy with the Runge-Kutta FD shell formulation model predictions. The lower predictions of the displacement by the FE model can again be attributed to the convergence tolerance together with the overestimation of the stiffness of the numerical method. The Runge-Kutta solution using the solid formulation rendered the greatest magnitudes in displacement fields; this can be attributed to the solid formulation stiffness being less than that of the shell analogues. It should also be noted that all modelling predictions are higher than the 10% threshold dictated by a number of design codes (e.g., [7, 8]), where the maximum values of the displacement are predicted as 11.2% (470 mm), 12% (504 mm), and 13% (564 mm) for the FE, FD-Shell, and FD-Solid models, respectively, noting that the absolute difference between the modelling predictions is as much as 94 mm, that is an overall 20% difference.
(a)
(b)
(a)
(b)
Figure 11(b) shows the displacement results from each of the models in the presence of a 16 mm thick stiffening web, where similar trends are observed as in Figure 10(a) with the largest magnitude being predicted by the Runge-Kutta solution to the solid formulation, followed, respectively, by solution to the FD and ANSYS-FE shell formulations. The web thickness is in fact an optimum value based on the aforementioned deflexion design criterion of 10% of the span obtained via successive bisection of the maximum magnitude of the deflexion of all models that is the Runge-Kutta solution to the solid formulation. It should be noted that the optimization routine established the required thickness. It took between three to five iterations to obtain and was very much more efficient than the more generalized ANSYS equivalent algorithm. Hence the modelling methods predict upper and lower expectations of the displacement field and therefore the maximum blade deflexion in the case of Figure 10(b) this is between 355 mm (8.46% of the span) and 420 mm (10% of the span).
The strain field and hence failure predictions of the idealized blade, without the stiffening web, are shown in Figure 11(a), while the equivalent field with the stiffening rib is shown in Figure 11(b). In both cases the maximum strain field values are predicted from the Runge-Kutta solution of the solid formulation, which is expected as this particular solution also predicted greater magnitude displacements as shown earlier in Figure 10.
Very close agreement is observed for the strain field predictions of both FE and FD solutions to the shell formulation, for the aerofoil model without the stiffening rib, neglecting the Saint-Venant effects evident in the ANSYS-FE model data. In this case the field values being less than the solid formulation by almost a constant amount of about 0.03% until maxima values are realized. As with the elliptical tapering beam, Figure 8, the maximum strain and hence possible failure position do not occur at the root. In fact the maximum strain value of 0.42% is predicted by the Runge-Kutta solution to the solid formulation and occurs at 54% of the span (2.27 m). The FE solution from the ANSYS code predicts the maximum strain value of 0.39 at 52% of the span (2.78 m), while identical values were predicted from the Runge-Kutta solution to the shell formulation. In the presence of the stiffening rib, Figure 11(b), as expected the magnitude of the strain field values is reduced by just under 0.1%. The maximum strain value of 0.32% is being predicted by the Runge-Kutta solution to the solid formulation which occurs at 56% of the span. The solution to the shell formulation models predict very similar maxima of 0.30 and 0.31 for, respectively, the FE and FD, occurring at 52% and 56% of the span. We note here that according to the design codes [8] these strain values are above the tensile 0.3% strain and the compressive 0.2% strain. This implies that failure of the structure will occur at about the halfway point not at root. From a structural design viewpoint this implies that the idealized blade is not fully stressed [11], or more correctly, since the fibre-glass fails though maximum strain, in this case fully strained. In order to minimize the weight of the structure it would be necessary to vary the dimensions of the cross section so as to give a blade of constant strength. However, this idealized blade is not in this favorable condition due to the aforementioned aerodynamic considerations coupled with the constant shell thickness. Nonetheless, especially the FE modelling goal driven design optimization methods employed throughout, the work presented in this paper could be used to drive the idealized model toward this optimum condition.
3.3. Experimental Verification
As a final verification of the methods used calculations were compared with empirical data from the reference [17]. Here a model was constructed which approximated the FX 66-S-196 aerofoil at the root and NACA 63-621 aerofoil at the tip. Here the chord and thickness values at the root were, respectively, 600 mm and 178 mm, with tip dimension analogues being 200 mm and 25 mm. The modelling method of this work then naturally matched the aerofoils; the average profile shell thickness and mechanical properties are being employed, which guaranteed the required smooth transition from the FX to NACA profile. The loading conditions were amended to model the proof load test detailed in the reference [17]. Here the blade is in the usual flap-wise position (Figure 6) from the clamped root flange to represent the vertical structural columns at the factory building and loaded around 2/3 of its length from the root. The experimental tip position was recorded first under the blade's own weight and then with the platform which carried the load. The load was applied perpendicular to the chord direction in increments of 270 N until the 6.4 kN load was reached making a total load of 7.4 kN with platform and ropes.
Comparisons between Runge-Kutta FD solutions and the solid and shell formulations with the experimental data of [17] are shown in Figure 12, where a linear static response is evident.
Considering the amount of engineering and geometrical assumptions the agreement between the modelling and experimental data is quite remarkable. At low forces the solution to the shell formulation underestimates the structural stiffness when compared to the empirical data as evidenced by larger displacement being predicted by the model for identical forces. Moreover, there appears to be evidence to support that this structural stiffness of the shell is load dependent since after a load of about 5.8 kN this modelling formulation begins to underestimate the experimental data for any given load. Though this possible phenomenon is also present, to less of an extent in the solid formulation data, this model appears to always overestimate the structural stiffness when compared to the experimental data presented.
All the empirical data fall within the 95% error bars of the modelling data. These error bars being evaluated from the average sample standard deviations of the Runge-Kutta and ANSYS-FE displacement-fields for each of the given loads and multiplying the result by the relevant -value (). The results therefore demonstrate that modelling methods can predict the static performance of aerodynamic components such as wind-turbine blades to within 5% of real manufactured components. This implies that the small discrepancies between geometrical and/or intrinsic material properties enforced by the engineering modelling assumptions are insignificant. This can be explained in part by inconsistencies in geometry (e.g., twist), lay-up and/or fibre/matrix volume fractions, and/or changes in intrinsic constituents of the composite materials used for construction.
4. Conclusions
Mathematical solid and shell formulations of composite aerodynamic components have been presented and solved using analytical, finite difference, and finite element methods. These models were shown to predict the displacement, strain-fields, and hence failure of real wind-turbine blades in service. Though the models presented throughout allow only for changes in cross-section, their robust formulations allow for the variation of many other relevant parameters. Noting that the FD Runge-Kutta solution codes are freely available (e.g., via the freeware SciLab: https://www.scilab.org/) the models presented throughout this work can be employed as the basis of an expert system type design code; indeed, with relatively minor amendments, even fatigue and modal analyses being possible. The more specific conclusions from the work presented in this paper being as follows.(i)Solution to an elliptical shell beam has been offered for the first time as a response to a lack of treatment in the literature. The solution has served as a salient benchmark in order to compare with well-established FD methods and commercial ANSYS and design evaluation codes.(ii)All model calculations of strain fields indicate that failure does not occur at the root of the beam/blade for all of the geometries investigated, meaning that the components under investigation were not fully strained (stressed) and therefore not yet optimized from a purely structural viewpoint.(iii)In order to conform to a 10% displacement design criterion the aerofoil section approximations investigated that NACA 63-621 with average shell thickness of 5 mm requires a stiffening rib that is required in all cases. Based on this criterion a successive bisection optimization method reveals that the thickness of this web should be not less than 16 mm.(iv)Solution and optimization protocols obtained via analytical method and/or Runge-Kutta methods presented throughout this work are generally orders of magnitude faster than commercial FE optimization routines employed by ANSYS. We note here that the ANSYS design optimization routines are general purpose while the one presented in this paper is specific to this problem and hence no generality is inferred.(v)Though modelling of idealized turbine blades can show up to 20% difference with regard to predictions of the maximum displacement values, these methods show remarkable agreement with empirical data considering the underlying assumptions.
Further work in now progressing to use the modelling methods presented here in order to investigate variations in composite material lay-up, hence mechanical properties, shell thickness, and loading conditions (from CFD and empirical data), with position.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.