1 Introduction

Truss design problems are often formulated based on the so-called ground structure approach (Dorn et al. 1964), in which a set of joints are distributed in the design space and are connected by potential bars. We are here concerned with a truss design problem where the goal is to optimize both the topology and geometry of the structures, i.e., when the design variables are the cross-sectional areas of the bars and the coordinates of the joints. These problems have been studied in many articles, for example by Dobbs and Felton (1969), Kirsch (1990b), Ben-Tal et al. (1993), Bendsøe et al. (1994), Pedersen (1972), Sergeyev and Pedersen (1996), Achtziger (1998), Tejani et al. (2018), and Miguel and Miguel (2012), to mention just a few. The problems are highly nonlinear, mainly due to the variation of the joint coordinates. However, the models are known to obtain optimal designs that are more practically useful as they require less post-processing and contain fewer joints connecting bars in the design space. Alternatively, one can also obtain efficient (least weight) structures by considering many more joints with fixed coordinates and solving large-scale linear or nonlinear topology optimization problems (Bendsøe and Sigmund 2003; Jarre et al. 1998; Gilbert and Tyas 2003; Sokół and Rozvany 2013). However, the resulting designs usually have many active joints and bars, and are thus often far from practical.

Due to the high nonlinearity of the geometry and topology optimization problem formulations, several solution techniques have been proposed, mainly to improve the computation tractability of the problems, e.g., Imai and Schmit (1981), Ringerts (1985), Svanberg (1981), Kočvara and Zowe (1996), Achtziger (2007), and He and Gilbert (2015). One of the most common techniques is to use the so-called alternating method, for example described in Ringerts (1985) and Kočvara and Zowe (1996), where the problems are solved for a fixed geometry but with a variation in topology and vice versa. Such a technique involves a block-coordinate type approach and obtains designs which satisfy the optimality conditions for the solved subproblems, which are often acceptable. In other studies, the topology and geometry of the structure are optimized simultaneously (see Achtziger (2007)), and it has been proved that the solution is (locally) optimal. In both approaches, some studies also use first-order information (e.g., Svanberg (1981)), while others attempt to use the second-order primal-dual method (e.g., Imai and Schmit (1981) and He and Gilbert (2015)) to improve convergence properties. For an overview of these approaches and other solution strategies, we refer the reader to Achtziger (2007).

Many studies on the optimization of truss structures incorporate far more constraints than the classical formulations, where the weight or volume is minimized subject to constraints on the compliance or the other way round, to improve the practicality of the optimal designs. These include constraints on stresses and/or on local buckling based on Euler’s formula (e.g., Kirsch (1990a), Guo et al. (2001a), Stolpe and Svanberg (2001), Stolpe and Svanberg (2003), Zhou (1996), Achtziger (1999), Rozvany (1996), Guo et al. (2001b), Guo et al. (2005), and Mela (2014)), the indirect imposition of nodal stability via nominal forces and associated constraints (e.g., Tyas et al. (2006) and Descamps and Coelho (2014)), and direct imposition of global stability constraints (e.g., Ben-Tal et al. (2000), Levy and Su (2004), Stingl (2006), Evgrafov (2005), and Tugilimana et al. (2018)). There are also variants of these formulations incorporating buckling constraints using frame structures (Torii et al. 2015; Mitjana et al. 2019), beam modeling (Madah and Amir 2017), and continuum structures (Ferrari and Sigmund 2019).

In this paper, we address truss design problems with global stability constraints using a linear buckling model that is formulated as a nonlinear semidefinite programming problem. Such problems have been extensively studied by Ben-Tal et al. (2000), Kanno et al. (2001), Levy and Su (2004), Kočvara (2002), Stingl (2006), and Evgrafov (2005) and solved, for example by Fiala et al. (2013) and Kočvara and Stingl (2003), but always for cases where the nodes or joints are assumed to be fixed. In Weldeyesus et al. (2019), an invariant large-scale problem formulation similar to those in Kočvara (2002), obtained by relaxing nonlinear kinematic compatibility constraints, but still for fixed joints, has been solved using a customized primal-dual interior point method by exploiting the sparsity and low-rank properties of the associated element stiffness matrices (Bendsøe et al. 1994; Ben-Tal 1993; Achtziger et al. 1992; Bendsøe and Sigmund 2003), the use of column generation (member adding) procedures (Gilbert and Tyas 2003; Sokół and Rozvany 2013; Weldeyesus and Gondzio 2018), and a warm-start strategy (Weldeyesus and Gondzio 2018).

The goal of this paper is to introduce geometry optimization to existing truss topology optimization with global stability constraints problems via nonlinear semidefinite programming. In particular, the goal is to extend the models proposed by Kočvara (2002) and Weldeyesus and Gondzio (2018), and to show that there are problems of this type that can be solved using a standard primal-dual interior point method. This is perhaps surprising, considering the severity of the nonlinearity and nonconvexity of the problem formulations. We refer the reader to Yamashita and Yabe (2015) for a recent survey of numerical methods for nonlinear semidefinite programming, and further discussions on ongoing challenges in the field.

