Investigating infill density and pattern effects in additive manufacturing by characterizing metamaterials along the strain-gradient theory

Abstract

Infill density used in additive manufacturing incorporates a structural response change in the structure. Infill pattern creates a microstructure that affects the mechanical performance as well. Whenever the length ratio of microstructure to geometry converges to one, metamaterials emerge and the strain-gradient theory is an adequate model to predict metamaterials response. All metamaterial parameters are determined by an asymptotic homogenization, and we investigate the effects of infill density and pattern on these parameters. In order to illuminate the role of infill characteristics on the strain-gradient parameters, an in-depth numerical investigation is presented for one, widely used case in three-dimensional (3D) printers, rectangular grid.

1. Introduction

The usage of additive manufacturing (AM) techniques in different industrial applications has recently gained remarkable momentum, trending in different fields like aerospace, automotive, medical, architecture, and construction. AM techniques fabricate three-dimensional (3D) designs prepared in a CAD software using a layer-wise approach, enabling to manufacture complex parts, which are challenging or even impossible to be realized by conventional manufacturing techniques such as CNC milling, turning, and casting. As many design constraints are eliminated, AM methods have revolutionized the perception of engineers and designers leveraging complexity, creativity, and multi-functionality. Therefore, with their inexorable progress and distinctive properties, AM techniques are currently preferred in many applications as primary fabrication technique [1–4] enabling to introduce new innovative design concepts [5–11].

One of the most significant advantages of AM is that lightweight and topologically optimized parts may be designed and fabricated per its layer-by-layer fabrication approach. In this way, less material is used while providing desired performance and a considerable reduction in waste. For this purpose, during the tool planning process, also referred to as “slicing,” the virtual model (i.e. CAD design) is divided into layers taking heed of infill pattern and infill ratio. Then, each layer is formed accordingly by following the generated tool path, and the part with desired internal structural pattern is fabricated. The infill pattern and infill ratio may vary depending on parts to be printed and expected mechanical performance. Therefore, as AM deals with internal structure without requiring any extra manufacturing steps or equipment, design and fabrication of mechanical metamaterials is emerging in basic research. Metamaterials exhibit tailored, “exotic” properties allowing to design multiscale approaches [12–17]. Therefore, in mechanical metamaterials such as auxetic [18,19] or pantographic [3] structures, the unique overall behavior is obtained by mechanical interaction of microconstituents composing the structure.

Ability of fabricating complex materials with new advanced manufacturing methods has also revealed the need for new mathematical models as the classical ones are not sufficient or practical [20–25]. Indeed, this is an ineluctable fact from a theoretical point of view in a more complex design and manufacturing environment. For this purpose, higher-gradient theories have been extensively developed and investigated in the literature. Briefly, higher-gradient theories assume that the deformation energy may also include higher gradients of displacement field as opposed to the classical theory [26]. Although such theoretical efforts have been made since 1950s [27–33], fabrication of metamaterials showing higher-gradient effects by AM has succeeded to ignite an applied research. Specifically, the recent studies on pantographic structures have epitomized the importance of higher-gradient theories, especially for numerical predictions of such a complex structure. In order to investigate the mechanical behavior of pantographic structures, various higher-gradient models have been proposed [34–38] and examined through experiments [39–44], considering different mechanical phenomena [45–50]. We refer interested readers to the previous works [51–54] on different aspects of higher-gradient modeling.

In this study, strain-gradient modeling of additively manufactured specimens is investigated considering different infill ratios and infill patterns. In order to elicit the strain-gradient model of an additively manufactured specimen, an asymptotic homogenization framework is utilized. The computational homogenization procedure, which is based on the balance of micro- and macro-energies, determines all the constitutive parameters emerging in the strain-gradient theory. Via the determined parameters, simulations based on the strain-gradient model are conducted to investigate the mechanical behavior of the specimen. In this study, we examine the effects of infill ratio and infill pattern on the strain-gradient model and its constitutive parameters. To this end, a specimen with grid infill pattern has been first investigated considering different infill ratios. Then, by keeping the infill ratio fixed, various infill patterns have been studied. The numerical study has been aimed to divulge how strain-gradient parameters and deformation energy related to the strain-gradient terms change due to infill ratio and infill pattern. As a result, insights have been unraveled for design and strain-gradient modeling of such materials.

The rest of the paper is as follows. In section 2, the asymptotic homogenization is presented identifying strain-gradient constitutive parameters in order to provide mathematical preliminaries for the subsequent numerical study. In section 3, the conducted numerical study is presented for the specimen under study. Then, the conclusions are drawn in section 4.

