Uncertainty Quantification of the Mechanical Properties of 2D Hexagonal Cellular Solid by Analytical and Finite Element Method Approach

TL;DR

This paper compares analytical and numerical methods to evaluate the mechanical properties of a 2D hexagonal cellular solid, highlighting discrepancies and uncertainty.

Contribution

The study introduces a systematic comparison of analytical and finite element methods for mechanical properties under geometric and material variations.

Findings

Analytical and numerical results for Young’s modulus showed percentage differences ranging from −8.28% to 10.87%.

Yield strength discrepancies varied between −2.86% to −5.53% for segmental length and 7.76% to 10.57% for thickness.

Entropy measures like Bhattacharyya and Kullback–Leibler were used to assess uncertainty in homogenized properties.

Abstract

The mechanical properties of cellular materials are critical to their performance and must be accurately determined through both analytical and numerical methods. These approaches are essential not only for understanding material behavior but also for evaluating the effects of parametric variations within the unit cell structure. This study focuses on the in-plane comparison of analytical and numerical evaluations of key mechanical properties, including Young’s modulus, yield strength, and Poisson’s ratio of a 2D hexagonal unit cell subjected to systematic geometric and material variations. Analytically, the mechanical properties were derived based on the geometric configuration of the hexagonal unit cell. Numerically, the finite element method (FEM) simulations employed three different meshing methods: quadrilateral, quad-dominated, and triangular elements, to ensure precision and consistency in the results. The elastic response (Young’s modulus) was examined through a parametric sweep involving segmental length variations (4.41 to 4.71 mm) and material modulus (66.5 to 71.5 GPa), revealing percentage differences between analytical and numerical results ranging from −8.28% to 10.87% and −10.58% to 11.95%, respectively. Similarly, yield strength was evaluated with respect to variations in segmental length (4.41 to 4.71 mm) and wall thickness (1.08 to 1.11 mm), showing discrepancies between −2.86% to −5.53% for segmental length and 7.76% to 10.57% for thickness. For Poisson’s ratio, variations in the same parameters led to differences ranging from −7.05% to −12.48% and −9.11% to −12.64%, respectively. Additionally, uncertainty was assessed through relative entropy measures—Bhattacharyya, Kullback–Leibler, Hellinger, and Jeffreys—to evaluate the sensitivity of homogenized properties to input variability. These entropy measures quantify the probabilistic distance between core material distributions and their effective counterparts, reinforcing the importance of precise modeling in the design and optimization of cellular structures.

Figures

Full text

1. Introduction

Cellular materials are widely observed in natural [1,2,3,4] and engineered structures [5,6,7] as a hexagonal structure that exhibits various mechanical properties, including high strength due to their light weight, high energy absorption, and thermal insulation due to their structural configuration or integrity. Various studies discussed cellular material as a hexagonal structure in the form of honeycombs, lattice structures, and foam, which play a significant role in multiple conditions [5,6,7,8,9,10,11,12,13,14,15].

Several studies have investigated the mechanical properties of hexagonal cellular structures. Yu et al. (2023) [16] identified topological defects, specifically five to seven defects, in hexagonal honeycombs through observations of heterogeneous structures in nature, such as graphene, insect eyes, and beehives; Swellam et al. (1997) [10] examined the mechanical behavior of foam cells within a series of hexagonal configurations; Habib et al. (2018) [8] analyzed the compressive behavior and strain energy of nine types of honeycomb structures, including the hexagonal configuration; Liu et al. (2021) [17] numerically investigated the in-plane crushing behavior of hexagonal honeycombs, focusing on the influence of different Poisson’s ratio values; and Carlsson et al. (2022) [18] provided a numerical solution for the elastic–plastic response of hexagonal unit cells subjected to uniaxial tensile and compressive loading. In light of these literature studies, the present work conducts analytical and numerical homogenization to evaluate the Young’s modulus and yield strength of the hexagonal unit structure.

Homogenization techniques play a crucial role in the complex geometry and boundary conditions of cellular materials. Various approaches to homogenization, such as analytical and numerical methods, have been used to determine the effective mechanical properties of cellular materials [19,20,21,22,23]. Analytical homogenization based on the assumption of periodic microstructure and idealization of the unit cell, such as square, hexagonal, tetrakaidekahedral, etc. The basic principles of mechanics, such as beam theory [24,25,26], strain energy method [27,28,29], Gibson–Ashby Model [30,31], asymptotic [21,32,33], and closed-form expressions for effective moduli [17,34], were used to derive a mathematical equation as a function of geometric parameters and base material properties. Most analytical homogenization of the effective mechanical properties of cellular materials yields linear results due to the following assumptions: uniform material distribution, linear elasticity, and a regular periodic structure. While numerical homogenization often employs the finite element method (FEM), providing a more general and accurate framework, this approach considers the representative volume element (RVE) and subjects it to specified boundary conditions and various meshing analyses, including quad, quad-dominated, and triangular meshes [23,35,36,37,38].

