Topologically robust CAD model generation for structural optimisation

TL;DR

This paper presents an automated method to convert topology-optimized finite element models into topologically accurate, parametric CAD models using skeletonisation, graph algorithms, and constructive solid geometry, streamlining design workflows.

Contribution

It introduces a fully automated approach for generating topologically accurate CAD models from topology-optimized structures, combining digital topology, graph processing, and CAD regeneration techniques.

Findings

Produces compact, editable CAD models with minimal parameters.

Ensures topological accuracy through skeletonisation and graph algorithms.

Facilitates seamless integration of topology optimization results into CAD systems.

Abstract

Computer-aided design (CAD) models play a crucial role in the design, manufacturing and maintenance of products. Therefore, the mesh-based finite element descriptions common in structural optimisation must be first translated into CAD models. Currently, this can at best be performed semi-manually. We propose a fully automated and topologically accurate approach to synthesise a structurally-sound parametric CAD model from topology optimised finite element models. Our solution is to first convert the topology optimised structure into a spatial frame structure and then to regenerate it in a CAD system using standard constructive solid geometry (CSG) operations. The obtained parametric CAD models are compact, that is, have as few as possible geometric parameters, which makes them ideal for editing and further processing within a CAD system. The critical task of converting the topology optimised structure into an optimal spatial frame structure is accomplished in several steps. We first generate from the topology optimised voxel model a one-voxel-wide voxel chain model using a topology-preserving skeletonisation algorithm from digital topology. The weighted undirected graph defined by the voxel chain model yields a spatial frame structure after processing it with standard graph algorithms. Subsequently, we optimise the cross-sections and layout of the frame members to recover its optimality, which may have been compromised during the conversion process. At last, we generate the obtained frame structure in a CAD system by repeatedly combining primitive solids, like cylinders and spheres, using boolean operations. The resulting solid model is a boundary representation (B-Rep) consisting of trimmed non-uniform rational B-spline (NURBS) curves and surfaces.

Tables

Table 1. Table 1: Efficiency and scaling of the skeletonisation for the quadcopter frame on three different voxel grids.

	Voxel grid # voxels	Voxel model # voxels	Skeleton # voxels	removal steps	time per step	total time
Grid 1	$139\times 16\times 139$	18222	1278	4	$4.43$ s	$17.73$ s
Grid 2	$172\times 19\times 172$	38076	1677	5	$9.33$ s	$46.65$ s
Grid 3	$197\times 21\times 197$	53981	1894	6	$13.17$ s	$78.99$ s

Equations

$$J(\mathbf{\rho })=\mathbf{f}\cdot \mathbf{u}(\mathbf{\rho })$$

$$\mathbf{K}(\mathbf{\rho })\mathbf{u}=\mathbf{f}$$

$$\frac{V(\mathbf{\rho })}{\overline{V}}\le V_f$$

$$0\le \mathbf{\rho }\le 1,$$

$$E({\rho }_i)=E_{\text{min}}+{\rho }_i^p(\overline{E}-E_{\text{min}}),$$

