Analysis of a Preliminary Design Approach for Conformal Lattice Structures

: An important issue when designing conformal lattice structures is the geometric modeling and prediction of mechanical properties. This paper presents suitable methods for obtaining optimized conformal lattice structures and validating them without the need for high computational power and time, enabling the designer to have quick feedback in the ﬁrst design phases. A wireframe modeling method based on non-uniform rational basis spline (NURBS) free-form deformation (FFD) that allows conforming a regular lattice structure inside a design space is presented. Next, a previously proposed size optimization method is adopted for optimizing the cross-sections of lattice structures. Finally, two different commercial ﬁnite element software are involved for the validation of the results, based on Euler–Bernoulli and Timoshenko beam theories. The ﬁndings highlight the adaptability of the NURBS-FFD modeling approach and the reliability of the size optimization method, especially in stretching-dominated cell topologies and load conditions. At the same time, the limitation of the structural beam analysis when dealing with thick beams is noted. Moreover, the behavior of different kinds of lattices was investigated.


Introduction
Nowadays, the development of additive manufacturing (AM) technologies increases the freedom in the production of parts with a complex shape, enabling the designer to overcome the limitations imposed by traditional manufacturing technologies. Indeed, the conventional product design evolves toward the design for additive manufacturing (DfAM) concept [1]. Thanks to the additive processes, which are based on a layer-bylayer manufacturing approach, in DfAM, it is possible to exploit four forms of design complexity: shape complexity, functional complexity, hierarchical complexity, and material complexity [2].
Several design methods are developed for these goals, mostly inspired by nature, where each approach is the result of millions of years of evolution. Among them, some material distribution strategies mimic the mineral structures and the animals' skeletons. They are called cellular structures and are obtained by the arrangement and the repetition of a unit cell in the space [3]. Lattice structures present interesting properties: lightweight, high strength, energy absorption, and vibration reduction [4][5][6][7][8]. These properties are directly where b is the number of struts, j is the number of locked joints, s is the number of states of self-stress, and m is the number of mechanisms. According to the Maxwell criterion, if M is lower than 12, the 3D cell is bending-dominated; otherwise, it is stretching-dominated, and it is redundant from a rigidity point of view. However, the mechanical behavior of the lattice also depends on the unit cell topology, orientation, periodicity, and the load direction; for these reasons, it is not possible to determine the deformation mechanism directly from the Maxwell criterion. For example, a lattice structure composed of simple cubic cells is defined as bending-dominated but, when loaded along with the directions of its main axes, it can exhibit stretching-dominated behavior [5]. Lastly, a fundamental parameter of lattice structures is the relative density, defined as the ratio between the volume of the cellular structure and the volume of the design space. With their model [18], Gibson and Ashby demonstrated that it is possible to correlate the lattice mechanical proprieties, such as Young's modulus, to the relative density. Indeed, material selection charts that are related to those properties are available in the literature for the designer as a tool for the selection of the desired relative density and cell [4,5]. Tobias et al. [19] compared a significant amount of data of experimental compression tests on 3D-printed lattices available in the literature with the predictions of the Gibson-Ashby model. The findings highlight a suitable correlation between the experimental tests and the mathematical model; anyway, a precise prediction of these properties is not possible, so the model is intended to be used in the early design stages. Kaur et al. [20] studied the deformation of 3D-printed lattices, obtained by various thermoplastic materials and tested by compression testing and FEA; they found that the material can influence the bending and stretching-dominated behavior of the lattice structures.
Moreover, the cell element size can be optimized to satisfy structural requirements. Size optimization, together with the shape and topology optimization, is part of the structural optimization and aims to iteratively find the best configuration of the size of the elements, both beams diameter and shell thickness, to enhance the structural properties of the structure [21][22][23]. Once the lattice structure has been optimized, the information is gathered, and the geometrical model has to be obtained. Even though some lattice design tools were integrated inside commercial computer-aided design (CAD) software [24][25][26], the lack of suitable modeling methods for such structures has been highlighted in the literature [27][28][29]. Traditional CAD software are usually based on boundary representation (B-Rep) modeling approaches through function-based methods, such as non-uniform rational B-splines (NURBS) or discrete mesh methods [30]. Modeling conformal lattice structures inside these environments is challenging. Boolean, offsetting, and filleting operations on complex shapes require high computational resources and time, and, even worse, they often fail; furthermore, lattice structures with a high number of elements are difficult to manage and visualize. One more issue comes from the file format when exchanging data between different software; especially with AM parts, the integration between CAD software, where the part is modeled, computer-aided engineering (CAE) software, where the part is optimized and validated, and computer-aided manufacturing (CAM) software, where the part is prepared for the production phase, is fundamental, and a suitable file format is required to describe the parts and to be able to interact and go "back and forth" at every time [31,32]. To avoid these limitations, the design of lattice structure can be performed via wireframe modeling, a suitable method for the design and optimization of lattice structure by describing the geometry using mono-dimensional elements [33]. The lattice structure, constructed as an array of unit cells, can be morphed along a direction or surfaces thanks to advanced CAD tools [32,34,35]. An alternative method for conforming a lattice inside a design space, already used in virtual reality and biomedical fields, is the NURBS free-form deformation (FFD) technique, able to morph any kind of geometry into a design domain [14]. It is based on the creation of a trivariate tensor volume around the object; the volume can be deformed through its vertices to deform the object. Several evolutions of the FFD approach were developed from the first proposed by Sederberg and Parry [36], such as improvements in the flexibility of the control volume, combining several tricubic Bezier volumes [37], and enhancing the control of deformations through vertex weighting via NURBS-based approach [38]. Some software [39,40] already implemented this method through the creation of a trivariate NURBS cage around the object and the translation of its control points.
Furthermore, performing FEAs on lattice structures is not trivial due to the complex polyhedral mesh creation and the high number of nodes and finite elements that have to be processed. Again, high computational resources and time are needed. At the same time, FEAs need to be performed both during the size optimization phase and during the final validation of the model and cannot be neglected. Due to these reasons, alternative structural analysis methods, such as the asymptotic homogenization (AH) and mono-dimensional element methods, were developed [41][42][43][44]. In the AH approach, the lattice structure is studied as a fully dense material that has equivalent mechanical properties; this allows reducing the complexity of the model, but as a drawback, AH precludes working with anisotropic structures, so it is not possible to analyze pseudo-random and random lattice structures. As an alternative, the mono-dimensional finite element methods (FEM) based on the beam theory can be used. Two types of beams theory are widely used in the commercial FEM software: Euler-Bernoulli and Timoshenko. The Euler-Bernoulli approach assumes that the shear deformation is negligible compared to the bending deformation, so the beam cross-section remains perpendicular to its neutral axis. The critical parameter of this theory is the slenderness ratio, defined as the ratio between the length and the size of the side of the cross-section of the beam that affects the flexural rigidity matrix. When this value is low, in bending conditions, the flexural rigidity increases exponentially, causing underestimation of the deflection result. Due to the typical small dimensions of the elements of the lattice structures in AM, this ratio becomes an applicability limit [45]. The Timoshenko theory includes the deformation contribution coming from the shear stress, allowing the rotation between the cross-section and the neutral axis of the bending curve; therefore, it is more accurate. The limitation of this theory is based on the beam thickness parameter; if it is too small, a phenomenon called shear locking can occur. The shear locking problem lies in a bad schematic of the shear part of the elastic energy. To avoid the issue, a convergence Appl. Sci. 2021, 11, 11449 4 of 19 study needs to be performed subdividing the beam into n elements until the results are comparable to the exact theoretical values of displacement [46].
Bacciaglia et al. [47] reconstructed a lattice structure wireframe from a B-Rep using a novel voxel approach. They then created an ad-hoc MATLAB algorithm to transfer the mono-dimensional geometry in the FEM environment. Finally, they compared the result of the Euler-Bernoulli beam structural analysis approach with the results of the AH and 3D mesh FEM analysis approaches. The findings highlight the beam structural analysis has the best time-saving approach, but with uncertainty in the displacement estimation of 15%. Ravari et al. [48] studied the effects of the variation of the beam diameter on the mechanical proprieties of a lattice structure produced by fused filament fabrication. They compared the result of the mechanical tests with the Timoshenko beam and the solid finite element analysis showing a suitable agreement between the stress and Young's modulus results of the solid model and the experimental tests. The mono-dimensional beam theory showed a better precision in the prediction of the elongation and was more computationally efficient. Campoli et al. [49] carried out an FEA to predict the mechanical proprieties of porous titanium produced using selective laser melting (SLM), taking into account irregularities caused by the manufacturing process, comparing it to the analytical model and the experimental result. They found that the beam approach can predict the proprieties of pours biomaterial manufactured by SLM more accurately than the analytical. A similar conclusion was found by Genovese et al. [50], who studied the compression proprieties of the body-centered cubic (BCC) cell and the diamond (DA) cell, produced by SLM technologies, adopting digital image correlation (DIC) measurements during the experiments. Their numerical results using the beam theory on the BBC framework are more accurate than the ones from the solid model analysis. They hypothesized that the reason is the additional contribution of the stiffness provided by the joints of the struts. The results of the DA are not accurate in predicting the mechanical response of the real structure. Dong et al. [51] proposed to modify the Timoshenko beam element to consider the influence of the joint on the stiffness properties of the lattice structures due to the concentration of the material in the nodal points of the lattices. The validation of their joint stiffening element, through experimental mechanical tests and numerical analyses, demonstrated better accuracy than the traditional beam elements; more, a suitable agreement with experimental results and FE analyses with tetrahedral elements was found, with a considerable reduction in the computational costs. Zhu et al. [12] used the Timoshenko beam theory to determine the effect of the cell irregularity on the mechanical properties of 3D random open-cell foams. Each strut was divided into five elements according to the length of the strut. They found that highly structural irregularity brings to higher elastic and shear modulus but smaller bulk modulus concerning a perfect foam. They also found that their analytical model over-estimated the elastic modulus of a perfect foam with a large relative density. They suggested using a slenderness ratio greater than 4 to use the beam theory correctly. Evangelos et al. [44] studied the elastic modulus and Poisson ratio of a BCC cellular structure by analytical, numerical, and AH methods. Their experimental results are in excellent correlation with the numerical results derived by Timoshenko and AH methods. They also found a deviation of the results when adopting the Euler theory with respect to other numerical approaches when the slenderness ratio is higher than 10.
Even though the research community has proposed several solutions for reducing the simulation times and for simplifying the modeling of uniform lattice structures, a gap found in the literature regards the conformal lattice structures; indeed, the lack of a methodological approach for design, optimization, and geometric modeling inside the same environment was highlighted. Moreover, a comparison among the behavior of different cells in conformal lattices needs to be investigated in more detail [15,32,34,52,53].
Based on the mentioned issues, the aim of this paper is manifold. First, a novel method for designing conformal lattice structures based on the NURBS-FFD approach and for performing a size optimization in the same environment is proposed. The method is im-plemented using Python programming language and Karamba3D solver in the Rhinoceros environment. The integrated approach enables the exploitation of the previously developed geometric modeling methods based on subdivision surfaces [30]. Then, the monodimensional beam algorithm used for the size optimization of the conformal lattice is tested on four cells (simple cubic, body-centered cubic, modified Wallach-Gibson, and diamond), considering the technological limits of the additive manufacturing technologies such as minimum beam diameter and escape holes. Finally, the method was validated by performing numerical analyses based on the mono-dimensional beam theory inside two commercial FEM software, namely ANSYS Mechanical and MSC Nastran.

