Micromechanics-based Material Properties of Inorganic Shale

The scope of this project encompasses finding the material properties of the inorganic shale matrix using Maxwell's scheme which is the most approximate scheme for anisotropic multiphase materials. The inorganic shale is considered to be consisting of four ingredients; illite dominated clay, quartz, tetrahedral and elliptical pores. The clay particles are taken as elliptical particles while the quartz is a crystal with tetrahedral shape. The shape of elliptical pores is affected mostly by the clay so, they have the same shape as elliptical clay particles. This elliptical shape depends on the compaction of the shale due to stresses, which is affected by the ratio of the minerals. A good accordance is found out between the reproduced results and the benchmark data, else along the axis of TI symmetry.

Introduction

There have been a lot of efforts to estimate elastic properties of the sedimentary rocks in terms of their mineral composition [@Louis2018; @Vernik2020]. Sevostianov and Vernik [@Sevostianov2021] presented a micromechanics-based rock-physics model for inorganic shale using the Maxwell's scheme in the interpretation of Sevostianov and Giraud [@Sevostianov2013]. In the presented model, the inorganic shale is taken to be consisted of a binary mixture of illite dominated clay and silt particles of quartz, along with two shapes of pores. The shape of pores strongly depends on the sedimentary compaction. It has been shown that compacted mudstones are typically characterized by the preferred clay platelet alignment in the bedding plane, which is also the place of transverse isotropic symmetry in these rocks. Clay platelet alignment functions are developed and calibrated by the experimental results on X-ray goniometry. The orientation functions can be defined as a single parameter function, relying on the aspect ratio of the strain ellipse related to the uniaxial sediment compaction. The Owens-March function is defined as

$$ODF(\theta_3) = \frac{MPD}{(cos^2 \theta_3 + MPD sin^2 \theta_3)}$$

where maximum plot density (MPD) is the aspect ration of the strain ellipse and $\theta_3$ is the pole orientation angle. $\theta_1$ is taken as the angle about the TI symmetry axis and the orientation of $\theta_2$ and $\theta_3$ can be taken arbitrarily because of the symmetry in these orientations. The unit for MPD is Multiples of Random Distribution (MRD) and lower number depict higher degree of misalignment. MPD = 1.0 means the composition minerals are fully random in their alignment and $MPD = \infty$ being fully aligned.

In this work, the results produced by [@Sevostianov2021], are being reproduced as most of the data was lost. The same scheme is being followed as presented by the authors. In this report, the readers will get more in-depth information about the practical implementation and computation of Maxwell's scheme of the subject research problem. For computations, MATLAB codes have been developed, and a good accordance is found between the bench-marked results, and those reproduced in this study.

Scheme of Approach

Composition of Inorganic Shale

The inorganic shale is composed of a binary porous mixture of two minerals clay and quartz. As an initial step, the pores are neglected in this study and different ratios of the minerals are studied in respect of the resulting elastic properties of the shale. Quartz are crystal particles and are considered to have purely isotropic properties. On the other hand, clay particles are transversely isotropic with axis of symmetry normal to the bedding plane. The shape of clay particles is oblate spheroids. When the rock is totally free of stresses the clay particles are fully randomly oriented, eventually giving an isotropic elastic properties. In natural compaction process of the rocks, the constituents are presses in the bedding plane, and the clay particles are somewhat aligned with each other giving over all TI elastic behavior.

Maxwell's Scheme

Effective media approaches and the effective field scheme are the most commonly used for inhomogeneity. [@Sevostianov2019] shows that the most appropriate scheme for anisotropic materials is the Maxwell's scheme in terms of property contribution tensors. A representative volume element (RVE) $V_\omega$ from the multiphase material is considered to be placed into a matrix material of know properties. The effective produced by this RVE is presented by the stiffness contribution tensor which is a volumetric weighted average of the stiffness contribution tensors of all the inhomogeneities of the RVE. For the case of the ellipsoidal shape of RVE, the effective elastic properties can be found using $$C^{eff} = C^0 + ((\frac{1}{V_\Omega} \sum_i V_i N^{(i)})^{-1}) - P_\Omega)^{-1} \label{maxwell}$$

