Rapid computation of effective conductivity of 2D composites by equivalent circuit and spectral methods

1.   Introduction

A key problem in materials science and engineering is representing a heterogeneous characteristic within a composite material with a single effective parameter. This has been variously called homogenization, coarse-graining, upscaling, or the effective medium approximation. Solutions have been developed in the context of a variety of linear response moduli: electronic, ionic, and thermal conductivity; mass diffusivity; dielectric constant; elastic modulus; hydraulic permeability, and viscosity. Early homogenization efforts were made by workers in a variety of fields including Maxwell, Faraday, Rayleigh, Poisson, Einstein, and Bruggeman ^[1]. Accurate methods of homogenization enable one to compute needed effective properties for the composite, but also lead to insights into the mechanisms controlling the effective property, thus aiding in material design. For instance, such methods can aid in the search for metamaterials or artificial composites that have unexpected properties ^[2,3].

Homogenization efforts make the common assumption that the property being studied is governed by the same physical laws at both the microscopic and macroscopic scales, which allows for averaging operations that generalize the known characteristics of the underlying homogeneous materials or phases ^[4,5]. This is a difficult problem because the approximation of an effective modulus depends on the moduli of the underlying pure materials, the relative amounts of those materials, and the physical arrangement or microstructure of the materials. Furthermore, these relationships can be nonlinear and depend on whether the underlying phases are isotropic or anisotropic. It is possible to develop methods for specific cases, such as a two-phase mixture of dilute inclusions in a matrix, but more general solutions are less forthcoming. One hindrance is that for the problems considered here that are governed by an elliptic partial differential equation (PDE), the solution in one region depends on conditions in every other region and thus the solution cannot be truly localized.

There is recent interest in homogenizing and simulating digitized media in which heterogeneous phases are specified on a grid ^[2,6,7,8,9]. Such can be tied to the growth in microstructure analysis by serial sectioning methods. Serial sectioning methods such as computed tomography (x-ray or neutron), nuclear magnetic resonance imaging, or focused-ion-beam milling allow one to experimentally determine 3D microstructures for composite materials in the form of a rectangular array of pixels or voxels. The microstructure array can be used as input for multi-physics computations, but some level of coarse-graining is typically required ^[10,11].

Direct numerical homogenization by solving the governing PDE for a conserved quantity is particularly challenging when phases have greatly contrasting moduli, sharply angled boundaries, or both, leading to singularities or rapid variations in currents or fluxes, as would be true for a rectangular grid. Efforts to solve this by decreasing the resolution of the mesh leads to matrix ill-conditioning and other slow-convergence problems. While there are specialized methods to address these problems ^{[7,9,12,13,14,15]}, there is no escaping the fact that substantial computational and user effort is required when there are rapid oscillations in conductivity in a medium, as with many important percolation problems.

When homogenizing large systems with more than two phases and without repeating symmetry, renormalization is an inexpensive numerical approach, provided a rapid calculation can be used for the coarse-grain procedure. Karim and Krabbenhoft ^[5] and others ^[16,17] suggest renormalization schemes as an alternative to direct numerical simulations that can be applied to square grid problems of multiple material phases. Renormalization in this case means a systematic recursive process of upscaling $2\times2$ blocks of squares, such that a grid of $2^{n}\times2^{n}$ squares (with $n$ a whole number) can be successively reduced to a single large square with a homogenized property. To upscale each block, Karim and Krabbenhoft apply a closed-form formula given by Mortola and Steffé ^[18] (and confirmed by others ^[4,19]) that is exact for assumed localized boundaries on the block. Such local confinement of the current or flux in each $2\times2$ block is not exact for a material with larger-scale heterogeneities or anisotropy, but nevertheless may be a good approximation.

In this paper electrical conductivity $\sigma$ is the representative modulus that will be homogenized, although obviously results can apply to other physical situations governed by a linear response modulus with the same tensorial character. We reinterpret and extend the homogenization formula first introduced by Mortola and Steffé ^[18] for four-phase square periodic grids (Section 2), to allow solutions for expanded geometries and anisotropic underlying phases, albeit approximately in some cases. This is done through equivalent electrical circuits applied to the $2\times2$ rectangular grid, leading to closed-form solutions (Section 3). In addition we provide a spectral method for solving the $2\times2$ grid problem that can likewise allow stretched, asymmetric, and rotated primitive cell geometries as illustrated in Figure 1 (Sections 4 and 5). The spectral method is a hybrid between fully analytic and numeric calculations and is more accurate than standard finite element methods for high contrast ratio of phase conductivities, but does suffer from some numerical limitations. Validations for the equivalent circuit and spectral methods are made by comparing to known exact solutions.

2.   Square lattice methods

