1 Introduction

The simple and effective 99- and 88-line Matlab implementation of the classical minimum compliance topology design problem (Bendsøe and Sigmund 2003; Andreassen et al. 2011) is well known in the topology optimization community. At the heart of the implementation is the optimality criteria (OC) design update (Bendsøe and Sigmund 2003). Given a design variable array \({\varvec{x}}^{k}\) with n elements, in some iteration k (also the zeroth iterate), the design is updated according to

$$\begin{aligned} x^{k+1}_i = x_i^k \left( \frac{ -{\textrm{d}}^k_{x_i} f_0}{ \lambda \, {{\textrm{d}}^k_{x_i} f_1}} \right) ^{\eta } , \end{aligned}$$
(1)

with bounds applied to each value in the new design iterate \(\varvec{x}^{k+1}\), such that

$$\begin{aligned} x^{k+1}_i \leftarrow \min ( \max ({\check{x}}^k_i, x^{k+1}_i ), {\hat{x}}^k_i ) , \end{aligned}$$
(2)

for \(i = 1,2,\ldots ,n\). The design variable bounds in each design iterate, \({\check{\varvec{x}}}^k\) and \({\hat{\varvec{x}}}^k\), are determined by a prescribed move limit and the global design variable bounds at 0 and 1. The \(\eta\) is a ‘tuning parameter’ (or ‘numerical damping’ coefficient), typically set at 0.5. The \({\textrm{d}}^k_{\varvec{x}} f_0\) and \({\textrm{d}}^k_{\varvec{x}} f_1\) are the sensitivitiesFootnote 1 of the objective and constraint function, evaluated at design iterate k. When applied to the classical minimum compliance topology design problem, \(f_0\) is the compliance objective and \(f_1\) is the resource constraint. The \(\lambda\) is the Lagrange multiplier (dual variable) associated with the constraint, computed with a bisectioning algorithm.

With fixed (constant) external loads, the minimum compliance objective function is monotonic: the first-order derivatives are negative \({{\textrm{d}}_{x_i}^k f_0} \le 0 \, \forall i\) at any design \(\varvec{x}^k\). The resource constraint has strictly positive sensitivities \({{\textrm{d}}_{x_i}^k f_0} > 0 \, \forall i\), and the Lagrange multiplier is positive \(\lambda > 0\) due to activity of the inequality constraint in negative null form. Therefore, the square-root implicated by \(\eta = 0.5\) in the OC update (1) is guaranteed to yield real numbers. Bendsøe and Sigmund (2003) address the issue which arises when a positive number under the square-root cannot be guaranteed. For example, in eigenvalue maximization problems, it is suggested that the OC update is modified to

$$\begin{aligned} x^{k+1}_i = x_i^k \left( \frac{ \max \left( -{\textrm{d}}^k_{x_i} f_0, 0\right) }{ \lambda \, {{\textrm{d}}^k_{x_i} f_1}} \right) ^{\eta } . \end{aligned}$$
(3)

The same is done in Bendsøe and Sigmund’s (2003) treatment of (linear) compliant mechanism design, with \(\eta\) ‘sometimes chosen a bit lower’ to promote stable convergence. A smaller value for \(\eta\) corresponds to reciprocal intervening variables with a larger (negative) exponent, increasing the effective curvature of the resulting objective function approximation (Groenwold and Etman 2008b). Generalized and adaptive intervening (or intermediate) variables form the basis of the very successful CONLIN (Fleury 1989) and MMA (Svanberg 1987) algorithms.

Arguably, the simplest topology optimization problem with a non-monotonic objective function is the classical minimum compliance design problem posed with design-dependent self-weight loads (i.e. body forces). See Bruyneel and Duysinx (2005). In classical minimum compliance design with fixed external loading, increasing the material at a point in a structure can only decrease the total elastic energy—proportional to compliance. However, if the load on the structure is the self-weight, then a non-trivial trade-off occurs: increasing the amount of material at a point in the structure increases its stiffness, but it also increases the load, and, depending on the configuration, it may increase the total elastic energy. Considering only self-weight loads, without a constraint on the minimum amount of material to be distributed, the global optimum is, strictly speaking, the trivial all-void design associated with no loading whatsoever.

Bruyneel and Duysinx (2005) summarize the difficulties in the self-weight problem, as follows: (i) non-monotonicity of the compliance objective function, (ii) the ‘unconstrained character of the optimum’, and (iii) the ‘parasitic effect for low densities’.Footnote 2 To address (iii), Bruyneel and Duysinx (2005) introduce a modification to the SIMP material law, making it linear at small design variable values. To address (ii), a constraint on the minimum amount of material to be distributed is introduced. And to address (i), the ‘structural [response function] approximations’ in the MMA are modified.

