Influence of microstructure on size effect for metamaterials applied in composite structures

Microstructure related deviation from elastic response is known as ‘‘size-effect.’’ Metamaterials – for example modeled by strain gradient elasticity – capture this effect adequately by means of additional parameters to be determined. We employ a methodology based on asymptotic homogenization in order to obtain metamaterials parameters and then present the influence of these additional parameters by using simulations. By means of the finite element method, we solve metamaterials deformation modeled by the strain gradient elasticity. The implementation is established by open-source packages (FEniCS) for a realistic, composite structure with round and oval inclusions.


Introduction
It is known that microstructure affects the mechanical response at a larger length-scale.A typical example is found in composite materials, where a high tensile glass or carbon fiber is used as a reinforcement in the epoxy matrix, the alignment of fibers create an anisotropic response.At microscale, fibers and matrix are isotropic; however, at macroscale the response is anisotropic.This microstructure related response is of importance, when porous materials design is aimed for [1][2][3][4] or mechanical response of such structures is searched for [5][6][7].An efficient method for this analysis is known as homogenization.As formally analyzed in [8][9][10][11], a microstructure exhibits higher order parameters in mechanical response at the macroscale.One of such models is strain gradient theory and its roots go back to a mathematical analysis for involving microstructure [12][13][14].In order to utilize this method, different assumptions are proposed [15][16][17][18], discussed [19][20][21], and verified experimentally [22][23][24].In the strain gradient theory, additional parameters are used to compensate the microstructure effect.We stress the analogy of isotropic materials at microscale for fiber reinforced composites, where ''additional'' anisotropic terms appear in the stiffness tensor.In the strain gradient theory, these parameters are not only in the stiffness tensor, but also in the tensor in front of strain gradients.It is difficult to differentiate between ''material'' and ''structural'' parameters such that we call this continuum body a ''metamaterial''.
In addition to a formal analysis perspective, recent developments in sintering and additive manufacturing also demonstrate the potential of metamaterials as shown in [25][26][27] as well as in [28][29][30].A design may involve on purpose an inner substructure or a chosen production method causes an anisotropic microstructure by adding texture or infill ratio in (especially metal) 3-D printing [31][32][33].
We briefly show an asymptotic homogenization method and determine metamaterials parameters.Asymptotic analysis [34][35][36][37] has already been applied in [38][39][40][41] for such problems.As a homogenization method, this approach is employed in one-dimensional problems for reinforced composites [42,43] and in two-dimensional continuum [44][45][46][47] by using numerical solutions.As studied in [48][49][50][51], this method is adequate for the strain gradient theory.Such multiscale structures exist in biological materials [52][53][54] as well as man-made designs [55][56][57].In this work, we employ the asymptotic homogenization method for determining all parameters such that a computation is possible in a homogenized continuum.In order to clarify the idea, we depict an SEM micrograph in Fig. 1 from recycled rubber-mortar, where different rubber particles are added to the cement paste in order to control effective properties.We aim at predicting the effective properties including higher gradient parameters as the main outcome in this work.

Parameter determination
In this work we follow the methodology introduced in [58] and the implementation verified in [59] with applications in [60].Metamaterial is composed of a microstructure expressed in  denoting space.Its corresponding homogenized continuum, macroscale, is expressed in .Macro-to-microscale transformation is handled by a so-called homothetic ratio, .In this way, we circumvent a scale separation and use https://doi.org/the same coordinate system for both length-scales.Furthermore, this approach allows to obtain the metamaterial parameters for a specific microstructure modeled in specific geometric dimensions.Then the same parameters are valid for any macroscale geometry composed of such microstructure.
The approach is based on the ''known'' microscale leading to the ''sought after'' parameters for the macroscale.We use linear elastic isotropic model with the known stiffness tensor,  m .For a porous material, the stiffness varies in space depending on void or material.In the case of a composite material, the stiffness tensor depends on space, , as well, where the value is either fiber or matrix related.At the microscale, the deformation energy is modeled as a quadratic one, Herein we use standard continuum mechanics formulation and understand summation convention over repeated indices.For the sake of simplicity, we use linearized displacement gradients as a strain measure at the microscale, By using Castigliano's theorem,  m ∕ m  , for stress, we obtain Hooke's law.A linear strain measure is used and the stored energy density is quadratic in displacement gradients.Therefore, we use small displacement for the homogenization process leading to metamaterials parameters.Once the parameters are found, at the macroscale, geometric nonlinearities may be involved by using a nonlinear strain measure.For material nonlinearities, the proposed method needs to be revisited.In the homogenized continuum at the macroscale, we model by using the strain gradient elasticity with the same type of (quadratic) stored energy density [61] A comma notation denotes a partial derivative in .The same strain measure, at the macroscale, results in Strain gradient theory involves additional metamaterials tensors,  M ,  M , and  M .The microstructure is reduced to a unit cell, where the repetition of this unit cell composes the topology.Therefore, a representative volume element (RVE) is at least one unit cell and by using periodic boundary conditions on the RVE, we capture the complete continuum body's geometry.The whole methodology uses one axiom that the deformation energy within the RVE, , is identical at micro-and macroscales For the complete derivation of the governing equations, once more, we refer to [58].We introduce a field  and acquire its numerical values by solving Once the field  is calculated, we obtain metamaterials parameters, In the classical homogenization method, used also for composite materials, this stiffness tensor is the outcome.However, we emphasize that we search for  M and  M as well.Therefore, we introduce another field  and solve its governing equation, where  m is the known mass density at microscale and  M = ∫   m d is the (constant within the RVE) mean value at macroscale.With this solution, we calculate and obtain all strain gradient parameters ) . ( We observe that  M ,  M , and  M depend on  0 ,  1 , and  2 , respectively.Thus, we comprehend that their role in deformation vanishes in first and second order for  M and  M in the case of larger geometries.Especially for  M , such a quadratic dependence on geometric characteristic length is known as the size effect.