In this equation $N^{(i)}$ is the stiffness contribution tensor for each of the inhomogeneity in the mixture. In this study, there are only two inhomogeneities; clay and quartz. $P_Omega$ is the Hill's tensor [@Hill1963] of the domain RVE ($\Omega$). The stiffness contribution tensor and the Hills Tensor are a function of elastic properties and shape of the inclusion and are given in detail in [@Kachanov2018]. Once, stiffness contribution tensors of the inclusions, $N_c$ for clay and $N_q$ for quartz are known, the above equation can be simply written as

$$C^{eff} = C^0 + (( V_q N_q + V_c N_c)^{-1} - P_\Omega)^{-1} \label{maxwell2}$$

The volume of the RVE is taken to be unity because there are only two inclusions in the mixture and both the volumes $V_c$ and $V_q$ combined makes the whole volume where these are volumetric ratios of the inclusions clay and quartz, respectively.

Orientation Distribution

Another parameter of study in the effective properties of the shale matrix, is the orientation distribution because of the compaction, with causes a change in the properties and property contribution tensor for the clay. As quartz is isotropic so, ODF does not affect material properties or property contribution tensors of quartz. The effective stiffness contribution tensor of clay is found out through the orientation distribution function. ODF gives the volume ratio for each orientation about the symmetric and any of the non-symmetric axis. The Owens-March distribution function is defined for one axis (or 2D distribution) but we need a 3D distribution. This will be elaborated in the methodology section later. In the methodology section, the implementation of this scheme will be discussed in implementation perspective. In MATLAB, the tensors are notated and calculated in matrix form using binary Voigt notation. So, we will first review the tensors in the interpretation of matrices.

Tensors

In mechanics stress, strain, stiffness and compliance are written in tensors notation essentially, because of their nature. The stress in 3D can be written in 2x2 tensor as follows

$$\sigma_{ij} = \begin{bmatrix} \sigma_{11} & \sigma_{12} & \sigma_{13}\\ \sigma_{21} & \sigma_{22} & \sigma_{23}\\ \sigma_{31} & \sigma_{32} & \sigma_{33} \end{bmatrix}$$

and the strain can be written as $$\sigma_{ij} = \begin{bmatrix} \epsilon_{11} & \epsilon_{12} & \epsilon_{13}\\ \epsilon_{21} & \epsilon_{22} & \epsilon_{23}\\ \epsilon_{31} & \epsilon_{32} & \epsilon_{33} \end{bmatrix}$$

The stress and strain are related by stiffness, in tensor notation as

$$\sigma = C:\epsilon \label{stiffness_strain_tensor}$$

The stiffness tensor is a 4th rank tensor which may not be very easily visualized on a paper. But here I will try to elaborate a little bit, using its resemblance with 4-D array, or a 2-D matrix of 2-D matrix. It gets even simpler to understand when these 2-D matrices are square with dimension 3x3. We would write

$$C_{ijkl} = \begin{bmatrix} \begin{bmatrix} C_{1111} & C_{1112} & C_{1113}\\ C_{1121} & C_{1122} & C_{1123}\\ C_{1131} & C_{1132} & C_{1133} \end{bmatrix} & \begin{bmatrix} C_{1211} & C_{1212} & C_{1213}\\ C_{1221} & C_{1222} & C_{1223}\\ C_{1231} & C_{1232} & C_{1233} \end{bmatrix} & \begin{bmatrix} C_{1311} & C_{1312} & C_{1313}\\ C_{1321} & C_{1322} & C_{1323}\\ C_{1331} & C_{1332} & C_{1333} \end{bmatrix} \\ \begin{bmatrix} C_{2111} & C_{2112} & C_{2113}\\ C_{2121} & C_{2122} & C_{2123}\\ C_{2131} & C_{2132} & C_{2133} \end{bmatrix} & \begin{bmatrix} C_{2211} & C_{2212} & C_{2213}\\ C_{2221} & C_{2222} & C_{2223}\\ C_{2231} & C_{2232} & C_{2233} \end{bmatrix} & \begin{bmatrix} C_{2311} & C_{2312} & C_{2313}\\ C_{2321} & C_{2322} & C_{2323}\\ C_{2331} & C_{2332} & C_{2333} \end{bmatrix} \\ \begin{bmatrix} C_{3111} & C_{3112} & C_{3113}\\ C_{3121} & C_{3122} & C_{3123}\\ C_{3131} & C_{3132} & C_{3133} \end{bmatrix} & \begin{bmatrix} C_{3211} & C_{3212} & C_{3213}\\ C_{3221} & C_{3222} & C_{3223}\\ C_{3231} & C_{3232} & C_{3233} \end{bmatrix} & \begin{bmatrix} C_{3311} & C_{3312} & C_{3313}\\ C_{3321} & C_{3322} & C_{3323}\\ C_{3331} & C_{3332} & C_{3333} \end{bmatrix} \end{bmatrix} \label{4tensor}$$

