8 Optimization of continuum structures

In the previous chapter, we learned how to transform a continuous elastic problem defined on a two-dimensional domain into a discrete problem using finite elements. Now we want to optimize such problems.

8.1 Size optimization

Let’s consider a two-dimensional domain discretized with finite elements, in which we can adjust the thickness of each element. We may seek the best distribution of thicknesses \(\mathbf{d}\) to minimize compliance \(C\) with a constrained amount of material volume \(V_0\). This problem can be formulated in analogy to the standard problem of sizing optimization stated in Equation [eq:size_optimization]

\[\begin{aligned} \min_{\mathbf{d}} \quad & C(\mathbf{d}) = \mathbf{f} \cdot \mathbf{u}(\mathbf{d})\\ \textrm{s.t.} \quad & \mathbf{d} \cdot \mathbf{A} - V_0 \le 0 \\ & \mathbf{d} \in \mathcal{D}\\ \end{aligned} \label{eq:sheet_optimization}\] where \(\mathbf{A}\) denotes the area of an element. Mathematically, this problem is equivalent to the truss problem and can be solved in the same way: We formulate the Lagrangian \[\mathcal{L}(\mathbf{d}, \mu) = C(\mathbf{d}) + \mu \left( \mathbf{d} \cdot \mathbf{A} - V_0 \right) \label{eq:lagrangian_sheet_optimization}\] and determine its derivative \[\frac{\partial \mathcal{L} (\mathbf{d}, \mu)}{\partial d_j} = \frac{\partial C(\mathbf{d})}{\partial d_j} + \mu A_j\] with \[\frac{\partial C}{\partial d_j} = - \mathbf{u}_j(\mathbf{d}) \cdot \mathbf{k}^0_j \cdot \mathbf{u}_j(\mathbf{d}) = -2w_j(\mathbf{d}). \label{eq:sensitivity_sheet}\] In comparison to Equation [eq:compliance_sensitivity] for trusses with four degrees of freedom, the element displacement vector contains eight degrees of freedom for these 2D problems with linear quad elements (\(\mathbf{u}_j \in \mathcal{R}^8, \mathbf{k}^0_j \in \mathcal{R}^{8\times 8}\)).

Just like the truss problem, the variable thickness sheet problem can be approximated using MMA with lower asymptotes only. Subsequently, it is solved with the dual method and a line search to find the Lagrange parameter \(\mu\). The entire procedure is identical to Section [sec:truss_topology], if we replace the truss cross sections \(\mathbf{a}\) with element thicknesses \(\mathbf{d}\) and the truss lengths \(\mathbf{l}\) with element areas \(\mathbf{A}\).

8.2 Topology optimization with MMA

In a topology optimization of continuum structures, we want to optimize the neighborhood relations of a structure. While the size or shape of an object could be changed by some deformation, we cannot change a topology to another by a "rubber-like" deformation. Essentially this means, in topology optimization we are seeking the fundamental structure of a component: Which parts are connected? How many holes are in it?

There are two main approaches to optimize the topology (Eschenauer and Olhoff 2001): In a geometry-based formulation, we describe the domain in a dynamic fashion ( e.g. by a Level-Set method) and modify this geometry towards the optimum (e.g. by empirical growth rules like SKO or CAIO). In a material-based formulation, we describe the entire domain with finite elements and then seek a material distribution where each element is either completely filled or completely unfilled with material. In this lecture, we focus entirely on the latter approach.

We can formulate such a problem by reinterpreting the thickness variable as a normalized material density \(\pmb{\rho}\), where \(\rho_j=1\) indicates presence of material and \(\rho_j=0\) indicates absence of material in element \(j\). The stiffness of each element becomes a function of \(\rho_j\), e.g. \[\mathbb{C}(\rho_j)= \begin{cases} \mathbb{C} & \text{if} \quad \rho_j = 1 \\ 0 & \text{if} \quad \rho_j = 0 \end{cases}\]

