Mechanical analysis of heterogeneous materials with higher-order parameters

Even though heterogeneous porous materials are widely used in a variety of engineering and scientific fields, such as aerospace, energy-storage technology, and bio-engineering, the relationship between effective material properties of porous materials and their underlying morphology is still not fully understood. To contribute to this knowledge gap, this paper adopts a higher-order asymptotic homogenization method to numerically investigate the effect of complex micropore morphology on the effective mechanical properties of a porous system. Specifically, we use the second-order scheme that is an extension of the first-order computational homogenization framework, where a generalized continuum enables us to introduce length scale into the material constitutive law and capture both pore size and pore distribution. Through several numerical case studies with different combinations of porosity, pore shapes, and distributions, we systematically studied the relationship between the underlying morphology and effective mechanical properties. The results highlight the necessity of higher-order homogenization in understanding the mechanical properties and reveal that higher-order parameters are required to capture the role of realistic pore morphologies on effective mechanical properties. Furthermore, for specific pore shapes, higher-order parameters exhibit dominant influence over the first-order continuum.


Introduction
Heterogeneous porous materials, such as concrete, metallic foams, and biomaterials are unique subsets of heterogeneous composites with a porous structure. Due to the diverse pore structure, they exhibit a wide range of mechanical, thermal, and electrical properties, making them appealing in the aerospace industry, energy-storage technologies, and bioengineering among other engineering fields. Over the years, different homogenization techniques were developed to analyze the effects of pore morphology on effective material properties. Regardless of the manufacturing techniques, controlling pore shapes and/or distribution remains a challenging task. Consequently, most researches in homogenization methodologies at the macroscopic length-scale (macroscale) have focused on analyzing an existing porous structure rather than designing a new one with controllable microscopic length-scale (microscale) to achieve a desirable macroscopic length-scale (macroscale) properties. However, due to the advancements in additive manufacturing, we have gained a better understanding of these parameters in designing and analyzing tailored porous materials [6,13,21,50,51,54,86].
Different homogenization methods can be classified into analytical and numerical methods. Most analytical methods for porous solids were derived from research on heterogeneous composites, assuming a porous solid as a special case of a two-phase material. Analytical methods can be further delineated into indirect methods based on work by Eshelby [32] and variational methods. Universally used indirect methods are self-consistent scheme by Hill [48], the generalized-self-consistent scheme introduced by Christensen and Lo [22], Benveniste [15], the differential methods of Hashin [45], and the Mori-Tanaka method by Mori and Tanaka [65]. While most widely used variational methods are Voigt by Voigt et al. [89], Reuss by Reuss [79], and Hashin-Shtrikman [80]. Although indirect methods and variational methods are extensively used to obtain effective material properties, different approaches are necessary for a particular shape of inclusion. Indirect methods encounter problems when dealing with arbitrary microstructures as such models lack deterministic studies. Furthermore, these models usually fail to accurately estimate the role of the structures, where phases have large differences in their respective material properties. Variational methods also have limitations as they only account for the volume fraction of different phases and they do not account for the geometry of the microstructure. Hence, analytical methods are not suitable for precise analysis of effective material properties, and thus are not useful in the development of new purpose-built porous materials.
Numerical methods are usually based on averaging different fields, such as stress, strain, and deformation energy density. The average and the calculation of local field quantities are carried out by solving the underlying physics problem within a so-called representative volume element (RVE) to model the microstructure. Based on continuum formulation, numerical methods can be separated into two groups: methods based on classical continuum theory, or first-order multiscale computational homogenization approach, and methods based on generalized continuum theory, or second-order continuum generally referred to as second-order multiscale computational homogenization [69,70]. Although first-order methods are well established in the computational homogenization community and are widely used by researchers to account for porosity and pore shape such as He et al. [46], the first-order approach has its disadvantages. Particularly, the first-order approach requires strict separation of scales, and also adherence to the continuum mechanics concept of local action. This in turn results in an inability to capture the absolute size of the pores and pore distribution or to deal with localization problems. To overcome these problems, the first-order method has been extended to a second-order approach, where the use of a generalized continuum enables us to introduce the length scale into the material constitutive law to capture the pore size and distribution [38,63]. The ability of second-order homogenization methods to accurately capture pore morphology (e.g., pore sizes, pore distribution, and pore shapes) significantly advances our understanding of the interaction between pore structure and material behavior. This allows us to numerically analyze purpose-built porous materials with complex microstructures.
For application purposes, additive manufacturing, also called 3-D printing is capable of constructing microstructure resulting in a significant change in macrostructure response, called "metamaterials." A well-known example in composites is an isotropic microstructure causing an anisotropic behavior at the macroscale. Herein, we consider a higherorder effect arising from the microstructure and indicating the importance of "metamaterials".
Additive manufacturing is capable of producing metamaterials by infill ratio settings, adding texture on the surface, or by creating multi-architecture systems by design [37,53,64,88,92]. Complex multiscale structures are typically designed by setting a lattice of porous material [52,63,64,77], where complexity of such structures and effects of local microstructure are best described by incorporating highergradient effects. This can be observed in biological materials [39,44,57,83] that are highly heterogeneous (e.g. bone tissue) and have hierarchical structures that demand highergradient effects to be included in modeling their mechanical responses [40,42] as well as in their multi-physics simulations [39,62]. We stress that the homogenization approach is necessary because of the computational cost. A porous structure may be simulated by considering all details of the microstructure, but the computational time and memory requirements for a successful computation make this choice impractical.
For metamaterials and their description, the scientific scrutiny is still ongoing [7,41,67]. The use of homogenization techniques for metamaterials is under development as well [56,66]. Thus, there exist different homogenization approaches [33,43,55,71,72], by means of direct approaches [3,36,78] and also computational homogenization methods [29,75,76,84]. These methods are computationally costly. Semi-inverse methods have been developed as well, based on gamma-convergence [4,5] and by means of asymptotic analysis [14,23,49,85]. Asymptotic analysis approach has been applied in one-dimensional problems for reinforced composites [12,18], in two-dimensional continuum [10,19,25,73], and for formal analysis [30,34,35,87]. In this paper, we apply an asymptotic homogenization technique developed in Abali and Barchiesi [1] to investigate the effects of pore size, shape, and distribution on effective material properties of heterogeneous porous media. By means of ample studies [8,11,58,59,74], we have shown that this method is applicable for metamaterials. Herein, we make the connection that porous materials are metamaterials such that we must consider higher-order terms in their analysis.
The rest of the paper is organized as follows. The higher-order asymptotic homogenization method and computational implementation are briefly introduced in the second section. Numerical results and discussion of Young's modulus, Poisson's ratio, and higher-order parameters are presented in the third section, followed by the conclusion.

