A NEW METHOD FOR NUMERICAL SOLUTION OF CHECKERBOARD FIELDS

for a standard checkerboard of conductivity λg and λw. The explicit solution of the corresponding temperature-field was later found by Berdichevskii [1]. In particular he found that the heat-flux is infinitely high in the corners of the squares. Subsequently, explicit solutions for rectangular and triangular checkerboards were found in [17, 19, 20]. Mortola and Steffé [18] presented in 1985 an interesting conjecture concerning the effective conductivity of four phase checkerboards. Many attempts were made to prove/disprove the conjecture, even by specialists in homogenization theory (see [22]), but the problem remained unsolved for the rest of the century. Very recently the conjecture has been proved by Craster and Obnosov [5] and independently by Milton [16] (using a completely different proof). There are many interesting works on other types of checkerboards. Concerning three-dimensional checkerboards, we refer to [12, 15] (see also [14]).


Introduction
Very few microstructures yield explicit formulae for their effective conductivity.One type of such structures is checkerboards.By using a duality argument Dychne [6] proved about 30 years ago the famous formula for the effective conductivity for a standard checkerboard of conductivity λ g and λ w .The explicit solution of the corresponding temperature-field was later found by Berdichevskii [1].
In particular he found that the heat-flux is infinitely high in the corners of the squares.Subsequently, explicit solutions for rectangular and triangular checkerboards were found in [17,19,20].Mortola and Steffe ´ [18] presented in 1985 an interesting conjecture concerning the effective conductivity of four phase checkerboards.Many attempts were made to prove/disprove the conjecture, even by specialists in homogenization theory (see [22]), but the problem remained unsolved for the rest of the century.Very recently the conjecture has been proved by Craster and Obnosov [5] and independently by Milton [16] (using a completely different proof).
There are many interesting works on other types of checkerboards.Concerning three-dimensional checkerboards, we refer to [12,15] (see also [14]).
Random checkerboards where studied in [2,13,7,21].A new type has recently been considered in [10] where the conductivity blows up at the corners in a diagonal direction and degenerates in the orthogonal direction.There appears to be two natural ways of defining the effective conductivity for this problem in terms of variational problems involving weighted Sobolev spaces.However, the corresponding effective conductivities are very different (Lavrentiev phenomenon).Similar observations have been done for nonlinear checkerboards in [23], which also contains a generalization of the Dykhne formula to power-law materials.
Due to the behavior of the solutions near the corner points it is difficult to solve the corresponding variational problems by usual numerical methods, even for the standard checkerboard.In this paper, we consider a generalized version of the standard checkerboard and focus on these difficulties, both theoretically and experimentally (numerical experiments).Moreover, we present a new numerical method for determining the corresponding field which converges in the energy norm independent of the local conductivities.

