Topology Optimization as a Conceptual Tool for Designing New Airframes

(1)

Topology Optimization as a Conceptual Tool

for Designing New Airframes

Joakim Torstensson

ISRN: LIU-IEI-TEK-A16/02590SE

Division of Solid Mechanics | Department of Management and Engineering Spring 2016 | Linköping University, SE-581 83 Linköping, Sweden

Linköping University | 013-28 10 00 | www.liu.se June 20, 2016

(2)

Topology Optimization as a Conceptual Tool

for Designing New Airframes

Joakim Torstensson

Supervisors : Erik Holmberg, Saab AB

Anders Klarbring, Linköping University Examiner : Carl-Johan Thore, Linköping University

ISRN: LIU-IEI-TEK-A16/02590SE

Division of Solid Mechanics | Department of Management and Engineering Spring 2016 | Linköping University, SE-581 83 Linköping, Sweden,

Linköping University | 013-28 10 00 | www.liu.se June 20, 2016

(3)

This masters thesis concludes my studies at Linköping University and is the nal piece in obtaining my Master of Science in Mechanical Engineering.

I would like to express my gratitude towards my supervisor at Saab, Erik Holmberg, for the support, guidance and interesting discussions all along. I would also like to thank the people supporting me at the Division of Solid Mechanics at Linköping University: my supervisor Anders Klarbring and examiner Carl-Johan Thore for help all along.

Additionally, the people at the Altair support deserve a mention for the help they have pro-vided related to the OptiStruct software.

Joakim Torstensson Linköping, June 2016

(4)

During the two last decades, topology optimization has grown to be an accepted and used method to produce conceptual designs. Topology optimization is traditionally carried out on a component level, but in this project, the possibility to apply it to airframe design on a full scale aeroplane model is evaluated.

The project features a conceptual ying-wing design on which the study is to be carried out. Inertia Relief is used to constrain the aeroplane instead of traditional single point constraints with rigid body motion being suppressed by the application of accelerations instead of tradi-tional forces and moments. The inertia relief method utilized the inertia of the aeroplane to achieve a state of quasi-equilibrium such that static nite element analysis can be carried out. Two load cases are used: a steep pitch-up manoeuvre and a landing scenario. Aerodynamic forces are calculated for the pitch-up load case via an in-house solver, with the pressure being mapped to the nite element mesh via a Matlab-script to account for dierent mesh sizes. Increased gravitational loads are used in the landing load case to simulate the dynamic loading caused in a real landing scenario, which is unable to be accounted for directly in the topology optimization.

It can be concluded that the optimization is unable to account for one of the major design limitations: buckling of the outer skin. Approaches to account for the buckling of the outer skin are introduced and analysed, with a focus on local compression constraints throughout the wing. The compression constraints produce some promising results but are not without major drawbacks and complications.

In general, a one-step topology optimization to produce a mature conceptual airframe design is not possible with optimization algorithms today. It may be possible to adopt a multiple-step optimization approach utilizing topology optimization with following size and shape optimiza-tion to achieve a design, which could be expanded on in a future project.

(5)

List of Figures

1 Simple topology optimization example . . . 2

2 Conceptual ying-wing unmanned aerial vehicle model . . . 3

3 Overview of internal cavities . . . 4

4 Finite element discretization . . . 6

5 Solution with intermediate densities . . . 9

6 Checkerboard visualization . . . 10

7 Mesh size dependency . . . 11

8 Gravitational eects on topology optimization results . . . 13

9 Gravity eects on the TO results using a minimum lower volume constraints and point masses . . . 14

10 Optimization problem using SPCs . . . 16

11 Showcasing the eect inertia relief on the optimal design . . . 17

12 Problematic geometries . . . 19

13 Geometry clean-up of the aeroplane nose . . . 20

14 Meshing aspects . . . 20

15 Pressure distribution (in [Pa]) on the upper part of the body . . . 22

16 Pressure distribution (in [Pa]) on the lower part of the body . . . 22

17 SPC locations . . . 23

18 Constrained nodes for the inertia relief analysis . . . 24

19 Dierent optimal designs by density cut-o thresholds . . . 26

20 Wing-only model . . . 27

21 Wing-only model with rod elements positions . . . 28

22 The simplied model used for proof-of-concept testing . . . 28

23 Non-symmetric result from an optimization . . . 31

24 Same problem as in Figure 23, solved with an added symmetry constraint . . . 32

25 Optimized topology of the tetrahedron model . . . 36

26 Inuence of full wing compression constraint . . . 37

27 Inuence of compression constraint. 26 constraints used throughout the wing . 38 28 Simplied model using SPCs, no displacement constraint . . . 39

29 Simplied model using only inertia relief . . . 39

30 Simplied model, no displacement constraint, 2mm shell thickness . . . 40

31 Simplied model with displacement constraint on the whole wing . . . 40

32 Displacement constraint, only in specic points . . . 41

List of Tables

1 Optimization data, full wing constraints . . . 37

(8)

Nomenclature

Abbreviations

TO Topology Optimization FE Finite Element

SIMP Solid Isotropic Material with Penalization QUAD4 4-node Quadrilateral element

CST 3-node Constant Strain Triangle element SPC Single Point Constraints

DOF Degrees of Freedom

Mathematical Notations

Ω Design Domain

x Set of all design variables

xe Design variable (element density) for element e

ˆ

Ke Non-penalized element Stiness Matrix

K(x) Global Stiness Matrix Ke(x) Element Stiness Matrix

M (x) Global Mass Matrix u(x) Nodal displacement vector

¨

u(x) Nodal acceleration vector f (x) Nodal force vector

fe(xe) Element force vector due to gravitational loading

¯

f Nodal force vector for the inertia relief case

I Identity matrix

m Number of design variables c Number of constraints

N Number of elements in the nite element mesh o Number of compression constraints

gi Constraint functions

C(x) Compliance

Lower element density limit Vi, Ve Element volume

V0 Maximum allowed total volume

smax Maximum allowed strain in rod elements

smin Minimum allowed strain in rod elements

p SIMP penalization factor rmin Filter radius

dmin Minimum member size

dist(i, k) Distance between center of element i and k Ωk Set of elements within distance rmin of element k

ρe Element specic material density

ρ Isotropic material density α Scalar load factor

(9)

g Gravitational acceleration vector cp Pressure coecient

P Pressure where cp is evaluated

P∞ Surrounding airstream pressure

ρ∞ Surrounding airstream density

v Aeroplane speed relative surrounding airstream

α Angle of attack

(10)

1 Introduction

With rapid technological advancements keeping product development cycles as short as pos-sible is of utmost importance for companies in a competitive market. Traditional product development is an iterative process where a design is created manually based on experience and evaluated to ensure that the design satises certain design criteria. This traditional ap-proach can be very time-consuming and will often result in a design that satises the design criteria, but still has lots of room for improvement.

Using structural optimization in the design process can increase performance and reduce de-velopment time signicantly. More so, the nal design can be specically optimized for a particular objective, such as weight minimization. This is particularly useful in the automo-tive and aeronautics industry, where lower weight leads to better performance regarding e.g. lower fuel consumption, longer range, increased manoeuvrability or the ability to carry more payload.