Materials and Methods
The methodology adopted for designing and validating optimized conformal lattice structures is explained in the following sub-sections. In the first subsection, the modeling methodology based on NURBS free-form deformation is described; in the second subsection, the size optimization method is briefly explained; in the third subsection, the validation methods through the beam linear structural analysis are described. The case study is a piston rod previously investigated by Rosso et al. [54,55]. Figure 1 shows the main dimensions of the test case, highlighting in green the design space that was replaced by lattices. is implemented using Python programming language and Karamba3D solver in the Rhinoceros environment. The integrated approach enables the exploitation of the previously developed geometric modeling methods based on subdivision surfaces [30]. Then, the mono-dimensional beam algorithm used for the size optimization of the conformal lattice is tested on four cells (simple cubic, body-centered cubic, modified Wallach-Gibson, and diamond), considering the technological limits of the additive manufacturing technologies such as minimum beam diameter and escape holes. Finally, the method was validated by performing numerical analyses based on the mono-dimensional beam theory inside two commercial FEM software, namely ANSYS Mechanical and MSC Nastran.

Materials and Methods
The methodology adopted for designing and validating optimized conformal lattice structures is explained in the following sub-sections. In the first subsection, the modeling methodology based on NURBS free-form deformation is described; in the second subsection, the size optimization method is briefly explained; in the third subsection, the validation methods through the beam linear structural analysis are described. The case study is a piston rod previously investigated by Rosso et al. [54,55]. Figure 1 shows the main dimensions of the test case, highlighting in green the design space that was replaced by lattices.

