Scaling in Anti-Plane Elasticity on Random Shear Modulus Fields with Fractal and Hurst Effects

: The scale dependence of the effective anti-plane shear modulus response in microstructures with statistical ergodicity and spatial wide-sense stationarity is investigated. In particular, Cauchy and Dagum autocorrelation functions which can decouple the fractal and the Hurst effects are used to describe the random shear modulus ﬁelds. The resulting stochastic boundary value problems (BVPs) are set up in line with the Hill–Mandel condition of elastostatics for different sizes of statistical volume elements (SVEs). These BVPs are solved using a physics-based cellular automaton (CA) method that is applicable for anti-plane elasticity to study the scaling of SVEs towards a representative volume element (RVE). This progression from SVE to RVE is described through a scaling function, which is best approximated by the same form as the Cauchy and Dagum autocorrelation functions. The scaling function is obtained by ﬁtting the scaling data from simulations conducted over a large number of random ﬁeld realizations. The numerical simulation results show that the scaling function is strongly dependent on the fractal dimension D , the Hurst parameter H , and the mesoscale δ , and is weakly dependent on the autocorrelation function. Speciﬁcally, it is found that a larger D and a smaller H results in a higher rate of convergence towards an RVE with respect to δ .


Introduction
The continuum mechanical properties of any natural or engineering material are dependent on the material's physical structure at various scales below the macroscopic (or continuum). Thus, any continuum model corresponds to having effective (or averaged) properties of the heterogeneous substructure in some manner. The pertinent number of length scales and the relevant characteristics of substructure at each scale for obtaining these effective properties vary greatly from one material to another. One of the major challenges lies in homogenization, whose goal is to develop the phenomenological constitutive equations for the macroscopic response of idealized continua with heterogeneous substructures. The problem of homogenization for heterogeneous materials with linear mechanical and thermal properties has been studied extensively, with the conventional focus being on the effective (or macroscopic) responses of two or multi-phase random composites. This is the case for a representative volume element (RVE), as opposed to the statistical volume element (SVE), for which one is additionally interested in establishing the scaling to RVE [1,2]. Key papers in research on the SVE are [3,4].
In recent years, extensive literature has been produced on the computation of effective properties based on Monte-Carlo type simulations in various fields [5][6][7][8][9][10][11][12][13][14]. In the context of elasticity, Murshed et al. [7] obtained rigorous bounds on randomly distributed elastic polycrystals comprised of lower symmetry single crystals. Kale et al. [9] studied the scale-dependent bounds on the effective thermal conductivity of 2d Gaussian correlated microstructures. Savvas et al. [12] coupled an extended finite element method with Monte-Carlo simulations to determine effective stiffness and compliance elasticity tensors in random composites with differences in local volume fraction. Recently, Karimi et al. [14] investigated electrostatic and magnetostatic properties of 2d and 3d two-phase media after adapting Hill-Mandel conditions to electric and magnetic fields.
The present investigation focuses on the effects of fractal and Hurst characteristics of random media [15] on the scaling of their elastic properties when going from SVE to RVE. This is beyond all the previous studies that assumed either white-noise or very short-range correlations. Our study is focused on the anti-plane elasticity of media modeled by random fields of local shear modulus. Let the SVE in 2d be characterized by a mesoscale: a ratio of L = √ volume of SVE relative to the smallest heterogeneity size d. We ask: how do the random fluctuations in response to the SVE diminish to zero with mesoscale L increasing? This diminishing signifies attaining the effective deterministic response. This leads to stochastic partial differential equations under either Dirichlet or Neumann boundary conditions. The random fields are specified through Cauchy and Dagum autocorrelation functions which allow independent choices of fractal and Hurst characteristics. We employ a physics-based cellular automaton, a numerical method in which it is convenient to introduce randomness on a discretized domain, to solve the stochastic BVPs. From this, we obtain scale-dependent bounds on the effective shear modulus response at multiple length scales for various combinations of fractal and Hurst parameters.
This article is organized as follows. In Section 2, we first establish the governing equations, and discuss the Hill-Mandel conditions for elastostatics. We then introduce the concepts of apparent and effective properties and that of the scaling function. Next, we discuss the Cauchy and Dagum random fields that define the microstructure in this study, followed by an overview of the cellular automaton method that is applicable to anti-plane elasticity. Section 3 presents and discusses the homogenization results over a wide range of fractal and Hurst parameters for various SVE sizes. The major findings and conclusions are outlined in Section 4.