Asymptotic homogenization method in generalized continua
Let us start by defining a domain Ω that represents a representative volume element (RVE). Then we assert that deformation energies within Ω at the macroscopic length scale (macroscale) and the microscopic length scale (microscale) are equal where Φ is deformation energy density (per volume), "M" and "m" represent macroscale and microscale, respectively, and dV is an infinitesimal volume element expressed in Cartesian coordinates as dV = dx dy dz . To satisfy Eq. (1), the size of the domain Ω must be large enough to effectively capture the underlying microstructure. It should be noted that a unit cell represents the microstructure and its characteristic length is lower than the threshold for a generalized continuum model. An RVE is at least one unit cell and periodically repeated RVEs lead to the entire structure in the macroscale. However, they are expressed in the same coordinate system without scale separation such that we may assert energy equivalence. Energy densities at macroscale and microscale are different, since we have distinct displacement fields, u m and u M , implying that the microscale deformed position x m is no longer equal to the macroscale deformed position x M . This condition, x m ≠ x M , is to be depicted in Fig. 1 and the displacement is defined as a deviation from its reference position, X , as follows: At the microscale, deformation energy is expressed by a linear elastic model (first-order continuum). By following [61], we know that the energy density will have higherorder terms at the macroscale (higher-order or generalized continuum). We start making simplifications and assume that the microscale is accurately defined by a linear strain measure, m . Moreover, the material relation between stress and strain is linear such that we obtain a quadratic microscale deformation energy density, Here and henceforth we use Einstein's summation convention over repeated indices and a usual comma notation for space derivative. Furthermore, by specifying minor symmetries of the stiffness matrix C m ijkl = C m jikl = C m ijlk , we obtain: Conversely, deformation energy density at the macroscale is described by strain gradient model (second-order continuum) where alongside ℂ M we additionally account for higher-order parameters: higher-order elastic constants M and coupling constant M as described by Liebold and Müller [60]. Considering minor symmetries

