Giulia Bertolino

Marco Montemurro
email: marco.montemurro@ensam.eu

Two-scale topology optimisation of cellular materials under mixed boundary conditions

Keywords:

This work proposes a theoretical/numerical framework for the topology optimisation of anisotropic architected cellular materials at different scales. In particular, the topological variable (i.e., the pseudo-density field) is defined at both the scale of the representative volume element (i.e., the unit cell) of the material and at the macroscopic scale of the structure. The two-scale topology optimisation problem is formulated in the most general sense, i.e., by considering non-zero Neumann-Dirichlet boundary conditions. The proposed method is based on: (a) non-uniform rational basis spline hyper-surfaces to represent the topological variable at each scale, (b) the solid isotropic material with penalisation approach, (c) a general numerical homogenisation scheme based on the strain energy to establish the link between scales. The proposed formulation exploits the properties of non-uniform rational basis spline entities to determine the relationships occurring among the topological variables defined at different scales to correctly state the optimisation problem and to satisfy the hypotheses at the basis of the homogenisation method. In particular, scales separation (a necessary condition to be met in order to apply the homogenisation method) and manufacturing requirements are implicitly ensured by introducing minimum length scale constraints on the topological variables defined at both macroscopic scale and unit cell scale, respectively. Furthermore, the sensitivity of the optimised topology (at each scale) to the applied boundary conditions and to the elastic symmetry group of the representative volume element is investigated by founding new and original results. The effectiveness of the approach is tested on 2D and 3D benchmark problems taken from the literature.

Introduction

Nowadays architected cellular materials (ACMs) are of great interest in different engineering fields, like aerospace, automotive, bio-mechanical industries or in the energetic and chemical sectors, thanks to the possibility of designing the micro/meso structure to get tailored performances, e.g., high stiffness-to-weight and strength-to-weight ratios, high permeability, energy absorption and thermal insulation.

To this end, the scientific community is developing different approaches to optimise ACMs: parametric optimisation of the geometrical variables of pre-defined representative volume element (RVE) topologies [START_REF] Tang | Multifunctional design of heterogeneous cellular structures[END_REF][START_REF] Benedetti | Architected cellular materials: A review on their mechanical properties towards fatigue-tolerant design and fabrication[END_REF], topology optimisation (TO) of the RVE to satisfy a given macroscopic elastic behaviour [START_REF] Andreassen | Topology optimization of periodic microstructures for enhanced dynamic properties of viscoelastic composite materials[END_REF][START_REF] Collet | Topology optimization for microstructural design under stress constraints[END_REF][START_REF] Borovinšek | Designing 2d auxetic structures using multi-objective topology optimization[END_REF][START_REF] Zheng | Evolutionary topology optimization for mechanical metamaterials with auxetic property[END_REF][START_REF] Cheng | Efficient design optimization of variable-density cellular structures for additive manufacturing: theory and experimental validation[END_REF][START_REF] Gao | Isogeometric topology optimization for computational design of re-entrant and chiral auxetic composites[END_REF][START_REF] Li | Design of architected materials for thermoelastic macrostructures using level set method[END_REF][START_REF] Wallin | Nonlinear homogenization for topology optimization[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] Wang | Topological shape optimization of microstructural metamaterials using a level set method[END_REF] and TO of the ACM at multiple scales [START_REF] Li | Design of architected materials for thermoelastic macrostructures using level set method[END_REF][START_REF] Hoang | Robust multiscale design of incompressible multi-materials under loading uncertainties[END_REF][START_REF] Chougrani | Parts internal structure definition using non-uniform patterned lattice optimization for mass reduction in additive manufacturing[END_REF][START_REF] Wu | The substructuring-based topology optimization for maximizing the first eigenvalue of hierarchical lattice structure[END_REF][START_REF] Hoang | Adaptive concurrent topology optimization of cellular composites for additive manufacturing[END_REF][START_REF] Jia | Multiscale topology optimization for non-uniform microstructures with hybrid cellular automata[END_REF][START_REF] Wang | Practical design optimization of cellular structures for additive manufacturing[END_REF][START_REF] Liu | Two-scale concurrent topology optimization of lattice structures with connectable microstructures[END_REF][START_REF] Xu | Multi-scale topology optimization with shell and interface layers for additive manufacturing[END_REF][START_REF] Wu | Design and optimization of conforming lattice structures[END_REF][START_REF] Da | Design of heterogeneous mesostructures for nonseparated scales and analysis of size effects[END_REF][START_REF] Li | Topology optimization design of multi-scale structures with alterable microstructural length-width ratios[END_REF][START_REF] Li | Topology optimization design of quasi-periodic cellular structures based on erode-dilate operators[END_REF][START_REF] Xia | Optimizing cellular structures through the m-vcut level set method with microstructure mapping and high order cutting[END_REF][START_REF] Kijanski | Two-scale shape optimisation based on numerical homogenisation techniques and variational sensitivity analysis[END_REF][START_REF] Li | Integrated design of cellular composites using a level-set topology optimization method[END_REF][START_REF] Li | Topology optimization for concurrent design of structures with multi-patch microstructures by level sets[END_REF]. TO is identified as the most promising approach to carry out the concurrent topology and material optimisation (from a macroscopic scale perspective) since it allows for a total freedom in the choice of the material properties at the macroscopic scale, which depend upon the topology of the RVE at the lower scale (microscopic or mesoscopic, depending on the problem at hand). Moreover, when the problem is stated in the most general way, the optimisation of the topology of the ACM at multiple scales allows for avoiding the introduction of pre-defined RVE geometries, thus a wider design space can be explored and solutions more efficient than conventional ones can be found. The basic idea behind the multi-scale TO of an ACM is that the material is iteratively removed from the design domain (at each scale) and redistributed in order to minimise a prescribed merit function by satisfying the set of design requirements (DRs). Different TO methods have been proposed in the literature to carry out the multi-scale TO of ACMs, like, for instance, the level-set method [START_REF] Sivapuram | Simultaneous material and structural optimization by multiscale topology optimization[END_REF][START_REF] Liu | Minimum length scale constraints in multi-scale topology optimisation for additive manufacturing[END_REF][START_REF] Gao | Topological shape optimization of 3d micro-structured materials using energy-based homogenization method[END_REF], the solid isotropic material with penalisation (SIMP) approach [START_REF] Collet | Topology optimization for microstructural design under stress constraints[END_REF][START_REF] Cheng | Efficient design optimization of variable-density cellular structures for additive manufacturing: theory and experimental validation[END_REF][START_REF] Zhang | Using strain energy-based prediction of effective elastic properties in topology optimization of material microstructures[END_REF][START_REF] Wang | Multiscale isogeometric topology optimization for lattice materials[END_REF][START_REF] Guest | Optimizing multifunctional materials: Design of microstructures for maximized stiffness and fluid permeability[END_REF] or the bidirectional evolutionary structural optimisation method [START_REF] Xia | Concurrent topology optimization design of material and structure within nonlinear multiscale analysis framework[END_REF]. The above strategies can be applied, simultaneously, at the macroscopic scale of the structure and at the RVE scale, or just at this latter, with the aim of finding the optimised RVE topology showing a prescribed macroscopic behaviour, e.g., prescribed stiffness with the least mass as proposed in [START_REF] Xu | Isogeometric topology optimization for rational design of ultra-lightweight architected materials[END_REF], maximum shear stiffness with a prescribed volume [START_REF] Du | Topological design optimization of lattice structures to maximize shear stiffness[END_REF], or unconventional properties, like negative Poisson's ratio with a prescribed volume [START_REF] Gao | Isogeometric topology optimization for computational design of re-entrant and chiral auxetic composites[END_REF].

As far as TO of ACMs at multiple scales is concerned, many research works are available in the literature: an exhaustive review on this topic can be found in [START_REF] Wu | Topology optimization of multi-scale structures: a review[END_REF]. For instance, Sivapuram et al. [START_REF] Sivapuram | Simultaneous material and structural optimization by multiscale topology optimization[END_REF] make use of both the level-set method and the shape sensitivity to optimise, at the same time, the topology of the ACM at both RVE and structure scales in order to minimise the macroscopic compliance subject to constraint on the volume fraction at each problem scale. The scale transition is ensured via the asymptotic homogenisation method. A different approach is proposed by Wang et al. [START_REF] Wang | Multiscale isogeometric topology optimization for lattice materials[END_REF]: firstly, the RVE topology is chosen from a database of pre-defined architectures; secondly, equivalent material properties at the macroscopic scale are assessed via asymptotic homogenisation method as a function of the RVE relative density; finally, the TO is performed at the macroscopic scale via the SIMP approach to minimise the macroscopic compliance subject to a constraint on the volume fraction. Another strategy, based on the lattice infill Bertolino and Montemurro: Preprint submitted to Elsevier technique, is proposed by Yu et al. [START_REF] Yu | Stress-constrained shell-lattice infill structural optimisation for additive manufacturing[END_REF]. It consists of replacing the fictitious macroscopic heterogeneous material with a graded lattice material of a pre-defined topology characterised by a relative density equal to the local value of the pseudo-density field used to penalise the macroscopic stiffness tensor. To reduce the discrepancy between the results of the macroscopic TO and the structural responses evaluated after performing the lattice infill phase, a surrogate model of the ACM is developed and the penalisation scheme to be used at the macroscopic scale (during TO) is assessed a priori via lattice unit sampling and polynomial curve fitting. Finally, the TO problem is formulated by considering the volume as cost function to be minimised subject to a constraint on the Tsai-Hill failure criterion assessed on the lattice core. Xia et al. [START_REF] Xia | Concurrent topology optimization design of material and structure within nonlinear multiscale analysis framework[END_REF] proposed an approach for the concurrent optimisation of the material properties and of the macroscopic topology of a structure based on the bidirectional evolutionary structural optimisation method and finite element square technique to set the strong coupling between problem scales accounting for the non-linear behaviour of the material at the RVE scale and for stress redistribution at the macroscopic one. Li et al. [START_REF] Li | Topology optimization design of multi-scale structures with alterable microstructural length-width ratios[END_REF] proposed a multi-scale TO of ACMs wherein the RVE characteristic size can change locally over the structure. Zhang et al. [START_REF] Zhang | Topological design of sandwich structures with graded cellular cores by multiscale optimization[END_REF] proposed a method to optimise simultaneously the macroscopic design variables representing the distribution of different lattice materials and the microscopic design variables defining the topologies of different RVEs taking into consideration for the possibility of designing graded microstructures.

In all the aforementioned works, the key step is represented by the homogenisation method, which is needed to set the link between the scales of the problem at hand, i.e., the microscopic (or mesoscopic) scale of the ACM RVE and the macroscopic scale of the structure. Thanks to the homogenisation method, at the macroscopic scale the ACM RVE is modelled as an equivalent homogeneous anisotropic material, allowing, in this way, the concurrent optimisation of the macroscopic topology and material properties (through the TO of the RVE) without increasing the computational effort.

In the majority of the existing works the homogenisation approach is based on the asymptotic homogenisation method: e.g., in [START_REF] Arabnejad | Mechanical properties of lattice materials via asymptotic homogenization and comparison with alternative homogenization methods[END_REF] the equivalent elastic properties obtained through asymptotic homogenisation method are compared to those resulting from alternative homogenisation methods. In [START_REF] Sivapuram | Simultaneous material and structural optimization by multiscale topology optimization[END_REF], the equivalent elastic properties at the macroscopic scale are computed starting from the approximation, via Taylor expansion, of the displacement field function often used in the context of the asymptotic homogenisation method. Similarly, in other works [START_REF] Gao | Isogeometric topology optimization for computational design of re-entrant and chiral auxetic composites[END_REF][START_REF] Wang | Multiscale isogeometric topology optimization for lattice materials[END_REF][START_REF] Deng | Connected morphable components-based multiscale topology optimization[END_REF], the asymptotic homogenisation method is coupled with the isogeometric analysis to carry out the multi-scale TO of ACMs.

