Robust topology optimisation of microstructural details without length scale separation - using a spectral coarse basis preconditioner
RRobust topology optimisation ofmicrostructural details without length scaleseparation - using a spectral coarse basispreconditioner
Joe Alexandersen ∗ , Boyan S. Lazarov Department of Mechanical Engineering, Solid Mechanics,Technical University of Denmark, Nils Koppels All´e 404,DK-2800 Kgs. Lyngby, Denmark
Abstract
This paper applies topology optimisation to the design of struc-tures with periodic microstructural details without length scale sepa-ration, i.e. considering the complete macroscopic structure and its re-sponse, while resolving all microstructural details, as compared to theoften used homogenisation approach. The approach takes boundaryconditions into account and ensures connected and macroscopicallyoptimised microstructures regardless of the difference in micro- andmacroscopic length scales. This results in microstructures tailored forspecific applications rather than specific properties.Dealing with the complete macroscopic structure and its responseis computationally challenging as very fine discretisations are neededin order to resolve all microstructural details. Therefore, this articleshows the benefits of applying a contrast-independent spectral pre-conditioner based on the multiscale finite element method (MsFEM)to large structures with fully-resolved microstructural details. ∗ E-mail: [email protected] a r X i v : . [ c s . C E ] N ov he density-based topology optimisation approach combined witha Heaviside projection filter and a stochastic robust formulation is usedon various problems, with both periodic and layered microstructures.The presented approach is shown to allow for the topology optimisa-tion of very large problems in Matlab , specifically a problem with26 million displacement degrees of freedom in 26 hours using a singlecomputational thread.
Keywords: topology optimisation, multiscale design, microstructure,robust design, multiscale FEM, spectral preconditioner
The systematic design of novel materials with extremal properties is possibleby applying topology optimisation to the design of material microstructures.The usual approach is to consider the homogenised properties of a single unitcell. This completely decouples the problem, as one does not consider themacroscopic response at all. One may then optimise a periodic microstruc-ture to achieve a certain effective property, such as directional stiffness [1, 2],negative Poisson’s ratio or thermal expansion coefficient [3], enhanced poroe-lastic pressure-coupling [4], among many others.Topology optimisation [5] is an iterative design process which distributesmaterial in a design domain by minimising a selected objective functional,e.g. compliance, material weight or Poisson’s ratio, while satisfying a set ofconstraints. The material distribution is represented by a design field whichtakes the value one if the point is occupied with material and zero if not.In order to utilise gradient-based optimisation techniques, the design field isallowed to take intermediate values. At each iteration step, the design field isupdated using the gradients of the objective and constraint functionals. Thediscretisation determines the resolution of the optimisation process, as well asits computational cost. Therefore, in order to avoid excessive computations,the design is often limited to a single mechanical element with homogeneousmaterial properties or a single periodic cell.An alternative multiscale approach to the topological design of materialunit cells can be taken by introducing the homogenised macroscopic responsein the optimisation process. This is actually the origin of the homogenisationapproach to topology optimisation presented in the seminal paper by Bendsøeand Kikuchi [6] and later by Bendsøe [7]. A hierarchical approach is taken in28], where the macroscopic density and microstructure is designed iteratively.This approach has been applied to bone modelling [9] and is well suited forparallel computations due to the decoupled homogenisation cells [10]. Similarmethodologies have recently been proposed for nonlinear elasticity [11, 12].The lack of connectivity between the varying microstructural cells is acommon problem of the above two-scale optimisation approaches. This maybe neglectable when the microstructures are infinitely small, but when con-sidering current manufacturability, a finite size must be attributed to the mi-crostructures. Ensuring connected microstructures has been sparsely treatedin the literature, but good examples include [13] in the context of free ma-terial optimisation (FMO) and [14] in the context of functionally gradedmaterials (FGM).The convergence of periodic structures towards homogenised microstruc-tures is investigated in [15, 16]. These papers show that an optimised periodicstructure, with finite geometric periodicity, converges to optimised materialunit cells with an infinite geometric periodicity, obtained using homogeni-sation for a single microstructure repeated throughout the computationaldomain.This paper restricts itself to a single-scale optimisation approach, wherethe periodic microstructural details filling a given domain are optimised forthe macroscale response. The complete macroscopic structure and its re-sponse is considered, while resolving all microstructural details as comparedto the homogenisation approach common to all of the above papers, except[15, 16]. The approach takes boundary conditions into account and ensuresconnected and macroscopically optimised microstructures regardless of thedifference in micro- and macroscopic length scales. Furthermore, the restric-tion to microstructural details controlled by the size of the coarse mesh anda global volume constraint, implicitly introduces a maximum length scalesimilar to the results presented in [17]. Such a feature can be utilised forproviding manufacturability of the design, e.g. in the structural optimisa-tion of high-rise buildings, where the unit cell models a single building unit(room, office or apartment), or for imposing properties which are not directlyrelated or not included in the physical model due to computational cost, e.g.,requiring specified porosity in a scaffold design.Dealing with the complete macroscopic structure and its response is com-putationally challenging as very fine discretisations are needed in order to re-solve all microstructural details. This work therefore explores the applicationof a contrast-independent spectral preconditioner based on the multiscale fi-3ite element method (MsFEM), as presented in [18], to large structures withfully-resolved microstructural details. The origins of MsFEM can be tracedback to [19, 20] where special multiscale basis functions are utilised for solv-ing elliptic problems. The method has been extended further to be applicableto a wide range of linear and non-linear multiscale problems, e.g. [21, 22].The main idea is to construct basis functions which provide a good approx-imation of the solution on a coarse grid [22]. For obtaining good accuracy,the coarse basis functions need to have similar oscillatory behaviour as thefine scale solution, which can be achieved by oversampling techniques [22].The application of MsFEM to linear elasticity and material jumps isolatedin the interior of the coarse elements is reported in [23]. This conditioncannot be guaranteed in the topology optimisation process, which limits theapplicability of the original MsFEM in topology optimisation. For oscilla-tory high-contrast material properties, an alternative approach is proposedin [24, 25]. The method constructs several basis functions per coarse node,which are capable of representing the important features of the solution andthe accuracy of the approximation is controlled by the dimension of the coarsespace.The layout of the paper is as follows: Section 2 covers the theoreticalbackground of the methods used; Section 3 demonstrates the capabilities ofthe MsFEM approach through numerical experiments; Section 4 covers theimplementation details of the presented approach; Section 5 details numer-ical experiments with the approach; Section 6 goes into details with severalmicrostructural design examples; and Section 7 covers a discussion and con-clusions of the presented work.
The considered problem is static linear elasticity, which is governed by theNavier-Cauchy equation: − div σ ( u ) = f (1a) σ ( u ) = C : ε ( u ) (1b)where σ is the stress tensor, ε is the linearised strain tensor, C is the linearelastic stiffness tensor, u denotes the displacement field and f is the input to4he system in the form of distributed or concentrated forces. The displace-ments and the input forces are vector quantities with number of elementsequal to the dimensionality of the considered physical space. The consideredproblems are two-dimensional, however, the presented approach is applicableto three dimensions without any significant amendments. The system occu-pies a bounded physical domain Ω ∈ R with the boundary Γ = Γ D i ∪ Γ N i being decomposed into two disjoint subsets for each component i = 1 ,
2. Γ D i is the part of the boundary with prescribed zero displacement, u i = 0, andΓ N i denotes the part with prescribed traction, t i . The stiffness tensor C isassumed to be isotropic and has the following form C ( x ) = E ( x ) C where C is the stiffness tensor for a predefined Poisson’s ratio, ν < .
5, and unitYoung’s modulus, E = 1. The Young’s modulus E ( x ) is spatially varyingand bounded E ( x ) ∈ [ E min , E max ].The variational formulation, see e.g. [26], of equation (1) is to find u ∈ V such that: a ( u , v ) = l ( v ) for all v ∈ V (2)where V = (cid:110) v ∈ [ H (Ω)] | v i = 0 on Γ D i , i = 1 , (cid:111) ⊂ V = [ H (Ω)] andthe bilinear form a and linear form l are defined as: a ( u , v ) = (cid:90) Ω ( C : ε ( u )) : ε ( v ) d x (3a) l ( v ) = (cid:90) Ω ( f · v ) d x + (cid:90) Γ N ( t · v ) d s (3b)The weak formulation is discretised using the finite element space V h ⊂ V with vector-valued shape functions defined on a uniform mesh T h . Thediscretization leads to a linear system of equations of the form: Ku = f (4)where K is the so-called stiffness matrix, the vector u consists of the nodaldisplacements and the vector f contains the nodal forces applied to the dis-crete problem. Topology optimisation is a material distribution method that seeks to findan optimised structure for a given physical problem with respect to a certain5bjective functional subject to design constraints [5]. The physical problem,in this work linear elasticity, is discretised using the finite element method, asdescribed in section 2.1. The design is parametrised by attributing each finiteelement with a relative density, ρ e , which determines whether the element issolid, ρ e = 1, or void, ρ e = 0. In order to solve the problem using gradient-based optimisation, the binary density variables are relaxed to continuousvariables, ρ e ∈ [0 , K e ( ρ e ) = ( E min + ( E max − E min ) ρ ep ) K (5)where E min is the minimum Young’s modulus attributed to void areas inorder to avoid ill-conditioning of the stiffness matrix, E max is the maximumYoung’s modulus attributed to solid areas, ρ e is the element relative density, p is a penalisation parameter and K is the stiffness matrix for an elementwith a Young’s modulus of 1.In this work, we consider the minimum compliance problem under a con-straint on the volume fraction of solid material. The optimisation problemis posed as follows: minimise: ρ ∈ R nd f ( ρ , u ) = f T u subject to: g ( ρ ) = ρ T v v f e T v − ≤ (cid:82) ( ρ , u ) = ≤ ρ e ≤ e = 1 , ..., n d where f is the compliance functional, g is the volume constraint functional, ρ is a vector of the n d -number of design variables, v is a vector of the n d -number of element volumes, e is a n d -entry long vector of ones and (cid:82) ( ρ , u ) = K ( ρ ) u − f is the residual of the discretised system of equations, .The optimisation problem is solved using the nested formulation, where thediscretised system of equations for the state field is solved separately fromthe design problem. 6 .3 Computational issues of topology optimisation Topology optimisation is an iterative process and requires the solution of alarge system of equations for the state field, equation (4), at every designiteration. Solving the linear systems of equations can often account for morethan 99% of the total computational time [27]. Therefore, the solution cost ofthe optimisation procedure depends strongly on the time necessary for solv-ing the state equations. Realistic mechanical problems with microstructuralscales comparable to the macroscale require very fine discretisations to cap-ture the microstructural scales, making direct solvers prohibitive even in twodimensions. An alternative is to utilise iterative solvers, where the conver-gence and the solution cost is controlled mainly by the applied preconditioner.Classical preconditioning techniques, such as incomplete factorisation, sparseapproximate inverse or stationary iterative methods, cannot provide a mesh-independent number of iterations. Furthermore, for topology optimisationproblems with high contrast between the material parameters, E max /E min ,the number of iterations increases with increasing contrast [27, 28]. Geo-metric multigrid techniques [29] can provide mesh-independent solution forsmooth material properties, however, similar to the classical precondition-ers, the number of iterations depends on the contrast when the mesh is notaligned with the different materials. An effective alternative to geometric multigrid and classical preconditionerswas demonstrated in [24, 25]. The method is a multiscale approach providingthe solution of the diffusion equation with strongly heterogeneous coefficients.The solution cost is independent of the contrast of the diffusion coefficientsand is significantly reduced by using coarse space approximations which con-tain important features of the solution. The approach has been extendedand applied in topology optimisation of linear elastic problems in [18] andthe main steps are outlined below in the context of structures with periodicmicrostructural details.The fine mesh T h , utilised for discretising equation (1) in section 2.1, isassumed to be obtained by a refinement of a coarser one T H = { K j } N c j =1 ,where K denotes a coarse mesh cell and N c is the number of coarse nodes.In the context of microstructural design, the coarse cells represent the unit7igure 1: Illustration and dimensions of a double-clamped beam subjectedto a concentrated load in the centre. The beam is made up of a singleperiodic microstructure with square unit cells. The fine mesh, coarse meshand agglomerates are illustrated.cells and the fine cells the discretisation of the unit cells, which is illustratedin figure 1.The nodes of the coarse mesh are denoted as { y i } N c i =1 and the neighbour-hood of node y i is defined as: ω i = (cid:91) (cid:8) K j ∈ T H ; y i ∈ K j (cid:9) (7)which is illustrated in figure 1 using dotted lines. These neighbourhoods willbe denoted as agglomerates, since they can be viewed as a group of coarseelements agglomerated together with an overlap δ [30].A set of coarse basis functions, (cid:8) φ i,j , j = 1 . . . N i,j (cid:9) , defined with respectto T h , is introduced for each node y i in the coarse mesh with support on ω i , where N i,j is the number of basis functions for the i ’th coarse node.The coarse basis functions are based on the modes obtained by solving thefollowing local generalised eigenvalue problem on each agglomerate ω i : − div ( C ( x ) : ε ( u )) = λ E ( x ) u , x ∈ ω i (8)The local eigenvalue problem is discretised using the fine mesh, which leads tothe following discrete generalised eigenvalue problem in matrix-vector form: K ω i ψ ω i j = λ ω i j M ω i ψ ω i j (9)where K ω i is the stiffness matrix, M ω i is a stiffness-based mass matrix, ψ ω i j is the j ’th eigenvector and λ ω i j is the j ’th eigenvalue, all for the i ’th agglom-erate, ω i . It is important to note that the local matrices take the physical8oundary conditions of the full problem into account and the coarse basisthus automatically fulfils the specified boundary conditions.The eigenvalues are ordered as λ ω i ≤ λ ω i ≤ · · · ≤ λ ω i j ≤ ... and thefirst several eigenvectors corresponding to eigenvalues smaller than a pre-scribed global threshold, λ Ω , are selected to form the coarse basis. The coarsebasis functions, represented on the fine grid, are defined as φ i,j = ψ ω i j χ i ;that is they are constructed by multiplying the eigenfunctions, ψ ω i j , witha partition of unity, { χ i } N c i =1 , subordinated to ω i such that χ i ∈ H (Ω)and |∇ χ i | ≤ /H, i = 1 , . . . , N c , and H is the characteristic length of acoarse element K . Therefore, the set of coarse basis functions (cid:8) φ i,j (cid:9) asso-ciated with node y i is defined as the fine space finite element interpolant of (cid:8) ψ ω i j ( x ) χ i ( x ) , j = 1 . . . N i (cid:9) , see e.g. [31]. For each ω i , N i is determined asthe number of eigenvalues smaller than a globally selected threshold value, λ Ω .The solution in the coarse space is sought as u a = (cid:80) i,j c i,j φ i,j , where c i,j is the coefficient of the j ’th coarse basis function for the i ’th agglomerate, ω i .A coarse discretisation of the variational formulation is given as K c u c = f c ,where the coarse stiffness matrix K c and the coarse right-hand side f c areobtained as: K c = R c KR c T , f c = R c f (10)The vector u c contains the coefficients c i,j and R c T = (cid:2) φ , φ , . . . , φ N t (cid:3) ,where N t = (cid:80) N c i =1 N i , is a matrix describing the mapping from the coarse tothe fine space and consisting of nodal values of the coarse basis functions inthe fine space. An approximation of the nodal solution in the fine space isobtained as u a = R c T u c . Building the MsFEM basis for general linear elastic systems requires solv-ing the eigenproblem, equation (9), for each agglomerate and as discussed in[18] this process is computationally very expensive. It is important to notethat when dealing with periodic microstructural details, the computation ofeigenfunctions is limited to a small amount of unique agglomerates. Figure 1shows the three unique agglomerates for a double-clamped beam with a pe-riodic microstructure. Due to the imposed periodicity of the microstructure,all agglomerates along the left-hand side are the same (constrained in both x -and y -direction at the left-hand side), all agglomerates along the right-hand9ide are the same (constrained in both x - and y -direction at the right-handside) and all internal agglomerates are the same (unconstrained and allowingrigid body motion). This means that significant time can be saved on thecomputation of the local spectral bases. MsFEM can be applied as a solver in topology optimisation, however, ifthe coarse solver does not provide an accurate enough approximation of thefine-scale solution, the optimisation will results in sub-optimal designs withdisconnected features [18]. An alternative is proposed in [24], where thebasis is utilised in a two-level additive Schwarz preconditioner for the iter-ative solution of the fine scale solver. The preconditioner is scalable andresults in contrast independent number of iterations. Here, the coarse basisis simply applied as a coarse grid solve in a multigrid-like preconditioner [30]for the iterative solution of the fine-scale system of equations, equation (4).The residual of the fine-scale system is restricted to the coarse basis andthe correction is approximated using a MsFEM coarse-scale solve and pro-jected back to the fine space. In the current work, the MsFEM coarse-scalesolve is combined with a single pre- and post-smoothing using symmetricGauss-Seidel and used as a preconditioner for the generalised minimal resid-ual (GMRES) iterative method [32]. The linear elasticity problem, as wellas the preconditioner, is symmetric and positive definite and, hence, the pre-conditioned conjugate gradient (PCG) method could be used. However, theauthors’ experience is that GMRES performs better compared to PCG forthe presented problems and thus GMRES is used throughout.
In order to demonstrate the performance and applicability of MsFEM fortopology optimisation problems, a test problem is formulated. The test prob-lem is a short beam clamped at both ends and subjected to a concentratedload in the centre, as shown in figure 1. The height of the beam is set to B = 1, the length to L = 2 and the size of the force is set to P = 0 . H = 1 /
8, containing the same microstructure. Each coarse cell is discretised10 ndeformed − Periodic cell design distribution 00.10.20.30.40.50.60.70.80.91 (a) “Final” (c) Macrostructure: “Intermediate”
Undeformed − Periodic cell design distribution 00.10.20.30.40.50.60.70.80.91 (b) “Inter-mediate”
Figure 2: The investigated cross structure with a unit cell discretised by20x20 elements. The macrostructure is made up of 8 by 16 unit cells, whichis shown in two versions: (a) final 0-1 design and (b) the corresponding“intermediate” design generated by filtering the first. Subfigure (c) showsthe macrostructure for the “intermediate” design.with 20 by 20 linear finite elements, leading to a full macroscopic mesh of160 by 320. The problem is investigated for varying values of the eigenvaluethreshold, λ Ω , and MsFEM is applied both as a solver and as a preconditionerin combination with GMRES.Figure 2 shows the structure under consideration for these numerical ex-periments. The microstructure is chosen as a cross structure and is consid-ered in two versions, the first of which is the clear 0-1 design shown in figure2a. This design is representative of a final design of a topology optimisationprocess using a Heaviside projection as will be described in section 4.3. Thestructure is also considered in a smeared out version, as can be seen in figure2b, which is representative of an intermediate design during the optimisationprocess where the design usually contains significant amounts of intermedi-ate stiffnesses. The smeared out design is obtained using a density filter, aswill be described in section 4.4, with a filter radius of 4 times the elementside length. The macrostructure for the intermediate design case is shown infigure 2c. The structure is analysed with a SIMP penalisation parameter ofp=3. Similar trends are observed for the two structures and thus only theresults for the filtered version are shown.Figure 3 shows the relative error in the energy norm, (cid:107) u − u a (cid:107) K / (cid:107) u (cid:107) K =11 −5 −4 −3 −2 −1 −2 −1 E rr o r o f c oa r s e a pp r o x i m a t i o n [ % ] × × × (a) Eigenvalue threshold, λ Ω −2 −1 E rr o r o f c oa r s e a pp r o x i m a t i o n [ % ] × × × (b) Number of coarse DOF Figure 3: Relative error in the energy norm for a single MsFEM solve ofthe filtered cross structure. The error is shown as a function of (a) varyingeigenvalue threshold and (b) the corresponding number of coarse degrees offreedom. Colours denote different agglomerate sizes, see legend, and symbolsdenote different contrasts: lines; E min = 10 − , circles; E min = 10 − , squares; E min = 10 − and crosses; E min = 10 − .( u − u a ) T K ( u − u a ) / u T Ku , of the solution obtained from a single coarse-scale MsFEM solve for the intermediate design case. The analysis is donefor a constant maximum Young’s modulus, E max = 1 and a set of minimumYoung’s moduli, E min ∈ { − , − , − , − } , as well as for differentsizes of the agglomerates, denoted by the numbers of unit cells covered byeach coarse cell, { × , × , × } . It can be seen that the approximatedsolution converges with respect to both the eigenvalue threshold and thecorresponding number of degrees of freedom. Even though small differencesare observed between E min = 10 − and the lower, it can clearly be seenthat the behaviour is contrast-independent for E min ∈ { − , − , − } .It can be seen that the energy norm error follows a common trend, in relationto the eigenvalue threshold, independent of the agglomerate size, similar tothe scalar case [25]. Clear differences can be seen in the lower end of thedatasets for the error in relation to the number of coarse degrees of freedom.The curves do, however, converge towards a common trend as the number ofdegrees of freedom is increased.As can be seen from the error plots, a very large set of coarse basisfunctions need to be used to obtain a good accuracy of a single MsFEMsolve. This results in a very large computational cost, both for solving the12 −5 −4 −3 −2 −1 N u m b e r o f G M R E S i t e r a t i o n s × × × (a) Eigenvalue threshold, λ Ω N u m b e r o f G M R E S i t e r a t i o n s × × × (b) Number of coarse DOF Figure 4: Number of GMRES iterations until convergence for the filteredcross structure. The number of iterations is shown as a function of (a) varyingeigenvalue threshold and (b) the corresponding number of coarse degrees offreedom. See figure 3 for explanation of symbols and colours.generalised eigenvalue problems but also for solving the coarse-scale system.This is why the MsFEM basis is applied as a preconditioner to the full systeminstead.Figure 4 shows the number of GMRES iterations needed to convergeto a precision of ε rel = 10 − relative to the preconditioned norm of theforcing vector. As expected, the number of GMRES iterations decreasesas the eigenvalue threshold is increased and the basis is enlarged. Besidesthe results showing clear contrast-independence, they also show that theconvergence of the GMRES solver is not significantly affected by the size ofthe agglomerate. More specifically, by looking at figure 4a it can be seen thatby choosing the same eigenvalue threshold, the number of GMRES iterationsneeded to converge is very close to the same. Instead of computing the stiffness-based mass matrix for a given agglomer-ate, as specified by equation (8), it is possible to simply use the diagonal ofthe agglomerate stiffness matrix as the weighting matrix of the generalisedeigenproblem, equation (9). This is because the two weighting matrices arespectrally equivalent and λ diag( K ) ≈ Ch λ M for some constant C [29, 33].This is interesting from a computational point of view, as the diagonal leads13o a cheaper solution procedure; both in that the stiffness-based mass ma-trix is not assembled, and in that the matrix-products involved in solvingthe eigenproblems become cheaper. The numerically observed convergencebehaviour is very similar to that seen in figure 3 and is thus not shown. Thediagonal of the agglomerate stiffness matrix is used as the weighting matrixthroughout the rest of this paper. Considering the total computational time taken for the projection operation,equation (10), and the iterative solution, there exists an optimal basis sizethat provides the fastest solution process. Many factors influence this time,such as the sparsity of the projection matrix and also the number of iterationstaken to reach convergence. It should also be noted that the times are not thesame for the three agglomerate sizes, because the larger the agglomerate, thelarger the support area and the denser the projection matrix becomes. Thus,it is clear that increasing the size of the basis does not necessarily lead to afaster solution even though it results in fewer GMRES iterations, as each ofthese become more expensive along with the projection operation. It can thusbe concluded that it cannot pay to increase the size of the agglomerate, asit does not appear to influence the error nor the convergence of the iterativesolver, figure 4, as well as leading to a more expensive solution proceduredue to larger eigenproblems and denser projection matrices.
The main assumption in standard homogenisation is that the macroscalemedium consists of an infinite number of periodic unit cells. The homogenisedsolution is an asymptotic solution, which is valid only when the cell character-istic length H is orders of magnitude smaller than the macroscopic problemscale, i.e., H << L . For cell characteristic lengths in the intermediate range,
H < L , the homogenised solution can differ significantly from the true so-lution. Also, the homogenised material properties cannot account for theboundary conditions and any load variations of the fine scale problem.Furthermore, in the design of functionally graded materials, i.e. ma-terials with spatially varying properties, the periodicity condition is oftenviolated. The transition from one cell type to another cell type is appliedabruptly within a distance comparable to H . In the above mentioned cases,14 −3 −2 −1 −3 Size of coarse cells, H C o m p li a n ce HomogenisedFull scale (a) Double-clamped beam with concentrated load −3 −2 −1 Size of coarse cells, H C o m p li a n ce HomogenisedFull scale (b) Cantilever beam with distributed load
Figure 5: Convergence of compliance with respect to the size of the coarse cellusing both homogenisation and full scale analysis using MsFEM-GMRES.the proposed MsFEM-GMRES approach will provide the physical responsewithout any simplifications for relatively low computational cost. Hence, itwill be superior, but more costly, compared to the standard homogenisationfor optimisation problems, where localised constraints (loading, boundaryconditions, surface effects) significantly affect the design topology.Figure 5 shows the change of compliance with respect to the size of thecoarse cell using both homogenisation and full-scale analysis using MsFEM-GMRES. The compliance is shown both for the double-clamped beam with aconcentrated load in the centre, shown in figure 1, as well as a cantilever witha distributed load at the end, as will be introduced in section 6.2. It should benoted that convergence is not expected for the case with a point-load, whichis also observed for both methods in figure 5a. It is important to note thatany observed convergence is slightly different for the two solution types. Forthe full-scale analysis, the size of the microstructure decreases when decreas-ing the size of the coarse cells. The convergence of compliance with respectto microstructure size is assumed to be dominant over the increased finiteelement resolution of the macroscopic domain. For the homogenised case,decreasing the size of the coarse cells solely increases the resolution of thefinite element discretisation of the homogenised problem. That is, the mi-crostructure size can be seen as constant and is assumed to be small enough15o fulfil the assumptions of homogenisation. Thus, figure 5 illustrates thatthe macroscopic compliance of structures with finite size periodic microstruc-tural details converges towards that estimated by homogenisation. However,it also emphasises that if one attributes a large finite size to a microstruc-ture designed using homogenisation, one may be introducing a substantialerror in the macroscopic response. This is especially pronounced for prob-lems with localised behaviour as the illustrated double-clamped beam witha concentrated load in the centre, figure 5a.
Referring to figure 1, the design domain is made up of repeated unit cellscovering the entire domain. The unit cells contain a periodic microstructurethat is discretised with a certain number of finite elements and each of theseare assigned a relative density which is controlled by a corresponding designvariable. Due to the imposed periodicity of the design, the same designvariable controls the relative density of several fine-scale elements in themacroscopic problem. The total sensitivity for a given design variable is thuscalculated by summing up local sensitivities of all fine-scale elements whichit controls.It is important to emphasise that the state field is solved on the macro-scopic level resolving all microstructural details. The objective functional isthus the fine-scale macroscopic compliance. By considering the full structure,one ensures that all microstructural details will be connected in order to pro-vide a beneficial macroscopic behaviour, as disconnected members would bedetrimental to the compliance objective.One of the aims of this paper is to investigate the performance of the spec-tral basis preconditioner, as described in section 2.4, and the performanceis thus compared to reference results. For calculating reference results, theequation system is solved using the standard backslash ( mldivide ) in
Mat-lab and these are then compared to the results obtained when using thespectral preconditioner based on MsFEM combined with the standard
Mat-lab implementation of the GMRES iterative method ( gmres ).The optimisation problem is solved using the method of moving asymp-totes (MMA) [34], of which the used
Matlab implementation is courtesy of16rister Svanberg.
In order to use a gradient-based optimisation method, such as MMA, oneneeds the sensitivities, or gradients, of the objective functional and any con-straint functionals. These sensitivities are easily found for the structuralcompliance using the discrete adjoint method [35]. Because the problem isself-adjoint, the sensitivity for the e ’th element simply becomes: ∂f∂ρ e = − u e ∂ K e ∂ρ e u e (11)where u e is the element nodal displacements and K e is the element stiffnessmatrix. The derivative of the element stiffness matrix with respect to theelement density is easily found as: ∂ K e ∂ρ e = p ( E max − E min ) ρ p − e K (12) Heaviside projection filtering [36, 37] is applied in order to facilitate the for-mation of crisply-defined topologies, or 0-1 topologies, without intermediatedensities. First, the design variables are filtered using density filtering [38]and subsequently the filtered density variables are passed through a smoothedHeaviside function in order to project the densities below or above a giventhreshold to 0 or 1, respectively:¯˜ ρ e = tanh( βη ) + tanh( β ( ˜ ρ e − η ))tanh( βη ) + tanh( β (1 − η )) (13)where η ∈ [0; 1] determines the projection threshold and β > ρ , have been filtered, ˜ ρ , and projected, ¯˜ ρ ,one has obtained the actual physical densities. These physical densities areused for the calculation of the element stiffness matrix, equation (5), andtherefore, the sensitivities, equation (11), are corrected using the chain rule,as detailed in e.g. [39]. 17 a) Periodic filtering (b) Multipleblocks Figure 6: Illustration showing the filtering procedures: (a) single periodicmicrostructure where the filtering is performed on a 3 × × All manufacturing methods introduce sources of uncertainty in the final re-alised structure, the most prevalent of which is an uncertainty in the manu-factured geometry. This is caused by either too little or too much materialbeing removed or added during the manufacturing process. The possiblevariability in the manufactured structures can be taken into account by con-sidering several design realisations during the optimisation process [37].The optimisation problem is cast in a stochastic robust formulation [41],in which the mean and variance of the macroscopic compliance is considered.The different realisations of the design are obtained by using several pro-jections using different threshold values simulating uniform over- and under-18aterial deposition or removal. The optimisation problem is posed as follows:minimise: ρ ∈ R nd f ( ρ , u ) = E (cid:2) f T u (cid:3) + κ (cid:112) Var [ f T u ]subject to: g ( ρ ) = E (cid:2) ¯˜ ρ T v (cid:3) v f e T v − ≤ (cid:82) ( ρ , u ) = ≤ ρ i ≤ i = 1 , ..., n d where E [ · ] and Var [ · ] denotes the expected value and variance of a givenquantity, respectively, and κ is a weighting factor for the inclusion of thevariance in the objective. The sensitivities are computed as described in[41].Throughout this paper, the weighting factor, κ , is set to 1.0 and threedesign realisations are used with equal weights, specifically the set η ∈{ . , . , . } . The final designs are verified using a large set of uniformlydistributed points, η ∈ [0 . , . β approach for moderate β A constant β approach as presented in [42] is adopted in order to eliminatechanges in the allowable minimum length scale during the optimisation pro-cess. By having a constant β , controlling the slope of the projection, in com-bination with a constant filter radius and constant thresholds, one ensures aconstant minimum length scale. This ensures that small features that cannotbe supported by the mesh and length scale do not appear during the opti-misation process. Using the classical β -continuation, these smaller featuresare progressively removed during the optimisation process and convergenceissues can be encountered along the way. The approach relies on simplechanges to the size of the initial asymptotes in the MMA algorithm, whichis set to . β . No external move limits are imposed and thus the full controlof the updates are left in the hands of MMA. Even though one removes theneed to use a continuation approach on the β -parameter by making thesemodifications, it is observed to be beneficial to perform continuation of theSIMP penalisation parameter in order to get high quality solutions, as alsonoted in the original paper [42]. If starting from a homogeneous intermediateguess, the constant beta approach retains the possibility to relatively freelyform topologies. However, if starting from a random guess, many areas will19e projected to solid and void and the constant beta approach then in manyways resembles some form of level set approach.For the single design results presented in this work, a continuation ap-proach is used on the SIMP penalisation parameter, increasing it in steps of p ∈ { , , , , } when the maximum relative change in objective functionreaches below 0 .
001 or when the number of iterations for the same parame-ter size reaches 50. After the last step, the optimisation is stopped when amaximum number of iterations is reached or the maximum relative changecriteria is met.Unfortunately, when combined with the robust topology optimisation for-mulation, one is somewhat restricted in the maximum size of β . The size of β needs to be moderate in order to ensure the ability of the optimiser toform a good initial topology from both homogeneous and random startingdistributions. The important thing is for the smooth parts of the shiftedapproximative Heaviside functions to overlap, or rather that the parts withnon-negligible sensitivities overlap, and this sets an upper limit for the ini-tial β for the optimisation to proceed in practice. Thus, the robust topologyoptimised designs in this paper have been obtained from a three step contin-uation approach. During the initial stage where the topology is forming, theSIMP penalisation parameter is set to p = 1 . β is set to a moderatevalue of β ∈ [8 ,
32] depending on the filter radius - the larger the filter radiusrelative to the element size, the larger β is needed to provide the same levelof intermediate densities [36]. During the next stage, the SIMP penalisationparameter is raised to p ∈ [3 . , .
0] depending on the problem - some mi-crostructural design problems require higher penalisation than others to form0-1 topologies. For the final stage of the optimisation, the β parameter isincreased to β ∈ [32 , β is already set to a moderately high value from the beginning reduces thelikelihood of unsupportable features forming during the optimisation process.The next step is introduced after 100 iterations or when the maximum rela-tive change in objective function reaches below 0 . .5 Multiple blocks Restricting the design to be made up of a single periodic microstructure isquite a severe constraint on the design freedom. The obtained designs per-formance will be far from the unconstrained optimal design. This is due tothe fact that the same microstructure is placed throughout all parts of thedesign domain and the optimised microstructure will thus be one that sat-isfies optimality in an average sense, but where the parts with the largestsensitivities will dominate the optimisation procedure. This restriction canbe partially alleviated by partitioning the design domain into multiple sub-domains, within which a unique locally periodic microstructure is designed.Smoothly connected microstructures, both within the subdomains andacross their interfaces, are obtained by considering the full macroscopic prob-lem and by imposing periodicity and connectivity constraints on the opti-mised design. This is achieved by using the PDE-based density filter [40],which easily allows for imposing such constraints on the filtered density fieldby collapsing the corresponding degrees of freedom in the discrete filter prob-lem. Figure 6b illustrates the periodicity and connectivity constraints im-posed on the filtered density field for a design with three layers of locallyperiodic microstructures. The colours indicate which edges are coupled to-gether explicitly by collapsing their degrees of freedom. The top and bottom(red) of all cells are to be connected ensuring both periodicity, in the verti-cal direction, within the block and connectivity across the interfaces. Theleft and right sides (blue, green and magenta) of each cell type are to beconnected ensuring periodicity within the blocks in the horisontal direction.The increase in the number of microstructures present in the design,expands the number of unique agglomerates to consider when calculatingthe eigenmodes for the spectral basis. For instance, considering a double-clamped beam where the design is made up of three layers with each theirmicrostructure, the number of agglomerates increases to 15. There is 3 setsto accommodate the different boundary conditions, as discussed for the singlemicrostructure case in section 2.4.1, within which there now exists 5 uniqueagglomerate types, more specifically one for each layer and one for each in-terface. 21 (a) Iteration 1 (b) Iteration10 (c) Iteration30 (d) Iteration170 (e) Iteration190 (f) Iteration200
Figure 7: The design distributions for the first and final stages of the opti-misation process; design iterations 1, 10, 30, 170, 190 and 200.
It is rather trivial to argue that while the design is evolving slowly, whichis prevalent during the first and final stages of the presented optimisationprocedure, it is possible to reuse the previously computed spectral basis forthe iterative solution. In order to test this hypothesis, the previously de-scribed double-clamped beam problem, as shown in figure 1, is optimisedusing different strategies for updating the MsFEM basis during the optimi-sation process. The beam consists of 4 by 8 coarse cells within which themicrostructure is discretised using 40 by 40 elements. The investigation iscarried out for a constant β = 64 using a single realisation with threshold η = 0 .
5. The maximum number of design iterations is set to 200 in order tokeep the total number of iterations the same for all the different strategiesin order to conduct a fair comparison.Figure 7 illustrates the relatively slow development of the design duringthe first and final stages of the optimisation. It can be seen that the designis evolving quite slowly during the first iterations, figures 7a to 7c, due tothe low SIMP penalisation value and this supports the arguments that theMsFEM basis does not need to be updated very frequently during this period.Furthermore, it is clearly shown that the design is hardly changing duringthe final iterations, figures 7d to 7f, where the topology has been determinedand the final minor adjustments are being carried out.22 .2 Simple update scheme based on GMRES-iterations
In [18] the spectral basis is updated for the agglomerates within which thedesign is changing more than a specified tolerance. Here a simple heuristicupdate scheme based on the number of GMRES iterations is used. Thespectral basis is updated when the number of GMRES iterations changesmore than a given percentage from the design iteration where the basis waspreviously computed.In order to further speed up the iterative solution of the linear systemsof equations, the GMRES solver can be initialised with the displacementsolution from the previous design iteration as the initial guess for the GMRESsolver for the new design iteration, rather than starting from a zero-vectorinitial guess. When initialisation of GMRES is used, the heuristic basisupdate rule is modified so that if the number of GMRES iterations decreasesafter a basis update, the update threshold is updated to be with respect tothe lowest number of GMRES iterations encountered after the previous basiscomputation.Furthermore, it has been shown in [43, 28] that the convergence criteriacan be relaxed significantly during topology optimisation when using an iter-ative solver. The reason for this is that one does not necessarily require thedisplacement solution to be solved to full precision in order for the optimisa-tion to progress correctly. This leads to significant reductions in computationtime compared to solving to full precision. For the investigations presentedin this section, the stopping tolerance is set to ε rel = 10 − relative to thepreconditioned norm of the forcing vector.A metric that is correlated to the development of the design is the non-discreteness measure, M nd , as introduced in [39]: M nd = (cid:80) n d i =1 ρ i (1 − ¯˜ ρ i ) n d (15) M nd is an indication of the amount of intermediate densities and thus anindication of the contrast of the stiffness distribution. When M nd is high,there is large amount of intermediate densities and when M nd is low, there issmall amount of intermediate densities. For volume-constrained complianceminimisation where the design will never contain all 0 or all 1 densities, thisis equivalent to saying: when M nd is high, there is low contrast and when M nd is low, there is high contrast. The top plot of figure 8 shows how M nd changes during the optimisation process. It is clearly seen that the level of23
20 40 60 80 100 120 140 160 180 200050100
Design iteration M n d CU Design iteration N u m b e r o f G M R E S i t e r a t i o n s ConstantHeuristicHeuristic - initialisedHeuristic - initialised, low λ Ω Figure 8: Development of the non-discreteness measure, M nd , and a compar-ison of the number of GMRES iterations for the different basis computationstrategies for a relative tolerance of ε rel = 10 − . The lines show the num-ber of GMRES iterations for the various strategies and the symbols showwhen the basis is updated. Crosses highlight forced updates due to changein continuation parameter.intermediate densities is relatively constant at the first and final stages ofthe optimisation process and that the contrast goes from low to high. Thebottom plot of figure 8 shows the number of GMRES iterations during theoptimisation process for the described strategies. It can be seen that thenumber of GMRES iterations is more or less constant throughout the op-timisation process when the MsFEM basis is recomputed at every designiteration (denoted by “constant” in the legend). Keeping the developmentof M nd in mind, this clearly supports the evidence that the MsFEM pre-conditioner is contrast-independent for a constant λ Ω . This is of course anexpensive approach and results in a total time for the entire design processof 768.5 seconds, which is significantly higher as compared to the referencetime of 306.5 seconds when using the direct solver.When adopting the heuristic update rule, the number of GMRES itera-tions can be seen to remain very close to constant throughout the optimi-sation process (denoted by “heuristic” in the legend). Here it can be seenthat the MsFEM basis is only computed a total of 7 times during the 200design iteration optimisation process, yielding a significantly lowered total24ime for the entire design process of 444.0 seconds. It is interesting that ofthe 7 updates, the first is the initial basis computation and 4 are when thepenalisation parameter, p , is updated which is a requirement made in thealgorithm.However, this is still significantly slower than the reference design processusing the direct solver. Thus, as suggested, the displacement solution fromthe previous design iteration is used as the initial guess for the GMRESsolver for the new design iteration. This significantly decreases the numberof iterations needed to obtain the required accuracy when the design is slowlyevolving. This is shown in figure 8 (denoted by “heuristic - initialised” inthe legend), where it can be seen that 4-10 GMRES iterations are needed forlarge parts of the optimisation process, namely the initial and final stages,whereas significantly more are needed in the middle stage where the topologyis forming and the design is changing rapidly.By lowering the eigenvalue threshold, from λ Ω = 6 . × − to λ Ω =4 . × − , and thus decreasing the basis size, it is possible to decreasethe computational work associated with the projection and a single GMRESiteration, as discussed in section 3.3. As can be seen in figure 8 (denoted by“heuristic - initialised, low λ Ω ” in legend), the required number of GMRESiterations generally increases as expected, however, despite this the resultingoptimisation process only takes 298.7 seconds, which is 7.8 seconds fasterthan the reference time of 306.5 seconds when using the direct solver. Thisis an interesting observation and shows that the choice of the eigenvaluethreshold, λ Ω , is hardly trivial. Further investigation into determining theeigenvalue threshold using error estimates is the subject of future research.All four final topologies are qualitatively very similar to the reference caseand thus only a selected one is shown. The final topology for the optimisationprocess using the heuristic update rule combined with an initialised GMRESsolver (“heuristic - initialised”) is shown in figure 9c and can be seen to bevery similar to the reference design obtained using the direct solver shownin figure 9a, except for hard corners as compared to rounded corners. Thefinal compliance values are as follows; reference using direct solver: C end =1 . × − , “constant”: C end = 1 . × − , “heuristic”: C end = 1 . × − , “heuristic - initialised”: C end = 1 . × − , “heuristic - initialised,low λ Ω ”: C end = 1 . × − . It is very interesting to note that all designsobtained using the MsFEM-GMRES approach actually performs better thanthe reference case, specifically 1 . . .
9% and 1 . (a) Directsolver -reference case (c) Deformed macrostructure: Direct solver - referencecase (b) Iterativesolver - relaxedconvergencecriteria andseeded initialguess Figure 9: Comparison of (a) the reference design by use of a direct solverand (b) the optimised design when using the proposed iterative method withrelaxed convergence criteria. The deformed macrostructure is shown in (c)and is made up of 4 by 8 unit cells. 26ptimisation problem and can in no way be guaranteed in general.
Requiring robustness of the design response with respect to uniform erosionand dilation usually results in designs with similar topologies across the de-sign realisations. This property can be utilised in the construction of theMsFEM preconditioner in order to further reduce the computational cost.Two strategies can easily be identified: 1) an offline building of the basis forthe possible variations of the design and using it in an online optimisationprocess [44], and 2) building a basis for a representative design realisationand re-utilising it in the solutions of all other realisations. The first strategyis computationally expensive in the preparatory part (offline computations)compared to the second one, however it is expected to provide faster onlinecomputations. The second strategy does not require any additional modi-fications of the algorithm for building the coarse basis and provides a fastsolution for the online optimisation process, as demonstrated later. There-fore, the second strategy is utilised in the examples presented in section 6.The mass matrix in the eigenproblem defined, equation (9), acts as a filterwhich selects modes dominant within the solid region. As the dilated designrealisation embeds all other design realisations, the basis for the dilated de-sign will provide a good basis for all other design realisations. Therefore,the representative design is selected to be the dilated one. For small differ-ence between the design realisations, which is the case here, the basis forthe dilated design ensures a similar number of iterations for all solutions. Ifthe difference between the design realisation is large, more rigorous selectioncriteria needs to be derived, which is left for future research.
Throughout the microstructural design examples, the contrast in stiffness isconstant and the Young’s moduli of the solid and void are E max = 1 and E min = 10 − , respectively. 27 .1 Double clamped beam with concentrated load The first numerical example is the double clamped beam with a concen-trated load, shown in figure 1, as considered earlier in section 3.1. Theproblem is first investigated using the standard single design topology opti-misation formulation and secondly using the robust formulation. The prob-lem is optimised for an increasing number of unit cells in order to inves-tigate the convergence of the compliance, as well as to compare the speedof the MsFEM-GMRES iterative solver, combined with the proposed basisre-utilisation scheme, to using the standard
Matlab direct solver. In orderto make a fair comparison between the two, all of the timing runs have beenexecuted on the same computer, where
Matlab has been restricted to usea single computational thread. The computer is a Dell Precision T7500 withtwo Intel Xeon X5650 CPUs and 96GB memory, running CentOS 6.4 and
Matlab
The problem is first investigated using the standard single design topologyoptimisation formulation. The following parameters are used: r min = 4 h , v f = 0 . β = 64, η = 0 . λ Ω = 6 . × − , ε rel = 10 − .Figure 10 shows the optimised periodic microstructures for the double-clamped beam with a concentrated load in the centre. The microstructureis shown for increasingly smaller unit cell sizes. The number of unit cellsacross the height of the beam is characterised by M x = B/H , which will beused throughout the paper. It is observed that the microstructure convergesvery fast and the overall topology does not qualitatively change for M x > M x ∈{ , , , , , , , } . This trend was also observed in [16] for the sameproblem, although the optimised designs differ. This difference in topology islikely due to a difference in the effective length scale of the current approachand that of [16]; the former using density filtering combined with a projectionat η = 0 . a) M x = 2 (b) M x = 4 (e) M x = 64(c) M x = 8 (d) M x = 16 Figure 10: Optimised microstructures for the double-clamped beam with aconcentrated load. The microstructure is shown for varying number of cellsacross the height, M x ∈ { , , , } , and the full macrostructure is shownfor M x = 64.using a direct solver, however, slight deviations are observed in the compli-ance values.From the numerical experiments of this section, it is also observed that interms of average time per design iteration, the direct solver approach is fasterfor smaller problems, but that the MsFEM-GMRES approach becomes fasterfor larger problems, M x >
6, as expected. For M x = 16 the observed averagetimes are 59.44 and 83.41 seconds for the MsFEM-GMRES and direct ap-proaches, respectively. The relative difference is 28 .
7% and thus, the savingsin time are quite substantial. Using a linear fit, the projected reduction intime for M x = 32 and M x = 64 is 41 .
0% and for 54 . M x = 64 it has been observed in practice that it takes approximately26 hours for 200 design iterations. This yields a significantly lower averagetime per design iteration than the projected, yielding an expected reductionin time of 79 . Matlab re-stricted to a single computational thread and taking into account that for M x = 64, the total number of elements is 13 , ,
200 and the total number ofdegrees of freedom is 26 , , The problem is secondly investigated using the robust topology optimisationformulation. The following parameters are used: r min = 4 h , v f = 0 . β =16, β = 64, p = 5, η = 0 . λ Ω = 6 . × − , ε rel = 10 − .Figure 11 shows the optimised robust periodic microstructures for thedouble-clamped beam with a concentrated load in the centre. As for thesingle design case, it can be seen that the microstructure converges veryfast and the overall topology does not qualitatively change for M x > M x = { , , , , , } . It is interesting to note, that the robust topologies are verysimilar to the optimised topologies presented in [16]. The main difference isthat the topologies presented in figure 11 have been optimised to be robustwith respect to erosion and dilation of the design, which can be seen in figure12 for M x = 16. Comparing the non-robust and robust topologies, figures 10and 11 respectively, it can be observed that the small horisontal and verticalbars in the left- and right-hand side of the unit cell of the non-robust designare not present in the robust designs. This is due to the length scale imposedby the robust formulation, as opposed to the lack of length scale for thenon-robust result.As for the non-robust case, it is observed that the compliance increases asthe number of unit cells is increased. It is also seen that in terms of averagetime (total for all three realisations) per design iteration, the direct solverapproach is faster for smaller problems, but the MsFEM-GMRES approachbecomes faster for larger problems, M x >
8. For M x = 16 the observed (a) M x = 2 (b) M x = 4 (c) M x = 16 Figure 11: Optimised robust microstructures for the double-clamped beamwith a concentrated load. The microstructures are shown for varying numberof cells across the height, M x = { , , } .30 a) Eroded (b) Intermedi-ate (c) Dilated Figure 12: The three realisations of the optimised robust microstructures forthe double-clamped beam with a concentrated load and M x = 16. Design iteration N u m b e r o f G M R E S i t e r a t i o n s DilatedIntermediateErodedBasis update
Figure 13: The number of GMRES iterations as a function of design iterationfor the double clamped beam with a concentrated load using the robustformulation and M x = 16.average times were 188.91 and 246.68 seconds for the MsFEM-GMRES anddirect approaches, respectively. The relative difference is 23 .
4% and thus,the savings in time are also quite substantial for the robust case. Using alinear fit, the projected reduction in time for M x = 32 and M x = 64 is 36 . . β , the three realisations are significantly31igure 14: Illustration and dimensions of a cantilever beam subjected to adistributed force along the right-hand side. The beam is made up of threelayers containing separate periodic microstructures with square unit cells.different in the initial phase of the optimisation process when the topology isbeing formed. It is in this phase, where large differences in iteration numbersare observed, especially for the eroded case. Figure 13 shows the number ofGMRES iterations during the optimisation process for M x = 16. It can beseen that around design iterations 25-35, the number of GMRES iterationsfor the eroded design peaks at over 100. It is observed that the peak inGMRES iterations occurs at the same time as the designs are changing veryquickly. When the topology has formed and the three realisations are similar,the performance can be seen to become much better and the difference inthe number of GMRES iterations is relatively small. One could argue thatit would be better to adopt the usual β -continuation, since these problemsare not observed for this approach, because the realisations are very similarthroughout the optimisation process. However, it has been observed thatthe presented approach needs significantly fewer design iterations to producea competitive design and the approach is thus still faster, even with theincreased computational time associated with the mismatch of the realisa-tions. Improvement of the spectral basis for robust topology optimisation isa subject of future research. In order to demonstrate the methodology for design domains with multiplemicrostructural layers, a cantilever, as shown in figure 14, is investigated.The cantilever beam is subjected to a distributed force in the positive x -direction along the right-hand side. The height of the beam is set to B = 1,the length to L = 2 and the size of the force is set to P = 0 .
01 per unit32 a) M x = 2 (b) M x = 4 (c) M x = 8 (d) M x = 12 (e) M x = 16(f) M x = 2 (g) M x = 16 Figure 15: Optimised robust microstructures for the cantilever beam with adistributed load. The microstructures are shown for varying number of cellsacross the height, M x = { , , , , } . The full macrostructures are shownfor M x = 2 and M x = 16.length. The beam is made up of three layers with equal thickness, eachsubdivided into square coarse cells, with an edge length of H , containing alocally periodic microstructure. Each coarse cell is discretised with 40 by40 linear finite elements. The following parameters are used: r min = 4 h , v f = 0 . β = 16, β = 64, p = 5, λ Ω = 6 . × − , ε rel = 10 − .The problem is first investigated for a single periodic microstructure andafterwards for a varying number of coarse cells across the thickness of thethree layers. Finally, the problem is investigated for varying thickness of theouter two layers. This problem was also investigated in [16], however, in the present work,the distributed load is in direct contact with the design domain, whereas asolid non-design region was included in [16] in order to allow for a directcomparison with results obtained through homogenisation. This is not donein the current work, in order to highlight the fact that by considering thefull macroscopic structure, one actually obtains a design that automaticallytakes the boundary and loading conditions into account.33 a) Eroded (b) Intermedi-ate (c) Dilated
Figure 16: The three realisations of the optimised robust microstructures forthe cantilever beam with a distributed load and M x = 4.Figure 15 shows the optimised robust microstructures for the cantileverbeam with a distributed load along the right-hand edge. It can be observedthat, unlike for the double clamped beam problem of section 6.1, the mi-crostructure does not converge within the observed range of unit cell refine-ment. It can be seen that the microstructure has the same overall topologyfrom M x = 4 −
16, but that the size of the holes are changing. By performinga cross-check analysis, this development in the topology appears to be dueto local minima. This is not considered as important for the current study,as the aim of this paper is not to investigate the convergence of the unit celldesign with respect to unit cell size .Comparing the microstructures in figure 15 to those presented in [16], animportant difference is the vertical bar at the right-hand side of the unit cell.This vertical bar is needed for the unit cells in contact with the distributedload and due to the imposed periodicity it exist in all unit cells within thebeam, as seen in figures 15f and 15g. Of course, this leads to an increasinglysub-optimal design within the design domain, as the vertical bar is onlyneeded at the very edge. But as already stated, it is seen as important toillustrate that the presented methodology, of taking the entire macrostructureinto account, ensures a physically valid design with respect to boundary andloading conditions. The optimised microstructures can be seen to exhibita combination of bending and shear stiffness, the first characterised by thehorisontal bars and the latter characterised by the cross-like members. Thisis as expected, as forcing the microstructure to be the same all over the shortbeam, ensures that the microstructure needs to exhibit both characteristics.Figure 16 shows the three realisations of the final optimised design for M x = 4. Here it can be clearly seen, that the design is robust with respectto erosion and dilation. It is important to point out that the vertical bar34igure 17: Optimised robust layered macrostructure for the cantilever beamwith a distributed load for M x = 18 and layer thickness (cid:8) , , (cid:9) .ensuring connection to the distributed load is ensured for the eroded designalso, as seen in figure 16a. One of the main benefits of analysing the full macrostructure with all mi-crostructural details resolved, is the ability to design structures with varyingmicrostructure and ensuring that these microstructures are connected. Thepresented implementation has the ability to optimise layered microstructures.The cantilever beam with a distributed load along the right-hand edge is splitinto three layers of equal thickness each made up of a number of periodic cells,as illustrated in figure 14.Figure 17 shows the full macrostructure with optimised robust layeredmicrostructures for M x = 18 with layer thickness (cid:8) , , (cid:9) , where the firstentry denotes the relative thickness of the top layer, the second denotes themiddle layer and the third denotes the bottom layer. It can clearly be seenthat different microstructures exist in the three layers. As expected, thetop and bottom layers exhibit microstructures with a high bending stiffness,characterised by the thick horisontal bars, whereas the middle layer exhibitsa microstructure with a high shear stiffness, characterised by the cross struc-ture. This is perfectly in line with the stress distribution of a short beamunder bending; high axial stresses near the top and bottom and high shearstresses in the middle. Furthermore, due to the skew-symmetry of the stressdistribution, the optimal topology should be symmetric [45] which is alsoobserved for the optimised designs. It is important to note that no symme-try has been imposed during the optimisation process. It is observed thatthe average compliance does not change significantly with respect to increas-35 , , (cid:9) (cid:8) , , (cid:9) (cid:8) , , (cid:9) (cid:8) , , (cid:9) (cid:8) , , (cid:9) t o p m i dd l e b o tt o m Figure 18: Optimised robust microstructures for the cantilever beam with adistributed load for three layers of varying outer layer thickness.ing the number of coarse cells. Also, the microstructural topology does notchange significantly and these are therefore not shown.
The three-layered cantilever beam with a distributed load along the right-hand edge is now investigated for varying thickness of the outer layers. Figure18 shows the optimised robust microstructures for the three-layered cantileverbeam with a distributed load for varying outer layer thickness. It can be seenthat the microstructure of the middle layer does not change significantly whenthe thickness of the outer layers is reduced, in turn leading to an increasein the thickness of the middle layer. It should be noted, that although thetopologies for (cid:8) , , (cid:9) and (cid:8) , , (cid:9) appear significantly different thanfor the others, they are merely the same topology as the others but shiftedvertically half a coarse cell. However, the topology of the outer layers doeschange quite significantly as the thickness of the outer layers is reduced. Itappears that the outer cells are becoming more dense and this is supportedby looking at the local volume fraction of the outer coarse cells. Figure 19bclearly shows that as the thickness of the outer layers is reduced, the localvolume fraction of material is increased in the outer layers. On the contrary,the local volume fraction of material in the middle layer is more or lessconstant. Figure 19a shows the average compliance and it can be seen that as36 .1 0.15 0.2 0.25 0.3 0.358910x 10 −3 Relative thickness of outer layers A v e r ag ec o m p li a n ce (a) Average compliance Relative thickness of outer layers V o l u m e f r a c t i o n OuterMiddleGlobal (b) Local volume fraction
Figure 19: Average compliance and local volume fractions as a function ofthe relative thickness of the outer layers for the three-layered cantilever beamwith a distributed load. (a) (cid:8) , , (cid:9) (b) (cid:8) , , (cid:9) Figure 20: Optimised robust layered macrostructure for the cantilever beamwith a distributed load for two selected layer thicknesses.the thickness of the outer layers is reduced, the compliance decreases. Despitethe compliance anomaly for (cid:8) , , (cid:9) , likely due to a local minimum ofthe optimisation problem, it is concluded that it is beneficial to have thinouter layers with very high volume fractions and a thick middle layer with arelatively low volume fraction. This is perfectly in line with sandwich theory,where very thin and axially-stiff plates are joined to a low-density shear-stiffcore in order to best resist the high axial stresses near the top and bottomand high shear stresses in the middle.In order to fully appreciate the adaptivity of the approach to the vary-ing outer layer thickness and the relation to sandwich theory, figure 20shows the full macrostructure for two select layer thicknesses: (cid:8) , , (cid:9) and (cid:8) , , (cid:9) . Here it can be clearly seen that the middle microstructuretopology is essentially the same and that the outer layers become increasingly37igure 21: Illustration and dimensions of a cantilever beam subjected to aconcentrated force at the bottom rightmost corner. Every layer has a uniquemicrostructure which is periodic in the horisontal direction only.dense. In closing, in order to showcase the capabilities of the proposed methodology,a cantilever beam with a concentrated load at the lower right-hand corneris optimised where every layer of microstructure is different. Figure 21shows the problem layout, where it can be seen that every layer has a uniquemicrostructure which is periodic in the horisontal direction only. A smallchange is made to the filtering procedure, described in section 4.5, which isthat internal periodicity in the vertical direction is no longer required. Thatis, the constraints applied on the horisontal edges of the microstructures, seefigure 6b, are no longer imposed and the filtering domain essentially becomesa single vertical slice of the design domain with periodic boundary conditionsin the horisontal direction. The following parameters are used: r min = 6 h , v f = 0 . β = 8, β = 32, η ∈ { . , . , . } , λ Ω = 5 × − , ε rel = 10 − .Figure 22 shows the optimised robust macrostructure for the cantileverbeam with a concentrated load and twenty-four layers, M x = 24. It can beseen that the microstructural details vary continuously throughout the thick-ness. Many layers are similar with small variations, due to gradual change ofthe principal strain directions. The top and bottom layers exhibits bendingstiffness, where most of the inner layers exhibit shear stiffness. It can benoticed that there is a horisontal bar of material in a layer just below themiddle. This is likely due to convergence to local minima, due to the ratheraggressive continuation approach used. It has been observed that these fea-tures are less likely to appear when robustness is not required, however, they38igure 22: Optimised robust layered macrostructure for the cantilever beamwith a concentrated load for M x = 24.do still appear at times and the problem increases with an increasing numberof layers. The problem appears to be due to the propagation of the designinformation during the optimisation procedure, where the topologies of thelocal microstructural details are determined at very different speeds. Thus,the topologies of some layers are uniquely determined by the surrounding,already formed, layers. The need for manufacturable microstructures of finite size is rather impor-tant. Current manufacturing processes, e.g. additive manufacturing, twophoton polymerisation and photolithography, do not allow for the manufac-ture of microstructural details with length scales far below that of the spec-imen size. Furthermore, if the microstructural details are spatially-varying,it is crucial that connectivity is ensured throughout the manufactured struc-ture in order to ensure the expected performance. It can be argued thatgoing from a structure with disconnected, but optimal, microstructural de-tails, as in [13] in the context of FMO, to modified connected microstructuraldetails ensures a manufacturable and well-performing structure. But when39uch modifications are made after the optimisation procedure itself, withno regard to the optimisation objective, the final structure may perform farfrom what was originally predicted. This may work to some extent for theforgiving compliance objective, but when moving to more advanced objec-tives or physics, the only rational approach is to include the requirement forconnectivity in the original optimisation itself.
It is important to note that small seemingly useless features can appear atthe interfaces of the microstructural cells when dealing with multiple blocks.These small features are artefacts resulting from the very tough connectivityconstraints imposed through the filtering procedure. In effect one imposesthat all horisontal interfaces between any two given coarse cells should be thesame. Generally it is observed that the microstructures that develop a cleartopology first, that is the ones with the largest sensitivities, have a dominantinfluence on the interfaces. This is as expected, but the filtering procedureensures connected microstructures that have as smoothly connected inter-faces as possible in order to ensure manufacturability. When the design isfurther relaxed into multiple layers with unique microstructures, these arte-facts do not occur because internal periodicity is not required in the verticaldirection.
In this article, the microstructural details have been forced to be periodic orin layers. The ideal situation is, of course, the ability of the microstructuraldetails to change in all directions. Forcing the design to have microstruc-tural details everywhere, essentially imposes a maximum length scale on theresulting members. Providing the possibility of unrestricted microstructuraldetails is trivial for the presented approach, as the optimisation problem re-verts back to the original unrestricted topology optimisation problems. Inthis case, all length scales above the minimum length scale imposed by the ro-bust approach are allowed in the design space and thus a significantly betterperforming structure is to be expected than when restricting the design spaceusing forced microstructural details. If one for instance wants microstructuraldetails everywhere, one can imposed a maximum length scale on the designusing local lower-bound constraints on the void volume fraction, similar to40he approach presented in [17].
Most multiscale finite element methods have been developed for use as anapproximate solver, rather than as a preconditioner. It is also possible to usethe presented spectral MsFEM as a solver for high-contrast linear elasticityproblems. As the error convergence studies presented in section 3.1 sug-gest, a lot of eigenmodes are needed to approximate the true solution well.Using only a few modes yields results similar to homogenisation and thusdesigns can end up being massively disconnected, as discussed in the origi-nal work [18]. By increasing the number of modes slightly, one can increasethe sensitivity of the structural response to disconnected members. How-ever, for microstructural design problems, these low-dimension bases requireeigenfunction derivatives for the optimisation to proceed correctly, which iscomputationally expensive. Hence, using the presented spectral MsFEM asa solver is left as a subject for future research.
For the presented problems, the addition of local constraints on the minimumvoid volume fraction in each unique microstructure cell is a cheap way toensure a maximum length scale, similar to [17]. Of course, the number of localconstraints increases with the amount of unique microstructures, that is thenumber of layers. This approach has been tested successfully, however, it doesnot show significant differences for the presented problems with relatively lowvolume fractions. The extension to other local constraints will be investigatedin future work.A full comparison of the presented approach to homogenisation for vary-ing microstructural details is part of current research, which will includeinvestigations on the possibility of using the presented methodology to findthe optimal design for cores of sandwich beam, plate and shell structures.The extension of the MsFEM-GMRES solver to three-dimensional generaltopology optimisation problems in a large scale parallel framework [27, 46] iscurrently being pursued. Here the computational performance of the methodis expected to excel due to the inherent parallisation properties of the decou-pled local problems. 41
Conclusion
The presented approach treats the optimisation of microstructural detailswith respect to the macroscopic response, while resolving all microstructuraldetails as compared to the homogenisation approach. This results in anoptimisation that takes localised behaviour, such as boundary conditionsand forces, into account and ensures microstructures tailored for a specificapplication rather than specific properties. The approach is fully flexible andcan treat microstructural details with and without clear separation of lengthscale compared to macrostructure. Finally, the approach yields robust andmanufacturable solutions, with smoothly connected and continuously varyingmicrostructural details in one direction. A further advantage of the proposedmethodology, is that the application of the MsFEM-GMRES solver allowsfor huge problems to be solved in
Matlab and that it exhibits contrast-independent behaviour ensuring fast computations for topology optimisationproblems.
The authors would like to thank Erik Andreassen for providing
Matlab subroutines for homogenisation and the TopOpt group for useful commentsand discussions during the preparation of the work.Both authors were funded by Villum Fonden through the NextTop project,as well as the EU FP7-MC-IAPP programme LaScISO.
References [1] O. Sigmund, Materials with prescribed constitutive parameters: Aninverse homogenization problem, International Journal of Solids andStructures 31 (17) (1994) 2313 – 2329. doi:10.1016/0020-7683(94)90154-6 .[2] O. Sigmund, Tailoring materials with prescribed elastic properties,Mechanics of Materials 20 (4) (1995) 351 – 368. doi:10.1016/0167-6636(94)00069-7 .[3] O. Sigmund, S. Torquato, Composites with extremal thermal expansioncoefficients, Applied Physics Letters 69:21 (1996) 3203–3205.424] C. S. Andreasen, O. Sigmund, Multiscale modeling and topology opti-mization of poroelastic actuators, Smart Materials and Structures 21 (6)(2012) 065005. doi:10.1088/0964-1726/21/6/065005 .[5] M. P. Bendse, O. Sigmund, Topology Optimization: Theory, Methodsand Applications, Springer, 2003, iSBN: 3-540-42992-1.[6] M. P. Bendsøe, N. Kikuchi, Generating optimal topologies in struc-tural design using a homogenization method, Computer Methods inApplied Mechanics and Engineering 71 (2) (1988) 197–224. doi:10.1016/0045-7825(88)90086-2 .[7] M. P. Bendsøe, Optimal shape design as a material distribution prob-lem, Structural optimization 1 (4) (1989) 193–202. doi:10.1007/BF01650949 .[8] H. Rodrigues, J. Guedes, M. Bendsoe, Hierarchical optimization of mate-rial and structure, Structural and Multidisciplinary Optimization 24 (1)(2002) 1–10. doi:10.1007/s00158-002-0209-z .[9] P. G. Coelho, P. R. Fernandes, H. C. Rodrigues, J. B. Cardoso, J. M.Guedes, Numerical modeling of bone tissue adaptation - a hierarchicalapproach for bone apparent density and trabecular structure, Journalof Biomechanics 42 (2009) 830–837.[10] P. G. Coelho, J. B. Cardoso, P. R. Fernandes, H. C. Rodrigues, Parallelcomputing techniques applied to the simultaneous design of structureand material, Advances in Engineering Software 42 (2011) 219–227.[11] P. B. Nakshatrala, D. A. Tortorelli, K. B. Nakshatrala, Nonlinear struc-tural design using multiscale topology optimization. Part I: Static for-mulation, Computer Methods in Applied Mechanics and Engineering261-263 (2013) 167–176.[12] L. Xia, P. Breitkopf, Concurrent topology optimization design of mate-rial and structure within nonlinear multiscale analysis framework, Com-puter Methods in Applied Mechanics and Engineering 278 (2014) 524–542. doi:http://dx.doi.org/10.1016/j.cma.2014.05.022 .4313] F. Schury, M. Stingl, F. Wein, Efficient two-scale optimization of man-ufacturable graded structures, SIAM Journal on Scientific Computing34 (6) (2012) B711–B733.[14] A. Radman, X. Huang, Y. Xie, Topology optimization of functionallygraded cellular materials, Journal of Materials Science 48 (4) (2013)1503–1510.[15] Y. Xie, Z. Zuo, X. Huang, J. Rong, Convergence of topological pat-terns of optimal periodic structures under multiple scales, Structuraland Multidisciplinary Optimization 46 (1) (2012) 41–50.[16] Z. H. Zuo, X. Huang, X. Yang, J. H. Rong, Y. M. Xie, Comparing opti-mal material microstructures with optimal periodic structures, Compu-tational Materials Science 69 (2013) 137–147.[17] J. Guest, Imposing maximum length scale in topology optimiza-tion, Structural and Multidisciplinary Optimization 37 (2009) 463–473,10.1007/s00158-008-0250-7.[18] B. S. Lazarov, Topology optimization using multiscale finite elementmethod for high-contrast media, in: I. Lirkov, S. Margenov, J. Was-niewski (Eds.), Large-Scale Scientific Computing, Lecture Notes in Com-puter Science, Springer Berlin Heidelberg, 2014, pp. 339–346. doi:10.1007/978-3-662-43880-0_38 .[19] I. Babu˘ska, G. Caloz, J. E. Osborn, Special finite element methods fora class of second order elliptic problems with rough coefficients, SIAMJournal on Numerical Analysis 31 (4) (1994) pp. 945–981.[20] I. Babu˘ska, J. E. Osborn, Generalized finite element methods: Theirperformance and their relation to mixed methods, SIAM Journal onNumerical Analysis 20 (3) (1983) pp. 510–536.[21] T. Hou, X. Wu, A multiscale finite element method for elliptic problemsin composite materials and porous media, Journal Of ComputationalPhysics 134 (1) (1997) 169–189.[22] Y. Efendiev, T. Y. Hou, Multiscale Finite Element Methods: Theoryand Applications, Springer, 2009.4423] M. Buck, O. Iliev, H. Andr¨a, Multiscale finite element coarse spaces forthe application to linear elasticity, Central European Journal of Math-ematics 11 (4) (2013) 680–701. doi:10.2478/s11533-012-0166-8 .[24] J. Galvis, Y. Efendiev, Domain decomposition preconditioners for mul-tiscale flows in high-contrast media, Multiscale Modeling & Simulation8 (4) (2010) 1461–1483. doi:10.1137/090751190 .URL http://link.aip.org/link/?MMS/8/1461/1 [25] Y. Efendiev, J. Galvis, X.-H. Wu, Multiscale finite element methods forhigh-contrast problems using local spectral basis functions, Journal ofComputational Physics 230 (4) (2011) 937 – 955. doi:10.1016/j.jcp.2010.09.026 .[26] D. Braess, Finite Elements: Theory, Fast Solvers, and Applications inSolid Mechanics, Cambridge University Press, 2007.[27] N. Aage, B. Lazarov, Parallel framework for topology optimization usingthe method of moving asymptotes, Structural and Multidisciplinary Op-timization 47 (4) (2013) 493–505. doi:10.1007/s00158-012-0869-2 .[28] O. Amir, N. Aage, B. Lazarov, On multigrid-CG for efficient topol-ogy optimization, Structural and Multidisciplinary Optimization 49 (5)(2014) 815–829. doi:10.1007/s00158-013-1015-5 .URL http://dx.doi.org/10.1007/s00158-013-1015-5 [29] P. S. Vassilevski, Multilevel Block Factorization Preconditioners:Matrix-based Analysis and Algorithms for Solving Finite Element Equa-tions, Springer, New York, 2008.[30] Y. Efendiev, J. Galvis, P. Vassilevski, Spectral element agglomerate al-gebraic multigrid methods for elliptic problems with high-contrast coef-ficients, in: Y. Huang, R. Kornhuber, O. Widlund, J. Xu (Eds.), DomainDecomposition Methods in Science and Engineering XIX, Vol. 78 of Lec-ture Notes in Computational Science and Engineering, Springer BerlinHeidelberg, 2011, pp. 407–414.[31] Y. Efendiev, J. Galvis, R. Lazarov, J. Willems, Robust domain decom-position preconditioners for abstract symmetric positive definite bilin-ear forms, ESAIM: Mathematical Modelling and Numerical Analysis 46(2012) 1175–1199. 4532] Y. Saad, M. Schultz, GMRES: A generalized minimal residual algorithmfor solving nonsymmetric linear systems, SIAM Journal on Scientific andStatistical Computing 7 (3) (1986) 856–869. doi:10.1137/0907058 .[33] M. Brezina, P. S. Vassilevski, Smoothed aggregation spectral element ag-glomeration AMG: SA- ρ AMGe, in: I. Lirkov, S. Margenov, J. Waniewski(Eds.), Large-Scale Scientific Computing, Vol. 7116 of Lecture Notes inComputer Science, Springer Berlin Heidelberg, 2012, pp. 3–15. doi:10.1007/978-3-642-29843-1\_1 .[34] K. Svanberg, The method of moving asymptotes - a new method forstructural optimization, International Journal for Numerical Methodsin Engineering 24 (2) (1987) 359–373. doi:10.1002/nme.1620240207 .[35] P. Michaleris, D. A. Tortorelli, C. A. Vidal, Tangent operators anddesign sensitivity formulations for transient non-linear coupled prob-lems with applications to elastoplasticity, International Journal forNumerical Methods in Engineering 37 (14) (1994) 2471–2499. doi:10.1002/nme.1620371408 .[36] J. K. Guest, J. H. Prvost, T. Belytschko, Achieving minimum lengthscale in topology optimization using nodal design variables and projec-tion functions, International Journal for Numerical Methods in Engi-neering 61 (2) (2004) 238–254. doi:10.1002/nme.1064 .[37] F. Wang, B. S. Lazarov, O. Sigmund, On projection methods, con-vergence and robust formulations in topology optimization, StructuralMultidisciplinary Optimization 43 (6) (2011) 767–784. doi:10.1007/s00158-010-0602-y .[38] B. Bourdin, Filters in topology optimization, International Journal forNumerical Methods in Engineering 50 (9) (2001) 2143–2158. doi:10.1002/nme.116 .[39] O. Sigmund, Morphology-based black and white filters for topology op-timization, Structural Multidisciplinary Optimization 33 (4-5) (2007)401–424. doi:10.1007/s00158-006-0087-x .[40] B. S. Lazarov, O. Sigmund, Filters in topology optimization based onhelmholtz-type differential equations, International Journal for Numer-ical Methods in Engineering 86 (2011) 765–781.4641] B. S. Lazarov, M. Schevenels, O. Sigmund, Topology optimization con-sidering material and geometric uncertainties using stochastic colloca-tion methods, Structural and Multidisciplinary Optimization 46 (2012)597–612.[42] J. K. Guest, A. Asadpoure, S.-H. Ha, Eliminating beta-continuationfrom heaviside projection and density filter algorithms, Structural andMultidisciplinary Optimization 44 (4) (2011) 443–453. doi:10.1007/s00158-011-0676-1 .[43] O. Amir, O. Sigmund, B. S. Lazarov, M. Schevenels, Efficient reanaly-sis techniques for robust topology optimization, Computer Methods inApplied Mechanics and Engineering 245/246 (2012) 217 – 231.[44] Y. Efendiev, J. Galvis, T. Y. Hou, Generalized multiscale finite elementmethods (gmsfem), Journal of Computational Physics 251 (0) (2013)116 – 135. doi:http://dx.doi.org/10.1016/j.jcp.2013.04.045 .[45] G. I. Rozvany, On symmetry and non-uniqueness in exact topology op-timization, Structural and Multidisciplinary Optimization 43 (3) (2011)297–317. doi:10.1007/s00158-010-0564-0doi:10.1007/s00158-010-0564-0