With the work of Bruyneel and Duysinx (2005) serving as backdrop, and prompted by contemporary works on the subject (Kumar 2022; Garaigordobil et al. 2022), we wish to show that the self-weight problem can also be solved with a design update based on a quadratic program (QP) with an adaptive curvature approximation, without modifying the material interpolation scheme. That is to say, we retain a straight-forward optimization problem formulation:

$$\begin{aligned} \begin{aligned} \underset{\varvec{x}}{\text {min}}{} & {} {}&{f}_0\left[ \varvec{x}\right] = \varvec{u}[\varvec{x}] \cdot \varvec{K}[\varvec{x}] \varvec{u}[\varvec{x}] \\ \text {subject to}{} & {} {}&{f}_1[\varvec{x}] = 1 - {\textstyle \sum _i} \tfrac{\rho _i}{v{n}} \le 0 \, , \\{} & {} {}&0 \le x_i \le 1, \quad i = 1,2,\ldots ,{n} \, , \end{aligned} \end{aligned}$$
(4)

wherein \(\varvec{u}\) is the displacement degrees-of-freedom dependent implicitly on the design variables \(\varvec{x}\), \(\varvec{u}[\varvec{x}]\), via solution of the linear, finite element discretized equilibrium equations \(\varvec{K}[\varvec{x}] \varvec{u}= \varvec{g}[\varvec{x}]\) in each design iterate. The design variables are subjected to density filtering \(\varvec{\rho } = \varvec{H}\varvec{x}\) and SIMP is applied (Bendsøe and Sigmund 2003) in the construction of the global stiffness matrix \(\varvec{K}\). The self-weight (body-force) per element \(\varvec{g}_e\) is taken to scale linearly with the corresponding (filtered) density variable \(\rho _e\), and the resource constraint is modified to enforce distribution of a minimum amount of material in terms of a volume fraction \(v\). See the Appendix for the treatment of the more general problem with the addition of a fixed external load and a second resource constraint on the maximum amount of material to be distributed—a variant also studied by Bruyneel and Duysinx (2005). The reader is invited to consult the cited literature and the Matlab implementation provided in the supplementary material for more detail.

In the following section theory is recited, after which the implementation and numerical demonstrations are presented.

2 An adaptive QP

Groenwold et al. (2007) generalize the notion of intervening variable functions with the introduction of incomplete Taylor series expansions for function approximation. Arguably, the simplest response function approximation is quadratic with a single (scalar) curvature term \(c^k\), written as

$$\begin{aligned} {\tilde{f}}_0[\varvec{x}] =&{f}^k_0 + {\textrm{d}}^k_{\varvec{x}} {f}_0 \cdot (\varvec{x}- \varvec{x}^{k}) \nonumber \\&+ \frac{1}{2} c^k (\varvec{x}- \varvec{x}^{k}) \cdot (\varvec{x}- \varvec{x}^{k}) \, . \end{aligned}$$
(5)

In general, the quadratic (or curvature) terms may be constructed in all manner of ways, including the approximation (5) to the response function linearized in terms of intervening variables (Groenwold et al. 2010). The so-called spherical quadratic approximation (Groenwold et al. 2007, 2010) is obtained with

$$\begin{aligned} c^k = 2 \frac{{f}_0^{k-1} - {f}_0^{k} - {\textrm{d}}^k_{\varvec{x}} {f}_0 \cdot (\varvec{x}^{k-1} - \varvec{x}^{k}) }{ (\varvec{x}^{k-1} - \varvec{x}^k)\cdot (\varvec{x}^{k-1} - \varvec{x}^k)} , \end{aligned}$$
(6)

a result of requiring the function approximation to render the same zero-order information \({f}_0^{k-1}\) at the previous design iterate \(\varvec{x}^{k-1}\). Strict convexity of the objective function approximation is enforced by

$$\begin{aligned} c^k \leftarrow \max (\epsilon , c^k) , \end{aligned}$$
(7)

with \(\epsilon\) a small positive number—e.g. 1E\(-9\). Limiting ourselves to a single linear constraint, herein, we can write an OC-like update of the form

$$\begin{aligned} x^{k+1}_i = x_i^k - \frac{{\textrm{d}}^{k}_{x_i}{f}_0 + \lambda {\textrm{d}}^{k}_{x_i} {f}_1}{ c^k } , \end{aligned}$$
(8)

which results from requiring the Lagrangian of the QP associated with the iterate to be stationary with respect to \(\varvec{x}\). In this case, the Lagrange multiplier (dual variable) \(\lambda\) can be determined with the same bisectioning algorithm employed in the OC update, and the bounds (2) are applied as before. This is a special case of the dual due to Falk (1967), as discussed in more detail in the Appendix.

