Table of Contents Author Guidelines Submit a Manuscript
Journal of Optimization
Volume 2015, Article ID 345120, 14 pages
Research Article

A Structural Optimization Framework for Multidisciplinary Design

Department of Engineering, The Pennsylvania State University, Altoona, PA 16601, USA

Received 29 September 2014; Revised 23 December 2014; Accepted 6 January 2015

Academic Editor: Qingsong Xu

Copyright © 2015 Mohammad Kurdi. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


This work describes the development of a structural optimization framework adept at accommodating diverse customer requirements. The purpose is to provide a framework accessible to the optimization research analyst. The framework integrates the method of moving asymptotes into the finite element analysis program (FEAP) by exploiting the direct interface capability in FEAP. Analytic sensitivities are incorporated to provide a robust and efficient optimization search. User macros are developed to interface the design algorithm and analytic sensitivity with the finite element analysis program. To test the optimization tool and sensitivity calculations, three sizing and one topology optimization problems are considered. In addition, flutter analysis of a heated panel is analyzed as an example of coupling to nonstructural discipline. In sizing optimization, the calculated semianalytic sensitivities match analytic and finite difference calculations. Differences between analytic designs and numerical ones are less than 2.0% and are attributed to discrete nature of finite elements. In the topology problem, quadratic elements are found robust at resolving checkerboard patterns.

1. Introduction

Increasingly engineering products are required to maintain high efficiency but yet perform well under conflicting customer requirements. Examples of these include thermal loads [1], combination of aeroelastic, thermal, and buckling loads, and acoustic and impact loads. Structural optimization strives to improve performance by interrogating multiple disciplines simultaneously. However, there is continuous need to consider new disciplines within existing optimization tools.

Industry recognized potential for structural optimization since the 1980s. Now there are multiple commercial structural optimization tools that address product design. These tools have emerged primarily from the academic research community and then commercialized to focus on a particular application. Early structural optimization implementations include adding sizing and shape optimization to finite element software such as NASTRAN, ANSYS, and ABAQUS [2, page 243]. Thomas et al. [3] indicates that early version of Altair OptiStruct is based on application of homogenization method for topology optimization [4]. GENESIS [5] focused on employing advanced approximation techniques [6, 7] and sensitivity tools to provide an efficient optimization search [8]. To obtain design sensitivities in ABAQUS, Yi et al. [9] application utilized response surface technique to approximate design derivatives for structural optimization of gear damper. The author utilized python to access the ABAQUS kernel and formulate the noise criterion.

Thomas et al. [3] highlights important requirements of structural optimization in commercial tools. Of these requirements, the ease of use of the interface and narrow application focus of the tool provide a disincentive to expand the tools to wider range of applications. For example, until now topology optimization is limited to structural problems and does not consider important physics [10] critical to the structure. On the other hand, research tools frequently utilize in-house codes to model a narrow physical problem and then couple to structural optimization. Extension of the methodology into commercial tools may require significant development.

It is desired to have a design optimization tool that offers development environment intermediate between research and commercial applications. The tool should be amenable to the use of robust optimization methods, analytic sensitivities, advanced approximate methods, and inclusion of all important physics. The development of analytic sensitivities for new criteria is particularly important for effectiveness of gradient-based optimization tools. In addition and as suggested by Thomas et al. [3], structural optimization software should have a broad analysis capability by containing different types of elements and are amenable to integration of new physics and material types. The finite element analysis program (FEAP) [11], developed by Professor Taylor at UC Berkeley, is a research finite-element code available in FORTRAN. It offers wide range of elements and encompasses structural analysis and thermal capabilities in linear and nonlinear and static and transient domains. The code is designed to incorporate new elements, physics, materials, meshes, and solution procedures by development of user macros. In addition, the macros lend an interfacing capability to an optimization program.

Several works exist on the utilization of FEAP in structural optimization considering aeroelastic constraints. Moon [12] conducted an aerostructural optimization of wings under divergence constraint. Alonso et al. [13, 14] demonstrated an optimization framework by coupling high-fidelity CFD code to FEAP. In particular in [13] finite difference tools are used to compute the structural sensitivities. In the above studies, a Python-based programming interface is used to wrap the fluid code to structural analysis tool. More recently, Schafer et al. [15] carried out a derivative-free shape optimization study by coupling a finite volume flow solver (FASTEST) to FEAP using the coupling interface (MpCCI).

