On the role of the microstructure in the deformation of porous solids

This study explores the role that the microstructure plays in determining the macroscopic static response of porous elastic continua and exposes the occurrence of position-dependent nonlocal effects that are strictly correlated to the configuration of the microstructure. Then, a nonlocal continuum theory based on variable-order fractional calculus is developed in order to accurately capture the complex spatially distributed nonlocal response. The remarkable potential of the fractional approach is illustrated by simulating the nonlinear thermoelastic response of porous beams. The performance, evaluated both in terms of accuracy and computational efficiency, is directly contrasted with high-fidelity finite element models that fully resolve the pores' geometry. Results indicate that the reduced-order representation of the porous microstructure, captured by the synthetic variable-order parameter, offers a robust and accurate representation of the multiscale material architecture that largely outperforms classical approaches based on the concept of average porosity.


Introduction
Porous materials have a long and distinguished history of applications spanning the most diverse fields of science and engineering [1][2][3]. In nature, the porous architecture is pervasive particularly in those applications that require the simultaneous achievement of low weight and high stiffness; a concept known as high specific stiffness. The most notable and broadly available example of such material is found in animal bones, although other compelling examples of naturally occurring porous materials include wood, marine shells, and rocks [4][5][6]. The pursuit of materials endowed with extreme specific stiffness has long been a major quest in structural engineering. Learning from nature, it would appear logical to make more extensive use of porous architectures, hence providing transformational alternatives to more traditional materials like alloys and composites. While the design of natureinspired engineering materials is not new in the scientific community and has seen many successful applications [7], the use of porous materials within the design cycle has not progressed accordingly. The main reasons contributing to this lag are rooted in either technical or theoretical aspects. At a technical level, porous materials are extremely challenging to fabricate and the few realizations typically rely on a random microstructure (e.g. metal foams [8]). At a theoretical level, the incomplete understanding of how size effects emerge from the underlying microstructure and affect the intrinsic multiscale response of the material limits significantly the ability to design and simulate their performance. The recent rapid advances in additive manufacturing have markedly eased the technical concerns, hence suggesting that the design of engineered structural materials with carefully crafted porosity could be finally in reach.
On the contrary, the theoretical gap is still to be vanquished, hence severely hindering the use of porosity as a design concept. Looking more closely at the theoretical aspects, several theoretical and experimental investigations have shown the central role that size effects play in determining the response of porous structures as well as nontrivial phenomena such as softening behavior [9,10], and power-law scaling of structural constitutive relations (e.g. dispersion relations, effective structural strength, and average permeability [4,6,[11][12][13]). However, existing models do not accurately capture the size effects because they adopt homogenized microstructural descriptions that oversimplify the effect of the porous geometry by capturing only the degree of porosity, while ignoring the size and distribution of the pores [10,14]. From a practical perspective, the homogenization process discards many physical attributes of the porous microstructure when modeling its size effects within the continuum. Indeed, existing models either fully ignore size effects (e.g. rule of mixture models [15,16]) or assume an ad hoc spatially-uniform distribution of the size effects (e.g. classical nonlocal and micro-mechanical models [17]). We will show that, at the continuum level, this latter assumption leads to an inaccurate physical representation of the effect of the porous microstructure. From a more general perspective, and independently of the specific approach, existing modeling strategies share a common set of critical limitations including: [L1] the demand for computationally prohibitive resources, when attempting to fully resolve the pores' geometry in pursue of high-fidelity approaches; and [L2] the lack of accuracy in numerical predictions, typically driven by the inability to account for the effect of the pores' configuration and their resulting impact on size effects. The ability to overcome these critical limitations would unlock the many opportunities offered by the broad class of porous engineered materials, hence enabling a fundamental leap forward in a variety of pivotal applications including, but not limited to, light-weight components for aerospace, terrestrial, and marine transportation [18,19], batteries and energy storage systems [20,21], and even robotics and biomedical implants [3].
In this study, we attempt to bridge the fundamental theoretical gaps highlighted by the above limitations L1 and L2. We start by embracing the intrinsically nonlocal nature of porous solids by directly linking the size effects to the specific characteristics of the porous microstructure. This approach provides a solid physical foundation to develop highly accurate and computationally efficient nonlocal continuum models and it opens the way to new generation theoretical models based on variable-order (VO) fractional calculus 1 . Notably, the use of VO operators is not a mere mathematical expedient to more conveniently model these complex media, but instead it is a natural consequence of the nature of the underlying microstructure. In other terms, we will show that the fractional-order mechanics emerges naturally from the physical characteristics of porous media so that the resulting model can accurately resolve the interplay between the porous geometry, the size effects, and the experimentally observed power-law characteristics [28,29].

Results and Discussion
The physics of deformation A porous medium can be approximated as a collection of hollow spherical inclusions dispersed within a bulk solid. While the pores are not necessarily spherical in nature, the total volume of a pore can be effectively segmented into a series of spheres using standard concepts of segmentation [30] [ Fig. 1(a)]. It follows that a porous structure can be modeled as a network of hollow spherical inclusions dispersed within a solid matrix. Henceforth, we refer to the 'hollow spherical inclusions' as 'spherical inclusions' for the sake of brevity. Further, the spherical inclusions interact with each other via long-range forces [ Fig. 1(a)], hence forming a network of nonlocal interacting particles in the bulk solid. In order to visualize this concept, consider a porous solid with only two inclusions. Any change in the shape of one inclusion (e.g. due to the application of external mechanical forces) will affect also the shape and response of the other inclusion. This example is a direct realization of the principle of action-at-a-distance which is at the foundation of atomistic, molecular, and nonlocal continuum theories [31][32][33][34][35]. Note that the magnitude of the nonlocal force between two inclusions decays monotonically with the increasing (relative) distance following a power-law distribution [36]. This behavior is closely related to the decay of the magnitude of the stress field, which decreases as a power-law function of the radial distance from the spherical inclusion [37]. Indeed, the presence of a second inclusion within the horizon of influence (that is the region affected by the stress concentration due to the spherical inclusion) will deform under the stress field of the first inclusion and produce, in turn, a stress field around itself that will affect the next inclusion [see Fig. 1(c)]. It is readily seen how, by following this cascading process, the deformation of the first inclusion can affect the third inclusion (and beyond), hence producing an evident nonlocal effect. Further, the strength of the interaction between inclusions is directly dependent on the concentration of inclusions in a target area as well as on their relative size [38]. More specifically, the presence of either a larger number of inclusions or of inclusions of larger size will enhance the stress concentration and the corresponding nonlocal effect. It follows that spatial variations in size and distribution of the inclusions will result in variations in the strength of the nonlocal interactions within the network of inclusions. From the above discussion, it clearly emerges that porous media exhibit nonlocal effects characterized by a position-dependent strength strongly influenced by the porous microstructure.

The VO nonlocal elasticity model
Building on the previous observations of the physical mechanisms taking place in a porous medium, we derive a nonlocal continuum theory capable of capturing the complex nonlocal interactions of the solid. Consider a 1D infinite lattice of inclusions with spatial period ∆x, as shown in Fig. 1(b) [39]. The location and displacement of the n th inclusion (where n ∈ Z) are denoted as x n and u n , respectively. The nonlocal interactions between the inclusions are modeled as lumped springs. The lumped springs capture both the exchange and the redistribution of elastic energy between the inclusions, and the analogous elastic interactions are indicated as dashed lines within the bulk in Fig. 1(a). The power-law decay in the interaction force between two inclusions i and j( = i) is embedded in the stiffness term characterizing the lumped spring k ij : where α(x i ) captures the position-dependent strength of the nonlocal interactions. To obtain a well-posed and causal theory, it is required that α(x i ) ∈ (0, 1] [29]. Notably, the continuum model with the power-law distribution of order 1 + α(x i ) ∈ (1, 2] directly reflects the power-law dispersion relations [40], which is widely supported by experimental observations [4,11,29]. Additionally, we also note that strength of long-range forces between spherical inclusions also exhibit power-law distributions with exponents in the interval (1, 2) [36,38]. c 0 is a book keeping parameter defined in Equation (3). The total internal force on the i th inclusion is given by: To facilitate continualization (that is, ∆x → 0), ψ 0 is defined as ψ 0 = EAl where E and A denote the Young's modulus and cross-sectional area of the equivalent 1D continuum, respectively. l * is a length-scale parameter that ensures dimensional consistency of the formulation [analogous to l ∓ in Equation (7)]. Note that, ψ 0 = EA∆x for α(x i ) = 1, indicating that it captures the energy stored in an interaction between adjacent inclusions in a local lattice. Finally, the stress in the 1D nonlocal solid is derived from the continuum limit of the internal force as (see Supplementary Note 1): where Γ(·) denotes the Gamma function, the stress can be expressed in a compact form via VO fractional derivatives: where ε denotes the strain in the 1D nonlocal solid, and ∞ u denote the VO left-handed and right-handed Caputo derivatives, respectively. D α(x) x u denotes the VO Riesz-Caputo derivative. The above analysis suggests that VO operators are naturally equipped to model porous structures because the differ-integral nature of the VO fractional operator allows capturing the nonlocal interactions, while the spatially varying exponent of the power-law kernel of the operator (that is, α(x)) captures the position-dependent strength of the interactions.
In order to leverage the potential of VO operators to model porous structures of practical interest, we develop a 3D VO nonlocal continuum theory by directly extending the 1D strain relations obtained in Equation (4): where u i denotes the displacement field alongx i direction at a given point x. Next, the stress tensor is derived from thermodynamic balance laws as: where λ and µ denote the Lamé constants, α θ is the coefficient of thermal expansion, and θ denotes the temperature difference between the solid and the environment. The detailed derivation of the strain and stress tensors can be found in Supplementary Note 2 and Note 5, respectively.
The VO Riesz-Caputo derivative D α(x) x j u i in Equation (5) is defined on the interval [x − j , x + j ] as: l − j and l + j are length scales along thex j direction and ensure both the dimensional consistency and the frameinvariance of the model. Frame-invariance also requires that l − j = x j − x − j and l + j = x + j − x j . Physically, the length scales determine the dimensions of the horizon of nonlocality on either side of a point, that is the distance beyond which nonlocal forces cease to act. Detailed discussions on the consistent behavior of the model at material boundaries and interfaces, and the frame-invariance of the model can be found in Supplementary Note 3 and Note 4, respectively.

Overview of the numerical experiments
In order to validate and demonstrate the advantages of the VO theory when applied to the modeling of porous structures, we simulate and analyze three configurations of porous beams with porosity p 0 ∈ {0.20, 0.22, 0.24} [ Fig. 2(a)]. All beams have length L = 1m, width b = 0.02m, height h = 0.02m, and length scales l − = l + = 0.2L [14]. The bulk solid has Young's modulus E 0 = 30MPa, Poisson's ratio ν = 1/3, and α θ = 2.3 × 10 −5 K −1 . Each porous core is enclosed in a rectangular solid outer shell to simplify the application of external loads and boundary conditions [ Fig. 2(b)]. The porous beams were simulated using the nonlinear VO beam model, that was obtained from the 3D theory by employing Euler-Bernoulli assumptions (see Supplementary Note 6). Note that, to apply the VO theory, it is necessary to first determine the VO α(x) that characterizes the underlying microstructural configuration and the resulting nonlocal interactions. In order to determine α(x) [Fig. 2(c)], we formulated an inverse problem by means of deep learning techniques [41]. The details concerning the network architecture and training phase are provided in Supplementary Note 7. With the VO distribution available, the models can be solved numerically using the fractional-order finite element method [42], to validate the VO model and assess its performance with respect to the previously identified limitations L1 and L2. Result sets R# ([R1] Analysis of computational performance, and [R2] Pore configuration, nonlocal effects, and role of α(x)) map directly to the limitations L#. The numerical accuracy of the VO model can be assessed by considering the nonlinear static response of the porous beams. All configurations are assumed clamped at both ends and subject to a uniformly distributed transverse load (UDTL) of magnitude P 0 . To provide meaningful comparisons with established methods widely adopted in literature, the beams' response was also simulated by using the 3D finite element (FE) approach and the rule of mixtures approach [15,16]. The 3D FE approach, implemented in this study via the commercial FE software COMSOL Multiphysics, is chosen because it can fully resolve the porous geometry of the beam and hence it provides a reliable reference solution for performance assessment. The rule of mixtures approach is an integer-order (IO) homogenization approach that converts the initial porous beam into an anisotropic (non-porous) beam with a spatially varying modulus of elasticity E(x) = [(1 − p 0 (x)]E 0 , where p 0 (x) denotes the average porosity of planes perpendicular to the mid-plane of the beam (see Supplementary Note 8). The IO approach allows highlighting the remarkable effect on accuracy that VO mechanics offers over classical IO mechanics.
The nonlinear static response of the porous beams are presented in Fig. 3(a) in terms of the average transverse displacement ( w 0 ) and the transverse displacement of the geometric center [x(L/2, 0, 0); see Fig. S2] of the beam (w 0 ). Note that w 0 presents a global comparison of the three approaches, while w 0 allows a point-wise comparison. To further explore the point-wise comparison, the deformed shapes of the mid-planes for P 0 = 4N are presented in Fig. 3(b). Additionally, for a more complete validation, we also simulated the thermoelastic response of the porous beams when subjected to a steady state uniform temperature θ = 100K, in addition to the UDTL P 0 . The results of the thermoelastic study are presented in Fig. 4.

Analysis of computational performance
The total number of degrees of freedom and the average run time (for three repeated simulations) are presented in Fig. 5 for each method. The number of degrees of freedom was chosen so to guarantee proper resolution of the displacement profile. This criterion was defined by setting an arbitrary threshold of 1% on the magnitude of the difference of the displacement field between successive spatial refinements. The VO and IO models were simulated on a personal computer equipped with an Intel(R) Core(TM) i7-1065G7 processor and 8 GB RAM, and required less than 2 minutes run time. The FE models were run on a cluster with 2×AMD EPYC 7662 processors with 20  cores and 256 GB RAM, and required at least 10 hours run time 2 . While the two platforms are different, the cluster is significantly more powerful than the personal computer, hence making the comparison of the computational performance even more conservative. It follows that, based on the above results, the VO approach delivers an unparalleled computational efficiency (compared with the 3D FE approach) while simultaneously maintaining same level of accuracy. In this regard, note that the VO approach can achieve superior accuracy when compared to the rule of mixtures approach. This should not be surprising since the VO model fully captures the position-dependent nonlocal interactions between the pores (as we will demonstrate further in R2), while the the rule of mixtures model is a local approach.

Pore configuration, nonlocal e ects, and role of α(x)
The different results (in R1) presented for the numerical validation and analysis of the performance, suggest that the VO approach provides a physically meaningful parameter (that is, the VO α(x)) which accurately captures the characteristics of the underlying porous microstructure, hence enabling an efficient continuum level formulation endowed with detailed microscale information. Indeed, as motivated via the lattice model, the VO α(x) directly links the pores' configuration to the strength of the nonlocal effects. This latter concept is further supported, here below, by the results of the VO model that consistently captures the impact of the nonlocal effects on the anomalous softening behavior of the porous beams.
At the macro scales, the most direct impact of the presence of nonlocal effects is observed as structural softening, wherein the strength of the nonlocal effects determine the degree of softening [40]. In this regard, for porous beams, the position-dependent strength of the nonlocal interactions (due to the spatially-varying porous microstructure) is expected to result in a spatial variation in the degree of softening, when compared to an isotropic non-porous beam (p 0 = 0). Indeed, the VO model successfully captures the position-dependent softening of the porous beam, as evident from the variation of S = w 0 (p 0 = 0)/w 0 (p 0 = 0) across the length of the beam [ Fig. 5(b)]. The parameter S captures the increase in displacement at a given point on the porous beam, when compared to an isotropic non-porous beam, on account of the nonlocal interactions resulting from the porous microstructure. This result establishes clearly that the simplifying assumption of spatially-uniform size-effects in classical nonlocal and micro-mechanical models introduces severe inaccuracies.
Further, an analysis of the maximum transverse displacement (w 0 ) [ Fig. 5(c)] shows that an increase in the average porosity, for a fixed P 0 , does not result in a consistent softening of the porous beam (that is, an increase inw 0 ). This observation is in contrast to the expected behavior based on a classical IO model. However, an increase in the strength of nonlocality (e.g. via a decrease in the average VO [α 0 = α(x) ; α 0 for beams with p 0 = {0, 0.20, 0.22, 0.24} are α 0 = {1, 0.88, 0.84, 0.86}]) correlates in a consistent manner with the softening phenomenon 3 . Thus, it emerges that at the continuum level the average nonlocal strength α 0 , and not the average porosity p 0 , provides a much more appropriate and physically consistent representation of the porous microstructure. On a more theoretical note, we emphasize that the softening behavior observed in the porous beams (in both the nonlinear elastic and thermoelastic response) was captured consistently (that is, without any display of paradoxical effects [17,43]) by the VO theory. This result is possible thanks to the underlying foundation provided by the displacement-driven formulation [based on differ-integral kinematic relations, see Equation (5)], which is intrinsically well-posed and thermodynamically consistent [40].
For completeness, we also highlight that the presence of nonlocal effects does not always lead to softening behavior at the continuum level. Particularly, in the case of nonlocal solids with multiscale effects, an increase in the strength of nonlocality could potentially lead to either stiffening or softening behavior [44,45]. In this case, a different approach should be used, for example, distributed-order fractional models [45] or fractional-order strain gradients in the nonlocal constitutive formulation [44]. Nevertheless, in the present study, the radii of the pores distributed within the bulk are within the same scale (or, order) when compared to the scale of the beam and hence, we do not account for multiscale effects. As also verified numerically (based on the fully resolved geometry via 3D finite element analysis), the porous beams exhibit a consistent softening behavior indicating the absence of any multiscale effects.

Outlook
In summary, this study has provided two key contributions to the physical understanding and modeling of porous media. First, it was shown that porous solids exhibit position-dependent nonlocal effects that cannot be neglected if accurate predictions are sought. Second, a VO nonlocal continuum theory capable of capturing these complex nonlocal effects was developed and numerically tested. Results highlighted that VO mechanics offers a uniquely powerful approach to develop efficient continuum models endowed with fine level details capturing the underlying microstructure. The VO approach provides a rigorous methodology to develop physically-consistent reduced-order models of multiscale systems with an accuracy comparable with fully resolved 3D models. While the results were presented in the context of porous materials, it is expected that the VO framework could be extended to a variety of applications characterized by multiscale features including, but not limited to, composites, architectured materials, seismology, biotechnology, and much more.

Methods
The key methods adopted in deriving the results include: [M1] the use of VO fractional calculus to develop the VO nonlocal elasticity model presented in the results; and [M2] the use of a deep learning method to obtain the VO laws that characterize the porous solids. These methods have been described in detail in the supplementary information document. The method M1 has been divided into specific analytical derivations which are presented in the Supplementary Notes 1 -6 and summarized here below: • Supplementary Note 1 presents the analytical continualization of the discrete lattice to obtain the stress and strain fields in the 1D continuum, that is essentially the derivation of Equation (2) from Equation (4).
• Supplementary Note 2 presents the derivation of the 3D VO nonlocal strain tensor.
• Supplementary Note 3 presents the analytical proof of the consistent behavior of the VO formulation at material boundaries and interfaces.
• Supplementary Note 4 presents the frame-invariance of the VO nonlocal continuum model.
• Supplementary Note 5 presents the derivation of the stress tensor and the proof of the thermodynamic consistency of the VO nonlocal continuum model.
• Supplementary Note 6 presents the derivation of a geometrically nonlinear model for nonlocal slender beams from the 3D VO nonlocal continuum model.
• Supplementary Note 7 presents the deep learning framework for the identification of the VO that characterizes the response of different nonlocal beams.
The method M2 is presented in Supplementary Note 7.

Data Availability
All the necessary data and information required to reproduce the results are available in the paper. Additional data in the form of equations supporting this article are provided in the supplementary information.

Code Availability
The files of the numerical simulation data corresponding to the results in the paper will also be shared. The actual physical code may be subject to restrictions from the sponsors.