Due to the simplicity of the QP, and because the problem is simply constrained, so-called global convergence (enforced convergence) of the GCMMA (Svanberg 2002) is readily applied by requiring a descent step in each design iterate \({f}_0^{k+1} < {f}_0^{k}\). If a descent step is not taken, the previous iterate is restored and the curvature approximation of the objective function is increased by some factor, e.g. \(c^k \leftarrow 2 c^k\). We return to this topic later. See Groenwold and Etman for a detailed description of enforced convergence by conditional acceptance of iterates (Groenwold and Etman 2010a).

3 Implementation and demonstration

We take as basis the 88-line Matlab implementation of the classical minimum compliance problem with fixed external loads and the OC solution procedure (Andreassen et al. 2011). Please see the modified code topqpsw.m provided in the Supplementary Material. The fixed external loading is removed, design-dependent (distributed) loads in the y direction are introduced (lines 52–56), the objective sensitivities are modified (line 64), density filtering is used exclusively, and the resource constraint and its sensitivities are modified to limit the minimum amount of material to be distributed (lines 65–66). Please see the Appendix for the treatment of the problem with additional external loading and a second resource constraint on the maximum amount of material to be distributed. Variables and logic are introduced to calculate the adaptive curvature term (6) (line 85) and optionally enforce convergence (lines 77–91) with the input flag ec. Other modifications include scaling of the objective function (lines 68–71), provision of additional output information—e.g. the ‘black-and-white’ fraction (B &W) of the design—and the convergence criteria (line 49).

The OC update is replaced by the QP update (8) (lines 93–102):

figure a

A smaller move limit than typical is set (0.1 vs 0.2) in the presence of this less well-behaved, non-monotonic objective function.

The first example concerns the well-known MMB beam design domain and kinematic boundary conditions (symmetry and roller support). It may be executed with topqpsw(180,60,0.1,3.,3.,0). The result—the filtered design variable field \(\varvec{\rho }\) at termination—is given in Fig. 1. The number of iterations k, the objective function value \({f}_0\), the volume fraction V, and the black-and-white fraction \(B \& W\) are given in the caption of the figure. The constraint on the minimum amount of material (volume) to be distributed is not active. Executed with a larger filter radius (e.g. 4), increasing the minimum feature size of the design, increases the volume fraction at solution further (e.g. 0.19). It is understood that a particular amount of material is required to construct a topology which connects the roller supports in some stiffness optimized manner with respect to its own weight.

Fig. 1
figure 1

Self-weight MBB beam minimum compliance topology design with QP update: k = 132, \({f}_0\) = 6.194, V = 0.14, \(B \& W\) = 0.96

Full size image

The second example is a self-weight arch. The roller is changed to a pin support, with line 20 replaced by

figure b

The example may then be executed with topqpsw(120,120,0.1,3,3,0). This case confirms the difficulties in the self-weight problem: iterates oscillate with low-density material and termination occurs at a design with a low black-and-white fraction (0.87). This is remedied by running the case with enforced convergence; topqpsw(120,120,0.1,3,3,1). The result is given in Fig. 2. In this case, the constraint on the minimum amount of volume is active. If executed with a lower volume fraction (e.g. 0.05), it is recommended to modify the starting point (line 45) to, e.g. x = repmat(0.5,nely,nelx); in order that, too much material is not zeroed in the initial iterates, causing the run to have a bad starting position.

Fig. 2
figure 2

Self-weight arch minimum compliance topology design: k = 171, \({f}_0\) = 0.1261, V = 0.10, \(B \& W\) = 0.99

Full size image

The same is achieved in 3D using a simple C/C++ implementation (provided on request). Figures 3 and 4 illustrate the result of a 3D version of the self-weight arch problem considered above, with volume fractions of 0.05 and 0.07, respectively. The density field is interpolated to the nodes and clipped at 0.25. The domain is 100\(\times\)100\(\times\)50 elements. In both cases, the volume constraint is active at the solution.

Fig. 3
figure 3

3D self-weight arch minimum compliance topology design: k = 263, \({f}_0\) = 0.413, V = 0.05, \(B \& W\) = 0.99

Full size image
Fig. 4
figure 4

3D self-weight arch minimum compliance topology design: k = 193, \({f}_0\) = 0.472, V = 0.07, \(B \& W\) = 0.99

Full size image

4 Concluding remarks

It is fairly remarkable that the self-weight problem, typically addressed with modification of material interpolation schemes and sophisticated algorithms, can be solved with only a simple QP design update with an adaptive curvature approximation. The QP update presented is a special case of the dual of Falk (1967), with which multiple constraints may be treated—see the Appendix. Curvature information (if present) of constraint functions may also be taken into account (Groenwold 2012).