Application of generalized boundary conditions for homogenization of thermal and filtration properties of soils

Random composites, including porous, multi-phase media, such as soils, are often characterized by a high variability of the filtration and thermal properties of the components. Solving problems of mechanics for these media on a macroscopic scale is therefore associated with the necessity of volumetric averaging of these properties. In numerical analysis, the problem of averaging (homogenization) is based on solving the appropriate boundary value problem (BVP) on a series of statistical volume elements (SVEs) using the finite element method. Special boundary conditions (BCs) are utilized which ensure that the requirement of energy equivalence between the macro and micro scales is fulfilled. In the case of deformation analysis, the uniform, kinematic or static, and periodic boundary conditions are used [3,4,5,6,7,8,9,10, 13, 15, 16]. However, these BCs have disadvantages: uniform ones generate only limit values of effective properties (upper limit for kinematic conditions and lower limit for static conditions) and periodic conditions require periodic SVE. In the case of random media, such as soil, ensuring the periodicity of the solution, therefore, means altering the microstructure of the material.

Recently, the so-called generalized BCs have been proposed. They allow for determination of the apparent composite properties, which reside between the solutions provided by the uniform BCs, and, at the same time, the periodicity of SVE is not required. In the work [22], the usefulness of these conditions for homogenization of the elastic properties of random matrix-inclusion composites and the pixelized 2D media was analyzed. In the present paper the same idea is used for coefficients of the Laplace's partial differential equation (PDE). General equations of the method and also the finite element formulation are provided. Utilization of the generalized boundary conditions is then demonstrated on two 3D examples: homogenization of the filtration coefficient and thermal conductivity of the soil. In this paper, only the microscopic BVP is considered, since the boundary conditions applied to the statistical volume elements are investigated. The macroscopic values are visible here only as the loading quantities and the averages over SVEs.

Let us consider anisotropic Laplace equation defined on SVE ω: (1) −(aijp,j),i=0, - {({a_{ij}}{p_{,j}})_{,i}} = 0, where a_ij is a tensor function of parameters (positive definite) and p is an unknown function. In the scope of computational homogenization, BCs for this PDE are determined by the averaging equation: (2) ∫ωp,i dω=∫∂ωpni ds=VP,I, \int\limits_\omega {{p_{,i}}\;{\rm{d}}\omega } = \int\limits_{\partial \omega } {p{n_i}\;{\rm{d}}s = V{P_{,I}},} where P_,I is a known quantity, which usually results from the upper scale computations, V is the SVE volume and ∂ω is the SVE boundary. Weak form of equation (1) reads as (3) −∫ω(aijp,j),i q dω=0, - \int\limits_\omega {{{({a_{ij}}{p_{,j}})}_{,i}}\;q\;{\rm{d}}\omega } = 0, with q being an arbitrary testing function belonging to the same function space as p. Applying Gauss-Ostrogradsky theorem the following is obtained: (4) ∫ωaijp,jq,i dω−∫∂ωrq ds=0, \int\limits_\omega {{a_{ij}}{p_{,j}}{q_{,i}}\;{\rm{d}}\omega } - \int\limits_{\partial \omega } {rq\;{\rm{d}}s} = 0, where r is interpreted as the boundary flux. In the spirit of generalized BCs [22], the following is assumed at boundary: (5) r=−κp˜, r = - \kappa \tilde p, where κ is called the boundary transfer function and (6) p˜=p−P,Ixi \tilde p = p - {P_{,I}}{x_i} is the perturbation of p at boundary. Function κ has to be guessed since it depends on the properties of SVE and its surrounding (see section 4). In addition, the averaging equation (2) is enforced by means of the Lagrange multipliers term. The considered weak form is then given as (7) ∫ωaijp,jq,i dω+∫∂ωκp˜q ds+ΛI(∫∂ωqni ds−VP,I)=0, \int\limits_\omega {{a_{ij}}{p_{,j}}{q_{,i}}\;{\rm{d}}\omega } + \int\limits_{\partial \omega } {\kappa \tilde pq\;{\rm{d}}s} + {\Lambda _I}\left( {\int\limits_{\partial \omega } {q{n_i}\;{\rm{d}}s - V{P_{,I}}} } \right) = 0, where Λ_I are the mentioned multipliers which are real numbers.

