Marco Montemurro
email: marco.montemurro@ensam.eu

Khalil Refai

Anita Catapano

Thermal design of graded architected cellular materials through a CAD-compatible topology optimisation method

Keywords: Architected cellular materials, Additive manufacturing, Homogenisation, Topology optimisation, Heat conduction, NURBS hyper-surfaces

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Laser powder bed fusion (LPBF) allows producing efficient solutions for heat transfer applications. Components operating at high temperatures require efficient dissipation of thermal energy to prevent failure. Whilst the demand on cooling systems is increasing in the pursuit of component performance gains, the package size available for heat-sinks and tools is decreasing due to lightness-oriented requirements. To this end, architected cellular materials (ACMs) with a periodic micro-structure can be employed to satisfy, concurrently, lightness and heat transfer-related requirements, including heat dissipation, heat shielding, thermal insulation and flame arresting.

A review on the design of a particular class of ACMs, i.e., the well-known lattice materials, can be found in [START_REF] Pan | Design and Optimization of Lattice Structures: A Review[END_REF]. As far as the thermal behaviour of ACMs is concerned, it depends on several factors such as the representative volume element (RVE) topology, material and manufacturing process characteristics, as well as loading conditions. By modifying the RVE topology and the relative density at the microscopic/mesoscopic scale (i.e., the scale of the RVE), one can fulfil the desired requirements at the macroscopic one (i.e., the scale of the structure). Many works on the optimal design of ACMs focuses mainly on the design and assessment of their macroscopic elastic properties [START_REF] Hu | Two-scale concurrent topology optimization method of hierarchical structures with self-connected multiple lattice-material domains[END_REF][START_REF] Montemurro | A General Multi-Scale Topology Optimisation Method for Lightweight Lattice Structures Obtained through Additive Manufacturing Technology[END_REF][START_REF] Xu | An isogeometric approach to topology optimization of spatially graded hierarchical structures[END_REF][START_REF] Behrou | Topology optimization of nonlinear periodically microstructured materials for tailored homogenized constitutive properties[END_REF][START_REF] Borovinšek | Designing 2d auxetic structures using multi-objective topology optimization[END_REF]. Nevertheless, to the best of the authors' knowledge, only a limited number of investigations on the design of ACMs for heat transfer is available in the literature [START_REF] Wong | Selective laser melting of heat transfer devices[END_REF][START_REF] Wong | Convective heat transfer and pressure losses across novel heat sinks fabricated by selective laser melting[END_REF][START_REF] Wadley | Thermal applications of cellular lattice structures[END_REF][START_REF] Catchpole-Smith | Thermal conductivity of tpms lattice structures manufactured via laser powder bed fusion[END_REF][START_REF] Fink | Fluid dynamics of flow through microscale lattice structures formed from self-propagating photopolymer waveguides[END_REF][START_REF] Do | Additive manufacturing of interpenetrating periodic open cellular structures (interpocs) with in operando adjustable flow characteristics[END_REF][START_REF] Dong | A survey of modeling of lattice structures fabricated by additive manufacturing[END_REF][START_REF] Han | A novel design method for nonuniform lattice structures based on topology optimization[END_REF][START_REF] Maloney | Multifunctional heat exchangers derived from three-dimensional micro-lattice structures[END_REF][START_REF] Lu | Active cooling by metallic sandwich structures with periodic cores[END_REF][START_REF] Valdevit | Optimal active cooling performance of metallic sandwich panels with prismatic cores[END_REF][START_REF] Han | Optimization design of nonuniform cellular structures for additive manufacturing[END_REF][START_REF] Cheng | The equivalent thermal conductivity of lattice core sandwich structure: A predictive model[END_REF], and most of them are of experimental nature [START_REF] Wong | Selective laser melting of heat transfer devices[END_REF][START_REF] Wong | Convective heat transfer and pressure losses across novel heat sinks fabricated by selective laser melting[END_REF][START_REF] Wadley | Thermal applications of cellular lattice structures[END_REF][START_REF] Catchpole-Smith | Thermal conductivity of tpms lattice structures manufactured via laser powder bed fusion[END_REF][START_REF] Fink | Fluid dynamics of flow through microscale lattice structures formed from self-propagating photopolymer waveguides[END_REF][START_REF] Do | Additive manufacturing of interpenetrating periodic open cellular structures (interpocs) with in operando adjustable flow characteristics[END_REF]. In particular, Wong et al. [START_REF] Wong | Selective laser melting of heat transfer devices[END_REF] used LPBF to manufacture pin-fin heat exchanger made of stainless steel (316L) and aluminium alloy (6061). The geometric shapes include circular and diamond cross-sections, as well as additional ACMs composed of elliptical cross-sections placed at the corners of the substrate. This research was extended in [START_REF] Wong | Convective heat transfer and pressure losses across novel heat sinks fabricated by selective laser melting[END_REF], using a strut-type lattice, which is based on a geometric model commonly used as a support structure for LPBF. Wadley et al. [START_REF] Wadley | Thermal applications of cellular lattice structures[END_REF] compared the effective thermal conductivity (ETC) of copper textile-based structures including random open-cell metal foam with the one of lattice materials. They found that the overall thermal efficiency of lattice materials is better than the one of the open-cell foams having the same relative density. Catchpole-Smith et al. [START_REF] Catchpole-Smith | Thermal conductivity of tpms lattice structures manufactured via laser powder bed fusion[END_REF] showed that the thermal conductivity of lattice materials mainly depends on material properties and relative density. They also showed that the Schwarz RVE topology always provide the highest ETC, while the one of diamond and gyroid RVEs topologies is slightly lower.

Nevertheless, performing experimental tests on ACMs to assess their ETC at the macroscopic scale is an expensive and time-consuming process. Moreover, as mentioned above, there are many factors influencing the macroscopic thermal behaviour of ACMs, thus experimental tests should be performed for each configuration of the RVE geometry. Therefore, before manufacturing and testing samples, it is very useful to develop a predictive multi-scale numerical model to assess the influence of each variable, characterising the RVE geometry of the AMC, on the macroscopic thermal behaviour of the structure and, subsequently, plan a reduced number of experimental tests for validation. A survey on modelling strategy for lattice materials obtained through various additive manufacturing technologies can be found in [START_REF] Dong | A survey of modeling of lattice structures fabricated by additive manufacturing[END_REF]. Generally speaking, in the framework of multi-scale modelling strategy, at the mesoscopic scale, the RVE of the ACM can be modelled, from a mechanical point of view, as a heterogeneous medium composed of two phases, i.e., the bulk material and the void. Conversely, at the macroscopic scale, it can be modelled as an equivalent homogeneous anisotropic continuum whose ETC is described by a set of equivalent material properties [START_REF] Montemurro | A General Multi-Scale Topology Optimisation Method for Lightweight Lattice Structures Obtained through Additive Manufacturing Technology[END_REF][START_REF] Delucia | Determination of the effective thermoelastic properties of cork-based agglomerates[END_REF][START_REF] Refai | Determination of the effective elastic properties of titanium lattice structures[END_REF][START_REF] Bertolino | Multi-scale shape optimisation of lattice structures : an evolutionary-based approach[END_REF].

Many analytical models for predicting the ETC of general composite materials have been proposed in the literature, each one being characterised by a different degree of complexity, such as those of Schapery [START_REF] Schapery | Thermal expansion coefficients of composite materials based on energy principles[END_REF], Chamberlain [START_REF] Chamberlain | Derivation of expansion coefficients for a fibre reinforced composites[END_REF], Chamis [START_REF] Chamis | Simplified composite micromechanics equations for hygral, thermal and mechanical properties[END_REF], Rosen and Hashin [START_REF] Rosen | Effective thermal expansion coefficients and specific heats of composite materials[END_REF]. However, all these analytical models have a common limitation: they do not consider the influence of local geometrical features of the RVE on its macroscopic behaviour. To go beyond the limitations of analytical models, numerical homogenisation techniques have been developed in the last five decades. Finite element (FE)-based models are very versatile, and, depending on the sophistication of the model, they can lead to a more accurate assessment of the equivalent macroscopic behaviour of the ACM [START_REF] Montemurro | A General Multi-Scale Topology Optimisation Method for Lightweight Lattice Structures Obtained through Additive Manufacturing Technology[END_REF][START_REF] Delucia | Determination of the effective thermoelastic properties of cork-based agglomerates[END_REF][START_REF] Refai | Determination of the effective elastic properties of titanium lattice structures[END_REF][START_REF] Bertolino | Multi-scale shape optimisation of lattice structures : an evolutionary-based approach[END_REF].