In this work a structural optimization and analytic sensitivity framework is developed. The framework is built from within the FEAP command structure by accessing the user interface. The advantage of this is direct access to finite element matrices necessary for formulation of analytic gradients, flexibility in the choice of the numerical optimization method, and incorporation of new physics. Analytic gradients are derived and provided to the method of moving asymptotes [16] to expedite the optimization search (this method is utilized here as an example, and the user can select other gradient-based optimization techniques depending on the application). To handle competing objectives, a Pareto approach [17] multiobjective is utilized to address the multiobjective problem. The Pareto front is computed via the -constraint method. The paper proceeds in Section 2 to provide the general form of nonlinear structural system of equations. In Section 3, the analytic sensitivity equations are summarized. In Section 4, the implementation of the structural optimization framework and sensitivity analysis is described. In Section 5, numerical studies are carried out to demonstrate application of the framework to structural problems with sizing and density design variables. Finally Section 6 demonstrates incorporation of new physics by analyzing the buckling and flutter of a heated panel.

2. Structural Analysis

In static nonlinear analysis of structures, the finite element method gives a set of equilibrium equations [2, page 274]: where is a design variable, is a global response deformation, is the external applied load, is the internal load in structure due to the deformation, and is a load continuation parameter.

Newton method is used to compute the solution to system of equations. A Taylor series approximation of the residual is where is the tangent stiffness matrix due to internal and external loads as follows: The tangent stiffness matrix is calculated at the element level and then assembled for the global system. The matrix is symmetric when external loads are independent of deformations. In this case, direct solver is used to solve the system of equations for . The deformation is updated with the new solution until convergence to specified accuracy.

3. Sensitivity Analysis

For performance functions with dependence on deformation, one can write the sensitivity of to design variable : In order to compute the sensitivity using finite differences, one needs to recompute the response at perturbed value of design variable . In addition to computational cost, this sensitivity is prone to inaccuracies in the finite element solution. There are two alternative options for computing the sensitivity semianalytically [2, page 274]. The direct method computes the sensitivity from equilibrium equations (1) at a converged response : Forward finite difference is used to compute sensitivity of the residual. Note that this only requires build-up of the residual at new design. The factored tangent stiffness matrix is available from final iteration in (2). Equation (5) is repeatedly solved for all design variables. The cost of computation increases for optimization problems with large number of design variables. The adjoint method avoids this repetition by premultiplying (5) with an adjoint vector and adding to (4). The adjoint vector is computed by setting the coefficient of to zero: where refers to transpose of matrix. For structural systems where the external load is independent of deformation, the tangent stiffness matrix is symmetric and this operation can be skipped. Sensitivity of response is then computed as Note that (6) is repeatedly solved for each deformation dependent performance function. In optimization problems with large number of such functions, the direct method is more appropriate.

4. FEAP Implementation

The optimization framework is applied in the problem solution stage of the FEAP input file [11, page 3]. Here a direct interface (referred to as user macro [18]) to the program is developed to carry out the optimization search, sensitivity analysis, and definitions of the objective function and constraints. The macro allows access to the FEAP program pointers and defines new user pointers for use in the framework. New user commands are interweaved with existing solution commands in FEAP to achieve the algorithm shown in Figure 1, where a user command calls the user macro; see Table 1.

Table 1: Description of FEAP input file for structural optimization.
Figure 1: Optimization algorithm.

The user macro for the optimization search integrates the Fortran source code of the method of moving asymptotes (MMA) [16]. In addition to its availability, this method is well suited for structural optimization.

Due to complexity of data management in the sensitivity analysis portion of the algorithm, the basic structure of direct sensitivity calculation is explained in the following steps.(i)Calculate response at a specified design variable . This is implemented using Newton method. FEAP command consists of “TANG,,1” for linear analysis or a loop of same command for nonlinear analysis.(ii)The user macro saves response and design variable variable vector and residual into user designated pointers.(iii)Loop on design variable .(iv)User macro assigns saved user pointers of and into FEAP program pointers.(v)User macro perturbs a design variable , .(vi)Build the residual for the perturbed design. This is implemented by the FEAP command “FORM.”(vii)User macro calculates sensitivity of residual using forward finite difference. The saved residual in the user pointer is recalled here.(viii)User macro calculates the sensitivity of response due to change in , (5).(ix)Go to step (iii) for another design variable.