Solution of problem (7) is equivalent to the unconstrained minimization: (8) minp,ΛI Π(p, ΛI), \underset{p,{{\Lambda }_{I}}}{\mathop{\text{min}}}\,\text{ }\!\!~\!\!\text{ }\Pi \left( p,~{{\Lambda }_{I}} \right), where potential Π is given by (9) Π(p, ΛI)=12∫ωaijp,jp,i dω+ ∫∂ωκp˜p ds+ΛI(∫∂ωpni ds−VP,I). \begin{array}{*{35}{l}} \Pi \left( p,~{{\Lambda }_{I}} \right) & =\frac{1}{2}\int\limits_{\omega }{{{a}_{ij}}{{p}_{,j}}{{p}_{,i}}\ \text{d}\omega } \\ {} & +\ \int\limits_{\partial \omega }{\kappa \tilde{p}p\ \text{d}s}+{{\Lambda }_{I}}\left( \int\limits_{\partial \omega }{p{{n}_{i}}\ \text{d}s-V{{P}_{,I}}} \right). \\ \end{array} Minimum is found by equating variations of Π in p and Λ_I to 0, that is, (10) δΠ(p, ΛI, δp)=∫ωaijp,jδp,i dω +∫∂ωκp˜δp ds+ ΛI∫∂ωδpni ds=0, \begin{array}{*{35}{l}} \delta \Pi \left( p,~{{\Lambda }_{I}},~\delta p \right) & =\int\limits_{\omega }{{{a}_{ij}}{{p}_{,j}}\delta {{p}_{,i}}\ \text{d}\omega }\ +\int\limits_{\partial \omega }{\kappa \tilde{p}\delta p\ \text{d}s} \\ {} & +\ {{\Lambda }_{I}}\int\limits_{\partial \omega }{\delta p{{n}_{i}}\ \text{d}s=0}, \\ \end{array} (11) δΠ(p, ΛI, δΛI)=δΛI(∫∂ωpni ds−VP,I)=0. \delta \Pi \left( p,~{{\Lambda }_{I}},~\delta {{\Lambda }_{I}} \right)=\delta {{\Lambda }_{I}}\left( \int\limits_{\partial \omega }{p{{n}_{i}}\ \text{d}s-V{{P}_{,I}}} \right)=0. The above equations must be fulfilled for any admissible variations δp and δΛ_I. So, let us consider δp = x_i1_i, and thus, δp_,i = 1_i. From the first equation, the following is obtained: (12) −ΛI=1V∫ωaijp,j dω+1V∫∂ωκp˜xi ds. - {\Lambda _I} = \frac{1}{V}\int\limits_\omega {{a_{ij}}{p_{,j}}\;{\rm{d}}\omega } + \frac{1}{V}\int\limits_{\partial \omega } {\kappa \tilde p{x_i}\;{\rm{d}}s} . Recognizing that (13) aijp,j=fi {a_{ij}}{p_{,j}} = {f_i} is the flux quantity defined on SVE, then the first component in the right-hand side of equation (12) is its volumetric average, interpreted as the macroscopic flux F_I, and: (14) −ΛI=FI+1V∫∂ωκp˜xi ds. - {\Lambda _I} = {F_I} + \frac{1}{V}\int\limits_{\partial \omega } {\kappa \tilde p{x_i}\;{\rm{d}}s} . Thus, Lagrange multipliers stand for the macroscopic flux (with minus sign) modified by the component related to the boundary perturbation p˜ \tilde p . Usually this second component is enforced to vanish by applying p˜=0 \tilde p = 0 or by assuming its periodicity on ∂ω. However, this is not the case in the current framework of generalized BCs, where perturbation is allowed, but controlled by κ transfer coefficient.