1.1 Project Description

The project has been carried out at the section Structural Dynamics, Aeroelasticity and Stores Separation at Saab AB (2016) in Linköping during the rst half of 2016.

The goal of this master thesis has been to evaluate the possibility to use Topology Optimization (TO) on a full aeroplane model in an early stage of the design process. The assignment was, given an aerodynamic outer structure, to nd an optimized inner structure given several load cases, geometric, and manufacturing constraints and computational cost restrictions. Things to consider included studying the level of required detail in the Finite Element (FE) model, evaluation of which structural requirements the design has to satisfy, which optimization formulation, parameters and settings to consider, boundary conditions for the model and evaluation of the optimization results.

1.2 Saab AB

Saab AB (2016) was founded as Svenska Aeroplan AB in 1937 in Trollhättan, Sweden. Since then the company has branched out in several markets and provides products such as missile systems, surveillance systems, submarines and most famously known, the Gripen-system. Saab is also active in the civilian market, having previously produced civilian aircraft as the Saab 340 and Saab 2000, and now serving as a subcontractor to leading aircraft manufacturers such as Boeing and Airbus.

1.3 Structural Optimization

The eld of structural optimization can be divided into three major categories: size, shape and topology optimization. Size and shape optimization is usually performed on more mature designs and does not introduce any major design changes. For example, size optimization can be used to determine the optimal cross-section area of a beam or the radius of a hole, whereas

(11)

shape optimization can be used to determine the optimum shape of a cross section or whether a hole should be circular or elliptical.

TO is seldom performed on an already existing design. TO begins with a predened design space, denoted Ω, in which a certain amount of material can freely be placed to nd an optimized structure given specic load cases and boundary conditions. This makes TO a great tool early in the design process where it can be used to produce a conceptual design. Traditional TO is usually done by minimizing the compliance of the structure while con-straining the volume fraction allowed in the optimization. For a given load case, compliance is dened as the inverse of the stiness so minimizing the compliance maximizes the stiness of the structure.

Illustrated in Figure 1 is an example of a simple TO solution. The objective is to minimize the compliance while only being allowed to ll 40% of the original design space. The design space is subject to a force on the right-hand side while the left-hand side is xed, as seen in Figure 1a. Displayed in Figure 1b is the optimized structure obtained via TO, given the specic optimization settings and boundary conditions.

Ω

(a) Design space with boundary conditions (b) Result of the TO

Figure 1: Example of a simple TO problem: minimize compliance formulation using 40% of the design space

The minimum compliance formulation is great for nding the most ecient load paths in the design space but might lead to a design which could experience problems such as stress concentrations, a risk of buckling, dangerous eigenfrequencies and fatigue related problems. These problems can be considered in the TO but due to increased computational time and the introduction of even more diculties and uncertainties associated with each respective area, this is usually not done. Instead, size and shape optimization can be performed on the optimal design acquired from TO in order to eliminate the previously mentioned problems. This thesis will, however, only cover the TO part of the design process.

1.4 Historical Background and Previous Research

Structural TO as it is used today was introduced by Bendsøe (1989). TO is usually conducted on a component level and not on a global scale, which is to be attempted in this project. While TO with large FE-models has been carried out, no examples of TO on the scale of a full aeroplane model has been found in the studied literature.

(12)

Since its introduction TO has been used in many industrial applications (Bendøse and Sig-mund, 2003). A well known industrial example is the use of topology, size and shape opti-mization when designing conceptual wing box ribs for the Airbus A380 aircraft (Krog et al., 2002). The use of this optimization focused methodology has reportedly reduced the total weight of the aeroplane by over a thousand kilos (Krog et al., 2004).

More closely related to this thesis; research has been done by Quinn (2010) on large-scale TO problem using a combination of several weighed load cases and inertia relief to avoid rigid body modes. While the work carried out by Quinn (2010) was done on a car chassis, many similarities exist with this project. Modelling was done on a full vehicle model with critical components modelled as point masses, which basically is the same approach used in this thesis. Cavazzuti et al. (2011) also studied a similar problem without the use of point masses. Also studied by Luo et al. (2006) is the interaction between TO and inertia relief methods when performing TO on a missile body. Several load cases were considered which combined to a single optimized body. All these studies show promising behaviour when coupling inertia relief and TO.

1.5 Software

The Altair HyperWorks package (Altair, 2016a), including the pre-processor HyperMesh, the solver OptiStruct (Altair, 2016b) and the post-processor HyperView have been used through-out this project. Matlab has been used to simplify and automate certain processes.

1.6 The Aeroplane Model

The aeroplane model which was to be used during the project is a conceptual ying-wing unmanned aerial vehicle design, showcased in Figure 2. The model contains approximately 140 internal components which have to be taken into consideration in the TO.

Figure 2: Conceptual ying-wing unmanned aerial vehicle model

The conceptual design also features a number of internal cavities, seen in Figure 3, such as air intake to the engines, fuel tanks and landing gears, with some of those internal cavities accommodating internal components.

(13)

Figure 3: Overview of internal cavities

1.7 Other Considerations

No gender or environmental questions have been considered. The project has been carried out at Saab AB which main area of operation is in the defence industry, but no ethical considerations regarding this have been made.

(14)

2 Theoretical background

The following chapter serves as an introduction to the eld of TO and inertia relief, as well a short introduction to aerodynamic panel methods and aerodynamic forces.

2.1 The General Structural Optimization Problem

A general structural optimization problem (SO) can be described as in equation (2.1) Chris-tensen and Klarbring (2008). The problem consists of a objective function go(x, y), e.g. mass,

displacement etc., which is minimized. Maximization problems can also be considered, but are usually rewritten as minimization problems.

(SO)                min g0(x, y) s.t.          behavioral constraints on y design constraints on x equilibrium constraints. (2.1)

Design constraints constitute constraints on the design variables x, which could represent a minimum allowed thickness of a plate or the maximum length of a rod. The variables y are state variables which represent the response of a structure. Considering structural optimization, this response is usually mass, displacement, force, stress etc. The equilibrium constraint reads

K(x)u = f (x) (2.2)

where K(x) is the global stiness matrix, u is the nodal displacement vector and f(x) is the nodal force vector. The force vector may be design variable independent, i.e. f, in cases where no internal body forces such as gravity, centrifugal forces or inertia forces are taken into account (Zheng et al., 2009). More regarding body forces will be discussed later in section 2.7. Concerning structural optimization it is very frequent that the equilibrium constraint uniquely denes the state variable y. For this case the state variable, the displacement u, is implicitly dependent on the equilibrium constraint such that u = u(x) = K(x)−1_{f (x)}_{. The problem}

can then be formulated, on nested form, as

(SO)nf          min x∈Rmg0(x) s.t.    gi(x) ≤ ¯gi, i = 1, . . . , c ¯ xe ≤ xe≤ ¯xe, e = 1, . . . , m (2.3)

