Open access peer-reviewed chapter

Optimization of Functionally Graded Material Structures: Some Case Studies

Written By

Karam Maalawi

Submitted: May 28th, 2018 Reviewed: November 6th, 2018 Published: November 26th, 2018

DOI: 10.5772/intechopen.82411

Chapter metrics overview

1,624 Chapter Downloads

View Full Metrics


This chapter focuses on some of the most advances made in the field of stability, dynamic, and aeroelastic optimization of functionally graded composite structures. Practical realistic optimization models using different strategies for measuring structural performance are presented and discussed. The selected design variables include the volume fractions of the composite material constituents as well as geometrical and cross-sectional parameters. The mathematical formulation is based on dimensionless quantities; therefore, the analysis can be valid for different configurations and sizes. Such normalization has led to a naturally scaled optimization model, which is favorable for most optimization techniques. Case studies include structural dynamic optimization of thin-walled beams in bending motion, optimization of drive shafts against torsional buckling and whirling, and aeroelastic optimization of subsonic aircraft wings. Other stability problems concerning fluid-structure interaction has also been addressed. Several design charts that are useful for direct determination of the optimal values of the design variables are introduced. The proposed mathematical models have succeeded in reaching the required optimum solutions, within reasonable computational time, showing significant improvements in the overall structural performance as compared with reference or known baseline designs.


  • functionally graded materials
  • composite structures
  • optimum design
  • buckling stability
  • structural dynamics
  • fluid-structure interaction

1. Introduction

Functionally graded materials (FGMs) are new generation of advanced composites that have gained interest in several engineering applications such as spacecraft heat shields, high-performance structural elements, and critical engine components [1, 2]. They are formed of two or more constituent phases with a continuously variable composition producing properties that change spatially within a structure. FGMs possess a number of advantages that make them attractive in improving structural performance, such as maximized torsional rigidity of composite shafts [3], improved residual stress distribution and enhanced thermal properties [4], higher natural frequencies of composite beams [5], and broader aeroelastic stability boundaries of aircraft wings [6]. Actually, the concept of FGMs was originated in Japan in 1984 during a space project, in the form of proposed thermal barrier material capable of withstanding high-temperature gradients.

Figure 1 shows variation of the volume fraction through the thickness of a plate fabricated from ceramic and metal. Ceramic provides high-temperature resistance because of its low thermal conductivity, while metal secures the necessary strength and stiffness. FGMs may also be developed using fiber-reinforced layers with a volume fraction of fibers that is coordinate dependent, rather than constant, producing favorable properties or response [6]. In this chapter, much attention is given to fibrous-type and their constitutive relationships.

Figure 1.

FGM ceramic/metal particulate composite with volume fraction graded in the vertical direction [1].

An excellent review paper dealing with the basic knowledge and various aspects on the use of FGMs and their wide applications is given by Birman and Byrd [7], who presented comprehensive discussions of the development related to stability and dynamic of FGM structures. Closed-form expressions for calculating the natural frequencies of an axially graded beam were derived in [8], where the modulus of elasticity was taken as a polynomial of the axial coordinate along the beam’s length. An inverse problem was solved to find the stiffness and mass distributions so that the chosen polynomial serves as an exact mode shape. Qian and Batra [9] considered frequency optimization of a cantilevered plate with variable volume fraction according to simple power laws. Genetic algorithm was implemented to find the optimum values of the power exponents, which maximize the natural frequencies. They concluded that the volume fraction needs to be varied in the longitudinal direction of the plate rather than in the thickness direction. Another work presented an analytical approach for designing efficient patterns of FGM bars having maximized frequencies while maintaining the total mass at a constant value [10]. The distribution of the volume fractions of the material constituents was optimized using either discrete or continuous variations along the bar length.

In the context of structural stability, Elishakoff and Endres [11] considered buckling of an axially graded cantilevered column and derived closed form solution for the mode shape and the critical load. A semi-inverse method was employed to obtain the spatial distribution of the elastic modulus in the axial direction. In Ref. [12], the buckling of simply supported three-layer circular cylindrical shell under axial compressive load was analyzed. The middle layer sandwiched with two isotropic layers was made of an isotropic FGM whose Young’s modulus varies parabolically in the thickness direction. Classical shell theory was implemented under the assumption of very small thickness/radius and very large length/radius ratios. Numerical results showed that the buckling load increases with an increase in the average value of Young’s modulus of the middle layer. An exact method was given in [13] for obtaining column’s designs with the maximum possible critical buckling load under equality mass constraint. Both material and wall thickness grading in the axial direction have been applied to determine the required optimal solutions. Case studies and detailed results were given for the cases of simply supported and cantilevered columns. Another work by Maalawi [14] presented a mathematical model for enhancing the buckling stability of composite, thin-walled rings/long cylinders under external pressure using radial material grading. The main structure is made of multiangle fibrous laminated layups having different volume fractions within the individual plies. This produced a piecewise grading of the material and thickness in the radial direction. The critical buckling contours are plotted for different types of materials, showing significant improvement in the overall stability limits of the structure under the imposed mass constraint.

Considering dynamic aeroelasticity of FGM structures, Shin-Yao [15] investigated the effect of variable fiber spacing on the supersonic flutter of composite laminates using the finite element method and quasi-steady aerodynamic theory. The formulation of the location-dependent stiffness and mass matrices due to nonhomogeneous material properties was derived. This study first demonstrates the flutter analysis of composite laminates with variable fiber spacing. Numerical results show that the sequence of the natural mode may be altered, and the two natural frequencies may be close to each other because the fiber distribution may change the distributed stiffness and mass of the plate. Therefore, it may change the flutter coalescent modes, and the flutter boundary may be increased or decreased due to the variable fiber spacing. More detailed discussions on stability, dynamic, and aeroelasticity of FGM structures are outlined in Ref. [16]. The attained optimal solutions were determined by applying the global search techniques [17, 18], which construct a number of starting points and use a local solver, such as “fmincon” routine in the MatLab optimization toolbox [19]. Global search technique is distinguished with fast converging to the global optima even if it starts with a design point far from the optimum. The local solver “fmincon” uses the method of sequential quadratic programming (SQP), which has a theoretical basis related to the solution of a set of nonlinear equations using Newton’s method and applies Kuhn-Tucker conditions to the Lagrangian of the constrained optimization problem.

It is the main intend of this chapter to present some fundamental issues concerning design optimization of different types of functionally graded composite structures. Practical realistic optimization models using different strategies for enhancing structural dynamics, stability, as well as aeroelastic performance are presented and discussed. Case studies include frequency optimization of thin-walled box beams, optimal design of drive shafts against torsional buckling and whirling, and aeroelastic optimization of subsonic aircraft wings. Design of pipelines against flow-induced flutter and/or divergence has also been addressed. Several design charts that are useful for direct determination of the optimal values of the design variables are introduced. In all, the given mathematical models can be regarded as useful design tools, which may save designers from having to choose the values of some of their variables arbitrarily.


2. Mathematical modeling of material grading

There are different scenarios in modeling the spatial variation of material properties of a functionally graded structure. For example, Chen and Gibson [20] considered distributions represented by polynomial functions and applied Galerkin’s method to calculate the required polynomial coefficients from the resulting algebraic equations. They found that the linear variation of the volume fraction is a best fit with that predicted experimentally for selected composite beam specimens. Chi and Chung [21] studied the mechanical behavior of FGM plates under transverse loading, where a constant Poisson’s ratio and variable moduli of elasticity throughout the plate thickness were assumed. The volume fraction of the constituent materials was defined by simple power laws, and closed form solutions using Fourier series were given for the case of simply supported plates. In general, the distribution of the material properties in functionally graded structures may be designed by either continuous or piecewise variation of the volume fraction in a specified direction. The most commonly utilized distributions are summarized in what follows.