Let us also consider δp = p, and thus, δp_,i = p_,i. Again, from equation (10), the following is obtained: (15) ∫ωaijp,jp,i dω+∫∂ωκp˜p ds+ΛI∫∂ωpni ds=0. \int\limits_\omega {{a_{ij}}{p_{,j}}{p_{,i}}\;{\rm{d}}\omega } + \int\limits_{\partial \omega } {\kappa \tilde pp\;{\rm{d}}s} + {\Lambda _I}\int\limits_{\partial \omega } {p{n_i}\;{\rm{d}}s = 0} . Inserting (2), (13), and (14) into (15), after some transformations, the following is derived: (16) FIP,I=1V∫ωfip,i dω+1V∫∂ωκp˜p˜ ds. {F_I}{P_{,I}} = \frac{1}{V}\int\limits_\omega {{f_i}{p_{,i}}\;{\rm{d}}\omega } + \frac{1}{V}\int\limits_{\partial \omega } {\kappa \tilde p\tilde p\;{\rm{d}}s} . This formula is identified as the Hill-Mandel macrohomogeneity principle, which relates macroscopic energy (left-hand side) to the average of the microscopic energy (right-hand side). Again, the second component on the right will vanish with vanishing or periodic perturbation p˜ \tilde p .

It has been shown [22] that solution of the system of equations (10) and (11) is unique for the given SVE and for the specific boundary transfer function κ chosen. This function is the key ingredient of the generalized BCs. In general, the only requirement is that it must be positively real-valued, that is, κ : ∂ω → ℝ⁺. Specifically, if the values of κ tend to 0⁺ all over the boundary ∂ω, then the solution becomes governed by the averaging, Lagrange multipliers term. This results in (almost) unconstrained perturbation p˜ \tilde p , so this is equivalent in applying the so-called static or minimal kinematic BCs (see, e.g., [20]). On the other hand, if κ → ∞, then it acts as the penalty parameter and the solution is such that perturbation p˜ \tilde p vanishes. This is equivalent in applying the kinematic uniform BCs, where p = P_,Ix_i is enforced on ∂ω. If the macroscopic fluxes obtained with these two κ selections and by application of the macroscopic gradient P_,I are denoted as FI0 F_I^0 and FI∞ F_I^\infty , then the following holds: (17) FI0P,I≤FIPI≤FI∞P,I, F_I^0{P_{,I}} \le {F_I}{P_I} \le F_I^\infty {P_{,I}}, where F_I is the flux obtained with any other admissible κ choice. In other words, solutions of the considered homogenization problem are properly bounded by the limiting cases.

In this paper, by analogy with the approach presented in [22] function κ is taken as: (18) κ=1ℓaijnjni, \kappa = \frac{1}{\ell }{a_{ij}}{n_j}{n_i}, where ℓ is interpreted as the microscopic length parameter, a_ij describes the material crossing the boundary, and n_i is a unitary vector, normal to ∂ω. Basic assumption which stands behind this choice is that at the boundary, the following holds: (19) p˜nj=p˜,jℓ. \tilde p{n_j} = {\tilde p_{,j}}\ell . If the SVE neighbourhood is taken exactly the same as at the boundary, then this assumption leads to the following external loading term: (20) r=−aijp˜,jni=−1ℓaijnjnip˜=−κp˜. r = - {a_{ij}}{\tilde p_{,j}}{n_i} = - \frac{1}{\ell }{a_{ij}}{n_j}{n_i}\tilde p = - \kappa \tilde p.

