1 Introduction

Today’s state of the art within industrial applications requires the usage of efficient and high-performance materials, which are optimally designed in terms of production costs, material savings, fuel consumption or their mechanical behaviour. In many cases, a useful choice of materials of components contributes to the overall performance significantly. Materials can not be classified as homogeneous, but are composed of various ingredients. Often, considerations at the structural (macroscale) level are no longer sufficient and investigations on a lower length scale, the so-called material- or microscale, become necessary. Based on experimental data, practical experience and available methods for computer-aided simulation and analysis, the responsible design engineer has to decide which ingredients are suitable for a goal-oriented composition of materials. The evaluation of different designs can be performed based on so-called multiscale methods for structural analysis. These also include the FE\(^2\) method, which allows a combined analysis of the macroscopic structural behaviour and the microscopic material behaviour by a numerical homogenisation scheme applied to the microscopic material structure.

Overall, a broad range of multiscale optimisation problems are considered intensively in literature. Particularly, solution strategies based on topology optimisation methods are of major interest. An extensive introduction on topology optimisation methods in general is given in [13]. Several problem formulations for the choice of optimal thermoelastic properties with the target of maximum, zero or negative thermal expansion as well as of weight minimisation and manufacturing constraints can be found in [62,63,64,65]. The design of periodic linear elastic microstructures of cellular materials in terms of maximisation of weighted sum of equivalent strain energy density or linear combinations of mechanical properties under volume and material symmetry constraints are discussed in [51] as well as for minimisation of mean compliance subject to volume fractions of constituents in an underlying RVE in [22]. The key idea of extreme materials is often represented by applications with negative Poisson’s ratios. Several authors [36, 49, 59, 62, 63] tackle optimisation problems and provide results and material representations for this kind of structures, which in general are characterised by their ability for high energy absorption and fracture resistance. Two-phase composite materials with targeted properties such as piezoelectric properties in combination with thermoelectric coefficients can be found in [30]. General investigations for multiphase composites and their design process for optimal performance under certain constraints are proposed in [24, 70, 71] and especially in [69] for high stiffness or high strength composites. Procedures which tackle orthotropic materials, thickness of components and especially the optimisation of the angle of material rotation or material orientation in general are derived in [54, 55] and [56]. A computational model to design bi-material composite laminates to minimise structural compliance with mixed sets of micro and macro design variables is introduced in [19]. The target is to find optimal composite microstructures and optimal fibre orientation on the macroscale. Aspects on lightweight materials for non-periodic topologies within minimum compliant mechanism problems can be found in [37].

Strategies based on evolutionary optimisation algorithms, like the bi-directional evolutionary structural optimisation (BESO), for the concurrent design of material and structures within frameworks for topology optimisation and non-linear \(\hbox {FE}^{2}\) analysis with design variables on both scales are investigated in [82,83,84,85]. Aspects on efficient topology optimisation schemes for the design of multiscale non-linear heterogeneous structures and high performance materials are summarised in [81]. The authors in [31] bridge topology and shape optimisation schemes to design three-dimensional microstructured materials with extreme properties using energy-based homogenisation and parameteric level set methods. Patterns for engineered materials in terms of size, shape and layout of inclusion-like phases in continuum domains are obtained in [26].

Considering shape optimisation methods, theory and numerical methods for homogenisation techniques are formulated in [1]. Classical (fixed connectivity) shape optimisation has been outlined over the last decades in detail including the corresponding sensitivity analysis, see e.g. [16, 17, 29, 32]. In this manuscript, special emphasis is given to methods for (fixed-connectivity) shape optimisation applied to numerical homogenisation between two adjacent scales. Since the numerical effort for multiscale structural analysis raises with the complexity of referred problems, it is useful to focus on computational efficient methods for structural optimisation. Thus, so-called gradient based optimisation methods utilising the variational approach for sensitivity analysis are chosen. The focus is on relations for design (e.g. shape) modifications in terms of general geometrical parameters on both scales. But modifications like number, location and shapes of holes and inclusions are postponed. Similarly, characteristic material properties are not tackled here. The proposed method extends the formulation of a single-scale optimisation task and allows choices of objective functions or goals, constraints and design parameters on two scales.

In accordance to the well-known principal framework for structural optimisation pictured in Fig. 1, which is valid for single- and two-scale optimisation problems, the paper is organised as follows. Section 2 elaborates the advocated enhanced formulation of variational sensitivity analysis for shape optimisation problems on a single scale. Section 3 outlines a brief review of relations used for the structural analysis block, i.e. relations for numerical homogenisation and \(\hbox {FE}^{2}\) methods in terms of the Lagrange formalism are outlined. Section 4 provides required two-scale sensitivity information with respect to defined design parameters (DP), which is necessary for the evaluation of objective functions (OF), constraints (CON) and side-constraints (SCON) within the block design sensitivity analysis. The setup of two-scale optimisation problems for the block mathematical optimisation is presented in Sect. 5. Finally, Sect. 6 demonstrates the applicability of the derived formulation investigating an illustrative example.

Fig. 1
figure 1

Principle framework for structural optimisation

Full size image

The following notation is used. Second order tensors are depicted by Roman, upright, boldface letters, e.g. \(\mathbf{A }\). Vectors are printed in Roman, slanted, boldface letters, i.e. \(\varvec{X}\). Scalars are written, as usual, using Roman, slanted letters, i.e. \(\alpha \). All matrices are denoted by sans serif, slanted, boldface letters, i.e. \({{{\varvec{\mathsf{{ A}}}}}}, {{{\varvec{\mathsf{{ X}}}}}}\).

Furthermore, the notation of central quantities are mentioned at this early stage to enhance readability, especially in Sect. 2. Continuous mappings for geometry, motion, deformation and displacements are denoted by \({\varvec{\kappa }}, \varvec{\mu }, {\varvec{\varphi }}\) and \(\varvec{\upsilon }\). The corresponding vectors are written as \(\varvec{X}, \varvec{x}\) and \(\varvec{u}\). The (global) matrices in the finite element context read \({{{\varvec{\mathsf{{ \textsf {x}}}}}}}\) (nodal coordinates) and \({{{\varvec{\mathsf{{ u}}}}}}\) (nodal displacements), respectively. The dependency on time is indicated, as usual, by the parameter t. Similarly, the scalar variable s highlights the dependency on design.

Variations are written using the operator \(\delta \). A more compact representation of variations is chosen using the notation \( \delta {(\cdot )} = (\cdot )^\prime \). The total variation of a quantity of interest appears in the following compact way

$$\begin{aligned} \begin{aligned} \delta {(\cdot )}(\varvec{\upsilon }, \varvec{\kappa })&= (\cdot )^{\prime }_{\upsilon } + (\cdot )^{\prime }_{\kappa } = (\cdot )^\prime (\varvec{\upsilon }, \varvec{\kappa }). \end{aligned} \end{aligned}$$

Partial variations with respect to state and geometry functions can be identified by \(\delta _{\upsilon }{(\cdot )} = (\cdot )^{\prime }_{\upsilon }\) and \(\delta _{\kappa }{(\cdot )} = (\cdot )^{\prime }_{\kappa }\).

2 Variational sensitivity analysis on a single scale

The essential statements for answering the central question

“How do resulting quantities change

if input parameters are modified?”

can be deduced from sensitivity analysis, which plays the key role within frameworks for structural optimisation and is an important prerequisite for the success of mathematical optimisation algorithms. The design (e.g. shape) sensitivity analysis leads to insights about impacts and significant effects of design (e.g. shape) modifications on the physical properties and the resulting physical behaviour. Different approaches are well-known from literature including continuous variational, discrete (semi-) analytical as well as numerical methods, to name a few. A detailed discussion is beyond the scope of this paper and the readers are referred to the vast literature on structural optimisation, see e.g. [16, 17, 29, 34, 75].

2.1 Advantages of variational sensitivity analysis

The benefits of variational sensitivity analysis (VSA) compared to discrete approaches are emphasised in more detail. Within this class, especially the material derivative approach (MDA), see [2, 16, 17, 35], and the domain parametrization approach (DPA), see [15, 27, 57, 72], gained a high degree of respect based on numerous valuable theoretical and computational results for a wide range of challenging problems. Both approaches yield mutally equivalent results. Again, the reader is referred to literature for a detailed characterisation of the approaches. Importantly, all results are derived on the continuous level using variations of continuous functions. A subsequent discretisation yield the matrix equations needed for implementation and subsequent computation.

A central purpose of our contribution is to show, that a variational sensitivity analysis of complex problems such as the subject of this paper can be performed and presented in an understandable manor. We emphasise, that this statement is true for both mentioned variants as well as for the variational approach to (shape) sensitivity analysis outlined in Sect. 2.3.

2.2 Challenge of an enhanced variational sensitivity analysis

The advantages of an enhanced theoretical framework for variational (shape) sensitivity analysis are argued briefly. Most of the mentioned approaches use already worked out and widely recognised standards of continuum mechanics and numerical methods as starting point for performing sensitivity analysis. But historically, the development of continuum mechanics, computer aided geometric design, finite element method as well as structural and multidisciplinary optimisation neither took place simultaneously nor in a coordinated way. Similarly to shape optimisation, neighbouring research fields such as mechanics of growth and configurational mechanics matured after fundamental concepts of continuum mechanics had been consolidated. An integrated presentation of these fields is desirable. Thus, the concepts of variations of the body and its shape should be clarified within an enhanced formulation of continuum mechanics. Additionally, continuum mechanics and the other mentioned areas would benefit significantly from (i) rigorous mathematical concepts and (ii) their understandable translation to engineering language. Thus, a demanding and still unsolved optimisation problem for its own is the question on how to enhance both theory and numerics in a minimal form to achieve maximal insight, unifications and improvements for the above-mentioned research areas.

In a series of papers, Barthold and coworkers contributed to this problem, see e.g. [5, 7,8,9, 11] on the theory and numerics of variational sensitivity analysis, [38,39,40,41,42] on the interaction of variational sensitivity analysis and configurational mechanics, and [9] on the interaction of structural optimisation and growth. A thorough outline of these ideas is beyond the scope of this paper. Thus, only a brief sketch is given highlighting the central ideas.

2.3 Details of an enhanced variational sensitivity analysis

2.3.1 Remarks on an intrinsic, frame-free formulation

The material body \(\mathcal {B}\) described in [73], see [74, Sect. 15], is considered as given continuum consisting of material points \(\mathfrak {X}\). It is homogenized from the real matter composed of a large but finite integer number of elementary entities measured in the SI unit mole. Mathematically speaking, continuum mechanics is based on differentiable manifolds [74, Sect. 15]. It is important for understanding our approach to refer to a revised formulation by Noll [52], i.e. an intrinsic, frame-free viewpoint, which is completely formulated on the material body \(\mathcal {B}\) using intrinsic coordinates \(\varvec{\theta }\) of material points \(\mathfrak {X}\). Again, a detailed outline of Noll’s ideas is beyond the scope of this presentation and the readers are referred to e.g. [14, 52, 53] for details. The most central ingredients are highlighted.