Governing Equations
Consider a random mesoscale body, B δ = {B δ (ω); ω ∈ Ω}, where B δ (ω) is a specific realization corresponding to the sample event ω in the Ω space. The dimensionless parameter δ = L/d characterizes the mesoscale, where L is the domain (window) size. Focusing on the quasi-static loading with negligible body forces, we have the following governing equation of local equilibrium: where σ is the the Cauchy stress field in the body, B δ (ω). Let the strain field in B δ (ω) be . The present study focuses on anti-plane linear elasticity for B δ ⊂ R 2 , so the equilibrium equation becomes: In general, C(x, ω) is the stiffness tensor random field (TRF) [16] and u ≡ u 3 (x, y) is the anti-plane displacement. Equation (2) reduces to a stochastic generalization of the Navier-Cauchy equation of elastostatics. We consider the TRF C(x, ω) = c(x, ω)I to be locally isotropic, where I is a second-rank identity tensor and c(x, ω) is henceforth referred to as the random shear modulus. Thus, the TRF model is modeled directly by a scalar random field. The central issue addressed here concerns the scale-dependence of C as it is smoothed on ever larger scales, while assuming the fractal and Hurst properties for the finest scale resolution.

Hill-Mandel Macrohomogeneity Condition
The transition from micro to mesoscale is studied in terms of the so-called Hill-Mandel macrohomogeneity condition [17,18]. This condition requires the equivalence of the energetic and mechanical interpretations of spatially-averaged work performed on a material domain. If these interpretations give the same result for different uniform boundary conditions imposed on the domain, then we have an RVE, and if the material domain is too small, we have an SVE. To study this mesoscale (L/d) dependence, consider σ and to be the volume averages of σ and over the SVE. The Hill-Mandel macrohomogeneity condition reads [1]: where t and n are traction and normal vectors, respectively. The notation a : b = tr(a T b) denotes the tensor scalar product between two second order tensors a and b. Equation (3) is satisfied by three types of uniform boundary conditions [2,7,19]: • Uniform traction (Neumann): • Uniform displacement-traction (mixed-orthogonal): where 0 and σ 0 are constant strain and stress tensors, respectively. Note that 0 = and σ 0 = σ on account of the average strain and stress theorems [2].

Apparent and Effective Properties
At any particular mesoscale δ, upon solving the stochastic boundary value problem under the boundary conditions (4) or (5), and consequent ensemble averaging, results in a different apparent stiffness (or compliance) tensor. In the current work, the uniform displacement and the uniform traction boundary conditions are considered. Equation (4) results in the apparent stiffness tensor C d δ (ω) , where d stands for displacement boundary conditions. On the other hand, Equation (5) results in the apparent compliance tensor S t δ (ω) , where t stands for traction boundary conditions and note that C t As δ → ∞, the randomness vanishes and the tensor represent the macroscale effective elastic response C eff of an RVE. In general, C d δ (ω) , which gives an upper estimate for the effective stiffness, is different from C t δ (ω) , which gives a lower estimate for the effective stiffness [10]. Indeed, upon using the principles of minimum potential energy and complementary energy with the assumptions of spatial homogeneity and ergodicity of the random microstructure, one can obtain the heirarchy of δ-dependent bounds on C eff as follows [3,4]: Upon using the definition of isotropic stiffness tensor along with Equation (9), we get the δ-dependent bounds on the shear modulus c eff as follows: where c R and c V are the Reuss (harmonic mean estimate of the shear modulus) and Voigt (arithmetic mean estimate of the shear modulus) estimates of the shear modulus.

Scaling
Consider the contraction between the apparent stiffness tensor and the apparent compliance tensor with the assumption of local isotropy for anti-plane elasticity in 2d: Following the previous studies, we postulate the following relation between the left-hand side of Equation (11) and the left-hand side of Equation (12): where g(δ) is the scaling function. Besides δ, it is also a function of the random microstructure. In general, g(δ) has two exact properties: (i) in the RVE limit: lim δ→∞ g(δ) = 0; (ii) for a homogeneous microstructure: g(δ) = 0. The next step is to introduce the random fields that specify the random medium.