Numerical implementation
All aforementioned governing equations will be solved numerically.We use open-source packages known as FEniCS [62,63] and solve the governing equations by generating a weak form for each of them.This procedure is often called a variational formulation, we refer to [64] for engineering applications and their implementation.The numerical calculation is established by the finite element method [65], where the field functions are approximated by their discrete couterparts with a compact support.Technically, we span a finite dimensional Hilbertian Sobolev space for trial functions.Owing to Galerkin approach, the same space is used for test functions.Since micro-and macroscales are both expressed in the same coordinate system, we use the same discretized space within the domain.The computational domain is the RVE.The domain is discretized by nodes and their connectivity by elements.This standard discretization is called a triangulation that is generated by NetGen algorithms in Salome CAD software.
We solve the discrete problem by minimizing the weak form obtained after multiplying by a corresponding test function and integrating by parts.In 3-D (2-D) case, for  from Eq. ( 6), we obtain 6 (3) Trial functions,  and , as well as their test functions,  and , are all approximated by linear form functions.Periodic boundary conditions are applied on all boundaries.Hence,  and  solutions are equivalent on corresponding boundaries.For example, on the left boundary the solution of  is mapped to the solution at the right boundary, in other words, the degree of freedom (DOF) on each node on the left is set equivalent to the DOF on the right with the same coordinates.This restriction is applied on the higher-level directly for creating the space for trial functions.This restriction is necessary to obtain -periodic  and  fields leading to a -periodic displacement that is the condition to be fulfilled for an RVE.Since all weak forms are linear, we solve them directly by using Mumps solver from PetSc packages [66].
We use two different RVEs with spherical and oblate inclusions.In 2-D continuum, these inclusions are circles (round shape) and ellipses (oval shape), for a more general approach, we distribute them randomly.The matrix is out of mortar approximated by mass density,  = 1900 kg/m 3 , Young's modulus,  = 20 000 MPa, and Poisson's ratio,  = 0.2.The inclusions are rubber with  = 840 kg/m 3 ,  = 2 MPa, and  = 0.45.We obtain two set of parameters  M ,  M ,  M for both RVEs, their CAD models and corresponding mesh are shown in Figs. 2 and 3.
Determined parameters are given by using Voigt notation, where ,  indicate {11, 22, 12} and ,  denote {111, 221, 121, 112, 222, 122}.Comparison of parameters for round and oval inclusions reveal that along the horizontal,  1 -axis, the anisotropic character is more pronounced in oval inclusions than round.This result is expected as the shape of inclusion is more dominant than the random distribution.It is challenging to examine the role of strain gradient parameters,  M ,  M .Therefore, we utilize these parameters and solve one simple example of beam bending for analyzing the influence of these parameters on the mechanical response.We minimize the energy involving second gradient of displacement at macroscale,  M .Hence, for the computation, we use quadratic form functions.In order to enforce a continuous and differentiable field across the elements, we use a mixed formulation and solve derivative of displacement as an additional unknown, ∇.For ensuring consistency,  M , = ∇  , we use a Lagrange multiplier, , as an additional unknown.In order to discretize the multiplier, we use constant (discontinuous) form functions.In this way, we employ strain and its gradient from two unknown fields, and generate the weak form, as follows: For technical details, we refer to an analogous implementation described in [67].This weak form solves displacement, gradient of displacement, and multiplier in a coupled manner.The unknowns are tensor rank 1, 2, and 2, respectively.Therefore, for a 2-D problem 2 + 4 + 4 = 10 degrees of freedom (DOFs) and for a 3-D problem 3 + 9 + 9 = 21 DOFs are used in each node.Despite the computational cost, this implementation allows to introduce  1 continuity, which ensures a monotonous convergence in the strain gradient theory solved Table 1 Beam bending and the size effect.by the finite element method [68,69].On Dirichlet boundaries, the displacement is given such that the test function vanishes.Hence, only on Neumann boundaries,  N , traction vector, t, is defined.We use FEniCS packages for solving this mixed spaces finite element problem.
As a solver, we use SNES line search algorithm.We stress that the use of a mixed space with different types of elements creates a numerical difficulty in solvers based on Krylov subspaces, hence, we use a direct solver umfpack for the linear problem without making use of Krylov subspaces.