2.1. Thickness grading

The first early model of volume fraction variation through the thickness of a plate fabricated from ceramic and metal was considered in [20]. This volume fraction is based on the mixture of metal and ceramic and is an indicator of the material composition (volumetric wise) at any given location in the thickness. If the volume fraction of ceramic is defined as “v” then the volume fraction of metal is the remainder of the material, or (1 − v), assuming no voids are present. A typical example, which was considered by numerous researchers in the field [1- 4] assumed that the volume fraction “v” can be varied through the thickness coordinate z by the power law (refer to Figure 1):


where h is the plate thickness and p is a volume fraction exponent, which dictates the amount and distribution of ceramic in the plate. With higher values of p, the plate tends toward metal while lower values tend toward ceramic (p = 0.0: fully ceramic, p = ∞: fully metal). Accordingly, the distribution of the mechanical and physical properties of FGM can be defined in terms of the material constants of the constituent phases based on a selected power-law model. Designers can vary the p-value to tailor the FGM to specific applications. In case of fibrous composites, Eq. (1) ought to be modified to account for the limits imposed on the fiber volume fractions at ẑ=zlh=±1/2 for consideration of other strength requirements and/or manufacturing restrictions. The modified form can be expressed as follows [5]:


Another type of the power-law expression was utilized by Bedjilili et al. [22], who considered vibration of fibrous composite beams with a variable volume fraction through the thickness of the cross section, as shown in Figure 2. It was concluded that by varying the fiber volume fraction within the beam thickness to create a FGM, certain vibration characteristics are significantly affected. The utilized formula was given as:


Figure 2.

Thickness distribution of the fiber volume fraction in FGM beam, vf (0) = 40%, vf (1/2) = 60% [22].

2.2. Spanwise grading

Some researchers considered grading of the fiber volume fraction in the spanwise (longitudinal) direction of a composite plate. Librescu and Maalawi [6] investigated optimization of composite wings using the concept of material grading in the spanwise direction. Both continuous and discrete distributions of the fiber volume fractions were considered in the developed optimization models. The following power-law expression was implemented:


where vft and vfr are the fiber volume fractions at wing tip and root, respectively. Δf is called the tapering ratio of the fiber volume fraction. Figure 3 shows the different patterns of the fiber volume fraction distribution for different values of the power exponent p. Both configurations of fibers aligned in the transverse (chordwise) and in the longitudinal (spanwise) directions are shown. The volume fraction is constrained to lie between 25% and 75% in order not to violate other strength and manufacturing requirements.

Figure 3.

Spanwise grading of fibers in a fibrous composite plate [6]. (a) Fibers aligned in chordwise direction, (b) fibers aligned in spanwise direction, and (c) fiber volume fraction distribution for different power exponents.

A more general distribution, given in Eq. (5), was tried by Shih-Yao [15], who applied it successfully to investigate the effect of grading on the supersonic flutter of rectangular composite plates.


2.3. Determination of mechanical properties

A variety of approaches have been developed to predict the mechanical properties of fibrous composite materials [23]. The common approaches fall into the following general categories: mechanics of materials; numerical methods; variational approach; semiempirical formulas; experimental measurements. Mechanics of materials approach is based on simplifying assumptions of either uniform strain or uniform stress in the constituents. Its predictions can be adequate only for longitudinal properties of unidirectional continuous fibrous composites. Numerical methods using finite difference, finite element, or boundary element methods yield the best predictions; however, they are time-consuming and do not yield closed-form expressions. Variational methods based on energy principles have been developed to establish bounds (inequality relations) on the effective properties. The bounds are close to each other in the case of longitudinal properties, but they can be far apart in the case of transverse and shear properties. Semiempirical relationships have been developed to avoid the difficulties with the above theoretical approaches and to facilitate computations. The so-called Halpin-Tsai relationships have consistent forms for all properties of fibrous composite materials and can be used to predict the effects of a large number of system variables. Table 1 summarizes the mathematical formulas for determining the equivalent mechanical and physical properties for known type and volume fractions of the fiber (Vf) and matrix (Vm) materials [23]. The 1 and 2 subscripts denote the principal directions of an orthotropic lamina, defined as follows: direction (1) principal fiber direction, also called fiber longitudinal direction; direction (2) in-plane direction perpendicular to fibers, transversal direction. The factor ξ is called the reinforcing efficiency and can be determined experimentally for specified types of fiber and matrix materials. Whitney [24] suggested the range 1 < ξ < 2 depending on the fiber array type, for example, hexagonal, square, etc. Usually, ξ is taken equal to 1.0 for theoretical analysis procedures in the case of carbon or glass fibrous composite laminates.

PropertyMathematical formula*
Young’s modulus in direction (1) E11Em Vm + E1f Vf
Young’s modulus in direction (2) E22Em (1 + ξηVf)/(1-ηVf); η = (E2f –Em)/(E2f + ξEm)
Shear modulus G12Gm (1 + ξηVf)/(1-ηVf); η = (G12f –Gm)/(G12f + ξGm)
Poisson’s ratioϑ12ϑmVm+ϑ12fVf
Mass density ρρm Vm + ρf Vf

Table 1.

Halpin-Tsai semiempirical relations for calculating composite properties [23].

Subscripts “m” and “f” refer to properties of matrix and fiber materials, respectively.

Assuming no voids are present, then Vm + Vf = 1.


3. Frequency optimization of FGM thin-walled box beams

This section presents a mathematical model for optimizing the dynamic performance of thin-walled FGM box beams with closed cross sections. The objective function is to maximize the natural frequencies and place them at their target values to avoid the occurrence of large amplitudes of vibration. The variables considered include fiber volume fraction, fiber orientation angle, and ply thickness distributions. Various power-law expressions describing the distribution of the fiber volume fraction have been implemented, where the power exponent was taken as the main optimization variable [25]. The mass of the beam is kept equal to that of a known reference beam. Side constraints are also imposed on the design variables in order to avoid having unacceptable optimal solutions. A case study on the optimization of a cantilevered, single-cell spar beam made of carbon/epoxy composite is considered. The results for the basic case of uncoupled bending motion are given.

3.1. Structural dynamic analysis

Figure 4 shows a slender, composite thin-walled beam constructed from uniform segments, each of which has different cross-sectional dimensions, material properties, and length. Tapered shapes of an actual blade or wing spar can be adequately approximated by such a piecewise structural model with a sufficient number of segments. The various parameters and variables are normalized with respect to a reference beam, which is constructed from just one segment with single unidirectional lamina having equal fiber and matrix volume fractions, that is, Vfo = Vmo = 50%. The different quantities are defined in the following:

Figure 4.

General configuration of a multisegment, composite box beam [25]. (a) Circular cross section, (b) rectangular cross section.

Ns = number of segments (panels).

j = subscript for the j-th segment, j = 1, 2,…….Ns.

NL(j) = number of layers in the j-th segment.

k = subscript for the k-th layer, k = 1, 2,NL(j).

L̂j=(Lj/Lo) = normalized length of the j-th segment.

L̂=L/Lo=j=1NsL̂j = normalized total beam length.

Ĥj = Hj/H0 = Ĥj=k=1NLjĥkj= normalized total wall thickness of the j-th segment.

ĥkj=hkj/H0 = normalized thickness of the k-th layer in the j-th segment.

θkj = fiber orientation angle in the k-th layer in the j-th segment.

Γ̂j=Γj/Γo=Jds/ods = normalized circumference of the j-th segment cross section.