The intrinsic geometry mapping \({\varvec{\kappa }}\) and the intrinsic motion mapping \(\varvec{\mu }\) define the reference placement \(\varvec{\mathcal {K}}\) and the current placement \(\varvec{\mathcal {M}}\), respectively. The dependency of motion on time t is obvious but hidden at most places. Thus, the referential deformation mapping \({\varvec{\varphi }}\) is composed in form of

$$\begin{aligned} {\varvec{\varphi }}= \varvec{\mu }\circ {\varvec{\kappa }}^{\mathrm {-1}}. \end{aligned}$$
(1)

Please observe the distinction between placement and configuration argued by Noll [52]. In the novel formulation, a configuration means the metric at the considered point, i.e. the central information which enters the definition of strains.

The referential deformation gradient \(\mathbf{F }\) split up into two tangent mappings with respect to intrinsic coordinates \(\varvec{\theta }\), i.e.

$$\begin{aligned} \mathbf{F }= \mathbf{M }\mathbf{K }^{{-\mathsf {1}}}, \end{aligned}$$
(2)

with the intrinsic motion gradient \(\mathbf{M }= {\nabla }_{\theta }\varvec{\mu }= {{\text {GRAD}}} \varvec{\mu }\) and the intrinsic geometry gradient \(\mathbf{K }= {\nabla }_{\theta }{\varvec{\kappa }}= {{\text {GRAD}}} {\varvec{\kappa }}\), see [52, Eq. 5.3], [53, Eq. 2.7] as well as [5, 7, 8]. The absolute tensor presentation can be transferred into a matrix notation. Let \( \varvec{X}= {\varvec{\kappa }}(\varvec{\theta }) \in \varvec{\mathcal {K}}\) and \(\varvec{x}= {\varvec{\mu }}(\varvec{\theta }) \in \varvec{\mathcal {M}}\) be the reference and current coordinates of a material point \(\mathfrak {X} \in \mathcal {B}\) with intrinsic coordinates \(\varvec{\theta }\in \varvec{\mathcal {P}}\). Let \(\varvec{E}_i =\varvec{e}_i\) and \(\varvec{Z}_i\) be Cartesian base vectors of reference and current placements as well as of the set \(\varvec{\mathcal {P}}\) of intrinsic coordinates, respectively. Thus,

(3)

are used to compute the deformation gradient

(4)

The (still continuous) matrix equation reads

(5)

Thus, the tensor \(\mathbf{K }\) results into the Jacobian matrix \({{{\varvec{\mathsf{{ K}}}}}}\) of the coordinate transformation between (Cartesian) reference and (Cartesian) intrinsic coordinates. These results are valid for any choice of intrinsic coordinates due to the underlying manifold structure.

2.3.2 Remarks on evolving material bodies

In contrast to all investigations for a single and fixed material body based on [52, 73, 74], growth and structural optimisation consider evolving material bodies. The modification of the material body is done either by nature over time t or by engineers within an optimisation cycle. Common to both challenges is the observation, that an enhanced continuum mechanical theory should incorporate a family of material bodies which are parametrised by \(s \in \mathbb {R}\), say a real valued design variable. In case of growth, the design s is a function of time t.

A suitable extension based on a thorough analysis of the manifold concept is outlined in [5, 7, 8]. In detail, any ball around any point of the manifold, i.e. the material body, can be mapped to a (local) coordinate system using invertible mappings \(\phi \). Consequently, the inverse mappings \(\phi ^{-1}\) defines the material body locally, i.e. for a small portion of the body. The choice of (local) coordinate systems is arbitrary and does not change the properties of the material body. And the change of coordinate systems is continuously differentiable. Without loss of generality, a design independent coordinate domain can be chosen. Altogether, the concept of a reservoir \(\varvec{\mathcal {R}}\) is deduced from the manifold property. Thus, material bodies are generated from the set \(\mathcal {R}\) by a function of local coordinates \(\theta \) and design s, see [9] for details. Here, the symbols of mappings and coordinates for either the intrinsic or the local description coincide to limit the notational complexity. Thus, the above-mentioned matrix equations are similarly valid for any local coordinate system used on the reservoir.

The outlined extension of Noll’s framework of mechanics on differentiable manifolds [52, 53] is a direct consequence of the manifold properties. Importantly, a concept to define a material body in a continuous differentiable way is available. Now, structural optimisation and growth can be investigated in an unified theoretical framework based on variations of functions [9]. To our best knowledge, no comparable formulation is available elaborating the evolution of material bodies.

2.3.3 Continuum mechanics as theory of two mappings

As a consequence, continuum mechanics can be interpreted as theory of two independent mappings, i.e. the (design s dependent) geometry mapping \({\varvec{\kappa }}\) (from reservoir \(\varvec{\mathcal {R}}\) to reference placement \(\varvec{\mathcal {K}}\)) and the (time t dependent) motion mapping \(\varvec{\mu }\) (from reservoir \(\varvec{\mathcal {R}}\) to current placement \(\varvec{\mathcal {M}}\)) [7]. This concept eases the investigations for sensitivity analysis and structural optimisation because all known functional dependencies are disassembled into parts, i.e. \(\, \alpha (s, t) = \beta (s) \circ \gamma (t)\). The known implicit dependency of motion \(\varvec{\mu }(t)\) on the (design s dependent) geometry \({\varvec{\kappa }}(s)\) is achieved once any kind of equilibrium condition is formulated.

Furthermore, all mathematical structures concerning either geometry (\({\varvec{\kappa }}, \mathbf{K }\)) or motion (\(\varvec{\mu }, \mathbf{M }\)) are very similar, i.e. any quantity in the physical world (parametrised by time t) has a twin in the material (configurational) world (parametrised by design s), see e.g. [42]. Again, to our best knowledge, no comparable conclusion has been published so far.

2.3.4 Fundamental variations of continuous functions

The (local) mappings for geometry \({\varvec{\kappa }}\) and motion \(\varvec{\mu }\) with

$$\begin{aligned} \varvec{X}= \varvec{X}(s) = {\varvec{\kappa }}(\varvec{\theta },s) \quad \text {and}\quad \varvec{x}= \varvec{x}(t) = \varvec{\mu }(\varvec{\theta }, t) \end{aligned}$$
(6)

are starting point of the advocated variational formulation of (shape) sensitivity analysis. Thus, the general expressions of perturbations (variations) of either geometry or motion read

$$\begin{aligned} \delta \varvec{X}= \delta \varvec{X}(s) = \delta {\varvec{\kappa }}(\varvec{\theta }, s) \quad \text {and}\quad \delta \varvec{x}= \delta \varvec{x}(t) = \delta \varvec{\mu }(\varvec{\theta }, t) . \end{aligned}$$
(7)

The (local) displacement mapping \(\varvec{\upsilon }= \varvec{\mu }- \varvec{\mu }_\circ \) is defined as difference between the motion mapping \(\varvec{\mu }\) at current time t with current placement \(\varvec{\mathcal {M}}\) and the motion mapping \(\varvec{\mu }_\circ \) at initial time \(t_\circ \) with initial placement \(\varvec{\mathcal {M}}_\circ \), i.e.

$$\begin{aligned} \varvec{u}= \varvec{u}(t) = \varvec{\upsilon }(\varvec{\theta },t) = \varvec{\mu }(\varvec{\theta },t) - \varvec{\mu }(\varvec{\theta },t_\circ ) = \varvec{x}- \varvec{x}_\circ . \end{aligned}$$
(8)

Usually, the initial placement \(\varvec{\mathcal {M}}_\circ \) is assumed to coincide with the reference placement \(\varvec{\mathcal {K}}\), i.e. \(\varvec{x}_\circ = \varvec{X}\). Thus, the relationship and its variation read

$$\begin{aligned} \varvec{x}= \varvec{X}+ \varvec{u}\quad \text {and}\quad \delta \varvec{x}= \delta \varvec{X}+ \delta \varvec{u}. \end{aligned}$$
(9)

As the consequence, continuum mechanics can be presented as theory of two local mappings, i.e. using geometry \({\varvec{\kappa }}\) and either motion \(\varvec{\mu }\) or displacement \(\varvec{\upsilon }\). The dependencies on time t and design s have been argued but are hidden at most places. The linkage of mappings is outlined in Sect. 2.3.5.

Next, the variation of tangent mappings is considered, see [5, 7,8,9] for details. Based on the variation of local geometry and motion mappings, i.e.

$$\begin{aligned} \delta \mathbf{K }= {{\text {GRAD}}} \delta {\varvec{\kappa }}\quad \text {and}\quad \delta \mathbf{M }= {{\text {GRAD}}} \delta \varvec{\mu }, \end{aligned}$$
(10)

the variation of the (referential) deformation gradient \(\mathbf{F }\) reads

$$\begin{aligned} \begin{aligned} \delta \mathbf{F }= \delta {{\text {Grad}}} \varvec{x}&= {{\text {Grad}}}\, \delta \varvec{x}- {{\text {Grad}}} \varvec{x}\, {{\text {Grad}}} \delta \varvec{X}\\&= {{\text {Grad}}} \delta \varvec{u}- {{\text {Grad}}} \varvec{u}\, {{\text {Grad}}}\, \delta \varvec{X}. \end{aligned} \end{aligned}$$
(11)

Three different gradient operators \({{\text {GRAD}}}, {{\text {Grad}}}\) and \({{\text {grad}}}\) are used, which are defined on the intrinsic, referential and current domains, respectively.

These preparations are the foundation of an enhanced formulation of variational sensitivity analysis. In detail, a consistent linearisation technique for nonlinear mechanics, cf. [80], is utilised. A detailed report of all findings is beyond the scope of this paper, see [5, 7, 8]. Central results needed in this publication are mentioned next.

2.3.5 Variational sensitivity of continuous functions

The linkage between the (independent) geometry and motion mappings is given by an equilibrium condition. On the continuous level, strong and weak forms of equilibrium are equivalent if sufficient smoothness is assumed. With regard to the discrete formulations, only the weak form

$$\begin{aligned} R (\varvec{\upsilon }, \varvec{\kappa }; \varvec{\eta }) = 0 \end{aligned}$$
(12)

is considered at this place, where \( \varvec{\eta }\) denotes the test function. Thus, the motion and the displacements depend on the chosen geometry. The introduced notation using commas and semicolons allows to distinguish between linear and non-linear arguments in functionals and tangent forms, when needed. For instance, a functional \((\cdot )(\cdot ; \cdot , \cdot )\) is non-linear in arguments on the left hand side of the semicolon, and linear in arguments on the right hand side of the semicolon, respectively.

Equilibrium should never be violated. Therefore, the total variation of \(R = 0\) must vanish, i.e.

$$\begin{aligned} \delta {R} = \delta _{\upsilon }{R} + \delta _{\kappa }{R} = k (\varvec{\eta }, \delta \varvec{\upsilon }) + p (\varvec{\eta }, \delta \varvec{\kappa }) = 0 . \end{aligned}$$
(13)

Two tangent operators occur, i.e. a tangent stiffness operator k and a tangent pseudo load operator p. The variation of the (continuous) geometry mapping \(\delta \varvec{X}= \delta {\varvec{\kappa }}(\varvec{\theta }, s)\), is answered by a corresponding variation of the (continuous) state function, i.e. of the displacements \(\delta \varvec{u}= \delta \varvec{\upsilon }(\varvec{\theta }, t)\). From Eq. (13), the (continuous) sensitivity operator \(\mathfrak {s}\) is implicitely defined