As a final remark, it should be pointed out that the selected κ depends on the boundary orientation, and thus, it is variable even if a_ij is constant all over the boundary. However, for the case of isotropic materials, where a_ij = aδ_ij, with a being a constant scalar function and δ_ij being Kroenecker delta, the orientation becomes unimportant and κ=aℓ \kappa = \frac{a}{\ell } .

Solutions p of the problem defined by equations (9), (10), and (11) have to belong to the Sobolev functional space H¹(ω) of differentiable functions (at least once). This is ensured by the piece-wise polynomial base functions used in standard FEM discretization. Rearranging equation (9) to the form (21) Π(p, ΛI)=12∫ωaijp,jp,i dω+∫∂ωκpp ds+−P,I∫∂ωκxip ds+ΛI(∫∂ωpni ds−VP,I) \begin{array}{*{35}{l}} \Pi \left( p,~{{\Lambda }_{I}} \right)= & \frac{1}{2}\int\limits_{\omega }{{{a}_{ij}}{{p}_{,j}}{{p}_{,i}}\ \text{d}\omega }+\int\limits_{\partial \omega }{\kappa pp\ \text{d}s}+ \\ {} & -\,{{P}_{,I}}\int\limits_{\partial \omega }{\kappa {{x}_{i}}p\ \text{d}s}+{{\Lambda }_{I}}\left( \int\limits_{\partial \omega }{p{{n}_{i\ }}\text{d}s-V{{P}_{,I}}} \right) \\ \end{array} and skipping details of the finite element approximation (see [20, 22] for that) one ends up immediately in discreet version of this potential: (22) Π^(p, Λ)=12pT(K+D)p−pTGP+Λ(pTB−PTV), \hat{\Pi }\left( \mathbf{p},~\mathbf{\Lambda } \right)=\frac{1}{2}{{\mathbf{p}}^{T}}\left( \mathbf{K}+\mathbf{D} \right)\mathbf{p}-{{\mathbf{p}}^{T}}\mathbf{G}{\boldsymbol {P}}+\mathbf{\Lambda }\left( {{\mathbf{p}}^{T}}\mathbf{B}-{{\mathbf{P}}^{T}}V \right), where p is global vector of unknowns of length M (M – total number of nodes in discretization), K is a global linear operator of size M × M, D is the boundary constraint matrix, also of size M × M, G is a boundary loading matrix of size M × D (D – space dimension: 2 or 3), B is a matrix resulting from the application of averaging constraint, also of size M × D, Λ is a vector of unknown Lagrange multipliers of length D, and P is a vector of known macroscopic gradient (also of length D). Minimum of potential Π^ \hat \Pi is found by equating its variations in p and Λ to 0, which is denoted as (23) ∂Π^(p,Λ)∂p=(K+D)p+BΛ−GP=0, \frac{{\partial \hat \Pi \left( {{\bf{p}},{\bf{\Lambda }}} \right)}}{{\partial {\bf{p}}}} = \left( {{\bf{K}} + {\bf{D}}} \right){\bf{p}} + {\bf{B\Lambda }} - {\bf{G}}{\boldsymbol P} = {\bf{0}}, (24) ∂Π^(u,Λ)∂Λ=BTp−PV=0. \frac{{\partial \hat \Pi \left( {{\bf{u}},{\bf{\Lambda }}} \right)}}{{\partial {\bf{\Lambda }}}} = {{\bf{B}}^T}{\bf{p}} - {\boldsymbol {P}}V = {\bf{0}}. This is a system of linear equations (25) [K+DBBT0][pΛ]=[GPPV], \left[ {\matrix{{{\bf{K}} + {\bf{D}}} & {\bf{B}} \cr {{{\bf{B}}^T}} & {\bf{0}} \cr } } \right]\left[ {\matrix{{\bf{p}} \cr {\bf{\Lambda }} \cr } } \right] = \left[ {\matrix{{{\bf{G}}{\boldsymbol P}} \cr {{{\boldsymbol P}V}} \cr } } \right], which is solved by standard numerical methods.