In case one desires employing the adjoint sensitivity, then first the adjoint vector is calculated for a performance function using (6). For this purpose, a user macro is constructed to calculate sensitivity performance to response . This sensitivity can be assigned as a fictitious load in the input file or in the corresponding pointer of FEAP memory. The latter is used here. The adjoint vector is then computed by (1) forming tangent stiffness matrix (“TANG”); (2) constructing residual (“FORM”); and (3) solving system of equations (“SOLVE”). Now the direct sensitivity algorithm is carried out but with step (viii) above being replaced with a user macro that implements (7). Note that a new adjoint vector is computed for each new design.

5. Numerical Studies

Three problems are studied regarding the sizing design of a cantilever beam with different loading conditions. Two problems have analytic solutions and serve as validation of the optimization framework [19, 20]. Then, a topology design problem is demonstrated.

The sizing optimization problems are arranged with regard to variation of beam cross-section along its axis. First the area is restricted to vary uniformly, maintaining a constant aspect ratio; second the beam depth is designed, maintaining constant width [2, page 38].

5.1. Uniformly Tapered Beam

An Euler-Bernoulli beam element is used to represent the deformation. The deformation vector consists of axial, lateral, and slope at each node. The design variables are the cross-sectional area of the beam elements. Because separate variables define the cross-sectional area and moment of inertia in FEAP, the moment of inertia is recalculated after each design update. Table 2 provides properties common to the uniformly tapered beam.

Table 2: Length and material of beam studied in uniformly tapered beams.

5.1.1. Tip Deflection Objective and Tip Load

Consider a cantilever beam subjected to a static tip load; see Figure 2(a). The objective problem of minimizing the tip deflection and weight is formulated as subject to where is the cross-section area of each beam element, and are the lower and upper limits on the area design variables, is the lateral deflection at the free end of the beam, is the volume of beam, is the axial stress, and subscript refers to a scaling value for either objective function or constraints. Scaling value is computed at initial uniform beam cross-section. are set of limits imposed on volume constraint to generate the efficiency curve and is a factor for setting the maximum stress. Scaling is employed here to achieve an effective optimization search. Due to this scaling, the reported results are generally dimensionless.

Figure 2: Loading schematic for uniform tapered beam design.

Sensitivity of the objective is computed from (4), where first explicit term is zero and Therefore, the objective sensitivity consists only of the derivative of the response at tip in lateral direction. Computed sensitivity is checked by analytic derivation for one or two elements [2, page 265]. Sensitivity of the constraints is computed explicitly.

The beam has an initial uniform cross-sectional area of m2 and is allowed to change in the interval of  m2. The stress factor in (8c) is set to 3.0. Ten design variables are assigned to ten elements in the beam. For , the optimization search converges to within 8 iterations. Only the volume constraint is found active; see Figure 3(a). The efficiency curve of displacement and volume is reported in Figure 3(b). At the right end of the curve, the displacement is minimized with the volume constrained to be lower than initial one. At left end the displacement is minimized with volume constrained to 30% of the initial volume. At same volume as initial design the displacement is reduced by 19%. Additionally a 20% reduction in volume is achieved for tip displacement equal to initial design.

Figure 3: Tip displacement objective for beam with tip load.
5.1.2. Compliance Objective and Distributed Load

The cantilever beam considered here is subjected to a uniform distributed load; see Figure 2(b). The objective is to decrease total area under deformation curve without increasing the volume: subject to The compliance objective is defined as The integral is evaluated by approximating with cubic Hermite shape functions and summing over all elements: where is the length of each element. The sensitivity of the objective to the deformation for element is Sensitivity to design variables is evaluated by summing (13) over all elements then applying (4). The efficiency curve is computed starting from beam of uniform cross-section area of  m2. Lower and upper limits of the area design variables are 0.001 m2 to 0.05 m2, respectively. Ten elements model the deformation and one design variable is assigned to each element. Constraint on the volume is varied between and  in. increments of 0.1; see Figure 4. For the same weight, the stiffness of the beam improved by 2.5 times (). This is within 2.5% of the analytic result reported in [19], where the authors calculated . The difference decreases with increase in the number of design variables.

Figure 4: Efficiency curve of area under lateral deformation and total volume for beam with uniform distributed load.
5.2. Depth Tapered Beam

