Enhanced Cellular Materials through Multiscale, Variable-Section Inner Designs: Mechanical Attributes and Neural Network Modeling

In the current work, the mechanical response of multiscale cellular materials with hollow variable-section inner elements is analyzed, combining experimental, numerical and machine learning techniques. At first, the effect of multiscale designs on the macroscale material attributes is quantified as a function of their inner structure. To that scope, analytical, closed-form expressions for the axial and bending inner element-scale stiffness are elaborated. The multiscale metamaterial performance is numerically probed for variable-section, multiscale honeycomb, square and re-entrant star-shaped lattice architectures. It is observed that a substantial normal, bulk and shear specific stiffness increase can be achieved, which differs depending on the upper-scale lattice pattern. Subsequently, extended mechanical datasets are created for the training of machine learning models of the metamaterial performance. Thereupon, neural network (NN) architectures and modeling parameters that can robustly capture the multiscale material response are identified. It is demonstrated that rather low-numerical-cost NN models can assess the complete set of elastic properties with substantial accuracy, providing a direct link between the underlying design parameters and the macroscale metamaterial performance. Moreover, inverse, multi-objective engineering tasks become feasible. It is shown that unified machine-learning-based representation allows for the inverse identification of the inner multiscale structural topology and base material parameters that optimally meet multiple macroscale performance objectives, coupling the NN metamaterial models with genetic algorithm-based optimization schemes.


Introduction
Progress in additive manufacturing has allowed for the engineering of customized structures with highly refined inner designs, well below the micrometer scale [1]. Mechanical parts have been fabricated with characteristics that would have been infeasible using traditional manufacturing processes [2]. Moreover, a wide range of advanced materialsnamed metamaterials-have been developed, with tailorable effective mechanical properties [3][4][5] and design flexibility that has opened new frontiers in the control of the functional response of structural components [6][7][8].
Metamaterials have typically been based on the design of periodic inner structural patterns that can yield macroscopic mechanical properties with fundamentally different attributes from the ones observed for the base material used [9,10]. This quest has fostered the engineering of advanced materials with extraordinary mechanical behaviors such as the development of auxetics [11,12], and thus, materials which laterally expand instead of contracting upon the application of tensile loads [13][14][15][16]. The aforementioned nonconventional volumetric response has been primarily materialized by re-entrant patterns, with indicative examples being the star-shaped or the re-entrant honeycomb periodic cell designs [17][18][19]. The effective metamaterial normal E, bulk K, shear G and Poisson's ratio ν attributes are a function of the S2 and S1 scale design.
One of the major challenges in the development of high-performing machine learning methods of this kind is the generation of databases that are informative enough for the relevant solid mechanics response to be captured. Finite element models have been used to generate datasets for high-contrast [43,57] and tessellate composites [58]. Moreover, machine learning models have been developed to identify the relationship between salient structural features and the uniaxial compressive response of advanced materials. In particular, the geometric design of cellular materials has been used as an input for deep neural networks to be trained, so as to predict the macroscale uniaxial response of nonuniform cellular solids [59].
While it is currently well-established that the insertion of inner-material scales beyond the primal cellular pattern provides additional degrees of freedom for design, the corresponding material spaces remain, to a great extent, unexplored. Existing contributions have mainly emphasized on hollow inner designs [60,61], with the form of their inner structural components remaining unchanged. In particular, multiscale cellular metamaterials with hollow, variable-section inner designs have not been yet analyzed, their macroscale material performance remaining unquantified, both numerically and experimentally. Moreover, relevant data-based machine learning models that can associate multiple inner-scale input design parameters with the complete set of effective macroscale cellular material attributes have not been developed. Therefore, the modeling specifications and complexity required for such a task remain uncharacterized. Moreover, the potential use of surrogate neural network models of multiscale cellular material performance in the inverse identification of structural patterns has not been investigated.
In the current work, multiscale cellular material designs with hollow, variable-section inner designs are investigated. In particular, metamaterial designs with inner elements that are both hollow and follow a variable cross-sectional profile are analyzed for the first time. For the analysis, extensive numerical finite element simulations are performed, corroborated by selective experimental testing and analytical modeling results. It is shown that appropriate tuning of the second innermost scale can allow for the creation of lightweight, multiscale metamaterials, with enhanced normal, shear and bulk properties (Section 4). Furthermore, the necessary specifications for a neural-network-based association of the multi-scale input design parameters, with a complete set of constitutive and relative density effective properties, are identified. The high fidelity and low computation of the neural network model are highlighted, along with its potential coupling with inverse engineering analysis methods (Section 5).