The same can be written in components form, as $$\sigma_{ij} = C_{ijkl}\epsilon_{kl}$$ This components form (called as indices notation) means
$$\begin{gathered} \sigma_{11} = C_{1111}\epsilon_{11} + C_{1112} \epsilon_{12} + C_{1113}\epsilon_{13} + C_{1121} \epsilon_{21} + C_{1122} \epsilon_{22} + C_{1123}\epsilon_{11} + C_{1131} \epsilon_{31} + C_{1132} \epsilon_{32} + C_{1133}\epsilon_{33} \label{comp_long_eq} \end{gathered}$$

or it can be thought of a component wise multiplication of 3x3 stiffness matrix (2nd rank tensor) and 3x3 strain matrix (2nd rank tensor), as follows
$$ \sigma_{11} = \begin{bmatrix} C_{1111} & C_{1112} & C_{1113}\\ C_{1121} & C_{1122} & C_{1123}\\ C_{1131} & C_{1132} & C_{1133} \end{bmatrix}: \begin{bmatrix} \epsilon_{11} & \epsilon_{12} & \epsilon_{13}\\ \epsilon_{21} & \epsilon_{22} & \epsilon_{23}\\ \epsilon_{31} & \epsilon_{32} & \epsilon_{33} \end{bmatrix} $$ and the other components of stress matrix (or 2nd order tensor) can be written as
$$\sigma_{ij} = \begin{bmatrix} C_{ij11} & C_{ij12} & C_{ij13}\\ C_{ij21} & C_{ij22} & C_{ij23}\\ C_{ij31} & C_{ij32} & C_{ij33} \end{bmatrix}: \begin{bmatrix} \epsilon_{11} & \epsilon_{12} & \epsilon_{13}\\ \epsilon_{21} & \epsilon_{22} & \epsilon_{23}\\ \epsilon_{31} & \epsilon_{32} & \epsilon_{33} \end{bmatrix}$$ Remember, each of these colon (:) signed operation is summation over component wise multiplication, and not the traditional matrix rotation. Of course, this type of multiplication is not very new for all of us, so an analogy can be created, by writing this 3x3 2nd order tensor in the form of a 1D array (vector) or 1x9 matrix. Writing stiffness tensor for one component $\sigma_{11}$ in a row vector and the strain tensor as column vector, the above relations can be written in the form of simple vector notation as follows $$ \sigma_{ij} = \begin{bmatrix} C_{ij11} & C_{ij12} & C_{ij13}& C_{ij21} & C_{ij22} & C_{ij23}& C_{ij31} & C_{ij32} & C_{ij33} \end{bmatrix} \begin{bmatrix} \epsilon_{11} \\ \epsilon_{12} \\ \epsilon_{13}\\ \epsilon_{21} \\ \epsilon_{22} \\ \epsilon_{23}\\ \epsilon_{31} \\ \epsilon_{32} \\ \epsilon_{33} \end{bmatrix} $$