In cellular material, the geometric parameters of the unit cell, such as length [24,39,40,41], nodal connection [41,42,43], angle between nodes [4,25,44], thickness of the cell wall [31,39], and its width, play an essential role in the mechanical properties of the entire cellular solid. Additionally, the mechanical properties vary with the base material, which can be homogeneous or non-homogeneous [8,25]. In both cases, geometric topology and base material properties are crucial for optimizing and designing the cellular material; in this way, analytical and numerical homogenization are essential for predicting the effective mechanical properties as factors vary.

Cellular materials are widely studied for their excellent energy absorption in vehicle crashworthiness applications [45], but challenges remain in optimizing geometry and material combinations for consistent crash performance. Fatigue failure of platform screen doors under subway aerodynamic loads further highlights the need for structural systems that can dissipate energy and resist cyclic stresses [46]. Recently, hybrid friction energy-absorbing structures (HFEAS), comprising a bearing bush, anti-creeper device, pretension bolt, friction plate, and frictional metal/CFRP tubes, have emerged as a promising solution by combining frictional dissipation with structural compliance [47]. However, the integration of frictional mechanisms with cellular architectures remains unexplored, leaving a clear gap in developing multifunctional systems that combine crashworthiness, fatigue resistance, and reliable energy absorption for enhanced transportation safety. Within this context, 2D hexagonal cellular solids are particularly attractive because of their simple geometry, tunable mechanical properties, and resemblance to natural honeycomb structures. However, in practical applications, uncertainties inevitably arise from variability in material properties, manufacturing imperfections, and geometric irregularities, all of which can significantly influence their performance and reliability [48]. Accurate uncertainty quantification of the mechanical properties of hexagonal cellular solids is therefore crucial for robust design and safe deployment in engineering systems. Despite extensive studies on their deterministic mechanical behavior, limited research has systematically addressed stochastic aspects, especially through combined analytical and finite element method (FEM) approaches [49,50,51]. This gap motivates the present work, which utilizes Monte Carlo simulation, stochastic perturbation technique, and semi-analytical methods to perform an uncertainty analysis of both analytical and numerical models of cellular solids.

This study presents a systematic framework for the uncertainty quantification of the mechanical properties of 2D hexagonal cellular solids using a combined analytical and finite element method (FEM) approach. Variability in parameters such as cell wall thickness, elastic modulus, and geometric irregularities is modeled and propagated to quantify the stochastic response of the effective elastic and strength properties. Analytical homogenization methods are employed for computational efficiency, while FEM simulations provide high-fidelity verification and benchmarking of the probabilistic predictions. The originality of this work lies in bridging deterministic homogenization with stochastic modeling, highlighting the accuracy and limitations of analytical theories under uncertain inputs.

2. Analytical Model

The analytical approach for calculating material parameters for the 2D hexagonal unit cell is based on the elastic bar model. It is divided into three independent testing phases, each devoted to determining the elastic modulus, yield strength, and deriving the effective Poisson’s ratio. The mathematical model for each phase is based on the assumption that the elastic phase of the skeleton characterizes the geometric ability to return to its original shape after deformation, the plastic phase describes the yield strength, and that Poisson’s ratio quantifies the material’s lateral strain relative to axial strain under applied stress. The von Mises metal plasticity model is applicable in this study to model nonlinear deformations of such a solid. The geometric model of the skeletal material under consideration is shown in Figure 1. Let us adopt the following notation: $[eqn]$ is the analytical Young’s modulus, $[eqn]$ is the Young’s modulus of the base material, $[eqn]$ is the analytical yield strength, $[eqn]$ is the yield strength of the base material, $[eqn]$ is the tensile yield strength, $[eqn]$ is the compressive yield strength, $[eqn]$ is the analytical Poisson’s ratio, $[eqn]$ is the thickness of the cell wall, $[eqn]$ is the segmental length of cell, $[eqn]$ is the external length of unit rectangle, θ is the angle b/w applied force the segment or beam of the hexagon, $[eqn]$ is the hexagon internal height in vertical direction, $[eqn]$ is the width of hexagon on horizontal direction, $[eqn]$ is the height of unit rectangle, $[eqn]$ is the applied unit force, $[eqn]$ is the horizontal component of applied force, $[eqn]$ is the vertical component of applied force, $[eqn]$ is the longitudinal strain and $[eqn]$ is the transverse strain.

The homogenized Young’s modulus of the two-dimensional hexagonal lattice structure was determined using the analytical geometric homogenization method (also known as the micromechanical analytical approach) [21,52]. In this approach, the hexagonal cell was considered as a periodic representative volume element (RVE), and the overall stiffness was derived based on the geometric configuration and load transfer mechanisms of the cell struts. Let us begin with the mathematical expression describing the horizontal and vertical elements of the cell with applied deformation in the segmental length of hexagon ( $[eqn]$ ) as

$$[eqn]$$

Therefore, the effective Young’s modulus in terms of the homogenized Young’s modulus ( $[eqn]$ ) and effective area ( $[eqn]$ ) is

$$[eqn]$$

