
This work presents models for homogenizing or finding the effective transport or mechanical properties of microscale composites formed from highly contrasting phases described on a grid. The methods developed here are intended for engineering applications where speed and geometrical flexibility are a premium. A canonical case that is mathematically challenging and yet can be applied to many realistic materials is a 4-phase 2-dimensional periodic checkerboard or tiling. While analytic solutions for calculating effective properties exist for some cases, versatile methods are needed to handle anisotropic and non-square grids. A reinterpretation and extension of an existing analytic solution that utilizes equivalent circuits is developed. The resulting closed-form expressions for effective conductivity are shown to be accurate within a few percent or better for multiple cases of interest. Secondly a versatile and efficient spectral method is presented as a solution to the 4-phase primitive cell with a variety of external boundaries. The spectral method expresses the solution to effective conductivity in terms of analytically derived eigenfunctions and numerically determined spectral coefficients. The method is validated by comparing to known solutions and can allow extensions to cases with no current analytic solution.
Citation: Ben J. Ransom, Dean R. Wheeler. Rapid computation of effective conductivity of 2D composites by equivalent circuit and spectral methods[J]. Mathematics in Engineering, 2022, 4(3): 1-24. doi: 10.3934/mine.2022020
[1] | Zhiwen Zhao . Exact solutions for the insulated and perfect conductivity problems with concentric balls. Mathematics in Engineering, 2023, 5(3): 1-11. doi: 10.3934/mine.2023060 |
[2] | Giacomo Canevari, Arghir Zarnescu . Polydispersity and surface energy strength in nematic colloids. Mathematics in Engineering, 2020, 2(2): 290-312. doi: 10.3934/mine.2020015 |
[3] | E. Artioli, G. Elefante, A. Sommariva, M. Vianello . Homogenization of composite materials reinforced with unidirectional fibres with complex curvilinear cross section: a virtual element approach. Mathematics in Engineering, 2024, 6(4): 510-535. doi: 10.3934/mine.2024021 |
[4] | Raimondo Penta, Ariel Ramírez-Torres, José Merodio, Reinaldo Rodríguez-Ramos . Effective governing equations for heterogenous porous media subject to inhomogeneous body forces. Mathematics in Engineering, 2021, 3(4): 1-17. doi: 10.3934/mine.2021033 |
[5] | Valentina Volpini, Lorenzo Bardella . Asymptotic analysis of compression sensing in ionic polymer metal composites: The role of interphase regions with variable properties. Mathematics in Engineering, 2021, 3(2): 1-31. doi: 10.3934/mine.2021014 |
[6] | Caterina Ida Zeppieri . Homogenisation of high-contrast brittle materials. Mathematics in Engineering, 2020, 2(1): 174-202. doi: 10.3934/mine.2020009 |
[7] | Patrick Dondl, Martin Jesenko, Martin Kružík, Jan Valdman . Linearization and computation for large-strain visco-elasticity. Mathematics in Engineering, 2023, 5(2): 1-15. doi: 10.3934/mine.2023030 |
[8] | Elena Beretta, M. Cristina Cerutti, Luca Ratti . Lipschitz stable determination of small conductivity inclusions in a semilinear equation from boundary data. Mathematics in Engineering, 2021, 3(1): 1-10. doi: 10.3934/mine.2021003 |
[9] | Ko-Shin Chen, Cyrill Muratov, Xiaodong Yan . Layered solutions for a nonlocal Ginzburg-Landau model with periodic modulation. Mathematics in Engineering, 2023, 5(5): 1-52. doi: 10.3934/mine.2023090 |
[10] | Laurent Baratchart, Sylvain Chevillard, Adam Cooman, Martine Olivi, Fabien Seyfert . Linearized active circuits: transfer functions and stability. Mathematics in Engineering, 2022, 4(5): 1-18. doi: 10.3934/mine.2022039 |
This work presents models for homogenizing or finding the effective transport or mechanical properties of microscale composites formed from highly contrasting phases described on a grid. The methods developed here are intended for engineering applications where speed and geometrical flexibility are a premium. A canonical case that is mathematically challenging and yet can be applied to many realistic materials is a 4-phase 2-dimensional periodic checkerboard or tiling. While analytic solutions for calculating effective properties exist for some cases, versatile methods are needed to handle anisotropic and non-square grids. A reinterpretation and extension of an existing analytic solution that utilizes equivalent circuits is developed. The resulting closed-form expressions for effective conductivity are shown to be accurate within a few percent or better for multiple cases of interest. Secondly a versatile and efficient spectral method is presented as a solution to the 4-phase primitive cell with a variety of external boundaries. The spectral method expresses the solution to effective conductivity in terms of analytically derived eigenfunctions and numerically determined spectral coefficients. The method is validated by comparing to known solutions and can allow extensions to cases with no current analytic solution.
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×2 blocks of squares, such that a grid of 2n×2n 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×2 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 σ 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×2 rectangular grid, leading to closed-form solutions (Section 3). In addition we provide a spectral method for solving the 2×2 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.
The primitive unit cell considered here is the case of a 2×2 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×2 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×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 σ1=σ3≠σ2=σ4. Such a two-phase alternating checkerboard can be homogenized by the exact solution σ◻x=σ◻y=√σ1σ2, where σ◻i 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 (σLx) and upper (σUx) bounds that are themselves composed of alternating arithmetic and harmonic means of the quadrant conductivities:
σLx=σ1σ2σ3σ4(σ−11+σ−12+σ−13+σ−14)(σ1+σ2)(σ3+σ4)=(1σ1+1σ2)−1+(1σ3+1σ4)−1 | (2.1) |
σUx=(σ2+σ3)(σ4+σ1)σ1+σ2+σ3+σ4=(1σ2+σ3+1σ4+σ1)−1 | (2.2) |
σ◻x=√σLxσUx | (2.3) |
The results for σ◻y can likewise be obtained by transposing σ2 and σ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].
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 ∞. 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) K1y=σ1L/H and (b) K1x=σ1H/L where L is the quadrant length in x and H is the height in y and an isotropic σ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, ϕ12 and ϕ34, assuming an overall potential difference of 1 across the circuit. Conservation of current at the interior nodes leads to the matrix equation
[K1x+K2x+K∗y−K∗y−K∗yK∗y+K3x+K4x][ϕ12ϕ34]=[K1xK4x] | (3.1) |
The total current through the circuit is I=Kx=K2xϕ12+K3xϕ34. Combining these relations leads to an expression for the effective conductance:
Kx=[K2xK3x]T[K1x+K2x+K∗y−K∗y−K∗yK∗y+K3x+K4x]−1[K1xK4x] | (3.2) |
The same process is repeated for Figure 3(a) to give
Ky=[K4yK3y]T[K1y+K4y+K∗x−K∗x−K∗xK∗x+K2y+K3y]−1[K1yK2y] | (3.3) |
We now assert that K∗x=Kx and K∗y=Ky, 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
Kx=gx+√g2x+hxKy=gy+√g2y+hy | (3.4) |
where
gx=12(−axby+aybx+cxdy−cydx)/(axdy+aydx)gy=12(−axby+aybx−cxdy+cydx)/(bxdy+bydx)hx=(bxcy+bycx)/(axdy+aydx)hy=(axcy+aycx)/(bxdy+bydx)ai=(K1i+K2i)(K3i+K4i)bi=(K1i+K4i)(K2i+K3i)ci=K1iK2iK3iK4i(K−11i+K−12i+K−13i+K−14i)di=K1i+K2i+K3i+K4i | (3.5) |
and i represents the x or y direction. Considerable simplification is possible for cases with isotropic phases, geometric symmetry, or both. With Kx and Ky determined, they can be converted to effective σx and σ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 (σ1=σ3≠σ2=σ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 χ=σ1/σ2 is the contrast ratio.
σECsy√σ1σ2=√σ1σ2σECsx=(1−R−2√χ+√χ−1)+√(1−R−2√χ+√χ−1)2+R−2 | (3.6) |
When R→∞ this formula produces σy→2/(σ−11+σ−12), and when R→0 then σy→12(σ1+σ2), demonstrating the exact Maxwell lamellar limits. Similarly, when R=1 then σx=σy=√σ1σ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×2 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 Lmin. 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
σECny√σ1σ2=√σ1σ2σECnx={2(1−R−1)√χ+√χ−1+R−1ifR≥1(2(1−R)√χ+√χ−1+R)−1ifR<1 | (3.7) |
which can be compared to Eq 3.6 and similarly produces the correct limits for R=0,1,∞. 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 π4-rotated system (see Section 5.3) with the following simultaneous conductivity substitutions: σ1←√σ1σ2, σ2←√σ2σ3, σ3←√σ3σ4, σ4←√σ4σ1, σ◻x←σECru, and σ◻y←σECrv 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:
σECru=(1√σ2σ3+√σ3σ4+1√σ4σ1+√σ1σ2)−1σECrv=(1√σ1σ2+√σ2σ3+1√σ3σ4+√σ4σ1)−1 | (3.8) |
This prediction is compared to other methods in Section 5.3.
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 π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.
The spectral method starts with a proposed solution to the Laplace equation for potential ϕ in polar coordinates (Figure 6) by the standard technique of separation of variables:
ϕα(r,θ)=rλ[Aαcos(λθ)+Bαsin(λθ)] | (4.1) |
The subscript α indicates a solution for quadrant 1, 2, 3, or 4. Note that all physical quantities are taken to be dimensionless. Coefficients Aα and Bα depend on quadrant while λ is a common eigenvalue. The current density in the θ direction is given by
iθ,α=−σαr∂ϕα∂θ=λrλ−1[σαAαsin(λθ)−σαBαcos(λθ)] | (4.2) |
We enforce the internal boundaries between contiguous quadrants in order to determine Aα, Bα, and λ in terms of conductivities σα. At each boundary both ϕ and iθ are continuous, independent of r value. (Note that instead of potential we could have equivalently specified that the tangential electric field, in this case −∂ϕα/∂r, is continuous across the boundary.) The resulting 8 equations can be expressed in matrix form:
G[AB]=0 | (4.3) |
where A and B are respective column vectors of Aα and Bα values, stacked to become a single vector (i.e., a block matrix),
G=[c1−c100s1−s100σ1s1−σ2s100−σ1c1σ2c1000c2−c200s2−s200σ2s2−σ3s200−σ2c2σ3c2000c3−c300s3−s300σ3s3−σ4s300−σ3c3σ4c3−100c4000s4000σ4s4σ100−σ4c4] | (4.4) |
and for compactness cα=cos(απ2λ) and sα=sin(απ2λ).
Equation 4.3 requires that
detG=0 | (4.5) |
and this relationship is used to solve for λ. But first we reduce the number of unknowns by the substitutions cα=Tα(c1) and sα=s1Uα−1(c1), where Tα and Uα are respectively Chebyshev polynomials of the first and second kind and order α. These relationships come from the multiple angle formulas for sine and cosine. The identity c21=1−s21 is also used.
After these substitutions and considerable algebra we obtain from Eq 4.5 that s21=ς2 (the primary solution) and s1=0 (the secondary solution), where ς is a function of the quadrant conductivities:
ς=√1−(σ1σ3−σ2σ4)2(σ1+σ2)(σ2+σ3)(σ3+σ4)(σ4+σ1) | (4.6) |
From ς we can determine an infinite series of eigenvalues for the primary solution using
λ1=2πasin(ς)λn=(n−12)+(−1)n+1(λ1−12);n=1,2,…,∞ | (4.7) |
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 (s1=0) produces a series of paired eigenvalues λn, namely repeated even numbers (i.e., 2, 2, 4, 4, 6, 6, …). The eigenvalue λ=0 is omitted from the solution due to our intention to keep ϕ=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
ϕα(r,θ)=∑nCnrλn[Aαcos(λnθ)+(−1)n+1Bαsin(λnθ)] | (4.8) |
where for generality and compactness we have combined the primary and secondary solutions into a single set of spectral coefficients Cn and eigenvalues λ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 c1 relative to s1 for the sequence of eigenvalues, which in turn changes the sign of B relative to A.
In order to calculate coefficient vectors A and B, we consider below three cases: (1) s21=0, (2) s21=1, and (3) 0<s21<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 σ1=σ3=0 or σ2=σ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 A and B respectively and allow us one degree of freedom for each vector. For instance, one could choose
A1=A2=A3=A4=1σ1B1=σ2B2=σ3B3=σ4B4=1 | (4.9) |
The Aα and Bα terms in Eq 4.8 are further scaled and maintain their respective independence by virtue of the Cn coefficients, the pairs of duplicate eigenvalues, and the alternating sign of the (−1)n term.
Case (2) occurs when σ1σ3=σ2σ4. It is similar to case (1) because, with s21=1, Eq 4.7 produces a series of duplicated eigenvalues, namely the odd values 1,1,3,3,5,5,…. The internal boundary conditions again require certain ratios among elements of A and B respectively and allow independent scaling of A and B. For instance, one could choose
A1σ2+σ3=A2σ1+σ4=A3σ1+σ4=A4σ2+σ3=1B1σ3+σ4=B2σ3+σ4=B3σ1+σ2=B4σ1+σ2=1 | (4.10) |
Case (3) is the most general where 0<s21<1. Matrix G (Eq 4.4) then has rank 7, meaning there is one degree of freedom and that A and B are coupled. We therefore stipulate that A1=1. This can be incorporated into our system of equations by use of vector g=[10000000] T that can augment G to create an invertible matrix. A and B can therefore be obtained from
[AB]=(G+ggT)−1g | (4.11) |
where G is evaluated at s1=ς and c1=√1−ς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 c1 and s1 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 ς and λ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(πλ)=1−2Δ2. It can be shown that their eigenvalue λ is related to our λ1 by λ1=1−λ and that ς2=1−Δ2. The analytic solution to the effective conductivities σ◻x and σ◻y of the four-phase square checkerboard (i.e., Eq 2.3) can be expressed in terms of ς:
ς=(1σ2+σ3+1σ4+σ1)σ◻x=(1σ1+σ2+1σ3+σ4)σ◻y | (4.12) |
Furthermore, in the general case where the four-phase system is stretched and aspect ratio R ranges from 0 to ∞, 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 ς and its inverse:
ς≤(σxσ◻x,σyσ◻y)≤ς−1 | (4.13) |
or equivalently σLi/σUi=ς2 for either direction, where σLi and σUi are the Maxwell lamellar limits (Eqs 2.1 and 2.2). Eq 4.13 demonstrates for the case ς=1 (i.e., σ1σ3=σ2σ4) that effective conductivity of the system is independent of aspect ratio and is isotropic (i.e., σx=σy).
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 Cn 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 Cn 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 ξ, L−ξ, ψ, and H−ψ. Angles at the rectangle vertices are given by θ1=acot(ξ/ψ), θ2=acot[(ξ−L)/ψ], θ3=π+acot[(L−ξ)/(H−ψ)], and θ4=π+acot[ξ/(ψ−H)]. The bounding rectangle with its associated uv coordinates are rotated by angle β relative to the xy coordinates. Coordinate θ is zero at the positive u semi-axis, meaning θ 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 β, phase boundaries may or may not align with the uv axes or the four θi values, though in this work only multiples of π/4 were used for β.
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√2.
Our intention is to define a boundary equation completely in terms of θ. Therefore we use the following piecewise function for the radius at the perimeter
rper(θ)={ψcsc(θ)ifθ1≤θ<θ2(ξ−L)sec(θ)ifθ2≤θ<θ3(ψ−H)csc(θ)ifθ3≤θ<θ4ξsec(θ)if0≤θ<θ1orθ4≤θ<2π | (4.14) |
A second useful function tracks the accumulated linear distance along the perimeter from θ=0 to 2π:
ℓper(θ)={ξtan(θ)if0≤θ<θ1ξ+ψ−ψcot(θ)ifθ1≤θ<θ2L+2ψ+(L−ξ)tan(θ)ifθ2≤θ<θ32L+H−ξ+ψ−(H−ψ)cot(θ)ifθ3≤θ<θ42L+2H+ξtan(θ)ifθ4≤θ<2π | (4.15) |
The potential along the rotated perimeter is given by an application of Eq 4.8 to get
ϕper(θ)=∑nCnϕpern(θ) | (4.16) |
Perimeter eigenpotentials are given by
ϕpern(θ)=[rper(θ)]λn[A(θ∙)cos(λnθ∙)+(−1)n+1B(θ∙)sin(λnθ∙)] | (4.17) |
where θ∙=(θ+β)mod2π 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.,
A(θ∙)={A1if0≤θ∙<π2A2ifπ2≤θ∙<πA3ifπ≤θ∙<3π2A4if3π2≤θ∙<2π | (4.18) |
To properly formulate external boundary conditions in terms of fields and fluxes, we need to define the gradient of ϕ with respect to the uv coordinates and evaluate it along the perimeter:
∇ϕper(θ)=∑nCn∇ϕpern(θ) | (4.19) |
where
∇ϕpern(θ)=λn[rper(θ)]λn−1[cos(λnθ∙−θ)sin(λnθ∙−θ)−sin(λnθ∙−θ)cos(λnθ∙−θ)][A(θ∙)(−1)n+1B(θ∙)] | (4.20) |
With gradients defined we can also calculate eigencurrents on the perimeter in the u and v directions:
[ipern,uipern,v]=−σ(θ∙)∇ϕpern(θ) | (4.21) |
where σ(θ∙) is defined in the same manner as Eq 4.18.
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 ∂ϕper/∂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 F(θ)=Lϕ be a linear integro-differential function of the potential ϕ that is chosen for the boundary condition. Operator L can vary for different edges around the perimeter (i.e., mixed boundaries). However, for the standard case F=∂ϕper/∂u for all edges as discussed above. The homogeneous component of the boundary equation is to force Fk=F(θk)≈0 for multiple θk values distributed around the boundary. In accord with Eqs 4.16 and 4.19:
Fk=∑nCnMkn≈0 | (4.22) |
where Mkn=Lϕpern(θk) generally and in the standard case Mkn=∂ϕpern(θk)/∂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
F=MC≈0 | (4.23) |
Obviously C=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:
fTC=1 | (4.24) |
where vector element fn (given in Eq 4.30 below) is the average difference in eigenpotential ϕpern between opposing edges.
Our least squares procedure minimizes a squared residual, in this case a weighted norm of vector F=MC, subject to constraint Eq 4.24. The error function is
E=12FTWF+Λ(1−fTC) | (4.25) |
where Λ is a Lagrangian multiplier and W is a diagonal weighting matrix. The following operation will minimize E with respect to spectral coefficient vector C.
∂E∂C=0=(MTWM)C−Λf | (4.26) |
The solution takes the form C=ΛC(0). Combining this with Eq 4.24 gives
C=C(0)fTC(0) | (4.27) |
where a numerical matrix equation solver is used to determine C(0) from Eq 4.26:
C(0)=(MTWM)−1f | (4.28) |
Elements of the symmetric square matrix in Eq 4.28 are
(MTWM)mn=∑kLϕperm(θk)WkkLϕpern(θk):generally=∑k∂ϕperm(θk)∂uWkk∂ϕpern(θk)∂u:standardcase | (4.29) |
Let us further assume current is driven in the positive v direction. The corresponding potential difference vector f is given by
fn=∑k(ebotk−etopk)ϕpern(θk) | (4.30) |
where (ebotk−etopk) selects the difference between the bottom and top edges of the rotated primitive cell, incorporates integration weight Wkk, and normalizes by the edge length:
etopk=WkkL{1ifθ1≤θk<θ20otherwiseebotk=WkkL{1ifθ3≤θk<θ40otherwise | (4.31) |
With a solution for Cn 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:
σv=∑kH2(etopk+ebotk)∑nCnipern,v(θk) | (4.32) |
This averages the current for the top and bottom edges in order to get effective conductivity. One can calculate σ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 π2 to β to rotate the external boundaries relative to the phases.
There are multiple possibilities for choice of boundary nodes θk and weights Wkk. 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 θ. If a node spans an angular segment defined by [θk−,θk+] then we let 2ℓper(θk)=ℓper(θk+)+ℓper(θk−) and Wkk=ℓper(θk+)−ℓper(θk−), where function ℓper(θ) 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 ξ−L≤u≤ξ, v=ψ−H/2.
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.
Figure 8 shows comparative results for normalized effective conductivity of a stretched two-phase alternating checkerboard where σ1=σ3=100σ2=100σ4, meaning contrast ratio χ=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,∞ 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<χ<∞).
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).
Figure 9 gives results for a square grid of isolated square inclusions (σ1) in a continuous background medium (σ2=σ3=σ4) with contrast ratio χ=100. p1 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
σECnxσ2={(√3+χ3χ+1+1−2√p12√p1)−1+1−2√p1if0≤√p1≤12(√3+χ3χ+1+2√p1−1(1−√p1)(χ+1))−1+2√p1−1√p1(χ−1−1)+1if12<√p1≤1 | (5.1) |
Also included is the result for Maxwell Garnett (MG) effective medium theory [26], adapted for 2D inclusions:
σMGxσ2=(χ+1)+p1(χ−1)(χ+1)−p1(χ−1) | (5.2) |
The spectral, ECs, and ECn methods each reproduce the exact result (Eq 2.3) at p1=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 p1 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.
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 π4-rotated four-phase square checkerboard as shown in Figure 10. The diagonal lines of symmetry effectively generate a 4×4 doubly periodic system. Craster and Obnosov presented analytic solutions for two nontrivial cases, namely σ1σ2=σ3σ4 and σ1σ3=σ2σ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., σ1≠σ2=σ3=σ4, p1=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, σ1≪σ2, that effective conductivity is σu=2√σ1σ2 and that in the limit of insulating background medium, σ1≫σ2, that effective conductivity is σv=12√σ1σ2. However they incorrectly apply the duality principle to suggest that σv=12√σ1σ2 for σ1≪σ2 and that σu=2√σ1σ2 for σ1≫σ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 σ2, σ3, and σ4 in Figure 10), so that if σ1≪σ2 the effective conductivity σv is linear in σ2 and not in √σ1σ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:
σSEu√σ1σ2=[12+(12−w)χ1−w+w√χ]−1σSEv√σ1σ2=12+(12−w)χw−1+w√χ−1 | (5.3) |
A reasonably good fit of the spectral model data is given by w=0.66. If w=12 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
σu(σ1,σ2)σv(σ−11,σ−12)=1σu(σ1,σ2)σv(σ2,σ1)=σ1σ2 | (5.4) |
The problem with using the duality principle to derive a missing effective conductivity (e.g., σu from σv) in the high-contrast limit is that there can be a vital term that vanishes in the limit, such as w√χ−1 in Eq 5.3 with χ≫1, 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) σ1σ4=σ2σ3 and (2) σ1σ3=σ2σ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.
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×2 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.
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.
The authors declare no conflicts of interest.
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.
[1] | G. Milton, The theory of composites, Cambridge: Cambridge University Press, 2002. |
[2] |
S. Torquato, Optimal design of heterogeneous materials, Annu. Rev. Mater. Res., 40 (2010), 101-129. doi: 10.1146/annurev-matsci-070909-104517
![]() |
[3] |
D. Felbacq, G. Bouchitte, Homogenization of a set of parallel fibres, Waves in Random Media, 7 (1997), 245-256. doi: 10.1088/0959-7174/7/2/006
![]() |
[4] |
G. Milton, Proof of a conjecture on the conductivity of checkerboards, J. Math Phys., 42 (2001), 4873-4882. doi: 10.1063/1.1385564
![]() |
[5] |
M. Karim, K. Krabbenhoft, New renormalization schemes for conductivity upscaling in heterogeneous media, Transport Porous Med., 85 (2010), 677-690. doi: 10.1007/s11242-010-9585-9
![]() |
[6] |
W. Suen, S. Wong, K. Young, The lattice model of heat conduction in a composite material, J. Phys. D Appl. Phys., 12 (1979), 1325-1338. doi: 10.1088/0022-3727/12/8/013
![]() |
[7] | S. A. Berggren, D. Lukkassen, A. Meidell, L. Simula, A new method for numerical solution of checkerboard fields, J. Appl. Math., 1 (2001), 615091. |
[8] |
J. Bigalke, Derivation of an equation to calculate the average conductivity of random networks, Physica A, 285 (2000), 295-305. doi: 10.1016/S0378-4371(00)00262-4
![]() |
[9] |
D. Yushu, K. Matouš, The image-based multiscale multigrid solver, preconditioner, and reduced order model, J. Comput. Phys., 406 (2020), 109165. doi: 10.1016/j.jcp.2019.109165
![]() |
[10] |
L. Zielke, T. Hutzenlaub, D. Wheeler, C. W. Chao, I. Manke, A. Hilger, et al., Three-phase multiscale modeling of a LiCoO2 cathode: Combining the advantages of FIB-SEM imaging and X-ray tomography, Adv. Energy Mater., 5 (2015), 1401612. doi: 10.1002/aenm.201401612
![]() |
[11] | D. E. Stephenson, B. C. Walker, C. B. Skelton, E. Gorzkowski, Modeling 3D microstructure and ion transport in porous Li-ion battery electrodes, J. Electrochem. Soc., 158 (2011): A781-A789. |
[12] |
L. Jylha, A. Sihvola, Approximations and full numerical simulations for the conductivity of three dimensional checkerboard geometries, IEEE T. Dielect. El. In., 13 (2006), 760-764. doi: 10.1109/TDEI.2006.1667733
![]() |
[13] |
J. Helsing, The effective conductivity of random checkerboards, J. Comput. Phys., 230 (2011), 1171-1181. doi: 10.1016/j.jcp.2010.10.033
![]() |
[14] |
R. Edwards, J. Goodman, A. Sokal, Multigrid method for the random-resistor problem, Phys. Rev. Lett., 61 (1988), 1333-1335. doi: 10.1103/PhysRevLett.61.1333
![]() |
[15] |
S. Torquato, I. Kim, D. Cule, Effective conductivity, dielectric constant, and diffusion coefficient of digitized composite media via first-passage-time equations, J. Appl. Phys., 85 (1999), 1560-1571. doi: 10.1063/1.369287
![]() |
[16] |
S. S. Wei, J. S. Shen, W. Y. Yang, Z. L. Li, S. H. Di, C. Ma, Application of the renormalization group approach for permeability estimation in digital rocks, J. Petrol. Sci. Eng., 179 (2019), 631-644. doi: 10.1016/j.petrol.2019.04.057
![]() |
[17] | P. King, Upscaling permeability: Error analysis for renormalization, Transport Porous Med., 23 (1996), 337-354. |
[18] | S. Mortola, S. Steffé, A two-dimensional homogenization problem, Atti Accad. Naz. Lincei Rend. Cl. Sci. Fis. Mat. Natur., 78 (1985), 77-82. |
[19] |
R. Craster, Y. Obnosov, Four-phase checkerboard composites, SIAM J. Appl. Math., 61 (2001), 1839-1856. doi: 10.1137/S0036139900371825
![]() |
[20] |
J. Keller, A theorem on the conductivity of a composite medium, J. Math. Phys., 5 (1964), 548-549. doi: 10.1063/1.1704146
![]() |
[21] | J. Maxwell, A treatise on electricity and magnetism, 2 Eds., Oxford: Cambridge, 1881,398-411. |
[22] |
J. Keller, Effective conductivity of periodic composites composed of two very unequal conductors, J. Math. Phys., 28 (1987), 2516-2520. doi: 10.1063/1.527741
![]() |
[23] | R. Craster, Y. Obnosov, A model four-phase square checkerboard structure, Q. J. Mech. Appl. Math., 59 (2006), 20-23. |
[24] |
R. Craster, R. Obnosov, A three-phase tessellation: Solution and effective properties, Proc. R. Soc. Lond. A, 460 (2004), 1017-1037. doi: 10.1098/rspa.2003.1196
![]() |
[25] |
J. Helsing, Corner singularities for elliptic problems: special basis functions versus 'brute force', Commun. Numer. Meth. Eng., 16 (2000), 37-46. doi: 10.1002/(SICI)1099-0887(200001)16:1<37::AID-CNM307>3.0.CO;2-1
![]() |
[26] | T. Choy, Effective medium theory, Oxford: Clarendon Press, 1999. |
[27] | A. Dykhne, Conductivity of a two-dimensional two-phase system, Soviet Journal of Experimental and Theoretical Physics, 32 (1971), 63. |