The primitive unit cell considered here is the case of a $2\times2$ block with 4 distinct phases shown in Figure 2. This simple model can represent bulk material by use of periodic boundary conditions, which is equivalent to imposing unbalanced fixed potentials on two opposing edges and insulating conditions on the other two opposing edges of the block. Even if the larger system does not have this periodic symmetry, using the $2\times2$ primitive cell as the basis for renormalization can be used as a rapid means of homogenization as mentioned above ^[5]. However, the usefulness of such a renormalization scheme is expanded if it could work for any $n\times m$ size grid, rectangular grid elements, and anisotropic materials. The methods in this paper enable this.

The sharp corners and potentially highly contrasting phases make the grid system ill-suited for direct numerical calculations. Nevertheless, there are analytic solutions available for certain cases. A well-known case is when the four phases are squares and $\sigma_{1} = \sigma_{3}\ne\sigma_{2} = \sigma_{4}$. Such a two-phase alternating checkerboard can be homogenized by the exact solution $\sigma_{x}^{{ \square}} = \sigma_{y}^{{ \square}} = \sqrt{\sigma_{1}\sigma_{2}}$, where $\sigma_{i}^{{ \square}}$ indicates an effective conductivity in direction $i$ for a Cartesian grid system with periodic or mixed boundaries ^[20].

As Mortola and Steffé ^[18] first reasoned, there is an exact solution for the general four-phase square system of Figure 2. They proposed that effective conductivity is the geometric mean of lower ($\sigma_{x}^{L}$) and upper ($\sigma_{x}^{U}$) bounds that are themselves composed of alternating arithmetic and harmonic means of the quadrant conductivities:

The results for $\sigma_{y}^{{ \square}}$ can likewise be obtained by transposing $\sigma_{2}$ and $\sigma_{4}$ in the above formulas.

This conjecture was proven by Craster and Obnosov ^[19] and shortly thereafter by Milton ^[4] for square grids. Craster and Obnosov concurrently produced the exact solution for arbitrarily stretched rectangular grids. On further examination it becomes apparent that the lower and upper bounds given above are actually the limiting effective conductivities demonstrated much earlier by Maxwell for "stratified conductors" or laminates in which length in one direction greatly exceeds length in another, i.e., at infinite aspect ratio ^[21].

3.   Equivalent circuit methods

The lower and upper limit solutions for the $x$ and $y$ directions can be represented with the simple circuits shown in Figure 3. To generate the lower conductivity bounds (Eq 2.1) the central transverse element labeled $K_{x}^{*}$ or $K_{y}^{*}$ is given zero conductance; to generate the upper conductivity bounds (Eq 2.2) the transverse element is given infinite conductance. To complete the model, a closure relation is required to determine $K_{x}^{*}$ and $K_{y}^{*}$ for the intermediate cases when they are not 0 or $\infty$. These transverse elements must account for the 2D nature of current flow through the primitive cell.

Utilizing the equivalent circuit of Figure 3 allows for an extrapolation of the original exact formula to stretched geometries, unequal quadrant sizes, and anisotropic media. This flows naturally from using conductances, which contain geometrical information, rather than just phase conductivities. For instance quadrant 1 of Figure 3 is represented by the conductances (a) $K_{1y} = \sigma_{1}L/H$ and (b) $K_{1x} = \sigma_{1}H/L$ where $L$ is the quadrant length in $x$ and $H$ is the height in $y$ and an isotropic $\sigma_{1}$ is assumed for simplicity.

Effective or overall conductances for the ECs method are determined as follows. For Figure 3(b) we solve for the potential of the two interior nodes, $\phi_{12}$ and $\phi_{34}$, assuming an overall potential difference of 1 across the circuit. Conservation of current at the interior nodes leads to the matrix equation

The total current through the circuit is $I = K_{x} = K_{2x}\phi_{12}+K_{3x}\phi_{34}$. Combining these relations leads to an expression for the effective conductance:

The same process is repeated for Figure 3(a) to give

We now assert that $K_{x}^{*} = K_{x}$ and $K_{y}^{*} = K_{y}$, making it is possible to solve the two circuits simultaneously. This closure relation, in addition to being simple and physically intuitive, generates exact results in some cases and approximate though reasonably accurate results in other cases. The resulting closed-form formulas are

where

and $i$ represents the $x$ or $y$ direction. Considerable simplification is possible for cases with isotropic phases, geometric symmetry, or both. With $K_{x}$ and $K_{y}$ determined, they can be converted to effective $\sigma_{x}$ and $\sigma_{y}$ values through the geometry of the primitive cell.

The above model is exact (i.e., generates Eq 2.3) in the specific case of equally sized, isotropic square phases. The question is its accuracy in other instances. One canonical case is the stretched isotropic two-phase alternating checkerboard ($\sigma_{1} = \sigma_{3}\ne\sigma_{2} = \sigma_{4}$). Application of the ECs model produces the following effective conductivities in the $x$ and $y$ directions, where $R = L/H$ is the aspect ratio of individual quadrants and of the cell as a whole and $\chi = \sigma_{1}/\sigma_{2}$ is the contrast ratio.