Then, the problem statement reads \[\begin{aligned} \min_{\pmb{\rho}} \quad & C(\pmb{\rho}) = \mathbf{f} \cdot \mathbf{u}(\pmb{\rho})\\ \textrm{s.t.} \quad & \pmb{\rho} \cdot \mathbf{V} - V_0 \le 0 \\ & \rho_j \in \{0, 1\}\\ \end{aligned} \label{eq:topology_optimization}\]

where the only change compared to Equation [eq:sheet_optimization] is the discrete nature of the design variables \(\pmb{\rho}\) and the name of the design variables.

As mentioned during truss optimization, the binary problem is notoriously hard to solve. Hence we use Solid Isotropic Material with Penalization (SIMP) again, which is denoted as \[\mathbb{C}(\rho_j)= \rho_j^p \mathbb{C}\] for the generalized FEM problem.

Incorporating the SIMP method to the optimization is straight-forward: The element stiffness for a 2D problem with element thicknesses \(d_j\) is now \[\mathbf{k}_j(\rho_j) = \rho_j^p d_j \mathbf{k}^0_j\] and consequently, the gradient becomes \[\frac{\partial\mathbf{k}_j(\rho_j)}{\partial \rho_j} = p \rho_j^{p-1} d_j\mathbf{k}^0_j.\] Then, we just need to adjust the sensitivity (see Equation [eq:sensitivity_sheet]) slightly towards \[\frac{\partial C}{\partial \rho_j} (\pmb{\rho}) = - p \rho_j^{p-1} d_j\mathbf{u}_j(\pmb{\rho}) \cdot \mathbf{k}^0_j \cdot \mathbf{u}_j(\pmb{\rho}) = -2 p \rho_j^{p-1} d_j w_j(\pmb{\rho}). \label{eq:sensitivity_topology}\]

Obviously, the SIMP approach helps us to find a binary solution of the discretized problem stated in Equation [eq:topology_optimization]. However, we may ask ourselves if there is a unique solution to this problem independent of discretization. And it turns out, the answer is no: You can improve any given design by replacing it with finer structures in a process that goes on indefinitely (Christensen and Klarbring 2008). In addition to the theoretical non-existence of a solution, this also means that our solution is mesh-dependent: If we want to refine the solution, we may end up with a totally different solution.

Even if we were to accept mesh dependence and the theoretical problem of ill-posedness, we are still left with challenges in the resulting designs. First of all, fine structures may be hard to manufacture. Even additive manufacturing methods are limited to some minimal structure size. Secondly, we may observe checkerboard-like patterns on the structure. These are numerical artifacts from the fact that we use linear shape functions which lead to an overestimation of the stiffness for that configuration.

We can address all these problems with the introduction of filters as a regularization (Harzheim 2014; Lazarov and Sigmund 2011). There are two traditional approaches for mesh-independent filtering: density filtering and sensitivity filtering. In density filtering, we replace the density \(\pmb{\rho}\) with a weighted sum of neighboring densities, and use this filtered density \(\tilde{\pmb{\rho}}\) in the stiffness evaluation \[\mathbb{C}(\rho_j)= \tilde{\rho}_j^p \mathbb{C}.\] The filtered density is computed according to \[\tilde{\rho}_j (\pmb{\rho}) = \frac{\sum_i H_{ji} A_i \rho_i}{\sum_i H_{ji} A_i}\] with a distance weighting matrix \[H_{ji} = \begin{cases} R-\textrm{dist}(i,j) & \text{if} \quad \textrm{dist}(i,j) \le R \\ 0 & \text{else}, \end{cases}\] where \(R\) describes the filter radius and \(\textrm{dist}(i,j)\) the distance between element centers of elements \(i\) and \(j\). This filter results in a structural weakening of structures finer than \(R\), as they get smeared with neighboring empty elements. Hence, the filter solves the mesh-dependency problem, prevents structures that cannot be manufactured and prevents checkerboard patterns. However, we need to account for the filter during sensitivity computation by \[\frac{\partial C}{\partial \rho_j} = \frac{\partial C}{\partial \tilde{\rho}_m} \frac{\partial \tilde{\rho}_m}{\partial \rho_j} = \frac{\partial C}{\partial \tilde{\rho}_m} \frac{H_{jm} A_m }{\sum_i H_{ji} A_i}\] and in the volume constraint.