macroscale deformation energy reads as
If microscale is represented by a centro-symmetric RVE then M = 0 , see dell'Isola et al. [26]. Furthermore, if we assume M = 0 , then the solution reverts to the classical linear elastic macroscale model (first-order continuum) without higher-order terms [2]. Inserting expressions for microscale and macroscale deformation energy from Eqs. 4 and 5 back to Eq. (1), energy equivalence becomes: The rest of this section will be devoted to the solution of Eq. (6). First, we will solve the right-hand side integral (macroscale energy equation) to extract ℂ M , M , and M parameters. Next, we will apply a two-scale expansion of the microscale displacement field (asymptotic expansion) and solve the left-hand side integral (microscale energy equation). Lastly, by comparing the left-hand side and the right-hand side solution, we will link ℂ M , M , and M to the solution of the microscale energy equation and obtain the homogenized values.

Solution for macroscale energy of RVE
We emphasize that ℂ M , M and M are constant in strain and strain gradient, To solve Eq. (7), a geometric center of the RVE, c X , has been employed, Assuming displacement field u M is continuous over the microscale, we can approximate the macroscale displacement by a Taylor expansion around the value at the geometric center by truncating terms with orders higher than quadratic as shown in Yang et al. [91]. Macroscopic displacement field and its displacement gradients are evaluated as: In Eq. (9), the first and second derivatives of macroscopic deformation field are the unknowns. They are obtained by spatial averaging and exploiting the fact that terms evaluated at c X are constant within Ω , thus, Going back to the Eq. (9), we replace the displacement gradients with spatially averaged values from Eq. (10), as follows: Finally, we import the above averages into the macroscale energy definition on the right-hand side of Eq. (6) and take spatial averaged terms out of the integral.
where In the following section, we solve microscale deformation energy from Eq. (6) and compare it to the Eq. (12) for extracting the values for ℂ M , M and M .

Solution for microscale energy of RVE
Asymptotic homogenization methods as described in Pinhoda Cruz et al. [24] and Abali and Barchiesi [1] are used to approximate the microscale deformation energy of the RVE. The basic idea behind the two-scale expansion is to formulate microscale deformation energy in terms of the gradients of the macroscopic deformation field. To do that, we first introduce homothetic ratio (see Fig. 2) as a transformation between RVE and unit cell in a so-called local coordinate.
Microscale displacement field for the RVE is then expanded with regard to as: where n u(X, y) (n = 0, 1, 2, ...) are y-periodic as an assertion from the unit cell. In order to solve the left-hand side of Eq. (6), we impose equilibrium condition ∇ ⋅ + m f = 0 on the RVE as: where m and f i are microscale mass density and external forces, respectively. By substituting local coordinate y into Eq. (15) and using the chain rule, we obtain the first derivative of microscale displacement field, Substituting Eq. (17) into Eq. (16) and again using chain rule, we obtain an asymptotically expanded governing equation: Comparing coefficients in Eq. (18) of the same order of leads to: ijkl depends only on local variable y . This leads to a straightforward conclusion: After substituting Eq. (22) in Eq. (20), and introducing y-periodic function abc = abc (y) , we define the following partial differential equation: General solution of Eq. (20) with the respect to Eq. (23) is given as: Repeating the same procedure for Eq. (21) and introducing y-periodic function abci , we define the following partial differential equation: General solution of Eq. (21) with the respect to Eq. (25) is given as: Therefore, the microscale displacement field can be rewritten as: From the left-hand side integral in Eq. (6), we observe that we need to obtain displacement gradient.
Substituting macroscale gradients from Eq. (11) into Eq. (28) and neglecting terms higher than the second gradient, as in connection to the simplification introduced in Eq. (17); we acquire: Now we can substitute microscale displacement gradient in Eq. (6) with Eq. (29) and derive the solution of microscale deformation energy: As it was said before, we can now determine macroscale (homogenized) material parameters ℂ M , M and M by comparing Eq. (12) to Eq. (31) : Eq. (32) is focal output from higher-order asymptotic homogenization that gives us homogenized material properties. Numerical solution of Eq. (32) will be briefly explained in the next section.