One drawback of the solutions obtained by applying geometry optimization is that the resulting optimal designs depend hugely on the initial positions of, and the number of, joints. As an attempt to overcome this, and another challenge associated with numerical instability, namely when some joints come too close to each other leading to singularity, we perform adaptive geometry and topology optimization, inspired by He and Gilbert (2015). This is an iterative procedure where the problems are initially solved by restricting the movement of the joints to smaller regions, and then progressively updating these. With this approach, the joints can ultimately navigate larger regions, which could be far beyond the design space defined by the initial joint configuration. During the procedure, inactive joints, i.e., the joints connected entirely to thin bars that have cross-sectional areas below a prescribed threshold, are removed. Moreover, we perform some of the common techniques in geometry optimization such as node merging when joints are too close, and node melting when a joint just connects two active collinear bars. The overall procedure can amount to a certain extent to post-processing or rationalization of an optimal truss design (He and Gilbert 2015).

The paper is organized as follows. In Section 2, we present the essential mathematical background for geometry optimization of trusses. In Section 3, the truss geometry and topology optimization with global stability constraints problem formulation is presented, modelled via nonlinear semidefinite programming. We describe the general framework of the primal-dual interior point method in Section 4. The numerical experiments, and an adaptive geometry and topology optimization scheme are described in Section 5. Finally, conclusions and future research directions are presented in Section 6.

2 Background

In this section, we describe the essential mathematical concepts that are useful for modelling truss geometry and topology optimization with global stability constraints problems. Part of this section closely follows Weldeyesus et al. (2019) and Kočvara (2002). We adopt the ground structure approach (Dorn et al. 1964), i.e., a finite set of joints are distributed in a given N-dimensional design space, where N ∈ {2, 3}, and the joints are connected by some potential bars. Let d be the number of the joints with \(\bar {v}_{j}, j=1,...,d\) corresponding to the coordinates of the positions of the joints. Hence, \(\bar {v}_{j}=(x_{j},y_{j})\) if N = 2, and \(\bar {v}_{j}=(x_{j},y_{j},z_{j})\) if N = 3. Note that these coordinates will be considered as the initial positions of the joints in this paper. Let m be the number of bars with the cross-sectional areas ai, i = 1,..., m. In geometry optimization, the joints are allowed to move within certain limits; we refer the reader to Achtziger (2007) and He and Gilbert (2015) for a brief discussion of various types of admissible move limits. Throughout the course of the paper, we assume that the positions of supported and loaded joints are always fixed, and the other joints, say d0, are allowed to move within a move limit defined by the neighborhood:

$$ \mathcal{V}=\mathcal{V}_{1} \cap \mathcal{V}_{2} $$
(1)

where \(\mathcal {V}_{1}\) is a region defined by balls of radii r around the joints and is given by:

$$ \mathcal{V}_{1}=\{v \in\mathbb{R}^{d_{0}N}| ||v_{j}-\bar{v}_{j}|^{2}\leq {r^{2}_{j}}, j=1,...,d_{0}\}, $$
(2)

where ||⋅|| is the Euclidean norm, and \(\mathcal {V}_{2}\) is a region described by a set of linear constraints defined by:

$$ \begin{aligned} \mathcal{V}_{2}=\{v \in\mathbb{R}^{d_{0}N}| \bar{v}_{j,k}-v_{j,k}^{min}\leq v_{j,k}\leq \bar{v}_{j,k}+ v_{j,k}^{max}, \\ j=1,...,d_{0}, k=1,...,N\}. \end{aligned} $$
(3)

Note that, we have d0 < d.

Remark 1

For most joints, particularly those inside the design domain, we set \(\min \limits \{v_{j,k}^{min},v_{j,k}^{max}\}\geq r_{j}\) making the interval constraints (3) inactive. In some cases, for example, when the joints are required to remain in the design domain, we set \(v_{j,k}^{min}\) or \( v_{j,k}^{max}\) to values that are less than rj, or even to zero.

Let \(v^{(2)}_{i}\) and \(v^{(1)}_{i}\) denote the coordinates of the start and end joints of the bar i, i = 1,⋯ , m. In order to avoid singularity (or non-differentiability, for example, when \(v^{(2)}_{i}=v^{(1)}_{i}\) in (5)), the radius rj of each of the balls in (2) is chosen to satisfy:

$$ 0<r_{j}=\frac{1}{2}\min\{||\bar{v}_{j}-\bar{v}_{p}||, p\in I\}-\varepsilon, $$
(4)

where I is the set of indices of the joints connected to joint j, and ε > 0. For every bar i, i = 1, ⋯ , m, its length li(v) is given by

