• Keine Ergebnisse gefunden

Gerges Edwar Mehanny Beshay and Karam Yousef Maalawi

C. Applied loads and stress analysis

The distributed load vectors are expressed in the undeformed coordinates (see Figure C.1), as follows:

Distributed forces:P¼PAþPIþPGþPD (C.1) Distributed moments:q¼qAþqIþqGþqD (C.2) where the subscripts A, I, G, and D refer to the aerodynamic, inertial, gravita-tional, and damping contributions, respectively. The aerodynamic forces PAand moments qAcan be obtained using the quasi-steady blade-element strip

theory [1, 2].

Argument Description

fun The function to be minimized x0, lb., and

ub

Starting, lower, and upper boundaries of design variables

nonlcon The function that computes the nonlinear inequality constraints “c(x)0” and the nonlinear equality constraints “ceq(x) = 0”

fval Value of the function “fun” at x

exitflag Integer identifying the reasons causing the algorithm terminated grad Gradient at x

hessian Hessian at x

lambda Structure containing the Lagrangian multipliers at the solution x

output Structure containing information about the optimization process such as number of function evaluations and iterations taken and the used optimization technique

Table B.1.

Arguments related to fmincon function [35].

Figure C.2 shows the velocity triangle and the coordinate system used for computing aerodynamic loads. Vyrand Vzrare the tangential and axial velocity components, respectively. The resultant velocity Vrcan be calculated from:

Vr≃ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V2yrþV2zr

q ¼ �Vyr

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1þ Vzr

Vyr

� �2

s

(C.3) The distributed lift and drag force vectors on an arbitrary airfoil section are given by the aerodynamic formulas (see Figure C-2):

Lift L¼1

2 ρaV2r C CL�sin φi^j0þcos φi^k0� DragD¼1

2 ρaV2rC CD��cos φi^j0þsin φi^k0� (C.4) whereρa= air density; C = local chord at spanwise location r; CL= CL(α) = CLαα, lift coefficient; CD= CD(α), drag coefficient; CLα= lift – curve slope;α= angle of attack;φi¼inflow angle¼θBþα.

Figure C.1.

Deformed and undeformed configurations of a blade segment.

optimization problems is called the optimization toolbox [35]. MATLAB optimiza-tion toolbox contains a library of programs or m-files, which can be used for the solution of different optimization problems. A commonly implemented function denoted by fmincon is applied to most constrained objective functions. It is the most suitable function to the optimization problem of this investigation, since it uses sequential quadratic programming (SQP) method.

fmincon starts at “x0” and attempts to find a minimizer “x” of the objective function described in the m-file named “fun” subject to the linear inequalities

Axb” and the linear equalities “Aeqx¼beq.” If there are no linear equali-ties or inequaliequali-ties, A, b, Aeq, and beq are replaced with “[].” fmincon can define lower and upper bounds on the design variables in “x” so that the solution is always in the range “lb≤x≤ub.” fmincon can subject the minimization process to the nonlinear inequalities “c(x)” or equalities “ceq(x)” defined in the m-file constraint function named “nonlcon.” fmincon optimizes such that “c(x)≤0” and “ceq(x) = 0.” fmincon minimizes with the optimization options involved in the structure

options.” Table B.1 defines all input and output arguments related to fmincon optimization function.

C. Applied loads and stress analysis

The distributed load vectors are expressed in the undeformed coordinates (see Figure C.1), as follows:

Distributed forces:P¼PAþPIþPGþPD (C.1) Distributed moments:q¼qAþqIþqGþqD (C.2) where the subscripts A, I, G, and D refer to the aerodynamic, inertial, gravita-tional, and damping contributions, respectively. The aerodynamic forces PAand moments qAcan be obtained using the quasi-steady blade-element strip

theory [1, 2].

Argument Description

fun The function to be minimized x0, lb., and

ub

Starting, lower, and upper boundaries of design variables

nonlcon The function that computes the nonlinear inequality constraints “c(x)0” and the nonlinear equality constraints “ceq(x) = 0”

fval Value of the function “fun” at x

exitflag Integer identifying the reasons causing the algorithm terminated grad Gradient at x

hessian Hessian at x

lambda Structure containing the Lagrangian multipliers at the solution x

output Structure containing information about the optimization process such as number of function evaluations and iterations taken and the used optimization technique