Thus, the final analytical equation of the effective Young’s modulus proceeds using the following equivalence:

$$[eqn]$$

Hence, the analytical geometric method provides a simplified yet accurate estimation of the equivalent Young’s modulus for cellular structures.

The effective yield strength [7] of the 2D hexagonal lattice was also derived using the analytical geometric method combined with the yield criterion approach. The cell walls were modeled as slender beam elements subjected to combined axial and shear stresses. By applying the principle of static equivalence and the von Mises yield criterion, the effective analytical yield strength of the lattice can be expressed as

$$[eqn]$$

Thus, the maximum force at which the cellular structure exhibits irreversible yielding ( $[eqn]$ ) is can be expressed as

$$[eqn]$$

Substituting the value of force ( $[eqn]$ ) from (5) into (4), one can develop the final equation for yield strength as [7]

$$[eqn]$$

The homogenized Poisson’s ratio for the hexagonal lattice structure was also obtained through the analytical geometric method by analyzing the deformation kinematics of the hexagonal cell. Under uniaxial loading, the lateral deformation was derived from the rotation and bending of the inclined cell walls. The axial deformation $[eqn]$ with applied force $[eqn]$ , bending moment on the $[eqn]$ with segmental length $[eqn]$ , and end rotation $[eqn]$ is expressed as

$$[eqn]$$

The transverse deformation $[eqn]$ and the macroscopic strain with the regular hexagonal height $[eqn]$ is expressed as

$$[eqn]$$

The simplification of Equation (8) in terms of axial $[eqn]$ and transversal $[eqn]$ strain is

$$[eqn]$$

Since the term $[eqn]$ , we neglect this term and Equation (9) is reframed as

$$[eqn]$$

The final equation for Poisson’s ratio is given below. According to Ongaro, F. (2018) [49], the value of this coefficient C for a hexagonal structure is 0.33; therefore, Equation (10) becomes

$$[eqn]$$

These closed-form equations for the effective material parameters of the elastic–plastic homogeneous and isotropic equivalent solid enable a relatively easy probabilistic model, where the solid skeleton material and geometrical parameters become random variables and/or fields.

3. Numerical Analysis

Consider a static equilibrium of a statistically homogeneous and bounded continuum $[eqn]$ . Let us assume further that there are non-empty subsets of the external boundaries of this region $[eqn]$ , namely $[eqn]$ and $[eqn]$ , where the Dirichlet and von Neumann boundary conditions are defined, correspondingly. They result in a displacement field and stress tensor components that satisfy the linear elasticity elliptic boundary value problem given by Equations (12)–(16). It should be emphasized that the proposed homogenization method does not require any Dirichlet boundary conditions in contrast to the approach presented in (Kamiński, 2005; Kamiński, 2013) [53,54]; they are introduced for a brevity of presentation. Considering the proposed numerical technique, the entire set of homogenization boundary value problems with the same boundary conditions and with additionally modified input variables (Kamiński, 2005) [53] is to be solved, where n represents the number of discrete realizations of the mean value. The upper index α indicates the different structural responses associated with these input values. The solution to the particular boundary differential equation systems that describe the static equilibrium near the average value of this parameter is sought:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The elasticity tensor of the core material satisfies the symmetry, boundedness, and ellipticity conditions, and is defined as

$$[eqn]$$

where $[eqn]$ and $[eqn]$ denote the Young’s modulus and Poisson’s ratio, respectively, with $[eqn]$ . The following series variational formulation is introduced and discretized to obtain an appropriate numerical solution with the finite element method:

$$[eqn]$$

The left-hand side of Equation (18) corresponds to the elastic behavior of the structure, and the right-hand side is equivalent to the applied stress boundary conditions. It should be noted that the indexing with α, in addition to the Response Function Method (RFM), should be incorporated into the computational domain Ω for modeling shape sensitivity; the corresponding extension and additional conditions may reflect fluctuations in this external boundary structure. The determination of the effective material tensor requires the introduction of the strain energy for a heterogeneous medium, i.e.,

$$[eqn]$$

The homogenized medium accumulates the same amount of energy as a medium with a series of effective elastic characteristics $[eqn]$ , so that we can compare this with the strain energy stored in the original medium:

$$[eqn]$$

where $[eqn]$ are the macro strains in the equivalent macro-homogeneous body. The following boundary conditions are applicable in this case:

$$[eqn]$$

and, independently, along the x2 axis:

$$[eqn]$$

The macro strains within the RVE can be determined as

$$[eqn]$$

and reduced to x1 and x2 by only, while assuming additionally uniform size l in both directions.

Then, the following system of linear algebraic equations is solved symbolically in the MAPLE system to compute the effective characteristics:

$$[eqn]$$

where the strain energy is equal to

$$[eqn]$$

in the case of biaxial tension of the RVE. The computational implementation in the framework of the finite element method is relatively straightforward and follows directly from traditional approaches (see, e.g., Kamiński, 2013) [54,55,56,57,58,59].