$$ l_{i}(v)=||v^{(2)}_{i}-v^{(1)}_{i}||, $$
(5)

and the associated vector of direction cosines \({\gamma ^{e}_{i}}(v)\) is defined by

$$ {\gamma^{e}_{i}}(v)=\frac{1}{l_{i}(v)}(v^{(2)}_{i}-v^{(1)}_{i})^{T}. $$
(6)

Let n = Nd − n0, where n0 is the number of fixed degrees of freedom. Then, the corresponding global vector \(\gamma _{i}(v) \in \mathbb {R}^{n}\) can be appropriately constructed by embedding \((-{\gamma ^{e}_{i}}(v),\)\({\gamma ^{e}_{i}}(v) )^{T}\) and setting all of the remaining entries to zero.

Now, given an external load \(f\in \mathbb {R}^{n}\), the resulting displacement \(u\in \mathbb {R}^{n}\) satisfies the (linear) elastic equilibrium equation:

$$ K(a,v)u = f, $$
(7)

where the stiffness matrix K(a, v) is computed as:

$$ K(a,v)= \sum\limits^{m}_{i=1} a_{i} \frac{E}{l_{i}(v)} {\gamma}_{i} (v){\gamma}_{i}^{T}(v), $$
(8)

with E being Young’s modulus of the material.

Next, we define the so-called geometrical stiffness matrix G(a, v, u) as given by:

$$ G(a,v,u)=\sum\limits^{m}_{i=1} \frac{a_{i}E\gamma_{i} (v)^{T}u}{{l^{2}_{i}}(v)} ({\delta_{i} (v)}{\delta_{i} (v)}^{T}+{\eta_{i} (v)} {\eta_{i} (v)}^{T} ). $$
(9)

The vectors δi(v) and ηi(v) are determined so that γi(v), δi(v), and ηi(v) are mutually orthogonal (Kočvara 2002). Hence, there are many possible ways of choosing these vectors. When all of the joints are fixed, i.e., when solving topology optimization, these have been computed as the orthogonal basis of the null space of \({\gamma _{i}^{T}}\) in Kočvara (2002). We follow a similar approach but additionally derive the vectors δi(v) and ηi(v) explicitly as described in Section 2.1 since we need to compute the derivatives of these vectors with respect to the coordinates of the joints during the optimization process.

2.1 Computing the vectors δi(v) and ηi(v)

When N = 2 (in this case, there is no ηi(v)), then we have

$$ {\delta^{e}_{i}}(v)=A_{2\times2}{\gamma^{e}_{i}}(v), $$
(10)

where the rotation matrix A2×2 is given by

$$ A_{2\times2}= \left( \begin{array}{ll} 0 &-1\\ 1& \phantom{-}0\end{array} \right). $$
(11)

In this case, the set \(\{{\gamma ^{e}_{i}}(v),{\delta ^{e}_{i}}(v)\}\) is an orthonormal basis for \(\mathbb {R}^{2}\).

When N = 3, first let us define the following three rotation matrices:

$$ \begin{aligned} A^{(1)}_{3\times3}&=\left( \begin{array}{lll} 0 &0 & \phantom{-}0\\0&0&-1\\ 0 &1 &\phantom{-}0\end{array}\right), A^{(2)}_{3\times3}=\left( \begin{array}{lll} 0 &0 & -1\\0&0&\phantom{-}0\\ 1 &0 &\phantom{-}0\end{array}\right),\\ A^{(3)}_{3\times3}&=\left( \begin{array}{lll} 0 &-1 & 0\\1&\phantom{-}0&0\\ 0 &\phantom{-}0 &0\end{array}\right). \end{aligned} $$
(12)

Then, we determine the non-zero entries of δi(v) and ηi(v) using the following steps.

  1. Step 1

    Compute \({\delta ^{e}_{i}}(v)\) as:

    $$ {\delta^{e}_{i}}(v) =\frac{A^{(j)}_{3\times3} {\gamma^{e}_{i}}(v)}{||A^{(j)}_{3\times3} {\gamma^{e}_{i}}(v)||}, $$
    (13)

    where j is the index of \(|({\gamma ^{e}_{i}}(v))_{j}|\) with the smallest magnitude.

  2. Step 2

    Compute \({\eta ^{e}_{i}}(v)\) using the vector product:

    $$ {\eta^{e}_{i}}(v) ={\gamma^{e}_{i}}(v){\times\delta^{e}_{i}}(v). $$
    (14)

Remark 2

In the implementation, if vector \({\gamma ^{e}_{i}}(v)\) contains two entries with the same smallest magnitude, then we set j in 2.1 to the smallest index of the corresponding entries. If all entries are equal, then we use j = 1.

Remark 3

We compute \({\delta ^{e}_{i}}(v)\) as given in (13), i.e., keep the two largest absolute value entries, to ensure robustness.