NURBS Free-Form Deformation Approach for Conformal Wireframe
The modeling of the conformal wireframe is based on the morphing of a lattice structure inside a design space and is performed in Rhinoceros 3D CAD software (release 7.8.2) and the relative visual programming environment Grasshopper.
First, a regular lattice structure is created by repeating the unit cell with a linear array along with the X-, Y-, and Z-axes. Four different beam-based unit cells are used: simple cubic (SC), body-centered cubic (BCC), modified Wallach-Gibson (WG), and diamond (DA) unit cell [56] (Figure 2). These unit cells are selected to show different behavior of bending and stretching-dominated lattices, according to the Maxwell rule [8], and to consider symmetry and the number of elements. The SC and the DA cells have a bendingdominated behavior, whereas the BCC and the WG cells a stretching-dominated one. SC and BCC belong to the same cubic cell family, sharing the beams at the edges of the cube; the main difference is based on the additional beams of the BCC, from the vertices to the center, which reinforce the cell. As opposed to the SC and BCC cells, the WG and DA ones

NURBS Free-Form Deformation Approach for Conformal Wireframe
The modeling of the conformal wireframe is based on the morphing of a lattice structure inside a design space and is performed in Rhinoceros 3D CAD software (release 7.8.2) and the relative visual programming environment Grasshopper.
First, a regular lattice structure is created by repeating the unit cell with a linear array along with the X-, Y-, and Z-axes. Four different beam-based unit cells are used: simple cubic (SC), body-centered cubic (BCC), modified Wallach-Gibson (WG), and diamond (DA) unit cell [56] (Figure 2). These unit cells are selected to show different behavior of bending and stretching-dominated lattices, according to the Maxwell rule [8], and to consider symmetry and the number of elements. The SC and the DA cells have a bendingdominated behavior, whereas the BCC and the WG cells a stretching-dominated one. SC and BCC belong to the same cubic cell family, sharing the beams at the edges of the cube; the main difference is based on the additional beams of the BCC, from the vertices to the center, which reinforce the cell. As opposed to the SC and BCC cells, the WG and DA ones are asymmetric, so the cell orientation affects the results related to the load condition. Finally, BCC and WG have more elements in the unit cell. are asymmetric, so the cell orientation affects the results related to the load condition. Finally, BCC and WG have more elements in the unit cell. Then, each array is morphed by the NURBS-FFD method to fit inside the design space of the piston road, while the straight elements of the cells are transformed into curve elements. More in detail, the NURBS-FFD technique is available in several commercial CAD software such as Rhinoceros, Blender, Autodesk Maya, and Autodesk 3ds Max and provides a user interface that simplifies the transformation of the regular lattice into a conformal one. Indeed, it is sufficient to move the vertices of a trivariate NURBS to adapt a box to the design space, simultaneously transforming the rods contained within it.
The deformation aims to maintain the cell aspect ratio similar to the undeformed cells in the central zone of the design space and, at the same time, cells that conform to the borders at the boundaries. In detail, this is performed by creating a trivariate NURBS volume around the object and moving its control points. The detailed procedure is as follows:  A regular wireframe based on the repetition along the X-, Y-, and Z-axis is generated. The number of instances is 14 cells along the Z-direction, 10 cells along the X-direction, and 4 cells along Y-direction for the BC, BCC, and WG. However, for the DA, the instances along the X-direction are 6 to maintain the aspect ratio ( Figure 3a);  A 3D NURBS volume is built around the wireframe using the Cage-Edit command.
In the case study, the cage is governed by 5 control points in X-direction, 2 control points in Y-direction, and 3 control points in Z-direction ( Figure 3b).
The control cage is adapted to the vertices of the design space and the curves of the boundaries by moving the control points ( Figure 3c).
This approach can be automated by defining a specific transformation of the cage, and it adapts to different geometric configurations of the design space. Moreover, it allows warping transformations if densification of the cells is needed.