2. Asymptotic homogenization to identify strain-gradient parameters

We follow the asymptotic homogenization approach introduced [55], applied [56,57], and verified [58]. Continuum with a microstructure (microscale) is modeled as a metamaterial equivalent to a homogenized continuum (macroscale) with the smeared out microstructure. At the microscale, expressed in $\mathit{y}$ , displacement, ${\mathit{u}}^{\mathrm{m}}$ , strain, and stress may be different than its corresponding homogenized continuum at the macroscale, expressed in $\mathit{x}$ . Between $\mathit{x}$ and $\mathit{y}$ , macro-to-microscale transformation is handled by a homothetic ratio, $ϵ$ . Technically, no scale separation is used, i.e., the same coordinate system is incorporated leading to a different material equation at micro- and macro-scales. By knowing one material model, we determine the parameters of the other material model. Herein, the approach is based on the given microscale with the material properties. Metamaterials’ parameters at the macroscale are aimed for. At the microscale, we use linear elastic isotropic model with the known stiffness tensor, ${\mathit{C}}^{\mathrm{m}}$ , the deformation energy density, $e^{\mathrm{m}}$ , is modeled as a quadratic function:

$$e^{\mathrm{m}}=\frac{1}{2}{\epsilon }_{\mathrm{ij}}^{\mathrm{m}}{C}_{\mathrm{ijkl}}^{\mathrm{m}}{\epsilon }_{\mathrm{kl}}^{\mathrm{m}}.$$

(1)

We use standard continuum mechanics formulation and understand Einstein’s summation convention over repeated indices. Linearized displacement gradients, ${\epsilon }^{\mathrm{m}}$ , create an objective energy formulation, they are given as,

$${\epsilon }_{\mathrm{ij}}^{\mathrm{m}}=\frac{1}{2}(\frac{\partial u_i}{\partial y_j}+\frac{\partial u_j}{\partial y_i}).$$

(2)

By using Castigliano’s theorem:

$$\frac{\partial e^{\mathrm{m}}}{\partial {\epsilon }_{\mathrm{ij}}^{\mathrm{m}}}={\sigma }_{\mathrm{ij}}^{\mathrm{m}}.$$

(3)

We obtain stress and the latter is called Hooke’s law. Since the energy is quadratic, the material law is linear. We use small displacement approach for the homogenization process leading to metamaterials’ parameters. Once the parameters are found, at the macroscale, geometric nonlinearities may be involved using a nonlinear strain measure [59,60]. At the macroscale, we model using the strain-gradient elasticity with the same type of (quadratic) stored energy density

$$e^{\mathrm{M}}=\frac{1}{2}{\epsilon }_{\mathrm{ij}}^{\mathrm{M}}{C}_{\mathrm{ijkl}}^{\mathrm{M}}{\epsilon }_{\mathrm{kl}}^{\mathrm{M}}+{\epsilon }_{\mathrm{ij}}^{\mathrm{M}}{G}_{\mathrm{ijklm}}^{\mathrm{M}}{\epsilon }_{\mathrm{kl},m}^{\mathrm{M}}+\frac{1}{2}{\epsilon }_{\mathrm{ij},k}^{\mathrm{M}}{D}_{\mathrm{ijklmn}}^{\mathrm{M}}{\epsilon }_{\mathrm{lm},n}^{\mathrm{M}},$$

(4)

with the usual restrictions as objectivity and positive definiteness [61,62]. A comma notation denotes a partial derivative in $\mathit{x}$ , which denotes the homogenized continuum. The same strain measure results in

$${\epsilon }_{\mathrm{ij}}^{\mathrm{M}}=\frac{1}{2}(u_{i,j}^{\mathrm{M}}+u_{j,i}^{\mathrm{M}}).$$

(5)

These additional metamaterial tensors, ${\mathit{C}}^{\mathrm{M}}$ , ${\mathit{G}}^{\mathrm{M}}$ , and ${\mathit{D}}^{\mathrm{M}}$ , need to be determined. In the case of the repeating microstructure, we characterize only one unit cell or a representative volume element (RVE) built by unit cells. Using periodic boundary conditions on the RVE boundaries, we capture the complete response of the body. By varying additional parameters, it is possible to determine the correct variables that fulfill one axiom: the deformation energy within the RVE, $\mathrm{\Omega }$ , is identical at micro- and macro-scales

$${\int }_{\mathrm{\Omega }}{e}^{\mathrm{m}}\mathrm{d}V={\int }_{\mathrm{\Omega }}{e}^{\mathrm{M}}\mathrm{d}V.$$

(6)

