Asymptotic and Numerical Modelingof Magnetic Field Profiles inSuperconductors with RoughBoundaries and Multi-Component GasTransport in PEM Fuel CellsbyMichael Robert LindstromB. Sc. (Hons.), University of British Columbia, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate Studies(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2010c© Michael Robert Lindstrom 2010AbstractThis thesis is a combination of two research projects in applied mathematics,which use the applied math techniques of numerical and asymptotic analysisto study real-world problems.The first problem is in superconductivity. This section is motivated byrecent experimental results at the Paul Sherrer Institute. Here, we needto determine how the surface roughness of a superconductor influences thepenetration properties of an externally applied magnetic field. We applyasymptotic analysis to study the influences, and then verify the accuracy -even going well-beyond the limits of the asymptotics - by means of compu-tational approximations. Through our analysis, we are able to offer insightsinto the experimental results, and we discover the influence of a few partic-ular surface geometries.The second problem is in gas diffusion. The application for this studyis in fuel cells. We compare two gas diffusion models in a particular fuelcell component, the gas diffusion layer, which allows transport of reactantgases from channels to reaction sites. These two models have very differentformulations and we explore the question of how they differ qualitatively incomputing concentration changes of gas species. We make use of asymptoticanalysis, but also use computational methods to verify the asymptotics andiiAbstractto study the models more deeply. Our work leads us to a deeper under-standing of the two models, both in how they differ and what similaritiesthey share.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction to Mathematical Modeling . . . . . . . . . . . 11.1 Introduction to Numerical and Asymptotic Analysis . . . . . 11.1.1 Numerical Analysis . . . . . . . . . . . . . . . . . . . 21.1.2 Asymptotic Analysis . . . . . . . . . . . . . . . . . . 51.2 Nondimensionalization . . . . . . . . . . . . . . . . . . . . . 82 Superconductivity Modeling . . . . . . . . . . . . . . . . . . . 112.1 Introduction to Superconductivity . . . . . . . . . . . . . . . 112.1.1 Superconducting Properties . . . . . . . . . . . . . . 112.1.2 Relevant Equations . . . . . . . . . . . . . . . . . . . 122.1.3 Physical Dead Layer . . . . . . . . . . . . . . . . . . . 15ivTable of Contents2.2 Asymptotic Analysis . . . . . . . . . . . . . . . . . . . . . . . 182.2.1 Asymptotic Formulation . . . . . . . . . . . . . . . . 182.2.2 General Procedure . . . . . . . . . . . . . . . . . . . . 272.2.3 Geometry One: Surface with Roughness in One Spa-tial Direction and Parallel Applied Magnetic Field . . 312.2.4 Geometry Two: Surface with Roughness in One Spa-tial Direction and Applied Field Not Uniformly Par-allel to Surface . . . . . . . . . . . . . . . . . . . . . . 432.2.5 Geometry Three: Surface with Roughness in Two Spa-tial Directions . . . . . . . . . . . . . . . . . . . . . . 492.3 Finite Difference Program for Geometry One . . . . . . . . . 572.3.1 Numerical Formulation . . . . . . . . . . . . . . . . . 572.3.2 Validation . . . . . . . . . . . . . . . . . . . . . . . . 642.3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . 682.4 Finite Difference Program for General Sinusoidal Surface . . 702.4.1 Numerical Formulation . . . . . . . . . . . . . . . . . 702.4.2 Validation . . . . . . . . . . . . . . . . . . . . . . . . 842.4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . 872.5 Conclusions of Superconductor Modeling . . . . . . . . . . . 902.5.1 More Complicated Field Behaviour . . . . . . . . . . 912.5.2 Orientation of the Roughness Affects the Profile . . . 922.5.3 Effective Dead Layer . . . . . . . . . . . . . . . . . . 923 Fuel Cell Modeling . . . . . . . . . . . . . . . . . . . . . . . . . 943.1 Modeling of Gas Diffusion in Fuel Cells . . . . . . . . . . . . 94vTable of Contents3.1.1 PEM Fuel Cell Overview . . . . . . . . . . . . . . . . 953.1.2 The Model . . . . . . . . . . . . . . . . . . . . . . . . 963.1.3 Standard Operating Conditions of a Fuel Cell . . . . 983.1.4 Diffusion Equations . . . . . . . . . . . . . . . . . . . 1003.2 Asymptotic Formulation . . . . . . . . . . . . . . . . . . . . 1043.2.1 Nondimensionalization . . . . . . . . . . . . . . . . . 1053.3 Asymptotic Analysis . . . . . . . . . . . . . . . . . . . . . . . 1083.3.1 Asymptotic Analysis of Fick Diffusion . . . . . . . . . 1083.3.2 Asymptotic Analysis of Maxwell-Stefan . . . . . . . . 1113.4 Numerical Analysis of Diffusion Models . . . . . . . . . . . . 1123.4.1 Discretizing the Diffusion Equations . . . . . . . . . . 1123.4.2 Verification of Program Results . . . . . . . . . . . . 1133.5 Exploring the Fundamental Differences between the Models . 1143.6 Conclusions of Gas Diffusion Modeling . . . . . . . . . . . . 1163.6.1 Formulations . . . . . . . . . . . . . . . . . . . . . . . 1163.6.2 Quantitative Differences . . . . . . . . . . . . . . . . 1174 Summary and Future Work . . . . . . . . . . . . . . . . . . . 1194.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1194.1.1 Superconductor Research Summary . . . . . . . . . . 1194.1.2 Gas Diffusion Research Summary . . . . . . . . . . . 1204.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 1214.2.1 Future Work for Superconductor Project . . . . . . . 1214.2.2 Future Work for Gas Diffusion Project . . . . . . . . 122Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124viTable of ContentsAppendicesA Eigenvalues of Finite Difference Matrix . . . . . . . . . . . . 127B The General Interface . . . . . . . . . . . . . . . . . . . . . . . 129C Dead Layers and Averages . . . . . . . . . . . . . . . . . . . . 132D Proof of Periodicity of ˜g . . . . . . . . . . . . . . . . . . . . . 133E Proof of Existence of a Null Vector . . . . . . . . . . . . . . 135viiList of Tables2.1 Estimated orders of convergence for (ωx,ωy) = (pi,0). . . . . . 872.2 Estimated orders of convergence for (ωx,ωy) = (pi,pi). . . . . 873.1 Physical constants in our gas diffusion model. . . . . . . . . . 993.2 Nondimensionalized and rescaled parameters and variablesfor Fick diffusion. . . . . . . . . . . . . . . . . . . . . . . . . . 1073.3 Nondimensionalized and rescaled parameters and variablesfor Maxwell-Stefan diffusion. . . . . . . . . . . . . . . . . . . 1083.4 Different modeling predictions for the relative changes in theconcentrations. . . . . . . . . . . . . . . . . . . . . . . . . . . 115viiiList of Figures2.1 The magnetic field is constant up to the interface and thendecays exponentially in magnitude. Left: a visual represen-tation. Right: Plot of field magnitude vs z, with z = 0 as theinterface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.2 Reprinted figure with permission from [6]. Copyright (2010)by the American Physical Society. This figure displays anexperimentally measured magnetic field profile. The a andb represent different magnetic field orientations, with slightlydifferent decay length scales. These scales are believed to bedue to an anisotropy in the YBCO superconducting material.In both cases, there is a lag in the exponential decay. . . . . . 162.3 A sketch of the effective dead layer, in this case δ. . . . . . . 172.4 The overall geometry, PDEs, and boundary conditions we aresolving. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.5 There is no control over the error at x = 0 but there is controlat x = epsilon1. To ensure O(epsilon1n) accuracy everywhere we use x = epsilon1instead of x = 0 as the point where we switch from lscriptleft tolscriptright. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30ixList of Figures2.6 The applied field is parallel to the surface. . . . . . . . . . . . 322.7 A profile of the field magnitude with epsilon1 = 0.05 and ω = 2piin the y−t. Here the perturbation is quite obvious, but as tgets large the perturbation gradually disappears. . . . . . . . 402.8 Peak and Valley Profiles. . . . . . . . . . . . . . . . . . . . . . 412.9 A profile of the field magnitude with epsilon1 = 0.05 and ω = 2pi inthe peak and valley cases, and the average field profile. . . . . 422.10 Left: the effective dead layer at fixed epsilon1 = 0.05. Right: theeffective dead layer at fixed ω = pi. . . . . . . . . . . . . . . . 432.11 The applied field has nonzero perpendicular components withrespect to the surface. . . . . . . . . . . . . . . . . . . . . . . 442.12 Left: a profile of b1 with epsilon1 = 0.05 and ω = 2pi in the peakand valley cases. Right: a profile of b3 with epsilon1 = 0.05, ω = 2piand x = −0.26 fixed. . . . . . . . . . . . . . . . . . . . . . . . 472.13 A profile of |b|avg with epsilon1 = 0.05 and ω = 2pi. . . . . . . . . . . 482.14 Left: the effective dead layer at fixed epsilon1 = 0.05. Right: theeffective dead layer at fixed ω = pi. . . . . . . . . . . . . . . . 482.15 Top left: a profile of b1 from peak and valley. Top right:profile of b2 with (x,y) = (−0.27,−0.27). Bottom left: profileof b3 with (x,y) = (−0.27,−0.50). Bottom right: differencebetween average field profile in perturbed geometry and flatgeometry. All figures computed with epsilon1 = 0.05, and (ωx,ωy) =(2pi,2pi). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 542.16 Left: the effective dead layer at fixed epsilon1 = 0.05 and ωx = ωy..Right: the effective dead layer at fixed (ωx,ωy) = (pi,pi). . . . 55xList of Figures2.17 The effective dead layer at fixed epsilon1 = 0.05 and (ω2x+ω2y)1/2 = 8pi. 552.18 Near t = 0 the grid is square but as t goes farther out thespacing in the t−direction increases. . . . . . . . . . . . . . . 602.19 A plot of τ = f(t) . . . . . . . . . . . . . . . . . . . . . . . . 612.20 A verification that the two-dimensional code has the rightbehaviour in the flat geometry. . . . . . . . . . . . . . . . . . 652.21 Checking second order convergence by observing the errorbehaviour where the exact solution exp(−z) is known. Hereepsilon1 = 0.1 and ω = 2pi. . . . . . . . . . . . . . . . . . . . . . . . 662.22 Checking the convergence order of the asymptotics. Here ω = pi. 682.23 Profiles of the field magnitude at epsilon1 = 0.05 and ω = 2pi. Fig-ure displays decay from a peak and valley, and the averagemagnitude. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 692.24 Profiles of the field magnitude at epsilon1 = 0.2 and ω = 100pi.Figure displays decay from a peak and valley and shows theflat geometry solution. . . . . . . . . . . . . . . . . . . . . . . 702.25 Formulation for three-dimensional code. . . . . . . . . . . . . 712.26 The mesh in the three-dimensional system. Note how ˜g and bmeshes interlock. The circles with crosses indicate the pointsoccur deeper into the page than the circles wit the dots. Ateach circle with a cross, three components are specified. Ateach circle with a dot, the scalar value of ˜g is specified. . . . . 772.27 Visual confirmation that the three-dimensional program cor-rectly handles the flat interface N = 11, and ωx = ωy = pi. . . 85xiList of Figures2.28 Left: resolution of first-order asymptotic term for first ge-ometry with ω = pi. Centre: resolution of first-order asymp-totic term for second geometry with ω = pi. Right: resolu-tion of first-order asymptotic term for third geometry with(ωx,ωy) = (pi,pi). . . . . . . . . . . . . . . . . . . . . . . . . . 862.29 Top left: a profile of b1 from peak and valley. Top right: aprofile of b2 with (x,y) = (−0.27,−0.27). Bottom left: profileof b3 with (x,y) = (−0.27,−0.50). Bottom right: differencebetween average field profile in perturbed geometry and flatgeometry. All figures computed with epsilon1 = 0.05, and (ωx,ωy) =(2pi,2pi). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 892.30 The field profiles for different roughness orientations with epsilon1 =0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 902.31 The asymptotic and numeric mean field profiles. The two arequite close given the large epsilon1 and ω values. . . . . . . . . . . . 913.1 A cross-section of a PEM fuel cell. Hydrogen and Oxygendiffuse primarily in the XY−plane, but there is a diffusion inthe Z−direction as well. . . . . . . . . . . . . . . . . . . . . 963.2 Our one-dimensional model of the gas diffusion layer. . . . . . 1063.3 Verifying second-order convergence with the Fick code. . . . . 1133.4 The concentration profile of Oxygen in the GDL for both gasdiffusion models. . . . . . . . . . . . . . . . . . . . . . . . . . 1163.5 The concentration profile of water vapor in the GDL for bothgas diffusion models. . . . . . . . . . . . . . . . . . . . . . . . 117xiiList of Figures3.6 The concentration profile of Nitrogen in the GDL for bothgas diffusion models. . . . . . . . . . . . . . . . . . . . . . . . 118xiiiAcknowledgementsA huge thanks to my thesis supervisor, Brian Wetton, for all his help. When-ever I had questions, he was willing to give guidance and direction. He reallyknows his stuff, and without his generosity in time and help, I’d likely stillbe stuck at the beginning.I would also like to thank Rob Kiefl for his supervision in the super-conductor work. This modeling work was initially part of a summer job,which eventually became part of my thesis. Also many thanks to Rob forhis insights on the problem, too.Thanks to Jon Chapman, who made a critical observation in our earlystruggles to understand the proper formulation of the superconductor prob-lem.I would also like to thank Michael Ward, Joerg Rottler, and Matt Chop-tuik for their help along the way.xivChapter 1Introduction toMathematical Modeling1.1 Introduction to Numerical and AsymptoticAnalysisThis thesis explores two physical problems. The first pertains to modelingthe effects of surface roughness on superconductors. In our work, we considerthe surface roughness as a small perturbation from a perfectly flat supercon-ductor with well known properties. The effect of such a small perturbationcan be studied by means of asymptotic analysis. Then, to reach perturba-tions that go beyond the region where the asymptotic work is reliable, weuse finite difference discretizations to numerically approximate the solutionwith a computer. The starting point for all of this is in nondimensionalizingthe system.Our second project involves modeling gas diffusion in a porous media.We have two diffusion models to work with, and by nondimensionalizingtheir equations, we are able to use asymptotic methods to compare them.This work is followed up by numerical computations to verify the accuracy11.1. Introduction to Numerical and Asymptotic Analysisand explore the models from various angles.With this basic context, we begin by providing a brief overview of thepremise of computational approximations to the solutions of differentialequations using finite difference discretizations, the principles of asymptoticanalysis, and the method of nondimensionalization.1.1.1 Numerical AnalysisNumerical analysis is a powerful tool to solve real-world problems that be-come too complicated to solve (or even approximate) analytically. Much ofthis thesis involves the approximation of derivatives by finite differencing.A good reference to which our reader might like to refer is reference [1].The idea is that a differential equation (ordinary or partial) can be writ-ten in a discretized form, with a “small” error.Let’s consider the real-valued function f : R −→ R i.e. f is a scalarfunction that takes a single real number as an input and returns a single realnumber. If we assume that f ∈ C2 (f is twice continuously differentiable)then the limit limh→0 f(x+h)−2f(x)+f(x−h)h2 exists everywhere and is in factequal to fprimeprime(x).The above limit is exact. However, if we only know f on a discreteset of points, we could only come up with an approximation to the secondderivative. Let’s assume that we have values for f at [x0,x1,...,xn] withxj = x0 + jh, where h, the spacing between grid points, is a small positivenumber. Denote these values of f by f0,...,fn.An approximation to fprimeprime(xj), denoted by D2hfj, where j ∈ {1,...,n−1}21.1. Introduction to Numerical and Asymptotic AnalysisisD2hfj = fj+1 −2fj +fj−1h2 . (1.1)It turns out this is actually a very good approximation and the erroris bounded by a constant times h2, as long as f is smooth enough. If weassume f ∈ C4 then we can write fj±1 = fj±hfprime(xj)+h22 fprimeprime(xj)±h36 fprimeprimeprime(xj)+h424f(4)(t) with t being some number between xj and xj ±h. If we substitutethis into (1.1) we get:D2hfj = fprimeprime(xj) + h212f(4)(t).Thus, the error is bounded by Ch2 for some constant C. We write D2hfj =fprimeprime(xj)+O(h2). So if h decreases by a factor of two then the (bound on the)error in the approximation should decrease by a factor of 4. As we makeh smaller (tending to zero) then the error in our approximation also getssmaller (and tends to zero).Being able to discretize in this way means that differential equationscan be approximated by systems of equations. We will consider the simpleexample in solving yprimeprime(x) = x with y(0) = y(1) = 1 on the interval [0,1].We consider the points [x0 = 0,x1 = h,...,xn = nh = 1] and we let yjbe an approximation to y(xj). The boundary conditions tell us that y0 =yn = 1. And in the interior (1 ≤ j ≤ n − 1) we can impose yprimeprime(xj) ≈yj+1−2yj+yj−1h2 = xj.This allows us to write the (n+ 1)×(n+ 1) system31.1. Introduction to Numerical and Asymptotic Analysis1 0 0 0 ··· 0 0 0 01h2−2h21h2 0 ··· 0 0 0 00 1h2 −2h2 1h2 ··· 0 0 0 0· · · · ... · · · ·0 0 0 0 ··· 1h2 −2h2 1h2 00 0 0 0 ··· 0 1h2 −2h2 1h20 0 0 0 ··· 0 0 0 1bracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightLhy0y1y2...yn−2yn−1yn=1h2h...(n−2)h(n−1)h1.Solving this matrix system will produce a set of yprimejs that approximatethe exact solution y(xj) to within an error bounded by a constant timesh2. To see this, we note that we have written the equation Lhˆy = ˆf, whereˆy = (y0,...,yn)T and ˆf = (1,h,...,(n − 1)h,1)T is the load vector. Hereand throughout this thesis, the superscript T will represent the transposeoperator.We know the exact solution, y = (y(x0),...,y(xn))T, will satisfy Lhy =ˆf + O(h2) (O(h2) is a vector with each component bounded by a constanttimes h2) and if we subtract the equation Lhˆy = ˆf then Lhδ = O(h2) whereδ = y − ˆy is the error in the approximation. Multiplying by the inversematrix yields δ = L−1h O(h2). As long as the matrix norm ||L−1h ||∞ staysbounded for all h (or equivalently for all n) then the error δ is of order h2.The proof that the matrix norm is bounded is generally not trivial. How-ever, for a symmetric matrix, if all its eigenvalues are bounded away fromzero then its inverse has a bounded norm. The eigenvalues in this symmetric41.1. Introduction to Numerical and Asymptotic Analysismatrix are all bounded away from zero, and in fact as n gets large the eigen-values of smallest magnitude can be approximated by −pi2,−4pi2,−9pi2,....Therefore, ||L−1h ||∞ is bounded. We include a proof for the values of theeigenvalues in appendix A.1.1.2 Asymptotic AnalysisAsymptotic analysis is a means to study the behaviour of solutions to equa-tions in response to important parameters. These parameters arise in anequation after it has been nondimensionalized i.e. rewritten so that thereare as few constants as possible in the equation (see next section).In asymptotic analysis, we generally know the solution of a particularequation, but are interested to know how the solution changes if the equa-tion is perturbed a little away from the case we know how to solve. Goodbackground material can be found in reference [2].Supposing we had an equation for the variable/function u, f(u;epsilon1) = 0with the small parameter epsilon1, we could assume that u can be written as anasymptotic series in epsilon1. We write u = summationtext∞j=0 fj(epsilon1)uj where fj+1 = o(fj) i.e.as epsilon1 → 0, fj+1(epsilon1)fj(epsilon1) → 0.The most obvious example of an asymptotic series would be a Taylorseries, for example u = u0 +epsilon1u1 +epsilon12u2 +.... Asymptotic series are far moregeneral, however. Their terms can be functions, and even if their terms areconstants, in some cases, peculiar expansions such as u = u0 +epsilon1log(1epsilon1)u1 +sinepsilon1u2 +epsilon13/2u3 +... could be necessary.Here, by means of a simple example we will show how asymptotic seriescan be generated. We will find one that converges and one that diverges.51.1. Introduction to Numerical and Asymptotic AnalysisWe consider the error function, which is important in statistics and prob-ability, erf(z) = 2√pi integraltextz0 exp(−t2)dt.For small |z| i.e. |z| lessmuch 1 we can use Taylor series to approximate theerror function. In this case we are perturbing the integrand away from z = 0.erf(z) = 2√piintegraldisplay z0(∞summationdisplayn=0(−1)nt2nn! )dt =2√pi(∞summationdisplayn=0(−1)n z2n+1(2n+ 1)n!).This Taylor series (which is also an asymptotic series) converges for all zas it comes from term-by-term integration of a power series that convergesuniformly on every compact subset ofC(the series for exp(−t2)).This is not, however, a useful way of approximating erf(z) if |z| is large.As erf(z) is an odd function we will consider z greatermuch 1 here (in this case if wewanted a small parameter we could take epsilon1 = 1/z). Given that erf(∞) = 1we can consider a “small perturbation” from z = ∞ and write:erf(z) = 1− 2√piintegraldisplay ∞zexp(−t2)dt.We will write the integrand in a form conducive to integration by parts,(−2texp(−t2))(−12t ).Focusing on the integral, we will integrate by parts (a couple of times):−12t e−t2|∞z −integraldisplay ∞z12t2e−t2bracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipupright(−2texp(−t2))( −14t3 )dt = −12t e−t2|∞z + 14t3e−t2|∞z +integraldisplay ∞z34t4 exp(−t2)dt.The boundary terms are zero when evaluated at t = ∞. We thus furnish61.1. Introduction to Numerical and Asymptotic Analysisan asymptotic expansion for the error function for large z:erf(z) = 1− 2exp(−z2)√pi ( 12z − 18z3 +...)By induction we could go further to write:erf(z) = 1− exp(−z2)√piz (1 +∞summationdisplayn=1(−1)n(1)(3)...(2n−1)2nz2n ).This asymptotic series diverges. The coefficients cn of (z−2)n are givenby cn = (1)(3)...(2n−1)2n and we use the ratio test to compute limn cncn+1 =limn 22n+1 = 0. The radius of convergence is zero.Although the series diverges, it still obeys the properties of being asymp-totic. The ratio of term n + 1 in the series to the term n in absolute valueis 2n+12 z−2 which tends to 0 as z →∞ (or as epsilon1 = 1/z tends to 0).Also, this is an excellent means of approximating the error function forlarge z. If z = 2 then erf(z) ≈ 0.9953 and if we had used the asymptoticexpansion taking only the first two nonzero terms we would have obtained0.9948. By taking 5 terms we arrive at the best asymptotic approximation,0.9954. After 5 terms the error increases.To get the same precision with the Taylor series we would have needed14 terms to obtain the precision found with only 2 terms in the asymptoticseries, but only 15 terms to obtain the precision found with the 5 terms inthe asymptotic series for large z. The Taylor series will converge to the exactvalue, however.Asymptotic analysis is a reliable tool to give insight into the behaviour ofperturbed systems. Not all asymptotic series converge, but even when they71.2. Nondimensionalizationdo not, they can still provide valuable insights into a perturbed system.In the cases where the series diverge, there is an optimal number of termsneeded to obtain a “best” approximation.Finding an asymptotic series for a problem without proving its conver-gence is known as formal asymptotics. In this thesis we are primarily con-cerned with finding only the first two or three nonzero terms in an asymp-totic series so we will not be concerned with proving convergence and allour expansions will be formal. In these formal expansions, our test as towhether we have found a suitable series is whether we can actually computethe terms in the series.1.2 NondimensionalizationNondimensionalization is used in the study of differential equations to reducethe number of constants appearing in the equation and to gain insight intothe quantities that govern the system’s behaviour.The procedure consists of replacing the variables with dimensionlessnumbers (i.e. getting rid of the units). More details can be found in reference[3].As an example, we consider the heat equation withut = σuxx +f(t)for x ∈ [0,L] and t ≥ 0 with u(0,t) = a, and u(L,t) = b and u(x,0) = u0(x).The heat equation models the temperature, u inside of a material asa function of time t and the spatial position x. The coefficient σ here is81.2. Nondimensionalizationthe thermal diffusivity - describing how efficiently the energy can move inthe system - with units of length squared over time. The function f is theexternal heating: here, we could think of it as a uniform heating that variesin time. The units of f are temperature divided by time. Physically, thetemperatures a and b represent a fixed temperature on either end of thematerial, and u0 is an initial temperature distribution.We begin by writing f = ¯fF, where F has no dimensions and where ¯fis a dimensional quantity holding a representative value for f. Later we willselect the size of ¯f to make the equation as simple as possible.We can rewrite the PDE as:1¯fut =σ¯fuxx +F(t).We will similarly rescale the independent variables now with x = ¯xX,and t = ¯tT along with the dependent variable u = ¯uU. This yields:¯u¯f¯tUT =σ¯u¯f¯x2UXX +G(T),where G(T) = F(¯tT).By equating the coefficients of the U terms, we can select ¯t = ¯x2σ . Wecan also make one final selection in letting ¯u = ¯f¯x2σ .The resulting PDE is to solve:UT = UXX +G(T)with U(0,T) = a/¯u, U(L/¯x,T) = b/¯u, and U(X,0) = ˜u0/¯u where ˜u0(X) =91.2. Nondimensionalizationu0(¯xX) and with X ∈ [0,L/¯x].This equation is in nondimensionalized form, but it can also be nice touse the characteristic scales in the problem. We can choose ¯x = L, anddivide the equation by an upper bound on |G|, say, M. Replacing U byU/M givesUT = UXX +H(T)for some H(T) that takes on values between [−1,1], and with the boundaryconditions U(0,T) = A, U(1,T) = B, initial condition U(X,0) = U0, withA = a¯uM, B = b¯uM, U0 = ˜u0/M, x ∈ [0,1], and T ≥ 0.The problem is a lot simpler now. Our spatial variable ranges between0 and 1. We no longer need to deal with the constant σ. The function H isnicely bounded between −1 and 1.Had the problem been more complicated, we could have had some di-mensionless parameters appearing in the PDE itself, which could be used asasymptotic parameters (this will be apparent in the next two chapters).10Chapter 2Superconductivity Modeling2.1 Introduction to Superconductivity2.1.1 Superconducting PropertiesA superconductor is a material that when cooled to a sufficiently low tem-perature (near absolute zero), exhibits a phase transition to a state with zeroelectrical resistance. This means that an electric current can run through asuperconductor without generating any heat, as long as the current does notexceed some critical value. Superconductors fit into two categories: type 1and type 2.Another defining property of a superconductor is the Miessner effectwhereby the material expels magnetic flux provided the external field doesnot exceed some critical value [4]. Experimentally this means that if there isan applied magnetic field outside of a superconductor, it can only penetratea little ways into the substance before becoming negligibly small. This isknown as the Meissner effect. A sketch of this behavior is given in figure2.1.In an experimental setting (such as in Muon Spin Rotation), a beam ofpolarized, low energy muons is sent through a vacuum region and enters112.1. Introduction to SuperconductivityFigure 2.1: The magnetic field is constant up to the interface and thendecays exponentially in magnitude. Left: a visual representation. Right:Plot of field magnitude vs z, with z = 0 as the interface.the sample, where the magnitude of the local magnetic field causes themto precess. Muons are unstable particles with a very short half-life andthe position and distribution of their decay products can be measured bydetectors. Based on the decay properties and an experimentally controlleddepth at which they decay, muons act as a probe to measure the averagemagnetic field at a given depth within a superconductor [5]. They cannotmeasure the exact field magnitude and direction at specific points.2.1.2 Relevant EquationsThe governing equations for these experimental settings are Maxwell’s equa-tions describing electromagnetic phenomena, which hold everywhere in space,and the London Equation (within the superconductor). In the region of the122.1. Introduction to Superconductivityvacuum, the magnetic field (in SI units) B will satisfy:∇•B = 0 (2.1)and∇×B = µ0J = 0 (2.2)where µ0 is the permeability of free space and where J the local currentdensity. In the vacuum, J = 0 as there are no currents. These experimen-tal settings do not require the consideration of the electric field E as it isconstant.Within the superconductor we have the London equation, which states:∇×∇×B = −B/λ2 (2.3)where λ is the London penetration depth [4]. This length scale λ describeshow far an applied field can penetrate into the sample. Typically it is of theorder of 100 nm [6].The London equation was discovered by the London brothers and itwas initially derived in the efforts to describe an experimentally verifiedrelation within superconductors that mne2 ˙J = E (where the dot signifies atime derivative, m is the mass of an electron, n is the number density ofelectrons and e is their charge) [7]. Their work found two equations, but(2.3) is often referred to as the London equation in many books [4].By using the vector identity that ∇ × ∇× = ∇∇ • −triangle where triangle isthe Laplacian operator, along with (2.1), and (2.3) we find that inside the132.1. Introduction to Superconductivitysuperconductor we can work with:triangleB = B/λ2. (2.4)In the superconducting region both (2.1) and (2.4) must be implemented.We don’t combine (2.1) and (2.2) to use the Laplacian in the vacuumregion, where we would rightfully be able to state that triangleB = 0, because itturns out that we lose information that way: all B for which ∇×B = 0 and∇•B = 0 satisfy triangleB = 0, but there are some B for which ∆B = 0 and∇•B = 0 where ∇×B negationslash= 0 (take B = (z,0,0)T for example).Physically, all components of the magnetic field must be continuousacross the vacuum-superconductor inferface.If we have a flat superconductor occupying the region z > 0 with vacuumin z ≤ 0, and an applied magnetic field B0 that is parallel to the supercon-ductor surface i.e. it has no z−component, then the field should be constantin the vacuum and inside the superconductor it is given by B0e−z/λ.This is a well-known result that’s easily derived in noting that for a flatsuperconductor all x− and y−derivatives must vanish (due to translationalinvariance in x and y). Then in the vacuum, (2.2) implies B1 and B2 areconstant, and (2.1) implies B3 is constant. As B3 = 0 at z = −∞ then B3is zero everywhere in the vacuum. As it must vanish at z = ∞, B3 is zeroeverywhere in the superconductor. So it must be zero everywhere.Looking at B1 and B2 they will have the same value on the interfaceas they did at z = −∞ and the PDE (∂zzB1 = B1λ2 ,∂zzB2 = B2λ2 ) withthe Dirichlet conditions at z = ∞ imply that they decay exponentially in142.1. Introduction to Superconductivitymagnitude.2.1.3 Physical Dead LayerUntil recently, it was assumed that in physical experiments the supercon-ductors were completely flat and hence the fields would decay exponentiallyonce within the superconductor. Experiments at the Paul Sherrer Institute(PSI) have found experimental profiles that differ quite significantly fromthe field profile predicted by this London theory with a flat interface [6].See figure 2.2.The research carried out in the first portion of this thesis has been di-rected to modeling surface roughness, what is believed to be a potentialexplanation to the field perturbations. By surface roughness, we mean asurface that has bumps or wiggles instead of being completely flat.The field seems to decay slower than an exponential near the surfaceand there is the notion of a dead layer - a distance over which the fieldmagnitude does not decay, even within the superconductor and after whichthe field decays exponentially. The cartoon function (in dimensionless units)describing this is that|b| =1 if t ≤ δexp(−(t−δ)) if t > δ.(2.5)The sketch of this cartoon function is given in figure 2.3.We remark that the modeling in this chapter is independent of time. Thevariable t will always represent a nondimensionalized depth and not time.152.1. Introduction to SuperconductivityFigure 2.2: Reprinted figure with permission from [6]. Copyright (2010)by the American Physical Society. This figure displays an experimentallymeasured magnetic field profile. The a and b represent different magneticfield orientations, with slightly different decay length scales. These scales arebelieved to be due to an anisotropy in the YBCO superconducting material.In both cases, there is a lag in the exponential decay.162.1. Introduction to SuperconductivityFigure 2.3: A sketch of the effective dead layer, in this case δ.One way of defining a dead layer, δeff, would be to compute the areaunder the curve |b| from t = 0 to t = ∞. In (2.5), the area is 1 +δ.In our asymptotic work to come, after finding the average field magnitudeat depth t, |b|avg, we can use this to find an effective deadlayer:δeff =integraldisplay ∞0|b|avg(t)dt−1. (2.6)As we will see soon, the definition we give in (2.6) observes many of theproperties we would like in a dead layer (in particular becoming zero if theroughness approaches zero).For any surface that is not flat, it is extremely difficult if not impossi-ble to obtain an exact solution to the problem. However, by considering asmall perturbation (small with respect to λ) then it is possible to obtain172.2. Asymptotic Analysisan approximate analytic solution by means of an asymptotic series. Thismeans it is possible to have reasonably accurate analytic results for cer-tain geometries. In the next section we consider the procedure of findingan asymptotic solution. In subsequent parts of this paper, we numericallycompute approximate solutions for particular surface geometries.2.2 Asymptotic Analysis2.2.1 Asymptotic FormulationNondimensionalization of the ProblemIn general, the vacuum-superconductor interface could be parametrized byz = ah(x,y) where|h|≤ 1 is dimensionless and a is a length (the amplitude).This describes a perturbed flat boundary. However, working with such ageneral surface generally will not lead to analytic formulas that give insightinto the problem and by considering a surface given by acos(ωxx)cos(ωyy),we are still working in quite general terms. Indeed, as long as h and itspowers can be Fourier transformed (where we allow the space of distribu-tions), this is still general. When we refer to h, we refer to this Fourierform cos(ωxx)cos(ωyy). See appendix B for a demonstration of computinga general surface.Vacuum is found in the region z ≤ ah(x,y) and superconductor in theregion z > ah(x,y). We will require that the field approach an appliedfield far away from the superconductor, B(z = −∞) = B0 = |B0|ˆθ =|B0|(θ1,θ2,0)T, and that the Meissner effect is observed, B(z = ∞) = 0.182.2. Asymptotic AnalysisWe remark that ˆθ has no z−component as we would otherwise not be ableto obtain a solution in the flat geometry. The asymptotics rely upon beingable to use the flat geometry solution as a base solution, so our applied fieldcannot have a z−component. We will explain this in more detail shortly.Another condition we require is that the magnetic field is continuous acrossthe superconductor-vacuum interface, [B] = 0 where here [•] denotes thejump of •.We now nondimensionalize as in section 1.2. We select the representativeunits of the magnetic field to be |B0| and the representative units of lengthto be λ.We write b = B|B0| and (xprime,yprime,zprime) = 1λ(x,y,z). Then ∂x = 1λ∂xprime andsimilarly for the other spatial coordinates. So we have ∇• = 1λ∇prime•, ∇× =1λ∇prime× and triangle = 1λ2triangleprime where the prime signifies we’re looking in the prime-coordinates.This allows us to writetriangleb = 1λ2triangleprime = 1λ2b,∇•b = 1λ∇prime•b = 0, and∇×b =1λ∇primeb = 0. We define epsilon1 = a/λ, and if we work in these nondimensionalizedprime coordinates and get rid of the primes we have∇•b = 0 and ∇×b = 0, if z ≤ epsilon1h∇•b = 0 and ∆b = b, if z > epsilon1h(2.7)with boundary conditions b(x,y,−∞) = ˆθ, b(x,y,∞) = 0, [b(x,y,epsilon1h(x,y)] =0.See figure 2.4 for a visual summary of the problem we are solving.To explain the no z−component requirement, we suppose there were an192.2. Asymptotic AnalysisFigure 2.4: The overall geometry, PDEs, and boundary conditions we aresolving.applied field (0,0,1)T and we try to solve the flat geometry. Then all x−and y−derivatives vanish (translational invariance) and we only work withz− derivatives.If the field were b = (b1,b2,b3) then by (2.7) in this reduced case with∂x = ∂y = 0 by using the curl and divergence conditions we have that∂zb1 = ∂zb2 = ∂zb3 = 0 in the vacuum. Thus the field is constant there. Itmust be the value (0,0,1)T because of the boundary condition at z = −∞.In the superconductor from the divergence condition we also have ∂zb3 =0 so that b3 is constant. We have b3(z = ∞) = 0 and so b3 = 0 in thesuperconductor. But b3 = 1 in the vacuum and so b is not continuous, andwe don’t have a solution.Physically, given the divergenceless condition, b can be thought of as202.2. Asymptotic Analysisthe flux of a fluid. By mass conservation, the fluid must go somewhere. Ifthere is some nonzero b3 at z = −∞ and there can be no flux in the x or ydirections there must be some b3 at z = ∞ which cannot come about dueto the Meissner effect. The apparent paradox comes about because we areconsidering a superconductor of infinite size. If it had a finite size then thefield could bend around the superconductor.Finding a Basis in Fourier ComponentsA solution can be expressed as a sum (or an integral) of functions withparticular Fourier frequencies in the x− and y− directions. What we dohere is establish some properties of a general term in the sum. If ˆθ has noz−component then there is a base solution for a flat interface given asb(0) =ˆθ if z ≤ 0ˆθe−z if z > 0(2.8)which satisfies the boundary conditions at ±∞. So the far-field boundaryconditions for these extra Fourier components (which would be added tob(0) to yield a solution) would be to vanish at z = ±∞ and to bring aboutcontinuity along the interface (the base solution is not continuous if z ≤ 0,>0 are replaced by z ≤ epsilon1h,> epsilon1h).If a Fourier component of the field in the vacuum took on the formb = f(z)eiαxxeiαyy for a vector-valued f having value 0 at z = −∞ thenthere would be certain restrictions upon f.212.2. Asymptotic AnalysisThe divergence-free condition of (2.7) yields(iαxf1 +iαyf2 +∂zf3)eiαxxeiαyy = 0so we require:iαxf1 +iαyf2 +∂zf3 = 0. (2.9)And the irrotional condition (the fact that there is a vanishing curl) of(2.7) yields(ˆx(iαyf3 −∂zf2) + ˆy(∂zf1 −iαxf3) + ˆz(iαxf2 −iαyf1))eiαxxeiαyy = 0so we can read off 3 more equations:∂zf2 = iαyf3 (2.10a)∂zf1 = iαxf3 (2.10b)iαxf2 = iαyf1 (2.10c)From (2.9), and (2.10a) and (2.10b) we have the linear equations:∂zf = ∂zf1f2f3=0 0 iαx0 0 iαy−iαx −iαy 0bracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightFf1f2f3.We have a first order linear differential equation in matrix form. In222.2. Asymptotic Analysissolving this problem, we first need to find the eigenvalues of F. They arefound by setting the characteristic polynomial det(F −λI) to 0 :−λ(λ2 −(α2x +α2y)) = 0.We read the eigenvalues off as λ = {0,±α} where we have α =radicalBigα2x +α2y.If vλ is the eigenvector with eigenvalue λ then:f = c0v0 +c−αv−αe−αz +cαvαeαzfor cprimes being constants.To ensure the that f goes to 0 at z = −∞ only cα can be nonzero andthus we only seek the eigenvector with this eigenvalue.If vλ is an eigenvector of F with eigenvalue λ then (F −λI)vλ = 0. Sowe find the null space of−α 0 iαx0 −α iαy−iαx −iαy −α.If β is in the null space then row 1 tells us that αβ1 = iαxβ3 and row2 tells us that αβ2 = iαyβ3. Then we can parametrize β with one degree offreedom, β3: iαxiαyαβ3232.2. Asymptotic AnalysisFor ease of notation later, we’ll rename β3 as β1 and conclude that solu-tions on the vacuum side take on the form:b = β1iαxiαyradicalBigα2x +α2ye√α2x+α2yzeiαxxeiαyy.It’s interesting to note that (2.10c) holds even though we didn’t directlysolve it. Indeed we see that if f = (iαx,iαy,radicalBigα2x +α2y)Te√α2x+α2yz then(iαx)f2 = −αxαye√α2x+α2yz = (iαy)f1.Now we consider a solution on the superconducting side g(z)eiαxxeiαyyvanishing at z = ∞. By having zero divergence we arrive at an equationsimilar to (2.9):iαxg1 +iαyg2 +∂zg3 = 0 (2.11)and by imposing the Laplacian part of (2.7) we have−α2xgi −α2ygi +∂2zzgi = gi for i = 1,2,3after cancelling out the exponential terms.Then we find that ∂2zzgi = α2gi where α ≡radicalBig1 +α2x +α2y. And this tellsus thatg = w−αe−αz +wαeαz.242.2. Asymptotic AnalysisFor g to be 0 at ∞ we set wα to 0. To satisfy (2.11) we requireiαxw−α1 +iαyw−α2 −αw−α3 = 0.We can parametrize g with two degrees of freedom, β2 and β3:g = (α0iαxβ2 +0αiαyβ3)e−αz.So solutions on the superconducting side have the form:b = (radicalBig1 +α2x +α2y0iαxβ2 +0radicalBig1 +α2x +α2yiαyβ3)e−√1+α2x+α2yzeiαxxeiαyy.In the asymptotic version of the problem, we will later need to consider252.2. Asymptotic Analysisa jump condition at z = 0 for a function b of form:b = eiαxxeiαyyβ1iαxiαyradicalBigα2x +α2ye√α2x+α2yz if z ≤ 0(β2radicalBig1 +α2x +α2y0iαx+β30radicalBig1 +α2x +α2yiαy)e−√1+α2x+α2yz if z > 0.(2.12)We can find the jump for given αx and αy frequencies:[b] =−iαxradicalBig1 +α2x +α2y 0−iαy 0radicalBig1 +α2x +α2y−radicalBigα2x +α2y iαx iαybracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightMαx,αyβ1β2β3eiαxxeiαyy.(2.13)This is useful for us because it means if we know the jump at z = 0 for apiecewise Fourier component (of the form expressed in (2.12)), then we canquickly solve for the βprimes and determine its exact form. All we really need isthe inverse matrix. It’s worth noting that when (αx,αy) = (0,0) thenMαxαyis not invertible. However, in such a case, as long as there is no z−componentto the jump, we can still find a solution as β2(1,0,0)T+β3(0,1,0)T still spansR2 ×{0}.262.2. Asymptotic AnalysisLater on we will use the notationM−10,0c in which case we mean the vector(0,β2,β3) where the jump condition [b] = (β2,β3,0)T = c is satisfied.2.2.2 General ProcedureIn this section we present a general asymptotic procedure, and illustrate thesolution method up to second order.We will take the interface to have form z = epsilon1h(x,y). We will applyformal asymptotics (as discussed in section 1.1.2) and express b as a regularasymptotic expansion in epsilon1:b = b(0) +epsilon1b(1) +epsilon12b(2) +O(epsilon13) (2.14)We make use of a Taylor expansion to evaluate b along the surface:b(x,y,epsilon1h) = b(x,y,0) +epsilon1h∂zb(x,y,0) + 12epsilon12h2∂2zzb(x,y,0) +O(epsilon13). (2.15)By substituting (2.14) into (2.15) we can approximate the field to secondorder in epsilon1 along the surface:b(x,y,epsilon1h) = (b(0) +epsilon1(b(1) +h∂zb(0)) +epsilon12(b(2) +h∂zb(1) + 12h2∂2zzb(0)))|z=0+O(epsilon13).Our solution for a flat superconductor is given below:272.2. Asymptotic Analysisb(0) =ˆθ if z ≤ 0ˆθe−z if z > 0.By demanding continuity of b along the interface we arrive at:[b(x,y,epsilon1h)] = ([b(0)]+epsilon1[b(1)+h∂zb(0)]+epsilon12[b(2)+h∂zb(1)+12h2∂2zzb(0)])|z=0+O(epsilon13) = 0which gives us a set of equations at each order in epsilon1:O(1) : [b(0)]|z=0 = 0 (2.16a)O(epsilon1) : [b(1)]|z=0 = −h[∂zb(0)]|z=0 (2.16b)O(epsilon12) : [b(2)]|z=0 = −h[∂zb(1)]|z=0 − 12h2[∂2zzb(0)]|z=0 (2.16c)Of course (2.16a) already holds. From (2.16b) we can find the jump inb(1) and from (2.16c) we can find the jump in b(2).The most general case of a jump is something of the form[b(j)]|z=0 =summationdisplayαx,αyc(j)αx,αyeiαxxeiαyyfor cprimes being some constant vectors (which can be found for different per-turbed geometries as in the following sections). But from (2.13) we can282.2. Asymptotic Analysisactually find the form of b(j) :b(j) =summationdisplayαx,αyiαxiαyαpi1M−1αx,αyc(j)αxαyeαzeiαxxeiαyy if z ≤ 0(α0iαxpi2 +0αiαypi3)(M−1αx,αyc(j)αxαy)e−αzeiαxxeiαyy if z > 0(2.17)where for each (αx,αy) we have that α =radicalBigα2x +α2y and α =radicalBig1 +α2x +α2yand where pij = (ej,•) for j = 1,2,3 is the jth coordinate (pi2(3,4,1)T = 4for example).For shorthand we may later write β(j)αx,αy in place of M−1αx,αyc(j)αx,αy andrefer to the vectors (iαx,iαy,α)T,(α,0,iαx)T,(0,α,iαy)T as γ(1)αx,αy, γ(2)αx,αy,and γ(3)αx,αy respectively.As a final comment on the solutions, we replace the conditions z ≤ 0and z > 0 with z ≤ epsilon1h(x,y) and z > epsilon1h(x,y) respectively. This is becausethe partial differential equations that need to be satisfied hold for z ≤ epsilon1hand z > epsilon1h and not for z ≤ 0 and z > 0.This also ensures continuity to within O(epsilon13). The sketch in figure 2.5should explain this a bit better in a one-dimensional sense. The figuredepicts the use of this convention: for x ≤ 0 there is a function lscriptleft (thestraight horizontal line) and for x > 0 there is a function lscriptright (the othercurve) and we assume they are chosen so that at x = epsilon1 there is continuity292.2. Asymptotic Analysisto within O(epsilon1n) for some n > 0.Figure 2.5: There is no control over the error at x = 0 but there is control atx = epsilon1. To ensure O(epsilon1n) accuracy everywhere we use x = epsilon1 instead of x = 0as the point where we switch from lscriptleft to lscriptright.If we chose to take lscriptleft for x ≤ 0 and lscriptright for x > 0 then we wouldn’thave control over the jump at x = 0. But by taking lscriptleft for x ≤ epsilon1 and lscriptrightfor x > epsilon1 there is only one jump and it is zero to within O(epsilon1n).In terms of applying these computations to the understanding of physicalexperiments, the vector magnetic field at coordinates (x,y,z) is not useful(see sections 2.1.1 and 2.1.3). Instead, the average field magnitude at agiven depth t = z −epsilon1h(x,y) beyond the interface is useful, because this iswhat is actually measured. Therefore in later sections, after obtaining theasymptotic solution in terms of (x,y,z) we will need to convert to a fieldmagnitude expressed in terms of (x,y,t) and average over x and y.302.2. Asymptotic AnalysisIn this last portion of the section, we introduce a few final equations andpieces of notation that will be useful later. The inverse of Mjωx,kωy (definedin (2.13)) is given as:Mjωx,kωy−1 = Γj,kjiωj,kωx kiωj,kωy −ωj,k2k2ω2y +ωj,kωj,k −jkωxωy −jiωj,kωx−jkωxωy j2ω2x +ωj,kωj,k −kiωj,kωywhere we define ωj,k =radicalBig1 +j2ω2x +k2ω2y, ωj,k =radicalBigj2ω2x +k2ω2y, and Γj,k =(ωj,kωj,k2 +j2ω2xωj,k +k2ω2yωj,k)−1.In the following subsections, we consider some specific geometries withincreasing complexity. It may seem peculiar that the analysis is presentedin this order - instead of immediately going to the most general (althoughhighly complex) geometry - but it turns out each geometry needs to beanalyzed independently. Due to a lack of uniform convergence, the resultsof the simpler geometries cannot be obtained from the general geometry.This will become apparent very soon.2.2.3 Geometry One: Surface with Roughness in OneSpatial Direction and Parallel Applied Magnetic FieldHere we explore a simple geometry in the asymptotic regime. We consideran applied magnetic field that is parallel to the surface. Due to the simplic-ity, much of the algebra and work can easily be shown. In more complexgeometries the same techniques are applied but the computations are toolengthy to include in full detail and we also need to make heavy use of the312.2. Asymptotic Analysisnotational shorthand.First-Order TermFigure 2.6: The applied field is parallel to the surface.In our non-dimensionalized regime, we consider a surface of form z =epsilon1cos(ωy) = epsilon12(eiωy + e−iωy) with applied field ˆθ = (1,0,0)T. See figure 2.6.In this case we have that (ωx,ωy) = (0,ω).Thenb(0) =(1,0,0)T if z ≤ 0(1,0,0)Te−z if z > 0is the solution at zeroth order in epsilon1.Then (2.16b) tells us:322.2. Asymptotic Analysis[b(1)]|z=0 = −cos(ωy)[∂zb(0)]|z=0 = (cos(ωy),0,0)T =12(1,0,0)Tbracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightc(1)0,−ωe−iωy + 12(1,0,0)Tbracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightc(1)0,ωeiωy.There is a term with e−iωy and we will needβ(1)0,−ω = M−10,−ω1200=012√1+ω20.It turns out that β(1)0,ω = β(1)0,−ω. Then we can arrive at a first-order pertur-bation by using (2.17):332.2. Asymptotic Analysisb(1) =0−iω|ω|bracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightγ(1)0,−ω(0)bracehtipupleftbracehtipdownrightbracehtipdownleftbracehtipuprightpi1β(1)0,−ωe|ω|ze−iωy +0iω|ω|bracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightγ(1)0,ω(0)bracehtipupleftbracehtipdownrightbracehtipdownleftbracehtipuprightpi1β(1)0,ωe|ω|ze−iωy if z ≤ 0(√1 +ω200bracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightγ(2)0,−ω( 12√1 +ω2)bracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightpi2β(1)0,−ω+0√1 +ω2−iωbracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightγ(3)0,−ω(0)bracehtipupleftbracehtipdownrightbracehtipdownleftbracehtipuprightpi3β(1)0,−ω)e−√1+ω2ze−iωy+(√1 +ω200bracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightγ(2)0,ω( 12√1 +ω2)bracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightpi2β(1)0,ω+0√1 +ω2iωbracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightγ(3)0,ω(0)bracehtipupleftbracehtipdownrightbracehtipdownleftbracehtipuprightpi3β(1)0,ω)e−√1+ω2zeiωy if z > 0.This expression simplifies greatly when we use Euler’s identity to get ourfinal answer:b(1) =0 if z ≤ 0cos(ωy)00e−√1+ω2zif z > 0.342.2. Asymptotic AnalysisSecond-Order TermWe apply the same procedure at second order.[b(2)]|z=0 = −cos(ωy)[∂zb(1)]|z=0−12 cos2(ωy)[∂2zzb(0)]|z=0 = (−cos2(ωy)(−radicalbig1 +ω2+12),0,0)T= 14(radicalbig1 +ω2−12)(1,0,0)Te−2iωy+12(radicalbig1 +ω2−12)(1,0,0)T+14(radicalbig1 +ω2−12)(1,0,0)Te2iωywhere we again expressed the jump as a sum of complex exponentials. Ofcourse the middle term has ω = 0.In looking at (2.17) we need the terms M−10,−2ωc(2)0,−2ω = M−10,2ωc(2)0,2ω =14√1+4ω2(√1 +ω2−12)(0,1,0)T,andM−10,0c(2)0,0 = 12√1+4ω2(√1 +ω2−12)(0,1,0)T.After we apply (2.17) we conclude that:b(2) =0 if z ≤ 012(√1 +ω2 − 12)(e−z + cos(2ωy)e−√1+4ω2z)00if z > 0.We can now write the asymptotic expansion for the magnetic field tosecond order (after replacing z ≤ 0,z > 0 with z ≤ epsilon1cos(ωy),z > epsilon1cos(ωy)):352.2. Asymptotic Analysisb =100if z ≤ epsilon1cos(ωy)100e−z +epsilon1cos(ωy)00e−√1+ω2z+epsilon12(12(√1 +ω2 − 12))(cos(2ωy)00e−√1+4ω2z+e−z00) if z > epsilon1cos(ωy)+O(epsilon13). (2.18)Asymptotic Results for First GeometryThe depth past the surface, t = z −epsilon1h(x,y) is of more physical relevancethan the z−coordinate. We therefore expand (2.18) to second order in epsilon1 (sothe zeroth-order term is expanded to second order, the first-order to firstorder, and the second-order to zeroth order) using Taylor expansions underthe substitution that z = t+epsilon1cos(ωy).Here we’ll illustrate the procedure for the zeroth-order term, but thesame idea holds for all the terms.362.2. Asymptotic Analysis100if t+epsilon1cos(ωy) ≤ epsilon1cos(ωy)100e−t−epsilon1cos(ωy) if t+epsilon1cos(ωy) > epsilon1cos(ωy)=100if t ≤ 0100e−t(1−epsilon1cos(ωy) + 12epsilon12 cos2(ωy)) if t > 0=100if t ≤ 0100e−t if t > 0+epsilon10 if t ≤ 0−cos(ωy)00e−t if t > 0+epsilon120 if t ≤ 012cos2(ωy)00e−t if t > 0+O(epsilon13).In the end, after combining and simplifying the terms up to second order372.2. Asymptotic Analysisin epsilon1, we arrive at:b1(y,t) =1 if t ≤ 0e−t if t > 0+epsilon10 if t ≤ 0(e−√1+ω2t−e−t)cos(ωy) if t > 0+epsilon120 if t ≤ 012(√1 +ω2 − 12)(e−t +e−√1+4ω2t cos(2ωy))+(12e−t −√1 +ω2e−√1+ω2t)cos2(ωy) if t > 0+O(epsilon13)(2.19)with b2 and b3 identically 0.We would finally like to take (2.19) and find the average field magnitude(averaged in this case over y as x doesn’t play a role) at a given depth t. Aswe only need to consider one component and it is positive, we already havethe field magnitude: it is just b1.The average field is given here as|b|avg = ω2piintegraldisplay pi/ω−pi/ω|b|(y,t)dy.Recalling that the cosine function when integrated over a period yieldszero and that the average value of cos2 over a period is 12, this immediatelyyields:|b|avg(t) =1 if t ≤ 0e−t +epsilon1212√1 +ω2(e−t −e−√1+ω2t) if t > 0+O(epsilon13). (2.20)382.2. Asymptotic AnalysisWe can now determine the profile of the magnetic field for small epsilon1 from(2.19) and (2.20). As far as what small epsilon1 means, we need to consider thiscarefully.For a fixed ω, the asymptotic results for the field and average field areaccurate to within O(epsilon13) and the accuracy only approaches zero as epsilon1 → 0.As |ω| → ∞, unless epsilon1 decreases, the accuracy will not remain. This isbecause the second-order term goes like epsilon12√1 +ω2. But if |ω|∼ 1epsilon1 then thesecond-order becomes first-order and the asymptotics break down. Thus theasymptotics for epsilon1 ↓ 0 is not uniformly valid for |ω|→∞.However, the breakdown isn’t necessarily so severe. For large |ω|, theFourier components decay very quickly, so the region over which these larger|ω| values would cause problems would be very small.If epsilon1 lessmuch 1 then we generally seek |ω| lessmuch 1epsilon1. From now on we will alwaysassume ω > 0.We now consider the deadlayer. We compute integraltext∞0 |b|avg(t)dt = 1 +epsilon122 (√1 +ω2 −1) so that δeff = epsilon122 (√1 +ω2 −1) by (2.6).If epsilon1 = 0 or if the spatial frequency goes to zero (so that the interfacelooks flat) we should recover the flat solution and the dead layer should bezero. This formula is in agreement.We have selected a value of epsilon1 that seems physically reasonable. In speak-ing with the experimentalists, there is reasonable confidence the amplitudeof the perturbation, epsilon1 is at most 0.05 [8]. We then chose ω with the conditionthat epsilon1ω is small.Our first plot depicts the field magnitude profile predicted by the asymp-totics - showing how the magnitude varies both with y and with t. This is392.2. Asymptotic AnalysisFigure 2.7: A profile of the field magnitude with epsilon1 = 0.05 and ω = 2pi inthe y − t. Here the perturbation is quite obvious, but as t gets large theperturbation gradually disappears.402.2. Asymptotic AnalysisFigure 2.8: Peak and Valley Profiles.seen in figure 2.7. We can see some effects due to the surface roughness.We now consider the field profile for fixed values of y. When we refer toa peak profile, we will refer to a profile of the field magnitude as a functionof t where cos(ωy) takes on a maximum value and when we refer to a valleywe refer to a point where cos(ωy) takes on a minimum. See figure 2.8.Intuitively, the field should decay slower (with respect to t) from a valleythan a peak. This is because in this geometry the field will not decay inthe vacuum and points located at t > 0 in a valley profile are closer to thevacuum region (where the first component is always 1) than points at t > 0in a peak profile.Our next plot, figure 2.9 shows the peak and valley profiles. It alsoincludes the plot of the average field magnitude after averaging over y. The412.2. Asymptotic AnalysisFigure 2.9: A profile of the field magnitude with epsilon1 = 0.05 and ω = 2pi in thepeak and valley cases, and the average field profile.average lies right in the middle of the two extreme profiles.Given our definition of the dead layer δeff = epsilon122 (√1 +ω2−1), we plottedit. We held epsilon1 fixed at 0.05 and varied ω from 0 to 8pi. We also held ω fixedat pi and varied epsilon1 from 0 to 0.15. The effective dead layer is plotted in figure2.10.This may seem like an extreme range of values given that the asymp-totics are valid for small epsilon1 and ω of order one. However, in computing themaximum difference between the asymptotic and numeric solutions (whichwe explain in section 2.3) at the upper end of these ranges, the sup-norm ofthe difference was 0.0027 (for epsilon1 = 0.15 and ω = pi) and 0.0137 (for epsilon1 = 0.05and ω = 8pi), demonstrating that the asymptotics are reasonably accurateeven in these more extreme limits.422.2. Asymptotic AnalysisFigure 2.10: Left: the effective dead layer at fixed epsilon1 = 0.05. Right: theeffective dead layer at fixed ω = pi.2.2.4 Geometry Two: Surface with Roughness in OneSpatial Direction and Applied Field Not UniformlyParallel to SurfaceFirst-Order TermWe now consider a surface of form z = epsilon1cos(ωx) = epsilon12(eiωx + e−iωx) withapplied field ˆθ = (1,0,0)T. We note that in this case the applied field is nolonger parallel to the surface as can be seen in figure 2.11. In this regime(ωx,ωy) = (ω,0).Then our zeroth-order term isb(0) =(1,0,0)T if z ≤ 0(1,0,0)Te−z if z > 0.By (2.16b), [b(1)]z=0 = −cos(ωx)[∂zb(0)]z=0 = 12(1,0,0)Te−iωx+12(1,0,0)Teiωx432.2. Asymptotic AnalysisFigure 2.11: The applied field has nonzero perpendicular components withrespect to the surface.so we computeβ(1)±ω,0 = 12Γ1,0±iω1,0ωω1,1ω1,10.Using these terms as per (2.17), our resulting first order perturbation isb(1) = Γ1,0−ω2ω1,0 cos(ωx)0−ω1,0ω1,0ωsin(ωx)eω1,0z if z ≤ 0ω1,0ω1,02 cos(ωx)0−ω1,0ω1,0ωsin(ωx)e−ω1,0z if z > 0.442.2. Asymptotic AnalysisAlready we see new behavior emerging: the mixing of field components(b3 is no longer zero). We also note the magnetic field in the vacuum isperturbed.Second-Order TermBy carrying out similar computations as for the case of the perturbationz = epsilon1cos(ωy) now that we know b(1) and b(0) we find:b(2) =4ωE2,0,1 cos(2ωx)02ω2,0E2,0,1 sin(2ωx)eω2,0z if z ≤ 02ω2,0E2,0,2 cos(2ωx)0−4ωE2,0,2 sin(2ωx)e−ω2,0z +E0,0,200e−z if z > 0where E2,0,1 = Γ2,0(−12ω2,0ωσ1 − 14ω2,02σ3), E2,0,2 = Γ2,0(14ω2,0ω2,0σ1 −12ω2,0ωσ3), E0,0,2 =12σ1 and where σ1 = −Γ1,0(−ω1,0ω1,03 +ω1,0ω1,0ω2)− 12and σ3 = −Γ1,0(ω1,0ω1,02ω +ω1,02ω1,0ω).Asymptotic Results for Second GeometryWe have expressed b to second order in epsilon1 with the variable t = z−epsilon1cos(ωx).Although we do not write the expression here, it is of the form b = tildewidestb(0)(t)+epsilon1tildewidestb(1)(t) +epsilon12tildewidestb(2)(t) +O(epsilon13).452.2. Asymptotic AnalysisTo find the field magnitude to second order in epsilon1 we first compute|b|2 = (tildewidestb(0)21 +tildewidestb(0)22 +tildewidestb(0)23) +epsilon1(2tildewidestb(0)1tildewidestb(1)1 + 2tildewidestb(0)2tildewidestb(1)2 + 2tildewidestb(0)3tildewidestb(1)3)+epsilon12(2tildewidestb(0)1tildewidestb(2)1 + 2tildewidestb(0)2tildewidestb(2)2 + 2tildewidestb(0)3tildewidestb(2)3 +tildewidestb(1)21 +tildewidestb(1)22 +tildewidestb(1)23) +O(epsilon13).This is in full generality. Many of the terms are zero. We label the zerothorder term as u0, the first order as u1 and the second order as u2 so that|b| =radicalbig|b|2 =radicalbigu0 +epsilon1u1 +epsilon12u2 +O(epsilon13) = √u0(1+epsilon1 u12u0+epsilon12( u22u0− u218u20))+O(epsilon13).(2.21)Given this, we could compute the average field magnitude by averagingover x. The result is given below:|b|avg(t) =1 +epsilon12(ρ1e2ω1,0t −ρ2eω1,1t) if t ≤ 0e−t +epsilon12(ρ3e−t −ρ4e−ω1,0t +ρ5e(1−2ω1,0t)) if t > 0where we define ρ1 = 14Γ21,0ω1,02ω1,02ω2, ρ2 = 12Γ1,0ω1,0ω1,0ω2, ρ3 = (E0,0,2+14), ρ4 =12Γ1,0ω1,0ω1,03, and ρ5 = 14Γ21,0ω1,02ω1,02ω2.The effective dead layer is computed by (2.6) and we obtainδeff = epsilon12(ρ3 − ρ4ω1,0+ ρ52ω1,0 −1).In this second geometry, we plot the first component of the field in a peakand valley profile. At both the peak and valley, the third component hasvalue zero so we need to select another value of x (in this case −0.26). Here,462.2. Asymptotic AnalysisFigure 2.12: Left: a profile of b1 with epsilon1 = 0.05 and ω = 2pi in the peak andvalley cases. Right: a profile of b3 with epsilon1 = 0.05, ω = 2pi and x = −0.26fixed.peaks and valleys refer to slices of constant x which maximize or minimizecos(ωx) respectively. The plot of both the first and third component is givenin figure 2.12. The profile of the average field is given in figure 2.13.We also show the effective dead layer in figure 2.14, over the same rangesof epsilon1 and ω as in the previous geometry. The respective sup-norm errorsverified numerically in section 2.3 for epsilon1 = 0.15 with ω = pi and for epsilon1 = 0.05with ω = 8pi are 0.00462 and 0.00641 respectively.This geometry is a lot more interesting. We see that with respect tothe first component of the field, the decay begins even before the interface(when looking at the peak profile).We also note that the individual components can take on values thatexceed their value at z = −∞. It might be surprising, but there’s nothingthat says this cannot happen. Indeed, we have that triangleb = 0 in the vacuumregion and therefore each component should obey the maximum/minimum472.2. Asymptotic AnalysisFigure 2.13: A profile of |b|avg with epsilon1 = 0.05 and ω = 2pi.Figure 2.14: Left: the effective dead layer at fixed epsilon1 = 0.05. Right: theeffective dead layer at fixed ω = pi.482.2. Asymptotic Analysisprinciple (because it is harmonic). The maximum principle states that if f isa harmonic function in a domain D then its maxima and minima lie on ∂D(the boundary). In our case, the maxima of b1 and b3 lie on the boundaryof the vacuum region - right on the interface.The dead layer is actually negative here. By adding roughness in thedirection of the field, the field is no longer constant on the vacuum side, andthe net effect is to diminish the average field magnitude as compared to theflat geometry (so that δeff < 0). We also note that the magnitude of thedead layer in this geometry is much smaller than in the first geometry.2.2.5 Geometry Three: Surface with Roughness in TwoSpatial DirectionsFirst-Order TermAgain we take the applied field as ˆθ = (1,0,0)T and we have the samebase solution as in the previous geometries. By (2.16b), after decomposingthe surface z = epsilon1cos(ωxx)cos(ωyy) into components z = 14(ei(−ωx−ωy) +ei(−ωx+ωy) +ei(ωx−ωy) +ei(ωx+ωy)) we require[b(1)]|z=0 = 14(ei(−ωx−ωy) +ei(−ωx+ωy) +ei(ωx−ωy) +ei(ωx+ωy))(1,0,0)T.We then compute four terms such terms asM−1ωx,ωy1400and use (2.17)to get the first-order perturbation:492.2. Asymptotic Analysisb(1) =Γ1,1−ω11ω2x cos(ωxx)cos(ωyy)ω1,1ωxωy sin(ωxx)sin(ωyy)−ω1,1ω1,1ωx sin(ωxx)cos(ωyy)eω1,1z if z ≤ 0Γ1,1ω1,1(ω1,1ω1,1 +ω2y)cos(ωxx)cos(ωyy)ω1,1ωxωy sin(ωxx)sin(ωyy)−ω1,1ω1,1ωx sin(ωxx)cos(ωyy)e−ω1,1z if z > 0.Second-Order TermThe second order term can be worked out using the same algorithms as inprevious sections and we find502.2. Asymptotic Analysisb(2) =−8E2,2,1ωx cos(2ωxx)cos(2ωyy)8E2,2,1ωy sin(2ωxx)sin(2ωyy)−4E2,2,1ω2,2 sin(2ωxx)cos(2ωyy)eω2,2z+−4E2,0,1ωx cos(2ωxx)0−2E2,0,1ω2,0 sin(2ωxx)eω2,0z if z ≤ 04E2,2,2ω2,2 cos(2ωxx)cos(2ωyy)−4E2,2,3ω2,2 sin(2ωxx)sin(2ωyy)(−8E2,2,2ωx −8E2,2,3ωy)sin(2ωxx)cos(2ωyy)e−ω2,2z+2E2,0,2ω2,0 cos(2ωxx)0−4E2,0,2ωx sin(2ωxx)e−ω2,0z +2E0,2,2ω0,2 cos(2ωyy)00e−ω0,2z+σ1400e−z if z > 0where we define the E’s by: E2,2,1 = Γ2,2(18ω2,2ωxσ1−18ω2,2ωyσ2+ 116ω2,22σ3),E2,2,2 = Γ2,2( 116(4ω2y+ω2,2ω2,2)σ1+14ωxωyσ2−18ω2,2ωxσ3),E2,2,3 = Γ2,2(−14ωxωyσ1−116(4ω2x + ω2,2ω2,2)σ2 − 18ω2,2ωyσ3), E2,0,1 = Γ2,0(14ω2,0ωxσ1 +18ω2,02σ3),E2,0,2 = Γ2,0(18ω2,0ω2,0σ1−14ω2,0ωxσ3),andE0,2,2 = Γ0,2(18(4ω2y+ω0,2ω0,2)σ1).We also have that σ1 = −12 −J1, σ2 = Γ1,1(ω1,12ωxωy + ω1,1ω1,1ωxωy),σ3 = Γ1,1(−ω1,1ω1,12ωx−ω1,12ω1,1ωx), and J1 = Γ1,1(−ω1,12(ω2y+ω1,1ω1,1)+ω1,1ω1,1ω2x).512.2. Asymptotic AnalysisAsymptotic Results for Third GeometryWe can now expand the solution in powers of epsilon1 with the variable t =z−epsilon1cos(ωxx)cos(ωyy). We don’t show the result here as it is very complex,but we emphasize that our result has been verified thoroughly with Maple.We first checked that before re-expressing in terms of t, the terms on boththe vacuum and superconducting sides satisfied the PDEs. Taking these ex-pressions, we expanded to second order in epsilon1 with z = t+epsilon1cos(ωxx)cos(ωyy).Then we substituted t = 0 into the asymptotic solution and verified that forrandomly chosen x, y, ωx, and ωy the solution was continuous at each orderin epsilon1.Using maple, we can also compute|b|avg(t) = ωxωy4pi2integraldisplay pi/ωy−pi/ωyintegraldisplay pi/ωx−pi/ωx|b|(x,y,t)dxdywhere |b| has been expressed in powers of epsilon1 up to second order (by means of(2.21)). The result is given below:|b|avg(t) =1 +epsilon12(ρ1e2ω1,1t −ρ2eω1,1t) if t ≤ 0e−t +epsilon12(−ρ3e−ω1,1t +ρ4e(1−2ω1,1)t −ρ5e−t) if t > 0(2.22)where we have the definitions ρ1 = 18Γ21,1(ω1,12ω2xω2y + ω1,12ω1,12ω2x), ρ2 =14Γ1,1ω1,1ω1,1ω2x, ρ3 = 14Γ1,1(ω1,12ω2y + ω1,1ω1,13), ρ4 = 18Γ21,1(ω1,12ω1,12ω2x +ω1,12ω2xω2y), and ρ5 = J14 .522.2. Asymptotic AnalysisBy (2.6) we find in this model thatδeff = epsilon12(− ρ3ω1,1+ ρ41−2ω1,1−ρ5).We now compute profiles for the field components and average field withepsilon1 = 0.05 and (ωx,ωy) = (2pi,2pi). In this geometry, peaks and valleys arenot uniquely defined in terms of fixed x and y values. We define peaks asoccurring at x = −pi/ωx and y = −pi/ωy and valleys as occurring at x = 0and y = −pi/ωy. Both the second and third components vanished at peaksand valleys so to obtain any profile, we selected pairs of (x,y) where they didnot vanish. Figure 2.15 displays these results. In this regime, the averagefield still looks very much like the flat solution and we do not include itsplot.Again we see b1 exceeding its value of b1(z = −∞) = 1. The dead layeris another interesting point to investigate. We explore it from three angles,displayed in figures 2.16 and 2.17.Initially, we fix epsilon1 = 0.05 and set ωx = ωy = ω and plot how the deadlayer varies with what we could call the net spatial frequencyradicalBigω2x +ω2y.We also fix (ωx,ωy) = (pi,pi) and plot the dead layer as a function of epsilon1. Asverified numerically, the sup-norm of the errors at the upper end of theseplots were 0.0174 and 0.00662 respectively.Then, for fixedradicalBigω2x +ω2y = 8pi and epsilon1 = 0.05, we explore how the deadlayer varies with the ratio ωxωy. As ωx → 0, we would expect recover the firstgeometry, and as ωy → 0, we would expect to recover the second geometry.The plot depicts qualitatively what we expect. As the ratio gets small, the532.2. Asymptotic AnalysisFigure 2.15: Top left: a profile of b1 from peak and valley. Top right:profile of b2 with (x,y) = (−0.27,−0.27). Bottom left: profile of b3 with(x,y) = (−0.27,−0.50). Bottom right: difference between average fieldprofile in perturbed geometry and flat geometry. All figures computed withepsilon1 = 0.05, and (ωx,ωy) = (2pi,2pi).542.2. Asymptotic AnalysisFigure 2.16: Left: the effective dead layer at fixed epsilon1 = 0.05 and ωx = ωy..Right: the effective dead layer at fixed (ωx,ωy) = (pi,pi).Figure 2.17: The effective dead layer at fixed epsilon1 = 0.05 and (ω2x+ω2y)1/2 = 8pi.552.2. Asymptotic Analysisdead layer increases. As it gets big then the roughness acts more and morein the direction of the field and the dead layer decreases.More interestingly, as the ratio reaches zero, the dead layer of (ωx,ωy) =(8pi,0) is approximately half of the dead layer in the first geometry withω = 8pi. Also, as the ratio reaches its maximum value, the dead layer reachesa value of −2.95×10−4 (this value is off the plot range of the graph), whichis also about half of the dead layer size in the second geometry with (0,8pi).This may seem very surprising, but we show exactly where this factor of12 comes from in appendix C. Ultimately, taking a limit as a spatial frequencyapproaches zero in the third geometry is unphysical experimentally becausephysicists can only average over a finite range (if a frequency goes to zero,the region over which averaging needs to take place grows without bound).In the next section we turn to numerical computations.562.3. Finite Difference Program for Geometry One2.3 Finite Difference Program for Geometry OneHere we are considering an interfacez = epsilon1cos(ωy),with applied field (1,0,0)Twhich we studied asymptotically in section 2.2.3.2.3.1 Numerical FormulationConstant Field on the Vacuum SideWe shall take the field to be constant on the vacuum side. Here we showthat there is a solution where the field is constant on the vacuum side.Uniqueness is not proven here, but given the three-dimensional numericalwork that follows (section 2.4) it seems highly likely.If we takeb =100if z ≤ epsilon1cos(ωy)b(y,z)00if z > epsilon1cos(ωy)with the condition [b]|z=epsilon1cos(ωy) = 0, then as long as triangleb(y,z) = b(y,z),b(y,epsilon1cos(ωy)) = 1 and b(y,z) → 0 as z → ∞ (which we will solve fornumerically) all PDEs and boundary conditions are satisfied.572.3. Finite Difference Program for Geometry OneCoordinate TransformationsTrying to implement a finite difference mesh on a sinusoidal surface is ex-tremely difficult. In addition, the main variable of interest is the depth pastthe sample surface - not the z−coordinate. We can solve both problems bya transformation of coordinates.Here we will define σ = y and t = z −epsilon1cos(ωy). We now need to findthe Laplacian (∂yy +∂zz) in the new coordinates.∂y = ∂yσ∂σ +∂yt∂t = ∂σ +epsilon1ωsin(ωy)∂t = ∂σ +epsilon1ωsin(ωσ)∂tThen∂yy = ∂σ(∂σ +epsilon1ωsin(ωσ)∂t) +epsilon1ωsin(ωσ)∂t(∂σ +epsilon1ωsin(ωσ)∂t)= ∂σσ +epsilon1ω2 cos(ωσ)∂t +epsilon1ωsin(ωσ)∂tσ +epsilon1ωsin(ωσ)∂σt +epsilon12ω2 sin2(ωσ)∂2tt.And∂z = ∂zσ∂σ +∂zt∂t = ∂t.So clearly∂zz = ∂tt.We arrive at the Laplacian in the new coordinates (after assuming allfunctions are C2 and making use of the equality of mixed partial derivatives):triangle = ∂σσ + (1 +epsilon12ω2 sin2(ωσ))∂tt + 2epsilon1ωsin(ωσ)∂σt +epsilon1ω2 cos(ωσ)∂t.582.3. Finite Difference Program for Geometry OneThere’s one more thing that needs to be considered. Given that we knowthe field is zero at t = ∞, it must decay far from the interface. Near theinterface itself is where the most interesting phenomena will occur. In termsof a mesh, we need more points closer to the interface to numerically resolveall the intricacies in the solution and fewer points farther away. This will saveus computation time. It is also helpful as numerical results can be unreliableif two dimensions of a mesh differ significantly. This aspect of minimizinggrid points is of paramount importance later on in our three-dimensionalcode.To achieve these extra grid properties, one more coordinate transforma-tion is needed. We will do this transformation in the next section as wedescribe how we are building our mesh.The GridWe will have a mesh in (σ,τ)−coordinates where τ is a transformation oft. We will enforce periodicity in σ with spatial period 2pi/ω. We need onlyconsider the single nonzero component of b in this analysis. Instead of callingit b1 we can just call it b.On the vacuum side (up to t = 0) we have b = 1. We will focus on thesuperconducting side in our mesh.Numerically we will impose t = ∞ as t = M greatermuch 1.We now pick a natural number N and we discretize in σ. We choosehσ = 2piNω592.3. Finite Difference Program for Geometry OneFigure 2.18: Near t = 0 the grid is square but as t goes farther out thespacing in the t−direction increases.and chooseσi = −piω + (i−1)hσ for i = 1,2,...N.Because we want both grid dimensions to scale with N we will choose asimilar discretization for τ :τj = α+ (j −1)hσ for all j ∈Nfor which τj ≤ β.We call the number of τ points Nτ.We have some freedom over the choice of α and β and will select themlater.Now we consider what we want in terms of a mesh. We would like thatthe t−spacing is roughly equal to σ−spacing near the interface (to resolvedetails), but then far away we would like the t−spacing to scale roughlywith M/N. See figure 2.18.Near the interface if ht denotes the t−spacing then we want ht ≈ hσ.Thus, given that τ is incremented in units of hσ near the interface dtdτ ≈ 1.602.3. Finite Difference Program for Geometry OneFigure 2.19: A plot of τ = f(t)Far from the interface, we want ht ≈ MN = Mω2pi hσ. So dtdτ ≈ Mω2pi . We willdefine η = 2piMω to be the value of dτdt near t = M.The algebra is a little easier if we consider τ = f(t) where f is invertibleand satisfies f(0) = α < β = f(M). We will consider f to have the formalog(t−b). See figure 2.19 for the general shape of f.From the derivative conditions we have a−b = 1 and aM−b = η whichshows us that a = ηM1−η > 0 and b = −ηM1−η < 0.Knowing a and b, we now calculate α = alog(−b) and β = alog(M −b).Under this transformation, ∂t = fprime(t)∂τ = at−b∂τ = ae−τ/a∂τ (this comesfrom observing if τ = alog(t−b) then 1/(t−b) = 1/eτ/a).We also have ∂tt = e−2τ/a(−a∂τ +a2∂ττ).Substituting the τ− derivatives in place of t−derivatives gives:612.3. Finite Difference Program for Geometry One∆ = ∂σσ +a2(1 +epsilon12ω2 sin2(ωσ))e−2τ/abracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightξ(1)∂ττ + 2aepsilon1ωsin(ωσ)e−τ/abracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightξ(2)∂στ+a(epsilon1ω2 cos2(ωσ)e−τ/a −(1 +epsilon12ω2 sin2(ωσ))e−2τ/a)bracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightξ(3)∂τ.So a generic point on the mesh could be described by (σi,τj). It is alsonecessary to map (i,j) coordinates to just a single number, K(i,j) = j+(i−1)Nτ. This way we can turn the whole PDE system into a matrix equationof form Mu = f where M is a matrix that comes from the discretizedequations, u is our unknown, and f is the load vector (see section 1.1.1).The Discretized EquationsConditions along the Interface and at z = ∞ The interface corre-sponds toj = 1.The field being 1 on the interface corresponse toMK(i,1),K(i,1) =1 with fK(i,1) = 1 for all i.The far-field equation is for j = Nτ. The field being 0 at τNτ correspondsto MK(i,Nτ),K(i,Nτ) = 1 and fK(i,Nτ) = 0 for all i.Laplacian Here we consider 2 ≤ j ≤ Nτ −1.In this regime, we have the following second order approximations:∂σσb(σi,τj) = uK(i+1,j) −2uK(i,j) +uK(i−1,j)h2σ+O(h2σ)∂ττb(si,tj) = uK(i,j+1) −2uK(i,j) +uK(i,j−1)h2σ+O(h2σ)622.3. Finite Difference Program for Geometry One∂στb(si,tj) = uK(i+1,j+1) +uK(i−1,j−1) −uK(i+1,j−1) −uK(i−1,j+1)4h2σ+O(h2σ)∂τb(σi,tj) = uK(i,j+1) −uK(i,j−1)2hσ+O(h2σ)for all i, where to impose periodicity, when i = 1 we replace i−1 with Nsand when i = Ns we replace i+ 1 with 1.Inside the superconductor we impose the discretized version of ∆b = bby choosing the appropriate rows for M. Defining discretized ξ variables (sothat ξ(2)i,j = 2aepsilon1ωsin(ωsi)e−τj/a, for example) our matrix M is assembled as:MK(i,j),K(i,j) = −2ξ(1)i,jh2σ −2h2σ −1MK(i,j),K(i±1,j) = 1h2σMK(i,j),K(i,j±1) = ξ(1)i,jh2σ ±ξ(3)i,j2hσMK(i,j),K(i+1,j+1) = ξ(2)i,j4h2σMK(i,j),K(i+1,j−1) = −ξ(2)i,j4h2σMK(i,j),K(i−1,j+1) = −ξ(2)i,j4h2σMK(i,j),K(i−1,j−1) = ξ(2)i,j4h2σfor all i. We also set fK(i,j) = 0.632.3. Finite Difference Program for Geometry OneFinding the Results We define a discretization of t :tj = eτj/a +b for j = 1,...,Nτ.The approximate magnetic field at (y = σi,t = tj) is given by uK(i,j). Wherenecessary, an average field at depth tj can be found by averaging uK(i,j) overi with j fixed.2.3.2 ValidationThe program can handle a larger range values of ω and epsilon1 than the asymp-totics. However, we test the program first to see it is providing reasonableresults.We found that M = 9.25 and N = 50 provided results that were accurateto within 10−4 - an error we feel is very acceptable. Generally these are theparameters we set.Order of Convergence with Respect to Exact SolutionIf epsilon1 = 0 then for any arbitrary ω we have a flat superconductor so the solutionshould decay exponentially past the superconducting surface at exp(−t). Weselected ω = pi with N = 50 and plotted result. The numerical results areshown in figure 2.20.642.3. Finite Difference Program for Geometry OneFigure 2.20: A verification that the two-dimensional code has the rightbehaviour in the flat geometry.We can also choose boundary conditions so that the exact solution isexp(−z) even for nonzero epsilon1. We observe that if b(y,z) = exp(−z) thentriangleb = b and b(y,epsilon1cos(ωy)) = b(σ,t = 0) = exp(−epsilon1cos(ωσ)). That’s what wechoose as our boundary condition.If the error E = ||bnum −bex||∞ is of the order h2σ (which is of the orderN−2) then on a log-log plot, we would expect that logE is a linear functionof logN with slope−2. The subscripts “num” and “ex” denote the numericaland exact solutions respectively. Later on “asy,k” will signify the kth-orderasymptotic solution.The plot verifying second order convergence is given in figure 2.21.From these results, the error at N = 60 (not plotted) corresponded to1.42 × 10−4 which is very near the error in setting b(M) = exp(−9.25) ≈652.3. Finite Difference Program for Geometry OneFigure 2.21: Checking second order convergence by observing the error be-haviour where the exact solution exp(−z) is known. Here epsilon1 = 0.1 andω = 2pi.662.3. Finite Difference Program for Geometry One9.61×10−5 to 0.With the same amplitude and frequency, we set M = 11 and increasedN until the error plateaued. The error was 2.26 × 10−5 and exp(−11) ≈1.67×10−5. These results give us estimates on the errors coming from theapproximation of the far field condition at a finite length from the interface;it roughly scales with exp(−M).Comparison with Asymptotic SolutionBy selecting large N, in this case N = 50 and M = 9.25, the numericalsolution should be very close to the exact solution. We would thereforeexpect the asymptotics when taken to the second order term to convergeat a rate of epsilon13 to the numerical solution. However, this is a very delicatecomputation (and it is even more delicate in the three-dimensional code).We have bnum = bex + O(h2σ) + O(λ) (where λ represents the far-fielderror). Also, we have basy,2 = bex + O(epsilon13). Subtracting the numerical andsecond-order asymptotic solutions yields a difference O(h2σ)+O(λ)+O(epsilon13).In general, the error terms O(h2σ) and O(λ) can depend upon epsilon1. If epsilon1 istoo small then the O(epsilon13) term cannot be detected. Once epsilon1 gets too large,the fourth-order term in the asymptotics could exceed the third-order termand once again we wouldn’t detect O(epsilon13).For this code, by choosing epsilon1 in increments of 0.05 from 0.05 to 0.3, we canverify the order of convergence is O(epsilon13) by computing the sup-norm of thedifference between the asymptotic and numerical solution at the differentvalues of the amplitude. See figure 2.22.672.3. Finite Difference Program for Geometry OneFigure 2.22: Checking the convergence order of the asymptotics. Here ω =pi.2.3.3 ResultsHaving established the validity of the simulation and having some under-standing of the level of discretization needed to obtain a desired accuracy,we are now in a position to examine some physical problems.We start off by picking a case we have also studied in the asymptoticanalysis of the first geometry. We choose epsilon1 = 0.05 and ω = 2pi. A plotdepicting the field magnitude on average and from a peak and valley canbe found in figure 2.26. The results are nearly identical to the asymptoticvalues.We can also select values of epsilon1 and ω that are quite a bit larger. Inparticular, we choose a case with epsilon1 = 0.2 and ω = 100pi. This plot is found in682.3. Finite Difference Program for Geometry OneFigure 2.23: Profiles of the field magnitude at epsilon1 = 0.05 and ω = 2pi. Figuredisplays decay from a peak and valley, and the average magnitude.figure 2.24. Figure 2.24 depicts what we would expect intuitively. Given thevery high spatial frequency, we would expect the field to be nearly constantuntil past the deepest peaks (because it must be continuous and have theconstant value of 1 in the vacuum). As a result, the field decay should bedelayed when looking along a valley profile. In the plot, the field seems tobe almost constant in the valley profile for a scaled depth of nearly 0.4. Eventhe peak profile is significantly larger than the flat profile.692.4. Finite Difference Program for General Sinusoidal SurfaceFigure 2.24: Profiles of the field magnitude at epsilon1 = 0.2 and ω = 100pi. Figuredisplays decay from a peak and valley and shows the flat geometry solution.2.4 Finite Difference Program for GeneralSinusoidal SurfaceHere we consider the surface z = epsilon1cos(ωxx)cos(ωyy) with applied field(1,0,0)T.2.4.1 Numerical FormulationIn this three-dimensional setting, all three components of the field need tobe solved for both on the vacuum and superconducting sides. This higher-dimensional setting along with having two regions with different propertiesto consider makes this work considerably more difficult.We consider the formulation in parts. Our overall formulation (based on702.4. Finite Difference Program for General Sinusoidal SurfaceFigure 2.25: Formulation for three-dimensional code.the subsequent subsections) is displayed in figure 2.4.1.Vacuum Side InteriorWithin the vacuum we need to simultaneously satisfy∇×b = 0 and∇•b = 0.We make use of the vanishing curl to write b = ∇g where g : R3 →R is ascalar function (which for our purposes we will assume is at least C2).This then allows us to replace ∇•b = ∇•(∇g) with triangleg = 0.Superconducting Side InteriorIn the superconductor we require that both∇•b = 0712.4. Finite Difference Program for General Sinusoidal Surfaceandtriangleb = b. (2.23)In the interior we choose to only impose the Laplacian condition andimpose the divergenceless condition on the boundary. This is noted below.As the divergence and vector Laplacian operators commute i.e. ∇•triangle =triangle∇•then we can take the divergence of (2.23) to get∇•triangleb = ∇•b = triangle∇•bso that triangleq = q where q = ∇•b.If q = 0 on the boundary of the superconducting region and at z = +∞with triangleq = q, then q = 0 everywhere in the interior. Thus, we aim tonumerically implement ∇•b = 0 on the interface and at z = +∞.Boundary Conditions at z = ±∞Given that b = ∇g = (∂xg,∂yg,∂zg)T, along z = −∞ it is only possibleto specify ∂xg and ∂yg. We recall, however, that in the asymptotic analysisthat no solutions were possible if the third component of b were nonzero atz = −∞. As it turns out, we don’t need to impose the value of the thirdcomponent at z = −∞ and our results are still consistent.To impose b(x,y,−∞) = (1,0,0)T, it would at first seem appropriate tochoose g(x,y,−∞) = x; however, given that our numerical formulation isbased upon periodic boundary conditions this cannot work because x is nota periodic function (of x).We introduce ˜g = g−x. Then b = ∇˜g + (1,0,0)T and triangleg = triangle˜g.This b still satisfies all the PDEs and boundary conditions. If triangleg = 0then ∇•b = ∇•∇(˜g + (1,0,0)T) = triangleg = 0. And given that b = ∇g we722.4. Finite Difference Program for General Sinusoidal Surfacemust have that ∇×b = 0.There is a potential concern that this ˜g variable may not be periodic.However, given that the b must be periodic, we can show that ˜g is, too. Wegive this proof in appendix D.Our far-field vacuum boundary conditions are ˜g = 0 at z = −∞. At z =∞, we need to impose both b = 0 and ∇•b = 0. We impose the zero condi-tions on the first two components directly i.e. b1(x,y,∞) = b2(x,y,∞) = 0.This also forces ∂xb1(x,y,∞) = ∂yb2(x,y,∞) = 0.We now focus on the divergence condition and demand that∂zb3(x,y,∞) =0. Given the analysis in the asymptotic work, we see that the Fourier com-ponents would involve exponentials of multiples of z and ∂z could only goto zero if all the exponentials decayed. Thus ∂zb3(x,y,z) → 0 as z → ∞also weakly imposes b3(x,y,∞) = 0.Interface ConditionsAt the interface we need to specify a scalar equation (representing thevacuum-side interface condition) and a vector equation (representing thesuperconducting side interface condition). We impose ∇•b = 0 on the su-perconducting side of the interface, and b = ∇˜g +(1,0,0)T on the interfaceas discussed in section 2.4.1.Coordinate TransformationAgain, we note that the most interesting phenomena occur near the interfaceand so, along with considering the depth relative to the interface, t = z −epsilon1cos(ωxx)cos(ωyy), we also stretch our coordinates as in the two dimensional732.4. Finite Difference Program for General Sinusoidal Surfacecode.We define our coordinates asρ = x,σ = y,andτ =av log(−t−bv) if t ≤ 0as log(t−bs) if t > 0.with av, bv, as, bs being determined by the numerical approximations to±∞.This provides us with a new set of differential operators in the (ρ,σ,τ)parameters. Much of the work is very similar as in the two dimensional codeso we summarize the results here:∂x = ∂ρ±epsilon1ωx sin(ωxρ)cos(ωyσ)ae−τ/abracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightξ(1)∂τ∂y = ∂σ ±epsilon1ωy cos(ωxρ)sin(ωyσ)ae−τ/abracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightξ(2)∂τ∂z = ±ae−τ/abracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightξ(3)∂τandtriangle = ∂ρρ +∂σσ ±2epsilon1ωx sin(ωxρ)cos(ωyσ)ae−τ/abracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightξ(4)∂τρ±2epsilon1ωy cos(ωxρ)sin(ωyσ)ae−τ/abracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightξ(5)∂τσ+[(epsilon1(ω2x +ω2y)ω2y cos(ωxρ)cos(ωyσ))(±ae−τ/a)+(epsilon12ω2x sin2(ωxρ)cos2(ωyσ) +epsilon12ω2y cos2(ωxρ)sin2(ωyσ) + 1)(−ae−2τ/a)]bracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightξ(6)∂τ+742.4. Finite Difference Program for General Sinusoidal Surface[(epsilon12ω2x sin2(ωxρ)cos2(ωyσ) +epsilon12ω2y cos2(ωxρ)sin2(ωyσ) + 1)(a2e−2τ/a)]bracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightξ(7)∂ττwhere a = av on the vacuum side and as on the superconducting side andwe choose − for the vacuum side and + for the superconducting side.In this transformation we have a flat interface in the new coordinates.The GridAnother difficulty in this numerical work is in setting up a grid. Because themagnetic field is given in terms of the gradient of the vacuum ˜g function, itis necessary to consider two grids that interlock at the interface.We are given ωx, ωy, and epsilon1.We define a parameter N describing the number of points in the ρ−and σ− directions and Mm and Mp describing how far the grid spans inthe vacuum and superconducing sides respectively. Given N, we proceed todefine:Our spacings, hρ = 2piNxωx, hσ = 2piNyωy, hτ = min{hρ,hσ}; and our stretch-ing parameters, ηv = NhτMm , ηs = NhτMp , bv = ηgMm+hτ2ηg−1 , hI =MpN , bs =ηsMp+hI2ηs−1 , av = βg+hτ2 , as = −βs, αv = av log(Mm−bv), βv = av log(hτ2 −bv),αs = as log(−bs), βs = as log(Mp −bs).The careful selection of these α and β parameters allows for the gridsto interlock in the right way. When ωx or ωy are zero we defined hρ or hσrespectively to be 1.We setρi = − piωx+ (i−1)hρ for 1 ≤ i ≤ N,752.4. Finite Difference Program for General Sinusoidal Surfaceσj = − piωy+ (j −1)hσ for 1 ≤ j ≤ N,ρi = ρi + hρ2 ,σj = σj + hσ2and if ωx and/or ωy are zero then we only choose three coordinates −1,0,1for the ρprimes and/or σprimes. We set Nρ to be the number of ρ− points and Nσ tobe the number of σ−points.We set Nv = [βv−αvhτ ]+1 and Ns = [βs−αshτ ]+1 with [ ] being the integerfloor function and define τNv−k+1 = βv − (k − 1)hτ for k = 1,...,Nv andτk = αs + (k−1)hτ for k = 1,...,Ns.In the end we have NxNy(Nv + 3Ns) variables and equations to define.We define the mapping K that takes four coordinates and maps them toa single number:K(i,j,k,lscript) =(k−1)NxNy + (j −1)Nx +i if lscript = 0NxNyNv + (lscript−1)NxNyNs + (k−1)NxNy + (j −1)Nx +i if lscript negationslash= 0where to impose the periodic boundary conditions, if i = 0 then it mustbe replaced by i = Nρ and if i = Nρ +1 it is replaced by 1 and similarly forj with Nσ.Based on this numbering, thei,j,andk describe theρ,σ,andτ positions.The number lscript = 0 describes the ˜g variable (the variable in the vacuumincluding its ghost points), and lscript = 1,2,3 describes the first, second, andthird components of the magnetic field respectively in the superconducting762.4. Finite Difference Program for General Sinusoidal SurfaceFigure 2.26: The mesh in the three-dimensional system. Note how ˜g and bmeshes interlock. The circles with crosses indicate the points occur deeperinto the page than the circles wit the dots. At each circle with a cross, threecomponents are specified. At each circle with a dot, the scalar value of ˜g isspecified.coordinates.To discretize t, the depth past the surface, we define:tk =tk = −e(τk+τk+1)/(2av) −bv for 1 ≤ k ≤ Nv −2tk = e(τk+2−Nv)/as +bs for Nv −1 ≤ k ≤ Nv +Ns −3.This imposes t1 ≈ −Mm + hI2 , tNv−1 = 0, and tNv+Ns−3 ≈ Mp. Figure2.26 illustrates the mesh.772.4. Finite Difference Program for General Sinusoidal SurfaceThe Discretized EquationsHere we present how the equations were handled. We will define u as theunknown vector, M as the matrix, and f as the load vector. The equationschange based upon the τ position so in each slice the equations hold for alli and j.We also discretize the ξ variables, like in the two-dimensional code, withindices (i,j,k,lscript). For example, ξ(2)i,j,k,0 = −epsilon1ωy cos(ωxρi)sin(ωyσj)ave−τk/avConditions at −∞ This corresponds to k = 1 and lscript = 0.We need to impose ˜g(ρ,σ,−∞) = 0. To do so we set fK(i,j,1,0) = 0 withMK(i,j,1,0),K(i,j,1,0) = 1.Conditions at +∞ This corrresponds to k = Ns and lscript = 1,2,3.We need to impose b1(ρ,σ,∞) = b2(ρ,σ,∞) = 0 and ∂zb3(ρ,σ,∞) =ξ(3)∂τb3 = 0.This is an equation for the points with τ coordinate number Ns. Weimpose fK(i,j,Ns,lscript) = 0.The point at ∞ is between regular grid points and ghost-points. Weimpose an average condition for the first two components (imposing theiraverage be zero) and a regular second-order derivative condition for the thirdcomponent:MK(i,j,Ns,lscript),K(i,j,Ns−1,lscript) = 12MK(i,j,Ns,lscript),K(i,j,Ns,lscript) = 12782.4. Finite Difference Program for General Sinusoidal Surfacefor lscript = 1,2, andMK(i,j,Ns,3),K(i,j,Ns−1,3) = − 1hτξ(3)i,j,Ns,3MK(i,j,Ns,3),K(i,j,Ns,3) = 1hτξ(3)i,j,Ns,3.Vacuum Laplacian Here, k = 2,...,Nv −1 and lscript = 0.The Laplacian is zero in the vacuum region so fK(i,j,k,0) = 0.The discretization methods for mixed partial derivatives are the samehere as for the two-dimensional code. For the matrix values, we have:MK(i,j,k,0),K(i±1,j,k,0) = 1h2ρMK(i,j,k,0),K(i,j±,k,0) = 1h2σMK(i,j,k,0),K(i,j,k±1,0) = ±ξ(6)i,j,k,02hτMK(i,j,k,0),K(i,j,k,0) = −2h2ρ− 2h2σ− 2ξ(7)i,j,k,0h2τMK(i,j,k,0),K(i+1,j,k+1,0) = ξ(4)i,j,k,04hρhτMK(i,j,k,0),K(i−1,j,k−1,0) = ξ(4)i,j,k,04hρhτ792.4. Finite Difference Program for General Sinusoidal SurfaceMK(i,j,k,0),K(i+1,j,k−1,0) = −ξ(4)i,j,k,04hρhτMK(i,j,k,0),K(i−1,j,k+1,0) = −ξ(4)i,j,k,04hρhτMK(i,j,k,0),K(i,j+1,k+1,0) = ξ(5)i,j,k,04hσhτMK(i,j,k,0),K(i,j−1,k−1,0) = ξ(5)i,j,k,04hσhτMK(i,j,k,0),K(i,j+1,k−1,0) = −ξ(5)i,j,k,04hσhτMK(i,j,k,0),K(i,j−1,k+1,0) = −ξ(5)i,j,k,04hσhτSuperconductor Laplacian Here, k = 2,...,Ns −1 and lscript = 1,2,3.This is identical to the discretization of the vacuum Laplacian aboveexcept that the 0 is replaced by lscript = 1,2,3 and the slight change in nowhaving and extra −1 in these entries:MK(i,j,k,0),K(i,j,k,0) = −2h2ρ− 2h2σ− 2ξ(7)i,j,k,0h2τ −1.Zero Divergence on the Interface This is the scalar equation corre-sponding to the ˜g-ghost points. We have k = Nv and lscript = 0. We chose toimpose zero divergence at the ˜g−ghost points.These points, however, are defined between the b−points. To resolve thisissue, we define the divergence at ghost points by taking averages. Consider-ing each spatial direction ˆx, ˆy, ˆz individually, every ghost point is surrounded802.4. Finite Difference Program for General Sinusoidal Surfaceby four edges along which a derivative can be taken. Taking the average ofthese four discrete derivatives gives an approximation to the derivative (inone of the directions) at the ghost point.For example, to compute the z− derivative of b1 at (ρi,σj,τNv,0) wewould compute the average of the discrete derivatives ξ(3)i,j,3/2,1hτ (uK(i,j,2,1) −uK(i,j,1,1)), ξ(3)i,j,3/2,0hτ (uK(i−1,j,2,1)−uK(i−1,j,1,1)),ξ(3)i,j,3/2,0hτ (uK(i,j−1,2,1)−uK(i,j−1,1,1)),ξ(3)i,j,3/2,0hτ (uK(i−1,j−1,2,1) −uK(i−1,j−1,1,1)) where ξ(3)i,j,3/2,1 is defined by the av-erage 12(ξ(3)i,j,1,1 +ξ(3)i,j,2,1). Referring back to the figure 2.26 should clarify.The divergence is zero and we are looking at the equations for the vacuumghost points. Therefore we set fK(i,j,Nv,0) = 0.The entries of M corresponding to the terms in ∂xb1 are shown belowand the other terms are similar. For the matrix M, we impose:M(K(i,j,Nv,0),K(i,j,2,1)) = 14hρ+ξ(1)i,j,3/2,14hτM(K(i,j,Nv,0),K(i,j,1,1)) = −14hρ+ξ(1)i,j,3/2,14hτM(K(i,j,Nv,0),K(i−1,j,2,1)) = 14hρ+ξ(1)i−1,j,3/2,14hτM(K(i,j,Nv,0),K(i−1,j,1,1)) = −14hρ+ξ(1)i−1,j,3/2,14hτM(K(i,j,Nv,0),K(i,j −1,2,1)) = 14hρ+ξ(1)i,j−1,3/2,14hτ812.4. Finite Difference Program for General Sinusoidal SurfaceM(K(i,j,Nv,0),K(i,j −1,1,1)) = −14hρ+ξ(1)i,j−1,3/2,14hτM(K(i,j,Nv,0),K(i−1,j −1,2,1)) = 14hρ+ξ(1)i−1,j−1,3/2,14hτM(K(i,j,Nv,0),K(i−1,j −1,1,1)) = −14hρ+ξ(1)i−1,j−1,3/2,14hτwhere the factor of 14 comes from taking averages.The Gradient Condition This corresponds to k = 1 and lscript = 1,2,3.The condition is that b = ∇˜g + (1,0,0)T or that ∇˜g −b = (−1,0,0)T.These are equations for b on the interface.Therefore we set fK(i,j,1,1) = −1 and fK(i,j,1,2) = fK(i,j,1,3) = 0.As with the divergence, we need to take averages in finding the deriva-tives. For an illustration, we will show the equations for the third componentof b. They correspond to setting the entries of M as follows:MK(i,j,1,3),K(i,j,Nv,3) =ξ(3)i,j,Nv−1/2,04hτMK(i,j,1,3),K(i,j,Nv−1,3) =−ξ(3)i,j,Nv−1/2,04hτMK(i,j,1,3),K(i−1,j,Nv,3) =ξ(3)i−1,j,Nv−1/2,04hτMK(i,j,1,3),K(i−1,j,Nv−1,3) =−ξ(3)i−1,j,Nv−1/2,04hτ822.4. Finite Difference Program for General Sinusoidal SurfaceMK(i,j,1,3),K(i,j−1,Nv,3) =ξ(3)i,j−1,Nv−1/2,04hτMK(i,j,1,3),K(i,j−1,Nv−1,3) =−ξ(3)i,j−1,Nv−1/2,04hτMK(i,j,1,3),K(i−1,j−1,Nv,3) =ξ(3)i−1,j−1,Nv−1/2,04hτMK(i,j,1,3),K(i−1,j−1,Nv−1,3) =−ξ(3)i−1,j−1,Nv−1/2,04hτMK(i,j,1,3),K(i,j,1,3) = −1.We similarly use the notation ξ(3)i,j,Nv−1/2,0 = 12(ξ(3)i,j,Nv−1,0 +ξ(3)i,j,Nv,0).Finding the Results This averaging actually brings about a null vectorfor M when N is even. As a result, the program is restricted to odd valuesof N. We explain how this null vector comes about in appendix E for theflat geometry. The same ill-conditioning also occurs for rough geometries.Given the system Mu = f we found u = M−1f (by implementing adirect solve routine for sparse matrices in Matlab). This gives us a magneticfield on the superconducting side, but only gives us ˜g on the vacuum side.To get b on the vacuum side, we took the discretized gradient of ˜g justas above, thereby finding the magnetic field at values at the centres of cubeswith ˜g values. We also added 1 to the first component.Where necessary, an average field magnitude at depth tk can be foundby averaging the field magnitude over i and j with k fixed.832.4. Finite Difference Program for General Sinusoidal Surface2.4.2 ValidationTo test that the code results are accurate, we compare the results to theasymptotically computed solutions.As our problem of interest is a three-dimensional vector problem, thesimulation is extremely difficult to run, due to memory constraints. Thislimitation stems from the use of a direct-solver for the linear system.As a result, in the fully three-dimensional problem, the resulting systemof equations can only be solved with values of N up to about 21. When thiscode has either ωx or ωy zero, we can take N much larger, but still nowherenear the values possible in the two-dimensional code. In all computations,we tried to take N as large as possible.In terms of Mp and Mm, we generally choose Mp = 9.25 and set Mm =Mp/radicalBigω2x +ω2y. This choice of Mm ensures that even if epsilon1 were as large as1 that the decay of the first order perturbation (which, in the asymptoticcase, decays with exp(radicalBigω2x +ω2yt) would be of the same order as the errorat Mp.The most basic check is how the program behaves with a flat interface.For epsilon1 = 0, the field magnitude should be identically 1 on the vacuum side anddecay exponentially with e−t on the superconducing side. Taking N = 11we see the plot in figure 2.27Comparison with Asymptotic SolutionsWe now compare the numerical approximations to the asymptotic solutionfor different values of epsilon1.842.4. Finite Difference Program for General Sinusoidal SurfaceFigure 2.27: Visual confirmation that the three-dimensional program cor-rectly handles the flat interface N = 11, and ωx = ωy = pi.Our first-order asymptotic solution satisfies basy,1 = b(0) + epsilon1b(1) = bex +O(epsilon12) where b(0) and b(1) are the zeroth and first-order terms of the asymp-totic expansion and bex is the exact solution.Our numerical solution bnum satisfies bnum = bex + O(h2τ) + O(λ) whereλ is the far-field error. Taking the difference of basy,1 and bnum and dividingby epsilon1maxi||b(1)i ||∞ and taking norms yields an error:E = || bex −basy,1epsilon1maxi||b(1)i ||∞||∞ = O(epsilon1) +O(h2τ/epsilon1) +O(λ/epsilon1).Provided h2τ is small enough (by taking N large) and λ is small, thenthis should be a linear function of epsilon.To verify the linear nature, we plot E vs epsilon1 to see that E varies roughly852.4. Finite Difference Program for General Sinusoidal SurfaceFigure 2.28: Left: resolution of first-order asymptotic term for first geometrywith ω = pi. Centre: resolution of first-order asymptotic term for secondgeometry with ω = pi. Right: resolution of first-order asymptotic term forthird geometry with (ωx,ωy) = (pi,pi).linearly with epsilon1 and that the E ↓ 0 as epsilon1 ↓ 0. From this, we see the programis remarkably accurate. Even for epsilon1 ≈ 0.3 the program can discern the first-order term to within 15%. See figure 2.28.As in section 2.3.2, the second-order asymptotic solution basy,2 should benear the numerical solution to within an error of order O(h2τ)+O(epsilon13)+O(λ).Being able to test this higher order of convergence is extremely delicate (dueto memory constraints). When ωx,ωy negationslash= 0, we can only take N as large as21. Another difficulty is that the magnetic field is essentially comprised ofsix parts: three components on two sides of the interface. Each part has itsown terms and could presumably have optimal ranges over which the O(epsilon13)term could be detected by using second-order asymptotics.We begin this test by taking the three-dimensional code with the condi-tions (ωx,ωy) = (pi,0) so that we can test how it converges with respect tothe results of section 2.2.4. This also means we can take N much larger (inthis case we 41). We allow epsilon1 to range from 0.05 to 0.40 with increments of0.05. From this we obtain the asymptotic errors for each side of the inter-862.4. Finite Difference Program for General Sinusoidal SurfaceTable 2.1: Estimated orders of convergence for (ωx,ωy) = (pi,0).Field Part Range of epsilon1 Slopeb1 vacuum [0.1,0.25] 2.77b1 superconductor [0.2,0.4] 3.08b3 vacuum [0.05,0.2] 2.71b3 superconductor [0.1,0.25] 2.86Table 2.2: Estimated orders of convergence for (ωx,ωy) = (pi,pi).Field Part Range of epsilon1 Slopeb1 vacuum [0.2,0.4] 2.94b1 superconductor [0.15,0.4] 3.03b2 vacuum [0.25,0.4] 2.98b2 superconductor [0.15,0.25] 2.71b3 vacuum [0.15,0.35] 2.74b3 superconductor [0.1,0.2] 2.50face, E = ||bnumi−basy,2i||∞ (i = 1,2,3), and we can perform a least-squaresregression for E = c1 logepsilon1+c0. We do our fitting on the range of data thatseems to best detect O(epsilon13). We tabulate the results of c1 in table 2.1.Seeing these very convincing results, we also test the case where (ωx,ωy) =(pi,pi), keeping the same range of epsilon1 values as above. Again, we perform aleast-squares linear fit to the best-behaved subset of the data. Results aretabulated in table 2.2.These results show a convincing agreement between the computationaland asymptotic results in the three-dimensional setting, which validates bothapproaches.2.4.3 ResultsChoosing N, we show the plots of the field profile.We begin by selecting (ωx,ωy) = (2pi,2pi) with epsilon1 = 0.05 to see the profiles872.4. Finite Difference Program for General Sinusoidal Surface(shown in figure 2.29) are similar to that of the asymptotics.We also have the freedom to choose much larger values of epsilon1 and the ω’s.To examine how the profile changes when the frequencies are out of pro-portion, we fix epsilon1 = 0.1 and select (ωx,ωy) = (pi,8pi) and (ωx,ωy) = (pi,8pi).The results are plotted together along with the flat geometry solution infigure 2.30. We see the mean field magnitudes in both perturbed geometriesexceed that of flat geometry after a certain depth.With ωx lessmuch ωy we see a similar shape to that of the experimental results(see figure 2.2) although the dead layer is considerably smaller. With ωx greatermuchωy we see a shape that resembles the experimental plot quite closely, butthe decay in magnitude begins even before the superconducting interface.The experimental results in figure 2.30 only measure the field profile inthe superconductor. However, it is possible to coat the superconductorswith silver to measure the field outside [8].From this, we observe that by rotating the applied magnetic field (as-suming that we have a single value of λ that does not depend on orientation),the field profile changes purely because the roughness differs in the two di-rections. However, the amount by which the field profile differs from the flatgeometry is not large enough to account for the experimental results unlessepsilon1 could be larger i.e. the surface is rougher than experimentalists believe.As final plot for this section, we show how well the asymptotics andnumerics agree. With epsilon1 = 0.1 and (ωx,ωy) = (8pi,8pi), their predicted meanfield profiles are plotted together in figure 2.31.882.4. Finite Difference Program for General Sinusoidal SurfaceFigure 2.29: Top left: a profile of b1 from peak and valley. Top right: aprofile of b2 with (x,y) = (−0.27,−0.27). Bottom left: profile of b3 with(x,y) = (−0.27,−0.50). Bottom right: difference between average fieldprofile in perturbed geometry and flat geometry. All figures computed withepsilon1 = 0.05, and (ωx,ωy) = (2pi,2pi).892.5. Conclusions of Superconductor ModelingFigure 2.30: The field profiles for different roughness orientations with epsilon1 =0.1.2.5 Conclusions of Superconductor ModelingRecent experiments have suggested the need for more sophisticated modelsof a superconductor. Our work examined the notion of surface roughness andthe effects this would have on the magnetic field profile in an experimentalsetting.Through analytic and numerical techniques, we have shown that surfaceroughness does indeed play a role in perturbing the field magnitude. Notonly does it perturb the field, but the orientation of the geometry itself playsa big role in the nature of the perturbation.Here we briefly summarize the most interesting of these results.902.5. Conclusions of Superconductor ModelingFigure 2.31: The asymptotic and numeric mean field profiles. The two arequite close given the large epsilon1 and ω values.2.5.1 More Complicated Field BehaviourThe standard model with a flat superconductor and parallel applied mag-netic field yields a constant solution within the vacuum region and a fieldpointing only in one direction decaying exponentially once within the super-conductor.Our results show that in general the field magnitude would not be con-stant in the vacuum and that once within the superconducting region themagnitude could decay differently than a purely exponential decay - decayswith multiple length scales and field components with sinusoidal variationsin the longitudinal directions.The individual components can rise above or fall below their value at912.5. Conclusions of Superconductor Modeling−∞. In addition, the field components themselves become mixed, and thefield is no longer pointing in a single direction throughout the experimentalregion.2.5.2 Orientation of the Roughness Affects the ProfileFor a given roughness that’s not equal in both spatial dimensions, dependingon whether the applied field is pointing in a rough or smooth direction, theamount by which it decays in the vacuum region could be significant ornearly negligible.If the experiments were precise enough, this would mean that if a su-perconductor did have a rough surface (which was not equal in both spatialdirections), then by changing the field orientation by 90◦ (assuming the pen-etration depth remained nearly the same) then the field profile should bedifferent. Having to take into account the anisotropy in the sample wouldmake this more difficult.2.5.3 Effective Dead LayerBased on an asymptotic regime, the effective dead layer (a length over whichwe can think of the field magnitude remaining constant beyond the interface)is of the order epsilon12. Recall, epsilon1 is the scaled size of the surface roughness definedin section 2.2.1.The actual dependence on the spatial frequencies is quite complicated,and in some cases the effective dead layer is actually negative.Experimentally the dead layer appears to be on the order of 0.05. Ourresults show that a surface roughness of 0.05 cannot account for the dead922.5. Conclusions of Superconductor Modelinglayers of the size shown in figure 2.30. However, results for larger roughnessesare qualitatively similar to the experimental profiles.93Chapter 3Fuel Cell Modeling3.1 Modeling of Gas Diffusion in Fuel CellsHydrogen fuel cells hold a promising future for being able to produce cleanenergy. If Hydrogen and Oxygen are used in the fuel cell, they react electro-chemically producing water (as the only by-product) and energy. For thisto take place, these gases must diffuse from channels to reaction sites. Be-ing able to set ideal running conditions for hydrogen fuel cells to maximizeefficiency requires an in-depth understanding of the efficiency of the gas dif-fusion among other things. Running experiments can be costly and in recentyears a lot of emphasis has been placed on modeling fuel cells numerically.The question then arises as to how to model this gas diffusion.Fickian diffusion, wherein a species diffuses in proportion to the gradi-ent of its mole fraction, has traditionally been chosen. However it is be-lieved that a more sophisticated (and more complex) formulation knownas Maxwell-Stefan diffusion yields more physically accurate results. Thisdiffusion model takes into account the interspecies competition in diffusion.Although the theory behind these models is different, they are used to modelthe same phenomena.Due to the interaction terms and mathematical subtleties, it is not clear943.1. Modeling of Gas Diffusion in Fuel Cellshow to best pose a Maxwell-Stefan diffusion problem, and exactly how muchit differs from Fickian diffusion in situations relevant to fuel cells is unknown.These issues are explored in this chapter.We begin this portion of the paper with a brief overview of the releventPhysics involved in PEM fuel cells and proceed to develop our model fromthere.3.1.1 PEM Fuel Cell OverviewA Polymer Electrolyte Membrane (PEM) fuel cell generates power through areaction between Hydrogen (H2) and Oxygen (O2) producing water (H2O).The simplest cartoon is that the two gases enter through two separate chan-nels and diffuse through a concentration gradient to arrive at catalyst sitesthat facilitate the reaction [9]. See figure 3.1.1.Hydrogen gas present in the Hydrogen flow channel enters the gas diffu-sion layer (GDL) anode where it diffuses through a concentration gradientuntil it reaches a catalyst layer (typically comprised of Platinum or a Plat-inum alloy). Here, it oxidizes with the reaction H2 → 2H+ + 2e−. Theelectrons become part of the electric current drawn from the fuel cell andthe protons diffuse through the polymer electrolyte membrane (PEM) untilreaching the cathode Oxygen catalyst layer.The other reacting gas is Oxygen which diffuses from the Oxygen flowchannel through the GDL cathode. Once reaching the catalyst, an Oxygenmolecule combines with four protons (from the Hydrogen side) and 4 elec-trons (generated through the current) to produce water under the reactionO2 +4e− +4H+ → 2H2O. This reaction is also enhanced by Platinum or a953.1. Modeling of Gas Diffusion in Fuel CellsFigure 3.1: A cross-section of a PEM fuel cell. Hydrogen and Oxygen diffuseprimarily in the XY−plane, but there is a diffusion in the Z−direction aswell.Platinum alloy.The properties of the GDL have a large impact on the diffusion. TheGDL is not an open channel. Instead, it is a thin layer of teflonated carbonfibre paper, which can be considered a porous media, through which thegases are conducted. The efficiency of this transport depends on how “open”these pathways are. These pathways make for a greater overall distance forthe molecules to travel. They cannot travel straight; instead, they must passthrough tortuous pathways to reach the catalyst site.3.1.2 The ModelIn our work here we will look at the GDL of the cathode. We will model theconcentration profile of the gases and explore how Fick and Maxwell-Stefan963.1. Modeling of Gas Diffusion in Fuel Cellsdiffusion laws differ.We will consider the presence of three gas species, O2, H2Ovap (watervapor), and N2 in the GDL at steady state transport. We approximatethe system as isothermal (the temperature is assumed to be constant). Wealso approximate the system as one-dimensional: the gases have prescribedconcentrations at the channels and travel in the X−direction to reach thecatalyst layer where the fluxes are known.We consider the different concentrations of gas species Ci, i = 1,2,3,where 1 indicates O2, 2 indicates H2Ovap, and 3 indicates N2. Together weconsider the vector of concentrations C = (C1,C2,C3)T.The total molar concentration is the sum of the concentrations of theindividual speciesC1+C2+C3 which we will denote by||C||.Unless otherwisestated, the norm will always represent the L1-norm of a vector. The mass-density ρi of a species i is found by multiplying its molar mass Mi by itsmolar concentration Ci.The temperature (assumed constant) will be denoted by T, the totalpressure by P, and the mass-averaged velocity by U.Important physical constants in this setting are the ideal gas constant,R, the porosity, epsilon1, the viscosity, µ, and the permeability, κ. In Fickian dif-fusion, D will represent the diffusion coefficient, and in the Maxwell-Stefanformulation, Dij, (i,j) ∈{1,2,3}2 will be the binary diffusivities.The important fluxes are the molar flux with respect to the molar-averaged velocity J∗ and the molar flux with respect to the mass-averagedvelocity J. Both J and J∗ are in general rank 2 tensors (in one dimensionthey are vectors).973.1. Modeling of Gas Diffusion in Fuel CellsIn our model, the channel concentrations are specified at position X = 0so that C(0) = C∗ is known. A known current is drawn from the reac-tions and based upon this current, under standard operating conditions, thefluxes, Ni of the individual species at the catalyst layer (at X = L) areknown. We will now determine our parameters of interest.3.1.3 Standard Operating Conditions of a Fuel CellHere, based upon standard operating conditions (see references [9], [10], [11],and [12]), we compute our channel concentrations and fluxes.The main equation we need here is the ideal gas law ((3.2) in the nextsection). We take the pressure within the fuel cell channels to be 3 atm(approximately 3×105 Pa). The temperature is taken to be 350 K. Basedon this we can compute the total molar concentration of the gases as ap-proximately 100 mol m−3. These are typical conditions for fuel cells beingdeveloped for the automotive sector [11].At this same temperature, the saturation pressure of water vapor is3.8×104 Pa. This corresponds to a saturation concentration, Csat of 13 molm−3. If we assume 75% humidity then within the fuel cell, the molar densityof water vapor is approximately 10 mol m−3.We then assume the remaining 90 mol m−3 are comprised of Oxygen at21% and Nitrogen at 79%. Our channel concentrations are: C1 = 19 molm−3, C2 = 10 mol m−3, and C3 = 71 mol m−3.A typical current density drawn from a fuel cell is 1 A cm−2 which is 104C s−1 m−2. Based on the reaction in the GDL, we see that for every fourelectrons consumed, one Oxygen molecule is consumed. Taking the current983.1. Modeling of Gas Diffusion in Fuel CellsTable 3.1: Physical constants in our gas diffusion model.constant valueC1 19 mol m−3C2 10 mol m−3C3 71 mol m−3D 8×10−6 m2 s−1D1,2 1.24×10−5 m2 s−1D1,3 1.04×10−5 m2 s−1D2,3 1.23×10−5 m2 s−1epsilon1 0.74F 9.6×104 C mol−1κ 10−12 m2L 10−4 mM1 32 g mol−1M2 18 g mol−1M3 28 g mol−1N1 2.6×10−2 mol m−2 s−1N2 −5.2×10−2 mol m−2 s−1N3 0 mol m−2 s−1µ 2.24×10−5 kg m−1 s−1R 8.314 kg m2 mol−1 s−2 K−1T 350 Kdensity and dividing by 4 F (where F is Faraday’s constant, 9.6 × 104 Cmol−1) we find the Oxygen flux at the catalyst layer to be N1 = 0.026 molm−2 s−1.The water vapor flux should be in the opposite direction and have twicethe magnitude. So N2 = −0.052 mol m−2 s−1. Nitrogen does not react soits flux is zero.Our particular choice of parameters are summarized in table 3.1. Asidefrom the constants we used or obtained above, we obtained the other con-stants from references [9] and [11].In the next section we will review the relevant gas diffusion equations.993.1. Modeling of Gas Diffusion in Fuel Cells3.1.4 Diffusion EquationsHere we give the formulations of Fick and Maxwell-Stefan diffusion. Al-though our model is one-dimensional, the equations here are given in fullgenerality. For further reference, we suggest reference [13].Einstein summation convention is used throughout in that repeated in-dices in products are summed over (aijbj =summationtextj aijbj for example).Basic EquationsIn both the Fick and Maxwell-Stefan settings, we will assume Darcy’s law,which applies to porous media [9], and the Ideal Gas law. Respectively, theyare stated below:U = − κepsilon1µ∇P (3.1)P = ||C||RT. (3.2)These equations can be combined nicely by substituting (3.2) into (3.1)to yield a modified Darcy Law:U = −κRTepsilon1µ ∇||C||. (3.3)In later sections we will define the constant σ = KRTepsilon1µ .The mass-averaged velocity is given by U = ρiVi||ρ||. The molar-averagedvelocity is given by U∗ = CiVi||C||. Here the Vi are the velocities of the individualspecies.1003.1. Modeling of Gas Diffusion in Fuel CellsFickian DiffusionFick’s law describes diffusive fluxes. The formulation we adopt is the sameas that implemented by Berg et al. [11]. For each species of gas, we have:∂tCi +∇•(CiU −D||C||∇( Ci||C||)) = 0. (3.4)The term CiU is a convective term - the gas species move with the mass-averaged velocity. The term D||C||∇( Ci||C||) is a diffusive term, where diffu-sion is assumed to be driven by the gradient of the mole fractions (Ci/||C||).Fickian diffusion cannot actually be reconciled with Darcy’s law as stated[14]. By summing (3.4) over i we find:∂t||C||+∇•(||C||U −D||C||∇(||C||||C||)) = 0and the latter term in the diffusive term is actually zero (∇(1) = 0.) There-fore we see that ∂t||C|| = −∇•(||C||U) which only makes sense if U is themolar-averaged velocity U∗ (different to how Darcy’s law is posed where Uis the mass-averaged velocity).In spite of this, many researchers continue to use Fickian diffusion inconjunction with Darcy’s law [9]. As we see next, Maxwell-Stefan diffusionis compatible with Darcy’s law.1013.1. Modeling of Gas Diffusion in Fuel CellsMaxwell-Stefan DiffusionThe Maxwell-Stefan diffusion law is expressed below:∂tCi +∇•(CiU +Ji) = 0, (3.5)where we define J through a series of transformations. To start,Aij ˜Jj = ∇( Ci||C||) (3.6)gives a relationship between the molar diffusive flux ˜J (with respect to anarbitrary velocity) and the gradients of the mole fractions, whereAij(C) = 1||C||2summationtextlscriptnegationslash=iClscriptDilscript if i = j− CiDij if i negationslash= j.The Dij for (i,j) ∈ {1,2,3}2 are binary diffusivities. We will denoteˆD = maxij Dij.The matrix A is not invertible. The vector (C1,C2,C3)T spans its nullspace. To find the general solution of (3.6) we will pick a base point,(J∗1,J∗2,J∗3), by using a projection and replacing one of the rows of A bya multiple of (1,1,1) and then add an arbitrary multiple of ζ(C1,C2,C3)Tto the solution. We find:( ˜J1, ˜J2, ˜J3)T = (PA+B)−1P∇( C||C||)bracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipupright(J∗1,J∗2,−J∗1−J∗2)T+ζ(C1,C2,C3)T1023.1. Modeling of Gas Diffusion in Fuel CellswhereP =1 0 00 1 00 0 0andB = 1ˆD||C||0 0 00 0 01 1 1.The factor of 1ˆD||C|| within B may seem strange but it is needed to keepconsistency with the dimension and scaling of A. The usefulness of this willbecome apparent in the asymptotic and numerical work to come (section3.2.1).Also note that the −J∗1 −J∗2 for the third component of the particularsolution comes from our requirement that J∗1 +J∗2 +J∗3 = 0 (i.e. that J∗ isthe molar flux with respect to the molar-averaged velocity).We have found the molar-diffusive fluxes (which have an arbitrary pa-rameter). We would like (3.5) to involve a mass-averaged velocity [14]. Inthis case, if we multiply (3.5) by Mi and add we find:∂tρ+∇•(ρU +MiJi) = 0where ρ is the total density, CiMi (the sum).For a mass-averaged U we require that MiJi = 0. For a particular choiceof ζ, the ˜J we found is the molar flux with respect to the mass-averagedvelocity, J. So we now find ζ by solving1033.2. Asymptotic FormulationM1(J∗1 +C1ζ) +M2(J∗2 +C2ζ) +M3(−J∗1 −J∗2 +C3ζ) = 0.We obtain ζ = (M3−M1)J∗1+(M3−M2)J∗2ρ . Using this ζ allows us to write:J1J2J3=S11 S12 0S21 S22 0−M1M3 −M2M3 0J∗1J∗2J∗3where Sij = δij − CiMjρ (1− M3Mj ) as noted by Stockie et. al. [9].Overall we can find J as shown below:J1J2J3=S11 S12 0S21 S22 0−M1M3 −M2M3 0(PA+B)−1Pbracehtipupleft bracehtipdownrightbracehtipdownleft bracehtipuprightQ∇( C||C||).The Maxwell-Stefan system is now expressed as:∂tCi +∇•(CiU +Qij(C)∇( Cj||C||)) = 0. (3.7)3.2 Asymptotic FormulationOur analysis will be in a one-dimensional setting so that ∇(•) = •prime with primebeing the spatial derivative. We begin by nondimensionalizing the systemsof equations as we explained in section 1.2.1043.2. Asymptotic Formulation3.2.1 NondimensionalizationFickian DiffusionLet our spatial coordinate be X. By combining (3.3) and (3.4), we obtainthe equation∂tCi + (−σCi||C||prime −D||C||( Ci||C||)prime)prime = 0 (3.8)which, in steady state (i.e. with ∂tC = 0), and after integrating with respectto X yields−σCi||C||prime −D||C||( Ci||C||)prime = Ni (3.9)where the N’s are fluxes.Defining the bars as dimensional quantities with a representative scaleas in section 1.2, we write X = ¯Xx,C = ¯Cc, and N = ¯Nn. Substitutingthis into equation and dividing by ¯N yields−σ¯C2¯X ¯Nci||c||prime −D ¯C¯X ¯N||c||(ci||c||)prime = n (3.10)where prime now represents the derivative with respect to x.By equating the dimensions of the terms on the left side of (3.10) wefind σ ¯C2¯X ¯N = D ¯C¯N ¯X which can be satisfied by setting ¯C = Dσ .We will choose the length scale ¯X = L as the length of the diffusionlayer. Requiring that there be no dimension on the left side of (3.10) meansthat D ¯CL ¯N = 1, which can be satisfied in selecting ¯N = D2σL.1053.2. Asymptotic FormulationFigure 3.2: Our one-dimensional model of the gas diffusion layer.With these substitutions, we find−ci||c||prime −||c||( ci||c||)prime = n.We define γ = 1||c∗|| (where c∗ = c(0)) and rescale with ˜c = γc. We alsodefine the number epsilon1 = γ2||n|| and the vector ˜n = n||n|| so that γ2n = epsilon1˜n with||˜n|| = 1. Writing (3.10) with the tilde variables and removing the tildesyields.ci||c||prime −γ||c||( ci||c||)prime = epsilon1ni (3.11)for x ∈ [0,1] and with |c(0)| = |c∗| = 1. Note that epsilon1 is not the porosity here.Our problem is displayed in figure 3.2.It turns out epsilon1 is very small, and so is γ, although γ greatermuch epsilon1. In this case,the Fickian term is a small perturbation from the Darcy term (which alonecould not be solved).Our parameters are given in table 3.2.1063.2. Asymptotic FormulationTable 3.2: Nondimensionalized and rescaled parameters and variables forFick diffusion.Paramater Valuec∗1 0.19c∗2 0.10c∗3 0.71||c∗|| 1n1 0.333n2 −0.667n3 0γ 4.56×10−4epsilon1 4.44×10−6Maxwell-Stefan DiffusionAlthough the problems have very different structures, the nondimensional-ization procedure is very similar. We begin by combining (3.1) and (3.5) toobtain∂tCi + (−σCi||C||prime +Qij(C)( Ck||C||)prime)prime = 0 (3.12)which becomes−σCi||C||prime +Qij( Cj||C||)prime = Ni (3.13)in steady-state.We again use the ¯• notation to refer to a dimensional parameter with arepresentative scale for • and write C = ¯Cc, X = ¯Xx, and N = ¯Nn. Notingthat (PA+B) ∝ ( ˆD||C||)−1 so that (PA+B)−1(C) ∝ ˆD||C|| and applying1073.3. Asymptotic AnalysisTable 3.3: Nondimensionalized and rescaled parameters and variables forMaxwell-Stefan diffusion.Paramater Valuec∗1 0.19c∗2 0.10c∗3 0.71||c∗|| 1n1 0.333n2 −0.667n3 0γ 7.06×10−4epsilon1 4.44×10−6the same techniques (with ˆD replacing D) that led to (3.10) we find−ci||c||prime +Qij(c)( cj||c||)prime = ni. (3.14)Rescaling as in the Fickian case we can write equation 3.14 as:−ci||c||prime +γQij(c)( cj||c||)prime = epsilon1ni (3.15)where |c(0)| = |c∗| = 1 and where we consider x ∈ [0,1]. Again, epsilon1 andγ are small with γ greatermuch epsilon1. Our parameters are given in table 3.3. Figure 3.2applies here as well.3.3 Asymptotic Analysis3.3.1 Asymptotic Analysis of Fick DiffusionBy expanding the derivatives of (3.11) we have−ci||c||prime−||c||( ci||c||)prime = −ci(summationtextj cprimej)−γcprimei + γ ci||c||(summationtextj cprimej) = epsilon1ni. It is possible to rearrange this in the form of1083.3. Asymptotic AnalysisMFcprime = epsilon1n for a matrix MF.We see−cisummationtextj cprimej = Dijcprimej whereDij = −ci,−γcprimei = −δijcprimej,andγ ci||c||(summationtextj cprimej) =˜Fijcprimej where ˜Fij = γ ci||c||.Thus, if we define MF = D+F withDij = −ciandFij = γ( ci||c|| −δij)thenMF(c)cprime = epsilon1n. (3.16)Now we expand c = c(0) +epsilon1c(1) +O(epsilon12) and substitute this into (3.16) toget MF|c(0)+epsilon1c(1)+O(epsilon12)(c(0) + epsilon1c(1) + O(epsilon12)) = epsilon1n which yields an equationat each order of epsilon:O(1) : MF|c(0)c(0)prime = 0 (3.17a)O(epsilon1) : MF|c(0)c(1)prime + (∂MF∂ci|c∗c(1)i )c(0)prime = n. (3.17b)As MF(c(0)) is invertible we obtain that c(0)prime = 0 so c(0) is a constantvector. Given the boundary conditions here that c(0) = c∗ we have c(0) = c∗.This makes (3.17b) much simpler because all the terms with c(0)prime arezero (as shown by Promislow et al. [10]). So the order epsilon1 equation is simplyMF|c∗c(1)prime = n which can be solved to find c(1)prime :1093.3. Asymptotic Analysisc(1)prime = MF|−1c∗ nwhich is constant.As the boundary conditions are already satisfied with c(0) we need thatc(1)(0) = 0 and the first-order approximation to the concentration is:c(x) = c∗ +epsilon1c(1)primex+O(epsilon12). (3.18)Substituting x = 1 into (3.18), we obtain an approximation to the con-centrations at the catalyst site. The most interesting information for us incomputing the difference between the two diffusion models is the relativechange in concentrations from channel values. So we compute ri = ci(1)−c∗ic∗ifor the three species.The computations yieldr1 = −0.0203r2 = 0.0617r3 = −0.00325.The signs of these quantities are generally what we would expect. As wa-ter vapor (component 2) is being produced on the catalyst layer, its concen-tration should be higher than in the channel. Oxygen, the first component,which is reacting on the catalyst layer to produce water would decrease inconcentration as we move along the GDL from the channel. Interestingly1103.3. Asymptotic Analysisenough, although Nitrogen does not react, because of the diffusion and massconservation, it still has a concentration gradient.3.3.2 Asymptotic Analysis of Maxwell-StefanThe analysis here is nearly identical to Fick. Instead of γ we have γQij.By expanding the derivatives in (3.7) we can write:MMScprime = epsilon1n (3.19)with MMS = D+G whereDij = −ciandGij = γQik(− ck||c||2 + δkj||c||).By expanding c asymptotically with c = c(0) +epsilon1c(1) +O(epsilon12) the identicalanalysis holds as in the equations for Fickian diffusion.We again have that c(0) is a constant equalling c∗ and c(1)prime = MMS|c∗−1nso that our first-order solution is identical to (3.18) but with the differentvalue for c(1)prime.We again compute the ri = ci(1)−c∗ic∗ivalues and find this time:r1 = −0.0140r2 = 0.04061113.4. Numerical Analysis of Diffusion Modelsr3 = −0.00196The signs are again what we expect, but the results are quite differentfrom Fick’s law. The magnitude of these relative changes is smaller than inthe Fickian model.3.4 Numerical Analysis of Diffusion ModelsHere we present some numerical computations used to verify the validity ofthe asymptotic results.3.4.1 Discretizing the Diffusion EquationsWe use (3.16) and (3.19) to express the derivative at each point as the inverseof a matrix (dependent upon the local concentration) times the fluxes.Given N, the number of grid points, we define h = 1/N and set xi =(i−1)h for i = 1,...,N +1. We define ci to be the approximation to c at xi.We now use a modified Euler stepping scheme [15] beginning with c1 =c∗.To get ci+1 from ci we define:k1 = epsilon1M|−1ci nk2 = epsilon1M|−1ci+hk1nandci+1 = ci + h2(k1 +k2).1123.4. Numerical Analysis of Diffusion ModelsFigure 3.3: Verifying second-order convergence with the Fick code.The approximate concentration vector at x = 1 is cN+1. From here wecan also approximate ri = cN+1i−c1ic1i .3.4.2 Verification of Program ResultsThe modified Euler method is second-order accurate. To test the approxi-mation, we compute the ri values at N = 10,20,40,80 and N = 200.Given the simplicity of the problem, we take the result at N = 200 asthe “exact” answer and compute the maximum magnitude of the error foreach component of r for N = 10,20,40,80. We can plot the error trend ona log-log plot to verify the second-order convergence. The plot is given forthe Fickian code in figure 3.3.We can also test the programs against the asymptotic solutions. Running1133.5. Exploring the Fundamental Differences between the Modelsthe program (at N = 200) in the Fickian and Maxwell-Stefan models yieldsconcentrations that are identical to the asymptotically computed solutionsto within 1.00 × 10−4 and 3.29 × 10−5 respectively. We now have a highdegree of confidence in both the asymptotic work and the coding.3.5 Exploring the Fundamental Differencesbetween the ModelsThere are numerous questions that arise from the analysis and we will nowcarry out some investigations.In speaking with fuel cell engineers, modern fuel cells actually have apermeability of 10−15 m2 [8]. This would make the Darcy’s law less efficientand it could impact the way the models behave.There is a question of what happens to the Maxwell-Stefan predictionswhen we use the “wrong” fluxes (the fluxes that Fick’s law inevitably uses).We can redo the calculations for the Maxwell-Stefan model by replacing J,the molar flux with respect to the mass-averaged velocity, by J∗, the molarflux with respect to the molar-averaged velocity. In this regime, Q (as in(3.15)) is simply (PA+B)−1P.We also remark that the diffusivity D as used by Promislow et al. [10] isnoticeably smaller than the binary diffusivities used in the Maxwell-Stefanformulation. We wish to see what happens if we replace D by the arithmeticmean of the binary diffusivities, 1.17×10−5 m2 s−1 and compute the Fickianrelative changes.Our results are tabulated in table 3.4 and they were achieved by means1143.5. Exploring the Fundamental Differences between the ModelsTable 3.4: Different modeling predictions for the relative changes in theconcentrations.i (M,κ) (M,κ) ( ˜M,κ) ( ˜M,κ) (F,κ) (F,κ) (˜F,κ) (˜F,κ)1 −0.0140 −0.0142 −0.0143 −0.0137 −0.0203 −0.0188 −0.0139 −0.01242 0.0406 0.0404 0.0403 0.0410 0.0617 0.0632 0.0422 0.04373 −0.00196 −0.00211 −0.00185 −0.00118 −0.00325 −0.00177 −0.00222 −0.000743of the numerical programs in section 3.4. Notationally, (X,Y) is used whereX can be M (the standard Maxwell-Stefan formulation), ˜M (the Maxwell-Stefan with the incorrect flux), F (Fick’s law with D as initially stated), ˜F(Fick’s law with the updated value of D); and Y can be κ (the permeabil-ity initially stated) or κ (the smaller permeability that is possible in moremodern fuel cells). The i indicates relative change ri.From table 3.4 it is very clear that by choosing the updated value ofD, the Fickian model is in better agreement with Maxwell-Stefan. With asmaller permeability, the updated Fickian model fares worse. This is likelydue to the fact that Darcy’s law is weaker and simply updating the valueof D isn’t good enough because we start to see the different interactionscontained in the Maxwell-Stefan law.We also note that the changes in using the wrong flux in the Maxwell-Stefan setting are negligible.We include the plots (figures 3.4 to 3.6) of the different concentrationprofiles for the gas species in both the Fick and Maxwell-Stefan settings,where we have chosen an improved value of the diffusivity and the smallervalue of the permeability.We see in both cases, the plots are nearly linear and we see the asymp-totics are a good match.1153.6. Conclusions of Gas Diffusion ModelingFigure 3.4: The concentration profile of Oxygen in the GDL for both gasdiffusion models.3.6 Conclusions of Gas Diffusion Modeling3.6.1 FormulationsMaxwell-Stefan diffusion in its purest form, is a highly nonlinear system ofequations with many intricacies. Fick diffusion is nonlinear, but has fewerdifficulties involved in its computation.Fick’s law has formulation inconsistencies when coupled with Darcy’slaw, and it does not take into account the individual competitions betweendiffusing species. Instead, it lumps all the binary diffusivities together into asingle constant (which can agree quite well with Maxwell-Stefan if properlyselected).Although they differ, the dominant driving force in both models is the1163.6. Conclusions of Gas Diffusion ModelingFigure 3.5: The concentration profile of water vapor in the GDL for bothgas diffusion models.same: the bulk diffusion of gas species given in Darcy’s law in porous media.The system with Darcy’s law alone would not be solvable; both Fickian andMaxwell-Stefan diffusion are perturbations to a singular system.3.6.2 Quantitative DifferencesWe saw that in a very simple setting with one dimension and an isothermal,steady state set of conditions that there are very small differences betweenthe two diffusion models when the Fick diffusivity constant is judiciouslyselected and for a permeability 10−12 m2. For the smaller permeability of10−15 m2 the differences are no longer so small (although much smaller thanif the diffusivity had not been modified).It’s important to note that the differences we obtained could be larger1173.6. Conclusions of Gas Diffusion ModelingFigure 3.6: The concentration profile of Nitrogen in the GDL for both gasdiffusion models.in more realistic settings. Our one dimensional model has many limitations,one being that it doesn’t take into account the difficult diffusion pathwaysof the gas species. The effective length each species must diffuse could be asmany as four or five times larger than what our computations were basedupon. It’s also possible the differences could be more significant in a higherdimensional setting or with more physically realistic conditions (temperaturevariations, etc.).118Chapter 4Summary and Future Work4.1 SummaryHere we provide a summary of the research results and discuss future ex-tensions to the work done.4.1.1 Superconductor Research SummaryIn superconductor systems governed by the London equation, many physi-cal properties are well-understood with certain assumptions on the surfacegeometry of the superconductor. In particular, if a magnetic field is ap-plied parallel to a superconducting interface, and the interface is flat, thenthe field magnitude should decay exponentially with the distance into thesuperconductor (when in the Meissner state).Recent experiments have measured magnetic field profiles that deviatefrom this exponential decay. Our research project in chapter 2 explored apossible explanation for the deviations: the notion of surface roughness.To reach our conclusions, careful asymptotics were done that accuratelydescribe field profiles in superconductors with sinusoidal surface roughness.Novel numerical methods with a carefully chosen mesh were used to verifythe asymptotic results and then provide results beyond what the asymptotics1194.1. Summaryallow.Two new phenomena are predicted by the study. The first is that theindividual field components (and even the field magnitude), on particularregions on the vacuum-superconducting interface can their exceed valuesgiven by the applied field (experimentally this cannot be detected becausethe experiments only measure the average field). The second is that if thereis a rough interface, then in general the field in the vacuum region would alsobe perturbed - and none of the components will be identically zero either inthe vacuum or the superconductor.4.1.2 Gas Diffusion Research SummaryIn the modeling of gas transport, there exist a number of diffusion models.Fick diffusion is often implemented in modeling due to its relative simplic-ity. However, another diffusion model known as Maxwell-Stefan diffusion isbelieved to be more accurate due to various interspecies interaction termsthat appear in its formulation. This model is a great deal more complicated.We sought to determine the qualitative and quantitative differences be-tween the two models, to find out if in fact there is a significant differencewhen computations are done with one model or the other in situations rel-evant to fuel cell operation.Using a simple, one-dimensional model of a gas diffusion layer in a fuelcell, we computed the relative changes in concentration for three gas specieswith both sets of equations. Our work here (chapter 3) was carried out inan asymptotic regime and verified with numerical programs. The numericalprograms were also used to further study the predictions of the models.1204.2. Future WorkWe showed that the differences between relative concentration changespredicted by the two models are negligible if the Fickian diffusivity is chosencorrectly and if the permeability is not too small. The models do, however,show small differences in their predictions for smaller (more modern) per-meability values.We saw that both models are dominated by the bulk diffusion of theentire system and not the diffusion of individual species (results of Stockieet al. [9] are in agreement).4.2 Future Work4.2.1 Future Work for Superconductor ProjectThrough various validations, we have seen that the asymptotic and numer-ical computations are in agreement and agree with physical intuition. Evenwith a small number of grid points, the three-dimensional code we wroteshows a high level of accuracy.However, in the future if this code needed to be even more accurate, wewould need to maximize the number of grid points in the three-dimensionalfinite difference code. However, the program and environment both poselimitations.We started with a Matlab implementation to test the formulation, butrewriting the code in a more computationally friendly language such as Cwould be a next step. Matlab doesn’t have enough memory or speed todeal with systems of the size we desire. In addition, we need the compu-tation to be as fast as possible. By making use of the known solution for1214.2. Future Worka flat-interface, we could use Fast Fourier Transformations to come up anefficient preconditioner for the matrix equations which could then be solvediteratively using a Krylov subspace method such as GMRES [16].Exploring symmetries could also help reduce the number of unknowns.Given the surfaces we examined were sinusoidal, there are be planes of sym-metry for the solution. This could reduce the number of grid points neededby up to a factor of 4.Analytically, another consideration is the asymptotic analysis. One phys-ical limit that was inaccessible in the asymptotics was if the spatial fre-quencies became very large. To analyse these limits, it’s possible that ho-mogenization theory holds the answers. Research has been done in solvingMaxwell’s equations with very rough boundaries [17]. This could be possiblein our case as well.4.2.2 Future Work for Gas Diffusion ProjectHaving seen there is such a small differences between Fick Diffusion andMaxwell-Stefan diffusion in our model, some serious questions arise. At whatpoint do the two models differ significantly (to the point where Maxwell-Stefan diffusion could be necessary for modeling)? Also, what is the bestchoice of the Fickian diffusivity when used in conjunction with Darcy’s law,given the set of binary diffusivities?A next course of action would be to model different systems of gases anddetermine where the two models differ, and how to select the diffusivity.It would also be nice to consider more realistic conditions. Not onlywould we ideally take the model here to a higher dimensional setting, but1224.2. Future Workwe would need to add temperature variations, as well as the presence ofliquid water, to the calculations.123Bibliography[1] Morton, K. and Mayers, D., Numerical Solutions of Partial DifferentialEquations. Cambridge University Press: 1994.[2] Bender, C. and Orzag, A., Advanced Mathematical Methods for Sci-entists and Engineers: Asymptotic Methods and Perturbation Theory,Springer-Verlag New York, Inc.: 1999.[3] Folwer, A., Mathematical Models in the Applied Sciences, CambridgeTexts in Applied Mathematics: 1997.[4] Ashcroft, N. and Mermin, N., Solid State Physics, Saunders CollegePublishing: 1976.[5] Sonier, J., “Muon Spin Rotation/Relaxation/Resonance (µSR)”(http://musr.org/ jess/musr/muSRBrochure.pdf).[6] Kiefl, R.; Hossain, M.; Wojek, B.; Dunsiger, S.; Morris, G.; Prokscha,T.; Salman, Z.; Baglo, J.; Bonn, D.; Liang, R.; Hardy, W.; Suter,A.; and Morenzoli, E., “Direct Measurement of the London Pen-etration Depth in YBa2Cu3O6.92 Using Low-Energy µSR”, Phys-ical Review B 2010, 81:18, 180502. Abstract can be viewed at:http://link.aps.org/doi/10.1103/PhysRevB.81.180502124Bibliography[7] London, F. and London, H., “The Electromagnetic Equations of theSupraconductor”, Series A, Mathematical and Physical Sciences 1935,149:886, 71-88.[8] Private communication[9] Stockie, J.; Promislow, K.; and Wetton, B., “A Finite Volume Methodfor Multicomponent Gas Transport in a Porous Fuel Cell Electrode”,International Journal for Numerical Methods in Fluids 2003, 41, 577-599.[10] Promislow, K.; Stockie, J.; and Wetton, B., “A Sharp Interface Re-duction for Multiphase Transport in a Porous Fuel Cell Electrode”,Proceedings of the Royal Society A: Mathematical, Physical and Engi-neering Science 2006, 462, 789-186.[11] Berg, P.; Promislow, K.; St. Pierre, J.; Stumper, J.; and Wetton, B.,“Water Management in PEM Fuel Cells”, Journal of the Electrochem-ical Society 2004, 151:3, A341-A353.[12] CRC Handbook of Chemistry and Physics, 90th ed., CRC Press: 2009.[13] Taylor, R. and Krishna, R., Multicomponent Mass Transfer, Wiley Se-ries in Chemical Engineering: 1993.[14] Bear, J. and Bachmat, Y., Introduction to Modelling Transport Phe-nomena in Porous Media, Kluwer Academic: 1990.[15] Hamming, R., Numerical Methods for Scientists and Engineers, GeneralPublishing Company: 1973.125[16] Quarteroni, A. and Sacco, R.; and Saler, F., Numerical Mathematics,Springer-Verlag, Inc.: 2000.[17] Nevard, J. and Keller, J., “Homogenization of Rough Boundaries andInterfaces”, SIAM Journal on Applied Mathematics 1997, 57:6, 1660-1686.126Appendix AEigenvalues of FiniteDifference MatrixWe seek the eigenvalues of the (n+ 1)×(n+ 1) matrix M in section 1.1.1.The matrix is clearly invertible (by simple row reduction it can be seen tohave n+ 1 pivots).If Mu = λu then u1 = λu1, un+1 = λun+1 and uj+1−2uj+uj−1 = h2λujfor 2 ≤ j ≤ n−1.The endpoints tell us that u1 = un+1 = 0 since λ negationslash= 0 (M is invertible).On the interior we will look for a solution of the form uj = θj. We theneasily findθj−1(θ2 −βθ + 1) = 0 (A.1)where β = (2 +h2λ).If β2 − 4 ≥ 0 then uj = Aθj+ + Bθj− (where θ± satisfies the quadraticequation for β2 − 4 > 0) or uj = Aj + B (when β2 − 4 = 0) neither ofwhich can lead to anything other than the zero vector with the restrictionsu1 = un+1 = 0.On the other hand, if β2−4 < 0 then we have two roots that are complexconjugates, each of which has modulus 1 (since their product must be 1).127Appendix A. Eigenvalues of Finite Difference MatrixWe denote these roots as exp(±iα). We require that Aeiα + Be−iα = 0so that A = −Be−2iα. Also we need Aei(n+1)α + Be−i(n+1)α = 0 so A =−Be−i(2n+2)α.By dividing the two equations for A we see this is only possible if e2niα =1 so that 2niα = 2pii. Thus α = pikn , k = 1,2,...,n−1. The constant α cannotbe zero (or n) because we require two distinct complex roots.Solving the quadratic for θ in (A.1) and taking the real part, we findθ = β/2. The real part of eikpi/n is cos(kpi/n). This allows us to find theeigenvalues:β/2 = 1 +h2λk/2 = cos(kpi/n)which meansλk = 2cos(kpi/n)−2h2 = 2n2(1− 12pi2k2n2 +O(n−4)−1) = −pi2k2 +O(n−2).The eigenvalue of smallest magnitude occurs at k = 1 and is given by−pi2 as n →∞. Therefore the matrix M−1 is has a bounded sup-norm.128Appendix BThe General InterfaceHere we show how to asymptotically compute the solution for a generalinterface z = epsilon1h(x,y) with the assumption that h and its powers h2,h3,...can be Fourier transformed. We will use the symmetric Fourier transformthroughout.For all surfaces, the zeroth order term b(0) will be the base solution givenin (2.8). We can write the interfacez = 1√2piintegraldisplayR2epsilon1ˆh(ωx,ωy)eiωxx+iωyydωxdωy,where ˆh is the Fourier transform of h. Then by (2.16b) we have[b(1)]|z=0 = 1√2piintegraldisplayR2ˆh(ωx,ωy)ei(ωxx+ωyy)dωxdωy.We could also go to the equations at higher powers of epsilon1.129Appendix B. The General InterfaceBy modifying (2.17) for its Fourier integral analogue, we can write:b(j) = 1√2piintegraldisplayR2dωxdωyei(ωxx+ωyy)γ(1)pi1β(j)eωz if z ≤ 0(γ(2)pi2 +γ(3)pi3)β(j)e−ωz if z > 0(B.1)whereω(ωx,ωy) =radicalBigω2x +ω2y, ω(ωx,ωy) =radicalBig1 +ω2x +ω2y, γ(1)ωx,ωy = (iωx,iωy,ω)T,γ(2)ωx,ωy = (ω,0,iωx)T, γ(3)ωx,ωy = (0,ω,iωy)T, and β(j)ωx,ωy = M−1ωx,ωy hatwider[b(j)]|z=0.We can thus find the perturbation terms at every order.For an easy example we will consider a surface with roughness in onlyone direction, z = epsilon1 sinc y, with applied field (1,0,0)T and we will computethe first-order perturbation. Recall sinc y = (siny)/y.The symmetric Fourier transform of sinc y is radicalbigpi2χ[−1,1](ω), where χ isthe characteristic function. Thereforehatwider[b(1)]|z=0 = −ˆh[∂zb(0)]z=0 = (radicalbiggpi2χ[−1,1],0,0)T.By carrying out the computation, we findγ(j)ωx,ωypijβ(1)ωx,ωy =radicalbigpi2χ[−1,1]δ2,j.Then (by (B.1)) we conclude, for b(1), all components in the vacuum arezero and all components except for b(1)1 are zero in the supercondor. We findb(1)1 = 1√2piintegraldisplay ∞−∞radicalbiggpi2χ[−1,1]eiωye−√1+ω2zdω = 12integraldisplay 1−1e−√1+ω2zeiωydω.If we expand the complex exponential as cos(ωy) + isin(ωy) and take notethat the imaginary part of the integrand is an odd function of ω integrated130Appendix B. The General Interfaceover a symmetric range (and hence zero) then we obtainb(1)1 = 12integraldisplay 1−1e−√1+ω2zcos(ωy)dω =integraldisplay 10e−√1+ω2zcos(ωy)dω.To first-order (neglecting b2 and b3 which remain identically zero) wecan write:b1 =1 if z ≤ epsilon1 sinc ye−z +epsilon1integraltext10 e−√1+ω2zcos(ωy)dω if z > epsilon1 sinc yThis process could go on to obtain higher-order approximations.131Appendix CDead Layers and AveragesThe terms of importance in taking the averages are the constants (withrespect to x and y) and those involving the squares of sine and cosine (oth-erwise they would average to zero).Let us suppose a term had a factor of cos2(ωxx)cos2(ωyy) present in itin the third geometry (i.e. neither ωx nor ωy is zero). Then, taking theaverage of this over x and y is trivially 14.However, if ωx were zero then we would be computing the average valueof cos2(ωyy) which is actually 12.One of the issues here is that cos2(ωxx) does not approach cos2(0) = 1uniformly as ωx → 0.132Appendix DProof of Periodicity of ˜gWe wish to prove that if b : R3 → R3 (where b = b(x,y,z)) is periodic inx and y with periods Tx and Ty and it is defined by b = ∇g for some C∞scalar function g satisfying g(x,y,−∞) = x, then the function ˜g = g−x isperiodic in x and y satisfying the same periodicity conditions as b.The boundary condition at z = −∞ is clearly ˜g(x,y,−∞) = 0.As b is periodic (and C∞) we have the following equations for any x0,y0 and z0 fixed:integraldisplay x0+Txx0∂xb(x,y0,z0)dx =integraldisplay y0+Txy0∂yb(x0,y,z0)dy = 0.Replacing b with ∇˜g +(1,0,0)T and looking at the third component weobtain:integraldisplay x0+Txx0∂zx˜g(x,y0,z0)dx =integraldisplay y0+Tyy0∂zy˜g(x0,y,z0)dy = 0. (D.1)To show that ˜g is periodic, it would suffice to prove that bothintegraltextx0+Txx0 ∂x˜g(x,y0,z0)dxand integraltexty0+Tyy0 ∂y˜g(x0,y,z0)dy = 0.We prove the first is zero, as the second is done by the same methodology.133Appendix D. Proof of Periodicity of ˜gWe writeintegraldisplay x0+Txx0∂x˜g(x,y0,z0)dx =integraldisplay x0+Txx0(∂x˜g(x,y0,−∞)+integraldisplay z0−∞∂xz˜g(x,y0,z)dz)dxand by Fubini’s theorem (to re-order the integrations) and Claurant’s theo-rem (in the equality of the mixed partial derivatives) along with noting that˜g(x,y0,−∞) = 0 we findintegraldisplay x0+Txx0∂x˜g(x,y0,z0)dx =integraldisplay z0−∞integraldisplay x0+Txx0∂zx˜g(x,y0,z)dxdz = 0where the last equality comes from the fact that the inner integral is zeroby equation D.1.We have shown that ˜g is periodic with the same periodicity as b.134Appendix EProof of Existence of a NullVectorWe show here that for even N, the three dimensional code matrix (in theflat geometry) admits a null vector, making the system unsolvable. We willuse the same notation as in section 2.4.1.Here we consider an unstretched coordinate system with x,y and z as thevariables. The transformed coordinate system would still be approximatingthe same system (so if there is a null vector in Cartesian coordinates, thereis also a null vector in the stretched coordinates to within an error imposedby the numerical system).We wish to find a vector u so that Mu = 0. Such an equation would imply(in a discretized sense) the following: firstly, ˜g(x,y,−∞) = 0; secondly,triangle˜g = 0; thirdly, ∇˜g−b = 0 along the interface; fourthly, ∇•b = 0 along theinterface on the superconducting side; fifthly, triangleb = b in the superconductingregion; and sixthly, b1(x,y,∞) = b2(x,y,∞) = ∂zb3(x,y,∞) = 0.On the vacuum side, we begin by finding a solution vector with zeroLaplacian (in the discretized setting) of form uK(i,j,k,0) = (−1)i+jλk fork = 2,...,Nv −1 with a uniform mesh in all directions. Here, a point on the135Appendix E. Proof of Existence of a Null Vectormesh can be described by (xi,yj,zk) with spacing h.Computing the discretized Laplacian and setting it equal to zero we haveui+1,j,k +ui−1,j,k +ui,j+1,k +ui,j−1,k +ui,j,k+1 +ui,j,k−1 −6ui,j,kh2 = 0.Thus,(−1)i+jλk−1(λ2 −6λ + 1) = 0, which has roots in reciprocal pairs, λ =3±√32/2.If α is small (close to machine epsilon) then α(−1)i+j(3 + √32/2)k isnearly zero at k = 1 and grows exponentially in magnitude with k.Because of how we discretized∇˜g−b by taking averages, at the interface,the discretized ∇˜g is zero along the interface making all components of bzero along the interface. Then if the discretized b is zero everywhere on thesuperconducting side then the remaining conditions are upheld.Note the periodicity implies that (−1)1+j = (−1)N+1+j for fixed j (andsimilarly switching i and j) and this could only happen if N is even. There-fore, uK(i,j,k,lscript) = α(−1)i+jλkδlscript0 is a null vector for even values of N.136
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Asymptotic and numerical modeling of magnetic field...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Asymptotic and numerical modeling of magnetic field profiles in superconductors with rough boundaries.. Lindstrom, Michael Robert 2010-08-31
pdf
Page Metadata
Item Metadata
Title | Asymptotic and numerical modeling of magnetic field profiles in superconductors with rough boundaries and multi-component gas transport in PEM fuel cells |
Creator |
Lindstrom, Michael Robert |
Publisher | University of British Columbia |
Date Issued | 2010 |
Description | This thesis is a combination of two research projects in applied mathematics, which use the applied math techniques of numerical and asymptotic analysis to study real-world problems. The first problem is in superconductivity. This section is motivated by recent experimental results at the Paul Sherrer Institute. Here, we need to determine how the surface roughness of a superconductor influences the penetration properties of an externally applied magnetic field. We apply asymptotic analysis to study the influences, and then verify the accuracy - even going well-beyond the limits of the asymptotics - by means of computational approximations. Through our analysis, we are able to offer insights into the experimental results, and we discover the influence of a few particular surface geometries. The second problem is in gas diffusion. The application for this study is in fuel cells. We compare two gas diffusion models in a particular fuel cell component, the gas diffusion layer, which allows transport of reactant gases from channels to reaction sites. These two models have very different formulations and we explore the question of how they differ qualitatively in computing concentration changes of gas species. We make use of asymptotic analysis, but also use computational methods to verify the asymptotics and to study the models more deeply. Our work leads us to a deeper understanding of the two models, both in how they differ and what similarities they share. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-08-31 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NoDerivs 3.0 Unported |
DOI | 10.14288/1.0071198 |
URI | http://hdl.handle.net/2429/28000 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2010-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nd/3.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2010_fall_lindstrom_michael.pdf [ 1.88MB ]
- Metadata
- JSON: 24-1.0071198.json
- JSON-LD: 24-1.0071198-ld.json
- RDF/XML (Pretty): 24-1.0071198-rdf.xml
- RDF/JSON: 24-1.0071198-rdf.json
- Turtle: 24-1.0071198-turtle.txt
- N-Triples: 24-1.0071198-rdf-ntriples.txt
- Original Record: 24-1.0071198-source.json
- Full Text
- 24-1.0071198-fulltext.txt
- Citation
- 24-1.0071198.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0071198/manifest