Multiscale structural optimization with concurrent coupling between scales

Abstract

A robust three-dimensional multiscale structural optimization framework with concurrent coupling between scales is presented. Concurrent coupling ensures that only the microscale data required to evaluate the macroscale model during each iteration of optimization is collected and results in considerable computational savings. This represents the principal novelty of this framework and permits a previously intractable number of design variables to be used in the parametrization of the microscale geometry, which in turn enables accessibility to a greater range of extremal point properties during optimization. Additionally, the microscale data collected during optimization is stored in a reusable database, further reducing the computational expense of optimization. Application of this methodology enables structures with precise functionally graded mechanical properties over two scales to be derived, which satisfy one or multiple functional objectives. Two classical compliance minimization problems are solved within this paper and benchmarked against a Solid Isotropic Material with Penalization (SIMP)–based topology optimization. Only a small fraction of the microstructure database is required to derive the optimized multiscale solutions, which demonstrates a significant reduction in the computational expense of optimization in comparison to contemporary sequential frameworks. In addition, both cases demonstrate a significant reduction in the compliance functional in comparison to the equivalent SIMP-based optimizations.

1 Introduction

The pursuit of optimum structural performance has led to the development of numerous structural optimization methodologies in recent decades. At present, structural optimization frameworks which seek to combine the potential of bespoke periodic microstructures (Huang et al. 2011; Andreassen et al. 2014; Xia and Breitkopf 2015; Xu et al. 2016) and the design flexibility of additive manufacturing (Vaezi et al. 2013) are attracting increasing attention, as they enable manufacturable structures to be optimized over two scales. Typically, frameworks of this nature assume multiscale approximations (Weinan et al. 2007) governed by homogenization theory (Bendsøe and Kikuchi 1988), to represent the mechanical properties of a periodic microstructure as a single point within the macroscale domain, as this reduces the computational expense of otherwise expensive analysis.

The extremal point properties afforded by bespoke periodic microstructures within multiscale optimization frameworks (Rodrigues et al. 2002; Schumacher et al. 2015; Sivapuram et al. 2016; Zhu et al. 2017; Li et al. 2018; Watts et al. 2019; Garner et al. 2019; Gao et al. 2019) represent an advantage over topology optimization (Bendsøe and Sigmund 2004), which is constrained to isotropic point material properties. This is an advantage shared by free material optimization (Kočvara et al. 2002, 2008), which uses all 21 unique elements of the elasticity tensor as design variables to derive theoretically optimal but non-manufacturable structures. The accessibility to extremal point properties should not be understated, as achieving equivalent local properties using topology optimization requires the heterogeneous assembly of isotropic material, which is heavily dependent on initialization conditions, optimization algorithm, mesh resolution, and regularization.

Many multiscale optimization frameworks opt for a parameterization-based strategy, where the design variables controlling the microstructural geometry are optimized to fulfil functional objectives at the macroscale (Cheng et al. 2018; Imediegwu et al. 2019; Colabella et al. 2019; Wang et al. 2020; Thillaithevan et al. 2020). Efficient parameterizations can provide clear accessibility to extremal point properties, as individual parameters can be configured to couple tightly with specific elements of the elasticity tensor, enabling straightforward control of point properties during optimization.

The parameterization-based multiscale optimization frameworks referenced within this section correspond to the state-of-the-art methodologies within this field. It is noted, however, that these frameworks all require a precomputed database of microstructures before optimization can be performed; so-called sequential frameworks. Precomputing a database is prohibitive, as it places a limit on either the number of design variables available to parameterize the microscale geometry or the fidelity of microscale simulations. The work presented within this paper seeks to circumvent the limitations imposed by precomputing a database, by coupling the micro- and macroscale models concurrently (Weinan et al. 2007). This permits a previously intractable number of design variables to be optimized at the microscale, enabling enhanced control over the point properties during optimization, by not requiring a precomputed database and instead assembling a partial database on-the-fly.

2 Overview

The objective of this paper is to present a three-dimensional parameterization-based multiscale optimization framework, embedded with a novel concurrent coupling strategy between scales, and demonstrate its efficiency in comparison to sequential parameterization-based methods (Imediegwu et al. 2019; Colabella et al. 2019). Concurrent coupling ensures that only the microscale data required to evaluate the macroscale model during each iteration of optimization is collected and represents the principal novelty of this framework. Additionally, the microscale data collected during optimization is stored in a reusable database, further reducing the computational expense of optimization. Through the application of this methodology, it is permissible to derive structures with precise functionally graded mechanical properties over two scales, satisfying one or multiple functional objectives. In addition, this framework prioritizes the generation of structures within the space of manufacturable designs.

The present framework is founded upon an adaptation of the fully anisotropic microscale model developed by Imediegwu et al. (2019), which has been expanded to accommodate three additional parameters in the microscale geometry definition. The additional design variables constitute the availability of 10 unique parameters per microstructure within the macroscale domain, enabling greater flexibility in the microscale geometry in comparison to the base microscale model. The microscale parameters serve as continuous design variables during optimization, directed by adjoint sensitivities (Plessix 2006), which enable the satisfaction of functional objectives.

To demonstrate the efficacy of the concurrent multiscale optimization framework detailed within this paper, two classical problems are presented and solved. To highlight the efficiency, the total computational resource expended in the derivation of all optimized solutions presented within this paper is compared to the requirements of sequential parameterization-based methods (Imediegwu et al. 2019; Colabella et al. 2019). Objective values derived from optimized multiscale solutions are also compared to equivalent values derived through monoscale SIMP-based topology optimizations, to understand the relative merit of the present framework in comparison to the de facto standard method for structural optimization.

3 Methods

Outlined within this section is a detailed description of the methodology underpinning the microscale model, macroscale model, concurrent coupling strategy, and optimization formulation. Structural considerations for the remainder of this paper are limited to the theory of linear elasticity, so small deformations are assumed at both scales. As a consequence, constituent nodes used in the construction of both the microscale and macroscale structural models are required to satisfy the constitutive law for linear elasticity

$$ \sigma_{ij} = E_{ijkl} \epsilon_{kl} $$

(1)

where σ_ij represents the Cauchy stress tensor, E_ijkl denotes the fourth-order elasticity tensor, and 𝜖_kl represents the infinitesimal strain tensor. Strain is calculated using the strain-displacement equations

$$ \epsilon_{ij} = \frac{1}{2} \left( \frac{\partial{u_{i}}}{\partial{x_{j}}} + \frac{\partial{u_{j}}}{\partial{x_{i}}} \right) $$

(2)

where u represents displacement and x denotes coordinate axes. In total, there are 6 independent equations relating strains and displacements with 9 independent unknowns. Further to these initial assumptions, static equilibrium is also assumed in the evaluation of both structural models

$$ \sigma_{ji,j}+ F_{i} = 0 $$

(3)

where F_i denotes the sum of the external forces acting on the design domain.

3.1 Microscale model

The microscale model serves as a callable function within the present framework, which accepts an array of input parameters describing a single microstructure geometry and after computation returns the corresponding homogenized mechanical properties. The periodic microscale geometry (illustrated in Fig. 1) is constructed using 7 parameterized structural members, each with a unique but fixed orientation and located within the boundaries of a unit cell.

Fig. 1

The parameterized microscale geometry used within the present framework, with all 10 unique parameters depicted