Displacements or stress may differ, but the energy is the same. This axiom is the starting point for the asymptotic homogenization. For the complete derivation of the governing equations, we refer to Abali and Barchiesi [55], and summarize the algorithm herein briefly. The core idea is to expand displacement field up to the second order, and then find out additional governing equations by coefficients comparison in this expansion in ${ϵ}^0$ , ${ϵ}^1$ , and ${ϵ}^2$ . Three equations emerge and the first equation explains that the leading term with ${ϵ}^0$ is equivalent to the displacement at the macroscale. Second equation delivers the newly introduced field $\phi$ to be calculated by

$$\frac{\partial }{\partial y_j}\left(C_{\mathrm{ijkl}}^{\mathrm{m}}{L}_{\mathrm{abkl}}\right)=0,L_{\mathrm{abkl}}={\delta }_{\mathrm{ak}}{\delta }_{\mathrm{bl}}+\frac{\partial {\phi }_{\mathrm{abk}}}{\partial y_l}.$$

(7)

Once the field $\phi$ is calculated, we obtain ${\mathit{C}}^{\mathrm{M}}$ as follows:

$$C_{\mathrm{abcd}}^{\mathrm{M}}=\frac{1}{V}{\int }_{\mathrm{\Omega }}{C}_{\mathrm{ijkl}}^{\mathrm{m}}{L}_{\mathrm{abij}}{L}_{\mathrm{cdkl}}\mathrm{d}V.$$

(8)

Third equation delivers $\psi$ by calculating,

$$\begin{array}{c}\multicolumn{1}{c}{-\frac{{\rho }^{\mathrm{m}}}{{\rho }^{\mathrm{M}}}{C}_{\mathrm{icab}}^{\mathrm{M}}+C_{\mathrm{ickl}}^{\mathrm{m}}{L}_{\mathrm{abkl}}+\frac{\partial }{\partial y_j}\left(C_{\mathrm{ijkl}}^{\mathrm{m}}{N}_{\mathrm{abckl}}\right)=0,}\\ \multicolumn{1}{c}{N_{\mathrm{abckl}}={\phi }_{\mathrm{abk}}{\delta }_{\mathrm{cl}}+\frac{\partial {\psi }_{\mathrm{abck}}}{\partial y_l}.}\end{array}$$

(9)

By means of $\phi$ and $\psi$ , we obtain

$$\begin{array}{c}{M}_{\mathrm{abcij}}=y_c{L}_{\mathrm{abij}}+N_{\mathrm{abcij}},\hfill \end{array}$$

(10)

which is used to determine all strain-gradient parameters:

$$\begin{array}{cc}\multicolumn{1}{c}{G_{\mathrm{abcde}}^{\mathrm{M}}=}& \frac{ϵ}{V}{\int }_{\mathrm{\Omega }}{C}_{\mathrm{ijkl}}^{\mathrm{m}}{L}_{\mathrm{abij}}{M}_{\mathrm{cdekl}}\mathrm{d}V,\hfill \\ \multicolumn{1}{c}{D_{\mathrm{abcdef}}^{\mathrm{M}}=}& \frac{{ϵ}^2}{V}({\int }_{\mathrm{\Omega }}{C}_{\mathrm{ijkl}}^{\mathrm{m}}{M}_{\mathrm{abcij}}{M}_{\mathrm{defkl}}\mathrm{d}V-C_{\mathrm{abde}}^{\mathrm{M}}{\int }_{\mathrm{\Omega }}{y}_c{y}_f\mathrm{d}V).\hfill \end{array}$$

(11)

We observe that ${\mathit{C}}^{\mathrm{M}}$ , ${\mathit{G}}^{\mathrm{M}}$ , and ${\mathit{D}}^{\mathrm{M}}$ depend on ${ϵ}^0$ , ${ϵ}^1$ , and ${ϵ}^2$ , respectively. Hence, stiffness tensor is constant for different homothetic ratios. The term ${\mathit{G}}^{\mathrm{M}}$ vanishes in case of centro-symmetric microstructure, which is the case herein. However, we compute their values in order to make sure that the code is working as expected. Obviously, ${\mathit{D}}^{\mathrm{M}}$ depends on $ϵ$ quadratically, this dependence on geometric characteristic length is known as the size effect in the literature.

3. Numerical implementation and results