Γj = πDj for circular C.S., Γj = 2(aj + bj) for rectangular C.S.

ρ̂kj=ρkj/ρo = normalized density of the k-th layer in the j-th segment.

Vf,kj = fiber volume fraction in the k-th layer in the j-th segment.

ρkj = ρfVf,kj+ρm1Vf,kj, ρo= 0.5(ρf + ρm).

ρf = fiber mass density, ρm = matrix density.

m̂j=mj/mo= normalized mass per unit length of the j-th segment.

mj=Γjk=1NLjρkjhkj = mass per unit length of the j-th segment, mo=ΓoρoHo.

Ij = mass polar moment of inertia per unit length of the j-th segment.

= kNLjρkjhkjy2+z2 ds.

The normalized total structural mass is given by the expression:


where Mo=moLo=ΓoρoHoLo is the total mass of the uniform baseline design. A quantity with subscript “o” refers to a reference beam parameter.

3.1.1. Constitutive relationships

The displacement field of anisotropic thin-walled closed cross-sectional beams was derived by Dancila and Armanios [26], who used a variational asymptotic approach to obtain the following constitutive equations:


where Fx, Mx, My , and Mz stand for the axial force, torsional, and bending moments, respectively, and Cmn are the cross-sectional stiffness coefficients derived in terms of closed-form integrals of the geometry and material constants. The notations U1, U2, U3, and ϕ are the kinematic variables representing the average displacements and rotation of the cross section. The primes denote differentiation with respect to x.

3.1.2. Equations of motion

The general equations of motion for the free vibration analysis are derived using Hamilton’s principle and expressed in terms of the kinematic variables, where it was shown that a closed form solution is not available [25]. However, particular choices of cross-sectional shape and layup can produce zero coupling coefficients in the equations of motion. Two special layup configurations can be considered, namely circumferentially uniform stiffness (CUS) and circumferentially asymmetric stiffness (CAS). The equations of the CUS type consist of two coupled equations for extension-twist and two uncoupled bending equations. For the CAS type, the extension displacement (U1) is uncoupled, as well as the edgewise bending (U2), while the flapping displacement (U3) is coupled with twist (φ). The general solution can be obtained by separating the space and time variables, where the time dependence is assumed to be harmonic with circular frequency, ω. The solutions for the uncoupled axial and bending equations are straightforward, while those for the coupled equations involve much mathematics [26].

3.1.3. Solution procedure of uncoupled bending motion

The basic important case to be considered first is the uncoupled bending response, which exists in both CUS and CAS layup configurations. Using the multisegment model depicted in Figure 4 and considering flapping motion (U3), the associated eigenvalue problem can be written directly in the form:


which must be satisfied over the length Lj of any segment composing the beam structure. Normalizing with respect to the reference beam, we get:


where β̂j=ω̂m̂j/Ĉ33,j1/4, Ĉ33,j=C33,j/C33,0, and ω̂=ωL02m0C33,01/2. Eq. (9) must be satisfied in the interval 0x¯L̂j, where x¯=x̂x̂j is a local coordinate of the j-th segment and x̂=x/Lo. The general solution is well known to be:


Expressing the constants ai,i=1,2,3,4 in terms of the state variables vector {S}T = {U3U3C33U3C33U3}T at both ends of the j-th segment, we get


where [T(j)] is called the transfer matrix of the j-th segment with its elements given in detail in Ref. [25]. The state variable vectors can be computed progressively along the length of the beam by applying continuity among the interconnecting joints of the different segments composing the beam structure. An overall transfer matrix denoted by [T], which relates the state variables at both ends of the beam, can be obtained from the following matrix multiplication:


The required frequency equation for determining the natural frequencies can then be obtained by applying the associated boundary conditions and considering only the nontrivial solution of the resulting matrix equation.

3.2. Formulation of the optimization problem

Several design objectives can exist in structural optimization including minimum mass, maximum natural frequencies, minimum manufacturing cost, etc. [17]. Considering the reduction of vibration level, two optimization alternatives can be formulated, namely, frequency placement by separating the natural frequencies from the harmonics of the excitations or direct maximization of the natural frequencies. The latter can ensure a simultaneous balanced improvement in both stiffness and mass distributions of the vibrating structure. The related optimization problems are usually formulated as nonlinear mathematical programming problem where iterative techniques are implemented for finding the optimal solution in the selected design space. Numerous computer programs [18] are available to solve nonlinear optimization models, which can be interacted with structural and eigenvalue analyses routines. The MATLAB toolbox optimization routines can be useful in solving some types of unconstrained and constrained optimization problems. One of the most commonly applied routines that find the constrained optima of a nonlinear merit function of many variables is named “fmincon” [19].

3.2.1. Basic optimization problem

Before performing the necessary mathematics, it is essential to recognize that design optimization is only as meaningful as its core model of structural analysis. Any deficiencies therein will absolutely be affected in the optimization process. Consider the basic problem of a uniform cantilevered, thin-walled, single-cell spar constructed from just one segment with one unidirectional lamina (Ns = 1, NL = 1). The total length and outer cross-sectional dimensions are given preassigned values equal to those of the baseline design. The remaining set of variables is, therefore, X = VfĤθ. The associated frequency equation for such a basic case is:


It is seen that ω̂ is an implicit function of the design variables and can be calculated numerically by any suitable method such as Newton-Raphson or the Bisection method. However, the frequency equation can be solved directly for the whole term ω̂m̂/Ĉ3314L̂ without regard to the specific values of the design variables. The computed roots are:


In Eq. (14), the frequency parameter ω̂i can be imagined as an explicit function of the design variables. So, for prescribed values of the design variables within the domain of side constraints, ω̂i can be obtained directly from the above equation. Therefore, it is possible to place the frequency at its desired value and obtain the corresponding value of any one of the design variables directly from Eq. (14). The selected composite material of construction is made of epoxy-3501-6 and carbon-AS4 (see Table 2), which has favorable properties and is highly recommended in many applications of civil, aerospace, and mechanical engineering [23].

PropertyFiber: Carbon-AS4Matrix: Epoxy-3501-6
Density (g/cm3)ρf = 1.81ρm = 1.27
Modulus of elasticity (GPa)E1f = 235.0
E2f = 15.0
Em = 4.30
Modulus of rigidity (GPa)G12f = 27.0Gm = 1.60
Poisson’s ratioν12f = 0.20νm = 0.35

Table 2.

Material properties of fiber and matrix materials [23].

Figure 5 depicts the functional behavior of the dimensionless fundamental frequency parameter ω̂1 combined with the structural mass constraint (M̂s=1). The imposed side constraints are:


Figure 5.

Level curves of ω̂1 function augmented with the constraint M̂s=1 in Vfθ design space (Ns = 1, NL = 1, L̂=1).

It is remarked that the function is continuous and well behaved everywhere in (Vf – θ) design space. The contours are symmetric about the horizontal line θ = 0 where the constrained global maxima occurs when the fiber volume fraction reaches its upper limiting value. It can then be concluded that the unidirectional lamina is favorable when considering beam designs with maximum bending frequency. The optimal design point was found to be (Vf, Ĥ, θ) = (0.75, 0.92, 0) at which ω1max=2.02589. This corresponds to an optimization gain of about 8.04% as measured from the reference value 1.8751. Before ending this section, it is interesting to address here the dual optimization problem of minimizing the total structural mass under preserved frequency (ω̂1=1.8751). The optimal solution was calculated to be (Vf, Ĥ, θ) = (0.50, 0.915, 0) and M̂s,min = 0.915, which corresponds to a mass saving of 8.5% as compared to the baseline design.