Table B.1.

Arguments related to fmincon function [35].

Figure C.2 shows the velocity triangle and the coordinate system used for computing aerodynamic loads. Vyrand Vzrare the tangential and axial velocity components, respectively. The resultant velocity Vrcan be calculated from:

Vr≃ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V2yrþV2zr

q ¼ �Vyr

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1þ Vzr

Vyr

� �2

s

(C.3) The distributed lift and drag force vectors on an arbitrary airfoil section are given by the aerodynamic formulas (see Figure C-2):

Lift L¼1

2 ρaV2rC CL�sin φi^j0þ cosφi^k0� DragD¼1

2 ρaV2rC CD��cosφi^j0þ sinφi^k0� (C.4) whereρa= air density; C = local chord at spanwise location r; CL= CL(α) = CLαα, lift coefficient; CD= CD(α), drag coefficient; CLα= lift – curve slope;α= angle of attack;φi¼inflow angle¼θBþα.

Figure C.1.

Deformed and undeformed configurations of a blade segment.

Expressed in the undeformed system, the aerodynamic force vector per unit length of the blade is therefore given as:

PA¼LþD¼PxA^iþPyA^jþPzA^k (C.5) The distributed inertia loads per unit blade span are obtained by applying D’Alembert’s principle as follows:

Inertial force vector:PIB¼ � ð ð

ρBrBdy dz Inertial moment vector:qIB¼ �

ð ð

ρBr0Br0E:C

�€rB

� �

dy dz

(C.6)

whereρBis the blade material mass density,€rBis the acceleration vector, and r0B and r0E:Care the position vectors of an arbitrary point and the elastic center of the blade cross section, respectively.

Next, considering the equilibrium of a differential element of the deformed blade, the equilibrium equations expressed in the rotating (xyz) axes are:

Force equilibrium: ∂F

∂xþP¼0 Moment equilibrium: ∂M

∂x þ^i0Fþq¼0

(C.7)

Internal force vector: F¼Fx^iþFy^jþFz^k

Internal moment vector: M¼Mx^iþMy^jþMz^k (C.8) A wind turbine blade cross section is composed of thin-walled, closed, single or multicellular composite beams. The periphery of each beam cross section is assumed to be constructed of flat composite laminates as shown in Figure C.5. The stress-strain relations of the composite laminates which discretize each cross section are computed using the classical laminate theory CLT. Although each laminate is

Figure C.2.

Velocity triangle of an arbitrary airfoil section in the deformed state.

actually an assembly of multiple layers with different constitutive properties, CLT is used to calculate a set of effective stiffness coefficients that allows a composite laminate to be treated as a single structural element [12]. Therefore the blade cross section may be considered to be composed of discrete sections of homogenous materials.

Generally, for a heterogeneous composite section shown in Figure C.3 and C.4, the modulus-weighted cross-sectional properties are defined as [12]:

Figure C.4.

Cross section for a heterogeneous composite beam.

Figure C.5.

Resultant loads and moments applying to a composite laminate.

Figure C.3.

Blade cross section with discretized composite laminated plates.

Expressed in the undeformed system, the aerodynamic force vector per unit length of the blade is therefore given as:

PA¼LþD¼PxA^iþPyA^jþPzA^k (C.5) The distributed inertia loads per unit blade span are obtained by applying D’Alembert’s principle as follows:

Inertial force vector:PIB¼ � ð ð

ρBrBdy dz Inertial moment vector:qIB¼ �

ð ð

ρBr0Br0E:C

�€rB

� �

dy dz

(C.6)

whereρBis the blade material mass density,€rBis the acceleration vector, and r0B and r0E:Care the position vectors of an arbitrary point and the elastic center of the blade cross section, respectively.

Next, considering the equilibrium of a differential element of the deformed blade, the equilibrium equations expressed in the rotating (xyz) axes are:

Force equilibrium: ∂F

∂xþP¼0 Moment equilibrium: ∂M

∂x þ^i0Fþq¼0

(C.7)

Internal force vector: F¼Fx^iþFy^jþFz^k

Internal moment vector: M¼Mx^iþMy^jþMz^k (C.8) A wind turbine blade cross section is composed of thin-walled, closed, single or multicellular composite beams. The periphery of each beam cross section is assumed to be constructed of flat composite laminates as shown in Figure C.5. The stress-strain relations of the composite laminates which discretize each cross section are computed using the classical laminate theory CLT. Although each laminate is