A good alternative is represented by the strain energy-based homogenisation method (SEHM), which is a general numerical homogenisation method that has been used in different works [START_REF] Collet | Topology optimization for microstructural design under stress constraints[END_REF][START_REF] Zhang | Using strain energy-based prediction of effective elastic properties in topology optimization of material microstructures[END_REF][START_REF] Guest | Optimizing multifunctional materials: Design of microstructures for maximized stiffness and fluid permeability[END_REF][START_REF] Du | Topological design optimization of lattice structures to maximize shear stiffness[END_REF][START_REF] Gao | Topology optimization of micro-structured materials featured with the specific mechanical properties[END_REF] due to its straightforward numerical implementation and direct coupling with the TO algorithms based on pseudo-density field. The SEHM is based on the equivalence between the strain energy of the heterogeneous ACM RVE and that of the corresponding volume of the equivalent homogenised anisotropic medium. The difference between the variants of the SEHM available in the literature is in the post-processing of the results to assess the macroscopic elasticity tensor of the ACM. Indeed, Bertolino and Montemurro: Preprint submitted to Elsevier one can retrieve either the averaged elements stresses inside the RVE [START_REF] Guest | Optimizing multifunctional materials: Design of microstructures for maximized stiffness and fluid permeability[END_REF][START_REF] Xia | Concurrent topology optimization design of material and structure within nonlinear multiscale analysis framework[END_REF] (in this case the resulting homogenisation scheme is called volume-averaged stress-based SEHM) or the elements strain energy [START_REF] Montemurro | A general multi-scale topology optimisation method for lightweight lattice structures obtained through additive manufacturing technology[END_REF][START_REF] Gao | Topological shape optimization of 3d micro-structured materials using energy-based homogenization method[END_REF].

The goal of the present study is to propose a theoretical/numerical framework for the two-scale TO of ACMs by using two different topological descriptors defined at each characteristic scale of the problem at hand. Each topological descriptor is represented through a pseudo-density field in the context of density-based TO approaches. The first pseudo-density field describes the topology at the scale of the RVE of the ACM (this scale is also referred to as lower scale in the following). The second pseudo-density field describes the topology at the macroscopic scale of the structure (also called upper scale in the following). It is noteworthy that, in this study, the RVE architecture is the same at each point of the design domain at the upper scale. This approach is different from those presented in [START_REF] Li | Integrated design of cellular composites using a level-set topology optimization method[END_REF][START_REF] Li | Topology optimization for concurrent design of structures with multi-patch microstructures by level sets[END_REF], where the optimised topology at the upper scale results from the assembly of RVEs of different topology and (eventually) different material properties.

In this work, the coupling between lower and upper scales is established by means of the SEHM, whose goal is to compute the elasticity tensor of the equivalent homogeneous anisotropic material at the upper scale. Of course, at the upper scale, the structural responses, like compliance, displacements, etc. depend upon the two pseudo-density fields defined at both lower and upper scales. The dependency upon the density field at the upper scale is explicit because it intervenes in the penalisation of the stiffness matrix of the elements composing the finite element (FE) model at the macroscopic scale, whilst the dependency upon the topological descriptor at the lower scale is implicit because it is involved in the calculation of the equivalent elasticity tensor of the material at the upper scale through the homogenisation technique. As widely known in the field of multi-scale analysis of composite materials and structures [START_REF] Izzi | A multi-scale two-level optimisation strategy integrating a global/local modelling approach for composite structures[END_REF][START_REF] Scardaoni | A General Global-Local Modelling Framework for the Deterministic Optimisation of Composite Structures[END_REF][START_REF] Delucia | A stochastic approach for predicting the temperature-dependent elastic properties of cork-based composites[END_REF][START_REF] Montemurro | Least-weight composite plates with unconventional stacking sequences : design, analysis and experiments[END_REF][START_REF] Montemurro | A general B-Spline surfaces theoretical framework for optimisation of variable angle-tow laminates[END_REF][START_REF] Montemurro | On the effective integration of manufacturability constraints within the multi-scale methodology for designing variable angle-tow laminates[END_REF], this kind of coupling is referred to as weak coupling because there is only a one-way relationship among scales (i.e., bottom-up scale transition), instead of a two-way relationship in terms of the effect of the macroscopic displacement, strain and stress fields on the lower scale responses.

Nevertheless, when dealing with the multi-scale TO of ACMs a special attention must be put on the problem formulation in order to integrate pertinent DRs dealing with scale transition, manufacturing aspects, anisotropy and boundary conditions (BCs) type. As it can be inferred from the works cited above, the approaches available in the literature to perform the multi-scale TO of ACMs suffer from some drawbacks. The first issue is related to the hypotheses at the basis of the homogenisation procedure. Of course, regardless of the adopted technique, the accuracy of the homogenisation technique depends upon the satisfaction of the following conditions: a) the ACM is obtained by distributing the RVE according to a periodical (or quasi-periodical) pattern within the macroscopic topology; b) the characteristic lengths of the problem scales must differ (at least) of about two/three order of magnitude to fulfil the well-known scale separation condition. Regarding the latter condition, in the literature, only few works integrate the scale separation condition in the problem formulation [START_REF] Zhang | Scale-related topology optimization of cellular materials and structures[END_REF][START_REF] Alexandersen | Topology optimisation of manufacturable microstructural details without length scale separation using a spectral coarse basis preconditioner[END_REF]. However, the approaches proposed in [START_REF] Zhang | Scale-related topology optimization of cellular materials and structures[END_REF][START_REF] Alexandersen | Topology optimisation of manufacturable microstructural details without length scale separation using a spectral coarse basis preconditioner[END_REF] are substantially different from the one presented in this study because in [START_REF] Zhang | Scale-related topology optimization of cellular materials and structures[END_REF][START_REF] Alexandersen | Topology optimisation of manufacturable microstructural details without length scale separation using a spectral coarse basis preconditioner[END_REF] there is only one pseudo-density field describing the topology of the continuum at the RVE scale. Indeed, the RVE can continuously vary, in terms of topology, between the "layers" composing the design domain at the upper scale, whilst there is no descriptor defining the macroscopic topology of the structure. As a matter of fact, in [START_REF] Zhang | Scale-related topology optimization of cellular materials and structures[END_REF][START_REF] Alexandersen | Topology optimisation of manufacturable microstructural details without length scale separation using a spectral coarse basis preconditioner[END_REF] the macroscopic topology is the result of the assembling of the RVE topology of each macroscopic layer according to predefined periodic patterns.

Conversely, when the TO is performed simultaneously at both the macroscopic scale and the RVE one, it is difficult to ensure scale separation because the characteristic length of the macroscopic topological branches could become of the same order of the characteristic length of the RVE scale (recall that the topology of the continuum at both upper and lower scales is continuously changing during the optimisation process). Therefore, a suitable constraint must be introduced in the problem formulation to fulfil the scale separation condition at each iteration of the TO.

The second issue is related to the integration of suitable manufacturing constraints in the problem formulation. In particular, the minimum length scale requirement should be considered at the RVE scale to ensure that small topological branches could be manufactured by means of the considered process.

Moreover, in none of the aforementioned works, mixed non-zero Neumann-Dirichlet BCs are considered and only in few of them (e.g., in [START_REF] Gao | Concurrent topology optimization of multiscale composite structures in matlab[END_REF]) the influence of the RVE anisotropic behaviour on the optimised topology at the macroscopic scale is taken into account, but a systematic analysis considering the influence of the RVE elastic symmetry groups on the macroscopic topology is still lacking.