All aforementioned governing equations are solved numerically by means of the finite element method (FEM). We use open-source packages known as FEniCS [63,64] and solve the governing equations by generating a weak form for each of them by following the standard variational formulation. We refer to Abali [65] for engineering applications and their implementation by means of FEniCS packages. Standard FEM properties are used, where the field functions are approximated by their discrete counterparts with a compact support. Technically, we span a finite dimensional Hilbertian Sobolev space for trial functions [66]. As known as the Galerkin approach, the same space is used for test functions as well. We emphasize that micro- and macro-scales are both expressed in the same coordinate system, which is an important benefit in the computational implementation. The computational domain is the RVE that is discretized by nodes and their connectivity by elements. This discretization procedure is handled by NetGen algorithms in Salome CAD software.

In two-dimensional (2D) case, for $\phi$ from equation (7), we obtain three different weak forms and for $\psi$ from equation (9), we obtain six different weak forms, by separately setting $a$ , $b$ , and $c$ indices (no sum), we may write them all together,

$$\begin{array}{c}{\int }_{\mathrm{\Omega }}{C}_{\mathrm{ijkl}}^{\mathrm{m}}{L}_{\mathrm{abkl}}\frac{\partial \delta {\phi }_{\mathrm{abi}}}{\partial y_j}\mathrm{d}V=0,\hfill \\ {\int }_{\mathrm{\Omega }}(-C_{\mathrm{icab}}^{\mathrm{M}}\delta {\psi }_{\mathrm{abci}}+C_{\mathrm{ickl}}^{\mathrm{m}}{L}_{\mathrm{abkl}}\delta {\psi }_{\mathrm{abci}}-\left(C_{\mathrm{ijkl}}^{\mathrm{m}}{N}_{\mathrm{abckl}}\right)\frac{\partial \delta {\psi }_{\mathrm{abci}}}{\partial y_j})\mathrm{d}V=0.\hfill \end{array}$$

(12)

Trial functions, $\mathbf{\phi }$ and $\mathbf{\psi }$ , as well as their test functions, $\delta \mathbf{\phi }$ and $\delta \mathbf{\psi }$ , are all approximated using standard Lagrange elements. We use continuous and linear form functions. Periodic boundary conditions are applied on all boundaries by mapping the solution at the corresponding surfaces. Simply stated, the solution is restricted in such a way that values of $\mathbf{\phi }$ and $\mathbf{\psi }$ are equal on corresponding surfaces. Therefore, the node positions have to match, which is established using projection method in meshing in Salome. Since the weak form is linear, we solve them using a direct method, specifically, using Mumps solver from PETSc package [67].

A specimen of length $L=110\mathrm{mm}$ and height $H=22\mathrm{mm}$ as shown in Figure 1 has been selected to investigate in this study, where apparently the microstructure is of the same order of magnitude. Additional support has been added to the left and right sides of the grid structure with the purpose of applying boundary conditions to the grid structure equally distributed.

Figure 1.

Geometric properties of specimen.

The specimen was first assumed to be designed with grid infill pattern and 40% infill ratio. RVE of the specimen is given in Figure 2, and the asymptotic homogenization procedure is applied to this RVE to extract strain-gradient constitutive parameters.

Figure 2.

RVE of initial design: grid infill pattern and 40% infill ratio.

We consider this specimen as built by PLA, which is one often used filament type in fused deposition modeling (FDM)-based 3D printers. We utilize Young’s modulus of PLA as $E=3500\mathrm{MPa}$ and Poisson’s ratio as $\nu =0.3$ . For the “voids,” we use Young’s modulus of $0.01\mathrm{MPa}$ .

Asymptotic homogenization results in ${\mathit{C}}^{\mathrm{M}}$ , ${\mathit{G}}^{\mathrm{M}}$ , and ${\mathit{D}}^{\mathrm{M}}$ values, representing the same deformation energy for the metamaterial as a homogenized continuum as shown in Figure 3. Therefore, strain-gradient simulations are conducted on the homogenized continuum to predict the mechanical behavior of specimen under study.

Figure 3.

Homogenized continuum to conduct strain-gradient simulations.

In order to present the parameters, we use Voigt-like notation. Concretely, we use $A$ , $B$ indicating $\{11,22,12\}$ and $\alpha$ , $\beta$ denoting $\{111,221,121,112,222,122\}$ , such that we obtain

$$C_{\mathrm{AB}}^{\mathrm{M}}=\left(\begin{array}{ccc}{C}_{1111}^{\mathrm{M}}\hfill & C_{1122}^{\mathrm{M}}\hfill & C_{1112}^{\mathrm{M}}\hfill \\ \multicolumn{1}{c}{C_{2211}^{\mathrm{M}}}& C_{2222}^{\mathrm{M}}\hfill & C_{2212}^{\mathrm{M}}\hfill \\ \multicolumn{1}{c}{C_{1211}^{\mathrm{M}}}& C_{1222}^{\mathrm{M}}\hfill & C_{1212}^{\mathrm{M}}\hfill \end{array}\right),$$