Numerical solution and implementation in FEniCS
Calculation of macroscale parameters in Eq. (32) requires the solution of abc and abcd tensors from Eqs. (23) and (25). This calculation is done by following the standard procedure of the finite element method, where we define a weak formulation for Eqs. (23) and (25). Those weak forms are solved with the FEniCS platform [17], as shown in Figure 3.
We have defined a CAD model in the open-source software SALOME platform [16] and utilized Netgen algorithms for generating the mesh. Models in 2D were discretized with triangle elements. Standard periodic boundary conditions were imposed on the model, where nodes on opposite boundaries have the same solution for tensors and , see Figure 4.

Model setup
In this section, we are going to verify and demonstrate the predictive capability of the higher-order asymptotic homogenization. For this purpose, several case studies with diverse concrete morphologies, such as different combinations of inclusion/pore volume fractions, shapes and distributions, were generated and their mechanical responses were investigated. First, validation was preformed by comparing the numerical results to experimental data for porous concrete reported by Yaman et al. [90] and Du et al. [28]. Next, two sets of problems were solved to investigate the influence of structural heterogeneity on homogenized material properties, such as Young's modulus, ℂ M , M , and M .

Problem description
To study the behavior of homogenized material parameters, two types of problems were designed. The first problem was designed so that inclusions were treated as pores and surrounding material was defined as intrinsic concrete (zero porosity), which enables us to study the influence of porosity, pore shape, and pore distribution on homogenized material properties. In the second problem, inclusion properties were defined as a fraction of matrix material properties, where matrix properties remained the same as the previous case (intrinsic concrete). Thus, we analyze three distinct cases: (i) inclusions have the same material properties as a matrix-homogeneous case, (ii) inclusions are defined as a fraction of matrix material properties, where they simulate different mixtures of cement, aggregates (coarse & fine), air, and waterheterogeneous composite case, (iii) inclusions are defined as pores.
Additionally, the influence of RVE size on the homogenized material properties was investigated. This systematic study  Table 1 shows the properties obtained from the literature [90] that were used in our numerical models. It should be noted that plane stress models were studied in all cases.

Mesh convergence
To determine mesh convergence, the Richardson Exploration method [81] was utilized to obtain a higher-order estimate ("exact solution") of the continuum value from a series of lower-order discrete values. Furthermore, gird (mesh) convergence index GCI was calculated to consistently report the results of convergence studies and provide an error band on the convergence of the solution. The GCI is based upon a mesh refinement error estimator derived from the theory of generalized Richardson Extrapolation and can be computed using two levels of the mesh; however, three levels are recommended to estimate the order of convergence accurately and to check that the solutions are within the asymptotic range of convergence.
Single circular pore RVEs with coarse, medium, and fine mesh were used for Richardson exploration and convergence study was done on three RVE sizes: 1 × 1, 2 × 2, and 3 × 3. As shown in Table 2 the convergence study was explicitly focused on the higher-order parameters of tensor, specifically D 111111 . Values for the exact solution from Richardson exploration, CGI fine error for the fine mesh, and asymptotic range of convergence index ARCi are also listed in the Table 2. It should be noted that ARCi value close to one indicates that the solution is well within the asymptotic range of convergence. Moreover, CGI error value less than 0.1% and ARCi index of ∼ 1, shows that solutions obtained on the fine mesh are mesh independent and well within the asymptotic range of convergence.

