Abstract
A novel effective medium theory for homogenized transport coefficients of anisotropic mixtures of possibly anisotropic materials is developed. Existing theories for isotropic systems cannot be easily extended, because that would require geometric characterizations of anisotropic connectivity. In this work anisotropic connectivity is characterized by introducing a tensor that is constructed from a histogram of local percolating directions. The construction is inspired by local porosity theory. A large number of known and unknown generalized effective medium approximations for anisotropic media are obtained as limiting special cases from the new theory. Among these limiting cases the limit of strong cylindrical anisotropy is of particular interest. The parameter space of the generalized theory is explored, and the advanced results are applied to experiment.
Graphic abstract
Similar content being viewed by others
1 Introduction
Accessible and applicable characterizations of the irregular geometry and connectivity of disordered media are currently enjoying a surge of attention [1,2,3,4,5,6,7,8,9,10,11]. Disordered systems as diverse as nano-structured composites [12,13,14], randomly doped crystals [15], metallic alloys [16], sandstones [17], sedimentary rocks [18, 19], chalk [20], wood [21], soil [22], zeolites [23], polymer blends [24], porous silicon [25], or ceramics [26] all require quantitative characterizations of their connectivity as a prerequisite for predicting the transport of charge, volume, mass, momentum or energy in these materials.
Much attention was paid in recent years [1, 11] to the study of certain additive quantities [27] known as intrinsic volumes, quermaßintegrals [28] (p.240), valuations [29] (p.173) or Minkowski functionals [30, 31]. Additivity means that adding the values for two samples equals the sum of the value for the union plus that for the intersection [27], Eq. (42). In particular the additive Euler characteristic [27] (p.31), [28] (p.237) is often believed to determine transport in disordered media [32, 33] or to “account for connectivity in a macroscale way” [1](p.325). On the other hand it is widely appreciated, that transport in random systems depends upon the existence of percolating paths from inlet to outlet [34,35,36,37,38]. Rigorous relations between the Euler characteristic and percolation were established in Eqs. (48)–(52) in Ref. [39], p. 386. Except for certain special geometries these relations show that in general the Euler characteristic cannot possibly suffice to characterize connectivity or transport. Most importantly, the Euler characteristic is additive, while connectivity, connectedness or “percolativity” are not.
Disordered systems often exhibit anisotropic geometry and/or connectivity. Examples include metamaterials [40], composites [41,42,43] metallic filters [44], shales and clay [45, 46], wood [21, 47], layered porous media [48, 49], fibrous media [50, 51], piezoelectric polymer foams, photonic crystals [52] or liquid crystals. It is therefore important to generalize theories and geometric characterizations to the case of anisotropic media.
Geometric quantities for anisotropic media, especially vectorial [53] and tensorial [54, 55] extensions of additive functionals are currently under active investigation [56]. Let us emphasize that higher order Minkowski functionals as well as their tensorial extensions are studied primarily as descriptors of geometry and have, to the best of our knowledge, not been used as predictors for physical transport. Our intention in this work is to introduce novel effective medium approximations to predict percolation and conduction in disordered materials. Recall that classical effective medium theories [35, 57,58,59,60] contain only volume, more precisely volume density i.e. volume fraction, as a geometrical descriptor. In self-consistent effective medium theories for mixtures with infinite material contrast a percolation transition emerges at volume fraction 1/3. Anisotropic variants of effective medium theory based on volume fraction [41, 61, 62] exhibit the same percolation transition and they continue to attract much interest [42, 52, 63]. Minkowski functionals other than volume, however, have, as far as we know, not been used to generalize effective medium theories to anisotropic media, and the same applies to pair connectivity [38], the pair connectedness [64], or n-connectedness [65].
The self-consistent effective medium approximation for the effective dc-conductivity of homogeneous and isotropic media with porosity \({\overline{\phi }}\) reads [35, 57, 58]
where \(\upsigma _\mathbbm {p}^\prime \) and \(\upsigma _\mathbbm {m}^\prime \) are the dc-conductivities of two isotropic materials filling pores \(\mathbbm {p}\) and matrix \(\mathbbm {m}\) (they will be introduced more precisely in Eq. (7) below). For infinite contrast \(\upsigma _\mathbbm {p}^\prime /\upsigma _\mathbbm {m}^\prime =\infty \) this equation has a singularity at \({\overline{\phi }}={\overline{\phi }}_c=1/3\) which marks the percolation transition. Nonzero solutions for \({\overline{\upsigma }}^\prime \) exist only for \({\overline{\phi }}>{\overline{\phi }}_c\), while for \({{\overline{\phi }}\le 1/3}\) the effective conductivity \({\overline{\upsigma }}^\prime =0\) vanishes.
Equation (1) has been generalized to anisotropic media in the special case where \(\mathbbm {p}\) is a homogeneous mixture of aligned oblate ellipsoids \((x/b)^2+(y/b)^2+(z/c)^2\le 1\) with half axes \(b\ge c>0\). In this case the effective conductivity \({\overline{\upsigma }}^\prime \) becomes tensorial. Its eigenvalues \({\overline{\upsigma }}^\prime _x,{\overline{\upsigma }}^\prime _y,{\overline{\upsigma }}^\prime _z\) obey \({\overline{\upsigma }}^\prime _x={\overline{\upsigma }}^\prime _y\) and Eq. (1) becomes essentially two coupled equations [41], Eq. (8)
for \({\overline{\upsigma }}^\prime _x={\overline{\upsigma }}^\prime _y\) and \({\overline{\upsigma }}^\prime _z\). The depolarization factors \(N_a,N_c\) depend on a combined ratio R of the two unknowns \({\overline{\upsigma }}^\prime _x,{\overline{\upsigma }}^\prime _z\) and the two half axes b, c through the relations
that couple the two Eqs. (2a) and (2b). As in the isotropic case, these two coupled equations exhibit a percolation singularity at \({\overline{\phi }}={\overline{\phi }}_c=1/3\).
Local porosity theory [36] is a generalized self-consistent effective medium theory that contains local porosity distributions and local percolation probabilities above and beyond volume fractions to characterize the geometry and connectedness of a disordered medium. An important consequence are percolation thresholds \({\overline{\phi }}_c\ne 1/3\). Predictions from local porosity theory agree with observations such as the empirically known Archie’s law [36]. Local porosity distributions and local percolation probabilities have become standard tools for image analysis of porous media [66, 67]. They were applied to clay rocks [68], fractured coal [69], Opalinus clay [70, 71], chalk [72], reticulate porous ceramics [73], macroporous reticulated silicon oxycarbide [74], dielectric characterization of wood [21], highly porous polymeric matrices [75], metallic filters [44], beryllium pebbles and arctic firns [64], and recently even to Ibuprofen tablets in pharmaceutical research [76]. Numerous publications [17, 32, 77,78,79,80,81] have demonstrated the predictivity of percolation probabilities and the total fraction of percolating cells for physical transport.
2 Problem and objective
A shortcoming of local porosity theory for electrical transport in heterogeneous media (introduced in Ref. [36] and reviewed in Refs. [37, 82]) is that the local percolation probability \(\uplambda \), which characterizes the geometrical connectivity of the medium, is a scalar but not a tensorial quantity. Disordered or heterogeneous media, however, generally have tensorial transport properties even if the constituent materials are isotropic.
An attempt to extend LPT to anisotropic systems was made in Ref. [83]. The anisotropy is described using eigenvalues of the “orientation matrix”. The orientation matrix is constructed from the position vectors of those voxels within a measurement cell, that belong to the largest percolating cluster in the cell. The position vectors are computed relative to the “cluster center” defined as an arithmetic average. Two eigenvalue ratios of the orientation matrix are collected into local anisotropy distribution. The relation of the local anisotropy distribution to local percolation probabilities is not discussed in Ref. [83]. As a consequence the generalized effective medium equation is unchanged and transport quantities such as conductivities or permeabilities remain scalar quantities.
The objective of this work is to extend local porosity theory and effective medium theories to anisotropic media in such a way, that the transport quantities become tensorial. For brevity and simplicity we consider only transport of charge, i.e. the case of electrical transport.
3 Model
Let \(\mathbbm {p},\mathbbm {m}\subset \mathbbm {R}^d\) denote the domains (pore space, matrix) occupied by the two components of a heterogeneous or porous sample \(\mathbbm {S}=\mathbbm {p}\cup \mathbbm {m}\subset \mathbbm {R}^d\). The Maxwell equations [84] read
where \({\varvec{\mathrm {\mathrm {D}}}}({\varvec{\mathrm {r}}},t)\) is the electric displacement, \({\varvec{\mathrm {\mathrm {B}}}}({\varvec{\mathrm {r}}},t)\) is the magnetic induction, \({\varvec{\mathrm {\mathrm {E}}}}({\varvec{\mathrm {r}}},t)\) is the electric field, and \({\varvec{\mathrm {\mathrm {H}}}}({\varvec{\mathrm {r}}},t)\) is the magnetic field. Taking the divergence of Eq. (3d) gives the equation of continuity
It relates the macroscopic charge density \(Q({\varvec{\mathrm {r}}},t)\) and the macroscopic current density \({\varvec{\mathrm {\mathrm {J}}}}({\varvec{\mathrm {r}}},t)\).
Causality and locality in space restrict the constitutive equations to convolutions in time written as
for \({\varvec{\mathrm {r}}}\in \mathbbm {p}\cup \mathbbm {m}\), where \({\varvec{\mathsf {\upsigma }}}^\prime \) is the electrical conductivity tensor. The symbols \(\upepsilon _0,{\varvec{\mathsf {\upepsilon }}}_r\) denote the of dielectric permittivity of the vacuum resp. the permittivity tensor of the medium, and \(\upmu _0,{\varvec{\mathsf {\upmu }}}_r\) are the scalar resp. tensorial magnetic permeabilities of the vacuum resp. the medium with \(\upmu _0=4\pi \times 10^{-7}\)N/A\(^2\) and \(\upepsilon _0=1/(\upmu _0 c^2)\approx 8.8542\times 10^{-12}\)F/m. The conduction current due to free charges \({\varvec{\mathrm {\mathrm {J}}}}_c\) appearing in Ohm’s law, Eq. (4a), is related to the total current via the current of bound charges defined as \({\varvec{\mathrm {\mathrm {J}}}}_b({\varvec{\mathrm {r}}},t)={\varvec{\mathrm {\mathrm {J}}}}({\varvec{\mathrm {r}}},t)-{\varvec{\mathrm {\mathrm {J}}}}_c({\varvec{\mathrm {r}}},t)\). Fourier transformation with respect to t gives
where \(\omega =2\pi \nu \) is the angular frequency and \(\nu \) denotes frequency. One has
for the complex conductivity and dielectric tensors.
The material parameter functions for a composite material consisting of two homogeneous but anisotropic components are assumed to have the form
where
is the characteristic (or indicator) function of a subset \(\mathbbm {X}\subset \mathbbm {R}^d\). The relaxation frequencies, defined as
are assumed to exist. They can be used to make \(\omega \) dimensionless.
The discontinuity of the material parameters at the interface requires one to specify boundary conditions at \(\partial \mathbbm {p}=\partial \mathbbm {m}\). Mathematically the equations are interpreted in a weak sense as equations for distributions. Depending on the boundary conditions a suitable domain could be a Sobolev space for the potentials resp. a space of potential fields for the electromagnetic fields [85, 86]. The boundary conditions at the interface are
for \({\varvec{\mathrm {r}}}\in \partial \mathbbm {p}=\partial \mathbbm {m}\), where \(Q_{\partial }({\varvec{\mathrm {r}}},\omega )\) (resp. \({\varvec{\mathrm {\mathrm {J}}}}_{\partial }({\varvec{\mathrm {r}}},\omega )\)) are the Fourier transforms of (possibly time dependent) surface charge (resp. surface current) densities with support in \(\partial \mathbbm {p}=\partial \mathbbm {m}\). The notation \({\varvec{\mathrm {\mathrm {B}}}}_{\partial \mathbbm {p}}({\varvec{\mathrm {r}}},\omega )\) (resp. \({\varvec{\mathrm {\mathrm {B}}}}_{\partial \mathbbm {m}}({\varvec{\mathrm {r}}},\omega )\)) is the limiting value of the vector field as the point \({\varvec{\mathrm {r}}}\in \partial \mathbbm {p}=\partial \mathbbm {m}\) is approached from within \(\mathbbm {p}\) (resp. from within \(\mathbbm {m}\)).
It is assumed that there are neither volume nor surface charges or currents inside the medium so that \(Q=0\), \(Q_{\partial }=0\), \({\varvec{\mathrm {\mathrm {J}}}}=0\) and \({\varvec{\mathrm {\mathrm {J}}}}_{\partial }=0\) from now on. This assumption is commonly made to exclude phenomena and effects arising from electrical double layers and other physico-chemical processes in the interfacial region.
When the wavelength is large compared to the scale of heterogeneities the time derivatives in in Faraday’s law, Eq. (3), are small compared to spatial derivatives. Setting them to zero results in the quasistatic approximation valid for small frequencies
where \(\ell \) denotes the scale of heterogeneities and c is the speed of light. For heterogeneities with \(\ell \approx 10...100\mu \)m one finds \(\omega \ll 3...30\)THz approaching the infrared. The quasistatic approximation decouples the full equations (3). Fourier transformation of eqs. (3) with respect to t results in the quasistatic equations
where \({\varvec{\mathrm {r}}}\in \mathbbm {p}\cup \mathbbm {m}\) in Eqs. (12a)–(12c), \({\varvec{\mathrm {r}}}\in \partial \mathbbm {p}=\partial \mathbbm {m}\) in Eqs. (12d)–(12e), and where \(\upepsilon ({\varvec{\mathrm {r}}},\omega )\) denotes the complex frequency dependent local dielectric function.
The effective macroscopic dielectric tensor \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\) of self-consistent effective medium approximations [58] is usually defined by averaging Eq. (12c)
where the angular brackets \(\langle \cdot \rangle \) denote (spatial or ensemble) averaging and \({\varvec{\mathrm {\Pi }}}\) are parameters of physical importance (such as porosity or connectivity) that arise from averaging the smaller scale heterogeneities.
4 Local geometry
4.1 Disordered media
A disordered (or porous) medium \({\mathbbm {m}}\) (or \({\mathbbm {p}}\)) is assumed to be a random set. A random set is a set-valued random variable [87]. The probability space of random disordered geometries is denoted \((\Omega ,{\mathcal {A}},P)\). The set \(\Omega \) is the set of all closed subsets of \({\mathbbm {R}^d}\) including \(\emptyset \). The \(\sigma \)-algebra is generated by the system of closed sets having nonempty intersection with compact sets and the random geometry distribution P is the image measure of the measure on the probability space underlying the set-valued random variable. For mathematical details on random sets the interested reader is referred to Ref. [82].
4.2 Local porosity
Let \({\mathbbm {K}\subset \mathbbm {R}^d}\) be a convex and compact set with centroid at the origin \({{\varvec{\mathrm {0}}}\in \mathbbm {R}^d}\). The centroid is the geometrical center of \({\mathbbm {K}}\). Then
denotes its translate by a vector \({{\varvec{\mathrm {r}}}\in \mathbbm {R}^d}\). The local porosity in a measurement cell \({\mathbbm {K}({\varvec{\mathrm {r}}})}\) placed at position \({\varvec{\mathrm {r}}}\) is defined as
where
is the volume of a set \({\mathbbm {K}}\). The local matrix volume fraction is \({\phi _{\mathbbm {m}}=1-\phi }\).
4.3 Anisotropic local matrix fraction
Consider from now on \(d=3\) and let \({\mathbbm {m}_i}\) denote the local matrix space clusters defined as (path-)connectedness components of \({\mathbbm {m}\cap \mathbbm {K}({{\varvec{\mathrm {r}}}})}\). Then
where \(\mathbbm {m}_i\cap \mathbbm {m}_j=\emptyset \) for \(i\ne j\) and \(n_\mathbbm {m}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))\) is the total number of local clusters within \({\mathbbm {K}({{\varvec{\mathrm {r}}}})}\). The barycenters of these local clusters are
and
is their local volume fraction. The matrix is assumed to have a constant density taken to be unity. Averaging over the local clusters the tensor
is introduced, where \({\varvec{{\mathsf {X}}}}({{\varvec{\mathrm {r}}}})={{\varvec{\mathrm {r}}}}\otimes {{\varvec{\mathrm {r}}}}\), is the tensor product. The tensor is real, symmetric and non-negative definite. In the following the tensor is normalized as
by dividing with its trace. The non-negative eigenvalues \({J_{\mathbbm {m}a}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))}\), \({J_{\mathbbm {m}b}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))}\), \({J_{\mathbbm {m}c}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))}\) obey \(J_{\mathbbm {m}a}(\mathbbm {K}({{\varvec{\mathrm {r}}}})) +J_{\mathbbm {m}b}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))+J_{\mathbbm {m}c}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))=1\), and are assumed to be ordered as \({1\ge J_{\mathbbm {m}a}\ge J_{\mathbbm {m}b}\ge J_{\mathbbm {m}c}\ge 0}\).
4.4 Local percolation
Let \({\partial \mathbbm {K}({{\varvec{\mathrm {r}}}})}\) denote the boundary of a measurement cell with centroid \({{\varvec{\mathrm {r}}}}\). A measurement cell \({\mathbbm {K}({{\varvec{\mathrm {r}}}})={{\varvec{\mathrm {r}}}}+\mathbbm {K}}\) is called percolating if there exists a continuous path
with
where \({\mathbbm {B}({\varvec{\mathrm {x}}},a)=\{{{\varvec{\mathrm {r}}}}\in \mathbbm {R}^3\!:\!|{{\varvec{\mathrm {r}}}}-{\varvec{\mathrm {x}}}|\le a\}}\) denotes a ball of radius a centered at \({\varvec{\mathrm {x}}}\). Equation (21) is illustrated in Figure 1.
The set of local percolating directions for a measurement cell \({\mathbbm {K}({{\varvec{\mathrm {r}}}})}\) centered at \({{\varvec{\mathrm {r}}}}\) is defined as
where \({\mathbbm {S}^2\!=\!\partial \mathbbm {B}({\varvec{\mathrm {0}}},1)}\) is the unit sphere and \({\varvec{\mathrm {p}}}_0\rightsquigarrow {\varvec{\mathrm {p}}}_1\) means that \({\varvec{\mathrm {p}}}_0\) is the starting point and \({\varvec{\mathrm {p}}}_1\) is the end point of a percolation path in the sense of Eq. (21). The local percolation indicators, defined as
are a family of \(\{0,1\}\)-valued random variables
indexed by \({\mathbbm {K}}\). The image of the measure P under the map
defines the total percolation probability \({p(\mathbbm {K})}\) of the measurement cell \({\mathbbm {K}}\).
4.5 Anisotropic local percolation
The set \({\mathbbm {L}}\) of local percolating directions defines a tensor
with respect to the cell center \({{\varvec{\mathrm {r}}}}\). It is again normalized
by dividing with the trace. The non-negative eigenvalues \({J_{\mathbbm {p}a}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))}\), \({J_{\mathbbm {p}b}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))}\), \({J_{\mathbbm {p}c}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))}\) obey \({J_{\mathbbm {p}a}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))}\) \({+J_{\mathbbm {p}b}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))+J_{\mathbbm {p}c}(\mathbbm {K}({{\varvec{\mathrm {r}}}}))=1}\), and are assumed to be ordered as \({1\ge J_{\mathbbm {p}a}\ge J_{\mathbbm {p}b}\ge J_{\mathbbm {p}c}\ge 0}\).
4.6 Local geometry distribution
The local geometry of a random set \({\mathbbm {p}}\) inside a measurement cell \({\mathbbm {K}}\) is characterized by the geometric parameters \({g(\mathbbm {p};\mathbbm {K})=(\Lambda ,\phi ,d,e)}\) where
and
parametrize the anisotropy of the measurement cell \({\mathbbm {K}}\) with local porosity \(\phi \). These quantities form a family of \({\mathbbm {G}}\)-valued random variables
indexed by \({\mathbbm {K}}\) where \({\mathbbm {G}=\{0,1\}\times [0,1]^3 }\) and \({\mathcal {B}}\) is the standard Borel \(\sigma \)-algebra of the unit hypercube in \({\mathbbm {R}^4}\). The family of image measures
of P under g is the family of local geometry distributions with family index \({\mathbbm {K}}\).
In practical applications the distribution P of the random geometries is unknown. Instead one is given a sample assumed to be a realization of a random geometry. The local geometry distribution is then estimated as the joint empirical probability measure
where \({\varvec{\mathrm {r}}}_i\) is a sequence of cell centers such that all \({\mathbbm {K}({\varvec{\mathrm {r}}}_i)}\) are mutually disjoint, and N denotes the total number of different placements of \({\mathbbm {K}}\). Mathematically this empirical density is an estimator for a Young measure (see Refs. [39, 82, 89] for details).
4.7 Limiting local geometry distributions
The dependence of \({\mu (g;\mathbbm {K})}\) on \({\mathbbm {K}}\) disappears in the limit where \({\mathbbm {K}}\) shrinks to a point or expands to become \({\mathbbm {R}^3}\). In the limit \({|\mathbbm {K}|\rightarrow 0}\) of small pointlike measurement cells one has with Eq. (27)
where \({{\overline{\phi }}=\phi (\mathbbm {S})}\) is the total porosity of the sample.
In the opposite \({|\mathbbm {K}|\rightarrow \infty }\) limit the limiting local geometry distributions are also expected to become independent [89] of \({\mathbbm {K}}\). Then Eq. (27) gives
where \({\mu _1(\phi ,d_\mathbbm {p},e_\mathbbm {p})}\) is the conditional probability density that a cell is percolating in the scaling limit, has porosity \(\phi \) and percolation anisotropy parameters \({d_\mathbbm {p},e_\mathbbm {p}}\), while \({\mu _0(\phi ,d_\mathbbm {m},e_\mathbbm {m})}\) is the conditional probability density that a cell is non-percolating in the scaling limit, has porosity \(\phi \) and percolation anisotropy parameters \({d_\mathbbm {m},e_\mathbbm {m}}\). The quantity \(\uplambda (\phi )\) is the conditional probability that a very large measurement cell is percolating given that it has local porosity \(\phi \). It may be called limiting local percolation probability (cf. Ref. [36]) in the sense that
is expected to hold from Eq. (25).
The possible macroscopic limits have been classified in Ref. [89], Sec.VII.b, into four cases. In the following the macroscopically heterogeneous cases will be neglected and the porous medium is assumed to be macroscopically homogeneous and anisotropic. If it is macroscopically non-random one has
where \({\overline{\phi }}\) is the macroscopic porosity and \({{{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m}}\) are the macroscopic anisotropy parameters for pore and matrix. The simple choice
will be used below. Expectation values of measurable functions \({f:\mathbbm {G}\rightarrow \mathbbm {R}}\) are computed as
5 Coated ellipsoidal inclusion
5.1 Formulation of the local problem
Let \(d=3\) and let \({{\varvec{\mathrm {r}}}^{\!\mathrm {T}}=(x,y,z)\in \mathbbm {R}^3}\) denote a vector in cartesian coordinates. The potential problem for the unknown electrostatic potentials \({U_\mathbbm {c},U_\mathbbm {s},U_\mathbbm {u}}\) in the domains \({\mathbbm {c},\mathbbm {s},\mathbbm {u}\subset \mathbbm {R}^3}\) reads
where \({{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}}\) are the dielectric tensors of the homogeneous but anisotropic dielectric materials filling the domains \({\mathbbm {c},\mathbbm {s},\mathbbm {u}}\) and \({\varvec{\mathrm {\mathrm {E}}}}_{\mathrm {ext}}\) denotes a constant applied field. The union of \({\mathbbm {c}}\) and \({\mathbbm {s}}\) forms an ellipsoidal inclusion
consisting of an ellipsoidal core
surrounded by a shell
that is embedded into an environment (or background)
defined as the complement of \({\mathbbm {e}}\). It is assumed that the quadratic forms \({{\varvec{{\mathsf {Q}}}}_\mathbbm {e},{\varvec{{\mathsf {Q}}}}_\mathbbm {c}}\) are non-degenerate, and that \({\mathbbm {c}\subset \mathbbm {s}}\) and \({\mathbbm {s}\setminus \mathbbm {c}\ne \emptyset }\) holds throughout. The finite two-dimensional boundaries are denoted
and the vector fields \({{\varvec{\mathrm {\mathrm {n}}}}_{\partial \mathbbm {e}},{\varvec{\mathrm {\mathrm {n}}}}_{\partial \mathbbm {c}}}\) in Eq. (34) are unit normal vector fields on \({\mathbbm {e}}\) and \({\mathbbm {c}}\). The boundary conditions (34d), (34e) at \({\partial \mathbbm {u}\!\cap \!\partial \mathbbm {s}\!=\!\partial \mathbbm {e}}\) and (34f), (34g) at \({\partial \mathbbm {s}\cap \partial \mathbbm {c}\!=\!\partial \mathbbm {c}}\) express continuity of the potentials and the normal component of the dielectric displacement. The boundary condition (34h) demands that the potential vanishes at infinity \(|{{\varvec{\mathrm {r}}}}|\rightarrow \infty \).
The non-degenerate quadratic form matrices \({{\varvec{{\mathsf {Q}}}}_\mathbbm {e},{\varvec{{\mathsf {Q}}}}_\mathbbm {c}}\) (with \({\det {\varvec{{\mathsf {Q}}}}_\mathbbm {e}\ne 0, \det {\varvec{{\mathsf {Q}}}}_\mathbbm {c}\ne 0}\)) depend on the material properties \({{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}}\) of the shell as follows. Assume that \({\det {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}\ne 0}\) and define the coordinate transformation \({{{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {s}}}\) by the formula
as the square root of \({{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}}\). The linear, non-orthogonal and non-singular coordinate transformation (36b) relates the vector \({{\varvec{\mathrm {r}}}}=(x,y,z)^{\!\mathrm {T}}\) in the original coordinate system with the vector \({{{{\varvec{\mathrm {r}}}}_\mathbbm {s}}=({x_\mathbbm {s}},{y_\mathbbm {s}},{z_\mathbbm {s}})^{\!\mathrm {T}}}\) in the transformed coordinates. To specify the transformation uniquely it is assumed that
with
holds true. Here \({a_{\mathbbm {e},\mathbbm {s}}, b_{\mathbbm {e},\mathbbm {s}}, c_{\mathbbm {e},\mathbbm {s}}}\) denote the semiaxes of the outer shell ellipsoid \({\mathbbm {e}}\) and \({a_{\mathbbm {c},\mathbbm {s}}, b_{\mathbbm {c},\mathbbm {s}}, c_{\mathbbm {c},\mathbbm {s}}}\) for the core ellipsoid \({\mathbbm {c}}\) in \({({x_\mathbbm {s}},{y_\mathbbm {s}},{z_\mathbbm {s}})}\)-coordinates.
Finally, it is assumed that the ellipsoids \({\mathbbm {e}_{\mathbbm {s}}={{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {s}}^{-1}\mathbbm {e}}\) and \({\mathbbm {c}_{\mathbbm {s}}={{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {s}}^{-1}\mathbbm {c}}\) are confocal, i.e. it is assumed there exists a number \({0<C<c_{\mathbbm {e},\mathbbm {s}}^2}\) such that
holds. Then \({a_{\mathbbm {c},\mathbbm {s}}<a_{\mathbbm {e},\mathbbm {s}}}\), \({b_{\mathbbm {c},\mathbbm {s}}<b_{\mathbbm {e},\mathbbm {s}}}\) and \({c_{\mathbbm {c},\mathbbm {s}}<c_{\mathbbm {e},\mathbbm {s}}}\) and the coated ellipsoid is parameterized by four parameters. If the half axes \({a_{\mathbbm {e},\mathbbm {s}},b_{\mathbbm {e},\mathbbm {s}},c_{\mathbbm {e},\mathbbm {s}}}\) are normalized by division with their sum \({a_{\mathbbm {e},\mathbbm {s}}+b_{\mathbbm {e},\mathbbm {s}}+c_{\mathbbm {e},\mathbbm {s}}}\), then three parameters suffice to characterize the coated ellipsoid. They could be chosen as the volume fraction \({\phi _\mathbbm {c}}\), and two normalized half axes, or as \({\phi _\mathbbm {c}}\) and two of the axis ratios
for shell and core. Here the axis ratio \({\varphi _{\mathbbm {c}}}\) is a core to shell length fraction for the intermediate axes. The confocality condition (38) implies the relations
between these axis ratios. The core \({\mathbbm {c}}\) within the coated ellipsoidal inclusion \({\mathbbm {e}}\) has a local volume fraction
that can be used to eliminate \({\varphi _{\mathbbm {c}}}\) in exchange for \({\phi _\mathbbm {c}}\). Then with \({\phi _\mathbbm {c}\rightarrow 1}\) for \(C\rightarrow 0\), while \({\phi _\mathbbm {c}\rightarrow 0}\) for \({C\rightarrow c_{\mathbbm {e},\mathbbm {s}}^2}\). Given \({d_\mathbbm {s},e_\mathbbm {s}}\) and \({\phi _\mathbbm {c}}\) these equations are solved for \({\varphi _{\mathbbm {c}}}\) to determine
as functions of \({d_\mathbbm {s},e_\mathbbm {s}}\) and \({\phi _\mathbbm {c}}\). The functions \(D_\mathbbm {c},E_\mathbbm {c}:[0,1]^3\rightarrow [0,1]\) show clearly, that \({d_\mathbbm {c},e_\mathbbm {c}}\) are not independent parameters. Instead, they are fixed as soon as \({\phi _\mathbbm {c},d_\mathbbm {s},e_\mathbbm {s}}\) are given.
5.2 Solution of the local problem
To solve Eq. (34), one first transforms Eq. (34c) into the Laplace equation in vacuum using the transformation
where it is assumed that \({\det {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}\ne 0}\). To specify the linear, non-orthogonal and non-singular coordinate transformation \({{{{\varvec{\mathsf {\sqrt{\upepsilon }}}}}_\mathbbm {u}}}\) uniquely it is assumed that
holds true in \(({x_\mathbbm {u}},{y_\mathbbm {u}},{z_\mathbbm {u}})\)-coordinates with \(a_{\mathbbm {e},\mathbbm {u}} \ge b_{\mathbbm {e},\mathbbm {u}} \ge c_{\mathbbm {e},\mathbbm {u}} > 0\).
Following Refs. [90] and [91], Eq. (44), the solution of the boundary value problem (34) can be written as
where \({{\xi _\mathbbm {u}}={\xi _\mathbbm {u}}({{\varvec{\mathrm {r}}}})}\) (resp. \({{\xi _\mathbbm {s}}={\xi _\mathbbm {s}}({{\varvec{\mathrm {r}}}})}\)) is the first ellipsoidal coordinate of the vector \({{\varvec{\mathrm {r}}}}\) in \({({x_\mathbbm {u}},{y_\mathbbm {u}},{z_\mathbbm {u}})}\)-coordinates (resp. \({({x_\mathbbm {s}},{y_\mathbbm {s}},{z_\mathbbm {s}})}\)-coordinates) and \(C>0\) is the confocal constant in Eq. (38). The ellipsoidal coordinates \({({\xi _\mathbbm {u}},{\eta _\mathbbm {u}},{\zeta _\mathbbm {u}})}\) of a point \({({x_\mathbbm {u}},{y_\mathbbm {u}},{z_\mathbbm {u}})}\) are defined as the three roots \({\xi _\mathbbm {u}},{\eta _\mathbbm {u}},{\zeta _\mathbbm {u}}\) of the equation [92] (p. 19)
which obey
Similarly, the ellipsoidal coordinates \(({\xi _\mathbbm {s}},{\eta _\mathbbm {s}},{\zeta _\mathbbm {s}})\) of a point \(({x_\mathbbm {s}},{y_\mathbbm {s}},{z_\mathbbm {s}})\) are defined as the three roots of the equation
which obey
To formulate expressions for \({\varvec{{\mathsf {G}}}}_\mathbbm {c},{\varvec{{\mathsf {G}}}}_\mathbbm {s}({\xi _\mathbbm {s}}),{\varvec{{\mathsf {G}}}}_\mathbbm {u}({\xi _\mathbbm {u}})\) three diagonal tensors, \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {u}}}({\xi _\mathbbm {u}})\), \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}({\xi _\mathbbm {s}})\) and \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}}({\xi _\mathbbm {s}})\) of depolarizing factors are needed. All of them have unit trace, i.e. are normalized. The diagonal entries \(N_{\mathbbm {e},\mathbbm {u}}^{({x_\mathbbm {u}})}\), \(N_{\mathbbm {e},\mathbbm {u}}^{({y_\mathbbm {u}})}\), and \(N_{\mathbbm {e},\mathbbm {u}}^{({z_\mathbbm {u}})}\) of \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {u}}}({\xi _\mathbbm {u}})\) are
where \(w=a_{\mathbbm {e},\mathbbm {u}},b_{\mathbbm {e},\mathbbm {u}},c_{\mathbbm {e},\mathbbm {u}}\) for \(i={x_\mathbbm {u}},{y_\mathbbm {u}},{z_\mathbbm {u}}\) and \(|\mathbbm {e}|\) is the volume of \(\mathbbm {e}\). The diagonal entries \(N_{\mathbbm {c},\mathbbm {s}}^{({x_\mathbbm {s}})},N_{\mathbbm {c},\mathbbm {s}}^{({y_\mathbbm {s}})},N_{\mathbbm {c},\mathbbm {s}}^{({z_\mathbbm {s}})}\) of \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}({\xi _\mathbbm {s}})\) are
where \(|\mathbbm {c}|\) is the volume of \(\mathbbm {c}\) and \(w=a_{\mathbbm {c},\mathbbm {s}},b_{\mathbbm {c},\mathbbm {s}},c_{\mathbbm {c},\mathbbm {s}}\) for \(i={x_\mathbbm {s}}_\mathbbm {u},{y_\mathbbm {s}}_\mathbbm {u},{z_\mathbbm {s}}_\mathbbm {u}\). The diagonal entries \(N_{\mathbbm {e},\mathbbm {s}}^{({x_\mathbbm {s}})},N_{\mathbbm {e},\mathbbm {s}}^{({y_\mathbbm {s}})},N_{\mathbbm {e},\mathbbm {s}}^{({z_\mathbbm {s}})}\) of \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}}({\xi _\mathbbm {s}})\) are
where \(w=a_{\mathbbm {e},\mathbbm {s}},b_{\mathbbm {e},\mathbbm {s}},c_{\mathbbm {e},\mathbbm {s}}\) for \(i={x_\mathbbm {s}}_\mathbbm {u},{y_\mathbbm {s}}_\mathbbm {u},{z_\mathbbm {s}}_\mathbbm {u}\) and \(|\mathbbm {e}_\mathbbm {s}|=|\mathbbm {e}|\) is the volume of \(\mathbbm {e}\). Then
are the depolarization tensors in (x, y, z)-coordinates. With the abbreviations \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {u}}}={{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {u}}}(0)\), \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}={{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}(0)\), \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}={{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}(0)\), \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {c},\mathbbm {s}}}={{{\varvec{{\mathsf {L}}}}}_{\mathbbm {c},\mathbbm {s}}}(0)\) etc., one has
where \(\phi _\mathbbm {c}=|\mathbbm {c}|/|\mathbbm {e}|\) is the local core volume fraction and
is the dielectric difference between the core and the shell material.
6 Self-consistent approximation
The potentials \(U_\mathbbm {c}({{\varvec{\mathrm {r}}}};{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}})\), \(U_\mathbbm {s}({{\varvec{\mathrm {r}}}};{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}})\), and \(U_\mathbbm {u}({{\varvec{\mathrm {r}}}};{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}})\) in Eq. (45) depend on the data \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}\). The dielectric permittivity \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}\) of the coated ellipsoidal inclusion is determined as the solution of the equation
for all \({{\varvec{\mathrm {r}}}}\in \mathbbm {u}\) by demanding that the potential outside of the inclusion is identical to that of the full solution. Setting \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}}={{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}={{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}\) in Eq. (52) gives for the potential
and
for the fields. Equating Eq. (54b) with Eq. (52c) yields
where
and
Solving Eq. (56) for \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}\) gives
which simplifies to become
for the dielectric permittivity of the ellipsoidal inclusion in terms of its geometric and material parameters. The geometric parameters are restricted to \(\phi _\mathbbm {c},d_\mathbbm {s},e_\mathbbm {s}\) because \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {c},\mathbbm {s}}}\) follows from \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {s}}}\) and the latter depends only on \(d_\mathbbm {s},e_\mathbbm {s}\).
Following Ref. [58] and Eq. (13), the effective dielectric permittivity \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\) is defined by averaging Eq. (55) over the disorder to get
where \(\left\langle \dots \right\rangle \) is the disorder average. The implicit dependencies of \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}\) are clear from Eq. (60). The implicit dependencies of \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}\) emerge from noting that the quadratic form matrix \({\varvec{{\mathsf {Q}}}}_\mathbbm {e}\) of the ellipsoidal inclusion from Eq. (35a) can be expressed equivalently as
in \({{{\varvec{\mathrm {r}}}}_\mathbbm {u}}\)-coordinates or \({{{\varvec{\mathrm {r}}}}_\mathbbm {s}}\)-coordinates. Noting that Eq. (48) can be written as
allows to rewrite \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}\) as
and this shows that \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}={{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}},{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}})\) depends on \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}\) and \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}\).
Self-consistency determines \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\) and \({{\overline{{{\varvec{{\mathsf {L}}}}}}}}\) such that the background equals the effective medium. Setting
in Eq. (61b) yields
the self-consistent effective medium approximation for \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\). The effective depolarisation tensor \({{\overline{{{\varvec{{\mathsf {L}}}}}}}}\) defined as
with \({\varvec{{\mathsf {Q}}}}_\mathbbm {e}^{-1}\) given in Eq. (62b) depends on \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\). In addition to \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\) also the shape of the ellipsoidal inclusion \(\mathbbm {e}\) has to be determined self-consistently. This differs from conventional Bruggemann theory.
7 Local porosity theory
The self-consistent equation of local porosity theory is obtained by expressing \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}\) and \({{\overline{{{\varvec{{\mathsf {L}}}}}}}}\) in terms of local geometric quantities and replacing the disorder average with an average over the local geometry distribution \(\mu \). This gives
with \(\langle \dots \rangle _\mu \) given in Eq. (33). The dependencies
have to be related to the materials filling \(\mathbbm {p}\) and \(\mathbbm {m}\). The cases of \(\Lambda =0\) (nonpercolating local geometry) and \({\Lambda =1}\) (percolating local geometry) need to be distinguished, because for percolating local geometries the shell \(\mathbbm {s}\) should contain \(\mathbbm {p}\)-material, while it should contain \(\mathbbm {m}\)-material in the nonpercolating case. The core and shell properties are specified as
Inserting Eq. (69) into Eq. (68) and using (33) gives
as the mixing law of local porosity theory. The unknown quantity in Eq. (70) is \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\). The functional \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {e}}\) is given in Eq. (60) with Eqs. (57) and (52d), while \({{\overline{{{\varvec{{\mathsf {L}}}}}}}}\) is from eq.(67).
8 Analytical results
8.1 Small measurement cells
In the limit \(|\mathbbm {K}|\rightarrow 0\) of pointlike measurement cells the mixing law of local porosity theory becomes
after inserting Eq. (30) into Eq. (70). For isotropic materials \({{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}}=\upepsilon _\mathbbm {p}{\varvec{{\mathsf {1}}}}}, {{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}}=\upepsilon _\mathbbm {m}{\varvec{{\mathsf {1}}}}\) the classic Bruggemann effective medium approximation (EMA) equation
ensues, which should be compared to Eq. (1). The known anisotropic effective medium formula (2) is recovered setting \(\mu (1,\phi ,d_\mathbbm {p},e_\mathbbm {p})\,{=}\, {\overline{\phi }}\delta \Big [(\phi ,d_\mathbbm {p}, e_\mathbbm {p})\!-\!(1,\!1,\!c/b)\Big ]\) and \(\mu (0,\phi ,d_\mathbbm {m},e_\mathbbm {m})=(1\!-\!{\overline{\phi }}) \delta [(\phi ,d_\mathbbm {m},e_\mathbbm {m}) -\!(0,\!1,\!c/b)]\) in Eq. (30), where b, c are the parameters in Eq. (2). The limit of small measurement cells is the limit in which the basic assumption of statistical independence underlying effective medium approximations is violated. It is also worth to remark that, in a limiting sense, Eq. (2) requires measurement cells to remain anisotropic when they shrink to a point.
8.2 Large measurement cells
For large measurement cells \(|\mathbbm {K}|\rightarrow \infty \) Eq. (32) gives
For geometrical isotropy \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}={{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}}={\varvec{{\mathsf {1}}}}/3\) and material isotropy \({{{\varvec{\mathsf {\upepsilon }}}_\mathbbm {p}}=\upepsilon _\mathbbm {p}{\varvec{{\mathsf {1}}}}}\), \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {m}}=\upepsilon _\mathbbm {m}{\varvec{{\mathsf {1}}}}\) this becomes
the generalized effective medium approximation from Ref. [36], Sec.V.A, where
are the local dielectric functions.
8.3 Geometrical isotropy
Geometrical isotropy is the special case where the local geometry is spherical instead of ellipsoidal. Computing the integrals gives
and hence
Inserting this into Eq. (57) and (60) gives
for coated spheres. In this case Eq. (70) becomes
where the shorthand notation \({{\overline{{{\varvec{{\mathsf {L}}}}}}}}({\overline{{\varvec{\mathsf {\upepsilon }}}}},\,\cdot \,)={{\overline{{{\varvec{{\mathsf {L}}}}}}}}({\overline{{\varvec{\mathsf {\upepsilon }}}}},\;\cdot \;;1,1)\), was used and \(\mu (\phi )=\mu _1(\phi )=\mu _0(\phi )\) is the local porosity distribution. This generalizes local porosity theory to anisotropic materials.
For \(|\mathbbm {K}|\rightarrow 0\) (pointlike measurement cells) Eqs. (30), (57), (60) and (77) for spheres give again Eq. (71) as the generalized effective medium equation for \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\) in the case of anisotropic media.
9 Material isotropy
A situation of great practical and applied interest is the case of material isotropy but geometrical anisotropy. In this case the effective material parameters the effective dielectric properties become anisotropic only due to geometric properties of the mixture, while the materials themselves are isotropic. General results for macroscopically random systems will be investigated in the first Sect. 9.1 of this section. The non-random case will be discussed in more detail in the remaining subsections.
9.1 Macroscopically random systems
Isotropic materials occupying \(\mathbbm {p}\) and \(\mathbbm {m}\) have dielectric tensors proportional to \({\varvec{{\mathsf {1}}}}\)
where the scalar prefactors
are in general frequency dependent. Then \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {c}}=\upepsilon _\mathbbm {c}{\varvec{{\mathsf {1}}}}\), \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {s}}=\upepsilon _\mathbbm {s}{\varvec{{\mathsf {1}}}}\) so that
where \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}},{{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}\) and hence \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {s}}},{{{\varvec{{\mathsf {L}}}}}_{\mathbbm {c},\mathbbm {s}}}\) are all diagonal. It follows now from Eq. (62b) that \({\varvec{{\mathsf {Q}}}}_\mathbbm {e}^{-1}\) is diagonal and hence from Eq. (62a) that \({{\varvec{\mathsf {\upepsilon }}}_\mathbbm {u}}\) is diagonal, and thus that \({{{\varvec{{\mathsf {L}}}}}_{\mathbbm {e},\mathbbm {u}}}\) is diagonal. As a result
and \({{\overline{{{\varvec{{\mathsf {L}}}}}}}}\) are diagonal. Equation (60) becomes
where the diagonal matrices \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}},{{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}}\), given in Eqs. (49) and (50), can be written as
with
and \(D_\mathbbm {c},E_\mathbbm {c}\) from Eq. (42). Inserting Eq. (82) into Eq. (67), using
and substituting \(s = (b_{\mathbbm {e},\mathbbm {s}}^2 \upepsilon _\mathbbm {s}u)/{\overline{\upepsilon }}_y\) in (67) gives
with
where
are self-consistently determined parameters that depend upon the unknown \({\overline{\upepsilon }}\). Equation (68) becomes diagonal
with \({{\overline{{{\varvec{{\mathsf {N}}}}}}}}={{\overline{{{\varvec{{\mathsf {N}}}}}}}}(d_{{\overline{\upepsilon }}\mathbbm {s}},e_{{\overline{\upepsilon }}\mathbbm {s}})\) depending on \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\). This coupled system of three equations needs to be solved for \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\) with material data from (69).
Substituting \(u=e^2s\) in Eq. (85) gives a form that is useful when \(d=0\) or \(e=0\)
This form is convenient, because it allows one to read off the limiting cases
and
for all \(0\le d,e\le 1\).
Compared to the elliptic integrals in Eqs. (48), (49), (50) the parameters \(d_{{\overline{\upepsilon }}\mathbbm {s}},e_{{\overline{\upepsilon }}\mathbbm {s}}\) in Eq. (88b) are neither from the unit interval nor real, but in general complex numbers. As a consequence the well known formulae for real numbers given in Ref. [93] do not apply. The results for \(N_a,N_b,N_c\) with complex parameters d, e have been computed below in Appendix A.
9.2 Macroscopically non-random systems
For macroscopically homogeneous but non-random systems Eq. (32) yields the diagonal equation
for \({\overline{{\varvec{\mathsf {\upepsilon }}}}}\) where the notation
is used again for local functions. Written in components the diagonal Eq. (92) becomes three coupled equations (\(j=x,y,z\))
where the abbreviations
were introduced. The local functions are obtained from Eqs. (68b), (93), (83) as
where
are the diagonal elements in \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}}\), \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}\) from Eq. (84), and \(\upepsilon _\mathbbm {p}\), \(\upepsilon _\mathbbm {m}\) are from Eq. (80). Here
from Eq. (42). The material parameters \(\upepsilon _\mathbbm {p}^\prime ,\upepsilon _\mathbbm {m}^\prime ,\upsigma _\mathbbm {p}^\prime ,\upsigma _\mathbbm {m}^\prime \) in \(\upepsilon _\mathbbm {p},\upepsilon _\mathbbm {m}\) and the geometrical parameters \({\overline{\phi }},{\overline{\uplambda }},{{\overline{d}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {m}\) in these equations are measurable or known from experiment.
The complex effective depolarization factors for \(\mathbbm {p}\) and \(\mathbbm {m}\) depend on \({\overline{\upepsilon }}\). They couple the three equations (94) via Eq. (88b). Clearly, \({\overline{\upepsilon }}_j=0\) solves Eq. (94) for any j. For \({\overline{\upepsilon }}_j\ne 0\) one can rewrite the system of equations as
where the normalized effective dielectric permittivity
is the unknown,
is the right hand side, and
with \(j=x,y,z\) are parameters. Equation (97) can also be formulated as
This is not a solution of Eqs. (94) but a reformulation, because \(A_j,B_j\) depend through \(R_j\) on the two ratios \({\widehat{\upepsilon }}_x/{\widehat{\upepsilon }}_y\) and \({\widehat{\upepsilon }}_y/{\widehat{\upepsilon }}_z\) .
To obtain analytical information about solutions it is useful to consider limiting and special cases. The identity
holds generally. It is obtained by summation of Eq. (97a) with the help of Eq. (A13).
9.3 Spherical (Isotropic) Limit
The spherically isotropic limit is the limit \({{\overline{d}}}_\mathbbm {p}\rightarrow 1\), and \({{\overline{d}}}_\mathbbm {m}\!\rightarrow 1\), and \({{\overline{e}}}_\mathbbm {p}\rightarrow 1\), and \({{\overline{e}}}_\mathbbm {m}\!\rightarrow 1\). In this limit \({{{\varvec{{\mathsf {N}}}}}_{\mathbbm {e},\mathbbm {s}}}={{{\varvec{{\mathsf {N}}}}}_{\mathbbm {c},\mathbbm {s}}}={\varvec{{\mathsf {1}}}}/3\) and \({{\overline{N}}}_{\mathbbm {p} j}\!=\!{{\overline{N}}}_{\mathbbm {m} j}\). Inserting this into Eq. (97) renders \(P_j\), \(M_j\) independent of j and \(R_j\!=\!-{{\overline{N}}}_{\mathbbm {p} j}\) follows from Eq. (95a) and (95b). Subtracting Eq. (97a) with \(j=x\) from (97a) with \(j=y\) gives
where Eq. (97c) was used. The integral on the right hand side depends on \({\widehat{\upepsilon }}_z\), while the left hand side does not. Thus the equation can only be fulfilled for \({\widehat{\upepsilon }}_x={\widehat{\upepsilon }}_y\). Repeating the same subtraction with \(j=y\) and \(j=z\) gives \({\widehat{\upepsilon }}_y={\widehat{\upepsilon }}_z\), now invoking that the integral depends on \({\widehat{\upepsilon }}_x\). In this way the equality
and the generalized effective medium equation from Ref. [36] is recovered.
9.4 Planar (flat) limits
The planar or flat limits arise when \({{\overline{e}}}_\mathbbm {p}\rightarrow 0\) or \({{\overline{e}}}_\mathbbm {m}\rightarrow 0\). Then Eq. (41) implies \(E_\mathbbm {p}\rightarrow 0\) and \(E_\mathbbm {m}\rightarrow 0\) and Eq. (40) gives \(\varphi _{\mathbbm {c}}=1\), \({{\overline{d}}}_\mathbbm {p}^2=D_\mathbbm {p}^2\) and \({{\overline{d}}}_\mathbbm {m}^2=D_\mathbbm {m}^2\). From Eq. (96) follows \({N_C}_x,{N_C}_y,{N_B}_x,{N_B}_y\rightarrow 0\), \({N_C}_z\!\rightarrow \!{\overline{\phi }}\), \({N_B}_z\!\rightarrow \! 1-{\overline{\phi }}\). With this (97d) and (97e) give
and eqs. (97a) simplify to become
for \(j=x,y,z\). Note that these equations do not decouple. For \({{\overline{d}}}_\mathbbm {p}={{\overline{d}}}_\mathbbm {m}\) and \({{\overline{e}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {m}\) Eq. (103) is independent of \({\overline{\uplambda }}\) and hence the solution \({\widehat{\upepsilon }}\) becomes independent of \({\overline{\uplambda }}\).
9.5 Percolation limit
For sufficiently small frequencies \(\omega \ll \omega _\mathbbm {p}\), \(\omega \ll \omega _\mathbbm {m}\) the effective dielectric function \({\widehat{\upepsilon }}_j(\omega )\rightarrow {\overline{\upsigma }}^\prime (0)\) approaches the effective d.c. conductivity and the same holds for the material parameter functions \(\upepsilon _\mathbbm {p}(\omega )\rightarrow \upsigma _\mathbbm {p}^\prime (0)\) and \(\upepsilon _\mathbbm {m}(\omega )\rightarrow \upsigma _\mathbbm {m}^\prime (0)\). The complex quantity K approaches a non-negative real value
which will be called conductivity contrast. The d.c. conductivities obey \(\upsigma _\mathbbm {p}^\prime (0),\upsigma _\mathbbm {m}^\prime (0)\ge 0\). Without loss of generality it can be assumed that
holds true, i.e. that the material of higher conductivity occupies the pore space \(\mathbbm {p}\). Vanishing contrast \(\kappa =0\) means \(\upsigma _\mathbbm {p}^\prime =\upsigma _\mathbbm {m}^\prime \).
In the limit \(\kappa \rightarrow \infty \) of infinite conductivity contrast one has for \(j=x,y,z\)
in Eq. (97d) and (97e). Then Eq. (99) becomes
where
is the normalized d.c.-conductivity. Because \({\widehat{\upsigma }}^\prime _j\le P_j^\infty \), the left hand side is non-negative. For \({\overline{\uplambda }}\le 1/3\) the right hand side becomes negative, and hence \({\overline{\uplambda }}=1/3\) is a percolation threshold.
Because \(M_j\!=\! 0\), Eq. (98c) gives \(B_j=0\). Thus \({\widehat{\upsigma }}^\prime _j\!=\! 2A_j\) and
in the percolation limit. These are three coupled equations, that can be reduced to two coupled equations
by taking ratios. They have to be solved for the two unknown ratios \(u_{xy}\), \(u_{yz}\). The functions \(R_x,R_y,R_z\) depend only on these two ratios. For \({\overline{\uplambda }}=1\) and for \({\overline{\uplambda }}=0\) one finds
because the expression \(P_j^\infty {{\overline{N}}}_{\mathbbm {p} j}/({{\overline{N}}}_{\mathbbm {p} j}-1\)) resulting from Eq. (109) becomes negative.
9.6 Degenerate cases of Eq. (97)
Neither equation (97) nor Eq. (98) can be solved for \({\widehat{\upepsilon }}\) in full generality. To simplify the equations consider the special case \({{\overline{d}}}_\mathbbm {p}={{\overline{d}}}_\mathbbm {m}={{\overline{d}}}\) and \({{\overline{e}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {m}={{\overline{e}}}\), with \(0<{{\overline{d}}}\), \({{\overline{e}}}\le 1\). Then \({{\overline{N}}}_{\mathbbm {p} j}={{\overline{N}}}_{\mathbbm {m} j}={{\overline{N}}}_j\) and \(R_j=-{{\overline{N}}}_j\) ensues.
9.6.1 Spherical LPT
The spherical case \({{\overline{d}}}=1={{\overline{e}}}\) implies \({\overline{\upsigma }}^\prime _x={\overline{\upsigma }}^\prime _y = {\overline{\upsigma }}^\prime _z\) and \(P_x^\infty =P_y^\infty =P_z^\infty \) from Section 9.3. Then Eq. (107) recovers
the known results from Ref. [36], eqs.(5.1),(6.3) and Section 9.3 above.
Note, that \({\widehat{\upsigma }}^\prime \) depends on two geometric parameters, namely porosity \({\overline{\phi }}\) and connectivity \({\overline{\uplambda }}\) of the medium. Here \({\overline{\uplambda }}=\uplambda ({\overline{\phi }})\) was assumed to be constant. In general \({\overline{\uplambda }}({\overline{\phi }})\) is a nonlinear or singular function of \({\overline{\phi }}\) (see Ref. [36], Sec. VI.C).
9.6.2 Planar (flat) LPT
In the planar case with \({{\overline{e}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {m}={{\overline{e}}}\rightarrow 0\) the ellipsoids degenerate into flat ellipses. As discussed in Section 9.4 one has \({{\overline{d}}}_\mathbbm {p}^2=D_\mathbbm {p}^2\) in this case. Inserting Eq. (C7) into Eq. (96) and (106a) yields
for the local functions. Then Eq. (109) implies that \({\widehat{\upsigma }}^\prime _z =0\) for all \({\overline{\uplambda }}\). Furthermore, \({\widehat{\upsigma }}^\prime _x={\widehat{\upsigma }}^\prime _y={\overline{\phi }}\) must hold for \({\overline{\uplambda }}=1\), because of Eq. (111).
9.6.3 Oblate LPT
Here \({{\overline{d}}}_\mathbbm {p}={{\overline{d}}}_\mathbbm {m}=1\) and \({{\overline{e}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {m}={{\overline{e}}}\) with \({0<{{\overline{e}}}<1}\). Evaluating the integrals in Eq. (90) for \(d=1\) and using Eqs. (B3), (B4) yields Eq. (C3). Inserting Eq. (C3a) into Eq. (97) gives \(P_x=P_y\) and \(M_x=M_y\). Subtracting Eq. (97a) with \(j=y\) from (97a) with \(j=x\) yields
by an argument analogous to that in Section 9.3 for the spherically isotropic limit.
With this result Eq. (109) reduces to a system of two coupled equations
where \(u=\sqrt{{\widehat{\upsigma }}^\prime _y/{\widehat{\upsigma }}^\prime _z}\) generates the coupling, and the functions \(N_b(1,e)\), \(N_c(1,e)\) are given in Eq. (C3). They obey \(2N_b+N_c=1\). From Eq. (40) one finds
9.6.4 Prolate LPT
In this case \({{\overline{e}}}=1\), i.e. one has \({{\overline{e}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {m}=1\) and \({{\overline{d}}}_\mathbbm {p}={{\overline{d}}}_\mathbbm {m}={{\overline{d}}}\) with \(0<{{\overline{d}}}<1\).
Evaluating the integrals in Eq. (90) for \(e=1\) using Eq. (B3) and Eq. (B4) yields Eq. (C1). Inserting Eq. (C1c) into Eq. (97) shows \(P_y=P_z\) and \(M_y=M_z\). Subtracting Eq. (97a) with \(j=z\) from (97a) with \(j=y\) yields
by an argument analogous to that in Section 9.3 for the spherically isotropic limit.
With this result Eq. (109) reduces to a system of two coupled equations
where \(u=\sqrt{{\widehat{\upsigma }}^\prime _x/{\widehat{\upsigma }}^\prime _y}\) generates the coupling, and the functions \(N_a(d,1)\), \(N_b(d,1)\) are given in Eq. (C1). They obey \(N_a+2N_b=1\). From Eq. (40) one finds
For \({{\overline{d}}}\approx 1\) an approximate solution of Eq. (119) can be obtained for \(1/3\le {\overline{\uplambda }}\le 1\) by setting
with \(N_a,N_b\) from Eq. (C1). This decouples Eqs. (119) and yields \({\widehat{\upsigma }}^\prime \) as a function of \({\overline{\phi }},{\overline{\uplambda }},{{\overline{d}}}\).
For \({{\overline{d}}}\approx 0\) the prolate ellipsoids become strongly elongated and Eq. (C1) can be approximated by Eq. (C2) to lowest order in d. Approximating \(N_a({{\overline{d}}} u,1)\approx 0\) and \(N_b({{\overline{d}}} u,1)\approx 1/2\) in Eq. (119) by their limiting values gives
for \(1/2\le {\overline{\uplambda }}\le 1\). Then \({\widehat{\upsigma }}^\prime _j\ge 0\) implies \({\widehat{\upsigma }}^\prime _y={\widehat{\upsigma }}^\prime _z=0\) for \({\overline{\uplambda }}\le 1/2\). An anisotropic percolation threshold arises at \({\overline{\uplambda }}=1/2\). Inserting \({\widehat{\upsigma }}^\prime _y={\widehat{\upsigma }}^\prime _z=0\) into Eq. (107) implies
for \(1/3\le {\overline{\uplambda }}\le 1/2\). In this case, i.e. for \({{\overline{d}}}\approx 0\) an approximate solution of Eq. (119) can be obtained for \(1/3\le {\overline{\uplambda }}\le 1\) by using u from Eq. (122) with \(N_a,N_b\) from Eq. (C2).
9.6.5 Circocylindrical LPT
Here \({{\overline{d}}}_\mathbbm {p}={{\overline{d}}}_\mathbbm {m}={{\overline{d}}}=0\) and \({{\overline{e}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {m}={{\overline{e}}}=1\). This case resembles the spherical case. Inserting Eq. (91b) decouples Eq. (119) and yields
as a function of \({\overline{\phi }},{\overline{\uplambda }}\). Here the case \(1/3\le {\overline{\uplambda }}\le 1/2\) in Eq. (125a) follows from inserting \({\widehat{\upsigma }}^\prime _y = {\widehat{\upsigma }}^\prime _z=0\) in Eq. (107).
9.6.6 Ellipsocylindrical LPT
In this case \({{\overline{d}}}=0\) and \({{\overline{e}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {m}={{\overline{e}}}\) and \(0<{{\overline{e}}}<1\).
Evaluating the integrals in Eq. (90) for \(d=0\) gives eqs. (91f) resp. (C5). From Eq. (40) one finds
or
so that Eqs. (106) and (96) yield
for the local functions, and
holds true for their ratio. Dividing Eq. (109) for \(j=y\) by Eq. (109) for \(j=z\) results in the simple expression
for the ratio of the unknown effective conductivity components. Inserting this into Eq. (95a) gives the surprising result
for all \({{\overline{e}}}\), i.e. independent of \({\overline{\upsigma }}^\prime _j\) or \({{\overline{e}}}\). The result decouples eqs. (109) to give with the help of Eq. (106)
for \(1/2\le {\overline{\uplambda }}\le 1\) where \(E_\mathbbm {p}\) is from Eq. (127). For \({\overline{\uplambda }}\le 1/2\) it follows that \({\widehat{\upsigma }}^\prime _y={\widehat{\upsigma }}^\prime _z=0\) and inserting this into Eq. (107) gives
for \(1/3\le {\overline{\uplambda }}\le 1/2\).
10 Numerical results
10.1 Percolation transitions
The Bruggemann effective medium approximation (1) exhibits a percolation transition at porosity \({\overline{\phi }}_c=1/3\) at which the dc-conductivity \({\widehat{\upsigma }}^\prime (0)\) vanishes. While qualitatively correct, the percolation threshold is generally too high for direct application to porous and disordered media. They often become percolating at much lower volume fractions. Indeed, Archie’s empirical law [36], which relates electrical conductivity to porosity, suggests that there is no lower limit for the percolation threshold.
The Bruggemann result (1) has been generalized from spheres to oblate ellipsoids in Ref. [41]. It is governed by the generalized Eq. (2) which reduces to Eq. (1) for \(b=c\). The numerical solution of the anisotropic or ellipsoidal Bruggemann equation (2) is shown as dash-dotted lines in Fig. 2 as a function of \({\overline{\phi }}\). The isotropic or spherical Bruggemann result is plotted as the solid line Fig. 2. The percolation threshold is \({\overline{\phi }}_c=1/3\) in both cases.
The numerical solution of Eq. (109) for a comparable oblate case with \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p}) = (1.0,0.5)\) is shown as two dashed lines in Fig. 2. The long dashed line is \({\widehat{\upsigma }}^\prime _x = {\widehat{\upsigma }}^\prime _y\) and \({\widehat{\upsigma }}^\prime _z\) is short dashed. Here, \(\uplambda ({\overline{\phi }}) = {\overline{\phi }}^{1/3}\) has been assumed for the local percolation probability, and the remaining two axis ratios are \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})=(1.0,0.01)\). The percolation threshold \({\overline{\phi }}_c\) is then found as the solution of the equation \(\uplambda ({\overline{\phi }}) = 1/3\), which gives \({\overline{\phi }}_c = 1/27 \approx 3.7\%\). This is confirmed in Fig. 2. It is much more realistic than \({\overline{\phi }}_c = 1/3\) for natural media. The value of \({\overline{\phi }}_c\) depends on the local percolation probability \(\uplambda ({\overline{\phi }})\), thereby highlighting the importance of this geometrical quantity.
The solutions of Eq. (2) and Eq. (109) both display anisotropic \({\widehat{\upsigma }}^\prime \) conductivity. In both cases \({\widehat{\upsigma }}^\prime _x = {\widehat{\upsigma }}^\prime _y\), because the anisotropy is of oblate type. However, several differences exist between Eq. (2) and Eq. (109).
While \({\widehat{\upsigma }}^\prime _{x,\mathrm {Sch}}\) and \({\widehat{\upsigma }}^\prime _{z,\mathrm {Sch}}\) in Fig. 2 do not differ much, the difference between \({\widehat{\upsigma }}^\prime _{x}\) and \({\widehat{\upsigma }}^\prime _{z}\) is significant for \({\overline{\phi }} < 1\). This is in part due to the small value of \({{\overline{e}}}_\mathbbm {m}\). In addition \({\widehat{\upsigma }}^\prime _{z}\) does not approach unity for \({\overline{\phi }}\rightarrow 1\). This behaviour is caused by the local functions and stems from the confocality condition. It appears in the limit \({\overline{\phi }}\rightarrow 1\) whenever \(e_\mathbbm {p}\ne 1\) and is absent in the prolate case when \(e_\mathbbm {p}= 1\). Of course, the two additional axis ratios \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) in Eq. (109) allow for a greater variety of solutions than for Eq. (2). Recall also that Eq. (109) is only one of the many special cases from the generalized LPT in Eq. (70).
10.2 dc-Conductivity at infinite contrast
This section gives numerical solutions of the anisotropic effective medium equations (109) for the unknown normalized dc-conductivity \({\widehat{\upsigma }}^\prime (0)\). Equations (109) hold for material isotropy and geometric anisotropy in the percolation limit, i.e. at frequency \(\omega =0\) and infinite conductivity contrast \(\kappa =\infty \). The solutions for the unknown real part of the normalized effective conductivity \({\widehat{\upsigma }}^\prime \) are displayed in twelve subfigures of Fig. 3.They depend on six parameters \(({\overline{\phi }},{\overline{\uplambda }},{{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) with values in the six dimensional unit cube \([0,1]^6\). The porosity is fixed at \({\overline{\phi }}=0.1\) in all cases.
The central subfigure in Fig. 3 visualizes the four dimensional parameter space \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\in [0,1]^4\). The anisotropy parameters \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) for \(\mathbbm {p}\) are indicated by open symbols, the parameters \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) for \(\mathbbm {m}\) by filled symbols in the central unit square. Every parameter quadruple \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\in [0,1]^4\) is represented by a pair consisting of an open and a closed symbol. The two symbols of such a pair are connected by a dashed line.
Twelve subfigures in Fig. 3 are arranged around the central subfigure. They show \({\widehat{\upsigma }}^\prime =({\widehat{\upsigma }}^\prime _x,{\widehat{\upsigma }}^\prime _y,{\widehat{\upsigma }}^\prime _z)\) as functions of \({\overline{\uplambda }}\in [0,1]\). Because \({\overline{\phi }}=0.1\) the range of values is bounded by \({\widehat{\upsigma }}^\prime _j\le 0.1\) for \(j=x,y,z\). Note also that \({\widehat{\upsigma }}^\prime _x({\overline{\uplambda }})\ge {\widehat{\upsigma }}^\prime _y({\overline{\uplambda }})\ge {\widehat{\upsigma }}^\prime _z({\overline{\uplambda }})\) holds true for all \({\overline{\uplambda }}\) in all subfigures.
The twelve subfigures are labeled from (a) through (l) as follows
-
(a)
spherical - - - spherical
-
(b)
oblate - - - prolate
-
(c)
prolate - - - oblate
-
(d)
circular - - - circular
-
(e)
(nearly) prolate - - - flat
-
(f)
flat - - - (nearly) prolate
-
(g)
flat - - - flat
-
(h)
flat - - - circocylindrical
-
(i)
cylindrical - - - flat
-
(j)
circocylindrical - - - circocylindrical
-
(k)
cylindrical - - - none
-
(l)
none - - - cylindrical
The first word in e.g. “oblate - - - prolate” characterizes the position of the pair \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) in the unit square, while the second gives the position of the pair \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\).
10.2.1 \({{\overline{d}}}_\mathbbm {p}\approx {{\overline{d}}}_\mathbbm {m}\) and \({{\overline{e}}}_\mathbbm {p}\approx {{\overline{e}}}_\mathbbm {m}\)
The subfigures (a), (d), (g) and (j) located on the four corners of Fig. 3 illustrate cases with \({{\overline{d}}}_\mathbbm {p}\approx {{\overline{d}}}_\mathbbm {m}\) and \({{\overline{e}}}_\mathbbm {p}\approx {{\overline{e}}}_\mathbbm {m}\). In the isotropic case (a) the three components \({\widehat{\upsigma }}^\prime _x,{\widehat{\upsigma }}^\prime _y,{\widehat{\upsigma }}^\prime _z\) are nearly identical as expected from Eq. (101). In the circular - - - circular case (d) one has \({\widehat{\upsigma }}^\prime _x\approx {\widehat{\upsigma }}^\prime _y\) but \({\widehat{\upsigma }}^\prime _z \approx 0\). In both these cases the percolation transition occurs at \({\overline{\uplambda }}=1/3\). In the flat - - - flat case (g) a second percolation transition starts to appear in \({\widehat{\upsigma }}^\prime _x\) and \({\widehat{\upsigma }}^\prime _y\) at \({\overline{\uplambda }}=1/2\) while \({\widehat{\upsigma }}^\prime _z\approx 0\) remains small. Note also that \({\widehat{\upsigma }}^\prime _x\) and \({\widehat{\upsigma }}^\prime _y\) are nearly linear above \({\overline{\uplambda }}\approx 1/2\). In the circocylindrical - - - circocylindrical case (j) the percolation transition in \({\widehat{\upsigma }}^\prime _x\) at \({\overline{\uplambda }}=1/3\) and at \({\overline{\uplambda }}=1/2\) in \({\widehat{\upsigma }}^\prime _y\) and \({\widehat{\upsigma }}^\prime _z\) remain clearly visible. Now all components are nearly linear above \({\overline{\uplambda }}\approx 1/2\). The main difference to the flat - - - flat case (g) is that \({\widehat{\upsigma }}^\prime _z\ne 0\) and \({\widehat{\upsigma }}^\prime _x\ne {\widehat{\upsigma }}^\prime _y\) at \({\overline{\uplambda }}=1\). Other cases with \({{\overline{d}}}_\mathbbm {p}\approx {{\overline{d}}}_\mathbbm {m}\) and \({{\overline{e}}}_\mathbbm {p}\approx {{\overline{e}}}_\mathbbm {m}\) interpolate between these four limiting cases.
10.2.2 Flat pores \({{\overline{e}}}_\mathbbm {p}\approx 0\)
In the case \({{\overline{e}}}_\mathbbm {p}\approx 0\) the pores approach flat ellipses and the open symbol falls near the lower boundary of the central square. The corresponding subfigures are (d), (f), (g) and (h). As discussed in Section 9.6.2 one expects \({\widehat{\upsigma }}^\prime _z\approx 0\) for all \({\overline{\uplambda }}\), because \(P_z^\infty \rightarrow 0\) for \({{\overline{e}}}_\mathbbm {p}\approx 0\). Secondly one expects \({\widehat{\upsigma }}^\prime _x={\widehat{\upsigma }}^\prime _y=0.1\) at \({\overline{\uplambda }}=1\), because \(P_x^\infty =P_y^\infty ={\overline{\phi }}\) for \({{\overline{e}}}_\mathbbm {p}\approx 0\). Indeed, these expectations are confirmed in all cases.
10.2.3 Cylindrical pores \(d_\mathbbm {p}\approx 0\)
For \({{\overline{d}}}_\mathbbm {p}\approx 0\) the pores tend to become cylindrical and the open symbol \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) falls near the left boundary of the central square. This case is seen in subfigures (g), (i), (j) and (k). In all these cases the percolation transition in \({\widehat{\upsigma }}^\prime _x\) occurs at threshold 1/3, while \({\widehat{\upsigma }}^\prime _y\) approaches a percolation transition with threshold \({\overline{\uplambda }}=1/2\) as \({{\overline{d}}}_\mathbbm {p}\rightarrow 0\).
10.2.4 Exchanging \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) with \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\)
The parameter pairs \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) and \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) have very different impact on the result. This is illustrated in an approximate way by exchanging them.
Subfigures (b) and (c) illustrate such an exchange for parameters near the prolate and oblate boundary. Note that the curvature of \({\widehat{\upsigma }}^\prime _y\) changes.
Subfigures (e) and (f) illustrate the exchange between positions close to the flat boundary and positions somewhat removed from the upper prolate boundary. Again the curvature of \({\widehat{\upsigma }}^\prime _y\) changes, but in addition the change in \({\widehat{\upsigma }}^\prime _z\) is significant. In the case (f) the component \({\widehat{\upsigma }}^\prime _z\) is nearly zero for all \({\overline{\uplambda }}\), while it becomes strongly curved upwards in subfigure (e). This indicates a percolation transition at \({\overline{\uplambda }} = 1\) in the limit \({{\overline{e}}}_\mathbbm {m} \rightarrow 0\). Keep in mind, that the condition \({\overline{\uplambda }} = 1\) can be fulfilled for any \(0\le {\overline{\phi }}\le 1\).
The subfigures (h) and (i) illustrate the exchange between points near the flat and cylindrical boundary. Subfigure (h) shows that \({\widehat{\upsigma }}^\prime _z \approx 0\) holds for all \({\overline{\uplambda }}\). There is no percolation transition in \({\widehat{\upsigma }}^\prime _y\) at \({\overline{\uplambda }} = 1/2\) in (h). In (i) a percolation transition in \({\widehat{\upsigma }}^\prime _y\) occurs at \({\overline{\uplambda }} = 1/2\), while the strong curvature of \({\widehat{\upsigma }}^\prime _z\) indicates again a percolation transition at \({\overline{\uplambda }} = 1\) for \({{\overline{e}}}_\mathbbm {m} \rightarrow 0\).
The cylindrical cases with \({{\overline{d}}}_\mathbbm {p} \approx 0\) or \({{\overline{d}}}_\mathbbm {m} \approx 0\) are shown in subfigures (k) and (l). In (k) the percolation transition for \({\widehat{\upsigma }}^\prime _y\) and \({\widehat{\upsigma }}^\prime _z\) at \({\overline{\uplambda }} = 1/2\) is clearly visible. This transition does not appear in subfigure (l).
The eight subfigures illustrate the difference between the parameter pairs \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) and \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\). It is mainly the pair \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) that determines whether or not an additional percolation transition emerges at \({\overline{\uplambda }}=1/2\). It is also responsible for the smallness of \({\widehat{\upsigma }}^\prime _z\) in the limit \({{\overline{d}}}_\mathbbm {p} \rightarrow 0\) or \({{\overline{e}}}_\mathbbm {p} \rightarrow 0\).
10.2.5 Interpolation
The twelve results for \({\widehat{\upsigma }}^\prime \) shown in Fig. 3 depend in many cases continuously on the parameters \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m}, {{\overline{e}}}_\mathbbm {m})\). The continuous parameter dependence can be used to interpolate between different cases and to envisage results for other cases. This will be illustrated with two examples.
Consider first the sequence (d) \(\rightarrow \) (f) \(\rightarrow \) (h) \(\rightarrow \) (j). In the first transition the pair \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) moves from the lower boundary (flat) to the upper boundary (prolate) while \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) moves slightly away from the oblate boundary. The change in \({\widehat{\upsigma }}^\prime \) is small and mainly seen in \({\widehat{\upsigma }}^\prime _y\). In the transition (f) \(\rightarrow \) (h) both parameter pairs move towards the cylindrical boundary. The result is that the curvature in \({\widehat{\upsigma }}^\prime _y\) becomes much more pronounced and also \({\widehat{\upsigma }}^\prime _x\) becomes curved. The curvature in \({\widehat{\upsigma }}^\prime _x\) is due to the percolation transition that emerges at \({\overline{\uplambda }}=1/2\), because the open symbol has moved to the left. This is confirmed in the transition (h) \(\rightarrow \) (j) as the parameter pair \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) moves even closer to the cylindrical boundary. The change in \({\widehat{\upsigma }}^\prime _y\) on the other hand is caused by the closed symbol being close to the cylindrical boundary.
Next, consider the sequence (c) \(\rightarrow \) (e) \(\rightarrow \) (i) \(\rightarrow \) (l). In the first transition (c) \(\rightarrow \) (e) the pair \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) remains essentially fixed, while \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) moves from prolate to flat. As a result \({\widehat{\upsigma }}^\prime _z\) becomes strongly curved upwards, indicating a percolation transition at \({\overline{\uplambda }} = 1\). Shifting \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) close to the cylindrical boundary and \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) only slightly keeps the upward curvature in \({\widehat{\upsigma }}^\prime _z\) and exposes the emerging percolation transition at \({\overline{\uplambda }}=1/2\) in \({\widehat{\upsigma }}^\prime _x\) and \({\widehat{\upsigma }}^\prime _y\). The transition (i) \(\rightarrow \) (l) then undoes most of these changes and suggests to close the circle with an additional transition (l) \(\rightarrow \) (c).
10.3 ac-Conductivity at infinite contrast
Figure 4 shows numerical solutions of Eq. (94) for the normalized frequency dependent conductivity tensor \({\widehat{\upsigma }}^\prime (\omega )\) at infinite conductivity contrast \(\kappa = \infty \). The line styles for \({\widehat{\upsigma }}^\prime _x,{\widehat{\upsigma }}^\prime _y\) and \({\widehat{\upsigma }}^\prime _z\) in Fig. 4 are the same as in Fig. 3. Each of the subfigures in Fig. 4 has a fixed parameter quadruple \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p},{{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) and is labeled (a), (d), (g), (j), (e), and (i), to indicate the analogous pair \(({{\overline{d}}}_\mathbbm {p},{{\overline{e}}}_\mathbbm {p})\) - - - - \(({{\overline{d}}}_\mathbbm {m},{{\overline{e}}}_\mathbbm {m})\) labelled identically in the central subfigure of Fig. 3. The six subfigures in Fig. 4 are further subdivided into four subfigures. The two columns are distinguished by \({\overline{\uplambda }} = 0.35\) (left column) and \({\overline{\uplambda }} = 0.55\) (right column). The values are chosen in this way, because they are slightly above the two percolation thresholds at \({\overline{\uplambda }}=1/3\) and \({\overline{\uplambda }}=1/2\). In each subplot the upper two plots show \(\log _{10}({\widehat{\upsigma }}^\prime )\) against \(\log _{10}(\omega )\), the lower row displays the behaviour of the derivative \({\mathrm d} \log _{10}({\widehat{\upsigma }}^\prime (\omega ))/{\mathrm d}\log _{10}(\omega )\) against \(\log _{10}(\omega )\). The numerical value of the derivatives may be viewed as a local scaling exponent. The abscissa in all plots is \(\log _{10}(\omega )\). The frequency is dimensionless and given in units of the relaxation frequency \(\omega _\mathbbm {p}\) defined in Eq. (9).
The porosity is \({\overline{\phi }} = 0.1\) in all cases.
10.3.1 General observations
Two different sources can be identified which cause dispersion of \({\widehat{\upsigma }}^\prime \). The first source is the dispersion of the components \(P_j\) and \(M_j\), defined in Eq. (97d) and Eq. (97e). It is the dispersion in the response of a single local geometry. The second source causing dispersion arises from disorder and the percolation limit.
The spherical case (a) is shown for reference in the upper left corner of Fig. 4. In this case the dispersion is weak and it has a limited range in frequency.
The ordering \({\widehat{\upsigma }}^\prime _x(\omega ) \ge {\widehat{\upsigma }}^\prime _y(\omega ) \ge {\widehat{\upsigma }}^\prime _z(\omega )\) holds in all subfigures as for dc-conductivities.
In general, the dispersion decreases for increasing values of \({\overline{\uplambda }}\). Contrary to the limit \(\omega = 0\), the values of \({\widehat{\upsigma }}^\prime \) at high frequencies are barely affected by \({\overline{\uplambda }}\).
10.3.2 Flat pores \({{\overline{e}}}_\mathbbm {p} \approx 0\)
The case (d) from Fig. 3 is shown in the upper right corner of Fig. 4. Clearly, \({\widehat{\upsigma }}^\prime _z \ll {\widehat{\upsigma }}^\prime _x\) and \({\widehat{\upsigma }}^\prime _z \ll {\widehat{\upsigma }}^\prime _y\) holds for all frequencies. As in the case \(\omega = 0\), this results from \(P_z \rightarrow 0\) for \({{\overline{e}}}_\mathbbm {p} \approx 0\). For \({\widehat{\upsigma }}^\prime _x\) and \({\widehat{\upsigma }}^\prime _y\), there is one dispersive frequency range, which arises from the percolation transition at \({\overline{\uplambda }} = 1/3\). For \({\widehat{\upsigma }}^\prime _z\), there is an additional dispersion at higher frequencies. It arises from the dispersion of \(P_z\).
10.3.3 Cylindrical pores \({{\overline{d}}}_\mathbbm {p} \approx 0\)
The case (j) from Fig. 3 is placed on the right side of the middle row in Fig. 4. As discussed in section 10.2.3, \({\widehat{\upsigma }}^\prime _y\) and \({\widehat{\upsigma }}^\prime _z\) approach a percolation transition at \({\overline{\uplambda }} = 1/2\). This yields two effects for the dispersion of these components at \({\overline{\uplambda }} = 0.35\). First, the dispersion is much larger than the dispersion of \({\widehat{\upsigma }}^\prime _x\), and second, it occurs in a much broader frequency range. Both observations can be explained by the fact that the percolation limits at \({\overline{\uplambda }} = 1/3\) and \({\overline{\uplambda }} = 1/2\) lead to dispersion in two different frequency ranges which overlap. This effect cannot be observed at \({\overline{\uplambda }} = 0.55 > 1/2\).
10.3.4 Completely flat case
The case (g) from Fig. 3 is shown on the middle left side of Fig. 4. The dispersion of \({\widehat{\upsigma }}^\prime _x\) and \({\widehat{\upsigma }}^\prime _y\) is similar to case (j), because \({\widehat{\upsigma }}^\prime _x\) approaches the percolation limit at \({\overline{\uplambda }} = 1/3\) and \({\widehat{\upsigma }}^\prime _y\) at \({\overline{\uplambda }} = 1/2\). The behavior of \({\widehat{\upsigma }}^\prime _z\) can be qualitatively understood as a combination of the cylindrical (j) and the circular case (d). Compared to case (j), there is an additional dispersion at high frequencies which arises from the dispersion of \(P_z\) for \({{\overline{e}}}_\mathbbm {p} \approx 0\).
10.3.5 Flat matrix \({{\overline{e}}}_\mathbbm {m} \approx 0\)
The case (e) from Fig. 3 is shown in the lower left corner of Fig. 4. At first sight the dispersion looks similar to the flat case (g). However, the additional dispersion at high frequencies for \({\widehat{\upsigma }}^\prime _z\) results from the percolation limit at \({\overline{\uplambda }} = 1\), and not from the dispersion of \(P_z\) for \({{\overline{e}}}_\mathbbm {p} \approx 0\). This can be concluded from the fact that this dispersion barely changes between \({\overline{\uplambda }} = 0.35\) and \({\overline{\uplambda }} = 0.55\).
10.3.6 Mixed case \({{\overline{d}}}_\mathbbm {p} \approx 0\), \({{\overline{e}}}_\mathbbm {m} \approx 0\)
The case (i) from Fig. 3 is shown in the lower right corner of Fig. 4. Because \({{\overline{d}}}_\mathbbm {p} \approx 0\), the conductivity eigenvalue \({\widehat{\upsigma }}^\prime _y\) approaches the percolation limit at \({\overline{\uplambda }} = 1/2\), and \({\widehat{\upsigma }}^\prime _z\) approaches the percolation limit at \({\overline{\uplambda }} = 1\) because \({{\overline{e}}}_\mathbbm {m} \approx 0\) holds. This means that the three different eigenvalues of \({\widehat{\upsigma }}^\prime \) approach three different percolation limits. From this follows that the components are dispersive over different broad frequency ranges. At \({\overline{\uplambda }} = 0.35\), the dispersion at low frequencies results from the transition at \({\overline{\uplambda }} = 1/3\), the dispersion at intermediate frequencies from the transition at \({\overline{\uplambda }} = 1/2\), and the dispersion at high frequencies from the transition at \({\overline{\uplambda }} = 1\).
11 Application to experiment
A simple application of the new approximations to experiment can be made for anisotropic shaly quartz sands. It has recently been suggested by Nguyen, Vu, and Vu [94] that the electrical conductivities of shaly sands measured by Waxman and Smits [95] and Hill and Milburn [96] are anisotropic. Anisotropy ratios as large as 5 between normal and transversal conductivities are suggested in Nguyen, Vu, and Vu [94], Table 1. For details on the samples and their description the reader is referred to Waxman and Smits [95], p. 113. The electrical conductivity of the anisotropic shaly sands was measured using an impedance bridge for only one particular orientation. The data with an accuracy of 0.1 percent will be viewed here as resulting from randomly orienting the sample. Figure 5 shows the logarithm of the measured d.c. electrical conductivities normalized according to Eq. (108). Square symbols are the data reported in Table 1, while circles are from Table 3 in Waxman and Smits [95].
Detailed sample specific information for the anisotropy of pore \(\mathbbm {p}\) or matrix \(\mathbbm {m}\) in the shaly sands are not available from the publication [95]. To apply the theory nevertheless, it will be assumed for simplicity that \({{\overline{d}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {p}={{\overline{d}}}_\mathbbm {m}={{\overline{e}}}_\mathbbm {m}\) holds true. For the local percolation probability \(\uplambda ({\overline{\phi }})\) the function underlying Fig. 2 is assumed again, so that the percolation threshold occurs again at porosity \({\overline{\phi }}_c=1/27\approx 3.7\%\). The numerical solution of Eq. (109) yields the three eigenvalues \({\widehat{\upsigma }}^\prime _x({\overline{\phi }}), {\widehat{\upsigma }}^\prime _y({\overline{\phi }})\) and \({\widehat{\upsigma }}^\prime _z({\overline{\phi }})\) of the conductivity tensor as functions of average porosity. The three curves for \({{\overline{d}}}_\mathbbm {p}={{\overline{e}}}_\mathbbm {p}={{\overline{d}}}_\mathbbm {m}={{\overline{e}}}_\mathbbm {m}=0.58\) are displayed in Fig. 5.
The figure suggests that the scatter in the experimental data may indeed arise from an anisotropy of the samples. In conclusion, the theoretical assumptions underlying Eq. (109) and the experimental indications, that the measured conductivities of shaly sands are anisotropic, can be consistently quantified with the help of the approximations developed in this paper.
12 Conclusion
Generalized effective medium approximations for anisotropic media have been derived from local porosity theory. A purely geometric and experimentally measurable characterization of the anisotropy in (local) connectivities was introduced as a prerequisite for the generalized theory. The resulting anisotropic local porosity theory is based on the familiar local geometry distributions from the isotropic case. These distributions characterize the complex pore-space geometry. The generalized theory contains isotropic local porosity theory as a special case. The local geometry distributions characterizing local anisotropic connectivity are obtained from local percolating directions. Local percolating directions were defined based on the generalized percolation criterion introduced in Eq. (21). The geometric quantities enter directly into the generalized anisotropic self-consistent effective medium equation. No adjustable parameters appear in the generalized anisotropic local porosity theory. Several hitherto unknown generalized effective medium approximation formulae are reported. Among them the analytical results (125) and (132) for the cylindrical cases are of greatest interest. The generalized theory recovers previously known anisotropic effective medium approximations in the literature as special cases. The practical applicability of the theoretical predictions has been tested by a simple comparison with experiment.
Data Availability Statement
Not applicable.
References
R. Armstrong, J. McClure, V. Robins, Z. Liu, C. Arns, S. Schlüter, S. Berg, Trans. Porous Med. 130, 305 (2019)
S. Aramideh, P. Vlachos, A. Ardekani, Phys. Rev. E 98, 013104 (2018)
Z. Liu, J. McClure, R. Armstrong, Phys. Rev. E 98, 043102 (2018)
M. Blunt, Multiphase flow in permeable media (Cambridge University Press, Cambridge, 2017)
P. Tahmasebi, M. Sahimi, Phys. Rev. E 91, 032401 (2015)
K. Ding, Q. Teng, Z. Wang, X. He, J. Feng, Phys. Rev. E 97, 063304 (2018)
L. Zhu, C. Zhang, C. Zhang, X. Zhou, Z. Zhang, X. Nie, W. Liu, B. Zhu, Geofluids 2019, 7814180 (2019)
K. Alim, S. Parsa, D. Weitz, M. Brenner, Phys. Rev. Lett. 119, 144501 (2017)
L. Zhao, H. Li, D. Zhang, Phys. Rev. E 102, 023308 (2020)
H. Chen, X. He, Q. Teng, R. Sheriff, J. Feng, S. Xiong, Phys. Rev. E 101, 023305 (2020)
Y. Wang, J. Arns, S. Rahman, C. Arns, Phys. Rev. E 98, 043310 (2018)
B. Nøst, B. Hansen, E. Haslund, Phys. Scrip. T44, 67 (1992)
V. Shalaev, Phys. Rep. 272, 61 (1996)
C. Siclen, Phys. Rev. E 59, 2804 (1999)
P. Evesque, J. Physique 44, 1217 (1983)
S. Cho, Y. Cubides, H. Castaneda, Electrochim. Acta 236, 82 (1996)
P. Øren, S. Bakke, Trans. Porous Med. 46, 311 (2002)
P. Øren, S. Bakke, R. Held, Water Resour. Res. 43, W12S04 (2007)
M. Sahimi, Flow and transport in porous media and fractured rock, 2nd edn. (VCH Verlagsgesellschaft mbH, Weinheim, 2011)
M. Talukdar, O. Torsaeter, J. Pet. Sci. Eng. 33, 265 (2002)
K. Pentos, D. Luczycka, T. Wysoczanski, Wood Res. 62, 727 (2017)
H. Vogel, Eur. J. Soil Sci. 48, 365 (1997)
J. Perez-Ramirez, C. Christensen, K. Egeblad, C. Christensen, J. Groene, Chem. Soc. Rev. 37, 2530 (2008)
R. Hilfer, B. Nøst, E. Haslund, Th. Kautzsch, B. Virgin, B.D. Hansen, Phys. A 207, 19 (1994)
H. Lee, Y. Seo, D. Oh, K. Nahm, E. Suh, Y. Lee, H. Lee, Y. Hwang, K. Park, S. Chang, E. Lee, Appl. Phys. Lett. 62, 855 (1993)
K. Yoshida, J. Phys. Soc. Jpn. 59, 4087 (1990)
H. Hadwiger, Altes und Neues über konvexe Körper (Birkhäuser, Basel, 1955)
H. Hadwiger, Vorlesungen über Inhalt, Oberfläche und Isoperimetrie, Die Grundlehren der mathematischen Wissenschaften, vol. XCIII (Springer, Berlin, 1957)
R. Schneider, Convex bodies: The Brunn–Minkowski theory (Cambridge University Press, Cambridge, 1993)
R. Schneider, W. Weil, Stochastische Geometrie (Teubner, Stuttgart, 2000)
K. Mecke, D. Stoyan (eds.), Statistical Physics and Spatial Statistics, Lecture Notes in Physics, vol. 554. (Springer, Berlin, 2000)
S. Schlüter, H. Vogel, Advances in Water Resources 34, 314 (2011)
C. Scholz, F. Wirner, J. Götz, U. Rüde, G. Schröder-Turk, K. Mecke, C. Bechinger, Phys. Rev. Lett. 109, 264504 (2012)
S. Broadbent, J. Hammersley, Math. Proc. Camb. Philos. Soc. 53, 629 (1957)
S. Kirkpatrick, Rev. Mod. Phys. 45, 574 (1973)
R. Hilfer, Phys. Rev. B 44, 60 (1991)
R. Hilfer, Adv. Chem. Phys. XCII, 299 (1996)
D. Stauffer, A. Aharony, Introduction to Percolation theory (Taylor and Francis, London, 1992)
R. Hilfer, Trans. Porous Med. 46, 373 (2002)
N. Liu, H. Guo, L. Fu, S. Kaiser, H. Schweizer, H. Giessen, Nat. Mater. 7, 31 (2008)
L. Schwartz, Phys. A 207, 131 (1994)
L. Gao, J. Huang, K. Yu, Eur. Phys. J. B 36, 475 (2003)
M. Aouaichia, N. McCullen, C. Bowen, D. Almond, C. Budd, R. Bouamane, Eur. Phys. J. B 90, 39 (2017)
V. Heijtmanek, M. Vesely, P. Capek, J. Phys. Conf. Ser. 410, 012110 (2013)
S. Hemes, G. Desbois, J. Urai, B. Schröppel, J. Schwarz, Microp. Mesop. Mater. 208, 1 (2015)
P. Cosenza, D. Pret, A. Giraud, S. Hedana, Mech. Mater. 84, 55 (2015)
S. Avramidis, P. Englezos, T. Papathanasiou, AIChE J. 38, 1279 (1992)
Y. Tang, M. Aral, Water Resour. Res. 28, 1399 (1992)
M. Quintard, S. Whitaker, Trans. Porous Med. 5, 429 (1990)
J. Higdon, G. Ford, J. Fluid Mech. 308, 341 (1996)
F. Phelan-Jr., G. Wise, Composites Part A 27A, 25 (1996)
J. Flores-Mendez, B. Zenteno-Mateo, M. Moreno-Moreno, A. Morales-Sanchez, G. Minquiz, H. Vazquez-Leal, I. Vivaldo, S. Cortes-Lopez, A. Pinon-Reyes, R. Ambrosio, Eur. Phys. J. B 93, 124 (2020)
H. Hadwiger, R. Schneider, Elemente der Mathematik 26, 49 (1971)
S. Alesker, Geom. Funct. Anal. 8, 402 (1998)
D. Hug, R. Schneider, R. Schuster, St. Petersburg Math. J. 19, 137 (2008)
E. Vedel-Jensen, M. Kiderlen (eds.), Tensor Valuations and Their Applications in Stochastic Geometry and Imaging, Lecture Notes in Mathematics, vol. 2177. (Springer Verlag, Berlin, 2017)
D. Bruggeman, Annalen der Physik, 5.Folge 24, 636 (1935)
D. Bruggeman, Annalen der Physik, 5.Folge 24, 666 (1935)
R. Elliott, J. Krumhansl, P. Leath, Rev. Mod. Phys. 46, 465 (1974)
R. Landauer, in Electrical Transport and Optical Properties of Inhomogeneous Media, edited by J. Garland and D. Tanner (American Institute of Physics, New York, 1978) p. 2
B. Hornby, L. Schwartz, J. Hudson, Geophysics 59, 1570 (1994)
S. Berthier, J. Phys. I France 4, 303 (1994)
V. Levin, M. Markov, A. Mousatov, E. Kazatchenko, E. Peravgo, Eur. Phys. J. B 90, 192 (2017)
J. Ohser, C. Ferrero, O. Wirjadi, A. Kuznetsova, J. Düll, A. Rack, Int. J. Mat. Res. 103, 184 (2012). ((Formerly Z. Metallkd))
J. Ohser, K. Schladitz, 3D images of materials structures (Wiley, Weinheim, 2009)
X. Li, Q. Teng, Y. Zhao, S. Xiong, J. Feng, Phys. Rev. E 101, 053308 (2020)
P. Øren, S. Bakke, H. Rueslatten, “Digital core laboratory: Rock and flow properties derived from computer generated rocks,” (2006), International Symposium of the Society of Core Analysts held in Trondheim, Norway 12-16 September (2006)
P. Cosenza, D. Pret, M. Zamora, J. Geophys. Res. Solid Earth 120, 145 (2015)
Y. Jing, R. Armstrong, H. Ramandi, P. Mostaghimi, J. Geophys. Res. Solid Earth 122, 9849 (2017)
L. Keller, L. Holzer, P. Schuetz, P. Gasser, J. Geophys. Res. Solid Earth 118, 2799 (2013)
L. Keller, L. Holzer, J. Geophys. Res. Solid Earth 123, 285 (2017)
M. Talukdar, O. Torsaeter, M. Ioannidis, J. Howard, Trans. Porous Med. 48, 101 (2002)
J. Petrasch, F. Meier, H. Friess, A. Steinfeld, Int. J. Heat Fluid Flow 29, 315 (2008)
A. Choudhary, A. Agrawal, S. Pratihar, B. Singh, S. Behera, Adv. Eng. Mater. 21, 1900172 (2019)
A. Barbetta, C. Cametti, G. Rizzitelli, M. Dentini, Soft Matter. 8, 1120 (2012)
M. Neumann, R. Cabiscol, M. Osenberg, H. Markkötter, I. Manke, J. Finke, V. Schmidt, J. Microsc. 274, 102 (2019)
B. Biswal, C. Manwart, R. Hilfer, Phys. A 255, 221 (1998)
B. Biswal, C. Manwart, R. Hilfer, S. Bakke, P. Øren, Phys. A 273, 452 (1999)
K. Okabe, M. Blunt, J. Pet. Sci. Eng. 46, 121 (2005)
B. Biswal, P. Øren, R. Held, S. Bakke, R. Hilfer, Phys. Rev. E 75, 061303 (2007)
B. Biswal, R. Held, V. Khanna, J. Wang, R. Hilfer, Phys. Rev. E 80, 041301 (2009)
R. Hilfer, in Räumliche Statistik und Statistische Physik, Vol. Lecture Notes in Physics, Vol. 554, edited by D. Stoyan and K. Mecke (Springer, Berlin, 2000) p. 203
J. Liu, K. Regenauer-Lieb, C. Hines, K. Liu, O. Gaede, A. Squelch, Geochem. Geophys. Geosyst. 10, 1 (2009)
J. van Vleck, The theory of electric and magnetic susceptibilities (Clarendon, Oxford, 1932)
E. Sanchez-Palencia, Non-Homogeneous Media and Vibration Theory, Lecture Notes in Physics, vol. 127. (Springer Verlag, Berlin, 1980)
V. Jikov, S. Kozlov, O. Oleinik, Homogenization of differential operators and integral functionals (Springer, Berlin, 1994)
G. Matheron, Elements pour une Theorie des Milieux Poreux (Masson, Paris, 1967)
R. Hilfer, J. Hauskrecht, Trans. Porous Med. (2022) (in print)
R. Hilfer, J. Math. Phys. 59, 103511 (2018)
V. Frank, in Electrmagnetic Theory and Antennas, edited by E. Jordan (Mac Millan, New York, 1963) pp. 615-623
I. Lavrov, V. Yakovlev, Tech. Phys. 63, 1435 (2018)
L. Landau, E. Lifshitz, Theoretical Physics vol VIII, Electrodynamics of Continuous Media (Pergamon Press, Oxford, 1984)
J. Osborn, Phys. Rev. 67, 351 (1945)
S. Nguyen, M. Vu, M. Vu, J. Appl. Geophys. 123, 211 (2015)
M. Waxman, L. Smits, SPE J. 8, 107 (1968)
H. Hill, J. Milburn, AIME Pet. Trans. 207, 65 (1956)
M. Abramowitz, I. Stegun (eds.), Pocketbook of mathematical functions (abridged edition of handbook of mathematical functions) (Verlag Harri Deutsch, Thun, Frankfurt, 1984)
A. Hurwitz, Vorlesungen über allgmeine Funktionentheorie und elliptische Funktionen (Springer, Berlin, 1964) 4. vermehrte und verbesserte Auflage
Acknowledgements
The authors thank Tillmann Kleiner for discussions.
Funding
Open Access funding enabled and organized by Projekt DEAL.
Author information
Authors and Affiliations
Corresponding author
Appendices
Appendix A: Analytic continuation of \(N_a,N_b,N_c\)
To continue Eqs. (85) and (90) analytically into the complex plane let \(d,e\in \mathbbm {C}\) with \(d\ne e\) and \(d,e\notin \{0,1\}\). Assume \(|d|,|e|<\infty \) and define \(a,b,c\in \mathbbm {C}\) by
so that with
Equation (90) can be expressed as
in terms of three elliptic integrals. These can be written as [97] (p. 569)
where \(\mathrm {sn}\), \(\mathrm {sd}=\mathrm {sn}/\mathrm {dn}\) and \(\mathrm {tn}=\mathrm {sn}/\mathrm {cn}\) are Jacobian elliptic functions in Gudermann notation [98] (p.212). These integrals can be expressed in terms of the basic elliptic integrals of the first, second and third kind. They are independent of the path in the complex plane as long as the path stays within the fundamental period-parallelogram [98] (p.241). For each crossing of the boundaries of the period-parallelogram the corresponding period is added. The incomplete elliptic integral of the first kind is
and
is the incomplete elliptic integral of the second kind with complex amplitude \(\varphi (a,b,c)\), modulus k(a, b, c), and complementary modulus \({k^\prime }(a,b,c)\) defined by
The primitive functions are [97] (p. 576, 16.26)
where \(\mathrm {cd}(u)=\mathrm {cn}(u)/\mathrm {dn}(u)\). Using the special values
and the defining relations
one finds
where u and E are given in Eq. (A5). Summation in Eq. (A4) yields
with the help of identities for elliptic functions [97] (p. 569).
Returning to the parameters from Eq. (A1) the result becomes
where \(d,e\in \mathbbm {C}\) with \(d\ne e\) and \(d,e\notin \{0,1\}\),
and \(|d|,|e|<\infty \). Inserting the parameters from Eq. (A1) into Eq. (A10) shows that the sum of the depolarization factors is unity,
also for complex parameters d, e. For real parameters the same result can be seen more directly from the substitution \(s^2=(u+(1/d^{2}))(u+1)(u+e^2)\) in Eq. (85).
Appendix B: Some useful integrals
Appendix C: Special cases
1.1 1 Prolate spheroid
For prolate spheroids with \(a\ge b=c\) and \(e=1\) the depolarization factors are
where
with \(0<d<1\). For strongly elongated prolate ellipsoids with \(a\gg b=c\) and \(d\approx 0\) one finds
to lowest order.
1.2 2 Oblate spheroid
For oblate spheroids where \(a=b\ge c\) one has \(d=1\). In this case
where
with \(0<e<1\). When \(e\approx 0\) this becomes
as an approximation for very flat oblate spheroids.
1.3 3 Elliptical cylinder
The case \(d=0\) or \(a=\infty \) has depolarization factors
and corresponds to elliptical cylinders. For strongly elongated general ellipsoids with \(a\gg b\ge c\) one has \(d\approx 0\) and
where \(0<e<1\).
1.4 4 Flat ellipse
For flat ellipses where \(c=0\) one has \(e=0\) and
for all \(0<d<1\). For almost flat ellipses with \(e\approx 0\) the result is
where
are complete elliptic integrals of the first and second kind.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Hilfer, R., Hauskrecht, J. Effective transport coefficients of anisotropic disordered materials. Eur. Phys. J. B 95, 117 (2022). https://doi.org/10.1140/epjb/s10051-022-00338-5
Received:
Accepted:
Published:
DOI: https://doi.org/10.1140/epjb/s10051-022-00338-5