Once the vectors \({\delta ^{e}_{i}}(v)\) and \({\eta ^{e}_{i}}(v)\) are determined as in (13) and (14), then it is straightforward to verify that the set \(\{{\gamma ^{e}_{i}}(v),{\delta ^{e}_{i}}(v),{\eta ^{e}_{i}}(v)\}\) forms an orthonormal basis for \(\mathbb {R}^{3}\).

3 The problem formulation

In this section, we present the formulation for the minimum weight truss geometry and topology optimization with global stability constraints problem, which is:

$$ \begin{aligned} & \underset{a,v,u}{\text{minimize}} & & l(v)^{T}a \\ & \text{subject to} &&f^{T}u\leq \zeta\\ & & & K(a,v)u=f\\ & & & K(a,v)+\tau G(a,v,u)\succeq 0 \\ & & & v\in \mathcal{V}\\ & & & a \geq 0, \end{aligned} $$
(15)

where ς is a given bound on the compliance, and the set \(\mathcal {V}\) is as in (1). The matrix inequality constraint addresses global instability via a linearized buckling analysis in which the design load factor τ should be at least 1. In that case, the obtained optimal design is stable for the load τf.

If we additionally assume that a > 0 and using (4), we have K(a, v) ≻ 0. Then, we can solve for u in the elastic equilibrium equation (7) and remove it from (15) to obtain other equivalent nested formulations. Furthermore, by applying the so-called Schur complement method, problem (15) could be rewritten as:

$$ \begin{aligned} & \underset{a,v}{\text{minimize}} & & l(v)^{T}a \\ & \text{subject to} &&\left( \begin{array}{lll} \zeta&&f^{T} \\f&&K(a,v)\end{array}\right)\succeq 0\\ & & & K(a,v)+\tau \tilde{G}(a,v)\succeq 0 \\ & & & v\in \mathcal{V}\\ & & & a \geq 0, \end{aligned} $$
(16)

where

$$ \begin{aligned} \tilde{G}(a,v)=\sum\limits^{m}_{i=1} \frac{a_{i}E\gamma_{i} (v)^{T}K^{-1}(a,v)f}{{l^{2}_{i}}(v)} ({\delta_{i} (v)}{\delta_{i} (v)}^{T}+\\{\eta_{i} (v)} {\eta_{i} (v)}^{T} ). \end{aligned} $$

However, in the numerical experiments presented in Section 5, we only solve problem (15).

Remark 4

Problems (15) and (16) can be considered as a natural extension of the model described by Kočvara (2002), where it is formulated only for topology optimization with global stability constraints, i.e., all nodal coordinates (or joints) are fixed.

Remark 5

In Weldeyesus et al. (2019), a variant of formulation (15) written based on member forces was addressed for fixed joints. After relaxing nonlinear constraints, it is shown that the problem can be solved for a very large number of bars by applying an adaptive method, i.e., a technique in which a sequence of much smaller problems is solved to obtain the solution of the original large-scale problem. The computational gain is remarkable and the qualities of the solution of the relaxed problem are good. However, since the joints are fixed, a large number of nodes are required to obtain least weight structures. This can lead to optimal designs that contain many active joints and bars. In the present contribution, the focus is on solving problems with somewhat fewer nodes/joints, but on exploiting geometry optimization to realize enhanced solutions.

Problem (15) (and (16)) is a highly nonlinear and noncovex semidefinite program and is very difficult to solve (Yamashita and Yabe 2015). To the best of the authors’ knowledge, no solutions to the problem have been proposed in the current literature.

As mentioned in Section 1, it is generally a challenge to solve truss geometry and topology optimization with or without stability constraints. In many studies, the so-called alternating optimization method is used, in which the optimization problem is split into two subproblems, i.e., fixing joints and optimizing with respect to the cross-sectional areas of the bars, and vice versa (Ringerts 1985; Kočvara and Zowe 1996). This block-coordinate type approach delivers acceptable optimal designs that satisfy the optimality conditions of the subproblems. Another approach is to solve the problem simultaneously and obtain some local solutions (Achtziger 2007).

In this paper, we solve problem (15) simultaneously for both geometry and topology optimizations, using a second-order primal-dual interior point method for nonlinear semidefinite programming described below, i.e., in Section 4.

4 The primal-dual interior point framework

In this section, we present an overview of the interior point method and an algorithm we applied to solve problem (15), which is in some sense similar to that described in Section 3 of Weldeyesus and Stolpe (2015), except for some slight reformulation of the general nonlinear semidefinite programming and associated notations, now adopted from Yamashita et al. (2012). This is to make the flow of the presentation consistent with this paper.

Consider the nonlinear semidefinite programming problem of the form:

$$ \begin{aligned} & \underset{x \in\mathbb{R}^{m}}{\text{minimize}} & & f(x)\\ & \text{subject to} & & g(x)=0\\ & & & \mathcal{A}(x)\succeq0. \end{aligned} $$
(17)