We provide in turn a mathematical basis for the Least-Squares Method (LSM) adjacent to the fourth-order tensor in both the non-weighted (NLSM) and weighted (WLSM) versions of the homogenization method (Kamiński, 2013) [54]. We use a polynomial approximation of the sth order (indexed by β here) through n numerical tests of the homogenization problem solved near the mean value of the given design parameter $[eqn]$ . As a result, we obtain n different pairs $[eqn]$ $[eqn]$ We seek the following polynomial representation of the effective tensor versus the given design parameter:

$$[eqn]$$

We introduce for this purpose the residuals at each trial point and each component of the homogenized tensor, i.e.:

$$[eqn]$$

The coefficients $[eqn]$ are found by minimizing the weighted residual function:

$$[eqn]$$

Such a minimization can be performed as follows (Björck, 1996) [57,60]:

$$[eqn]$$

Adopting the following notation for the Jacobian:

$$[eqn]$$

one can get

$$[eqn]$$

Thus, we obtain the following matrix of normal equations:

$$[eqn]$$

This system of equations is solved symbolically in the system MAPLE for three components of the effective elasticity tensor, namely $[eqn]$ , $[eqn]$ , and $[eqn]$ , while determining the orthotropic effective medium for the given hexagonal material. The situations simplify in the case of an isotropic homogenized medium, which is investigated below, and where direct approximation of Young’s modulus and Poisson’s ratio is completed.

4. Computer Simulation

The analytical and numerical homogenization approaches for a hexagonal unit cell are employed here to determine the effective mechanical parameters, including Young’s modulus, yield strength, and Poisson’s ratio. This analysis was performed by introducing slight variations, both increases and decreases, in key geometric parameters such as segment length, cell wall thickness, and material properties of the unit cell. Analytical approaches have been implemented directly in the symbolic algebra package MAPLE 2025.1, whereas the FEM analysis was delivered in the system ABAQUS/CAE 6.14-5.

Some specific geometric dimensions of the hexagonal unit cell used in the numerical simulations were employed to calculate the corresponding values of Young’s modulus, yield strength, and Poisson’s ratio. The segmental length ranged from 4.41 mm to 4.71 mm, and the material’s Young’s modulus ranged from 66.5 GPa to 71.5 GPa. A similar procedure was followed for determining the yield strength and Poisson’s ratio, incorporating segment length variations of 0.3 mm (4.41 to 4.71 mm) and cell wall thickness variations of 0.01 mm (1.04 to 1.14 mm). Three types of meshing were employed in ABAQUS, as shown in Figure 2, namely quadrilateral (Quad, quad-dominated) and triangular meshes, with the principal aim of collecting reliable numerical results that closely correspond to the analytical predictions, thereby validating the accuracy and robustness of the homogenization methodology.

In numerical homogenization, a uniform mesh element size of 0.3 is used across all models. These models are categorized based on mesh element geometry and the number of nodes per element. Each case of the model has distinct mesh geometries: quadrilateral (Quad), quad-dominated, and triangular, with 4, 4, and 3 nodes per element, and for more precise numerical results, each finite element is further divided into 8, 8, and 6 nodes per element in the quadrilateral (Quad), Quad-dominated, and triangular cases, respectively. The meshing configurations and number of nodes for all three models are shown in Figure 2 and Figure 3, respectively, and also the total number of mesh elements and nodes of each model as shown in Appendix A.

In all three-phase testing of the numerical model, aluminum material was assigned with Young’s modulus of 69 GPa and a yield strength of 69 MPa in the case of the plastic test. In the numerical measurement of an elastic model, Young’s modulus is determined using the following expression:

$$[eqn]$$

where $[eqn]$ is numerical Young’s modulus, $[eqn]$ is stress on the right and left side, and $[eqn]$ is strain on the right and left side. Finally, Poisson’s ratio is quantified numerically using the following expression:

$$[eqn]$$

Hence,

$$[eqn]$$

where $[eqn]$ and $[eqn]$ represent axial and shear stress, respectively, and $[eqn]$ and $[eqn]$ represent axial and shear strain. Numerically, the yield strength is obtained from the stress–strain curve, where the point at which yielding begins is taken as the yield strength of the model. The stress–strain curves for the various models, generated using different mesh types and varying numbers of nodes, are shown in Figure 4. The axial and shear strains are determined by the distance between nodes 1 and 2 before applying the load, divided by the distance between nodes 1 and 2 after the load (deformed shape), as shown in Figure 5.