Figure C.2.

Velocity triangle of an arbitrary airfoil section in the deformed state.

actually an assembly of multiple layers with different constitutive properties, CLT is used to calculate a set of effective stiffness coefficients that allows a composite laminate to be treated as a single structural element [12]. Therefore the blade cross section may be considered to be composed of discrete sections of homogenous materials.

Generally, for a heterogeneous composite section shown in Figure C.3 and C.4, the modulus-weighted cross-sectional properties are defined as [12]:

Figure C.4.

Cross section for a heterogeneous composite beam.

Figure C.5.

Resultant loads and moments applying to a composite laminate.

Figure C.3.

Blade cross section with discretized composite laminated plates.

Ac¼ 1 Eref

Xn

i¼1

EiAi yc¼ 1

ErefAc

Xn

i¼1

EiAiyi zc¼ 1

ErefAc

Xn

i¼1

EiAizi

Iy¼ 1 Eref

Xn

i¼1

EiIvo,iþAiz2iIz¼ 1

Eref Xn

i¼1

EiIwo,iþAiy2iIyz¼ 1

Eref Xn

i¼1

EiIvwo,iþAiyizi

(C.9)

where Eref is a reference modulus of elasticity, (yi,zi) denotes the geometric centroid of each discrete element of the cross section, and (vo,i, wo,i) denotes the principal axes of each discrete element.

The parallel axes theorem can be applied to compute the second moments of area about the cross-sectional principal axes as follows:

Iv¼IyAcð Þzc 2

Iw¼IzAc� �yc 2 Ivw¼IyzAcyczc

(C.10)

Once the global cross-sectional properties are computed using the method of modulus-weighted properties, the effective axial stress applied to the blade cross section can be given by:

σxðy, zÞ ¼Fx

AcMzIvþMyIvw

IvIwIvw2yyc

þMyIwþMzIvw

IvIwIvw2 ðzzcÞ (C.11) Finally, by converting the distribution of the effective beam stresses into equiv-alent in-plane distributed loads on the flat laminates composing the cross section, as shown in Figure C.5, the lamina-level strains and stresses can be computed using the CLT.

Author details

Gerges Edwar Mehanny Beshay1* and Karam Yousef Maalawi2

1 Faculty of Engineering (Shoubra), Department of Mechanical Engineering, Benha University, Cairo, Egypt

2 Department of Mechanical Engineering, National Research Centre, Dokki, Giza, Egypt

*Address all correspondence to: gerges.beshay@feng.bu.edu.eg

© 2020 The Author(s). Licensee IntechOpen. Distributed under the terms of the Creative Commons Attribution - NonCommercial 4.0 License (https://creativecommons.org/

licenses/by-nc/4.0/), which permits use, distribution and reproduction for non-commercial purposes, provided the original is properly cited. –NC

Ac¼ 1 Eref

Xn

i¼1

EiAi yc¼ 1

ErefAc

Xn

i¼1

EiAiyi zc¼ 1

ErefAc

Xn

i¼1

EiAizi

Iy¼ 1 Eref

Xn

i¼1

EiIvo,iþAiz2iIz¼ 1

Eref Xn

i¼1

EiIwo,iþAiy2iIyz¼ 1

Eref Xn

i¼1

EiIvwo,iþAiyizi

(C.9)

where Eref is a reference modulus of elasticity, (yi,zi) denotes the geometric centroid of each discrete element of the cross section, and (vo,i, wo,i) denotes the principal axes of each discrete element.

The parallel axes theorem can be applied to compute the second moments of area about the cross-sectional principal axes as follows:

Iv¼IyAcð Þzc 2

Iw¼IzAc� �yc 2 Ivw¼IyzAcyczc

(C.10)

Once the global cross-sectional properties are computed using the method of modulus-weighted properties, the effective axial stress applied to the blade cross section can be given by:

σxðy, zÞ ¼Fx

AcMzIvþMyIvw

IvIwIvw2yyc

þMyIwþMzIvw

IvIwIvw2 ðzzcÞ (C.11) Finally, by converting the distribution of the effective beam stresses into equiv-alent in-plane distributed loads on the flat laminates composing the cross section, as shown in Figure C.5, the lamina-level strains and stresses can be computed using the CLT.

