Abstract
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.
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.
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.
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.
(a) Iterations of displacement objective for volume constraint set at 90% of initial volume
(b) Efficiency curve of tip displacement and total volume
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.
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.
(a) Loading and geometry of depth tapered beam
(b) Displacement contours of designed beam. The mesh shows rectangular elements. A design variable is assigned per two elements
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].
(a) 1 element
(b) 2 elements
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.
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
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.
(a) 4-node element.
(b) 4-node element.
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.
(a) 8-node element
(b) 9-node element
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.
(a) . Iterations 60
(b) . Restart. Iterations 81
(a) . Iterations 80
(b) . Restart. Iterations 80
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.
(a) Design iterations
(b) Topology
(a) Mesh 10 × 30
(b) Mesh 40 × 120
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.
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].
(a)
(b)
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.
(a)
(b)
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.
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.
Acknowledgment
The author would like to thank Professor Robert L. Taylor at University of California, Berkeley, for valuable insights on the use of FEAP.