Random Fields
Many natural material properties are usually described through fractal behavior and Hurst effects. We introduce the randomness into the shear modulus fields, c(x, ω), through two stochastic models, which can both capture and decouple fractal behavior and Hurst effects: Cauchy and Dagum [20][21][22][23]. The fractal dimension, D ∈ [n, n + 1), is a measure of roughness for a surface defined over R n . The larger is D, the rougher is the surface. The Hurst exponent or the Hurst parameter represents the "spatial memory." In the case of a 1d RF (i.e., parametrization by R), H ∈ [0, 0.5) represents an anti-persistent system, H ∈ (0.5, 1) represents a persistent system, and H = 0.5 corresponds to a true random walk. Both the Cauchy and Dagum shear modulus random fields are taken as Gaussian RFs, truncated to ensure the positive-definiteness of energy density.
The autocorrelation function for two random variables Z(x 1 ) and Z(x 2 ), R : We consider random fields, which are wide-sense-stationary (WSS) and isotropic. Define r = x 1 − x 2 , for the Cauchy class Gaussian random field, c(x, ω), x ∈ R 2 , the generalized Cauchy correlation function [20] is The special case with α = 2 corresponds to the well-known Cauchy model [22]. The generalized Cauchy random field will be hereafter referred to as Cauchy random field following the nomenclature in [23]. The correlation function corresponding to Dagum random field is given as follows [23]: The fractal dimension (D) and the Hurst parameter (H) are related to (n, α, β) when α ∈ (0, 2] and β ∈ (0, 2) as follows: where n = 2 for a 2d system. The special case when α = 2 − β corresponds to a self-affine fractal. The power spectral density of the Cauchy random field is given in [24], and that of the Dagum random field is given in [25]. In the current study, the Cauchy and Dagum random fields are generated for various combinations of α and β. We used the so-called circulant embedding method [26] to generate the random fields in R software.

Cellular Automata
The cellular automaton (CA) method is a computational approach first proposed by John von Neumann [27]. Leamy [28] studied inplane elastodynamics using a physics-based rectangular-celled CA in an idealized linear elastic medium. Later Hopman et al. [29] extended this method to study arbitrary-shaped two-dimensional geometries. The application of CA to the anti-plane shear elastodynamic problem is studied by Xian et al. [30]. The rectangular-celled CA yields discrete equations which are mathematically equivalent to that of the centered-difference finite difference method. In CA, the computational domain is discretized into cells and the material properties are assigned to each cell making it convenient to introduce randomness to the shear modulus field.
Following [28,30], the CA formulation for anti-plane shear forces at cell (i, j) is given below: where ∆x and ∆y are cell sizes in x and y directions, respectively; w is the width in the z direction; and c is the shear modulus value at cell (i, j). The balance of linear momentum equation for the cell (i, j) is given below: The cell (i, j) with its von Neumann neighbors and the corresponding forces acting on it are shown in Figure 1. Based on this method, we are concerned with solving BVPs under the uniform displacement and uniform traction boundary conditions (Equations (4) and (5)). The displacement values can be directly assigned to the boundary cells. The traction values can be assigned by applying equivalent forces at the centers of boundary cells.

Numerical Results and Discussion
The anti-plane elastostatic behavior resulting for various Cauchy and Dagum shear modulus random fields at successively increasing δ are reported in this section. The δ values considered are δ = 1, 5, 9, 31, 93, and 279. The combinations of α and β in Cauchy and Dagum random fields considered in this study are given in Table 1.  The apparent shear modulus responses, s t δ and c d δ , were obtained after applying CA for a large number of random field realizations corresponding to each (α, β, δ). The number of realizations was increased until the mean and the standard deviation of the apparent shear modulus responses did not change considerably with a further increase in the number of realizations. For all (α, β, δ), the ensemble average of all the realizations was set to 5, and the standard deviation was set to 1. After obtaining the apparent shear modulus responses, the scaling data were calculated for all (α, β, δ) using Equation (14). Then, a stretched exponential fit to the scaling data was employed following the studies conducted in [2,14]: Cauchy and Dagum-type functions, as defined below, are also postulated as possible fits to the scaling data.
All the fit coefficients in Equations (24)-(26) were obtained using the Curve Fitting Toolbox in MATLAB and are shown in Table 1, which uses the method of least squares.