The most important notice related to the effective mechanical properties (Young’s modulus, yield strength, and Poisson’s ratio) of the cellular solid exhibit distinct relationships with geometric parameters and material properties when subjected to small incremental variations. Analytically, as described by the mathematical formulation of Young’s modulus in Equation (3), a clear inverse relationship exists between Young’s modulus and the segmental length $[eqn]$ , while a direct relationship was observed with the material’s intrinsic Young’s modulus. Accordingly, as the segmental length increases, the effective Young’s modulus decreases, which is consistent with the trends observed in the numerical simulations as shown in Figure 6. Similarly, an increase in the material’s Young’s modulus $[eqn]$ leads to a corresponding increase in the effective Young’s modulus. Quite similarly, the relationship of yield strength and Poisson’s ratio with incremental variations in segmental length $[eqn]$ and cell wall thickness $[eqn]$ reveals distinct trends in the mechanical response of the hexagonal cellular structure as shown in Figure 7 and Figure 8, respectively. Analytically, these relationships are defined by Equations (6) and (11) for yield strength and Poisson’s ratio, respectively. As shown in Figure 6, Figure 7 and Figure 8, both properties exhibit an inverse correlation with segmental length, indicating that as $[eqn]$ increases, the corresponding values of yield strength and Poisson’s ratio decrease. Conversely, a direct relationship is observed with respect to the cell wall thickness $[eqn]$ , where increases in $[eqn]$ lead to an increase in both yield strength and Poisson’s ratio.

In our numerical simulations, we primarily selected quadrilateral (Quad) and Quad-dominated meshes for evaluating Young’s modulus and yield strength, since these meshes produced results closest to the analytical solutions. This is because Quad elements are generally more accurate for representing axial and shear responses under uniform loading conditions, leading to better convergence in stiffness-related properties such as Young’s modulus and yield strength. However, in the case of Poisson’s ratio, triangular meshes provided results closer to the analytical values. This can be attributed to the ability of triangular elements to better capture localized deformations and lateral strain distributions, which play a critical role in accurately estimating Poisson’s ratio.

The analytical results for Young’s modulus demonstrate a decreasing trend from 6.6725 GPa to 5.9036 GPa as the incremental variation in the segmental length $[eqn]$ from 4.41 mm to 4.71 mm. Accompanying numerical results show a slightly less pronounced decline, from 6.019 GPa to 5.392 GPa. This discrepancy may be attributed to the simplifications inherent in the analytical model, which assumes idealized boundary conditions, whereas the numerical model incorporates more realistic constraints and deformation behavior. Furthermore, when varying the base material’s Young’s modulus $[eqn]$ from 66.5 GPa to 71.5 GPa, the analytically determined effective Young’s modulus increases from 6.2275 GPa to 6.6957 GPa, while the numerical values range from 5.563 GPa to 5.996 GPa. This slight deviation is likely due to mesh sensitivity and stress concentration effects captured in the finite element analysis, which are not accounted for in the analytical expressions.

Although with an incremental variation in segmental length $[eqn]$ from 4.41 mm to 4.71 mm, the analytical results showed a reduction in yield strength from 0.831 MPa to 0.747 MPa, and analogously a decrease in Poisson’s ratio from 0.0808 to 0.0756. Correspondingly, the numerical results show a decline in yield strength from 0.8701 MPa to 0.776 MPa and in Poisson’s ratio from 0.0899 to 0.0810. In the case of incremental variation in cell wall thickness $[eqn]$ from 1.08 mm to 1.14 mm, the analytical yield strength increases from 0.867 MPa to 1.036 MPa, while Poisson’s ratio slightly decreases from 0.0757 to 0.830. And numerically, the yield strength rises from 0.785 MPa to 0.944 MPa, accompanied by an increase in Poisson’s ratio from 0.0866 to 0.0908.

The elastic response was examined through a parametric sweep involving variations in segmental length and material modulus, revealing percentage differences between analytical and numerical results ranging from −8.28% to 10.87% and from −10.58% to 11.95%, respectively. Similarly, yield strength was evaluated with respect to variations in segmental length and wall thickness, showing discrepancies between −2.86% to −5.53% and 7.76% to 10.57%, respectively. For Poisson’s ratio, variations in the same parameters led to differences ranging from −7.05% to −12.48% and −9.11% to −12.64%, respectively.

The discrepancies between analytical and numerical results of the effective mechanical parameters, such as Young’s modulus, yield strength, and Poisson’s ratio of hexagonal cellular material, can be attributed to several factors. In the case of an analytical model, we assumed a uniform material distribution, linear elastic behavior, and idealized geometry that may not fully capture the stress distribution and deformation percentage in the actual structure. While in numerical simulation (FEM), consider geometric non-linearities, different stress distributions, and more accurate localized plastic deformation, which leads to a realistic representation of effective mechanical behavior. Additionally, mesh density, type of mesh element, and boundary conditions in the numerical simulation (FEM) can influence the output result.

5. Development of Polynomial Model

Based on variations in material properties, segmental length, and cell wall thickness in both analytical and numerical models, we obtained multiple sets of results. From these results, six-order polynomial models for all effective parameters (Young’s modulus, yield strength, and Poisson’s ratio) were generated using trend analysis in MS Excel. The resulting sixth-order polynomial expressions are presented below, and the corresponding regression statistics are summarized in Table 1.

The following polynomials have been determined:

Polynomial for Young’s Modulus in the case of variation with material property $[eqn]$ :

$$[eqn]$$

Polynomial for Young’s Modulus in case of variation with segmental length $[eqn]$ :

$$[eqn]$$

Polynomial for Yield strength in case of variation with segmental length $[eqn]$ :