Structural Size Optimization
This work refers to the method firstly implemented by Savio et al. [56] using a Phyton script and exploiting the Karmaba3D libraries [57], a Grasshopper plug-in for performing FEAs, and it is based on the utilization factor described in the Eurocode 3 [58]. The utilization factor considers the axial force, biaxial bending, torsion, shear force, and buckling load in compression. As an input, the method requires the wireframe model, the material properties, the boundary conditions, and the size constraints. Finally, the iterative FEA is run, varying the beam size to obtain the desired utilization for each element inside a predefined range [54]. Then, each array is morphed by the NURBS-FFD method to fit inside the design space of the piston road, while the straight elements of the cells are transformed into curve elements. More in detail, the NURBS-FFD technique is available in several commercial CAD software such as Rhinoceros, Blender, Autodesk Maya, and Autodesk 3ds Max and provides a user interface that simplifies the transformation of the regular lattice into a conformal one. Indeed, it is sufficient to move the vertices of a trivariate NURBS to adapt a box to the design space, simultaneously transforming the rods contained within it.
The deformation aims to maintain the cell aspect ratio similar to the undeformed cells in the central zone of the design space and, at the same time, cells that conform to the borders at the boundaries. In detail, this is performed by creating a trivariate NURBS volume around the object and moving its control points. The detailed procedure is as follows: • A regular wireframe based on the repetition along the X-, Y-, and Z-axis is generated. The number of instances is 14 cells along the Z-direction, 10 cells along the X-direction, and 4 cells along Y-direction for the BC, BCC, and WG. However, for the DA, the instances along the X-direction are 6 to maintain the aspect ratio ( Figure 3a); • A 3D NURBS volume is built around the wireframe using the Cage-Edit command. In the case study, the cage is governed by 5 control points in X-direction, 2 control points in Y-direction, and 3 control points in Z-direction ( Figure 3b).  Two different load cases are applied, a traction load case of 7.5 kN for comparing the result with the previous work of Rosso et al. [54] and a compression load case of −10 kN for evaluating the structure's stability and buckling effects. The load is evenly applied to the nodes of the beams connected to the big rod's end (Table 1). All the displacements and rotations of the beams nodes connected to the small rod's end are locked; the displacements of the beams nodes connected to the big rod's end are locked along the X-and Ydirection, while they can move along the Z-direction and rotate around all the axes ( Figure  4). Even if the boundary conditions do not perfectly represent the actual conditions, they The control cage is adapted to the vertices of the design space and the curves of the boundaries by moving the control points (Figure 3c). This approach can be automated by defining a specific transformation of the cage, and it adapts to different geometric configurations of the design space. Moreover, it allows warping transformations if densification of the cells is needed.

Structural Size Optimization
This work refers to the method firstly implemented by Savio et al. [56] using a Phyton script and exploiting the Karmaba3D libraries [57], a Grasshopper plug-in for performing FEAs, and it is based on the utilization factor described in the Eurocode 3 [58]. The utilization factor considers the axial force, biaxial bending, torsion, shear force, and buckling load in compression. As an input, the method requires the wireframe model, the material properties, the boundary conditions, and the size constraints. Finally, the iterative FEA is run, varying the beam size to obtain the desired utilization for each element inside a predefined range [54].
Two different load cases are applied, a traction load case of 7.5 kN for comparing the result with the previous work of Rosso et al. [54] and a compression load case of −10 kN for evaluating the structure's stability and buckling effects. The load is evenly applied to the nodes of the beams connected to the big rod's end ( Table 1). All the displacements and rotations of the beams nodes connected to the small rod's end are locked; the displacements of the beams nodes connected to the big rod's end are locked along the X-and Y-direction, while they can move along the Z-direction and rotate around all the axes ( Figure 4). Even if the boundary conditions do not perfectly represent the actual conditions, they may be sufficient to test and validate the method on different configurations. Due to the impossibility of producing cellular structures with traditional manufacturing technologies, the piston rod is supposed to be produced by AM; in particular, the SLM process, which belongs to the powder bed fusion category [59], is chosen due to its consolidated capability in producing almost fully dense lattices and, more in general, components characterized by mechanical properties similar to the parts obtained by conventional technologies [60,61], and due to its diffusion in the industry [62]. The selected AM machine is the Renishaw AM 400, and the material is the aluminum AlSi10Mg-0403 powder produced by Renishaw [63]. The mechanical proprieties declared by the producer are summarized in Table 2. The material was assumed as a linear elastic isotropic material without considering the influence of the printing direction, and the deformations are assumed to be small. The utilization factor for the optimization is set to (90 ± 1)% of the yield strength. The beam crosssection is circular. The size constraints were defined to avoid interference between the beams, to avoid powder evacuation problems, and considering the manufacturing size constraints given by the SLM technology, especially the minimum manufacturable diameter of beams [64]: the allowed beam diameter ranges from 0.5 to 1.5 mm for the SC, BCC, and WG structures, whereas, due to the bigger cell dimension of the DA lattice, the diameter ranges from 0.95 to 4.75 mm. The FE models set up in the three software are presented in Figure 4. The structural analysis of the case study was simplified, and only the results regarding the displacement and direct stress were investigated. Direct stress is defined as the stress component due to the axial load encountered in a beam element, and the axial load is the force acting on an object parallel to its axis. For a better understanding of the behavior of the piston rod, a detailed investigation should be performed considering, for instance, failure criteria such as the Von Mises criterion [73,74].

Results and Discussions
The developed method allows designing conformal lattice structures based on the NURBS-FFD approach and performing size optimization. Due to the integration in the same environment, other kinds of lattices can be further explored, such as strut-TPMS