An alternative to filtering densities is filtering of sensitivities. We may compute the filtered sensitivity as \[\widetilde{\frac{\partial C}{\partial \rho_j}} = \frac{A_j}{\rho_j} \frac{\sum_i H_{ji} \frac{\rho_i}{A_i} \frac{\partial C}{\partial \rho_i}}{\sum_i H_{ji}}\] which is a purely heuristic concept, i.e. there is no mathematical proof for this filter to work (Sigmund and Petersson 1998; Sigmund 2007). However, we can implement this formulation simply by replacing the sensitivities in [eq:sensitivity_topology] with filtered sensitivities. This is very efficient, simple to implement and experience shows that this filter works quite well.

There is one more challenge in the SIMP formulation of the standard topology optimization problem: For \(p>1\), the problem is not convex anymore and different starting points may result in different local minima (Christensen and Klarbring 2008). A strategy to cure termination at non-global minima is a gradual increase of \(p\) starting from \(p=1\) or using multiple starting points.

8.3 Topology optimization with optimality criteria

The solution via MMA is in line with the previous chapters of this manuscript and is able to solve arbitrary problems beyond the compliance minimization stated in problem [eq:topology_optimization]. Historically, this specific problem was first solved with the optimality criteria method instead of convex approximations. The idea here is that we can use the gradient of the Lagrangian \[\mathcal{L}(\pmb{\rho}, \mu) = C(\pmb{\rho}) + \mu \left( \pmb{\rho} \cdot \mathbf{V} - V_0 \right),\] which is \[\frac{\partial \mathcal{L} (\pmb{\rho}, \mu)}{\partial\rho_j} = \frac{\partial C(\pmb{\rho})}{\partial \rho_j} + \mu V_j = 0,\] at the optimal point, to formulate an optimality condition \[G_j = \frac{-\frac{\partial C(\pmb{\rho})}{\partial \rho_j}}{\mu V_j} = \frac{p \rho_j^{p-1} d_j \mathbf{u}_j(\pmb{\rho}) \cdot \mathbf{k}^0_j \cdot \mathbf{u}_j(\pmb{\rho})}{\mu V_j} = 1\] that must be fulfilled at the optimum whenever a variable does not reach its bounds. A heuristic algorithm can now compute \(G_j\) and adjust \(\rho_j\) such that \(G_j\) approaches 1. We may realize that increasing \(\rho_j\) increases the stiffness and hence decreases the displacement for a given load. Consequently, increasing \(\rho_j\) decreases \(G_j\) and vice versa. Hence, an update rule \[\hat{\rho}_j^{k+1} = \left(G_j^k\right)^\xi \rho_j^k\] drives the solution towards the optimal point \(\rho_j^*\), where \(G^*_j=1\). The exponent \(\xi \in (0,1]\) stabilizes the algorithm with typical values being \(\xi=0.5\) (Harzheim 2014).

In addition, we employ move limits \[\begin{aligned} \rho_j^{-,k} &= \max \left(\rho_j^-, (1-m)\rho_j^k \right) \\ \rho_j^{+,k} &= \min \left(\rho_j^+, (1+m)\rho_j^k \right) \end{aligned}\] to compute the next value as \[\rho_j^{k+1} = \textrm{clamp} \left(\hat{\rho}_j^{k+1}, \rho_j^{-,k}, \rho_j^{+,k} \right).\] The move limits reduce the maximum step width and account for bounds on the design variables.