Figure 1.
A multiscale metamaterial design, with the effective material properties (S 0 ) to be determined by the first-scale unit-cell (UC, S 1 ) and inner-element scale design (S 2 ) of the periodic pattern. The effective metamaterial normal E, bulk K, shear G and Poisson's ratio ν attributes are a function of the S 2 and S 1 scale design.
One of the major challenges in the development of high-performing machine learning methods of this kind is the generation of databases that are informative enough for the relevant solid mechanics response to be captured. Finite element models have been used to generate datasets for high-contrast [43,57] and tessellate composites [58]. Moreover, machine learning models have been developed to identify the relationship between salient structural features and the uniaxial compressive response of advanced materials. In particular, the geometric design of cellular materials has been used as an input for deep neural networks to be trained, so as to predict the macroscale uniaxial response of non-uniform cellular solids [59].
While it is currently well-established that the insertion of inner-material scales beyond the primal cellular pattern provides additional degrees of freedom for design, the corresponding material spaces remain, to a great extent, unexplored. Existing contributions have mainly emphasized on hollow inner designs [60,61], with the form of their inner structural components remaining unchanged. In particular, multiscale cellular metamaterials with hollow, variable-section inner designs have not been yet analyzed, their macroscale material performance remaining unquantified, both numerically and experimentally. Moreover, relevant data-based machine learning models that can associate multiple inner-scale input design parameters with the complete set of effective macroscale cellular material attributes have not been developed. Therefore, the modeling specifications and complexity required for such a task remain uncharacterized. Moreover, the potential use of surrogate neural network models of multiscale cellular material performance in the inverse identification of structural patterns has not been investigated.
In the current work, multiscale cellular material designs with hollow, variable-section inner designs are investigated. In particular, metamaterial designs with inner elements that are both hollow and follow a variable cross-sectional profile are analyzed for the first time. For the analysis, extensive numerical finite element simulations are performed, corroborated by selective experimental testing and analytical modeling results. It is shown that appropriate tuning of the second innermost scale can allow for the creation of lightweight, multiscale metamaterials, with enhanced normal, shear and bulk properties (Section 4). Furthermore, the necessary specifications for a neural-network-based association of the multi-scale input design parameters, with a complete set of constitutive and relative density effective properties, are identified. The high fidelity and low computation of the neural network model are highlighted, along with its potential coupling with inverse engineering analysis methods (Section 5).

Analytical and Numerical Characterization
We consider multiscale lattice material architectures composed of elements with a hollow variable-section inner form. The elements have outer and inner cross-sectional thickness t o and t i , respectively, at their ends, as schematically depicted in Figure 2a. Both the outer-and inner-element thickness vary along the element length L, following a sinusoidal, half-period wave-form evolution of magnitude e-here named swelling-at the middle of the element (Figure 2a). The element's area A(x) and moment of inertia I e (x), thus, vary along the element length, and their evolution is defined as follows: where t opl and t v opl in Equation (1) stand for the out-of-plane total and hollow thickness parts, respectively (Figure 2a), while sin[πx/L] characterizes the sinusoidal strut crosssection evolution. The hollow variable-section element form leads to a modification of the inner mass distribution, which depends on the swelling e and the inner-to-outer element thickness t i /t o selection. Using as a reference the prismatic element case with a volume V p = t o t opl L, the element-scale density modification ρ * e is defined as: where ρ * e in Equation (2) simplifies to unity for the prismatic, non-hollow case with zero swelling e, while the expression is independent of the element length L. The effective relative density ρ * at the metamaterial macroscale is determined by the relative density of the primal lattice design ρ * S 1 , accounting for the second-scale density modification ρ * e . In Figure 2, different multiscale, hollow variable-section square (b), honeycomb (c) and starshaped re-entrant (d) designs are depicted. The first-scale star-shaped lattice patterns refer to the auxetic lattice variant [17] and not to petal-form architecture variants, as investigated in [62].  For the computation of the effective multiscale metamaterial properties, two equivalent but distinct methodologies are followed. In the first case, multi-scale periodic finite element Abaqus models are created with a relative density of 0.03, 0.04 and 0.05 and an e value of 0.5, with element length values of 24, 30 and 40 mm and a / i o t t of 0.5. The inner metamaterial structure is discretized with tetrahedral C3D10 solids using a fine mesh of more than 20 elements per inner constituent. For each of the lattice multiscale patterns (Figure 2e), a total of two independent loading cases, namely of a uniaxial tensile and shear load, suffice to compute the effective normal, shear and bulk properties. We note that all multiscale patterns investigated fall within the tetragonal symmetry space, with equal x E and y E moduli, as well as equal Poisson's ratio values xy  and yx  . The For the computation of the effective multiscale metamaterial properties, two equivalent but distinct methodologies are followed. In the first case, multi-scale periodic finite element Abaqus models are created with a relative density of 0.03, 0.04 and 0.05 and an e value of 0.5, with element length values of 24, 30 and 40 mm and a t i /t o of 0.5. The inner metamaterial structure is discretized with tetrahedral C3D10 solids using a fine mesh of more than 20 elements per inner constituent. For each of the lattice multiscale patterns (Figure 2e), a total of two independent loading cases, namely of a uniaxial tensile and shear load, suffice to compute the effective normal, shear and bulk properties. We note that all multiscale patterns investigated fall within the tetragonal symmetry space, with equal E x and E y moduli, as well as equal Poisson's ratio values ν xy and ν yx . The mechanical properties are obtained through the application of infinitesimal normal and shear strains ε x and γ xy of 0.001.
For periodic models with more than five unit-cell repetitions in each material direction, more than 60,000 elements are required for a single design case for convergent analysis results. The multiscale nature of the geometry poses practical analysis limitations in the parametric investigation of a wide range of inner multiscale topologies, and thus, different inner geometries and lattice configurations. This, in turn, limits the dataset size that can be created, an important parameter for the neural network analysis part elaborated in Section 3.
In the second methodological modeling approach, the previous limitation is surpassed, applying a two-scale homogenization process. Initially, the homogenized axial and bending stiffness of the inner variable-section element geometry are computed as a function of geometric and material properties (t o , t i , L, e, E s ) for the innermost scale S 2 (Equation (1), Figure 1). Thereafter, the macroscale lattice's effective attributes are obtained based on the homogenized second-scale S 2 properties, which are employed as inputs for the first-scale S 1 computations, making use of the homogenization algorithm elaborated upon in [32]. The effective metamaterial properties (E, K, G, ν) are a direct result of the asymptotic homogenization, which computes the complete flexibility and stiffness matrix tensor of a given multiscale periodic pattern, provided with the homogenized normal and bending stiffness attributes of its inner constituents. The reader is referred to [32] for a detailed description of the theoretical formulation, as well as for its numerical implementation. The previously explicated algorithmic process allows for the parametric computation of a wide range of multiscale cellular configurations.
For the computation of the homogenized S 2 element-scale stiffness attributes, analytical formulas for the effective axial k e N and bending k e B stiffness attributes are derived, using the unit-force method. The corresponding analytical mechanics expressions contain integral forms with varying area and moment-of-inertia properties (Equation (1)). In particular, the axial flexibility f e N is computed as M · M /(EI e (x))dx − M being the virtual bending moment developed within a clamped beam element, subsect to a shear type forces at its ends -. The corresponding stiffness terms are obtained by computing the inverse of the flexibility components, k e N = 1/ f e N and k e B = 1/ f e B accordingly. The reader is referred to [27,63] for a detailed description of the virtual force method.

