A multiscale topology optimisation framework for hollow spheres as cellular materials Bruno G. Christo ff, Jos´e Humberto S. Almeida Jr., Eduardo L. Cardoso, V olnei Tita Ω∂Ωu∂Ωf x1,y1x2,y2 x3,y3 A multiscale topology optimisation framework for hollow spheres as cellular materials Bruno G. Christo ffa, Jos´e Humberto S. Almeida Jr.b,∗, Eduardo L. Cardosoc, V olnei Titaa,d aDepartment of Aeronautical Engineering, Sao Carlos School of Engineering, University of Sao Paulo, Sao Carlos, SP , Brazil bSchool of Mechanical and Aerospace Engineering, Queen’s University Belfast, Belfast, UK cDepartment of Mechanical Engineering, State University of Santa Catarina, Joinville, Brazil dDepartment of Mechanical Engineering, Faculty of Engineering of the University of Porto, Porto, Portugal Abstract Cellular materials can be designed to achieve mechanical properties with low weight. Among them, hollow spheres have been investigated since modern manufacturing methods can accurately produce them. This work aims at developing an optimisation procedure to distribute hollow spheres to min- imise the compliance of a structure through a multiscale approach. The model has two geometrical parameters: internal and external diameters, in which the asymptotic homogenization method (AHM) is employed to predict the e ffective properties of the material. Equations relating the e ffective prop- erties to geometrical parameters are obtained through the least square method. Then, a topology optimisation is employed to minimise compliance constrained by an admissible relative density. The optimised structures are compared to homogeneous distributions of hollow spheres, and compliance reductions up to 79% are reached. The optimum distributions of hollow spheres are validated against the classical solid isotropic with material penalisation (SIMP) approach. Moreover, some structures are successfully 3D–printed to show the feasibility and capabilities of the proposed approach. Highlights •Cellular materials for improved specific mechanical properties; •A topology optimisation framework for Engineering Structures; •A multiscale approach to obtain the optimal distribution of cellular materials; •Compliance reduction of up to 79% when compared to non-optimised solution; •Feasibility of the proposed method demonstrated by 3D printing hollow spheres. Keywords: Cellular materials, Hollow spheres, Multiscale approach, Asymptotic homogenization method, Topology optimisation ∗Corresponding author Email addresses: brunochristoff@usp.br (Bruno G. Christo ff),humberto.almeida@qub.ac.uk (Jos´e Humberto S. Almeida Jr. ), eduardo.cardoso@udesc.br (Eduardo L. Cardoso), voltita@sc.usp.br (V olnei Tita) Preprint submitted to Engineering Structures September 28, 2022 1. Introduction Cellular materials are formed, necessarily, by two phases [1, 2]: the first is composed of solid and continuous material and the second lacks material, that is, void. These materials can roughly be divided into two distinct groups: open and closed cell materials [3]. In closed-cell materials, voids are discontinuous and dispersed into the solid phase. An open-cell structure presents a continuous void phase within the structure [1]. In practice, most cellular materials present characteristics of both open and closed–cell structures simultaneously. Another way of defining a cellular material is through the concept of a composite. Considering that a composite is a set of two or more di fferent materials combined into a macroscopic scale, cellular materials can be seen as a particular type of composite, in which one of the constituent materials is solid and the one is void [4]. The most important feature of cellular materials is their relative density, i.e., the density of the cel- lular material divided by the density of the solid of which the cellular material is made [5] since these materials can be designed to achieve very small relative densities. According to Gibson and Ashby [5], the four major areas of application of cellular materials are thermal insulation, packaging, struc- tural, and buoyancy. The cellular structure provides excellent properties to the structure, such as high energy absorption, good formability, and excellent insulation [6]. In addition, such cellular structures allow tailoring their mechanical, thermal, and acoustic properties, which aids in multi-functionality, such as in thermo-structural applications [7]. As these kinds of materials can be designed to achieve a broad range of properties, they can be used in several industrial applications. Cellular materials are widely used in thermal conductivity applications [8, 9], damping and vibration control [10, 11], structural applications [12, 13, 14, 15] and frequency shielding [16, 17]. A type of cellular material is hollow spheres [18, 19, 20]. These structures provide a specific mi- crostructure with the likelihood of enhanced properties for low-density lightweight structures [21, 22]. This type of structure is characterised by easily reproducible geometry, guaranteeing physical con- sistency in the properties of the material [23, 24, 25]. Additionally, hollow sphere structures can be manufactured from a broad range of materials and can be assembled in a periodic structure [26]. Au- gustin et al. [27] show the technological development of manufacturing processes for hollow spheres focusing on their production and prospective applications. Waag et al. [28] present several manu- facturing processes for hollow spheres made of porous metals, ceramics and polymers. Kalamar et al. [29] performed a structural analysis of full-scale hollow, laminated glass columns under combined 3 in-plane compression and impact loads. Liu et al. [30] investigated the residual mechanical properties of hollow spherical joints with H-beam (H-WHSJ) after loaded in elevated temperatures. Gasser et al. [31] predicted the elastic behaviour of a hollow sphere foam using a model derived from the thin shell elastic theory. Marcadon and Feyel [32] studied the influence localised plasticity and buckling of hollow spheres on the e ffective behaviour of their packing. Yan et al. [33] assessed the mechanical behaviour of H-shaped steel member (WHS-HSM) joints under axial compression. Hosseini et al. [34] investigated the thermal properties of a perforated hollow sphere structure. Solorzano et al. [35] analysed the thermal conductivity of hollow sphere structures analytically. Fiedler and Solorzano [36] used Monte Carlo and finite element methods to predict the e ffective thermal conductivity properties of randomly oriented hollow sphere structures. Marcadon [37] aimed at exploiting the potential of using hollow spheres in high-temperature applications, wherein elastic, plastic, and viscous properties on the e ffective behaviour of hollow sphere structure were investigated. Liu et al. [38] exploited the dynamic crushing of hollow spheres, Marcadon [39] evaluated the influence of geometrical defects on the mechanical response of hollow sphere structures, whereas Wu et al. [40] focused on tailoring the sound absorption properties. Shufrin et al. [41] investigated the negative Poisson’s ratio in auxetic cellular structures. After this non–extensive state-of-the-art, three main gaps are found: i) How to predict the ef- fective material properties of a media formed by hollow spheres? ii) How to optimally distribute hollow spheres in a domain to maximise their performance? iii) Is it possible to manufacture the structures with the design variables considered? To answer these research questions, this paper aims at proposing an original optimisation framework to obtain the optimal distribution of material for hollow spheres within a certain domain to minimise their compliance. The problem is tackled in a multiscale environment and it is constrained by a maximum allowable relative density. The global domain is locally represented by means of a Representative V olume Element (RVE). An idealised model for a single hollow sphere is also proposed, which is described by its geometrical parameters. The e ffective properties are obtained by the asymptotic homogenization method by using the Least Square Method, and they are written as a continuous and di fferentiable function of the geometrical parameters of the RVE. The framework is evaluated through two classical benchmark cases: i) can- tilever beam and ii) three–point bending. The maximum allowed relative densities considered are 5%, and 10%. 4 2. The approach An anisotropic and quasi-periodic elastic body is considered. Let Ωbe a domain in R3with boundary∂Ω(see Figure 1). Consider that Dirichlet and Neumann boundary conditions are applied on∂Ωuand∂Ωf, respectively, where ∂Ω =∂Ωui∪∂Ωfi;∂Ωui∩∂Ωfi=∅i=1,2,3, (1) such that ui=ug ion∂Ωuii=1,2,3, (2) where uiis the i-th component of the displacement and ug iis known, and σi jnj=tion∂Ωtii=1,2,3, (3) where σis the stress tensor, nis a normal vector, and tis a traction vector. The elastic body is assumed to be a cellular solid with a microstructure formed by a distribution of hollow spheres. Two scales are used to define the problem. On the macroscale, defined by the slow variable x, the mechanical properties are homogenised in the considered material point. On the microscale, represented by the fast variable y, the e ffective properties of the medium are calculated with a two-scale Asymptotic Homogenization Method (AHM), such that CH i jkl=CH i jkl(y). (4) Both coordinate systems are assumed to be coincident. In Figure 1, a small portion of the domain is highlighted, depicting the distribution of the hollow sphere material. An idealised material model is assumed, in which a regular distribution of hollow spheres is considered. A Representative V olume Element (RVE) is used to depict the medium. The geometrical param- etersρo, andρi(outer and inner diameters of the RVE, respectively), describe the entire microscale domain. In addition, two adjacent hollow spheres are connected by a flat area and it is considered that the thickness is constant and a function of the geometrical parameters. Therefore, the homogenised 5 elasticity tensor of the media can be written as CH=CH(ρo,ρi). (5) With the macroscopic properties of the domain written as a function of the RVE geometry, it is assumed that the distribution of cellular material within the macroscopic domain varies as a function of its position. This can be written as a function of the geometric parameters of the RVE as ρo=ρo(x)andρi=ρi(x). (6) An adaptation of the classical topology optimisation problem is used in the present study to find the optimum distribution of hollow spheres in the domain. Thus, every material point of the domain is described by a pair of geometrical parameters ρoandρi, describing the hollow sphere for that particular point. The optimisation problem is set as min ρo,ρif(ρo(x),ρi(x)) Subject to V(ρo(x),ρi(x))≤¯V,(7) where fis the compliance, Vis the relative density of the structure, and ¯Vis a given admissible relative density. It is also considered that the design variables are bounded, such that ρmin o≤ρo(x)≤ρmax o ρmin i≤ρi(x)≤ρmax i,∀x∈Ω.(8) Linear elastic behaviour is assumed, and Figure 2 shows the design domains considered for the optimisation cases. The AHM equilibrium relations and the entire optimisation procedure are carried out on in-house routines written in Julia Language [42]. 3. Mathematical background 3.1. Asymptotic homogenisation method In the context of linear elasticity, the e ffective properties of a material with a periodic (or quasi- periodic) microstructure can be determined by the AHM method [43]. This method is based on the assumption that the RVE, a small part of the domain which repeats itself with periodicity in the whole 6 Ω∂Ωu∂Ωf x1,y1x2,y2 x3,y3 Idealisation Figure 1: Depiction of the microscopic structure of a small portion of the macroscopic domain, idealised model and the RVE. domain, is considerably smaller than the design domain [44]. The notation, derivation of the equilibrium problems, and the finite element (FE) solution are based on [45, 46, 43, 47]. An anisotropic and quasi-periodic elastic body is considered. The problem is formulated in the bounded subset ΩofR3, and a slow variable xand a fast variable yare considered to describe the problem. The RVE is denoted by Y={y=(y1,y2,y3)∈R3: 0<yi<li,i=1,2,3} (9) where Yis the size of the RVE, with lias positive numbers. The macroscopic domain is denoted as Ω =εY={x=(x1,x2,x3)∈R3:ε−1xi∈Y,i=1,2,3}, (10) 7 Ω(a) Ω (b) Ω (c) Figure 2: Design domains considered for the optimisation problem: (a) cantilever beam, (b) 3–point bending beam with load application on the upper side, and (c) 3–point bending beam with load application on the bottom side. where y=x ε, (11) withε<< 1. The development of the method follows three basic considerations [48]. The first demands that the displacement field of the media can be written in an asymptotic expansion as uε(x)=u(0)(x)+εu(1)(x,y)+ε2u(2)(x,y)+... (12) where uεis the total displacement field and u(0),u(1)andu(2)are the contributions for the displacement of the macroscopic scale, microscopic scale, and eventual smaller scales, respectively. The second consideration is that the coordinates at each level are related by a small parameter ε, as shown in Equation 11. The third consideration is that the displacements on the boundaries of the RVE are periodic. This consideration implies that the local domain is periodic, thus describing the macroscopic domain. 8 A general field variable fε iis thus dependent on both fast and slow variables, such as fε i(x)=fi(x,y) (13) and its partial derivatives take the form ∂fε i ∂xj=∂fi ∂xj+ε−1∂fi ∂yj. (14) Assuming that no body forces are acting on the media, one can write the elastic equilibrium as ∂σε i j ∂xj=0 in Ωε, (15) where σε i j=Ci jkl(y)∂uε k ∂xl. (16) Using the asymptotic expansion, Equation 12, and the stress field, Equation 16, yields σε i j(x)=σ(0) i j(x,y)+εu(1) i(x,y)+ε2u(2) i(x,y)+..., (17) where σ(m) i j(x,y)=Ci jkl(y)∂u(m) k ∂xl+Ci jkl(y)∂u(m+1) k ∂yl,...m=0,1,2,... (18) Now, using Equations 12 and 17 into Equation 15, and rearranging and grouping like exponents ε, one has for ε−1andε0, respectively ∂σ(0) i j ∂yj=0 (19) and ∂σ(0) i j ∂xj+∂σ(1) i j ∂yj=0. (20) 9 Using Equation 18 for m=0 and 19 one has ∂ ∂yjCi jkl∂u(1) k ∂yl=−∂Ci jkl ∂yj∂u(0) k ∂xl(21) where the solution is represented in the form of u(1) i(x,y)=χjk i(y)∂u(0) j(x) ∂xk, (22) and obtained from the equilibrium state given by Z YCi jpq∂χkl p ∂yq∂vi ∂yjdY=−Z YCi jkl∂vi ∂yjdY, (23) where vis a virtual displacement field and χis a characteristic periodic displacement field. Assuming m=0 in Equation 18, and considering Equation 22, one obtains σ(0) i j=Ci jkl+Ci jpq∂χkl p ∂yq∂u(0) k ∂xl. (24) The macroscopic equilibrium is obtained by averaging Equation 20 as ∂D σ(0) i jE Y ∂xj=0 (25) where D σ(0) i jE Y=1 |Y|Z Yσ(0) i jdY=CH i jkl∂u(0) k ∂xl(26) with CH i jkl=* Ci jkl+Ci jpq∂χkl p ∂yq+ Y, (27) which is the homogenised elasticity tensor of the media. The cell average operator is defined as ⟨•⟩Y=1 |Y|Z Y•dY. (28) Equations 23 and 27 are su fficient to obtain the complete elasticity tensor of a periodic media as a function of its microstructure. Several analytically explicit relations were derived in the last years 10 (see [46, 47, 49, 50]). However, in all cases, only simple geometries of the RVE are considered. If a complex geometry is considered, the explicit relations for the elasticity tensor of the media might become tedious to obtain, and thus, numerical methods are typically used in this kind of problem. In this study, both the equilibrium in the microscopic scale and the relation for the e ffective prop- erties are solved by the FE method. Equation 23 can be written in a vector form as Z Y ∂v ∂y!T C∂χ ∂ydY=−Z Y ∂v ∂y!T CdY (29) Dividing the domain into finite elements and assuming a general isoparametric interpolation, one can write the derivatives of displacement fields as ∂v ∂y=Bˆv,and∂χ ∂y=Bˆχ, (30) where ˆvand ˆχare the nodal virtual and characteristic displacement vectors, respectively, and Bis the global strain-displacement relation. Also, ˆχis not a function of the local coordinates, and the equilibrium becomes Z Y(Bˆv)TCBdYˆχ=−Z Y(Bˆv)TCdY, (31) and since the virtual displacement vector is arbitrary, one can choose it in such a manner that Z YBTCBdYˆχ=−Z YBTCdY. (32) Equation 32 can be written in a standard finite element form as Kˆχ=−f. (33) One can deduce by dimensional analysis that 6 equilibrium problems are depicted in Equation 33 and each one is denoted by a load case kl. Using index notation, one gets the global equilibrium problem as Ki jˆχkl j=−fkl i. (34) 11 For a single element e, and assuming a quadrature scheme to evaluate the integrals, the local stiffness matrix, and the local load vector can be written, respectively, as K(e) i j=npX m=1h C(e) i jpqB(e)(m) pB(e)(m) q W(m)J(m)i (35) f(e)kl i=npX m=1h C(e) i jklB(e)(m) j W(m)J(m)i , (36) where npis the total number of Gauss points used for the numerical integration, and W(m)andJ(m) are, respectively, the quadrature weight, and determinant of the Jacobian matrix associated with the Gauss point m. The global arrays are assembled in the usual way. A similar procedure is used for the discretisation of the integral Equation 27, which can be written as CH=1 |Y|Z YC+CBˆχdY. (37) The integral equation 37 can be substituted by a sum of all elements of the FE mesh and solved by a quadrature scheme, as CH i jkl=neX e=1npX m=1 C(e) i jkl+C(e) i jpqB(e)(m) pˆχkl(e) q W(m)J(m), (38) where neis the total number of elements of the mesh. The complete homogenised elasticity tensor of the media is obtained by Equation 38 subjected to the equilibrium problem of Equation 34. 3.2. Optimisation formulation The domain Ωis divided in a FE mesh, with element subdomains donated as Ω(e), such that Ω =Ne] e=1Ω(e)={x=(x1,x2,x3)∈R3: 0<xi<Li,i=1,2,3}, (39) where Neis the number of elements on the mesh,Uis the assembly operator, and Liare positive numbers. Furthermore, the mechanical properties of each element of the mesh are described by the geomet- rical parameters of a hollow sphere as 12 C(e)=C(e) ρ(e) o,ρ(e) i , (40) whereρ(e) oandρ(e) iare the outer, and inner diameters, respectively, of a hollow sphere associated to element e. The vectors containing the design variables are written as ρo=h ρ(1) oρ(2) oρ(3) o... ρ(Ne) oiT, (41) and ρi=h ρ(1) iρ(2) iρ(3) i... ρ(Ne) iiT, (42) and for simplicity, the following notation is adopted: ρ=ρoρi. (43) A compliance minimisation problem is considered constrained by a known maximum allowed relative density of the structure. The optimisation takes the following form: min ρf(ρ)=FTU(ρ) Subject to V(ρ)≤¯V,(44) where fis the compliance of the structure, Fis the global external load vector, Uis the global displacement vector, Vrelative density of the structure, and ¯Vis a fixed known relative density (upper bound). Additionally, side constraints are applied to the design variables, as the minimum and maximum admissible values for the diameters, as ρmin o≤ρ(e) o≤ρmax o ρmin i≤ρ(e) i≤ρmax i e=1,..., Ne.(45) Linear elastic behaviour is considered, and thus the FE equilibrium can be written as K(ρ)U(ρ)=F, (46) 13 where Kis the global sti ffness matrix. The optimisation problem can also be stated as an unconstrained problem, as min ρL(ρ) (47) where L(ρ)=f(ρ)+λ V(ρ)−¯V is the Lagrangian function and λis the Kuhn-Tucker multiplier. The minimum is attained when ∇ρL=0 (48) such that ∇ρf(ρ)+λ∇ρV(ρ)=0, (49) which is rewritten as −∂f(ρ) ∂ρq λ∂V(ρ) ∂ρq=1, (50) whereρqis the q-th design variable of the problem. Equation 50 is also known as the optimum condition and states that in the optimum (minimum) point ρ∗, all elements which are not bounded by the side constraints should attain a constant unitary value. Motivated by such behaviour, Bendsøe and Sigmund [51] proposed the use of an optimum criteria method in which the update of the design variables at the k-th iteration can be done as ρk+1 q=ρk qβη q (51) whereηis a relaxation parameter, and βqis the update parameter, given by βq=−∂f(ρ) ∂ρq λ∂V(ρ) ∂ρq, (52) whereλis the Kuhn Tucker Multiplier, obtained by a simple bisection algorithm written to satisfy the volume constraint. For simplicity, hereinafter the dependency of the arrays on the design variables is suppressed. The derivatives of the compliance, with respect to the q-th design variable, are given by 14 ∂f ∂ρq=−UT∂K ∂ρqU (53) and the local sti ffness matrix is K(e)=Z Ω(e)BT(e)C(e)B(e)dΩ(e), (54) where B(e)andC(e)are, respectively, the strain-displacement matrix and the elasticity tensor of the e-th element of the mesh. The derivatives of the local sti ffness matrix are obtained through ∂K(e) ∂ρq=Z Ω(e)BT(e)∂C(e) ∂ρqB(e)dΩ(e), (55) and the derivatives of the global sti ffness matrix through ∂K ∂ρq=Ne] e=1∂K(e) ∂ρq=Ne] e=1Z Ω(e)BT(e)∂C(e) ∂ρqB(e)Ω(e). (56) As the constitutive tensor of a given element edepends only on the design variables associated with that element, it is possible to state that ∂C(e) ∂ρq=δeq∂C(e) ∂ρq, (57) whereδeqis the Kronecker delta. The derivatives of the compliance with respect to the q-th element of the mesh can thus be written in terms of the local arrays as ∂f ∂ρq=uT(q)"Z Ω(q)BT(q)∂C(q) ∂ρqB(q)dΩ(q)# u(q), (58) where u(q)is the displacement vector of the q-th element of the mesh. Similarly, the relative density v(e)of an element e, is also dependent on the internal and external diameters of the hollow sphere, such that v(e)=v(e) ρ(e) o,ρ(e) i , (59) and the global relative density is given by 15 V=NeX e=1v(e), (60) with derivatives given by ∂V ∂ρq=NeX e=1∂v(e) ∂ρq. (61) As the derivatives of a given element depend only on the design variables associated with it, one can write ∂V ∂ρq=NeX e=1δeq∂v(e) ∂ρq=∂v(q) ∂ρq. (62) Equations 58 and 62 are valid for both outer and inner diameters associated to the q−th element of the mesh. Since the derivatives are local information, a moving limit scheme is adopted. For both external and internal diameters, the moving limits constraints can be set as ρk+1 q∈h ρk q−dq,ρk q+dqi (63) where dqis a positive moving limit of the q-th design variable. The maximum and minimum values for the moving limits are chosen as a percentage for the admissible range of the design variables. To update these limits, three consecutive iterations are considered. If the considered design variables vary monotonically in three successive iterations, the associated moving limit is relaxed. Otherwise, the moving limit for the considered design variable is tightened. It is reasonable to assume that for the external diameter of a hollow sphere, the derivatives of the compliance and the derivatives of the volume are given, respectively, by ∂f ∂ρ(e) o≥0 and∂V ∂ρ(e) o≤0,∀Ω(e)∈Ω (64) and, for the internal diameter, by ∂f ∂ρ(e) i≤0 and∂V ∂ρ(e) i≥0,∀Ω(e)∈Ω, (65) which ensures that the value of the update parameter βq(Equation 52) is always positive. 16 The behaviour of the external diameter of a hollow sphere and the pseudo densities at the classical SIMP case [51] are equal, since that, in both cases, when the value of the design variable increases, the volume increases, and the compliance reduces. Hence, analogously to the update of the design variables to the SIMP case, the external diameter is updated with a positive value of the relaxation coefficient, and it is adopted that ηo=0.5. (66) Conversely, the behaviour of the internal diameter is opposite to the behaviour of the external diameter. Therefore, it is adopted a relaxation coe fficient of ηi=−0.5, (67) such that, in both cases, the update of the variables is consistent with the optimal criteria method. To reduce the complexity of the structure, a filtering method is used. The chosen method is the basic density filtering [52], which consists in relating the variables attached to an element e, to the variables of its neighbourhood, which is defined by N(e)=n i,∥c(i)−c(e)∥≤Ro , (68) where c(i)is the centroid of the element i,c(e)is the centroid of the element eandRis a predefined filtering radius. The dependency of the design variable of the element eon its neighbours can be written as ˜ρ(e)=˜ρ(e)(ρi∈N(e))=P i∈N(e)w c(i) v(i)ρ(i) P i∈N(e)wc(i)v(i), (69) where w c(i) is a weighting function [53] chosen as a linearly decaying one. The optimisation is solved as a nested problem and it is considered that the range for the inner diameter is dependent on the value of the outer diameter for every iteration. Thus, the outer diameter is updated first, and the inner diameter is updated next as a function of the outer diameter on the current iteration. Figure 3 shows a flowchart for the nested optimisation approach. The convergence is assumed to be reached when the compliance and both outer and inner diameters have variation inferior to 1% in 5 subsequent steps. 17 Data Reading Mesh GenerationStart Define k = 1;ρo(1),ρi(1) Convergence?No Yes Return f,ρo,ρiSolve KU=F K=K(ρo(k),ρi(k)) Optimum Criterion: Findρo(k+1)k = k+1Solve KU=F K=K(ρo(k+1),ρi(k))Outer Diameter UpdateInner Diameter UpdateNESTED OPTIMISATION APPROACH Optimum Criterion: Findρi(k+1)Find df/dρoDensity Filtering Find df/dρ~iDensity Filtering ~Figure 3: Flowchart of the optimisation framework depicting the update of outer and inner diameters. 4. Results and discussion 4.1. E ffective properties The homogenised elasticity tensor of the hollow sphere material is obtained by Equation 38 sub- jected to the periodic solution of the equilibrium problem given in Equation 34. Some considerations are made to simplify the solution. Two di fferent isotropic materials are considered as base materials in the analyses. As the distri- bution of the material within the RVE is complex, a regular cubic RVE is adopted. In the regions containing material, the properties of solid material are assigned, representing an aluminium alloy (Young’s modulus, E=70 GPa and Poisson’s ratio, ν=0.3), whilst in regions with no material, prop- erties of void elements are assigned ( E=1.0−9GPa andν=0.3 ). In this way, void elements do not need to be removed from the FE mesh, hence the connectivity 18 of the elements does not need to be updated when the material distribution changes. Additionally, the application of periodic boundary conditions does not change as well. Conversely, this approach is computationally ine fficient since void elements contribute to the assembly of global arrays, even though it does not influence the homogenised elasticity tensor of the media. Additionally, a regular FE mesh is considered in the simulations, in which all elements on the mesh are cubes of the same size. The implication is that the hollow sphere itself can not be faithfully represented by a coarse mesh. However, this also renders the mesh generation and application of periodic boundary conditions simpler. The normalised dimensions liof the RVE, and number of elements niin the direction yiare shown in Table 1. Table 1: Normalised dimension liand number of finite elements nialong direction yifor the determination of e ffective properties. y1 y2 y3 l1 n1 l2 n2 l3 n3 1.0100 1.0100 1.0100 Figure 4(a) shows the FE model used to discretize the RVE, in which both void and solid ele- ments are portrayed. In Figures 4(b) and 4(c), the complete FE mesh and a cross-section are shown, respectively, with void elements suppressed in the visualisation. The element used is the Trilinear Isoparametric Hexahedral element with incompatible modes. The FE mesh is chosen to be both convergence-independent and able to faithfully represent the domain. (a) (b) (c) Figure 4: The FE mesh for the hollow sphere model: (a) complete FE mesh; FE mesh with void elements suppressed for (b) the entire RVE and (c) a cross-section of the RVE. For this model, two adjacent spheres are connected by a flat area delimited by the hollow sphere and the boundaries of the RVE. The thickness of the hollow sphere is constant and it is a function of 19 the geometrical parameters ρo, andρi. By varying the geometrical parameters in a given range, one can obtain a discrete set of the e ffective properties of the media. The considered normalised range for the geometrical parameters of the hollow sphere is 1.09≤ρo≤ 1.2 ρo−0.2≤ρi≤ρ0−0.01.(70) The ranges are chosen so that the geometry of the sphere is maintained. If a large outer diameter is used, the sphere tends to become a cube, due to the area used to connect two adjacent spheres. On the other hand, if a small outer diameter is used, adjacent spheres have small contact areas, which lowers the overall sti ffness. Additionally, the minimum inner diameter is chosen so that the finite element mesh can represent the thickness of the hollow sphere. For all combinations of parameters, the homogenised elasticity tensor has the form of CH=CH 1111CH 1122CH 11330 0 0 CH 2222CH 22330 0 0 CH 33330 0 0 CH 12120 0 CH 23230 sym. CH 1313, (71) in which CH 1111=CH 2222=CH 3333, CH 1122=CH 1133=CH 2233, CH 1212=CH 2323=CH 1313,(72) indicating a cubic symmetry on the material. Thus, the homogenised elasticity tensor is described by 3 independent constants. For the optimisation procedure, the homogenised elasticity tensor has to be written as a continuous and di fferentiable function, in the following form CH=CH(ρo,ρi). (73) 20 Thus, the least square method is used to adjust a continuous and di fferentiable surface based on the data given by the simulations. The Least Square Function is given by Φ =vtnsX n=1 CH i jkl ρ(n) o,ρ(n) i −C(n) i jkl2, (74) where nsis the number of simulations, ρn oandρn iare, respectively, the external and internal diam- eters of the hollow sphere at the n-th simulation and Cn i jklis the value of the component i jklof the homogenised elasticity tensor at the n-th simulation. A polynomial function is assigned to each independent coe fficient and the optimum function is found when ∇aΦ =0, (75) where ais the vector containing the coe fficients of the function used to approximate the e ffective properties. For all components of the homogenised elasticity tensor, the ideal polynomial function is ob- tained when a Coe fficient of Determination R2≥0.999 is achieved. The continuous functions for the independent components of the homogenised elasticity tensor are given by CH 1111(ρi,ρo)= +216.82ρo3−20.83ρiρo2−348.38ρo2 −588.11ρi2ρo+847.58ρiρo−183.31ρo +326.72ρi3−266.65ρi2−92.39ρi +108.94[GPa ],with R2=0.9998,(76) CH 1122(ρi,ρo)= +1148.55ρiρo2−671.20ρo3+282.47ρi3 +280.89ρi2+263.22ρi−833.54ρiρo +1337.64ρo2−981.53ρi2ρo−1191.81ρo +366.79[GPa ]with R2=0.9995,(77) and 21 CH 1212(ρi,ρo)= −8.13ρo3−2.44ρiρo2+80.76ρo2 +22.56ρi2ρo−117.16ρiρo−26.55ρo +23.46ρi3−81.61ρi2+157.47ρi −48.40[GPa ],with R2=0.9999.(78) Furthermore, the optimisation is constrained by the relative density of the structure. Thus, the Least Square Method is used to adjust a continuous function for the relative density of a hollow sphere as a function of the geometric parameters and is given by v(ρi,ρo)= +0.80ρo2+0.73ρiρo−0.67ρo−1.68ρi2 +1.00ρi−0.19,with R2=0.9998.(79) Equations 76-79 are valid for the ranges defined in Equation 70, and are used in the optimisation problem. Recent works on AHM show a good agreement with experimental data as well as a better perfor- mance compared to other homogenisation methods (see [54, 55, 56, 57, 58]). It is important to note that the homogenisation method demands that the size of the RVE is much smaller than the size of the macroscopic domain. As a limit theory, one expects to achieve better results when the relation of Equation 11 is satisfied. In addition, the approach used to find e ffective properties can be extrapolated to any kind of het- erogeneous material. Thus, it is possible to derive a set of equations to fit manufacturing constraints, such as the resolution of the domain, the materials used in the manufacturing process, and the geom- etry of the RVE. The two main considerations made to simplify the AHM analyses have a couple of drawbacks. The first one considers void elements with degenerated properties in regions with no material. The second one is that all elements on the mesh are regular, and thus the geometry of the hollow sphere is depicted as a pixel-like structure. Both considerations significantly a ffect the computational e ffort of the analyses. V oid elements contribute to the assembly of the global arrays in the finite element equi- librium and a regular mesh requires a fine discretisation to truly represent the geometry, which is not 22 possible with a coarse mesh. It is noteworthy to mention that the dimensions of a three-dimensional finite element equilibrium problem increase in the third power with the number of elements in the mesh. This leads to an expensive computational e ffort and makes the solution to these problems the most computationally demanding part of the algorithm. On the other hand, for complex geometries, such as the ones considered, the mesh generation and the application of periodic boundary condi- tions are the most challenging parts to implement. Consequently, the considerations to simplify the procedure are used even if they significantly a ffect the processing time. An alternative to avoid this high computational e ffort is either using external software for mesh generation or algorithms implemented alongside commercial software (see [43, 59]). However, for both AHM and optimisation procedures, it is important to have full control and access to all variables in the finite element algorithms, mainly because the local arrays have to be manipulated, which is not always possible with external commercial software. 4.2. Optimisation results The optimum distribution of hollow spheres to minimise structural compliance is considered and the problem is constrained by a maximum admissible relative density of the structure. Four cases are chosen to investigate the proposed formulation. Distinct relative densities, boundary conditions, and geometrical parameters are considered and discussed. The nomenclature and definitions for each case are shown in Table 2. The design variables associated with each element on the mesh are the outer and inner diameters of a hollow sphere. The update of the design variables is made through Equation 51. The mechanical properties and the relative density, as functions of the design variables, are given in Equations 76-79, and their derivatives are obtained explicitly. The filtering radius, for all cases, comprises 4 adjacent elements on the mesh. No mesh depen- dency was observed for the discretisation used in all cases. The visualisation of the results is made using Gmsh software [60]. Table 2: Nomenclature and definitions for the optimisation cases. The normalised dimension Liand the number of finite elements Niare used in direction xi, domain and load definitions, and admissible relative density. a>0.0. Casex1 x2 x3DomainRelative L1N1L2N2L3N3 Density [%] cb1 10a200 a 20 a20 Figure 2(a) 5 & 10 cb2 8a160 5a100 a10 Figure 2(a) 5 3pb110a200 a 20 a20 Figure 2(b) 5 & 10 3pb28a160 5a100 a10 Figure 2(c) 5 23 A convergence analysis of the objective function is shown in Figure 5 for the cb 1 optimisation case. The compliance is normalised by its value in the first iteration. In this case, the compliance is calculated at the end of each global step of the optimisation procedure, meaning that both design variables are updated to obtain the value of the objective function. All optimisation cases follow the same convergence behaviour, therefore only one plot is presented as it is representative of all cases. Furthermore, the convergence is monotonic, which shows the consistency of the implemented algorithm. Iteration0 10 20 30 40 50Normalised Compliance 0.30.40.50.60.70.80.91 Figure 5: Convergence plot for the cb 1 case. For the visualisation of the optimal distribution of hollow spheres, both design variables are pre- sented separately and represented in a linear grayscale. Black regions represent the maximum admis- sible value for a design variable, whilst light grey regions represent the minimum allowed value. In addition, as the relative density is a combination of outer and inner diameters, it is easier to analyse the results with the volume fraction alone instead of analysing both external and internal diameters separately. The optimum distribution of the relative density for each element of the mesh is also shown in the same linear grayscale. The e fficiency of the proposed multiscale approach is made by quantitative and qualitative analy- ses. Initially, the distribution of design variables and relative densities are compared for the same case using di fferent admissible global relative densities for the structures. Thus, one can verify that the patterns of material distribution within the domain are the same, regardless of the admissible relative density. As a quantitative comparison, the compliance of the same structure, now with a homoge- neous distribution of hollow spheres with the same global relative density is presented. With this 24 comparison, it is possible to measure the reduction of the compliance for optimised structures with the structure made of a homogeneous distribution of hollow spheres. As a qualitative comparison, the classical 0−1 SIMP [51] approach is presented for the same cases. This comparison is made to verify that the patterns of material distribution are the same for both cases. As the same physical problem is addressed for both approaches, it is expected that the regions with hollow spheres with a higher relative density, consequently a higher sti ffness, are the same that those regions with material in the classical SIMP approach. Figures 6 and 7 show the distribution of design variables and relative densities for cases cb 1 and 3pb 1, respectively, for admissible global relative densities of 5% and 10%. The distribution of design variables is consistent with a maximum allowed relative density of 5% and 10% since the same patterns are observed. This indicates consistency in the results as the patterns of material distribution should not change regardless of the maximum assumed relative den- sity. The di fference between the results obtained for the di fferent admissible relative densities is the amount of material in the same regions of the domain. Also, the results are compatible with the beam theory, as the distribution of material, for both cases, occurs in the regions where the highest stresses are expected, i.e., in regions with the highest bending moments and distant from the neutral axis. x1x2 x3 (a) x1x2 x3 (b) Figure 6: Optimum distribution of hollow spheres for the case cb 1. External diameter, internal diameter and relative density for cases with (a) ¯V=5%, and (b) ¯V=10%. A linear scale is assumed between the minimum admissible value (light grey) and the maximum admissible value (black) for the considered variables. Table 3 shows quantitative comparisons for all cases, where the optimised results are compared to the compliance obtained by a homogeneous distribution of hollow spheres defined by ρh oandρh i. For all cases, a significant reduction in compliance is observed, with values ranging from 42% to 79%. For the cases cb 1 and 3pb 1, relative densities of 5% and 10% are considered. In those cases, the higher the admissible relative density, the higher the compliance reduction. Since all elements have the same relative density in the case of the homogeneous distribution, regions that barely a ffect 25 x1x2 x3 (a) x1x2 x3 (b) Figure 7: Optimum distribution of hollow spheres for the case 3pb 1. External diameter, internal diameter and relative density for cases with (a) ¯V=5%, and (b) ¯V=10%. A linear scale is considered between the minimum admissible value (light grey) and the maximum admissible value (black) for the considered variables. the compliance still have a high relative density. On the other hand, for the optimised results, regions that are not significant to the compliance are occupied by elements with a low relative density. This behaviour is intensified when higher admissible relative densities are considered. Moreover, the higher the outer diameter, the lower the compliance minimisation. This is expected since the flat area connecting two adjacent spheres is a function of the outer sphere, and for a higher outer diameter, this region is larger, thus there is more material connecting the hollow spheres, hence increasing the overall sti ffness of that particular region. This behaviour is observed in all cases. Figures 8 - 11 show, for all cases, a qualitative comparison for the optimum distribution of mate- rial using SIMP approach. The results show the optimum distribution of maximum allowed relative density, that is, the elements that contribute more to the minimisation of the compliance remain in the optimised topology towards minimising their compliance. In a cantilever beam, when bending dominates the problem, some considerations concerning the stiffness of the structure can be made. The bending moment is the highest at the cantilever and it has a null value on the opposite side, where the load is applied. Sections in which the bending moment is higher are more susceptible to higher displacements, rendering the beam in these sections to be more flexible. Considering the cross-section of a beam, the moment of inertia increases when more material is added to the regions away from the neutral axis of the section. The higher the moment of inertia of a section, the lesser its flexibility. Adding material to regions close to the neutral axis almost does not modify the moment of inertia of a given section. One can see that, for both cases, the optimisation procedure leads to topologies in which there is 26 Table 3: Compliance reduction for optimised and homogeneous distributions of hollow spheres of same relative density (¯V) for all cases. Homogeneous distributions are defined by ρh oandρh i. Caseρh oρh i¯V[%]Compliance Reduction [%] cb11.2 1.1767 545.56 1.15 1.1252 58.17 1.10 1.0730 70.61 cb11.2 1.1520 1062.00 1.15 1.0987 67.59 1.10 1.0443 73.55 cb21.2 1.1767 554.90 1.15 1.1252 64.56 1.10 1.0730 74.25 3pb11.2 1.1767 542.83 1.15 1.1252 55.94 1.10 1.0730 68.09 3pb11.2 1.1520 1057.30 1.15 1.0987 63.48 1.10 1.0443 70.12 3pb21.2 1.1767 564.24 1.15 1.1252 71.88 1.10 1.0730 79.19 (a) (b) x1x2 x3 (c) Figure 8: Qualitative comparison for case cb 1 between the optimised distribution of relative densities for (a) ¯V=5%, (b)¯V=10%, and (c) the optimum distribution of material using SIMP. a concentration of material close to the cantilever, alongside material in the regions away from the neutral axis. These facts show a consistency of the results with the physical aspect of the problem. The structure obtained with SIMP follows the same principles of hollow sphere structures. It can be seen that material is added to regions closer to the cantilever and, at the cross-section, in the regions far from the neutral axis. The main di fference between the two approaches is that, for the SIMP case, 27 x1x2 x3(a) x1x2 x3 (b) Figure 9: Qualitative comparison for case cb 2 between the optimised distribution of relative densities for cases (a) ¯V=5% and (b) the optimum distribution of material using SIMP. x1x2 x3 (a) x1x2 x3 (b) x1x2 x3 (c) Figure 10: Qualitative comparison for case 3pb 1 between the optimised distribution of relative densities for (a) ¯V=5%, (b)¯V=10%, and (c) the optimal distribution of material using SIMP. x1x2 x3 (a) x1x2 x3 (b) Figure 11: Qualitative comparison for case 3pb 2 between the optimised distribution of relative densities for (a) ¯V=5% and (b) the optimum distribution of material using SIMP. there are solid connections between the upper and lower regions of the beam since void regions do not provide sti ffness to the structure. Here, it is not necessary to connect the two regions with spheres with a high relative density since even spheres with low relative density can provide sti ffness to the structure. With this consideration, it is clear that the algorithm is concerned with adding material to the most significant regions of the domain. On the other hand, for the 3-point bending cases, the highest bending moment is at the central 28 section and the lowest is near the supports. Concerning the moment of inertia, the same considerations of the cantilever beam cases apply. Thus, the optimisation procedure leads to a material distribution in which there is more material close to the central part of the beam. Also, the hollow spheres with higher relative density, consequently higher sti ffness, are concentrated at the upper and lower parts of the beam, which complies with the beam theory for a beam under 3–point bending load. As the hollow spheres with lower relative densities also add sti ffness to the structure, they can be used to join the upper and lower parts of the beam, letting the material with higher sti ffness be placed only at places with higher bending moment and at regions in which the moment of inertia is increased, thus reducing the overall relative density of the optimised structure. The comparison of the proposed problem with the classical SIMP approach shows consistency. In both cases, an arc is formed to provide higher sti ffness for the structures. It is important to highlight that, in the SIMP case, the arc must be joined with the straight section, and radial reinforcements appear, otherwise, the structure would not present a proper sti ffness. For the hollow sphere case, the material in the interior of the arc still provides sti ffness to the structure, rendering it not necessary to connect the arc to the region where the load is applied with hollow spheres with high relative density. The reinforcement that appears in the hollow sphere case allows the load to be distributed homogeneously in the regions inside the arc, allowing material with a low relative density to be used in this region. 5. Manufacturing of hollow spheres Some hollow sphere structures with di fferent dimensions were manufactured to demonstrate the feasibility of constructing the models and the possibility of manufacturing the spheres when a range of different geometrical features is considered, and also by varying the dimensions of the structure. The topologies of the hollow spheres are obtained directly from the FE mesh of the homogenisation procedure (see Figure 4). The FE mesh contains information about the nodes and connectivities of every node on the mesh. As most commercial additive manufacturing (AM) processes use STL file format as default, an in-house code was developed to convert the FE data into an STL file. The void elements are removed during this process. The smoothing of the geometries is carried out on MeshLab [61] using the Laplacian Smooth Algorithm [62]. All spheres were produced by AM by using a Form3 3D printer from FormLabs, which utilises the stereolithography (SLA) principle. Here, Clear®resin, a photo-polymerising resin has been used, 29 which polishes to near optical transparency, ideal for showcasing internal features. The parts are printer layer-by-layer with a layer thickness of 100 µm. After the printing is finished, the parts are washed in a Form Wash chamber with isopropyl alcohol (IPA) bath for 15 min, followed by curing in a Form Cure oven at 65◦C for 60 min, in which the parts sit on a rotating turntable. In Figure 12, the design variables of the hollow spheres are varied, where both non-smooth and smooth geometries are shown. Figures 12(a)-(b) show a homogeneous distribution of the spheres. It can be seen that the smoothing of the geometry can maintain the geometry of the structure with high accuracy. In Figures 12(c)-(d), a variation on the internal diameters is considered. Four di fferent internal diameters are considered, with ρ(1) i>ρ(2) i>ρ(3) i>ρ(4) i, showing that this type of structure can be manufactured with thin walls. Figure 12(e)-(f) shows the variation of the external diameters of the hollow spheres. It can be seen that the higher the diameters considered, the higher the contact area between two adjacent hollow spheres. In Figure 13, samples with di fferent sizes and numbers of hollow spheres are shown. Structures with 16, 512, and 128 spheres were printed, which demonstrates the feasibility of using this type of material. Here it is demonstrated that it is possible to manufacture this type of structure strictly following the proposed design variables, that is, it is feasible to produce samples by varying both internal and external radii of the hollow spheres. Moreover, it is demonstrated that the number and size of the spheres can be varied in the manufacturing process without losing quality and resolution and that it is feasible to obtain structures comprised of small hollow spheres, which depends only on the resolution of the used equipment. In this way, it is possible to manufacture samples with an optimised distribution of hollow spheres without violating the conditions imposed by the AHM to determine the effective properties of the medium. 6. Conclusion In this work, a multiscale optimisation framework to find the distribution of hollow-sphere materi- als has been developed. The e ffective mechanical properties are obtained by the Asymptotic Homog- enization Method and written as a function of the geometrical parameters of a single hollow sphere. The topology optimisation procedure is an adaptation of the SIMP method. The objective function for all cases is compliance minimisation and the relative density of the structure is used as a constraint. A 30 (a) (b) (c) (d) (e) (f) Figure 12: Variation of geometrical parameters of the printed hollow sphere structures. Homogeneous distribution with a (a) non-smooth geometry and a (b) smooth geometry. Variation of the internal diameter for a (c) non-smooth geometry, and a (d) smooth geometry. Variation of the external diameter for a (e) non-smooth geometry, and a (f) smooth geometry. cantilever beam and a 3–point bending beam are considered as design domains, and relative densities of 5%, and 10% are considered. Two comparisons are performed to validate the proposed multiscale approach. As a quantitative comparison, the optimised structures are compared to structures formed by a homogeneous distribu- tion of hollow spheres with the same relative density. Compliance reductions between 42% and 79% are reached. Qualitatively, the optimised structures are compared to an optimal distribution of mate- rial using the SIMP approach. In both cases, the distribution of the material in the domain follows the same pattern, validating the originally proposed framework. In the regions where there is material for the SIMP material model, the developed algorithm adds hollow spheres with high volume fraction. In contrast, for regions with no material in the SIMP material model, the algorithm uses a material with 31 (a) (b) (c) Figure 13: Variation on the size and number of hollow spheres. The structures were printed with (a) 16, (b) 512, and (c) 128 spheres. low volume fraction. These comparisons show the e fficiency of the proposed multiscale approach to optimise the sti ffness of hollow sphere structures. The feasibility of the proposed approach is tested by additively manufacturing some hollow sphere structures. Hollow spheres are 3D printed by varying their design variables, as well as the number and size of hollow spheres, thus indicating that it is feasible to design, optimise, and manufacture such structures without losing quality and resolution. Acknowledgements HA is supported by the Royal Academy of Engineering under the Research Fellowship scheme [Grant No. RF /201920 /19/150]; ELC thanks CNPq (process number 303900 /2020-2) and FAPESC (process number 2021TR843); VT acknowledges CNPq (process No. 310656 /2018-4), CAPES-FCT (No. AUXPE 88881.467834 /2019-01 - Finance Code 001), and FAPESP-FAPERGS (process num- ber: 2019 /15179-2 and 19 /2551). 32