It is quite basic of mechanics, that the shear stains are symmetric as $\epsilon_{21} = \epsilon_{12}$, $\epsilon_{31} = \epsilon_{13}$ and $\epsilon_{23} = \epsilon_{32}$, and knowing the same symmetry for stiffness 2nd order tensors, Eq. (\ref{comp_long_eq}) can be rewritten as $$\sigma_{11} = C_{1111}\epsilon_{11} + C_{1122} \epsilon_{22} + C_{1133}\epsilon_{33} + 2C_{1123}\epsilon_{11} + 2C_{1113}\epsilon_{13} + 2C_{1112} \epsilon_{12} \label{comp_long2_eq}$$

With this equation the vector form can be contracted as $$\sigma_{ij} = \begin{bmatrix} C_{ij11} & C_{ij22} & C_{ij33}& 2C_{ij23} & 2C_{ij13} & 2C_{ij12} \end{bmatrix} \begin{bmatrix} \epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{33}\\ \epsilon_{23} \\ \epsilon_{13} \\ \epsilon_{12} \end{bmatrix}$$

As, the overall result for any component $\sigma_{ij}$ is just a summation, the order of these components in Eq. (\ref{comp_long_eq}) and Eq. (\ref{comp_long2_eq}) is arbitrary. But it is the most followed notation of order which is used here. And eventually, these indices are also contracted as per the vector notation as follows
11
22
33
23
13
12
Now, one last step is left to make this tensor notation of stress, strain, and stiffness, pretty simple $$\sigma_{i} = \begin{bmatrix} C_{i1} & C_{i2} & C_{i3}& 2C_{i4} & 2C_{i5} & 2C_{i6} \end{bmatrix} \begin{bmatrix} \epsilon_{1} \\ \epsilon_{2} \\ \epsilon_{3}\\ \epsilon_{4} \\ \epsilon_{5} \\ \epsilon_{6} \end{bmatrix}$$

Remember, we had the choice to write this equation as follows as well

$$\sigma_{i} = \begin{bmatrix} C_{i1} & C_{i2} & C_{i3}& C_{i4} & C_{i5} & C_{i6} \end{bmatrix} \begin{bmatrix} \epsilon_{1} \\ \epsilon_{2} \\ \epsilon_{3}\\ 2\epsilon_{4} \\ 2\epsilon_{5} \\ 2\epsilon_{6} \end{bmatrix}$$

Each one of this shear strain elements are sometime written as $2\epsilon_i = \gamma_i$
Note: In both of the above equations i = 1 to 6 and not 1 to 3 as in all previous equations and relations
Now we have everything sorted out so, we can define our stress strain relationship (for all $\sigma_i$) in vector and matrix form

$$\begin{bmatrix} \sigma_{1} \\ \sigma_{2} \\ \sigma_{3}\\ \sigma_{4} \\ \sigma_{5} \\ \sigma_{6} \end{bmatrix} = \begin{bmatrix} C_{11} & C_{12} & C_{13}& C_{14} & C_{15} & C_{16}\\ C_{21} & C_{22} & C_{23}& C_{24} & C_{25} & C_{26}\\ C_{31} & C_{32} & C_{33}& C_{34} & C_{35} & C_{36}\\ C_{41} & C_{42} & C_{43}& C_{44} & C_{45} & C_{46}\\ C_{51} & C_{52} & C_{53}& C_{54} & C_{55} & C_{56}\\ C_{61} & C_{62} & C_{63}& C_{64} & C_{65} & C_{66} \end{bmatrix} \begin{bmatrix} \epsilon_{1} \\ \epsilon_{2} \\ \epsilon_{3}\\ 2\epsilon_{4} \\ 2\epsilon_{5} \\ 2\epsilon_{6} \end{bmatrix}$$

Similarly, the compliance is defined as $$\epsilon = S:\sigma \label{stiffness_stress_tensor}$$

And the 4th rank tensor is written as