Beam Approach for Linear Structural Analysis
In this section, the implemented beam structural linear analysis approaches are presented. The FEM analyses were performed to validate the optimized lattice structure models.
Three displacement-based FEM software were employed: 3) plug-in with Rhinoceros 7 (7.8.2) [54]. Karamba3D is a parametric structural engineering tool that performs structural analyses of spatial trusses, frames, and shells, and it is fully embedded in the parametric environment of Grasshopper, the visual programming language inside Rhinoceros. This makes it easy to combine parameterized geometric models, finite element analyses, and optimization algorithms. Karamba3D has been mostly used in the early stage design [65]; • MSC Patran/Nastran [47]. Patran is a pre-/post-processing software for FEA that provides modeling, meshing, analysis setup, and post processing for solvers. Nastran is a multidisciplinary structural solver for accurate static, dynamic, and thermal analyses; • Ansys 2020 R1 Mechanical [44]. Ansys is an FEA software tool for multidisciplinary simulation with a common platform, ANSYS Workbench, that can connect the tools together. As for Nastran, ANSYS Mechanical is the dedicated tool for multidisciplinary structural analyses.
In Karamba3D, the structural beam analysis is based on the Timoshenko beam theory with a linear shape function [65]. In Nastran, the beam model is based on the Euler-Bernoulli beam element. In ANSYS Mechanical, two types of the beam element can be selected: beam 188 (or beam 189) based on the Timoshenko theory and beam 3 based on the Euler-Bernoulli theory; moreover, through the Mechanical APDL interface, it is possible to select the type of the shape function (linear, quadratic, cubic) to define more integration points along with the element for enhancing the results and visualization accuracy; in this work, the beam 188 element with a quadratic shape function was selected.
In the Euler-Bernoulli beam theory, it is assumed that the plane cross-sections remain plain and perpendicular to the axis of the beam after the deformation; this assumption involves that all the transverse shear strains are null. The Timoshenko beam theory is an extension of the Euler-Bernoulli one: the plane of the sections remains plain but not necessarily perpendicular to the longitudinal axis after the deformation, thus admitting a nonzero transverse shear strain [66]. In operational terms, for a static case, Timoshenko's beam model, if compared to Euler-Bernoulli's model, has an additional energetic contribution given by the shear. A number of refined theories have been proposed for the analysis of shear deformable beams in the last four decades [67][68][69][70][71][72], but they are not adopted by the used commercial FEM software and are beyond the scope of this work. Regardless of the software, a simplification of the geometry is mandatory since the conformal cells are described by curves, whereas the application of the beam theories requires lines. To do this, the start and the end of each transformed element are connected by a line.
The beam structural validation using Karamba3D is carried out in the same environment of the NURBS-FFD modeling approach and the size optimization, so the information about the structure, material, and boundary conditions are directly available and easy to modify. Patran/Nastran and Ansys software require instead to import the model information through text files and to develop ad-hoc scripts using the dedicated programming languages of the FEM software. The lattice structure wireframe information needed to reconstruct the model are the following:

•
The node coordinates list contains the cartesian coordinates of all the nodes of the structure; • The beam connection list contains the indexes of the starting and ending node of each beam of the structure; • The beam size list contains the diameter of the cross-section of each beam of the structure.  Table 2).
The FE models set up in the three software are presented in Figure 4. The structural analysis of the case study was simplified, and only the results regarding the displacement and direct stress were investigated. Direct stress is defined as the stress component due to the axial load encountered in a beam element, and the axial load is the force acting on an object parallel to its axis. For a better understanding of the behavior of the piston rod, a detailed investigation should be performed considering, for instance, failure criteria such as the Von Mises criterion [73,74].

Results and Discussions
The developed method allows designing conformal lattice structures based on the NURBS-FFD approach and performing size optimization. Due to the integration in the same environment, other kinds of lattices can be further explored, such as strut-TPMS (Triply periodic minimal surfaces) [75], easily evaluating other projects and boundary conditions. The detailed findings are presented below.

Conformal Wireframe Modeling
The front and side views of the wireframes obtained by the NURBS-FFD method are presented in Figure 5. The method allows obtaining wireframe lattice structures that conform to the boundaries in the same environment but, as a drawback, it can be complicated to reach the aimed deformation of the captive parts depending on the number of control points; more, in Rhinoceros software, the NURBS-FFD control volume is rectangular or cubical leading to a trial-and-error approach when moving the points to follow the external shape of the object.
The front and side views of the wireframes obtained by the NURBS-FFD method are presented in Figure 5. The method allows obtaining wireframe lattice structures that conform to the boundaries in the same environment but, as a drawback, it can be complicated to reach the aimed deformation of the captive parts depending on the number of control points; more, in Rhinoceros software, the NURBS-FFD control volume is rectangular or cubical leading to a trial-and-error approach when moving the points to follow the external shape of the object. In all the configurations, the cell aspect ratio remains similar to the undeformed structure in the central zone, whereas the cells are more stretched along the design space borders to conform to the boundary. The types of the structure resulting from the application of Maxwell's rule to the lattices are reported in Table 3, together with the number of nodes and beams and the total length of the beams of the models. The SC and DA structures are bending-dominated, whereas the BCC and WG structures are stretching-dominated. Even if the SC is a bending-dominated structure, considering the load direction, it behaves as a stretching-dominated, as explained below. The DA structure presents the lowest number of elements due to the bigger size of the unit cell. The WG unit cell is the most populated, so the lattice structure presents the highest number of elements. More, being the BCC an evolution of the SC unit cell, the number of elements increases as expected.  In all the configurations, the cell aspect ratio remains similar to the undeformed structure in the central zone, whereas the cells are more stretched along the design space borders to conform to the boundary. The types of the structure resulting from the application of Maxwell's rule to the lattices are reported in Table 3, together with the number of nodes and beams and the total length of the beams of the models. The SC and DA structures are bending-dominated, whereas the BCC and WG structures are stretching-dominated. Even if the SC is a bending-dominated structure, considering the load direction, it behaves as a stretching-dominated, as explained below. The DA structure presents the lowest number of elements due to the bigger size of the unit cell. The WG unit cell is the most populated, so the lattice structure presents the highest number of elements. More, being the BCC an evolution of the SC unit cell, the number of elements increases as expected.  Figure 6 shows the size-optimized lattice structures. In both the traction and compression load cases, the optimization algorithm converged after the first iteration for the SC, BCC, and WG lattice structures, whereas it required 20 iterations for the DA structure.