where m is the number of design variables contained in x ∈ Rm _{with c being the number of}

constraints. gi is a constraint function which must be below a dened upper limit ¯gi.

¯ xe and

¯

(15)

2.2 Finite Element Discretization

In order to nd the optimal placement of material in the design domain Ω via structural TO the domain is discretized into a FE-mesh consisting of N elements, illustrated in Figure 4. The optimization determines which elements to keep (representing material) and which elements to remove (representing void). Traditionally TO is carried out using the lowest order nite elements, i.e. 3-node CST-triangles and 4-node Quadrilaterals in two dimensions and either 6-node Tetrahedron or 8-node Hexahedron elements in three dimensions. The reason these simple elements are used, when a more accurate result probably could be an achieved using higher-order elements is simply the large increase in computational time by using higher-order elements.

Ω

(a) Design space Ω before nite element discretiza-tion

x

e

(b) FE mesh after discretization

Figure 4: FE discretization of the design space Ω with each element given a designated design variable xe

Associated with each element in the design domain is a design variable xe, which determines

whether the element should be included in the nal design. A value of 1 implies that the element should be kept in the nal design and 0 (or very close to 0, more on that later) implies that the element should not be included.

2.3 The Topology Optimization Problem Formulation

The most common TO problem formulation is the minimum compliance formulation. Com-pliance, here denoted C(x), is dened as the inverse of the stiness, as that minimizing the compliance returns a structure as sti as possible given the boundary conditions and opti-mization settings. The TO problem for minimizing the compliance can be formulated as

(TO)                min x∈RmC (x) = 1 2f T_(x)u(x) s.t.        m X e=1 Vexe ≤ V0 0 < ≤ xe≤ 1, e = 1, . . . , m. (2.4)

(16)

The design domain is subject to a volume constraint, where the currently lled volume of the design domain, calculated by the sum of all individual element volumes Ve times the design

variable xe has to be lower or equal to a predened maximum allowed volume V0. It should

be noted that this volume constraint also can be dened over the entire volume rather than just the design domain, making the constraint read PN

i=1Vixi ≤ V0. The design variables xe

are constrained to vary continuously between a small number and 1. The lower bound is introduced in order to avoid zero-stiness diagonal terms in the stiness matrix K(x), making it singular, and thus making the nite element equilibrium equation (2.2) non-solvable. Several other problem formulations can be formulated, including objective functions such as structural mass, maximum displacements (Sigmund, 1997), stress (Holmberg et al., 2013), eigenfrequency (Pedersen, 2000), fatigue life (Holmberg et al., 2014) and buckling (Gao and Ma, 2015). Constraints and objective functions are mostly interchangeable: functions that can be used as objective functions can be used as constraints and vice verse.

2.4 Solving the Optimization Problem

Due to the non-convex nature of most structural optimization problems and the diculties as-sociated, it is not feasible to solve large-scale optimization problems directly. Instead, explicit convex approximations of the general optimization problem are generated and subsequently solved instead (Christensen and Klarbring, 2008).

The optimization algorithms used to solve these subproblems are commonly of the rst order, meaning that the highest order of derivatives used is 1; implying that gradients of the objective and constraint functions have to be calculated. Common algorithms include The Method of Moving Asymptotes (MMA) (Svanberg, 1987), Convex Linearization (CONLIN) (Fleury, 1989) and Sequential Linear Programming (SLP), among others.

Solving the optimization problem is an iterative process and can be described by the following steps:

1. Start with an initial guess x0_{. Set iteration counter k = 0.}

