Optimization of Functionally Graded Material Structures: Some Case Studies

Abstract

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.

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):

$$v\left(\mathrm{z}\right)={\left(\frac{1}{2}+\frac{z}{h}\right)}^p-1/2\le \left(z/h\right)\le \mathit{1}/\mathit{2},p\ge \mathit{0}$$E1

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 $\widehat{z}=\left(\mathit{zlh}\right)=\pm 1/2$ for consideration of other strength requirements and/or manufacturing restrictions. The modified form can be expressed as follows [5]:

$$v_f\left(\widehat{z}\right)=v_f\left(-0.5\right)+\left[v_f\left(0.5\right)-v_f\left(-0.5\right)\right]{\left(\widehat{z}+0.5\right)}^P,$$E2

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:

$$\begin{array}{c}{v}_f\left(\widehat{z}\right)=v_f\left(0\right)+\left[v_f\left(0.5\right)-v_f\left(0\right)\right]{\left(2\left|\widehat{z}\right|\right)}^P,\\ -0.5\le \widehat{z}\left(=\frac{z}{h}\right)\le 0.5,p\ge 0\end{array}$$E3

Figure 2.

Thickness distribution of the fiber volume fraction in FGM beam, v_f (0) = 40%, v_f (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:

$$\begin{array}{cc}{v}_f\left(\widehat{y}\right)& =v_{\mathit{fr}}\left(1-{\beta }_f{\widehat{y}}^P\right),0\le \widehat{y}\left(=\frac{y}{L}\right)\le 1\\ {\beta }_f& =\left(1-{∆}_f\right),{∆}_f=v_{\mathit{\text{ft}}}/v_{\mathit{fr}}\end{array}$$E4

where v_ft and v_fr 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.

$$v_f\left(\widehat{y}\right)=v_{\mathit{fr}}\left[{\beta }_f{\left(1-{\widehat{y}}^n\right)}^p+{∆}_f\right],\mathit{n}=\mathit{1},\mathit{2},\mathit{3}\mathit{p}\ge \mathit{0}$$E5

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 (V_f) and matrix (V_m) 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.

Property	Mathematical formula^*
Young’s modulus in direction (1) E11	Em Vm + E1f Vf
Young’s modulus in direction (2) E22	Em (1 + ξηVf)/(1-ηVf); η = (E2f –Em)/(E2f + ξEm)
Shear modulus G12	Gm (1 + ξηVf)/(1-ηVf); η = (G12f –Gm)/(G12f + ξGm)
Poisson’s ratio${\vartheta }_{12}$	${\vartheta }_m{V}_m+{\vartheta }_{12f}{V}_f$
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 V_m + V_f = 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, V_fo = V_mo = 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.

N_s = number of segments (panels).

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

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

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

${\widehat{L}}_j=(L_j/L_o$) = normalized length of the j-th segment.

$\widehat{L}=\left(L/L_o\right)={\sum }_{j=1}^{N_s}{\widehat{L}}_j$ = normalized total beam length.

${\widehat{H}}_j$ = $\left(H_j/H_0\right)$ = ${\widehat{H}}_j={\sum }_{k=1}^{N_L\left(j\right)}{\widehat{h}}_{\mathit{kj}}$= normalized total wall thickness of the j-th segment.

${\widehat{h}}_{\mathit{kj}}=\left(h_{\mathit{kj}}/H_0\right)$ = normalized thickness of the k-th layer in the j-th segment.

${\theta }_{\mathit{kj}}$ = fiber orientation angle in the k-th layer in the j-th segment.

${\widehat{\mathrm{\Gamma }}}_j={\mathrm{\Gamma }}_j/{\mathrm{\Gamma }}_o={\oint }_J\mathit{ds}/{\oint }_o\mathit{ds}$ = normalized circumference of the j-th segment cross section.

Γ_j = πD_j for circular C.S., Γ_j = 2(a_j + b_j) for rectangular C.S.

${\widehat{\rho }}_{\mathit{kj}}={\rho }_{\mathit{kj}}/{\rho }_o$ = normalized density of the k-th layer in the j-th segment.

$V_{f,\mathit{kj}}$ = fiber volume fraction in the k-th layer in the j-th segment.

${\rho }_{\mathit{kj}}$ = ${\rho }_f{V}_{f,\mathit{kj}}+{\rho }_m\left(1-V_{f,\mathit{kj}}\right)$, ${\rho }_o=$ 0.5(ρ_f + ρ_m).

ρ_f = fiber mass density, ${\rho }_m$ = matrix density.

${\widehat{m}}_j=\left(m_j/m_o\right)=$ normalized mass per unit length of the j-th segment.

$m_j={\mathrm{\Gamma }}_j{\sum }_{k=1}^{N_L\left(j\right)}{\rho }_{\mathit{kj}}{h}_{\mathit{kj}}$ = mass per unit length of the j-th segment, $m_o={\mathrm{\Gamma }}_{\mathrm{o}}{\rho }_o{H}_o$.

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

= ${\sum }_k^{N_L\left(j\right)}\oint {\rho }_{\mathit{kj}}{h}_{\mathit{kj}}\left(y^2+z^2\right)$ ds.

The normalized total structural mass is given by the expression:

$${\widehat{M}}_s=M_s/M_0=\sum _{j=1}^{N_s}{\widehat{M}}_j=\sum _{j=1}^{N_s}{\widehat{m}}_j{\widehat{L}}_j=\sum _{j=1}^{N_s}{\widehat{\mathrm{\Gamma }}}_j{\widehat{L}}_j\sum _{k=1}^{N_L\left(j\right)}{\widehat{\rho }}_{\mathit{kj}}{\widehat{h}}_{\mathit{kj}}$$E6

where $M_o=m_o{L}_o={\mathrm{\Gamma }}_{\mathrm{o}}{\rho }_o{H}_o{L}_o$ 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:

$$\left\{\begin{array}{c}{F}_x\\ M_x\\ M_y\\ M_z\end{array}\right\}=\left[\begin{array}{c}\begin{array}{cc}{C}_{11}& C_{12}\\ C_{12}& C_{22}\end{array}\begin{array}{cc}{C}_{13}& C_{14}\\ C_{23}& C_{24}\end{array}\\ \begin{array}{cc}{C}_{13}& C_{23}\\ C_{14}& C_{24}\end{array}\begin{array}{cc}{C}_{33}& C_{34}\\ C_{34}& C_{44}\end{array}\end{array}\right]\left\{\begin{array}{c}{U}_1^{\prime }\\ {\phi }^{\prime }\\ U_3^{″}\\ U_2^{″}\end{array}\right\}$$E7

where F_x, M_x, M_y , and M_z stand for the axial force, torsional, and bending moments, respectively, and C_mn are the cross-sectional stiffness coefficients derived in terms of closed-form integrals of the geometry and material constants. The notations U₁, U₂, U₃, 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 (U₁) is uncoupled, as well as the edgewise bending (U₂), while the flapping displacement (U₃) 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 (U₃), the associated eigenvalue problem can be written directly in the form:

$$C_{33,j}{U}_3^{\prime \prime \prime \prime }-{\omega }^2{m}_j{U}_3=0$$E8

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

$${\widehat{U}}_3^{\prime \prime \prime \prime }-{\widehat{\beta }}_j^4{\widehat{U}}_3=0$$E9

where ${\widehat{\beta }}_j$=$\sqrt{\widehat{\omega }}{\left({\widehat{m}}_j/{\widehat{C}}_{33,j}\right)}^{1/4}$, ${\widehat{C}}_{33,j}=C_{33,j}/C_{33,0}$, and $\widehat{\omega }=\omega L_0^2{\left(\frac{m_0}{C_{33,0}}\right)}^{1/2}$. Eq. (9) must be satisfied in the interval $0\le \overline{x}\le {\widehat{L}}_j,$ where $\overline{x}=\widehat{x}-{\widehat{x}}_j$ is a local coordinate of the j-th segment and $\widehat{x}=\left(x/L_o\right).$ The general solution is well known to be:

$${\widehat{U}}_3\left(\overline{x}\right)=a_1\mathit{sin}{\widehat{\beta }}_j\overline{x}+a_2\mathit{cos}{\widehat{\beta }}_j\overline{x}+a_3\mathit{sinh}{\widehat{\beta }}_j\overline{x}+a_4\mathit{cosh}{\widehat{\beta }}_j\overline{x}$$E10

Expressing the constants $a_i,i=1,2,3,4$ in terms of the state variables vector {S}^T = {$U_3-U_3^{\prime }-C_{33}{U}_3^{″}-C_{33}{U}_3^{\prime ″}\}{}^T$ at both ends of the j-th segment, we get

$${\left\{S\right\}}_{j+\mathit{1}}=\left[T^{\left(j\right)}\right]{\left\{S\right\}}_j$$E11

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:

$$\left[T\right]=\left[T^{\left(\mathit{Ns}\right)}\right]\left[T^{\left(\mathit{Ns}-\mathit{1}\right)}\right]\dots \dots ..\left[T^{\left(\mathit{2}\right)}\right]\left[T^{\left(\mathit{1}\right)}\right]$$E12

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, N_L = 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 = $\left(V_f,\widehat{H},\theta \right)$. The associated frequency equation for such a basic case is:

$$\mathit{cos}\widehat{\beta }\widehat{L}\mathit{cosh}\widehat{\beta }\widehat{L}=-1,\text{or}cos\sqrt{\widehat{\omega }}{\left(\widehat{m}/{\widehat{C}}_{33}\right)}^{\frac{1}{4}}\widehat{L}\mathit{cosh}\sqrt{\widehat{\omega }}{\left(\widehat{m}/{\widehat{C}}_{33}\right)}^{\frac{1}{4}}\widehat{L}=-1$$E13

It is seen that $\sqrt{\widehat{\omega }}$ 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 $\sqrt{\widehat{\omega }}{\left(\widehat{m}/{\widehat{C}}_{33}\right)}^{\frac{1}{4}}\widehat{L}$ without regard to the specific values of the design variables. The computed roots are:

$$\sqrt{{\widehat{\omega }}_i}={\left({\widehat{C}}_{33}/\widehat{m}\right)}^{\frac{1}{4}}\left(\frac{1}{\widehat{L}}\right)(1.8751,4.6941,7.8548,\dots \pi \left(i-0.5\right)i\ge 4$$E14

In Eq. (14), the frequency parameter $\sqrt{{\widehat{\omega }}_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, $\sqrt{{\widehat{\omega }}_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].

Property	Fiber: Carbon-AS4	Matrix: Epoxy-3501-6
Density (g/cm^3)	ρ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.0	Gm = 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 $\sqrt{{\widehat{\omega }}_1}$ combined with the structural mass constraint (${\widehat{M}}_s=1$). The imposed side constraints are:

$$\left(0.25,0.75,-\pi /2\right)\le \left(V_f,\widehat{H},\theta \right)\le \left(0.75,1.25,\pi /2\right)$$E15

Figure 5.

Level curves of $\sqrt{{\widehat{\omega }}_1}$ function augmented with the constraint ${\widehat{M}}_s=1$ in $\left(V_f-\theta \right)$ design space (Ns = 1, N_L = 1, $\widehat{L}=1)$.

It is remarked that the function is continuous and well behaved everywhere in (V_f – θ) 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 (V_f, $\widehat{H}$, θ) = (0.75, 0.92, 0) at which ${\left(\sqrt{{\omega }_1}\right)}_{\mathit{max}}=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 ($\sqrt{{\widehat{\omega }}_1}$=1.8751). The optimal solution was calculated to be (V_f, $\widehat{H}$, θ) = (0.50, 0.915, 0) and ${\widehat{M}}_{s,\mathit{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, N_L = 1), the reduced optimization problem can be defined as follows:

$$\begin{array}{c}\mathit{Minimize}F\left(\underset{¯}{X}\right)=–\left(\sqrt{{\widehat{\omega }}_1}\right)\\ \mathit{\text{Subject to}}{\widehat{M}}_s=1\\ \sum _1^{\mathit{Ns}}{\widehat{L}}_j=1\\ \left(\mathit{0.25},\mathit{0.0}\right)\le {\left(\mathit{Vfj},{\widehat{L}}_j\right)}_{j=1,2}\le \left(\mathit{0.75},\mathit{1.0}\right)\end{array}$$E16

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 (V_A −$\widehat{L}$)₁. 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 V_f1 = 50%. The constrained optimum is to found to be ($V_{\mathit{fj}},{\widehat{L}}_j){}_{j=1,2}$ = (0.75, 0.50), (0.25, 0.50) corresponding to ($\sqrt{{\widehat{\omega }}_1}$)_max = 2.0645 with 10.10% optimization gain.

Figure 6.

Level curves of $\sqrt{{\widehat{\omega }}_1}$ function augmented with ${\widehat{M}}_s=1$ in $\left(V_{f1},{\widehat{L}}_1\right)$ 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:

$$F\left(\underset{¯}{X}\right)=-\sqrt{{\widehat{\omega }}_1}$$

subject to the constraints:

$$\begin{array}{c}{\widehat{M}}_s=1\\ 0.33\le {\mathrm{∆}}_f\le 3.0\\ P\ge 0\end{array}$$E17

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 $\sqrt{{\widehat{\omega }}_1}$ combined with the mass constraint in the design space $\left({\mathrm{∆}}_f,p\right)$. 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., V_f = 50% at root and tip) split the domain into two zones. The lower zone encompasses the constrained optimum solution: ($\sqrt{{\widehat{\omega }}_1}$)_max = 2.01875 at the design point (∆_f, P)_opt. = (0.34, 1.01).

Figure 7.

Level curves of $\sqrt{{\widehat{\omega }}_1}$ function augmented with ${\widehat{M}}_s=1$ in $\left({\mathrm{∆}}_f,p\right)$ design space (Ns = N_L = 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.,	${\left(\sqrt{{\omega }_1}\right)}_{\mathit{max}}$,	gain %
Thickness grading (Eq. (3))		(1.0, 0.0),	1.8751	0.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 (${\widehat{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:

$$\begin{array}{c}{N}_{x,x}+N_{\mathit{yx},y}-2T{u}_{,\mathit{xy}}=0\\ N_{\mathit{xy},x}+N_{y,y}+\left(M_{\mathit{xy},x}/R\right)+\left(M_{y,y}/R\right)-2T\left(v_{,y}+w_{,x}/R\right)=0\\ M_{x,\mathit{xx}}+{\left(M_{\mathit{xy}}+M_{\mathit{yx}}\right)}_{,\mathit{xy}}+M_{y,\mathit{yy}}-N_y/R+2T\left(v_{,x}/R-w_{,\mathit{xy}}\right)=0\end{array}$$E18

where N_x and N_y are the normal forces, N_xy and N_yx are shear forces, M_x and M_y are bending moments, and M_xy and M_yx 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]:

$$T_{\mathit{cr}}=\left(2\pi R^{2}H\right)\left(0.272\right){\left(E_x\right)}^{0.25}{\left(E_y\right)}^{0.75}{\left(H/R\right)}^{1.5}$$E19

where T_cr 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 (E_x) and hoop (E_y) 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 [0^o/90⁰]_N with equal volume fractions of the fibers and matrix materials, that is, V_f = V_m = 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 $\overrightarrow{X}={\left(V_f,\theta ,\widehat{h}\right)}_{k=1,2,\dots N_L},$ which minimizes the objective function:

$$\begin{array}{ccc}& Minimize\hfill & F=-{\widehat{T}}_{\mathit{cr}}\\ \text{subject to}& \text{mass limitation}:& \widehat{M}-1\le 0\end{array}$$E20

$$\begin{array}{cc}\text{Torsional strength}:& \left(\frac{{\tau }_{\mathit{max}}}{{\tau }_{\mathit{\text{allow}}}}\right)-1.0\le 0\\ \text{Whirling}:\hfill & {\widehat{\mathrm{\Omega }}}_{\mathit{max}}-{\widehat{\mathrm{\Omega }}}_{\mathit{cr}}\le 0\end{array}$$E21

$$\begin{array}{cc}\text{Side constraints}:& \left(0.30,-\frac{\mathrm{\pi }}{2},0.015\right)\le {\left(V_f,\mathrm{\theta },\widehat{h}\right)}_{k=1,2,\dots N_L}\le \left(0.70,\frac{\mathrm{\pi }}{2},0.20\right)\\ & 0.75\le \sum _{k=1}^{N_L}{\widehat{h}}_k\le 1.25\hfill \end{array}$$E22

where ${\widehat{T}}_{\mathit{cr}}=T_{\mathit{cr}}/T_{\mathit{cro}}$, $\widehat{\Omega }=\left(\Omega \ast 2\pi \right)/\left(60{{\omega }_1}_{,o}\right)$ are the dimensionless critical torque and rotational speed, respectively. The baseline design parameters are denoted by subscript “o.” ${\tau }_{\mathit{max}}=\left(T_{\mathit{max}}/2\pi R^{2}H\right)$ is the maximum shear stress, T_max 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 (3N_L) 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 (V_f1-θ) design space, which are, as seen, monotonic and symmetric about the zero ply angle. A local maximum of ${\widehat{T}}_{\mathit{cr}}$ can be observed near the design point (V_f1,θ) = (0.7, 90^o) with ${\widehat{T}}_{\mathit{cr}}$ = 1. This figure illustrates that the maximum critical buckling torque can be achieved when the fiber orientation angle is close to 90^o. 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, N_L = 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 N_L = 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.

${\widehat{T}}_{\mathit{cr}}$—contours in (V_f1-θ) design space under mass constraint $\widehat{M}=1$. (Case of drive shaft with eight symmetric, balanced, carbon/ epoxy layers)

The final attained optimal solution was a cross ply layup [90⁰/0⁰]₄ 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]:

$${\mathit{C}}_{\mathit{B}}\frac{{\mathit{\partial }}^{\mathbf{2}}\mathit{v}}{{\mathit{\partial X}}^{\mathbf{4}}}+\mathit{\rho A}\frac{{\mathit{\partial }}^{\mathbf{2}}\mathit{v}}{{\mathit{\partial t}}^{\mathbf{2}}}-\left(\mathit{\rho I}+\frac{{\mathit{C}}_{\mathit{B}}}{{\mathit{C}}_{\mathit{S}}}\mathit{\rho A}\right)\frac{{\mathit{\partial }}^{\mathbf{4}}\mathit{v}}{{\mathit{\partial X}}^{\mathbf{2}}{\mathit{\partial t}}^{\mathbf{2}}}-\mathbf{2}\mathit{\rho I}\mathbf{\Omega }\left(\frac{{\mathit{\partial }}^{\mathbf{3}}\mathit{w}}{{\mathit{\partial x}}^{\mathbf{2}}\mathit{\partial t}}-\frac{\mathit{\rho A}}{{\mathit{C}}_{\mathit{S}}}\frac{{\mathit{\partial }}^{\mathbf{3}}\mathit{w}}{\mathit{\partial }{\mathit{t}}^{\mathbf{3}}}\right)+\frac{\mathit{\text{ρIρA}}}{{\mathit{C}}_{\mathit{S}}}\frac{{\mathit{\partial }}^{\mathbf{4}}\mathit{v}}{{\mathit{\partial t}}^{\mathbf{4}}}+\left({\mathit{C}}_{\mathit{BT}}/\mathbf{2}\right)\frac{{\mathit{\partial }}^{\mathbf{3}}\mathit{\varphi }}{{\mathit{\partial x}}^{\mathbf{3}}}=\mathbf{0}$$E23

$${\mathit{C}}_{\mathit{B}}\frac{{\mathit{\partial }}^{\mathbf{4}}\mathit{w}}{{\mathit{\partial x}}^{\mathbf{4}}}+\mathit{\rho A}\frac{{\mathit{\partial }}^{\mathbf{2}}\mathit{w}}{{\mathit{\partial t}}^{\mathbf{2}}}-\left(\mathit{\rho I}+\frac{{\mathit{C}}_{\mathit{B}}}{{\mathit{C}}_{\mathit{S}}}\mathit{\rho A}\right)\frac{{\mathit{\partial }}^{\mathbf{4}}\mathit{w}}{{\mathit{\partial x}}^{\mathbf{2}}{\mathit{\partial t}}^{\mathbf{2}}}+\mathbf{2}\mathit{\rho I}\mathbf{\Omega }\left(\frac{{\mathit{\partial }}^{\mathbf{3}}\mathit{v}}{{\mathit{\partial x}}^{\mathbf{2}}\mathit{\partial t}}-\frac{\mathit{\rho A}}{{\mathit{C}}_{\mathit{S}}}\frac{{\mathit{\partial }}^{\mathbf{3}}\mathit{v}}{{\mathit{\partial t}}^{\mathbf{3}}}\right)+\frac{\mathit{\text{ρIρA}}}{{\mathit{C}}_{\mathit{S}}}\frac{{\mathit{\partial }}^{\mathbf{4}}\mathit{w}}{{\mathit{\partial t}}^{\mathbf{4}}}+\left({\mathit{C}}_{\mathit{BT}}/\mathbf{2}\right)\frac{{\mathit{\partial }}^{\mathbf{3}}\mathit{\varphi }}{{\mathit{\partial x}}^{\mathbf{3}}}=\mathbf{0}$$E24

$$\left({\mathit{C}}_{\mathit{BT}}/\mathbf{2}\right)\left(\frac{{\mathit{\partial }}^{\mathbf{3}}\mathit{v}}{{\mathit{\partial x}}^{\mathbf{3}}}+\frac{{\mathit{\partial }}^{\mathbf{3}}\mathit{w}}{{\mathit{\partial x}}^{\mathbf{3}}}\right)-\left({\mathit{C}}_{\mathit{BT}}/\mathbf{2}\right)\frac{\mathit{\rho A}}{{\mathit{C}}_{\mathit{S}}}\left(\frac{{\mathit{\partial }}^{\mathbf{3}}\mathit{v}}{{\mathit{\partial x\partial t}}^{\mathbf{2}}}+\frac{{\mathit{\partial }}^{\mathbf{3}}\mathit{w}}{\mathit{\partial x}{\mathit{\partial t}}^{\mathbf{2}}}\right)+{\mathit{C}}_{\mathit{T}}\frac{{\mathit{\partial }}^{\mathbf{2}}\mathit{\varphi }}{{\mathit{\partial x}}^{\mathbf{2}}}-\mathit{\rho J}\frac{{\mathit{\partial }}^{\mathbf{2}}\mathit{\varphi }}{{\mathit{\partial t}}^{\mathbf{2}}}=\mathbf{0}$$E25

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 (C_B and C_BT) associated with bending about the y and z axes are identical; likewise for the transverse shear stiffness (C_s) [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:

$$\begin{array}{c}\left(C_{11}^2-C_{12}^2\right)C_{33}-2C_{11}{C}_{13}{C}_{31}=0\\ C_{11}=C_B{\lambda }^4-\left(\mathit{\rho I}+\frac{C_B}{C_S}\mathit{\rho A}\right){\omega }^2{\lambda }^2+\left(\frac{\mathit{\text{ρIρA}}}{C_S}{\omega }^2-\mathit{\rho A}{\omega }^2\right)\end{array}$$E26

$$\begin{array}{c}{C}_{12}=2\mathit{\rho I}\mathrm{\Omega }\omega \left({\lambda }^2-\frac{\mathit{\rho A}}{C_S}{\omega }^2\right);C_{13}=\left(C_{\mathit{BT}}/2\right){\lambda }^3\\ C_{31}=C_{13}-\left(C_{\mathit{BT}}/2\right)\left(\mathit{\rho A}/C_S\right)\lambda {\omega }^2;C_{33}=C_T{\lambda }^2-\mathit{\rho J}{\omega }^2\end{array}$$E27

where $\lambda =\mathit{n\pi }/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 $\overrightarrow{X}={\left(V_f,\mathrm{\theta },\widehat{h}\right)}_{k=1,2,\dots N_L}$, which minimizes the objective function:

$$\begin{array}{c}MinimizeF=-{\widehat{\Omega }}_{\mathit{cr}}\\ \text{Subject to}\widehat{M}-1\le 0\end{array}$$E28

$$\begin{array}{c}\left(\frac{{\tau }_{\mathit{max}}}{{\tau }_{\mathit{\text{allow}}}}\right)-1.0\le 0\\ \left(\frac{T_{\mathit{max}}}{T_{\mathit{cr},o}}\right)-{\widehat{T}}_{\mathit{cr}}\le 0\end{array}$$E29

(b) Placement of the critical speed

The other alternative of the objective function is defined by:

$$MinimizeF={\left({\widehat{\Omega }}_{\mathit{cr}}-{\widehat{\Omega }}^{\ast }\right)}^2$$E30

The same set of constraints given in Eq. (23) is applied. The notation ${\widehat{\Omega }}^{\ast }$ 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:

$${\mathit{V}}_{\mathit{f}}\left(\widehat{\mathit{x}}\right)={\mathit{V}}_{\mathit{f}}\left(\mathbf{0.5}\right)+\left[{\mathit{V}}_{\mathit{f}}\left(\mathbf{0}\right)-{\mathit{V}}_{\mathit{f}}\left(\mathbf{0.5}\right)\right]{\left({\left(\mathbf{1}-\mathbf{2}\left|\widehat{\mathit{x}}\right|\right)}^{\mathit{n}}\right)}^{\mathit{P}},$$E31

where V_f (0) is the fiber volume fraction at the right or left end of the drive shaft, while V_f (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 ${\widehat{\mathit{\Omega }}}_{\mathit{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 ${\widehat{\mathit{\Omega }}}_{\mathit{cr}}=1.045,$ is located inside the zone where the fiber taper ratio $∆f=V_f\left(0\right)/V_f\left(0.5\right)$ 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 $\overrightarrow{X}=\left(V_f\left(0\right),V_f\left(\frac{1}{2}\right),p,\widehat{H}\right)$ with lower and upper limits ${\overrightarrow{X}}_L$ = (0.3, 0.3, 0, 0.75) and ${\overrightarrow{X}}_U$ = (0.7, 0.7, ∞, 1.25). The attained optimal design variable vector was calculated to be ${\overrightarrow{X}}_{\mathit{opt}}=\left(0.7,0.3,5.61,0.955\right)$ at which the maximum critical speed increased by 14% above that of the baseline design with active mass constraint.

Figure 10.

Normalized critical speed ${\widehat{\Omega }}_{\mathit{cr}}$ augmented with the mass constraint ($\widehat{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.

$$\begin{array}{cc}Minimize& F=-\left({\widehat{\Omega }}_{\mathit{cr}}+{\widehat{T}}_{\mathit{cr}}\right)\\ \text{Subject to}& \widehat{M}-1\le 0\hfill \\ & \left(\frac{{\tau }_{\mathit{max}}}{{\tau }_{\mathit{\text{allow}}}}\right)-1.0\le 0\end{array}$$E32

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 ${\widehat{\Omega }}_{\mathit{cr}}=1.135$ and ${\widehat{T}}_{\mathit{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:

$$C\left(x\right)=C_r\left(\mathit{1}-{\beta }_{c}x\right),{\beta }_c=\left(\mathit{1}-{\Delta }_c\right)$$E33

The symbol Δ_c denotes the chord taper ratio (= tip chord C_t/root chord C_r) and x (= x₁/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]:

$$G=f_{\mathit{1}}{G}_{\mathit{12}}$$E34

where f₁ is a function that depends on the geometry and thickness ratio of the cross section (h/C) and the ratio (G₁₂/ G₁₃), where G₁₂ and G₁₃ 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, G₁₂ ≈ G₁₃.

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

$${\left(\widehat{G}\widehat{J}{\alpha }^{\prime }\right)}^{\prime }+{\widehat{V}}^2\alpha \left(x\right)=0$$E35

?The associated boundary conditions of the elastic angle of attack, α, are α(0) = 0 and ${\alpha }^{\prime }\left(1\right)=0$. The symbol $\widehat{G}=G_{12}/G_{12,o}$ denotes the dimensionless shear modulus, $\widehat{J}=J/J_r$ denotes the dimensionless torsion constant, and the prime denotes differentiation with respect to the dimensionless coordinate x = x₁/L. The dimensionless flight speed is defined by $\widehat{V}=V{C}_{r}b\sqrt{\mathit{\rho ae}/2G{J}_r}$, where (GJ)_r is the torsional stiffness of the baseline design at root. The shear modulus G_12,0 of the baseline design can be calculated by taking V_fo = 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:

$$z{\alpha }^{″}+3{\alpha }^{\prime }+a_k^2\alpha =0,\left(1‐\beta {\mathrm{x}}_{\mathrm{k}+1}\right)\le y\le \left(1-\beta {\mathrm{x}}_{\mathrm{k}}\right)$$E36

where ${\mathrm{a}}_{\mathrm{k}}=\widehat{\mathrm{V}}/\mathrm{\beta }\sqrt{{\widehat{\mathrm{h}}}_{\mathrm{k}}{\widehat{\mathrm{G}}}_{\mathrm{k}}}$, ${\widehat{\mathrm{h}}}_{\mathrm{k}}$ and ${\widehat{\mathrm{G}}}_{\mathrm{k}}$ are the normalized wall thickness and shear modulus of the kth wing panel, respectively. The general solution of Eq. (26) is:

$$\alpha \left(y\right)=A_1\frac{J_2\left(2\sqrt{a_{k}y}\right)}{y}-A_2\frac{Y_2\left(2\sqrt{a_{k}y}\right)}{y}$$E37

where J₂ and Y₂ are Bessel’s function of the first and second kind with order 2, respectively [35], and A₁ and A₂ 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 A₁ and A₂ 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:

$$\left\{\begin{array}{c}{\mathrm{\alpha }}_{\mathrm{k}+1}\\ {\mathrm{T}}_{\mathrm{k}+1}\end{array}\right\}=\left[{\mathrm{E}}^{\left(\mathrm{k}\right)}\right]\left\{\begin{array}{c}{\mathrm{\alpha }}_{\mathrm{k}}\\ {\mathrm{T}}_{\mathrm{k}}\end{array}\right\}$$E38

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:

$$\begin{array}{c}Minimize-{\widehat{V}}_{\mathit{div}}\\ \text{Subject to}{\widehat{M}}_s-1.0\le 0\\ \left(0.25,0.5,0.0\right)\le {\left(V_f,\widehat{h},\widehat{b}\right)}_{k=1,2,..\mathit{Np}}\le \left(0.75,1.25,1.0\right)\\ \sum {\widehat{b}}_k=1.0\end{array}$$E39

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 (V_div) 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, V_fo = 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 V_div 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 (V_f1- b₁) design space for a two-panel wing model (h₁ = h₂ = 1.0, Δ_c = 0.5,${\widehat{M}}_s=1$).

Figure 13.

Isodivert in (V_f1- b₁) design space for a two-panel wing model (h₂ = 0.5, V_f2 = 0.3, Δ_c = 0.5,${\widehat{M}}_s=1$).

Np	Type of grading	${\overrightarrow{X}}_{\mathit{opt}}={\left(V_f,\widehat{h},\widehat{b}\right)}_{k=1,2,..\mathit{Np}}$	${\widehat{V}}_{\mathit{div},\mathit{max}}$	Optimization gain %
2	Material	(0.5906, 1.0, 0.6594)1; (0.2500, 1.0, 0.3406)2	2.0523	3.93%
	Material and thickness	(0.75, 1.0625, 0.7)1; (0.25, 0.5250, 0.3)2	3.2275	63.44%
3	Material	(0.6000, 1.0, 0.5750)1; (0.4125, 1.0, 0.1375)2; (0.2500, 1.0, 0.2875)3	2.05625	4.13%
	Material and thickness	(0.75, 1.10, 0.65)1; (0.30, 0.65, 0.05)2; (0.25, 0.50, 0.3)3	3.3475	69.52%

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]:

$${\left(\mathit{EI}\right)}_k{V}^{\prime \prime \prime \prime }+U^2{V}^{\prime \prime }+2\mathit{i\omega U}\sqrt{{\beta }_o}{V}^{\prime }-m_k{\omega }^{2}V=0$$E40

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:

$$V\left(x\right)=\sum _{j=1}^4{A}_j{e}^{i{p}_{j}x}$$E41

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

$$p^4-\left({\alpha }_k{U}^2\right)p^2-\left(2\mathit{\omega U}\sqrt{{\beta }_o}{\alpha }_k\right)p-m_k{\omega }^2{\alpha }_k=0,{\alpha }_k=1/\left(E_k{h}_k\right)$$E42

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 ${\mathit{\alpha }}_{\mathit{k}}$, β_o, ω, and U. The transfer matrix [T_k] of the kth pipe module can be obtained by performing the matrix multiplications

$$\left[T_k\right]=\left[P^{\left(k\right)}\right]\left[E^{\left(k\right)}\right]{\left[P^{\left(k\right)}\right]}^{-\mathit{1}}$$E43

$$\left[p^{\left(k\right)}\right]=\left[\begin{array}{cc}\begin{array}{cc}1& 1\\ -{\mathit{ip}}_1& -{\mathit{ip}}_2\end{array}& \begin{array}{cc}1& 1\\ -{\mathit{ip}}_3& -{\mathit{ip}}_4\end{array}\\ \begin{array}{cc}{\left(\mathit{Eh}\right)}_k{p}_1^2& {\left(\mathit{Eh}\right)}_k{p}_2^2\\ i{\left(\mathit{Eh}\right)}_k{p}_1^3& i{\left(\mathit{Eh}\right)}_k{p}_2^3\end{array}& \begin{array}{cc}{\left(\mathit{Eh}\right)}_k{p}_3^2& {\left(\mathit{Eh}\right)}_k{p}_4^2\\ i{\left(\mathit{Eh}\right)}_k{p}_3^3& i{\left(\mathit{Eh}\right)}_k{p}_4^3\end{array}\end{array}\right]$$E44

and [E] is a diagonal matrix with elements:

$$E_{\mathit{jj}}^{\left(k\right)}=e^{i{p}_j^{\left(k\right)}{L}_k}$$E45

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:

$$T_{\mathit{33}}{T}_{\mathit{44}}-T_{\mathit{34}}{T}_{\mathit{43}}=\mathit{0},$$E46

where T_ij 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 U_f-MR_o 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 MR_o = 2.0. The calculated values of the dimensionless flutter velocity and frequency are found to be U_f = 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 (V_f, h₁) = (0.70, 0.9345), which corresponds to U_f = 12.517 and ω_f = 31.8615.

Figure 15.

Variation of flutter speed and frequency with mass ratio for a uniform one-module cantilevered pipe (V_f1 = 50%, h₁ = 1, L₁ = 1).

6.1.1. Solutions for cantilevered two-module pipe with uniform thickness (h₁ = h₂ = 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 (V_f1, 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 (U_f) and its corresponding flutter frequency (ω_f) occur in the region colored with dark brown L₁ = 0.36 and V_f1 = 0.3. The maximum U_f 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 (V_f, 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 (V_f1-L₁) design space.

Design variables: (Vf, L)k = 1,2	Dimensionless Uflutter	Stability improvement %
(0.30, 0.20), (0.55, 0.80); (0.30, 0.50), (0.70, 0.50)	13.15; 12.07	21.99; 11.97
(0.35, 0.25), (0.55, 0.75); (0.35, 0.40), (0.60, 0.60)	13.04; 13.31	20.96; 23.47 (max.)
(0.40, 0.50), (0.60, 0.50); (0.40, 0.60), (0.65, 0.40)	12.70; 12.77	17.81; 18.46
(0.45, 0.50), (0.55, 0.50); (0.45, 0.75), (0.65, 0.25)	12.82; 11.73	18.92; 8.81
(0.50, 0.50), (0.50, 0.50)	10.78	0.00 (baseline)

Table 5.

Standard solutions for a cantilevered two-module pipeline with uniform thickness (h₁ = h₂ = 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.