Structural Size Optimization
based on Hooke's equation derived for the tensional stress of a uniform bar. The stretching-dominated structures size optimization process converges immediately, as the characteristic feature of these structures is the predisposition for axial loads. Consequently, the large number of iterations required for the DA is due to the algorithm that tends to increase the size of certain beams without increasing the diameters of the surrounding ones, leading to divergence, especially in the case of bending-dominated structures, in line with previous studies [56,76]. The bar chart in Figure 7 shows that the median value of the diameters of all the configurations equals the lower diameter constraint, i.e., 0.5 mm for the SC, BCC, and WG lattices and 0.95 mm for the DA lattice; this is explained by the presence of several oversized beams that, according to the size optimization, should have a diameter smaller than the minimum allowable size, but in the final configuration present the minimum crosssection to guarantee the printability of the part. The mean value of the diameter of the SC structure is higher than the BCC and WG ones due to the lower number of beams not aligned to the load direction. The DA lattice structure is the only one that requires beams with a cross-section as high as the upper manufacturing constraint, i.e., 4.75 mm. Figure 8 shows the values of the slenderness ratio and the relative density of the structures. The relative density is calculated through the Karamba3D plug-in, which does not consider the nodal connection between the beams, overestimating the value. The presence of a high number of oversized beams (beams with the minimum manufacturable diameter, i.e., 0.5 mm, even though the optimized size would be lower) affects the relative density of the structures with a stretching-dominated behavior, namely the SC, the BCC, and WG, which increases according to the total number of elements. Moreover, due to Maxwell's criterion, the BCC and WG structures are over-constrained, whereas the SC and DA structures are under-constrained; this reflects on the different ways the lattices react when changing the load case.
Moving from the traction to the compression load case, the relative density of the BCC and WG increases by about 13%, whereas the relative density of the SC and DA increases by about 24%. Even if the cell size of the DA structure is 67% bigger than the others, due to the number of cells in the X-direction, this structure presents the lowest minimum slenderness ratio, 0.85 and 0.67 for the traction and compression load case, due to The increase in the number of iterations given by the DA structure is due to the diameter variation mechanism of the algorithm, which calculates the optimized diameter based on Hooke's equation derived for the tensional stress of a uniform bar. The stretchingdominated structures size optimization process converges immediately, as the characteristic feature of these structures is the predisposition for axial loads. Consequently, the large number of iterations required for the DA is due to the algorithm that tends to increase the size of certain beams without increasing the diameters of the surrounding ones, leading to divergence, especially in the case of bending-dominated structures, in line with previous studies [56,76].
The bar chart in Figure 7 shows that the median value of the diameters of all the configurations equals the lower diameter constraint, i.e., 0.5 mm for the SC, BCC, and WG lattices and 0.95 mm for the DA lattice; this is explained by the presence of several oversized beams that, according to the size optimization, should have a diameter smaller than the minimum allowable size, but in the final configuration present the minimum cross-section to guarantee the printability of the part. The mean value of the diameter of the SC structure is higher than the BCC and WG ones due to the lower number of beams not aligned to the load direction. The DA lattice structure is the only one that requires beams with a cross-section as high as the upper manufacturing constraint, i.e., 4.75 mm.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 12 of 19 the bending-dominated behavior. As previously said, it is not advised to apply the beam theory when dealing with low values of slenderness ratio [77], and consequently, the results related to DA could be meaningless. This, together with the objective of maintaining the aspect ratio of the cell, led to keeping a low number of cells for the DA structure. On the opposite, since the overall maximum slenderness ratio is a low value of 13.16, the lattices are unaffected from the instability issues, having, therefore, a high Euler's critical buckling load [78].   Figure 8 shows the values of the slenderness ratio and the relative density of the structures. The relative density is calculated through the Karamba3D plug-in, which does not consider the nodal connection between the beams, overestimating the value. The presence of a high number of oversized beams (beams with the minimum manufacturable diameter, i.e., 0.5 mm, even though the optimized size would be lower) affects the relative density of the structures with a stretching-dominated behavior, namely the SC, the BCC, and WG, which increases according to the total number of elements. Moreover, due to Maxwell's criterion, the BCC and WG structures are over-constrained, whereas the SC and DA structures are under-constrained; this reflects on the different ways the lattices react when changing the load case.

Structural Analysis Validation
In this section, the structural analyses results are presented. Figure 9 shows the simulation environments of the three FEM software for the BBC lattice structure during the compression load case (−10 kN). As previously discussed, the boundary conditions consist Moving from the traction to the compression load case, the relative density of the BCC and WG increases by about 13%, whereas the relative density of the SC and DA increases by about 24%. Even if the cell size of the DA structure is 67% bigger than the others, due to the number of cells in the X-direction, this structure presents the lowest minimum slenderness ratio, 0.85 and 0.67 for the traction and compression load case, due to the bending-dominated behavior. As previously said, it is not advised to apply the beam theory when dealing with low values of slenderness ratio [77], and consequently, the results related to DA could be meaningless. This, together with the objective of maintaining the aspect ratio of the cell, led to keeping a low number of cells for the DA structure. On the opposite, since the overall maximum slenderness ratio is a low value of 13.16, the lattices are unaffected from the instability issues, having, therefore, a high Euler's critical buckling load [78].