The model problem
We consider the stationary heat conduction problem for the checkerboard structure in Figure 2.1, where the conductivity λ(x) in each quarter V of a period is given by λ(x) = k g l(r) on grey parts, k w (l(r)) −1 on white parts, ( where 0 < β ≤ l(r) ≤ γ < ∞ for some constants β and γ, r is the distance from the center of V, and l is continuous at r = 0. Due to symmetries it is enough to only consider the set V for the calculation of the effective conductivity λ eff .This value can be determined by the following variational formula: or equivalently by (cf.[9]) where where e 1 = [1, 0].By (2.2) and (2.3) we can obtain the formula Stein A. Berggren et al. 159 which was first proved by Dychne [6].If W is any subspace of W and λ + eff and λ − eff are the values defined by then it is clear that λ − eff ≤ λ eff ≤ λ + eff .

Solutions found in finite-dimensional spaces
The minimizer u + of (2.6) is an approximation of the minimizer u of (2.2).When the ratio k g /k w is large, it becomes almost impossible to obtain good approximations by using the finite element method (FEM).In fact we have the following result.
Remark 3.2.As a consequence we obtain that u − u + 2 λ /λ eff → ∞ as k g /k w → ∞, where v λ denotes the energy norm 160 A new method for numerical solution of checkerboard fields This follows since Proof.Let c = min{l(r)} and let Ω g be the grey area of V. Since ) we obtain that Hence, the theorem follows if we can prove that inf u∈W Ωg Du + e 1 2 dx > 0. (3.6) Let Z denote the subspace of L 2 (Ω g ,R 2 ) of the gradients of all functions in W . Since Z is finite-dimensional, it is closed, and, therefore the closest approximation to −e 1 exists in Z, that is, we have the existence of the minimum Du + e 1 2 dx = 0 (3.8) (this will lead to a contradiction).Then To show this, let F be the interior of the triangle defined by the three corners (1, 0), (0, 1) and (0, 0).Rotate the coordinate axis x 1 and x 2 by 45 • counterclockwise and denote the new coordinate axis by x1 and x2 .u has a representative (still denoted u) which is absolutely continuous on almost all line segments and whose partial derivatives belong to L 2 (F) (cf.[24, Theorem 2.1.4.]).Let I t be the line segment between the points (0, t) and (t, 0).By the Jensen inequality We observe that Thus because (3.9) and (3.10) give us the inequality Noting that the left side is ∞, we obtain that F |Du| 2 dx = ∞, which implies u / ∈ W 1,2 (V), and the proof is complete.
In order to illustrate the theorem presented above we have computed the effective conductivity by a standard FEM program in the classical case, l(r) = 1, k g = k and k w = 1/k (concerning the finite element method in general, (cf.[3,4]or [11]).In this case λ eff = 1 (see (2.5)).We have made an element mesh as shown in Figure 3.1, with increasing number of elements close to the midpoint of V.The total number of nodes and elements (quadrilateral 8-nodes elements) in this triangulation are 28207 and 9350, respectively.Even with this large number of elements we clearly see from Table 3.1 that λ + eff /λ eff diverge rapidly as k increases.This agrees with the theoretical result presented in Theorem 3.1.It is important to note that not only the number of elements is critical, but also the choice of the refinement.Our choice of triangulation is based on the fact that the gradients are large near the midpoint of V (see the calculated temperature distribution on V, Figure 3.2).
In Table 3.1 we present the results for some values of k.

An improved numerical method
We will, in this section, describe a method to overcome the difficulties which arise when using the minimizer of (2.6) as an approximation of the minimizer of (2.2) for a finite-dimensional subspace W .This method will guarantee uniform convergence in the energy norm with respect to k w /k g when the elementsize tends to zero.
Let 0 < α ≤ 1.We insert a disk O in V with centre in 0 and radius R (see Figure 4.1).The disk O is divided into equally sized sectors {S j } with angle φ 0 , on which we define the functions {t j } of the form where Moreover, on the triangular elements {K j } we define the polynomials {p j } of the form The coefficients of {t j } and {p j } are coupled in such a way that these functions agree on common traces.For a fixed triangulation with given R, h and φ 0 , where h is the longest side length of the triangular elements {K j }, we let W h denote the space of functions u such that u = v − x 1 ∈ W, where v is antisymmetric (v(x) = −v(−x)), of the form v = t j on S j and v = p j on K j .Now (2.6) takes the form 164 A new method for numerical solution of checkerboard fields For a given u where In the case when l(r) = 1 for 0 ≤ r ≤ R we obtain that K 1 = K 2 = 1.We note that the integrals in (4.6) are easily calculated.For example where each of the integrals can be found by using standard integration formulae for trigonometric functions.Thus becomes a quadratic form in the coefficients a j ,b j ,c j ,d j ,e j , and f j .Summing up becomes a quadratic form in the collection of all such coefficients.The minimizer of (4.4) is therefore found easily by standard numerical treatment.
In order to find an approximation of the minimizer we may solve the problem for a number of α-values, α ∈ 0, 1], and search for the α-value which gives the minimum value of However, the following theorem shows that the value is always a good choice.
Proof.Assume that k w /k g ≤ 1 (otherwise we just use the same arguments for k g /k w ).Without loss of generality, put k g = 1, and First we want to prove that there exist a > 0, and u ∈ W such that for every 0 < k ≤ a. Next, we will show that there exist R, φ, h, and u ∈ W h , such that (4.14) holds for k ∈ 0, a] and such that and let g(θ) = −g(θ − π) for π/2 < θ < 3π/2.Let Q R,g be the continuous function defined by Moreover, in the grey part of V\O we let Q R,g = 1 and −1 (on opposite quadrants), and let Q R,g be linear in x 1 in the white parts of V\O.For the function u = Q R,g − x 1 ∈ W, we have that