Validation
Any model is an approximation of the actual physical system, thus validation against verified experimental data is a necessary step. To the best of the authors' knowledge, there are no experimental data on higher-order parameters ( M and M ), which is mainly due to the lack of available experimental procedures to measure these parameters. In light of this, Young's modulus and Poisson's ratio are used to validate the accuracy of the higher-order homogenization approach, as they are readily available in the literature. Specifically, experimental data for porous concrete with porosity range, from 9% to 21%, has been utilized. In addition to the experimental data, results were also validated against conventional empirical models: Voigt, Hashin and Shtrikman (referred to as "Hashin model"), and Gibson and Ashby (referred to as "Ashby model"), Ashby and Medalist [9]. Voigt model shows the theoretical maximum for Young's modulus, while Hashin and Ashby models are able to characterize a two-phase material with relatively good accuracy. In order to mimic this material closely, a single circular pore RVE has been used in both examples. Voigt and Hashin models are defined as: where E m and E p are Young's modulus of matrix (e.g., intrinsic concrete) and pore, respectively, while v m and v p are their corresponding volume fractions. Unlike Voigt and Hashin models, Ashby model in addition to Young's modulus and volume fraction, has two additional constants, C and n, defined as follows: where C is a geometrical constant of proportionality, which is set equal to one in our case, C = 1 , Salvini et al. [82]. Constant n is an empirical exponent depending on pore morphology, where for open pores, n = 3∕2 , and n = 2 for closed ones [82]. In this study, we set n = 2 as we are dealing with closed pore in all our models.
Comparison between homogenized values of Young's modulus E h and Poisson's ratio h with experimental data and empirical models are shown in Fig. 5a, b, where one can see the obtained results from the homogenization model are in good agreement with the experimental data and empirical models. Moreover, in Table 3 we show (33)  macroscale form (first-order continuum), it is be reasonable to assume that data for higher-order parameters is also in good agreement with physical reality. In the subsequent sections, we are going to apply homogenization model to investigate the influence of structural heterogeneity. First, on material parameters coming from linear elastic model ( ℂ M ), and then material parameters additionally accounted by strain gradient model ( M and M ).