Let us assume that between the macroscopic quantities P_,J and F_I exists a constitutive relation analogous to the microscopic one, that is, (26) FI=AI JP,J, {F_I} = {A_{I\;J}}{P_{,J}}, where A_IJ is an apparent, homogenized material tensor. Using the discretized version of equation (14) and the linear system (25), this tensor can be derived as (27) AI J=A=1V[R−GTG⋆+(GTB⋆−IV)(BTB⋆)−1(BTG⋆−IV)], {A_{I\;J}} = {\bf{A}} = \frac{1}{V}\left[ {{\bf{R}} - {{\bf{G}}^{\bf{T}}}{{\bf{G}}^ \star } + \left( {{{\bf{G}}^{\bf{T}}}{{\bf{B}}^ \star } - {\bf{I}}V} \right){{\left( {{{\bf{B}}^{\bf{T}}}{{\bf{B}}^ \star }} \right)}^{ - 1}}\left( {{{\bf{B}}^{\bf{T}}}{{\bf{G}}^ \star } - {\bf{I}}V} \right)} \right], where I is a D × D identity matrix, R is a D × D matrix resulting from the following integration (computed by means of the finite element discretization) (28) R=RI J=∫∂ωκxixj ds, {\bf{R}} = {R_{I\;J}} = \int\limits_{\partial \omega } {\kappa {x_i}{x_j}\;{\rm{d}}s} , and G^⋆ and B^⋆ are the M × D - shaped solutions (29) G⋆=(K+D)−1G, {{\bf{G}}^ \star } = {({\bf{K}} + {\bf{D}})^{ - 1}}{\bf{G}}, (30) B⋆=(K+D)−1B, {{\bf{B}}^ \star } = {({\bf{K}} + {\bf{D}})^{ - 1}}{\bf{B}}, found using direct or iterative numerical solvers. It is clear that there is no dependence between the macroscopic tensor A_IJ and the applied macroscopic gradient P_,J. This is due to the linearity of the problem. In other words, explicit solution p of the system (25) is actually not needed for homogenization purposes. However, this solution can be always recovered for a specific P_,J = P via the linear operation (31) p=HP {\bf{p}} = {\bf{H}}{\boldsymbol P} where H is a M × D - shaped matrix given by: (32) H=G⋆+VB⋆−B⋆(BTB⋆)−1BTG⋆ {\bf{H}} = {{\bf{G}}^ \star } + V{{\bf{B}}^ \star } - {{\bf{B}}^ \star }{({{\bf{B}}^T}{{\bf{B}}^ \star })^{ - 1}}{{\bf{B}}^T}{{\bf{G}}^ \star }