Comprising the 7 members are 4 diagonal members, connecting opposing vertices of the unit cell and 3 face members, connecting the centroids of opposing faces. The diagonal members have a cylindrical shape, individually parameterized by a normalized radius ($\hat {r}_{1}, \hat {r}_{2}, \hat {r}_{3}, \hat {r}_{4}$). Normalized radii ($\hat {r}_{\eta }$) are calculated by dividing the dimensioned radii (r_η) by the characteristic unit cell length (l) and are used as they enable a clearer interpretation of the optimized results. In contrast to the diagonal members, the face members have configurable elliptical cross sections, parameterized by normalized radius ($\hat {r}_{5},\hat {r}_{6}, \hat {r}_{7}$) and aspect ratio (e₅, e₆, e₇). The normalized radii and aspect ratios are used in combination to derive the two orthogonal axes ($\hat {r}_{\eta ,1},\hat {r}_{\eta ,2}$) defining the elliptical cross sections, using a piecewise continuous relationship defined ∀η ∈{5,6,7}

$$ \hat{r}_{\eta,1}(\hat{r}_{\eta},e_{\eta})= \begin{cases} \hat{r}_{\eta} e_{\eta}, & AR_{max} \leq e_{\eta} \leq 1 \\ \hat{r}_{\eta}, & 1 < e_{\eta} \leq 2 - AR_{max} \end{cases} $$

(4)

$$ \hat{r}_{\eta,2}(\hat{r}_{\eta},e_{\eta})= \begin{cases} \hat{r}_{\eta}, & AR_{max} \leq e_{\eta} \leq 1 \\ \hat{r}_{\eta} \left( 2 - e_{\eta} \right) , & 1 < e_{\eta} \leq 2 - AR_{max} \end{cases} $$

(5)

where AR_max ∈ (0,1] denotes the maximum permissible aspect ratio of an elliptical cross section. This formulation permits either orthogonal axis (r_5,y, r_5,z, r_6,x, r_6,z, r_7,x, r_7,y) to become the dominant axis; as illustrated in Fig. 2. The combination of normalized radius and aspect ratio is preferred over directly using the two orthogonal axes as parameters, as it enables a limit to be placed on the maximum aspect ratio of an elliptical cross section.

Fig. 2

The permissible aspect ratios of an elliptical cross section, defined as a function of AR_max

This parametrization of the microscale geometry is selected as it provides clear accessibility to a range of mechanical properties, while satisfying the requirements of the heterogeneous multiscale method (Weinan et al. 2007), which requires the geometry to be periodic and the mechanical property space to vary continuously as a function of the microscale parameters. In addition, this parametrization prohibits the formation of problematic fabrication features such as locally enclosed volumes. Through the precise configuration of the microscale parameters, orthotropic microstructures with efficient stiffness characteristics can be configured to resist a range of load cases, which are ideal for compliance minimization-based problems. Additionally, the parameters controlling the precisely orientated diagonal members enable direct manipulation of the shear-extension terms within the elasticity matrix, providing a mechanism to couple shear and normal forces, and therefore facilitating clear accessibility to anisotropic mechanical properties. Consequentially, in addition to optimizing for minimum compliance, this microstructure is well suited at targeting other structural optimization problems, such as obtaining a target deformation profile. It is important to note, however, that the target of this paper is to demonstrate the potential of concurrent coupling within a multiscale optimization framework, and not to focus on the optimality of the microstructure used with respect to the Hashin-Shtrikman upper bounds for bulk and shear moduli (Hashin and Shtrikman 1963). Furthermore, it should be also noted that it is permissible to interchange the microstructure used within this framework with another, so long as it is satisfies the requirements of the heterogeneous multiscale method.

The microscale model is initialized for all parameter configurations by discretizing a cubic domain into a structured mesh using hexahedral elements. A fixed mesh is preferred, as it aids with the implementation of periodic boundary conditions, due to the consistent discretization on opposing faces. Microstructures are assembled within the fixed domain by assigning one of two possible material phases to each constituent node. The primary phase is represented by an isotropic material with properties E = 209 GPa and ν = 0.3, and the secondary phase is represented by a quasi-void material with isotropic properties E = 10^− 7 GPa and ν = 0.3. The primary phase is assigned to a node if it lies on the interior or boundary of one or more structural members. To determine which nodes are assigned this phase, a node set (N_η) is derived for each member by computing the set of nodes (P) which lie on the interior or boundary, using the perpendicular distance separating the nodes and the member centerlines

$$ \begin{array}{@{}rcl@{}} N_{\eta} &= &\left\{ P \in \mathbb{R}^{3} \left| \frac{\left\| P \circ {\Xi}_{\eta} \times (P \circ {\Xi}_{\eta} - L_{\eta}) \right\|}{\left\| L_{\eta} \right\| } \leq \hat{r}_{\eta} \right. \right\}, \\ && \forall \eta \in \{1,2,3,\hdots ,7\} \end{array} $$

(6)

where L_η denotes the direction vector of a member ($\hat {r}_{\eta }$) and