(13)

and

$$G_{A\alpha }^{\mathrm{M}}=\left(\begin{array}{cccccc}{G}_{11111}^{\mathrm{M}}\hfill & G_{11221}^{\mathrm{M}}\hfill & G_{11121}^{\mathrm{M}}\hfill & G_{11112}^{\mathrm{M}}\hfill & G_{11222}^{\mathrm{M}}\hfill & G_{11122}^{\mathrm{M}}\hfill \\ \multicolumn{1}{c}{G_{22111}^{\mathrm{M}}}& G_{22221}^{\mathrm{M}}\hfill & G_{22121}^{\mathrm{M}}\hfill & G_{22112}^{\mathrm{M}}\hfill & G_{22222}^{\mathrm{M}}\hfill & G_{22122}^{\mathrm{M}}\hfill \\ \multicolumn{1}{c}{G_{12111}^{\mathrm{M}}}& G_{12221}^{\mathrm{M}}\hfill & G_{12121}^{\mathrm{M}}\hfill & G_{12112}^{\mathrm{M}}\hfill & G_{12222}^{\mathrm{M}}\hfill & G_{12122}^{\mathrm{M}}\hfill \end{array}\right),$$

(14)

as well as

$$D_{\alpha \beta }^{\mathrm{M}}=\left(\begin{array}{cccccc}{D}_{111111}^{\mathrm{M}}\hfill & D_{111221}^{\mathrm{M}}\hfill & D_{111121}^{\mathrm{M}}\hfill & D_{111112}^{\mathrm{M}}\hfill & D_{111222}^{\mathrm{M}}\hfill & D_{111122}^{\mathrm{M}}\hfill \\ \multicolumn{1}{c}{D_{221111}^{\mathrm{M}}}& D_{221221}^{\mathrm{M}}\hfill & D_{221121}^{\mathrm{M}}\hfill & D_{221112}^{\mathrm{M}}\hfill & D_{221222}^{\mathrm{M}}\hfill & D_{221122}^{\mathrm{M}}\hfill \\ \multicolumn{1}{c}{D_{121111}^{\mathrm{M}}}& D_{121221}^{\mathrm{M}}\hfill & D_{121121}^{\mathrm{M}}\hfill & D_{121112}^{\mathrm{M}}\hfill & D_{121222}^{\mathrm{M}}\hfill & D_{121122}^{\mathrm{M}}\hfill \\ \multicolumn{1}{c}{D_{112111}^{\mathrm{M}}}& D_{112221}^{\mathrm{M}}\hfill & D_{112121}^{\mathrm{M}}\hfill & D_{112112}^{\mathrm{M}}\hfill & D_{112222}^{\mathrm{M}}\hfill & D_{112122}^{\mathrm{M}}\hfill \\ \multicolumn{1}{c}{D_{222111}^{\mathrm{M}}}& D_{222221}^{\mathrm{M}}\hfill & D_{222121}^{\mathrm{M}}\hfill & D_{222112}^{\mathrm{M}}\hfill & D_{222222}^{\mathrm{M}}\hfill & D_{222122}^{\mathrm{M}}\hfill \\ \multicolumn{1}{c}{D_{122111}^{\mathrm{M}}}& D_{122221}^{\mathrm{M}}\hfill & D_{122121}^{\mathrm{M}}\hfill & D_{122112}^{\mathrm{M}}\hfill & D_{122222}^{\mathrm{M}}\hfill & D_{122122}^{\mathrm{M}}\hfill \end{array}\right),$$

(15)

For the specific case of 40% infill ratio, parameters read as follows:

$$\begin{array}{c}{C}_{\mathrm{AB}}^{\mathrm{M}}=\left(\begin{array}{ccc}910\hfill & 105\hfill & 0\hfill \\ \multicolumn{1}{c}{105}& 910\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 28\hfill \end{array}\right)\mathrm{MPa},\hfill \\ G_{A\alpha }^{\mathrm{M}}=\left(\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right)\mathrm{N}/\mathrm{mm},\hfill \\ D_{\alpha \beta }^{\mathrm{M}}=\left(\begin{array}{cccccc}67\hfill & 213\hfill & 0\hfill & 0\hfill & 0\hfill & -145\hfill \\ \multicolumn{1}{c}{213}& 3150\hfill & 0\hfill & 0\hfill & 0\hfill & 435\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 252\hfill & 436\hfill & -146\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 436\hfill & 3150\hfill & 213\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & -146\hfill & 213\hfill & 67\hfill & 0\hfill \\ \multicolumn{1}{c}{-145}& 435\hfill & 0\hfill & 0\hfill & 0\hfill & 252\hfill \end{array}\right)\mathrm{N}.\hfill \end{array}$$