The functions \(f:\mathbb {R}^{m}\rightarrow \mathbb {R}\), \(g:\mathbb {R}^{m}\rightarrow \mathbb {R}^{k}\), and \(\mathcal {A}:\mathbb {R}^{m}\rightarrow \mathbb {S}_{+}^{n}\) are assumed to be sufficiently smooth, and \(\mathbb {S}_{+}^{n}\) is the cone of positive semidefinite matrices in the space \(\mathbb {S}^{n}\) of symmetric n × n matrices.

After introducing a barrier parameter μ > 0, the associated barrier problem is:

$$ \begin{aligned} & \underset{x\in\mathbb{R}^{m}}{\text{Minimize}} & & f(x)-\mu\ln(\det(\mathcal{A}(x))\\ & \text{Subject to} & & g(x)= 0. \end{aligned} $$
(18)

The Lagrangian to problem (18) is:

$$ \mathcal{L}_{\mu}(x,\lambda)=f(x)-\mu\ln(\det(\mathcal{A}(x))+\lambda^{T}g(x), $$

where \(\lambda \in \mathbb {R}^{k}\) is a Lagrangian multiplier. Introducing the additional matrix variable \(Z:=\mu (\mathcal {A}(x))^{-1}\), we have:

$$ \nabla_{x}\mathcal{L}_{\mu}(x,\lambda)=\nabla_{x}f(x)-\mathcal{G}(x)Z+\nabla_{x}g(x)^{T}\lambda $$
(19)

where:

$$ \mathcal{G}(x)Z=\left( \begin{array}{ll} \langle\frac{\partial \mathcal{A}(x)}{\partial x_{1}},Z\rangle\\ ~~~~~~~\vdots\\ \langle\frac{\partial \mathcal{A}(x)}{\partial x_{m}},Z\rangle \end{array}\right). $$
(20)

Then, first-order optimality conditions of the barrier problem (18) are:

$$ \begin{array}{@{}rcl@{}} \nabla_{x}\mathcal{L}_{\mu}(x,\lambda)=0 \end{array} $$
(21a)
$$ \begin{array}{@{}rcl@{}} g(x)= 0 \end{array} $$
(21b)
$$ \begin{array}{@{}rcl@{}} \mathcal{A}(x)Z-\mu I=0. \end{array} $$
(21c)

To apply the Newton method to solve the system of nonlinear equations (21a), first we need to maintain symmetry, which is done by replacing equation (21c) with

$$ H_{P}(\mathcal{A}(x)Z)=\mu I, $$
(22)

where the linear operator \(H_{P}:\mathbb {R}^{n\times n}\rightarrow \mathbb {S}^{n}\), introduced in Zhang (1998), is defined as follows:

$$ H_{P}(Q):=\frac{1}{2}\left( PQP^{-1}+(PQP^{-1})^{T}\right) $$

with \(P\in \mathbb {R}^{n\times n}\) being some non-singular matrix. There are various ways to choose the matrix P. In this paper, we use NT direction (Nesterov and Todd 1997; 1998).

Now, to get the search directions (Δx,Δλ,ΔZ), we can apply the Newton method to (21a) with the last equation replaced by (22), and solve the system:

$$ \begin{aligned} \nabla^{2}_{xx}\mathcal{L}_{\mu}(x,\lambda){\Delta} x - \mathcal{G}(x){\Delta} Z + \nabla_{x}{g(x)^{T}}{\Delta} \lambda = -\nabla_{x}\mathcal{L}_{\mu}(x,\lambda)\\ \nabla_{x}{g(x)}{\Delta} x=-g(x)\\ \mathcal{E}{\Delta} X+ \mathcal{F}{\Delta} Z=\mu I-H_{P}(\mathcal{A}(x)Z), \end{aligned} $$
(23)

where \(\mathcal {E}=\mathcal {E}(x,Z)\) and \(\mathcal {F}=\mathcal {F}(x,Z)\) are the derivatives of \(H_{P}(\mathcal {A}(x)Z)\) with respect to x and Z, and \({\Delta } X = {\sum }_{i=1}^{m} {\Delta } x_{i} \frac {\partial \mathcal {A}(x)} {\partial x_{i}}\).

An overview of the interior point method is summarized in Algorithm 1. It has two loops. The norm of the optimality conditions in (23) is the stopping criteria, with μ = 0 for the outer loop and μ = μk for the inner loop. The parameters α and σ control step-length of search directions and the centrality of the points, respectively.

figure a

5 Numerical results

The interior point method has been implemented in MATLAB (R2018a). All numerical experiments have been performed on a PC equipped with an Intel(R) Core(TM) i5-8350U CPU running at 1.90 GHz with 8 GB RAM. In all examples, we use input data without units. We use the values of Young’s modulus E = 1 and the design loading factor τ = 1. In the plots showing optimal designs, only bars with cross-sectional area \(\geq 0.001a_{\max \limits }\) are shown. Balls are used to show the active joints connecting these bars.

In Algorithm 1, unless stated otherwise, we use σ = 0.5 and α = 0.8αmax, where αmax is the maximum step length such that the current point is positive (definite).

5.1 Fixed versus moving joints

In this section, we present three examples to demonstrate the use of geometry optimization (allowing the joints to move) to help minimize the volume of the optimal design. This is done by comparing solutions obtained with fixed and moving joints. Moreover, even though it is not the primary purpose the paper, we also demonstrate the importance of global stability constraints in these examples. We do this by comparing designs obtained with and without stability constraints. For more details on such stability constraints, see Kočvara (2002).

Example 1

We start with the L-shaped truss problem shown in Fig. 1a that has 132 members and has previously been solved by Levy and Su (2004) assuming fixed joints. It has overall dimensions 1 × 3 × 4, including a cut-out region of dimensions 1 × 2 × 3, with two point loads (0,0,− 0.001) applied simultaneously. The bound on the compliance is ς = 0.0005. When solving the topology optimization problem for fixed geometry without stability constraints, we obtain the design shown in Fig. 1b, consisting of two disjoint parallel planar trusses of total volume 7.6880. Now, including the stability constraints and solving the problem again for fixed joints, we obtain the design shown in Fig. 1c, where connectivity between the parallel planar trusses is established. In this case, the optimal volume is 7.8340.

Fig. 1
figure 1

Example 1: a design domain, boundary conditions, and loads; b side (left) and 3D (right) views of the optimal design without stability constraints; c–e same views of the optimal designs with stability constraints

Full size image

Next, we allow the joints to move. In this case, we consider two scenarios as follows.

  1. (i)

    The joints are allowed to move in a ball of radius r = 0.4, but still remain within the L-shaped design space of Fig. 1a. We obtain the design shown in Fig. 1d that has volume 7.1557, which is approximately 9% lower than that of the fixed joint design shown in Fig. 1c. When the joints are allowed to move, a part of the external envelope of the optimized structure appears curved in form, as shown in Fig. 1d.

  2. (ii)

    The joints are allowed to move in a ball of radius r = 0.4, but only the joints in the inner most surfaces are required to remain within the design space (Fig. 1a). In this case, we obtain the design shown in Fig. 1e that has a volume of 6.1192, which is approximately 21% lower than that of the fixed joint solution shown in Fig. 1c. As can be seen in Fig. 1e, some of the joints on the outermost edges have left the original design space (Fig. 1a), such that the entire external envelope of the optimized structure now appears curved in form.

The computational statistics for all cases are presented in Table 1. The algorithm obtained the solution within a reasonable number of iterations and CPU time.

Table 1 Example 1 (L-shaped problem): numerical results
Full size table

Finally, we note that there can be seen some degree of resemblance of the optimal designs for the L-shaped problem in Fig. 1e to the solutions obtained by Schwarz et al. (2018), Descamps and Coelho (2014), and Ohsaki and Hayashi (2017), where solutions are reported for geometry and topology optimization problems, but where other types of stability constraints are employed.

Example 2

We now solve the bridge problem of dimensions 8 × 1 × 1 as shown in Fig. 2a. The nodal loads are applied simultaneously and have magnitude (0,0,− 0.001); the bound on the compliance is ς = 0.003. When we solve the problem with and without stability constraints for fixed joints, we obtain the designs shown in Fig. 2b and c, with volumes of 10.3253 and 10.8904, respectively. The disjoint planar trusses in Fig. 2b are now connected in Fig. 2c, which is obtained by solving the problem with stability constraints. Next, we solve the problem for moving joints. First, the joints are allowed to move in a neighborhood of radius r = 0.4 but are restricted to remain within the design space. In this case, we obtain the design shown in Fig. 2d, with volume 8.5030, which is approximately 21% lower than that of the fixed joint design shown in Fig. 2c. Moreover, the structure clearly becomes more arch-like in form. Finally, when we remove the requirement for the joints to remain in the design space, we obtain the solution shown in Fig. 2e, with volume 5.7345, which is 47% lower than that of the original fixed joint design shown in Fig. 2c. This large reduction is principally because the arch has been able to extend upwards, above the original extent of the design space. The computational statistics are presented in Table 2; once again, the solution was found within a reasonable number of iterations and in a reasonable CPU time.

Fig. 2
figure 2

Example 2: a design domain, boundary conditions, and loads; b side (top left), top (lower left), and 3D (right) views of the optimal design without stability constraints; c–e same views of the optimal designs with stability constraints

Full size image
Table 2 Example 2 (bridge problem): numerical results
Full size table

Example 3

We now consider the tower problem shown in Fig. 3a which utilizes a design domain of dimensions 2 × 2 × 4. The load is (0,0,− 0.001) and the bound on compliance is ς = 0.00025. The goal of this example is to consider a condition where compressive global buckling is dominant. When the problem is solved assuming fixed joints and without stability constraints, the resulting optimal design has a volume of 0.0640, and comprises a chain of vertical bars with no bracing to connecting nodes, as shown in Fig. 3b. Next, global stability constraints are included, but the joints are kept fixed. In this case, the optimal design has a volume of 0.1536, a 140% increase. This increase is much more significant than observed in the previous examples considered in the paper. The associated optimal design is shown in Fig. 3c; here, all nodes are braced to ensure the design is stable. Finally, we solve the same problem but now allowing the joints to move in a ball of radius r = 0.3. The resulting optimized design is shown in Fig. 3d; here, the on-plan area of the tower reduces with increasing height. The design has volume 0.1230, a 20% reduction compared with fixed joint structure shown in Fig. 3c. The computational statistics are reported in Table 3. The global buckling modes for these structures are shown in Fig. 4a–c, respectively.

Fig. 3
figure 3

Example 3: a design domain, boundary conditions, and load; b optimal design without stability constraints; c–d 3D (left), middle (top), and side (right) views of the optimal design with stability constraints

Full size image
Fig. 4
figure 4

Example 3: a–c buckling modes for the designs shown in Fig. 3b–d respectively

Full size image
Table 3 Example 3 (tower problem): numerical results
Full size table

5.2 Adaptive geometry and topology optimization

We now consider development and application of an adaptive geometry and topology optimization technique, illustrated via a number of supporting numerical examples. The proposed technique is iterative in nature.

As mentioned in Section 1, and as can be seen in Examples 1–2, the designs obtained by solving the geometry and topology optimization problem (15) strongly depend on the initial configuration of the joints. A natural observation is that the designs can be improved if we allow the joints to navigate much wider regions. However, due to imposing the requirement that the joints should not come too close to each other to avoid numerical instabilities, we had to restrict the movement to smaller neighborhoods. The approach now is to solve the problem over and over again, with a view to obtaining improved (lower volume) designs. This is achieved by an iterative procedure that involves many strategies, using an overall process motivated by He and Gilbert (2015).

First, any joints that are entirely connected to bars with cross-sectional areas below \(0.001a_{\max \limits }\), where \(a_{\max \limits }\) is the maximum attained cross-sectional area, are removed. We call these joints inactive nodes. Moreover, any joints connecting only two collinear bars of cross-sectional area \(\geq 0.001a_{\max \limits }\) are melted, i.e., set to vanish, resulting in the collinear bars being merged to form a single bar. Note that there are also other approaches described in the literature that can be used to melt nodes, including that recently proposed by Ohsaki and Hayashi (2017) based on the force density method. We consider two bars as collinear if the angle 𝜃 between them satisfies |π − 𝜃|≤ 0.01. The third strategy is to merge joints that are too close to each other. In our implementation, we perform node merging if the distance between them is less than or equal to 0.25. If the nodes constitute supported or loaded joints, then the nodes are merged to these nodes. Otherwise, the nodes are merged to the average of the coordinates of the merging nodes. By choosing small positive values for ε in (4), we determine joint dependent move limits, that is, balls of radius rj given by:

$$ r_{j}=\min\{k\min\{||\bar{v}_{j}-\bar{v}_{p}||, p\in I\},0.3\}, $$
(24)

where \(\bar {v}_{j}\) and \(\bar {v}_{p}\) are the coordinates of the joints, and I is the set of indices of the joints connected to joint j. Moreover, we use \(k=\frac {1}{3}\) except in Example 7, where we have used \(k=\frac {1}{4}\) due to issues with convergence.

The iterative adaptive optimization procedure stops when there is no significant change in the volume of the optimal designs, and there is no bar with length less than or equal to 0.25.

The overall use of adaptive geometry and topology optimization can also be considered as a post-processing or rationalization technique, since designs with a small number of joints and bars are achievable. The procedure is summarized in Algorithm 2.

figure b

Next, we present examples to demonstrate the benefits of the adaptive geometry and topology optimization technique.

Example 4

We first solve the simple two-dimensional bracket problem shown in Fig. 5a. The design domain has dimensions 3 × 6, the load is (0, − 0.001), and the bound on the compliance is ς = 0.00025. Initially, the radius of the move limits is set to r = 0.3. When we solve the problem without stability constraints, we obtain the optimal design shown in Fig. 5b. Next, when we solve the problem with stability constraints and with the initial configuration of the joints, we obtain the design shown in Fig. 5c. Now, we apply Algorithm 2 to obtain the design shown in Fig. 5d. Successively repeating the process allows us to find the designs shown in Fig. 5e–i.

Fig. 5
figure 5

Example 4: a design domain, boundary conditions, and load; b optimal design without stability constraints; c–i optimal designs with stability constraints obtained in seven stages of the adaptive optimization

Full size image

In this example, we have removed inactive nodes and performed node melting (see Fig. 5d), and merged nodes that are too close (see Fig. 5h). The numerical statistics are presented in Table 4. The CPU time for the last six problem instances are shorter compared with that of the first. This is because many inactive nodes and bars are removed while applying adaptive optimization to the problem shown in Fig. 5c. The designs shown in Fig. 5d–i look qualitatively similar to a design reported by Ferrari and Sigmund (2019), obtained by solving topology optimization with buckling constraints for continuum structures.

Table 4 Example 4 (2D bracket problem): numerical results
Full size table

Example 5

We now apply Algorithm 2 to the optimal L-shaped design shown in Fig. 1e. The results are reported in Fig. 6. Merged nodes can be observed in Fig. 6e and the final design can be observed to have a volume of 5.8911, which is 3.7% lower than that of the initial solution shown in Fig. 6a. The curvature of the outer envelope is even more pronounced, with even more bars connected to the joints in the re-entrant corner. The computational statistics are given in Table 5.

Fig. 6
figure 6

Example 5: optimal designs obtained in six stages of of the adaptive optimization for the L-shaped problem

Full size image
Table 5 Example 5 (L-shaped problem): adaptive numerical results
Full size table

Example 6

Here, we start with the solution of the bridge problem shown in Fig. 2e, where the top four corner nodes in the design space of Fig. 2a are no longer present. The entire process is demonstrated in Table 6 and Fig. 7. No nodes are melted or merged in this case. The final design shown in Fig. 7f has a volume of 4.2135, which is 26.5% lower than the initial solution shown in Fig. 7a. Looking at the final bridge design in Fig. 7f, the semicircular-like side view is similar to the results presented by Kočvara and Zowe (1996) and Achtziger (2007), which are reported for two-dimensional topology and geometry optimization problems. Moreover, the side and top views show that the structure is wider in the middle region and narrower near the two end points, somewhat similar to a form shown in Jiang et al. (2019).

Fig. 7
figure 7

Example 6: optimal designs obtained in the first six stages of the adaptive optimization for the revisited bridge problem

Full size image
Table 6 Example 6 (bridge problem): adaptive numerical results
Full size table

Example 7

We solve the two-span bridge problem, shown in Fig. 8a, supported at span end points at both edges. The design space has dimensions 12 × 1 × 1 and all of the loads have magnitude (0,0,− 0.001), applied simultaneously. The bound on the compliance is set to ς = 0.005. The problem has been solved in stages and the obtained designs are shown in Fig. 8b–f. Once again, we see that the final design (Fig. 8f) has a lower volume than that of the initial solution (Fig. 8b), in this case 14% less. Looking at the final optimal design (Fig. 8f), the side view can be compared to that described in Kočvara and Zowe (1996), where we have two large asymmetric semicircular-like geometries on both sides, and with a smaller symmetric curved structure in the middle. The computational results are reported in Table 7.

Fig. 8
figure 8

Example 7: a design domain, boundary conditions, and loads; b–f optimal designs obtained in five stages of the adaptive optimization

Full size image
Table 7 Example 7 (two-span bridge problem): adaptive numerical results
Full size table

5.3 General comments based on the numerical experiments

As expected, the nonlinear semidefinite programming truss geometry and topology optimization problem given in (15) is not easy to solve. However, in most cases, by setting the radius ri of the move restrictions to a reasonable value, the problems are solvable, mainly because short bars are avoided. The designs can then be improved by applying an adaptive optimization procedure. There were some instances, mainly in the final stages of the adaptive optimization, when problems became harder to solve, due to an inevitable existence of short bars (e.g., see Fig. 6d). However, this can be successfully resolved using a significantly less aggressive update strategy of the barrier parameter and the use of conservative step lengths at the cost of more iterations; see Table 5. As a last comment, it is also very important that the joints are constrained to always remain within the prescribed restricted regions (2) and (3), both initially and in subsequent iterations.

6 Conclusions

We have introduced geometry optimization to an existing truss topology optimization with global stability constraints formulation, posed as a nonlinear semidefinite programming problem. We have demonstrated that these problems can in fact be solved by a standard second-order primal-dual interior point method, despite the associated challenging mathematical properties, such as nonlinearity and nonconvexity. We have also presented an iterative adaptive geometry and topology optimization procedure to improve the quality of the optimized designs. The work is supported by several numerical experiments, showing that allowing the joints to move can lead to reduced material usage and to designs that are more practically relevant.

There seems to be a high degree of similarity between the solutions obtained in successive iterations near the final stage of the presented adaptive optimization procedure. Hence, this suggests that there could be an opportunity to use a warm-start strategy (Weldeyesus and Gondzio 2018; Weldeyesus et al. 2019) to improve the performance of the interior point method further.