$$S_{ijkl} = \begin{bmatrix} \begin{bmatrix} S_{1111} & S_{1112} & S_{1113}\\ S_{1121} & S_{1122} & S_{1123}\\ S_{1131} & S_{1132} & S_{1133} \end{bmatrix} & \begin{bmatrix} S_{1211} & S_{1212} & S_{1213}\\ S_{1221} & S_{1222} & S_{1223}\\ S_{1231} & S_{1232} & S_{1233} \end{bmatrix} & \begin{bmatrix} S_{1311} & S_{1312} & S_{1313}\\ S_{1321} & S_{1322} & S_{1323}\\ S_{1331} & S_{1332} & S_{1333} \end{bmatrix} \\ \begin{bmatrix} S_{2111} & S_{2112} & S_{2113}\\ S_{2121} & S_{2122} & S_{2123}\\ S_{2131} & S_{2132} & S_{2133} \end{bmatrix} & \begin{bmatrix} S_{2211} & S_{2212} & S_{2213}\\ S_{2221} & S_{2222} & S_{2223}\\ S_{2231} & S_{2232} & S_{2233} \end{bmatrix} & \begin{bmatrix} S_{2311} & S_{2312} & S_{2313}\\ S_{2321} & S_{2322} & S_{2323}\\ S_{2331} & S_{2332} & S_{2333} \end{bmatrix} \\ \begin{bmatrix} S_{3111} & S_{3112} & S_{3113}\\ S_{3121} & S_{3122} & S_{3123}\\ S_{3131} & S_{3132} & S_{3133} \end{bmatrix} & \begin{bmatrix} S_{3211} & S_{3212} & S_{3213}\\ S_{3221} & S_{3222} & S_{3223}\\ S_{3231} & S_{3232} & S_{3233} \end{bmatrix} & \begin{bmatrix} S_{3311} & S_{3312} & S_{3313}\\ S_{3321} & S_{3322} & S_{3323}\\ S_{3331} & S_{3332} & S_{3333} \end{bmatrix} \end{bmatrix} \label{4tensor_C}$$

$$\epsilon_{ij} = \begin{bmatrix} S_{ij11} & S_{ij12} & S_{ij13}\\ S_{ij21} & S_{ij22} & S_{ij23}\\ S_{ij31} & S_{ij32} & S_{ij33} \end{bmatrix}: \begin{bmatrix} \sigma_{11} & \sigma_{12} & \sigma_{13}\\ \sigma_{21} & \sigma_{22} & \sigma_{23}\\ \sigma_{31} & \sigma_{32} & \sigma_{33} \end{bmatrix}$$

$$\epsilon_{ij} = \begin{bmatrix} S_{ij11} & S_{ij12} & S_{ij13}& S_{ij21} & S_{ij22} & S_{ij23}& S_{ij31} & S_{ij32} & S_{ij33} \end{bmatrix} \begin{bmatrix} \sigma_{11} \\ \sigma_{12} \\ \sigma_{13}\\ \sigma_{21} \\ \sigma_{22} \\ \sigma_{23}\\ \sigma_{31} \\ \sigma_{32} \\ \sigma_{33} \end{bmatrix}$$

$$\epsilon_{ij} = \begin{bmatrix} S_{ij11} & S_{ij22} & S_{ij33}& 2S_{ij23} & 2S_{ij13} & 2S_{ij12} \end{bmatrix} \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{23} \\ \sigma_{13} \\ \sigma_{12} \end{bmatrix}$$

$$\begin{bmatrix} \epsilon_{1} \\ \epsilon_{2} \\ \epsilon_{3} \\ 2\epsilon_{4} \\ 2\epsilon_{5} \\ 2\epsilon_{6} \end{bmatrix} = \begin{bmatrix} S_{11} & S_{12} & S_{13}& 2S_{14} & 2S_{15} & 2S_{16}\\ S_{21} & S_{22} & S_{23}& 2S_{24} & 2S_{25} & 2S_{26}\\ S_{31} & S_{32} & S_{33}& 2S_{34} & 2S_{35} & 2S_{36}\\ 2S_{41} & 2S_{42} & 2S_{43}& 4S_{44} & 4S_{45} & 4S_{46}\\ 2S_{51} & 2S_{52} & 2S_{53}& 4S_{54} & 4S_{55} & 4S_{56}\\ 2S_{61} & 2S_{62} & 2S_{63}& 4S_{64} & 4S_{65} & 4S_{66} \end{bmatrix} \begin{bmatrix} \sigma_{1} \\ \sigma_{2} \\ \sigma_{3}\\ \sigma_{4} \\ \sigma_{5} \\ \sigma_{6} \end{bmatrix}$$