All the above aspects are investigated in this work. The proposed multi-scale TO approach makes use of: (a) the SIMP method reformulated in the framework of non-uniform rational basis spline (NURBS) hyper-surfaces [START_REF] Montemurro | A general multi-scale topology optimisation method for lightweight lattice structures obtained through additive manufacturing technology[END_REF][START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | NURBS hyper-surfaces for 3D topology optimization problems[END_REF][START_REF] Costa | Minimum length scale control in a NURBS-based SIMP method[END_REF][START_REF] Rodriguez | Structural displacement requirement in a topology optimization algorithm based on isogeometric entities[END_REF][START_REF] Costa | Eigen-frequencies and harmonic responses in topology optimisation: A CAD-compatible algorithm[END_REF][START_REF] Roiné | Stress-based topology optimization through non-uniform rational basis spline hyper-surfaces[END_REF] and (b) the SEHM to set the link between the problem scales. In the following, the method is referred to as NURBS-based SIMP method. Unlike the classical SIMP approach, the NURBS-based SIMP method separates the pseudo-density field, describing the topology of the continuum at each characteristic scale, from the mesh of the FE model. More precisely, if the dimension of the TO problem is , a NURBS hyper-surface of dimension + 1 is used as topology descriptor. This entity is used to describe the pseudo-density field at each scale, which is projected over the mesh of the FE model to penalise the element stiffness matrix according to the SIMP method. As discussed in [START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | NURBS hyper-surfaces for 3D topology optimization problems[END_REF], the use of NURBS entities to describe the topology makes the CAD reconstruction phase a straightforward task [START_REF] Bertolino | An efficient hybrid optimization strategy for surface reconstruction[END_REF][START_REF] Costa | A general hybrid optimization strategy for curve fitting in the non-uniform rational basis spline framework[END_REF]. More details on the NURBS-based SIMP method are available in [START_REF] Montemurro | A general multi-scale topology optimisation method for lightweight lattice structures obtained through additive manufacturing technology[END_REF][START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | NURBS hyper-surfaces for 3D topology optimization problems[END_REF][START_REF] Costa | Minimum length scale control in a NURBS-based SIMP method[END_REF][START_REF] Rodriguez | Structural displacement requirement in a topology optimization algorithm based on isogeometric entities[END_REF][START_REF] Costa | Eigen-frequencies and harmonic responses in topology optimisation: A CAD-compatible algorithm[END_REF][START_REF] Roiné | Stress-based topology optimization through non-uniform rational basis spline hyper-surfaces[END_REF].

In this paper, the multi-scale TO of ACMs is formulated in the most general case by considering mixed non-zero Neumann-Dirichlet BCs and by integrating multiple DRs defined at different scales in the problem formulation, like the constraint on the volume fraction at both scales, the constraint on the mass fraction at the macroscopic scale, the manufacturing constraint on the minimum member size at the RVE scale and the scale separation condition (imposed through the minimum member size at the macroscopic scale) to guarantee the validity of the homogenisation results.

Moreover, the influence of the RVE elastic symmetry group on the optimised topology at both scales is analysed. The effectiveness of the proposed method is tested on 2D and 3D benchmark problems taken from the literature.

The reminder of the paper is as follows. The fundamentals of the NURBS hyper-surfaces theory and of the SEHM are briefly recalled in Sec. 2. In Sec. 3, the mathematical formulation of the multi-scale TO problem of ACMs and the associated numerical strategy are introduced. The numerical results on 2D and 3D benchmark problems are illustrated and discussed in Sec. 4, whilst Sec. 5 ends the paper with meaningful conclusions and prospects.

Notation.

Theoretical background

In this section, the fundamentals of NURBS hyper-surfaces and of the numerical homogenisation procedure used in this work are briefly recalled for a fruitful understanding of the multi-scale TO problem of ACMs presented in Sec.

3.

NURBS hyper-surfaces

A NURBS hyper-surface is a polynomial-based function, defined as ∶ ℝ → ℝ , where is the dimension of the parametric space, whilst is the dimension of the co-domain. The formula of a NURBS hyper-surface reads:

( 1 , … , ) ∶= 1 ∑ 1 =0 ⋯ ∑ =0 1 … ( 1 , … , ) 1 ,…, , (1)
where 1 … ( 1 , … , ) are the piece-wise rational basis functions, which are related to the standard Bernstein's polynomials , ( ), ( = 1, … , ) by means of the relationship:

1 … ∶= 1 … ∏ =1 , ( ) ∑ 1 1 =0 ⋯ ∑ =0 1 … ∏ =1 , ( ) . ( 2
)
In Eqs. ( 1) and ( 2) , ( 1 , … , ) is a -dimension vector-valued rational function, ∈ [0, 1] is the -th dimensionless coordinate (or parametric coordinate), whilst is stored in the array ( ) ∈ ℝ ( 1 +1)×⋯×( +1) , = 1, … , . The explicit expression of CPs coordinates is:

T 1 ,…, = { (1) 1 ,…, , … , ( ) 1 ,…, }. (3)
The CPs layout is referred to as control hyper-net [START_REF] Piegl | The NURBS Book[END_REF]. The generic CP affects the shape of the NURBS entity by means of its coordinates. The overall number of CPs constituting the hyper-net is:

CP ∶= ∏ =1 ( + 1). (4)
In Eq. ( 2), a weight 1 … is associated to the generic CP. The higher the weight the more the NURBS entity is attracted towards the associated CP. For each parametric direction , the NURBS blending functions , ( ) appearing in Eq. ( 2) can be defined recursively as:

,0 ( ) ∶= ⎧ ⎪ ⎨ ⎪ ⎩ 1, if ( ) ≤ < ( ) +1 , 0, otherwise, (5) , ( ) = - ( ) ( ) + - ( )
, -1 ( ) + ( )

+ +1 - ( ) + +1 - ( ) +1 +1, -1 ( ), = 1, ..., , = 1, … , , (6)
where each constitutive blending function is defined on the knot vector:

T ( ) = {0, … , 0 ⏟⏟⏟ +1 , ( ) +1 , … , ( ) --1 , 1, … , 1 ⏟⏟⏟ +1 }, ( ) ∈ ℝ +1 , = 1, … , , (7)
with:

= + + 1. (8)
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 more details on the topic. Here, only the local support property is recalled because it is exploited in the Multi-scale topology optimisation of anisotropic ACMs context of 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 hyper-surfaces for 3D topology optimization problems[END_REF]:

1 … ( 1 , … , ) ≠ 0, if ( 1 , … , ) ∈ (1) 1 , (1) 1 + 1 +1 × ⋯ × ( ) , ( ) + +1 . (9)
Eq. ( 9) means that each CP (and the respective weight) affects only a precise zone of the parametric space, which is denoted as local support.

The strain energy-based homogenisation method

Since the goal of this study is to concurrently optimise the topology of the continuum both at the lower scale and at the upper one, a homogenisation procedure is needed to correctly set the transition among the scales of the problem at hand. Of course, at the lower 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 is modelled as an equivalent homogeneous anisotropic continuum whose mechanical response is described by a set of equivalent material properties. Therefore, the structural responses evaluated at the macroscopic scale will depend upon both the topological descriptor defined at the RVE scale (through the calculation of the equivalent material properties of the homogenised medium) and the one introduced at the macroscopic scale to describe the structure topology (see Sec. [START_REF] Andreassen | Topology optimization of periodic microstructures for enhanced dynamic properties of viscoelastic composite materials[END_REF] for more details on the definition of the topological variable at each scale).

To correctly assess the macroscopic structural responses and their dependency upon the topological variable introduced at the lower scale a homogenisation technique is, thus, necessary. In this work, the macroscopic elastic behaviour of the RVE is determined using the SEHM of periodic media [START_REF] Barbero | Finite Element Analysis of Composite Materials Using ANSYS®[END_REF]. It is noteworthy that two variants of the SEHM are available in the literature: the first one makes use of the elements averaged stresses to compute the components of the macroscopic elasticity matrix [START_REF] Barbero | Finite Element Analysis of Composite Materials Using ANSYS®[END_REF], whilst the second one [START_REF] Montemurro | A general multi-scale topology optimisation method for lightweight lattice structures obtained through additive manufacturing technology[END_REF] employs the elements strain energy to determine . As rigorously proven in [START_REF] Montemurro | A general multi-scale topology optimisation method for lightweight lattice structures obtained through additive manufacturing technology[END_REF], the SEHM based on elements strain energy requires the least number of FE analyses to compute the components of and the related gradient with respect to the topological variable defined at the RVE scale when coupled with a TO algorithm. Accordingly, only this variant of the SEHM is used in this paper to perform the scale transition and its main features are briefly introduced in the following.

The SEHM technique makes use of the repetitive unit of the periodic structure to evaluate the resulting properties at the macroscopic scale. The basic feature of the SEHM is the assumption that the strain energy of the RVE and the one of the corresponding volume of the homogeneous solid at the macroscopic scale are equal. This homogenisation scheme has proven to be an efficient numerical procedure able to determine the equivalent properties of different 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] Bertolino | Multi-scale shape optimisation of lattice structures: an evolutionary-based approach[END_REF][START_REF] Cappelli | Multi-scale identification of the elastic properties variability for composite materials through a hybrid optimisation strategy[END_REF][START_REF] Delucia | Determination of the effective thermoelastic properties of cork-based agglomerates[END_REF].

To evaluate the elastic response of the ACM at the macroscopic scale, two main hypotheses are considered: a) linear, elastic behaviour for the bulk material of the RVE; b) the buckling of the ACM topological branches is disregarded.

Multi-scale topology optimisation of anisotropic ACMs

To assess the components of the stiffness tensor of the ACM at the macroscopic scale, the RVE is submitted to an average strain field 0 , with , = 1, 2, 3 (tensor notation). The six independent components of the average strain tensor are applied through the following set of periodic BCs on the RVE faces [START_REF] Barbero | Finite Element Analysis of Composite Materials Using ANSYS®[END_REF]:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ (2 1 , 2 , 3 ) - (0, 2 , 3 ) = 2 1 0 1 , (0 ≤ 2 ≤ 2 2 , 0 ≤ 3 ≤ 2 3 ), ( 1 , 2 2 , 3 ) - ( 1 , 0, 3 ) = 2 2 0 2 , (0 ≤ 1 ≤ 2 1 , 0 ≤ 3 ≤ 2 3 ), ( 1 , 2 , 2 3 ) - ( 1 , 2 , 0) = 2 3 0 3 , (0 ≤ 1 ≤ 2 1 , 0 ≤ 2 ≤ 2 2 ), (10)
where = 1, 2, 3.

In the above equations, = 2 is the characteristic length of the 3D lattice RVE along the axis, while is the component of the displacement field along the same axis.

Consider, now, the equilibrium equation of the FE model of the RVE. In the most general case it reads:

̂ ̂ = ̂ ; ̂ , ̂ ∈ ℝ ̂ DOF , ̂ ∈ ℝ ̂ DOF × ̂ DOF , (11)
where ̂ DOF is the overall number of degree of freedoms (DOFs) of the structure before the application of the BCs, while ̂ is the non-reduced (singular) stiffness matrix of the RVE. ̂ is the non-reduced vector of generalised displacements containing both imposed and unknown DOFs of the FE model and ̂ is the non-reduced vector of generalised nodal forces (both known and unknown quantities). The expression of the above vectors and matrix is:

̂ ∶= ⎡ ⎢ ⎢ ⎣ BC T BC ̃ ⎤ ⎥ ⎥ ⎦ , ̂ ∶= ⎧ ⎪ ⎨ ⎪ ⎩ BC ⎫ ⎪ ⎬ ⎪ ⎭ , ̂ ∶= ⎧ ⎪ ⎨ ⎪ ⎩ ⎫ ⎪ ⎬ ⎪ ⎭ , , ∈ ℝ DOF , BC , ∈ ℝ BC , ∈ ℝ DOF × DOF , BC ∈ ℝ DOF × BC , ̃ ∈ ℝ BC × BC . ( 12
)
In Eq. ( 12), DOF is the number of unknown DOFs, while BC represents the number of DOFs where BCs on generalised displacements are applied (of course ̂ DOF = DOF + BC ). and BC are the unknown and known vectors of generalised displacements, respectively. is the vector of generalised external nodal forces, whilst is the vector of generalised nodal reactions where BCs on generalised displacements are imposed. , BC and ̃ are the stiffness matrices of the FE model of the RVE after applying BCs. Inasmuch as the periodic BCs of Eq. [START_REF] Wallin | Nonlinear homogenization for topology optimization[END_REF] are imposed in terms of displacements and no external forces are applied to the FE model of the RVE, i.e., = , the equilibrium problem of Eq. ( 12) is of the Dirichlet's type.

Multi-scale topology optimisation of anisotropic ACMs