When $R\rightarrow \infty$ this formula produces $\sigma_{y}\rightarrow2/(\sigma_{1}^{-1}+\sigma_{2}^{-1})$, and when $R\rightarrow0$ then $\sigma_{y}\rightarrow \tfrac{1}{2}(\sigma_{1}+\sigma_{2})$, demonstrating the exact Maxwell lamellar limits. Similarly, when $R = 1$ then $\sigma_{x} = \sigma_{y} = \sqrt{\sigma_{1}\sigma_{2}}$ as expected.

An alternative equivalent circuit network (ECn) was developed, as illustrated in Figure 4. The main idea is that in stretched geometries most of the current flow is 1D. Only in the vicinity of the intersection of the four phases is the current flow 2D. Therefore we handle the 2D flow region with the exact solution for a $2\times2$ square and append additional elements to the circuit to account for surrounding 1D current flow. Assembling the circuit is aided by first identifying the phase-quadrant of the stretched primitive cell with the minimum length or height given by $L_{\mathrm{min}}$. This defines the maximum size square which is centered on the four-phase intersection. This "divide and conquer" strategy allows the ECn method to be more accurate generally than the ECs method for stretched and asymmetric geometries of isotropic phases.

Applying the ECn model to the stretched isotropic two-phase alternating checkerboard leads to the following for effective conductivities

which can be compared to Eq 3.6 and similarly produces the correct limits for $R = 0, \, 1, \, \infty$. Additional comparison of the ECs and ECn models is given in Section 5.

Either of the ECs and ECn models, which for square geometry are exact, can be adapted to the $\tfrac{\pi}{4}$-rotated system (see Section 5.3) with the following simultaneous conductivity substitutions: $\sigma_{1}\leftarrow\sqrt{\sigma_{1}\sigma_{2}}$, $\sigma_{2}\leftarrow\sqrt{\sigma_{2}\sigma_{3}}$, $\sigma_{3}\leftarrow\sqrt{\sigma_{3}\sigma_{4}}$, $\sigma_{4}\leftarrow\sqrt{\sigma_{4}\sigma_{1}}$, $\sigma_{x}^{{ \square}}\leftarrow\sigma_{u}^{\mathrm{{ ECr}}}$, and $\sigma_{y}^{{ \square}}\leftarrow\sigma_{v}^{\mathrm{{ ECr}}}$ as illustrated in Figure 5. We call this the rotated equivalent circuit (ECr) model. These substitutions are by analogy with conductance through a corner region under high-contrast conditions ^[22]. This leads to a closed-form prediction of effective conductivities of the rotated square system:

This prediction is compared to other methods in Section 5.3.

4.   Spectral method

In addition to solving for the effective conductivity of the four-phase checkerboard, Craster and Obnosov solved for certain cases when the square checkerboard is rotated by $\tfrac{\pi}{4}$ radians before imposing the doubly periodic boundaries ^[23]. They additionally solved for a three-phase diamond geometry ^[24]. Here we develop a spectral 2D model, a hybrid of analytic and computational methods, which is able to accurately solve for the conductivity of many of these cases, as well as other types of tilings for which there are not yet analytic solutions. Thus this model can be used to evaluate proposed exact solutions and otherwise develop rapid and approximate solutions. In addition, we consider this method a stepping stone to the development of solutions to 3D checkerboard problems, and therefore explicitly avoid use of complex analysis to solve the 2D problem. Also note that our spectral method shares a few ideas with a method developed by Helsing for treating two-phase corner singularities numerically ^[25].

The basic procedure is to propose a general solution to the electric potential in the four-phase 2D composite, ensure that the solution satisfies a uniform set of internal boundary conditions, and then impose external boundary conditions, which can vary from problem to problem. Effective conductivity of the composite material is then derived from a ratio of average current density to average electric field.

4.1.   Solution to the internal boundaries

The spectral method starts with a proposed solution to the Laplace equation for potential $\phi$ in polar coordinates (Figure 6) by the standard technique of separation of variables:

The subscript $\alpha$ indicates a solution for quadrant 1, 2, 3, or 4. Note that all physical quantities are taken to be dimensionless. Coefficients $A_{\alpha}$ and $B_{\alpha}$ depend on quadrant while $\lambda$ is a common eigenvalue. The current density in the $\theta$ direction is given by

We enforce the internal boundaries between contiguous quadrants in order to determine $A_{\alpha}$, $B_{\alpha}$, and $\lambda$ in terms of conductivities $\sigma_{\alpha}$. At each boundary both $\phi$ and $i_{\theta}$ are continuous, independent of $r$ value. (Note that instead of potential we could have equivalently specified that the tangential electric field, in this case $-\partial\phi_{\alpha}/\partial r$, is continuous across the boundary.) The resulting 8 equations can be expressed in matrix form:

where $\mathbf{A}$ and $\mathbf{B}$ are respective column vectors of $A_{\alpha}$ and $B_{\alpha}$ values, stacked to become a single vector (i.e., a block matrix),