The cantilever beam considered has a tip moment, geometry, and dimensions as indicated in Figure 5(a). Use of U.S. customary units is favored here for direct comparison to the analytical results. Haug and Kirmser [20] derived analytic thickness design for minimum volume under stress and tip displacement constraints. A tip displacement objective is retained as subject to is the plate thickness design variables, is the lower limit, and is the upper limit. Plate elements are utilized to model the deformation [21]. The mesh of the plate is generated using block command in FEAP; see Figure 5(b). Each block corresponds to a distinct material number. In Figure 5(b), a total of twenty materials or design variables are generated. The materials are then merged together to form the complete structure. Coincident nodes are removed in this process.

Figure 5: Depth tapered beam example.

Sensitivity of the tip deflection with respect to thickness design variables is reported in Figure 6. The sensitivity is seen to be affected by fixed end boundary condition; see Figure 6(a) in the plate element. To improve accuracy of the sensitivity, more elements need to be used per design variable [22, 23]. The sensitivity is recomputed with more elements per design variable; see Figure 6(b). Smoother sensitivity is seen near the fixed end. Alternatively one can eliminate this by setting to retrieve the beam stiffness [24, page 533].

Figure 6: Comparison between analytic and finite difference sensitivities when number of elements per design variable is increased. in and .

The analytic minimum beam volume reported in [20] is specified for a maximum tip displacement of 0.5 in. (0.0127 m), maximum stress of 30,000 psi (207 MPa), and lower and upper limit on the design variables of 0.3 and 0.5 in., respectively. The material has  psi (69 GPa) and . A numerical design is also calculated in [6] with fewer optimization iterations by exploiting improved approximation techniques. The maximum volume is set to the minimum analytic one (). Comparison of results is provided in Table 3. Good agreement is found.

Table 3: Final design of cantilever beam with end moment.

5.3. Topology of Cantilever Structure

The topology design of the extensively studied MBB problem is analyzed [25]. The beam has a length-to-height ratio of 6 and is simply supported and loaded at its center. Due to symmetry, half of the beam is considered in the optimization; see Figure 7. The optimization problem is to minimize the compliance of the structure at specified limit of material volume: subject to The solid isotropic material with penalization (SIMP) is utilized [26] in the topology design. In this approach, the design variables correspond to the volume density of a material element. The volume densities are linked to mechanical properties of interest. The elastic modulus is normally chosen for compliance objective: where is the elastic modulus of the material and is minimum stiffness to avoid singularity. Or may be set as the modulus of a secondary material. is a penalty parameter selected to yield binary values of the densities (0 or 1). The compliance is calculated according to

Figure 7: Half of MBB beam discretized into 12 square elements. Node 42 restrained in vertical direction. Nodes 1, 3, and 7 restrained in horizontal direction. Load is applied at node 7.

In anticipation of large number of design variables in topology optimization, the adjoint method (7) is utilized to compute sensitivity of compliance objective. The adjoint vector is found from (6) by setting . Although the finite difference step used in computing the residual sensitivity must be selected carefully depending on the mesh size and penalty parameter , sensitivity of volume constraint is computed explicitly.

The design domain is discretized using two dimensional square blocks. Each block corresponds to a distinct material or design variable. An example of design domain is constructed with 12 blocks; see Figure 7. Displacement of the beam is computed under plane stress assumption. Parameters of the beam considered are , , , and .

The strict application of (16) in conjunction with 4-node elements per design variable often results in designs that are marred by checkerboard patterns. Díaz and Sigmund [27] showed that designed structures possess a fictitiously higher stiffness when the patterns are present. An example of this is demonstrated in Figure 8.

Figure 8: Designed topology with a volume fraction of 0.5. A mesh of design variables 20 × 60 is utilized. Each design variable corresponds to a 4-node element. 29 iterations of optimization are carried out.

Standard methods used to circumvent checkerboards are use of higher-order elements and filters [28]. Here advantage is made of elements library in FEAP to resolve the numerical difficulty. 8-node (Q8) and 9-node (Q9) elements are studied, providing extra degrees of freedom per design variable; see Figure 9.

Figure 9: Schematic of high order elements.

In contrast to the linear 4-node elements, these elements provide quadratic interpolation of the displacement field. The presence of additional degrees of freedom per design variable allows for better estimate of the displacement field and its sensitivity, especially for penalized formulation of the SIMP method. To arrive at designs with high efficiency, continuation method is employed on the penalty parameter. Topologies are computed for mesh of 26 × 78 design variables for Q8 and Q9 elements; see Figures 10 and 11, respectively.

