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

: In the paper, generalized boundary conditions were used for the homogenization of coefficients of the Laplace partial differential equation in the context of Darcy flow and heat diffusion phenomena. The mesoscopic boundary value problem was defined and analyzed from the variational perspective and the finite element formulation of the homogenization problem was provided. The matrix equation for the apparent macroscopic properties, resulting from FEM discretization, was derived and utilized in two illustrative examples: homogenization of the filtration coefficient of clay amended with expanded shale and thermal conductivity of the soil with multiple fractions. It is shown, that generalized boundary conditions can provide very good homogenization results without the assumption of the periodicity of the material. For best results, the microscopic length parameter has to be properly estimated.


Introduction
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-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.

Microscopic BVP -Laplace equation
Let us consider anisotropic Laplace equation defined on SVE : where   is a tensor function of parameters (positive definite) and  is an unknown function.In the scope of computational homogenization, BCs for this PDE are determined by the averaging equation: where  , is a known quantity, which usually results from the upper scale computations,  is the SVE volume and  is the SVE boundary.Weak form of equation ( 1) reads as with  being an arbitrary testing function belonging to the same function space as .Applying Gauss-Ostrogradsky theorem the following is obtained: where  is interpreted as the boundary flux.In the spirit of generalized BCs [22], the following is assumed at boundary: where  is called the boundary transfer function and is the perturbation of  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 where Λ  are the mentioned multipliers which are real numbers.

Variational considerations
Solution of problem ( 7) is equivalent to the unconstrained minimization: min where potential Π is given by Minimum is found by equating variations of Π in  and Λ  to 0, that is, The above equations must be fulfilled for any admissible variations  and Λ  .So, let us consider  =   1  , and thus,  , = 1  .From the first equation, the following is obtained: Recognizing that    , =   (13) 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   , and: Thus, Lagrange multipliers stand for the macroscopic flux (with minus sign) modified by the component related to the boundary perturbation p. Usually this second component is enforced to vanish by applying 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  = , and thus,  , =  , .Again, from equation (10), the following is obtained: Inserting ( 2), (13), and ( 14) into (15), after some transformations, the following is derived: 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.

Selection of 𝜅
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,  :  → R + .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, so this is equivalent in applying the socalled 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 vanishes.This is equivalent in applying the kinematic uniform BCs, where  =  ,   is enforced on .If the macroscopic fluxes obtained with these two  selections and by application of the macroscopic gradient  , are denoted as  0  and  ∞  , then the following holds: where   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: where ℓ is interpreted as the microscopic length parameter,   describes the material crossing the boundary, and   is a unitary vector, normal to .Basic assumption which stands behind this choice is that at the boundary, the following holds: If the SVE neighbourhood is taken exactly the same as at the boundary, then this assumption leads to the following external loading term: As a final remark, it should be pointed out that the selected  depends on the boundary orientation, and thus, it is variable even if   is constant all over the boundary.However, for the case of isotropic materials, where   =   , with  being a constant scalar function and   being Kroenecker delta, the orientation becomes unimportant and  =  ℓ .

Finite element implementation
Solutions  of the problem defined by equations ( 9), (10), and ( 11) have to belong to the Sobolev functional space  1 () 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 and skipping details of the finite element approximation (see [20,22] for that) one ends up immediately in discreet version of this potential: where p is global vector of unknowns of length  (total number of nodes in discretization), K is a global linear operator of size  ×  , D is the boundary constraint matrix, also of size  ×  , G is a boundary loading matrix of size  ×  ( -space dimension: 2 or 3), B is a matrix resulting from the application of averaging constraint, also of size  × , Λ is a vector of unknown Lagrange multipliers of length , and  is a vector of known macroscopic gradient (also of length

𝐷). Minimum of potential ̂︀
Π is found by equating its variations in p and Λ to 0, which is denoted as This is a system of linear equations which is solved by standard numerical methods.

Homogenized material tensor
Let us assume that between the macroscopic quantities  , and   exists a constitutive relation analogous to the microscopic one, that is, where   is an apparent, homogenized material tensor.
Using the discretized version of equation ( 14) and the linear system (25), this tensor can be derived as where I is a  ×  identity matrix, R is a  ×  matrix resulting from the following integration (computed by means of the finite element discretization) and G ⋆ and B ⋆ are the  ×  -shaped solutions found using direct or iterative numerical solvers.It is clear that there is no dependence between the macroscopic tensor   and the applied macroscopic gradient  , .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  , =  via the linear operation where H is a  ×  -shaped matrix given by:

