Introduction
Exploring novel material architectures with extremal properties has been a constant quest in the field of material design and lightweight engineering. These developments have been further promoted by advances in additive manufacturing facilitating fabrication of functional materials with unprecedented complexity (Meza et al., Reference Meza, Zelhofer, Clarke, Mateos, Kochmann and Greer2015; Zheng et al., Reference Zheng, Smith, Jackson, Moran, Cui, Chen, Ye, Fang, Rodriguez and Weisgraber2016). As substitutions to trial and error and limited human intuition, topology optimisation methods have been shown to be powerful tools in designing novel materials in various applications (Sigmund, Reference Sigmund2000; Bendsøe and Sigmund, Reference Bendsøe and Sigmund2003; Osanov and Guest, Reference Osanov and Guest2016). Up to now, however, the field of architected materials has focused on individual material properties-like stiffness, yield or buckling strength but less so on the intricate trade-off between these properties.
Material stiffness and strength are fundamental properties for determining material load-bearing capacity as they measure the material’s deformation resistance and ultimate load-carrying capacity, respectively. The quest for optimal stiff and strong materials depends on many aspects, including the loading conditions and constituent base materials. Stiffness-optimal materials meeting theoretical upper bounds have been obtained via systematic design approaches (Sigmund, Reference Sigmund1995; Guest and Prevost, Reference Guest and Prevost2006; Huang et al., Reference Huang, Radman and Xie2011; Andreassen et al., Reference Andreassen, Lazarov and Sigmund2014; Berger et al., Reference Berger, Wadley and McMeeking2017). It has been shown that plate microstructures reach the Hashin–Shtrikman bounds (Hashin, Reference Hashin1962) in the low-volume fraction limit and remain within 10% of the theoretical upper bounds at moderate volume fractions. A few studies have focused on improving microstructure strength, including stress minimisation and buckling strength optimisation. It has been shown that the maximum von Mises stress can be reduced when taking stress into account during the microstructure optimisation procedure either as constraints or objectives (Collet et al., Reference Collet, Noël, Bruggi and Duysinx2018; Coelho et al., Reference Coelho, Guedes and Cardoso2019; Alacoque et al., Reference Alacoque, Watkins and Tamijani2021; Ferrer et al., Reference Ferrer, Geoffroy-Donders and Allaire2021). The stress singularity issue, where the stress at a point approaches infinity as the density at that point approaches zero, was handled by stress relaxation methods (Cheng and Guo, Reference Cheng and Guo1997; Duysinx and Bendsøe, Reference Duysinx and Bendsøe1998; Bruggi, Reference Bruggi2008; Le et al., Reference Le, Norato, Bruns, Ha and Tortorelli2010). Large numbers of local stress constraints can be aggregated to a global quantity using the p-norm, Kresselmeier–Steinhauser (KS) (Kreisselmeier and Steinhauser, Reference Kreisselmeier and Steinhauser1980) methods or by the augmented Lagrangian method (da Silva et al., Reference da Silva, Aage, Beck and Sigmund2021). Microstructure buckling strength optimisation using topology optimisation was first studied by Neves et al. (Reference Neves, Sigmund and Bendsøe2002), focusing on cell-periodic buckling modes. More recently, 2D and 3D material microstructures were systematically designed to enhance buckling strength based on linear material analysis. This approach evaluates effective material properties using the homogenisation method (Hassani and Hinton, Reference Hassani and Hinton1998) and material buckling strength under a given macro stress using linear buckling analysis (LBA) with Bloch–Floquet boundary conditions to capture all the possible buckling modes (Triantafyllidis and Schnaidt, Reference Triantafyllidis and Schnaidt1993; Thomsen et al., Reference Thomsen, Wang and Sigmund2018; Wang and Sigmund, Reference Wang and Sigmund2021). Both studies showed that the optimised materials possess several times higher buckling strength than their references at the cost of some stiffness degradation. Subsequent 2D and 3D experimental verifications have further verified the buckling superiority of the optimised microstructures and validated the linear material evaluation (Bluhm et al., Reference Bluhm, Christensen, Poulios, Sigmund and Wang2022; Wang et al., Reference Wang, Brøns and Sigmund2023).
So far, optimised microstructures have been designed considering material buckling or yield failure separately. However, yield-optimised microstructures tend to be vulnerable to buckling failure. On the other hand, buckling-optimised microstructures assume constituent base materials with relatively high yield strength (yield-strength-to-Young’s-modulus ratio,
$\sigma _1/E_1$
), for example, elastomers. For other base materials with low relative yield strengths, for example, steel, the buckling-optimised microstructures fail due to localised yield. Hence, the optimised strength superiority in both cases may significantly degrade in real applications. Furthermore, a previous numerical study showed that for a given microstructure topology, the failure mechanism switches from buckling- to yield-dominated failure as the volume fraction increases (Andersen et al., Reference Andersen, Wang and Sigmund2021). Hence, it is crucial to consider both failure mechanisms in the design procedure considering different volume fractions. Furthermore, it is essential to provide microstructure candidates with programmable properties working for various base materials and volume fractions fitting various applications.
In this study, we extend the stiffness/buckling studies from Thomsen et al. (Reference Thomsen, Wang and Sigmund2018) and Wang and Sigmund (Reference Wang and Sigmund2021) to also include yield strength. Yield strength and Young’s modulus share a monotonic relation in the optimised designs, as shown in the result section. Hence, microstructure will be designed systematically by maximising buckling strength with different yield strength bounds and considering different volume fractions. The optimised microstructure are further evaluated for different practical base materials ranging from low relative yield strength steel to high relative yield strength thermoplastic polyurethane (TPU).
The paper is organised as follows: Section ‘Optimisation problem for 2D architected materials with enhanced stiffness and strength’ summarises basic formulations to evaluate microstructure stiffness and strength, and formulates the optimisation problem for designing 2D microstructures with extremal properties. Section ‘Results’ validates the proposed approach and presents topology-optimised microstructure sets for different volume fractions and corresponding performances considering different base materials. Finally, Section ‘Conclusions’ concludes the study.
Optimisation problem for 2D architected materials with enhanced stiffness and strength
This section summarises the basic formulations for designing 2D architected materials with enhanced stiffness and strength using topology optimisation. The finite element method is combined with homogenisation theory and LBA to evaluate material properties (Cook et al., Reference Cook, Malkus and Plesha2002). To accurately represent the stress situation, we employ the incompatible elements from Wilson et al. (Reference Wilson, Taylor, Doherty and Ghaboussi1973) and Wilson and Ibrahimbegovic (Reference Wilson and Ibrahimbegovic1990), that is, the so-called
$Q_{6}$
element. Two additional so-called incompatible modes are considered to represent bending deformations accurately. The reader is referred to the work by Wilson et al. (Reference Wilson, Taylor, Doherty and Ghaboussi1973) and Wilson and Ibrahimbegovic (Reference Wilson and Ibrahimbegovic1990) for additional formulations for the
$Q_{6}$
elements.
Material stiffness and strength evaluations
Under small strain assumptions, the material stiffness and strength are evaluated using the homogenisation approach and LBA together with Bloch–Floquet theory to account for buckling modes at different wavelengths based on the periodic microstructure. Figure 1 summarises the material stiffness and strength evaluation procedure.
Flowchart for material stiffness and strength evaluations based a stiffness-optimal orthotropic material with a volume fraction of 0.2. (a) Homogenisation of a periodic material using a representative volume element (RVE), that is, microstructure. (b) Calculation of yield strength and illustration of the corresponding irreducible Brillouin zone (IBZ) using uniaxial compression. (c) Buckling band structure, buckling strength, and mode.
The microstructure is assumed of unit size here. The symmetry properties of the effective elasticity matrix are exploited to represent the equations in a more compact form using the abbreviation
$kl\rightarrow \alpha $
:
$11\rightarrow 1$
,
$22\rightarrow 2$
,
$(12,21)\rightarrow 3$
. The effective elasticity matrix is calculated by an equivalent energy-based homogenisation formulation (Sigmund, Reference Sigmund1995; Hassani and Hinton, Reference Hassani and Hinton1998) via
$$ \begin{align} \bar{D}_{ \alpha\beta} = \frac{1}{|Y|} \sum_{e=1}^N \int_{Y^e} \left(\tilde{\boldsymbol{\varepsilon}}_{\alpha} - \hat{\mathbf{B}}^e\boldsymbol{\chi}^e_{\alpha} \right)^T \mathbf{D}^e \left(\tilde{\boldsymbol{\varepsilon}}_{\beta} - \hat{\mathbf{B}}^e\boldsymbol{\chi}^e_{\beta} \right) dY,\nonumber \\ \mathbf{K}_0 \boldsymbol{\chi}_{\kern-1.5pt\alpha} = \mathbf{f}_{\alpha}, \quad \alpha=1,2,3,\\ \boldsymbol{\chi}_{\kern-1.5pt\alpha}|_{x=1}=\boldsymbol{\chi}_{\kern-1.5pt\alpha}|_{x=0}, \quad \boldsymbol{\chi}_{\kern-1.5pt\alpha}|_{y=1}=\boldsymbol{\chi}_{\kern-1.5pt\alpha}|_{y=0}.\nonumber \end{align} $$
Here,
$|Y|$
denotes the volume of the microstructure,
$\sum $
represents a finite element assembly operation over all N elements, the superscript
$ \left ( \right ) ^T$
denotes the transpose,
$\hat {\mathbf {B}}^e$
with a size of
$3\times 8$
is the condensed strain-displacement matrix of element e in the
$Q_{6}$
element formulation,
$\mathbf {D}^e$
is the elasticity matrix of the material in element e,
$\tilde {\boldsymbol {\kern-1.5pt\varepsilon }}_{\alpha }=\left [ \delta _{\kern-1.5pt\alpha \beta }\right ] $
denotes the three independent unit strain fields,
$\mathbf {f}_{\alpha }$
is the condensed equivalent load vector induced by the
$\alpha $
th unit strain, and
$\mathbf {K}_0$
is the global condensed elastic stiffness matrix. The detailed formulation of
$ \mathbf {K}_0$
and
$\mathbf {f}_{\alpha }$
can be found in Wilson and Ibrahimbegovic (Reference Wilson and Ibrahimbegovic1990) and Thomsen et al. (Reference Thomsen, Wang and Sigmund2018).
The effective material compliance matrix is
$\bar {\mathbf {C}}=\bar {\mathbf {D}}^{-1} $
. Under plane stress assumptions, the effective Young’s and bulk moduli are calculated using the effective elasticity or compliance matrices, stated as
$$ \begin{align} \bar{E}=\frac{1}{\bar{C}_{11}}, \quad \bar{\kappa} =\frac{ \bar{D}_{11}+\bar{D}_{12}}{2}. \end{align} $$
For a prescribed macroscopic stress state,
$\boldsymbol {\sigma }_0$
, the corresponding stress state of element e is obtained by the superposition of the three perturbation fields induced by the three unit strain fields in Eq. (1), expressed as
$$ \begin{align} \boldsymbol{\sigma}^e = \mathbf{D}^e \Big( \mathbf{I} - \hat{\mathbf{B}}^e \boldsymbol{X}^e \Big) \boldsymbol{\varepsilon}_0, \ \mathrm{with} \ \boldsymbol{\varepsilon}_0 = \bar{\mathbf{C}} \boldsymbol{\sigma}_0, \end{align} $$
where
$\boldsymbol {\sigma }_0= \boldsymbol {n}$
with
$\boldsymbol {n}$
indicating the loading direction,
$\boldsymbol {n}=[-1,0,0]^T$
for the uniaxial loading case,
$\boldsymbol {X}^e= \left [\boldsymbol {\chi }^e_{1}\ \ \boldsymbol {\chi }^e_{2}\ \ \boldsymbol {\chi }^e_{3}\right ] $
is a
$8\times 3$
matrix containing the element perturbation fields.
The yield strength of the microstructure,
${\sigma }_y$
, is determined by the ratio between the base material yield strength (
$\sigma _1$
) and the maximal von Mises stress in the microstructure, given as
$$ \begin{align} {\sigma}_y =\frac{{\sigma}_1 }{\max\limits_e(\sigma^e_{vm})}. \end{align} $$
Here,
$\sigma ^e_{vm}$
is the elemental von Mises stress, calculated by
$$ \begin{align} \sigma^e_{vm}=\sqrt{ \sigma^e_1\sigma^e_1 -\sigma^e_1\sigma^e_2+\sigma^e_2\sigma^e_2 +3\sigma^e_3\sigma^e_3 } =\sqrt{ \left( \boldsymbol{\sigma}^e\right)^T \boldsymbol{M}\boldsymbol{\sigma}^e}, \end{align} $$
where
$\boldsymbol {M}= \left [ \begin {array} {l l l} 1 & -1/2 & 0 \\ -1/2 & 1 & 0 \\ 0& 0 &3\end {array} \right ]. $
Based on the stress distribution in the microstructure, subsequent LBA is performed to evaluate the material buckling strength. Both short- and long-wavelength buckling is captured by employing the Floquet–Bloch boundary conditions in the LBA (Triantafyllidis and Schnaidt, Reference Triantafyllidis and Schnaidt1993; Neves et al., Reference Neves, Sigmund and Bendsøe2002; Meza et al., Reference Meza, Phlipot, Portela, Maggi, Montemayor, Comella, Kochmann and Greer2017; Thomsen et al., Reference Thomsen, Wang and Sigmund2018), stated as
$$ \begin{align} \left[ \boldsymbol{K}_0 + \lambda_h \boldsymbol{K}_{\sigma} \right]\boldsymbol{\phi}_h = \boldsymbol{0}, \\ \boldsymbol{\phi}_h|_{x=1}=e^{Ik_1}\boldsymbol{\phi}_h|_{x=0},\quad \boldsymbol{\phi}_h|_{y=1}=e^{Ik_2}\boldsymbol{\phi}_h|_{y=0}. \nonumber \end{align} $$
Here,
$\boldsymbol {K}_{\sigma }=\sum _e \boldsymbol {K}^e_{\sigma } $
is the stress stiffness matrix with
$\boldsymbol {K}^e_{\sigma } $
being the elemental stress stiffness matrix,
$I=\sqrt {-1}$
is the imaginary unit, the smallest eigenvalue,
$\lambda _1$
, is the critical buckling strength for the given wave vector (
$\boldsymbol {k}=[k_1,k_2]^T$
), and
$\boldsymbol {\phi }_1$
is the associated eigenvector.
The material buckling strength,
$\sigma _{c}$
, is determined by the smallest eigenvalue for all the possible wave vectors in Eq. (6), located in the first Brillouin zone,
$\lambda _{min}$
, that is,
$\sigma _{c}=\min \limits _{h,\boldsymbol {k}} \lambda _h$
. The associated eigenvector is the critical buckling mode. The first Brillouin zone is the primitive cell in reciprocal space (Brillouin, Reference Brillouin1953), spanning over
$k_j \in [-\pi ,\ \pi ], \ j=1,2$
. Previous studies have shown that the critical buckling mode can be captured by sweeping
$\boldsymbol {k}$
-vectors along the boundaries of the irreducible Brillouin zone (IBZ) (Geymonat et al., Reference Geymonat, Müller and Triantafyllidis1993; Thomsen et al., Reference Thomsen, Wang and Sigmund2018) that is determined by the shared symmetries between the microstructure geometry and the macroscopic stress state. In this study, we focus on designing a square cell in Figure 1(a) under uniaxial stress and the corresponding IBZs is illustrated in Figure 1(b). Figure 1(c) shows the corresponding buckling band diagrams calculated as customary along the boundaries of IBZ and critical buckling modes for the volume fraction of 20%. It is seen that buckling strength is determined by a global shear mode with
$\sigma _c/E_1=0.00060$
with
$\bar {E}/E_1=0.1074$
. From the computational point of view, it is more convenient to perform LBA via
$\left [ -\tau _h \boldsymbol {K}_0 - \boldsymbol {K}_{\sigma } \right ]\boldsymbol {\phi }_h = \boldsymbol {0} $
, where
$\tau _h =1/\lambda _h$
.
Design parameterisation
Based on the FEM discretisation, an element-wise constant physical design variable,
$\bar {\rho }_e$
, is employed to represent the material distribution in element e, with
$\bar {\rho }_e=1$
and
$\bar {\rho }_e=0$
representing base material and void, respectively. As in Thomsen et al. (Reference Thomsen, Wang and Sigmund2018), to suppress the spurious buckling modes associated with the low-stiffness elements, different interpolation schemes are employed for the elastic stiffness and stress stiffness using the solid isotropic material with penalisation (SIMP) model (Bendsøe and Sigmund, Reference Bendsøe and Sigmund1999). The
$\varepsilon $
-relaxed approach from Cheng and Guo (Reference Cheng and Guo1997) is employed for the von Mises stress interpolation to avoid the stress singularity phenomenon. The interpolation schemes are written as
$$ \begin{align} {E}_e = \left\{ \begin{array}{ll} \bar{\rho}^p_e(E_1-E_0)+E_0, & \text{for } \boldsymbol{K}_0, \mathbf{f}_{\alpha}, \\ \bar{\rho}^p_e {E}_1, & \text{for } \boldsymbol{K}_{\sigma}, \ \\ \bar{\rho}_e/\left( \varepsilon (1-\bar{\rho}_e)+\bar{\rho}_e \right) {E}_1, & \text{for } {\sigma}^e_{vm}. \end{array} \right. \end{align} $$
Here,
$E_1$
is the base material Young’s modulus,
$E_0=10^{-5} E_1$
represents void regions to avoid spurious modes located at the low-stiffness region,
$p=3$
is chosen as the standard penalisation factor.
$\varepsilon =0.002$
is chosen for the von Mises interpolation.
A hyperbolic tangent threshold projection is employed to generate physical design variables from the design variables,
${\rho }_e$
, to enhance the discreteness of the optimised design (Wang et al., Reference Wang, Lazarov and Sigmund2011). This is given as
$$ \begin{align} \bar{\rho}_e=\frac{\tanh{\left(\beta_1\eta\right)}+\tanh{\left(\beta_1\left(\widetilde{\rho}_e-\eta\right)\right)}}{\tanh{\left(\beta_1\eta\right)}+\tanh{\left(\beta_1\left(1-\eta\right)\right)}}, \end{align} $$
where
$\widetilde {\rho }_e$
is the filtered design variable calculated from the design variables,
$\boldsymbol {\rho }$
, using a PDE filter presented in Lazarov and Sigmund (Reference Lazarov and Sigmund2011). When
$\beta _1$
is big,
$\bar {\rho }_e\approx 1$
if
$\widetilde {\rho }_e>\eta $
and
$\bar {\rho }_e\approx 0$
if
$\widetilde {\rho }_e<\eta $
. Hence, the projection in Eq. (8) suppresses grey element density regions induced by the PDE filter, when
$\beta _1$
is big and ensures black–white designs when the optimisation converges. Moreover, it mimics the manufacturing process, and is used in the robust design formulation context (Wang et al., Reference Wang, Lazarov and Sigmund2011), where manufacturing errors are taken into accounts by choosing different thresholds,
$\eta $
. In this study, the maximum value of
$\beta $
is chosen to
$\beta =8$
allowing a small amount of grey regions in the design to avoid the mesh dependency of the maximum von Mises stress as discussed by da Silva et al. (Reference da Silva, Aage, Beck and Sigmund2021).
Design problem formulation
The KS function (Kreisselmeier and Steinhauser, Reference Kreisselmeier and Steinhauser1980) is employed to aggregate the elemental von Mises stress (
$\sigma ^e_{vm}$
) to represent the maximum von Mises stress; to aggregate the considered eigenvalues for given
$\boldsymbol {k}$
-vectors (
$\tau _h\left (\boldsymbol {k}_l \right )$
); to represent the material buckling strength; or to aggregate both quantities to represent both yielding and buckling strength, stated as
$$ \begin{align} KS\left(\kappa_1 {\sigma^e_{vm}}/{\sigma_1}, \kappa_2 \tau_h\left(\boldsymbol{k}_l \right) \right) &= \frac{1}{\zeta} \ln\left( \kappa_1 \sum_{e=1}^{N} e^{\zeta \sigma^e_{vm}/\sigma_1}+ \kappa_2 \sum_{l=1}^{n_h} \sum_{h=1}^{m_l}e^{\zeta \left( \tau_h\left( \boldsymbol{k}_l\right)\right) }\right), \qquad \kappa_1,\kappa_2 \in \left\lbrace 0;1 \right\rbrace. \end{align} $$
Here,
$\kappa _1=1$
and
$\kappa _2=0$
aggregate the elemental von Mises stress,
$\kappa _1=0$
and
$\kappa _2=1$
aggregate the eigenvalues for given
$\boldsymbol {k}$
-vectors,
$\kappa _1=1$
and
$\kappa _2=1$
to aggregate both quantities.
The optimisation problem for enhancing material stiffness and strength can be formulated to minimise a weighted value of the KS function in Eq. (9) and
$\bar {E}^{-1}$
, stated as
$$ \begin{align} \min\limits_{\boldsymbol{\rho}} \quad & \qquad \left( \gamma_1 KS\left(\kappa_1 \sigma^e_{vm}/\sigma_1 ,\kappa_2 \tau_h\left(\boldsymbol{k}_l\right) \right) + \left( 1-\gamma_1 \right)\bar{E}^{-1} \right) {E_1}, \qquad \kappa_1,\kappa_2 \in \left\lbrace 0;1 \right\rbrace, \nonumber \\ s.t. & \qquad \left[ -\tau_h \boldsymbol{K}_0 \left(\boldsymbol{k}_l \right) - \boldsymbol{K}_{\sigma} \left( \boldsymbol{k}_l\right) \right]\boldsymbol{\phi}_h = \boldsymbol{0}, \nonumber \\ & \qquad \mathbf{K}_0 \boldsymbol{\chi}_{ \alpha} = \mathbf{f}_{\alpha}, \quad \alpha=1,2,3, \nonumber \\ & \qquad {KS \left(\sigma^e_{vm}/{\sigma_1} \right)} {E_1} \leq 1/{\sigma}^*.\ \\ & \qquad \bar{E}{/E_1}\geq E^*, \nonumber\\ & \qquad f={\sum_e v_e \bar{\rho}_e}/ {\sum_e v_e} \leq f^*, \nonumber \\ & \qquad \boldsymbol{0} \leq \boldsymbol{\rho} \leq \boldsymbol{1}. \nonumber \end{align} $$
Here,
$\gamma _1$
is the weight,
$\gamma _1=0$
and
$\gamma _1=1$
represent stiffness and strength optimisation, respectively.
$E_1$
in the objective normalises the microstructure strength and stiffness with respect to the base material Young’s modulus.
$ {\sigma }^*$
is the normalised yield strength lower bound,
$ E^*$
is the normalised Young’s modulus lower bound,
$v_e$
is the volume of element e, and
$f$
and
$f^*$
are the actual and the prescribed upper bound of the volume fraction in the microstructure.
The sensitivity of a component,
${\bar {D}}_{ \alpha \beta }$
, in the effective elastic matrix
${\bar {\boldsymbol {D}}}$
with respect to
$\bar {\rho }^e$
, is written as
$$ \begin{align} \frac{\partial {\bar{D}}_{\alpha\beta} } {\partial \bar{\rho}^e} = \frac{1}{|Y|} \int_{Y^e} \left(\tilde{\boldsymbol{\varepsilon}}_{\kern-1.5pt\alpha} - {\hat{\mathbf{B}}}^e\boldsymbol{\chi}^e_{\kern-1.5pt\alpha} \right)^T \frac{\partial {\mathbf{D}}^e } {\partial \bar{\rho}^e}\left(\tilde{\boldsymbol{\varepsilon}}_{\beta} - {\hat{\mathbf{B}}}^e\boldsymbol{\chi}^e_{\beta} \right). \end{align} $$
The sensitivities of
$\bar {E}$
can be analytical derived using Eqs. (2) and (11).
The sensitivities of elemental von Mises stress (
$\sigma ^e_{vm}$
) with the respect to
$\bar {\rho }^e$
is calculated using the adjoint sensitivity analysis by
$$ \begin{align} \frac{\partial \sigma^e_{vm} } {\partial \bar{\rho}^e} =\frac{\left( \boldsymbol{\sigma}^e\right)^T\boldsymbol{M}}{\sqrt{\left( \boldsymbol{\sigma}^e\right)^T \boldsymbol{M}\boldsymbol{\sigma}^e}}\left[ \frac{\partial {\mathbf{D}}^e } {\partial \bar{\rho}^e} \Big( \mathbf{I} - {\hat{\mathbf{B}}}^e \boldsymbol{X}^e \Big) \boldsymbol{\varepsilon}_0+ {\mathbf{D}}^e \Big( \mathbf{I} - {\hat{\mathbf{B}}}^e \boldsymbol{X}^e \Big) \frac{\partial \boldsymbol{\varepsilon}_0 } {\partial \bar{\rho}^e}\right] + \sum_{\alpha=1}^3 \left( \boldsymbol{ \varphi}_ \alpha^e\right)^T \left[ \frac{\partial \boldsymbol{K}^e_{0} }{\partial \bar{\rho}^e} \boldsymbol{\chi}^e_{ \alpha} - \frac{\partial \mathbf{f}^e_{\kern-1.5pt\alpha}}{\partial \bar{\rho}^e} \right]. \end{align} $$
Here,
${ \partial \boldsymbol {\varepsilon }_0}/{\partial \bar {\rho }^e}$
can be directly derived using Eqs. (3) and (11), the adjoint vectors
$ \boldsymbol {\varPhi }=\left [\boldsymbol {\varphi _1},\boldsymbol {\varphi _2},\boldsymbol {\varphi }_3 \right ] $
are obtained by solving
$$ \begin{align} \boldsymbol{K}_{0} \boldsymbol{\varPhi}=\sum_e \left[ \frac{\left( \boldsymbol{\sigma}^e\right)^T\boldsymbol{M}}{\sqrt{\left( \boldsymbol{\sigma}^e\right)^T \boldsymbol{M}\boldsymbol{\sigma}^e}} {\mathbf{D}}^e {\hat{\mathbf{B}}}^e\right] ^T \left[ \boldsymbol{\varepsilon}_0\right]^T. \end{align} $$
Assuming that the eigenvector is normalised, as
$\left ( \boldsymbol {\phi }_h \right )^H \boldsymbol {K}_0 \boldsymbol {\phi }_h = 1 $
, the sensitivity of eigenvalue
$\tau _h$
with respect to
$\bar {\rho }^e$
can be obtained via the adjoint sensitivity analysis as described below (Thomsen et al., Reference Thomsen, Wang and Sigmund2018),
$$ \begin{align} \frac{\partial \tau_h}{\partial \bar{\rho}^e}=& {\left( \boldsymbol{\phi}^e_h\right)}^H \left[ -\tau_h \frac{\partial \boldsymbol{K}^e_{0} }{\partial \bar{\rho}^e} - \frac{\partial \boldsymbol{K}^e_{\sigma} } {\partial \bar{\rho}^e} \right] \boldsymbol{\phi}^e_h + {\left( \boldsymbol{\phi}_h\right)}^H \left[ - \frac{\partial \boldsymbol{K}_{\sigma} }{ \partial \boldsymbol{\varepsilon}_0} \frac{ \partial \boldsymbol{\varepsilon}_0}{\partial \bar{\rho}^e} \right] \boldsymbol{\phi}_h +\sum_{\alpha=1}^3\left( \boldsymbol{\psi}_ \alpha^e\right)^H \left[ \frac{\partial \boldsymbol{K}^e_{0} }{\partial \bar{\rho}^e} \boldsymbol{\chi}^e_{ \alpha} - \frac{\partial \mathbf{f}^e_{\kern-1.5pt\alpha}}{\partial \bar{\rho}^e} \right]. \end{align} $$
Here,
$\left ( \right ) ^H$
denotes the complex conjugate,
$\boldsymbol {\psi }_ \alpha $
is the adjoint vector corresponding to
$\boldsymbol {\chi }_{\kern-1.5pt\alpha }$
, which is obtained by
$$ \begin{align} \boldsymbol{K}_{0} \boldsymbol{\psi}_\alpha=\sum_e \left(\boldsymbol{\phi}^e_h\right)^H \left[ \frac{\partial \boldsymbol{K}^e_{\sigma} }{ \partial \boldsymbol{\chi}^e_\alpha} \right] \boldsymbol{\phi}^e_h. \end{align} $$
The reader is referred to the work by Thomsen et al. (Reference Thomsen, Wang and Sigmund2018) for the detailed calculation of
$ \frac {\partial \boldsymbol {K}^e_{\sigma } }{ \partial \boldsymbol {\chi }^e_ \alpha }$
. One of the advantages of using aggregation functions is the uniqueness of the gradient of eigenvalues, even when eigenvalues are repeated as stated by Torii and Faria (Reference Torii and de Faria2017).
The sensitivities of the objective and constraints with respect to a design variable,
${\rho }_e$
, are obtained using the chain rule. The optimisation problem is implemented in a flexible framework for large-scale topology optimisation (Aage et al., Reference Aage, Andreassen and Lazarov2015) using the Portable Extensible Toolkit for Scientific computation (PETSc; Balay et al., Reference Balay, Abhyankar, Adams, Brown, Brune, Buschelman, Dalcin, Eijkhout, Gropp, Kaushik, Knepley, McInnes, Rupp, Smith, Zampini, Zhang and Zhang2016) and Scalable Library for Eigenvalue Problem Computations (SLEPc; Hernandez et al., Reference Hernandez, Roman and Vidal2005). The design is iteratively updated using the Method of Moving Asymptotes (MMA; Svanberg, Reference Svanberg1987) based on the gradients of the objective and constraints. One case robust formulation in Wang and Sigmund (Reference Wang and Sigmund2021) is employed in this study, where the objective, yield, and stiffness constrains are evaluated on an eroded microstructure generated using
$0.5+ \Delta \eta $
with the volume constraint working on the dilated microstructure generated with a threshold of
$0.5- \Delta \eta $
. The volume constraint is updated every 20 iterations such that the volume constraint of the intermediate microstructure of
$\eta = 0.5$
is satisfied. In this study, we choose
$\Delta \eta =0.05$
.
Results
The proposed optimisation formulation is employed to design 2D square microstructures with 45
$^{\circ }$
-symmetry to achieve tunable stiffness and buckling strength response while considering three different volume fractions, that is,
$f^*=0.2,0.1, 0.05$
. The unit cell is discretised by
$512\times 512\ Q_{6}$
elements for
$f^*=0.2$
and the resolution is doubled in order to enable the evolution of structural hierarchy for
$f^*=0.1$
and
$f^*=0.05$
. The filter radius is chosen to
$r=0.03$
for
$f^*=0.2$
and it is reduced to
$r=0.01$
for the two lower-volume fractions. A square cell with 45
$^{\circ }$
-symmetry is designed. The initial microstructures (starting guesses) are chosen as the Young’s modulus optimal microstructures with the prescribed volume fractions as shown in Figure 1. The presented designs are the blueprints with
$\eta =0.5$
unless otherwise stated. Smooth and sharp boundaries are extracted as contour lines from the optimised greyscale designs using the in-house code mentioned above, and subsequently evaluated in COMSOL using body-fitted meshes. The considered based material in the optimisation is pyrolytic carbon (PC) with a Poisson’s ratio of
$\nu =1/3$
and a relative yield strength of
$\sigma _1/E_1=0.044$
(Crook et al., Reference Crook, Bauer, Guell Izard, Santos de Oliveira, Martins de Souzae Silva, Berger and Valdevit2020).
Optimised microstructures with maximised Young’s modulus without (a) and with a yield strength constraint (b). Left: Optimised microstructures; Middle: von Mises stress distribution evaluated at the element centre using the in-house code with the optimised greyscale designs; Right: von Mises stress distribution evaluated using the extracted designs in COMSOL with body-fitted meshes. (c) Differences between the optimised designs. (d) Performance summary of the two designs.
Validation of the proposed approach
In the first case, we maximise microstructure Young’s modulus without or with the yield strength constraint by setting
$\gamma _1 =0$
. The corresponding normalised yield strength lower bound is set
${\sigma ^*=1/400}$
. The optimised designs are extracted using the 0.5 isocontour and imported into COMSOL for post-evaluation using second-order body-fitted meshes. Figure 2(a,b) presents the optimised microstructures obtained without and with the yield constraint, respectively. The left panels show the optimised microstructures, whereas the middle and right panels present the von Mises stress distribution evaluated using the in-house code and COMSOL. Both in-house code and COMSOL predict similar von Mises profiles and the same high von Mises stress regions. The intermediate physical densities in the optimised designs in the in-house code make the microstructure’s perform slightly weaker (that is, lower effective Young’s modulus) than the ones from COMSOL. Hence, the microstructures undergo bigger pre-strain as indicated by Eq. (3) and exhibit larger maximum von Mises stress compared to the post-evaluation in COMSOL. These conclusions are verified by the von Mises stress distributions in both designs shown in Figure 2(a,b).
The rather insignificant difference between the two designs is presented in Figure 2(c) using the 0.5 isocontours of the 1/8 cell at the upper right corner. The optimised microstructures possess very similar shapes with small deviations at the inner corner. The inner corner in the optimised microstructure with the yield constraint exhibits slightly lower curvatures to fulfil the yield strength constraint. The properties of the optimised microstructures evaluated from the in-house code and COMSOL are summarised in Figure 2(d). In the post-evaluation, both designs possess the same Young’s modulus with 2% to the theoretical Young’s modulus upper bound,
$E_{HS}=f/(2-f)E_1$
. The material buckling evaluations of both designs show that global shear modes dominate buckling failure, and both designs possess higher buckling strength than the initial design in Figure 1. The buckling strength enhancement is attributed to the round inner corner that makes the effect length of the member shorter than the initial design. Comparison between yield and buckling strengths shows that optimised microstructures fail due to yield.
The microstructure is subsequently optimised by maximising buckling strength with the same yield strength constraint, which results in the same topology due to the limited design freedom provided by the yield strength bound.
In the rest of the paper, all the designs and corresponding performance will be presented using the well-defined post-evaluation in COMSOL, unless otherwise stated.
Strength-optimised microstructure
In the second case, a microstructure is optimised to maximise the material strength, that is, maximising the minimum value between the buckling and yield strengths by choosing
$\kappa _1=1$
,
$\kappa _2=1$
, and
$\gamma _1 =1$
. Figure 3 presents the optimised design and corresponding performance. Figure 3(a) shows the optimised
$2\times 2$
cells where the red box highlights the optimised microstructure. Figure 3(b,c) presents the corresponding critical buckling mode and von Mises stress distribution with titles showing the corresponding buckling and yield strengths. As in the previous designs, the critical buckling mode is a global shear mode. Compared to them, the optimised design develops hierarchy to enhance the microstructure buckling strength. The hierarchy in the design enhances the effective width-to-length ratio of the members and leads to a higher buckling strength. On the other hand, it reduces member stiffness and leads to a lower effective Young’s modulus and yield strength. The optimised microstructure exhibits an effective Young’s modulus of
$\bar {E}/E_1= 0.08533$
, which is 21.9% degradation compared to the stiffness-optimal design. For the considered base material, PC, the corresponding buckling and yield strengths manifest that the optimised microstructure fails simultaneously due to buckling and yield failures.
Optimised microstructure with maximised material strength. (a)
$2\times 2$
optimised cells, the red region highlights the optimised microstructure. (b) Buckling failure mode with title representing buckling strength. (c) von Mises stress distribution under uniaxial compression with title representing the yield strength.
Optimised microstructures with tunable stiffness/yield and buckling strength
We optimise microstructure sets by maximising material buckling strength with different yield strength bounds for different volume fractions by choosing
$\kappa _1=0$
,
$\kappa _2=1$
, and
$\gamma _1 =1$
. The optimised design set for
$f^* =0.2$
is presented in Figure 4(a). As the required yield strength increases, the optimised microstructure becomes less hierarchical and converges to the stiffness/yield-strength-optimal one as shown in Figure 2(b) for a large yield strength constraint. The optimised microstructures’ performance is evaluated using different base materials ranging from low relative yield strength steel to high relative yield strength TPU. Table 1 summarises the considered base materials and corresponding properties, including Young’s modulus
$E_1$
, density
$\rho _1$
, and relative yield strengths, that is, the yield-strength-to-Young’s-modulus ratio
$\sigma _1/E_1$
.
Considered constituent base materials and corresponding properties (Crook et al., Reference Crook, Bauer, Guell Izard, Santos de Oliveira, Martins de Souzae Silva, Berger and Valdevit2020; Andersen et al., Reference Alacoque, Watkins and Tamijani2021)
Optimised microstructures for
$f^*=0.2$
. (a) Optimised microstructures with increased yield strength bounds from left to right. (b) Yield/buckling strength versus Young’s modulus performances for different constituent base materials. Black line with circles: Buckling strength versus Young’s modulus curve. The other coloured lines: Yield strength versus Young’s modulus. Aqua: TPU; purple: PC-Nano; teal: PC; yellow: epoxy; red: steel. (c) Strength versus Young’s modulus performances.
Figure 4(b) shows the normalised buckling strength (
$\sigma _c$
, black line with circles) and yield strengths for the considered base materials (
$\sigma _y$
, coloured lines with crosses, each colour represents one base material) versus Young’s modulus performance. Figure 4(c) presents the optimised microstructures’ strength versus Young’s modulus performance. For the considered base material in the optimisation, that is, PC, represented by the teal lines, it is seen that the increase of the yield strength bound leads to optimised microstructures with increased yield strength and Young’s modulus and reduced buckling strengths. The monotonic relation between the effective Young’s modulus and yield strength is observed in the microstructure set. It is explained by Eq. (3) that the higher effective Young’s modulus leads to a lower equivalent global strain for a given stress situation hence resulting in higher yield strength. For PC, the optimised design for material strength predicted by the optimised microstructure set is represented by the intersection between the buckling (the teal line with crosses) and yield strength (the black line with circles) curves. It is seen from Figure 4(b) that the strength-optimised design in Figure 3 possesses a material strength very close to the intersection point between the buckling curve (black) and the yield curve (teal), hence, it provides the best combined yield and buckling strength for PC.
The yield strength of the optimisation set for TPU represented by the aqua-coloured line is much bigger than the covered strength range in the plot; hence, it does not appear in Figure 4(b). For TPU, all the optimised designs are dominated by buckling failure due to its high relative yield strength, and all the optimised designs are dominated by yield failure for steel due to its low relative yield strength. Failure mechanism transition between yield and buckling is observed for the rest of the considered base materials. As the base material’s relative yield strength increases, the transition between yield-dominated to buckling-dominated failure shifts to the left, that is, to lower Young’s modulus microstructures. The strength-optimal microstructure switches to the optimised microstructure with a lower Young’s modulus as the relative yield strength of the base material increases as seen in Figure 4(c).
To further explore the optimised microstructure failure at the lower-volume fractions, the microstructures are optimised for
$f^*=0.1$
and
$f^*=0.05$
by maximising microstructure buckling strength for different yield strength bounds.
Figure 5 presents the optimised microstructure set for
$f^*=0.1$
with increased yield strength bounds and the corresponding performances for different based materials. The geometrical evolutions in the design set show that the hierarchy vanishes as the required yield strength bound increases and the optimised microstructure configuration converges to the stiffness-optimal one for the largest yield strength bound. This observation is the same as in the previous case.
Optimised microstructures for
$f^*=0.1$
. (a) Optimised microstructures with increased yield strength bound from left to right. (b) Buckling/yield strength versus Young’s modulus performance considering different base materials. (c) Microstructure strength versus Young’s modulus performances.
Unlike the previous case, the optimised microstructure set fails due to buckling for PC-Nano and TPU. The failure mechanism transits from the yield-dominated to the buckling-dominated failure for the rest of the base materials when the optimised microstructure evolves from the buckling-optimal microstructure at the left to the stiffness/yield-optimal one at the right in Figure 5(a) . Compared to the design set for
$f^*=0.2$
, the transition point shifts further to a lower Young’s modulus range with more hierarchy for
$f^*=0.1$
for the same base material.
Figure 6 summarises the optimised microstructure set for
$f^*=0.05$
and corresponding performances with different base materials. A similar geometrical evolution is observed as in the previous two design sets. The optimised microstructure set fails due to buckling failure for PC, PC-Nano, and TPU. Only the buckling-optimised design with the lowest yield strength bound fails due to yield strength for epoxy, and the rest designs fail due to buckling. The failure mechanism transition is observed for steel. However, buckling failure dominates only the stiffness-optimal design, whereas yield failure dominates the rest designs. Compared with the previous two design sets for higher-volume fractions, the transition between yield- and buckling-dominated failure shifts further towards the buckling-optimised designs.
Optimised microstructures for
$f^*=0.05$
. (a) Optimised microstructures with increased yield strength bound from left to right. Pink boxes highlight the regions for enlarged views at the bottom. (b) Buckling/yield strength versus Young’s modulus performances considering different base materials. (c) Microstructure strength versus Young’s modulus performances.
The overall comparison between the three microstructure sets shows that the optimised microstructures are mainly dominated by buckling failure for the low-volume fraction region, and they are dominated by yield failure for the high-volume fraction region. It can be expected that all the optimised designs, from the buckling-optimal design to the stiffness/yield-optimal design, fail only due to buckling at an even lower-volume fraction, and all the optimised designs have yield failure at much higher-volume fractions. The low critical volume fraction is determined when the full hierarchical buckling optimal microstructure fails due to buckling. The high critical volume fraction is determined when the corresponding stiffness-optimal microstructure exhibits yield failure. The critical volume fractions are determined by the relative yield strength of the base material. For a higher relative yield strength, the critical low-volume fraction is higher than the one with a lower one, whereas the high critical volume fraction is also higher.
To provide additional insight into the role of the choice of different base materials, Figure 7 summarises the non-normalised material properties versus mass density for the optimised microstructures. PC-Nano possesses the highest Young’s modulus and high relative yield strength (see Table 1), hence, the optimised microstructures using PC-Nano exhibit the highest effective Young’s modulus and strength while TPU results in the lowest effective Young’s modulus and strength due to its low Young’s modulus. As expected, the maximum Young’s modulus curve of the optimised microstructures matches the Hashin–Shtrikman bound,
$E_{HS}=f/(2-f)E_1=\rho /(2\rho _1-\rho )E_1$
; see the left panel.
Non-normalised properties of the optimised microstructures. Left: Young’s modulus versus mass density with the dashed curve showing the Hashin–Shtrikman bound. Right: Strength versus mass density. The slopes indicate the polynomial order of the maximum strength lines to mass density.
The right panel presents the microstructure strength versus mass density, where
$\times $
’s indicate microstructures failed due to buckling, whereas
$\blacktriangle $
’s indicate microstructures failed due to yield. The dashed lines represent the maximum strength obtained from the optimised microstructures. By assuming a strength–mass–density relation of
$ \min (\sigma _c,\sigma _y) = c_0 f^{n_0}=c_0 \left (\frac {\rho }{\rho _1} \right ) ^{n_0}$
, the order
$n_0$
is obtained using the two points with lower densities and the orders for the different base materials are indicated in the strength plot as the line slopes. At low mass densities, the analytical yield strength upper bound follows the same relation as Young’s modulus, that is,
$\sigma _y=f/(2-f)\sigma _1=\rho /(2\rho _1-\rho )\sigma _1$
, represented by the solid lines in the plot and it possesses a linear relation with respect to mass density with
$n_0=1$
. The minimum strength for TPU is controlled by the simple Young’s modulus-optimal microstructures with first-order hierarchy (zeroth order indicates base material properties according to Lakes (Reference Lakes1993)) and shows a close to cubic mass density dependency (Haghpanah et al., Reference Haghpanah, Papadopoulos, Mousanezhad, Nayeb-Hashemi and Vaziri2014). The maximum strength line for TPU is controlled by the optimised microstructures with second-order hierarchy, and it has a polynomial order of
$n_0=2.3$
that is slightly lower than the predicted one, 2.5, assuming self-hierarchy (see Andersen et al., Reference Andersen, Wang and Sigmund2021), moreover with a higher coefficient
$c_0$
. As the relative yield strengths of the base materials decrease, the mass density dependency order reduces and converges to
$n_0=1$
for steel that is governed by the yield upper bound. The maximum strength line meets the yield strength upper bound at relatively high mass density for epoxy, PC, and PC-Nano, while the one for TPU is solely governed by the maximum strength line, and the one for steel solely governed by the yield bound.
Conclusions
A general topology optimisation formulation is proposed to tailor microstructure properties, including stiffness, yield, and buckling strength. The formulation is validated by designing square cell microstructures with maximal stiffness and/or strength. Numerical results for stiffness-optimal microstructures show that yield strength can be enhanced by a yield strength constraint without sacrificing the material stiffness and buckling strength. Compared to the stiffness-optimised microstructures, strength-optimised microstructures develop geometric hierarchy that enhances their buckling resistance. These strength-optimised microstructures have considerably higher strengths and fail simultaneously due to yield and buckling at a slight decrease in stiffness.
Microstructure sets with volume fractions of 20%, 10%, and 5% are systematically designed by maximising buckling strength with different yield strength bounds. The geometric evolution of the optimised designs show that a lower yield strength bound allows microstructures to develop geometrical hierarchies, leading to higher buckling strength with decreased stiffnesses. Moreover, the effective Young’s modulus and yield strength evolutions show a monotonic relation. The optimised microstructures are post-evaluated for base materials with different relative yield strengths ranging from steel to TPU. The overall performances show that base materials with higher relative yield strength allow the strength-optimal designs to form geometrical hierarchies, leading to higher optimised strength with slightly decreased Young’s modulus. Furthermore, as the design volume fraction decreases, the transition between yield- to buckling-dominated failures moves further towards the buckling-optimised designs. The overall performance indicates that buckling strength determines optimised strength at low-volume fractions and yield strength at high-volume fractions. The base material’s relative yield strength determines the critical volume fractions for pure buckling-dominated or pure yield-dominated failures. For a volume fraction in between, the strength-optimised microstructure is the optimised one with simultaneous yield and buckling failures. Meanwhile, the failure mechanism switches from yield- to buckling-dominated when the optimised microstructure configurations evolve from the buckling-strength to the stiffness/yield-strength optimised one. The non-normalised properties of the optimised microstructures show that base materials with higher Young’s modulus lead to both higher Young’s modulus and strength. The polynomial order of the maximum strength lines of the optimised microstructures on the mass density reduces as the base material relative yield strength decreases, reducing from 2.3 for TPU with buckling dominated to 1 for steel with yield dominated.
Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Funding statement
This study was funded by the Villum Fonden through the Villum Investigator Project InnoTop.
Competing interest
The authors declare no competing interests exist.
Authorship contribution
F.W. designed and conducted the study, analysed the data, and wrote and edited the article. O.S. acquired funding, designed the study, analysed the data, and edited the article.
Open access