and for compactness $c_{\alpha} = \cos\left(\alpha\frac{\pi}{2}\lambda\right)$ and $s_{\alpha} = \sin\left(\alpha\frac{\pi}{2}\lambda\right)$.

Equation 4.3 requires that

and this relationship is used to solve for $\lambda$. But first we reduce the number of unknowns by the substitutions $c_{\alpha} = T_{\alpha}(c_{1})$ and $s_{\alpha} = s_{1}U_{\alpha-1}(c_{1})$, where $T_{\alpha}$ and $U_{\alpha}$ are respectively Chebyshev polynomials of the first and second kind and order $\alpha$. These relationships come from the multiple angle formulas for sine and cosine. The identity $c_{1}^{2} = 1-s_{1}^{2}$ is also used.

After these substitutions and considerable algebra we obtain from Eq 4.5 that $s_{1}^{2} = \varsigma^{2}$ (the primary solution) and $s_{1} = 0$ (the secondary solution), where $\varsigma$ is a function of the quadrant conductivities:

From $\varsigma$ we can determine an infinite series of eigenvalues for the primary solution using

An important element of the solution is that the eigenvalues come in pairs as discussed below. For instance the secondary solution to Eq 4.5 ($s_{1} = 0$) produces a series of paired eigenvalues $\lambda_{n}$, namely repeated even numbers (i.e., 2, 2, 4, 4, 6, 6, $\dots$). The eigenvalue $\lambda = 0$ is omitted from the solution due to our intention to keep $\phi = 0$ at the intersection of the four quadrants.

An infinite series of eigenvalues means that we can use a superposition of the proposed solutions from Eq 4.1, namely

where for generality and compactness we have combined the primary and secondary solutions into a single set of spectral coefficients $C_{n}$ and eigenvalues $\lambda_{n}$. Nevertheless, for many systems of interest the coefficients corresponding to the secondary solution are zero (within numerical uncertainties) and therefore may be neglected. The term $(-1)^{n+1}$ in the equation accounts for the alternating sign of $c_{1}$ relative to $s_{1}$ for the sequence of eigenvalues, which in turn changes the sign of $\mathbf{B}$ relative to $\mathbf{A}$.

In order to calculate coefficient vectors $\mathbf{A}$ and $\mathbf{B}$, we consider below three cases: (1) $s_{1}^{2} = 0$, (2) $s_{1}^{2} = 1$, and (3) $0 < s_{1}^{2} < 1$.

Case (1) occurs for the secondary solution and for the sake of generality needs to be considered as coexisting with either case (2) or (3). While case (1) could also occur for the primary solution when $\sigma_{1} = \sigma_{3} = 0$ or $\sigma_{2} = \sigma_{4} = 0$, zero conductivities for diagonally opposite quadrants would result in zero effective conductivity for the system and need not be considered further. For case (1) the internal boundary conditions establish ratios among elements of $\mathbf{A}$ and $\mathbf{B}$ respectively and allow us one degree of freedom for each vector. For instance, one could choose

The $A_{\alpha}$ and $B_{\alpha}$ terms in Eq 4.8 are further scaled and maintain their respective independence by virtue of the $C_{n}$ coefficients, the pairs of duplicate eigenvalues, and the alternating sign of the $(-1)^{n}$ term.

Case (2) occurs when $\sigma_{1}\sigma_{3} = \sigma_{2}\sigma_{4}$. It is similar to case (1) because, with $s_{1}^{2} = 1$, Eq 4.7 produces a series of duplicated eigenvalues, namely the odd values $1, 1, 3, 3, 5, 5, \dots$. The internal boundary conditions again require certain ratios among elements of $\mathbf{A}$ and $\mathbf{B}$ respectively and allow independent scaling of $\mathbf{A}$ and $\mathbf{B}$. For instance, one could choose

Case (3) is the most general where $0 < s_{1}^{2} < 1$. Matrix $\mathbf{G}$ (Eq 4.4) then has rank 7, meaning there is one degree of freedom and that $\mathbf{A}$ and $\mathbf{B}$ are coupled. We therefore stipulate that $A_{1} = 1$. This can be incorporated into our system of equations by use of vector $\mathbf{g} = \left[1\, \, 0\, \, 0\, \, 0\, \, 0\, \, 0\, \, 0\, \, 0\right]~^{\mathrm{T}}$ that can augment $\mathbf{G}$ to create an invertible matrix. $\mathbf{A}$ and $\mathbf{B}$ can therefore be obtained from

where $\mathbf{G}$ is evaluated at $s_{1} = \varsigma$ and $c_{1} = \sqrt{1-\varsigma^{2}}$. There is an algebraic solution to this equation, but it is somewhat lengthy and is not given here. A numerical solution by a standard matrix equation solver is sufficiently accurate. Again, the treatment as if $c_{1}$ and $s_{1}$ always have the same sign is corrected by the $(-1)^{n}$ term in Eq 4.8.