Example 1 -hydraulic conductivity
Darcy flow in porous media is ruled by the equations where    is a fluid flux,  , is a pressure gradient,   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.81 g cm 3 ) 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   = 6.5e−9 m s , expanded shale permeability is given by   = 8e−5 m s and the permeability of interface   is variable in this example, ranging from 6.5e−9 m s to 8e−3 m s [11].The size of grains was taken as  = 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 5 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.
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 ℓ  ratios was considered in computations: [10 3 , 10 4 , 10 5 ].
Every SVE realization was discretized into the second order, 10-node tetrahedral finite elements of variable size, resulting in meshes with over 10 6 nodes.A system of equations (25) was then formed, for every ℓ  choice, and apparent homogenized permeability tensors   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   tensors obtained for 10 different SVEs, are presented in Figure 2.An isotropic permeability  =   3 and its standard error are shown since the deviator components are small.This was controlled by the ratio which ranged from 0.01 to 0.035, for all different   and ℓ  combinations1 1 .The results show that for increasing ℓ parameter the solutions approach the limit established by the uniform static BCs -the results obtained for  increasing the   value from   to   , the macroscopic permeability  also increases, approximately from the Hashin-Shtrikman (H-S) lower bound when   =   to the same bound when   =   .This coincidence between the numerically obtained values and the H-S bound is a bit surprising, but actually it should be expected, since the H-S bounds are derived exactly for this kind of microstructure, that is, for spherical inclusions in the matrix surrounded by equivalent material.Thus, this result can be treated as the confirmation of correctness of the considered homogenization approach.In case of ℓ  ratio equal to 10 3 , the results are located above the H-S lower bound.With further decrease of ℓ  (increase of ), the macroscopic  would eventually approach the upper H-S bound.This case is not considered here because of the reasons mentioned before (but it is considered in the next example).
As a final remark, let us note that when increasing the interphase permeability beyond the permeability of the expanded shale, the averaged  starts to increase at a faster rate (see the last points of the lines generated with ℓ  equal to 10 3 and 10 4 ).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.

Example 2 -heat conduction
Stationary heat conduction equation is given by: where    is a heat flux,  , is a temperature gradient, and   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) , and for water (w)  w = 0.6 W m•K [17].The volume fractions of the components are taken as  Sa = 0.3,  Si = 0.15,  Cl = 0.25,  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 3 , 17 3 , 33 3 , 65 3 , 129 3 ] 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 ℓ  ratio have also been considered, namely [0.001, 0.01, 0.1, 1, 10, 100, 1000].
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 Λ = Λ  3 obtained for all SVE sizes and for different ℓ  values are presented in Figure 4.The isotropy of the resulting averages Λ  was tested using the deviation measure defined by equation ( 35).For all SVE sizes and ℓ  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.43 W m•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.35 W m•K for 9 3 voxels) and it becomes increasingly narrow as the size of the SVE increases (up to 0.027 W m•K for 129 3 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 ℓ  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 ℓ  ratios.The line for ℓ  = 0.001 is indistinguishable from the results which are obtained with uniform kinematic BCs.Similarly, the line for ℓ  = 1000 is equivalent to the line generated with uniform static BCs.Other lines reside between these limits.Additional computations performed for ℓ  = 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 ℓ  = 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.
Summarizing this illustrative example, let us note that the obtained macroscopic thermal conductivities with values above 3 W m•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.

Conclusions
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  in the second example gives homogenization results comparable to the results obtained with periodic BCs.However, even better ℓ  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   .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.

Figure 1 :
Figure 1: Example of SVE used for homogenization of permeability coefficient.Approximately 32 spherical inclusions are distributed in a random uniform way and surrounded by the interphase layer.

Figure 2 :
Figure 2: Averaged apparent permeability  relative to the interphase filtration coefficient , for different values of the ℓ  ratio.Standard error, shown on the error bars, never exceeds 4%.

Figure 4 :
Figure 4: Averaged and isotropized macroscopic thermal conductivity Λ relative to the ℓ  ratio and for different SVE sizes.Standard error, shown on the error bars, never exceeds 1.5%.

Figure 5 :
Figure 5: Averaged and isotropized macroscopic thermal conductivity Λ relative to the SVE size and for different ℓ  ratios.In addition, the results obtained with periodic BCs are presented.Standard error, shown on the error bars, never exceeds 1.5%.