As far as the multi-scale design of uniform ACMs and graded ACMs (GACMs) is concerned, only few studies dealing with thermal applications are available in the literature [START_REF] Han | A novel design method for nonuniform lattice structures based on topology optimization[END_REF][START_REF] Maloney | Multifunctional heat exchangers derived from three-dimensional micro-lattice structures[END_REF][START_REF] Lu | Active cooling by metallic sandwich structures with periodic cores[END_REF][START_REF] Valdevit | Optimal active cooling performance of metallic sandwich panels with prismatic cores[END_REF][START_REF] Han | Optimization design of nonuniform cellular structures for additive manufacturing[END_REF][START_REF] Li | Design and optimization of graded cellular structures with triply periodic level surface-based topological shapes[END_REF][START_REF] Huang | Topology optimization of lattice support structures for heat conduction in selective laser melting[END_REF][START_REF] Cheng | Coupling lattice structure topology optimization with design-dependent feature evolution for additive manufactured heat conduction design[END_REF][START_REF] Cheng | Concurrent lattice infill with feature evolution optimization for additive manufactured heat conduction design[END_REF]. For instance, in [START_REF] Han | A novel design method for nonuniform lattice structures based on topology optimization[END_REF], a topology optimisation (TO) method based on the solid isotropic material with penalisation (SIMP) approach is proposed to determine the optimal topology of non-uniform lattice materials with a pre-defined RVE geometry. In this context, the pseudo-density field of the SIMP method, defined at the macroscopic scale, is mapped onto a non-uniform lattice material (at the lower scale). Of course, the mapping between the pseudo-density field at the macroscopic scale and the thickness of the struts of the lattice RVE at the lower scale is determined through a suitable FE-based homogenization scheme. A similar approach is used in [START_REF] Li | Design and optimization of graded cellular structures with triply periodic level surface-based topological shapes[END_REF] and extended to the case of GACMs whose RVE has a triply periodic minimal surface (TPMS) thin-walled topology. In particular, the approach proposed in [START_REF] Li | Design and optimization of graded cellular structures with triply periodic level surface-based topological shapes[END_REF] is based on two main features: 1) a TO algorithm based on a pseudo-density field and 2) a numerical homogenisation scheme to set the link between the macroscopic pseudo-density field and the geometrical feature of the TPMS topology of the RVE. Huang et al. [START_REF] Huang | Topology optimization of lattice support structures for heat conduction in selective laser melting[END_REF] make use of lattice materials to optimise the support structure required in selective laser melting process to support overhanging regions. The goal of the methodology proposed in [START_REF] Huang | Topology optimization of lattice support structures for heat conduction in selective laser melting[END_REF] is to maximise the macroscopic ETC of the support structure made of lattice material in order to improve heat dissipation and to reduce part distortion during fabrication. In [START_REF] Cheng | Coupling lattice structure topology optimization with design-dependent feature evolution for additive manufactured heat conduction design[END_REF][START_REF] Cheng | Concurrent lattice infill with feature evolution optimization for additive manufactured heat conduction design[END_REF], the TO method for lattice materials is enhanced with design-dependent movable features (over which design-dependent loading conditions are imposed) and applied to heat condution problems. The method-ology presented in [START_REF] Cheng | Coupling lattice structure topology optimization with design-dependent feature evolution for additive manufactured heat conduction design[END_REF][START_REF] Cheng | Concurrent lattice infill with feature evolution optimization for additive manufactured heat conduction design[END_REF] is based, on the one hand, on the asymptotic homogenisation to determine the equivalent macroscopic thermal conductivity of the lattice and, on the other hand, on the level-set method to obtain, at each iteration, an implicit geometric description of the topology boundary where thermal boundary conditions (BCs) are applied via an immersed boundary method.

It is noteworthy that most of the aforementioned works make use of deterministic algorithms to update the design variables defining the topological descriptor. A sound alternative, which does not require the computation of the gradient of the physical responses involved in the problem formulation, is the hybrid cellular automaton method introduced by Tovar et al. [START_REF] Tovar | Hybrid Cellular Automata: a biologically-inspired structural optimization technique[END_REF][START_REF] Tovar | Bone Structure Adaptation as a Cellular Automaton Optimization Process[END_REF]. This method makes use of a metaheuristic algorithm simulating the process of structural adaptation in bones and combines the cellular automaton paradigm with FE analysis. Hybrid cellular automaton method proved its effectiveness in dealing with a variety of TO problems [START_REF] Tovar | Topology Optimization Using a Hybrid Cellular Automaton Method With Local Control Rules[END_REF][START_REF] Tovar | Optimality Conditions of the Hybrid Cellular Automata for Structural Optimization[END_REF] including those dealing with the design of the equivalent macroscopic properties of ACMs [START_REF] Da | Design of materials using hybrid cellular automata[END_REF].