$$\begin{aligned} \delta \varvec{\upsilon }= \mathfrak {s} (\varvec{\upsilon }, \varvec{\kappa }; \delta \varvec{\kappa }) . \end{aligned}$$
(14)

The difference of the scalar variable s, introduced to model the functional dependencies on design, and the sensitivity operator \(\mathfrak {s}\) is highlighted by using different letter fonts.

The variation of any continuum mechanical function f can be derived based on the already gathered information, i.e.

$$\begin{aligned} \delta f = \delta _{\upsilon }{f} + \delta _{\kappa }{f} = f^{\prime }_{\upsilon } + f^{\prime }_{\kappa } = a (\delta \varvec{\upsilon }) + b (\delta \varvec{\kappa }) , \end{aligned}$$
(15)

with the linear operators \(a = \delta _{\upsilon }{f} = f^{\prime }_{\upsilon }\) and \(b = \delta _{\kappa }{f} = f^{\prime }_{\kappa }\). All operators, evaluated at \((\varvec{\upsilon },\varvec{\kappa })\), are elaborated in detail for the two-scale optimisation problem in Sect. 4.

2.3.6 General remarks on discrete approximations

Some remarks on discrete approximations of continuous mappings and continuous equations are added.

A differentiable manifold, i.e. the atlas \(\mathcal {A}\) of charts with sufficiently smooth transformation mappings between coordinate systems, is the continuous master concept to computational techniques. The choice of a specific atlas does not change the manifold, i.e. the derived continuous equations are valid for any choice of atlas \(\mathcal {A}\).

Any discrete technique using some approximation spaces of shape and Ansatz functions yield discrete atlases, e.g. \(\mathcal {A}_{\mathrm {cagd}}\) (computer aided geometric design) and \(\mathcal {A}_{\mathrm {fem}}\) (finite element method) as outlined in [7, Sect. 2.3]. Importantly, the quality of the approximations depend on the chosen discrete atlases.

The following three function spaces and their approximations are mentioned without discussing their mathematical details. All discrete (finite dimensional) function spaces and all elements of theses spaces are denoted by a subscript h. First, the continuous geometry mapping \({\varvec{\kappa }}\in \mathcal {G}\) is approximated by \({\varvec{\kappa }}_h \in \mathcal {G}_h \subset \mathcal {G}\). Second, the continuous displacement mapping \(\varvec{\upsilon }\in \mathcal {V}\) is approximated by \(\varvec{\upsilon }_h \in \mathcal {V}_h \subset \mathcal {V}\). Furthermore, the specific optimization problem introduces a continuous design space \(\mathcal {S}\subset \mathcal {G}\) in order to fulfill certain constraints. Thus, third, the variation of geometry \(\delta {\varvec{\kappa }}\in \mathcal {S}\subset \mathcal {G}\) is approximated by \(\delta {\varvec{\kappa }}_h \in \mathcal {S}_h \subset \mathcal {G}_h\).

2.3.7 Discretisation of continuous geometry mappings

The combination of CAGD and FEM is the most widely used approximation method for geometry. The continuous mapping \({\varvec{\kappa }}\) is approximated using a CAGD model, say a Bézier surface, with \({\varvec{\mathsf{{ y}}}} \in \mathbb {R}^{n_\textsf {y}}\) parameters, say the Bézier control point coordinates, see e.g. [21] for details. Furthermore, the nodal coordinates \({\varvec{\mathsf{{ x}}}} \in \mathbb {R}^{n_\textsf {x}}\) of any subsequent FEM model are linked to the CAGD parameters \({{{\varvec{\mathsf{{ \textsf {y}}}}}}}\), i.e. \({\varvec{\mathsf{{ x}}}} = {{{\varvec{\mathsf{{ g}}}}}} ({\varvec{\mathsf{{ y}}}})\) for example using a mapped meshing technique. Last but not least, the finite element technique adds another approximation on the finite element domain, see Sect. 2.3.8 for further remarks. In all cases, the discrete approximations are member of finite dimensional spaces with dimensions \(n_\textsf {x}\) and \( n_\textsf {y}\), respectively.

The linkage between CAGD and FEM models is well-known in literature on shape optimisation since decades and has also been investigated in [6] describing the effect of mapped meshing on sensitivity analysis. In [4], the mesh refinement error for optimisation problems has been investigated. In [20] hierarchical geometry models are investigated concerning their approximation properties.

The variation of the geometry mapping \(\delta {\varvec{\kappa }}\) substitutes the design velocity field introduced in the material derivative approach (MDA), see e.g. [87] and [16, 17]. Any detailed computation of \(\delta {\varvec{\kappa }}\) needs a precise knowledge of the mapping \({\varvec{\kappa }}\), which is only available in the context of specific computational methods.

Here, the approximation \(\delta {\varvec{\kappa }}_h \in \mathcal {S}_h \subset \mathcal {G}_h\) is exemplified in the context of CAGD and FEM coupling. In (shape) optimisation, the CAGD parameters \({\varvec{\mathsf{{ y}}}}\) (coordinates of \(\mathcal {G}_h\)) are linked to some finite number of design variables \({{{\varvec{\mathsf{{ s}}}}}}\in \mathbb {R}^{n_s}\) (coordinates of \(\mathcal {S}_h\)) via \({\varvec{\mathsf{{ y}}}} = {{{\varvec{\mathsf{{ h}}}}}} ({{{\varvec{\mathsf{{ s}}}}}})\). Thus, the variation \(\delta {\varvec{\kappa }}\) is approximated by \(\delta {\varvec{\kappa }}_h\) using the matrices

$$\begin{aligned} {{{\varvec{\mathsf{{ G}}}}}} = \frac{\text {d}{\varvec{\mathsf{{ x}}}}}{\text {d}{\varvec{\mathsf{{ y}}}}} \in \mathbb {R}^{n_\textsf {x}\times \, n_\textsf {y}} \quad \text{ and }\quad {{{\varvec{\mathsf{{ H}}}}}} = \frac{\text {d}{\varvec{\mathsf{{ y}}}}}{\text {d}{{{\varvec{\mathsf{{ s}}}}}}} \in \mathbb {R}^{n_\textsf {y}\times \, n_s} \, . \end{aligned}$$
(16)

The combination yields the mesh design matrix

$$\begin{aligned} {{{{\varvec{\mathsf{{ V}}}}}}}= \frac{\text {d}{\varvec{\mathsf{{ x}}}}}{\text {d}{{{\varvec{\mathsf{{ s}}}}}}} = \frac{\text {d}{\varvec{\mathsf{{ x}}}}}{\text {d}{\varvec{\mathsf{{ y}}}}} \, \frac{\text {d}{\varvec{\mathsf{{ y}}}}}{\text {d}{{{\varvec{\mathsf{{ s}}}}}}} = {{{\varvec{\mathsf{{ G}}}}}}\, {{{\varvec{\mathsf{{ H}}}}}} \in \mathbb {R}^{n_\textsf {x}\times \, n_s} \, , \end{aligned}$$
(17)

which is frequently used in computations, see Sect. 6.

2.3.8 Approximations on the finite element domain

The continuous displacement mapping \(\varvec{\upsilon }\) is approximated by FEM using e.g. nodal displacement parameters \({{{\varvec{\mathsf{{ u}}}}}}\in \mathbb {R}^{n_u}\). At this place, the isoparametric finite element technique using a fixed parent domain and n element nodes is exemplified, i.e.

$$\begin{aligned} u^i = \sum _{I=1}^{n} h_I (\xi _j)\, u^i_I\quad \text {and} \quad X^i = \sum _{I=1}^{n} h_I (\xi _j)\, X^i_I \, . \end{aligned}$$
(18)

The method combines shape functions \( h_I \) of coordinates \( \xi _j \) on the parent domain and nodal coordinates \(X^i_I\) and nodal displacements \(u^i_I\). The technique to compute a Cartesian derivative for shape functions is using the Jacobian matrix \({{{\varvec{\mathsf{{ J}}}}}}_e\) of the coordinate transformation between referential and local coordinates of the finite element parent domain, see e.g. [80].

Concerning sensitivity analysis, the continuous variations \(\delta \varvec{X}= \delta {\varvec{\kappa }}(\varvec{\theta },s)\) and \(\delta \varvec{u}= \delta \varvec{\upsilon }(\varvec{\theta },s)\) are linked to perturbations of nodal coordinates \(\delta X^i_I\) and nodal displacements \(\delta u^i_I\).

Next, the computation (on the finite element level) of the Cartesian derivatives of shape functions, the displacement approximation as well as of the corresponding variations can be directly deduced from Eqs. (5) to (11), i.e. by evaluating continuum mechanical results. As already mentioned above, the tangent mapping \(\mathbf{K }\) (tensor of 2nd order) is the origin of the Jacobian matrix \({{{\varvec{\mathsf{{ J}}}}}}_e\). Furthermore, the derivatives of shape functions \(h_I\) with respect to (Cartesian) referential coordinates \(X^i\) are a direct result of combining Eqs. (5) and (18), i.e.

(19)

Similarly, analytical derivatives of finite element quantities with respect to nodal coordinates can be linked to the continuous variations such as Eq. (11). These remarks highlight the fundamental concepts of the advocated approach to sensitivity analysis based on Noll’s intrinsic, frame-free formulation of continuum mechanics and the subsequent variations with respect to both fundamental mappings. The final results are similar to those obtained by DPA and MDA. A detailed comparison is beyond the scope of this paper and the readers are referred to the cited literature.

2.3.9 Nonlinear problems solved by the advocated method

The advocated approach has been successfully applied to shape sensitivity analysis of nonlinear behaviour such as hyperelasticity with isotropic damage [10], multiplicative elasto-plasticity [76], shake down analysis [77] and the theory of porous media [67], to name a few early applications.

3 Structural analysis on two scales: computational homogenisation and the \(\hbox {FE}^{2}\) method

Heterogeneities are part of most natural and non-natural materials and often cause anisotropic effects in stress and strain fields. The spatial distribution of individual and so-called microscale constituents on the so-called mesoscale, the specific properties of each constituent as well as their size and shape have a strong influence on the overall mechanical behaviour on the macroscopic structural level. To be able to postulate reliable predictions about the overall behaviour or possible failure scenarios, fundamental knowledge about the physical behaviour on lower scales is required. Homogenisation methods and so-called \(\hbox {FE}^{2}\) techniques act like a bridging approach for the coupling of structural and material scales and allow to include mentioned effects into overall investigations. A basic introduction on numerical two-scale homogenisation methods is given in [86] or [60]. A discussion on computational homogenisation, a general overview and a summary of developments is given in [68] and [43]. A variety of methods is presented in [50], and especially the article [23] reviews the state-of-the-art of computational homogenisation and discusses trends and upcoming challenges in this field.

3.1 Continuous formulation of homogenisation

In this section, some of the aforementioned aspects on mechanical two-scale problems are summarised. The notation of following descriptions is related to the formulations in [12, 45, 46, 48]. Throughout this work, all frequently used quantities connected to the macroscale are identified by overlines, i.e. \(\overline{(\cdot )}\), and quantities connected to the microscale are represented without additional markers, i.e. \((\cdot )\).