$$[eqn]$$

Polynomial for Yield strength in case of variation with thickness of cell wall $[eqn]$ :

$$[eqn]$$

Polynomial for Poisson’s ratio in case of variation with segmental length $[eqn]$ :

$$[eqn]$$

Polynomial for Poisson’s ratio in case of variation with thickness of cell wall $[eqn]$ :

$$[eqn]$$

The comparative analysis of the regression coefficients (see Appendix B) and residual plots (seen in Figure 9) between the numerical and analytical models reveals that the numerical approach provides a more realistic representation of the effective mechanical parameters. Although the analytical solutions exhibit mathematically perfect precision, as evidenced by their near-zero standard errors and extremely large t-statistics, such results often reflect an idealized system that lacks sufficient accommodation of variability. In contrast, the numerical regressions, while presenting slightly larger error margins, align more closely with physically expected trends and practical observations as shown in Appendix B. For instance, Young’s modulus in the numerical model has the coefficients that are physically interpretable and consistent with realistic material stiffness values. In contrast, the analytical outputs tend to overfit the data space. Similarly, for yield strength, the numerical slope captures the expected sensitivity to geometry while leaving minor but meaningful residuals that reflect the natural scatter of mechanical testing. Poisson’s ratio also shows stronger plausibility in the numerical results, with moderate coefficients that represent realistic lateral contraction effects. At the same time, the analytical regressions tend to exaggerate precision without accounting for system imperfections. Residual graphs reinforce this interpretation: the numerical residuals are randomly scattered around zero without systematic bias, a hallmark of well-fitted but physically constrained models. In contrast, the analytical residuals converge too closely to the regression line, suggesting over-idealization. Therefore, despite the mathematical elegance of the analytical solutions, the numerical regressions provide a more faithful depiction of the material response by accommodating inherent variability and experimental uncertainty, thereby making them more reliable for practical applications.

6. Uncertainty Quantification

In this study, the mechanical response of a 2D hexagonal unit cell was modeled using a sixth-order polynomial function representing the nonlinear behavior of an effective mechanical property under varying microstructural conditions (variation in the core material and geometric parameters). Three different computational strategies were employed: the Monte Carlo simulation (MCS), the stochastic perturbation technique (SPT), and a semi-analytical method (SAM). All uncertainty analyses were conducted across a range of randomness parameters, with alpha varying from 0.025 to 0.275 in increments of 0.025. The Monte Carlo simulation, as the reference approach, was implemented to capture the full probabilistic behavior of the unit cell response under uncertainty. A total of 1,000,000 simulation trials were conducted for each value of α, wherein the input random variables (e.g., base material, geometric parameters $[eqn]$ ) were generated from prescribed probability distributions. The sixth-order polynomial model was evaluated numerically for each realization, allowing for the statistical characterization (mean, variance, higher moments) of the output mechanical property. The method did not involve iterative convergence in the classical sense, but instead relied on statistical convergence across the ensemble, with a sufficient sample size ensuring the accurate estimation of probabilistic measures.

The stochastic perturbation technique was implemented to provide a computationally efficient alternative to the full-scale MCS, and mechanical response function was expanded as a Taylor series about the mean values of the input random variables, capturing the effects of stochastic fluctuations up to higher-order terms in this method. The perturbation order was systematically modified from 1 to 10, allowing for an investigation of the sensitivity of response statistics to the level of approximation. For each $[eqn]$ increment, the perturbation-based expansion yielded approximations to the first and second statistical moments (mean and variance), with higher-order perturbations offering increased accuracy at the cost of analytical complexity. This method required no iterative numerical solver, as derivatives were evaluated symbolically or via automatic differentiation using MAPLE 2025 symbolic engine. Finally, the semi-analytical method (SAM) was employed as a hybrid approach that integrates analytical differentiation with numerical evaluation. This method involved computing partial derivatives of the sixth-order polynomial with respect to input parameters such as material properties and microstructural variables. These derivatives were then embedded in an incremental-iterative framework to simulate the nonlinear response under varying α. A modified Newton–Raphson algorithm was used for the iterative solution of nonlinear equations, with a convergence tolerance $[eqn]$ based on the norm of residual forces.

Across all three methods, the variation in the parameter $[eqn]$ served as a control parameter to examine the influence of microstructural geometry on the global mechanical response. The combination of stochastic and semi-analytical approaches provided both computational efficiency and deep insight into the reliability and robustness of the material system under uncertainty. The comparative analysis of these methods not only validated the consistency of the results but also highlighted the trade-offs between computational cost and solution fidelity. Monte Carlo provided the most accurate statistical results, stochastic perturbation offered efficient moment estimation, and the semi-analytical method enabled deterministic sensitivity analysis with moderate computational effort.