A couple of words are stated here regarding the side constraints in Eq. (15). First of all, it is reminded that the main focus of the present study is to optimize the fiber volume fraction in order to achieve higher values of the natural frequencies without mass penalty. The optimization is performed with respect to a known baseline design, which is considered to be conservative having reserve strength to withstand severe dynamic loads. The imposed side constraint on the total wall thickness, normalized with respect to that of the baseline design, is included for consideration of strength and stability requirements, which are not considered in the present study. So, the imposed limits with a percentage of 25% below or above that of the baseline can be practically accepted for the given model formulation. On the other hand, appropriate values of the upper and lower bounds imposed on the fiber volume fraction are chosen to avoid having unacceptable designs from the manufacture point of view. For example, the filament winding is usually associated with the highest fiber volume fractions. With careful control of fiber tension and resin content, values of around 75% would be reasonable [27].

3.2.2. Optimization model for discrete grading

A comprehensive analysis and formulation of discrete optimization models for beam structures considering both stability and dynamic performance were formulated in [28], where mathematical programming coupled with finite element analysis procedures was implemented. For the case of a two-segment spar beam, (Ns = 2, NL = 1), the reduced optimization problem can be defined as follows:

MinimizeFX¯=ω̂1Subject toM̂s=11NsL̂j=10.250.0VfjL̂jj=1,20.751.0E16

Using the equality constraints, two of the design variables can be expressed in terms of the other two variables. Figure 6 shows the functional behavior of the dimensionless frequency combined with the structural mass constraint. It is remarked that the function is well defined in the feasible domain of the selected design space (VAL̂)1. Two empty regions can be observed at the upper left and right parts of the design space, where violation of the equality mass constraint is indicated. In the left one, the fiber volume fraction is equal to 100%, violating the imposed side constraint. The feasible domain is seen to be split into two distinct zones separated by the baseline contour, which is represented by the vertical line Vf1 = 50%. The constrained optimum is to found to be (Vfj,L̂j)j=1,2 = (0.75, 0.50), (0.25, 0.50) corresponding to (ω̂1)max = 2.0645 with 10.10% optimization gain.

Figure 6.

Level curves of ω̂1 function augmented with M̂s=1 in Vf1L̂1 design space (Ns = 2).

3.2.3. Optimization model for continuous grading

For continuous grading models, the associated optimization problem is cast as follows: find the design variables vector X = (Δf, p), which minimizes the objective function:


subject to the constraints:


Solutions obtained by applying the power-law model of Eq. (3) have shown that no improvements can be achieved using grading of the fiber volume fraction in the thickness direction. On the other hand, grading in spanwise direction has shown some interesting results. Considering spanwise grading according to Eq. (4), Figure 7 depicts the level curves of the fundamental frequency parameter ω̂1 combined with the mass constraint in the design space fp. It is observed that the feasible domain is bounded from below and above by the constraint curves corresponding to the upper and lower bounds imposed on the fiber volume fractions at tip and root. The horizontal line f = 1.0 (i.e., Vf = 50% at root and tip) split the domain into two zones. The lower zone encompasses the constrained optimum solution: (ω̂1)max = 2.01875 at the design point (f, P)opt. = (0.34, 1.01).

Figure 7.

Level curves of ω̂1 function augmented with M̂s=1 in fp design space (Ns = NL = 1) with spanwise grading Eq. 4.

Table 3 summarizes the attained optimal solutions for the different grading patterns. It is seen that the highest optimization gain is obtained by using spanwise grading of Eq. (5) with the coordinate exponent n = 3.

Vf—power-law model(Δf, p)opt.,ω1max,gain %
Thickness grading (Eq. (3))(1.0, 0.0),1.87510.0%
Spanwise grading (Eq. (4))(0.34, 1.01),2.01875,7.66%
Spanwise grading (Eq. 5)
n = 1(0.34, 1.02),2.01938,7.70%
n = 2(0.34, 2.425),2.04813,9.23%
n = 3(0.34, 5.175),2.06125,9.93%

Table 3.

Optimal solutions using different grading patterns (M̂s = 1).


4. Optimization of FGM drive shafts against torsional buckling and whirling

One of the important design issues in mechanical industries is the buckling and whirling instabilities that may arise from the loads applied to a power transmission shaft. These instabilities result in a reduced control of the vehicle, undesirable performance, and often cause damage, sometimes catastrophic, to the vehicle structure. Therefore, by incorporating such considerations into an early design optimization [29], the design space of a power transmission shaft will be reduced such that undesirable instability effects can be avoided during the range of the vehicle’s mission profile. Figure 8 shows an idealized structural model of a long, slender composite shaft having circular thin-walled cross section. The main structure is constructed solely of functionally graded, fibrous composite materials. The laminate coordinates are defined by x parallel to the shaft axis, y points to the tangential direction, and z points to the radial direction. Predictions of both torsional buckling and whirling instabilities are based on simplified analytical solutions of equivalent beam and shell structures. The coupling between bending and torsional deformations, introduced by the composite construction, and its influence on such instabilities is considered.

Figure 8.

Shaft model and definition of reference axes.

4.1. Torsional buckling optimization problem

Bert and Kim [30] derived the governing differential equations of torsional buckling in the form:


where Nx and Ny are the normal forces, Nxy and Nyx are shear forces, Mx and My are bending moments, and Mxy and Myx are torsional moments. All are applied to the midsurface and measured per unit wall thickness of the shaft. T is the applied torque, R is the mean radius, and (u, v, w) are the displacements of a generic point on the middle surface of the shaft wall. An iterative process is outlined in Ref. [30] for calculating the buckling torque for specified boundary conditions. There are other simple empirical equations based on experimental studies that can give a reasonable estimate of the buckling torque. The most commonly used formula for the case of simply supported shaft is [31]:


where Tcr is the critical buckling torque and H is the total wall thickness of the shaft. Expressions of the equivalent modulii of elasticity in the axial (Ex) and hoop (Ey) directions for symmetric and balanced laminates are given in Ref. [31]. The various parameters and variables are normalized with respect to known baseline design, which is constructed from cross-ply laminates [0o/900]N with equal volume fractions of the fibers and matrix materials, that is, Vf = Vm = 50%. Optimized shaft designs shall have the same transmitted power, length, outer diameter, boundary conditions, and material properties of those known for the baseline design. The different dimensionless quantities are defined in Ref. [31]. The optimal torsional buckling problem is to find the design variables vector X=Vfθĥk=1,2,NL, which minimizes the objective function:

MinimizeF=T̂crsubject tomass limitation:M̂10E20
Torsional strength:τmaxτallow1.00Whirling:Ω̂maxΩ̂cr0E21
Side constraints:0.30π20.015Vfθĥk=1,2,NL0.70π20.200.75k=1NLĥk1.25E22