Fig. 2
figure 2

Basic homogenisation scheme

Full size image

The basic concept of the homogenisation scheme formulated by using continuous quantities is pictured in Fig. 2. The underlying representative volume element (RVE) with the domain \(\varvec{\mathcal {K}}\) is associated with each macroscopic point \(\overline{\varvec{X}}\) in the macroscopic domain \(\overline{\varvec{\mathcal {K}}}\). The macroscopic deformation gradient \(\overline{\mathbf{F }}\) in \(\overline{\varvec{X}}\) is the necessary driving factor for the formulation of boundary conditions for the microscopic boundary value problem (BVP). The solution investigated within homogenisation approaches results in effective stress and material parameters, i.e. \(\overline{\mathbf{P }}\) and \(\overline{\mathbb {A}}\). The effective response \(\overline{\mathbf{P }}\) and \(\overline{\mathbb {A}}\) of the microscale influences the solution of the macroscopic BVP.

The characterisation of the aforementioned micro to macro transition as an energy minimisation problem for microstructures with constituents of standard materials is presented in [44, 45, 47]. For each \(\overline{\varvec{X}} \in \overline{\varvec{\mathcal {K}}}\), an exact geometrical resolution of the referred microscale is required and a non-linear elastic material behaviour is considered, see Eq. (20),

(20)

The energy depends on the microscopic deformation \({\varvec{\varphi }}\), which itself is induced by the macroscopic deformation gradient \(\overline{\mathbf{F }}\). The average energy of the microstructure is minimised while boundary conditions on the microscale are enforced, i.e. either linear (D) or periodic (P) displacements or uniform tractions (S), see Sect. 3.2. The boundary conditions \(\{\text {D},\ \text {P},\ \text {S}\}\) can be applied using the Lagrange formalism, which allows to recast Eq. (20) into the following saddle point problem

$$\begin{aligned} \overline{W}_{\text {I}}^{\lambda } (\overline{\mathbf{F }}) \! = \inf _{{\varvec{\varphi }}} \sup _{\lambda _\text {I}} \left\{ \frac{1}{V}\int _{\varvec{\mathcal {K}}} W({\nabla }{\varvec{\varphi }};\varvec{X})\; \text {d}V - c_\text {I}({\varvec{\varphi }}, \varvec{\lambda }_\text {I}; \overline{\mathbf{F }}) \right\} . \end{aligned}$$
(21)

Here, a general form of enforced constraints \(c_\text {I}\) and a general description of the necessary Lagrange multiplier \(\varvec{\lambda }_\text {I}\) for introduced classes of boundary conditions \(\text {I}= \{\text {D},\ \text {P},\ \text {S}\}\) is used. The effective stress and material parameters can be obtained from Eq. (21) directly by the evaluation of

$$\begin{aligned} \overline{\mathbf{P }}_\text {I}= \partial _{\overline{\mathbf{F }}}\ \overline{W}_{\text {I}}^{\lambda } \qquad \text {and} \qquad \overline{\mathbb {A}}_\text {I}= \partial _{\overline{\mathbf{F }}}\ \overline{\mathbf{P }}_\text {I}= \partial ^2_{\overline{\mathbf{F }} \overline{\mathbf{F }}}\ \overline{W}_{\text {I}}^{\lambda }. \end{aligned}$$
(22)

3.2 General form of boundary conditions at discrete points

The abstract formulation of general boundary conditions is omitted at this place, because the boundary conditions in the framework of the FE\(^2\) method are evaluated at discrete nodal points of the finite element mesh. Nevertheless, all quantities are still continuous mappings but evaluated at discrete points. The very first step is the partitioning of nodes of the underlying FE mesh of the discretised microscopic domain \(\varvec{\mathcal {K}}\). Two sets of nodes can be distinguished: the nodes in the interior of the domain \(\varvec{\mathcal {K}}\), i.e. nodes indicated by the index \((\text {i})\), and nodes on the boundary \(\partial \varvec{\mathcal {K}}\) of the RVE, i.e. nodes indicated by the index \((\text {b})\), see [12, 46] or [33] for further explanations.

The structure of introduced boundary conditions \(\{\text {D},\ \text {P},\ \text {S}\}\) is known from standard literature on homogenisation. The obvious similarity allows to find a general formulation for the representation of boundary conditions and to use this unique and compact notation for further purposes and investigations. Using the Voigt notation for the deformation gradient \(\overline{{{\underline{{{\varvec{\mathsf{{ F}}}}}}}}}\), this general form is defined by

$$\begin{aligned} \mathcal {S}_1\, {{{\varvec{\mathsf{{ u}}}}}}_\text {b}- \mathcal {S}_2\, (\overline{{{\underline{{{\varvec{\mathsf{{ F}}}}}}}}} - {\underline{{{\varvec{\mathsf{{ I}}}}}}}) = \mathbf{0} \end{aligned}$$
(23)

for introduced classes of boundary conditions \(\text {I}= \{\text {D},\ \text {P},\ \text {S}\}\). For instance, in the two-dimensional case the deformation gradient reads . The corresponding coefficient matrices \(\mathcal {S}_1\) and \(\mathcal {S}_2\) are of the following forms

