Inverse design method for tailoring the deflection of beams with functionally graded metamaterials

The convergence of Additive Manufacturing (AM) technologies and advanced simulation methodologies has propelled lattice structures and cellular metamaterials into the forefront of engineering materials [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib1" id="fnref-smsad6d1fbib1"="">1</a>, <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib2" id="fnref-smsad6d1fbib2"="">2</a>]. Mechanical metamaterials exhibit exceptional tunability and distinctive properties, catalyzing the development of tailor-made materials for specific mechanical requirements and applications [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib3" id="fnref-smsad6d1fbib3"="">3</a>].

Extensive research has focused on developing various topologies with distinct properties such as negative Poisson's ratio, energy absorption, and lightweight characteristics, among others [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib4" id="fnref-smsad6d1fbib4"="">4</a>]. This approach of crafting porous materials draws strong inspiration from nature, where the properties of structures like trees or bones stem not only from the materials they are made of but also from their micro-structures [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib5" id="fnref-smsad6d1fbib5"="">5</a>].

However, nature does not exhibit uniform distributions of these cellular materials; rather, it features non-periodic arrangements with segments possessing different compositions or properties [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib6" id="fnref-smsad6d1fbib6"="">6</a>, <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib7" id="fnref-smsad6d1fbib7"="">7</a>]. Similarly, ongoing efforts in metamaterials aim to combine different topologies or adjust topology parameters to achieve customized mechanical responses. Ramirez–Chavez <em="">et al</em> [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib8" id="fnref-smsad6d1fbib8"="">8</a>] recently classified aperiodic metamaterials into three distinct types: perturbations, hybridization, and gradation.

Perturbation involves introducing alterations to the periodicity at the nucleus or node level, thereby impacting the structure of the cellular material. Hybridization encompasses the fusion of two or more topologies within a single structure, achieved at either the unit cell or member level [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib9" id="fnref-smsad6d1fbib9"="">9</a>]. Lastly, gradation is a method involving the adjustment of parameters within the features of a cellular material [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib10" id="fnref-smsad6d1fbib10"="">10</a>]. Gradation can be implemented at the unit cell, within members like walls or beams, and at the nodes. When gradation is tailored to specific mechanical requirements or applications, the resultant structures are termed Functionally Graded Metamaterials (FGMM). They imply a notable advancement in achieving extensive material tunability. The approach of generating a design based on a targeted behavior is referred to as inverse design. In contrast to the direct design paradigm, where a material or structure is designed first, and then its properties are measured, the inverse design allows the generation of a design based on expected properties [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib11" id="fnref-smsad6d1fbib11"="">11</a>].

Various strategies have been pursued to create transitions between layers in a graded material. For instance, gradation can be classified as continuous, where the gradient seamlessly flows throughout the volume without distinct separations between layers, or discontinuous, where the gradient shifts incrementally, resulting in clear separations between layers [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib12" id="fnref-smsad6d1fbib12"="">12</a>]. Within the realm of metamaterials, smooth transitions between densities are highly desirable and recommended. Gradual transitions help reduce stress concentrations and prevent the risk of interfacial defects, such as cracking and delamination, that can occur from abrupt shifts in density. Kriging models have been employed to facilitate seamless transitions between different stages of two- and three-dimensional graded lattices [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib13" id="fnref-smsad6d1fbib13"="">13</a>]. At the same time, genetic algorithms have been utilized for transitions between unit cells or topologies [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib14" id="fnref-smsad6d1fbib14"="">14</a>]. Topology optimization-based methods have been proposed as well for transitioning between stages [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib15" id="fnref-smsad6d1fbib15"="">15</a>].

Significant efforts have been dedicated to the mechanical characterization of graded metamaterials by adjusting parameters like wall thickness or the orientation of strut-based topology elements [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib16" id="fnref-smsad6d1fbib16"="">16</a>, <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib17" id="fnref-smsad6d1fbib17"="">17</a>]. Studies focusing on graded hexagonal honeycombs have investigated the impact of variations in the size and composition of the unit cells on their energy absorption properties, yielding improved mechanical responses from the graded variations [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib18" id="fnref-smsad6d1fbib18"="">18</a>]. Similarly, the grading of wall thickness and the compressive behavior using varying grading stages in hexagonal honeycombs has also been explored, revealing a shift in energy-absorbing capabilities in graded structures [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib19" id="fnref-smsad6d1fbib19"="">19</a>]. Furthermore, hexagonal and reentrant honeycomb variations have been studied in graded structures under compression, resulting in an increase in their energy absorption capabilities through the progressive grading of wall thickness [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib20" id="fnref-smsad6d1fbib20"="">20</a>].

Extensive studies have investigated the grading of two-dimensional lattices based on rectangular perforations, varying parameters, and grading stages to examine their Poisson's ratio and stiffness, observing significant effect of the geometric features of the structures in their mechanical response [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib21" id="fnref-smsad6d1fbib21"="">21</a>, <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib22" id="fnref-smsad6d1fbib22"="">22</a>]. Research on graded metamaterials under compressive loads has extensively analyzed failure mechanisms associated with different density distributions within different topologies, [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib23" id="fnref-smsad6d1fbib23"="">23</a>], or various geometries distributed within a single structure under compressive [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib24" id="fnref-smsad6d1fbib24"="">24</a>, <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib25" id="fnref-smsad6d1fbib25"="">25</a>], and flexural loads [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib26" id="fnref-smsad6d1fbib26"="">26</a>]. These studies report that the graded structures presented superior mechanical stability, strain concentration reduction, and enhanced energy absorption capacities. However, the previously highlighted works focused on specific gradation modes without being driven by any particular functionality or requirement.

Numerous methods have been proposed for developing structures with FGMM: topology optimization-based, implicit-based, and deep learning-assisted methods. Topology optimization stands as the most commonly employed method for designing structures with FGMM [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib27" id="fnref-smsad6d1fbib27"="">27</a>], employing different strategies such as homogenization method [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib28" id="fnref-smsad6d1fbib28"="">28</a>], solid isotropic material with penalization method [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib29" id="fnref-smsad6d1fbib29"="">29</a>], evolutionary structural optimization method [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib30" id="fnref-smsad6d1fbib30"="">30</a>], and level-set methods [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib31" id="fnref-smsad6d1fbib31"="">31</a>].

For instance, studies have analyzed the use of cell-wise graded distributions in sandwich structures to optimize their response to flexural loads [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib32" id="fnref-smsad6d1fbib32"="">32</a>], and different loading conditions for segment-wise graded distributions [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib33" id="fnref-smsad6d1fbib33"="">33</a>]. The asymptotic homogenization method has also been utilized for topology optimization under stress constraints to create lightweight structures with predictable mechanical responses [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib34" id="fnref-smsad6d1fbib34"="">34</a>]. Nguyen <em="">et al</em> [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib35" id="fnref-smsad6d1fbib35"="">35</a>] developed a robust implicit-based design framework for the aided design of structures with FGMM. However, the high computational cost in multiscale topology optimization remains a significant drawback [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib36" id="fnref-smsad6d1fbib36"="">36</a>]. Additional topology optimization-based methods aiming for low computational requirements have been proposed as alternative strategies for generating structures with FGMM [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib37" id="fnref-smsad6d1fbib37"="">37</a>, <a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib38" id="fnref-smsad6d1fbib38"="">38</a>].

Deep learning-assisted methods have become another trend in FGMM structures. For example, neural networks have been employed to obtain region-wise graded structures with specific responses by combining two variations of truss-based unit cells [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib39" id="fnref-smsad6d1fbib39"="">39</a>]. Additionally, a neural network with physical constraint inputs was developed to create truss structures with tailored mechanical responses [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib40" id="fnref-smsad6d1fbib40"="">40</a>]. Wilt <em="">et al</em> [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib41" id="fnref-smsad6d1fbib41"="">41</a>] generated region-wise two-dimensional graded structures with customized local Poisson's ratio using a neural network. Nonetheless, current artificial neural network-assisted methods have not yet shown effective results and require further exploration [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib42" id="fnref-smsad6d1fbib42"="">42</a>]. All the works mentioned above leverage AM technologies to validate their findings, benefiting from their capacity to manufacture highly complex and detailed geometries [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib43" id="fnref-smsad6d1fbib43"="">43</a>].

Alternative deep learning approaches have been explored for solving this type of problem. For instance, the deep energy method proposed by Samaniego <em="">et al</em> [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib44" id="fnref-smsad6d1fbib44"="">44</a>], where deep neural networks were used to approximate the solution to partial differential equations by using the energy of a mechanical system as the loss function. Even though analytical models for functionally graded materials have been proposed [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib45" id="fnref-smsad6d1fbib45"="">45</a>], limited research has focused on analytical design approaches for graded cellular metamaterials. Zhang <em="">et al</em> [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib46" id="fnref-smsad6d1fbib46"="">46</a>] developed different strategies to achieve tailored mechanical responses of tubular structures with graded mechanical metamaterial reinforcements. Good agreement between the analytical model and the experimental data was obtained, however, numerical results were not reported. For their optimization-based design method, high computational requirements were reported. Similarly, inverse design methods for phononic topological insulators based on hexagonal unit cells have been developed, employing level set topology optimization [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib47" id="fnref-smsad6d1fbib47"="">47</a>], and genetic algorithms [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib48" id="fnref-smsad6d1fbib48"="">48</a>].

The existing models, though effective, often come with high computational demands for inverse design, hindering rapid iteration of desired mechanical behaviors when required. Situations demanding rapid iteration of desired mechanical behaviors necessitate low computational cost and time. This work is part of an ongoing effort towards designing bending elements for their application in soft robotics. For this, a graded distribution of metamaterials under simple loading conditions is required, and employing high computational and robust models results unnecessary and inefficient.

Motivated by these needs, this work introduces a low computational cost, boundary conditions-specific inverse design method for segment-wise functionally graded beams, tailored to achieve specific transverse deflections under identical loading conditions. The design framework relies on a semi-analytical method, validated through Finite Element (FE) models and experimental tests conducted on additively manufactured samples. This approach serves as a proof-of-concept of using case-specific models for rapid design of solids with FGMM, that can be extrapolated to different requirements and structures.

The model's primary objective is to derive a beam featuring a functionally graded distribution of a chosen topology based on targeted deflection criteria, particularly focusing on cantilever beams subjected to bending. The method was assessed using three distinct honeycomb topologies: rectangular, reentrant, and hexagonal. Four different relative densities per topology were explored by modifying the wall thickness of the unit cells. Thermoplastic polyurethane (TPU) was chosen for its flexible properties, which serve as the material used. Mechanical properties of both the source material and the topologies were obtained through mechanical experimentation, serving as inputs for the FE models and the semi-analytical model respectively. The performed validation demonstrated agreement between the semi-analytical method and both the FE analysis and laboratory experimentation.

Our work introduces a novel inverse design method based on analytically synthesizing a segment-wise graded metamaterial beam. The results demonstrate the effectiveness and low computational cost of the method. The same proposed framework can be applied to any desired structure under specific loading and boundary conditions. The rapid prototyping nature of current AM technologies can leverage the proposed method to generate end-to-end fast design and manufacturing cycles.

The remainder of this document is structured as follows: section <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#smsad6d1fs2"="">2</a> describes the materials and methods used in this work. Section <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#smsad6d1fs3"="">3</a> outlines the semi-analytical model used for beam deflection calculations, along with the proposed inverse design methodology. Section <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#smsad6d1fs4"="">4</a> comprehensively describes the obtained results. Section <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#smsad6d1fs5"="">5</a> comprises the discussion of the results. Finally, section <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#smsad6d1fs6"="">6</a> presents the conclusions drawn from the work, alongside future work and potential applications.

The methodology followed in this work comprises 14 stages. These are depicted in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff1"="">1</a> and detailed in the upcoming subsections.

Zoom In Zoom Out Reset image size

The methodology can be briefly summarized as follows. The initial steps involve characterizing the source material utilized for the AM process and FE analysis. This is achieved through tension tests conducted on additively manufactured dogbone samples (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff1"="">1</a>(a)). Following this, parameterization of the topologies used (top-left step) was carried out. Subsequently, Computer-Aided Design (CAD) models of uniform beams were generated (middle-left) for FE analysis and laboratory experimentation using additively manufactured samples subjected to three-point bending tests (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff1"="">1</a>(b)). The comparison of these results serves to validate the FE models and determine the effective Young's Modulus of the topologies. Moreover, effective Young's Moduli are employed as input for the proposed semi-analytical model (bottom-middle). Additionally, from the earlier parameterization of the topologies, CAD models of randomly generated graded beams are derived (top-middle). These beams are then used to validate the semi-analytical model by conducting cantilever laboratory experiments on additively manufactured samples, as well as performing FE analysis under the same boundary conditions (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff1"="">1</a>(c)). Further steps for validating the developed inverse design model will be detailed in section <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#smsad6d1fs4-5"="">4.5</a>.

2.1. Parameterization of used topologies

The topologies studied here are illustrated in figure 2(a): rectangular, reentrant, and hexagonal honeycomb structures. These cellular structures consist of periodic distributions of the unit cells. Each unit cell is defined by the parameters h and t; h defines the dimensions of the unit cell, uniformly set to 2h on the x-axis and 3h on the z-axis. In the case of the reentrant and hexagonal honeycombs, the distance between inclined struts is also determined by h. For this study, a standardized value of 2.5 mm was assigned to h across all cases. Parameter t regulates the variation in wall thickness within the unit cells, thereby influencing their relative density ($\bar{\rho}$).

Zoom In Zoom Out Reset image size

These cellular structures represent two-dimensional extruded patterns, with their relative density defined as the ratio of the volume of the material making up the unit cell to the total volume of the unit cell's bounding box.

The equations required to calculate the relative densities of the topologies based on the parameters <em="">h</em> and <em="">t</em> can be found in the supporting information S1. Four relative density variations for all topologies were studied and their corresponding t values are listed in table <a xmlns:xlink="http://www.w3.org/1999/xlink" class="tabref" href="#smsad6d1ft1"="">1</a>.

<b="">Table 1.</b> Wall thickness values for each topology's four relative density variations.

	Rectangular	<em="">t</em> (mm) Reentrant	Hexagonal
$\bar{\rho}_{24} \approx 0.24$	0.8	0.4	0.6
$\bar{\rho}_{36} \approx 0.36$	1.2	0.8	0.8
$\bar{\rho}_{49} \approx 0.49$	1.6	1	1.2
$\bar{\rho}_{70} \approx 0.70$	2.6	1.6	2

2.2. Design of beams with graded metamaterials

The semi-analytical method, FE analysis, and laboratory experimentation were conducted on graded beams consisting of six segments (also known as layers), each with the same topology but different values of $\bar{\rho}$. An example using the hexagonal honeycomb is shown in figure 2(b). The CAD models for the graded beams were created using SOLIDWORKS^® 2023 (Dassault Systems, Vèlizy, France). The designated value of $\bar{\rho}$ for each segment was specified, and using programmed macros, the solid models were generated and exported. This tool facilitated an automated method for generating and exporting the CAD models intended for subsequent AM and FE analysis.

Each segment comprises four unit cells in both the <em="">x</em>- and <em="">z</em>-axes, resulting in a beam measuring 120 mm in the <em="">x</em>-axis and 30 mm in the <em="">z</em>-axis. An additional 25 mm long extension for fixture purposes was incorporated on the left-most side of the beam, as highlighted in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff2"="">2</a>(b). This extension will be secured using three screws in a fixture designed for the subsequent cantilever tests detailed in section <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#smsad6d1fs2-5"="">2.5</a>. The beam is extruded 4 mm in the <em="">y</em>-axis.

As highlighted in section 1, achieving seamless transitions between unit cells with differences in design parameters is crucial in the development of FGMM. Here, transition regions have been incorporated between the graded segments. These transitions between values of $\bar{\rho}$ involve modifying the thickness of neighboring elements. Hence, half of the unit cells in contact with the following segment will be modified to meet a transitional thickness. Figure 3 depicts the transitioning between $\bar{\rho}_{24}$ to $\bar{\rho}_{49}$ for the three topologies studied. Note that while four unit cells are depicted, the adjacent elements of the two unit cells in contact with the transition region continuously alter their wall thickness to match the thickness of the subsequent segment. This automated process was executed for any desired shift in $\bar{\rho}$.

Zoom In Zoom Out Reset image size

Considering the four selected values of $\bar{\rho}$ and the three studied topologies, a total of 12 different uniform beams were designed for the characterization of their effective Young's modulus. On the other hand, for these graded beams with six segments where different values of $\bar{\rho}$ can be assigned, a total of 4096 possible combinations exist. For validation purposes, 20 combinations were randomly chosen for FE analysis per topology, and from these, five were selected to be manufactured and subjected to experimental testing.

2.3. Additive manufacturing of cellular beams

2.4. Finite element models, meshing, and boundary conditions

Three-point bending and cantilever tests were used in both FE analysis and laboratory experiments to evaluate the mechanical properties of the topologies and to measure the transverse deflection of the beams under flexural loads. Structural static analyses were performed with ANSYS<sup="">®</sup> Workbench 2021 (ANSYS Inc. Canonsburg, United States).

All models were meshed using tetrahedral elements of a maximum of 2 mm. The equivalent stress mesh refinement convergence tool was used, with a 5% of allowable change. Computational models were fed with the material properties of Ultimaker<sup="">®</sup> TPU 95A filament. The Young's modulus was obtained from 7 tensile tests of dogbone ASTM d638 additively manufactured samples. The tensile tests were performed in a Shimadzu<sup="">®</sup> AGX-V2 universal testing machine (Kioto, Japan) with a strain rate of 0.1 mm/s, obtaining a Young's Modulus of 57.4762 ± 4.5 Mpa.

Three-point bending simulations were performed with the boundary conditions depicted in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff4"="">4</a>(a). Dimensions of the supporting and loading pins were set to meet the setup used for experimentation, which will be described in section <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#smsad6d1fs2-5"="">2.5</a>. Material properties of structural steel were used for the supporting and loading pins. Between the beam, and the supporting and loading pins, a friction contact with a coefficient of 0.2. A remote displacement in the <em="">y</em>-axis of 10 mm was applied to the loading pin, and the force reaction generated in the supporting pins was measured.

Zoom In Zoom Out Reset image size

Cellular beams under cantilever loading conditions were also simulated (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff4"="">4</a>(b)). The fixture extension was fully fixed, restricting all degrees of freedom. The loading pin was placed at the free end of the beam, and 0.1 N in the negative <em="">y</em>-axis was applied. Deflection of the graded beams was measured with the displacement probe tool, measuring <em="">y</em>-axis displacements in 6 different positions as indicated in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff4"="">4</a>(b).

2.5. Characterization via flexural tests on additively manufactured samples

2.6. Statistical analysis

As discussed in section <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#smsad6d1fs1"="">1</a>, existing design tools for the inverse design of FGMM are constrained by their computationally demanding resources. In response to the need for a low-computation-cost model to facilitate rapid generation and analysis of structures with FGMM, a semi-analytical model was derived based on Castigliano's second theorem.

The theorem states that <em="">when forces act on a system subject to small displacements, the displacement corresponding to the acting force is equal to the partial derivative of the total strain energy with respect to the force</em> [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib49" id="fnref-smsad6d1fbib49"="">49</a>]. This theorem can be used to find the deflection at points where no force is applied; the theorem used for beams under bending forces can be expressed as,

$$\begin{equation} \delta = \frac{\partial U}{\partial F} = \int \frac{1}{E I} \left(M \frac{\partial M}{\partial F}\right){\mathrm{d}}x\end{equation} \tag{ 1 }$$

where <em="">δ</em> represents the deflection of the beam, <em="">U</em> is the total energy strain, <em="">F</em> is the applied force, <em="">I</em> represents the inertial moment, <em="">M</em> is the internal moment, and <em="">E</em> is the Young's modulus of the constituent material.

Consider now a beam as the one shown in figure 5(a), it has been discretized into s<sub="">i</sub> segments of length l<sub="">i</sub> (i going from 1 to 6), each with different stiffness E<sub="">i</sub>. In the case of FGMM beams, the stiffness of each segment will be defined by adjusting its corresponding $\bar{\rho}$.

Zoom In Zoom Out Reset image size

To analyze the deflection of the beam at any given i point, we will consider a F force applied at the free end of the beam, that corresponds to the case scenario studied here. Additionally, a fictitious force labeled Q equivalent to 0 is located at the end of each segment i in the position l<sub="">i</sub> (see figure 5(a). Hence, the value of M can be obtain with $M = -(F+Q)x$. Considering Q = 0, then

$$\begin{equation} M\frac{\partial M}{\partial Q} = -\left(F+Q\right)\left(x\right)\left(-x\right) = {F}{x}^{2}.\end{equation} \tag{ 2 }$$

Substituting equation (<a xmlns:xlink="http://www.w3.org/1999/xlink" class="eqnref" href="#smsad6d1feqn2"="">2</a>) in equation (<a xmlns:xlink="http://www.w3.org/1999/xlink" class="eqnref" href="#smsad6d1feqn1"="">1</a>), we can obtain the deflection <em="">δ</em><sub=""><em="">i</em></sub> of the <em="">i</em>th segment as

$$\begin{equation} \delta_i = {\int_{l_i}^{l_{i+1}}{\frac{1}{E_i\ I}\left(Fx^2\right)}{\mathrm{d}}x}.\end{equation} \tag{ 3 }$$

By solving equation (3), we can obtain the deflection at any jth location of the beam if $0 \unicode{x2A7D} j \unicode{x2A7D} n$ as

$$\begin{equation} \delta_j = \sum_{m\ = \ 1}^{j}{\frac{F\left({l_{m+1}}^3-{l_m}^3\right)}{3{\ E}_m\ I}\ }.\end{equation} \tag{ 4 }$$

However, equation (<a xmlns:xlink="http://www.w3.org/1999/xlink" class="eqnref" href="#smsad6d1feqn4"="">4</a>) does not consider the slope continuity between segments. Beams periodically loaded or resulting in periodic displacements (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff5"="">5</a>(b)), should also consider the slope-continuity [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib50" id="fnref-smsad6d1fbib50"="">50</a>]. Hence, applying the slope-deflection condition, the slope <em="">ψ</em><sub=""><em="">i</em></sub> can be determined as

$$\begin{equation} \psi_i = \frac{{\mathrm{d}}\delta_j}{{\mathrm{d}}x}\ \ = \frac{\delta_j-\ \delta_{j-1}}{x_j-x_{j-1}},\end{equation} \tag{ 5 }$$

where x<sub="">j</sub> is a distance from the origin to a j position, considering that j is within the ith segment. By applying equations (4) and (5), the total deflection $\Delta_j$ of the beams in the j position considering the slope continuity can be calculated as

$$\begin{equation} \Delta_j\ = \delta_j + {\mathrm{d}}\psi_i\ \left(s_i\right),\end{equation} \tag{ 6 }$$

where $d \psi_i$ represents the change in slope in the ith segment, and s<sub="">i</sub> represents the length of the ith segment. By applying equation (6) at any desired point, the localized deflection of a beam with segments of different material properties can be obtained. To achieve an inverse design for beams with functionally graded metamaterials, we need a method that enables the generation of a design based on the desired beam deflection, as outlined in section 2.2. This design method will be described in the upcoming subsection.

3.1. Inverse design method for beams with functionally graded metamaterials

The design method proposed here is based on the workflow shown in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff6"="">6</a>. The 4096 possible combinations for each topology were subjected to an <em="">F</em> load of 0.10<em="">N</em> and their transverse deflections were estimated with a MATLAB<sup="">®</sup> R2023b (MathWorks, Massachusetts, United States) script. The results were stored in a data structure for each topology containing the deflection of the graded beams at any <em="">j</em> position for all possible distributions.

Zoom In Zoom Out Reset image size

To compare laboratory experimentation, FE analysis, and the semi-analytical method, the six values of Δ indicated in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff5"="">5</a>(b) (which correspond to the marker indicated in figures <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff4"="">4</a>(b) and (d) were used as points for measuring deflection in all the methods. These points coincide with the transition regions between the grading segments.

Figure 6 summarizes the inverse design framework used to determine the adequate $\bar{\rho}$ distribution for a desired transverse deflection. The desired deflections are used as inputs of the design framework. The values are then compared with the stored deflections of the 4096 possible combinations. To determine the best fit between the grading variations, the mean absolute error (MAE) of each variation is obtained as

$$\begin{equation} {\mathrm{MAE}} = \frac{1}{6} \sum_{i\ = \ 1}^6 |\Delta_i - \hat{\Delta}_i|,\end{equation} \tag{ 7 }$$

where $\Delta_i$ represents the deflection at the six selected positions along the beam. The grading variation with the least MAE when compared with the targeted deflections is considered the best fit. The grading variation output consists of six $\bar{\rho}$ values that represent the functionally graded beam design. The MATLAB^® script output will be then input to the previously described SOLIDWORKS^® macro to generate and export the functionally graded beam CAD model.

4.1. Additively manufactured samples: variations and microscopic inspection

Manufacturing defects can lead to significant discrepancies between CAD models and manufactured samples, resulting in deviations between FE and experimentally obtained data [51]. To mitigate this, the additively manufactured samples underwent microscopic inspection to assess manufacturing defects. The inspection was conducted using a Zeiss^® Axio Vert.A1 (Oberkochen, Germany) microscope. Micrographs from the inspection can be found in figure S2 of the supporting information, and measurements of wall thickness were compared with the expected values from table 1. The average results of the measurements are reported in table S2, supporting information. Six different measurements per topology and relative density variation were performed. The standard deviations range from 4.6% in the reentrant unit cells with $\bar{\rho_{24}}$ to 8.4% in the hexagonal honeycomb samples of $\bar{\rho_{70}}$, with an average standard deviation of 6.4% for all possible variations. Overall, larger deviations were found in samples with higher relative densities.

Figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff7"="">7</a> shows images of the additively manufactured samples compared with their CAD versions. These comparisons prompt a correction of the CAD models used for FE analysis to have a closer representation of the physical samples when performing computational experimentation.

Zoom In Zoom Out Reset image size

As described in section 2.2, 20 out of the 4096 possible combinations of the graded beams were randomly selected per topology for FE analysis. From these, 5 combinations were selected for AM and laboratory experimentation. Figure 8 depicts the 5 selected combinations of $\bar{\rho}$ for beams with rectangular (figure 8(a)), reentrant (figure 8(b)), and hexagonal (figure 8(c)) honeycombs. These illustrations only contain one of the three replicates manufactured per each variation.

Zoom In Zoom Out Reset image size

4.2. Effective stiffness: computational and experimental results

The effective Young's modulus of the three analyzed topologies and their four $\bar{\rho}$ variations were studied to feed the semi-analytical model. These were measured through computational and laboratory experimentation as described in sections 2.4 and 2.5. Figure 9 includes Young's Moduli obtained from the three-point bending tests. The stiffness of the structures was measured from the experimental stress–strain graphs, with a maximum strain of 10%. Error bars for the three replicates indicating the +/- standard deviation of the additively manufactured beams were included. Low variability between experimental samples was observed. The stiffness measured from the computational and laboratory experimentation can be seen in table S3, supporting information.

Zoom In Zoom Out Reset image size

A cubic proportional relation between relative density and stiffness ($E \propto \bar{\rho}^{\ 3}$) is observed for the reentrant and hexagonal honeycomb topologies. Whereas the rectangular honeycomb exhibits a linear proportional relation ($E \propto \bar{\rho}$), clearly noticed in the $\bar{\rho} \approx 0.7$ variation. A higher deviation between computational and laboratory experimentation was found in samples with higher relative densities. This is attributed to a higher number of manufacturing defects in denser samples, as reported in section 4.1.

The stiffness results were used to feed the semi-analytical model as described in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff1"="">1</a>. This is a key step in the construction of the semi-analytical model and in the validation of the laboratory and computational methods implemented.

4.3. Experimental cantilever tests results

The flexural cantilever tests described in section <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#smsad6d1fs2-5"="">2.5</a> were performed to validate the semi-analytical model described in section <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#smsad6d1fs3"="">3</a>. The three additively manufactured replicates of the variations of graded beams were subjected to these loading conditions. Figures <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff10"="">10</a>(a)–(c) include the force-displacement graphs of the performed flexural tests, the standard deviations of every combination are marked with a shading. Low variability between manufactured samples was found. The grading combinations results correspond to the beams shown in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff8"="">8</a>.

Zoom In Zoom Out Reset image size

For the variations of the rectangular honeycomb (figure 10(a)), the beams presented stiffness ranging from 0.005 N mm⁻¹ to 0.01 N mm⁻¹. The average relative density of all the segments of these samples is 0.338 and 0.587, respectively. For the grading variations with the reentrant honeycomb (figure 10(b)), combinations four and five achieved stiffness of up to 0.0133 N mm⁻¹, while the remaining resulted with stiffness between 0.006 N mm⁻¹ and 0.0073 N mm⁻¹. Once again, the combination with the lowest average relative density ($\bar{\rho} = 0.33$) presented the lowest stiffness, and the highest one was produced by the combination with an average relative density of 0.482.

In the beams composed of graded hexagonal honeycombs (figure 10(c)), four variations had similar stiffness between 0.004 N mm⁻¹ and 0.005 N mm⁻¹, with the least stiffness being obtained from the samples with an average relative density of $\bar{\rho} = 0.40$. The fifth combination had the highest stiffness with 0.0073 N mm⁻¹, obtained from the samples with average relative densities of 0.42. The difference in stiffness between these two groups is attributed to the positioning of the layer with the highest densities, as the first segment of the combination with the highest stiffness has a $\bar{\rho}_{70}$, while the less stiff variation has a segment with $\bar{\rho}_{24}$.

Note that these displacement measurements are obtained at the free end of the beams. For the desired analysis, all beams are studied when exposed to 0.10 N to keep within the linear behavior, and their deflections were measured in the six positions indicated in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff4"="">4</a>(d). These measurements were obtained through the computer vision-based analysis described in section <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#smsad6d1fs2-5"="">2.5</a>.

Figures <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff10"="">10</a>(d)–(f) illustrate the transverse deflection of additively manufactured graded beams with rectangular (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff10"="">10</a>(d)), reentrant (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff10"="">10</a>(e)), and hexagonal (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff10"="">10</a>(f)) honeycombs. Minimal deviation was observed between samples of the same variation.

The distribution of relative densities along the beams influences the deflection modes. To illustrate, hexagonal honeycomb variation 5 ($\bar{\rho}_{70}, \bar{\rho}_{49}, \bar{\rho}_{24}, \bar{\rho}_{49}, \bar{\rho}_{36}, \bar{\rho}_{24}$) from figure 10(f) shows deflections further along the beam, at 0.03 m. In contrast, variation 2 ($\bar{\rho}_{24}, \bar{\rho}_{36}, \bar{\rho}_{36}, \bar{\rho}_{70}, \bar{\rho}_{49}, \bar{\rho}_{36}$) deforms rapidly from the origin. These deformation mechanisms are attributed to the distribution of relative densities in the diverse segments, shifting the deflection curve depending on the $\bar{\rho}$ of each segment. This tunability of deflection mechanisms is desirable for tailor-made deflection of beams.

4.4. Transverse deflections of graded beams: reconciliation between computational, laboratory and analytically obtained data

Computational and laboratory experimentation were conducted to validate and compare transverse deflection values obtained through the semi-analytical method. Following the procedures outlined in section <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#smsad6d1fs2-4"="">2.4</a> for computational experimentation, the same boundary conditions and loads as the experimental setup were applied. Additionally, the semi-analytical method (described in section <a xmlns:xlink="http://www.w3.org/1999/xlink" class="secref" href="#smsad6d1fs3"="">3</a>) was used to estimate the transverse deflections. Transverse deflection data from the three approaches is plotted in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff11"="">11</a>. Two grading variations per topology are displayed to exemplify the deflection comparison. The deflection measurements of the three approaches are depicted, and the three replicates of experimental laboratory tests are shown. In all cases, experimental tests exhibit low deviation between samples.

Zoom In Zoom Out Reset image size

For the rectangular honeycomb, the deflections of variations ($\bar{\rho}_{36}, \bar{\rho}_{36}, \bar{\rho}_{36}, \bar{\rho}_{70}, \bar{\rho}_{24}, \bar{\rho}_{24}$) and ($\bar{\rho}_{49}, \bar{\rho}_{70}, \bar{\rho}_{24}, \bar{\rho}_{24}, \bar{\rho}_{24}, \bar{\rho}_{36}$) are shown in figures 11(a) and (d), respectively. Agreement between the semi-analytical and computational methods is observed, with a higher deviation from the experimental tests. For the reentrant honeycombs, the variations ($\bar{\rho}_{24}, \bar{\rho}_{70}, \bar{\rho}_{24}, \bar{\rho}_{49}, \bar{\rho}_{24}, \bar{\rho}_{49}$) and ($\bar{\rho}_{70}, \bar{\rho}_{49}, \bar{\rho}_{36}, \bar{\rho}_{49}, \bar{\rho}_{36}, \bar{\rho}_{49}$) are displayed in figures 11(b) and (e), accordingly. Slight deviations can be observed between the three approaches, mainly between the computational and semi-analytically obtained data. The hexagonal honeycomb variations ($\bar{\rho}_{36}, \bar{\rho}_{49}, \bar{\rho}_{24}, \bar{\rho}_{24}, \bar{\rho}_{49}, \bar{\rho}_{24}$) and ($\bar{\rho}_{24}, \bar{\rho}_{49}, \bar{\rho}_{49}, \bar{\rho}_{24}, \bar{\rho}_{24}, \bar{\rho}_{36}$) and their deflections with the different methods are depicted in figures 11(c) and (f), respectively, showing low deviations between approaches. Larger deviations are found mostly between the semi-analytical method and the computational and laboratory results. When comparing the variations of the three topologies, it can be observed that higher differences in the deflection measurements at the free end are found. Additionally, these deviations are greater in the beams composed of rectangular honeycombs, followed by the reentrant and the hexagonal honeycombs. Note that the study focuses on assessing transverse deflections, therefore, displacements along the x-axis were not measured in any of the methods.

Considering the validation depicted in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff1"="">1</a>(c), semi-analytical, computational, and laboratory experimentation was conducted for the selected variations. Transverse deflection measurements were obtained for all the variations and methods. Figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff12"="">12</a> illustrates the deflection discrepancies between the semi-analytical method and the computational and laboratory experimentation in the six positions where deflections were measured. Recall that per topology, five variations were experimentally tested, and 20 were subjected to FE analysis. Figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff12"="">12</a>(a) includes the deflection differences for the rectangular honeycomb variations. Laboratory tests resulted in more differences than the computational data, increasing as the measurements were taken further in the <em="">x</em>- positive direction. Negative differences increase as they approach the free end of the beams for the variations composed of reentrant honeycombs (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff12"="">12</a>(b)). Rectangular honeycomb beams show the highest deviations among the three topologies. The deflections of hexagonal honeycomb beams are depicted in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff12"="">12</a>(c), where fewer deviations can be observed compared to the other topologies.

Zoom In Zoom Out Reset image size

For the three topologies, deflection variations increase proportionally along the <em="">x</em>-axis. This is attributed to angle changes in the early segments of the beams, causing high increments in the later segments, thereby increasing the deflection deviations between the methods. Higher deviations can be observed when comparing the semi-analytical method with the laboratory experiments.

The MATLAB^® script that estimates the deflection calculations reports an elapsed time per variation of $1.29 \times 10^{-4}$ s when running on a quad-core laptop. However, to calculate the deflection of the 4096 possible combinations of the graded segments, an average of 0.52 s was required per topology. Considering the requirement of performing the calculation for the three topologies, a total of approximately 2 s was required, taking into account additional calculations performed for organizing the arrays with the resulting data. There are no specific computational requirements for calculating the deflections since only algebraic formulas are needed when numerically solving the integrals from the model.

4.5. Inverse design method: case studies

The proposed inverse design method allows obtaining the $\bar{\rho}$ distributions based on targeted deflections of beams with graded topologies. These desired deflections must fall within the flexural capacities of the beams. To validate the inverse design approach, a targeted deflection for each topology was fed as input into the method. As depicted in figure 6, the output was a combination of densities that could achieve the targeted deflections estimated semi-analytically. The targeted deflections were chosen by assessing the design space for beams composed of each topology and selecting an allowable set of deflections per segment. The proposed beam variations were then exported as CAD models to perform FE analysis and compare the deflections of these approaches.

The targeted deflections used for the case studies and the inverse design approaches are displayed in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff13"="">13</a>. The numerical results are found in table S4, supporting information. Additionally, deflections from the FE models were measured to assess the semi-analytical data. Note that the grading combinations were proposed by the semi-analytical inverse method, and not randomly chosen combinations like the variations studied in the previous section. Slight deviations between the semi-analytical data and computational experimentation can be observed. Overall, higher deflections were obtained from the FE models. The beam composed of rectangular honeycombs (figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff13"="">13</a>(a)) presents the lowest deviations between targeted deflection, semi-analytical, and computational results.

Zoom In Zoom Out Reset image size

The targeted deflections represent performance requirements, while the outputs of the inverse design method are the material arrangements needed in the beams to achieve them. For instance, in figure 13(a), the targeted deflection required a higher deflection at the free end of the beam, resulting in the variation ($\bar{\rho}_{36}, \bar{\rho}_{36}, \bar{\rho}_{70}, \bar{\rho}_{70}, \bar{\rho}_{24}, \bar{\rho}_{24}$), where the segments with less relative densities are found at the free end, thus, making this segment more flexible. For the deflection in the case study depicted in figure 13(b), high deflections were required in the first segments of the beam, with the remaining segments deflecting in a straight line. This was accomplished by the inverse design approach by arranging low $\bar{\rho}$ values in the first segments and higher ones in the second half of the beam, resulting in the combination ($\bar{\rho}_{24}, \bar{\rho}_{23}, \bar{\rho}_{24}, \bar{\rho}_{70}, \bar{\rho}_{70}, \bar{\rho}_{70}$). In the case study illustrated in figure 13(c), the slope of the deflection increments proportionally to the x-axis. This was achieved with the variation ($\bar{\rho}_{24}, \bar{\rho}_{24}, \bar{\rho}_{36}, \bar{\rho}_{36}, \bar{\rho}_{49}, \bar{\rho}_{70}$) by gradually incrementing the relative density of the segments.

Equivalent Von Misses stress contours were obtained from the analyzed beams in cantilever tests in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff14"="">14</a>. It can be observed that stress is concentrated in less dense segments, with almost no stress located in segments with high densities. The different stress distributions are notable when comparing the different topologies. For instance, in figure <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff14"="">14</a>(a) a beam composed of rectangular honeycombs is displayed. Here, the stresses are only exhibited in the longitudinal struts. In contrast, the beam displayed in figures <a xmlns:xlink="http://www.w3.org/1999/xlink" href="#smsad6d1ff14"="">14</a>(b) and (c) are composed of reentrant and hexagonal honeycombs, and stresses are located in the angled struts along the <em="">x</em>-axis, as well as some regions of the transversal elements.

Zoom In Zoom Out Reset image size

Data for the effective stiffness obtained from computational and laboratory experimentation resulted in agreement. Figure 9 shows that the effective Young's modulus scales with $\bar{\rho}$. Two distinct effective Young's modulus—relative density relations can be obtained, which are attributed to the dominating deformation mechanisms of the cellular materials [52]. A cubic proportional relationship between $\bar{\rho}$ and Young's modulus was observed for the reentrant and hexagonal honeycombs, consistent with the Gibson–Ashby model [53] for bending-dominated cellular materials. In contrast, the rectangular lattice demonstrated a linear response, as expected in stretch-dominated structures.

Additionally, for some beams with certain relative densities, the Young's modulus of experimental and computational results shows significant statistical differences, as noted in table S3. This is attributed to two main factors. First, the influence of topology on the effective properties of cellular materials diminishes as relative density increases, transitioning towards less porous media. In these cases, structures resemble those with incidental pores rather than engineered arrangements with well-defined geometric topologies. Consequently, such structures are more prone to manufacturing defects, such as over-extrusion and inner voids, due to the limitations of FFF. These defects notably impact the overall mechanical performance of high-density samples.

When using the semi-analytical method, the topology selection will depend on the mechanical performance requirements and the design space that the beams of each topology allow. The selection can be assisted with the guideline proposed by Arredondo–Soto <em="">et al</em> [<a xmlns:xlink="http://www.w3.org/1999/xlink" class="cite" href="#smsad6d1fbib54" id="fnref-smsad6d1fbib54"="">54</a>] for tailoring the stiffness of compliant systems through the use of metamaterials. This is based on iteratively modifying the selected topology or its design parameters based on stiffness and auxetic requirements.

The stiffness of the beams can be modulated segment-wise, achieving tailored deflections based on the $\bar{\rho}$ values of each segment. Consequently, beams with the same average $\bar{\rho}$ but distinct distributions yield different deflection curves, emphasizing the importance of graded distributions in achieving specific responses. This modulation is crucial for crafting customized beams as integral components of systems that demand precise deformation control in flexible mechanisms. Fields such as soft robotics and compliant mechanisms require this level of refined adjustment. For example, Cáceres et al [55] explored beams composed of metamaterials as components in compliant joints, utilizing neural networks for the parameter selection. Similarly, Mohammadi et al [56] investigated graded distributions of reentrant honeycombs as compliant joints in a bio-inspired soft robotic hand. These examples demonstrate the potential of using graded distributions of metamaterials to modulate stiffness, and resulting deformed shapes, in soft robotic elements and components.

The validation of the semi-analytical model exhibited differences between the deflections measured by the model, the computational, and laboratory experimentation. The differences were found to increase closer to the free end of the beams; this is attributed to the incremental sum of errors from previous segments, creating an angular error that causes larger differences. The mean deviations between the methods in the different measurement positions range from 4% to 16% with an average of 11.6% for the beams composed of rectangular honeycombs, from 1% to 15% with an average of 9% for the reentrant honeycombs, and for the hexagonal honeycombs, a range from 3% to 16% with an average of 10%.

The results suggest that the approached discretization of the structure as a continuous beam with six distinct material properties affects the accuracy of predicted deflections. A factor not accounted for in the semi-analytical model is the transitioning regions between segments with different relative densities, a relevant feature in graded structures. Furthermore, as the material used is a flexible polymer, gravity could play a role in the deformations observed in the performed experimentation, while gravity was not considered in the other explored methods. These omissions lead to variations in mechanical behavior, further exacerbating differences between the analyzed experimental and computational methods. Additionally, manufacturing defects in additively manufactured samples impact the mechanical responses of the beams. Significantly, fewer discrepancies were observed between the semi-analytical model and computational experimentation, revealing the repercussions of manufacturing errors.

Experimental measurements may also be influenced by the vision-based tracking method, with factors like camera positioning, lighting, and calibration constraining data accuracy. The observed deviations among deflection calculation methods impact the model's reliability when high accuracies are required. Nonetheless, the semi-analytical model exhibits a consistent direction, curvature, and similar response despite the reported discrepancies, underscoring its efficiency among the various studied cases.

The beams composed of rectangular honeycombs exhibited higher deviations between laboratory experimentation and the semi-analytical method. This is attributed to the composing struts of the additively manufactured samples, which include higher wall thickness values. The inclusion of thicker additively manufactured elements increases the presence of defects such as gaps inside the walls and rasters that affect the geometries. Furthermore, the stress distributions of the rectangular honeycombs focus exclusively on the longitudinal struts, therefore, having a higher stress concentration that intensifies the deviations when the manufactured samples do not match the exact structures. It is noteworthy that only two struts compose each unit cell of the rectangular honeycombs, in contrast to the reentrant or hexagonal honeycombs, where 6 struts with less wall thickness are found.

The application of Castigliano's second theorem facilitates the use of the same segment-wise grading solution for deflection calculations in scenarios involving alternative loading or boundary conditions applied to beams. As such, the investigated cantilever setting serves as a proof-of-concept, and diverse configurations, such as coupled forces and moments, can be explored. It is crucial to emphasize that Castigliano's second theorem is applicable only for small deflections, thereby restricting the model's utility to linear responses and small deflections. The method under evaluation is specifically focused on deflection calculation and cannot be employed for stress concentration or other additional mechanical properties.

Beams with targeted deflections could be achieved using the proposed inverse design method. The explored case studies presented low variation between the desired behavior and that obtained by a proposed graded beam with segment-wise $\bar{\rho}$ variations. This functionally graded strategy enables the user to generate a design based on the mechanical requirements. A similar strategy can be implemented by utilizing low computational models to obtain tailor-made distributions of cellular materials in different structures or loading conditions.

For developing the inverse design method, the <em="">a priori</em> calculation of deflections of the 4096 possible combinations is needed. For this, an elapsed time of 2 s is needed when running the MATLAB<sup="">®</sup> script on a quad-core laptop. This allows the user to rapidly iterate and modify not only the dimensions and design requirements but also the mechanical properties of the metamaterials used, in case <em="">a prior</em> mechanical characterization has been performed.

Beams with segment-wise graded metamaterials were proposed to analyze their transverse deflection under cantilever tests using FE analysis, laboratory experimentation on additively manufactured samples, and a proposed low computational cost semi-analytical method. The semi-analytical model successfully predicted the transverse deflection of the graded beams under cantilever flexural tests with some differences between the measurements. These differences are attributed mainly to manufacturing defects in the additively manufactured samples and the non-consideration of transitioning regions between graded segments in the model. Additionally, the model demonstrated very low computational cost, achieving the computational efficiency sought.

The inverse design approach demonstrated an accurate design of beams based on deflection performance requirements. This allows the user to rapidly iterate the grading based on the design requirements. Additionally, by implementing the described pipeline, solids with the desired topologies can be swiftly generated for FE modeling and manufacturing. This model can be used as a design method for elements of compliant mechanisms.

This study opens avenues for future research in computational low-cost models for FGMM. Non-linear analytical models can enhance the analysis of a broader range of deformations, and exploring different boundary conditions and loading cases can extend the model's utility. Incorporating additional equations for calculating various mechanical responses can strengthen analytical inverse design models. The proposed semi-analytical model allows for the inclusion of a wider range of lattice types, such as metamaterials with curved elements or chiral structures. Including more parameters to modify the topologies could also allow for finer tuning of the mechanical properties of structures. Evaluating diverse materials and additive manufacturing techniques can further broaden the design scope of compliant structures, achieving more complex structures, or encompassing diverse materials such as rubbers or even metals. Grading in other directions could lead to widening the possibilities of tuning the desired behavior of engineered structures. For instance, by varying the thickness along the <em="">y</em>-axis, diverse deformation modes can be triggered. Similarly, more refined gradings, such as region-wise or cell-wise distributions, could achieve more complex deformation modes, albeit with a proportional increase in the complexity of the required inverse design models.

This work is part of an ongoing project focused on the implementation of mechanical metamaterials in the reinforcement of soft actuators. The proposed inverse design method allows us to generate metamaterial elements that meet a designated mechanical response. Rapid design iteration is required due to the targeted modularity sought from these compliant systems, where the metamaterial structure can be substituted depending on the deformation needed under controlled loading conditions.

We acknowledge the support from the School of Engineering and Sciences at Tecnologico de Monterrey Campus Querétaro and CONAHCYT for the MSc scholarship of the leading author (CVU: 1317862). The assistance provided by Dr Tejada-Ortigoza in the statistical analysis is also acknowledged.

The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors.

The authors have no conflicts of interest to declare that are relevant to the content of this article.