where T̂cr=Tcr/Tcro, Ω̂=Ω2π/60ω1,o are the dimensionless critical torque and rotational speed, respectively. The baseline design parameters are denoted by subscript “o.” τmax=Tmax/2πR2H is the maximum shear stress, Tmax is the maximum applied torque, and τallow is the allowable shear stress that can be calculated according to the embedded material properties and volume fraction of the fiber [23]. This optimization problem may be thought as a search in an (3NL) dimensional space for a point corresponding to the minimum value of the objective function and such that it lies within the region bounded by subspaces representing the constraint functions. It must be noted that the outside dimensions (outer diameter and length) of the shaft are restricted by the available interior space of the vehicle and will be considered as preassigned parameters in the present model formulation. The first case study to be examined herein is a shaft with discrete thickness grading constructed from eight plies (±θ ±θ)s with the same properties of carbon/epoxy composites (see Table 2) and same thicknesses. This sequence is applied in filament wound circular shells, as such a process demands adjacent (±θ) layers. Figure 9 depicts the obtained contours in the (Vf1-θ) design space, which are, as seen, monotonic and symmetric about the zero ply angle. A local maximum of T̂cr can be observed near the design point (Vf1,θ) = (0.7, 90o) with T̂cr = 1. This figure illustrates that the maximum critical buckling torque can be achieved when the fiber orientation angle is close to 90o. Other case studies including both discrete and continuous grading with several optimal solutions can be found in Ref. [31]. At the start of the optimization process in each case, the shaft wall was divided into a large number of layers with equal thicknesses, for example, NL = 32. It has been found that the optimization algorithm treats the number of layers as an additional implicit variable. Sometimes the computer discards one or more layers by letting their thicknesses sink to the lower limits and sometimes makes some consecutive layers identical, that is, having the same fiber orientation and volume fraction. Such a situation was repeated for many cases of study. It was found that the appropriate number to be taken for the shaft problem under consideration is NL = 8. This would eliminate much of the numerical effort necessary for performing structural analysis in each optimization cycle and, consequently, reduces the computational time considerably.

Figure 9.

T̂cr—contours in (Vf1-θ) design space under mass constraint M̂=1. (Case of drive shaft with eight symmetric, balanced, carbon/ epoxy layers)

The final attained optimal solution was a cross ply layup [900/00]4 with the fiber volume fraction in the eight layers reached its upper value of 70%. The optimal dimensionless ply thickness was found to be [0.1994, 0.0967, 0.152, 0.019]s at which the shaft torsional buckling capacity was increased by 32.1% above that of the baseline design. However, the total structural mass has reached its baseline value, and whirling constraint became active at the achieved optimum design point.

4.2. Whirling optimization problem

The calculation of the critical speed, also referred to as whirl instability of a rotating shaft, is based on the work given in Refs. [32, 33]. The critical speed is defined as the point at which the spinning shaft reaches its first natural frequency. The shaft is modeled as a Timoshenko beam, which implies that first-order shear deformation theory with rotatory inertia and gyroscopic action was used. The shaft is assumed to be pinned at both ends with a Cartesian coordinate system (x, y, z), where x is measured along the longitudinal axis of the shaft. The displacements in the y and z directions are denoted by v and w, respectively, and Φ is the angle of twist. The cross-sectional area, second moment of area, and polar moment of area are denoted by A, I, and J, respectively. The equations of motion are derived by invoking Hamilton’s principle, with the following results [32]:


The symbol t denotes time and Ω the rotational speed. The ρA, ρI, and ρJ terms account for translational, rotary, and torsional inertias, respectively, while the 2ρI terms account for the gyroscopic inertia effects. It is assumed that the flexural and bending-twisting coupling rigidities (CB and CBT) associated with bending about the y and z axes are identical; likewise for the transverse shear stiffness (Cs) [32]. Bert and Kim [33] considered the case of simply supported shaft and assumed separable solution in space and time to solve the associated eigenvalue problem. The derived frequency equation is given by:


where λ=/L, ω = circular natural frequency, and n = mode number. For each natural frequency of the nonrotating shaft, the rotational speed (Ω) develops gyroscopic moments, which cause the natural frequency to bifurcate into two. The higher of the two increases with Ω and is associated with forward precession, while the lower one decreases with Ω and is associated with backward precession. A critical instability occurs when the rotational speed coincides with the first backward-precision natural frequency, which is termed as the first critical speed. Two alternatives may be considered regarding the whirling optimization problem [34]:

(a) Direct maximization of the critical rotational speed

Find the design variables vector X=Vfθĥk=1,2,NL, which minimizes the objective function:

MinimizeF=Ω̂crSubject toM̂10E28

(b) Placement of the critical speed

The other alternative of the objective function is defined by:


The same set of constraints given in Eq. (23) is applied. The notation Ω̂ is a dimensionless target rotational speed, which should be greater than the maximum permissible rotational speed by a reasonable margin (e.g., 10–20%). As a case study, a drive shaft with continuous material grading along the shaft axis is optimized considering the following power-law model:


where Vf (0) is the fiber volume fraction at the right or left end of the drive shaft, while Vf (0.5) is the fiber volume fraction at the middle of the shaft length. Figure 10 illustrates the level curves of the normalized critical speed augmented with the mass equality constraint. It is seen that there are four distinct zones separated by the contour lines Ω̂cr=1.0. The upper left zone and lower right zone contain local maximum solutions. The best point (p,f) = (4.53, 0.3), corresponding to Ω̂cr=1.045, is located inside the zone where the fiber taper ratio f=Vf0/Vf0.5 is less than one. The upper empty zone contains infeasible solutions that violate the imposed constraints. Another case study considers the through-thickness grading pattern given by Eq. (3). The corresponding design variable vector is defined by X=Vf0Vf12pĤ with lower and upper limits XL = (0.3, 0.3, 0, 0.75) and XU = (0.7, 0.7, ∞, 1.25). The attained optimal design variable vector was calculated to be Xopt= at which the maximum critical speed increased by 14% above that of the baseline design with active mass constraint.

Figure 10.

Normalized critical speed Ω̂cr augmented with the mass constraint (M̂=1.0) in (p-f) design space.

A last optimization strategy to be addressed here is to combine the two criteria in a single objective function subject to the mass, strength, and side constraints.

MinimizeF=Ω̂cr+T̂crSubject toM̂10τmaxτallow1.00E32

Eq. (22) assumes that whirling and torsional buckling instabilities are of equal relative importance. This model resulted in a balanced improvement in both stabilities with active mass constraint. The attained optimal solution was found to have a uniform distribution of the fiber volume fraction with its upper limiting value of 70% and wall thickness = 0.935. The corresponding optimal values of the design objectives were Ω̂cr=1.135 and T̂cr=1.161, representing optimization gains 13.5 and 16.1%, respectively, as measured from the baseline design.


5. Optimization of FGM wings against divergence

The use of the in-plane grading in aeroelastic design was first exploited by Librescu and Maalawi [6], who introduced the underlying concepts of using material grading in optimizing subsonic rectangular wings against torsional instability. Exact mathematical models were developed allowing the material physical and mechanical properties to change in the wing spanwise direction, where both continuous and piecewise structural models were successfully implemented. In this section, analytical solutions are developed for slender tapered composite wings through optimal grading of the material volume fraction in the spanwise direction. The enhancement of the wing torsional stability is measured by maximization of the critical flight speed at which aeroelastic divergence occurs. The total structural mass is maintained at a value equals to that of a known baseline design in order not to violate other performance requirements. Figure 11 depicts a slender wing constructed from Np panels with trapezoidal planform and known airfoil cross section. The wing is considered to be made of unidirectional fiber-reinforced composites with variable fiber volume fraction in the spanwise direction. The flow is taken to be steady and incompressible, and the aspect ratio is assumed to be sufficiently large so that the classical engineering theory of torsion can be applicable and the state of deformation described in terms of one space coordinate.

Figure 11.

Trapezoidal wing planform and cross section geometry. (a) Multipanel, piecewise wing model, (b) airfoil section and applied airloads.

The chord distribution is assumed to have the form:


The symbol Δc denotes the chord taper ratio (= tip chord Ct/root chord Cr) and x (= x1/L) denotes the dimensionless spanwise coordinate. The equivalent shear modulus G of a unidirectional reinforced composite, thin-walled cross section can be determined from the relation [35]:


where f1 is a function that depends on the geometry and thickness ratio of the cross section (h/C) and the ratio (G12/ G13), where G12 and G13 are the in-plane and out-of-plane shear moduli, respectively (refer to Table 1). C is the chord length and h is the maximum thickness of the cross section. For many types of fibrous composites that are commonly utilized in aerospace industry [23], such as carbon/epoxy and graphite epoxy, both moduli are approximately equal, G12 ≈ G13.

Using the classical elasticity and aerodynamic strip theories, the governing differential equation of torsional stability in dimensionless form is [35]:


?The associated boundary conditions of the elastic angle of attack, α, are α(0) = 0 and α1=0. The symbol Ĝ=G12/G12,o denotes the dimensionless shear modulus, Ĵ=J/Jr denotes the dimensionless torsion constant, and the prime denotes differentiation with respect to the dimensionless coordinate x = x1/L. The dimensionless flight speed is defined by V̂=VCrbρae/2GJr, where (GJ)r is the torsional stiffness of the baseline design at root. The shear modulus G12,0 of the baseline design can be calculated by taking Vfo = 50%.

Considering the K-th panel of the wing as shown in Figure 11a, and using the transformation y = (1-βx), Eq. (25) takes the form:


where ak=V̂/βĥkĜk, ĥk and Ĝk are the normalized wall thickness and shear modulus of the kth wing panel, respectively. The general solution of Eq. (26) is:


where J2 and Y2 are Bessel’s function of the first and second kind with order 2, respectively [35], and A1 and A2 are the constants of integration. The dimensionless internal torsional moment, T, can be obtained by differentiating Eq. (27) and multiplying by the dimensionless shear rigidity. Applying the boundary conditions at stations (k) and (k + 1), the constants A1 and A2 can be expressed in terms of the state variables at station (k), which can be related to those at station (k + 1) by the transfer matrix relation:


It is now possible to compute the state variables progressively along the wing span by applying continuity requirements of the variables (α, T) among the interconnecting boundaries of the various wing panels. The divergence speed can be calculated by applying the boundary conditions and considering the nontrivial solution of the resulting equations (similar to the procedure outlined in Section 3.1.3). The associated optimization problem may be cast in the following:

MinimizeV̂divSubject toM̂s1.ĥb̂k=1,2,..Np0.751.251.0b̂k=1.0E39

The preassigned parameters that do not change during the optimization process include the wing semispan (b), the chord taper ration (c), airfoil type and geometry, and fiber and resin material types. This model has been applied to obtain wing designs with improved torsional stability by maximizing the divergence speed (Vdiv) without weight penalty. The selected material is carbon-AS4/epoxy-3501-6 composite (see Table 2), which has favorable characteristics and is highly desirable in manufacturing aircraft structures. The baseline design has uniform mass and stiffness distributions and is made of uniform unidirectional fibrous composite with equal volume fraction of the matrix and fiber materials, that is, Vfo = 50%. Figures 12 and 13 show the developed level curves of constant divergence speed (also named isodiverts) for two-panel wings with chord tapering ratio, c = 0.5. Actually, these curves represent the dimensionless critical speed, augmented with the equality mass constraint. Examining Figure 12, it is seen that the Vdiv function is well behaved and continuous everywhere in the selected design space except in the empty regions to the upper left and right regions, where the equality mass constraint is violated. The feasible domain is bounded from above by the two curved lines representing the upper and lower limiting constraints imposed on the volume fraction of the outboard blade panel. The contours inside the feasible domain are not allowed to penetrate these borderlines and obliged to turn sharply to be asymptotes to them, in order not to violate the mass constraint. The final attained optimal solutions are summarized in Table 4. It can be observed that good wing patterns shall have the lower limit of the fiber volume fraction at the tip and the upper limit at root. Using material and wall thickness grading together results in a considerable enhancement of the wing torsional stability.

Figure 12.

Isodivert in (Vf1- b1) design space for a two-panel wing model (h1 = h2 = 1.0, Δc = 0.5,M̂s=1).

Figure 13.

Isodivert in (Vf1- b1) design space for a two-panel wing model (h2 = 0.5, Vf2 = 0.3, Δc = 0.5,M̂s=1).

NpType of gradingXopt=Vfĥb̂k=1,2,..NpV̂div,maxOptimization gain %
2Material(0.5906, 1.0, 0.6594)1
(0.2500, 1.0, 0.3406)2
Material and thickness(0.75, 1.0625, 0.7)1
(0.25, 0.5250, 0.3)2
3Material(0.6000, 1.0, 0.5750)1
(0.4125, 1.0, 0.1375)2
(0.2500, 1.0, 0.2875)3
Material and thickness(0.75, 1.10, 0.65)1
(0.30, 0.65, 0.05)2
(0.25, 0.50, 0.3)3

Table 4.

Optimal piecewise wing designs with chord taper = 0.5 (baseline value = 1.9747).


6. Optimization of composite thin-walled pipes conveying fluid

The subject of vibration and stability of thin pipes conveying flowing fluids is of a considerable practical interest. An advanced textbook by Paїdoussis [36] gives an excellent review of the several developments made in this research area. Practical models for enhancing static and dynamic stability characteristics of pipelines constructed from uniform modules were addressed by Maalawi et al. [37, 38], where the relevant design variables were selected to be the mean diameter, wall thickness, and length of each module composing the pipeline. The general case of an elastically supported pipe, covering a variety of boundary conditions, was also investigated. Distinct domains of the flutter instability boundaries were presented for different ratios of the fluid-to-pipe mass, and the variation of the critical flow velocity with support flexibility was examined and discussed. Concerning pipelines made of advanced FGMs, this section presents a mathematical model for enhancing the overall system stability against flutter and/or divergence under mass constraint. Figure 14 shows a FGM pipe conveying flowing fluid with the coordinate system chosen such that the x-axis coincides with the longitudinal centroidal axis in its undeformed position, while the y- and z-axes coincide with the cross section principal axes. The pipe model consists of rigidly connected thin-walled circular tubes made of unidirectional fibrous composite material. Each pipe module has different material properties, wall thickness, and length. Such a configuration results in a piecewise axial grading of either the material of construction or the wall thickness in the direction of the pipe axis. Assuming no voids are present, the distributions of the mass density (ρ) and modulus of elasticity (E) can be determined using the formulas of Table 1.

Figure 14.

Multimodule composite pipe conveying flowing fluid.

The associated eigenvalue problem is described by the fourth-order ordinary differential equation [39]:


where V(x) is the dimensionless mode shape satisfying boundary conditions, and ω is the corresponding dimensionless frequency of oscillation, which will be, in general, a complex number to be determined by the requirement of nontrivial solutions, V(x) ≠ 0. More details for the definition of the various parameters and dimensionless design variables are given in Ref. [39]. Both static and dynamic instability phenomena can be involved for the physical model described by Eq. (30), depending on the type of boundary conditions at the pipe ends. The notation ()′ means total spatial derivative. The general solution can be obtained using standard power series methods:


pj, j = 1,…4 are the four roots of the fourth-order polynomial:


An efficient method [38] is successfully implemented to find the complex roots of Eq. (32) by formulating a special companion matrix and finding the associated eigenvalues for any desired values of the variables αk, βo, ω, and U. The transfer matrix [Tk] of the kth pipe module can be obtained by performing the matrix multiplications


and [E] is a diagonal matrix with elements:


Applying the appropriate boundary conditions and considering only the nontrivial solution, the resulting characteristic equation can be solved numerically for the frequency and the critical flow velocity. The system is stable or unstable according to whether the imaginary component of the frequency, ω, is positive or negative, respectively. In case of neutral stability, ω is wholly real. As the flow velocity increases, the system may become unstable in one or more of its normal modes. The critical flow velocity is the greatest velocity for which the system is stable in all its modes. The characteristic equation for cantilevered pipe is:


where Tij are the elements of the overall transfer matrix. For a specific mode number, the proper starting frequency is determined by solving a subsidiary eigenvalue problem corresponding to the case of stationary fluid inside the vibrating pipe. A globally convergent optimization algorithm, known as Levenberg-Marquardt algorithm [38], has been applied to solve the resulting nonlinear equations derived from the consideration of nontrivial solution of characteristic equation. The effect of flow for small velocity is to damp the system in all modes. At higher velocities, some of the modes become less damped and the corresponding branches cross the Re(ω)-axis, indicating the existence of unstable oscillations of the system. If a branch passes through the origin, that i, ω = 0, the case of static instability (called divergence) is reached.

6.1. Flutter solutions

To verify the developed formulation, the classical problem dealing with one-module cantilevered pipe is considered first. The dimensionless wall thickness and length of the pipe are assigned at a value of 1.0, while the volume fraction at 50%. Figure 15 depicts variation of the critical flutter velocity and frequency with the mass density ratio MRo, covering a wide range of pipe and fluid mass densities. It is seen that there are four subdomains with the associated flutter modes defined in the specified intervals of the mass ratio. The upper and lower bounds determine the critical values of the mass ratios at which some of the frequency branches cross each other at the same value of the flutter velocity. The overall flutter mode may be regarded as composed of different quasimodes separated at the shown “jumps” in the Uf-MRo curve. The calculated mass density ratios at the three indicated frequency jumps are 0.4225, 2.29, and 12.33, respectively. They correspond to multiple points of neutral stability, where for a finite incremental increase in the flow velocity, the system becomes unstable, then regains stability, and then once again becomes unstable with a noticeable abrupt increase in the flutter frequency. Next, we consider a baseline design made of carbon/epoxy composites (see Table 2) with mass ratio MRo = 2.0. The calculated values of the dimensionless flutter velocity and frequency are found to be Uf = 10.78 and ωf = 26, respectively. Keeping the total dimensionless mass constant at the value corresponding to the baseline design, the best solution having the greatest flutter velocity was found to be (Vf, h1) = (0.70, 0.9345), which corresponds to Uf = 12.517 and ωf = 31.8615.

Figure 15.

Variation of flutter speed and frequency with mass ratio for a uniform one-module cantilevered pipe (Vf1 = 50%, h1 = 1, L1 = 1).

6.1.1. Solutions for cantilevered two-module pipe with uniform thickness (h1 = h2 = 1.0)

Considering the case of two-module pipe, a direct and fast way for checking out system stability for any desired set of the dimensionless design variables (Vf1, L)k = 1,2 is given here. Lower and upper bounds are imposed on the design variables in order not to violate other strength and manufacturing requirements. The fiber volume fraction is constrained to be within the range 30% up to 70%, while the dimensionless length is between 0.0 and 1.0. The mass ratio MRo is taken to be 2.0.

Dimensionless flutter velocity and flutter frequency are obtained from the frequency and velocity branches at the four modes. The lowest frequency and velocity among the four modes at which Imag(ω) = 0.0 are considered the flutter velocity and frequency. These computed values at different conditions are employed in constructing the flutter velocity and frequency contours as shown in Figure 16. The white regions shown in both figures indicate that the fiber volume fraction of the material of the second segment does not fall in range between 0.3 and 0.7. The maximum flutter velocity (Uf) and its corresponding flutter frequency (ωf) occur in the region colored with dark brown L1 = 0.36 and Vf1 = 0.3. The maximum Uf and its corresponding ωf are 13.67 and 58.4, respectively. Table 5 gives several standard optimal solutions for the two-module case study. The global optimum design point is seen to be (Vf, L)k = 1,2 = (0.35, 0.40), (0.60, 0.60), at which the normalized flutter velocity reached a value of 13.31 corresponding to 23.47% optimization gain.

Figure 16.

Contour plots of flutter velocity and frequency in (Vf1-L1) design space.

Design variables: (Vf, L)k = 1,2Dimensionless UflutterStability improvement %
(0.30, 0.20), (0.55, 0.80)
(0.30, 0.50), (0.70, 0.50)
(0.35, 0.25), (0.55, 0.75)
(0.35, 0.40), (0.60, 0.60)
23.47 (max.)
(0.40, 0.50), (0.60, 0.50)
(0.40, 0.60), (0.65, 0.40)
(0.45, 0.50), (0.55, 0.50)
(0.45, 0.75), (0.65, 0.25)
(0.50, 0.50), (0.50, 0.50)10.780.00 (baseline)

Table 5.

Standard solutions for a cantilevered two-module pipeline with uniform thickness (h1 = h2 = 1.0, M = 1.0).


7. Conclusions

As a major concern in producing efficient structures with enhanced properties and tailored response, this chapter presents appropriate design optimization models for improving performance and operational efficiency of different types of composite structural members. The concept of material grading has been successfully applied by incorporating the distribution of the volume fractions of the composite material constituents in the mathematical model formulation. Various scenarios in modeling the spatial variation of material properties of functionally graded structures are addressed. The associated optimization strategies include frequency maximization of thin-walled composite beams, optimization of drive shafts against torsional buckling and whirling instability, and maximization of the critical flight speed of subsonic aircraft wings. Design variables encompass the distribution of volume fraction, ply angle, and wall thickness as well. Detailed optimization models have been formulated and presented for improving the dynamic performance and increasing the overall stiffness-to-mass level of thin-walled composite beams. The objective functions have been measured by maximizing the natural frequencies and place them far away from the excitation frequencies, while maintaining the total structural mass at a constant value. For discrete models, the optimized beams can be constructed from any arbitrary number of uniform segments where the length of each segment has shown to be an important variable in the optimization process. It has also been proved that expressing all parameters in dimensionless forms results in naturally scaled design variables, constraints, and objective functions, which are favored by a variety of optimization algorithms. The attained optimal solutions using continuous grading depend entirely upon the prescribed power-law expression, which represents additional constraint on the optimization problem. Results show that material grading in the spanwise direction is much more better than grading through the wall thickness of the cross section. Regarding optimization of FGM drive shafts, it was shown that the best model is to combine torsional buckling and whirling in a single objective function subject to mass constraint. This has produced a balanced improvement in both stabilities with active mass constraint at the attained optimal design point.

In the context of aeroelastic stability of aircraft structures, an analytical model has been formulated to optimize subsonic trapezoidal wings against divergence. It was shown that by using material and thickness grading simultaneously, the aeroelastic stability boundary can be broaden by more than 50% above that of a known baseline design having the same total structural mass. Other stability problems concerning fluid-structure interaction have also been addressed. Both flutter and divergence optimization have been considered, and several design charts that are useful for direct determination of the optimal values of the design variables are given. It has been confirmed that the segment length is the most significant design variable in the whole optimization process. Some investigators who apply finite elements have not recognized that the length of each element can be taken as a main design variable in the whole set of optimization variables. The results from the present approach reveal that piecewise grading of the material can be promising in producing truly efficient structural designs with enhanced stability, dynamic, and aeroelastic performance. It is the author’s wish that the results presented in this chapter will be compared and validated through other optimization techniques such as genetic algorithms or any appropriate global optimization algorithm.