Before proceeding with the spectral method, it is worth making a few connections to prior work. The quantities $\varsigma$ and $\lambda_{1}$ defined above play a central role in the effective conductivity of four-phase checkerboards, as shown by comparing this work to the results of Craster and Obnosov derived from complex analysis ^[19]. They use an intermediate formula of the form $\cos(\pi\lambda) = 1-2\Delta^{2}$. It can be shown that their eigenvalue $\lambda$ is related to our $\lambda_{1}$ by $\lambda_{1} = 1-\lambda$ and that $\varsigma^{2} = 1-\Delta^{2}$. The analytic solution to the effective conductivities $\sigma_{x}^{{ \square}}$ and $\sigma_{y}^{{ \square}}$ of the four-phase square checkerboard (i.e., Eq 2.3) can be expressed in terms of $\varsigma$:

Furthermore, in the general case where the four-phase system is stretched and aspect ratio $R$ ranges from 0 to $\infty$, the ratio of effective conductivity for the stretched system relative to the square system (for either direction $x$ or $y$) is bounded by the values of $\varsigma$ and its inverse:

or equivalently $\sigma_{i}^{L}/\sigma_{i}^{U} = \varsigma^{2}$ for either direction, where $\sigma_{i}^{L}$ and $\sigma_{i}^{U}$ are the Maxwell lamellar limits (Eqs 2.1 and 2.2). Eq 4.13 demonstrates for the case $\varsigma = 1$ (i.e., $\sigma_{1}\sigma_{3} = \sigma_{2}\sigma_{4}$) that effective conductivity of the system is independent of aspect ratio and is isotropic (i.e., $\sigma_{x} = \sigma_{y}$).

4.2.   Geometry of the external boundary

Up to this point in the development of the spectral method, no external boundary was employed, meaning the solution is general and can in principle be a basis for a periodic or nonperiodic system; a square or rectangular (stretched) system; a non-rotated or rotated system; or any other system with a regular or irregular external boundary. Spectral coefficients $C_{n}$ in Eq 4.8 are determined numerically to satisfy the arbitrary external boundary with minimal error. The analytic solution up to this point ensures that regardless of the finite number of terms used in the series or the $C_{n}$ values chosen, exact continuities in potential and normal current at the internal boundaries are maintained. This allows the solution to retain accuracy under extreme values of contrast ratio (ratio of conductivities for adjacent phases) far beyond standard or even some specialized finite element formulations ^[7].

Figure 7 shows the generalized geometry we assume for the external boundary. A bounding rectangle is constructed relative to the origin from the distances $\xi$, $L-\xi$, $\psi$, and $H-\psi$. Angles at the rectangle vertices are given by $\theta_{1} = \mathrm{acot}(\xi/\psi)$, $\theta_{2} = \mathrm{acot}[(\xi-L)/\psi]$, $\theta_{3} = \pi+\mathrm{acot}[(L-\xi)/(H-\psi)]$, and $\theta_{4} = \pi+\mathrm{acot}[\xi/(\psi-H)]$. The bounding rectangle with its associated $uv$ coordinates are rotated by angle $\beta$ relative to the $xy$ coordinates. Coordinate $\theta$ is zero at the positive $u$ semi-axis, meaning $\theta$ is tied to the rotated external boundary. Internal phase boundaries are not rotated and remain aligned with the $xy$ axes. Thus, depending on the value of $\beta$, phase boundaries may or may not align with the $uv$ axes or the four $\theta_{i}$ values, though in this work only multiples of $\pi/4$ were used for $\beta$.

Effective conductivity depends on relative geometry of the phases, rather than on absolute size, therefore we can uniformly scale the distances as needed. For the sake of reasonable accuracy under finite machine arithmetic, we take the above distances to be normalized so that $L+H = 2\sqrt{2}$.

Our intention is to define a boundary equation completely in terms of $\theta.$ Therefore we use the following piecewise function for the radius at the perimeter

A second useful function tracks the accumulated linear distance along the perimeter from $\theta = 0$ to $2\pi$:

The potential along the rotated perimeter is given by an application of Eq 4.8 to get

Perimeter eigenpotentials are given by

where $\theta^{\bullet} = (\theta+\beta)\mod\mathrm{2\pi}$ is the angle expressed in the unrotated $xy$ coordinates and where for convenience we have turned the $A$ and $B$ coefficients into piecewise functions of angle, e.g.,

To properly formulate external boundary conditions in terms of fields and fluxes, we need to define the gradient of $\phi$ with respect to the $uv$ coordinates and evaluate it along the perimeter:

where

With gradients defined we can also calculate eigencurrents on the perimeter in the $u$ and $v$ directions:

where $\sigma(\theta^{\bullet})$ is defined in the same manner as Eq 4.18.

4.3.   Enforcement of the external boundary