Author details

Gerges Edwar Mehanny Beshay1* and Karam Yousef Maalawi2

1 Faculty of Engineering (Shoubra), Department of Mechanical Engineering, Benha University, Cairo, Egypt

2 Department of Mechanical Engineering, National Research Centre, Dokki, Giza, Egypt

*Address all correspondence to: gerges.beshay@feng.bu.edu.eg

© 2020 The Author(s). Licensee IntechOpen. Distributed under the terms of the Creative Commons Attribution - NonCommercial 4.0 License (https://creativecommons.org/

licenses/by-nc/4.0/), which permits use, distribution and reproduction for non-commercial purposes, provided the original is properly cited. –NC

References

[1]Al-Bahadly I. Wind Turbines.

Croatia: IntechOpen; 2011. Available from: http://www.intechopen.com/

books/wind-turbines. ISBN: 978-953-307-221-0

[2]Gasch R, Twele J. Wind Power Plants. New York: Springer; 2012. ISBN:

978-3-642-22938-1

[3]Gay D, Hoa SV. Composite Materials:

Design and Applications. New York:

CRC Press; 2007

[4]Daniel I, Ishai O. Engineering Mechanics of Composite Materials. 2nd ed. New York: Oxford University Press;

2006. pp. 432. ISBN: 978-0195150971 [5]Suresh S, Mortensen A.

Fundamentals of Functionally Graded Materials. Cambridge University Press;

1998

[6]Birman V, Byrd WL. Modeling and analysis of functionally graded materials and structures. Applied Mechanics Reviews. 2007;60(5):195-216 [7]Maalawi KY, Negm HM. Optimal frequency design of wind turbine blades. Journal of Wind Engineering and Industrial Aerodynamics. 2002;

90(8):961-986

[8]Maalawi KY. A model for yawing dynamic optimization of a wind turbine structure. International Journal of Mechanical Sciences. 2007;49:1130-1138 [9]Librescu L, Maalawi K. Material grading for improved aeroelastic stability in composite wings. Journal of Mechanics of Materials and Structures.

2007;2(7):1381-1394

[10]Chen K-N, Chen P-Y. Structural optimization of 3 MW wind turbine blades using a two- step procedure.

International Journal for Simulation and

Multidisciplinary Design Optimization.

2010;4:159-165. DOI: 10.1051/ijsmdo/

2010020

[11]Maalawi KY, Badr MA. Frequency optimization of a wind turbine blade in pitching motion. Journal of Power and Energy, Proceedings of the Institution of Mechanical Engineers, Part A. 2010;

224:545-554. DOI: 10.1243/

09576509JPE907

[12]Sale D, Aliseda A, Motley M, Li Y.

Structural optimization of composite blades for wind and hydrokinetic turbines. In: Proceeding of 1st Marine Energy Technology Symposium, METS 2013. Washington, DC, April 10–11.

2013

[13]Chehouri A, Younes R, Ilinca A, Perron J, Lakiss H. Optimal design for a composite wind turbine blade with fatigue and failure constraints.

Transactions of the Canadian Society for Mechanical Engineering, CSME-16.

2015;39(2):171-186

[14]Maalawi KY. Dynamic optimization of functionally graded thin-walled box beams. International Journal of

Structural Stability and Dynamics. 2017;

17(9):1-24

[15]Armanios EA, Badir AM. Free vibration analysis of anisotropic thin-walled closed-sections beams. American Institute of Aeronautics and

Astronautics (AIAA). 1995;33:1905-1910 [16]Dancila DS, Armanios EA. The influence of coupling on the free vibration of anisotropic thin-walled closed-section beams. International Journal of Solids and Structures. 1998;

35:3105-3119

[17]Durmaz S, Kaya MO. Free vibration of an anisotropic thin-walled box beam under bending-torsion coupling.

In: Proceedings of 3rd International Conference On Integrity, Reliability and Failure. 2009

[18]Shadmehri F, Haddadpour H, Kouchakzageh M. Flexural-torsional behaviour of thin-walled composite beams with closed cross-section. Thin Walled Structures. 2007;45:699-705 [19]Phuong T, Lee J. Flexural-torsional behavior of thin-walled composite box beams using shear-deformable beam theory. Engineering Structures. 2008;