Structural Analysis Validation
In this section, the structural analyses results are presented. Figure 9 shows the simulation environments of the three FEM software for the BBC lattice structure during the compression load case (−10 kN). As previously discussed, the boundary conditions consist in locking the displacements along the X-and Y-direction of the nodes of the beams connected to the big rod's end and locking the displacements and rotations of the nodes of the beams connected to the small rod's end. in locking the displacements along the X-and Y-direction of the nodes of the beams connected to the big rod's end and locking the displacements and rotations of the nodes of the beams connected to the small rod's end.  Figure 10 shows the maximum displacement along the Z-direction. The displacement results given by the three software are similar because, in the FE displacement-based type, the displacement is more accurate than the strain field [77]. A bending-dominated structure experiences the bending of its struts, whereas a stretching-dominated structure deforms primarily by the stretching or compressing of individual struts under loading. So, as can be seen from the results, the BC, BCC, and WG structures present a stretchingdominated behavior, whereas the DA has a bending-dominated behavior. According to the connectivity of the structure and Maxwell's law, the BC structure should be bendingdominated but appears to be the stiffest structure between the four. The BC structure shows a stretching-dominated behavior because it is mainly stressed along the main axes; different boundary conditions would lead to a different behavior [5]. The BCC and DA structures confirm Maxwell's rule prediction; indeed, even if they have a similar apparent density, the displacement of the DA structure is more than 10 times higher than the displacement of the BC one. Since the SC and the BCC structures come from the same unit cell, they exhibit similar mechanical behavior; the BCC is more constrained and therefore has a smaller deflection in both the load cases.   Figure 10 shows the maximum displacement along the Z-direction. The displacement results given by the three software are similar because, in the FE displacement-based type, the displacement is more accurate than the strain field [77]. A bending-dominated structure experiences the bending of its struts, whereas a stretching-dominated structure deforms primarily by the stretching or compressing of individual struts under loading. So, as can be seen from the results, the BC, BCC, and WG structures present a stretchingdominated behavior, whereas the DA has a bending-dominated behavior. According to the connectivity of the structure and Maxwell's law, the BC structure should be bendingdominated but appears to be the stiffest structure between the four. The BC structure shows a stretching-dominated behavior because it is mainly stressed along the main axes; different boundary conditions would lead to a different behavior [5]. The BCC and DA structures confirm Maxwell's rule prediction; indeed, even if they have a similar apparent density, the displacement of the DA structure is more than 10 times higher than the displacement of the BC one. Since the SC and the BCC structures come from the same unit cell, they exhibit similar mechanical behavior; the BCC is more constrained and therefore has a smaller deflection in both the load cases.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 13 of 19 in locking the displacements along the X-and Y-direction of the nodes of the beams connected to the big rod's end and locking the displacements and rotations of the nodes of the beams connected to the small rod's end.  Figure 10 shows the maximum displacement along the Z-direction. The displacement results given by the three software are similar because, in the FE displacement-based type, the displacement is more accurate than the strain field [77]. A bending-dominated structure experiences the bending of its struts, whereas a stretching-dominated structure deforms primarily by the stretching or compressing of individual struts under loading. So, as can be seen from the results, the BC, BCC, and WG structures present a stretchingdominated behavior, whereas the DA has a bending-dominated behavior. According to the connectivity of the structure and Maxwell's law, the BC structure should be bendingdominated but appears to be the stiffest structure between the four. The BC structure shows a stretching-dominated behavior because it is mainly stressed along the main axes; different boundary conditions would lead to a different behavior [5]. The BCC and DA structures confirm Maxwell's rule prediction; indeed, even if they have a similar apparent density, the displacement of the DA structure is more than 10 times higher than the displacement of the BC one. Since the SC and the BCC structures come from the same unit cell, they exhibit similar mechanical behavior; the BCC is more constrained and therefore has a smaller deflection in both the load cases.   Figure 11 represents the maximum direct stress of the structures. The values of the WG and DA lattice structures from the three FE software disagree, most likely due to the very small size of the elements for the WG structure and the very small slenderness ratio for the DA one. A divergence of the results due to the small ratios of slenderness was also observed by Ptochos and Ahmadi [44,79]. For thicker beams, the shear and rotational inertia effects, which are taken into account in the Timoshenko theory but are neglected in the Euler-Bernoulli theory, start to play an important role; the shear and rotational inertia effects cause the beam to be less stiff, explaining a higher elastic modulus calculated using the Euler-Bernoulli rather than the solution based on the Timoshenko beam theory [43,46,73].
Appl. Sci. 2021, 11, x FOR PEER REVIEW 14 of 19 Figure 11 represents the maximum direct stress of the structures. The values of the WG and DA lattice structures from the three FE software disagree, most likely due to the very small size of the elements for the WG structure and the very small slenderness ratio for the DA one. A divergence of the results due to the small ratios of slenderness was also observed by Ptochos and Ahmadi [44,79]. For thicker beams, the shear and rotational inertia effects, which are taken into account in the Timoshenko theory but are neglected in the Euler-Bernoulli theory, start to play an important role; the shear and rotational inertia effects cause the beam to be less stiff, explaining a higher elastic modulus calculated using the Euler-Bernoulli rather than the solution based on the Timoshenko beam theory [43,46,73]. The difference between the results of Ansys and Karamba3D could be explained by the shear locking phenomenon that occurs where the application of the linear shape function leads to an over-stiffened element. To avoid this issue, the element can be subdivided into more instances, or the shape function can be changed with a quadratic or cubic one [46]. The element used in Ansys has a quadratic shape function, with an additional interpolation point in the center of the beam with respect to the Karamba3D linear element, and this could explain the difference between the direct stress result.
Indeed, from the traction to the compression configuration, an increase in the diameter is observed with a consequent decrease in the minimum slenderness ratio, causing higher divergence between the FEM results. Therefore, a limitation of the application of the beam structural analysis inside FEM software is the assumption that the elements have a minimum slenderness ratio; the results of these analyses show that for slenderness ratios lower than 2.4, the solutions between the various structural analysis tools begin to diverge. Otherwise, mechanical validation methods for lattice structures with thick beams require experimental analyses or FEM simulations involving solid elements. As a drawback, those alternatives are expected to require high costs for manufacturing the specimens, difficulties to test components with complex shapes by standard testing protocols, high computational power, and difficulties in the meshing process. Furthermore, the homogenization method cannot be performed on conformal lattice structures since they are anisotropic and have different mechanical properties along with different directions. Moreover, the actual beam joints' stiffness due to the beams overlapping was not evaluated; this lack could lead to the underestimation of the stiffness of the model, giving an under-evaluation of the stiffness if compared to the experimental results and FE analysis based on solid elements [34,51,80]. Both the mono-dimensional beam approaches and the solid element approaches do not include manufacturing irregularities of the beam size, The difference between the results of Ansys and Karamba3D could be explained by the shear locking phenomenon that occurs where the application of the linear shape function leads to an over-stiffened element. To avoid this issue, the element can be subdivided into more instances, or the shape function can be changed with a quadratic or cubic one [46]. The element used in Ansys has a quadratic shape function, with an additional interpolation point in the center of the beam with respect to the Karamba3D linear element, and this could explain the difference between the direct stress result.
Indeed, from the traction to the compression configuration, an increase in the diameter is observed with a consequent decrease in the minimum slenderness ratio, causing higher divergence between the FEM results. Therefore, a limitation of the application of the beam structural analysis inside FEM software is the assumption that the elements have a minimum slenderness ratio; the results of these analyses show that for slenderness ratios lower than 2.4, the solutions between the various structural analysis tools begin to diverge. Otherwise, mechanical validation methods for lattice structures with thick beams require experimental analyses or FEM simulations involving solid elements. As a drawback, those alternatives are expected to require high costs for manufacturing the specimens, difficulties to test components with complex shapes by standard testing protocols, high computational power, and difficulties in the meshing process. Furthermore, the homogenization method cannot be performed on conformal lattice structures since they are anisotropic and have different mechanical properties along with different directions. Moreover, the actual beam joints' stiffness due to the beams overlapping was not evaluated; this lack could lead to the underestimation of the stiffness of the model, giving an under-evaluation of the stiffness if compared to the experimental results and FE analysis based on solid elements [34,51,80]. Both the mono-dimensional beam approaches and the solid element approaches do not include manufacturing irregularities of the beam size, which negatively affect the mechanical proprieties of the entire printed lattice structure [43,48].
To prevent manufacturing issues, the beam inclination angle relative to the printing direction, also referred to as the overhang angle, and the minimum and maximum diameter of the beams are factors that must be considered. The overhang angle is suggested to be greater than 30 • for avoiding loss of powder, which could lead to a partially melted zone compared to a vertical beam causing different thermal conditions that may create unique microstructures and highly irregular surface morphology [11,81]. To avoid the need for supports for the lattice structure, which would result difficult to be removed, the lattice structures could be angled for production, but it is not possible to define an ideal printing orientation due to the variable orientation of the conformal lattice. Nevertheless, recent studies have shown that small sizes of the beam may not require the use of supports [82,83], so considering the size of the beams of the studied lattice structures, it might be possible to produce them vertically by SLM technologies.