Taking inspiration from the aforementioned works dealing with the multi-scale design of ACMs, this study presents a computer-aided design (CAD)-compatible TO method for designing GACMs with an enhanced macroscopic thermal behaviour. The proposed approach makes use of: (1) a density-based TO method reformulated in the framework of the non-uniform rational basis spline (NURBS) hyper-surfaces [START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | NURBS Hypersurfaces for 3D Topology Optimisation Problems[END_REF] and (2) the strain energy-based homogenisation method (SEHM) of periodic media applied to thermal problems [START_REF] Montemurro | A General Multi-Scale Topology Optimisation Method for Lightweight Lattice Structures Obtained through Additive Manufacturing Technology[END_REF][START_REF] Delucia | Determination of the effective thermoelastic properties of cork-based agglomerates[END_REF][START_REF] Refai | Determination of the effective elastic properties of titanium lattice structures[END_REF][START_REF] Bertolino | Multi-scale shape optimisation of lattice structures : an evolutionary-based approach[END_REF] to determine the ETC of the RVE of the ACM. In particular, for a given RVE topology, the SEHM is used to derive the trend of the ETC at the macroscopic scale as a function of the RVE relative density. This function is then integrated in the NURBSdensity-based TO algorithm to penalise the thermal conductivity tensor of the continuum at the macroscopic scale. Unlike classical density-based TO approaches [START_REF] Bendsoe | Topology Optimization -Theory, Methods and Applications[END_REF], the NURBS-density-based TO method separates the pseudo-density field, describing the topology of the continuum, from the mesh of the FE model. More precisely, for general 3D problems, a 4D NURBS hypersurface is used as a topology descriptor, whilst for 2D problems a standard 3D NURBS surface is employed. In this way, the topological descriptor, i.e., the pseudo-density field, relies on a purely CAD-compatible geometric entity, i.e., the NURBS hyper-surface, and the optimisation variables are the pseudo-density and the associated weight at each control point (CP) . Regarding the application illustrated in this study, the thermal conductivity matrix of the generic element is penalised through a physically-based penalty function derived as result of the homogenisation process (for a given RVE topology). As discussed in [START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | NURBS Hypersurfaces for 3D Topology Optimisation Problems[END_REF], this approach is characterised by some advantages over conventional density-based TO methods. Firstly, since the topological descriptor consists in a high-level geometric parametrisation of the pseudo-density field in the form of a NURBS entity, the optimised topology does not depend upon the quality of the mesh of the finite element model. Secondly, unlike the classical SIMP approach, there is no need to define a further filter zone because the NURBS local support property establishes an implicit relationship among contiguous elements. Thirdly, when compared to the classical SIMP approach, the number of design variables is reduced (because the number of CPs is lower than the number of elements constituting the FE model). Fourthly, since the topology is described through a NURBS entity, the boundary of the topology is available at each iteration of the optimisation process, thus, the integration of constraints of geometric nature [START_REF] Costa | Maximum length scale requirement in a topology optimisation method based on NURBS hyper-surfaces[END_REF][START_REF] Costa | Minimum Length Scale Control in a NURBSbased SIMP Method[END_REF] (e.g., on the local curvature of the boundary, on the local direction of the tangent vector, maximum member size, minimum member size, etc.) in the problem formulation and the CAD reconstruction phase of the boundary of the optimised topology become easy tasks [START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Bertolino | An Efficient Hybrid Optimisation Strategy for Surface Reconstruction[END_REF]. Finally, some fundamental properties of the NURBS basis functions, like the local support property, can be conveniently exploited to derive the gradient of the physical responses with respect to the topological variables, i.e., pseudo-density at CPs and weights.

In this study, a wide campaign of numerical analyses is conducted on 31 RVE topologies of ACMs to determine, for each configuration, the equivalent thermal conductivity tensor at the macroscopic scale. Among them, seven topologies have been chosen and used in the framework of the NURBS-density-based TO method. The effectiveness of the proposed approach is tested on both 2D and 3D benchmark problems taken from the literature.

The paper is organised as follows: Sec. 2 introduces the problem description, the ACMs classes and the general work-flow of the design methodology proposed in this work. Sec. 3 briefly recalls the fundamental of the SEHM and the details of the numerical model of the RVE at the lower scale. The numerical results of the SEHM are presented and discussed in Sec. 4. Sec. 5 introduces the fundamentals of the NURBS-density-based TO algorithm used in this work as well as the numerical strategy to recover GACMs from the results of the TO process. Sec. 6 presents the results of the TO. Finally, Sec. 7 ends the paper with some concluding remarks and prospects.

Notation. Upper-case bold letters and symbols are used to indicate tensors and matrices, while lower-case bold letters and symbols indicate column vectors. Subscripts m and M denote quantities evaluated at RVE scale and macroscopic scale, respectively.

Problem description

Material and geometries of the architected cellular materials

The first part of this study aims to assess and compare the macroscopic ETC of various ACMs, which can be grouped in three main classes: the strut-based ACMs class, as shown in Fig. 1, collecting 17 RVE geometries, the thin-walled TPMS ACMs class including eight RVE geometries, as illustrated in Fig. 2(a), and the skeletal TPMS including six RVE geometries, as shown in Fig. 2(b). The RVE of thin-walled TPMS ACMs is characterised by a constant thickness and is generated by means of an implicit surface as follows:

cos(k m1 x m1 ) + cos(k m2 x m2 ) + cos(k m3 x m3 ) -t m = 0, (Schwarz), (1) cos
(k m1 x m1 ) sin(k m2 x m2 ) + cos(k m2 x m2 ) sin(k m3 x m3 )+ + cos(k m3 x m3 ) sin(k m1 m1) -t m = 0, (Gyroid), (2) 2
[cos(k m2 x m2 ) cos(k m3 x m3 ) + cos(k m1 x m1 ) cos(k m2 x m2 )+ + cos(k m1 x m1 ) cos(k m3 x m3 )] -(cos(2k m1 x m1 ) + cos(2k m2 x m2 )+ + cos(2k m3 x m3 )) -t m = 0, (IWP), (3) 3
[cos(k m1 x m1 ) + cos(k m2 x m2 ) + cos(k m3 x m3 )] + + 4 cos(k m1 x m1 ) cos(k m2 x m2 ) cos(k m3 x m3 ) -t m = 0, (Neovius). (4)
In the above formulae, x mi , (i = 1, 2, 3) are the Cartesian coordinate at the RVE scale, whilst the parameter t m is a threshold value (or isovalue) representing the offset from the zero level-set. Parameters k mi are defined as:

k mi := 2π n mi L mi , i = 1, 2, 3, (5)
where L mi := 2a mi and n mi represent the size of the RVE and number of half-waves along the x mi axis, respectively. The skeletal TPMS ACMs are generated by using boolean operations on RVE topologies belonging to the second class.

All the ACMs considered in this work are constituted of 6061 aluminium alloy whose thermal conductivity and density are K m = 170 W/mK and ρ m = 2700 kgm -3 . Wong et al. [START_REF] Wong | Selective laser melting of heat transfer devices[END_REF] used the same material to create multiple heat sinks with different geometries through the selective laser melting process. Regardless of its topology, each RVE is embedded in a parallelepiped whose overall volume is V RVE = 8a m1 a m2 a m3 , as illustrated in Fig. 3. In particular, in the following of the paper, a cubic embedding domain is considered for each RVE topology, i.e., 2a m1 = 2a m2 = 2a m3 = L m .

General work-flow of the design methodology

The main goal of this work is to propose a general design methodology for GACMs obtained through AM technology. The proposed method is based on the use of: a) the SEHM to determine the macroscopic behaviour of the RVE of the ACM; b) standard approximation techniques (i.e., curve fitting) to build the physically-based penalty schemes for the characteristic tensors involved in the TO problem at the macroscopic scale; c) NURBS hyper-surfaces as a descriptor of the continuum topology at the macroscopic scale; d) a simple reconstruction technique to pass from the optimised topology to the corresponding graded structure of the ACM at the lower scale.

The work-flow of the proposed approach is illustrated in Fig. 4. In particular, in this work, the proposed methodology is applied to a representative design case: the minimisation of the thermal compliance of a GACM subject to a constraint on its overall volume.

The design strategy is articulated in the following six steps:

1. For each RVE topology, the FE model is generated and its geometry and mesh are opportunely parametrised. In particular, the geometric parameters, which depend on the RVE topology, are those illustrated in Fig. 3: in the case of strut-based geometries, the element size is set as a fraction of the strut diameter, whilst in the case of TPMS the element size is chosen as a fraction of the wall thickness.

The number of elements composing the FE model is chosen as a result of a mesh sensitivity analysis. An example of mesh for octet and gyroid topologies is shown in Fig. 5 2. The SEHM (presented in Sec. 3) is applied to the RVE and the numerical results, obtained for different values of the relative density ρ M of the RVE, are collected in a database. 3. A standard curve fitting technique is applied to derive the physically-based penalisation scheme φ(ρ M ) to be used during TO at the macroscopic scale of the structure. 4. The penalisation function φ(ρ M e ) is affected to the thermal conductivity matrix of each element composing the FE model at the macroscopic scale. In particular, the pseudo-density field at this scale, i.e., ρ M , is described by means of a NURBS hyper-surface, as discussed in Sec. 5. 5. The TO problem is solved and the optimised solutions are available in a CADcompatible form as discussed in [START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | NURBS Hypersurfaces for 3D Topology Optimisation Problems[END_REF][START_REF] Bertolino | An Efficient Hybrid Optimisation Strategy for Surface Reconstruction[END_REF]]. 6. The optimised pseudo-density field at the macroscopic scale is mapped onto a GACM at the lower scale; in particular, an automatic reconstruction strategy is used to recover the parameters of the RVE corresponding to the local optimal value of ρ M .

It is noteworthy that, ACMs with an open-cell RVE can be easily manufactured through LPBF without adding support material because of the small RVE size. Moreover, the intermediate density values resulting from the TO problem can be mapped to a physicallyconsistent RVE. Of course, in the regions of the structure where heat must be efficiently transferred, high relative density ACMs are expected, whilst low relative density ACMs should occur in those regions where high insulation capability is required. This result can be easily achieved through a TO at the macroscopic scale. The corresponding GACMs can be recovered through a proper parametrisation of the RVE of the ACM to completely control the distribution of the material in terms of relative density for a prescribed RVE size (L m ) and a given topology. In particular, for ACMs composed of struts, the variable is the diameter of the strut cross-section d m , whilst for thin-walled TPMS ACMs the parameter varying over the structure is the thickness e m of the walls.

The strain energy-based homogenisation method

Homogenisation technique and associated boundary conditions

At the lower (i.e., mesoscopic or microscopic) scale, the RVE of the ACM can be interpreted, from a mechanical perspective, as a heterogeneous medium composed of two phases, i.e., the bulk material and the void. Conversely, at the macroscopic scale it can be modelled as an equivalent homogeneous anisotropic continuum whose mechanical response is described by a set of effective (or equivalent) material properties.