In this study, we also investigate the uncertainty quantification of probabilistic coefficients associated with effective properties. These coefficients are described in terms of statistical parameters (expectation, coefficient of variation, skewness, and kurtosis) and relative entropy measures, specifically, the Bhattacharyya distance, Kullback–Leibler divergence, Jeffreys divergence, and Hellinger distance in a two-dimensional hexagonal cellular unit cell subject to small variations in both base material $[eqn]$ properties and geometric parameters $[eqn]$ . Uncertainty quantification was achieved through the application of a sixth-order polynomial approximation derived from both analytical and numerical homogenization tests.

Let us consider for this purpose $[eqn]$ and $[eqn]$ as the analytical and numerical polynomials of 6th order determined by the generalized Taylor expansion, respectively, as follows:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the coefficients of the analytical and numerical polynomials, respectively, and are determined by curve fitting. According to the Bhattacharyya distance $[eqn]$ , the polynomial approximation is given as

$$[eqn]$$

According to the Kullback–Leibler divergence $[eqn]$ , the polynomial approximation can be proposed in the following way:

$$[eqn]$$

The Hellinger distance $[eqn]$ can be expressed using the following polynomial approximation:

$$[eqn]$$

According to Jeffrey’s divergence $[eqn]$ , the polynomial approximation is given as

$$[eqn]$$

where $[eqn]$ represents the expected value of the analytically derived polynomial, and $[eqn]$ represents the expected value of the numerically derived polynomial. Further numerical analysis of these distances has been delivered under the assumption that probability distributions of both properties of the cellular solid are Gaussian.

The effective Young’s modulus of the cellular material exhibited a clear nonlinear response to increasing microstructural randomness, as captured by a sixth-order polynomial expansion concerning input variability. The randomness parameter α, which ranged from 0 to 0.25 in increments of 0.05, was applied to both material stiffness $[eqn]$ and geometric parameters such as cell length $[eqn]$ . The expectation of the effective Young’s modulus showed a mild decreasing trend with increasing α, indicating that elevated randomness leads to a slight softening of the structure’s average stiffness, as shown in Figure 10.

This effect is attributable to the increased occurrence of weaker configurations within the stochastic ensemble. The coefficient of variation (COV), on the other hand, displayed a clearly nonlinear (near-quadratic) increase, confirming that even small input variability can lead to significant amplification in output variance. Higher-order statistical measures such as skewness and kurtosis also revealed the growing asymmetry and a peak of the modulus distribution with increasing randomness. Specifically, the skewness transitioned from near-zero to positive values, while kurtosis increased steadily, suggesting that the distribution becomes more right-tailed and sharply peaked. These effects were more accurately captured by Monte Carlo simulations and semi-analytical methods compared to stochastic perturbation techniques, which tended to underpredict higher moments at $[eqn]$ This highlights the importance of using comprehensive uncertainty quantification techniques when modeling effective stiffness in architected materials under manufacturing or material variability.

The response of the effective yield strength to increasing input randomness exhibited significantly greater sensitivity compared to elastic stiffness. The yield strength was modeled as a sixth-order polynomial function of both cell length $[eqn]$ and wall thickness, $[eqn]$ with randomness parameterized by the coefficient of variation α. As α increased, the expectation of the yield strength declined sharply, especially when wall thickness variations dominated, due to the pronounced influence of local instabilities and stress concentration effects in thinner or uneven regions. The COV increased more rapidly than for Young’s modulus, exhibiting a pronounced nonlinear trend as shown in Figure 11. This reflects the inherently more unstable nature of plastic deformation, where localized yielding can initiate prematurely under geometrically perturbed conditions. Skewness values became markedly positive at higher levels of randomness, indicating a right-skewed distribution in which a significant number of realizations exhibited much lower yield strength than the mean. Similarly, kurtosis rose steeply with α, indicating a leptokurtic distribution with a heavy tail and an increased risk of extreme failure modes. The stochastic perturbation approach provided reasonable estimates for expectation and variance at low $[eqn]$ but failed to capture the nonlinear growth in skewness and kurtosis. In contrast, both Monte Carlo simulations and semi-analytical methods were able to track these changes effectively, making them essential tools for the reliable prediction of yield behavior in the presence of structural uncertainty.

The effective Poisson’s ratio demonstrated a relatively muted response to increasing randomness compared to elastic modulus (seen in Figure 10 and Figure 11) and yield strength (seen in Figure 12 and Figure 13). Polynomial modeling of Poisson’s ratio under variability in cell length $[eqn]$ and wall thickness $[eqn]$ showed that the expectation remained nearly constant across the range of α as shown in Figure 14 and Figure 15 respectively, especially when cell length was the dominant uncertain parameter. This suggests that the global lateral-to-axial deformation behavior of the cellular material is less sensitive to small to moderate variations in microstructural geometry. The coefficient of variation increased slightly but linearly with α, implying a modest and predictable growth in dispersion. Skewness and kurtosis remained negligible for $[eqn]$ only starting to rise marginally at higher levels of randomness, particularly when wall thickness was perturbed. These higher-order effects suggest that Poisson’s ratio, though generally stable, may exhibit some stochastic sensitivity when cell geometry becomes highly irregular. Unlike the cases of Young’s modulus and yield strength, stochastic perturbation techniques remained sufficiently accurate for all four statistical measures across the full range of $[eqn]$ . This robustness, combined with the low sensitivity of Poisson’s ratio to geometric variability, makes it a relatively deterministic parameter in the mechanical characterization of cellular materials. Nevertheless, Monte Carlo and semi-analytical methods confirmed subtle nonlinear trends that may be critical in high-precision applications or when lateral compliance is a key design criterion.