More than one type of external boundary is possible in order to generate different periodic repeat structures or to compute effective conductivity in orthogonal directions. In every case the external boundary equation has both a homogeneous and inhomogeneous component. The homogeneous component for the standard case is $\partial\phi^{\mathrm{per}}/\partial u = 0$ on all edges, which is equivalent to zero normal current (insulation) on left and right edges of the boundary rectangle and zero tangential field (constant potential) on the top and bottom edges. This corresponds to reflection symmetry on the edges as shown in Figures 2 and 10. The homogeneous component is enforced at discrete boundary locations by a weighted least squares minimization. The inhomogeneous component expresses an average electric field in a particular direction—the positive $v$ direction in the standard case---in order to drive current flow. It is enforced as a Lagrangian constraint within the minimization procedure and at the same discrete locations. The general procedure is as follows.

Let $\mathcal{F}(\theta) = \mathcal{L}\phi$ be a linear integro-differential function of the potential $\phi$ that is chosen for the boundary condition. Operator $\mathcal{L}$ can vary for different edges around the perimeter (i.e., mixed boundaries). However, for the standard case $\mathcal{F} = \partial\phi^{\mathrm{per}}/\partial u$ for all edges as discussed above. The homogeneous component of the boundary equation is to force $F_{k} = \mathcal{F}(\theta_{k})\approx0$ for multiple $\theta_{k}$ values distributed around the boundary. In accord with Eqs 4.16 and 4.19:

where $M_{kn} = \mathcal{L}\phi_{n}^{\mathrm{per}}(\theta_{k})$ generally and in the standard case $M_{kn} = \partial\phi_{n}^{\mathrm{per}}(\theta_{k})/\partial u$ (obtained from Eq 4.20). A finite number of spectral coefficients is used in the sum, generally from 20 to 50, depending on accuracy needs. In matrix form Eq 4.22 is

Obviously $\mathbf{C} = \mathbf{0}$ would satisfy this equation. However, we seek a nontrivial solution that also satisfies the inhomogeneous driving-force component of the boundary equation, here given in matrix form:

where vector element $f_{n}$ (given in Eq 4.30 below) is the average difference in eigenpotential $\phi_{n}^{\mathrm{per}}$ between opposing edges.

Our least squares procedure minimizes a squared residual, in this case a weighted norm of vector $\mathbf{F = \mathbf{M\, C}}$, subject to constraint Eq 4.24. The error function is

where $\Lambda$ is a Lagrangian multiplier and $\mathbf{W}$ is a diagonal weighting matrix. The following operation will minimize $\mathcal{E}$ with respect to spectral coefficient vector $\mathbf{C}$.

The solution takes the form $\mathbf{C} = \Lambda\mathbf{C}^{(0)}$. Combining this with Eq 4.24 gives

where a numerical matrix equation solver is used to determine $\mathbf{C}^{(0)}$ from Eq 4.26:

Elements of the symmetric square matrix in Eq 4.28 are

Let us further assume current is driven in the positive $v$ direction. The corresponding potential difference vector $\mathbf{f}$ is given by

where $(e_{k}^{\mathrm{bot}}-e_{k}^{\mathrm{top}})$ selects the difference between the bottom and top edges of the rotated primitive cell, incorporates integration weight $W_{kk}$, and normalizes by the edge length:

With a solution for $C_{n}$ in hand, we can calculate the current through our primitive cell, which can immediately provide the effective conductivity of the composite periodic medium. Taking the standard case, noting that potential difference was constrained to be unity, and combining previously defined quantities gives:

This averages the current for the top and bottom edges in order to get effective conductivity. One can calculate $\sigma_{u}$ by a trivial adaptation of the above method to change the direction of the current flow. Alternatively for a symmetric square geometry one can just add $\tfrac{\pi}{2}$ to $\beta$ to rotate the external boundaries relative to the phases.

There are multiple possibilities for choice of boundary nodes $\theta_{k}$ and weights $W_{kk}$. Essentially we are choosing a quadrature routine for the implied integrals in Eqs 4.29, 4.30 and 4.32. A routine that proved effective was to use the midpoint rule with either uniform or nonuniform spacing in angle $\theta$. If a node spans an angular segment defined by $[\theta_{k-}, \theta_{k+}]$ then we let $2\ell_{\mathrm{per}}(\theta_{k}) = \ell_{\mathrm{per}}(\theta_{k+})+\ell_{\mathrm{per}}(\theta_{k-})$ and $W_{kk} = \ell_{\mathrm{per}}(\theta_{k+})-\ell_{\mathrm{per}}(\theta_{k-})$, where function $\ell_{\mathrm{per}}(\theta)$ is defined by Eq 4.15. Generally 240 nodes were used. Uniform spacing of angles worked well for many cases, but for more challenging rotated or stretched cases, nonuniform distribution of the boundary nodes proved more effective, with smaller node spacing where currents or potentials abruptly change. Other quadrature routines including adaptive refinement of nodes are possible but were not attempted.