This work focuses only on the thermal behaviour of the ACM at the macroscopic scale, thus, the macroscopic thermal conductivity tensor (represented as a matrix K M ∈ R 3×3 through the Voigt's notation) of the ACM is determined by means of the SEHM. It is notewoerthy that the SEHM has been initially introduced to determine the equivalent elastic properties of heterogeneous media [START_REF] Barbero | Finite element analysis of composite materials[END_REF]. However, this technique can be easily extended to assess the macroscopic thermal behaviour of ACMs. For thermal problems, the basic assumption of the SEHM is that the thermal energy of the RVE is equal to the counterpart of the corresponding "envelope volume" of the homogeneous anisotropic medium replacing the ACM at the macroscopic scale. The effectiveness of the SEHM has been shown in different applications dealing with the determination of the macroscopic properties of heterogeneous materials characterised by complex RVE topologies [START_REF] Montemurro | A General Multi-Scale Topology Optimisation Method for Lightweight Lattice Structures Obtained through Additive Manufacturing Technology[END_REF][START_REF] Delucia | Determination of the effective thermoelastic properties of cork-based agglomerates[END_REF][START_REF] Refai | Determination of the effective elastic properties of titanium lattice structures[END_REF][START_REF] Bertolino | Multi-scale shape optimisation of lattice structures : an evolutionary-based approach[END_REF][START_REF] Barbero | Finite element analysis of composite materials[END_REF]. To compute K M , three further hypotheses have been considered: (a) the bulk material of the ACM has a linear thermal behaviour; (b) the thermal response of the RVE is evaluated at steady-state regime; (c) the only heat transfer phenomenon occurring in the RVE is the heat conduction and adiabatic BCs are applied to the internal boundary of the RVE (i.e., on those nodes which are not located on the RVE external faces).

To assess the components of K M , the RVE is submitted to a uniform temperature gradient ∇T mβ (that can be set by the user), (β = 1, 2, 3). The three analyses, corresponding to the three different forms of ∇T mβ , are executed one at time by considering the following set of periodic boundary conditions (PBCs) [START_REF] Barbero | Finite element analysis of composite materials[END_REF]:

T m (-a m1 , x m2 , x m3 ) -T m (a m1 , x m2 , x m3 ) = 2a m1 ∇T T mβ e m1 , -a mj ≤ x mj ≤ a mj , j = 2, 3, T m (x m1 , -a m2 , x m3 ) -T m (x m1 , a m2 , x m3 ) = 2a m2 ∇T T mβ e m2 , -a mj ≤ x mj ≤ a mj , j = 1, 3, T m (x m1 , x m2 , -a m3 ) -T m (x m1 , x m2 , a m3 ) = 2a m3 ∇T T mβ e m3 , -a mj ≤ x mj ≤ a mj , j = 1, 2, β = 1, 2, 3, (6)
where e mj is the unit vector along the x mj axis. Then, three steady-state thermal analyses are performed on the FE model of the RVE, by considering the PBCs of Eq. ( 6) and the expression of the user-defined temperature gradient reported in Tab. 1 (the value of ∆T m appearing in this table is arbitrary and should be set by the user).

1 st load case 2 nd load case 3 rd load case ∇T T m1 = ( ∆Tm 2a m1 , 0, 0) ∇T T m2 = (0, ∆Tm 2a m2 , 0) ∇T T m3 = (0, 0, ∆Tm 2a m3 )
Table 1: Temperature gradient to be used in the PBCs of Eq. ( 6).

Once the steady-state thermal analyses are solved, the components of the average heat flux per unit surface can be determined as:

q mi = 1 V RVE V EFF q mi (x m1 , x m2 , x m3 ) dV, i = 1, 2, 3, (7)
where V EFF is the actual volume of the RVE. The ETC tensor components K M ij can be calculated column-wise as follows:

K M ij = q mi ∇T T mj e mj , i, j = 1, 2, 3. (8)
Finally, the relative density of the ACM is defined as:

ρ M := V EFF V RVE . (9)

Lower and upper bounds on the effective thermal conductivity

When considering a homogenisation technique, the ETC of the ACM must always fall within appropriate bounds determined via a suitable variational method. These aspects are briefly discussed in the following.

Reuss-Voigt bounds

Under the hypothesis of linear thermal behaviour, the Reuss-Voigt bounds [START_REF] Voigt | Ueber die beziehung zwischen den beiden elasticitätsconstanten isotroper körper[END_REF][START_REF] Reuss | Berechnung der fließgrenze von mischkristallen auf grund der plastizitätsbedingung für einkristalle[END_REF] require only the knowledge of the volume fraction and of the thermal conductivity of each phase composing the heterogeneous medium. The lower and upper bounds, i.e., K RV M LB and K RV M UB , respectively, on the ETC can be derived by using the simple law of mixtures:

K RV M LB = ( n i=1 ν mi K -1 mi ) -1 , (10)
and

K RV M UB = n i=1 ν mi K mi . (11)
In Eqs. ( 10) and [START_REF] Fink | Fluid dynamics of flow through microscale lattice structures formed from self-propagating photopolymer waveguides[END_REF], n is the number of constitutive phases, while ν mi and K mi are the volume fraction and the thermal conductivity of the i-th phase. Of course, the following relationship holds:

n i=1 ν mi = 1. ( 12
)
For a material composed of two phases (e.g., void and bulk material as in the case of ACMs) Eq. ( 10) and ( 11) read:

K RV M LB = K m1 K m2 ν m1 K m2 + ν m2 K m1 , K RV M UB = ν m1 K m1 + ν m2 K m2 . (13)

Hashin-Shtrikman bounds

Calvo-Jurado and Parnell [START_REF] Calvo-Jurado | Hashin-Shtrikman bounds on the effective thermal conductivity of a transversely isotropic two-phase composite material[END_REF] derived Hashin-Shtrikman upper and lower bounds on the ETC of composite materials. The original derivation was carried out by means of a suitable variational method based on the minimisation of the total thermal energy of the material. The only hypothesis at the basis of the approach used in [START_REF] Calvo-Jurado | Hashin-Shtrikman bounds on the effective thermal conductivity of a transversely isotropic two-phase composite material[END_REF] are: 1) the material needs to be homogeneous on a scale larger than the one of any inhomogeneity (i.e., the generic inclusion); 2) the temperature must be continuous at the interface between two adjacent phases. For two-phases materials wherein the component with higher thermal conductivity is indicated as component 2, i.e., K m2 > K m1 , the Hashin-Shtrikman lower and upper bounds on the ETC read:

K HS M LB = K m1 K m2 +2K m2 (K m1 ν m2 +K m2 ν m1 ) 2K m2 +K m1 ν m1 +K m2 ν m2 , K HS M UB = K m1 K m2 +2K m1 (K m1 ν m2 +K m2 ν m1 ) 2K m1 +K m1 ν m1 +K m2 ν m2 . (14)

The finite element model

The FE model of each RVE topology has been generated within the commercial FE code ANSYS. Each FE model is built through an ad-hoc script wherein the mesh is properly related to the relevant geometric parameters of the RVE topology at hand. The generic RVE is modelled by means of 20-nodes solid elements (SOLID278), with one degree of freedom (DOF) per node. A FE model made of solid elements is needed to provide both a realistic representation of the RVE geometry and an accurate assessment of the 3D temperature field inside the RVE.

A convergence study in terms of the average element size (not reported here for the sake of brevity) has been conducted in order to assess its influence on the ETC of the different RVE topologies. The mesh size has been chosen equal to 0.03 mm for a cell size L m = 4 mm.

As stated above, three steady-state thermal analyses are conducted on each RVE topology with the PBCs of Eq. ( 6) and the imposed thermal gradient provided in Tab. 1 in order to determine the components of the macroscopic thermal conductivity tensor K M according to Eq. ( 8). As reported in Tab. 2, ACMs belonging to the first class, i.e., strut-based ACMs, show almost the same behaviour in terms of ETC, except for the Diamond topology, which shows poor thermal conductivity. In particular, in agreement with the results related to the macroscopic elastic behaviour presented in [START_REF] Refai | Determination of the effective elastic properties of titanium lattice structures[END_REF], all RVE topologies are characterised by a cubic syngony behaviour, in terms of the conductivity tensor at the macroscopic scale, which means K M ii = K M (i = 1, 2, 3) and K M ij = 0 (i, j = 1, 2, 3 with i = j). Moreover, one can notice that for some strut-based RVE topologies (e.g., the CBCC configuration which is a combination of CC and BCC topologies), the ETC may experience a drop. This is due to the slight increase in the overall heat surface exchange area of the RVE vs. the important decrease in volume (i.e., through the reduction of the strut diameter d m in order to fulfil the condition on the imposed relative density): these details are provided in Tab. 3.

Strut

Influence of the relative density on the effective thermal conductivity

The influence of the relative density of the RVE ρ M on the macroscopic thermal conductivity K M is investigated in the following. For the sake of brevity, the analysis is limited only to the best and worst RVE topologies, in terms of ETC, of each class. More precisely, the following topologies are considered: diamond, truncated cuboctahedron, S-IWP, S-iBCC, T-IWP, T-Gyroid, T-Manta. The ETC vs. the relative density for each RVE topology is illustrated in Fig. 6. These curves have been obtained by launching the homogenisation procedure discussed in Sec. 3 for ten different values of ρ M in the range [0, 1]. It is noteworthy that, numerical results provided by the FE-based homogenisation scheme always fall within both Reuss-Voigt (VR-LB and VR-UB) and Hashin-Shtrikman bounds (HS-LB and HS-UB). The FE analyses have been carried out for each RVE topology in order to build a database of results. Subsequently, results have been fitted through a polynomial of degree 7 as follows:

φ(ρ M ) := K M K m = 7 i=1 c i ρ i M , c i ∈ R, i = 1, . . . , 7, (15)
where c i are suitable coefficients calculated to minimise the distance between the approximating curve and the database of results (for a given RVE topology). The values of these coefficients are reported in Tab.

Design of graded architected cellular material through topology optimisation

Fundamentals of NURBS hyper-surfaces

A NURBS hyper-surface is a polynomial-based function, defined as h : R N → R M . The mathematical formula of a generic NURBS hyper-surface reads

h(ζ 1 , . . . , ζ N ) := n 1 i 1 =0 • • • n N i N =0 R i 1 ,...,i N (ζ 1 , . . . , ζ N )P i 1 ,...,i N , (16)
where n j + 1 (j = 1, . . . , N ) is the number of CPs along the ζ j parametric direction, R i 1 ,...,i N (ζ 1 , . . . , ζ N ) are the piece-wise rational basis functions, which are related to the standard Bernstein's polynomials

N i k ,p k (ζ k ), k = 1, . . . , N through the following formula R i 1 ,...,i N (ζ 1 , . . . , ζ N ) = ω i 1 ,...,i N N k=1 N i k ,p k (ζ k ) n 1 j 1 =0 • • • n N j N =0 ω j 1 ,...,j N N k=1 N j k ,p k (ζ k ) . (17)
In Eqs. ( 16) and ( 17), h(ζ 1 , . . . , ζ N ) is a M -dimension vector-valued rational function, (ζ 1 , . . . , ζ N ) are scalar dimensionless parameters defined in the interval [0, 1], whilst P i 1 ,...,i N are the CPs coordinates. The j-th CP coordinate X (j)

i 1 ,...,i N is stored in the array X (j) ∈ R (n 1 +1)ו••×(n N +1)
. The explicit expression of CPs coordinates in R M is:

P i 1 ,...,i N = {X (1) i 1 ,...,i N , . . . , X (M ) i 1 ,...,i N }, X (j) ∈ R (n 1 +1)ו••×(n N +1) , j = 1, . . . , M. (18)
Curves and surfaces formulae can be easily deduced from Eq. ( 16). The CPs layout is referred to as control hyper-net [START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF]. The total number of CPs constituting the hyper-net is:

n CP := N i=1 (n i + 1). ( 19
)
The generic CP affects the shape of the NURBS entity through its coordinates. Moreover, a weight w i 1 ,...,i N is associated to the generic CP: the higher the weight w i 1 ,...,i N , the more the NURBS entity is attracted towards the related CP. For each parametric direction ζ k , k = 1, . . . , N , the Bernstein's polynomials are of degree p k and can be generated in a recursive way as

N i k ,0 (ζ k ) := 1, if v (k) i k ≤ ζ k < v (k) i k +1 , 0, otherwise, (20)
N i k ,q (ζ k ) := ζ k -v (k) i k v (k) i k +q -v (k) i k N i k ,q-1 (ζ k ) + v (k) i k +q+1 -ζ k v (k) i k +q+1 -v (k) i k +1 N i k +1,q-1 (ζ k ), q = 1, . . . , p k , (21)
where each blending function is defined on the knot vector v (k) := {0, . . . , 0

p k +1 , v (k) p k +1 , . . . , v (k) m k -p k -1 , 1, . . . , 1 p k +1 }, (22)
whose dimension is m k + 1, with