Cauchy Random Field Responses
The mesoscale apparent shear modulus responses for the case of Cauchy random fields are shown in Figure 2. It can be observed that for all (α, β, δ), the inequalities of Equation (10)  . For the cases with α = 1 and varying β, the δ at which convergence occurred increased as the β value decreased. For the case with β = 1 and varying α, the δ values at which convergence occurred for the cases with α = 1 and α = 1.8 are similar and are far greater than that for the case with α = 0.2. Note that a higher α indicates a lower fractal dimension and a higher β indicates a lower Hurst parameter and vice versa. These results qualitatively suggest that the δ at which convergence occurs increases as the Hurst parameter increases. Similarly, the δ at which convergence occurs increases as the fractal dimension decreases. In order to quantitatively understand these results, the scaling data are plotted for all combinations of (α, β), along with the corresponding stretched exponential fits and the Cauchy function fits in Figures 3 and 4, respectively. The coefficients, A, B, and C, of the stretched exponential fits for the scaling data, are given in Table 1a. The higher B's value is, the quicker is the convergence, and the above-mentioned qualitative inferences are in line with the change in B value with (α, β). The coefficients of Cauchy function fits, A C , B C , and C C , are also given in Table 1a for reference. The residual error comparison between the stretched exponential fit and the Cauchy function fit shows that the Cauchy function fit is better than the stretched exponential fit in five out of the eight cases.

Dagum Random Field Responses
The mesoscale apparent shear modulus responses for the case of Dagum random fields are shown in Figure 5. It can be observed that for all (α, β, δ), the inequalities of Equation (10) are satisfied. Convergence occurs within the range of δ values considered for all (α, β). For the cases with α = 0.2 and varying β, the rate of convergence with respect to δ is higher for the case with a larger β value. For the case with β = 0.8 and varying α, the rate of convergence with respect to δ is higher for the case with a smaller α value. Analogously to the case of Cauchy random fields, the Dagum random field results qualitatively suggest that a larger β and a smaller α indicate a higher rate of convergence with respect to δ. The scaling data are plotted for all combinations of (α, β), along with the corresponding stretched exponential fits and the Dagum function fits in Figures 6 and 7, respectively. The coefficients, A, B, and C, of the stretched exponential fits for the scaling data, are given in Table 1b. The above-mentioned qualitative inferences are in line with the change in B value with (α, β). The coefficients of Dagum function fits, A D , B D , and C D are also given in Table 1b for reference. The residual error comparison between the stretched exponential fit and the Dagum function fit shows that the Dagum function fit is better than the stretched exponential fit in five out of the eight cases. Figure 2. Apparent shear modulus responses corresponding to Cauchy random fields for various combinations of (α, β, δ).  Five out of the eleven combinations of (α, β) are common in Cauchy and Dagum random fields considered in the current study (see Table 1). It can be observed that for each (α, β), the convergence is a bit quicker for the corresponding Cauchy random field response than that of the Dagum random field response, which is evident from the B values given in Table 1. However, note that for the same α, β, the differences between the Cauchy random field response and the Dagum random field response are not very different. These results suggest that the random shear modulus fields having a particular (D, H) may be modeled either by Cauchy or Dagum random fields without significant difference in the anti-plane response. Figure 5. Apparent shear modulus responses corresponding to Dagum random fields for various combinations of (α, β, δ).

Conclusions
In this article, the length scale effect on the apparent shear modulus values within the context dictated by the Hill-Mandel macrohomogeneity condition was studied. The effective shear modulus for a random microstructure was bounded above and below by the apparent responses of stochastic boundary value problems under uniform displacement and the uniform traction boundary conditions, respectively. The microstructural randomness was introduced through Cauchy and Dagum stochastic models which can both capture and separate fractal and Hurst effects. A varied set of (α, β) in both Cauchy and Dagum random fields were studied at each mesoscale δ. A physics-based cellular automaton was used to obtain the apparent shear modulus responses of the stochastic boundary value problems for all (α, β, δ). Major results include: (i) The obtained scaling data show that the mesoscale δ at which the effective shear modulus response is strongly dependent on the choice of parameters varies with α and β.
(ii) It was found that for a fixed α, the rate of convergence with respect to δ increases with an increase in β. For a fixed β, the rate of convergence with respect to δ increases with a decrease in α.
(iii) The scaling data of the corresponding Cauchy and Dagum random field responses were found to be approximated better by Cauchy-type and Dagum-type functions, respectively, than by a stretched exponential function which works very well in white-noise-type random media.
The methodology of the current work can be extended to study the scale-dependent bounds in media represented by tensor random fields without assuming local isotropy. Data Availability Statement: All data are shown clearly in the plots and tables, and will be made available up on reasonable request.