Darcy flow in porous media is ruled by the equations (33) qi,iD=0, q_{i,i}^D = 0, (34) qiD=−kijγp,j, q_i^D = - \frac{{{k_{ij}}}}{\gamma }{p_{,j}}, where qiD q_i^D is a fluid flux, p_,j is a pressure gradient, k_ij is a permeability tensor depending on the position in SVE (in velocity units), and γ is a specific weight of the fluid. The minus sign signifies the problem is negatively defined, that is, the flow direction is from higher to lower pressure. Clearly, this is formally equivalent to the Laplace equation described in previous sections. As a homogenization example let us consider water flow (γ=9.81gcm3) \left( {\gamma = 9.81\frac{{\rm{g}}}{{{\rm{c}}{{\rm{m}}^3}}}} \right) in a porous material consisting of a clay matrix and the expanded shale grains surrounded by the interface layer. Isotropy of the constituents is assumed in this example. Clay permeability is given by km=6.5e−9ms {k_m} = 6.5{\rm{e}} - 9\frac{{\rm{m}}}{{\rm{s}}} , expanded shale permeability is given by ks=8e−5ms {k_s} = 8{\rm{e}} - 5\frac{{\rm{m}}}{{\rm{s}}} and the permeability of interface k_t is variable in this example, ranging from 6.5e−9ms 6.5{\rm{e}} - 9\frac{{\rm{m}}}{{\rm{s}}} to 8e−3ms 8{\rm{e}} - 3\frac{{\rm{m}}}{{\rm{s}}} [11]. The size of grains was taken as d = 0.005 m, and the expected grain volume ratio was taken as 20.5%. Ten spherical SVEs were considered which contained, on average, 32 randomly distributed spherical inclusions (see Figure 1). The interphase thickness has been taken as 0.00025 m, and this determines its overall volume ratio around 3.5%. SVEs were drawn carefully from the larger volume, containing 10⁵ inclusions distributed in a random uniform way, with the volume ratio equal to 24% (20.5% for grains and 3.5% for the interphase layer). The number of SVEs was selected in such a way that the mean value of the grain volume ratio does not differ from the expected value by more than 5% with the confidence level equal to 95%, as described in [21]. Detailed considerations related to the interdependence between the size and the required number of SVEs for obtaining meaningful homogenization results were not included in this example. Readers interested in this topic are referred to the literature [9, 12, 22]. Here, it was assumed, that the generated set of 10 relatively large SVEs can provide results which are not accidental. This assumption has been confirmed by the statistical properties of the homogenization results.

Figure 1:

In this example, the relatively large values of ℓ parameter have been considered for establishing the boundary transfer function κ. This is dictated by the fact, that in case of stiff (or highly permeable) inclusions crossing the SVE boundary, which is the case here, the Neumann kind, uniform static, boundary conditions are much more suitable. The uniform kinematic BCs would generate highly over-stiffened response, even for large SVEs [14, 18, 22]. Relatively large values of ℓ assure the solution is close to the results generated using the uniform static (uniform flux) BCs, as it was described in section 4. For comparisons, the following sequence of ℓd \frac{\ell }{d} ratios was considered in computations: [10³, 10⁴, 10⁵].

Every SVE realization was discretized into the second order, 10-node tetrahedral finite elements of variable size, resulting in meshes with over 10⁶ nodes. A system of equations (25) was then formed, for every ℓd \frac{\ell }{d} choice, and apparent homogenized permeability tensors K_IJ were found using equation (27). All calculations were performed using the fempy Python package [19] supported by the gmsh software [2] for mesh generation and umfpack [1] for fast sparse direct solver. A server equipped with 2x Intel Xeon Gold 5220R processor, 684 GB RAM, and Nvidia Tesla V100S was utilized for fast computations.

The results of homogenization, averaged over 10 apparent K_IJ tensors obtained for 10 different SVEs, are presented in Figure 2. An isotropic permeability K=KI I3 K = \frac{{{K_{II}}}}{3} and its standard error are shown since the deviator components are small. This was controlled by the ratio (35) ε=‖KI J−KδI J‖FK, \varepsilon = \frac{{{{\left\| {{K_{I\;J}} - K{\delta _{I\;J}}} \right\|}_F}}}{K}, which ranged from 0.01 to 0.035, for all different k_t and ℓd \frac{\ell }{d} combinations1

‖•‖_F in equation (35) denotes Frobenius matrix norm.

Figure 2:

As a final remark, let us note that when increasing the interphase permeability beyond the permeability of the expanded shale, the averaged K starts to increase at a faster rate (see the last points of the lines generated with ℓd \frac{\ell }{d} equal to 10³ and 10⁴). This cannot be explained physically and is considered as the limitation of the method. Namely, for increasing contrasts in material constants, accompanied with high (or low) values of ℓ parameter, the numerical errors related to the FEM solution may become important. This problem needs an attention and further investigations.