m k = n k + p k + 1. ( 23
)
Each knot vector v (k) is a non-decreasing sequence of real numbers. The NURBS blending functions are characterised by several interesting properties: the interested reader is addressed to [START_REF] Piegl | The NURBS book[END_REF] for a deeper insight into the matter. Here, only the local support property is recalled because it is of paramount importance for the NURBS-based SIMP method [START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | NURBS Hypersurfaces for 3D Topology Optimisation Problems[END_REF]:

R i 1 ,...,i N (ζ 1 , . . . , ζ N ) = 0, if (ζ 1 , . . . , ζ N ) ∈ v (1) i 1 , v (1)
i 1 +p 1 +1 × • • • × v (N ) i N , U (N ) i N +p N +1 . ( 24
)
According to the above formula, the generic CP (and the associated weight) affects only a precise zone of the NURBS parametric space, which is indicated as local support.

The NURBS-density-based topology optimisation method

The details of the formulation of the density-based TO method in the NURBS hypersurfaces framework are given in [START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | NURBS Hypersurfaces for 3D Topology Optimisation Problems[END_REF]. The main features of the proposed approach are briefly recalled here. The formulation is presented for 3D steady-state heat conduction problems under the hypothesis that the applied thermal loads and BCs do not depend upon the pseudo-density field (of course, this hypothesis can be easily relaxed).

Design variables

Consider the compact space D M ⊂ R 3 in a Cartesian orthonormal frame O(x M 1 , x M 2 , x M 3 ) defined at the macroscopic scale:

D M := {x T M = {x M 1 , x M 2 , x M 3 } ∈ R 3 : x M 1 ∈ [0, a M 1 ], x M 2 ∈ [0, a M 2 ], x M 3 ∈ [0, a M 1 ]}, ( 25
)
where a M j (j = 1, 2, 3) is a reference length defined along x M j axis. The goal of TO is to search for the best distribution of a fictitious material satisfying the requirements of the design problem.

In density-based TO methods, the material domain at the macroscopic scale Ω M ⊆ D M is identified by means of a pseudo-density function ρ M (x M ) ∈ [0, 1] for x M ∈ D M : ρ M (x M ) = 0 denotes absence of material, whilst ρ M (x M ) = 1 indicates presence of material. In the context of the NURBS-density-based TO method, the pseudo-density field for a TO problem of dimension D is represented through a NURBS hyper-surface of dimension D + 1. Therefore, for a 3D problem a 4D entity is needed and the pseudo-density field is defined as:

ρ M (ζ M 1 , ζ M 2 , ζ M 3 ) = n M 1 i 1 =0 n M 2 i 2 =0 n M 3 i 3 =0 R i 1 ,i 2 ,i 3 (ζ M 1 , ζ M 2 , ζ M 3 )ρ M i 1 ,i 2 ,i 3 . ( 26
)
In Eq. ( 26),

n M CP = (n M 1 + 1)(n M 2 + 1)(n M 3 + 1) is the total number of CPs, ρ M (ζ M 1 , ζ M 2 , ζ M 3
) constitutes the fourth coordinate of the array h of Eq. ( 16), while

R i 1 ,i 2 ,i 3 (ζ M 1 , ζ M 2 , ζ M 3
) are the basis functions of Eq. ( 17). The dimensionless parameter ζ M j is defined as:

ζ M j := x M j a M j , j = 1, 2, 3. ( 27
)
It is noteworthy that, according to the methodology discussed in Sec. 2 and illustrated in Fig. 4, the pseudo-density ρ M of Eq. ( 26) describes the GACM distribution at the macroscopic scale. In the following, only the pseudo-density at CPs and the associated weights are selected as design variables and are grouped in the vectors ξ 1 and ξ 2 , respectively, defined as:

ξ T 1 := (ρ M 0,0,0 , . . . , ρ M n M 1 ,n M 2 ,n M 3 ), ξ T 2 := (w M 0,0,0 , . . . , w M n M 1 ,n M 2 ,n M 3 ), ξ 1 , ξ 2 ∈ R n M CP . ( 28
)
Therefore, in the most general case, the overall number of design variables is n var = 2n M CP .

Objective function and optimisation constraints

Consider the steady-state equation of the FE model in the most general case:

KM ûM = fM ; ûM , fM ∈ R NMDOF , KM ∈ R NMDOF × NMDOF , (29)
where NMDOF represents the overall number of DOFs of the FE model at the macroscopic scale before applying the BCs, KM is the non-reduced (singular) conductivity matrix of the FE model, while fM and ûM are the non-reduced vectors of the external generalised nodal thermal forces and temperatures, respectively. Using standard FE notation [START_REF] Zienkiewicz | The Finite Element Method: Its Basis and Fundamentals[END_REF], Eq. ( 29) can be rewritten as:

K M K M BC K T M BC KM u M u M BC = f M r M , (30)
with:

u M , f M ∈ R N M DOF , u M BC , r M ∈ R N M BC , K M ∈ R N M DOF ×N M DOF , K M BC ∈ R N M DOF ×N M BC , KM ∈ R N N BC ×N M BC , NMDOF = N M DOF + N M BC , (31)
where N M BC represents the number of DOFs where temperature is imposed (Dirichlet's BC), while N M DOF is the number of unknown DOFs . In Eq. ( 30), u M and u M BC are the unknown and imposed vectors of nodal temperatures, respectively. f M is the vector of generalised external nodal thermal forces, whilst r M is the vector of (unknown) generalised nodal thermal reactions at nodes where BCs are imposed. K M , K M BC and KM are the conductivity matrices of the FE model after applying BCs. Consider, now, the case of zero Dirichlet's BCs and non-zero Neumann's BCs: the thermal compliance of the structure is defined as

C M := f T M u M . (32)
According to Eq. ( 15), the pseudo-density field ρ M affects the element conductivity matrix and, accordingly, the global conductivity matrix of the FE model as

K M := N M e e=1 φ e L T M e K 0 M e L M e = N M e e=1 L T M e K M e L M e , K 0 M e , K M e ∈ R N e M DOF ×N e M DOF , L M e ∈ R N e M DOF ×N M DOF (33)
where φ e = φ(ρ M e ) is the function of Eq. ( 15) computed at the centroid of the generic element e. N M e is the total number of elements of the FE model at the macroscopic scale and N e M DOF is the number of DOFs of the generic element. In Eq. ( 33), K 0 M e and K M e are the non-penalised and the penalised conductivity matrices of element e, expressed in the global reference frame of the FE model, whilst L M e is the connectivity matrix of element e. The non-penalised conductivity matrix of the element can be related to the thermal conductivity matrix of the bulk material of the RVE K m through the following relationship:

K 0 M e := V M e B T M e K m B M e dV, K m ∈ R 3×3 , B M e ∈ R 3×N e M DOF , (34)
where B M e is the matrix representing the product between the linear differential operator and the shape function matrices of the generic element.

The TO problem of thermal compliance minimisation subject to an inequality constraint on the volume can be formulated as a constrained non-linear programming problem (CNLPP) as follows:

min ξ 1 ,ξ 2 C M C M ref , s.t. :                K M u M = f M , V M V M ref -γ ≤ 0, ξ 1k ∈ [ρ M min , ρ M max ], ξ 2k ∈ [w M min , w M max ], ∀k = 1, . . . , n M CP . (35)
In Eq. ( 35), V M ref is a reference volume, V M is the volume of the material domain Ω at the macroscopic scale, while γ is the fixed volume fraction; ρ min represents the lower bound, imposed to the density field to prevent any singularity for the solution of the equilibrium problem. The objective function is divided by a reference compliance, C Mref , to obtain a dimensionless value. The volume of the material domain appearing in Eq. ( 35) is defined as:

V M := N M e e=1 ρ M e V M e , (36)
where V M e is the volume of element e. Moreover, in Eq. ( 35), the linear index k has been introduced for the sake of compactness. The relationship between k and i j , (j = 1, 2, 3) is:

k := 1 + i 1 + i 2 (n M 1 + 1) + i 3 (n M 1 + 1)(n M 2 + 1). ( 37
)
The other parameters involved in the definition of the NURBS entity (i.e., degrees, knot-vector components and number of CPs) are set a-priori at the beginning of the TO analysis and are not optimised. The computation of the gradient of both objective function and optimisation constraints with respect to the design variables is needed to perform the solution search of the CNLPP of Eq. ( 35) through a deterministic algorithm. This task is achieved by exploiting the local support property of Eq. [START_REF] Chamberlain | Derivation of expansion coefficients for a fibre reinforced composites[END_REF]. For instance, the general expressions of the gradient of both the thermal compliance (in the case u M BC = 0) and the volume read

∂C M ∂ξ ik = - e∈S k C M e φ e ∂φ e ∂ρ M e ∂ρ M e ∂ξ ik , i = 1, 2, k = 1, . . . , n M CP , (38)
∂V M ∂ξ ik = e∈S k V M e ∂ρ M e ∂ξ ik , i = 1, 2, k = 1, . . . , n M CP , (39)
with

C M e := u T M L T M e K M e L M e u M . (40)
In Eqs. ( 38) and ( 39), S k is the discretised version of the local support of Eq. ( 24), while

∂ρ M e ∂ξ ik reads ∂ρ M e ∂ξ ik =      R ke , if i = 1, R ke ξ 2k (ξ 1k -ρ M e ) , if i = 2. (41)
The scalar quantity R ke , appearing in Eq. ( 41), is the NURBS rational basis function of Eq. ( 17) evaluated at the element centroid. More details on the analytical passages to derive the gradient of each response function are available in [START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | NURBS Hypersurfaces for 3D Topology Optimisation Problems[END_REF].

Reconstruction of graded architected cellular materials topologies

Once the result of the TO are available, a post-processing phase is needed to recover point-wise the RVE geometry corresponding to the local optimal value of the pseudodensity field at the macroscopic scale. Although this process is quite easy for strut-based ACMs (for which simple formulae are available to determine the effective volume at the RVE scale), this task reveals anything but trivial in the case of thin-walled or skeletal TPMS ACMs. Therefore, the algorithm 1 has been implemented to recover the RVE geometry for GACMs made of thin-walled (or skeletal) TPMS RVEs.

Algorithm 1 Recover the GACM from TO results (case of thin-walled TPMS RVE). determine e m corresponding to ρ M e from the database of results

4:

generate the iso-surface of the thin-walled TPMS RVE for t m = ±e m /2 in Eqs.

(1)-( 4)

5:
fill the volume between the boundary surfaces 6: end for 7: assemble the RVEs to generate the GACM 8: generate the Standard Tessellation Language file of the GACM Of course, a similar algorithm has been developed for skeletal TPMS RVE topologies. An example of the result of the reconstruction phase described in algorithm 1 applied to the case of the Gyroid topology for different values of ρ M is illustrated in Fig. 7. Moreover, two examples of GACMs with a linear variation of ρ M from 0.1 to 0.5 along z axis are shown in Fig. 8 for Schwartz and Gyroid topologies, respectively. Remark 5.1. As illustrated in algorithm 1, the recovery operation is performed for each element composing the macroscopic domain (in this case it is tacitly assumed that one element at the macroscopic scale corresponds to one RVE). In particular, the value of the pseudo-density field at the end of the optimisation process is available in the form of a NURBS entity which can be projected over the mesh (and thus over each element composing the design domain) at the macroscopic scale. Remark 5.2. Depending on the characteristic size of both the RVE and of the design domain at the macroscopic scale, each element constituting the mesh of the macroscopic FE model can correspond to one or more RVEs. Indeed, this aspect can be easily handled by the proposed approach because the pseudo-entity field is available in the form of a NURBS entity at the end of the optimisation process. Therefore, the mesh used to recover the GACM during the post-processing phase can be different from the mesh of the FE model used to calculate the physical responses. density field at the macroscopic scale on the optimised topology is studied; (2) the influence of the RVE geometry at the lower scale on the optimised topology at the upper scale is investigated.

As far as the design space of the CNLPP of Eq. ( 35) is concerned, lower and upper bounds of design variables are set as: ρ M min = 10 -3 , ρ M max = 1; w M min = 0.5, w M max = 10. It is noteworthy that the non-trivial components of the knot-vectors in Eq. ( 22) are evenly distributed in the interval [0, 1] for each benchmark problem. Finally, according to [START_REF] Alexandersen | Topology optimization for convection problems[END_REF][START_REF] Montemurro | A topology optimization method based on non-uniform rational basis spline hyper-surfaces for heat conduction problems[END_REF] the thermal conductivity of the material used for all the considered test cases is K m = 1 Wm -1 K -1 .

2D Benchmark problem

The 2D benchmark problem (BK-2D), illustrated in Fig. 9, is characterised by inplane dimensions a M 1 = a M 2 = 20 m, and a thickness t M = 1 m. Regarding the BCs, the temperature is set equal to zero for nodes located at The optimised topology, for each RVE geometry, is illustrated in Figs. 10 and11, for B-spline and NURBS solutions, respectively. In each figure, results are provided in terms of macroscopic thermal compliance C M and number of iterations N iter to achieve convergence. For each solution, the requirement on the volume fraction is always satisfied and the solution is located on the boundary of the feasible domain. A synthesis of the results is illustrated in Fig. 12 where the thermal compliance characterising each optimised topology is plotted for each RVE geometry and for both B-spline and NURBS solutions. From the analysis of these results, one can infer the following remarks:

x M 1 ∈[ a M 1 -a M t 2 , a M 1 +a M t 2 ], x M 2 = a
For each RVE geometry, the optimised topology at the macroscopic scale obtained through a NURBS entity is characterised by a thermal compliance lower than the one characterising the B-spline counterpart; moreover, in agreement with the results found in [START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Montemurro | A topology optimization method based on non-uniform rational basis spline hyper-surfaces for heat conduction problems[END_REF] each NURBS solution is characterised by a boundary smoother than the B-spline counterpart.

In agreement with the results of the homogenization presented in Sec. 4, better performances (in terms of thermal compliance) can be achieved when using thinwalled TPMS RVE topologies; particularly, the best result is obtained in the case of the Manta configuration, whilst the worst one is obtained in the case of the Diamond topology.

The RVE topology at the lower scale affects not only the value of the objective function but also the number of branches (and their curvature) of the optimised topology at the macroscopic scale.

An example of GACM recovered from the results of the TO is illustrated in Fig. 13. In particular, the optimised topology obtained with a B-spline entity related to the Gyroid RVE geometry shown in Fig. 10 has been reconstructed by following the procedure described in algorithm 1. It is noteworthy that the regions coloured in blue in Fig. 13 represent zones wherein the relative density ρ M = 1, thus the RVE at the lower scale is completely filled by the bulk material.

3D Benchmark problem

The 3D benchmark problem (BK-3D) deals with the TO of a cubic domain illustrated in Fig. 14; the geometrical parameters of BK-3D are: a M 1 = a M 2 = a M 3 = 20 m, a M t = 2 m. The material properties are the same as those used in the test case BK-2D.

The FE model is made of N M e = 64000 SOLID279 elements (eight nodes, one DOF per node). The temperature is equal to zero on the nodes belonging to the heat skin located at

(x M 1 ∈[ a M 1 -a M t 2 , a M 1 +a M t 2 ], x M 2 ∈[ a M 2 -a M t 2 , a M 2 +a M t 2 ],
x M 3 = a M 3 ). A heating source s h = 0.001 Wm -3 is evenly distributed over the design domain. Regarding problem [START_REF] Da | Design of materials using hybrid cellular automata[END_REF], the reference volume and the volume fraction are V M ref = a M 1 a M 2 a M 3 and γ = 0.3, respectively. An initial guess characterised by a uniform pseudodensity field ρ M (ξ 1 , ξ 2 , ξ 3 ) = γ is considered. The reference thermal compliance C M ref associated to the starting point is reported in Tab. 9 for each RVE considered in the analysis. Moreover, the optimisation has been carried out by considering both B-spline and NURBS entities characterised by the following integer parameters: p j = 2 (j = 1, 2, 3), n M CP = 30×30×30. The CNLPP formulation of Eq. ( 35) has been enhanced by adding a double orthogonal symmetry constraint with respect to planes x M 1 = a M 1 2 and x M 2 = a M 2 2 . The optimised topologies are illustrated in Figs. 15 and16, for both B-spline and NURBS entities: numerical results are provided in terms of thermal compliance C M of the optimised solution and number of iterations N iter to achieve convergence. For each solution the requirement on the volume fraction is always satisfied. A synthesis of results is reported in Fig. 17. The same remarks already done for benchmark problem BK-2D can be repeated here.

RVE

Conclusions

In this work, a general approach for the optimisation of GACMs has been presented. In particular, the proposed methodology relies on: a) the SEHM to perform the scale transition and determine the macroscopic properties of the RVE of the ACM; b) NURBS hyper-surfaces to represent the topological variable at the macroscopic scale; c) a densitybased TO algorithm reformulated in the framework of NURBS hyper-surfaces; d) a general strategy to recover the GACM pattern from the results of the TO at the macroscopic scale of the structure. This methodology has been applied to the problem of minimising the thermal compliance of GACMs subject to a constraint on the volume and its effectiveness has been shown on both 2D and 3D test cases. In this context, the properties of NURBS entities are efficiently exploited to derive the analytical expression of the gradient of the physical responses integrated in the problem formulation. Some features of the proposed methodology need to be highlighted.

1. The results of the homogenisation procedure show that the ETC depends upon the relative density of the ACM. In particular, thin-walled TPMS ACMs are characterised by an ETC higher than strut-based and skeletal TPMS configurations. The T-Manta topology shows the highest ETC, while the one of the T-Nevious geometry is slightly lower. 2. TO based on NURBS hyper-surfaces is characterised by three advantages: (a) unlike the classical SIMP approach, the NURBS local support establishes an implicit relationship among the pseudo-density of contiguous mesh elements, thus an explicit filter is not required; (b) the number of design variables is reduced with respect to the classical SIMP approach; (c) the CAD reconstruction of the boundary of the optimised topology is an easy task. 3. A sensitivity analysis of the optimised topology at the macroscopic scale to the RVE geometry at the lower scale has been performed. The RVE topology at the lower scale influences not only the optimal distribution of the pseudo-density at the macroscopic scale but also the structural response in terms of thermal compliance.

Among the considered RVE topologies, the best results can be achieved by the T-Manta configuration followed by the T-Gyroid and T-IWP topologies. 4. The role of NURBS weights has been assessed: by considering the same number of CPs and the same degrees, the thermal compliance of the NURBS solution is lower than the one of the B-spline counterpart. 5. An efficient and automatic post-processing strategy to recover GACMs patterns from the results of the TO at the macroscopic scale has been proposed. The algorithm makes use not only of the optimal distribution of the pseudo-density field resulting from the TO process but also of the database of results obtained through the SEHM before the optimisation process.

Regarding the prospects on this topic, several challenges still need to be faced. Firstly, the transient regime should be integrated within the TO problem formulation. Secondly, the proposed strategy should be extended to thermomechanical problems (by considering both weak and strong couplings). Thirdly, design-dependent BCs should be integrated into the problem formulation to deal with more realistic engineering applications, like heat exchangers made of GACMs. Finally, a dedicated experimental campaign of tests to measure the thermal conductivity of optimised heat exchangers obtained with the proposed methodology is in progress. Research is ongoing on all the above aspects.

Figure 1 :

 1 Figure 1: Strut-based ACMs [42].

Figure 2 :Figure 3 :

 23 Figure 2: (a) Thin-walled and (b) skeletal TPMS ACMs [42].

Figure 4 :

 4 Figure 4: Work-flow of the design methodology for GACMs.

Figure 5 :

 5 Figure 5: Examples of mesh quality of the FE model of the RVE for two different topologies.

4 . 4 . 1 .

 441 Numerical results: homogenisation of the architected cellular materials and database generation Effective thermal conductivity calculation A first campaign of numerical analyses has been carried out to compare the ETC of different RVE topologies, which are summarised in Tabs. 2, 4 and 5, for strut-based, thin-walled and skeletal ACMs, respectively, by considering a relative density ρ M = 0.1. It is noteworthy that, for each RVE topology, the characteristic length of the cubic RVE has been set as L m = 4 mm, while the strut diameter d m (for strut-based ACMs) and the wall thickness e m (for thin-walled ACMs) have been computed to meet the condition on the relative density of the RVE, i.e., ρ M = 0.1.

Figure 6 :

 6 Figure 6: Effective thermal conductivity vs. the relative density for the best and worst RVE topologies for each class.

1 :

 1 for e ∈ [1, N M e ] do 2:recover ρ M e from the results of the TO 3:

Remark 5 . 3 .

 53 The discontinuities shown in Figs.7 and 8are only apparent since they are due to the periodicity of the trigonometric functions used to generate TPMS geometries (which are obtained after a threshold operation). After the assembly of all RVEs, these discontinuities occur only on the external boundary of the macroscopic domain (due to the periodicity of the micro-structure of the material) and can be easily removed by applying standard Boolean operations.

6 Figure 7 :

 67 Figure 7: Examples of reconstruction of the T-Gyroid RVE topology for different values of relative density.

  M 2 with a M t = 2 m, while a heating source s h = 0.001 Wm -2 is uniformly distributed over the whole design domain. The FE model is made of N M e = 80 × 80 PLANE55 elements (4 nodes with a single DOF per node). The reference volume V M ref in Eq. (35) is the overall volume of the design domain, i.e., V M ref = a M 1 a M 2 t M . The volume fraction has been set as γ = 0.3. An initial guess characterised by a uniform pseudo-density field ρ M (ξ 1 , ξ 2 ) = γ has been considered for each analysis. An extensive campaign of numerical tests has been conducted on BK-2D: the goal is to assess the influence of the RVE topology and of the geometric entity (B-spline or NURBS surface) on the optimised topology at the macroscopic scale. In particular, the CNLPP of Eq. (35) is solved by considering the following RVE topologies: Diamond, Truncatedcuboctahedron, S-iBCC, S-IWP, T-Gyroid, T-IWP and T-Manta. For each case, B-spline and NURBS surfaces with blending functions degrees p j = 2 (j = 1, 2) and number of CPs n M CP = 68 × 68 are employed. Moreover, a symmetry constraint with respect to axis x M 1 = a M 1 2 has been considered in the CNLPP formulation of Eq. (35). The reference value of the thermal compliance C M ref appearing in Eq. (35) corresponds to the thermal compliance of the initial guess and is listed in Tab. 8 for each RVE topology.

Figure 9 :

 9 Figure 9: BK-2D: geometry and boundary conditions.

  (a) Diamond: C M = 0.26897 WK, N iter = 236 (b) Truncated-cuboctahedron: C M = 0.25154 WK, N iter = 222 (c) S-iBCC: C M = 0.24410 WK, N iter = 246 (d) S-IWP: C M = 0.25785 WK, N iter = 240 (e) T-Gyroid: C M = 0.23613 WK, N iter = 288 (f) T-IWP: C M = 0.24291 WK, N iter = 243 (g) T-Manta: C M = 0.23211 WK, N iter = 248

Figure 10 :

 10 Figure 10: BK-2D: sensitivity of the optimised topology to the RVE geometry using B-spline entities with pj = 2 (j = 1, 2) and nMCP = 68 × 68.

Figure 11 :

 11 Figure 11: BK-2D: sensitivity of the optimised topology to the RVE geometry using NURBS entities with pj = 2 (j = 1, 2) and nMCP = 68 × 68.

Figure 12 :

 12 Figure 12: BK-2D: synthesis of results.

Figure 13 :

 13 Figure 13: BK-2D: optimised GACM made of gyroid RVE reconstructed from B-spline solution of Fig. 10.

Figure 14 :

 14 Figure 14: BK-3D: geometry and boundary conditions.

  (a) Diamond: C M = 14.69802 WK, N iter = 244 (b) Truncated-cuboctahedron: C M = 14.27048 WK, N iter = 275 (c) S-iBCC: C M = 14.10750 WK, N iter = 166 (d) S-IWP: C M = 14.43974 WK, N iter = 228 (e) T-Gyroid: C M = 13.93192 WK, N iter = 150 (f) T-IWP: C M = 14.08170 WK, N iter = 165 (g) T-Manta: C M = 13.83650 WK, N iter = 181

Figure 15 :

 15 Figure 15: BK-3D: sensitivity of the optimised topology to the RVE geometry using B-spline entities with pj = 2 (j = 1, 2, 3) and nMCP = 30 × 30 × 30.

Figure 16 :

 16 Figure 16: BK-3D: sensitivity of the optimised topology to the RVE geometry using NURBS entities with pj = 2 (j = 1, 2, 3) and nMCP = 30 × 30 × 30.

Figure 17 :

 17 Figure 17: BK-3D: synthesis of results.

Table 2 :

 2 Effective thermal conductivity of strut-based RVE topologies for ρM = 0.1 and Lm = 4 mm.The ETC of thin-walled and skeletal RVE topologies is listed in Tabs. 4 and 5, respectively. As it can be inferred from these results, strut-based and skeletal ACMs are characterised by an ETC lower than the one of thin-walled ACMs. In particular, among the

	-based RVE	K M [W/mK]
	ACC	7.12
	BCC	7.09
	BFCC	7.03
	CBCC	6.43
	CBFCC	7.13
	CC	7.16
	CIC	6.65
	CFCC	7.11
	Diamond	4.93
	FCC	7.07
	FCC2+	7.39
	Hexahedral	7.19
	Octet	7.09
	Rhombic	7.28
	Truncated-cube	6.99
	Truncated-cuboctahedron	7.42
	Truncated-cuboctahedron2+	6.91

Table 3 :

 3 Value of the strut diameter and heat surface exchange area of BCC, CC and CBCC RVE topologies for ρM = 0.1 and Lm = 4 mm.

	investigated thin-walled geometries, the highest ETC occurs for T-Manta and T-Neovius
	configurations.
	Thin-walled RVE K M [W/mK]
	T-Batwing	10.66
	T-BTC	10.22
	T-F-RD	10.63
	T-Gyroid	10.29
	T-IWP	8.89
	T-Manta	11.32
	T-Neovius	11.21
	T-Schwarz	10.43

Table 4 :

 4 Effective thermal conductivity of thin-walled RVE topologies for ρM = 0.1 and Lm = 4 mm.

	Skeletal RVE K M [W/mK]
	S-BTC	7.19
	S-Gyroid	8.04
	S-iBCC	8.59
	S-IWP	6.35
	S-OrthoCircle	6.72
	S-Schwarz	7.21

Table 5 :

 5 Effective thermal conductivity of skeletal RVE topologies for ρM = 0.1 and Lm = 4 mm.

Table 6 :

 6 [START_REF] Borovinšek | Designing 2d auxetic structures using multi-objective topology optimization[END_REF] for the best and worst RVE topologies of each ACM class. Values of the coefficients ci (i = 1, . . . , 7) of Eq. (15) for the best and worst RVE topologies for each ACM class.

	RVE topology	c 1	c 2	c 3	c 4	c 5	c 6	c 7
	Diamond	0.2753 0.1338 0.1194 0.1180 0.1178 0.1178 0.1178
	Truncated-cuboctahedron 0.4231 0.1236 0.0933 0.0902 0.0899 0.0899 0.0899
	S-iBCC	0.4925 0.1188 0.0810 0.0772 0.0768 0.0768 0.0768
	S-IWP	0.3596 0.1280 0.1045 0.1021 0.1019 0.1019 0.1019
	T-Gyroid	0.5934 0.1119 0.0631 0.0583 0.0578 0.0577 0.0577
	T-IWP	0.5103 0.1176 0.0778 0.0739 0.0735 0.0734 0.0734
	T-Manta	0.6545 0.1077 0.0523 0.0468 0.0462 0.0462 0.0462

Table 8 :

 8 BK-2D: reference thermal compliance for each RVE topology.

	RVE	C M ref [WK]
	Diamond	1.52642
	Truncated-cuboctahedron	1.06614
	S-iBCC	0.93382
	S-IWP	1.22485
	T-Gyroid	0.79116
	T-IWP	0.90502
	T-Manta	0.72413

Table 9 :

 9 BK-3D: reference thermal compliance for each RVE topology.

		C M ref [WK]
	Diamond	128.98658
	Truncated-cuboctahedron	90.09139
	S-iBCC	78.91058
	S-IWP	103.50328
	T-Gyroid	66.85503
	T-IWP	76.47695
	T-Manta	61.19098

Acknowledgements

Part of this work has been funded by Nouvelle-Aquitaine region through OCEAN-ALM project.

Data availability

The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

Numerical results: topology optimisation

The effectiveness of the proposed method is illustrated on 2D and 3D benchmark problems taken from the literature [START_REF] Alexandersen | Topology optimization for convection problems[END_REF][START_REF] Montemurro | A topology optimization method based on non-uniform rational basis spline hyper-surfaces for heat conduction problems[END_REF]. For each case, the pseudo-density field and the optimum topology are shown. The results presented here are obtained through the code SANTO (SIMP and NURBS for topology optimisation) developed at the I2M laboratory in Bordeaux [START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | NURBS Hypersurfaces for 3D Topology Optimisation Problems[END_REF]. SANTO is coded in the Python ® environment and can be interfaced with any FE code. In this study, the commercial code ANSYS ® is used to build the FE model at the macroscopic scale and to assess the structural responses, i.e., nodal temperature and thermal compliance.

Moreover, the globally-convergent method of moving asymptotes (GCMMA) algorithm [START_REF] Svanberg | A class of globally convergent optimization methods based on conservative convex separable approximations[END_REF] is employed to carry out the solution search for the CNLPP of Eq. ( 35). The parameters governing its behaviour are listed in Table 7. Regarding the numerical analyses, the following aspects are taken into account: (1) the influence of the geometric entity, i.e., B-spline or NURBS, used to describe the pseudo-