Department of Mechanical & Aerospace Engineering
New Mexico State University
Las Cruces, New Mexico 88003
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^ = C0 + ((\frac{1}{V_\Omega} \sum_i V_i N{(i)}){-1}) - P_\Omega){-1}
\[
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}
\]
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_ = \begin\sigma_{11} & \sigma_{12} & \sigma_{13}\\
\sigma_{21} & \sigma_{22} & \sigma_{23}\\
\sigma_{31} & \sigma_{32} & \sigma_{33}
\end$$
and the strain can be written as $$\sigma_ = \begin\epsilon_{11} & \epsilon_{12} & \epsilon_{13}\\
\epsilon_{21} & \epsilon_{22} & \epsilon_{23}\\
\epsilon_{31} & \epsilon_{32} & \epsilon_{33}
\end$$
The stress and strain are related by stiffness, in tensor notation as
$$\sigma = C:\epsilon
\[
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}
\]
The same can be written in components form, as
\(\sigma_{ij} = C_{ijkl}\epsilon_{kl}\) This components form (called as
indices notation) means
$$\begin\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}
\end$$
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_ = \beginC_ & C_ & C_\\
C_ & C_ & C_\\
C_ & C_ & C_\end:
\begin\epsilon_{11} & \epsilon_{12} & \epsilon_{13}\\
\epsilon_{21} & \epsilon_{22} & \epsilon_{23}\\
\epsilon_{31} & \epsilon_{32} & \epsilon_{33}
\end$$ 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. (5) 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}
\[
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.
\(5\) and Eq.
\(6\) 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
\]
And the 4th rank tensor is written as
$$S_ = \begin\beginS_{1111} & S_{1112} & S_{1113}\\
S_{1121} & S_{1122} & S_{1123}\\
S_{1131} & S_{1132} & S_{1133}
\end &
\beginS_{1211} & S_{1212} & S_{1213}\\
S_{1221} & S_{1222} & S_{1223}\\
S_{1231} & S_{1232} & S_{1233}
\end &
\beginS_{1311} & S_{1312} & S_{1313}\\
S_{1321} & S_{1322} & S_{1323}\\
S_{1331} & S_{1332} & S_{1333}
\end \\
\beginS_{2111} & S_{2112} & S_{2113}\\
S_{2121} & S_{2122} & S_{2123}\\
S_{2131} & S_{2132} & S_{2133}
\end &
\beginS_{2211} & S_{2212} & S_{2213}\\
S_{2221} & S_{2222} & S_{2223}\\
S_{2231} & S_{2232} & S_{2233}
\end &
\beginS_{2311} & S_{2312} & S_{2313}\\
S_{2321} & S_{2322} & S_{2323}\\
S_{2331} & S_{2332} & S_{2333}
\end \\
\beginS_{3111} & S_{3112} & S_{3113}\\
S_{3121} & S_{3122} & S_{3123}\\
S_{3131} & S_{3132} & S_{3133}
\end &
\beginS_{3211} & S_{3212} & S_{3213}\\
S_{3221} & S_{3222} & S_{3223}\\
S_{3231} & S_{3232} & S_{3233}
\end &
\beginS_{3311} & S_{3312} & S_{3313}\\
S_{3321} & S_{3322} & S_{3323}\\
S_{3331} & S_{3332} & S_{3333}
\end
\end
\[
$$\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)*(4*Cc(1,1)+Cc(3,3) +4*Cc(1,3) - 4*Cc(6,6));
Gc_v = (1/30)*(2*Cc(1,1) + 2*Cc(3,3) - 4*Cc(1,3) + 10*Cc(6,6) + 12*Cc(4,4));
Cc_v_11 = Kc_v + (4/3)*Gc_v;
Cc_v_13 = Cc_v_11 - 2*Gc_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].
``` {.octave language="Octave"}
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.
``` {.octave language="Octave"}
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.
``` {.octave language="Octave"}
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\]](#maxwell){reference-type="ref" reference="maxwell"} are
ready, the effective stiffness can be found
``` {.c++ language="c++"}
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](#results){reference-type="ref" reference="results"} as a
comparison to the bench-marked data.
<figure id="results">
<figure>
<span class="image placeholder" data-original-image-src="MPD2.png"
data-original-image-title="" width="\textwidth"></span>
<figcaption>MPD = 2 MRD</figcaption>
</figure>
<figure>
<span class="image placeholder" data-original-image-src="MPD6.png"
data-original-image-title="" width="\textwidth"></span>
<figcaption>MPD = 6 MRD</figcaption>
</figure>
<figure>
<span class="image placeholder" data-original-image-src="MPD10.png"
data-original-image-title="" width="\textwidth"></span>
<figcaption>MPD = 10 MRD</figcaption>
</figure>
<figure>
<span class="image placeholder" data-original-image-src="MPD50.png"
data-original-image-title="" width="\textwidth"></span>
<figcaption>MPD = 50 MRD</figcaption>
</figure>
<figcaption>Comparison of reproduced results (dashed line) and the
bench-marked data (solid line)</figcaption>
</figure>
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.
\]