2. Calculate the displacement vector u(xk₎_{by the nite element equilibrium equation (2.2).}

3. Calculate objective and constraint function values and respective gradients given the current design xk_.

4. Formulate and solve an explicit convex approximation of the optimization problem to achieve a new design xk+1_.

5. Update k = k + 1. Control if a convergence criteria is satised; if not, return to step 2. Convergence criteria are most commonly formulated in either of two ways: a maximal allowed change of the objective function value between two iterations or a maximum dierence in design variables between iterations.

(17)

2.4.1 Sensitivity analysis

Calculation of gradients is referred to as sensitivity analysis. Two types of methods exist to calculate the sensitivities: numerical and analytic methods. Numerical methods are based on nite dierences which may be inaccurate and computationally expensive and are therefore seldom used in TO. These methods still have their uses due to the ease of implementation and can give an estimate whether an analytic implementation has been done correctly or not, but for pure calculations, analytic methods are used almost exclusively.

Analytic methods are based on the dierentiation of the objective and constraint function with respect to the design variable xe. Due to the nature of the problems solved in structural

TO the gradients are almost exclusively solved using an adjoint method where a second set of linear equations is introduced and solved each iteration. Refer to Appendix A for more details regarding this.

2.5 Intermediate Element Densities and the SIMP-Method

As previously mentioned, a design variable xe is allowed to vary continuously between a small

number and 1. This is, however, contradictory to the fact that the nal solution is only supposed to contain elements with design variable values of either or 1, representing void or material. If a normal nite element problem is considered, the global stiness matrix K is assembled from all element stiness matrices Ke, according to Cook et al. (2007), as

K =

N

X

e=1

Ke. (2.5)

This standard formulation considers all elements in the FE-mesh to have a stiness which is not governed by the design variable. As such, even elements with a design variable value of will give a full stiness contribution to the global stiness matrix. Because of this, the stiness of each design element is said to vary linearly with the design variable xe, such that

the global stiness matrix is assembled by K(x) =

N

X

e=1

Ke(x), (2.6)

where K(x) is the global stiness matrix and Ke(x) is the (penalized) element stiness

matrix. The element stiness matrix for design elements is calculated as

Ke(x) = xeKˆe, e = 1, . . . , m, (2.7)

where ˆKe is the non-penalized stiness matrix. Note that for problems with only design

elements m = N holds.

No physical representation of the design variable has been done up until this point other than representing void and 1 representing material. While a physical representation is not necessary since the nal design only should consist of void and solid elements, this can be done by thinking of the design variable as a thickness in a two-dimensional problem or, as it will be referenced throughout this report, a density in a three-dimensional problem.

(18)

The solution of an optimization problem using the global stiness matrix formulation in (2.6) may show large areas consisting of elements with intermediate element densities, as showcased in Figure 5. While a physical representation of these intermediate design variables has been introduced, it is rarely feasible to manufacture products with a continuously varying thickness or density and thus the existence of intermediate element densities should be minimized.

Ω

(a) The design problem (b) Solution with intermediate element densities

Figure 5: Showcase of a solution with intermediate element densities, represented by grey areas in the nal design. Minimize compliance formulation with mass constraint of 50%. In order to eliminate intermediate element densities in the nal design a penalization scheme is used. The Solid Isotropic Material with Penalization (SIMP) method was rst introduced by Bendsøe (1989) and is still the most commonly used approach. Instead of the design element stiness varying linearly with the density as in (2.7) a power-law penalization scheme is used:

Ke(x) = xpeKˆe, e = 1, . . . , m. (2.8)

With a value of p > 1 (usually p = 3 (Bendøse and Sigmund, 2003)) the optimization will favor non-intermediate density elements. Intermediate densities gain an uneconomically low stiness compared to the amount of volume they occupy in the design space and therefore tend towards a density of or 1. A nal solution consisting of only void and solid elements is referred to as a black-and-white solution.

2.6 Numerical Instabilities

A number of numerical problems arise during the TO process. Sigmund and Petersson (1998) divide these problem into three categories; checkerboard patterns, mesh dependency and the existence of local minima. These problems are discussed in the following chapter.

2.6.1 Checkerboard-patterns

The appearance of checkerboard-like patterns (as shown in Figure 6) is common when the SIMP-method is used. These patches of altering void and solid elements appear due to numer-ical approximations introduced in the FE formulation which causes the checkerboard pattern to have a higher stiness than a solid block of material (Diaz and Sigmund, 1995).

(19)

Ω

(a) The design problem (b) Solution showing checkerboard behavior

Figure 6: Showcase of checkerboard pattern. Minimize compliance formulation with mass constraint of 50%.

The problem with checkerboard pattern can partially be avoided by using higher order nite elements in combination with a lower SIMP-penalization factor (Diaz and Sigmund, 1995), but the use of higher order elements drastically increase computation time. With TO runs already being computationally demanding, the use of higher order nite elements is usually avoided.

In order to circumvent these problems dierent ltering techniques are used. First introduced by Sigmund (1994) is a sensitivity lter which modies the design sensitivity of a specic element based on a weighted average of nearby elements. The method is purely heuristic but has shown to produce good results with no major increase in computation time, while being simple and thus easy to implement. Another approach is a density slope algorithm introduced by Petersson and Sigmund (1998). The algorithm enforces a certain local gradient on the slope of element densities such that element densities of adjacent elements are not allowed to vary largely.

Since many of these techniques also deal with problems associated with the detail of the FE-mesh, they will be covered in detail later on in section 2.6.2.

2.6.2 Mesh-dependency

With an increasingly ne FE-mesh more thin structural members will exist in the optimal solution. This behaviour can be seen in Figure 7 where two FE discretizations of 1800 and 20000 quadrilateral elements respectively are used to solve the same problem as in Figure 6a. As clearly seen, more detail is emerging in the solution using the ner mesh.

(20)

(a) Result using a 60×30 mesh (b) Result using a 200×100 mesh

Figure 7: Minimize compliance formulation with mass constraint of 50%. Showcasing the dierences between dierent mesh resolutions

Therefore, for an innite detailed mesh an innite number of structural members will exist. This is of course not a wanted behavior and instead, an optimal structure of the same shape should be produced regardless of the detail of the nite element model. The two ltering techniques introduced in 2.6.1 both deal with this problem. The sensitivity lter modies the sensitivity of the objective function in (2.3) according to

∂ ˆf ∂xk = (xk)−1 1 PN i=1Hi N X i=1 Hixi ∂f ∂xi , (2.9)

in which the weight factor Hi is equal to

Hi = rmin−dist(i, k), {i ∈ N | dist(i, k) ≤ rmin}.

A minimum member size dmin= rmin/2is dened such that the formation of increasingly thin

members is suppressed. The lter modied the sensitivity with respect to design variable xk

based on a weighted average of design variables associated with elements within the ltering radius rmin. Ωkdenote the set of elements inside this radius with dist(i, k) being the distance

between the center of element k and i. ∂ ˆf /∂xkis the updated sensitivity of element k, which

is used in the optimization instead of the original sensitivity ∂f/∂xk. It should be noted that

the weight operator Hi is equal to zero outside the ltering radius as it decays linearly with

the distance between the elements.

While the sensitivity lter has shown to be very ecient, numerical diculties might arise when the method is used in conjunction with multiple constraints in the optimization problem. The quality of the constraint approximation might be heavily inuenced by the ltering and as such satisfying behavioural constraints might be dicult. (Zhou et al., 2001)

Zhou et al. (2001) also suggest a new algorithm - a modication of the previous mentioned density slope method. This algorithm is implemented in OptiStruct and works by constraining the lower bound of the element density by

xi≥ max[, xj− (1 − )dist(i, j)/rmin] (2.10)

where element j is the highest density element adjacent (inside the radius rmin) to element i

at the previous iteration, as xj = max(xk | k ∈ Ωi).

(21)

Both these presented approaches for eliminating checkerboards and mesh dependency intro-duce a new problem: transition areas consisting of intermediate density elements appear at the border of created structural members. The most common way to combat this behaviour is a several-step optimization scheme. For the modied density slope algorithm, a three-step scheme is suggested:

Step I: Perform an initial optimization with density slope control according to (2.10) using a default SIMP-penalization parameter p

Step II: Once Step 1 has converged, increase the SIMP-penalty exponent to p = p + 1 and run the optimization again

Step III: The density slope constraints on the boundary between void and solid elements are relaxed so the solution can reach a black-and-white design

This method might cause slight checkerboard-pattern while solving three-dimensional prob-lems. The method is, regardless of this, used in commercial TO software with great results. 2.6.3 Local Minima

Local minima refer to the fact that the same optimization problem might produce dierent solutions depending on initial optimization parameters such as x0, optimization settings, etc.

This is due to the non-convexity introduced with the SIMP-penalization, such that minima found in the optimization are stationary points and not necessarily global minima.

This problem can partially be avoided by using so-called continuation methods, where the initial convex problem gradually is transformed to the non-convex problem. Sigmund and Torquato (1997) suggest when using a sensitivity lter, a ltering radius rmin that is started

at a large value and gradually decreased throughout the optimization. For the modied density slope algorithm, the three-step optimization scheme described earlier reduces the problem.

2.7 Inuence of Body Forces

The introduction of body forces, i.e., gravitational, centrifugal or inertial forces, introduces new diculties into the optimization problem. Since these forces are dependent on the mass of specic elements, the resulting force will change throughout the optimization (Bruyneel and Duysinx, 2005). Only gravitational forces will be studied in this project, so the modied FE-equilibrium equation reads (Lee et al., 2012):

K(x)u = f (x),

with the contributing element force fe(xe)being given by

fe(xe) = xeρeVeαg

where ρe is the element density, α is a load factor used to scale the gravitational acceleration

g. By only considering homogeneous isotropic material throughout the design domain the density ρecan be written without element index, i.e. only ρ, giving the element force equation

(22)

A lower bound on the volume is needed when only gravitational forces are considered since the optimal solution will be the removal of almost all material. With element density ap-proaching the lower bound the ratio between the body force and the SIMP-penalized sti-ness approaches innity. This will lead to a solution that converges to a design with lots of intermediate density elements (Holmberg et al., 2015). As such, a modication of the SIMP-penalization scheme has to be used. In what way this is done in OptiStruct is kept a company secret, and thus can not be covered here.

To exemplify this, consider the problem seen in Figure 8a. Acting on the body is a gravitational acceleration of 9810 mm/s2 _{in negative z-direction. The compliance is minimized with a}

maximum allowed usable volume fraction of 25% of the design space volume and minimum member size control via the modied density slope algorithm. Visualized in Figure 8b is the problem associated with a lack of lower volume constraint; the optimization simply removes the majority of the material.

(a) Problem visualization. Gravity acting in

nega-tive z-direction (b) No minimum volume constraint

Figure 8: Gravity eects on the TO result

With an added lower bound of 20% on the volume constraint the optimized design seen in Figure 9a is found. OptiStruct experiences problems converging to a black-and-white solution, most likely due to the nature of the modied density slope algorithm. If a point mass with a weight equal to 25% of the total weight of the design domain is added, the solution in Figure 9b is found. The optimization now opts to use the full allowed volume fraction of 25%.

(23)

(a) Added minimum volume constraint of 20% (b) Added point mass to the top of the design do-main

Figure 9: Gravity eects on the TO results using a minimum lower volume constraints and point masses

2.8 The Inertia Relief Method

In order to perform static FE analysis on a body, the body is required to be in (quasi)-equilibrium. Traditionally, this is done using Single Point Constraints (SPC) which restrict translational and rotational movement by applying reactional forces and moments in specic nodes, causing the body to be in static equilibrium.

Using normal SPC all externally applied forces will ultimately be reacted in the constrained nodes, with load paths between the external load and the SPC. In reality, these load paths do not always naturally exist. Rather, it is possible that in some cases the external forces often are reacted via local deformations and accelerations (Christensen et al., 2012). This makes normal SPCs not always applicable, and special techniques are needed when studying unconstrained bodies subject to external loads such as aeroplanes and satellites. Stress concentrations due to SPCs are removed when inertia relief is used, which is exemplied in Appendix B.

The Inertia Relief Method works by restricting rigid body motion by application of calculated translational and rotational accelerations (Liao, 2011). With the applied inertia loads the body is found to be in a state of quasi-equilibrium and static analysis can be performed without the creation of non-naturally existing load paths (Quinn, 2010).

To start, the inertia relief accelerations are calculated from the dynamic equilibrium equation without damping terms

M ¨u + Ku = f , (2.12)

which is an extension of the static equilibrium equation (2.2). M is the mass matrix and ¨u the nodal acceleration vector. The displacement vector u may be written as

u = ur+ ud (2.13)

where uris the rigid body displacement and udthe local deformation of the body. Combining

(2.12) and (2.13) gives

(24)

Inertia from local deformations may be neglected, so that (2.14) can be reduced to

M ¨ur+ Kud= f , (2.15)

where ¨urhas to be eliminated in order to allow static FE-analysis to be carried out. The rigid

body displacement satises Kur= 0

and is an eigenvector associated with the eigenvalue of K. By collecting the m zero-eigenvectors in an orthonormal matrix R ∈ Rn×m _{it is possible to express the rigid body}

displacements as

ur= Rc (2.16)

in which c(t) ∈ Rm_{. With R being independent of time and by combining (2.15) and (2.16):}

M R¨c + Kud= f .

Pre-multiplication with RT _gives

RTM R¨c + RTKud= RTf . (2.17)

The symmetry of K allows the second term on the left hand side of (2.17) to be re-written as

RTKud= (KR)Tud= 0, (2.18)

and (2.15) can thus be simplied as RTM R¨c = RTf ,

from which ¨c can be solved as ¨

c = (RTM R)−1RTf . (2.19)

Substituting (2.19) in (2.16) gives the rigid body accelerations according to ¨

ur= R(RTM R)−1RTf . (2.20)

Combining (2.15) and (2.20) gives the nal FE equilibrium equation as

Kud= f − M R(RTM R)−1RTf , (2.21)

which for simplicity is expressed as

Kud= ¯f (2.22)

where ¯f is the inertia relief force vector. The system still has to be constrained to remove rigid body motion, i.e by locking six degrees of freedom for the three-dimensional case. This will however, in contrast to SPCs, not introduce any stresses caused by the locked DOFs.

(25)

2.8.1 Interaction with Topology Optimization

The major dierence when considering TO with inertia relief in contrast with SPC is the inuence of mass. The mass will have no eect on a optimization assuming no body forces, but the inertia associated with the mass will impact the result signicantly when inertia relief is used.

Furthermore, body forces are not compatible with inertia relief analysis. The (for example) gravitational accelerations are simply nullied by the inertia relief acceleration acting in oppo-site direction. As such, a model analyzed using inertia relief under only gravitational loading will not be deformed or be subject to any stresses, other than what is caused by numerical errors. Mathematically, the sensitives of the displacement and force vectors are dierent using inertia relief; Pagaldipti and Shyy (2004) and the MSC Nastran (2012) Design Sensitivity and Optimization User's Guide cover this more thoroughly and it will as such not be discussed further here.

To visualize the eects of inertia relief in TO lets study an example. As a reference, consider the problem in Figure 10a with a load acting on top of the design space as well as xed supports in each corner, with the corresponding optimized topology in Figure 10b. Optimization settings include a compliance minimization with a volume fraction constraint of 25% and minimum member size control via the modied density slope algorithm.

(a) Standard problem with SPCs (b) Solution of 10a

Figure 10: TO problem and solution using normal SPCs

Figure 11a shows the solution of a slight modication of the problem in Figure 10a; the SPC in Figure 10a are removed and replaced with inertia relief. In this case there are no load paths to the bottom corners of the model and as such no material is distributed towards to corners. If instead masses are added to the bottom corner nodes the program can utilize the added inertia by connecting the masses to the point where the external load is applied. The result of this is seen in Figure 11b.

(26)

(a) Automatic inertia relief control (b) Automatic inertia relief control with point masses added to each lower corner

Figure 11: Showcasing the eect of inertia relief on the TO result, both with and without added point masses

2.9 Short on Aerodynamic Panel Methods and Aerodynamic Forces

Aerodynamic potential-ow methods, also called Panel methods, are methods to evaluate uid velocity around an object. Knowing the uid velocity pressure coecients can be calculated and subsequently the pressure distribution. No eort will be made throughout this report to describe the theory behind these methods, but the interested reader can consult an introduc-tion made by Erickson (1990).

However, some aerodynamic relations which are used are worth mentioning. The relationship between pressure coecient cp and the surrounding pressure can be described by

cp =

P − P∞

q , (2.23)

where P is the pressure at the point where the pressure coecient is being evaluated, P∞ is

the far-away pressure and q is the dynamic pressure calculated by q = 1

2ρ∞v

2_, _(2.24)

given the density of the surrounding uid ρ∞ and the velocity v of the object relative to the

uid.

Given specic ight conditions (such as cruise) the aeroplane is supposed to be in static equilibrium. In this state forces and moments acting on the plane are zero with rudders being used to balance the moments. The dominating moment, and the only moment taken into account here is the pitch moment. Other ying-wing aeroplanes use a pitching so-called "beaver-tail" to balance this moment, but for simplicity the moment will be countered by modelling aperons (a combination of aps and ailerons) at the trailing edge of the wing. To nd the angle of attack α needed to achieve the required lift force a sweep is done over the angle of attack. For example, the aerodynamic problem is solved for an incremental increase

(27)

of angle of attack from 0◦ _{to 10}◦_{. From this the angle of attack needed to achieve the required}

lift force can be determined, and subsequently, the pressure distribution given the specic angle of attack.

(28)

3 Method

The following chapter describes a suggested work ow when considering structural TO of inter-nal aeroplane structures. The suggested approach utilizes an aerodynamic solver to calculate pressure distributions over the wing given specic ight manoeuvres. An interpolation scheme is then created based on the pressure distribution such that the pressure can be mapped to dierent FE-meshes easily. With the pressure mapped to the nite element mesh, classical TO can be carried out.

3.1 Geometry Clean-up

Before the model can be meshed it has to be simplied so that the model can be meshed in an ecient way. It is desirable to have a mesh that is as homogeneous as possible, which an original model might not produce. Problematic areas include duplicate geometry lines or geometry lines very close to each other (Figure 12a), thin edges (Figure 12b), intersecting surfaces (Figure 12c) or unnecessarily complex geometries (Figure 13).

(a) Duplicate lines (b) Thin edges (c) Intersecting surfaces

Figure 12: Visualization of common areas that needs to be xed before meshing Geometry lines in HyperMesh determine the way the FE mesh is created. The FE mesh will follow all geometry lines but lines can be suppressed such that they are not taken into consideration when a mesh is created. Suppression of lines should not necessarily be restricted to suppressing duplicate or adjacent lines as an excess of unnecessary lines will make it harder for HyperMesh to create a good mesh.

Other problematic areas require more complex ways to x and are therefore more easily dealt with by modifying the model in a CAD-software instead of in HyperMesh. It is possible to simplify the model by clever use of functions found in HyperMesh, but it would be preferred to have a detailed ready-to-mesh model directly from the CAD-software.

Due to the absence of CAD-software, all geometry clean-up throughout the project has been carried out in HyperMesh. The thin edge seen in Figure 12b have been given a thickness and the intersecting surfaces in Figure 12c have been removed by slightly moving surfaces around. Also modied is the complex nose cone seen in Figure 13, where the "bubbles" have been replaced with a at surface.

(29)

(a) Nose before geometry clean-up (b) Nose after geometry clean-up

Figure 13: Close up shot of the geometry clean-up done on the aeroplane nose. It is now possible to mesh the nose without incorporating lots of ill-conditioned elements

3.2 Meshing

While it is desired to have a homogeneous element size throughout mesh it is not always possible. By geometry clean-up most areas can be modied such that a homogeneous mesh can be created, but modifying the geometry too much might result in information being lost and thus give a misleading result. Therefore, dierent mesh sizes can be necessary. For this model, this is done along the engine exhaust as seen in Figure 14a. Also note the added thickness along the edge, which was added during the geometry clean-up.

Also, it might not be necessary to include the whole model in the design domain. Areas where a very ne mesh is needed to keep element quality high are preferably removed from the design domain. This is done along the trailing edge of the aeroplane model, as seen in Figure 14b, where the area in red is removed from the design domain and thus not included in the mesh. If previous optimization runs have shown large areas where no material is placed, the empty areas can also be removed from the design domain beforehand to reduce computational time.

(a) Dierent mesh size around the engine exhaust (b) Trailing edge (in red) not being meshed

Figure 14: Dierent aspects of mesh creation 3.2.1 Meshing Approach

The general approach to meshing the FE-models has been with an outer layer of shell elements and an internal solid mesh. The outer shell layers, consisting of QUAD4 and CST elements,

(30)

serves as non-design space and creates the outer shape of the aeroplane. It is on this shell mesh that aerodynamic pressures will be applied. The internal space of the aeroplane is meshed with tetrahedron elements. The shell and solid layers share FE-nodes, so forces applied on the outer shell are transferred to the solid mesh.

Additionally, to help HyperMesh create a good tetrahedron mesh internal cavities can also be meshed with shell elements. The shell mesh is used to guide the creation of the internal tetrahedron mesh so specically in areas with complex geometries this is very favourable to do. Also, the tetrahedron mesh can be optimized for factors such as element size, aspect ratio and similar to give a better result.

It is possible that the optimization is heavily governed by the thickness of the outer shell layer. If the outer shell layer is too thin, no load will be taken by the shell layer, and if it is too thick, all load will be taken by the outer shell. How much load the shell takes directly translates to how much load the internal structure has to absorb, which directly will aect the optimization result.

3.3 Component Modelling

Internal components such as engines, landing gears etc. are modelled as point masses con-nected to the FE-mesh. In cases where the components are far away from the mesh, rigid elements are created to connect the mass to the FE-mesh. The program described in section 6.1 is used to speed up this process.

Flaperons are modelled with rigid elements connecting to a node in which forces are applied. The point where the aps are localized on the aeroplane has not yet been specied in the design process and therefore their position in the FE-mesh is somewhat ctional.

While in a real scenario the internal components would occupy volume, this is mostly not considered as it would make the model overly complicated as it would introduce more uncer-tainties and restrictions on the design space. For example, which attachment points to use would have to be evaluated for each component. Some major components, such as the engines with respective air intakes and landing gears have predened cavities.

3.4 Flight Manoeuvres and Load Cases

Several load cases have to be considered in the optimization to achieve a robust design. If only a single load case is considered, a design that is well suited for the specic load case will be found, but the design might perform poorly in other load cases.

The aerodynamic forces due to the dierent manoeuvres are calculated using an in-house Matlab-code at Saab. In the program, an approximation of the ying wing geometry is modelled, and aerodynamic panel methods are used to calculate output data such as pressure coecient distribution, lift force and pitch moment. This pressure is then mapped to the FE-model via a Matlab-script, described in section 6.2.

Two load cases are used during the optimization. These are obtained from a steep pitch-up manoeuvre scenario and a landing scenario, both of which are described in the following

(31)

sections. The reasoning behind choosing those two load cases is pretty simple; the aeroplane is simply not designed to perform any advanced manoeuvres and because of this the landing and pitch-up manoeuvre load cases are among the most demanding. Further, the two load cases will cause wing bending in dierent directions, upwards during the pitch-up manoeuvre and downwards during the landing.

3.4.1 Pitch-up Manoeuvre

During the pitch-up manoeuvre load case only aerodynamic forces are considered due to the incompatibility between the inertia relief analysis and the gravitational load. The pressure distribution over the aeroplane can be seen in Figures 15 and 16.

Components in this load case are modelled as forces. The forces are supposed to act in a global z-direction directly towards the centre of the earth to simulate gravitational pull. The aeroplane, however, is not perpendicular to the ground as it ies with an angle of attack to increase lift, and therefore the forces are angled to be perpendicular to the ground. Flight data used for the pitch-up manoeuvre case is a ight height of 3 km and a speed of 400 km/h.

Figure 15: Pressure distribution (in [Pa]) on the upper part of the body

(32)

3.4.2 Landing

Landing simulations are normally time-dependent problems which are far to complicated to be considered in the TO. To simulate landing, increased gravitational loads are used, i.e. the load factor α in (2.11) is larger than 1. Further, SPCs are used with no aerodynamic loads being considered.

3.5 Boundary Conditions

Boundary conditions are set up for the two load cases in either of two ways: Landing, normal SPCs: Constrained xyz-directions in landing gear positions

Pitch-up, Inertia Relief: Constrained nodes along the symmetry plane of the aeroplane to avoid introducing non-symmetric loading

The SPC locations are chosen based on the location of components associated with the landing gears, and can be seen in Figure 17. While in a real landing the plane usually lands with the two back landing gears touching the ground rst it is assumed for simplicity that all landing gears touch the ground at the same time. Further, in a real scenario the plane will obviously roll forward while being able to move sideways, thus only really being constrained in the z-direction. For simplicity this is not taken into consideration either and all DOF are constrained.

Figure 17: SPC locations

For the inertia relief method to work all six degrees of freedom of the model have to be constrained. It is also important not to introduce any non-symmetric boundary conditions which would introduce any non-realistic loads into the structure. Because of this, the inertia relief reference points are chosen along the aeroplane symmetry plane as seen in Figure 18.

(33)

Figure 18: Constrained nodes for the inertia relief analysis

Constraint number 1 in Figure 18 represents a constraint in x-, y- and z-directions, number 2 is a constraint in y- and z-directions and number 3 in only y. This set-up of constraints will remove all three translational DOF and all three rotational DOF of the system without introducing any unwanted stresses or unsymmetrical behaviour.

3.6 The Optimization Problem, Constraints and Settings

The compliance minimization problem (2.4) with a volume fraction constraint is slightly mod-ied to account for the two load cases, the pitch-up manoeuvre and landing load cases. The modied optimization problem reads

(TO)                min x∈RmC (x) = 1 2 w1f T 1 ud,1(x) + w2f2T(x)u2(x) s.t.        m X e=1 Vexe ≤ V0 0 < ≤ xe≤ 1, e = 1, . . . , m, (3.1)

where w1 and w2 are weight factors, u1 solves the inertia relief equilibrium equation

K1(x)ud,1(x) = ¯f1,

and u2 solves the static FE equilibrium equation with body forces

K2(x)u2(x) = f2(x).

The index on the stiness matrices, the force and displacement vectors represent the two load cases, with 1 representing the pitch-up manoeuvre case and 2 the landing load case.

(34)

Constraints to reduce compression of the wing has also been tested, making the optimization problem read (TO)                        min x∈RmC (x) = 1 2 w1f T 1 ud,1(x) + w2f2T(x)u2(x) s.t.                m X e=1 Vexe ≤ V0 0 < ≤ xe≤ 1, e = 1, . . . , m smin< so(x) < smax, o = 1, . . . , n, (3.2)

in which so(x) is the strain in rod element o, smin and smax are bounds on minimum and

maximum allowed strain and n is the number of strain constraints. The reason why the rod elements are needed is discussed in 5.1 and how the constraints are implemented and applied is covered more thoroughly in 6.3. Symmetry constraints on the design variable are also used along the symmetry plane of the aeroplane to achieve a symmetric design, more on this in 5.2. Apart from objective function and constraints, there is an incredible amount of optimization settings in OptiStruct and it has not been feasible to try out every combination available. The settings most commonly used are:

Minimum Member Size dmin: Three times the average element size

Penalization factor p: 3. The variable in OptiStruct, DISCRETE, corresponds to p − 1 Initial density fraction: Equal to the allowed mass fraction

Minimum element density : 0.001

Optimization convergence tolerance: 0.005

Max allowed iterations: 200, which should never be reached assuming standard compli-ance minimization

3.7 Result Interpretation

Optimal designs acquired via TO are often complex and need to be manually rened. The designs may contain geometries which are not manufacturable or geometries which would be too expensive to produce. Also, some creativity from the designer might be needed when interpreting the TO result.

An acquired design from TO will, in cases considering compliance minimization, be a rather non-mature design regarding e.g. stress evaluation and buckling. The design will have to undergo testing to make sure that it satises all design requirements, either via other forms of optimization (size or shape) or by static analyses.

When interpreting the optimization results in OptiStruct, it is possible to get dierent levels of details depending on which element density threshold is used to illustrate which elements to keep in a nal design, as shown in Figure 19.

(35)

(a) Element densities >0.5 shown (b) Element densities >0.25 shown

(c) Element densities >0.1 shown (d) Element densities >0.01 shown

Figure 19: Dierent topologies depending on density threshold

Altair claim this is a feature in the software, allowing the user to 'choose his design'. Due to the penalization elements with a low density will contribute with almost zero stiness to the nal design, so it is uncertain why these elements are kept. The problem with kept intermediate element densities is further discussed in 5.4.

(36)

4 Models

A number of simplied FE models have mostly been used throughout the project. The models have been used to test how the optimization behaves with dierent settings, including, but not limited to: inertia relief, manufacturing constraints, dierent load cases and boundary conditions, and optimization settings.

4.1 Wing Only

The main use of the wing model, seen in Figure 20, has been testing of compression constraints on the wing. The wing model consists of an outer shell layer mesh with 2350 elements, an internal tetrahedron mesh of 9011 elements and 1166 rod elements connecting the top and bottom shell layer.

Figure 20: Wing-only model

The wing is subject to a pressure acting on the outer shell layer taken from the pitch-up load case. The pressure should, therefore, be the same as the pressure which would be acting on the wing had it been attached to the aeroplane. Further, components in the wing are also the same as those that should be included in the full model. To simulate the wing connection to the aeroplane the root of the wing is constrained with SPCs. Gravitational loads are also applied to account for the landing load case. In short, the wing is subject to the same forces as it would have been when connected to the aeroplane.

In some optimization runs, several or the majority of the rod elements are removed to give a representation of how the displacement constraints associated with these rod elements aect computational times. Some rod elements are kept along the thick line shown in Figure 21 with remaining elements being equally distributed along the thick lines to get compression constraints evenly distributed along the wing.

(37)

Figure 21: Wing-only model, rod element positions along the thick lines and corresponding compression constraints

4.2 Simplied Aeroplane, Outer Shell Mesh

The simplied aeroplane model in Figure 22 has mostly been used throughout the project. The idea is that this model can be used for proof-of-concept testing and when done, the mesh density can be increased to achieve a better, more detailed result. The model features a mesh that gives a good trade-o between accuracy and run times, making it great for testing dierent optimization settings.

The mesh consists of an outer shell layer of 12616 combined CST and QUAD4 elements and internal tetrahedron mesh consisting of 73616 elements. Due to meshing diculties the frontal nose and the tail section are not included. These areas include structures that are far to detailed to model using a rough FE-mesh. Including these areas into the FE-mesh would require ner elements and therefore the-the number of elements would be increased by tens of thousands.

Figure 22: The simplied model used for proof-of-concept testing

The outer shell layer is subject to a pressure load calculated for the take-o-loading case. The pressure loading is balanced by inertia relief with gravitational loads being considered in a separate load case balanced by SPC. Components are modelled as point masses and thereby taken into consideration via the gravitational load.

(38)

4.3 Simplied Aeroplane, Tetra-only Mesh

An even simpler model was used to test the inuence of the modelled components. The model only consists of an internal tetrahedron mesh and does not contain any non-design elements where pressure loading can not be applied in a good way and thus only the loading due to internal components is considered. The model uses xed SPC in the positions of the landing gear, located at the bottom of the plane.

(39)

5 Complications

Some complications have been reoccurring during the project. These include problems re-lated to poor material placement in the outer wing of the aeroplane, i.e. no safety against wing compression, non-symmetric results, neglected component loads, intermediate element densities and poor convergence. These problems are introduced in the following sections.

5.1 Wing Compression

A major problem throughout the project has been the lack of material in the outer part of the wing. This is troublesome as one of the main tasks of the airframe of the wing is to cancel buckling of the outer skin and due to the absence of material in the outer wing, it is safe to assume that the structure will buckle.

The reason for the lack of material in the outer wing is thought to be due to small compression of the wing relative to the bending of the whole wing. As such, the optimization favours reducing the amount of wing ex rather than limiting the compression of the wing. While mathematically correct due to the global compliance dependency of the displacements, the behaviour is not wanted.

Ideally, buckling could be taken into consideration in the optimization problem, but current versions of OptiStruct does not support this for meshes containing solid elements. Instead, the amount of compression of the wing, i.e. how much the upper shell layer and the bottom shell layer are allowed to move towards each other, can be measured and constrained.

These constraints are implemented in problem 3.2 by creating one-dimensional rod elements connecting nodes on the top shell layer to the bottom shell layer of the wing. The strain of these elements can be used as a design response which in turn can be constrained. For example, constraining the strain of the rod elements to 1% allows the wing to compress 1%. The rod elements do not introduce any additional DOF to the system so the FE-equilibrium equation should take approximately the same time to solve as before. The assembly process of the stiness matrix will become slightly longer, but the increase in time should be negligible. However, this method introduces a lot of constraints to the optimization problem. These constraints have to be evaluated, including evaluation of the sensitivities for each constraint, which will result in a massive increase of computational time.

The easiest way to avoid this massive increase in computational time is to introduce compres-sion constraints at specic points of the wing, instead of the whole wing. While this reduces the time it takes to solve the problem, the constraints are very local, and several questions have to be considered:

Where to place these constraints?

How many constraints should be included?

What distance is required between the constraints? What percentage of compression is allowed?

(40)

Should all compression constraints be of the same magnitude, i.e. same values of smin

and smax throughout the wing?

How does the number of constraints needed correspond to the detail in the FE-mesh? It has to be noted that the optimization algorithms used to solve TO problems are specically designed to solve problems with a large amount of design variables and few constraints. The increase in computational time is therefore expected, and not much can be done about it.

5.2 Non-symmetric Results

It has been seen throughout the project that a symmetric model and symmetric load might not always produce a symmetric optimized design. The dierence in those cases is often not major and can often be neglected, but if a fully symmetric problem is solved a fully symmetric solution is probably wanted. These small changes in the topology may be caused by numerical errors in the FE-calculations which may grow as the optimization proceeds.

To counteract this symmetry planes can be dened on the design variables in OptiStruct. By forcing the optimization to balance element densities on both sides of the symmetry plane, an optimal, symmetric, design is found.

Symmetry constraints can also be used to reduce the complexity of the optimization problem by reducing the number of constraints and design variables. Figure 23 and 24 show an example of this. The model features constraints on the wing compression on the right wing seen in the gures. The application of this constraint on one wing obviously causes a non-symmetric design, as Figure 23 shows.

Figure 23: Non-symmetric result from an optimization

However, if the symmetry constraint is added to the aeroplane symmetry plane the topology in Figure 24 is acquired. So, instead of introducing tens or hundreds of additional constraints, the same eect can be achieved by symmetry constraints.

(41)

Figure 24: Same problem as in Figure 23, solved with an added symmetry constraint It would be, given that the load cases used are symmetric, possible to only include half the aeroplane in the FE-mesh. While the load cases used throughout the project are symmetric, other real load cases are likely to be anti-symmetric, such as the loads caused by a turn.

5.3 Neglected Component Loads

Another recurring problem has been components "oating" in the design space, without any material being distributed towards them to connect them to the shell layer. There are some possible causes to this:

The loads associated with the components are small relative to the other loads causing the optimization to focus on countering the pressure or gravitational loading

The loads generated by the components are so small that they are carried by the -density elements

The rough mesh used does not allow the optimization to create small enough structural members to carry the component load

It might be possible to avoid this problem by scaling up the forces due to the components to the same magnitude as the external pressure, but it would not be a good representation of reality. While not being optimal, it might be necessary to do so to get material distribution to the components. A better approach would probably be a rened FE-mesh but with increased run times as a result.

5.4 Intermediate Element Densities

While intermediate element densities are to be expected at the boundary layer between solid and void elements due to the density slope algorithm, the optimization often favours large areas of intermediate element densities and does not converge to a black-and-white solution. These elements will, after penalization, provide close to zero stiness to the model.

(42)

The reason for this behaviour is not clear. It might be caused by a conservative mass fraction constraint leaving the optimization with too little material to account for all loading in the model. Simple tests have however shown this not be the case. Another consideration could be a penalization factor that is too low, but tests have also been conducted which show no signicant dierence in the number of intermediate density elements with an increased penalization factor.

Another cure to the problem could be a stricter convergence tolerance. This has shown to have some promising behaviour, but with increased computational time as a side eect. The increase of black-and-white structure is as best minimal whereas the amount of iterations usually is at least double, so overall this approach is not eective either. The upper limit on number of iterations might also have to be increased to reach convergence with the stricter tolerance.

5.5 Poor Convergence

Problems related to poor optimization convergence or failure to converge are, in this case, mostly connected to the use of compression constraints on the wing. The cause of this could be two-fold: too restrictive constraints will make it impossible for the optimization to nd a feasible design and OptiStruct will run until a maximum number of iterations is reached without converging. Changing the value of the constraint could be a simple x if it is allowed from a design point of view.

The other reason is related to a large number of constraints in the optimization problem. Op-tiStruct only evaluates a maximum number of constrains in each iteration. Which constraints that are to be evaluated is based on how much they violate the constraint or how close they are to being active. Constraints which violate the boundary the most will be evaluated rst, with constraints being the closest to the constraint boundary will be evaluated second-hand. By only evaluating a set amount of constraints each iteration while having to satisfy a large number of constraints, the optimization is forced to do a larger amount of iterations to evaluate and satisfy all constraints.

Topology Optimization as a Conceptual Tool for Designing New Airframes

Topology Optimization as a Conceptual Tool

for Designing New Airframes

Joakim Torstensson

Topology Optimization as a Conceptual Tool

for Designing New Airframes

Joakim Torstensson

Contents

List of Figures

List of Tables

Nomenclature

1 Introduction

Ω

2 Theoretical background

Ω

x

Ω

Ω

3 Method

4 Models

5 Complications