Influence of pore shape and distribution on homogenized mechanical properties
To explore the role of pore shape and distribution on homogenized mechanical properties, we first explicitly look at the role of pore shape (i.e., circle, square, ellipse, and triangle), while increasing the porosity from 0 to 50%. We then explore how the presence of multiple pores with different distributions (i.e., uniform and random) will impact the overall mechanical response, while keeping the porosity constant among all cases. Increase in porosity, results in a reduction of homogenized Young's modulus, see Choren et al. [21] and Herakovich and Baxter [47]. This is a well known effect, and thus we avoid additional discussion in this section. Furthermore, according to the first-order homogenization approach, RVE size does not affect values of ℂ M , therefore RVE size was kept constant at 1 × 1. To have a fair comparison between uniform and random distribution in multi-pore models, porosity was kept constant at 20%. Figure 6 shows the influence of porosity and pore shape on homogenized Young's modulus, where clear distinction between each cases have been highlighted. More an in-depth analysis of these results indicates that circular and square pores produce cubic behaviour ( C M 1111 = C M 2222 ), while elliptical and triangular pores produce orthotropic material behaviour ( C M 1111 ≠ C M 2222 ). Furthermore, square and circular pores have a similar influence on the reducing the value of the Young's modulus in comparison with triangular and elliptical pores. This difference in material behaviour is a product of different symmetries that each shape possesses. Since porosity values higher than 30% violate the boundary conditions by intersecting the edges, there is no values for triangular pores after such threshold.
In Table 4, we show the difference between homogenized values of Young's modulus for a single pore and multi pore distributions at 20% porosity. From the presented data, it is evident that pore distribution for the same porosity has little to no effect on the homogenized value of Young's modulus, a similar conclusion was also observed by other researchers, Carr et al. [20], where the authors also reported that pore geometry has no influence on Young's modulus, which is different from our observations. The authors came to this conclusion because they use a 3D SEM (scanning electron microscope) image for their FEM model, where underlying morphology is  defined by random size and spatial distributions of pores, thus obscuring the influence of pore shape. In our work, both uniform and random distributions have the same pore shapes and we have also noticed that uniform distribution does not change material behavior (e.g. uniformly distributed elliptical pores still produce orthotropic behavior), while random distribution creates anisotropic behavior independent from the pore shape.
Observations made in this case study, and the ability of the asymptotic homogenization method to obtain them could potentially have a wide industrial applications. Pore influence analysis allows for accurate prediction of mechanical properties from knowledge or characterization of the pore structure. This, in turn, allows for the fabrication of tailored components with a defined porous structure that satisfy a variety of applications, from biomedical implants to lightweight structures for the space industry. With rapid advances in additive manufacturing that allow for the design and manufacturing of ordered porous materials, such needs become even more apparent.

Influence of inclusion's material property on homogenized mechanical properties
Now that we have systematically explored the impact of the structural features on overall mechanical behavior, it is time to understand how material properties of inclusion would impact such properties. Therefore, we borrow the same models from the previous section where the underlying structures remain the same while we are varying the modulus of inclusion. This becomes specifically interesting when inclusion material is fluctuating between matrix material properties and pore. Figure 7 shows the influence of the inclusion shape and material properties on Young's modulus. Results reveal two distinct regions: one is highly influenced by both structural and material properties, while the other is only influenced by the material properties. Specifically, when the property of inclusion is larger than 50% of the property of the matrix, the homogenized value of Young's modulus E h is only impacted by the material's property. However, as the property of inclusion becomes smaller than 50% of the matrix's property, the homogenized value of Young's modulus, E h , is not only influenced by the property of the inclusion but also by the structural property of the inclusion (i.e., inclusion shape). Note that the volume fraction of the inclusion in all these cases remains the same. Figure 7 also shows that shape will have maximum impact on E h once inclusion transitions to a pore. Furthermore, from the case study on the pore shape/ distribution, we have learned that pore distribution has little to no influence on Young's modulus, which also true once pore transitions to inclusion.
Another noteworthy finding from this study is that orthotropic and cubic material behavior caused by inclusion/pore shape starts to shift to isotropic behavior as Young's modulus ratio goes from 0 to 1 (e.g. porous to homogeneous material). To illustrate this shift in material's behavior, we have used anisotropy measure introduced by Zener [93], which quantifies the anisotropy of cubic crystals. Zener anisotropy index is defined as: where C 11 , C 12 , and C 44 are the three independent components that collectively represent the elastic modulus tensor. Fig. 7 Decreasing E inclusion ∕E matrix ratio introduces shape influence on the homogenized value of Young's modulus, E h . Symmetric RVEs (circle and square) have a similar influence on the degradation of Young's modulus as E inclusion ∕E matrix is decreasing, but there are significant differences when comparing them to triangular and elliptical pores Figure 8 shows Zener anisotropy index for circular pore/ inclusion case, where the shift from initial cubic behavior is easily noticeable.
Summarizing the results from Figs. 7 and 8, two basic conclusions are derived: • E inclusion ∕E matrix < 0.5 inclusion shape starts to influence homogenized value of Young's modulus and overall material behavior (culminating once inclusion becomes a pore). • Anisotropy is introduced by random distribution, yet it rapidly diminishes once E inclusion ∕E matrix > 0.5.
This observation is of importance for explaining the lack of anisotropicity in the literature. Namely, a random distribution is often overseen since we know from polycrystals that this case is still isotropic in the homogenized continuum. Yet the polycrystal structure has moduli ratio equals to 1. Herein, we observe the anisotropic behavior for distinct materials, where pore and matrix have significantly different stiffness values. This observation can be potentially used in tailoring composite materials where fiber direction plays the primary role in determining material properties. From observations obtained from this case study, one can conclude that alongside fiber direction, fibers cross-sectional shape can also plays a role in dictating the homogenized material properties. This influence becomes even more pronounced when dealing with porous composite materials as pore shape has a strong influence on material properties and material's behavior.