Results
As it is evident, the chosen periodic structure is representing a random distribution of inclusions.The oval inclusions are longer along the horizontal axis leading to a higher stiffness in this direction.For a direct comparison, we construct a simple example of a bending beam with a length to thickness ratio of 10.In all simulations, this ratio remains the same.The beam is clamped on left and right.In the middle, we use a traction of 10 MPa applied in the middle on 10% of the length of the beam.By varying the thickness, and consequently length, we observe two experimentally known phenomena.First, the anisotropy leads to a stiffer beam along the horizontal (neutral) axis along  1 in oval inclusions such that the bending rigidity is, as expected, higher.Second, smaller the beam more dominant is the effect of second gradient.This stiffening behavior is known as size effect.Quantitative results are compiled in Table 1, where the maximum vertical displacement (deflection in the middle of the natural axis) is normalized by the thickness, ūmax .The stiffness tensor,  M , steers the maximum deflection.As  M  11 in Voigt notation is greater in oval inclusions than in round, normal stress is less leading to a smaller deflection.Higher order terms,  M ,  M , manipulates the size effect.As obvious, size effect is captured in both examples; however, their values are different.In order to quantify this difference, we simply use the normalized values and compare results in thicknesses 100 mm with 3 mm.For round inclusions, we obtain (7.37 − 4.78)∕7.37= 35% and for oval inclusions (6.57− 5.06)∕6.57= 23%.This decrease depends directly on the higher order terms determined by the inner substructure.By using examples, it seems unfeasible to determine all parameters.Yet by using such a simple beam bending analysis for different specimen sizes, it is possible to validate results obtained herein.We emphasize that the first order theory leads to the same (normalized) deflection in all cases.
Not apparent in the first sight, but we obtain an additional insight to bounds of this theory.All parameters are obtained for the unit cell of 1 mm characteristic length.As seen in Table 1, for larger than 100 mm sizes, the higher order effects are insignificant.We may call 100 mm as the upper threshold for the strain gradient theory.Approximately 3 mm size is the lower threshold.We demonstrate the deformation in ParaView by scaling it 10 times for the sake of visualization in Fig. 4. By using the same mesh of 200 × 20, we observe that the smooth displacement field is experiencing numerical disturbances at the lower threshold.This observation is not justified by the element size, but rather the geometric size.An easier interpretation is trying to simulate the microstructure without periodic boundaries, which is equal to the choice of 1 mm thickness.Naturally, one expects that 1 mm would be the absolute minimum possible where the homogenization would lead to wrong results.In this suggested way, as we know that the weak form is adequate for a smooth displacement field, we may determine the lower threshold for the underlying microstructure.

Conclusion
Strain gradient theory is understood as a necessary tool for metamaterials, i.e. structures with an inner substructure.We demonstrate direct results made possible by the recently developed computational framework based on asymptotic homogenization.By using randomly distributed round and oval inclusions, we determine parameters,  M ,  M ,  M , and use them in simulations.By this methodology, we acquire an understanding on how to distinguish between first order parameters  M and second order parameters,  M ,  M .Namely, for the beam bending, the former controls the bending rigidity (as it is known in any textbook in mechanics of materials).The latter is of importance on how the size effect changes the stiffness.Furthermore, quite naturally, we estimate bounds for the strain gradient theory.Despite the fact that the strain gradient parameters are different for different types of inclusions, for the chosen unit cell (herein also RVE) of 1 mm, we have obtained the lower threshold around 3 mm and upper threshold around 100 mm.In other words, for a design with such a metamaterial in sizes within 3 − 100 mm, strain gradient theory must be taken into account.In the worst case, up to 23% error in round inclusions and up to 35% error in oval inclusions will occur in beam bending.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Scanning electron microscopy (SEM) image of a concrete (cement with aggregates) structure.Light gray denotes the mortar with approximately 20 GPa of modulus, whereas the dark gray ''inclusions'' are rubber with an approximately 2 MPa of modulus.

Fig. 2 .Fig. 3 .
Fig. 2. Round inclusions, denoted by orange (representing mortar) and green (rubber), RVE is 1 mm × 1 mm.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig. 4. Beam bending shown by scaling the deformation 10 times, top: thickness of 100 mm, bottom: thickness of 3 mm, colors denote the magnitude of displacement.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 10.1016/j.mechrescom.2022.103877Received 5 November 2021; Received in revised form 28 February 2022; Accepted 4 March 2022