The accuracy of the homogenized conductivity is essentially determined by how well the external boundary equations are enforced including currents and fields being zero or nonzero where they are supposed to. Therefore a check on the spectral solution was to compare the result of Eq 4.32 to a similarly calculated current orthogonal to a line crossing the middle of the primitive cell, namely for the segment $\xi-L\le u\le\xi$, $v = \psi-H/2$.

5.   Comparison and application of methods

Here we compare and apply the methods proposed above to show capabilities and in some cases limitations for grids of interest, including stretched rectangular checkerboards, isolated square inclusions, and rotated checkerboards. Additional illustrations of the results are given in the Appendix.

5.1.   Stretched two-phase alternating checkerboard

Figure 8 shows comparative results for normalized effective conductivity of a stretched two-phase alternating checkerboard where $\sigma_{1} = \sigma_{3} = 100\sigma_{2} = 100\sigma_{4}$, meaning contrast ratio $\chi = 100$. Aspect ratio $R$ is varied over 8 orders of magnitude. The exact results from Craster and Obnosov are generated from a combination of Eq 2.3 and a hypergeometric-function-based aspect-ratio correction ^[19]. The curves for the equivalent circuit models come from closed-form Eqs 3.6 and 3.7. The ECs and ECn methods each reproduce the exact result at $R = 0, 1, \, \infty$ and obey the duality principle (Eq 5.4) as shown by the symmetry in the plots. The ECn model differs from the exact answer by no more than 1.4% for this contrast ratio and no more than 1.7% for any contrast ratio ($0 < \chi < \infty$).

As shown in plot (b), for aspect ratios close to 1 the spectral model most accurately reproduces the exact result. The spectral model as implemented here suffers from numerical problems for aspect ratios less than 0.1 or greater than 5, though the errors appear bounded. This is due to the difficulty in getting a sufficient number of eigenfunctions to accurately describe the elongated external boundary while also having a well-conditioned matrix for determining the spectral coefficients. The spectral result for $R = 3$ is illustrated in Figure 1 (upper).

5.2.   Isolated square inclusions

Figure 9 gives results for a square grid of isolated square inclusions $(\sigma_{1}$) in a continuous background medium ($\sigma_{2} = \sigma_{3} = \sigma_{4}$) with contrast ratio $\chi = 100$. $p_{1}$ is the area fraction of quadrant 1 in the primitive cell and hence of phase 1 in the periodic system. The effective conductivity for the ECn method is given by

Also included is the result for Maxwell Garnett (MG) effective medium theory ^[26], adapted for 2D inclusions:

The spectral, ECs, and ECn methods each reproduce the exact result (Eq 2.3) at $p_{1} = 0, \, 0.25, \, 1$. For other cases the exact result is not known, but a smoothed fit of the spectral method was used as an estimate. Figure 9(b) shows how the conductivity curves depart from that estimate. The ECn solution differs by no more than 1.1% and is therefore considered the more reliable of the equivalent circuit methods. The simple 2D MG result is surprisingly accurate given that it was derived for dilute circular inclusions, differing from the best estimate by no more than 4.2%.

As shown in Figure 9(b), the spectral method has numerical accuracy problems when $p_{1}$ approaches (but is not exactly at) 1 and 0, though the errors are reasonably well bounded. As with the corresponding case in Section 5.1, this is due to the difficulty in describing a rapid oscillation in potential along the external boundaries with no more than 50 spectral coefficients, a limit beyond which matrix ill-conditioning became a problem in our implementation.

5.3.   Rotated cells

Periodic grids generated from diagonal or rotated primitive cells are another interesting comparison case. Craster and Obnosov ^[23] extended their earlier analytic work to the geometry of a $\tfrac{\pi}{4}$-rotated four-phase square checkerboard as shown in Figure 10. The diagonal lines of symmetry effectively generate a $4\times4$ doubly periodic system. Craster and Obnosov presented analytic solutions for two nontrivial cases, namely $\sigma_{1}\sigma_{2} = \sigma_{3}\sigma_{4}$ and $\sigma_{1}\sigma_{3} = \sigma_{2}\sigma_{4}$. Despite considerable progress, they were not able to solve the general four-phase rotated case. One particular case, for an isolated inclusion (i.e., $\sigma_{1}\ne\sigma_{2} = \sigma_{3} = \sigma_{4}$, $p_{1} = 0.25$), was not solved generally but was tackled in the high-contrast limit with the widely used principle of duality ^[20,27]. That problem is re-examined below using the spectral and equivalent circuit methods. The discussion is translated into our notation and geometry, consistent with Figure 10. This rotated system is also illustrated in Figure 1(lower).