(16)

In order to conduct computational asymptotic homogenization procedure, discretization of the target RVE may play an important role, especially for $D_{\alpha \beta }^{\mathrm{M}}$ . Therefore, global mesh size of RVE is identified by h-convergence, we emphasize the use of the FEM providing the monotonous convergence. In Appendix 1, the mesh convergence is shown for one of the cases considered in this study (see Figure 13).

3.1. Infill density variation

In the first part of the numerical study, the effect of infill density on second-gradient parameters has been investigated, varying infill density from 10% to 100% for the same type of infill pattern to be depicted in Figure 1. For a better visualization, how the infill density manipulates the microstructure, we demonstrate in Figure 4 four different infill ratios for the same RVE as used in the analysis.

Figure 4.

RVE samples in the study of infill density variation. Light gray denotes the void and dark gray symbolizes the PLA material. Obviously, (a) very thin structures in 10% infill ratio are impracticable in engineering design; however, (b) 40% is in common use. Because of reduction in weight, (c) 80% and (d) 100% are seen as “conservative” choices among the 3D printing practice.

With the aforementioned approach, all variations are calculated and we provide an insight in Figure 5, how the constitutive parameters, $D_{11}^{\mathrm{M}},D_{33}^{\mathrm{M}},D_{12}^{\mathrm{M}},$ and $D_{34}^{\mathrm{M}}$ , vary with respect to infill ratio. These parameters are two on-diagonal and two off-diagonal terms in $D_{\alpha \beta }^{\mathrm{M}}$ providing a general tendency. Obviously, all the presented parameters initially increase with increasing infill density, reaching their peak values around infill density of 80%. Then, parameters decrease and all vanish when the infill density is 100%. It is of importance to emphasize that ${\mathit{D}}^{\mathrm{M}}$ parameters vanish in the case of homogeneous material of 100% infill (without any microstructure). Therefore, we obtain a non-monotonous dependency on the infill ratio. This property is indeed valuable since a tailored mechanism may be initiated using a varying microstructure with high and low parameters within the body.

Figure 5.

Constitutive parameters $D_{11}^{\mathrm{M}},D_{33}^{\mathrm{M}},D_{12}^{\mathrm{M}},$ and $D_{34}^{\mathrm{M}}$ as functions of infill density.

We stress that ${\mathit{D}}^{\mathrm{M}}$ parameters are relating strain-gradient terms’ significance in the energy description. In other words, one may predict for on-diagonal terms a “stiffness” against corresponding strain-gradient value. For example, $D_{11}^{\mathrm{M}}$ provides normal strain gradient along the same axis as the normal strain. In a beam bending problem, we expect a linear increase in the normal strain and this takes some of the deformation leading to a stiffening effect in smaller beams. This so-called size effect is apparent depending on these values. Interestingly, we realize that one specific choice of infill ratio may trigger an increase in such effects.

In addition, energy due to strain-gradient parameters has been elucidated by conducting simulations with the determined constitutive parameters considering the specimen under shear and tensile tests. We refer to [59,60] for implementation details. For both cases, we divide the energy to a first-order and second-order term and compare them. Specifically, we use equation (4) and write

$$e_1^{\mathrm{M}}=\frac{1}{2}{\epsilon }_{\mathrm{ij}}^{\mathrm{M}}{C}_{\mathrm{ijkl}}^{\mathrm{M}}{\epsilon }_{\mathrm{kl}}^{\mathrm{M}},e_2^{\mathrm{M}}=\frac{1}{2}{\epsilon }_{\mathrm{ij},k}^{\mathrm{M}}{D}_{\mathrm{ijklmn}}^{\mathrm{M}}{\epsilon }_{\mathrm{lm},n}^{\mathrm{M}},$$

(17)

that reads

$$\begin{array}{c}{e}^{\mathrm{M}}=e_1^{\mathrm{M}}+e_2^{\mathrm{M}},\hfill \end{array}$$

(18)

since ${\mathit{G}}^{\mathrm{M}}=0$ effected by the centro-symmetric microstructure. In this way, we obtain a ratio, $e_2^{\mathrm{M}}/e^{\mathrm{M}}$ providing the role of second gradient term in the overall deformation. This ratio depends on the microstructure as well as on the loading scenario as provided in Figure 6 as functions of infill density, along with the change in constitutive parameter $D_{22}^{\mathrm{M}}$ .