Additive Manufacturing and Experimental Characterization
The elastic stiffness of two-scale honeycomb metamaterial designs was experimentally investigated using 3D-printed periodic specimens (Figure 3a). Two distinct, multiscale honeycomb designs were fabricated with an element length L of 7.1 mm, an out-of-plane thickness of 10 mm (Figure 3b,c) and different inner-element hollow profiles. In particular, element profiles with an e value of 0.17 and 0.36 mm were fabricated, with the corresponding honeycomb periodic designs named D1 and D2. The L/t o ratio was set, in all design cases, as equal to 18, while the t opl /t o parameter was set to 0.5. For the fabrication, a BMF 3D printer was employed with a high-temperature resin. For each honeycomb design, three repeat samples were conducted. A total fabrication time of several tenths of hours was required for the 3D printing of the specimens, noting the substantially high number of slices necessary for a sufficiently accurate geometric representation. The multiscale honeycomb structures consisted of a total of 20 periodic cells stacked in the out-of-plane thickness direction ( Figure 3c); the stacking was employed to ensure sufficient stiffness and surface area along the specimen thickness direction for the load application upon mechanical testing. For the testing, an Instron machine with a 50 N loadcell was employed. A rather low testing speed of 1 mm/min was used, corresponding to a strain rate below 0.01/s to ensure static loading conditions for all specimens.