30:1958-1968

[20]Piovan MT, Filipicha CP,

Cortínez VH. Coupled free vibration of tapered box beams made of composite materials. Mecánica Computacional.

2006;25:1767-1779

[21]Kargarnovin MH, Hashemi M. Free vibration analysis of multilayered composite cylinder consisting fibers with variable volume fraction.

Composite Structures. 2012;94:931-944 [22]Liu Y, Shu DW. Free vibration analysis of exponential functionally graded beams with a single

delamination. Composites Part B Engineering. 2014;59:166-172 [23]Whitney JM. Elastic moduli of unidirectional composites with anisotropic filaments. Journal of Composite Materials. 1967;1:188-194 [24]Meirovitch L. Principles and Techniques of Vibrations. Englewood Cliffs, NJ: Prentice-Hall; 1997

[25]Mihail B, Valentin C, Costin D. A transfer matrix method for free vibration analysis of Euler-Bernoulli beams with variable cross section.

Journal of Vibration and Control. 2014:

1-12. DOI: 10.1177/1077546314550699 [26]Halpin JC, Tsai SW. Effect of environmental factors on composite

materials. In: AFML-TR. Vol. 67. 1969. pp. 1-53

[27]Rao S, Singiresu. Engineering Optimization: Theory and Practice. John Wiley & Sons Inc; 2009

[28]Maalawi K, Badr M. Design optimization of mechanical elements and structures: A review with

application. Journal of Applied Sciences Research. 2009;5(2):221-231

[29]Fletcher R. The sequential quadratic programming method. In: Nonlinear Optimization. Springer; 2010. pp. 165-214

[30]Boggs PT, Tolle JW. Sequential quadratic programming for large-scale nonlinear optimization. Journal of Computational and Applied Mathematics. 2000;124(1):123-137 [31]Bir GS, Oyague F. Estimation of blade and tower properties for the gearbox research collaborative wind turbine. In: Technical Report NREL/TP-500-42250. U.S. Department of Energy: National Renewable Energy Laboratory; November 2007

[32] Khennane A. Introduction to Finite Element Analysis Using MATLAB and Abaqus. Boca Raton: CRC Press; 2013 [33]Available from: http://www.plm. automation.siemens.com/en_us/ products/velocity/femap/features/ [34]Bryan H. Engineering Composite Materials. The Institute of Materials; 1999

[35]Grace A, Branch MA, Coleman T. Optimization Toolbox for use with MATLAB. The Math Works Inc.; 1999

References

[1]Al-Bahadly I. Wind Turbines.

Croatia: IntechOpen; 2011. Available from: http://www.intechopen.com/

books/wind-turbines. ISBN: 978-953-307-221-0

[2]Gasch R, Twele J. Wind Power Plants. New York: Springer; 2012. ISBN:

978-3-642-22938-1

[3]Gay D, Hoa SV. Composite Materials:

Design and Applications. New York:

CRC Press; 2007

[4]Daniel I, Ishai O. Engineering Mechanics of Composite Materials. 2nd ed. New York: Oxford University Press;

2006. pp. 432. ISBN: 978-0195150971 [5]Suresh S, Mortensen A.

Fundamentals of Functionally Graded Materials. Cambridge University Press;

1998

[6]Birman V, Byrd WL. Modeling and analysis of functionally graded materials and structures. Applied Mechanics Reviews. 2007;60(5):195-216 [7]Maalawi KY, Negm HM. Optimal frequency design of wind turbine blades. Journal of Wind Engineering and Industrial Aerodynamics. 2002;

90(8):961-986

[8]Maalawi KY. A model for yawing dynamic optimization of a wind turbine structure. International Journal of Mechanical Sciences. 2007;49:1130-1138 [9]Librescu L, Maalawi K. Material grading for improved aeroelastic stability in composite wings. Journal of Mechanics of Materials and Structures.

2007;2(7):1381-1394

[10]Chen K-N, Chen P-Y. Structural optimization of 3 MW wind turbine blades using a two- step procedure.

International Journal for Simulation and

Multidisciplinary Design Optimization.

2010;4:159-165. DOI: 10.1051/ijsmdo/

2010020

[11]Maalawi KY, Badr MA. Frequency optimization of a wind turbine blade in pitching motion. Journal of Power and Energy, Proceedings of the Institution of Mechanical Engineers, Part A. 2010;