Methodology

In this section, the implementation of Maxwell's scheme will be focused and the procedure in detail will be discussed.

Isotropic equivalent of Transversely Isotropic Material

If a TI material is distributed randomly in all orientations, the effective properties are purely isotropic, and an effective bulk and shear modulii can be approximated, which are required in the calculations of Hill's tensors. [@Vernik2020] mentioned a direct formulation to interpret TI material as an approximated isotropic material. With this formulation isotropic equivalent of the clay material can be found in terms of $C_{11}$ and $C_{66}$. Following MATLAB code was written

`` {.pre code language="MATLAB"} Kc_v = (1/9)(4Cc(1,1)+Cc(3,3) +4Cc(1,3) - 4Cc(6,6)); Gc_v = (1/30)(2Cc(1,1) + 2Cc(3,3) - 4Cc(1,3) + 10Cc(6,6) + 12Cc(4,4));

Cc_v_11 = Kc_v + (4/3)Gc_v; Cc_v_13 = Cc_v_11 - 2Gc_v; Cc_v = ToMatrix([Cc_v_11, Cc_v_11, Gc_v, Gc_v, Cc_v_13]); ``

Another way of doing the same is using ODF with an MPD = 1. Both gives numerically the same results.

$$C_c_v = \begin{bmatrix} 66.0 & 29.3 & 29.3 & 0.0 & 0.0 & 0.0 \\ 29.3 & 66.0 & 29.3 & 0.0 & 0.0 & 0.0 \\ 29.3 & 29.3 & 66.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & 18.4 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & 0.0 & 18.4 & 0.0 \\ 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 18.4 \\

\end{bmatrix}$$

Hill's Tensor

The stiffness contribution tensor are also found using the expressions given in [@Kachanov2018].

if (gamma<1)
    g = (1/(gamma*sqrt(1-gamma^2)))*atan(sqrt(1-gamma^2)/gamma);
    f0 = (gamma^2*(1-g))/(2*(gamma^2-1));
    f1 = (gamma^2*((2*gamma^2 +1)*g -3))/(4*(gamma^2 -1)^2);
    
else
    g = 1;
    f0 = 1/3;
    f1 = 1/15;   
end

P11 = (f0*(4-3*ko)+3*ko*f1)/(4*uo);
P33 = ((1-2*f0)*(1-ko) + 2*ko*f1)/(uo);
P66 = (f0*(2-ko) + ko*f1)/(4*uo);
P55 = (1-f0 - 4*ko*f1)/(4*uo);
P13 = -(ko*f1)/uo;
P12 = ko*(-f0+f1)/(4*uo);

P = [P11,P12,P13,0,0,0;
    P12,P11,P13,0,0,0;
    P13,P13,P33,0,0,0;
    0,0,0,4*P55,0,0;
    0,0,0,0,4*P55,0; 
    0,0,0,0,0,4*P66];
    

where in the above equation is the aspect ratio of the ellipsoid. For clay, it is experimentally found that the aspect ratio is 0.1 and for quartz is 1, as it is fully spherical in shape. $K_0$ and $G_0$ are the bulk and shear modulii of the respective inclusion.

Stiffness Contribution Tensor

Stiffness contribution tensor is defined by the stiffness properties of the matrix material and the inclusion and Hill's tensor of the material. As defined in [@Sevostianov2014], $$N_i = ((C_i - C_0)^{-1} + P_i )^{-1}$$

where the $C_0$ is the stiffness of the fictitious matrix material in which the inclusions (RVE) is place. The material properties of this fictitious matrix material are volumetric average of the material properties of the two inclusions.