Higher-order parameters
In this subsection, we exclusively discuss the role of higherorder parameters on previous case studies. Unlike first-order parameters where RVE size does not influence the material parameters, higher-order parameters, ( M and M ), are influenced by the RVE size, inclusion/pore shape, and inclusion/pore distribution. Such influence is easily explained by the homothetic ratio [see Eq. (31)], where values of will change linearly and values of will change quadratically by an increase in the RVE size as shown in Table 5. It should be noted that 5 one can observe negative values in M , which may be comprehended as a treat to the positive definiteness. This condition is automatically fulfilled in the presented homogenization method since we have started with the assertion that the strain energy is equivalent at the microscale to its corresponding model at the macroscale. By using a positive definite microscale model, herein positive values in ℂ m tensor, we assure that the positive definiteness is fulfilled. Thus, we emphasize that the energy at the macroscale with all parts related to ℂ M , M , and M must be positive definite [31,68].
Since values of for centro-symmetric RVEs become zero, D 111111 component of tensor has been selected as an example to illustrate the effects of porosity, pore shape, and pore distribution on higher-order parameters. Other components of tensor show similar impacts, which can be easily examined. Figure 9 highlights how the porosity impacts the D M 111111 for a single circular pore case. Based on Eq. (32, the homogenized value of D M 111111 is a function of D 111111 and homogenized value of ℂ M . Similar to C M 1111 , value of D 111111 also decreases as the porosity increases. However, this reduction initially occurs at a faster rate compared with C M 1111 . Thus, value of D M 111111 is increasing until approximately 20% porosity and then decreasing to 0 as porosity reaches to 100%. Similar behavior is observed for other components. This trend indicates that there is a critical porosity value at which M components have Fig. 8 The higher-order asymptotic homogenization method is capable of capturing the shift in material behavior. Cubic material behavior shifts to isotropic material behavior as E inclusion ∕E matrix goes from 0 to 1. Data is for single circular inclusion/pore RVE largest absolute value. Figure 10 shows similar analysis as Fig. 9 but for a square pores. As one can see, the change in the pore shape results in completely different behavior in D M 111111 as porosity reaches 100%. Moreover, C M 1111 shows similar behavior for both circular and square shape pores, which is expected based on the results and discussion from Sect. 4.1. However, with an increase in porosity, D 111111 first decreases, then starts increasing exponentially as porosity reaches higher values. Such exponential growth of D 111111 consequentially increases homogenized value of D M 111111 . Similar behaviour is observed in pantographic metamaterials, whose first-gradient energy is nearly zero. Thus, deformation is dominated by the strain gradient part of the energy, which calls for pantographic metamaterials to be described by a second-gradient continuum model, see dell'Isola et al. [27]. Figure 11 illustrates the variation in the D M 111111 for constant porosity value of 20%, the change in number of pores, pore morphology and orientation for three different RVE sizes are highlighted. Figure 11a, b and c show that depending on the shape and distribution of pores, values of D M 111111 are different, where uniform distribution has the smallest value and single pore has the largest value. This difference comes from the size influence, as it was shown in Table 5. Looking at random distribution, it can be concluded that higher-order parameters will oscillate between those two extremes, as random distribution statistically can become uniform distribution with a single pore due to coalesce of the pores. From random distribution values in Fig. 11d, one can see that circular pore has the biggest influence on the higher-order parameters for the specified porosity.
Comparing the influence of inclusion/pore shape on ℂ M and M values for different cases, several interesting observations can be obtained. First of all, there is a critical porosity at which M has the biggest influence, that can be observed from the analysis of circular and square pores in Figs. 9 and 10. Specifically, as one can see from Fig. 9, the parameters in M for the circular pore case shows an optimum value around 20% porosity, whereas values in ℂ M monotonously decreases with an increase in the porosity. This optimum value is of interest to maximize higher-order effects through simple design choices. Analogous study shown in Fig. 10 reveals that changing circular to rectangular pores influenced this optimum value and instead of 20% porosity, the maximum higher-order effects are achieved as we are approaching theoretical porosity of 100%. Moreover, taking a closer look at symmetric shapes, such as circles and squares, one can see that they have very similar Young's modulus values. In contrast, their respective M parameters are significantly different. This indicated that first-order homogenization is unable to capture the shape effects. The aforementioned observation becomes important when analyzing metamaterials for bending rigidity, where pore shape has a significant influence. Furthermore, Young's modulus is not sensitive to the distribution of inclusions/pores, while M parameters are highly sensitive. The ability of higherorder homogenization becomes important in designing and analyzing complex problems where understanding damage and fracture initiation and propagation is influenced by the inclusion/pore distribution.
Overall, with a broad interest in porous materials and metamaterials in various scientific and engineering applications, there is an absolute need to study the influence of these higher-order parameters that go beyond just a stiffness tensor. There is also a need to be able to probe these impacts, which demands advancements in existing experimental capabilities.

Conclusion
In this paper, we investigated the effects of pore morphology on the effective mechanical properties of heterogeneous porous materials using higher-order asymptotic homogenization. Specifically, the effects of pore morphology are assessed through macroscale stiffness matrix ℂ M and higherorder parameters M and M that are explicitly computed by assuming linear elastic material model at the microscale.
In the adopted framework, the higher-order asymptotic homogenization is implemented in the FEniCS platform, and was used to solve the partial differential equations generated from the homogenization procedure. Mesh convergence demonstrates the stability of the methodology, and numerical results from the validation problem show a good agreement with experimental data.
To demonstrate the predictive capability of the method, several case studies with diverse concrete morphologies and different combinations of inclusion/pore volume fractions, shapes, and distributions were generated, and their mechanical response was investigated. The study of the pore influence revealed that pore shape significantly influences Young's modulus and material's behavior. In contrast, pore distribution has little to no effect on Young's modulus although introduces anisotropy. Furthermore, problems where the pore is replaced by inclusion showed that the inclusion shape influences the homogenized value of Young's modulus once there is a sufficient difference between the material properties of inclusion and matrix. At first glance, case studies on the interaction between Young's modulus and pore morphology seem straightforward, but that was not actually the case. Taking a closer look at symmetric shapes, such as circle and square, and pore distribution showed that they have very similar Young's modulus values, while their respective higher-order parameters have a significant difference. This indicates the lack of first-order homogenization and the need for higher-order parameters in some application. For instance, analyzing metamaterials for bending rigidity, stress concentration, damage, and fracture, where pore shape and distribution play a significant role, need a better understanding of these parameters. Importance of higher-order parameters can be further emphasized when looking into biomechanics, especially when dealing with mechanics of bones and processes of bone tissue growth where multi-scale organization and high heterogeneity requires application of generalized continuum theories together with higher-order homogenization. Even though these numerical results for higher-order parameters still need experimental validation, these findings provide invaluable insights to study the interaction between underlying pore morphology and material properties. In summary, the results of this study can bring us one step closer to designing hierarchical heterogeneous porous materials/metamaterials with desirable macroscopic properties.
Funding Open access funding provided by Uppsala University.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.