Figure 6.

Energy ratio, $e_2^{\mathrm{M}}/e^{\mathrm{M}}$ , denoting the role of strain gradient in the overall deformation in shear and tensile tests with respect to the infill ratio. As a reference, $D_{22}^{\mathrm{M}}$ equally as in Figure 5 is given as well.

Furthermore, the effect of homothetic ratio, $ϵ$ , has been studied. The parameters ${\mathit{D}}^{\mathrm{M}}$ depend directly on this value as apparent in equation (11). For $ϵ\to 0$ , the microstructure vanishes. For comprehending the effect of the homothetic ratio in a beam bending, the mechanical behavior of the specimen has been simulated with $ϵ=1$ and $ϵ=2$ . Both solutions are plotted in Figure 7 for 40% infill density, where a more dominant microstructure, $ϵ=2$ shows a higher stiffening effect than $ϵ=1$ . From an engineering perspective, the homothetic ratio being equal to 1 means exactly the given RVE for 40% in Figure 2. In the case of homothetic ratio being equal to 2, it is scaled by 2, in other words, all dimensions are doubled. Infill density (infill ratio) is the same, yet the microstructure is greater so the effects are more dominant.

Figure 7.

Variation of homothetic ratio and its stiffening effect for a larger microstructure, $ϵ=$ 1 (red, front figure) has four holes along the thickness of the beam and $ϵ=2$ (green, back figure) has two holes along the thickness as the RVE is scaled by 2.

3.2. Infill pattern variation

In the second part of the numerical study, the effect of infill pattern on the strain-gradient constitutive parameters has been investigated. To this end, infill pattern has been changed while keeping infill density fixed as ratio 80%. In order to control the form of voids, the following generalized ellipsoid equation [68]:

$${\left(\frac{|x-x_0|}{r}\right)}^s+{\left(\frac{|y-y_0|}{r}\right)}^s=1,$$

(19)

has been used, where $s$ is the exponent, which controls the convexity of the generated shape, and $r$ is the positive real number. Non-convex shapes are generated by values $s<1$ while block-like shapes are generated by values $s>2$ . Using different values of $s$ , different patterns have been generated, and corresponding $r$ values providing the same infill density have been computed. In Table 1, corresponding $r$ values are listed for the considered $s$ values in this study. The eight different void shapes used in this part are presented in Figure 8. As seen from generated shapes, grid infill pattern is obtained for higher values of exponent $s$ . For each pattern, we applied the asymptotic homogenization procedure and investigated how constitutive parameters ( $C_{\mathrm{AB}}^{\mathrm{M}}$ and $D_{\alpha \beta }^{\mathrm{M}}$ ) change due to infill pattern.

Table 1.

Parameters characterizing infill patterns.

$s$	0.8	1	1.5	2	3	5	7	10
$r(\mathrm{mm})$	1.895	1.895	1.512	1.387	1.265	1.262	1.255	1.238

Figure 8.

RVEs used in the study of infill pattern variation, variation of one parameter, $s$ , incorporates convex and non-convex patterns as well as anisotropic and isotropic inclusions.

In Figure 9, constitutive parameters, namely, $C_{11}^{\mathrm{M}},C_{12}^{\mathrm{M}},$ and $C_{33}^{\mathrm{M}}$ , are presented as functions of $s$ . Roughly speaking, we observe that they remain the same without a clear trend in altering. This fact is justified since the infill density, equivalently the porosity of the material has been kept constant. The stiffness is not dependent on the microstructure, as in equation (11), it is independent of the homothetic ratio, $ϵ$ . Yet there is slight variation in the values, specifically the terms $C_{11}^{\mathrm{M}}$ and $C_{12}^{\mathrm{M}}$ tend to increase with increasing $s$ values up to around $s=2$ and then slightly decrease for values $s>2$ . We stress that the microstructure is introducing some level of anisotropic properties that vanish as $s=2$ is reached (see Figure 8).

Figure 9.

Constitutive parameters $C_{11}^{\mathrm{M}}=C_{1111}^{\mathrm{M}}$ , $C_{12}^{\mathrm{M}}=C_{1122}^{\mathrm{M}}$ , and $C_{33}^{\mathrm{M}}=C_{1212}^{\mathrm{M}}$ as functions of $s$ .