By considering elementary (i.e., uni-axial and bi-axial) strain components in Eq. [START_REF] Wallin | Nonlinear homogenization for topology optimization[END_REF], and by imposing the equivalence between the strain energy of the equivalent homogeneous anisotropic continuum and that of the ACM RVE, the components of tensor (Voigt's notation) can be determined as [START_REF] Montemurro | A general multi-scale topology optimisation method for lightweight lattice structures obtained through additive manufacturing technology[END_REF]:

=  0 RVE 0 2 , = 1, … , 6, (13)
=  0 , 0 2 RVE 0 0 - 0 2 0 - 0 2 0 , , = 1, … , 6, ≠ . ( 14
)
Eq. ( 13) is used to assess the terms belonging to the main diagonal of tensor , whilst Eq. ( 14) allows for determining the terms outside the main diagonal.  0 and  0 , 0 represent the work of internal forces, evaluated for uni-axial and bi-axial strain fields, respectively, while RVE = 8 1 2 3 is the volume of the 3D domain. Of course, Eq. ( 13) must be solved before Eq. ( 14): in the most general case, i.e., for a macroscopic anisotropic behaviour, six static analyses must be solved for 3D problems in order to assess the diagonal terms of the stiffness tensor and to deduce the extra diagonal terms via the superposition principle (the interested reader is addressed to [START_REF] Montemurro | A general multi-scale topology optimisation method for lightweight lattice structures obtained through additive manufacturing technology[END_REF] for a deeper insight in the matter).

Multi-scale topology optimisation of architected cellular materials in the NURBS-based

SIMP framework

The main goal of this work consists of determining the optimised topology of the ACM at both RVE scale (or lower scale) and macroscopic scale (or upper scale). This problem can be interpreted as a design problem of a complex anisotropic medium whose goal is the concurrent optimisation of the material properties (performed through the optimisation of the RVE topology at the lower scale) and of the macroscopic topology (i.e., the way wherein the equivalent homogeneous anisotropic material is distributed at the upper scale). Of course, since the RVE topology is the same at all points of the upper scale domain, the properties of the homogenised material are uniform at the structure-level.

Unlike the approaches presented in [START_REF] Andreassen | Topology optimization of periodic microstructures for enhanced dynamic properties of viscoelastic composite materials[END_REF][START_REF] Collet | Topology optimization for microstructural design under stress constraints[END_REF][START_REF] Borovinšek | Designing 2d auxetic structures using multi-objective topology optimization[END_REF][START_REF] Zheng | Evolutionary topology optimization for mechanical metamaterials with auxetic property[END_REF][START_REF] Cheng | Efficient design optimization of variable-density cellular structures for additive manufacturing: theory and experimental validation[END_REF][START_REF] Gao | Isogeometric topology optimization for computational design of re-entrant and chiral auxetic composites[END_REF][START_REF] Li | Design of architected materials for thermoelastic macrostructures using level set method[END_REF][START_REF] Wallin | Nonlinear homogenization for topology optimization[END_REF][START_REF] Wang | Topological shape optimization of microstructural metamaterials using a level set method[END_REF], where the goal is to obtain an RVE topology showing a given (possibly optimised) macroscopic elastic behaviour, in this study, the problem characteristic scales (i.e., the lower scale and the upper one) are weakly coupled in the sense that the DRs calculated at the macroscopic scale (and involved in the problem formulation) depend upon the topological variables defined at both scales. On the one hand, the structural

Bertolino and Montemurro: Preprint submitted to Elsevier

Multi-scale topology optimisation of anisotropic ACMs responses calculated at the lower scale, like the equivalent elastic properties of the homogeneous anisotropic material replacing the RVE at the upper scale, the volume fraction (or the mass fraction) of the RVE depend solely upon the topological descriptor defined at the RVE scale. On the other hand, the structural responses defined at the macroscopic scale, like the compliance, the displacement field, etc. depend upon the topological descriptors defined at both scales.

Of course, the dependency of the macroscopic structural responses upon the topological descriptor defined at the lower scale is implicit because the RVE topology affects the components of the macroscopic elasticity matrix .

The goal of the two-scale TO presented in this study is to minimise the generalised compliance of the continuum at the upper scale (under mixed non-zero BCs of the Neumann-Dirichlet type) by considering the following DRs:

1. A constraint on the volume fraction (or, equivalently, on the mass fraction) imposed on the topological descriptor defined at each scale.

2. A constraint on the minimum thickness that can be fabricated through the selected manufacturing process. This constraint is formulated as a minimum member size constraint at the RVE scale.

3. A constraint on the scale separation condition to ensure the validity of the results of the homogenisation technique in calculating equivalent elastic properties of the material at the upper scale. Since the topology of the RVE and the one of the structure are continuously changing during the optimisation process, introducing this type of constraint reveals of paramount importance to avoid the occurrence of too small topological branches at the upper scale whose size could become of the same order of magnitude of the RVE characteristic length (in such circumstances the results of the homogenisation technique are not correct). The requirement on the scale separation condition is formulated in terms of the minimum member size constraint on the topological descriptor defined at the upper scale. In particular, the smallest thickness of the topological branches occurring at the upper scale must be greater than or equal to a multiple of the characteristic length of the design domain at the lower scale.

It is noteworthy that the notion of scale separation introduced in this work is the one ensuring the validity of the results of the homogenisation process at the upper scale within a continuously changing topology. This concept should not be confused with the notion of scales coupling. The scales are separated only in terms of characteristic lengths, but they are coupled because the structural responses at the upper scale depend upon the geometrical and material parameters defined at the lower scale. In summary, imposing the scale separation conditions does not mean that the scales are uncoupled during the optimisation process.

The main features of the approach are briefly described here only for 3D multi-scale TO problems. The characteristic problem scales (and the relative geometrical features) are illustrated in Fig. 1.

Bertolino and Montemurro: Preprint submitted to Elsevier

 ∶= { T = 1 , 2 , 3 ∈ ℝ 3 ∶ ∈ [0, ], = 1, 2, 3}, = , , (15)
where , is the characteristic length of the domain defined along axis, as shown in Fig. 1. In the SIMP approach the material domain Ω ⊆  at the generic scale is identified by means of the pseudo-density function

( ) ∈ [0, 1]
for ∈  : ( ) = 0 means absence of material, whilst ( ) = 1 implies presence of material.

In the framework of the NURBS-based SIMP method, the topological variable (at each scale) is represented by a NURBS entity. More specifically, a NURBS entity of dimension + 1 is used to describe the pseudo-density field, i.e., the topological variable, of a problem of dimension . Therefore, if a 3D TO problem is considered, a 4D NURBS hyper-surface is needed to describe the part topology [START_REF] Costa | NURBS hyper-surfaces for 3D topology optimization problems[END_REF]. The first three coordinates of the NURBS entity correspond

Bertolino and Montemurro: Preprint submitted to Elsevier to the Cartesian coordinates defining the domain, while the last coordinate corresponds to the pseudo-density field that reads:

( 1 , 2 , 3 ) = 1 ∑ 1 =0 2 ∑ 2 =0 3 ∑ 3 =0 1 2 3 ( 1 , 2 , 1 ) 1 2 3 , = , . (16)
In Eq. ( 16), 1 2 3 is the pseudo-density value at the generic CP, i.e., the fourth coordinate of the vector in Eq. ( 1), while 1 2 3 are the rational basis functions of Eq. ( 2). The dimensionless parameters are obtained as:

= , = 1, 2, 3. (17)
Among the parameters governing the shape of the NURBS entity, only the pseudo-density at CPs and the associated weights are included in the design variable vectors 1 and 2 defined as:

T 1 ∶= 000 , ⋯ , 1 2 3 , T 2 ∶= 000 , ⋯ , 1 2 3 , 1 , 2 ∈ ℝ CP , ( 18
)
accordingly, the number of design variables for each scale is, at most, var = 2 CP , = , .

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 and are not optimised: for more details the reader is addressed to [START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | NURBS hyper-surfaces for 3D topology optimization problems[END_REF].

Objective function

At the lower scale, the stiffness matrix of the FE model of the RVE is penalised according to the following formula:

̂ = ∑ =1 1 , 2 ̂ T ̂ , ( 19
)
where is the fictitious density of Eq. ( 16) computed at the centroid of the generic element and is the total number of elements composing the FE model of the RVE. ̂ ∈ ℝ DOF × ̂ DOF is the connectivity matrix of element (whose number of DOFs is DOF ), whilst ∈ ℝ DOF × DOF is the non-penalised element stiffness matrix expressed in the global reference frame of the model. In Eq. ( 19), ≥ 1 is a suitable parameter that aims at penalising

Multi-scale topology optimisation of anisotropic ACMs

all the meaningless densities between 0 and 1: in agreement with the classic SIMP approach this parameter has been set as = 3.

The pseudo-density function defined at the RVE scale affects also the mechanical performances of the continuum at the macroscopic scale, via the elasticity tensor , as indicated in Eqs. ( 13) and [START_REF] Chougrani | Parts internal structure definition using non-uniform patterned lattice optimization for mass reduction in additive manufacturing[END_REF]. Of course, the elasticity tensor of the equivalent homogeneous anisotropic material, replacing the true RVE architecture at the macroscopic scale, is involved in the definition of the global stiffness matrix of the FE model at the macroscopic scale ̂ . Accordingly, matrix ̂ depends upon the pseudo-density fields defined at both scales as follows:

̂ ∶= ∑ =1 1 , 2 ̂ T ∫ Ω T 1 , 2 Ω ̂ , ( 20
)
where is the number of elements constituting the FE model at the macroscopic scale, is the penalised fictitious density of Eq. ( 16) computed at the centroid of the generic element of the FE model at macroscopic scale, ̂ ∈ ℝ DOF × ̂ DOF is the connectivity matrix of element (whose number of DOFs is DOF ), while ∈ ℝ 6× DOF is the matrix representing the product between the linear differential operator and the shape function matrices of the generic element.

The formulation of the equilibrium problem for a linear static analysis, considering the most general case of nonzero mixed Neumann-Dirichlet BCs reads:

̂ ̂ = ̂ , ̂ , ̂ ∈ ℝ ̂ DOF , ̂ ∈ ℝ ̂ DOF × ̂ DOF , ( 21
)
where ̂ DOF is the overall number of DOFs of the FE model at the macroscopic scale, whilst ̂ and ̂ are the non-reduced vectors of generalised displacements and forces, respectively. Similarly to Eq. ( 12), the above formula can be rewritten as follows:

⎡ ⎢ ⎢ ⎣ BC BC T ̃ ⎤ ⎥ ⎥ ⎦ ⎧ ⎪ ⎨ ⎪ ⎩ BC ⎫ ⎪ ⎬ ⎪ ⎭ = ⎧ ⎪ ⎨ ⎪ ⎩ ⎫ ⎪ ⎬ ⎪ ⎭ , ∈ ℝ DOF , BC , ∈ ℝ BC , ∈ ℝ DOF × DOF , BC ∈ ℝ DOF × BC , ̃ ∈ ℝ BC × BC , ( 22
)
where the physical meaning of the different quantities is the same as the counterparts defined at the lower scale, see Sec. 2.2.

Multi-scale topology optimisation of anisotropic ACMs

Under mixed non-zero Neumann-Dirichlet BCs the macroscopic compliance  is defined as:

 ∶= T -T BC . ( 23
)
Inasmuch as the solution search for the multi-scale TO problem is carried out by means of a suitable deterministic algorithm, the derivation of the formal expression of the gradient of the objective function with respect to the topological variables at each scale (and of the constraint functions too) is needed to speed up the iterations. Such expressions have already been derived in previous works [START_REF] Montemurro | A general multi-scale topology optimisation method for lightweight lattice structures obtained through additive manufacturing technology[END_REF][START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | NURBS hyper-surfaces for 3D topology optimization problems[END_REF] and are reported here below for the sake of completeness. In particular, when differentiating  with respect to ( = 1, 2) one obtains:

 = - ∑ =1 6 ∑ =1 6 ∑ =1 , = 1, 2, = 1, ⋯ , CP , (24) with
= ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 1 RVE 0 2 ∑ ∈  0 , if = , 1 2 RVE 0 0 ∑ ∈  0 , 0 - 0 2 0 - 0 2 0 , if ≠ . ( 25
)
In Eq. ( 24), ( = 1, … , 6) is the -th component of the strain vector (Voigt's notation) of the generic element of the FE model at the macroscopic scale, whilst is its volume. In Eq. ( 25), 0 is the -th elementary strain imposed on the RVE through the periodic BCs of Eq. ( 10) and  is the work of internal forces of the generic element of the FE model of the RVE.

The gradient of the macroscopic compliance with respect to the topological variable at the upper scale reads:

 = - ∑ ∈  , = 1, 2, = 1, ⋯ , CP , ( 26
)
where  is the internal work of the generic element of the FE model at the macroscopic scale. In Eqs. ( 24)-( 26), the linear index ( = , ) has been introduced for the sake of compactness. The relation between and , ( = 1, 2, 3) is:

∶= 1 + 1 + 2 ( 1 + 1) + 3 ( 1 + 1)( 2 + 1), = , . (27)
while reads:

= ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ , if = 1, 2 1 - , if = 2. ( 28
)
The scalar quantity appearing in Eq. ( 28) is the NURBS rational basis function of Eq. ( 2) evaluated at the element centroid.

Constraint functions

The physical responses functions related to the design requirements considered in this study are presented in the following. The first requirement is related to the lightness of the structure. It is formulated either in terms of a constraint on the overall mass of the structure or in terms of a constraint on its overall volume . These quantities read:

= ∑ =1 , (29)
= . ( 30
)
In Eq. ( 30) is the density of the equivalent homogeneous anisotropic medium which is defined as:

∶= RVE ∑ =1 , (31)
where is the density of the bulk material composing the lattice RVE and is the volume of the generic element composing the FE model of the RVE.

By differentiating and one obtains:

= 0, = 1, 2, = 1, ⋯ , CP , = ∑ ∈ , = 1, 2, = 1, ⋯ , CP , (32)
Bertolino and Montemurro:

Preprint submitted to Elsevier = = RVE ∑ ∈ , = 1, 2, = 1, ⋯ , CP , = , = 1, 2, = 1, ⋯ , CP . (33)
Therefore, the lightness requirement can be formulated either as:

1 ( 1 , 2 ) ∶= ref - ≤ 0, (34)
or as

2 ( 1 , 2 , 1 , 2 ) ∶= ref - ≤ 0, (35)
In the above equations, ref and ref are the reference values of the macroscopic volume and mass, respectively, whilst and are the imposed fraction for the macroscopic volume and mass, respectively. A further requirement is considered in terms of the volume fraction of the solid phase at the RVE scale. To this end, the RVE volume must be introduced as follows:

= ∑ =1 , ( 36
)
whose gradient can be expressed as:

= ∑ ∈ , = 1, 2, = 1, ⋯ , CP , = 0, = 1, 2, = 1, ⋯ , CP . (37)
Accordingly, the requirement on the RVE solid phase volume fraction can be expressed as:

3 ( 1 , 2 ) ∶= RVE - ≤ 0, (38)
where is the imposed volume fraction.

As stated in the introduction, two further requirements will be included in the problem formulation: the scale separation requirement and the manufacturing constraint on the minimum allowable dimensions at the RVE scale, which corresponds to the minimum thickness that can be fabricated through the selected manufacturing process. Both

Multi-scale topology optimisation of anisotropic ACMs requirements are formulated as minimum length scale constraints at the respective scales. On the one hand, the scale separation requirement must be introduced to guarantee the presence of a given number of RVEs within the thinner branches of the optimised topology at the macroscopic scale: this ensures the accuracy of the results of the SEHM. On the other hand, the minimum member size constraint at the RVE scale must be introduced for manufacturing purposes.

Accordingly, these requirements are formulated as:

4 ( 1 , 2 ) ∶= 1 - ,min RVE max =1,2,3 ≤ 0, = 1, 2, 3, (39)
5 ( 1 , 2 ) ∶= 1 - ,min MP ≤ 0. (40)
In Eqs. ( 39) and ( 40), ,min ( = , ) is the minimum length scale of the topology at the -th scale, while RVE is the number of RVEs to be foreseen within the thinner branches of the topology at the macroscopic scale in order to ensure the accuracy of the results of the SEHM, max is the characteristic length at the RVE scale and MP is the minimum dimension that can be obtained through the considered manufacturing process.

As discussed in [START_REF] Costa | Minimum length scale control in a NURBS-based SIMP method[END_REF], the main advantage of the NURBS-based SIMP method is in the handling of the geometric constraints imposed on the topological variable at the generic scale. In particular, since the pseudo-density field describing the topology of the continuum, at both lower and upper scales, is described by means of a pure geometric entity, i.e., a NURBS hyper-surface, it is possible to properly set the integer parameters (number of CPs and basis functions degree along each parametric direction) governing its shape to automatically satisfy the minimum length scale requirement, without introducing an explicit optimisation constraint in the problem formulation. Therefore, in this study, the scale separation requirement of Eq. ( 39) and the manufacturing requirement of Eq. ( 40) are controlled by means of this feature.

Moreover, as discussed in [START_REF] Costa | Minimum length scale control in a NURBS-based SIMP method[END_REF], unlike classical density-based topology optimisation methods, in the NURBS-based SIMP method the minimum member size requirement does not depend upon the size of the elements composing the mesh, but only on the integer parameters involved in the definition of the NURBS entity. This means that an eventual mesh refinement has an impact only on the value of the structural responses (displacements, strains, stresses, etc.), but not on the minimum member size of the topology.

Problem formulation

The multi-scale TO problem considered in this study is formulated as a constrained non-linear programming problem (CNLPP) as follows:

min 1 , 2   ref , subject to ∶ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ̂ ̂ = ̂ , 1 ≤ 0 or 2 ≤ 0, 3 ≤ 0, 1 ∈ [ min , max ], 2 ∈ [ min , max ], = , , ∀ = 1, ..., CP . (41)
In Eq. ( 41),  ref is the reference value of the macroscopic compliance, whilst min and max are lower and upper bounds on the pseudo-density at each CP, and min and max are the bounds on the generic weight (the bounds on the design variables are introduced at each characteristic scale). It is noteworthy that the lower bound of the pseudo-density must be strictly positive to prevent any singularity for the solution of the equilibrium problem. The overall number of design variables of problem ( 41) is equal to var = 2 CP + CP .

Results

The proposed methodology is tested on both 2D and 3D problems: all calculations are carried out by means of 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 hyper-surfaces for 3D topology optimization problems[END_REF].

The software, coded in python language, is interfaced with the FE code ANSYS ® to compute the mechanical responses of the structure at each scale. Moreover, the CNLPP of Eq. ( 41) has been solved through 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], whose parameters are listed in Tab. 1.

The design variables bounds of Problem ( 41) are set as follows: min = 10 -3 , max = 1; min = 0.5, max = 10. Regarding the other continuous parameters involved in the NURBS entity definition, the non-trivial knot vectors components in Eq. ( 7) have been evenly distributed in the interval [0, 1], for both 2D and 3D cases.

Furthermore, symmetry constraints on the pseudo-density field describing the RVE topology are imposed in all the numerical tests, except those wherein the effect of anisotropy is investigated: double symmetry for 2D problems (with respect to axes = , = 1, 2) and three planes of symmetry ( = , = 1, 2, 3) for 3D problems, in order to have an optimised topology characterised, at most, by an orthotropic behaviour. Of course, the presence of symmetry axes/planes implies a reduction in the design variables count, at the -th scale, as follows:

var = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∏ =1 , for B -spline entity, 2 ∏ =1 , for NURBS entity, (42)
with = 2 and = 3 for 2D and 3D problems, respectively, and

= ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ + 1 2 , if is odd, ⌊ + 1 2 ⌋ + 1, otherwise, (43)
where ⌊⋯⌋ is the floor operator. Regarding the campaign of numerical tests, the following aspects are investigated:

1. The influence of the geometric entity, i.e., B-spline or NURBS, used to describe the TV at each scale on the optimised solution (only for 2D problems);

2. The influence of mixed BCs applied at the macroscopic scale on the optimised topology at each scale (for 2D and 3D problems);

3. The influence of the minimum member size at the RVE scale (related to the minimum printable dimension) on the optimised solution (for both 2D and 3D problems);

4. The influence of the scale separation constraint on the optimised solution (for both 2D and 3D problems);

5. The influence of the constraint type (mass or volume) imposed at the macroscopic scale on the optimised solution (only for 3D problems).

Bertolino and Montemurro: Preprint submitted to Elsevier

Unless explicitly stated, for all benchmarks, an aluminium alloy, characterised by = 71 GPa, = 0.33, = 2.8 × 10 -6 kg mm 3 , is used as bulk material of the RVE. The reference macroscopic mass of the structure and the reference macroscopic compliance are those characterising the starting solution. Moreover, the reference volume (at the generic scale) is the volume of the overall design domain of dimension .

2D benchmark problems

The design domain of the RVE is shown in Fig. 2 The first 2D benchmark problem (denoted as BK1-2D in the following), taken from [START_REF] Gao | Concurrent topology optimization of multiscale composite structures in matlab[END_REF], is shown in Fig. 3a and deals with a 2D cantilever beam. The geometrical parameters of the design domain at the macroscopic scale are:

1 = 600 mm and 2 = 300 mm. The set of nodes located at 1 = 0 mm is clamped, while a point load

2 = -1 N is applied at ( 1 , 2 ) = ( 1 , 2 ∕2).
A static analysis is conducted on the macroscopic FE model whose mesh is made of = 7200 PLANE182 elements (plane stress hypothesis with unit thickness).

The second 2D benchmark problem (denoted as BK2-2D in the following) is characterised by the same geometrical parameters and the same mesh of BK1-2D. Also in this case the nodes located at 1 = 0 are clamped, while the unit point force is applied at node B, as illustrated in Fig. 3b. The third 2D benchmark problem (denoted as BK3-2D in the following), taken from [START_REF] Niu | A general formulation of structural topology optimization for maximizing structural stiffness[END_REF] and shown in Fig. 4, is characterised by the same geometrical parameters and the same mesh of benchmarks BK1-2D and BK2-2D. BK3-2D is used to assess the influence of non-zero mixed Neumann-Dirichlet BCs on the optimised topology. To this purpose, two sets of BCs are considered, as shown in Fig. 4. Each set of BCs corresponds to a different test case, i.e., BK3-2Da and BK3-2Db. In particular, both configurations share the following BCs:

a M 1 a M 2 A B C D x M 2 F M 2 E x M 1 (a) BK1-2D a M 1 a M 2 x M 2 x M 1 A B C D F M 2 (b) BK2-2D
1 = 2 = 0 at node A, 2 = 0 at node B, 2 = -1 N at node E (which is located at 1 = 1 ∕2, 2 = 0).
The difference between BK3-2Da and BK3-2Db is in the location where a variable displacement 1 mm is imposed: as shown in Fig. 4, in the case of BK3-2Da this displacement is imposed at node B and takes value in the interval [-0.05, 0.05] mm, while in the case of BK3-2Db it is imposed at node C and varies in the range [-0.2, 0.2] mm.

a M 1 a M 2 F M 2 u M 1 x M 2 x M 1 A B C D E (a) BK3-2Da a M 1 a M 2 F M 2 u M 1 x M 2 x M 1 A B C D E (b) BK3-2Db

Sensitivity of the optimised solution to the mass/volume constraint fraction

A campaign of numerical analyses is performed on BK1-2D, with the aim of showing the influence of the RVE volume fraction and of the macroscopic mass fraction on the optimised topology at both scales. The topology of the RVE is enforced to be symmetric with respect to axes = ( = 1, 2) resulting, thus, in a ACM with an orthotropic behaviour at the macroscopic scale. Moreover, as discussed in [START_REF] Gao | Topological shape optimization of 3d micro-structured materials using energy-based homogenization method[END_REF][START_REF] Xia | Concurrent topology optimization design of material and structure within nonlinear multiscale analysis framework[END_REF], the initial guess at the RVE scale is characterised by a central hole, which is generated by setting to zero the value of the central CPs of the NURBS/B-spline entity. In particular, the number of CPs whose pseudo-density is set equal to zero is the result of an iterative process that stops when the constraint on the RVE volume fraction of Eq. ( 36) is met. The CNLPP of Eq. ( 41), for which only constraints

At the RVE scale a B-spline surface with CP = 25 × 25 and = 3, ( = 1, 2) is employed. According to the methodology discussed in [START_REF] Costa | Minimum length scale control in a NURBS-based SIMP method[END_REF], this configuration ensures a minimum member size equal to ,min = 0.13 mm within the design domain and equal to ,min = 0.12 mm at the boundary of the domain. This value is greater than the minimum thickness that can be obtained through the direct metal laser sintering process.

At the macroscopic scale two B-spline surfaces are used: the first one has CP = 29×15 CPs, while the second one has CP = 23 × 14 CPs and both configurations are characterised by blending functions degrees = 3, ( = 1, 2).

Both solutions are characterised by a minimum length scale min = 30 mm within the domain and ,min = 6.5 mm along its boundary, but a different number of CPs is required. In fact, as explained in [START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | Minimum length scale control in a NURBS-based SIMP method[END_REF], in the post-processing phase the 3D pseudo-density descriptor, i.e., the B-spline surface, is cut by a plane placed at a given value of threshold, which is cut = 0.35 for the first surface and cut = 0.45 for the second one (this value depends upon the constraint imposed on the mass fraction as discussed in the following). Of course, for both B-spline surfaces the minimum number of RVEs within the thinnest topological branch is equal to 10.

The optimised solutions are reported in Tabs. 2 and 3 by considering two different values of the mass fraction at the macroscopic scale, i.e., = 0.35 and = 0.45, respectively, and different values of the volume fraction of the RVE. In each table, results are provided in terms of the number of iterations to achieve convergence ( iter ), of the value of reference and optimised compliance as well as of the minimum member size measured at the end of the optimisation process at each scale, i.e., ̂ ,min ( = , ). The following remarks can be inferred from the analysis of these results.

Firstly, inasmuch as the constraint on the mass fraction is imposed at the macroscopic scale, and since this constraint is proportional to the product of the volume fractions at both RVE and macroscopic scales (see Eq. ( 35)), the greater the volume fraction of the RVE the lower the volume fraction of the structure at the macroscopic scale for a given mass fraction.

Secondly, for a given mass fraction at the macroscopic scale, the greater the volume fraction of the RVE the better the performances of the optimised solution in terms of compliance. In particular, for BK1-2D, the optimised solution is characterised by a dense isotropic material filling the whole RVE domain and by the least volume fraction at the macroscopic scale.

Finally, for each solution, the minimum length scale requirement at both the RVE scale and the macroscopic scale is systematically fulfilled thanks to the local support property of the NURBS entities, which establishes an implicit filter according to Eq. ( 9). In particular, the scale separation requirement is satisfied in both the configurations of the B-spline entities. This is in agreement with the methodology presented in [START_REF] Costa | Minimum length scale control in a NURBS-based SIMP method[END_REF] whose aim is to set the integer parameters of the NURBS entity (i.e., number of CPs and basis functions degrees) in order to assess the lowest value of the member size achievable during the optimisation process. Of course, the minimum member size of the topology at the end of the optimisation process can be greater than this value because of the presence of further design requirements.

Sensitivity of the optimised solution to the starting guess at the lower scale

Due to the non-convexity of the CNLPP of Eq. ( 41), another sensitivity analysis is performed on BK1-2D to take into account the effect of the starting guess at the RVE scale on the optimised solution at both scales. This (3) a topology with four holes generated symmetrically with respect to axes = ( = 1, 2). The generation of the holes is carried out according to the iterative process described in the above subsection in order to fulfil the design requirement on the volume fraction of the RVE. Conversely, the initial guess at the macroscopic scale is characterised by a uniform density field whose value is determined in order to met the constraint on the macroscopic mass fraction.

For these analyses, an orthotropic behaviour of the RVE is set by imposing an orthogonal symmetry with respect to axes = ( = 1, 2) to its topology during the optimisation process. Moreover, the volume fraction of the RVE and the mass fraction of the whole structure are set to = 0.6 and = 0.35, respectively.

The optimised topologies are listed in Tabs. 4 and 5, for B-spline and NURBS solutions, respectively. From the analysis of these results some remarks can be drawn.

1. Problem ( 41) is non-convex because, for both B-spline and NURBS solutions, the optimised topologies obtained by considering an initial guess with one or two holes, at the RVE scale, correspond to equivalent local minima (the value of the macroscopic compliance is the same), although the components of the macroscopic elasticity matrix are different (this difference is more evident in the case of NURBS solutions). This behaviour is analogous to the one occurring when dealing with the multi-scale optimisation of composite (or anisotropic) structures [START_REF] Montemurro | An extension of the polar method to the first-order shear deformation theory of laminates[END_REF][START_REF] Montemurro | The polar analysis of the third-order shear deformation theory of laminates[END_REF][START_REF] Catapano | Strength optimisation of variable angle-tow composites through a laminate-level failure criterion[END_REF][START_REF] Montemurro | A general multi-scale two-level optimisation strategy for designing composite stiffened panels[END_REF] wherein different stacking sequences can correspond to the same macroscopic behaviour.

2. The optimised topologies obtained by considering a starting RVE domain with one or two holes show performances better (in terms of macroscopic compliance) than those characterising the optimised topologies obtained when considering an initial RVE domain with four holes (for both B-spline and NURBS solutions).

3. Due to the geometric properties of NURBS entities, minimum length scale requirements are always fulfilled at both scales ensuring, thus, that the optimised topologies can be actually manufactured and the results of the homogenisation process are accurate.

4. For each case, the NURBS solution outperforms the B-spline counterpart either in terms of the objective function or in terms of smoothness of the boundary of the optimised topology at both scales (in agreement with the results presented in [START_REF] Montemurro | A general multi-scale topology optimisation method for lightweight lattice structures obtained through additive manufacturing technology[END_REF][START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | NURBS hyper-surfaces for 3D topology optimization problems[END_REF]). This behaviour is related to the ability of NURBS entities to better approximate quadric hyper-surfaces as widely known in the CAD community.

5.

It is noteworthy that B-spline solutions are characterised by transition regions larger than NURBS solutions. This is, indeed, an expected result in agreement with the results obtained for single-scale TO problems presented in [START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | NURBS hyper-surfaces for 3D topology optimization problems[END_REF]. The reason behind the difference in the quality of the solution (in terms of the extent of the transition regions) is twofold. On the one hand, NURBS entities allow for a freedom greater than B-spline entities in describing the details of the topology. Indeed, this increased freedom is related to the number of design variables involved in the definition of the NURBS entity (i.e., density value at each CP and the associated weight), which is higher than the one characterising the B-spline entity (only density value at each CP). This aspect has already been deeply discussed in previous works on the NURBS-based SIMP method, see [START_REF] Montemurro | A general multi-scale topology optimisation method for lightweight lattice structures obtained through additive manufacturing technology[END_REF][START_REF] Costa | A 2D topology optimisation algorithm in NURBS framework with geometric constraints[END_REF][START_REF] Costa | NURBS hyper-surfaces for 3D topology optimization problems[END_REF]. On the other

RVE starting topology RVE optimised topology

Structure optimised topology Optimisation results

RVE starting topology RVE optimised topology

Structure optimised topology Optimisation results

Sensitivity of the optimised solution to the anisotropy of the RVE and to the macroscopic loads

The effect of the elastic symmetry group of the RVE and of the macroscopic loads is investigated on benchmark problems BK1-2D and BK2-2D. In particular, two configurations of the RVE are considered for both benchmarks.

In the first case, the RVE is forced to show an orthotropic behaviour at the macroscopic scale by imposing two symmetry constraints (with respect to axes = , = 1, 2) to the RVE topological variable.

Multi-scale topology optimisation of anisotropic ACMs

In the second case, no symmetry constraint is imposed on the topological variable of the RVE, which can get an arbitrary asymmetric topology during the optimisation showing, thus, a completely anisotropic behaviour at the macroscopic scale.

The RVE volume fraction and the macroscopic mass fraction are set to = 0.6 and = 0.35, respectively, while the NURBS parameters are chosen as follows: CP = 25 × 25, = 3, ( = 1, 2) at the RVE scale and CP = 29 × 15, = 3, ( = 1, 2) at the macroscopic one. The minimum length scale (within the domain and along its boundary) corresponding to this set of integer parameters of the NURBS surface is the same as the discussed in the previous section (for each scale). For each case, the initial guess at the RVE scale is characterised by a single hole in the centre of the RVE, while the initial topology at the macroscopic scale is uniform.

RVE optimised topology

Macroscopic optimised topology Optimisation results

RVE optimised topology

Structure optimised topology Optimisation results

MPa

The optimised solutions for both BK1-2D and BK2-2D are reported in Tabs. 6 and 7, respectively, for both the orthotropic and the anisotropic cases. As it can be inferred from these results, the influence of the RVE behaviour, in terms of elastic symmetry type, is strictly related to the macroscopic loading conditions.

In particular, for BK1-2D the optimised topology is practically unchanged, at both scales, when passing from an orthotropic RVE to an anisotropic one. In fact, in the latter case the RVE topology converges towards a configuration characterised by (almost) a double orthogonal symmetry and by (almost) the same macroscopic behaviour and performances of the first case.

Conversely, for BK2-2D, the elastic symmetry type of the RVE has a strong influence on the optimised topology at both RVE and macroscopic scales. As it can be seen from Tab. 7, although the topologies of the RVE and of the structure are completely different and the macroscopic elasticity matrix has different components, both solutions show equivalent performances in terms of the macroscopic compliance. This is a further evidence of the strong non-convexity of the CNLPP of Eq. ( 41).

Sensitivity of the optimised topology to mixed non-zero boundary conditions

The influence of mixed non-zero Neaumann-Dirichlet BCs on the optimised topologies at both RVE and macroscopic scales is analysed on benchmark problems BK3-2Da and BK3-2Db. In all the analyses, no symmetry constraint is imposed to the topology at the RVE scale, which can show, eventually, an anisotropic behaviour, thus adapting its elastic response to the BCs imposed at the macroscopic scale. The RVE volume fraction and the macroscopic mass fraction are set to = 0.4 and = 0.2, respectively, while the NURBS parameters are chosen as follows: CP = 25 × 25, = 3 for the lower scale and CP = 31 × 17, = 3, ( = 1, 2) for the macroscopic one. These configurations are characterised by a minimum member size (within the design domain) of ,min = 0.13 mm at the RVE scale and ,min = 30 mm at the macroscopic one. Conversely, the minimum member size at the boundary of the domain is equal to ,min = 0.12 mm and ,min = 6.5 mm at RVE scale and macroscopic scale, respectively. For each analysis, the starting guess at the lower scale is characterised by a central hole, while at the upper scale the initial topology is obtained by considering an uniform pseudo-density field. As in the other 2D test cases, the initial topologies satisfy the constraints on the RVE volume fraction and on the macroscopic mass fraction.

The influence of mixed BCs on the optimised topologies at both scales is firstly investigated on BK3-2Da: the applied horizontal displacement 1 can take values in the set [-0.05, -0.025, 0, +0.025, +0.05] mm. Unlike the previous analyses, the isotropic bulk material of the RVE is characterised by = 2525 MPa, = 0.25 and = 2.8 × 10 -6 kgmm -3 . The optimised solutions are listed in Tab. 8. Like the single-scale optimised topologies obtained by using an isotropic material presented in [START_REF] Niu | A general formulation of structural topology optimization for maximizing structural stiffness[END_REF], the optimised solutions at both the RVE scale and the macroscopic scale change when varying the entity of the applied displacement. This means that the result of the Optimised topology (isotropy) Optimised topology (transverse isotropy) The following remarks can be inferred from the analysis of these results:

Firstly, as expected, the optimised topology (at both scales) is different when considering either the constraint on the macroscopic volume fraction 1 or the one on the macroscopic mass fraction 2 . In particular, the solution of the CNLPP of Eq. ( 41) when integrating a constraint on the macroscopic mass fraction is characterised by a value of the objective function lower than the one characterising the solution obtained by considering constraint 1 . This is mainly because a macroscopic mass fraction = 0.2 corresponds to a value of the macroscopic volume fraction higher than = 0.2, thus the optimised topology in this case is stiffer than the solution obtained by considering constraint 1 . Therefore, the designer should pay attention to the problem formulation because the constraint on the macroscopic mass fraction 2 influences the topology at both scales (unlike the constraint on the macroscopic volume fraction 1 ).

Secondly, regardless of the problem formulation, the minimum length scale requirement is systematically fulfilled at both scales (for both test cases BK1-3D and BK2-3D). In particular, when constraint 2 is included in the CNLPP formulation, the optimised topology at the macroscopic scale is characterised by a more compact distribution of the material.

Bertolino and Montemurro: Preprint submitted to Elsevier

Sensitivity of the optimised solution to the scale separation requirement

The effect of the implicit scale separation constraint on the optimised topologies at both scales is analysed on BK1-3D and BK2-3D. In particular, two configurations of the macroscopic NURBS hyper-surface parameters are considered, for both benchmarks, while the same set of integer parameters is used for the topology at the RVE scale, i.e., CP = 11 × 11 × 11 and = 3, ( = 1, 2, 3), to ensure a minimum member size within the domain and along its boundary equal to ,min = 0.13 mm and ,min = 0.12 mm, respectively. Therefore, for each benchmark structure, two cases are considered.

In the first case the scale separation requirement is disregarded and the NURBS integer parameters are chosen as As shown in Tabs. 13 and 14, the implicit constraint on the minimum member size of the macroscopic topology has an impact on the optimised topologies at both scales. In terms of objective function value, for BK1-3D the optimised solution obtained without considering the scale separation requirement outperforms the optimised solution satisfying this constraint (see Tab. 13), the difference between the compliance of the two configurations being 7%.

Conversely, for BK2-3D an unexpected result is found: the optimised topology obtained by considering the scale separation requirement outperforms the one wherein it is disregarded (the relative difference between the compliance of the two solution is equal to 28%). This result is mainly due to the optimised topology at the RVE scale which is strongly different between the two optimised solutions reported in Tab. 14.

Sensitivity of the optimised topology to mixed non-zero boundary conditions

The influence of the mixed non-zero Neumann-Dirichlet BCs on the topologies at both scales is analysed on BK3-3D. In all the analyses, the topological variable at the RVE scale is characterised by the following NURBS Three values of displacement 1 are considered: 1 = {-1, 0, 1} mm. The optimised solutions are listed in Tab. 15. As expected, the RVE optimised topology is characterised by an orthotropic behaviour and, similarly to the results of the test case BK3-2Da, the optimised topology at both scales change when varying the entity of the applied displacement. This means that the result of the concurrent optimisation of the elastic symmetry of the material (related to the RVE topology) and of the macroscopic topology strongly depends upon the intensity of the applied displacement.

Finally, the same remarks provided for benchmark problem BK3-2Da can be repeated here too.

Bertolino and Montemurro: Preprint submitted to Elsevier

Conclusions

In this work, an approach for the concurrent optimisation of the material properties and the structure topology of ACMs has been presented. More precisely, the problem of designing of ACMs is formulated as a two-scale TO problem. In this context, two topological variables have been introduced: the first one is defined at the RVEs scale of the ACM (which can be the microscopic scale or the mesoscopic one, depending on the problem characteristic lengths), while the second one is defined at the macroscopic scale. A weak coupling among lower and upper scales has been considered in the problem formulation. In fact, the structural responses at the upper scale depend both on the topological descriptor defined at this scale and to the one introduced at the lower scale through the calculation of the equivalent elasticity matrix of the fictitious homogeneous material, which replaces the RVE at the upper scale. In this context, the topological variable at each scale is represented by the pseudo-density field of the SIMP method, which is described by means of a NURBS hyper-surface and the scale transition is correctly set through the SEHM based on elements strain energy. NURBS entities are efficiently exploited to derive the relationship among topological variables defined at different scales and also to determine the analytical expression of the gradient of the physical responses involved at all scales.

Some features of the proposed methodology need to be highlighted.

Firstly, some advantages of the NURBS formalism can be clearly identified: (a) 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, unlike the classical SIMP method;

(b) unlike the classical SIMP approach, there is no need to define a further filter zone, since the NURBS local support property establishes an implicit relationship among contiguous mesh elements; (c) when compared to the classical SIMP approach, the number of design variables is reduced; (d) 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 (e.g., on the local curvature of the boundary, on the local direction of the tangent vector, maximum member size, etc.) in the problem formulation and the CAD reconstruction phase of the boundary of the optimised topology become easy tasks. Of course, as in the case of the results of single-scale problems found through 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 hyper-surfaces for 3D topology optimization problems[END_REF], the optimised topology depends upon the NURBS integer parameters, i.e., number of control points and degrees of Bernstein's polynomials, which have a direct impact on the size of the local support of the blending functions.

Secondly, it has been shown that, depending on the problem formulation, the values of the macroscopic volume/mass fraction and of the RVE volume fraction, as expected, have a strong influence at each scale, thus a particular care must be put in the choice of the combination of such constraints to avoid "weird" topologies.

Thirdly, the influence of the RVE elastic symmetry group on the final topology, at both scales, has been considered.

Results highlight that, depending on the problem at hand, forcing a pre-defined elastic symmetry (e.g., orthotropy, transverse isotropy, etc.) could not reveal a good choice because completely anisotropic solutions characterised by equivalent or better performances could exist. In analogy with the multi-scale optimisation of composite materials and structures, this is due to the non-convexity of the optimisation problem and to the existence of equivalent local minima characterised by very different properties. Fourthly, the influence of the initial guess on the optimised topologies at both scales has been taken into account. As expected, the obtained results represent a confirmation of a well-established result: under the same working hypotheses and BCs, multi-scale TO problems of ACMs are non-convex; therefore, the choice of the starting guess at the RVE scale is of paramount importance and has a strong influence on the RVE optimised topology.

Furthermore, the influence of mixed non-zero Neumann-Dirichlet BCs on the optimised solution has been investigated.

The obtained results highlight that, to satisfy the requirements of the problem at hand and to withstand the applied loads, both the RVE topology and the macroscopic one evolve towards a configuration optimising the macroscopic elastic response of the continuum, as well as the material distribution at the upper scale.

Finally, the effect of the minimum-length scale requirement at each problem scale on the optimised solutions is correctly taken into account, without the need of introducing an explicit optimisation constraint. Indeed, this requirement can be easily fulfilled by properly setting the integer parameters of the NURBS entity representing the topological variable at each scale. In particular, the technological constraint on the minimum manufacturable dimension is handled by controlling the value of the NURBS blending functions degrees and number of CPs at the lower scale, whilst the scale separation requirement is ensured by controlling the same quantities related to the NURBS entity representing the topological variable at the upper scale.

Regarding the prospects on this topic, several challenges still need to be faced. Firstly, a set of pertinent technological requirements, related to the additive manufacturing process, should be integrated into the multi-scale TO problem formulation, especially in terms of the RVE geometrical features, e.g., overhang angle, admissible curvature, etc. Secondly, in order to design ACMs against failure, a suitable homogenisation procedure with strong coupling between scales should be developed and integrated in the optimisation process. Lastly, suitable failure criteria at upper and lower scales should be derived and integrated in the multi-scale TO problem formulation for the homogeneous anisotropic material at the macroscopic scale and for the bulk material at the RVE scale.

1

  ,…, ∈ ℝ are the control points (CPs) coordinates , while ∈ ℕ and ∈ ℕ ∪ 0 ( = 1, … , ) are the number of CPs and the basis functions degree, respectively, along the parametric Multi-scale topology optimisation of anisotropic ACMs direction. The -th CP coordinate ( ) 1 ,…,

Figure 1 :

 1 Figure 1: Characteristic scales of the two-scale topology optimisation problem: (a) the macroscopic (upper) scale of the structure and (b) the lower (i.e., mesoscopic or microscopic) scale of the representative volume element of the architected cellular material

Figure 2 :

 2 Figure 2: Finite element model of the representative volume element of the architected cellular material with its characteristic size for 2D problems.

Figure 3 :

 3 Figure 3: Finite element model, geometrical parameters and applied boundary conditions of benchmark problems (a) BK1-2D and (b) BK2-2D.

Figure 4 :

 4 Figure 4: Finite element model, geometrical parameters and applied boundary conditions of benchmark problems (a) BK3-2Da and (b) BK3-2Db.

Bertolino and Montemurro :

 Montemurro Preprint submitted to Elsevier Table 2 BK1-2D: optimised topologies at both scales, when the mass fraction is set as = 0.35. Solutions obtained by employing B-spline entities, with CP = 25 × 25 and CP = 29 × 15 and = 3, ( = 1, 2, = , ).  ref = 0.0126 Nmm;  opt = 0.0030 Nmm; ̂ ,min = 0.38 mm; ̂ ,min = 45.38 mm; ̂ ,min = 0.44 mm;  ref = 0.0136 Nmm;  opt = 0.0026 Nmm; ̂ ,min = 0.62 mm; ̂ ,min = 34.03 mm; ̂ ,min = 0.56 mm; ̂ ,min = 32.4 mm

  analysis has been conducted by considering both B-spline and NURBS entities (at each problem scale) characterised Bertolino and Montemurro: Preprint submitted to Elsevier Table 3 BK1-2D: optimised topologies at both scales, when the mass fraction is set as = 0.45. Solutions obtained by employing B-spline entities, with CP = 25 × 25 and CP = 23 × 14 and = 3, ( = 1, 2, = , ).  ref = 0.0059 Nmm;  opt = 0.0022 Nmm; ̂ ,min = 0.38 mm; ̂ ,min = 65.54 mm; ̂ ,min = 0.53 mm;  ref = 0.0063 Nmm;  opt = 0.0018 Nmm; ̂ ,min = 0.21 mm; ̂ ,min = 39.67 mm; ̂ ,min = 0.51 mm; ̂ ,min = 57.82 mm by CP = 25 × 25, CP = 29 × 15 and = 3 ( = , , = 1, 2). These configurations are characterised by a minimum member size equal to ,min = 0.13 mm ,min = 30 mm, within the domain, while its value at the domain boundary is ,min = 0.12 mm and ,min = 6.5 mm at RVE scale and macroscopic scale, respectively. For each case, three starting guesses are considered at the RVE scale: (1) a topology with a central hole; (2) a topology with two holes generated symmetrically with respect to axis 1 = 1 ;

Bertolino and Montemurro :

 Montemurro Preprint submitted to Elsevier Table 4 BK1-2D: optimised topologies at both scales for dierent starting topologies at the RVE scale. Solutions obtained by employing B-spline entities, with CP = 25 × 25 and CP = 29 × 15 and = 3, ( = 1, 2, = , ).

  ;  ref = 0.0628 Nmm;  opt = 0.0031 Nmm; ̂ ,min = 0.38 mm; ̂ ,min = 50.42 mm; ̂ ,min = 0.57 mm;  ref = 0.018 Nmm;  opt = 0.0035 Nmm; ̂ ,min = 0.17 mm; ̂ ,min = 41.48 mm; ̂ ,min = 0.28 mm; using B-spline entities as topological descriptors at both scales, it was observed that, very often, the GCMMA algorithm stops because the criterion on the tolerance on the (normalised) objective function is satisfied, which could explain the presence of the greater transition region (i.e., the algorithm converges towards a point falling in the neighbourhood of the local minimum, but not exactly towards the local minimum). Bertolino and Montemurro: Preprint submitted to Elsevier Table 5 BK1-2D: optimised topologies at both scales for dierent starting topologies at the RVE scale. Solutions obtained by employing NURBS entities, with CP = 25 × 25 and CP = 29 × 15 and = 3, ( = 1, 2, = , ).

  ;  ref = 0.0628 Nmm;  opt = 0.0025 Nmm; ̂ ,min = 0.34 mm; ̂ ,min = 39.04 mm; ̂ ,min = 0.7 mm; ̂ ,min = 47.5 mm;  ref = 0.018 Nmm;  opt = 0.0028 Nmm; ̂ ,min = 0.2 mm; ̂ ,min = 39.04 mm; ̂ ,min = 0.3 mm;

  : optimised topologies at both scales by considering an RVE with orthotropic and anisotropic behaviour. Solutions obtained by employing NURBS entities, with CP = 25 × 25 and CP = 29 × 15 and = 3, ( = 1, 2, = , ).

  ;  ref = 0.0126 Nmm;  opt = 0.0025 Nmm; ̂ ,min = 0.35 mm; ̂ ,min = 45.20 mm; ̂ ,min = 0.36 mm; ̂ ,min = 47.60 mm;  ref = 0.0127 Nmm;  opt = 0.0025 Nmm; ̂ ,min = 0.35 mm; ̂ ,min = 45.20 mm; ̂ ,min = 0.36 mm; ̂ ,min = 47.60 mm; : optimised topologies at both scales by considering an RVE with orthotropic and anisotropic behaviour. Solutions obtained by employing NURBS entities, with CP = 25 × 25 and CP = 29 × 15 and = 3, ( = 1, 2, = , ).

  ;  ref = 0.0142 Nmm;  opt = 0.0031 Nmm; ̂ ,min = 0.36 mm; ̂ ,min = 28.00 mm; ̂ ,min = 0.46 mm; ̂ ,min = 18.00 mm;  ref = 0.0145 Nmm;  opt = 0.0030 Nmm; ̂ ,min = 0.27 mm; ̂ ,min = 25.00 mm; ̂ ,min = 0.33 mm; ̂ ,min = 35.71 mm;

  : optimised topologies at both scales under mixed non-zero BCs. Solutions obtained by employing NURBS entities, with CP = 25 × 25 and CP = 31 × 17 and = 3, ( = 1, 2, = , ). ;  ref = 0.5742 Nmm;  opt = 0.0230 Nmm; ̂ ,min = 0.22 mm; ̂ ,min = 148.68 mm; ̂ ,min = 0.33 mm;  ref = 0.5860 Nmm;  opt = 0.0408 Nmm; ̂ ,min = 0.2143 mm; ̂ ,min = 39.26 mm; ̂ ,min = 0.3513 mm;  ref = 0.5919 Nmm;  opt = 0.0484 Nmm; ̂ ,min = 0.28 mm; ̂ ,min = 39.02 mm; ̂ ,min = 0.19 mm; ̂ ,min = 23.23 mmTable 9BK3-2Da: optimised topologies at the macroscopic scale under mixed non-zero BCs when the bulk material is considered isotropic or transversely isotropic. Solutions obtained by employing NURBS entities, with CP = 25×25 and CP = 31×17 and = 3, ( = 1, 2, = , ).

Figure 5 :

 5 Figure 5: BK3-2Da: Macroscopic compliance vs. applied displacement for optimised topologies resulting from two-scale topology optimisation. Solutions obtained by employing NURBS entities, with CP = 25 × 25 and CP = 31 × 17 and = 3, ( = 1, 2, = , ).

Figure 6 : 1 a M 3 a M 2 x M 1 x M 2 x M 3 AFigure 9 :

 61321239 Figure 6: BK3-2Da: Macroscopic compliance vs. applied displacement for optimised topologies, resulting from single-scale topology optimisation, made of a) isotropic material and b) transversely isotropic material. .Solutions obtained by employing NURBS entities, with CP = 25 × 25 and CP = 31 × 17 and = 3, ( = 1, 2, = , ).

4 ̂

 4 ,min ( ( )) ≥ 25 mm iter = 157;  ref = 0.4030 Nmm;  opt = 0.0332 Nmm; ̂ ,min = 0.202 mm; ̂ ,min = 45.51 mm; ̂ ,min = 0.2 mm; ̂ ,min = 17.5 mm For each analysis, the starting guess at the lower scale is characterised by a central hole, while at the upper scale the initial topology is obtained by considering an uniform pseudo-density field. Of course, the initial topologies satisfy the constraints on the RVE volume fraction and on the macroscopic mass/volume fraction. The optimised solutions for BK1-3D and BK2-3D are reported in Tabs. 11 and 12, respectively, by considering the macroscopic volume constraint (with volume fraction = 0.2) and the macroscopic mass one (with mass fraction = 0.2).

  one, the scale separation requirement at the macroscopic scale is taken into account and the NURBS parameters are set as follows: CP = 27 × 15 × 10 and = 3, ( = 1, 2, 3) for BK1-3D; while CP = 15 × 12 × 15 and = 3, ( = 1, 2, 3) for BK3-3D. This configuration is characterised by a minimum member size inside the design domain equal to ,min = 25 mm, which guarantees a minimum number of 10 RVEs along each axes within the thinnest topological branch. Moreover, the minimum member size at the boundary of the domain is ,min = 5.8 mm. The RVE volume fraction and the macroscopic volume fraction are set to = 0.6 and = = 0.2, respectively. Also in this case, the starting guess at the RVE scale is characterised by a central hole, while at the upper Bertolino and Montemurro: Preprint submitted to Elsevier Table 13 BK1-3D: inuence of the scales separation requirement on the optimised topology (at both scales). Solutions obtained by employing NURBS entities, with CP = 11 × 11 × 11, CP = 35 × 25 × 15 (rst case) or CP = 27 × 15 × 10 (second case) and = 3, ( = 1, 2, 3, = , )  ref = 0.4030 Nmm;  opt = 0.0310 Nmm; ̂ ,min = 0.07 mm; ̂ ,min = 23.8 mm; ̂ ,min = 0.07 mm; ̂ ,min = 29 mm ref ≤ 0.2, ref ≤ 0.4, ̂ ,min ( ( )) ≥ 25 mm iter = 157;  ref = 0.4030 Nmm;  opt = 0.0332 Nmm; ̂ ,min = 0.19 mm; ̂ ,min = 31.16 mm; ̂ ,min = 0.2 mm; ̂ ,min = 20.2 mm scale the initial topology is obtained by considering an uniform pseudo-density field; the initial topologies are set such that the constraints on the RVE volume fraction and on the macroscopic volume fraction are met.

parameters CP = 11 ×

 11 11 × 11 and = 3, ( = 1, 2, 3); while the parameters characterising the topological variable at the macroscopic scale are chosen as follows: CP = 27 × 15 × 10 and = 3, ( = 1, 2, 3). These configurations ensure a minimum member size within the design domain equal to ,min = 0.13 mm at the RVE scale and ,min = 25 mm at the structure scale. On the other hand, the minimum member size at the boundary of the design domain is ,min = 0.12 mm and ,min = 5.8 mm at the RVE and structure scales, respectively. The Table 14 BK2-3D: inuence of the scales separation requirement on the optimised topology (at both scales). Solutions obtained by employing NURBS entities, with CP = 11 × 11 × 11, CP = 23 × 17 × 23 (rst case) or CP = 15 × 12 × 15 (second case) and = 3, ( = 1, 2, 3, = , )  ref = 172.71 Nmm;  opt = 0.2746 Nmm; ̂ ,min = 0.232 mm; ̂ ,min = 34.65 mm; ̂ ,min = 0.232 mm; ̂ ,min = 26.5 mm ref ≤ 0.2, ref ≤ 0.4, ̂ ,min ( ( )) ≥ 25 mm iter = 466;  ref = 172.71 Nmm;  opt = 0.2142 Nmm; ̂ ,min = 0.2 mm; ̂ ,min = 43.14 mm; ̂ ,min = 0.2 mm; ̂ ,min = 28.7 mm RVE volume fraction and the macroscopic mass fraction are = 0.6 and = 0.35, respectively. Moreover, the starting topology at the RVE scale is characterised by a central hole, while at the upper scale it is obtained by considering an uniform pseudo-density field; the initial solutions satisfy the constraints on the RVE volume fraction and on the macroscopic volume fraction. The material properties used in this examples are = 2525 MPa, = 0.25 and = 2.8 × 10 -6 kgmm -3 .

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

Table 1

 1

	GCMMA algorithm parameters
	Parameter	Value
	move	0.1
	albefa	0.1
	Stop Criterion	Value
	Maximum n. of function evaluations	10000
	Maximum n. of iterations	20 × var
	Tolerance on objective function	10 -6
	Tolerance on constraints	10 -6
	Tolerance on input variables change	10 -6
	Tolerance on Karush Kuhn Tucker norm	10 -6

Table 11 BK1

 11 -3D: optimised topologies at both scales for dierent macroscopic constraints. Solutions obtained by employing NURBS entities, with CP = 11 × 11 × 11, CP = 27 × 15 × 10 and = 3, ( = 1, 2, 3, = , ).

Table 12 BK2

 12 -3D: optimised topologies at both scales for dierent macroscopic constraints. Solutions obtained by employing NURBS entities, with CP = 11 × 11 × 11, CP = 15 × 12 × 15 and = 3, ( = 1, 2, 3, = , )

	Constraints imposed	RVE optimised topology	Structure optimised topology	Optimisation results
					iter = 466;
					 ref = 172.71 Nmm;
		ref	≤ 0.2		 opt = 0.2142 Nmm;
		ref	≤ 0.4		̂	̂ ,min = 43.14 mm; ,min = 0.2 mm;
	̂	,min ( ( )) ≥ 25 mm		̂	,min = 0.2 mm;
					̂	,min = 28.7 mm

ref ≤ 0.2 ref ≤ 0.4 ̂ ,min ( ( )) ≥ 25 mm iter = 125;  ref = 0.4392 Nmm;  opt = 0.1061 Nmm; ̂ ,min = 0.17 mm; ̂ ,min = 92.3 mm; ̂ ,min = 0.17 mm; ̂ ,min = 70 mm

Table 15 BK3

 15 -3D: optimised topologies at both scales under mixed non-zero BCs. Solutions obtained by employing NURBS hypersurfaces, with

and

are included in the problem formulation, is solved by considering the following combination of blending functions degrees and CPs numbers at each scale.

Acknowledgements

This work benefited from the support of the project COFFA ANR-17-CE10-0008 of the French National Research Agency (ANR).

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.

concurrent optimisation of the elastic symmetry of the material (related to the RVE topology) and of the macroscopic topology strongly depends upon the intensity of the applied displacement.

A similar analysis has been conducted, by considering a single-scale optimisation and by modelling the macroscopic domain with two different bulk materials characterised by the following mechanical properties: an isotropic material with = 2525 MPa, = 0.25 and = 2.8 × 10 -6 kgmm -3 and a transversely isotropic material with The trend of the macroscopic compliance vs. the applied displacement for the optimised topologies obtained in the case of double-scale TO optimisation and single-scale TO (for the isotropic material and the transversely isotropic material) for BK3-2Da are illustrated in Figs. 5 and6. As it can be inferred from this results, the curve is monotonic (in the considered interval wherein 1 takes value) in the case of the multi-scale TO problem, whilst this is not the case for the single-scale TO problem (regardless the elastic symmetry of the material at the macroscopic scale).

Further analyses have been conducted on BK3-2Db, where the applied horizontal displacement 1 can take values in the set [-0.2, -0.1, -0.025, 0, +0.1, +0.2] mm. The properties of the bulk material at the RVE scale are the same as in the case of BK3-2Da. Also in this case the topology at both scales evolves when changing the applied displacement value 1 . Tab. 10 reports the optimised topologies for BK3-2Da, whilst the trend of the macroscopic compliance vs. the applied displacement is illustrated in Fig. 7. It is noteworthy that, except the case 1 = 0, the optimised topology at both the RVE scale and the macroscopic one is completely asymmetric (thus the optimised solution is characterised by an RVE with an anisotropic behaviour). Moreover, unlike the optimised topologies resulting from the multi-scale TO in the case BK3-2Da, the trend of the macroscopic compliance vs. the applied displacement is non-monotonic.

Bertolino and Montemurro: Preprint submitted to Elsevier

3D benchmark problems

The design domain of the RVE, shown in Fig. 8, is a cube of size = 2 1 = 2 2 = 2 3 = 2.5 mm. Six linear static analysis are performed on the FE model of the RVE to assess the mechanical properties of the equivalent homogeneous anisotropic material via the SEHM discussed in Sec. 2.2. For each analysis, the periodic BCs of Eq. [START_REF] Wallin | Nonlinear homogenization for topology optimization[END_REF] have been applied through constraint equations among homologous nodes belonging to the opposite faces of the RVE, by considering elementary unit strains. The FE model of the RVE, which remains unchanged for all the 3D benchmark problems, is implemented in the Ansys automatic parametric design language environment with a mesh of = 8000

SOLID185 elements (8 nodes, 3 DOFs per node). For all the 3D benchmark problems, the topology of the RVE is forced to be symmetric with respect to planes = , = 1, 2, 3) showing, thus, an orthotropic behaviour at the macroscopic scale.

The first 3D benchmark problem (denoted as BK1-3D in the following), taken from [START_REF] Gao | Concurrent topology optimization of multiscale composite structures in matlab[END_REF] and shown in Fig. 9a, is the well-known 3D cantilever beam. The geometrical parameters of the design domain are: 1 = 500 mm, 2 = 250 mm and 3 = 150 mm. The nodes belonging to the plane 1 = 0 are clamped, while a vertical negative force 2 = -10 N is applied to the nodes placed at ( 1 , 2 , 3 ) = ( 1 , 2 ∕2 ± size , 3 ∕2 ± size ), with size = 10 mm. A static analysis is conducted on the macroscopic FE model, discretised through = 18750 SOLID185 elements.

The second 3D benchmark problem (denoted as BK2-3D in the following), taken from [START_REF] Gao | Concurrent topology optimization of multiscale composite structures in matlab[END_REF] and shown in Fig. 9c, is a 3D supported structure of size mesh made of = 12500 SOLID185 elements. As shown in Fig. 9c, the domain is submitted to the following BCs: the nodes E, F, G, H are clamped, while a negative force 2 = -10 N is applied to the nodes located at

3 ) = ( 1 ∕2 ± size , 2 , 3 ∕2 ± size ). The third 3D benchmark problem (denoted as BK3-3D in the following) is the 3D Messerschmitt Bölkow Blohm beam and it is characterised by the same geometrical dimensions and mesh of BK1-3D. The nodes located at 1 ∈ [0, 3 size ], 2 = 0 and 3 ∈ [0, 3 ] are connected rigidly, via MPC184 beam elements, to a master node located at ( 1 , 2 , 3 ) = (0, 0, 3 ∕2), which is clamped. Moreover, nodes placed at 1 ∈ [ 1 -3 size , 1 ], 2 = 0 and 3 ∈ [0, 3 ] are connected to the master node located at ( 1 , 2 , 3 ) = ( 1 , 0, 3 ∕2) via MPC184 elements. To this latter a displacement 1 is applied. In addition, the negative force 2 = -10 N is applied to the master node located at ( 1 , 2 , 3 ) = ( 1 ∕2, 2 , 3 ∕2), which is rigidly connected via MPC184 to the nodes located at 1 = 1 ∕2 ± size , 2 = 2 and 3 ∈ [0, 3 ], as illustrated in Fig. 9b.

Sensitivity of the optimised solution to the constraint type imposed at the macroscale

Numerical tests are performed on BK1-3D and BK2-3D, with the objective of investigating the influence of the constraint nature (i.e., macroscopic volume or mass) on the final topology at both scales. The CNLPP of Eq. ( 41), for which 1 or 2 and 3 are included in the problem formulation, is solved by considering at the RVE scale a NURBS hyper-surface characterised by CP = 11 × 11 × 11 and = 3, ( = 1, 2, 3), for which a minimum member size within the domain and along its boundary are ,min = 0.13 mm and ,min = 0.12 mm, respectively.

Bertolino and Montemurro: Preprint submitted to Elsevier