UFDC Home  myUFDC Home  Help 



Full Text  
MULTISCALE ANALYSIS OF PARTIAL DIFFERENTIAL EQUATIONS MODELING VOLTAGE POTENTIAL By YERMAL SUJEET BHAT A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2006 Copyright 2006 by Yermal Sujeet Bhat To Y.L. Bhat, 19361994. ACKNOWLEDGMENTS My deep gratitude to Y.L.B., L.L.B., V.L.B., and S.L. Moskow and the University of Florida Mathematics Dept. I would also like to thank all my family and friends for their help and support. TABLE OF CONTENTS page ACKNOWLEDGMENTS ................... ....... iv LIST OF TABLES ........... ...... ........ ...... vii LIST OF FIGURES ............. .... ............ viii ABSTRACT ........... .. ............... .... x CHAPTER 1 INTRODUCTION .............................. 1 2 THE ELECTROLYTIC VOLTAGE POTENTIAL MODEL ....... 5 2.1 ButlerVolmer Boundary Conditions ................ 5 2.2 Existence and Uniqueness .......... ............ 8 2.3 Regularity . . . . . . . 15 3 A CORRECTOR BASED ON NEUMANN BOUNDARY DATA ..... 18 3.1 A Neumann Boundary Condition ........ .......... 18 3.2 Finite Element Method Implementation ..... ... ........ 24 3.3 Numerical Results .. .............. ........ 31 4 A CORRECTOR BASED ON ROBIN BOUNDARY DATA ....... 42 4.1 A Robin Boundary Condition ........ ......... .. 42 4.2 Numerical Results ........................... 46 5 SHIFTING MATERIAL BOUNDARIES ............ .... 58 5.1 The Electrostatic Conductivity Model ............ 58 5.2 Estimating the H'(Q) norm of Uo uh ............ .. 59 5.3 Formal Asymptotics ........... ............. 63 6 CONCLUSION ...................... ........ 69 REFERENCES ...... ...................... ...... 71 BIOGRAPHICAL SKETCH ....... LIST OF TABLES Table page 31 Table of estimates over 0 and convergence rates . . .... 32 32 Table of estimates over F and estimates of the gradient over .... 32 41 Table of estimates .................. .......... 47 LIST OF FIGURES Zinc loses electrons to silver . ............ The base is a heterogeneous surface . ........ 22 Perimeter increases while anodic area fraction stays constant. Twodimensional analogue .. ........ Limiting behaviour of u, on F as e approaches (b) = 1/11, (c) = 1/25, (d) = 1/40 u,, = 1/5 . . . . . uo + EU('), e = 1/5 ........... u, E = 1/11..................... Uo + C', E = 1/11 .. .......... The potential on the boundary F, e = 1/5 The potential on the boundary F, e = 1/11 The potential on the boundary F, e = 1/25 The potential on the boundary F, e = 1/40 L norm of Vu, on F as e approaches 0.. u e = 1/5 . . . . . uo + e U), e = 1/5 .. .......... u 6, e = 1/11 ...... 0 + E ', e 1/11 . . .... Graph of u, (above) and uo + u1) (below) on The approximation and the original, e = 1/5 47 Graph of u, (above) and uo + u1) (below) on zero for: (a) e the boundary the boundary 48 The approximation and the original, e = 1/11 . . Figure 11 21 page 3 6 1/5, = 1/5 = 1/11 r, P, < . . 53 49 Graph of u, on the boundary F, e = 1/25 410 Graph of uo + u1) on the boundary F, e = 411 Graph of u, on the boundary F, e = 1/40 412 Graph of uo + ui) on the boundary P, e = 51 Perturbation due to shifting between two 1/25......... = 1/40 . . . dielectrics a1 and 2..... Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy MULTISCALE ANALYSIS OF PARTIAL DIFFERENTIAL EQUATIONS MODELING VOLTAGE POTENTIAL By Yermal Sujeet Bhat May 2006 Chair: Shari Moskow Major Department: Mathematics We study a nonlinear elliptic boundary value problem arising from electro chemistry. The boundary condition is of an exponential type. We examine the questions of existence and uniqueness of solutions to this boundary value problem. We then treat the problem from the point of view of homogenization theory. The boundary condition has a periodic structure. We find a limiting or effective prob lem as the period approaches zero, along with a correction term and convergence estimates. This correction term satisfies a boundary value problem with Neumann boundary conditions. We do numerical experiments to investigate the behavior of galvanic currents near the boundary as the period approaches zero. We then consider a correction term which satisfies a boundary value problem with a Robin boundary condition. We do numerical experiments to investigate our approxi mation based on this corrector term. We then use asymptotics to analyze the behaviour of the steady state voltage potential of a conductor with a region of in homogeneity. The boundary of the inhomogeneity shifts slightly, we do asymptotics utilizing the small shift. CHAPTER 1 INTRODUCTION Modern perturbation theory is primarily concerned with constructing ap proximations of solutions to mathematical models that have a parameter which is approaching zero. One such class of models are boundary value problems in which the domain has a periodic structure. In this case the period size is the small scale parameter. These types of boundary value problems often arise in the study of, for example, composite materials, macroscopic parameters of crystalline struc tures, fluid mechanics and aerodynamics. Perturbation theory is an example of an analytical approximation as opposed to a numerical approximation. There are many techniques to formulate these analytical approximations. One fundamental technique is through the use of a multiple scales asymptotic expansion. The foun dational ideas for this approach appear in the early 1800s. In 1812 Laplace used asymptotic series to analyze some special functions. In 1823, Poisson constructed an expansion of a Bessel function. In 1886 Poincar6 used asymptotic expansions to study solutions of differential equations. The idea behind the asymptotic approx imations is that the solution can expressed as a sum of terms of different orders of magnitude with respect to the small scale. For example if e,(x) is the solution to some boundary value problem with small scale parameter E then we begin by assuming the solution 0, has an asymptotic expansion of the form o 0o + E0!1) + E20!2) + . The general procedure is to then substitute the above expansion back into the original boundary value problem to determine associated boundary value problems for o, 01), 02 ), ... We try to find simpler equations that describe the behaviour of the solution on various orders of e. Here we wish to use multiple scales analysis to develop an asymptotic expansion of the solution to some phenomena related to conductivity. We first use perturbation theory to approximate the solution of a nonlinear boundary value problem which models galvanic corrosion. As a second application of asymptotics we wish to describe the conductivity properties of a material with a shifting dielectric boundary In the electrochemistry community there is much interest in the study of galvanic interactions on heterogeneous surfaces [12], [13]. When two different metals in electrical contact, referred to as anode and cathode, are immersed in an electrolytic solution, the difference in rest potential generates an electron flow. This electron flow is called a galvanic current and may lead to a deterioration (corrosion) of the anode. In Figure 11 a strip of silver (Ag) and a strip of zinc (Zn) have been im mersed in a saltwater solution. The zinc strip gives up electrons to the silver strip. The silver strip is said to be cathodic, and reduction takes place (Ag gains elec trons.) Simultaneously oxidation takes place at the zinc strip, zinc loses electrons, and is said to be anodic. Zinc dissolves into the solution, the zinc electrode is being corroded and the electron flow is known as galvanic current. The driving force of the electron transport process is the difference in potential of the two metals involved. See Newman [15] for a complete introduction to the subject. Here we study the electrostatic problem on a surface where anodes are arranged periodically in a cathodic matrix. Mathematically the potential is modeled as a function, 0, over a Euclidean domain Q. Part of the boundary of Q is electrochemically active while the rest of the boundary is inert. It is the active region of the boundary that is made up of anodic and cathodic portions. The potential over both of these regions satisfies an exponential boundary condition of Butler and Volmer, but with different material parameters on each portion. Ag Zn hydrogen gas NaCI ^H ^ e ^ Figure 11: Zinc loses electrons to silver In Morris and Smyrl [12] the authors study such a problem numerically, using finite elements. Additionally various interesting aspects of the twodimensional, homogeneous model with the ButlerVolmer condition have been analyzed in Bryan and Vogelius [5], Turner and Hou [9], and Vogelius and Xu [19]. To the best of our knowledge, however, studies coming from the applied mathematics community have been restricted to two dimensions. The main reason for this is that one can bound exponentials of the twodimensional weak solution on the boundary by using an Orlicz estimate [18], [19]. Such an estimate would require more than H1 regularity in higher dimensions. In this paper, we attempt to treat a periodically heterogeneous problem, in two and three dimensions, from the point of view of homogenization theory. Our second application of asymptotics is to a problem pertaining to electro static conductivity. We consider a PDE which models the steady state voltage potential of a metal with a small inhomogeneity. There is a jump in the conduc tivity across the boundary of the inhomogeneous region. The boundary of the inhomogeneous region shifts by some small amount due to some type of physical stress. We now wish to describe the voltage potential of the conductor with shifted boundary as a perturbation of the voltage potential of the original conductor. We do asymptotics to establish an estimate. In Chapter 2 we formally present the electrolytic voltage potential model. We tackle the issue of existence and uniqueness of the solution to the model and then discuss the issue of regularity. In Chapter 3 we construct an asymptotic ap proximation of the solution to the original problem. Here the second term of the approximation satisfies a linear boundary value problem with Neumann data. We then establish some convergence estimates and do numerical implementation. Using a finite element method approach we implement and test our asymptotic approx imation and convergence estimates. In Chapter 4 we propose an approxiamation in which the second term satisfies a boundary value problem with Robin boundary data. We then do numerical implementation and testing of this approximation. In Chapter 5 we perform asymptotics on the "shifting boundary" problem. We establish some convergence estimates and do some formal asymptotics. In Chapter 6 we discuss future work to be done. CHAPTER 2 THE ELECTROLYTIC VOLTAGE POTENTIAL MODEL 2.1 ButlerVolmer Boundary Conditions Now we formally present the threedimensional model for electrolytic voltage potential on a heterogeneous surface. The domain 0 is of cylindrical shape with base some twodimensional domain. The bottom base is assumed to contain a periodic arrangement of islands (anodes). We call this collection of islands 8QA and the remainder of the bottom of the base tQc cathodicc plane). The electrolytic voltage potential, 0, satisfies the following nonlinear elliptic boundary value problem, A0 = 0 in Q JA [eas(4V) V eaac(4A)] on aA (2.1) 8n S= Jc[ea"(vc) e a vc)] on aQc on = 0 on 9 \ {tQA U 8Qc} 8n where aaai, c a, ce, are the transfer coefficients and it is assumed that the sums (aaa + aac) and (ac. + ace) are equal to one. The positive constants JA, JC are the anodic and cathodic polarization parameters and VA, Vc are the anodic and cathodic rest potentials respectively. Note that V4 represents galvanic current. These boundary conditions are the socalled the ButlerVolmer exponential boundary conditions. In the numerical studies of [12], the authors observed that for fixed ratios of anodic to cathodic areas on the bottom base, the resulting current increased ap proximately linearly with the length of the perimeter between the two regions, and Figure 21: The base is a heterogeneous surface they hypothesized that it is the ratio of anodic area to perimeter that determines the size of the resulting current. Figure 22: Perimeter increases while anodic area fraction stays constant. As a special case of increasing perimeter with approximately fixed area fraction, we consider a periodic structure with period approaching zero. Our goal is to expand the solution asymptotically with respect to the period size. Convergence results involving these approximations could provide insight into the behavior of the current for small period size; and possibly lead to techniques for computing approximate solutions to (2.1). We model the periodic structure by letting f (y, v) = A(y)[e"(y)(vv(y)) e(1"(Y))(vv(y))] for any v R and y E Y, the boundary period cell, which for simplicity we take to be the unit square; Y = [0, 1] x [0, 1]. Here A, a, and V are all piecewise smooth Yperiodic functions. We also assume there exist constants Ao, A0, ao, Ao and Vo such that: 0 < Ao < A(y) < Ao, (2.2) 0 < ao < a(y) < Ao < 1, (2.3) and V(y)I < Vo. (2.4) See [5] and [19] for an analysis of when A < 0. Consider the problem Au = 0 in Q S= f(x/,,u,) on r (2.5) On = 0 on 0Q \r. On As is typical in homogenization problems, one expects that as c  0, the solutions will converge in some sense to a solution of a problem with an averaged boundary condition. Define fo(v) to be a cell average of f(y,v), that is, fo(v) = f(y,v)dy. Consider the candidate for the homogenized problem Auo = 0 in Q Buo =u fo(uo) onr (2.6) an Buo = 0 on dQ \r. On Remark If, as is the case in [12], Y = Y UY2 and the functions A,a,V are piecewise constant, each taking on the values A,, ai, V respectively in Y5, then fo(v) = I1 jI[eal(1) e(l")(vvi)] + YA2 [ 2(VV2) e(1 2)(V2)]. That is, the above homogenized boundary condition would depend on the volume fraction of anodic to cathodic regions. 2.2 Existence and Uniqueness In this section we show that the energy minimization forms of the nonlinear problem (2.5) and (2.6) have unique solutions in H1(Q) in any dimension. Some elements of the proof are similar to those in [9] and [19]. For a given E, define the following energy functional, ,(v) = Vv 2 dx + F( v)da (2.7) where, F(y, v) = A(y) eY)(vV(y)) + ) e(1a(y))(vV(y)) a(y) 1 a(y) We show the existence and uniqueness of a minimizer of (2.7). Formally, we show the existence of a function u, E H1(V) such that E,(u) = min E(u). (2.8) Note that E, is not necessarily bounded on all of H1'() (unless n = 2 for which we can use an Orlicz estimate). However this does not pose a problem. We set E, equal to (2.7) where it is well defined and to +oo where it is not, as in [7], p.444. In the twodimensional case of the model, due to the boundedness of E, on Hl(S), direct calculation shows u, satisfies the variational form of (2.5), SVu Vv dx = f(x/e, u)v dc for any v c H1'(). (2.9) In the threedimensional case, if u, is an energy minimizer, we will have that / F(x/e,u,) da < oo, (2.10) and hence by the positivity of each term of F(x/e, u,), we have that each term is separately in L1(Q). Therefore, E,(ue + tv) < 00, (2.11) for any t E R and for any v which is smooth on P. Standard arguments then show that u, satisfies, / Vu. Vv dx f(x/e, ue)v doa for any v E C0(Q). Additionally, if we knew that u, E CO(Q) then f(x/e,ue) is bounded and hence clearly in H1/2(P). So by the density of C"(0) functions in HI(Q), u, in this case would satisfy / Vue Vv dx f(x/e, u)v doa for any v E H (). (2.12) Consider also the functional Eo(v) = V12 dx + Fo(v) do. (2.13) where, Fo(v) = Fj F(y,v)dy. Here again the energy Eo is not necessarily bounded but as before, we set Eo equal to (2.13) where it is well defined and to +oo where it is not. Direct calculations show that a minimizer uo of (2.13) will satisfy, / Vu o Vv dx = f(uo)v dao for any v E H1'(), (2.14) assuming u0 is continuous (actually we will see that uo is a constant.) Theorem 2.2.1 (Existence and Uniqueness of the Minimizer). Let E, be defined by (2.7), where X,a, and V satisfy (2.2)(2.4). Then there exists a unique function u, E H1(Q) satisfying E,(ue) = min E,(u). neH' () Proof. Note that 2 F(y,v) = A(y)a(y)ea(Y)(vV(y)) + A(y)(l a(y)) e(1a())(vV(y)) (9v2 since A > 0, a > 0, and 1 a > 0 we have that 92F > 0. Clearly the partial derivative is bounded below. That is, there exists a constant co, independent of y and v such that, 92 2 F(y, v) > co > 0. Since F is smooth in the second variable, for any v,w E H1(Q) and for any y, there exists some ( between v + w and v w such that a02 F(y, v + w) + F(y, v w) 2F(y, v) = Fly, J)w2 which from the lower bound yields F(, v + w) + F(, v w) 2F(X,v) > c0w2 6 E E whence E,(v + w) + E,(v w) 2E,(v) > +f Iw 2 dz + co0 w2cdr > Wco2H111 ) (2.15) where the last inequality follows by a variant of Poincar6. Now let {u'}'l be a minimizing sequence, that is E,(u) > inf E,(u) as n oo. ueHi(n) Since all the terms of (2.7) are nonnegative, clearly inf E,(u) > oo. ueHl(2) Note that without loss of generality we can choose the minimizing sequence so that all terms have finite energy (since infu Hi() E,(u) < E(0) and E(0) is bounded independently of E.) Let un + u2 V 2 and n U W u 2 Then v + w = un and v w = u, so (2.15) implies E(v + w) + E,(v w) 2E,(v) > 4u U7n1 which implies, E(t) + E,(u) 2 inf E,(v) > u u"Jn) veH1 (c) 4 Now if we let m, n  oc, we see that {uf}n is a Cauchy sequence in the Hilbert Space H'(Q). Define u, to be its limit in H1'(). Then we have S u, in H1(Q) which by the Trace Theorem implies, u u, in L2(P) which implies ([17],p.68) there exists a subsequence {u~k}k, which we label {zu}k such that u > us a.e. in P. Since F is smooth in the second variable, and u. > u a.e. in P we have that F(X,u6) = lim F( ,u,) a.e.. E ko0o E Now note that clearly F(Q, u) > 0 for any k. So, by Fatou's Lemma we have, SF(, u~)da, < lim inf F( zu)du,. Thus from this and the fact that u.  u, in H1(Q), we can conclude that, E,(u) < liminfke, EE(u) Hence, E,(u)= inf E (u). So we have shown the existence of a minimizer. Suppose u, and w, are both minimizers of the energy functional, i.e. E,(u) = inf E,(u)= ueHI(n) E, (w,). Now if we let v = (u, + w,)/2 w = (Ue w,)/2 then substituting v and w into (2.15) yields, Ue + w, jo 2 E,(uE) + E(w) 2E,( 2 >) W> H 'I i(n). However, this implies, u wll (n) < E,(ue) + E(w) 2 infuEHi(a) E,(u) = 0. Hence u, = wE in H'(Q). Thus we have shown the uniqueness of the minimizer. O Note that this argument can be generalized to address the ndimensional problem, i.e. the case in which we have f C R",r C R"1 with boundary period cell Y = [0, 1]"1. The existence and uniqueness of a minimizer uo of Eo follows from the same proof. limk_,, E(uk) inf.ueH(n) E,(u). Corollary 2.2.2. There exists a constant C, depending on Ao, ao,Ao and Vo but independent of e such that, IIU'IIHl() < c where u, is a weak solution to (2.5). Proof. Consider the function v 0. Then E,(v) == E(0) = F(, o)d < M for M independent of e (but depending on Ao, o, Ao and Vo ). Then since u, is a minimizer, E,(u,) < E,(0) < M. Since both terms in E, are positive, IIVuE 2( ) < M. We also have that SF(x,u,)da, By examining the form of F(y,v), we see that there exists some constant d, depending on Ao, ao, and Ao but independent of e and x such that duE V(x ) Hence, Su V(x) d\x < M/d, which by the boundedness of V implies that SIU, do, < M where M is independent of E. One variant of the Poincar6 inequality says that there exists C such that \\u, jueda\x\L2(Q) Finally the reverse triangle inequality yields, IIIIL2() < C IIVUEIL2(Q) + M, which proves the corollary. W We conclude this section by establishing the fact that the solution to the homogenized problem (2.6) is in a fact a constant. This fact follows easily once we have established the following lemma. Lemma 2.2.3. There exists a constant K such that fo(K) = 0. Proof. Recall that fo() = j f(y,v)dy (2.16) where f(y, v) = A(y)[e v()) e"))v())] for any v E R and y E Y, where Y = [0, 1] x [0, 1]. Also recall that A, a, and V are all piecewise smooth Yperiodic functions for which there exist constants Ao, Ao, ao, Ao and Vo such that: 0 < Ao < A(y) < Ao, (2.17) 0 < co < a(y) < Ao < 1, (2.18) and IV(y)I < Vo. (2.19) Now if v > Vo then (2.19) and (2.18) imply (1 a(y))(v V(y)) < 0 for all y Y. a(y)(v V(y)) > 0 and So we have e(a(y))(vV(y)) < ea(y)(vV(y)) for all y E Y. Thus (2.17),(2.18), and (2.19) imply that if v > Vo then f(y,v) > 0 for all y E Y. Similarly if v < Vo then f(y,v) < 0 for ally E Y. Now note that (2.16) implies that if v > Vo then fo(v) > 0 and if v < Vo then fo(v) < 0. Now since fo(v) is continuous by the Intermediate Value Theorem there exists a constant K E (Vo, Vo) such that fo(K) = 0. Thus the lemma is established. W Theorem 2.2.4. Let uo be a minimizer of (2.6) then uo is a constant. Proof. Suppose K is a constant such that fo(K) = 0. Such a constant exists by Lemma 2.2.3. Then clearly u0 = K is a strong solution of (2.6). Note that this argument and Lemma 2.2.3 hold in any dimension. O 2.3 Regularity We conclude this chapter with a short discussion of the regularity of the solu tions u, and u0 for the two and threedimensional case. For the twodimensional case of this problem, i.e. when the medium is layered as in Turner and Hou [9], and Vogelius and Xu [19] (see Figure 23), using imbeddings of Sobolev Spaces (W"'p) into Orlicz Spaces (L') we can show that f(x2/e, U) and fo(uo) are bounded in L2(F) independently of c. We first give a general definition of an Orlicz Space before presenting the imbedding result. See [2] or [16] for a thorough discussion of Orlicz Spaces. Let Q2 be a bounded domain in R". Let 4 be a Young function, i.e. a realvalued, convex function such that 4(x) = 4(x), 4(0) = 0, and 4 / oo as x / oo. Then the Orlicz space LQ(0) is the set of all measurable functions f on 2 such that 4 44(af) dx < o for some a > 0. oo soeo X2 1 Qr x, 1 Figure 23: Twodimensional analogue. Now define the norm I f L() =inf k> : j )( ) dx < l then L2(Q) becomes a Banach space. For example when 4(x) = IxlP then L' = LP. We now state the imbedding result established by [18] and [2]. Theorem 2.3.1 (Trudinger's Theorem). Let Q be a bounded domain in R" satisfying the cone condition. Let mp = n and p > 1. Set a(x) = eCPp1) 1. Then there exists the imbedding Wm,P()  Ln(Q). In the twodimensional case of (2.5) we can conclude from Trudinger's Theo rem [18], [19] that there exists a constant C such that for any v H'I() and any real / we have, Jefll~dx2 < ecf2(jlVil(nl')( r + 1). Then from standard elliptic regularity theory this implies that u, and uo are in H3/2(Q), with the norm bounded independently of e. By the trace theorem we then obtain bounds for u, and uo in H'(F). Since P is onedimensional it follows that u, and u0 are continuous on F and bounded pointwise, and their tangential derivatives are bounded in L2(F) For the homogenized solution we have much more regularity, u0 is in fact the constant that satisfies fo(uo) = 0. For nonzero boundary conditions on the inactive region, uo would still be a smooth bounded function. So for the twodimensional version of this problem we have the following lemma: Lemma 2.3.2. If Q C R2 is a rectangle and F is an edge, then u, e C(Q) where u, is a weak solution of (2.5). Furthermore, there exists a constant D, the value of which does not depend on e, such that, I6U(x) lc(n) < D. In the threedimensional case of (2.5) since u, E H'I() implies m = 1, so mp $ n and thus Trudinger's Theorem does not apply. In fact there is a large class of Sobolev imbedding theorems pertaining to the case mp = n, however there are no imbeddings of HI(Q) for a cyindrically shaped domain 2 C R3 that are applicable to (2.5). In general, there has been no rigorous analysis of the regularity of threedimensional solutions of elliptic boundary value problems with L1 Neumann boundary data. Now if we assume that (2.5) models electrolytic voltage potential then it is physically reasonable to make some assumptions about the regularity of u, in the threedimensional case. In the next chapter we show that if these physically justifiable assumptions are made we can establish convergence estimates for our asymptotic approximation for the threedimensional case. CHAPTER 3 A CORRECTOR BASED ON NEUMANN BOUNDARY DATA 3.1 A Neumann Boundary Condition To show u, converges to uo we will add a correction term and prove estimates in terms of powers of E. The convergence of u, to uo when n = 2 will then easily follow from this. The same estimate holds when n = 3 if we know that the solutions are continuous and uniformly bounded. When n = 2 we will see that the convergence is strong in HI(Q) and of the order of vE. Let u0 be a minimizer of (2.13) and define the correction ut4' to satisfy, Au)) = 0 in 0 a() 1 I u = (f(,uo) fo(uo)) + e, on (3.1) An 6 E uc 1) 6 = 0on 80 \r On I U1) d3 = 0 (3.2) where, ee = (fo(uo) f(x/,uo))dax. Hence e, is chosen such that the solution always exists and the condition (3.2) guarantees this solution is unique. We note that given uo, this is a linear problem. Now if u, and uo are in L(P), let De =max l IUL (), '0 o L(r)} (3.3) and let, M= sup (y, w). (3.4) (y,w)EYX[De,De] ov The next estimate holds for dimension n = 2 or 3 but depends on the constant Me. We do not know a prior that D, is finite in general when n = 3. However, such an assumption seems to be physically reasonable and known to be the case when the medium is layered. Proposition 3.1.1. Let n = 2 or 3 and let u,, Uo be minimizers of (2.7),(2.13) respectively, and let u(') be the solution to (3.1). Assume also that u, E Co(Q). Then there exists constants C and D independent of e such that U6 Uo eu Hi() < Ce(M, + D), where Me is defined by (3.4). Furthermore, there exists constants C1, and C2 independent of E such that, IU)L2:(r) < C1 and e6j < C2. Proof. Let ZC = UE UO ECU since u, is continuous, by (2.12), we have that for any v E H1(Q), Vze Vv dx= Vu Vv dx V Vv dx e Vu Vv dx j f(,uI)vdux + f( ,uo)vdorx +e evdox. So, Vze Vv dx + [ ( , e) f (x,o)]vdox c ecvdo =0. Now note that Uo and u, are defined pointwise on P. So, by the Mean Value Theorem, for each fixed c and x E F there exists between uo(x) and u,(x) such that, xx f x . f ( ,u) f( ,uo) = (u, uo) ( ,(X ). E E a2) E By subtracting and adding Eu() we have, f (,) f o) = Z (, + E ( af which, if we pick v = z,, yields, SVzVI dx+jf z29f(x, j ()dd d=c M 4)( C)dX + zd . Since f > co, this implies C~I J() < 4 z62dx+jzfi(,da = Me af (8 (X)z dax + ee zcd&e . So by applying Hilders Inequality and then the Trace Theorem we have, C 2 f< 1 )1/2 _< ll (t ,( e ) ll ( ) l (l ) L2() + l f i l/2) ~ L2 0 Io  I r, 11/2) 1 H (9~I~v E~~~lO~) U' L() e r/)Zlx~, Thus, we can write, tif x Hll'Z( ) < C ( vL() E 1). Now recall for any v we have, (f(y,v) fo(v))dy= 0 so there exists a continuous Yperiodic function r(y,v) such that Ary(y,v) = f(y,v) fo(v),Vv R. So we have, (3.5) (3.6) ee = (fo(uo) f(,uo))dar = Ayr(X, uo)dox p Jr 6 = Vyr(,uo) v dsx ar C where the last equality is arrived at using integration by parts and the fact that the 8r dr chain rule implies (x/, Uo) = (x/e, Uo). Note that the differential operators ay ax Vy and A. are with respect to y E Y, that is, they are surface operators. Now since u0 is bounded pointwise on F and since r(y, v) is a continuously differentiable Yperiodic function we have, ee < C (3.7) where C is bounded independent of c. Now we show that I.I')lIL2(r) is similarly bounded. Let w, E H1(Q) satisfy, Aw = 0 in 0 = uld( on r (3.8) On Owo = 0 on 9Q \ r 8n SWe do, = 0 then, (1) wd j=z = V Vw dx wdax fE (t 9n C2 e an C where the last two qualities follow from integration by parts. Now, since u(1' satisfies (3.1), we have f au ', ff(x/e, uo) fo(Uo) d weda = f + e[ w6d aQ tn E' I X l= Ayr(x/e,uo)wedcx e, w'da 1 Ayr(x/e, Uo)wdax e Jr where the second equality follows from (3.6) and the last equality holds since / wda, = 0. Now using the chain rule we can write Ayr(x/e, uo) = 2 Ar(x/E, o), where As is a surface Laplacian on F. Thus, we have, j(u())2dox = Ac Ar(x/e,uo)w da r Jr = e V7r(x/C,Uo)VwduE WEdsx, (3.9) where v is the outward unit normal to &F. Note that when n = 2, we use the last integral to represent endpoint evaluation. So, by Hilders Inequality, ef Vjr(xa/e,Uo)Vwcedax JE r wds, < Ell Vc r L2(r) VWE L2(r) +E IL2(ar)WllWll2(ar). (3.10) Then by the Trace Theorem we have, IIW6e 2(ar) < C1 IWe IH (r) < C2w~ e H3/2(). (3.11) Similarly, IIVW IL2() < C311 W6 Hi(rF) < CllEllH3/2(). (3.12) Then (3.9), (3.10), (3.11) and (3.12) imply, u 2() < c ( VXrL2(P) +  L2(P)) W H /2(). Now since we satisfies (3.8) we have from standard elliptic regularity theory [10], I HWll3,/2(n) < CU 1) llL2(r) where C is independent of e and so we can write, I LJ() JI2(r) < Ce  V r(x/e,U o) 2(rr) /E,~O IL2( = ( Vr(x/e, uo) L2() + 11r uo) L2( , where the last equality follows from the chain rule. Consequently, since we have that u0 is continuous on P and bounded pointwise and since r(y,v) is a continu ously differentiable Yperiodic function we can conclude that I4u) L2(r) < D (3.13) where D is bounded independently of e. Then (3.5), (3.7) and (3.13) imply the main result of the proposition: Of x KZ H'U() CE(Il e LP) l U L2(Plr) + e1 ) < Ce(M+ D), where M, is defined by (3.4). O Note that in light of Lemma 2.3.2, we can easily establish the following corollaries: Corollary 3.1.2. When n = 2, i.e. for the case in which f c R2, F c R with boundary period cell Y = [0, 1] there exists a constant C independent of e such that, Iu, Uo eu,1) Hi(n) < Ce. (3.14) Corollary 3.1.3. When n = 2, for u, the weak solution of (2.5), and uo the weak solution of (2.6), there exists a constant C independent of E such that, IUE uo \Hi(n) < C (3.15) Estimate (3.15) follows from the fact that, U("n) < C ll HH/2() < CC1/2, where the last inequality follows by interpolating between L2(P) and Hl(P) (see [11], Section 11.5) and then using duality (as in [14]). Finally note that estimate (3.14) also holds for n = 3 if we know that D, defined by (3.3) is uniformly bounded. 3.2 Finite Element Method Implementation We wish to numerically observe the behaviour of the homogenized boundary value problems as a way to describe the behaviour of the current near the bound ary. We use a finite element method approach to the 2D problem. For the 2D problem the domain 0 is a unit square and the boundary P is the left side of the unit square, that is r { (x21,2) E 1 =: 1} (see Figure 23). In this case we impose a grid of points(called nodes) on the unit square and triangulate the domain, then introduce a finite set of piecewise continuous basis functions. Here we use standard finite element method "tent" functions. We impose a grid of node points on the unit square which are evenly spaced both on the x and yaxes. We label the nodes starting at the origin and moving to the right. We label the nodes Pi, i = 1,..., m, where P1 = (0, 0) and P, = (1,1). Furthermore, if there are N nodes on the axis (i.e. if the axis is divided into N 1 pieces) then m = N2 and PN1 = (1 1/N, 0), PN = (1, 0), PN+I = (0,1/N), PN+2 = (1/N, 1/N) etc. Once the node points have been established we triangulate the domain in a predetermined fashion. Note that the node points P1, P2, PN+I and PN+2 form a square. We form one triangle by using as vertices the node points P1, P2 and PN+I and another triangle by the node points P2, PN+I and PN+2 Proceeding in this fashion we can form yet one more triangle by using as vertices the node points P2, P3 and PN+2 and yet another by using P3, PN+2 and PN+s. The entire domain can be triangulated in this fashion using the node points as vertices. Once the domain has been triangulated we introduce basis functions. Note that in the original problem (2.5) we attempt to minimize the energy functional E,(v)= j Vv2dx + F(X, v)dOx over the space H1(Q). However, numerically we minimize the energy functional over the space Vh where, Vh = {v E C(Q) : v is a linear function when restricted to each triangle in Q}. Now we introduce a set of basis functions bi(xi,X2), i = 1,... ,m. We use the standard finite element method "tent" function, that is for each i define bi E Vh by, bi(Pj)= I ifi= for j= 1,...,m. 0 if i j m So any function v E Vh can be written as v(x) = vjjbj(x), where q = v(Pj). j=1 Furthermore, we can reformulate the energy minimization problem in the following way, if we denote the minimizer of the energy functional as uh then m E,(uh) = min E,(v) = min E(E rjbj(x)). veVh i1,,/rm L j=1 m So if (1,..., m minimizes the energy functional then we can write uh Z= ( jbj(x). j=1 Thus, by triangulating the domain and introducing basis functions of Vh we are able to discretize the problem. Now we wish to solve the energy minimization m problem. Note that if we let v(x) = y rl2jb(x) then the energy functional has the j=1 form, E( T m EE() = A J+ F(x, 7bC(x))riT. j=1 where A = (aij), 1 < i, j < m such that aij = f Vbi Vbj di and ry is the m x 1 vector 7r = (qi,... ,rm)T. Hence we wish to minimize the energy EE(r7) over all ] E R". Since we plan to use a gradient descent based optimization method to minimize the energy we are required to calculate the gradient of the energy as well. For the homogenized problem (2.6) we wish to minimize the functional Eo(v) = I Vv2dx + Fo(v)d 2 JQ Jr over the space H1(Q). Thus numerically we wish to minimize m Eo(0) = A + Fo( yjbj(x))dax. j=1 over all 97 E R". To conclude this section, as an example, we include some detailed calculations to determine the boundary integral I Fo(uh)dx2 explicitly as a function of (. Then we find the gradient of the boundary integral with respect to (. Note that SNl 1 Fo(uh)dx2 = f Fo(uh)dx where i=1 r P i < x2 < N and, J Fo(uh)dx2 = ji F(v, h)dy dx2 In the twodimensional case Y = Y1 U Y2 = A1, (2 = A2, if Ye if Y2 , a(y) { 1, a2, S ea(hV) + (1a)(V)dy dx2 [, k] U (k, 1] and furthermore, [0, 1/k] U (I/k, 1] and furthermore, if y E Y if y Y2 ,and V(y)= , V2, if ye Yi if E Y2 Hence, F Fo(uh)dx2 . r A lea,(uhVi) eA (1ac)(uhV1) kr ka k(1 a ) (k 1)A2 e"2((2) + (k 1)A2 e 2)(Uhv2) dx2. ka2 k(1 a2) Now, note that Uh N2 S(p and uhF j=i N SiNYPiN. Also note that i=1 (N 1)x2 + on Fi (N r = IN(x2)= { and for i = 2,..., N 1 we have PiN r = OivN(2)= (N 1)x2 (i 2), (N 1)x2 + i, on Fr_1 on ri otherwise 1)2 (N 2), on rN1 otherwise This implies uh ri (iNV)iN + (i+l)N z(i+l)N for i 1,...,N 1. Whence, Uh rI = p iN1iN + (i+1)NN)(i+1)N = {iN((N 1)Z2 + i) + i(i+l)N((N 1)x2 = ((i+I)N iN)((N 1)x2 i) + (i+ )N = Gi(x2). (Q(i+l)N iN)((N 1)2 i) + (i+l)N, if iN # C(iS+)N ' (+i)N, if JiN = '(iS+)N Thus, SFo(uh)dx2 Jr ir kVl k(1 ) (k 1)A2ea2(Gi(X2)2) + (k 1)A2 (12)(Gi(X2)V2) dx2 ka2 k(l a2) otherwise and finally (i 1)) So, if iN ~(i+t)N then, F (uh)dx2 = l i ea1([(( +t)N CN)((N1)x2i)+ i+tN]vl ) 1 A e(1a )([((i+ )N IN)((N1)x2i)+2 ( i+)oN]V) i J 1h ka2 k (1 o1) ) + (k ka1)A ea2([((i+l,),N n )((N1)2i)+(E+l)N] V 2) (k 1) A2 e_(I )([(12 (i+ )INg )((N_I) 2i)+ (i+I)N]V2) dx2, k(1 a2) and, if iN = (i+t)N then, SFo(uh)dX2 ( ) A1 (1at)((i+1)NVi)}d2 rr ka, k(1 l) S (kr )A2 e2 + 2) ( 1)A2 e(1a2)((j+)Nv2) }d2 r ka2 k(1 a2) Recall Pi = {x rl i < X2 < }, so if iN # (i+I)N then for i = 1,...,N 1 we have, f A1 ead+tw ) eatiW v) SFo(u)dZ2= 2 AT   c r ka(N 1) (i+1) iN A1 e(1at)((i+1)Nvi) e(1a ()(iN vi) k(l ai)2(N 1) (+il)N iN (k 1)A2 e2((i+l)NVV2) ea2(iNV2) ka2(N 1) (i+i)N (iN (k 1)A2 e(102)(C( +)V2) e(1a2)(jNV2) k(1 a2)2(N 1) (i+1)N iN and, if iNv = (i+)N then for i = 1,..., N 1 we have, fI Ui2 17 A eA((i+(1)V) A (a)(i+) i) F, kl(N 1) k(1 a1)(N 1) (k A 1))A2) ( i1A2, e 1 a)(g e) V2) ka(N 1) k(1 2)(N 1) N2 Recall that Uh = (jp so that Eo(uh) = E0(1,... ,N2). But uhl j=1 iNV)iN + (i+i)NVb(i+1)N for i = 1,..., N 1 implies that the boundary integral is a function of N,(2N,*... (N1)N,(N2 only, i.e. we can write j Fo(uh) = g9(N, 2N, >(N1)N,N2). So we see that the gradient of the boundary integral with respect to ( is VC Fo(Uh) = (0,...,0, 0,...,0, 0,......,0, ). Ir N t 2N adN2 Now note that, g ar f 8g a (9iN 9 iN L _ 89 8 iN{ Jr. and g { (9&2 ( N9 & Thus if &N # (2N then, 89 { Fo(( A1 ea 1(2N ka (N 1) Fo(&N, 2N)dX2)}, ;N i) Fo({(i1)N, iN)dX2 + I Fo(iN, (i+I)N)dx2} for i = 2,..., N 1, Fo(%(N1)N, AN2)d2 }. , 2N)d2 =  eal(ENV) al(2N &N)eal(ENVi) (b2N N)2 A1 g(1a)(g2NV) (1 a)(NVt) + (1 a )(2N N)e(1at)(NV1) k(l a)2(N 1) (2N N)2 (k 1)A2 e"2(2NV2) e2(CV2) 2(N N)e2(1V2) ka2(N 1) (J2N (N)2 (k 1)A2 e(12)(E2NV2) e(1a2)(NV2) + (1 _ )(2N N)e1 2)(CNV2) k(1 aa2(N 1) (2N 2N)2 and if (N1)N # 4N2 then, N zag N(9 Fo((N1)N',N2)dx2} = A1 eaI(t(N1)N V) eai(N 2Vi) + a 1(N2 (N_1))eat(N2Vt) k(N 1) (N2 (N 1))2 AI e(1c)(e(N NV) e(lac)(E2 V) (1 ac)(CN2 i(N1)N)e(1a1)(12V1) k(1 )2(N 1) (NN2 (N_)N)2 1 V (k 1)A2 e"2((N 1NV2) ea62(e2V) + 02(N2 (N_)N)ea2(NV2) ka2(N 1) (N2 (Ni))2 (k 1)A2 e(1"2)((1)V2) e(102)(jN2V2) (1 Q2)(N2 $(N1)N) (1 2)(e22i) k(1 a2)2(N 1) (N2 (N1)N)2 Now if JN = (2N then we have, g9 A1 (ea(2nVi) (1a)(2NVi)) (k 1)A2 (ea2(2yV2) (1a)(21V2)) aN 2k(N 1) 2k(N 1) and if (N1)N = N2 then, ag A1 9 A ,l (eci*(w vi) e(1a1)(E(jv1)Nvv) 8(92 2k(N 1) (k ( 1)A2 (ea2( (N i)NV2) (102)((Ni)N_V2) 2k(N 1) Now recall that for i = 2,..., N 1 we have, ag a Fo((il)N, iN)d2 + Fo(iN, (i+l)N)dX2}. aiN iN i_ Ji So if (ii)N iN then for i= 2,..., N 1 we have, a{ f Fo(4(i1)N,$ iN)dx} = A, e i((I)NV i) ei (g Vi) + a1(iN (i1)N)ea(NVi) kac(N 1) (S (_i1)N)2 1 e(la)(( _)NVi) _ e(1ai)(aENV) _ (1 _ aQ)(SiN (i1)N)e(la2 )(lNV2) k(1 I)2( 1) iN (i1)N )2 (k 1)A2 ea2((_).V2) ea2(NV2) + 2(iN (i1)N)e'2(eNV2) +kac(N 1) (6N (i1)N)2 (k 1)A2 e(1ag)()NV2) e(1a)(N (1 V2) 2) (2iN (i1)N)e(1a2)(CNV2) k(1 2)2 (N 1) (N (1)N)2 and, if 'JN $(iS+i)N then for i = 2,..., N 1 we have, { Fo (iN, (i+1)N)dx2 } A1 eai((.+i)N Vi) eai(~t Vi) a ((i+i)N iN)eaL(EN Vi) ka(N ) ((i+l)N iN) 1 e(1a)((+i)N V) e(1i )(_ Vi) + a1 (i+)N iN)e(1ai)( Vi) k(1 i)2(N 1) ((+l)N iN )2 (k 1)A2 ea2(f(c+l,)v2) ea2(EN V2) 02((i+l)N $iN)ea2('N V2) ka(N 1) ((i)+)N Si)2 (k 1)A2 e(1a2)((+1)NV2) e(1a2)(NV2) + (1 Q)((i+l)N iN)e(1a2)( V2) k(1 c2)2(N 1) ($(i+l)N iN)2 Finally, if (i1)N = tiN then for i = 2,..., N 1 we have a Fo((i1)N, AiN)dx2 1 (e XNi. 2k(N 1) (k 1)A2 (ei2(@i1)NV2) 2k(N 1) and if AiN = k(i+1)N then for i = 2,..., N 1 we have, o iN FO(iN, (i+l)N)dx2} (9tiN 'i  A, (eai(g(i+i)NV1) 2k(N 1) (e + (k 1) A2 (ea2((i+1)N2) 2k(N 1) Once we are able to determine the energy and the gradient of the energy we can implement the optimization strategy developed in [8]. Coding was done in FORTRAN and the results of our implementation are discussed in the next section. 3.3 Numerical Results Here we will both test the accuracy of our asymptotic expansion and observe the behavior of the current by performing numerical experiments in two dimen sions. Note that for the twodimensional problem the domain Q is a unit square and the boundary F is the right side of the unit square, that is r = {(Xl,X2) E : Xi = 1}, (see Figure 23). To compute solutions u uo, and u0 we use piecewise linear finite elements on a regular mesh. To avoid singularities within elements, we chose a grid which conforms to the medium. To perform the nonlinear minimization e (1 a2)(C(~frt)N 1/2)). (when solving for u,), we use a conjugate gradient descent based algorithm devel oped by Hager and Zhang, [8]. Note that the homogenized solution uo is simply a constant value here, which we can find by Newton's Method. The correction, u(1', we compute using standard finite elements for a linear problem, again conforming to the media. We perform these computations for c = 1/5, e = 1/11, e = 1/25 and c = 1/40. We use the following parameter values for our simulation: JA = 1, Jc = 10, VA = 0.5, Vc = 1.0, a, = 0.5, a, = 0.85, and Y = YA U Yc where YA = [0, 1/3] and Yc = [1/3, 1]. Note that for the parameter values used in this implementation, we have u0 = 0.9758. We have analytically shown that the estimates below hold for the case of layered media and wish to numerically verify these estimates: ilUE Uo EU 1) ll(n) < C1E  UO H'o (Q) < C2v The results are summarized in Table 31. The estimates above are all bounded by a term of the form CCe. We estimate this exponent a in the table below. Note that the numerical results in Table 31 are in compliance with the given estimates. Table 31: Table of estimates over Q and convergence rates E 1/5 1/11 1/25 1/40 ce Iur (uo + eu'))lHi(n) .0189 .0090 .0040 .0025 .9699 .9843 .9913 IU, uollHH(n) .0537 .0360 .0238 .0188 .5057 .5061 .5060 IUe uo LL2(n .0063 .0027 .0011 .0007 1.0808 1.0722 1.0676 Table 32: Table of estimates over F and estimates of the gradient over F E 1/5 1/11 1/25 1/40 Iu (uo + i())lL2(r) 0.0108 0.0050 0.0022 0.0014 Iu Uol\L2(r) 0.0128 0.0057 0.0025 0.0015 IIVU V(uo + ,)) L2(r) 0.1027 0.0710 0.0475 0.0377 IIVu, VuoL2(r) 0.1235 0.0817 0.0536 0.0422 In Figure 32 and Figure 33 we plot the "correct" and asymptotic approxi mation of the potential on Q when e = 1/5. We see that the macroscopic behavior is captured by the expansion. Figure 34 and Figure 35 show the same for e = 1/11. In Figures 31(a)31(d) we can view the limiting behavior of u, on F as E approaches 0. To examine the influence of the corrector term more closely, in Fig ures 3639 we graph both the "correct" solution and the asymptotic expansion over P with material regions indicated. Note that the asymptotic approximation is not exact and in fact is slightly skewed. This is probably due to the linearization of the corrector term. In Figure 310 we graph the L"norm of Vu, on the boundary for various values of e. We see that according to our simulations of the layered me dia case, the current remains bounded as the perimeter becomes arbitrarily large, suggesting that the linear relation between current and perimeter observed in [12] may not hold for all geometries. Our results, however, do not directly contradict the observations made in [12], where the computations were done for a fixed num ber of anodes with a varying geometry. Furthermore, since the estimates here are merely in H1(Q), pointwise estimates for the gradient (current) on the boundary do not follow. 0.99  0.98 C S0.97 . 0.96 5 0 0.2 0.4 0.6 0.8 1 x2 1 0.99 0.98 0.97 0.96 0.95 0.94 0.93 0 ..... uO  u s 0.2 0.4 0.6 0.8 1 o 0.97 CL 01 S0.96 0.95 0.94 0.93  0 1 0.99 0.98 o 0.97 0. S0.96 0.95 0.94 0.93 0 .... UO SU 0.2 0.4 0.6 0.8 1 U 0.2 0.4 0.6 0 u.8 1 0.2 0.4 0.6 0.8 1 Figure 31: Limiting behaviour of u, on r e = 1/11, (c) e = 1/25, (d) e = 1/40 as c approaches zero for: (a) e = 1/5, (b)  1 0.99 0.98 0.97 0.96 0.95 0.94 0.93 0.8 .. . 0.6 ... ... " 0.4 0.6 0.4 0.2 0.2 x2 0 0 x2 xl Figure 32: u,, = 1/5 1.02 1. 0.96 ... 0.94 0.92 1 1 S0 0 x1 X2 _ Figure 33: Uo + Eu ', = 1/5 0.9856 0.955,,. , 0.975 . 0 .4 .. .." :.:. .B 0 W,8 097 S039, 0 97. 2 oe Figure 34:0, u = 1/11 0.99 0.95 ....... Figure 35: Uo + cu), e = 1/11 1.02 U 1.01 .. u UO+EUr / U+ E u(I1) 1 Region A / Region C I 0.99 / 50uo 0.97578 ...  .. . ....      0.97 t \t t \', \ / / 0.96 0.95 0.94 0.93 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 X2 Figure 36: The potential on the boundary F, e = 1/5 U 0.995 0.99 Region A 0.985 Region C 0.98 0.975 0.97 0.965 0.96 0.955 I 0.95 0 0.1 0.2 0.3 0.4 0.5 x2 0.6 0.7 0.8 0.9 1 Figure 37: The potential on the boundary F, e = 1/11 S 098 I I I\ 1 ' 0.97 0.965 0.96 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 2 Figure 38: The potential on the boundary F, E = 1/25 us U SRegion A Region C ** n J (, 1 1' 0.1 0.2 0.4 0.5 0.6 0.7 0.8 0.9 x2 Figure 39: The potential on the boundary F, e 0.984 0.982 0.98 H 0.978 , 0.976 > 0.974 0.9721 III. 0 97d 0 968 0.966  0 I j M 1.I.I m . m  I r, i i i r i 1 1 r l n l l   r  I I I I. . 1/40 0.8 0.78 0.76 0.74 e=1/5 e=1/11 0 "o 0 0.72 E=1/40 S0.7 o 0.68 0.66 0.64 0.62 0.6 I 0 0.05 0.1 0.15 0.2 0.25 Figure 310: L" norm of Vu, on P as E approaches 0. CHAPTER 4 A CORRECTOR BASED ON ROBIN BOUNDARY DATA 4.1 A Robin Boundary Condition Note that from Figures 36 through 39 we see that the corrector developed in the Chapter 3 is slightly shifted away from the original. While the shape of this approximation is good, the overall location of the approximation is poor. In this section we attempt to improve the model used to determine the corrector term u'). Thus instead of using (3.1), suppose the correction u') satisfies the Robin boundary condition problem Au' = 0 in 0 _() 1 x ()f  (f( ,uo) fo(uo)) + v (x/e,uo) on r (4.1) 8n E v 0 on 8Q \r 8n Before we develop rigorous estimates or provide numerical data to verify that this approximation will be more accurate let us first intuitively motivate our reason for proposing (4.1). The original model utilized in Chapter 3 to determine u(1) requires only Neumann boundary data. By adding Dirichlet boundary data to the Neumann boundary data we hope that the resulting approximation will have good shape as well as location. Note that it is not a priori obvious that the term u$'f(x/e, o) should be added to the boundary condition. Assuming u, u0o + cu'1 implies u( ) (u, Uo)/e which motivates the boundary condition (f(x/e, u) fo(uo))/c. Using f(x/e, u,) f(x/e, uo) motivates the Neumann boundary conditioned used in Chapter 3. Here using the Taylor Series expansion of f(y,v) in the variable v about uo, i.e. using af f(x/e,u,) f(x/e,uo) + a(x/e,Uo)(U, uo) with Eu') z (u, u,) then yields the Robin boundary condition (4.1). Proposition 4.1.1. Let n = 2 and let u,, Uo be minimizers of (2.7),(2.13) respectively, and let uit) be the solution to (4.1). Then there exists a constant C independent of e such that IIUC Uo CM) II(D) C2. \\uUoeu4\\ H'i(Q) < CE2. Proof. Let ZC = U UO CUM since u, is continuous, by (2.12), we have that for any v E H1(Q), / Vz Vv dx = '. Vv dx Vuo Vv dx e Vuz Vv dx j f(x,u)yvdc, + f( ,uo)vdca, + af So, Vze Vv dx + [f ( ,ue) f (,Uo)]vdox, c (x/c,Uo)U'vd = 0. Now note that uo and u, are defined pointwise on P. So, by Taylors Theorem, using Lagrange's form of the remainder term we have that for each fixed E and x E r there exists (f between uo(x) and u,(x) such that, af x I a2f (x f(X, U') f(X, uo) = f(l,o)(U, Uo) + lJ2f (,9 )( uo)2. E E 9v e 2V2 v2 E By subtracting and adding eiu4 within the parentheses of the first term on the right hand side we have, (X,) x fo(x 8 2(fx o 1 82f xo f(X,ue) f(,uo)= ( ,Uo)z, + ( ,UO)U1 + (, )(U UO)2 c 6 (9v E av E E 2 (9V2 E Thus making the above substitution yields z f^ ^af 1 / 2 jVz Vv dx j+ (x/e,uo)zevda= = j ({(x/, )(U )2vdaU . Now if we pick v = z, this yields /I f 1x af : d g 1 X2f .x SVzdx + d7;u U (2 dux I + v C 2 8v2 E Since tf > co, using a variant of Poincar6 yields Co <1() j< IVZ,2dx+ j ,uo)zd (a 02f (X, (U 2.xo)2 ,dor = 2 j Ea C Now note that 2 y, ) = A(y) [(y)2e(y)(vV(y)) )2e(a(y))(vV()) < A(y)[a(y)2e"(y)(v"(y)) + (1 a(y))2e(1(y))(vV(Y))] af < (y,v) 8v where the last inequality follows from the fact that 0 < a(y) < 1 for all y E Y. So collzE(ll ) < 2 ) I U uO)21z da. (4.2) Now note that when n = 2 IiU Uo11 L(r) < C l1 U 1o 11(,) < C2 1'E Uo1H/2 ) (4.3) where the first inequality follows from the Sobolev Imbedding Theorem [2] and the second inequality follows from the Trace Theorem. Now from standard elliptic regularity theory [10] we have lue Uo11x3/2( ) <_C a (u UO) L2(r). (4.4) By the Mean Value Theorem, for each fixed e and x E P there exists rvf between uo(x) and u,(x) such that, 8f x f(X, UE) f (, uo) = (E UO) (f, Now recall for any v we have, /y(f(y, v) fo(v))dy = 0 so there exists a continuous Yperiodic function g(y, v) such that (y, v) = f(y, v) fo(v), Vv E R. aR f(x/,, ) f(x/E, o) + f(x/,,uo) fo(uo) (f Xg (U o) (x ,)+ Ea (x/E, Uo). This implies f(x/e, U) fo(Uo) an (u, Uno) L =  f(x/, U') fo(UO) L2(r) an fx 8g = (u no) av + 6 ax (x/Ce,uo) L2(r) af x g < \\(ue Uo) ( ,x)\ L2(r) + (x/e, uo) L2(r) 8g < Me6IIUE UIL2(r) + 'll (x/6co) IL2(r) where Me is defined by (3.4). Note that if u' is a weak solution to (3.1) then IIE U0 L2(r) = IIUU CU + UI2(r) < iU UO 6UIL2(r) + EIU L2(r) < C(M, + D) + EC (4.5) where the last inequality follows from Proposition 3.1.1. So in the twodimensional case we have ilue 01L2(,r) < C (4.6) for some constant C independent of e. Thus in the twodimensional case (4.3), (4.4), (4.5) and (4.6) imply that IUe uoILs(r) < Ce (4.7) for some constant C independent of e. Then clearly (uE uo)2 < IU, U I 0 o() < (CE)2. So applying the above estimate to (4.2) yields, Co 2 ((Ce)2 x , C0 o Z, ) < J2 11lP) L1() < c8, lz6L2 P) < C0E2 11 !^:)I where the last inequality follows from the Trace Theorem. Thus, we can write, Izi where C is independent of c and so the proposition is proved. [ 4.2 Numerical Results We use the same parameter values as before and utilize the same finite element based approach as outline in Chapter 3 to discretize the model. The linear problem (4.1) can be solved directly in MATLAB. We provide graphs of the approximation as a threedimensional function in Figures 42 and 44. In Figures 45, 47, 410 and 412 we graph the new approximation on the active boundary. In Figures 46, and 48 we graph both the approximation and the 47 "correct" solution on the boundary. As we see from the graphs and Table 41 this approximation is much more accurate than the approximation used in Chapter 3. Using a corrector based on Robin boundary data yields a substantially more accurate approximation without becoming numerically cumbersome. The new corrector is a substantial improvement over the corrector used in Chapter 3. Table 41: Table of estimates 6 1/2 1/3 1/4 1/5 1/6 Ilut (uo + eut llIoo ) .0019 .0013 .0009 .0007 .0004 IIu (no + U7 )1ll2 r) .5. 29e4 .2908e4 .1i91e4 .1212e4 .0736e4 IIu, (uo + Eu1 )IIL2 ( .5254e4 .2427e4 .1290e4 .0720e4 .0417e4 lu 0 (uo + eu1))\Hi(q) .8406e4 .4420e4 .2877e4 2124e4 ,2104e4 u ,e=15 0.99 0.98 0.97 0.96 0.95 0.94 . 0.93 ..... .. .. .. . 0.6 0.2 0.8 0.9. 0.2 x2 0 0 Figure 41: u, 6 = 1/5 u e (1), =1/5 0.99 : 0.98 0.97 * 0.96 0.95. : '" " " 0.94. 0.93 . 0. 6 8" ... . .. 1 0.2 0 1 x2 0 0 0 Figure 42: uo + cu('), E = 1/5 49 u ,e.1/11 0.99 0.985 20 . 0. 975 0.965 0.96 0.955 x2 0 0 x Figure 43: ua, c 1/11 uo+e u(1), Coll1 0 .98 .. : ... . 0.985 0.97  0.965 0.96. 0.4 0:.6 0.2 2 0.4 X2 0 0x Figure 44: uo + eu) e = 1/11 1  0.99 0.98 0.97 0.96 0.95 0.94 0.93 0 1 0.99 0.98 0.97 0.96 0.95 0.94 0.93 0 0.2 0.4 0.6 0.8 1 Figure 45: Graph of u, (above) and u0 + u1) (below) on the boundary F, E = 1/5 I I I I I I 1 I I I I I I I I I I I 44444 4" *4 e c 4 4 * 44 4 4 * 4 4 4 4 4 44 S I I "BI 0 0.1 02 03 0.4 05 B. 0.7 08 X2 Figure 46: The approximation and the original, E = 1/5 os09 094 0. 1 I I I I I I I I I 0.99 0.985 0.98 0.975 0.97 0.965 0.96 0.955 0 0.99 0.985 0.98 0.975 0.97 0.965 0.96 0.955 0 0.2 0.4 0.6 0.8 1 Figure 47: Graph of u, (above) and uo + u'l) (below) on the boundary F, e = 1/11 0.2 0.4 0.6 0.8 1 448b 4 4 4 4 4 4 4 44 44 44 Vo 44 4 4 4 4 4 4 4 4 4 4 V V 0.1 0.2 0.3 0.4 0.5 X2 0.8 0.7 0.8 09 Figure 48: The approximation and the original, c = 1/11 0 OP 0.s95 0.975 0.97 e 4 4 4 4g 44 4 4 4 f 4 4 4 4 44 4 4 4 4 a 4 4r 4 4 4 0965s . 44p 44 44 Vg I Ilo 098 0 972 0966 0904  0 01 02 03 04 05 06 07 08 09 X2 Figure 49: Graph of u, on the boundary F, e = 1/25 0984 SRegion C 098 0978 0976 09 74 0972 097 0968 0966 0 01 02 03 04 05 06 07 08 09 X2 Figure 410: Graph of Uo+ u1) on the boundary r, ( = 1/25 0982 098 0978 0976 10974 0972 097 01 02 03 04 05 X2 06 07 08 09 Figure 411: Graph of u, on the boundary F, e = 1/40 0982 H Region A 09 Region C 0978 0976 0974 0972 097 01 02 03 04 05 X2 06 07 08 09 Figure 412: Graph of Uo + u1) on the boundary F, c = 1/40 CHAPTER 5 SHIFTING MATERIAL BOUNDARIES 5.1 The Electrostatic Conductivity Model We wish to utilize asymptotic expansions to approximate the solution of a linear elliptic boundary value problem over a twodimensional domain with shifting material boundaries. This is a problem that pertains to electrostatic conductivity and has applications to photonic bandgap (PBG) optical materials. We consider a PDE which models the steady state voltage potential of a conductor with a small inhomogeneity in which there is a discontinuity in the conductivity across the boundary of the inhomogeneity. The boundary of the inhomogeneity shifts by some small amount h (Figure 51). The shift results in a new steady state voltage potential for the conductor. Let uo be the solution to the boundary value problem V ((oVuo) = 0 on Q Buo an '0 o = g on 0 (5.1) where 2 C R2, and g E L2(Q) and ao(x)= { a(x), if x ED oa2(x), if x G \D where D is the region of the small inhomogeneity and a~, a2 > 0. Let uh be the solution to the problem with shifted boundaries, that is uh is a solution to V (ahVuh) = 0 on Q tuh ah = g on 90. (5.2) Here, Here, ah(x) l(x), ifx c Dh a2(x), ifx cf Q\Dh where Dh is the region of the small inhomogeneity but with shifted boundary. Here we assume D C Dh. The steady state voltage potential of the conductor with shifted inhomogeneity is viewed as a perturbation of the original steady state voltage potential. We wish to establish an estimate and do formal asymptotics. X, r,, ,D 1X Figure 51: Perturbation due to shifting between two dielectrics al and a2. 5.2 Estimating the H'I() norm of uo uh We conclude this chapter with an estimate. We use energy methods to rigorously develop an estimate characterizing the limiting behaviour of uh. First we must establish a lemma. Lemma 5.2.1. Let C = 1/ min{il, a2} then for e > C/4 we have 4Ce2(eg a1)2  V( Uh) L2() < e ) VU IIL2(Dh\D)* Proof. For v E H'l(), let Eh(v) be the energy defined by Eh(v) = f I Vv 2 dx gv dx, 2 C2 Jan then Eh(uo) = Vul 2 dx j 9 dx. L Jn Jac Thus we have /sh V(oUo Uh)2 dx = h VUo 02 dx 2 j h (VuLh Vuo)dxI + ) jh VUh2 dx = h  Vuo 2 dx 2 9uo dx+ f hVUh2 d J Jan Jn = 2Eh(O)+ gIh Vuhu2 dx (5.3) where we used the fact that the variational form of (5.2) implies 4 ah(VUh. Vo)dx = I guo dx. Now note that the variational form of (5.1) implies that f gud da = o f VUa 12 dx and the definition of Ch implies / s"hlVuol2 dx =\D U VoU2 dx + fD O Vuo2 dx + fh\D aO Vuo 2 dx J JO\Dh JD JDh\D and similarly the definition of ao implies / ao VIuo 2 dx= 2 Vuo 2 dx + Vu 2 dx + \ 2 V 2 dx. (5.4) JQ Jn\Dh JD JDh\D S ohlVUol2 dx )O 4 Vuol2 dx ( \h 2Vaol2 dx + fD (Uo2 dx+ f\D aUo2 dx) (f\Dh IVUo2 d + DVo2 dx + Dh\D 2IVU 2 dx) 2 2\1VUo12 dx 1 7lVuol2 dx + aC2) Vu0o2 dx. SOJn\Df z D hDh\D z Eh (Uo) Equation (5.4) then implies Eh(,u) = o o Vo 2 dx + 2 Vo0 2 dx + 2) VU0a12 dx Souo2dx+ h\D aia uodx. I romVUO 2 d + 2) dx. 2 Jn f Dh\D 2 So 2Eh(uo) + h VUh 2 dx = o7h VUh 2 dx ao Vuo0 2 dx nlo Jla Jn + fD\ (i  )uo 2 dx. JDh\D Now the variational form of (5.1) and (5.2) imply C/ao VU02dx = ah(Vuh Vuo) dx, and / ah Vuh2dx 0o (VUo VUh)dx J O JOi (5.5) so we have / (ah Vuh 2 _o I0 VU2)d = f (2 V oi)VUo V J O JDh\D So by (5.3), (5.5) and (5.6) we have that j ah V(uo uh)2 dx 2Eh(uo) + h VUha 2 dx = (02 Or)Vuo Vuhdx + f (\ J D \D JfDh\D Thus / ohIV(UO h)d2 dx JO JDh\D 2 JDh\D 'Uhdx. (5.6) 01 2) IVU1o2 dx. a1) Vu Vuah dx h\(2 ) V1)0 Uo2 dx JD(\D O1) (VUh Vuo) V80 dx. So, we have = / ((2 (1i) V(Uh Un ) Vt < IDKD4 u.)12dx+j J4c fDh\D 0o dx (\D E ( 6(72 1,)21Vuo02 dx 2 I)21Vuo0 2 dx thus C = 1/ min{ai, T2} implies I V(uo u 2 dx < C h V(uo Uh)2 dx < C V(U uo) dx+ Ce(x2 a 1l)2 VU02 dx 4c D h\D Hence (1 ) V(uo U) 1 dx< C'e(2 ,T)2 Vol2 dx 46 f Dh\D which implies IV(o h)12 dx < (4 Ce(a2 1)2f u 2 dx 2n (4E G C)Dh\D and so V(U4 Uh) C2(2 a1)2 V(o U)Ln) 4 L2(Dh\D) and thus the lemma is proved. E We now establish an estimate for the H1'() norm of uo uh. Proposition 5.2.2. Let C = 1/ min{(a,a2} and suppose e > C/4 then there exists a constant K(T1, (2), independent of h, such that ( U U 12 + IV(Uo ut)12) dx < K(7i1,a2)h. Proof. By a variant of Poincar6 we have So Uh12 dx < Ci( V(uo U,)2 dx +  Uoh 2 d) Ja/ hV(uo Uh) I dx = 1 IV(u, u)l2 dx (5.7) where the last equality follows from the fact that u, = uh on 80. Thus it suffices to show SIV(uo u)2 dx< 2h. Now by Lemma 5.2.1 we have that 4CG2(a:2 01)2 IIV(uo UI)IL2() L < 4CC (2 l u)2 V L2(D\D). (5.8) Note that due to elliptic regularity IVuo is uniformly bounded on Dh \ D thus IVUo IL2(Dh\D) < C!3Dh \ D 1/2 and note that there exists some positive real number a, independent of h, such that IDh \ D < ah for all h as h 0. So it follows that IIVolIL2(Dh\D) Define 2(^4\E2 ,( 42 U2 71)2 K(a,02) = (C' + 1)30 4c' 1 4E C then (5.7), (5.8), and (5.9) imply that / (ou u2 + I V(0 h)12) dx < K(a, 2)h and thus the proposition is proved. O 5.3 Formal Asymptotics We conclude this chapter with some formal asymptotics. We attempt to characterize the new steady state voltage potential by developing an asymptotic expansion in terms of the shift h using integral equation asymptotics and Green's function. Note that the shifting boundary also causes a perturbation in the conductivity of the conductor. The asymptotic expansion must be developed in such a way as to address this discontinuity. Let Jz(x) be the Dirac delta function centered at z, and let N(x, z) be the solution to V, *.go(x)VNN(x, z) = z(x) in 0 8N(x,z) 1 2 =on 8Q. Then uo(z) = uo(x)6r(x) dx = o() (Vx. ao(x)VxN(x, z)) dx = Vauo(x) go(o(x) VN(x,z) dx 2uo 8 dga 2 fan anx = Vxuo(x) *.o(x)VxN(x,z) dx where the last inequality follows from the fact that uo E H'I() and we assume / no d(x = 0. Now using integration by parts yields Jan uo(z) = (V~ .c o(x)Vuo(x))N(x,z) dx+ j a N(x,z) dgU. So (5.1) implies uo(z) = ~~2 auoN(x, z) da JOfan xn = 29g(x)N(x,z) dau. (5.10) Jan Similarly we have Uh(Z) = Uh(x)6z(x) dx = j (x)(Vx .o(x)VxN(x,z)) dx = VxUh(x) ao(x)VxN(x,z) dx 2Uh dar /n 1 2dnx where the last inequality follows by integration by parts. Now assuming / Uh dox = 0 an yields Uh(z) = f (o(x)VxUh(x) VN(x, z) dx = j ao(x)Vuh(x) VxN(x, z) dx + ao(x)VUh(x) VN(x,z) dx. Jn\Dh J Dh Then D C Dh and integration by parts in the first integral implies Uh(z) = oo(x)N(x,z) do, + o ,(x)Vuh(x) VN(x,z) dx J(n\Dh) anx D + f a2(x)VUh(x) VN(x,z) dx. (5.11) J D \D Note that Dh C 2 implies oo(x) N(x,) d N(x,z) d .a(n\Dh) an, fan anx X) a2x) N(x, z) dox faDh an, and using integration by parts on the last two integrals appearing in the right hand side of (5.11) yields Uh(z) = (x)l N(x, z) dao, a2(x) ( 1 N(x,z) dcx, L an n, faDh S i(x) N(x, z) do + 2(x) N(x, z) da. 'D an, 9(Dh\D) an. Now recall that D C Dh and let 3D = F1 and let r2 and P3 be such that r3 C F1 and aDh = (rl \ r3) U2 (Figure 51). Equation (5.10) implies Uh() = uo(Z) \ 2(X) (2 ) N xz) d U(r (\r3)Urh,2 ( )) + 1(x) UhN(x, z) dua, + a2 (x)aUhN(x,z) dat. ri x 13rsU F2 an, Note we have that /1a '2 (7Ua on F2 S8n, }n. and since duh/dn, is continuous across F3 we have a h = (a0 + on F3. K ian,,) K Z on an. Thus Uh(Z) = uo(z) (i(X) 2(x))(ah ) N(x,z) dao + (Ci() U2(X)) alh N(x, z, ) da, and so Uh(z) uO(z) (Ol(X) "2(X)) 2( J }N(zx, z) dc / Uh N(x, z) dux] Let (s) = (x1(s), 2(s)), s be a parametric equation for the curve OD and let q(s) = (x1(s), i2(s)), sel be a parametric equation for the perturbed boundary 8Dh where (i1(S),C2(S)) = (X1(S), 2(S)) + hO(S)v(S) and where 0(s) is a positive real, function and v is the outward pointing normal vector field on D. In particular let f(s) = (x1(s),X(s)), s E [a,b] be a parametriza tion of F3 and let q(s) = (1l(s), 2(s)), s E [a,b] be a parametrization of F2. Then 2 (an N(x, z) deo = (i(s), Z2(s))N((i,22 i), (z, 2))q'(s) ds JIdKa 2 V7 d a, p Uh< and (ah N(x, z) dc = (2X(S), 2(s))N(( 2) (Zl,Z2)) ds. NUh  Now we use a first order Taylor Series expansion of (1, i2) and N((1,2), 2 (1, 2)) about the point (xl, x2) E F3. Recall that the first order Taylor Series expansion of f(xl, x2) about the point (xl, 2) is given by 8 f (,, Z2 ,xf(xl, x2) f/(ix, 2) = f((l,2) + (l ( 1 l) + ()2 X2). To approximate Iq'(s)l we use a first order Taylor series expansion of g(yl,y2) v1+ y about the point (xl, 2) since I'(s)l = g(1,2). Note that for any s G [a,b] we have 2(s) x(s) = h4(s)v(s) and thus Uh(z) = uo(z) + h(2 Ca1)I where I,= [ (s) (V( i(s))N(x l(s), X2(S), )'(S) + (s) (VN(x (s),x2(s),z) V(S)) l'(s) a( )Uh)n,/i)( v')\1 + ) N(xl(s), X2(S),Z) ) ds. \ V(i) + (X2 2) Now if we assume OD and aDh are parametrized by arclength then I'(s)l = 1 and thus Uh(z) = uo(z) + h(2 71)12 where I2,= [b )(v(U) (s))N(xi(s),2),) + Uh () (VN(xi(s), 2(S), Z) V()) +( ) N(xil(s),x2(s),z)(P'(s) ('/u + Ow/))] ds. Note that to complete the formal asymptotics we need to show 8uh Quo as h 0. On on We leave this as the topic of future work. CHAPTER 6 CONCLUSION We have analyzed a ButlerVolmer type model which describes the potential distribution in a system of anodic islands in a coplanar cathodic matrix with a periodic structure. By using a mutliscale approach we have determined the limiting problem for the boundary value problem (2.5) as the period approaches zero. Furthermore, by introducing a linear correction, we have developed an asymptotic expansion which closely estimates the solution of the original boundary value problem. Essentially, we have taken a nonlinear heterogeneous problem and decomposed it, in a sense, into a nonlinear homogeneous problem and a linear heterogeneous problem. Hence the homogenization approach to this problem gives insight into the behaviour of the solution while also providing an efficient computational technique. The corrector term, although inhomogeneous, solves a linear problem, and was therefore not difficult to compute in our experiments. However, in higher dimen sions or for very small scale problems, one may want to homogenize the corrector term itself. This could perhaps be done by solving a cell problem or looking at the tail behaviour, as in Achdou et al. [1] or Allaire and Amar [3]. In this paper we have used the language and terminology of galvanic corrosion but this analysis could also carry over to a more general class of elliptic problems with nonlinear boundary conditions having periodic structure (assuming the appropriate convexity conditions.) Future work must address the continuity and boundedness issues of the threedimensional problem, i.e. the lack of an applicable Orlicz estimate must be resolved. We wish to do threedimensional numerics and we also wish to consider the model for the case A < 0. With respect to the corrector introduced in Chapter 4 we wish to develop a convergence estimate for the threedimensional case. In addition to implementing thrreedimensional numerical simulations we wish to use the multiscale analysis developed in Chapter 2 and Chapter 4 to construct a first order approximation of solutions to other nonlinear PDE. With regards to the electrostatic voltage potential model of Chapter 5, we wish to complete the asymptotic analysis presented there. We also wish to do numerics simulating electrostatic voltage potential. The end goal is to work up to a three dimensional time harmonic Maxwell's equation so that we may model propagation phenomena and apply this research to PBG structures. REFERENCES [1] Y. ACHDOU, O. PIRONNEAU, AND F. VALENTIN, Effective boundary conditions for laminar flows over periodic rough boundaries, J. Comput. Phys., 147, No.l, (1998), pp.187218. [2] R.A. ADAMS AND J.J. FOURNIER, Sobolev Spaces, 2nd edition, Elsevier, Oxford, U.K., 2003. [3] G. ALLAIRE AND M. AMAR, Boundary layer tails in periodic homogenization, ESAIM Control Optim. Cale. Var., Vol.4, (1999), pp.209243. [4] Y.S. BHAT AND S. MOSKOW, Homogenization of a nonlinear elliptic bound ary value problem modeling galvanic currents, Multiscale Model. Simul., 5 (2006), no.l, pp.149169. [5] K. BRYAN AND M. VOGELIUS, Singular solutions to a nonlinear elliptic boundary value problem originating from corrosion modeling, Quart. Appl. Math, Vol.60, No.4, (2002), pp.675694. [6] D.J. CEDIOFENGYA, S. MOSKOW AND M. VOGELIUS, Identification of conductivity imperfections of small diameter by boundary measurements, University of Minnesota, IMA Preprint Series, No.1502, (1997). [7] L.C. EVANS, Partial Differential Equations, American Mathematical Society, Providence, Rhode Island, 1998. [8] W.W. HAGER AND H. ZHANG CG DESCENT, A conjugate gradient method with guaranteed descent, ACM Transactions on Mathematical Software, (2005), accepted. [9] L.S. HOU AND J.C. TURNER, Analysis and finite element approximation of an optimal control problem in electrochemistry with current density controls, Numer. Math., Vol.71, No.3, (1995), pp.289315. [10] D. JERISON AND C.E. KENIG, The Neumann problem on Lipschitz domains, Bull. Amer. Math. Soc., Vol.4, (1981), pp.203207. [11] J.L. LIONS AND E. MAGENES, NonHomogeneous Boundary Value Problems and Applications, Vol.1, SpringerVerlag, New York, 1972. [12] R. MORRIS AND W. SMYRL, Galvanic Interactions on periodically regular heterogeneous surfaces, AIChE Journal, Vol.34, (1988), pp.723732. [13] R. MORRIS AND W. SMYRL, Galvanic Interactions on random heterogeneous surfaces, J. Electrochem. Soc, No.ll, (1989), pp.32373248. [14] S. MOSKOW AND M. VOGELIUS, First order corrections to the homogenized eigenvalues of a periodic composite medium: The case of Neumann boundary conditions, Indiana Journal, accepted. [15] J.S. NEWMAN, Electrochemical Systems, PrenticeHall, Englewood Cliffs, New Jersey, 1973. [16] M.M. RAO AND Z.D. REN, Theory of Orlicz Spaces, Marcel Dekker, New York, 1991. [17] W. RUDIN, Real and Complex Analysis, McGrawHill, New York, 1966. [18] N.S. TRUDINGER, On imbeddings into Orlicz Spaces and some applications, J. Math. Mech., Vol.17, (1967), pp.473483. [19] M. VOGELIUS AND J.M. Xu, A nonlinear elliptic boundary value problem related to corrosion modeling, Quart. Appl. Math., Vol.56, No.3, (1998), pp.479505. BIOGRAPHICAL SKETCH Sujeet Bhat was born in Bangalore, India in 1972. He lived in Malaysia, In donesia and the Philippines from 1973 to 1990. In 1990 he graduated from the International School, Manila. Sujeet obtained a bachelors degree in mathematics from the University of Florida in 1995, and a master's degree in applied mathe matics from the University of Texas at Dallas in May, 1998. He began his study towards a doctoral degree under the supervision of Professor Shari Moskow in 2003 and received his doctorate in mathematics from the University of Florida in May, 2006. He accepted a two year Industrial Postdoctoral Fellowship from the Institute for Mathematics and its Applications (IMA) at the University of Minnesota in April, 2006. 