To compute \(G_j\), we need to determine the value of the Lagrange multiplier \(\mu\). For the compliance problem, we can intuitively assume that the stiffest structure will use all material and consequently, that the volume constraint will be active. The Lagrange multiplier may be found with the bisection method on the feasibility condition \[\frac{\partial \mathcal{L}}{\partial \mu} = \pmb{\rho}(\mu) \cdot \mathbf{V} - V_0 = 0.\]

8.4 Shape optimization

Shape optimization of continua seeks the optimal shape of the outer perimeter of a structure without adjusting the topology, i.e. without adding new perimeters like holes. The shape optimization problem for minimum compliance for a given volume constraint can be formulated analogously to trusses (see Equation [eq:shape_optimization]) as

\[\begin{aligned} \min_{\mathbf{x}} \quad & C(\mathbf{x}) = \mathbf{f} \cdot \mathbf{u}(\mathbf{x})\\ \textrm{s.t.} \quad & g(\mathbf{x}) = \mathbf{d} \cdot \mathbf{A}(\mathbf{x}) - V_0 \le 0 \\ & \mathbf{x} \in \mathcal{X} \end{aligned} \label{eq:shape_optimization_continuum}\]

A natural choice for design variables are the coordinates of nodes bounding the meshed domain, just as it has been done for trusses. However, this leads to jagged meshes with inaccurate stress calculations and limited deformation capabilities, as illustrated in the following example.

We may prevent the disadvantages of the direct nodal optimization by describing the deformation of the domain by a linear combination of basic shape modifications. This could be achieved by describing the perimeter with a spline representation (Christensen and Klarbring 2008) or by morphing operations on the mesh. Here, we decide to morph the mesh using radial basis functions (Biancolini et al. 2020). A radial basis function (RBF) \(\varphi \in \mathcal{R}\) is a function, that depends on the Euclidean distance \(r = \lVert \mathbf{x}-\hat{\mathbf{x}} \rVert\) to some fixed point \(\hat{\mathbf{x}}\). We can use radial basis functions to interpolate the deformation of the mesh for a set of given control points \(\hat{\mathbf{x}}_i\) as \[\delta(\mathbf{x}) = \sum_{l=1}^P \gamma_i \varphi(\lVert \mathbf{x}-\hat{\mathbf{x}}_i \rVert).\]

The interpolation should coincide with prescribed deformations at the control points \(\delta_i\), i.e. \[\delta(\mathbf{x}_i) = \delta_i \quad \forall i \in [1, P]\] must hold. To find weights fulfilling this condition, we can solve the linear system of equations \[\mathbf{M} \cdot \pmb{\gamma} = \pmb{\delta}\] with the interpolation matrix \[M_{ij} = \varphi(\lVert \mathbf{x}_i-\hat{\mathbf{x}}_j \rVert) \quad i,j \in [1, P].\] Using the control points as design variables, we can modify the shape smoothly with a linear combination of basic shape modifications. The type of RBF constraints the solution space of shapes that can be achieved by shape optimization.

We may reformulate the optimization problem now in terms of the design variables \(\delta_i\) as \[\begin{aligned} \min_{\pmb{\delta}} \quad & C(\pmb{\delta}) = \mathbf{f} \cdot \mathbf{u}(\pmb{\delta})\\ \textrm{s.t.} \quad & g(\pmb{\delta}) = \mathbf{d} \cdot \mathbf{A}(\pmb{\delta}) - V_0 \le 0 \\ & \pmb{\delta} \in \mathcal{D} \end{aligned} \label{eq:shape_optimization_morph}\]

Just like in shape optimization of trusses, target function and constraint function are generally non-linear and non-separable. Hence, we use MMA to approximate the problem and solve it equivalent to Equations [eq:shape_lagrangian] to [eq:shape_dual_solution]. We do not elaborate the analytical derivation of the gradients (this can be found e.g. in Chapter 6.3.2.2 in (Christensen and Klarbring 2008)), but simply rely on automatic differentiation to compute \(\frac{\partial C}{\partial \delta_i}\) and \(\frac{\partial g}{\partial \delta_i}\).