Actually, the most economic structural design that will perform its intended function with adequate safety and durability requires much more than the procedures that have been described in this chapter. Further optimization studies must depend on a more accurate analysis of constructional cost. This combined with probability studies of load applications and materials variations should contribute to further efficiency achievement. Much improved and economical designs for the main structural components may be obtained by considering multidisciplinary design optimization, which allows designers to incorporate all relevant design objectives simultaneously. Finally, it is important to mention that, while FGM may serve as an excellent optimization and material tailoring tool, the ability to incorporate optimization techniques and solutions in practical design depend on the capacity to manufacture these materials to required specifications. Conventional techniques are often incapable of adequately addressing this issue. In conclusion, FGMs represent a rapidly developing area of science and engineering with numerous practical applications. The research needs in this area are uniquely numerous and diverse, but FGMs promise significant potential benefits that fully justify the necessary effort.


  1. 1. Suresh S, Mortensen A. Fundamentals of Functionally Graded Materials. London: IOM Communications Ltd; 1998. 165 p. ISBN: 9781861250636
  2. 2. Miyamoto Y, Kaysser WA, Rabin BH, Kawasaki A, Ford RG. Functionally Graded Materials: Design, Processing and Applications. Springer Science & Business Media; 2013. 330 p. ISBN: 9781461553014
  3. 3. Lipton R. Design of functionally graded composite structures in the presence of stress constraints. International Journal of Solids and Structures. 2002;39:2575-2586
  4. 4. Cho JR, Shin SW. Material composition optimization for heat-resisting FGM by artificial neural network. Composites Part A Applied Science and Manufacturing. 2004;35:585-594
  5. 5. Goupee AJ, Vel SS. Optimization of natural frequencies of bidirectional functionally graded beams. Journal of Structural and Multidisciplinary Optimization. 2006;32(6):73-484
  6. 6. Librescu L, Maalawi K. Material grading for improved aeroelastic stability of composite wings. Journal of Mechanics of Materials and Structures. 2007;2(7):1381-1394
  7. 7. Birman V, Byrd WL. Modeling and analysis of functionally graded materials and structures. Applied Mechanics Reviews. 2007;60(5):195-216
  8. 8. Elishakoff I, Guede Z. Analytical polynomial solutions for vibrating axially graded beams. Journal of Mechanics and Advanced Materials and Structures. 2004;11:517-533
  9. 9. Qian L, Batra R. Design of bidirectional functionally graded plate for optimal natural frequencies. Journal of Sound and Vibration. 2005;280:415-424
  10. 10. Maalawi KY. Functionally graded bars with enhanced dynamic performance. Journal of Mechanics of Materials and Structures. 2011;6(1-4):377-393
  11. 11. Elishakoff I, Endres J. Extension of Euler’s problem to axially graded columns. Journal of Intelligent Material Systems and Structures. 2005;16(1):77-83
  12. 12. Li S-R, Batra R. Buckling of axially compressed thin cylindrical shells with functionally graded middle layer. Journal of Thin-Walled Structures. 2006;44:1039-1047
  13. 13. Maalawi KY. Optimization of elastic columns using axial grading concept. Engineering Structures. 2009;31(12):2922-2929
  14. 14. Maalawi KY. Use of material grading for enhanced buckling design of thin-walled composite rings/long cylinders under external pressure. Composite Structures. 2011;93(2):351-359
  15. 15. Shin-Yao K. Flutter of rectangular composite plates with variable fiber spacing. Composite Structures. 2011;93:2533-2540
  16. 16. Maalawi KY. Stability dynamic and aeroelastic optimization of functionally graded composite structures. In: Coskun SB, editor. Advances in Computational Stability Analysis. Rijeka: InTech; 2012. pp. 17-42. DOI: 10.5772/45878.ch2
  17. 17. Rao S. Engineering Optimization: Theory and Practice. 4th ed. New York: John Wiley & Sons; 2009. ISBN: 978-0470183526
  18. 18. 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
  19. 19. Venkataraman P. Applied Optimization with MATLAB Programming. 2nd ed. New York: John Wiley & Sons; 2009. ISBN: 978-0470084885
  20. 20. Chen WH, Gibson RF. Property distribution determination of non-uniform composite beams from vibration response measurements and Galerkin’s method. Journal of Applied Mechanics. 1998;65:127-133
  21. 21. Chi SH, Chung YL. Mechanical behavior of functionally graded material plates under transverse load, I: Analysis. International Journal of Solids and Structures. 2006;43:3657-3674
  22. 22. Bedjilili Y, Tounsi A, Berrabah HM, Mechab I. Natural frequencies of composite beams with a variable fiber volume fraction including rotary inertia and shear deformation. Applied Mathematics and Mechanics. 2009;30(6):717-726
  23. 23. Daniel IM, Ishai O. Engineering Mechanics of Composite Materials. 2nd ed. New York: Oxford University Press; 2006
  24. 24. Whitney JM. Elastic moduli of unidirectional composites with anisotropic filaments. Journal of Composite Materials. 1967;1:188-194
  25. 25. Maalawi KY. Dynamic optimization of functionally graded thin-walled box beams. International Journal of Structural Stability and Dynamics. 2017;17(9):1-24
  26. 26. 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(23):3105-3119
  27. 27. Kieback B, Neubrand A, Riedal H. Processing techniques for functionally graded materials. Materials Science and Engineering. 2003;362:81-105
  28. 28. Maalawi KY, El-Chazly NM. Global optimization of multi-element beam-type structures. In: Proceedings of the Second International Conference on Advances in Structural Engineering and Mechanics (ASEM'02). Busan, South Korea; August 2002. pp. 21-23
  29. 29. Bijagare Anup A, Mehar PG, Mujbaile VN. Design optimization & analysis of drive shaft. VSRD International Journal of Mechanical, Automobile & Production Engineering. 2012;2(6):210-215
  30. 30. Bert CW, Chun-Do K. Analysis of buckling of hollow laminated composite drive shafts. Composites Science and Technology. 1995;53:343-351
  31. 31. Nasr MF, El-Zoghby AA, Maalawi KY, Azzam BS, Badr MA. Torsional buckling optimization of composite drive shafts. World Applied Sciences Journal. 2015;33(3):517-524
  32. 32. Bauchau OA. Optimal design of high speed rotating graphite/epoxy shafts. Journal of Composite Materials. 1983;17(3):170-180
  33. 33. Bert CW, Chun-Do K. Whirling of composite-material drive shafts including bending-twisting coupling and transverse shear deformation. Journal of Vibration and Acoustics, ASME. 1995;117:17-21
  34. 34. Nasr MF, Maalawi KY, Badr MA. Optimization of composite drive shafts against whirling instability. In: Proceedings of the 8th International Conference on Mechanical and Aerospace Engineering (ICMAE); 22-25 July 2017; Prague. Czech Republic: IEEE; 2017. pp. 816-821. DOI: 10.1109/ICMAE.2017.8038756
  35. 35. Librescu L, Maalawi KY. Aeroelastic design optimization of thin-walled subsonic wings against divergence. Thin-Walled Structures. 2009;47:89-97
  36. 36. Paїdoussis MP. Fluid Structure Interactions: Volume I; Slender Structures and Axial Flow. 2nd ed. New York: Academic Press; 2013. 880 p. ISBN: 9780123973122
  37. 37. Maalawi KY, Ziada MA. On the static instability of flexible pipes conveying fluid. Journal of Fluids and Structures. 2002;16(5):685-690
  38. 38. Maalawi KY. Abouel Fotouh AM: Investigation of instability of flexible pipes transporting fluids. Journal of Engineering and Applied Science. 2005;52(6):1199-1217
  39. 39. Maalawi KY, Abouel-Fotouh AM, El Bayoumi M, Ali Yehia KA: Design of composite pipes conveying fluid for improved stability characteristics. International Journal of Applied Engineering Research. 2016;11(12):7633-7639. ISSN 0973-4562

Written By

Karam Maalawi

Submitted: May 28th, 2018 Reviewed: November 6th, 2018 Published: November 26th, 2018