Stationary heat conduction equation is given by: (36) qi,iT=0, q_{i,i}^T = 0, (37) qiT=−λijT,j, q_i^T = - {\lambda _{ij}}{T_{,j,}} where qiT q_i^T is a heat flux, T_,j is a temperature gradient, and λ_ij is a thermal conductivity tensor. The minus sign signifies again that the problem is negatively defined, that is, the heat flow direction is from higher to lower temperature. Let us consider then a multiphase soil and the following thermal conductivities attributed to its components: for sand (Sa) λSa=8Wm⋅K {\lambda _{{\rm{Sa}}}} = 8\frac{{\rm{W}}}{{{\rm{m}} \cdot {\rm{K}}}} , for silt (Si) λSi=4Wm⋅K {\lambda _{{\rm{Si}}}} = 4\frac{{\rm{W}}}{{{\rm{m}} \cdot {\rm{K}}}} , for clay (Cl) λCl=2Wm⋅K {\lambda _{{\rm{Cl}}}} = 2\frac{{\rm{W}}}{{{\rm{m}} \cdot {\rm{K}}}} , and for water (w) λw=0.6Wm⋅K {\lambda _{\rm{w}}} = 0.6\frac{{\rm{W}}}{{{\rm{m}} \cdot {\rm{K}}}} [17]. The volume fractions of the components are taken as f_Sa = 0.3, f_Si = 0.15, f_Cl = 0.25, f_w = 0.3. The fully saturated soil is considered here and also no water flux is assumed. SVEs for this material are created using the randomly filled voxelized boxes (Figure 3). The voxel size is assumed to be constant and equal to 0.001 m. A number of sizes of SVEs containing [9³, 17³, 33³, 65³, 129³] voxels has been considered. For these sizes, a number of random material distributions – SVE realizations – were generated, that is, [81, 27, 9, 3, 1], respectively. Again, these numbers are somewhat arbitrary in this example, but the homogenization results show that they are chosen reasonably. To investigate the influence of ℓ parameter on the homogenization results, a number of choices for ℓd \frac{\ell }{d} ratio have also been considered, namely [0.001, 0.01, 0.1, 1, 10, 100, 1000].

Figure 3:

The same tools and devices as in the case of Darcy flow example were used for solving this homogenization problem. The only exception is that the finite elements are now 8-noded hexahedra with linear shape functions, which is a natural choice for the voxelized domains. The results, that is, the averaged and isotropized macroscopic conductivities Λ=ΛI I3 \Lambda = \frac{{{\Lambda _{I\;I}}}}{3} obtained for all SVE sizes and for different ℓd \frac{\ell }{d} values are presented in Figure 4. The isotropy of the resulting averages Λ_IJ was tested using the deviation measure defined by equation (35). For all SVE sizes and ℓd \frac{\ell }{d} choices, the value ɛ < 0.01 was obtained. One note that the standard errors visible in Figure 4 are decreasing with the SVE size. This justifies the decreasing number of SVEs being analyzed. All the results fall, as expected, between the upper and lower bounds given by the rule of mixtures, that is, 3.68 and 1.43Wm⋅K 1.43\frac{{\rm{W}}}{{{\rm{m}} \cdot {\rm{K}}}} , respectively (not shown in the figures for readability). Moreover, the presented lines approach asymptotically, from both sides, the limits which are identified as the upper and lower bounds of the macroscopic conductivity, obtained with uniform kinematic and uniform static BCs, respectively. The distance between the limiting values is relatively wide for small SVEs (0.35Wm⋅K for 93 voxels) (0.35\frac{{\rm{W}}}{{{\rm{m}} \cdot {\rm{K}}}}\;{\rm{for}}\;{9^3}\;{\rm{voxels}}) and it becomes increasingly narrow as the size of the SVE increases (up to 0.027Wm⋅K for 1293 voxels) ({\rm{up}}\;{\rm{to}}\;0.027\frac{{\rm{W}}}{{{\rm{m}} \cdot {\rm{K}}}}\;{\rm{for}}\;{129^3}\;{\rm{voxels}}) . Eventually, the largest SVE used in this example could be treated as a representative volume element (RVE) since it provides homogenization results which are (almost) independent on the kind of BCs applied. Another important observation is that all lines in Figure 4 cross each other at approximately the same point located at the ℓd \frac{\ell }{d} ratio, somewhere between 0.3 and 0.4. This suggests that there exists such a value of this ratio, for which the computed macroscopic conductivity is (almost) independent on the size of SVEs. In Figure 5, exactly the same results are presented, but now, the horizontal axis indicates the SVE size and the lines are plotted for different ℓd \frac{\ell }{d} ratios. The line for ℓd=0.001 \frac{\ell }{d} = 0.001 is indistinguishable from the results which are obtained with uniform kinematic BCs. Similarly, the line for ℓd=1000 \frac{\ell }{d} = 1000 is equivalent to the line generated with uniform static BCs. Other lines reside between these limits. Additional computations performed for ℓd=0.33 \frac{\ell }{d} = 0.33 show that indeed, the average macroscopic conductivity becomes almost constant for this ratio (black line), so it does not depend on the SVE size, at least for the sizes considered in this example. For comparisons, the results generated with periodic BCs have also been presented in Figure 5. These results closely follow the ℓd=1 \frac{\ell }{d} = 1 line, with slightly higher standard errors. However, this observation cannot be taken as a general rule; rather, it applies to the specific material composition considered in this example.