Figure 10: Volume fraction of 0.5. Mesh 26 × 78. Q8 elements. .
Figure 11: Volume fraction of 0.5. Mesh 26 × 78. Q9 elements. .

The higher-order elements are found to almost eliminate checkerboards from the design. A topology with finer design variable mesh is computed; see Figure 12. The objective reaches a minimum in around 20 iterations with remaining iterations devoted to removing intermediate densities. The use of sensitivity filter [26] also eliminates the checkerboard problem and provides a design independent of refinement in the mesh size; see Figure 13. Intermediate densities are now transferred to near the boundary of the material and can pose interpretation difficulties for manufacturing.

Figure 12: Volume fraction of 0.5. Mesh 33 × 99. Q9 elements. . Iterations 168. .
Figure 13: Volume fraction of 0.5. Q4 elements. .

6. Buckling and Flutter Analysis of a Heated Panel

Panel flutter is studied by many researchers in aerospace applications. See, for example, Stanford and Beran [29] and references therein. A heated panel is studied here that is fixed at both ends and subjected to aerodynamic flow on the upper surface; see Figure 14.

Figure 14: Panel geometry showing applied loads.

The upper and lower surfaces are heated from an initially unstressed reference temperature. The heating results in thermal stresses that may cause buckling and/or decrease resistance to lateral aerodynamic load. To analyze the panel, it is necessary to account for the aerodynamic and thermal loads in the structural model. Once these are developed, an eigenvalue analysis of the dynamic system is utilized to calculate the flutter and buckling metrics.

6.1. Panel Model

Equations of motion are derived based on principle of virtual work by balancing external work with elastic, inertial, and structural dissipation [24, page 375]. The work balance for an element reduces to Here , , and are the element mass, damping, and stiffness matrices, is the stress stiffness matrix resulting from thermal stress [24, page 642], and is the nodal deformation consisting of axial , lateral , and rotational degrees of freedom. is the external nodal load due to aerodynamic flow. In the following, the aerodynamic load and stress stiffness matrix are derived for an Euler-Bernoulli beam element.

6.2. Aerodynamic Load

The structural panel is exposed to quasi-steady supersonic flow. Piston theory is used to model the pressure on the top surface of the panel. The pressure is given by Dowell [30]: where is the flow velocity, is Mach number, is the dynamic pressure, is the air density, and . The surface pressure is resolved into a work-equivalent nodal load via the finite element shape functions [24, page 111]: Here are composed of linear shape functions that describe axial displacements and Hermitian shape functions that describe lateral deflection . The scalar pressure is derived by inserting the approximate lateral deformation into (19). Decomposing the aerodynamic load into elastic and damping components, The pressure acts normally to panel in opposite direction to the lateral deformation. For linear analysis, this direction can be assumed constant: The equivalent elastic nodal load becomes and aerodynamic damping load becomes The aerodynamic stiffness and damping matrices are found by taking derivative of aerodynamic forces to deformations [31]: Using MATLAB symbolic toolbox, the aerodynamic matrices are computed: Here zeros are added to the axial degree of freedom to match size of matrices to the element stiffness matrix.

6.3. Thermal Load

The effect of the thermal load is to reduce the overall stiffness of the panel by developing axial compressive stresses. This has adverse effect on buckling and flutter constraints. The net compressive stress in panel is computed by accounting for thermal strain panel [24, page 52]. The process of computing the thermal stress is summarized by (1) calculating a consistent nodal thermal load for an element; (2) assembling thermal nodal loads and conventional elastic stiffness matrix; (3) calculating thermal nodal deformations by solving where is the global stiffness matrix and is thermal load vector; (4) updating stresses in each element with thermal deformations; and (5) finally computing an element stress stiffness matrix.

The thermal load is computed assuming a prescribed steady state thermal boundary condition [29]. Then by applying Fourier law in one dimension, the temperature profile in the panel is computed: Here refers to the increase in temperature from a reference temperature with no thermal stress. The thermal axial force is found by integrating the initial thermal stress over the cross-sectional area [32]: where is the coefficient of thermal expansion, and the thermal bending moment is found by integrating the bending stress over the cross-sectional area [33, page 406]: A consistent nodal force and moment vector are computed by inserting the axial and bending loads into [24, page 90] where is the strain displacement matrix, giving Thermal deformations are computed from (27) and the axial thermal stress in an element with nodes is then computed as A stress stiffness matrix is calculated based on the in-plane membrane stresses [24, page 642]:

6.4. Equations of Motion

The element equation of motion (18) can be written in terms of aerodynamic damping and stiffness matrices: Assembling (35) for all elements, the equation of motion for the panel becomes To facilitate stability analysis, (36) is written in first order form: where and are the total stiffness and damping matrices, respectively. One may obtain an eigenvalue problem of (37) by assuming an exponential solution: where is the right eigenvector corresponding to the eigenvalue . Solution of the eigenvalue problem of this system determines its stability characteristics [34, page 248]. The eigenvalue solver must be selected to handle nonsymmetry of due to the aerodynamic stiffness matrix. The nonsymmetric eigenvalue driver DGGEVX.F is utilized here [35].

For the unheated panel and zero flow condition, solution of (38) for the uniform panel gives the natural modes of the system with the first mode given as where is the flexural rigidity of the panel. For the heated panel and zero flow condition, solution of (38) results in the critical buckling temperature . This temperature is nondimensionlized as For a thermal condition of , the first buckling temperature yields .

6.5. Flutter and Buckling Results

The onset of flutter is typified by emergence of a Hopf bifurcation [36] as a flow speed parameter is increased. Defining the parameter as in [30]: At the Hopf point, a pair of complex conjugate eigenvalues of the system cross the imaginary axis resulting in a destabilizing positive damping. To locate the flutter boundary, the damping of the least stable mode, is tracked to identify when .

An example of flutter computation is shown in Figure 15. The panel has a mass ratio . Five elements are used to analyze the structure and apply the aerodynamic load. No structural damping is considered. For the unheated panel and at , the eigenvalues yield the natural frequencies of the panel. As the flow parameter is increased, the frequencies of the first and second modes are seen to come closer. Near , their frequencies coalesce into one frequency at which point the damping of first mode becomes positive. This flutter point is slightly different from reported in [37] as the author discounted the aerodynamic damping term reducing model to static strip theory [38].

Figure 15: Depiction of system eigenvalues as flow parameter is increased for the unheated panel. (a) Frequencies. (b) Damping.

Tracing the eigenvalues of a heated panel at a uniform temperature of larger than , one can observe that the panel buckles due to the thermal stress at winds off. This is identified by zero buckling frequency and the positive damping of first mode; see Figure 16. As the aerodynamic load is increased beyond , the panel is blown stable as in [29], until a critical is reached, where the panel loses stability to flutter.

Figure 16: Depiction of buckling and flutter of a heated panel. (a) Frequencies. (b) Damping.

In a similar manner, the flutter and buckling boundaries can be traced in the domain of ; see Figure 17. At each heating load the flutter and buckling points are located. The boundary studied in [29] is calculated here.

Figure 17: Flutter and buckling boundaries.

The process of identifying the flutter point can be made faster with use of a root finding technique. The false position is a root bracketing method that can be used to facilitate locating of the flutter point by iterating [39, page 129] where and are the lower and upper values of the flow parameter, respectively. The method is based on locating positive and negative values of a function at opposite sides of the root. For zero structural damping, the systems damping is either positive or zero. Therefore, a small structural damping must be introduced for the method to work. In this work, Rayleigh type structural damping is used: where constants and are evaluated by assuming a constant damping factor of 0.01% for the first and 4th modes [40, page 455]. Recalculation of the flutter point of the unheated panel gives a nearly identical result as above (0.3% difference).

7. Conclusions

A structural optimization framework is developed that combines finite element analysis capability of FEAP with a gradient-based optimization search. Analytic sensitivities are incorporated into the framework to provide design gradients. The optimization algorithm is integrated directly with FEAP memory management system through the use of user macros, providing efficient execution of the optimization search.

Four test problems are analyzed using the framework with regard to sizing and topology optimization. In area of sizing design variables, the problems are found to give solutions identical to analytic designs with 1-2% difference due to discrete nature of finite element. In topology optimization, it is found that 8-node and 9-node square elements are robust in handling checkerboard patterns in the design, with slight advantage of the latter element. Finally buckling and flutter analysis of a heated panel is developed to demonstrate capability of introducing new physics into the structural analysis.

There are many advantages to the framework. First the framework provides easy access to all capabilities of FEAP program, which include thermal and mechanical analysis with a variety of elements, in addition to capability of interfacing with a parallel sparse matrix solver. Second because of the direct integration, the structural models can be extended to large problems. Third the user may expand on existing tools in FEAP by deriving his own elements and macros to analyze new physics.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.