Conclusions
In this manuscript, the modeling, size optimization, and structural linear analysis of conformal lattice structures were presented. In particular, the work aimed to provide design tools that require a low process time and low computational power in a preliminary design phase, thanks to the implementation of the design, optimization, and geometric modeling in a single environment and to the 1D analysis that is faster than solid FEA. A geometric design approach based on the deformation of the wireframe of a regular lattice structure through the NURBS-FFD method was proposed. The possibility of performing size optimization and structural linear analysis based on the beam theory in the same environment was explored. An automated process was then developed for importing the geometrical information and constructing FE models inside the Patran/Nastran and Ansys commercial FEA software. This method promotes the design, optimization, and structural analysis of conformal cellular structures and facilitates further real-time geometrical modification through the control points of the NURBS-FFD control volume inside the same environment. At the same time, limitations are still present; the NURBS-FFD approach offers a fixed topology of the control volume, and the graphic interface can appear confusing when dealing with volumes described by many control points. Furthermore, the proposed size optimization approach and the validation methods are limited to beam-based lattice structures with a restriction on the minimum slenderness ratio.
As a case study, the central part of a piston rod was filled with four different lattice structures, based on the SC, BCC, WG, and DA unit cells, two of them with thick struts to highlight the limitation of the beam structural size optimization and structural linear analysis. It was shown how bending behavior in conformal lattice structures not only depends on the Maxwell criterion but also on the load direction with respect to the main axes of the cell, as was seen for the SC structure. It was highlighted that the constraints of the size optimization algorithm should include not only the size limitations of the manufacturing process but also the minimum slenderness ratio to use the beam theory correctly and the maximum slenderness ratio to prevent buckling effects. In conclusion, it is expected that for beam-based lattice structures, the beam theory tends to be more conservative than experimental results in terms of stiffness, and therefore, the solutions have an advantage in terms of safety, which can be useful in the AM field since irregularities and defects of the parts negatively affect the mechanical performance.
As future works, new case studies will be used to validate the NURBS-FFD method, and other unit cells will be introduced, such as strut-TPMS. Dealing with size optimization, the algorithm could be modified to interact with field approaches for structural optimization and other size optimization strategies. Concerning the modeling method, an enhanced version able to provide a lattice structure with smooth connections between the beams and between the lattice and the solid model is going to be developed. Finally, experimental tests for validating the numerical results are needed. Funding: This work was partially funded by the Regional Operational Programme F.S.E. 2014-2020 Veneto Region and the European Regional Development Fund POR 2014-2020, project code 2105-0036-1463-2019, "Design and manufacturing by additive technologies of innovative components using functionally graded materials", and a Ph.D. grant of Fondazione Cassa di Risparmio di Padova e Rovigo (CARIPARO).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The raw/processed data required to reproduce these findings can be partially shared upon request.