The nonlinear response of the effective mechanical properties under uncertainty was modeled using sixth-order polynomial approximations derived from both analytical and numerical simulations. While this approach captures the essential trends in expectation, variance, skewness, and kurtosis of the output responses, the SPT exhibits limitations when applied beyond its rigorous domain of validity. In particular, SPT is mathematically reliable only for small perturbations, corresponding to low COV. In our study, the randomness parameter was varied up to COV = 0.275, and it was observed that SPT predictions begin to deviate at higher levels of input variability. Specifically, while SPT provides accurate estimates for the mean and variance at low to moderate COV (≤0.10–0.15), it systematically underpredicts skewness and kurtosis at larger perturbations, thereby failing to capture the emergence of heavy-tailed distributions and nonlinear amplification of fluctuations. MCS and SAM, by contrast, successfully reproduced these higher-order effects across the full range of variability. This comparative evidence demonstrates that SPT should be interpreted as reliable only in the low-COV regime while its application at higher variability levels is exploratory, with accuracy contingent upon benchmarking against more robust techniques such as MCS.

The entropy-based uncertainty quantification (as shown in Figure 16, Figure 17 and Figure 18) provides a rigorous framework for evaluating the divergence between analytical and numerical predictions of the effective mechanical properties of the hexagonal unit cell. Four relative entropy measures Kullback–Leibler, Jeffreys, Bhattacharyya, and Hellinger were employed to probe the sensitivity of Young’s modulus (see Figure 16), yield strength (see Figure 17), and Poisson’s ratio (cf. Figure 18) under perturbations in base material modulus $[eqn]$ , segmental length $[eqn]$ , and cell wall thickness $[eqn]$ . The results demonstrate that Young’s modulus is the most sensitive property, with Kullback–Leibler, Jeffreys, and Bhattacharyya measures exhibiting sharp increases in divergence as perturbation magnitudes increase, reflecting asymmetry and heavy-tail amplification in the probability distributions. By contrast, the Hellinger distance remains comparatively small, indicating a stronger emphasis on distributional overlap and thereby highlighting regions of robustness. Yield strength displays moderate sensitivity, primarily governed by variations in wall thickness, while Poisson’s ratio shows consistently low divergence across all measures, signifying its stability under parametric uncertainty. Importantly, each entropy metric provides complementary insight: Kullback–Leibler and Jeffreys are most effective in capturing asymmetry and parameter-driven instabilities, Bhattacharyya is suited to quantify overall separability between distributions, and Hellinger highlights zones of stability. From an engineering perspective, these divergences can be interpreted as indicators of performance reliability. Low divergence values (<0.05 across KL/Jeffreys, <0.02 in Bhattacharyya/Hellinger) suggest negligible degradation, moderate values (0.05–0.2 in KL/Jeffreys, 0.02–0.1 in Bhattacharyya/Hellinger) point to diminishing reliability margins, and high divergences (>0.2 in KL/Jeffreys, >0.1 in Bhattacharyya/Hellinger) signal critical thresholds where predictive fidelity is compromised. Safety factors or design adjustments become necessary. This comparative interpretation establishes the practical utility of the entropy measures, linking statistical divergences directly to performance degradation in mechanical properties.

Overall, the comparison reveals that Bhattacharyya, KL, and Jeffreys relative entropies are more sensitive to deviations in tail behavior and skewness, while Hellinger offers more conservative, overlap-based divergence measures. This entropy-based analysis confirms that the accuracy of homogenization under uncertainty is highly property- and parameter-dependent, with Young’s modulus being the most prone to divergence and Poisson’s ratio the most robust.

7. Concluding Remarks

In this study, we determined the mechanical properties of 2D hexagonal cellular solids-Young’s modulus, yield strength, and Poisson’s ratio under uncertainty using Monte Carlo simulation (MCS), stochastic perturbation technique (SPT), and the semi-analytical method (SAM). Comparative results between analytical and FEM-based homogenization approaches show strong agreement in predicting expected values, with FEM providing superior accuracy in capturing localized nonlinear effects. The findings demonstrate that material modulus strongly influences stiffness-related properties, whereas geometric variations primarily affect deformation metrics. The coefficient of variation and higher-order statistics highlight the growing impact of randomness, particularly on Young’s modulus and yield strength, underscoring the necessity of probabilistic modeling for reliability in design. Among these methods, SPT is reliable only in the low-randomness regime, whereas SAM and MCS are required to capture the full distributional behavior at higher variability levels. Overall, the framework developed here provides a validated basis for uncertainty-aware design of the lightweight and multifunctional materials, with direct implications for biomedical, aerospace, and structural applications at various industrial levels.