.16)
166 A new method for numerical solution of checkerboard fields We observe that (4.17) We obtain Du + e 1 2 rdrdθ, ( Du + e 1 2 rdrdθ where Using the formula we see that the right-hand side of (4.18) and (4.19), denoted F(g) and H(g), can be written as Stein A. Berggren et al. 167 respectively.By some further calculations we obtain that Observe that α → 0 as k → 0, which by the continuity of l at 0 implies that c + α and c − α → l(0).Since α α → 1 as k → 0, and we can conclude that The function Q R,g can be approximated by Q R,gn in W h , where Here, θ i = πi/(2n) for i = 0,1,...,n − 1 and g n (θ) = −g n (θ − π) for π/2 < θ < 3π/2 (note that sin(θ − θ i ) = sin θ cos θ i − sin θ i cos θ).Let φ = θ i+1 − θ i = π/(2n), then it is easy to see that (4.28) Using (4.28) in (4.22) and (4.23) we obtain from (4.25) that for every ε > 0, there is an N such that where We want to show that for any closed subset P ⊂ 0, 1], of the positive real numbers, there exist R, φ, and h such that,

.30)
Put k − = min k∈P {k} and put Choose a finite subset J of P such that for every k ∈ P there is some y ∈ J, such that k ≥ y and k − y = ∆t < δ/(x − A + ).Choose R, φ and h so small (using standard theory for error estimates in the finite element method, compare, e.g., with [11, page 97]) that Let k be some arbitrary point in P and let y ∈ J be such that, It follows that Hence, G(k, v k ) − G(y, v y ) < δ.Since the function k → 1/ √ k is uniformly continuous in P, there is a δ > 0 such that holds for every k ∈ P, hence from (4.31) and (4.36), we conclude that (4.30) holds. 1].Again by (4.29), there is a number k 2 , 0 < k 2 < k 1 , and (5.2) It is easily seen that (5.4) Setting we obtain from (5.5) that (5.7) Remark 5.1.We recall that the exact value λ eff = 1.In our example, we only use 4 elements outside the disk O. From (5.7) we obtain that λ + eff,h ≤ 1.12783 when k = 100.Note that the corresponding value from the standard numerical treatment (with 9350 elements) in Section 3 is 5.57028 (see Table 3.1).This illustrates the power of our new method.Remark 5.2.In the special case of standard checkerboards we have only compared our method with calculations obtained from using the finite element method (FEM).However, there exist methods based on solving integral equations numerically which are significantly better than standard FEM in case of difficult geometries.Here we want to refer to a recent work of Helsing [8] which considers checkerboards in which the white squares (with conductivity = 100) are a little bit smaller than the darker ones (with Stein A. Berggren et al. 171 conductivity = 1).In particular, by his method he has calculated the effective conductivity to be 9.89299 in the case when the white squares occupy 49,9999999% of the whole structure.According to Helsing it is also possible to apply a modified version of his method in case of standard checkerboards.

Some final comments
We think that the results obtained in this paper will be useful for mathematicians and engineers dealing with checkerboard composites.Our numerical experiments show that it is difficult to obtain good approximations by using the finite element method, and Theorem 3.1 (see also Remark 3.2) even states that there exists no finite-dimensional function space which can approximate the actual solution (independently of the material properties in the checkerboard).
By the new method we overcome this obstacle (see Theorem 4.1).It seems to be possible to extend the method to 4-phase checkerboards.We aim to develop these ideas in a forthcoming paper.

Figure 2 . 1 .
Figure 2.1.A part of a checkerboard with a quarter of a period, denoted V , shown to the right.

Figure 3 . 2 .
Figure 3.2.Distribution on V of the function v(x) = u(x) + x 1 , where u is the (FEM) solution (v = −1 on the left boundary and v = 1 on the right boundary).