Machine-Learning-Based Modeling and Design of Multiscale Metamaterial Architectures
Effective macroscale metamaterial properties were determined by a substantial number of inner design parameters. Even for a given base material ( s E ) and first-scale 1 S unit-cell pattern (e.g., honeycomb, square, Figure 2), the metamaterial attributes were a function of the element thickness o t and length L , as well as the inner-element hollow thickness i t and swelling parameter e . For a single combination of the previously reported parameters, a set of a minimum of two independent computational tasks needed to be performed for the effective normal moduli E , bulk K , shear G and Poisson's value  to be obtained. More importantly, no explicit functions associating the effective macroscale properties with the inner first 1 S and second 2 S scale design parameters ( Figure 1) were available, so that multiscale inverse engineering tasks could not be performed.
For the modeling of the mechanical response of a wide range of inner-material designs, a machine-learning-based approach, without prior assumptions with respect to the observed macroscale constitutive metamaterial response, was followed. Five primal inner design parameters were used as input features (I), namely two first-scale (  For the testing, an Instron machine with a 50 N loadcell was employed. A rather low testing speed of 1 mm/min was used, corresponding to a strain rate below 0.01/s to ensure static loading conditions for all specimens.

Machine-Learning-Based Modeling and Design of Multiscale Metamaterial Architectures
Effective macroscale metamaterial properties were determined by a substantial number of inner design parameters. Even for a given base material (E s ) and first-scale S 1 unit-cell pattern (e.g., honeycomb, square, Figure 2), the metamaterial attributes were a function of the element thickness t o and length L, as well as the inner-element hollow thickness t i and swelling parameter e. For a single combination of the previously reported parameters, a set of a minimum of two independent computational tasks needed to be performed for the effective normal moduli E, bulk K, shear G and Poisson's value ν to be obtained. More importantly, no explicit functions associating the effective macroscale properties with the inner first S 1 and second S 2 scale design parameters ( Figure 1) were available, so that multiscale inverse engineering tasks could not be performed.
For the modeling of the mechanical response of a wide range of inner-material designs, a machine-learning-based approach, without prior assumptions with respect to the observed macroscale constitutive metamaterial response, was followed. Five primal inner design parameters were used as input features (I), namely two first-scale (S 1 : t o , L) and two second-scale (S 2 : t i , e) attributes, as well as the base material modulus E s . For the current analysis, a unit out-of-plane thickness t opl , along with a half unit t v opl out-of-plane void thickness part, were considered (Figure 2a), without loss of generality. The neural network output parameters (O) included the complete set of the elastic E, bulk K and shear G values, as well as the Poisson's ratio value ν ( Figure 1) and relative density ρ * data, computed using multiscale homogenization (Section 2.1). The neural network model can be viewed as a multivariate regressor of the mechanical performance. We note that all multiscale patterns investigated fell within the tetragonal symmetry space, so that equal E x and E y moduli and Poisson's ratio values ν xy and ν yx applied, allowing for the simplification of the corresponding notations.
For each unit-cell primal design case, several effective metamaterial performance data in each feature direction were created. In particular, eleven data points along the element thickness t o , seven data points along the inner hollow element thickness t i and element length L, and eleven data points along the swelling e feature direction were employed. The data points were created in a hierarchical manner, so that for a given element thickness t o and base material moduli E s at the first cellular scale S 1 , normalized slenderness and second-scale feature attributes were created (t i /t o , t o /L, e/t o ), covering all possible parameter combinations (summarized in Table 1). For each feature, uniform spacing among the indicated bounds was applied for simplicity. It is noted that the bounds prescribed covered a wide range of inner cellular designs, albeit non-exhaustive with respect to the possible design space of multiscale metamaterial patterns. The number of sampling points for each parameter was selected based on initial reduced input-dimensionalityfitting studies; they were defined so as to allow for high-accuracy results, retaining an overall low data-creation computational cost. Moreover, it should be underlined that the multiscale metamaterial performance is affine with respect to the base material modulus E s , so that it could be omitted from the modelling process, if all the results were to be obtained in a normalized, non-dimensional form. The effective metamaterial attributes for each multiscale design were computed using the two-scale homogenization process explicated at the end of Section 2.1. Table 1. First-scale S 1 and second-scale S 2 metamaterial design parameter range used for the creation of training datasets for multiscale square, honeycomb and star re-entrant lattice patterns.

Input Features Range Sampling Points
The parameter range of Table 1, along with the discretization introduced for each feature, leads to more than 53,000 metamaterial design cases for a given multiscale lattice pattern (Figure 1). The different unit-cell multiscale designs were created in a parametric manner, following the above elaborated process. The input features (I) of Table 1 are related to the macroscale effective metamaterial attributes (O) through a neural network model ( Figure 4). The network design and computational complexity, as characterized by the number of layers n and neurons m per layer (Figure 4a) are parameters to be determined throughout the training process [64]. Neural network architectures with a minimum of two and up to five hidden layers were considered, with a minimum of 5 and up to 20 neurons per layer. A schematic of the machine modeling architecture is provided in Figure 4a.  For the training of the model, different activation functions were investigated, including the tanh, sigmoid, and deep learning common-purpose activation functions, such as the relu activation function [43]. For the training process, the mean squared error function was used (MSE), along with the Levenberg-Marquardt default optimization algorithm. The data were shuffled before training. The mean squared error was defined through the difference of the neural network (NN) from the multiscale homogenization For the training of the model, different activation functions were investigated, including the tanh, sigmoid, and deep learning common-purpose activation functions, such as the relu activation function [43]. For the training process, the mean squared error function was used (MSE), along with the Levenberg-Marquardt default optimization algorithm. The data were shuffled before training. The mean squared error was defined through the difference of the neural network (NN) from the multiscale homogenization (MH) mechanical . Throughout the training process, 30% of the available data created with the process previously explicated were used for testing. The trained NN models were independently controlled with respect to their accuracy, using an additional validation dataset of 1000 multiscale cellular designs, created through random input feature generation in a design space that exceeded the bounds set for each feature in Table 1 by up to 20%. The weighting matrix coefficients W i , the layer constants W i 0 , and the layer activation function f and network depth n were parameters to be determined throughout the training process [64].
The trained model was used as a surrogate for its coupling with the genetic, multiobjective optimization algorithm (Non-dominated Sorting Genetic Algorithm, Figure 4b) elaborated in [65] for the identification of Pareto set solutions that satisfy both base material moduli and density targets at the metamaterial macroscale. For the inverse analysis, the neural network modeling parameters of Figure 4a were used as inputs, while the effective elastic modulus and relative density value were concurrently used as macroscale optimization objectives. For the generic algorithm computation, a probability of crossover and mutation of 0.9 and 0.5 were used accordingly, along with a mutation parameter of 0.05 [65]. The population size and number of generations required for convergence are discussed in Section 5.

Effective Mechanical Attributes of Multiscale, Variable Inner Section Cellular Materials
The variable-section, doubly sinusoidal variation of the inner geometry result in a non-linear area and moment of inertia distribution along the element length, as indicated by the function definitions of Equation (1). An analytical computation of the effective axial k e N and bending k e B stiffness at the element scale using the classical mechanics formulations summarized in Section 2.1 is infeasible with the use of commercial integration routines (Mathematica 11.3, Maple 2020), retaining the sinusoidal geometric definitions of Equation (1). However, by employing the Bhaskara approximation of the sin(x) function [66] (Appendix A), closed-form expressions of the element's effective axial and bending stiffness attributes are obtained, as follows: where αcoth in Equation (5)  section design on the specific bending-element stiffness is illustrated in Figure 5b as a function of the swelling-to-length ratio e/L. The modified inner (S 2 , Figure 1) normal and bending stiffness attributes can lead to a different density scaling of the effective elastic macroscale cellular properties. In Figure 5c, the scaling law of the effective modulus with respect to the relative density (E/E s = Aρ * n ) of variable-section, multiscale honeycombs (MH), with a swelling parameter value e = 0.5, is provided, along with the case of prismatic, single-scale honeycombs (PH). The two-scale MH mechanical properties are computed using the multiscale homogenization process from Section 2.1. The periodic finite-element modeling results (P-FEM, Figure 2e) are depicted by blue cross symbols. The relative density of the hollow variable-section geometry is below unity over a wide range of variable-section hollow element forms (Figure 5a), allowing for lightweight material designs. More importantly, the specific bending stiffness of the variable-section inner geometry is enhanced (Figure 5b), with rather small e values to yield substantial stiffness improvements compared to the prismatic design case. It is noted that the inner stiffness enhancement is strongly non-linear. Doubling the maximum swelling e increases the specific bending stiffness more than two times (Figure 5b).
The improved second-scale material performance can modify the scaling of the effective macroscale elastic modulus E with respect to the metamaterial relative density *  . In particular, the least square fitting for the PH case yields A and n coefficients of 1.3 and 3, respectively, which is in agreement with the theoretical predictions [20] (e = 0, data from Table 1). For the MH case with e = 0.5, the A coefficient more than doubles to 2.8 compared to the theoretical prediction of 1.3 for the PH case, while the corresponding density exponent n is reduced to 2.75.
In Figure 6a, the experimentally obtained stress-strain results for the D1 and D2 design cases of Section 2.2 are provided, along with the prismatic reference case. The exper- Figure 5. Dependence of the relative density of hollow variable-section elements on the swelling-toend-element-thickness ratio e/t o and on the t i /t o ratio (a). The enhancement of the inner specific bending stiffness (normalized to ρ * e of Equation (2) and to the prismatic bending stiffness element case) as a function of the swelling-to-length e/L and inner-to-outer-element-thickness t i /t o is provided in (b). A length L = 20 and a t o , t opl and t v opl of 1, 1 and 0.5 mm are used, respectively. The effective elastic modulus scaling with relative density for a multiscale honeycomb lattice (e = 0.5) and a prismatic honeycomb (e = 0.0) is depicted in (c).
The relative density of the hollow variable-section geometry is below unity over a wide range of variable-section hollow element forms (Figure 5a), allowing for lightweight material designs. More importantly, the specific bending stiffness of the variable-section inner geometry is enhanced (Figure 5b), with rather small e values to yield substantial stiffness improvements compared to the prismatic design case. It is noted that the inner stiffness enhancement is strongly non-linear. Doubling the maximum swelling e increases the specific bending stiffness more than two times (Figure 5b).
The improved second-scale material performance can modify the scaling of the effective macroscale elastic modulus E with respect to the metamaterial relative density ρ * . In particular, the least square fitting for the PH case yields A and n coefficients of 1.3 and 3, respectively, which is in agreement with the theoretical predictions [20] (e = 0, data from Table 1). For the MH case with e = 0.5, the A coefficient more than doubles to 2.8 compared to the theoretical prediction of 1.3 for the PH case, while the corresponding density exponent n is reduced to 2.75.
In Figure 6a, the experimentally obtained stress-strain results for the D1 and D2 design cases of Section 2.2 are provided, along with the prismatic reference case. The experimentally obtained (E) elastic stiffness for each design case is depicted in Figure 6b-computed in the low-strain range (<1%)-along with the multiscale homogenization numerical analysis results (N). For each design case, two distinct experimental repeats are presented (R1, R2 in Figure 6). The elastic stiffness of the moderate e design case D1 is more than three times higher the prismatic, single-scale design case, a performance-improvement scaling factor that exceeds 7 for the high e D2 design case (Figure 6b). The experimental data lie in good accordance with the numerical predictions for both density cases tested, verifying the substantial difference of the MH case with respect to the PH single-scale designs. It should be noted that the stiffness increase effect observed upon the use of variable-section hollow inner-strut designs does not apply to all first-scale cellular patterns, as explicated in the sequel.
In Figure 7, the effect of the inner second-scale on the macroscale elastic E , shear G and bulk K material properties of multiscale square (s)-, honeycomb (h)-and star re-entrant (sr)-shaped cellular patterns is analyzed in a comparative manner. In particular, not only the specific elastic normal modulus, but also the specific shear and bulk modulus  The experimental results of Figure 6a indicate a clear elastic stiffness increase with higher swelling parameter e values (Figure 6a), with the highest swelling D2 design case considerably outperforming the stiffness of the reference prismatic design. The stiffer response is retained throughout the entire strain range investigated. It should be noted that the enhanced constitutive performance can be obtained at a practically invariant relative density value of the effective metamaterial, noting the near unity ρ * e value corresponding to the design specifications of the multiscale honeycomb designs D1 (Figure 6a). The elastic stiffness of the moderate e design case D1 is more than three times higher the prismatic, single-scale design case, a performance-improvement scaling factor that exceeds 7 for the high e D2 design case (Figure 6b). The experimental data lie in good accordance with the numerical predictions for both density cases tested, verifying the substantial difference of the MH case with respect to the PH single-scale designs. It should be noted that the stiffness increase effect observed upon the use of variable-section hollow inner-strut designs does not apply to all first-scale cellular patterns, as explicated in the sequel.
In Figure 7, the effect of the inner second-scale on the macroscale elastic E, shear G and bulk K material properties of multiscale square (s)-, honeycomb (h)-and star re-entrant (sr)-shaped cellular patterns is analyzed in a comparative manner. In particular, not only the specific elastic normal modulus, but also the specific shear and bulk modulus (Figure 7a-c) are given for different swelling-to-length e/L values. The mechanical properties are normalized with respect to the prismatic design reference case. All computations pertain to an L/t o of 20 with a t i /t o of 0.8 and a t o of 1 mm. Moreover, a unit out-of-plane thickness t opl , along with a half-unit t v opl out-of-plane void thickness part, is considered. The computations are performed following the multiscale homogenization process elaborated upon in Section 2.1. In Figure 7d, the evolutions of the normal-to-shear E/G and bulk-to-shear K/G moduli ratios over the same range of e/L values are provided.  (Figure 7d), contrary to the square lattice case, where a shear-stiffer response is observed (Figure 7d). Accordingly, the bulk-to-shear effective metamaterial performance ( / K G ) becomes shear-stiffer for the square and honeycomb lattice case, remaining unaffected for the re-entrant square lattice pattern. The employment of hollow, variable-section inner designs at the second inner scale (S 2 ) of the cellular material leads to an increased specific metamaterial shear stiffness for all lattice patterns in Figure 2. More specifically, the specific shear stiffness enhancement is up to approximately seven times the shear resistance of the reference hollow lattice pattern for an e/L value of 5%, with the increase being non-linear with respect to the swelling-to-length parameter. An analogous improvement is recorded for the axial modulus (E/E e=0 )/(ρ * /ρ * e=0 ) of the honeycomb and square re-entrant lattice patterns, with the specific axial modulus of the square lattice remaining practically unaffected (Figure 4b). Accordingly, the specific bulk modulus (K/K e=0 )/(ρ * /ρ * e=0 ) of the star re-entrant metamaterial design increases, contrary to the bulk resistance of the square and honeycomb lattice patterns, which remain practically invariant.
The multiscale designs yield a macroscale cellular material response (Figure 7a-c) that differs for each of the first-scale lattice patterns (S 1 ), so that the material design space extension induced by the innermost architectural scale (S 2 ) is non-unique. More specifically, the relative performance of the multiscale, variable-section metamaterial designs to normal loads and shear loads (E/G) remains practically unaffected for the honeycomb and reentrant square lattice patterns (Figure 7d), contrary to the square lattice case, where a shear-stiffer response is observed (Figure 7d). Accordingly, the bulk-to-shear effective metamaterial performance (K/G) becomes shear-stiffer for the square and honeycomb lattice case, remaining unaffected for the re-entrant square lattice pattern.

Neural-Network-Based Multiscale Metamaterial Forward Modeling and Inverse Design
Multiscale metamaterial performance is determined by several first-and second-scale parameters, as well as by the base material modulus E s . To predict the complete set of effective metamaterial properties, as summarized in Figure 4a, different depth neural network architectures were probed. It was observed that a highly accurate mapping of the input design features (I) with the output effective material properties (O) is feasible for all activation functions listed in Section 3. However, networks using the ReLu activation function require more than double the number of training parameters and deeper network architecture than the corresponding sigmoid or hyperbolic tangent counterparts for the same accuracy level to be achieved. Networks with a total number of four layers and less than one thousand training parameters were observed to provide the highest modeling accuracies, as quantified by the training process. The minimum testing loss performance was obtained for networks with 20 × 15 × 5 × 5, 16 × 10 × 5 × 4 and 16 × 12 × 6 × 5 neurons for the multiscale honeycomb, square and re-entrant square cellular patterns. All network architectures are provided in the form of supplementary material. Networks with higher computing cost and comparable loss performance were excluded from the analysis. In the independent validation sets used for each of the multiscale cellular patterns here investigated, a relative error among the multiscale homogenization and the NN predictions below 1% for all QoI of Figure 4b was recorded. In Figure 8a, insights in the training and testing loss curves for the case of a multiscale honeycomb lattice over a training period of 5000 epochs are provided. Figure 8b-d depict frequency error bars for the bulk (K) and normal (E) modulus and relative density neural network predictions for the validation dataset of the 1000 randomly generated multiscale cellular designs (Section 3). For the relative error estimation, the multiscale homogenization (MH) results are used as a reference for the NN prediction comparisons.

Neural-Network-Based Multiscale Metamaterial Forward Modeling and Inverse Design
Multiscale metamaterial performance is determined by several first-and secondscale parameters, as well as by the base material modulus s E . To predict the complete set of effective metamaterial properties, as summarized in Figure 4a, different depth neural network architectures were probed. It was observed that a highly accurate mapping of the input design features (I) with the output effective material properties (O) is feasible for all activation functions listed in Section 3. However, networks using the ReLu activation function require more than double the number of training parameters and deeper network architecture than the corresponding sigmoid or hyperbolic tangent counterparts for the same accuracy level to be achieved. Networks with a total number of four layers and less than one thousand training parameters were observed to provide the highest modeling accuracies, as quantified by the training process. The minimum testing loss performance was obtained for networks with 20 × 15 × 5 × 5, 16 × 10 × 5 × 4 and 16 × 12 × 6 × 5 neurons for the multiscale honeycomb, square and re-entrant square cellular patterns. All network architectures are provided in the form of supplementary material. Networks with higher computing cost and comparable loss performance were excluded from the analysis. In the independent validation sets used for each of the multiscale cellular patterns here investigated, a relative error among the multiscale homogenization and the NN predictions below 1% for all QoI of Figure 4b was recorded. In Figure 8a, insights in the training and testing loss curves for the case of a multiscale honeycomb lattice over a training period of 5000 epochs are provided. Figure 8b-d depict frequency error bars for the bulk ( K ) and normal ( E ) modulus and relative density neural network predictions for the validation dataset of the 1000 randomly generated multiscale cellular designs (Section 3). For the relative error estimation, the multiscale homogenization (MH) results are used as a reference for the NN prediction comparisons.  The loss curves of Figure 8a reveal a multi-level training process, which spans several orders of magnitude of performance improvement over 5000 epochs. The trained neural network model can predict the effective multiscale properties with substantial accuracy comparable to state-of-the-art relevant modeling contributions [40] over the entire space, as indicated from the loss curve in Figure 8a. The error frequency distribution for the bulk, elastic and relative density predictions denotes a relative material performance prediction error that is substantially lower than 1% for all NN model predictions (Figure 8b-d).
It should be noted that the sharp variation in the loss function in Figure 8a has to be associated-to a certain extent-with the incorporation of the base material modulus E s in a dimensional form in the input neural network model parameters (Section 3). All trained neural network architectures are provided in the form of complementary material for completeness purposes.
The elaborated NN-based metamaterial models provide a direct link of the inner scales (Figure 4a) with the complete set of elastic effective metamaterial attributes. As a result, inverse multiparametric identification tasks become feasible. In particular, given a set of macroscale material performance objectives, the inner design specifications that best match the prescribed criteria can be probed. We note that the inverse structural identification analysis conducted herein is restrained within the mechanical performance limits of the multiscale cellular patterns investigated. As such, its scope and functionality need to be clearly separated from free-morphology or full-topology optimization methods [49,51], which are beyond the analysis range and context of the current contribution. Using the genetic algorithm elaborated upon in [65], the base material modulus E s , along with the second S 2 and first S 1 material scale features required for a multiscale lattice pattern, to yield a desirable macroscale performance, can be identified. In Figure 9a,b, the evolution of the Pareto front for a target elastic modulus O 1 and relative density O 2 objective are provided. A target modulus E t of 2 GPa and 0.1 GPa, with relative densities ρ * of 0.1 (O 2 , Figure 9a,b), respectively, are employed for a multiscale square and square reentrant design. Accordingly, the evolution of the Pareto front solution for the case of a multiscale honeycomb lattice pattern is provided in Figure 9c. For simplicity, the hollow element thickness is set to be fifty percent of the outer-end element thickness t o in all computations. One of the best sorting structural patterns is provided in each case, along with the evolution of the optimal Pareto front solution (using a population size of 30 and 150 generations). In Figure 9d, the identified optimal effective material modulus is given in the form of an Ashby diagram, while the base material modulus required for each target metamaterial performance is provided in Figure 9a-c.
The results of Figure 9 provide insights into the inverse engineering potential of the neural network multiscale metamaterial models developed. More specifically, given the high accuracy and the low computational cost of a single model evaluation, multiobjective, genetic-algorithm-based computations, requiring thousands of model evaluations are feasible in a few seconds. By that means, the base material modulus E s and the geometric features of the second (S 2 ) and first (S 1 ) design scale can be probed to optimally meet the macroscale material objectives (Figure 9). Moreover, the potential of a given multiscale cellular material to meet a set of macroscale objectives (O) can be robustly assessed, providing quantitative estimates of its optimality in each material performance direction. This becomes explicitly evident in Figure 9b, where none of the objectives can be fully satisfied. However, the Pareto front solutions allow for the inverse identification of inner multiscale cellular patterns that either fully satisfy one of the objectives, or partially satisfy each, quantifying the relative performance discrepancies. What is more, the inferred base-material moduli values E s can be directly used to assess the manufacturing feasibility of the metamaterial part, confining the base-material selection range and 3D printing technologies that could be employed in the fabrication process. Figure 9. Optimal Pareto front for a multiscale square lattice with a target modulus of E t = 2 GPa and a relative density of 0.1 (a). The optimal pareto front solution for a multiscale square re-entrant and multiscale honeycomb cellular with a target modulus of E t = 0.1 GPa and E t = 0.05 GPa and a relative density of 0.1 are provided in (b) and (c), respectively. Inversely identified structural multiscale cellular patterns are depicted in each case. The metamaterial (MM) effective modulus is presented in an Ashby diagram comparative form in (d). For completeness, the modulus-density ranges for foams (F), elastomers (E), polymers (P), metals (M), composites (Co) and ceramics (Ce) are included.
The results of Figure 9 provide insights into the inverse engineering potential of the neural network multiscale metamaterial models developed. More specifically, given the high accuracy and the low computational cost of a single model evaluation, multi-objective, genetic-algorithm-based computations, requiring thousands of model evaluations are feasible in a few seconds. By that means, the base material modulus s E and the geometric features of the second ( 2 S ) and first ( 1 S ) design scale can be probed to optimally meet the macroscale material objectives (Figure 9). Moreover, the potential of a given multiscale cellular material to meet a set of macroscale objectives (O) can be robustly assessed, providing quantitative estimates of its optimality in each material performance direction. This becomes explicitly evident in Figure 9b, where none of the objectives can be fully satisfied. However, the Pareto front solutions allow for the inverse identification of inner multiscale cellular patterns that either fully satisfy one of the objectives, or partially satisfy each, quantifying the relative performance discrepancies. What is more, the inferred basematerial moduli values s E can be directly used to assess the manufacturing feasibility of the metamaterial part, confining the base-material selection range and 3D printing technologies that could be employed in the fabrication process.

Discussion
While it is well established that multiscale metamaterial architectures can provide additional degrees of freedom for design that can extend the performance limits of singlescale patterns, the corresponding bounds remain, to a great extent, unquantified. The employment of hollow variable-section elements provides the possibility to modulate the normal-and bending-stiffness attributes of the inner lattice constituents. The stiffness alteration is controlled by the variable-section inner geometry. Element profiles with a rather small variation in their inner-thickness-to-length ratio (e/L) allow for considerable bending stiffness enhancements (Section 2.1). The elaborated analytical stiffness forms indicate that the bending performance can be tuned depending on the inner length L and on the element profile geometric swelling parameter e. However, it can alter the reference bending stiffness performance by up to one order of magnitude for rather low e/L values, in a controlled manner (Figure 5b).
The multiscale designs modify the macroscale metamaterial performance in a nonlinear manner, which depends on both the first-scale S 1 pattern and on the second-scale design. In particular, while a specific shear stiffness increase is recorded for all the multiscale lattice patterns in Figure 7a (rs, s, h), denoting a shear strengthening per unit weight for all design cases, different behaviors with respect to the normal and bulk metamaterial attributes are observed. More specifically, variable-section multiscale honeycombs allow for increased normal moduli with an invariable specific bulk modulus (Figure 7), while multiscale square patterns retain invariance both in their specific bulk and normal stiffness characteristics. This observation indicates that the resistance of certain multiscale metamaterial designs to shape changes can be controlled through appropriate modifications of the second, innermost scale, with volume changes favored over shape changes. Moreover, there exist first-scale (S 1 ) lattice designs for which a simultaneous increase in the specific normal, shear and bulk stiffness is induced; this is the case for multiscale square-re-entrant patterns (Figure 7), with the relative volumetric and shape metamaterial resistance remaining practically invariable (Figure 7d).
The decoupled second-and first-scale geometry allows for the preservation of the connectivity of the upper-scale S 1 lattice pattern. All the multiscale lattice cases investigated above have a nodal connectivity that is less than 6, so that the bending stiffness of their inner constituents is significant [67,68], for the metamaterial's macroscale performance. Multiscale designs of this kind can yield enhanced density-scaling laws, which substantially differ from the single-scale metamaterial architecture, as demonstrated for the multiscale honeycomb case (Figure 5c).
However, the estimation of a complete set of effective mechanical properties can be a substantially cumbersome process due to the different scales and the number of parameters involved. The results of Section 5 indicate that a direct link between the inner multiscale geometric and material design attributes and the macroscale metamaterial performance is feasible, with the use of appropriately architected neural network models. In particular, rather low-computational-cost neural network architectures have been identified that can robustly predict the complete set of effective elastic properties, including the bulk and shear metamaterial response. The accuracy of the networks is in the third decimal order for all quantities of interest ( Figure 8), with the model execution time amounting to some fractions of a second on an ordinary personal computing system.
The low computing cost of the neural network models elaborated allows for their coupling with inverse, genetic-algorithm-based, multi-objective optimization schemes ( Figure 9) in a single-step process. By that means, a wide range of optimization engineering tasks can be performed. In particular, metamaterial designs that can optimally satisfy a combination of macroscale performance objectives can be inversely identified in the form of Pareto front solutions. This allows for the rigorous assessment of the engineering feasibility of a certain design; moreover, it allows for the identification of solution sets that can concurrently best meet performance requests within a given multiscale metamaterial class, as demonstrated in Figure 9b. In the inverse identification process, both the geometric and the base material modulus requirements are identified with the necessary moduli specifications and can be directly probed.

Conclusions
Overall, the current work investigated multiscale variable-element-section cellular metamaterial designs, combining analytical, numerical and experimental testing methods. Analytical expressions for the effect of variable-section, hollow element designs on the inner-axial and bending strut-scale stiffness were elaborated. The mechanical performance of a wide range of multiscale cellular designs were evaluated. It was observed that:

1.
Hollow, variable-section inner structural designs allow for an enhanced, specific normal, shear and bulk metamaterial response, well beyond the range of single-scale metamaterial architectures.

2.
The insertion of a second inner scale affects the macroscale metamaterial performance in a non-unique manner, which depends on the uppermost-scale cellular pattern design.

3.
Multiscale designs can modify the stiffness-to-density scaling of cellular materials from a bending-dominated towards a stretching-dominated response.

4.
Low-numerical-cost neural network models can derive a robust link between the different inner scales and the complete set of effective elastic cellular material properties.

5.
Inverse multi-objective engineering tasks can, therefore, be performed, identifying the optimal multiscale cellular patterns that best satisfy the macroscale performance requests.
We hope that the methodology elaborated and the results provided act as a framework in the analysis and design of multiscale cellular materials beyond the design cases investigated herein.