224:545-554. DOI: 10.1243/

09576509JPE907

[12]Sale D, Aliseda A, Motley M, Li Y.

Structural optimization of composite blades for wind and hydrokinetic turbines. In: Proceeding of 1st Marine Energy Technology Symposium, METS 2013. Washington, DC, April 10–11.

2013

[13]Chehouri A, Younes R, Ilinca A, Perron J, Lakiss H. Optimal design for a composite wind turbine blade with fatigue and failure constraints.

Transactions of the Canadian Society for Mechanical Engineering, CSME-16.

2015;39(2):171-186

[14]Maalawi KY. Dynamic optimization of functionally graded thin-walled box beams. International Journal of

Structural Stability and Dynamics. 2017;

17(9):1-24

[15]Armanios EA, Badir AM. Free vibration analysis of anisotropic thin-walled closed-sections beams. American Institute of Aeronautics and

Astronautics (AIAA). 1995;33:1905-1910 [16]Dancila DS, Armanios EA. The influence of coupling on the free vibration of anisotropic thin-walled closed-section beams. International Journal of Solids and Structures. 1998;

35:3105-3119

[17]Durmaz S, Kaya MO. Free vibration of an anisotropic thin-walled box beam under bending-torsion coupling.

In: Proceedings of 3rd International Conference On Integrity, Reliability and Failure. 2009

[18]Shadmehri F, Haddadpour H, Kouchakzageh M. Flexural-torsional behaviour of thin-walled composite beams with closed cross-section. Thin Walled Structures. 2007;45:699-705 [19]Phuong T, Lee J. Flexural-torsional behavior of thin-walled composite box beams using shear-deformable beam theory. Engineering Structures. 2008;

30:1958-1968

[20]Piovan MT, Filipicha CP,

Cortínez VH. Coupled free vibration of tapered box beams made of composite materials. Mecánica Computacional.

2006;25:1767-1779

[21]Kargarnovin MH, Hashemi M. Free vibration analysis of multilayered composite cylinder consisting fibers with variable volume fraction.

Composite Structures. 2012;94:931-944 [22]Liu Y, Shu DW. Free vibration analysis of exponential functionally graded beams with a single

delamination. Composites Part B Engineering. 2014;59:166-172 [23]Whitney JM. Elastic moduli of unidirectional composites with anisotropic filaments. Journal of Composite Materials. 1967;1:188-194 [24]Meirovitch L. Principles and Techniques of Vibrations. Englewood Cliffs, NJ: Prentice-Hall; 1997

[25]Mihail B, Valentin C, Costin D. A transfer matrix method for free vibration analysis of Euler-Bernoulli beams with variable cross section.

Journal of Vibration and Control. 2014:

1-12. DOI: 10.1177/1077546314550699 [26]Halpin JC, Tsai SW. Effect of environmental factors on composite

materials. In: AFML-TR. Vol. 67. 1969.

pp. 1-53

[27]Rao S, Singiresu. Engineering Optimization: Theory and Practice. John Wiley & Sons Inc; 2009

[28]Maalawi K, Badr M. Design optimization of mechanical elements and structures: A review with

application. Journal of Applied Sciences Research. 2009;5(2):221-231

[29]Fletcher R. The sequential quadratic programming method. In: Nonlinear Optimization. Springer; 2010.

pp. 165-214

[30]Boggs PT, Tolle JW. Sequential quadratic programming for large-scale nonlinear optimization. Journal of Computational and Applied Mathematics. 2000;124(1):123-137 [31]Bir GS, Oyague F. Estimation of blade and tower properties for the gearbox research collaborative wind turbine. In: Technical Report NREL/TP-500-42250. U.S. Department of Energy:

National Renewable Energy Laboratory;

November 2007

[32] Khennane A. Introduction to Finite Element Analysis Using MATLAB and Abaqus. Boca Raton: CRC Press; 2013 [33]Available from: http://www.plm.

automation.siemens.com/en_us/

products/velocity/femap/features/

[34]Bryan H. Engineering Composite Materials. The Institute of Materials;

1999

[35]Grace A, Branch MA, Coleman T.

Optimization Toolbox for use with MATLAB. The Math Works Inc.; 1999

Aerodynamic, Structural and