In Figure 10, the constitutive parameters $D_{\alpha \beta }^{\mathrm{M}}$ are given as functions of $s$ . Here, the terms $D_{11}^{\mathrm{M}},D_{45}^{\mathrm{M}},$ and $D_{55}^{\mathrm{M}}$ decrease initially with increasing $s$ values, reaching a minimum value at $s=2$ (circular void shape). Then, they tend to increase rather moderately with increasing $s$ values. This observation may be explained by the anisotropic to isotropic and then again to anisotropic change. However, a similar trend is observed for the terms $D_{22}^{\mathrm{M}}$ and $D_{44}^{\mathrm{M}}$ , having minima at $s=3$ . Hence, it is challenging to find a definite answer how the microstructure shape is altering the response. Moreover, other terms (i.e., $D_{16}^{\mathrm{M}},D_{26}^{\mathrm{M}},D_{33}^{\mathrm{M}},$ and $D_{35}^{\mathrm{M}}$ ) increase with increasing $s$ values.

Figure 10.

Constitutive parameters $D_{\alpha \beta }^{\mathrm{M}}$ variations effected by the microstructure shape change given by $s$ : (a) $D_{22}^{\mathrm{M}}=D_{221221}^{\mathrm{M}}$ and $D_{44}^{\mathrm{M}}=D_{112112}^{\mathrm{M}}$ , (b) $D_{16}^{\mathrm{M}}=D_{111122}^{\mathrm{M}}$ and $D_{35}^{\mathrm{M}}=D_{121222}^{\mathrm{M}}$ , and (c) $D_{11}^{\mathrm{M}}=D_{111111}^{\mathrm{M}}$ , $D_{26}^{\mathrm{M}}=D_{221122}^{\mathrm{M}}$ , $D_{33}^{\mathrm{M}}=D_{121121}^{\mathrm{M}}$ , $D_{45}^{\mathrm{M}}=D_{112222}^{\mathrm{M}}$ , and $D_{55}^{\mathrm{M}}=D_{222222}^{\mathrm{M}}$ .

Based on the extended Voigt notation, $D_{22}$ and $D_{44}$ matrix elements represent the variation in the orthogonal direction to normal strain. The fact that these two elements of constitutive parameter are similar, results from the orthotropic nature of the macro structure Clearly, it is possible to manipulate the macroscopic response, yet the steering mechanism is counter-intuitive. In a design, we have to rely to computation models to find an optimum topology.

In addition, energy due to strain-gradient terms has been calculated for each infill pattern, deforming the specimen under a shear force as before. In Table 2, energy percentages are given for each $s$ value considered in the analysis. Considering the change in energy values, it is seen that the stationary point is around $s=3$ just like in the variation of $D_{22}$ and $D_{44}$ parameters. Therefore, at least, we may claim that these values have a major effect on the strain-gradient energy.

Table 2.

Energy ratio, $e_2^{\mathrm{M}}/e^{\mathrm{M}}$ , showing the significance of the strain-gradient energy in the total deformation energy through variation of the microstructure shape parameter $s$ in a shear test.

$s$	0.8	1	1.5	2	3	5	7	10
Energy (in %)	4.57	4.06	3.59	3.38	3.17	3.37	3.43	3.43

Moreover, the mechanical behavior of the specimen designed with three different $s$ values, namely, $0.8$ , $1.5$ , and $3$ , has been studied, deforming the specimen in a shear test (Figure 11). The $s$ values have been chosen to emphasize the transition from convexity to concavity, considering the geometry of microstructure. According to the bending behavior of each specimen, it may be observed that an increase in concavity of the microstructure positively affects the stiffness of the structure.

Figure 11.

Deformed shapes obtained in shear test for different infill patterns showing the size effect variation regarding infill patterns.

4. Conclusion

Strain-gradient modeling in additively manufactured metamaterials is key for an effective modeling of mechanical response. We focus on different infill characteristics and their relevance in constitutive parameters. In order to identify strain-gradient constitutive parameters, an asymptotic homogenization procedure has been adopted and applied through RVE of each considered infill density and pattern. In the first part of the numerical study, for grid infill pattern, the effect of infill density on the constitutive parameters has been analyzed. Then, in the next part, the effect of infill pattern has been investigated, designing different void shapes parametrized by the exponent term of the adopted generalized ellipsoid equation.

Such studies shed light on the possibilities how to alter the macroscale mechanical response by varying microstructure properties like shape and porosity. We emphasize that the analysis has been carried out by envisioning that a varying infill density and pattern is possible with nowadays manufacturing possibilities. We give a visionary example in Figure 12 cartooning a semi-gradual change in $s$ parameters. By developing adequate optimization algorithms for such problems, we may claim that it is apparent why we need to determine a priori the strain-gradient parameters.

Figure 12.

Functionally graded metamaterial concept as possible tailored material properties to be determined by an optimization algorithm.