Figure 4:

Figure 5:

Summarizing this illustrative example, let us note that the obtained macroscopic thermal conductivities with values above 3Wm⋅K 3\frac{{\rm{W}}}{{{\rm{m}} \cdot {\rm{K}}}} are significantly higher than usually observed in laboratory for sandy clays (or clayey sands). This is because of the idealized material being analyzed, where no air or organic matter with low thermal conductivities is taken into account in the composition. The results can also be influenced by the finite element method accuracy, since the 8-noded hexahedra elements with linear shape functions are known to provide overestimated results. Again, this problem will require further investigations.

In this paper, generalized boundary conditions were used for homogenization of coefficients of the Laplace's PDE. The mesoscopic BVP was defined and analyzed from the variational perspective, and the finite element formulation of the problem was provided. The matrix equation for computing apparent properties, resulting from FEM discretization, was also provided and utilized in two examples: homogenization of the filtration coefficient of the clay amended with expanded shale and homogenization of the thermal conductivity of the soil with multiple fractions. It is shown that the computed macroscopic properties are bounded by the uniform kinematic and uniform static solutions. For the properly chosen microscopic length parameter ℓ, the generalized BCs provide reliable homogenization results, without the assumption of periodicity of the material. Thus, the generalized BCs offer an attractive unified homogenization framework for analyzing random composites.

The main challenge when using generalized BCs is the proper estimation of the ℓ parameter, which determines a value of the boundary transfer coefficient κ, via relation (18). As a rule of thumb, this parameter can often be identified to be close to the characteristic length of the microscopic inhomogeneity. For example, assuming ℓ equal to the voxel size d in the second example gives homogenization results comparable to the results obtained with periodic BCs. However, even better ℓd \frac{\ell }{d} ratio, equal to 0.33, has been identified in this example, for which generalized BCs provide the same macroscopic thermal conductivity, independently of the SVE size. On the other hand, in case of matrix-inclusion microstructure with high contrasts in parameters of the constituents crossing SVE boundary, the limiting values of ℓ must be used to obtain reliable homogenization results – as shown and explained in the first example. So, it seems, there is no simple answer to the question, what ℓ should be actually used for the specific microstructure. One of the possible solutions to this problem is implementation of the self-consistent approach, where the boundary transfer coefficient κ is a variable quantity depending on the current apparent solution A_IJ. System of equations (25) will then become non-linear, but it should be still fast and easy to be solved by means of iterative solvers with simple Jacobi preconditioning. This topic is currently a subject of the author's investigations.