$$\begin{aligned} \mathcal {S}_1:= {\left\{ \begin{array}{ll} \mathcal {I}\in \mathbb {R}^{n_{\mathrm {B}}\times n_{\mathrm {B}}} &{} (\text {D}) \\ \mathcal {P}\in \mathbb {R}^{n_{\mathrm {P}}\times n_{\mathrm {B}}} &{} (\text {P}) \\ \mathcal {A}\in \mathbb {R}^{d^2 \times n_{\mathrm {B}}} &{} (\text {S}) \end{array}\right. } \end{aligned}$$
(24)

and

$$\begin{aligned} \mathcal {S}_2:= \left\{ \begin{aligned} \mathcal {D}^T&\in \ \mathbb {R}^{n_{\mathrm {B}}\times d^2} (\text {D}) \\ {\mathcal {P}}\, \mathcal {D}^T&\in \ \mathbb {R}^{n_{\mathrm {P}}\times d^2} (\text {P}) \\ V\mathcal {I}\,&\in \ \mathbb {R}^{d^2 \times d^2} (\text {S}) \end{aligned}\right. \end{aligned}$$
(25)

with \(\mathcal {I}\) being an identity matrix, \(V\) the volume of the RVE, \(n_{\mathrm {B}}\) number of nodes on the boundary, \(n_{\mathrm {P}}\) number of periodicity conditions and \(d\) the dimension of the stated problem. For better clarity in terms of comprehensivness of the general form of boundary conditions, the matrices \(\mathcal {D}, \mathcal {P}\) and \(\mathcal {A}\) are depicted shortly. The matrix \(\mathcal {D}\in \mathbb {R}^{d^2 \times n_{\mathrm {B}}}\) is a global coordinate matrix

$$\begin{aligned} \mathcal {D}= \begin{bmatrix} \mathcal {D}_1&\mathcal {D}_2&\ldots&\mathcal {D}_{n_{\mathrm {B}}} \end{bmatrix} \end{aligned}$$
(26)

and contains coordinates of each node \({{{\varvec{\mathsf{{ \textsf {x}}}}}}}_q\) with \(q= 1,\ldots ,n_{\mathrm {B}}\) on the surface of the underlying RVE in terms of matrices \(\mathcal {D}_q\in \mathbb {R}^{d^2 \times d}\). Each matrix \(\mathcal {D}_q\) has the following explicit form

(27)

The periodicity matrix \(\mathcal {P}\in \mathbb {R}^{n_{\mathrm {P}}\times n_{\mathrm {B}}}\) is a topology matrix including only values \(\{0,1,-1\}\). The entries of the matrix \(\mathcal {P}\) can be constructed using following cases

$$\begin{aligned} \mathcal {P}_{ij} = \left\{ \begin{aligned} 1 \quad&\text {for } i= j, \\ -1 \quad&\text {for coupled nodes on surface } (i\ne j), \\ 0 \quad&\text {otherwise}. \end{aligned}\right. \end{aligned}$$
(28)

The matrix \(\mathcal {A}\in \mathbb {R}^{d^2 \times n_{\mathrm {B}}}\) contains components \(\mathcal {A}_q\in \mathbb {R}^{d^2 \times d}\)

$$\begin{aligned} \mathcal {A}= \begin{bmatrix} \mathcal {A}_1&\mathcal {A}_2&\ldots&\mathcal {A}_{n_{\mathrm {B}}} \end{bmatrix} \quad \text {and} \quad \mathcal {A}_q:= \begin{bmatrix} {{{\textsf {\textit{A}}}}}_1 &{} 0 \\ 0 &{} {{{\textsf {\textit{A}}}}}_2 \\ {{{\textsf {\textit{A}}}}}_2 &{} 0 \\ 0 &{} {{{\textsf {\textit{A}}}}}_1 \end{bmatrix}_q, \end{aligned}$$
(29)

where \( ({{{\textsf {\textit{A}}}}}_1,{{{\textsf {\textit{A}}}}}_2)\) in \(\mathcal {A}_{q}\) represent entries of the discrete normal area vector \({{{\varvec{\mathsf{{ A}}}}}}_q\) at node \(q\) in the referential placement. Further derivations, extended explanations and details in terms of the chosen notation can be found in [12, 46].

Fig. 3
figure 3

Enhanced kinematics based on two independent geometry and motion mappings on two scales with the following ingredients. (i) Point spaces: referential placements \(\overline{\varvec{\mathcal {K}}}\) and \(\varvec{\mathcal {K}}\), current placements \(\overline{\varvec{\mathcal {M}}}\) and \(\varvec{\mathcal {M}}\), and reservoirs \(\overline{\varvec{\mathcal {R}}}\) and \(\varvec{\mathcal {R}}\). (ii) Point mappings: deformation mappings \(\overline{{\varvec{\varphi }}}\) and \({\varvec{\varphi }}\), geometry mappings \(\overline{{\varvec{\kappa }}}\) and \({\varvec{\kappa }}\), and motion mappings \(\overline{{\varvec{\mu }}}\) and \({\varvec{\mu }}\). (iii) Tangent spaces not depicted (iv) Tangent mappings: deformation gradients \(\overline{\mathbf{F }}\) and \(\mathbf{F }\), geometry gradients \(\overline{\mathbf{K }}\) and \(\mathbf{K }\), and motion gradients \(\overline{\mathbf{M }}\) and \(\mathbf{M }\). (v) Decompositions: point mappings see Eq. (32) tangent mappings see Eq. (33)

Full size image

3.3 Discrete minimisation problem of homogenisation

Again, the continuous Lagrange functional is omitted to reduce the complexity of the presentation. Based on the specification of boundary conditions in Eq. (23) in their discrete form, the discrete Lagrange functional for the solution procedure of the discrete minimisation problem of homogenisation can also be stated in the general form

$$\begin{aligned} \overline{W}_{\text {I}}^{\lambda } (\overline{\mathbf{F }}) = \inf _{{{{\varvec{\mathsf{{ u}}}}}}} \sup _{{{\varvec{\uplambda }}}_{\text {I}}} \left\{ \tilde{W}({{{\varvec{\mathsf{{ u}}}}}}) - {{\varvec{\uplambda }}}^{T}_{\text {I}} \left[ \mathcal {S}_1\, {{{\varvec{\mathsf{{ u}}}}}}_{\text {b}} - \mathcal {S}_2\, (\overline{{{\underline{{{\varvec{\mathsf{{ F}}}}}}}}} - {\underline{{{\varvec{\mathsf{{ I}}}}}}}) \right] \right\} . \end{aligned}$$
(30)

To solve this problem in terms of minimisation, variations with respect to partitioned state functions \({{{\varvec{\mathsf{{ u}}}}}}_\text {i}\) and \({{{\varvec{\mathsf{{ u}}}}}}_\text {b}\) (evaluated at discrete points) and with respect to the (finite number of) Lagrange multipliers \({{\varvec{\uplambda }}}_\text {I}\) are necessary.

Based on the solution of Eq. (30), homogenised stresses and the overall tangent moduli can be obtained by the evaluation of Eq. (22) and, in the discrete sense, they result to

$$\begin{aligned} \begin{aligned} \overline{{{\underline{{{\varvec{\mathsf{{ P}}}}}}}}}_{\text {I}}&:= \partial _{\overline{\mathbf{F }}}\ \overline{W}_{\text {I}}^{\lambda } = \mathcal {S}_2^T {{\varvec{\uplambda }}}_\text {I}\quad \\ \overline{{\varvec{\mathsf{{ A}}}}}_{\text {I}}&:= \partial ^2_{\overline{\mathbf{F }} \overline{\mathbf{F }}} \overline{W}_{\text {I}}^{\lambda } = \partial _{\overline{\mathbf{F }}}\ \overline{{{\underline{{{\varvec{\mathsf{{ P}}}}}}}}}_{\text {I}} = \mathcal {S}_2^T \frac{\partial {{\varvec{\uplambda }}}_\text {I}}{\partial \overline{{{\underline{{{\varvec{\mathsf{{ F}}}}}}}}}}= \mathcal {S}_2^T \overline{{{{{\varvec{\mathsf{{ K}}}}}}}}_\text {I}\mathcal {S}_2. \end{aligned} \end{aligned}$$
(31)

It can be observed that effective stresses depend only on the Lagrange multiplier \({{\varvec{\uplambda }}}_\text {I}\), i.e. on terms defined on the surface of the RVE. The homogenised tangent moduli can be evaluated in terms of the condensed matrix \(\overline{{{{{\varvec{\mathsf{{ K}}}}}}}}_\text {I}= \frac{1}{V}(\mathcal {S}_1\tilde{{{{{\varvec{\mathsf{{ K}}}}}}}}_{\text {b}\text {b}}^{-1} \mathcal {S}_1^T)^{-1}\) with \(\tilde{{{{{\varvec{\mathsf{{ K}}}}}}}}_{\text {b}\text {b}} = {{{{\varvec{\mathsf{{ K}}}}}}}_{\text {b}\text {b}} - {{{{\varvec{\mathsf{{ K}}}}}}}_{\text {b}\text {b}} {{{{\varvec{\mathsf{{ K}}}}}}}_{\text {i}\text {i}}^{-1} {{{{\varvec{\mathsf{{ K}}}}}}}_{\text {i}\text {b}}\), cf. [46] and [12].

4 Two-scale structural design based on variational two-scale sensitivity analysis

Understanding two-scale homogenisation techniques and \(\hbox {FE}^{2}\) methods gives the opportunity to analyse complex and heterogeneous materials on different length scales. This outstanding and deep insight into the physical behaviour of individual constituents of heterogeneous material compositions comes along with the allure of improvements in terms of the overall performance or adjustments of the physical response. The presented framework for two-scale analysis is predestinated for purposes within structural optimisation, especially due to its variety of possible applications. The major target is to obtain a closed formulation, which contains the sensitivity information of the overall macroscopic BVP and therefore, the sensitivity information of all underlying microscopic BVPs.

4.1 Basic sensitivity relations on two scales

The layout for sensitivity analysis on single scales, see Sect. 2, is extended for two-scale problems, see Fig. 3 for an overview of mappings. In detail, the macroscopic deformation mapping \(\overline{{\varvec{\varphi }}}(\overline{\varvec{X}})\), depending on the macroscopic geometry \(\overline{{\varvec{\kappa }}}\), and the microscopic deformation mapping \({\varvec{\varphi }}(\varvec{X})\), depending on the microscopic geometry \({\varvec{\kappa }}\), are introduced. Both mappings can be equivalently decomposed in so-called local geometry and local deformation mappings with the distinction in macroscopic and microscopic contributions, i.e. using mappings \(\{\overline{{\varvec{\mu }}}, \overline{{\varvec{\kappa }}}\}\) for the macroscale and mappings \(\{{\varvec{\mu }}, {\varvec{\kappa }}\}\) for the microscale, see Eq. (1) for single scales and now Eq. (32). Thus, one obtains the decompositions

$$\begin{aligned} \overline{{\varvec{\varphi }}} = \overline{{\varvec{\mu }}} \circ \overline{{\varvec{\kappa }}}^{-1} \qquad \text {and} \qquad {\varvec{\varphi }}= {\varvec{\mu }}\circ \varvec{\kappa }^{-1}. \end{aligned}$$
(32)

Similarly and referring to Eq. (2), the decomposition of both deformations gradients into two tangent mappings read

$$\begin{aligned} \overline{\mathbf{F }} = \overline{\mathbf{M }} \, \overline{\mathbf{K }}^{-1} \qquad \text {and} \qquad \mathbf{F }= \mathbf{M }\, \mathbf{K }^{-1}. \end{aligned}$$
(33)

The full set of point mappings and affiliated tangent mappings in terms of macro- and microscopic notation is summarised in Table 1.

Table 1 Overview of macro- and microscopic point and tangent mappings
Full size table

All introduced tangent mappings can be used to perform pull-back and push-forward operations between introduced tangent spaces, i.e. those of the placements \(\overline{\varvec{\mathcal {K}}}, \overline{\varvec{\mathcal {M}}}\) and the reservoir \(\overline{\varvec{\mathcal {R}}}\) on the macroscale and between the tangent spaces of the placements \(\varvec{\mathcal {K}}, \varvec{\mathcal {M}}\) and the reservoir \(\varvec{\mathcal {R}}\) on the microscale for several variations of quantities from continuum mechanics. The advantage of the approach is, that the methodological procedure for obtaining analytical variations of physical quantities can be proceeded independent of the scale under investigation.

4.2 Variation of the macroscopic weak form of equilibrium

The full list of dependencies for the macroscopic weak form of equilibrium for all upcoming investigations within design sensitivity analysis, and later structural optimisation, reads

$$\begin{aligned} \overline{R} (\overline{\varvec{\upsilon }}, \overline{\varvec{\kappa }}, \varvec{\upsilon }, \varvec{\kappa }; \overline{\varvec{\eta }}) = 0. \end{aligned}$$
(34)

Considering these dependencies, the total variation of the macroscopic physical residual reads

(35)

Here, the variations of the macroscopic physical residual with respect to state functions \(\overline{\varvec{\upsilon }}\) and \(\varvec{\upsilon }\) are introduced by tangent stiffness operators \(\overline{k}\) and and variations with respect to geometry functions \(\overline{\varvec{\kappa }}\) and \(\varvec{\kappa }\) are introduced by tangent pseudo operators \(\overline{p}\) and , respectively. In the following, each tangent operator is stated

(36)

Throughout the entire work, the term “multilevel” indicates variations of the macroscopic residual \(\overline{R}\) with respect to functions related to the microscale, i.e. quantities like the microscopic state \(\varvec{\upsilon }\) and the microscopic geometry \(\varvec{\kappa }\).

4.3 Explicit weak form of equilibrium and its variation

The explicit form of the macroscopic residual in terms of the macroscopic deformation gradient \(\overline{\mathbf{F }}\) and now, in terms of effective or homogenised stresses \(\overline{\mathbf{P }}_\text {I}= \displaystyle \frac{1}{V}\int _{\varvec{\mathcal {K}}} \mathbf{P }\; \text {d}V\) according to relations for numerical homogenisation techniques presented in Sect. 3, reads

$$\begin{aligned} \overline{R}(\overline{\varvec{\upsilon }},\overline{\varvec{\kappa }}, \varvec{\upsilon }, \varvec{\kappa }; \overline{\varvec{\eta }}) = \int _{\overline{\varvec{\mathcal {K}}}} \overline{\mathbf{P }}_\text {I}:\overline{\mathbf{F }}^{\prime }_{\overline{\upsilon }}(\overline{\varvec{\upsilon }},\overline{\varvec{\eta }})\; \text {d}V - \overline{F}(\overline{\varvec{\kappa }},\overline{\varvec{\eta }}). \end{aligned}$$
(37)

The total variation of the residual form in Eq. (35) contains partial variations with respect to all introduced parameters in the argument list, i.e. with respect to \(\{ \overline{\varvec{\upsilon }}, \overline{\varvec{\kappa }}, \varvec{\upsilon }, \varvec{\kappa }\}\) and therefore, it reads \(\overline{R}^\prime = \overline{R}^{\prime }_{\overline{\upsilon }} + \overline{R}^{\prime }_{\overline{\kappa }} + \overline{R}^{\prime }_{\upsilon } + \overline{R}^{\prime }_{\kappa }\). All four contributions can be extracted from the total variation of Eq. (37)

$$\begin{aligned} \begin{aligned} \delta {\overline{R}}&(\overline{\varvec{\upsilon }}, \overline{\varvec{\kappa }}, \varvec{\upsilon }, \varvec{\kappa }; \overline{\varvec{\eta }}) = \delta { \int _{\overline{\varvec{\mathcal {K}}}} \overline{\mathbf{P }}_\text {I}:\overline{\mathbf{F }}^{\prime }_{\overline{\upsilon }}\; \text {d}V } - \overline{F}^{\prime }_{\overline{\kappa }}(\overline{\varvec{\kappa }}; \overline{\varvec{\eta }}, \delta {\varvec{\kappa }}) \\&= \int _{\varvec{\mathcal {R}}} \delta {\left[ \overline{\mathbf{P }}_\text {I}:\overline{\mathbf{F }}^{\prime }_{\overline{\upsilon }} J_{\overline{K}}\right] } \; \text {d}V_{\theta } - \overline{F}^{\prime }_{\overline{\kappa }}(\overline{\varvec{\kappa }}; \overline{\varvec{\eta }}, \delta {\varvec{\kappa }}) \\&= \int _{\varvec{\mathcal {R}}} \left[ \left( \delta {\overline{\mathbf{P }}_\text {I}} :\overline{\mathbf{F }}^{\prime }_{\overline{\upsilon }} + \overline{\mathbf{P }}_\text {I}:\delta {\overline{\mathbf{F }}^{\prime }_{\overline{\upsilon }}} \right) \, J_{\overline{K}}+ \ \overline{\mathbf{P }}_\text {I}:\overline{\mathbf{F }}^{\prime }_{\overline{\upsilon }} \delta {J_{\overline{K}}} \right] \; \text {d}V_{\theta } \\&\qquad - \overline{F}^{\prime }_{\overline{\kappa }}(\overline{\varvec{\kappa }}; \overline{\varvec{\eta }}, \delta {\varvec{\kappa }}). \\ \end{aligned} \end{aligned}$$
(38)

Here, the total variation (or mixed variation) of the variation of the deformation gradient \(\delta {\overline{\mathbf{F }}^{\prime }_{\overline{\upsilon }}}\) and the variation of the Jacobian \(\delta {J_{\overline{K}}}\) are provided from sensitivity analysis on single scales [7]. At this point, the variation of the contribution of external loads on the macroscale \(\overline{F}^{\prime }_{\overline{\kappa }}\) is neglected. The major difference in comparison to sensitivity analysis on single scales is the total variation of stresses, especially here of the homogenised stresses \(\overline{\mathbf{P }}_\text {I}\), which contains the following partial variations

$$\begin{aligned} \begin{aligned} \delta {\overline{\mathbf{P }}_\text {I}}&= \left( \overline{\mathbf{P }}_\text {I}\right) ^{\prime }_{\overline{\upsilon }} + \left( \overline{\mathbf{P }}_\text {I}\right) ^{\prime }_{\overline{\kappa }} + \left( \overline{\mathbf{P }}_\text {I}\right) ^{\prime }_{\upsilon } + \left( \overline{\mathbf{P }}_\text {I}\right) ^{\prime }_{\kappa } \\&= \frac{\partial \overline{\mathbf{P }}_\text {I}}{\partial \overline{\varvec{\upsilon }}} \delta {\overline{\varvec{\upsilon }}} + \frac{\partial \overline{\mathbf{P }}_\text {I}}{\partial \overline{\varvec{\kappa }}} \delta {\overline{\varvec{\kappa }}} + \frac{\partial \overline{\mathbf{P }}_\text {I}}{\partial \varvec{\upsilon }} \delta {\varvec{\upsilon }} + \frac{\partial \overline{\mathbf{P }}_\text {I}}{\partial \varvec{\kappa }} \delta {\varvec{\kappa }}. \end{aligned} \end{aligned}$$
(39)

Inserting the obtained relations in Eq. (38), an appropriate rearranging leads to

$$\begin{aligned}&\delta {\overline{R}} (\overline{\varvec{\upsilon }}, \overline{\varvec{\kappa }}, \varvec{\upsilon }, \varvec{\kappa }; \overline{\varvec{\eta }}) = \int _{\overline{\varvec{\mathcal {K}}}} \left( \overline{\mathbf{P }}_\text {I}\right) ^{\prime }_{\overline{\upsilon }} :\overline{\mathbf{F }}^{\prime }_{\overline{\upsilon }}\; \text {d}V \nonumber \\&\qquad + \int _{\overline{\varvec{\mathcal {K}}}} \left( \overline{\mathbf{P }}_\text {I}\right) ^{\prime }_{\overline{\kappa }} :\overline{\mathbf{F }}^{\prime }_{\overline{\upsilon }} + \overline{\mathbf{P }}_\text {I}:\overline{\mathbf{F }}^{\prime \prime }_{\overline{\upsilon }\, \overline{\kappa }} + \overline{\mathbf{P }}_\text {I}:\overline{\mathbf{F }}^{\prime }_{\overline{\upsilon }}\, {{\text {Div}}} \delta {\overline{\varvec{X}}} \; \text {d}V - \overline{F}^{\prime }_{\overline{\kappa }} \nonumber \\&\qquad + \int _{\overline{\varvec{\mathcal {K}}}} \left( \overline{\mathbf{P }}_\text {I}\right) ^{\prime }_{\upsilon } :\overline{\mathbf{F }}^{\prime }_{\overline{\upsilon }}\; \text {d}V \quad + \int _{\overline{\varvec{\mathcal {K}}}} \left( \overline{\mathbf{P }}_\text {I}\right) ^{\prime }_{\kappa } :\overline{\mathbf{F }}^{\prime }_{\overline{\upsilon }}\; \text {d}V. \end{aligned}$$
(40)

The first integral in Eq. (40) can be identified as the macroscopic tangent stiffness operator \(\overline{k}\) from Eq. (36). The second integral can be directly connected to the macroscopic tangent pseudo load operator \(\overline{p}\) from Eq. (36) and has a similar structure to the pseudo load operator within the sensitivity analysis on single scales. The last two terms in Eq. (40) contain partial variations of the macroscopic physical residual with respect to microscopic state and design (geometry) parameters, i.e. with respect to \(\varvec{\upsilon }\) and \(\varvec{\kappa }\) and can be identified with the multilevel stiffness operator and the multilevel tangent pseudo load operator , both introduced in Eq. (36).

4.4 Discretisation of variation of macroscopic weak form

Discretisation schemes for the state and geometry functions are employed to generate the discrete counterpart of the total variation of the weak form, see Sects. 2.3.7 and 2.3.8 for the results on a single scale. Here, \(\overline{{\varvec{\mathsf{{ u}}}}}\) and \({\varvec{\mathsf{{ u}}}}\) denote the matrix of all discrete state parameters of the macro- and microscale, respectively. Similarly, \(\overline{{\varvec{\mathsf{{ x}}}}}\) and \({\varvec{\mathsf{{ x}}}}\) are the matrices of discrete geometry parameters of the macro- and microscale, respectively. Thus, the discretisation of Eq. (35) reads

(41)

In case of the finite element method (FEM), the state parameters are nodal displacements and the geometry parameters are the nodal coordinates of the macro- and microscopic FE mesh, respectively.

4.5 Discrete sensitivity of the macroscopic physical state

To avoid confusion, prior to the evaluation of the sensitivity relation the following sensitivity operators are introduced

(42)

The sensitivity matrix is the discretisation of the continuous sensitivity operator, cf. Eq. (14). After a standard discretisation and the solution of the stated macroscopic BVP, each sensitivity matrix can be computed. Both, the macroscopic sensitivity matrix \(\overline{{{{{\varvec{\mathsf{{ S}}}}}}}}\) and the multilevel sensitivity matrix are compiled to a global sensitivity matrix . According to the list of sensitivity operators in Eq. (42), the sensitivity of the macroscopic state parameters reads

(43)

Here, the sensitivity relation \(\delta {{\varvec{\mathsf{{ u}}}}} = {{{{\varvec{\mathsf{{ S}}}}}}}\, \delta {{{{\varvec{\mathsf{{ \textsf {x}}}}}}}}\) for the discrete microscopic state parameters, which is evaluated in the equilibrium point of each stated microscopic BVP, is used. Furthermore, the multilevel terms are stored in a total multilevel pseudo load matrix , which represents the total derivative of the macroscopic residual \(\overline{{\varvec{\mathsf{{ R}}}}}\) with respect to microscopic geometry parameters \({{{\varvec{\mathsf{{ \textsf {x}}}}}}}\), i.e. \({\text {d}\overline{{\varvec{\mathsf{{ R}}}}}}/{\text {d}{{{\varvec{\mathsf{{ \textsf {x}}}}}}}}\). Importantly, this quantitiy is not computed by discrete derivatives but by the discretisation of continuous variations.

4.6 Variational sensitivity of arbitrary functionals

The total variation of an arbitrary objective or constraint functional \(f(\overline{\varvec{\upsilon }}, \overline{\varvec{\kappa }}, \varvec{\upsilon }, \varvec{\kappa })\), i.e. which is formulated on two scales, can be derived based on the preliminary remarks, cf. Eq. (15). Similarly to Eq. (36), the linear operators

(44)

are derived and composed to the total variation of the function

(45)

4.7 Discrete form of the variation of arbitrary functionals

The continuous formulation of the variation of an arbitrary objective or constraint functional \(f(\overline{\varvec{\upsilon }}, \overline{\varvec{\kappa }}, \varvec{\upsilon }, \varvec{\kappa })\) formulated on two scales, i.e. Eq. (45), can be transferred to its discrete counterpart using discretised sensitivity operators from Eqs. (43) to (42), respectively. Dependencies on quantities from two different scales have to be taken into account. The total geometry variation with respect to macroscopic and microscopic geometry, i.e. with respect to \(\overline{{{{\varvec{\mathsf{{ \textsf {x}}}}}}}}\) and \({{{\varvec{\mathsf{{ \textsf {x}}}}}}}\), is outlined using the notation of discrete derivatives of f, i.e. it reads

(46)

This notation is used to highlight the contents of the different continuous operators, i.e. the variations are performed and discretised. This relationship can be used for derivations and evaluations of several objective or constraint functionals on two scales within frameworks for two-scale structural optimisation. Especially their gradients or general sensitivity information are available.

4.8 Variation of effective parameters

To obtain the full formulation of the total variation of the macroscopic weak form of equilibrium from Eq. (40), variations of homogenised stresses are required. In accordance to the minimisation problem of homogenisation in terms of the Lagrange formalism from Eq. (21) and deduced homogenised stresses in Eq. (22) or in Eq. (31), effective stresses \(\overline{\mathbf{P }}_\text {I}\) or \(\overline{{{\underline{{{\varvec{\mathsf{{ P}}}}}}}}}_{\text {I}}\), respectively, explicitly depend on the Lagrange multiplier \({{\varvec{\uplambda }}}_\text {I}\) and their variation reads

$$\begin{aligned} \begin{aligned} \left( \overline{{{\underline{{{\varvec{\mathsf{{ P}}}}}}}}}_{\text {I}}\right) ^{\prime }_{\upsilon } + \left( \overline{{{\underline{{{\varvec{\mathsf{{ P}}}}}}}}}_{\text {I}}\right) ^{\prime }_{\kappa }&= \frac{\partial \overline{{{\underline{{{\varvec{\mathsf{{ P}}}}}}}}}_{\text {I}}}{\partial {{\varvec{\uplambda }}}_{\text {I}}} ({{\varvec{\uplambda }}}_\text {I})^\prime = \frac{\partial \overline{{{\underline{{{\varvec{\mathsf{{ P}}}}}}}}}_{\text {I}}}{\partial {{\varvec{\uplambda }}}_{\text {I}}} \left[ \, ({{\varvec{\uplambda }}}_{\text {I}})^{\prime }_{\upsilon } + ({{\varvec{\uplambda }}}_{\text {I}})^{\prime }_{\kappa } \, \right] \\&= \frac{\partial \overline{{{\underline{{{\varvec{\mathsf{{ P}}}}}}}}}_{\text {I}}}{\partial {{\varvec{\uplambda }}}_{\text {I}}} \left[ \, \frac{\partial {{\varvec{\uplambda }}}_{\text {I}}}{\partial {{{\varvec{\mathsf{{ u}}}}}}}\, \delta {{{{\varvec{\mathsf{{ u}}}}}}}+ \frac{\partial {{\varvec{\uplambda }}}_{\text {I}}}{\partial {{{\varvec{\mathsf{{ \textsf {x}}}}}}}}\, \delta {{{{\varvec{\mathsf{{ \textsf {x}}}}}}}}\, \right] . \end{aligned} \end{aligned}$$
(47)

This relation can be connected to the sensitivity of the micro-macro coupling condition. The sensitivity of the Lagrange multiplier \({{\varvec{\uplambda }}}_\text {I}\) has to be investigated. Referring the solution scheme on the lower scale, in the solution point of stated microscopic BVP the obtained Lagrange multiplier \({{\varvec{\uplambda }}}_\text {I}\), especially directly in the case of linear displacements (D) for boundary conditions, corresponds to resulting reaction forces or the external part of the residual on the surface of the underlying RVE scaled by the value of the microscopic volume. For different types of boundary conditions, like periodic displacements (P) or uniform tractions (S), the resulting sensitivity relations have to be adapted. Overall, the total design variation of the Lagrange multiplier for alternative boundary conditions \(\text {I}= \{\text {D},\ \text {P},\ \text {S}\}\) in discrete form is obtained from

$$\begin{aligned} V\mathcal {S}_1^T ({{\varvec{\uplambda }}}_\text {I})^\prime = \left[ {{{{\varvec{\mathsf{{ K}}}}}}}^{\textsf {int}}_{\text {b}\text {a}}\ {{{{\varvec{\mathsf{{ S}}}}}}}_{\text {a}} + {{{{\varvec{\mathsf{{ P}}}}}}}^{\textsf {int}}_{\text {b}} \right] \delta {{\varvec{\mathsf{{ x}}}}}. \end{aligned}$$
(48)

Compared to the set \((\text {i},\, \text {b})\) introduced in Sect. 3.2 for inner and boundary nodes of an RVE, the nodes on the Dirichlet boundary given in Eq. (48) are also indicated by the index \((\text {b})\). Whereas the set identified by the index \((\text {a})\) contains all nodes except boundary nodes, i.e. it contains nodes in the inner domain as well as nodes on the Neumann boundary. The superscript \((\cdot )^{\textsf {int}}\) indicates the contribution of the internal part of the residual and therefore, it indicates the contribution of the resulting stiffness \({{{{\varvec{\mathsf{{ K}}}}}}}^{\textsf {int}}\) and pseudo load \({{{{\varvec{\mathsf{{ P}}}}}}}^{\textsf {int}}\).

The chosen homogenisation approach requires the incorporation of the obtained sensitivity information into the context of homogenised stresses in Eq. (31) and formulated boundary conditions presented in Sect. 3. The variation of effective stresses is presented in Eq. (47). In the computational sense, they are adapted for different boundary conditions and are computed using the general formulation of the boundary conditions matrix \(\mathcal {S}_2\), i.e. \(\overline{{{\underline{{{\varvec{\mathsf{{ P}}}}}}}}}_{\text {I}} ({{\varvec{\uplambda }}}_\text {I}) = \mathcal {S}_2^T {{\varvec{\uplambda }}}_\text {I}\). In comparison to Eq. (47), the boundary conditions matrix \(\mathcal {S}_2\) can be identified with the partial derivative of effective stresses with respect to the Lagrange multiplier \({{\varvec{\uplambda }}}_\text {I}\). Due to the fact that \(\mathcal {S}_2\) is constant for all introduced types of boundary conditions, the variation of homogenised stresses contains only the variation of \({{\varvec{\uplambda }}}_\text {I}\) and therefore, its final form results to

$$\begin{aligned} \begin{aligned} \left( \overline{{{\underline{{{\varvec{\mathsf{{ P}}}}}}}}}_{\text {I}}\right) ^{\prime }_{\upsilon } \!+\! \left( \overline{{{\underline{{{\varvec{\mathsf{{ P}}}}}}}}}_{\text {I}}\right) ^{\prime }_{\kappa }&= \mathcal {S}_2^T ({{\varvec{\uplambda }}}_\text {I})^\prime = \frac{1}{V}\, \mathcal {S}_2^T\, \mathcal {S}_1^{-T} \left[ {{{{\varvec{\mathsf{{ K}}}}}}}^{\textsf {int}}_{\text {b}\text {a}}\ {{{{\varvec{\mathsf{{ S}}}}}}}_{\text {a}} + {{{{\varvec{\mathsf{{ P}}}}}}}^{\textsf {int}}_{\text {b}} \right] \delta {{{{\varvec{\mathsf{{ \textsf {x}}}}}}}}. \end{aligned} \end{aligned}$$
(49)

5 Optimisation setup for problems on two scales

After the solution of the stated BVP in terms of structural analysis, the subsequently performed design sensitivity analysis gives the opportunity to deal with following three questions and investigations: How will physical responses in \(\overline{\varvec{\mathcal {K}}}\) react,

  • if design parameters in \(\overline{\varvec{\mathcal {K}}}\) change?

  • if design parameters in \(\varvec{\mathcal {K}}\) change?

  • if design parameters in \(\overline{\varvec{\mathcal {K}}}\) and \(\varvec{\mathcal {K}}\) change?

These questions represent a central investigation and are evaluated in [33] within several numerical examples. Some possible model problems with different microstructures are schematically pictured in Fig. 4. Each individual microscale RVE can be embedded into the macroscale domain and investigated within several optimisation procedures.

Fig. 4
figure 4

Arbitrary subset of model problems with different microstructures for structural optimisation on two scales

Full size image

Beside the necessity of sensitivity relations for the coupling conditions of different scales, the choice of appropriate objective functionals and constraints in combination with useful design parameters is essential. Table 2 lists possible sets for objectives, constraints and design parameters, which are common and often used in standard publications on structural optimisation, see [3, 13, 18, 28, 61] to name a few.

Table 2 Possible sets for objectives (OF) in terms of goal setting, constraints (CON) in terms of restrictions and design parameters (DP) in the design space
Full size table

A complete formulation and the numerical treatment of a two-scale optimisation problem requires an accurate specification of the two-scale optimisation problem in its discrete form, which is articulated as follows.

figure a

Here, the following matrix notations are introduced: a matrix notation \(\overline{{{{\varvec{\mathsf{{ h}}}}}}}(\overline{{{{\varvec{\mathsf{{ u}}}}}}},\, \overline{{{{\varvec{\mathsf{{ s}}}}}}},\, {{{\varvec{\mathsf{{ u}}}}}},\, {{{\varvec{\mathsf{{ s}}}}}})\) of introduced macroscopic equality constraints \(\overline{h}_i(\overline{{{{\varvec{\mathsf{{ u}}}}}}},\, \overline{{{{\varvec{\mathsf{{ s}}}}}}},\, {{{\varvec{\mathsf{{ u}}}}}},\, {{{\varvec{\mathsf{{ s}}}}}}) = 0\), \(i \in \overline{\mathcal {E}}_h\), a matrix representation \(\overline{{{{\varvec{\mathsf{{ g}}}}}}}(\overline{{{{\varvec{\mathsf{{ u}}}}}}},\, \overline{{{{\varvec{\mathsf{{ s}}}}}}},\, {{{\varvec{\mathsf{{ u}}}}}},\, {{{\varvec{\mathsf{{ s}}}}}})\) of introduced macroscopic inequality constraints \(\overline{g}_j(\overline{{{{\varvec{\mathsf{{ u}}}}}}},\, \overline{{{{\varvec{\mathsf{{ s}}}}}}},\, {{{\varvec{\mathsf{{ u}}}}}},\, {{{\varvec{\mathsf{{ s}}}}}}) \le 0\), \(j \in \overline{\mathcal {I}}_h\), a matrix notation \({{{\varvec{\mathsf{{ h}}}}}}({{{\varvec{\mathsf{{ u}}}}}},\, {{{\varvec{\mathsf{{ s}}}}}})\) of introduced microscopic equality constraints \(h_k({{{\varvec{\mathsf{{ u}}}}}},\, {{{\varvec{\mathsf{{ s}}}}}}) = 0\), \(k \in \mathcal {E}_h\), and a matrix representation \( {{{\varvec{\mathsf{{ g}}}}}}({{{\varvec{\mathsf{{ u}}}}}},\, {{{\varvec{\mathsf{{ s}}}}}}) \) of introduced microscopic inequality constraints \(g_l({{{\varvec{\mathsf{{ u}}}}}},\, {{{\varvec{\mathsf{{ s}}}}}}) \le 0\), \(l \in \mathcal {I}_h\). The sets of indices for equality and inequality constraints are denoted by \(\overline{\mathcal {E}}_h,\, \mathcal {E}_h,\, \overline{\mathcal {I}}_h\) and \(\mathcal {I}_h\). The application of the general SQP method as a solver for the stated optimisation problem, which guarantees solutions in feasible regions with respect to the stated constraints, requires a compact and summarising representation of all introduced macro- and microscopic constraints and design parameters. Therefore, the following representation is introduced

(50)

The Lagrange function used within SQP as well as its gradient are modified to

(51)

For the solution of the optimality criteria using the SQP method, a sequential solution of quadratic subproblems with a quadratic approximation of the defined objective functional \(J(\overline{{{{\varvec{\mathsf{{ u}}}}}}},\, \overline{{{{\varvec{\mathsf{{ s}}}}}}},\, {{{\varvec{\mathsf{{ u}}}}}},\, {{{\varvec{\mathsf{{ s}}}}}})\) and the linearisation of equality and inequality constraints is required.

Fig. 5
figure 5

Two-scale optimisation of a bracket: mechanical system and FE mesh for the macro- and microscopic initial design (top line). Optimisation model for macroscale with two morphing boxes and for microscale with one morphing box (bottom line)

Full size image

Choice of design parameters In the previous section, all discrete shape sensitivity information is provided with respect to nodal coordinates \((\overline{{{{\varvec{\mathsf{{ \textsf {x}}}}}}}}, {{{\varvec{\mathsf{{ \textsf {x}}}}}}})\) of the finite element mesh. In the abstract setting of the two-scale optimisation problem in Eqs. (50) and (51), respectively, the design is represented by the vectors \(\overline{{{{\varvec{\mathsf{{ s}}}}}}}\) and \({{{\varvec{\mathsf{{ s}}}}}}\), which are a general placeholders for any arbitrary finite set of design parameters. For instance, the latter example in Sect. 6 contains coordinates of the Bézier patch of morphing boxes as design parameters of choice.

Verification of sensitivity information in numerical studies Prior to the discussion on numerical examples and obtained optimisation results, it is mentioned that all analytically obtained quantities and values within the design sensitivity analysis are verified using the finite difference method (FDM). These numerical tests guarantee correct gradient information and therefore, they guarantee an efficient optimisation process. The verifications are performed in accordance to the appropriate quantity either on the element (e.g. pseudo load matrices), on the global system level (e.g. sensitivities of the state) or after the transformation with the mesh design matrices (e.g. final objectives and constraints).

6 Example: two-scale optimisation of a bracket

The following example illustrates the applicability of the formulated sensitivity relations and corresponding numerical aspects from previous chapters. The stated optimisation problem focuses on the minimisation of the overall macroscopic compliance, which is directly related to the maximisation of the overall macroscopic stiffness. In addition to some constraints for defined design parameters, volume constraints on both scales are incorporated into the optimisation process. Additionally, a morphing based design parametrisation on the macro- and the microscale \(\overline{\varvec{\mathcal {K}}}\) and \(\varvec{\mathcal {K}}\) is used to perform shape optimisation on both scales simultaneously. Design parametrisation based on morphing is useful to handle forms and geometries without explicit geometrical properties, like diameters, positions and angles. Although morphing boxes have to be selected manualy, they provide a convenient means to specify which geometry degrees of freedom are subject to optimization. Furthermore, the morphing parametrisation technique provides the opportunity to reduce the dimension of the discrete geometry design space to make the numerical optimisation problem more tractable. It is basically established and known from image processing and is based on smooth and continuous transformations of target objects into other objects of interest. A review on different aspects and approaches can be found in [78] and extensions to three-dimensional investigations are presented in [66, 79]. Overall, mesh morphing techniques within mechanical problems and structural optimisation procedures on single scales are applied in [58] or in [25] using commercially-available software tools for instance. In the context of a two-scale shape optimisation scheme, its first use is presented in the work at hand. Extended explanations can also be found in [33]. In general, mesh morphing techniques can be classified as a combination of parameter free and CAGD based optimisation techniques. Consequently, transformations using mesh design matrices are necessary to obtain accurate information in the design space of interest.

Due to the wide range of applicability and practical relevance, the example is evaluated using the formulation for periodic boundary conditions (P) on the microscale.

The iterative solution in terms of structural analysis is provided by the application of the Newton’s method on the macro- and the microscale. The computation is based on the hyperelastic Neo-Hookean material law for the constituents on the microscale and are stopped after the criterion of \(\varepsilon _\upsilon = 10^{-8}\) for the residuals is reached. The iterative solution within the structural optimisation is carried out by the SQP method. The change in the objective function is defined as stop criterion and is set to \(\varepsilon _s= 10^{-4}\).

For the finite element analysis a half of a bracket on the macroscale with main dimensions A and B and a microscale with the main dimension a consisting of a stiff matrix and a softer kernel material is modeled, cf. Fig. 5. The illustrated microscale representation is connected with each macroscopic integration point and allows the computation of effective stress and material parameters in terms of computational homogenisation. Characteristic system parameters for both scales, like dimensions and material properties, are compiled in Table 3. The used finite element technology is a pure displacement formulation with linear shape functions and two degrees of freedom per node for triangular element types on the macroscale and rectangular element types on the microscale, respectively.

Table 3 Optimisation of a bracket: model parameters
Full size table

The optimisation setup with all design parameters on the macro- and the microscale also can be found in Fig. 5. On both scales selected coordinates of control points of defined morphing boxes are chosen as design parameters. On the microscale one morphing box with 25 control points is defined. The coordinates of the inner control points of the morphing box (indicated by numbers 1–8 in Fig. 5) are chosen as design variables. They can move nearly to the defined boundary of the domain. The control point in the center of the domain is neglected as design parameter due to symmetry properties of the system. On the macroscale two morphing boxes are defined. The first one is around the loading area and consists of 9 control points. Here, coordinates of only one control point are chosen as design variables (indicated by number 9 in Fig. 5). The second morphing box is defined over the main and regular part of the system. It contains 20 control points. For this morphing box, coordinates of 11 control points (indicated by the numbers 10–20 in Fig. 5) are chosen as design parameters within the optimisation procedure. Due to the characteristics of morphing based design parametrisation, see Sect. 2.3.7 for general remarks on functional dependencies, the following dependencies hold true

Fig. 6
figure 6

Two-scale optimisation of a bracket: objective (top) corresponds to macroscopic compliance \(\overline{C}\) and constraints (bottom) correspond to volume constraints \(\overline{V},V\) on the macro- and the microscale. Comparison of initial and optimised values of objectives for performed example results to: \(\overline{C}\, ^\text {opt}/\ \overline{C}\, ^\text {ini}= {0.329}\ /\ {0.398} = 0.83\). The volumes of referred domains remain constant, i.e. \(\overline{V}\, ^\text {ini}/ \overline{V}\, ^\text {opt}= V\, ^\text {ini}/ V\, ^\text {opt}= 1\)

Full size image
Fig. 7
figure 7

Two-scale optimisation of a bracket: distribution of final design parameters for the macroscopic optimisation model with two morphing boxes (left) and the microscopic optimisation model with one morphing box (right)

Full size image
$$\begin{aligned} \begin{aligned} \text {macroscale: }&\quad \overline{(\cdot )}(\overline{{{{\varvec{\mathsf{{ u}}}}}}}(\, \overline{{{{\varvec{\mathsf{{ \textsf {x}}}}}}}}(\, \overline{{{{\varvec{\mathsf{{ \textsf {y}}}}}}}} (\overline{{{{\varvec{\mathsf{{ s}}}}}}} ))), \overline{{{{\varvec{\mathsf{{ \textsf {x}}}}}}}}(\, \overline{{{{\varvec{\mathsf{{ \textsf {y}}}}}}}} (\overline{{{{\varvec{\mathsf{{ s}}}}}}}))) \quad \text {and} \\ \text {microscale: }&\quad (\cdot )({{{\varvec{\mathsf{{ u}}}}}}({{{\varvec{\mathsf{{ \textsf {x}}}}}}}( {{{\varvec{\mathsf{{ \textsf {y}}}}}}}({{{\varvec{\mathsf{{ s}}}}}})) ), {{{\varvec{\mathsf{{ \textsf {x}}}}}}}({{{\varvec{\mathsf{{ \textsf {y}}}}}}}({{{\varvec{\mathsf{{ s}}}}}}))) \, . \end{aligned} \end{aligned}$$
(52)

In detail, \((\overline{{{{\varvec{\mathsf{{ u}}}}}}}, {{{\varvec{\mathsf{{ u}}}}}})\) are the macro- and microscopic state, \((\overline{{{{\varvec{\mathsf{{ \textsf {x}}}}}}}}, {{{\varvec{\mathsf{{ \textsf {x}}}}}}})\) are the macro- and the microscopic coordinates of the FE mesh, \((\overline{{{{\varvec{\mathsf{{ \textsf {y}}}}}}}}, {{{\varvec{\mathsf{{ \textsf {y}}}}}}})\) are the coordinates of the macro- and the microscopic CAGD parameters and \((\overline{{{{\varvec{\mathsf{{ s}}}}}}}, {{{\varvec{\mathsf{{ s}}}}}})\) are the finite subsets of design variables. Therefore, the following transformations of sensitivity information for arbitrary functionals \(\overline{f}\) on the macroscale and \(f\) on the microscale are necessary

Fig. 8
figure 8

Two-scale optimisation of a bracket: macroscopic displacements (large) for initial and optimised supplemented designs. Microscopic von Mises stress distribution in selected macroscopic integration points (1–3) for initial and optimised design. Due to symmetry properties, microstructures in points (1–3) can be copied to design the full system

Full size image
(53)
$$\begin{aligned} f^\prime = \left( \frac{\partial f}{\partial {{{\varvec{\mathsf{{ u}}}}}}}\ {{{{\varvec{\mathsf{{ S}}}}}}}+ \frac{\partial f}{\partial {{{\varvec{\mathsf{{ \textsf {x}}}}}}}} \right) \frac{\partial {{{\varvec{\mathsf{{ \textsf {x}}}}}}}}{\partial {{{\varvec{\mathsf{{ s}}}}}}} \delta {{{{\varvec{\mathsf{{ s}}}}}}} = \left( \frac{\partial f}{\partial {{{\varvec{\mathsf{{ u}}}}}}}\ {{{{\varvec{\mathsf{{ S}}}}}}}+ \frac{\partial f}{\partial {{{\varvec{\mathsf{{ \textsf {x}}}}}}}} \right) {{{{\varvec{\mathsf{{ V}}}}}}}\delta {{{{\varvec{\mathsf{{ s}}}}}}}. \end{aligned}$$
(54)

Here, \( \overline{{{{{\varvec{\mathsf{{ V}}}}}}}}\) and \({{{{\varvec{\mathsf{{ V}}}}}}}\) are the mesh design matrices of both scales introduced in Sect. 2.3.7. In both cases the appropriate sensitivity matrix for the macro- and \({{{{\varvec{\mathsf{{ S}}}}}}}\) for the microscale is required. Equations (53) and (54) provide the transformation of sensitivity information calculated with respect to nodal coordinates of the finite element mesh, i.e. of \({\partial (\cdot )}/{\partial \overline{{{{\varvec{\mathsf{{ s}}}}}}}}\) on the macro- or \({\partial (\cdot )}/{\partial {{{\varvec{\mathsf{{ s}}}}}}}\) on the microscale, into the design space of interest.

The stated optimisation problem introduced in Sect. 5 contains the minimisation of an objective, which is here the macroscopic compliance \(\overline{C}\). Its minimisation is performed with respect to constant volume constraints \((\overline{V}, V)\) and lower and upper bounds \((\overline{s}^\text {l}_i,\ \overline{s}^\text {u}_i,\ s^\text {l}_i,\ s^\text {u}_i)\) for the design parameters on the macro- and the microscale. The solution is obtained by the application of the SQP method with an BFGS update for second derivatives.

To provide a well posed optimisation problem and to force a modification of the design, the contribution of each individual constituent is weighted, i.e. the overall volume is computed by \(V= \sum _{i=1}^k w_i V_i\). Otherwise, difficulties during the evaluation of the constraints and of their derivatives arise. The weight factors \(w_i\) can be related to the importance, overall costs or valency of each constituent in an abstract sense. The weight for the matrix material is set to \(w_1 = 2.0\) and for the softer kernel material to \(w_2=1.0\).

The optimisation results for this two-scale optimisation problem of a macroscopic bracket with a heterogeneous microstructure are plotted in Fig. 6. The objective can be minimised to a value that is \(17 \%\) lower than the initial value and the volume constraints on both scales are fulfilled. This means that the overall volume on the macroscale as well as the overall volume on the microscale remain constant. The optimisation algorithm takes 23 iterations to find a feasible solution.

The distribution of design parameters on both scales for the initial and optimised shape can be found in Fig. 7. All design parameters fulfill their constraints and lie in the feasible domain. The comparison of the initial and optimised designs in terms of overall performance is illustrated in Fig. 8. The results obtained for one half of the system can be mirrored due symmetry properties to design the overall system. It is shown that minimisation of the macroscopic compliance leads to a lower maximum value of the macroscopic displacements, i.e. \({\overline{u}^\text {opt}_\text {max}}\ /\ {\overline{u}^\text {ini}_\text {max}} = {0.0593}\ /\ {0.0733} = 0.81\), which is in the range of the reduction of overall compliance. Furthermore, the design of the obtained microscale representation leads to a lower distribution of von Mises stresses in the selected macroscopic integration points. For instance, the maximum von Mises stress in point 1 can be reduced by 6 %, in point 2 by 18 % and in point 3 even by 43 %, cf. Fig. 8. This example proves the applicability of the presented relations for the sensitivity analysis over two scales and that it comes along with a large gain in terms of overall behaviour of investigated systems.

7 Conclusion

The paper at hand outlines a framework for the design of structures and materials over two length scales. The formulation of the two-scale optimisation problem is compiled based on relations for numerical homogenisation and \(\hbox {FE}^{2}\) methods, methods for two-scale sensitivity analysis and methods for mathematical optimisation. Special emphasize is given to the determination of the sensitivity information of effective parameters, which in fact bridge underlying scales. An example for the minimisation of the overall macroscopic compliance with constraints and design parameters on different scales proves, that it is worthwhile to incorporate the material scale into the overall design process.

Detailed explanations on mentioned approaches and formulations as well as further illustrative examples can be found in the PhD thesis of the author Wojciech Kijanski, see [33].