$$H({\mathbf{x}}_i,{\mathbf{x}}_j)=\{\begin{array}{ll}R-\text{dist}({\mathbf{x}}_i,{\mathbf{x}}_j)& \text{if}\text{dist }({\mathbf{x}}_i,{\mathbf{x}}_j)\le R\\ 0& \text{else}\end{array},$$

$${\hat{\rho }}_i=\frac{{\sum }_{j}H({\mathbf{x}}_i,{\mathbf{x}}_j){\rho }_j}{{\sum }_{j}H({\mathbf{x}}_i,{\mathbf{x}}_j)}$$

$$\frac{\partial J(\hat{\mathbf{\rho }})}{\partial {\rho }_i}=\sum _j\mathbf{f}\cdot \frac{\partial \mathbf{u}(\hat{\mathbf{\rho }})}{\partial {\hat{\rho }}_j}\frac{\partial {\hat{\rho }}_j}{\partial {\rho }_i}.$$

$$\frac{\partial \mathbf{u}(\hat{\mathbf{\rho }})}{\partial {\hat{\rho }}_j}=-{\mathbf{K}}^{-1}(\hat{\mathbf{\rho }})\frac{\partial \mathbf{K}(\hat{\mathbf{\rho }})}{\partial {\hat{\rho }}_j}\mathbf{u}(\hat{\mathbf{\rho }})$$

$$\frac{\partial {\hat{\rho }}_j}{\partial {\rho }_i}=\frac{H({\mathbf{x}}_j,{\mathbf{x}}_i)}{{\sum }_{k}H({\mathbf{x}}_j,{\mathbf{x}}_k)}.$$

$$\frac{\partial J(\hat{\mathbf{\rho }})}{\partial {\rho }_i}=-\sum _j\mathbf{u}(\hat{\mathbf{\rho }})\cdot \frac{\partial \mathbf{K}(\hat{\mathbf{\rho }})}{\partial {\hat{\rho }}_j}\mathbf{u}(\hat{\mathbf{\rho }})\frac{H({\mathbf{x}}_j,{\mathbf{x}}_i)}{{\sum }_{k}H({\mathbf{x}}_j,{\mathbf{x}}_k)},$$

$$\frac{\partial V(\hat{\mathbf{\rho }})}{\partial {\rho }_i}=\sum _j{\hat{\rho }}_j\frac{V}{n_e}\frac{H({\mathbf{x}}_j,{\mathbf{x}}_i)}{{\sum }_{k}H({\mathbf{x}}_j,{\mathbf{x}}_k)}.$$

$$J(\mathbf{s})=\mathbf{u}(\mathbf{s})\cdot \mathbf{K}(\mathbf{s})\mathbf{u}(\mathbf{s})=\mathbf{f}\cdot \mathbf{u}(\mathbf{s})$$

$$\mathbf{K}(\mathbf{s})\mathbf{u}(\mathbf{s})=\mathbf{f}$$

$$\frac{V(\mathbf{s})}{\overline{V}}\le V_f$$

$${\mathbf{s}}_l\le \mathbf{s}\le {\mathbf{s}}_u.$$

$${\mathbf{K}}_i={\mathbf{\Lambda }}_i{\mathbf{K}}_i^l{\mathbf{\Lambda }}_i^{\mathsf{T}}$$

$${\mathbf{\Lambda }}_i=\mathrm{diag}({\mathbf{\lambda }}_i,{\mathbf{\lambda }}_i,{\mathbf{\lambda }}_i,{\mathbf{\lambda }}_i).$$

$${\mathbf{\lambda }}_i=\left(\begin{array}{ccc}\cos {\alpha }_i\cos {\beta }_i& -\sin {\alpha }_i& \cos {\alpha }_i\sin {\beta }_i\\ \sin {\alpha }_i\cos {\beta }_i& \cos {\alpha }_i& \sin {\alpha }_i\sin {\beta }_i\\ -\sin {\beta }_i& 0& \cos {\beta }_i\end{array}\right),$$

$$\frac{\partial J(\mathbf{s})}{\partial s_i}=-\mathbf{u}(\mathbf{s})\cdot \frac{\partial \mathbf{K}(\mathbf{s})}{\partial s_i}\mathbf{u}(\mathbf{s})=-\sum _{j=1}^{n_e}{\mathbf{u}}_j(\mathbf{s})\cdot \frac{\partial {\mathbf{K}}_j(\mathbf{s})}{\partial s_i}{\mathbf{u}}_j(\mathbf{s}),$$

$$\frac{\partial {\mathbf{K}}_i(\mathbf{s})}{\partial s_j}={\mathbf{\Lambda }}_i\frac{\partial {\mathbf{K}}_i^l}{\partial s_j}{\mathbf{\Lambda }}_i^{\mathsf{T}},$$

$$\frac{\partial {\mathbf{K}}_i(\mathbf{s})}{\partial s_j}=\frac{\partial {\mathbf{\Lambda }}_i}{\partial s_j}{\mathbf{K}}_i^l{\mathbf{\Lambda }}_i^{\mathsf{T}}+{\mathbf{\Lambda }}_i\frac{\partial {\mathbf{K}}_i^l}{\partial s_j}{\mathbf{\Lambda }}_i^{\mathsf{T}}+{\mathbf{\Lambda }}_i{\mathbf{K}}_i^l\frac{\partial {\mathbf{\Lambda }}_i^{\mathsf{T}}}{\partial s_j}.$$

$$\chi ({\mathcal{M}}_s)=n_0({\mathcal{M}}_s)-n_1({\mathcal{M}}_s)+n_2({\mathcal{M}}_s)-n_3({\mathcal{M}}_s),$$

$$\chi (\partial {\mathcal{M}}_s)=n_0(\partial {\mathcal{M}}_s)-n_1(\partial {\mathcal{M}}_s)+n_2(\partial {\mathcal{M}}_s),$$

$$\chi (\partial {\mathcal{M}}_s)=2\chi ({\mathcal{M}}_s)$$

$$\chi ({\mathcal{M}}_s)=o({\mathcal{M}}_s)-h({\mathcal{M}}_s)+c({\mathcal{M}}_s),$$

$$\chi ({\mathcal{M}}_s)=\sum _q\chi (\mathcal{O}(q))=\sum _q\left(1-\frac{n_1^{(q)}}{2}+\frac{n_2^{(q)}}{4}-\frac{n_3^{(q)}}{8}\right),$$

$$\Delta \chi (\mathcal{M}\left({\mathcal{N}}_{26}(v))\right)$$

$$\Delta h\left(\mathcal{M}({\mathcal{N}}_{26}(v))\right)$$

$$\Delta o\left(\mathcal{M}({\mathcal{N}}_{26}(v))\right)$$

$$\frac{EA}{L}\left(\begin{array}{cc}[r]1& -1\\ -1& 1\end{array}\right),$$

$$\left(\begin{array}{cccccccc}[c]\frac{12E{I}_z}{(1+b_z)L^3}& 0& 0& \frac{6E{I}_z}{(1+b_z)L^2}& -\frac{12E{I}_z}{(1+b_z)L^3}& 0& 0& \frac{6E{I}_z}{(1+b_z)L^2}\\ & \frac{12E{I}_y}{(1+b_y)L^3}& -\frac{6E{I}_y}{(1+b_y)L^2}& 0& 0& -\frac{12E{I}_y}{(1+b_y)L^3}& -\frac{6E{I}_y}{(1+b_y)L^2}& 0\\ & & \frac{(4+b_y)E{I}_y}{(1+b_y)L}& 0& 0& \frac{6E{I}_y}{(1+b_y)L^2}& \frac{(2-b_y)E{I}_y}{(1+b_y)L}& 0\\ & & & \frac{(4+b_z)E{I}_z}{(1+b_z)L}& -\frac{6E{I}_z}{(1+b_z)L^2}& 0& 0& \frac{(2-b_z)E{I}_z}{(1+b_z)L}\\ & & & & \frac{12E{I}_z}{(1+b_z)L^3}& 0& 0& -\frac{6E{I}_z}{(1+b_z)L^2}\\ & & \text{sym.}& & & \frac{12E{I}_y}{(1+b_y)L^3}& \frac{6E{I}_y}{(1+b_y)L^2}& 0\\ & & & & & & \frac{(4+b_y)E{I}_y}{(1+b_y)L}& 0\\ & & & & & & & \frac{(4+b_z)E{I}_z}{(1+b_z)L}\end{array}\right),$$

$$\frac{GJ}{L}\left(\begin{array}{cc}[r]1& -1\\ -1& 1\end{array}\right),$$

Full text

1 Introduction

The application of structural optimisation in industrial design requires the optimised geometries to be converted into CAD models. This is necessary because CAD systems are today an integral part of most product development processes. The prevalent parametric CAD systems are based on boundary representation (B-rep) techniques and trimmed NURBS (non-uniform rational B-splines). Geometries are created using constructive solid geometry (CSG) by recursively combining primitive shapes, such as cylinders, cubes etc. using boolean operations. This recursive process is stored in the form of a CSG tree. Although the input to most structural optimisation is a CAD geometry, their output is usually a finite element mesh of the optimised geometry. The finite element meshes consist of too many elements to be processed and edited in a CAD system. For a geometry to be editable in a CAD system, a compact representation in the form of a CSG tree is needed. The need to edit the optimised part geometry arises when additional geometric features have to be added or when the part is to be combined with other components into a functional product. Moreover, in industrial practice there are usually, in addition to structural efficiency and robustness, many equally important, often not explicitly quantified, design requirements, which require the geometry to be edited after optimisation.

The fully automated process we propose can robustly generate a compact CSG tree representation of topology optimised geometries and convert them to CAD models. The entire workflow is illustrated with the help of the topology optimisation of the rocker arm shown in Figure 1. For topology optimisation we use a standard density based SIMP (solid isotropic material with penalisation) approach on a structured hexahedral finite element mesh referred to as a voxel mesh, see e.g. [1, 2, 3]. The optimised density field splits the voxels after thresholding into the subsets solid and void, which gives a three-dimensional binary image with two grey levels. Subsequently, we extract from the binary image a voxel chain skeleton using a homotopic, i.e. topology-preserving, skeletonisation algorithm from digital topology [4]. The skeleton has the same number of connected components, handles and cavities as the three-dimensional voxel model. Evidently, topology preservation is critical in structural applications because a change of topology may disrupt the load paths discovered during optimisation. The voxel chain model defines a weighted undirected graph which is processed with graph algorithms to obtain a spatial frame structure. Each of the members of the frame structure is a beam and the beams are rigidly connected at the joints. During the conversion of the voxel model to a spatial frame structure the optimality of the structure is usually compromised. This can, however, be recovered after applying a few steps of size and layout optimisation to the spatial frame structure. Size optimisation adjusts the cross-sections of the members and layout optimisation adjusts the coordinates of the joints. The very compact representation which we generate from the original voxel model consists of the connectivity of the frame structure, its joint coordinates and the member cross-sections. These are used to create a binary CSG tree and a solid model of the spatial frame structure in a CAD system. Figures LABEL:fig:rockerArmE and LABEL:fig:rockerArmF show two different solid models generated in a CAD system. In the first model the cylindrical members are connected to spherical joints. In the second model the members are still cylinders, which are, however, smoothly blended at the joints. With the process presented in this paper, the first model is generated in a fully automated fashion while the second requires some user intervention. Note although we used a density based approach for topology optimisation, it is straightforward to apply the proposed approach to the more recent level-set based methods [5, 6] or the historical homogenisation method [7], see also [8] for an overview on optimisation approaches.

There are only a very limited number of structural optimisation approaches which arrive at a geometry in the form of a compact parametric CAD model. In more recent topology optimisation formulations the placement of a finite number of geometric primitives, such as rectangles or cuboids, embedded in a larger design domain is considered [9, 10, 11, 12]. The primitives are allowed to freely move within the design domain usually discretised with a fixed-grid finite element approach. The ensuing optimisation problem is inherently discrete and is converted into a continuous problem by careful regularisation. Although these techniques can provide a compact CAD model they are usually incompatible with the commercially prevalent density based topology optimisation implementations and have so far been applied mostly to 2D problems. There are also other approaches primarily intended for shape optimisation. For instance, it is quite common to use the parameters of a CAD or a reconstructed CAD-like spline model as geometric design variables in shape optimisation, see e.g. [13, 14, 15, 16]. Such techniques are especially appealing in an isogeometric analysis context when the same basis functions are used for geometry description and finite element analysis [17, 18, 19, 20, 21, 22]. The compact spline representation of the optimised geometry can in principle be imported into a CAD system. In the past, there has been some work in topology optimisation on generating a spline representation based on a given optimised finite element geometry. A common approach is to fit the density isocontours obtained with topology optimisation with a spline curve in 2D or a spline surface in 3D [23, 24, 25]. Especially in 3D, these techniques are not robust as they may require the solution of a nonlinear least squares problem and are hard to automate due to the need to manually position the control vertices of the splines. Possibly, therefore, they had no noteworthy impact on commercial software despite being developed around two decades ago. The use of skeletonisation, or thinning, for compact geometric representations has earlier been pioneered in Bremicker et al. [26] for 2D and very recently in Cuillière [27] for 3D. While in [26] the aim was not to arrive at a parametric CAD model, skeletonisation algorithms from digital topology were used to extract a truss structure, which was subsequently size and layout optimised. Homotopic skeletonisation in 3D is substantially more challenging than in 2D, and elimination of possible mechanisms in a spatial truss structure is far from trivial [28]. And so, in our approach, we convert the voxel chain model into a spatial frame structure (with rigidly connected joints) which intrinsically does not exhibit mechanisms. In [27] a curve skeleton is extracted from the surface mesh representing the density isocontour of the optimised geometry using the skeleton extraction method proposed by Au et al.  [29]. This specific extraction technique is rather elaborate as it builds on the successive contraction of a given surface mesh using Laplacian smoothing. The obtained skeleton and estimated cross-sections are used to generate a CAD model with no further postprocessing. It is not taken into account that skeletonisation may have impaired the optimality of the structure. A more manufacturing oriented overview of compact parametric representations for topology optimised geometries can be found in the recent review [30].

An essential component of the proposed CAD model generation process is skeletonisation, which is an active research area with applications in computer graphics, animation and volumetric image processing, see the comprehensive reviews [31, 32, 33, 34, 35]. The curve-skeleton, briefly the skeleton, of a 3D object is closely related to its medial surface, or surface-skeleton. In finite elements medial surfaces are known from mesh generation applications, see e.g. [36]. Informally, the medial surface is the set of all points that have two or more closest points on a 3D object’s surface. That is, a sphere centred on the medial surface touches the object’s surface at two or more points. In contrast to the medial surface, the skeleton of a 3D object has no rigorous definition. It is expected that the skeleton captures the essential topology of the object and is centred, i.e. lies on or close to the medial surface. There are many algorithms for determining the skeleton of an object starting from different types of geometry representations, such as implicit signed distances, polygonal surfaces etc. In topology optimisation, it is expedient to assume that the object is given as a voxelised binary image consisting of solid and void voxels. Digital topology provides a rigorous basis to study the topological properties of binary images consisting of pixels or voxels. An accessible introduction to digital topology can be found in the review [37] and the textbook [38]. Homotopic skeletonisation algorithms based on digital topology can extract a one-voxel-wide skeleton with the same topology as the binary image, i.e. the same number of connected components, handles and cavities. See the excellent reviews [37, 39] for early papers and [33, 35] for more recent papers on skeletonisation. Most of these methods start from the domain boundaries and proceed by iteratively removing voxels one by one that do not alter the topology of the object. When the algorithms terminate, the entire structure is represented by a network of one-voxel-wide chains. It is sufficient to inspect only a voxel’s immediate neighbourhood consisting of 3$\times$3$\times$3 voxels to decide whether that voxel can be removed. The possible states of such a small cluster can effectively be encoded in a look-up table [40, 4]. Crucially, the skeletonisation process does not rely on any floating-point operations, which makes it exceedingly robust. Recent research on skeletonisation focuses on parallelisation [41, 42] and improving the memory usage [43] of the algorithms. Skeletonisation is however far less computing intensive than finite element analysis so that advanced skeletonisation algorithms are of limited relevance for the present work. Our approach is based on Lee et al. [4] which is one of the first provably homotopy preserving skeletonisation algorithms. Its implementation is also discussed in the more recent papers [44, 45].

The outline of this paper is as follows. In Section 2, the standard density based topology optimisation and the size and layout optimisation of spatial frames are briefly reviewed. Subsequently, in Section 3 relevant aspects of digital topology are introduced and the implemented specific skeletonisation algorithm is described. The robustness and runtime of our implementation are studied with a relatively complex quadcopter frame geometry. The postprocessing of the obtained voxel chain model with graph algorithms and the generation of first a frame structure and then a CAD model are discussed in Section 4. The entire optimisation and geometry conversion process is summarised in Section 5. The application of the proposed approach to a standard 2D benchmark example from topology optimisation and more complex 3D examples is illustrated in Section 6. We study in particular the change of the compliance cost function during the conversion of the optimised voxel model to a CAD model.

2 Review of structural optimisation

In this Section, we briefly review the standard density based topology optimisation using the SIMP method and the size and layout optimisation of spatial frame structures. In all cases, the cost function is the compliance, and the discussion is restricted to aspects relevant to this paper. For further details see, e.g., the monographs [3, 46].

2.1 Topology optimisation of solids

The topology optimisation problem for finite element discretised solids reads

[TABLE]

where $J(\bm{\rho})$ is the compliance cost function, $\bm{\rho}$ is the vector of relative element densities, $\bm{u}$ is the displacement vector, $\bm{K}$ is the global stiffness matrix, $\bm{f}$ is the global external force vector, $V(\bm{\rho})$ is the material volume, $\overline{V}$ is the design domain volume and $V_{f}$ is the scalar prescribed volume fraction. The relative density of each element is constrained to be $0\leq\rho_{i}\leq 1$ with $i\in\{1,\dotsc,n_{e}\}$ and the number of elements is $n_{e}$.

As usual, the global stiffness matrix $\bm{K}$ and vector $\bm{f}$ are assembled from the $n_{e}$ element contributions $\bm{K}_{i}$ and $\bm{f}_{i}$ in the mesh. In this paper, we discretise the design domain always with a structured grid and hexahedral linear elements. In each element $i$ the material is isotropic and homogeneous, and the Young’s modulus $E$ is penalised depending on the relative density $\rho_{i}$ with

[TABLE]

where $\overline{E}$ is the prescribed Young’s modulus of the solid material, $p$ and $E_{\text{min}}$ are two algorithmic parameters. The penalisation parameter $p\geq 3$ ensures that elements with densities close to $\rho_{i}=0$ (void) and $\rho_{i}=1$ (solid) are preferred. The small Young’s modulus $E_{\text{min}}\approx 10^{-9}$ of the void material prevents ill-conditioning of the global stiffness matrix when $\rho_{i}=0$. Each element stiffness matrix $\bm{K}_{i}$ is computed using the corresponding Young’s modulus $E(\rho_{i})$ with the relative density $\rho_{i}$, which is constant within an element.

Furthermore, the topology optimisation problem (1) is regularised by filtering the element densities $\rho_{i}$ to prevent checker-board instabilities and mesh dependency of the solution. This is accomplished by convolving the element densities with the kernel function

[TABLE]

which depends on the coordinates of the centroids $\bm{x}_{i}$ and $\bm{x}_{j}$ of the elements $i$ and $j$, and the prescribed filter length $R$. With $H(\bm{x}_{i},\,\bm{x}_{j})$ at hand the filtered densities are given by

[TABLE]

assuming that all elements have the same volume.

For gradient-based optimisation the derivatives of the cost and constraint functions in (1) with respect to the relative densities $\rho_{i}$ are needed. In the following the relative densities $\rho_{i}$ in the topology optimisation problem (1) are replaced with the filtered relative densities $\hat{\rho}_{i}$. The derivative, or sensitivity, of the compliance cost function (1a) reads

[TABLE]

Differentiating and rearranging the equilibrium constraint (1b) gives

[TABLE]

and the derivative of the filtered densities (4) is

[TABLE]

Introducing both derivatives into (5) yields

[TABLE]

where the stiffness matrix derivative is to be assembled from element contributions considering the penalised Young’s modulus (2).

The derivative of the volume constraint (1c) is obtained analogously with

[TABLE]

It is assumed here that all the elements are of the same size, namely the design volume $\overline{V}$ divided by the total number of elements $n_{e}$.

Finally, with the obtained derivative of the compliance cost function (8) and the volume constraint (9) the optimised density distribution is determined iteratively with the optimality criteria method, see e.g. [3].

2.2 Size and layout optimisation of frames

We consider spatial frame structures consisting of straight beam members connected by joints that can transfer forces and moments. The members can deform by stretching, bending and torsion and are modelled as classical Timoshenko beams, see A. Without loss of generality, size and layout optimisation can be posed as an iterative sequential optimisation problem. In each optimisation step, either the size, i.e. cross-section area, or layout, i.e. joint coordinates, of the members is optimised.

Both the size and layout optimisation problems have the same structure as topology optimisation (1), namely,

[TABLE]

The vector of design variables $\bm{s}$ refers either to cross-section areas in case of size optimisation or joint coordinates in case of layout optimisation. The two bounds $\bm{s}_{l}$ and $\bm{s}_{u}$ have different interpretations in the two cases.

The frame consists of $n_{e}$ beam elements with the element stiffness matrices $\bm{K}_{i}$ and vectors $\bm{f}_{i}$. Each stiffness matrix $\bm{K}_{i}$ is obtained from a local stiffness matrix $\bm{K}_{i}^{l}$ formulated in a coordinate system attached to the element with the index $i$. The local matrix $\bm{K}_{i}^{l}\in\mathbb{R}^{12\times 12}$ with the displacements and rotations of the two end nodes of the beam element as degrees of freedom is given in A. In the local $x^{l}_{i}$, $y^{l}_{i}$ and $z^{l}_{i}$ coordinate system the beam axis is assumed to be aligned with the $x^{l}_{i}$ axis. The matrix $K_{i}^{l}$ is transformed to the global $x$, $y$ and $z$ coordinate system according to

[TABLE]

with the block-diagonal transformation matrix

[TABLE]

The rotation matrix $\bm{\lambda}_{i}\in\mathbb{R}^{3\times 3}$ can be chosen with

[TABLE]

which is composed of the two elemental rotations $\alpha_{i}$ around the $z_{i}^{l}$ axis and $\beta_{i}$ around $y_{i}^{l}$ axis.

Similar as in topology optimisation, the derivative, or sensitivity, of the compliance cost function (10a) reads

[TABLE]

where $\bm{u}_{j}$ is the vector of nodal displacements and rotations of the element $j$. In size optimisation the derivative of the element stiffness matrix (11) is given by

[TABLE]

whereas in layout optimisation it is

[TABLE]

3 Digital topology and skeletonisation

The input to skeletonisation is a structured hexahedral finite element grid with all the vertices within the domain having eight incident hexahedra. For digital topology only the connectivity of the mesh is relevant whereas the coordinates of the grid nodes can be arbitrary so that in the following hexahedron and voxel are used interchangeably. The hexahedral grid encodes the binary voxel model of the topology optimised structure. The binary voxel model is obtained by denoting voxels with a density $\rho_{i}$ above the threshold $\eta$ as solid and the remaining ones as empty (void). Hence, the set of all voxels $v\in\mathcal{V}=\mathcal{V}_{s}\cup\mathcal{V}_{e}$ is composed of the solid voxels $\mathcal{V}_{s}$ and empty voxels $\mathcal{V}_{e}$. Skeletonisation aims to reassign voxels one by one from $\mathcal{V}_{s}$ to $\mathcal{V}_{e}$ until the solid is represented by a network of one-voxel-wide chains. The decision about which voxels to reassign is guided by digital topology and proceeds starting from the voxels at the border between $\mathcal{V}_{s}$ and $\mathcal{V}_{e}$.

3.1 Review of digital topology

To retain the topology of the voxel model, its Euler characteristic is considered. The solid voxels in $\mathcal{V}_{s}$ and empty voxels in $\mathcal{V}_{e}$ yield two hexahedral volume meshes, or more abstractly cell complexes, which will be in the following denoted with $\mathcal{M}_{s}=\mathcal{M}(\mathcal{V}_{s})$ and $\mathcal{M}_{e}=\mathcal{M}(\mathcal{V}_{e})$. The volume mesh of a voxel set consists of the corresponding set of the vertices, edges, faces and voxels. It is also useful to define quadrilateral surface meshes $\partial\mathcal{M}_{s}$ and $\partial\mathcal{M}_{e}$ formed by the boundaries of $\mathcal{M}_{s}$ and $\mathcal{M}_{e}$, respectively. The Euler characteristic of the voxel model can be determined either from the volume mesh $\mathcal{M}_{s}$ with

[TABLE]

or from the surface mesh $\partial\mathcal{M}_{s}$ with

[TABLE]

where $n_{0}$ denotes the number of vertices, $n_{1}$ the number of edges, $n_{2}$ the number of the faces and $n_{3}$ the number of hexahedrons in $\mathcal{M}_{s}$ or $\partial\mathcal{M}_{s}$. In Figure 2 the Euler characteristics of three voxel meshes are given. Notice that there is also the relation

[TABLE]

between the Euler characteristics of the surface and volume meshes. The Euler characteristic of a voxel model is related to the total number of separate objects (connected components) $o(\mathcal{M}_{s})$, handles (tunnels) $h(\mathcal{M}_{s})$ and cavities $c(\mathcal{M}_{s})$ in the entire model according to the formula

[TABLE]

see, e.g., textbooks [47, 48] on topology. For instance, the Euler characteristic $\chi(\mathcal{M}_{s})$ of a sphere is $1$, of a torus is [math], of a double torus is $-1$ and so on, compare also Figure 2. Hence, the total number of different entities in the mesh and the described shape are intrinsically connected.

In determining the Euler characteristic of a voxel model, two different types of neighbourship definitions are possible, see Figure 3. These definitions define how pairs of touching voxels such as shown in Figure 3c and  3d are treated, whether connected and forming one surface or disconnected and forming two surfaces. That is, the number of entities in (17) and (18) depends on the chosen definition of the neighbourship. For a voxel $v$ its 6-neighbours consist of all the cells which share a face with $v$ and its 26-neighbours consist of all the cells which share a common vertex with $v$. The voxels in the 6-neighbourhood of $v$ are denoted with $\mathcal{N}_{6}(v)$ and the ones in the 26-neighbourhood with $\mathcal{N}_{26}(v)$. In our implementation, two solid voxels are considered to be neighbours when they are 26-neighbours. In contrast, two void voxels are considered to be neighbours when they are 6-neighbours. As pointed out in Kong et al. [37], it is essential not to use the same neighbourship definitions for both the solid and void voxels. The use of one single neighbourhood definition leads to ambiguities, for instance, in the planar case to the violation of the Jordan curve theorem. The Jordan curve theorem states that a simple closed curve divides the plane into an interior and an exterior region.

The Euler characteristic (17) and (18) for the entire voxel model can be determined by looping over all vertices in the grid and only considering in each step the voxels attached to one vertex. The 2$\times$2$\times$2 voxels attached to a vertex $q$ are defined as its octant $\mathcal{O}(q)$.111Note that it is also possible to associate the octants with the voxels instead of the vertices in the grid. Both viewpoints are equivalent [4]. There are as many octants as vertices in the volume grid (assuming that the voxel domain is padded with empty ghost voxels). The octants are overlapping and the edges, faces and voxels in the hexahedral mesh $\mathcal{M}_{s}$ appear in several octants. Hence, in computing the Euler characteristic (17) by iterating over the vertices the contribution of each octant has to be suitably weighted such that

[TABLE]

where $n_{1}^{(q)}$ is the number of edges and $n_{2}^{(q)}$ is the number of faces neighbouring to the $n_{3}^{(q)}$ solid voxels in the mesh $\mathcal{M}(\mathcal{O}(q))\setminus\partial\mathcal{M}(\mathcal{O}(q))$, see Figure 4. The contribution of an octant $\chi(\mathcal{O}(q))$ with different solid-empty voxel states can be precomputed and stored in a look-up table. The eight voxels in an octant have $2^{8}=256$ possible solid-empty states and considering symmetries this reduces to only 22 distinct cases. Tables with contributions of octants to the Euler characteristic can be found, e.g., in [4, 44, 45]

3.2 Skeletonisation by successive voxel removal

The Euler characteristic $\chi(\mathcal{M}_{s})$ is critical in determining the voxels at the boundary of the solid mesh $\mathcal{M}_{s}$ that can be deleted without changing the topology of the voxel model. A solid border voxel is defined as a solid voxel with at least one void voxel amongst its 6-neighbours. There is a large amount of research on voxels, also referred to as simple points, which can be removed without changing topology [37]. As evident from (20) simply conserving the Euler characteristic $\chi(\mathcal{M}_{s})$ of the entire or a portion of the mesh does not ensure that the topology is not changed, i.e. the number of objects $o(\mathcal{M}_{s})$, handles $h(\mathcal{M}_{s})$ and cavities $c(\mathcal{M}_{s})$ all remain the same. For a solid border voxel $v$ to be classified as simple its removal must not change, in addition to the Euler characteristic, the number of objects and handles for both $\mathcal{M}_{s}$ and $\mathcal{M}_{e}$ [37]. According to Lee et al. [4] these conditions can be checked by examining the state of the voxels in eight octants overlapping voxel $v$, or the voxels in $v$’s 26-neighbourhood $\mathcal{N}_{26}(v)$. That is, a border voxel $v$ is a simple point if and only if

[TABLE]

where $\Delta$ denotes the change of the respective quantity with and without voxel $v$ present. It is straightforward to determine the Euler characteristic (22a) and it is relatively easy to determine the number of objects (22c) in $\mathcal{N}_{26}(v)$. The change of the Euler characteristic is determined by summing the tabulated contributions of the eight octants belonging to the eight corners of the voxel $v$ from [4], see  (21) and Figure 5. To determine the number of objects (22c), we generate an undirected graph with the solid voxels as nodes and introducing edges between voxels that are 26-neighbours. Subsequently, the number of connected components is determined with a depth first search (DFS) algorithm [49].

The skeletonisation proceeds by removing in each step one layer of border voxels which are simple points according to (22), see Algorithm 1 in B. The removed solid voxels are reclassified as void. One removal step is split into six sub-steps and in each sub-step only the voxels approaching from one of the six grid directions are removed. The geometry of the obtained skeleton depends somewhat on the sequencing of these sub-steps, which is not critical for our purpose. It is straightforward to tag some voxels as non-removable irrespective of (22). In skeletonising topology optimisation results, for instance, the voxels at Dirichlet and non-zero Neumann boundaries are tagged as non-removable. Also, end points with only one voxel in their 26-neighbourhood have to be tagged as non-removable, in order to avoid a complete elimination of voxel chains. The skeletonisation terminates when none of the remaining voxels is removable without violating (22). The final result is a curve skeleton consisting of a network of one-voxel-wide voxel chains connected by joint voxels. We refer to this curve skeleton as the voxel chain skeleton.

3.3 Illustrative example and timing

We consider the skeletonisation of the quadcopter frame shown in Figure 6a obtained from GrabCAD. The frame has a non-trivial topology and has been designed in a CAD system. To generate a voxel model the implicit, or level set, representation of the frame on a uniform voxel grid is computed first. This is achieved by determining the distance of each voxel in the grid to the frame surface with the open source openVDB library [50]. In doing so, the frame geometry is approximated with the STL mesh exported from the CAD system depicted in Figure 6b. After thresholding the voxel grid the voxel model of the frame in Figure 6c is obtained. Its skeletonisation with the introduced digital topology algorithm yields the voxel chain skeleton in Figure 6d. It can be visually confirmed that the skeleton appears to have the same topology as the voxel model.

To investigate the efficiency and scaling of our C++ implementation of the introduced skeletonisation algorithm, we consider three different voxel grids for computing the implicit representation and skeletonisation. The size of the three grids, the number of voxels in the different models and the runtimes are given in Table 1. The STL mesh has in all cases 1086791 triangles. The number of solid voxels in the voxel model and the number of removal steps increase with decreasing voxel size because more and more of them cover the frame. The number of skeleton voxels increases because the smaller voxels can capture the topology of the frame better. The reported extremely short runtimes include only skeletonisation (no implicitisation) and confirm the efficiency of the skeletonisation algorithm. Moreover, notice that the average time per removal step is approximately linear with respect to the number of the voxels in the voxel grid. The experiments were performed on a Macbook Pro with an Intel Core CPU i7-4750HQ $@$ 2.0GHz and 16GB RAM.

4 Model extraction and conversion

The voxel chain skeleton is used to define a structural frame model. The skeleton provides both the connectivity and the geometry of the frame. Although it is feasible to obtain also the member cross-sections from the voxel model, in our implementation cross sections are obtained from size optimisation of the frame, as will be discussed in Section 5. In converting the voxel chain skeleton to a frame we express it as a weighted undirected graph and make use of its incidence matrix. The frame model provides a very compact representation of the optimised structure, which can be used to generate a parametric solid model using the scripting interface or API (application programming interface) of a CAD system.

4.1 Graph model

The voxel chain obtained from skeletonisation contains two different types of voxels, see Figure 7. We refer to the voxels with only two voxels in their 26-neighbourhood as regular voxels. The remaining voxels are the joint voxels, which have other than two voxels in their 26-neighbourhood. As illustrated in Figure 7b, it is not expedient to deduce from every joint voxel a joint for the frame model because this will lead to acutely short members and impractical joint designs.

To robustly merge joint voxels in close proximity to a single joint we resort to a weighted undirected graph representation of the voxel chain skeleton. The nodes of the graph are the joint voxels, and the edges are the connections between the joint voxels. Each of the nodes has an associated coordinate vector which is initially equal to the centroid of the corresponding joint voxel. The edge weights are proportional to the number of voxels between the two joints attached to an edge. The graph model is generated with the Algorithm 2 given in B. In a first step, all the solid voxels in the voxel model are categorised either as regular or joint. Subsequently, starting from the joint voxels we determine the voxel chains and their lengths by marching along the $\mathcal{N}_{26}(v)$ neighbours until a joint voxel is reached. The obtained graph model is stored in the form of an incidence matrix. In Figure 8a and 8d the graph model and its incidence matrix for the voxel chain skeleton in Figure 7b are given. The rows of the incidence matrix correspond to the nodes of the graph and the columns to its edges. Each edge has two identical non-zero entries corresponding to the length of the edge in voxel units.

The joint voxels in close proximity of each other are connected by very short edges. The corresponding graph nodes are merged by successively collapsing the edges of the graph. The edge collapse is implemented with the help of the incidence matrix. Figure 8 illustrates how the incidence matrix evolves during the collapse of the edge $e_{3}$ and the merging of its attached nodes $v_{3}$ and $v_{4}$. The coordinate vector associated to the node $v_{3}$ is the average of the node coordinate vectors before merging. After collapsing, the duplicate edge $e_{5}$ is detected and merged with $e_{4}$. The weight of $e_{4}$ is updated to be the closest integer to the mean of the weights of $e_{4}$ and $e_{5}$.

The voxel chain model usually contains branches which lead to structural frame members with zero stress. These are the chains directly connected to a zero Neumann boundary. Their pruning is again implemented with the incidence matrix. As illustrated in Figure 9, after identifying from the incidence matrix the nodes with only one attached edge, nodes and edges are deleted by simply removing the corresponding rows and columns.

4.2 Structural frame and CAD models

The graph model supplemented by the joint coordinates and member cross-sections provides sufficient information to generate a compact structural frame model. We assume that there is a one-to-one correspondence between nodes and edges of the graph model and joints and members of the frame model. Besides, it is assumed that members are only subjected to end forces and moments, i.e. no distributed loads, so that they can be straight and have uniform cross-sections along their lengths. Note that these assumptions are fully justified if there is no distributed loading in topology optimisation. Furthermore, for simplicity, we assume in the following that all cross-sections are circular. These restrictions can be mostly relaxed if necessary.

It is straightforward to create a parametric solid CAD model from the structural frame model. To this end, we use either the commercial CAD system Rhinoceros (Rhino) or the open source FreeCAD. However, the proposed approach can be realised in any CAD system which provides a scripting interface or API (application programming interface). In the solid model, the members are represented by cylinders and the joints by spheres, see Figure 10. They are combined with boolean operations. The radius of the sphere at a joint is chosen slightly larger than the largest attached member radius. In this paper, we use a factor of $1.05$. A binary CSG tree of the frame model is generated by starting with one of the members and by successively adding spheres and cylinders to the model. Different variations of this approach can be devised and implemented. For instance, it suggests itself to connect two of the most stressed members continuously across a joint or to add fillets to reduce stress concentrations. See also [51, 52] in which a somewhat similar approach was used for generating CAD models for optimised truss structures. After the model is generated it can be edited and postprocessed in the CAD system and be exported, for instance, in STL format for machining or IGES or STEP formats for mesh generation software.

5 Overall optimisation and conversion process

In this Section, we review the sequence of steps from the definition of the topology optimisation problem to obtaining the structurally-sound parametric CAD geometry. We provide additional details for each of the steps and focus on the interplay between the different steps. Each step corresponds to one Section in this paper.

Step 1: Topology optimisation

We assume that the finite element mesh for the optimisation problem (1) is a structured hexahedral grid. In most optimisation problems, the design domain is a parallelepiped which can be discretised with a structured grid. Other more complex design domains can be considered by computing their implicit, or level set, representation and embedding them in a structured hexahedral grid, see the quadcopter frame example in Figure 6. From the outset, the voxels outside the design domain are chosen as void by choosing their Young’s moduli with $E_{\text{min}}$.

Step 2: Skeletonisation

The input to the homotopic skeletonisation algorithm is a binary image defined on a hexahedral structured grid. The binary image is obtained by thresholding the topology optimised geometry. The threshold value $\eta$ is chosen such that the prescribed material volume constraint (1c) in topology optimisation is retained. The output of the skeletonisation is a curve skeleton in the form of a voxel chain with the same topology as the topology optimised geometry.

Step 3: Structural frame model generation

The topology and joint coordinates of the frame model are extracted from the voxel chain model. It is assumed that all the members are straight, have a circular cross-section with the same diameter and their total volume is equal to the prescribed material volume $V_{f}\overline{V}$ in topology optimisation. It is possible to obtain frame models with more complex member geometries and non-uniform cross-sections. This appears to be, however, usually not desired from an ease and cost of manufacturability viewpoint.

Step 4: Size and layout optimisation

The structural frame model extracted from the skeleton is usually suboptimal. That is, the compliance of the structural frame is larger than the compliance of the topology optimised geometry. To recover the optimality of the frame structure, we apply several steps of sequential size and layout optimisation (10). In size optimisation the diameters $d_{i}$ of the circular member cross-sections and in layout optimisation the coordinates $\bm{x}_{i}$ of the joints are updated. Both optimisation problems are solved with the SQP (sequential quadratic programming) method [53]. Throughout optimisation the volume of the frame is constrained to be equal to the prescribed material volume in topology optimisation $V_{f}\overline{V}$. It is straightforward to consider additional constraints pertaining to the member cross-sections, the positions of the joints, or positions and orientations of the members. The first step is always size optimisation, which is followed by as many as necessary alternating layout and size optimisation steps until the compliance cost function is converged. If the length of any member reduces to zero during shape optimisation, the member is removed, its end nodes are merged, and the iteration continues.

Step 5: CAD model generation

A compact CAD model of the structural frame can be generated essentially in any parametric CAD system using a fully automated process. The members are represented by cylinders and the joints by spheres which are combined by boolean operations. The underlying binary CSG tree representation makes it easy to edit further and to refine the optimised design.

6 Examples

We consider three examples of increasing complexity to demonstrate the application of the proposed approach. The minimised cost function is the compliance $J(\hat{\bm{\rho}})$ in topology optimisation and $J(\bm{d},\,\bm{x})$ in size and layout optimisation, where $\hat{\bm{\rho}}$ are the filtered element densities, $\bm{d}$ are the member diameters, and $\bm{x}$ are the joint coordinates. In all examples, the Young’s modulus and Poisson’s ratio of the solid material are $\overline{E}=2.1\cdot 10^{5}$ and $\nu=0.3$. In topology optimisation the penalisation factor is chosen as $p=3$ and the minimum Young’s modulus as $E_{min}=10^{-9}$. During layout optimisation, a member shorter than $1/20$ of the total length of its immediate neighbouring members is considered as too short, and its two end nodes are merged.

6.1 Cantilever

Cantilevered plates as shown in Figure 11 are one of the most widely studied benchmark examples in topology optimisation, see e.g. [3]. The length, height and width of the chosen design domain are $150\times 50\times 4$. The left face of the domain is clamped while all the other faces are free, and a distributed force with a total value of $F=100$ is applied along the centre of the right face. The finite element discretisation consists of $150\times 50\times 4$ linear hexahedral elements. The maximum volume fraction of the optimised structure is prescribed to be $V_{f}=0.3$ and the density filter radius is chosen as $R=3$. Figure LABEL:fig:cantilever_top depicts the optimised cantilever with only the voxels above a relative density $\eta=0.3$. The compliance of this voxel model is $J(\hat{\bm{\rho}})=3.40900$, which has been determined after topology optimisation by resetting $p=1$ in (2).

As a first step in obtaining the structural frame model the homotopic skeletonisation algorithm is employed. The skeletonisation yields after $5\times 6$ voxel removal steps the voxel chain skeleton shown in Figure LABEL:fig:cantilever_skeleton. As visually apparent, the voxel and voxel chain models have the same topology. The skeleton is converted into the structure depicted in Figure LABEL:fig:cantilever_frame by first identifying the 16 joints and then connecting them with 21 members. As evident from Figure LABEL:fig:cantilever_frame during the conversion to a frame model the optimality of the voxel model is compromised; note, for instance, the non-straight top and bottom members. Optimality is recovered by updating the joint positions, i.e. layout optimisation, and the cross-sectional areas, i.e. size optimisation, see Figure LABEL:fig:cantilever_result. Initially, all members are assumed to have the same diameter of  $d=4.547$, giving the prescribed total volume of $V=0.3\overline{V}$. In determining the volume only the member cross-sections and member lengths between the joints are considered.

According to Figure 14b, the frame consisting of beams with uniform diameter has a compliance $J(\bm{d},\bm{x})=4.59841$, which is larger than the voxel model. The stretch and bending strain energies of the frame are 2.02202 and 0.247108, respectively. The minimum and maximum axial stresses are 0 and 20.9274. The average axial stress is 8.04859 and its standard deviation is 6.26102. After the first size optimisation step the compliance is reduced to $J(\bm{d},\,\bm{x})=3.59948$ and the subsequent layout optimisation step to $J(\bm{d},\,\bm{x})=2.95264$. Several more steps of size and layout optimisation do not lead to a significant reduction in the compliance. Notice that the obtained final compliance $J(\bm{d},\,\bm{x})=2.90323$ is significantly lower than the compliance $J(\hat{\bm{\rho}})=3.40900$ of the topology optimised voxel model. The stretch and bending strain energies of the final frame are 1.42114 and 0.0190822, respectively. The minimum and maximum axial stresses are 7.00854 and 8.46234. The average axial stress is 8.09638 and its standard deviation is 0.44709. In the final optimised cantilever in Figure LABEL:fig:cantilever_result the axial stresses become more evenly distributed and the bending stresses are significantly reduced. The final optimised cantilever consists of 11 joints and 16 members. Finally, in Figure 15 the solid CAD model of the frame and a faceted triangular STL mesh exported from FreeCAD are shown. The lines in the CAD model show the layout of the trimmed NURBS patches.

6.2 Pipe bracket

As a structure with a truly three-dimensional load path, the pipe bracket shown in Figure 16 is considered. The design domain has a length, height and width of $120\times 40\times 60$ and contains two openings with each of radius of 18 for two pipes passing through the domain. Within the openings each pipe is supported at four points, applying at each support a vertical force of $F=100$. The four vertical outer edges of the design domain are chosen as fixed. The finite element discretisation consists of $120\times 40\times 60$ linear hexahedral elements. To represent the two openings in the finite element model the Young’s modulus of the elements within the void region is prescribed with $E_{\text{min}}$.

The maximum volume fraction of the optimised structure is prescribed to be $V_{f}=0.1$ and the density filter radius is chosen as $R=3$. Figure LABEL:fig:pipe_top depicts the optimised bracket with only the voxels above a relative density $\eta=0.5$. The compliance of the optimised voxel model is $J(\hat{\bm{\rho}})=2.47602$. As easily recognisable from Figure LABEL:fig:pipe_top, the optimised structure is a combination of an arch-like and a cable-like structure with vertical tension members between the two openings.

The skeletonisation algorithm yields after $6\times 6$ voxel removal steps the voxel chain skeleton in Figure LABEL:fig:pipe_skeleton. The skeleton is converted into the structure depicted in Figure LABEL:fig:pipe_frame by first identifying the 42 joints and then connecting them with 50 straight members. Alternatively, the 50 members can be represented with curved Bézier segments, which would lead to a structural frame with curved beam members. This has not been pursued here because such members may lead to an increase in manufacturing costs and are usually not preferred. Initially, all frame members are assumed to have the same diameter $d=6.296$, giving a total volume of $V=0.1\overline{V}$.

According to Figure 19b, the frame with beams of uniform diameter has a compliance $J(\bm{d},\,\bm{x})=4.69409$, which is significantly higher than $J(\hat{\bm{\rho}})=2.47602$ of the voxel model. After the first size optimisation step the compliance is reduced to $J(\bm{d},\,\bm{x})=3.03554$ and the subsequent layout optimisation step to $J(\bm{d},\,\bm{x})=2.26858$. In layout optimisation the beams are constrained not to move inside the volume occupied by the two cylindrical pipes. The intersection between a beam and a cylinder is determined by integrating the (positive) signed distance of the cylinder along the beam. It is straightforward to differentiate this constraint integral with respect to nodal coordinates of the beam. The intersection constraints for all the beams are aggregated with the Kreisselmeier-Steinhauser (K-S) function [54]. Several more steps of size and layout optimisation do not lead to a significant reduction in the cost function. Notice that the obtained final compliance $J(\bm{d},\,\bm{x})=2.26683$ is lower than the compliance $J(\hat{\bm{\rho}})=2.47602$ of the topology optimised voxel model. The final optimised frame structure depicted in Figure LABEL:fig:pipe_result consists of 36 joints and 44 members. In Figure 20 the solid CAD model of the pipe bracket and a faceted triangular STL mesh exported from FreeCAD are shown.

6.3 Frame-supported plate

In structural design it is very common to combine a plate- or shell-like skin structure with a frame-like support structure. In such composite structures the skin can be modelled as a standard plate or shell structure, and topology optimisation is applied only in the remaining part of the design domain. The design domain shown in Figure 21 combines a thin plate-like skin with a solid bottom part to be optimised. The entire domain is of size $120\times 80\times 20$ and contains a small opening of size $40\times 40\times 6$. The plate is of thickness $3$ and is subjected to a uniform pressure load of $1$. The four corners of the design domain are vertically supported. Due to symmetry only a quarter of the design domain is considered in topology optimisation. The finite element discretisation consists of $60\times 40\times 23$ linear hexahedral finite elements. Across the height of the design domain there are $6$ elements with a size of $1\times 1\times 0.5$ in the plate-like top part and $17$ unit sized elements in the bottom part. To represent the small opening in the finite element model the Young’s modulus of the relevant elements is prescribed as $E_{\text{min}}$.

In topology optimisation, the maximum volume fraction is prescribed to be $V_{f}=0.25$ and the density filter radius is chosen as $R=1.5$. The optimised structure with only the voxels above a relative density $\eta=0.55$ is shown in Figure LABEL:fig:deck_top and its voxel chain skeleton in Figure LABEL:fig:deck_skeleton. Although the voxels in the top part of the design domain are present during optimisation and skeletonisation, they are tagged as non-removable. The voxel chain skeleton is converted to the frame structure shown in Figure LABEL:fig:deck_frame. The structural model also contains a top plate not shown in Figure 25. The plate is modelled as a Kirchhoff–Love plate and is discretised with $4\times 4$ quadrilateral elements and Catmull-Clark subdivision basis functions. The frame has 124 joints and 136 members. Initially, all the members are assumed to have the same diameter of $d=5.010$ which gives a total frame volume of $V=4800$. The coupling between the plate and the frame is achieved with Lagrange multipliers as described in [55]. The joints between the plate and frame can transfer forces but not moments.

According to Figure LABEL:fig:deck_conv_frame, the compliance of the frame-supported plate with members of uniform diameter is $J(\bm{d},\,\bm{x})=1210.01$, which reduces after several sequential size and layout optimisation steps to $J(\bm{d},\,\bm{x})=322.967$. During the layout optimisation, the positions of the joints between the frame and plate are fixed. For comparison, the compliance of the topology optimised voxel model was $J(\hat{\bm{\rho}})=364.596$. The final optimised frame structure contains 56 joints and 64 members, see Figure LABEL:fig:deck_result. Its solid CAD model and faceted triangular STL mesh exported from FreeCAD are depicted in Figure 25.

7 Conclusions

The automated conversion of topology optimised geometries into compact parametric CAD models is a missing link in the wide-scale industrial adoption of optimisation. We introduced a fully automated workflow for synthesising a structurally-sound parametric CAD model from a voxelised solid-void binary image obtained with topology optimisation. Currently, there is no commercial design system, such as Altair Inspire, Abaqus CAE or similar, which offers such a fully automated CAD model generation capability. Homotopic skeletonisation of three-dimensional images is a key component of our proposed workflow, which is an extensively well-studied topic in digital topology. The resulting voxel chain skeleton is interpreted as a weighted undirected graph and processed with standard graph algorithms to yield a frame model. Both the skeletonisation and graph algorithms operate on topology rather than geometry and are not subject to floating-point errors, which makes them exceedingly robust. The so-obtained structural frame models are not optimal because the skeletonisation and graph algorithms do not take into account any mechanical considerations. However, as shown in the numerical examples, after several alternating steps of size and layout optimisation, the structural frames can achieve a cost function value similar to the original topology optimised geometry. The frame model is converted into a compact CAD solid model by recursively combining primitive shapes using boolean operations. Although we used in this paper only cylinders and spheres as shape primitives, it is straightforward to use other primitives. There is also scope in improving the CAD solid model itself, expressly, by using the recently proposed quadratic of revolution (quador) representation for lattices and frames [56, 57, 58]. A quador CAD model consists only of conic curves and surfaces, which greatly simplifies and accelerates frequently required geometry queries, like closest point computations or slicing, critical to meshing and manufacturing process planning.

The generated parametric CAD solid model and its compact CSG tree representation can be edited to refine the optimised design or to combine it with other components into a product. It is worth reminding that in industrial practice, as supported by user studies [59], geometries obtained from optimisation are often the starting point for design exploration and not the endpoint. This usually requires the ability to edit the optimised geometries in a CAD system. A high-level compact representation is also useful for taking into account geometric manufacturing constraints, like minimum thickness or slope or length of the members [30]. Beyond geometry editing, the frame model may also have advantages in structural analysis of the CAD model. The frame model makes it, for instance, easier to check the structure for buckling and inelastic deformations. As a low-order model, it can be analysed orders of magnitude faster than a three-dimensional solid model. This opens up the possibility of instant analysis of CAD models as it has become recently available in, e.g., Ansys Discovery Live or Creo Simulation Live.

In closing, one particular extension of the presented approach is worth mentioning. We assumed that the topology optimised voxel model can be approximated with tubular members. However, in three-dimensional optimisation problems, large flat or curved surfaces may also appear during optimisation. The introduced skeletonisation algorithm in Section 3 will reduce even those into a network of one-voxel-wide chains. As suggested in Lee et al. [4], the algorithm can be modified to preserve one-voxel-wide surfaces which appear during skeletonisation. This is achieved by not deleting the so-called surface points/voxels, which are defined based on the arrangement of solid voxels in an octant. In Figure 26 the topology optimised pipe bracket has been skeletonised with the modified skeletonisation algorithm. The curved surface at the bottom is preserved. In the structural model, this surface can be modelled as a thin-shell and in the CAD model it can be easily converted into a solid by providing a thickness. Currently, we do not have an automated process to extract a shell surface from the skeletonised voxel model, but this seems to be possible and is a promising direction for future research.

Appendix A Member stiffness matrix

We model the members of the structural frame as Timoshenko beams. As mentioned, the members are assumed to be straight and have uniform cross-sections along their lengths. Hence, each member can be approximated with a single beam finite element and its stiffness matrix can be derived with standard structural analysis techniques from undergraduate textbooks, see e.g. [60]. The stiffness matrix we use takes into account axial stretch, transversal shear, bending and torsion effects. The local stiffness matrix $\bm{K}_{i}^{l}$ of a beam $i$ is derived in a local coordinate system in which the beam axis is aligned with the $x_{i}^{l}$ axis. Each of the two end nodes of the beam have three displacement and three rotation degrees of freedom, see Figure 27. To ease the notation, we write in the following for the local coordinate axes only $x$, $y$ and $z$.

The element stiffness matrix of the beam is composed out of independent axial, bending and torsion stiffness matrices. The axial stiffness matrix corresponding to the displacement degrees of freedom $u_{x}^{(1)}$ and $u_{x}^{(2)}$ is given by

[TABLE]

where $E$, $A$ and $L$ denote the Young’s modulus, cross-sectional area and length.

The bending stiffness matrix corresponding to the displacement degrees of freedom $u_{y}^{(1)}$, $u_{z}^{(1)}$, $u_{y}^{(2)}$ and $u_{z}^{(2)}$ and the rotational degrees of freedom $\theta_{y}^{(1)}$, $\theta_{z}^{(1)}$, $\theta_{y}^{(2)}$ and $\theta_{z}^{(2)}$ is given by

[TABLE]

where $I_{y}$ and $I_{z}$ are the second moments of area about the $y$ and $z$ axis. The two dimensionless factors $b_{y}$ and $b_{z}$ take into account the transversal shear and are defined according to $b=12EI/(\kappa GAL^{2})$, with $G$ the shear modulus and $\kappa$ the shear correction factor. In case of a circular cross section we have $I_{y}=I_{z}=I$, $b_{y}=b_{z}=b$ and $\kappa=0.9$. Note that setting $b_{y}=b_{z}=0$ in (24) gives the stiffness matrix of an Euler-Bernoulli beam.

Finally, the torsional stiffness matrix corresponding to the rotational degrees of freedom $\theta_{x}^{(1)}$ and $\theta_{x}^{(2)}$ is given by

[TABLE]

where $J$ is the torsion constant of the cross section.

The local stiffness matrix $\bm{K}_{i}^{l}$ of an element $i$ is obtained by suitably assembling the three stiffness matrices (23), (24) and (25) into a $12\times 12$ matrix. This matrix in the local coordinate system $x_{i}^{l}$, $y_{i}^{l}$ and $z_{i}^{l}$ is transformed into a stiffness matrix $\bm{K}_{i}$ in the global coordinate system $x$, $y$ and $z$ according to (11).

Appendix B Algorithms

The proposed approach essentially relies on the Algorithms 1 and 2 mentioned in Sections 3.2 and 4.1, respectively. The input of Algorithm 1 for skeletonisation consists of the binary voxel model with $\mathcal{V}=\mathcal{V}_{s}\cup\mathcal{V}_{e}$ and the tagged non-removable voxels $\mathcal{V}_{t}$. In each iteration (loop starting line 1) a layer of border voxels in one of the six coordinate directions is removed. Within each iteration first all potentially removable voxels are determined and then removed (loops starting lines 1 and 1 respectively). Only voxels which have an empty neighbour, are not the end point of a chain, are a simple point and are not tagged as non-removable can be potentially removed. According to (22) a voxel $v$ is a simple point if its removal does not change the Euler characteristic and the number of connected components of its 26-neighbourhood $\mathcal{N}_{26}(v)$. We determine the change of the Euler characteristic with a look-up table [4, 44, 45] and the number of connected components using the Boost Graph Library [61]. The inner loop starting line 1 is necessary because each voxel is present in 27 26-neighbourhoods and the removal of one voxel can potentially render any one of the neighbouring 26 voxels non-removable.

The obtained voxel skeleton is converted with Algorithm 2 into a weighted undirected graph. Each joint voxel with $|\mathcal{N}_{26}(v)|\neq 2$ or tagged voxel $v\in\mathcal{V}_{t}$ becomes a node of the graph. The edges are obtained by starting from a joint voxel $v$ and marching along all its neighbours $\mathcal{N}_{26}(v)$, in turn, until a second joint or tagged voxel is reached (loop starting line 2). The edge weight is obtained by incrementing the variable $\ell$ during the marching (line 2). In Algorithm 2 each edge is traced twice which is redundant and could be improved.