The author would like to thank Professor Robert L. Taylor at University of California, Berkeley, for valuable insights on the use of FEAP.


  1. B. Chen and L. Tong, “Thermomechanically coupled sensitivity analysis and design optimization of functionally graded materials,” Computer Methods in Applied Mechanics and Engineering, vol. 194, no. 18-20, pp. 1891–1911, 2005. View at Publisher · View at Google Scholar · View at Scopus
  2. R. T. Haftka and Z. Gürdal, Elements of Structural Optimization, Kluwer Academic Publishers, Dordrecht, The Netherlands, 1992. View at MathSciNet
  3. H. Thomas, M. Zhou, and U. Schramm, “Issues of commercial optimization software development,” Structural and Multidisciplinary Optimization, vol. 23, no. 2, pp. 97–110, 2002. View at Google Scholar
  4. M. P. Bendsøe and N. Kikuchi, “Generating optimal topologies in structural design using a homogenization method,” Computer Methods in Applied Mechanics and Engineering, vol. 71, no. 2, pp. 197–224, 1988. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  5. H. L. Thomas, Y. K. Shyy, G. N. Vanderplaats, and H. Miura, “Genesis—a modern shape and sizing structural design tool,” in Proceedings of the 33rd Structures, Structural Dynamics and Materials Conference, AIAA Paper, 1992, View at Publisher · View at Google Scholar
  6. G. N. Vanderplaats and H. L. Thomas, “An improved approximation for stress constraints in plate structures,” Structural Optimization, vol. 6, no. 1, pp. 1–6, 1993. View at Publisher · View at Google Scholar · View at Scopus
  7. G. N. Vanderplaats and E. Salajegheh, “An efficient approximation technique for frequency constraints in frame optimization.International journal for numerical methods in engineering,” International Journal for Numerical Methods in Engineering, vol. 26, no. 5, pp. 1057–1069, 1988. View at Publisher · View at Google Scholar · View at Scopus
  8. R. Grandhi, “Structural optimization with frequency constraints—a review,” AIAA Journal, vol. 31, no. 12, pp. 2296–2303, 1993. View at Publisher · View at Google Scholar · View at Scopus
  9. G. Yi, Y. Sui, and J. Du, “Application of python-based abaqus preprocess and postprocess technique in analysis of gearbox vibration and noise reduction,” Frontiers of Mechanical Engineering, vol. 6, no. 2, pp. 229–234, 2011. View at Google Scholar
  10. J. D. Deaton and R. V. Grandhi, “A survey of structural and multidisciplinary continuum topology optimization: post 2000,” Structural and Multidisciplinary Optimization, vol. 49, no. 1, pp. 1–38, 2014. View at Google Scholar · View at MathSciNet · View at Scopus
  11. R. L. Taylor, FEAP—Finite Element Analysis Program: User Manual, University of California, Berkeley, Calif, USA, 2014,
  12. S. Moon, Aero-structural optimization of divergence-critical wings [Ph.D. thesis], University of Toronto, 2009.
  13. J. J. Alonso, J. R. R. A. Martins, J. J. Reuther, R. Haimes, and C. A. Crawford, “High-fidelity aero-structural design using a parametric cad-based model,” AIAA Paper 3429:2003, AIAA, 2003. View at Google Scholar
  14. J. J. Alonso, P. LeGresley, E. van der Weide, J. R. R. A. Martins, and J. J. Reuther, “pyMDO: a framework for high-fidelity multi-disciplinary optimization,” AIAA Paper 4480:2004, AIAA, 2004. View at Google Scholar
  15. M. Schafer, D. C. Sternel, G. Becker, and P. Pironkov, “Efficient numerical simulation and optimization of fluid-structure interaction,” in Fluid Structure Interaction II, pp. 131–158, Springer, New York, NY, USA, 2010. View at Google Scholar
  16. K. Svanberg, “The method of moving asymptotes—a new method for structural optimization,” International Journal for Numerical Methods in Engineering, vol. 24, no. 2, pp. 359–373, 1987. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  17. K. Sebaa, A. Tlemçani, M. Bouhedda, and N. Henini, “Multiobjective optimization using cross-entropy approach,” Journal of Optimization, vol. 2013, Article ID 270623, 9 pages, 2013. View at Publisher · View at Google Scholar
  18. R. L. Taylor, FEAP—A Finite Element Analysis Program, University of California, Berkeley, Calif, USA, 2014,
  19. S. M. Makky and M. A. Ghalib, “Design for minimum deection,” Engineering Optimization, vol. 4, no. 1, pp. 9–13, 1979. View at Publisher · View at Google Scholar · View at Scopus
  20. J. E. Haug and P. G. Kirmser, “Minimum weight design of beams with inequality constraints on stress and deflection,” Journal of Applied Mechanics, vol. 34, no. 4, pp. 999–1004, 1967. View at Publisher · View at Google Scholar
  21. F. Auricchio and R. L. Taylor, “A shear deformable plate element with an exact thin limit,” Computer Methods in Applied Mechanics and Engineering, vol. 118, no. 3-4, pp. 393–412, 1994. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  22. B. Prasad and R. T. Haftka, “Optimal structural design with plate finite elements,” Journal of the Structural Division, vol. 105, no. 11, pp. 2367–2382, 1979. View at Google Scholar · View at Scopus
  23. R. T. Haftka and R. V. Grandhi, “Structural shape optimization—a survey,” Computer Methods in Applied Mechanics and Engineering, vol. 57, no. 1, pp. 91–106, 1986. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  24. R. D. Cook, D. S. Malkus, M. E. Plesha, and R. J. Witt, Concepts and Applications of Finite Element Analysis, Wiley, New York, NY, USA, 2002.
  25. N. M. Patel, D. Tillotson, J. E. Renaud, A. Tovar, and K. Izui, “Comparative study of topology optimization techniques,” AIAA Journal, vol. 46, no. 8, pp. 1963–1975, 2008. View at Publisher · View at Google Scholar · View at Scopus
  26. M. P. Bendsoe and O. Sigmund, Topology Optimization: Theory, Methods and Applications, Springer, Berlin, Germany, 2003.
  27. A. Díaz and O. Sigmund, “Checkerboard patterns in layout optimization,” Structural Optimization, vol. 10, no. 1, pp. 40–45, 1995. View at Publisher · View at Google Scholar
  28. O. Sigmund and J. Petersson, “Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima,” Structural Optimization, vol. 16, no. 1, pp. 68–75, 1998. View at Publisher · View at Google Scholar · View at Scopus
  29. B. Stanford and P. Beran, “Aerothermoelastic topology optimization with flutter and buckling metrics,” Structural and Multidisciplinary Optimization, vol. 48, no. 1, pp. 149–171, 2013. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  30. E. H. Dowell, “Nonlinear oscillations of a fluttering plate,” AIAA journal, vol. 4, no. 7, pp. 1267–1275, 1966. View at Publisher · View at Google Scholar
  31. D. N. Herting and G. W. Haggenmacher, “Pressure follower matrix for geometric nonlinear finite element analysis,” in Proceedings of the MSC/Nastran User's Conference, 1987.
  32. R. B. Hetnarski and M. R. Eslami, Thermal Stresses—Advanced Theory and Applications, vol. 158, Springer, 2008.
  33. J. M. Gere and S. P. Timoshenko, Mechanics of Materials, PWS Engineering, Boston, Mass, USA, 1987.
  34. S. H. Strogatz, Nonlinear Dynamics and Chaos: with Applications to Physics, Biology, Chemistry, and Engineering, Westview Press, Boulder, Colo, USA, 2000.
  35. E. Anderson, Z. Bai, C. Bischof et al., LAPACK Users' Guide, vol. 9, SIAM, Philadelphia, Pa, USA, 1999.
  36. R. L. Bisplinghoff, H. Ashley, and R. L. Halfman, Aeroelasticity, Addison-Wesley, 1955.
  37. L. L. Erickson, Supersonic flutter of at Rectangular Orthotropic Panels Elastically Restrained against Edge Rotationn, National Aeronautics and Space Administration, Washington, DC, USA, 1966.
  38. C. Mei, K. Abdel-Motagaly, and R. Chen, “Review of nonlinear panel flutter at supersonic and hypersonic speeds,” Applied Mechanics Reviews, vol. 52, no. 10, pp. 321–332, 1999. View at Publisher · View at Google Scholar
  39. S. C. Chapra and R. P. Canale, Numerical Methods for Engineers: With Programming and Software Applications, McGraw-Hill, 2002.
  40. A. K. Chopra, Dynamics of Structures, vol. 3, Prentice Hall, Upper Saddle River, NJ, USA, 1995.