Craster and Obnosov ^[23] suggest correctly that in the limit of an insulating diagonal inclusion, $\sigma_{1}\ll\sigma_{2}$, that effective conductivity is $\sigma_{u} = 2\sqrt{\sigma_{1}\sigma_{2}}$ and that in the limit of insulating background medium, $\sigma_{1}\gg\sigma_{2}$, that effective conductivity is $\sigma_{v} = \frac{1}{2}\sqrt{\sigma_{1}\sigma_{2}}$. However they incorrectly apply the duality principle to suggest that $\sigma_{v} = \frac{1}{2}\sqrt{\sigma_{1}\sigma_{2}}$ for $\sigma_{1}\ll\sigma_{2}$ and that $\sigma_{u} = 2\sqrt{\sigma_{1}\sigma_{2}}$ for $\sigma_{1}\gg\sigma_{2}$. The error is revealed when considering, for instance, that in the $v$ direction there is a continuous diagonal pathway of phase 2 (colored squares associated with $\sigma_{2}$, $\sigma_{3}$, and $\sigma_{4}$ in Figure 10), so that if $\sigma_{1}\ll\sigma_{2}$ the effective conductivity $\sigma_{v}$ is linear in $\sigma_{2}$ and not in $\sqrt{\sigma_{1}\sigma_{2}}$.

Figure 11 shows numerical results from the spectral method and a corresponding semi-empirical (SE) fit for the diagonal inclusion system over 12 orders of magnitude in contrast ratio. Also shown are the ECr model results (Eq 3.8). The semi-empirical fit has analytically correct asymptotic behavior and obeys the duality principle, but nevertheless uses an empirically derived constant:

A reasonably good fit of the spectral model data is given by $w = 0.66$. If $w = \tfrac{1}{2}$ were used in the SE fit then it would correspond exactly to the ECr model. The ECr model produces a reasonably good approximation over a wide range of contrast ratios, but (based on the value of $w$) errors can be as high as 32% in high-contrast cases.

The duality principle ^[20,27] in this context means that

The problem with using the duality principle to derive a missing effective conductivity (e.g., $\sigma_{u}$ from $\sigma_{v}$) in the high-contrast limit is that there can be a vital term that vanishes in the limit, such as $w\sqrt{\chi}^{-1}$ in Eq 5.3 with $\chi\gg1$, and is therefore not recovered under an application of duality.

Before concluding this section, we note that the simple ECr model (Eq 3.8) reproduces the Craster and Obnosov rotated results for their two principal analytic cases: (1) $\sigma_{1}\sigma_{4} = \sigma_{2}\sigma_{3}$ and (2) $\sigma_{1}\sigma_{3} = \sigma_{2}\sigma_{4}$ ^[23]. For case (1) the result is exact. For case (2) the solution generated by the ECr model should be multiplied by a hypergeometric-based scaling function that, as Craster and Obnosov have shown, changes the effective value of conductivity by 5% or less. This and the analysis above suggests that the ECr model can with reasonable accuracy rapidly homogenize the rotated four-phase checkerboard for arbitrary conductivities.

6.   Conclusions

There are many difficulties involved in determining effective properties across varied microstructure, and a rich history of approaches to the problem. Effective conductivity, for instance, depends on the path that current takes through the heterogeneous material, which in turn depends on the boundary conditions employed. This is made particularly obvious by how effective conductivity can change with direction of the applied external field, even when underlying phases are isotropic. Only when the length scale of the heterogeneities is small compared to the length scale of the outer boundaries can one approach the idea of a single effective modulus for a composite.

In many cases of practical interest, one does not need an exact result for effective modulus, but rather one that is reasonably accurate, provided it can be computed rapidly. Such can be used for renormalization, multigrid solvers, and reduced-order models ^[5,9]. Although renormalization itself is not the focus of this paper, this work was motivated by a recognition that renormalization schemes are more useful if they can be applied to primitive cells (e.g., a $2\times2$ block of phases) in which the underlying phases are not necessarily square, symmetric, and isotropic.

The equivalent circuit models provide closed-form solutions easily applied to periodic grids or to renormalization schemes. Between the methods, ECn appears more accurate in general than ECs for nonsquare geometries, while ECs explicitly allows anisotropic underlying phases. The ECr method is also a useful tool for understanding effective conductance in primitive cells rotated with respect to the phase boundaries.

The spectral method expands the number of grid systems that can be accurately simulated, including those for which there is not yet an analytic solution and for which finite element methods are not well-suited. The method exactly handles traditionally difficult internal boundaries, namely large contrast ratios between neighboring phases. The external boundaries on the primitive cell have a larger impact on the stability of the spectral method. Numerical difficulty can arise in cases where current or potential fluctuate rapidly when moving along the perimeter of the primitive cell. The spectral method in its current form cannot match the accuracy of Helsing's numerical methods for square lattices ^[13], but programming simplicity and flexibility recommend it. We are currently exploring the extension of the spectral and equivalent circuit methods to 3D systems.

Acknowledgments

This work was partially supported by the Assistant Secretary for Energy Efficiency and Renewable Energy, Office of Vehicle Technologies of the U. S. Department of Energy.

Conflict of interest

The authors declare no conflicts of interest.

Appendix

Figures 12-14 shown below are wireframe potential maps generated from the spectral method for representative composites in Section 5. The potential from a single primitive cell has been replicated to illustrate the effect of the periodic phase boundaries.