$$ \begin{array}{@{}rcl@{}} {\Xi}_{\eta k} & = \left\{\begin{array}{ll} \hat{r}_{5,z} / \hat{r}_{5,y}, & \text{if } \hat{r}_{5,y} < \hat{r}_{5,z} \land \eta = 5 \land x_{k} = 2 \\ \hat{r}_{5,y} / \hat{r}_{5,z}, & \text{if } \hat{r}_{5,y} \geq \hat{r}_{5,z} \land \eta = 5 \land x_{k} = 3 \\ \hat{r}_{6,z} / \hat{r}_{6,x}, & \text{if } \hat{r}_{6,x} < \hat{r}_{6,z} \land \eta = 6 \land x_{k} = 1 \\ \hat{r}_{6,x} / \hat{r}_{6,z}, & \text{if } \hat{r}_{6,x} \geq \hat{r}_{6,z} \land \eta = 6 \land x_{k} = 3 \\ \hat{r}_{7,x} / \hat{r}_{7,y}, & \text{if } \hat{r}_{7,y} < \hat{r}_{7,x} \land \eta = 7 \land x_{k} = 2 \\ \hat{r}_{7,y} / \hat{r}_{7,x}, & \text{if } \hat{r}_{7,y} \geq \hat{r}_{7,x} \land \eta = 7 \land x_{k} = 1 \\ 1, & \text{otherwise } \end{array}\right. \\ & \forall k \in \{1,2,3\} \end{array} $$

(7)

denotes a nodal mapping matrix which enables the elliptical cross sections of the face members to be compensated for by applying transformations to the cubic domain, where x_k denotes a coordinate axis. The node sets are combined using a finite union, to establish which nodes are assigned the primary phase

$$ N_{P} = {\displaystyle \bigcup\limits_{\eta=1}^{7}N_{\eta}} $$

(8)

and conversely by symmetric difference, which nodes are assigned the secondary phase (N_S). Stiffness properties are allocated to each element within the cubic domain, as a function of the fraction of nodes assigned each phase

$$ E(n_{P},n_{S}) = n_{P} E_{P} +(1 - n_{S}) E_{S} $$

(9)

where E_P and E_S correspond to the stiffness properties of the primary and secondary phases respectively, and n_P and n_S correspond to the fraction of nodes assigned either the primary or secondary phase respectively. For a more in-depth presentation of nodal property allocation, the reader is directed toward Imediegwu et al. (2019).

The homogenized mechanical properties of assembled microstructures are derived through periodically constrained finite-element analyses, conducted using the open-source finite-element package FEniCS (Logg et al. 2012). Enforcing periodic boundary conditions (PBCs) enables the analysis of a single microstructure to be representative of an infinite repeating assembly. PBCs are implemented within the microscale model, by coupling the displacement of nodes on opposing faces and corners within the cubic domain. The displacement relations are derived by first considering an expression for the displacement field of a periodic structure (Xia et al. 2006)

$$ u_{i}({x}) = \bar{\epsilon}_{ik}{x}+ u^{*}_{i}({x}) $$

(10)

where $\bar {\epsilon }_{kl}$ denotes the macroscopic strain tensor, x denotes the position vector, and $u^{*}_{i}$ denotes a periodic function which modifies the linear displacement field ($\bar {\epsilon }_{ik}x_{k}$). Using (10), the displacement fields for two opposing faces within the cubic domain can be expressed as

$$ u^{+}_{i} = \bar{\epsilon}_{ik}x^{+}_{k} + u^{*}_{i} $$

(11)

$$ u^{-}_{i} = \bar{\epsilon}_{ik}x^{-}_{k} + u^{*}_{i} $$

(12)

where the + and − superscripts denote parallel opposing faces. By subtracting (12) from (11), it is possible to derive an equation relating the displacement fields on opposing faces

$$ u^{+}_{i} - u^{-}_{i} = \bar{\epsilon}_{ik}x^{+}_{k} - \bar{\epsilon}_{ik}x^{-}_{k} = \bar{\epsilon}_{ik}{\Delta} x_{k} $$

(13)

This relation is then used to generate constraint equations coupling node sets on opposing faces and corners during analysis. Satisfying (13) for all node sets during analysis ensures that the displacement field is continuous, with no overlaps or voids between adjacent microstructures.

For each microstructure, 6 periodically constrained finite-element analyses are conducted using a unique macroscopic unit-strain tensor for each case. Following the results of a mesh refinement study, the number of elements selected to discretize the cubic domain is 50³, as this enables a satisfactory balance between accuracy (< 5% error) and execution time to be obtained.

Homogenized stress ($\bar {\sigma }_{ij}$) and strain ($\bar {\epsilon }_{kl}$) quantities are derived from each analysis using volume-averaged homogenization

$$ \bar{\sigma}_{ij} = \frac{1}{|{\Omega}|}{\int}_{\Omega} \sigma_{ij} d{\Omega} $$

(14)

$$ \bar{\epsilon}_{kl} = \frac{1}{|{\Omega}|}{\int}_{\Omega} \epsilon_{kl} d{\Omega} $$

(15)

where Ω denotes the cubic analysis domain. In combination, the 6 cases provide sufficient information to derive a homogenized elasticity tensor, by applying the constitutive laws outlined at the beginning of Section 3 to assemble and solve a linear system. A single call of the microscale model takes on average 3 min to execute on a 3.5 GHz Dual-Core Intel Core i7 processor.

3.2 Macroscale model

The macroscale model is constructed using an assembly of N microstructures, individually defined by a row vector within the design variable matrix

$$ X_{\alpha \upbeta}= \left[\begin{array}{ccccccc} r_{1,1} & r_{2,1} & {\cdots} & r_{7,1} & e_{5,1} & e_{6,1} & e_{7,1} \\ r_{1,2} & r_{2,2} & {\cdots} & r_{7,2} & e_{5,2} & e_{6,2} & e_{7,2} \\ {\vdots} & {\vdots} & {\ddots} & {\vdots} & {\vdots} & {\vdots} & {\vdots} \\ r_{1,N} & r_{2,N} & {\cdots} & r_{7,N} & e_{5,N} & e_{6,N} & e_{7,N} \end{array}\right] $$

(16)

which has dimensions N × 10. Microstructures are represented within the macroscale finite-element model as singular points with homogenized mechanical properties. Similarly to the microscale model, the evaluation of the macroscale finite-element model is also computed using FEniCS. This approach is representative of the true macroscopic response of the system, so long as a separation of scales exists between the micro- and macroscales. To demonstrate this, an expansion of a two-scale displacement field in δ is considered (Hollister and Kikuchi 1992)

$$ \textbf{u}_{\delta }(\textbf{x}) = \textbf{u}(\textbf{x},\textbf{y}) = \textbf{u}_{0}(\textbf{x},\textbf{y})+\delta \textbf{u}_{1}(\textbf{x},\textbf{y}) + O(\delta^{2}) $$

(17)

where δ denotes the ratio of the micro- to macroscale characteristic lengths, u_δ denotes the true displacement field, u₀ denotes the macroscopic displacement field, u₁ and higher-order terms denote perturbations to the macroscopic displacement field, and $\vec x$ and $\vec y$ denote the macro- and microlevel coordinates respectively. This expansion implies that the microscopic displacement field varies δ^− 1 times quicker than the macroscopic displacement field. As a consequence, so long as the ratio of the characteristic lengths satisfies the following inequality constraint

$$ 0 < \delta \ll 1 $$

(18)

perturbations imposed by microstructures on the macroscopic displacement field can be considered negligible. A recent numerical study by Cheng et al. (2018) indicated that an order of magnitude is sufficient to achieve a separation of scales (δ < 0.1).

3.3 Concurrent coupling

The microscale model is coupled to the macroscale model concurrently during optimization, representing the principal novelty of this framework. This permits only the mechanical property data required to evaluate the macroscale model during each iteration of optimization to be collected. Additionally, the data collected during optimization is stored in a reusable database, for the benefit of successive iterations of optimization.

To facilitate this approach, in advance of the first optimization being undertaken, a full-factorial design of experiments (DOE) (Antony 2014) is generated for the microscale parameters, but not populated with simulation data. Nodes within the DOE represent unique permutations of the microscale parameters and in combination span the entire permissible range of design space. The DOE is constructed using the microscale parameter bounds and level intervals defined in Table 1.

Table 1 The individual microscale parameter bounds and level intervals used to generate the full-factorial DOE. In combination, this results in a total of 153,664,000 unique permutations of the microscale parameters

The bounds set for the normalized radii enable a broad range of mechanical properties to be obtained and permit both empty and fully dense unit cells to be configured, in contrast to ground-structure–based approaches (Hagishita and Ohsaki 2009). A higher normalized radius bound ($\hat {r}_{\eta ,5-7} = 0.49$) is afforded to each face member, in comparison to the diagonal members, to compensate for the loss in cross-sectional area resulting from the configurable elliptical cross sections. The maximum aspect ratio for an elliptical cross section is set to 5/3(AR_max = 0.6), to prevent highly eccentric members forming in either orthogonal axis.

The number of levels afforded to each dimension provides sufficient coverage of the mechanical property space using the minimum number of nodes in each dimension. In combination, the DOE contains 153,664,000 unique permutations of the microscale parameters, which is several orders of magnitude larger than the size of the databases assembled within contemporary sequential frameworks (Imediegwu et al. 2019; Colabella et al. 2019).

A full-factorial DOE is selected as the foundation of this method, as opposed to other less structured techniques, as the uniformly distributed nodes afford a number of benefits. The regular arrangement of nodes within the DOE permits rotational and reflective symmetries of the microstructure to be leveraged, enabling the mechanical properties of up to 24 geometrically similar microstructures to be derived from a single-parent microstructure simulation. Symmetries of the microscale geometry are derived by computing which combinations of rotations (90^∘ rotations about either the x, y or z axis) and reflections (defined about the X − Y, X − Z and Y − Z planes) uniquely permute the microscale geometry. Through the application of these symmetries, it is possible to derive the homogenized mechanical properties for all 153,664,000 microstructures within the database using only 6,823,320 calls of the microscale model. A further benefit of a full-factorial DOE is the ability to efficiently locate nodes within the DOE as a function of the microscale parameters. This offers an important reduction in the computational expense of an operation which is performed frequently during optimization and would otherwise require relatively expensive nearest-neighbor algorithms to be implemented.

In sequential parameterization based frameworks (Imediegwu et al. 2019; Colabella et al. 2019), the homogenized mechanical properties for all 6,823,320 parent microstructures would be collected prior to optimization, to facilitate the assembly of a database. In contrast, only the subset of microstructures which enable the evaluation of the macroscale model during each iteration of optimization is collected within the present framework. To derive this subset, microstructures within the DOE which form enclosing simplices (n-dimensional triangles) around the individual microscale parameter sets (X_α) are computed first

$$ {\Gamma} = \bigcup\limits_{\alpha=1}^{N} {\Phi}_{\alpha} $$

(19)

where Γ denotes the set of microstructures forming enclosing simplices, Φ_α denotes the n + 1 microstructures forming the enclosing simplex around the parameter set X_α and N denotes the number of microstructures within the macroscale domain. The algorithm used to determine which microstructures form the enclosing simplex is presented in detail within Appendix 1. A further subset (Λ) is subsequently derived, which contains the parent microstructures for all microstructures within Γ

$$ {\Lambda} = \bigcup\limits_{m=1}^{N_{\Gamma}} {\Psi} ({\Gamma}_{m}) $$

(20)

where N_Γ denotes the number of elements within the set Γ and the function Ψ denotes a child to parent mapping function, which is derived in advance of optimization by leveraging the reflective and rotational symmetries. The DOE is then examined to determine which microstructures within Λ have already been collected in previous iterations of optimization and which microstructures require simulation. After the mechanical properties for all microstructures within Λ have been collected, material property transformation matrices (Tsai 1966) are applied to the parent microstructures, to derive the mechanical properties for all of the microstructures within Γ. To reduce the storage requirements for the database, transformations are applied on-the-fly during optimization, so that only the mechanical properties of parent microstructures are stored.

A linear interpolant is then derived from each simplex for all 21 unique elements of the elasticity tensor (E_α) and also for volume-fraction (V_α), as a function of the continuous microscale parameters using simplex interpolation (Hemingway 2002)

$$ \begin{array}{@{}rcl@{}} E_{\alpha}(X_{\alpha}) =[ \text{vol}(\phi_{\alpha,2},\phi_{\alpha,3},\phi_{\alpha,4}, \hdots ,\phi_{\alpha,n+1}, X_{\alpha} ) E_{\phi_{\alpha,1}} & \\ + \text{vol}(\phi_{\alpha,1},\phi_{\alpha,3},\phi_{\alpha,4}, \hdots ,\phi_{\alpha,n+1} , X_{\alpha} ) E_{\phi_{\alpha,2}} & \\ + \text{vol}(\phi_{\alpha,1},\phi_{\alpha,2},\phi_{\alpha,4}, \hdots ,\phi_{\alpha,n+1} , X_{\alpha} ) E_{\phi_{\alpha,3}} & \\ {\vdots} & \\ + \text{vol}(\phi_{\alpha,1},\phi_{\alpha,2},\phi_{\alpha,3}, \hdots ,\phi_{\alpha,n} , X_{\alpha} ) E_{\phi_{\alpha,n+1}} ] & \\ / \text{vol}({\Phi}_{\alpha}) \end{array} $$

(21)

where ϕ_{α, n} denotes the n^th microstructure used in the construction of the enclosing simplex Φ_α and “vol” denotes the function

$$ \text{vol}({\Phi}_{\alpha} )= \frac{1}{n!} \left| \det \left[\begin{array}{cc} \phi_{\alpha,2} - \phi_{\alpha,1} \\ \phi_{\alpha,3} - \phi_{\alpha,1} \\ {\vdots} \\ \phi_{\alpha,n} - \phi_{\alpha,1} \\ \phi_{\alpha,n+1} - \phi_{\alpha,1} \\ \end{array}\right] \right| $$

(22)

which is used to calculate the volume enclosed by a simplex (Stein 1966). This approach is equivalent to modelling the response of the microscale model using a piecewise linear basis, which affords C⁰ continuity between adjacent simplicial interpolants. This is in contrast to other response surface methodologies such as radial or polynomial basis functions, which in general require fully populated databases prior to the generation of response surfaces. The limitation of using a linear interpolation-based approach is that the stiffnesses of microstructures with normalized radii close to the lower bound are overpredicted. However, the relative stiffnesses of these microstructures are negligible in comparison to microstructures with multiple larger normalized radii and therefore do not contribute significantly to the overall stiffness of the macroscale structure.

Simplex interpolation is favored over equivalent rectangular interpolation-based approaches, as fewer simulation nodes are required to derive an interpolant; n + 1 nodes are required to construct an enclosing simplex, whereas 2ⁿ nodes are required to construct an enclosing hypercube. This is due to the fact that a congruent simplex constructed from the nodes of a unit hypercube resolves 1/n! times the volume of the unit hypercube. However, there is no guarantee that the additional design space enclosed by the hypercube will be leveraged during optimization. Therefore, as microscale simulations are computationally expensive, it is more efficient to resolve the minimum amount of design space during optimization.

3.3.1 Concurrent coupling example (2D)

To demonstrate visually how a database could be developed over the course of an optimization, an example is presented and illustrated in Fig. 3. To permit clearer visualization, a two-dimensional microstructure with two configurable thickness parameters (t₁, t₂) is used within this example. Prior to optimization, a full-factorial DOE is generated using 6 uniformly distributed levels for both parameters. In total, the DOE contains 36 unique simulation nodes, for which the material properties can be derived using only 21 nodes, due to the existence of a 90^∘ rotational symmetry.

Fig. 3

A two-dimensional example is provided to demonstrate how a database may evolve over the course of an optimization, using the concurrent coupling technique detailed within Section 3.3. To permit clearer visualization, the path of only a single microstructure through the design space during optimization is illustrated

Again, to permit clearer visualization, the path of only a single microstructure through the design space is considered during optimization within this example. At the beginning of the optimization process, this microstructure is initialized with design variables [0.2,0.2]. An interpolant is then derived to determine the mechanical properties of this microstructure using the simulation data collected from the nodes [0.2,0.2], [0.4,0.2], and [0.4,0.4], as these microstructures form the enclosing simplex. By consequence, a further region of the design space is simultaneously resolved at no extra cost, as a result of the rotational symmetry which enables the mechanical properties for the microstructure [0.2,0.4] to be derived from the microstructure [0.4,0.2]. It is important to note that once a region of the design space has been resolved, it can be exhaustively searched within for a negligible computational expense.

Over the course of the next 3 iterations, the continuous design variables evolve by following sensitivities with respect to functional objectives. To facilitate the evaluation of the macroscale model during each iteration, additional microstructures are simulated to resolve the required regions of design space. By consequence, a reusable database is assembled during optimization, which can be leveraged to reduce the computational expense of successive iterations or entirely new optimization problems.

3.4 Optimization formulation

During the optimization process, the continuous microscale parameters contained with the design variable matrix (X) are varied to fulfil functional objectives at the macroscale. This approach enables structures with precise functionally graded mechanical properties over two scales to be derived.

Prior to the evaluation of the macroscale finite-element model during optimization, the design variables pass through a Helmholtz filter (Lazarov and Sigmund 2011) to promote smoothly varying microscale geometry and discourage “checkerboarded” geometries. The Helmholtz filter is implemented by solving the following partial differential equation

$$ -R^{2} \nabla^{2} \tilde{X} + \tilde{X} = X $$

(23)

where $\tilde {X}$ denotes the filtered design variable matrix and R denotes the filter length, which is set equal to two times the length of an element for the results presented within Section 4, as this filter length ensures continuous geometry is derived, which is a requirement of the heterogeneous multiscale method (Weinan et al. 2007).

Optimization problems within this framework are formulated as follows

$$ \begin{aligned} & \underset{X}{\text{Minimize:}} & & J(\tilde{X}) \\ & \text{Subject to:} & & K(\tilde{X})U = F \\ &&& V(\tilde{X}) - V_{f} \leq 0 \\ &&& X_{\alpha\upbeta} \geq X_{min,\upbeta}, \forall \upbeta \in \{1,2,3, \hdots , n \} \\ &&& X_{\alpha\upbeta} \leq X_{max,\upbeta}, \forall \upbeta \in \{1,2,3, \hdots , n \} \end{aligned} $$

(24)

where $J(\tilde {X})$ denotes the scalar valued objective function and “$K(\tilde {X})U = F$” denotes the linear elastic finite-element system constraining the physics of the macroscale model. A limit is placed on the availability of material within the macroscale domain through a volume-fraction constraint, where $V(\tilde {X})$ denotes the true volume-fraction and V_f denotes the target volume-fraction. Two further constraints are also applied to bound the design variables between the lowest (X_min,β) and highest (X_max,β) levels outlined within the full-factorial DOE (Table 1).

Sensitivities of the objective functional are calculated using a chain rule approach

$$ \frac{\partial J }{\partial X} = \frac{\partial \tilde{X}}{\partial X} \left[ \sum\limits_{i = 1}^{6} \sum\limits_{j = 1}^{6} \left( \frac{\partial J }{ \partial E_{ij}} \circ \frac{\partial E_{ij}}{\partial \tilde{X}} \right) \right] $$

(25)

where $\frac {\partial J}{\partial E_{ij}}$ is computed by solving an adjoint equation, derived using the open-source Python library dolfin-adjoint (Mitusch et al. 2019) which leverages high-level algorithmic differentiation to derive adjoint models. The accuracy of the sensitivities derived through adjoint methods is verified using Taylor convergence tests. The partial derivative $\frac {\partial E_{ij}}{\partial \tilde {X}}$ is computed by finite-differencing between the base nodes of the enclosing simplices and adjacent nodes in each dimension. By consequence of the piecewise linear nature of the simplicial interpolants, the derivative $\frac {\partial E_{ij}}{\partial \tilde {X}}$ may not be continuous across adjacent simplices.

A chain rule approach is used, as it places a limit on the number of adjoint equations to be solved per gradient evaluation to 21, irrespective of the number of design variables parameterizing the microscale geometry. As the design variables pass through a Helmholtz filter prior to the evaluation of the objective functional, the sensitivities are calculated with respect to the filtered design variables ($\tilde {X}$). It is therefore necessary to recover the unfiltered representation of the sensitivities to ensure compatibility with the unfiltered design variables (X) during optimization. As a result, the derivative of the Helmholtz operator, denoted by $\frac {\partial \tilde {X} }{\partial X}$, is applied as an additional step in the sensitivity derivation.

Optimization is performed using the primal-dual interior point algorithm with line search outlined within (Andrei 2017). This class of optimizer is preferred, as interior point algorithms typically outperform competing methods (sequential quadratic programming, method of moving asymptotes) when applied to structural optimization-based problems (Rojas-Labanda and Stolpe 2015). The second-order derivatives required to solve the reduced Karush-Kuhn-Tucker (KKT) system are calculated using L-BFGS (Liu and Nocedal 1989), as analytical derivatives of the Lagrangian are unavailable. The optimization is terminated when the objective function decreases by less than 0.01% for five successive iterations.

The methodology outlined within this section is distilled into a single flowchart and illustrated in Fig. 4.

Fig. 4

The flowchart outlining the concurrent multiscale structural optimization framework presented within this paper

4 Results

The framework presented within Section 3 is applied to solve two classical compliance minimization problems. The optimized results and computational expenditure are presented in detail for both problems. Additionally, the multiscale results are benchmarked against a SIMP-based topology optimization, to understand the relative merit of this approach to the de facto standard method for structural optimization. The initialization conditions are kept constant for both problems presented within this section and are detailed in Table 2.

Table 2 The initialization values used for all parameters at the beginning of optimization

4.1 Cantilever beam — compliance minimization

The first geometry considered is a cantilever beam with a uniformly distributed load applied to the top surface, as illustrated in Fig. 5. A 30% volume-fraction constraint is applied. Compliance in this instance is quantified by the homogenized mechanical properties

$$ J = {\int}_{{\Omega}_{TS}}{u^{T}_{\alpha}} K(\tilde{X}_{\alpha}) {u_{\alpha}} d{\Omega}_{TS} $$

(26)

where Ω_TS denotes the top surface of the cantilever geometry. The cantilever geometry is discretized into 24,000 elements using 60 × 20 × 20 divisions in the x, y, and z axes respectively. The optimized normalized radius distributions are illustrated in Fig. 6, with two thresholds applied to highlight the salient structures. One threshold $\tilde {\hat {r}}_{\alpha ,i} \geq 0.3$ is applied to highlight the dominant structures formed during optimization and the second (transparent) threshold $0 \leq \tilde {\hat {r}}_{\alpha ,1-7} < 0.3$ is applied to highlight where material is allocated within the design domain. The optimized aspect ratio distributions relating the orthogonal axes of the face members are illustrated in Fig. 7. To permit clearer visualization, these distributions are overlaid onto the primary normalized radius thresholds. Additionally, the convergence plot for this problem is illustrated in Fig. 8.

Fig. 5

The cantilever beam geometry targeted by compliance minimization, with the loads and boundary conditions highlighted

Fig. 6

The optimized normalized radius distributions for the cantilever compliance minimization problem. One threshold $\tilde {\hat {r}}_{\alpha ,i} \geq 0.3$ is applied to highlight the dominant structures formed during optimization and the second (transparent) threshold $0 \leq \tilde {\hat {r}}_{\alpha ,1-7} < 0.3$ is applied to highlight where material is allocated within the design domain

Fig. 7

The optimized aspect-ratio distributions for the cantilever compliance minimization problem. To permit clearer visualization, these distributions are overlaid onto the primary normalized radii thresholds ($\tilde {\hat {r}}_{\alpha ,1-7} \geq 0.3$)

Fig. 8

The convergence plot for the cantilever beam compliance minimization

As depicted in Fig. 6, material is primarily allocated near the root of the cantilever geometry. It is in this region of the cantilever that the largest principal stresses occur, as illustrated in Fig. 9. Therefore, by allocating material near the root, the resistance to bending is increased. Structural member 5 (Fig. 6e) receives the highest allocation of material and is distributed in the shape of an I-beam. The principal stresses in the flanges of the I-beam are oriented parallel to member 5, which enables this member to provide an efficient resistance to the defined load. This is an example of how a microscale parameter can be configured to couple tightly with specific elements of the elasticity tensor (E₁₁), enabling clear accessibility to efficient stiffness characteristics during optimization. Structural member 6 (Fig. 6f), which is orientated parallel to the force vector, forms an approximate I-beam shape also, but with smaller flanges and a thicker web section than member 5. The difference in shape is attributed to the orientation of member 6 aligning more closely to the principal stresses in the web, than at the flanges. The diagonal members (Fig. 6a–d) are primarily distributed at the root of the cantilever beam, as these members are able to provide an efficient resistance in regions where the principal stress directions transition between the web and flanges of the I-beam. Symmetries are noted about the X − Y plane at z = 0.5 for members 1 and 4 and members 2 and 3, as these member pairs are identically oriented in the X − Y plane. Structural member 7 (Fig. 6g) receives the lowest allocation of material, due to its orientation to the z-axis and therefore insensitivity to the compliance functional. As a result, this member is distributed in similar regions to the diagonal members, to also assist with the transition of the principal stress directions between the web and flanges.

Fig. 9

The maximum and minimum principal stresses for the optimized cantilever beam geometry, derived from the homogenized stress tensors. To permit clearer visualization, a threshold is applied to eliminate regions where the principal stresses are less than 4 kPa

It is apparent that the structural members are primarily distributed in regions where the principal stress directions are parallel to the respective member orientations. However, due to the finite number of members, it is not always possible to have a member directly oriented to the principal stresses. In this situation, a combination of members is used instead. By varying the two orthogonal axes defining each face member, it is possible to make these combinations more efficient. This is demonstrated in Fig. 7a for structural member 5, where the aspect ratio tends to the minimum permissible value (0.6) in the web section of I-beam, where the principal stresses are not directly aligned to any member, but in contrast remains unchanged from the initialization value (1.0) in regions where there is a direct alignment, for example in the flanges. The normalized radius parameters provide sufficient control over the point properties in regions where the principal stresses are directly aligned and therefore changes in the aspect ratios are not required. However, to compensate in regions where there is no direct alignment, changes to the aspect ratio enable the secondary component of stiffness to be efficiently increased. For structural member 5, this corresponds to increasing r_5,y relative to r_5,z, enabling the stiffness properties of this member to be enhanced in the y direction at the expense of the insensitive z direction for an equivalent cross-sectional area. A similar pattern is also noted for structural member 6 (Fig. 7b), where r_6,x is increased relative to r_6,z in the web section of the I-beam to increase stiffness in the x direction, while sacrificing stiffness in the insensitive z direction. Due to the insensitivity of structural member 7 (Fig. 7c) to the compliance functional, no significant changes in aspect ratio are noted from the initialization conditions (see Table 2).

The optimized parameter distributions, illustrated in Figs. 6 and 7, are combined to construct the physical representation of the macroscale structure, as illustrated in Fig. 10. The reconstruction is derived by populating the cantilever beam geometry with volumetric data describing the distance from each structural member surface. The volumetric data uses 25³ subdivisions for each unit cell. A marching cubes algorithm (Lewiner et al. 2003) is then used to reliably convert the volumetric data into a surface mesh. The number of unit cells per dimension is equal to 48 × 16 × 16 in the x, y, and z directions respectively and the corresponding manifold mesh is over 2.3 GB in size. The scale and resolution of the reconstructed geometry is selected to aid with visualization and is therefore coarser than the geometry which would be manufactured.

Fig. 10

Reconstruction of the optimized parameter distributions (see Figs. 6 and 7) into a single macroscale geometry

The computational expense incurred by microscale simulations over the course of this compliance minimization is illustrated in Fig. 11, by plotting the cumulative number of parent microscale simulations required after each iteration. In total, 32,490 parent microstructures are required to derive the optimized geometry. This corresponds to 0.48% of the total parent database, which represents a significant reduction in comparison to sequential frameworks which require a precomputed fully populated database. The number of parent microscale simulations required to perform the first 10 iterations of optimization is modest (364), as all microstructures are initialized with the same values and are therefore initially contained within a single simplex. This requirement then increases linearly at an average of 473 simulations per iteration until iteration 60, at which point the optimization begins to converge. As a result, a smaller proportion of the design space is explored in successive iterations; therefore, increasingly fewer parent microscale simulations are required.

Fig. 11

The cumulative number of parent microscale simulations required to derive the optimized cantilever beam geometry

Equivalent compliance minimizations are performed using the cantilever beam geometry, loads, boundary conditions, and discretization for the base 7-dimensional microstructure (Imediegwu et al. 2019) and a SIMP-based topology optimization with a penalization parameter of 3. As the smoothness requirements for topology optimizations are less stringent than for the present framework, a more typical Helmholtz filter length of 1 cell length (Lazarov and Sigmund 2011) is used instead. The compliance functionals for both frameworks are compared in Table 3. Through the application of the concurrent multiscale framework in combination with the present 10-variable microstructure parameterization (Fig. 1), a 27.4% reduction in the compliance functional is achieved in comparison to SIMP. This is a result of the clear accessibility afforded by the microscale parameterization to orthotropic microstructures with efficient stiffness characteristics, in comparison to the fixed isotropic material properties locally available within SIMP. A further reduction of 2.1% is also noted in comparison to the base 7-dimensional microstructure, demonstrating the increased accessibility to efficient stiffness characteristics permitted by the three additional design parameters. However, it should be reiterated that the principal objective of this paper is to demonstrate the computational savings afforded by a concurrent multiscale approach and not the optimality of the present microstructure with respect to the Hashin-Shtrikman upper bounds for bulk and shear moduli (Hashin and Shtrikman 1963).

Table 3 The optimized compliance functionals derived for the cantilever beam compliance minimizations, for both the present concurrent multiscale framework and a SIMP-based topology optimization. The 10 and 7 parameter optimizations correspond to the present and base microstructure respectively

4.2 Unit cube — compliance minimization

The second geometry targeted through compliance minimization using the derived framework is a unit-cube with four equal and uniformly distributed loads applied to four unique surfaces, constrained by a single Dirichlet boundary condition, as illustrated in Fig. 12. Compliance in this instance is quantified by the homogenized mechanical properties

$$ \begin{array}{@{}rcl@{}} J &=& {\int}_{{\Omega}_{1}}{u^{T}_{\alpha}} K(\tilde{X}_{\alpha}) {u_{\alpha}} d{\Omega}_{1} + {\int}_{{\Omega}_{2}}{u^{T}_{\alpha}} K(\tilde{X}_{\alpha}) {u_{\alpha}} d{\Omega}_{2}\\ &+& {\int}_{{\Omega}_{3}}{u^{T}_{\alpha}} K(\tilde{X}_{\alpha}) {u_{\alpha}} d{\Omega}_{3} + {\int}_{{\Omega}_{4}}{u^{T}_{\alpha}} K(\tilde{X}_{\alpha}) {u_{\alpha}} d{\Omega}_{4} \end{array} $$

(27)

where Ω_i denotes the 4 sub-domains illustrated in Fig. 12. The unit-cube geometry is discretized into 27,000 elements using 30 × 30 × 30 divisions in the x, y, and z axes respectively.

Fig. 12

The unit-cube geometry targeted by compliance minimization, with the loads and boundary conditions highlighted

To demonstrate the reusability of the microscale property data collected during optimization, a parametric study of volume-fraction is conducted. A total of 7 volume-fractions are investigated, starting at 20% and increasing linearly in 5% intervals before stopping at 50%. The optimized normalized radius distributions are illustrated for the 30% volume-fraction case only in Fig. 13, using equivalent normalized radius thresholds to those applied to the cantilever beam results, but with a lower bound for the primary threshold of 0.25. Similarly, the optimized aspect ratio distributions relating the orthogonal axes of the face members are illustrated for the 30% volume-fraction case also in Fig. 14. Additionally, the convergence plot for this problem is illustrated in Fig. 15.

Fig. 13

The optimized normalized radii distributions for the unit-cube compliance minimization problem, with a 30% volume-fraction constraint. One threshold $\tilde {\hat {r}}_{\alpha ,i} \geq 0.25$ is applied to highlight the dominant structures formed during optimization and the second (transparent) threshold $0 \leq \tilde {\hat {r}}_{\alpha ,1-7} < 0.25$ is applied to highlight where material is allocated within the design domain

Fig. 14

The optimized aspect-ratio distributions for the unit-cube compliance minimization problem, with a 30% volume-fraction constraint. To permit clearer visualization, these distributions are overlaid onto the primary normalized radii thresholds ($\tilde {\hat {r}}_{\alpha ,1-7} \geq 0.25$)

Fig. 15

The convergence plot for the unit-cube compliance minimization, with a 30% volume-fraction constraint

As depicted in Fig. 13, material is allocated for all members at the base of the unit-cube near the displacement constraint, as this is where the largest principal stresses occur, as illustrated in Fig. 16. Structural member 6 (Fig. 13f) receives the highest allocation of material, as it is oriented parallel to the two loads applied to the top surface of the unit-cube and therefore provides an efficient resistance and a load path to the boundary constraint. Structural member 5 (Fig. 13e) is also allocated a significant proportion of the available material, as it is oriented parallel to the two forces applied to the side surfaces (Y − Z) of the unit-cube and therefore also provides an efficient resistance. The four loads applied to the unit-cube are aligned to either the x or y axes and act in tension, which results in a compression of the unit-cube in the z axis, as demonstrated by the minimum principal stresses (Fig. 9a). As a result, structural member 7 (Fig. 13g), is distributed in regions where the minimum principal stresses are largest, as its orientation to the z-axis enables the compression forces to be efficiently resisted. The diagonal members (Fig. 13a–d) form structures in regions where both force vectors contribute significantly to the displacement field, as the principal stresses in these regions have components in all three coordinate axes. A 90^∘ rotational symmetry in the X − Z plane about the centroid of the base is also noted for the diagonal members, as the orientation of each member is suitable for a different combination of the two force vectors.

Fig. 16

The maximum and minimum principal stresses for the optimized unit-cube geometry, derived from the homogenized stress tensors. To permit clearer visualization, only the principal stresses in one half of the unit-cube domain are illustrated

The optimized aspect ratios defining the orthogonal axes of the face members display similar patterns to those noted for the cantilever beam compliance minimization. In regions where the principal stresses are parallel to a member’s orientation, the aspect ratio remains unchanged from the initialization value of 1 (see Table 2), but in regions where there is a misalignment, the aspect ratio is altered to increase stiffness in the secondary direction for an equivalent cross-sectional area. For example, the dominant axis for structural member 6 (Fig. 14b) near the top surface of the unit-cube is R_6,z, as this is where the strongest compression forces act. However, the dominant axis transitions to R_6,x nearer to the base of the unit-cube, as the loads applied to the side surfaces affect the displacement field more significantly within this region. Similarly, at the corners of the unit-cube, the dominant axis for structural member 5 (Fig. 14a) is R_5,y, as the secondary component of force acts in the y direction, as depicted by the maximum principal stresses (Fig. 16b). A similar pattern is also noted for structural member 7 (Fig. 14c), where the dominant axis is R_7,y at the top of the unit-cube, where the largest maximum principal stresses are aligned to the y axis.

The optimized parameter distributions, illustrated in Figs. 13 and 14, are combined to construct the physical representation of the macroscale structure, as illustrated in Fig. 17. The number of unit cells per dimension is equal to 24 × 24 × 24 in the x, y, and z directions, respectively, and the corresponding manifold mesh is over 2.6 GB in size. Again, the scale and resolution of the reconstructed geometry is selected to aid with visualization and is therefore coarser than the geometry which would be manufactured.

Fig. 17

Reconstruction of the optimized parameter distributions (see Figs. 13 and 14) into a single macroscale geometry

The individual computational expense incurred by microscale simulations for all 7 optimizations, each with a unique volume-fraction constraint, is illustrated in Fig. 18. The optimizations were conducted in ascending order of volume-fraction, starting at 20% and finishing at 50%. The microscale property data collected during the cantilever beam compliance minimization serves as the initial database for this parametric study.

Fig. 18

The number of parent microscale simulations required to derive the optimized geometries for each volume-fraction case within the parametric study. The optimizations were performed in ascending order of volume-fraction

The 20% volume-fraction case requires a total of 20,434 parent microscale simulations to derive the optimized geometry; however, 5,049 of those were already derived during the cantilever beam compliance minimization. As a result, the reusability of the database permitted a reduction in the computational expense incurred by microscale simulations of 24.7%. In general, excluding the 25% volume-fraction case, the number of parent microscale simulations required increased with volume-fraction, as all simulations used the same initialization conditions and therefore started at the same volume-fraction of 12%. Therefore, the optimizations with higher permitted volume-fractions have to cover more design space to reach the optimized solution. Over the course of the 7 optimizations, on average, 63.1% of the parent microscale data required is collected in a previous optimization. This represents a significant reduction in the computational expense of microscale simulations incurred during optimization and demonstrates the reusability of the collected data.

The results collected from the parametric study are used to derive a more coherent benchmark against a SIMP-based topology optimization. Equivalent optimizations are performed for all 7 volume-fractions, again using a SIMP algorithm with a penalization parameter of 3, identical discretization, and a Helmholtz filter length of 1 cell length. The compliance functionals derived by both the present framework and SIMP are compared and illustrated in Fig. 19. It is clearly demonstrated that the present concurrent multiscale optimization framework outperforms SIMP for all 7 volume-fractions, with a compliance functional on average 38.3% lower. Additionally, an approximate mass saving is also calculated for the present framework, by determining the volume-fraction required by SIMP to achieve an equivalent compliance for the 30% volume-fraction case. This is computed by intersecting a compliance isoline (see Fig. 19) with the compliance curve generated by the results from the SIMP parametric study. Approximately 25% more material is required for SIMP in comparison to the present framework. Similarly, to the cantilever beam compliance minimization, this advantage is a result of the clear accessibility to efficient orthotropic microstructures afforded by the microscale parameterization and the ability to perform optimization without penalizing the intermediate values.

Fig. 19

The compliance values derived from a parametric study of volume-fraction for the unit-cube, for both the present framework and a SIMP based topology optimization

4.3 Total computational expense incurred by microscale simulations

Over the course of the 7 unit-cube compliance minimizations and the single cantilever beam compliance minimization, a total of 126,122 parent microscale simulations were required to derive the optimized results. This corresponds to resolving 1.85% of the parent microstructure database and 0.09% of the total microstructure database. This represents a significant reduction in the computational expense in comparison to contemporary sequential methods, which require fully populated databases to be assembled prior to the first optimization. In addition, the average number of parent microscale simulations required per optimization for the present framework is only 38,456, which is lower than the requirement to assemble the databases in contemporary frameworks which have 3 fewer design variables, for example (Imediegwu et al. 2019) and (Colabella et al. 2019) require 40,817 and 41,990 simulations respectively. It should be noted, however, that the reduction in computational cost is a function of the sampling plan used to populate the microstructural design space (full-factorial DOE) and may scale differently for alternative sampling plans.

5 Conclusion

In this paper, a novel multiscale structural optimization framework with concurrent coupling between scales is successfully implemented. By embedding a concurrent coupling strategy, only the microscale data required to evaluate the macroscale model during each iteration of optimization is collected. The microscale data collected during optimization is stored in a reusable database to reduce the computational expense of successive iterations or entirely new optimization problems. This permits a previously intractable number of design variables (10) to be used in the parameterization of the microscale geometry, which enables accessibility to a greater range of extremal point properties during optimization.

To demonstrate the efficacy of this approach, two classical compliance minimization problems are solved and benchmarked against a Solid Isotropic Material with Penalization (SIMP)–based topology optimization, using a penalization parameter of 3 and an identical discretization. The first problem is a cantilever beam, subjected to a 30% volume-fraction constraint. A reduction in the compliance functional of 27.4% is noted in comparison to the equivalent SIMP case and only 0.48% of the total parent database is required to derive this result.

For the second compliance minimization problem, a parametric study of volume-fraction is undertaken using a unit-cube geometry, to demonstrate the reusability of the microscale data collected during the optimization process. A total of 7 volume-fractions are investigated between 20 and 50%. On average, 63.1% of the microscale data required to derive the optimized solution was collected in a previous optimization. In addition, the present framework outperformed SIMP for all 7 volume-fraction cases, with a compliance functional on average of 38.3% lower. The advantage of this framework over SIMP-based topology optimizations is the clear accessibility afforded by the microscale parameterization to orthotropic microstructures with efficient stiffness characteristics, in comparison to the fixed isotropic material properties locally available within a SIMP-based topology optimization.

In total, only 1.85% of the total parent database is required to derive all 8 optimized solutions presented within this paper. This represents a significant improvement in comparison to contemporary sequential methods, which require fully populated databases prior to undertaking optimization.

Appendix 1: Enclosing simplex function

To derive a simplex interpolant for the continuous microscale parameter set X_α, the microstructures within the full-factorial DOE which form an enclosing simplex around X_α need to be located first. This set is a subset of the microstructures which form the enclosing hypercube, which are trivial to locate, as a full-factorial DOE represents a collection of n-dimensional hypercubes. Initially, only two microstructures from the hypercube set are required, the microstructure closest to the origin (ζ_min) and conversely the microstructure furthest away from the origin (ζ_max). For example, to derive the microstructure, all elements within the microscale parameter set

$$ \begin{array}{@{}rcl@{}} X_{\alpha} = [0.22,0.3,0.11,0.41,0.03,0.42,0.19,0.72,1.31,0.87] \end{array} $$

are rounded down to the nearest DOE level (see Table 1)

$$ \begin{array}{@{}rcl@{}} \zeta_{min} = [0.21,0.28,0.07,0.35,0.00,0.42,0.14,0.60,1.20,0.80] \end{array} $$

and conversely the microstructure furthest away from the origin is derived by rounding up all elements to the nearest DOE level

$$ \begin{array}{@{}rcl@{}} \zeta_{max} = [0.28, 0.35,0.14,0.42,0.07,0.49,0.21,0.80,1.40,1.00] \end{array} $$

Then by performing an exhaustive binary search, all 2ⁿ microstructures within the hypercube set can then be found by permuting the minimum and maximum values derived for each microscale parameter.

The enclosing hypercube can be dissected into n! congruent simplices using n!/2!(n − 2)! bisecting hyperplanes, of which all pass through ζ_min and ζ_max by definition. As a consequence, all n! simplices within this hypercube are constructed using both ζ_min and ζ_max. To determine the remaining n − 1 microstructures defining the enclosing simplex, a hyperplane bisection algorithm is used. To simplify the bisection process, the relative position of each microscale parameter within the hypercube is first computed

$$ \begin{array}{@{}rcl@{}} \hat{X}_{\alpha} = (X_{\alpha} - \zeta_{min}) \oslash (\zeta_{max} - \zeta_{min}) \end{array} $$

(28)

where $\hat {X}_{\alpha }$ denotes the transformed microscale parameters and for the present example is equal to

$$ \begin{array}{@{}rcl@{}} [0.14,0.29,0.57,0.86,0.43,0.00,0.71,0.60,0.55,0.35] \end{array} $$

A bisecting hyperplane is defined by the equality of a unique combination of two microscale parameters. Depending on which side of the hyperplane the transformed parameters fall, a different set of microstructures can be eliminated from the hypercube set. An example is illustrated in Fig. 20, which depicts a hypercube, bisected by a hyperplane defined by the equality of the two transformed microscale parameters, $\hat {X}_{\alpha ,1}$ and $\hat {X}_{\alpha ,2}$.

Fig. 20

Hypercube bisection example, projected onto the $\hat {X}_{\alpha ,1}-\hat {X}_{\alpha ,2}$ plane for clearer visualization

The transformed vector of microscale parameters illustrated in Fig. 20 possesses an $\hat {X}_{\alpha ,2}$ value larger than its $\hat {X}_{\alpha ,1}$ value. As a result, the coordinate lies above the line $\hat {X}_{\alpha ,1} = \hat {X}_{\alpha ,2}$ and therefore the microstructures with transformed coordinates $(1, 0, ..., \hat {X}_{\alpha ,n})$ are eliminated from the hypercube set. This process is repeated for every hyperplane combination until only n + 1 microstructures remain in the hypercube set. The reverse transformation is then applied

$$ \begin{array}{@{}rcl@{}} X_{\alpha} = (\hat{X}_{\alpha} + \zeta_{min}) \circ (\zeta_{max} - \zeta_{min}) \end{array} $$

(29)

to recover the microscale parameters which form the enclosing simplex (Φ_α). Although this implementation enables the microstructures forming the enclosing simplex to be derived, it is computationally expensive in a high-dimensional space, as a large number of bisections are required.

To reduce the computational expense of this process, a novel method is presented which indirectly considers multiple bisecting hyperplanes simultaneously. The purpose of a bisecting hyperplane is to determine whether the transformed value of one parameter is larger than another. Instead of considering each hyperplane individually, by sorting the transformed vector of design variables into order of descending magnitude

$$ \begin{array}{@{}rcl@{}} [0.86,0.71,0.60,0.57,0.55,0.43,0.35,0.29,0.14,0.00] \end{array} $$

and recording the order of the sorted array

$$ \begin{array}{@{}rcl@{}} (\hat{X}_{\alpha,4}, \hat{X}_{\alpha,7}, \hat{X}_{\alpha,8}, \hat{X}_{\alpha,3}, \hat{X}_{\alpha,9}, \hat{X}_{\alpha,5}, \hat{X}_{\alpha,10}, \hat{X}_{\alpha,2}, \hat{X}_{\alpha,1}, \hat{X}_{\alpha,6}) \end{array} $$

all of the bisecting hyperplanes can be considered simultaneously. As the first design variable in the sorted array is larger than all other design variables by definition, it will always be incremented from the base point first, as all other permutations will be eliminated after bisection. This process is then repeated for all elements within the sorted array. In summary, to find the enclosing simplex, the base point of the enclosing hypercube (ζ_min) is incremented to the next DOE level in the order of the design variables within the sorted array, with each intermediate node serving as a microstructure within the enclosing simplex. For the microscale parameter set X_α, the microstructures forming the enclosing simplex are

$$ \begin{array}{@{}rcl@{}} {\Phi}_{\alpha} = \begin{array}{llllllllll} \{ {[0.21, 0.28, 0.07, 0.35, 0.00, 0.42, 0.14, 0.60 , 1.20 , 0.80 ]}, \\ {[0.21, 0.28, 0.07, 0.42, 0.00, 0.42, 0.14, 0.60 , 1.20 , 0.80 ]}, \\ {[0.21, 0.28, 0.07, 0.42, 0.00, 0.42, 0.21, 0.60 , 1.20 , 0.80 ]}, \\ {[0.21, 0.28, 0.07, 0.42, 0.00, 0.42, 0.21, 0.80 , 1.20 , 0.80 ]}, \\ {[0.21, 0.28, 0.14, 0.42, 0.00, 0.42, 0.21, 0.80 , 1.20 , 0.80 ]}, \\ {[0.21, 0.28, 0.14, 0.42, 0.00, 0.42, 0.21, 0.80 , 1.40 , 0.80 ]}, \\ {[0.21, 0.28, 0.14, 0.42, 0.07, 0.42, 0.21, 0.80 , 1.40 , 0.80 ]}, \\ {[0.21, 0.28, 0.14, 0.42, 0.07, 0.42, 0.21, 0.80 , 1.40 , 1.00 ]}, \\ {[0.21, 0.35, 0.14, 0.42, 0.07, 0.42, 0.21, 0.80 , 1.40 , 1.00 ]}, \\ {[0.28, 0.35, 0.14, 0.42, 0.07, 0.42, 0.21, 0.80 , 1.40 , 1.00 ]}, \\ {[0.28, 0.35, 0.14, 0.42, 0.07, 0.49, 0.21, 0.80 , 1.40 , 1.00 ]} \} \end{array} \end{array} $$

This method is extremely computationally efficient, as the only intensive process is the sorting of the array, which is in contrast to the large number of bisecting hyperplane calculations needed otherwise. During optimization, this calculation is required over one million times, so the reduction in execution time is significant.