$$C_0 = V_c C_c_v + V_q C_q$$ Note that we use the isotropic equivalent properties of the clay, though it is TI in nature.

Shape of the RVE

The shape of RVE depends upon the shape or Hill's tensors of the inclusions and there has been a lot of discussions about the choice of shape of RVE in [@Sevostianov2008; @Sevostianov2012; @Sevostianov2019; @Sevostianov2014]. The aspect ratio of the representative volume is $$\gamma_\Omega = \frac{V_c P_{1111}^{(c)} + V_q P_{1111}^{(q)} }{V_c P_{3333}^{(c)} + V_q P_{3333}^{(q)} }$$

Orientation Distribution Function

The given Owens-March ODF function define the distribution in 2D. The resulting stiffness matrix will average out the axis-1 (the axis of symmetry) and axis-3 (one of the axis with non-symmetry), but the material properties along the other non-symmetric axis, axis-3 will remain unchanged. For a 3D problem we need the ODF function to be defined in 3D.

vfc1 =(sind(theta1)*mpd)./((8*pi^2).*((cosd(theta1)).^2 + (mpd).*(sind(theta1)).^2).^(3/2));
vfc3 = zeros(1,Nd)+1;
vfc_2d = transpose(vfc3) * vfc1;
vfcn_2d = vfc_2d./sum(sum(vfc_2d));

The line 4 in the above script, normalizes the ODF function. vfcn_2d is a 2-D with an ODF and carrying the volume ratios of the rotated stiffness matrix by an amount $theta_1$ around axis-1 and $theta_3$ around axis-3.

Weighted average by their respective volumes of these rotated tensors, give the effective property (stiffness compliance or stiffness contribution tensor) tensor.

for j = 1:Nd
    for  k = 1:Nd
    M_ = M_ + vfcn_2d(k,j)*RotationTensorNotation(M,theta1(j),0,theta3(k));
    end
end

Depending upon the value of MPD input in the above function, the degree of anisotropy in the input tensor will be changed. Putting MPD = 1, will results in equivalent isotropic properties. RotationTensorNotation returns the rotated tensor about axis-1 and axis-3.

Implementation of Maxwell's Scheme

Once all the ingredients of the Eq. [maxwell]{reference-type="ref" reference="maxwell"} are ready, the effective stiffness can be found

P_q = ToHills(1,Go,Ko);
P_c = ToHills(0.1,Go,Ko);
    
    
    
Nc = inv(inv(Cc - Co)+ P_c);
Nq = inv(inv(Cq - Co)+ P_q);

Nc_ODF = ODF(Nc,MPD);
N = Vcl(i)*Nc_ODF +Vq(i)*Nq;
    
gamma(i) = (Vcl(i)*Q_c(3,3) + Vq(i)*Q_q(3,3))/(Vcl(i)*Q_c(1,1) + Vq(i)*Q_q(1,1));

P = ToHills(gamma(i),Go,Ko);
C_eff = Co + inv(inv(N) - P);

Results

The effective material properties as reproduced in this study, are shown in Fig. 1{reference-type="ref" reference="results"} as a comparison to the bench-marked data.

It can be seen that lower values of the MPD, the reproduced results are quite matching in all aspects. All the parameters of the stiffness tensor of the shale matrix (a mixture of clay and quartz) have been reproduced. As the we increase the MPD, and move towards some small degree of anisotropy, a mismatch can be seen specifically, $C11$ and $C33$, which define the degree of anisotropy. The reason behind this mismatch is mostly because of orientation distribution function, about which the implementation is not very clearly defined in the literature.

Conclusion

The material elastic properties of the inorganic shale matrix, which is considered to be composed of a binary mixture of clay and quartz crystal, under sedimentary compaction, has been reproduced for lower values of Maximum Pole Density (MPD), with an error of lower than 5 percent. This study provided a systemic implementation approach to Maxwell's scheme of homogeneity in the interpretation of stiffness contribution tensors with non-interactive approximation approach. It has been found that the developed implementation model in this work, does not produce good results for higher MPD values, thus this area needs to be explored, in future.