Table Of Content

Optimal design of a model energy conversion device Lincoln Collins and Kaushik Bhattacharya Division of Engineering and Applied Science California Institute of Technology Pasadena, CA 91125 7 Email: [email protected] and [email protected] 1 0 February 2, 2017 2 n a J Abstract 1 3 Fuel cells, batteries, thermochemical and other energy conversion devices involve the transport of a number of (electro-)chemical species through distinct materials so that they can meet ] and react at specified multi-material interfaces. Therefore, morphology or arrangement of these P differentmaterialscanbecriticalintheperformanceofanenergyconversiondevice. Inthispa- A per, we study a model problem motivated by a solar-driven thermochemical conversion device . h that splits water into hydrogen and oxygen. We formulate the problem as a system of cou- at pled multi-material reaction-diffusion equations where each species diffuses selectively through m a given material and where the reaction occurs at multi-material interfaces. We express the problem of optimal design of the material arrangement as a saddle point problem and obtain [ an effective functional which shows that regions with very fine phase mixtures of the material 1 arise naturally. To explore this further, we introduce a phase-field formulation of the optimal v design problem, and numerically study selected examples. 4 2 0 1 Introduction 0 0 . 2 The efficiency of fuel cells, batteries and thermochemical energy conversion devices depends on 0 inherent material characteristics that govern the complex chemistry and transport of multiple 7 species as well as the spatial arrangement of the various materials. Therefore, optimization of the 1 : spatial arrangement is a recurrent theme in energy conversion devices. Traditional methods of v i synthesis offer limited control of the microstructure and there has been much work in advanced X imaging for these uncontrolled microstructures (e.g., [35]) and optimizing gross features. However, r a the growing ability for directed synthesis [31, 23, 33, 26] allows us to ask the question of what microgeometriesareoptimalforparticularapplications. Inthissensewedirecttheproblemtooneof optimaldesignwherewearenotlimitedbytheimaginationindeterminingnewmicrostructuresbut insteadallowfortheunderlyingphysicalbehaviorandoptimizationtechniquestodirectarchitecture and microstructure, and eventually lead synthesis to unprecedented performance. The tailoring of material microstructure and nanostructure is not new to energy conversion and storage [26, 4]. The development of hierarchical structures and porosity affords balancing interfacial reactions and chemical transport to maximize efficiency. For example, in cathodes of solidoxidefuelcells, theefficiencyislargelydeterminedbyreactionsattriplephaseboundariesand the availability of transport pathways through each phase. These features are coupled with volume fractions, surface area densities, interfacial curvatures, and phase tortuosities to find the optimal balancebetweensurfacereactionandtransport[34,28]. Similarly,theimportanceofmicrostructure 1 onanodeperformance[5,30,15,13]hasalsobeenestablished. Themorphologyofmaterialsusedin lithium ion batteries is of interest from both the theoretical [29, 10, 9] and experimental standpoint [21, 22, 33]. Mass and ion transport and interface measure in battery electrodes directly impact the storage capacity and rate performance and is an ideal problem for optimization across many length scales. The application of metal oxides for solar-driven thermochemical conversion devices offers a promising new sustainable energy source [11]. Here, a porous, redox active oxide is cyclically exposedtoinertgasathightemperature,generatingoxygenvacanciesinthestructure,andreactant gas (H O, CO ), at moderate temperature, releasing fuel upon reoxidation the oxide [31]. The lack 2 2 of complex and expensive catalyst systems and full use of the entire solar spectrum separate these devices from many other photo-based energy sources. Recent advances made in the materials research community indicate many possible candidates for these applications, and lend themselves to advanced synthesis techniques facilitating directed architecture, where significant improvements can be made [32]. The thermodynamic and kinetic behavior of these materials are well-studied [25, 12, 20], fully describing the gas phase transport of reactant and product gases coupled to the solid state mass and charge transport occurring through the bulk. In this paper, we study a model system motivated by metal oxides in solar-driven thermochemical conversion devices. We have a two phase material (solid oxide and pore) where reactions at the surface create (gaseous) oxygen in the carrier gas in the pores and bound oxygen in the solid oxide; the oxygen diffuses through the carrier gas in the porous region and bound oxygen diffuses through the solid oxide. We seek to understand the arrangement of the solid and porous regions to maximize the transport given sources and sinks for the gaseous oxygen and vacancies. Thereisalargeliteratureinthestudyofoptimaldesignproblems,especiallyseekingtominimize compliance for a given weight as well as maximize conduction for a given mass. It is understood that the underlying problem is ill-posed in that the optimal designs often lie outside of the set of “classical admissible designs” and one has to either relax the problem by homogenization [19, 24] or regularize it by the introduction of perimeter constraints [3, 8]. It leads to two widely used methods, topology optimization (e.g. [7]) and shape optimization (e.g. [1]). The presence of two species lends a vectorial character to our problem, and the presence of the surface sources makes the problem at hand different from those in the literature. We begin with the formulation of physical problem, shown in Figure 2.1, in Section 2. We start with a sharp interface formulation. However, the optimal design of the sharp interface model is mathematically ill-posed, and therefore we study the analogous diffuse interface model. We also note that homogenization of the sharp interface model leads to equations of the same form as the diffuse interface equations. The transport of two chemical species with concentration u and u is 1 2 governed by the following reaction diffusion equations for i = 1,2:  k u = f , in Ω,  i i i ∇· ∇ k u nˆ = 0 on ∂Ω ∂ Ω (1) i i i  ∇ · \ u = u∗ on ∂ Ω, i i i where the isotropic conductivities are k (x) = k χ(x)+k (1 χ(x)), k (x) = k χ(x)+k (1 χ(x)); (2) 1 11 12 2 21 22 − − with k ,k >> k ,k > 0, and the sources are 11 22 12 21 f = f = χ(1 χ)k (u u ) (3) 1 2 s 1 2 − − − 2 @ ⌦ 2 ⌦ nˆ 2 S mˆ ⌦ 1 @ ⌦ 0 @ ⌦ 1 Figure 1: The physical setting: chemical species 1 enters through the source ∂ Ω, diffuses through 1 Ω , isconvertedtochemicalspecies2throughasurfacereactionattheinterface , chemicalspecies 1 S 2 diffuses through Ω and leaves through the sink ∂ Ω. 2 2 for k > 0 and χ : Ω [0,1]. Briefly, we have a two-material system and χ describes the s → volume fraction of material 1 (say solid phase). Chemical species 1 (say bound oxygen) diffuses preferentially in material 1 (χ = 1) while species 2 (say oxygen gas) diffuses preferentially in material 2 (say pore, χ = 0). The species react and therefore there is a source at the interface χ = 0,1. For future use, we write the source as f = χ(1 χ)Au where f = f ,f , u = u ,u 1 2 1 2 (cid:54) − { } { } and (cid:18) (cid:19) 1 1 A = ks − . 1 1 − This problem allows a variational formulation, and the direct method of the calculus of variation allows us to prove existence of a solution. We study the optimal design problem of maximizing the flux of species through the reactor over allpossiblearrangementsχinSection3. Weshowthatthisgivesrisetoasaddlepointproblem. We then obtain an explicit characterization which shows that the mixed phase regions arise naturally. To understand this further through particular examples, we introduce a phase field formulation in Section 5. Specifically, we add an Allen-Cahn type energy to that associated with the variational formulation of the transport problem and then solve the gradient flow associated with this energy. We solve this numerically in selected examples and conduct a parameter study. These show that the optimal design can be quite intricate as it seeks to balance transport and reaction. 2 Formulation 2.1 Sharp interface formulation Consider an open, bounded region Ω Rn with Lipschitz boundary separated into two regions Ω 1 ⊂ and Ω by an interface show in Figure 2.1. We consider the diffusion of one species with concen- 2 S tration u in region Ω with isotropic diffusivity K > 0, and a second species with concentration 1 1 1 u in region Ω with isotropic diffusivity K > 0. The two species meet at the interface and react 2 2 2 with reaction rate k > 0. The boundary of Ω is divided into three regions ∂Ω = ∂ Ω ∂ Ω ∂ Ω s 1 2 0 ∪ ∪ where ∂ Ω ∂Ω . The concentration of species i is held at a prescribed value u∗ on ∂ Ω while ∂ Ω i ⊂ i i i 0 3 is insulating. This is described by the following system of equations:  ∇·Ki∇ui = 0 in Ωi   K u mˆ = k (u u ) on i i s 1 2 − ∇ · − S (4) ui = u∗i on ∂iΩ   K u nˆ = 0 on ∂Ω ∂ Ω i i i ∇ · \ for i = 1,2 where mˆ represents the normal to pointing from Ω pointing to Ω , and nˆ represents 1 2 S the outward normal to ∂Ω. 2.2 Diffuse interface formulation It is often convenient to work with a smooth or diffuse interface formulation of the problem above. We now show formally that the diffuse interface formulation in (1)-(3) leads to the sharp interface formulation in (4) in an asymptotic limit. Let χ be the characteristic function of Ω as defined in 1 Section 2.1. Let χη be the mollification of χ with a standard mollifier at length-scale η: χη = ϕη χ where ϕη(x) = η−nϕ(x/η). Let kη,kη be as in (2) with χ = χη and kη = (1 exp( η))k∗ , η η 1 2 η 12 − − 12 k = (1 exp( η))k , and f (x) = f (x/η). Let u solve 21 − − 21 i i i kη uη = η2fη, i = 1,2 in Ω. (5) ∇· i∇ i i First consider the outer expansion η 0, and note that (5) formally gives (4) in Ω . Further, 1,3,4 1,2 → note that u (respectively u ) is indeterminate on Ω (respectively Ω ). However, this outer expan- 1 2 2 1 sion does not give any condition on the interface . To obtain this condition, denote the limiting S values on the interface to be u¯ ,u¯ . We seek to relate these to the flux as in (4) . 1 2 2 Now consider the inner expansion. Pick a point x and change variables x (x x )/η. 0 0 ∈ S (cid:55)→ − We obtain k u = f , i = 1,2 (6) i i i ∇· ∇ where k = k1. Further, as η 0, χ and hence the solution depend only on one dimension that is i i → normal to the interface. We take this direction to be x by changing variables if necessary. Let U 1 i solve (6) for the boundary conditions (u ,u )(x ) (1,0) as x and (u ,u )(x ) (0,1) 1 2 1 1 1 2 1 → → −∞ → as x , and V solve (6) for the boundary conditions (u ,u )(x ) (1,0) as x and 1 i 1 2 1 1 → ∞ → → −∞ (u ,u )(x ) (0, 1) as x . Note that 1 2 1 1 → − → ∞ u = αU +βV +γ i i i also solves (6) for any arbitrary scalars α,β,γ, and satisfies the boundary conditions u α+β+γ as x , u α β+γ as x . 1 1 2 1 → → −∞ → − → ∞ Further, by integrating (6), we find that the flux J = [[ k u e ]]∞ = K u(cid:48)( ) = [[k u e ]]∞ = K u(cid:48)(+ ) = αJ +βJ i 1 1 −∞ 1 1 i 2 1 −∞ 2 2 U V − ∇ · −∞ ∇ · − ∞ where J , J are the fluxes associated with the solutions U and V respectively. It is easy to verify U V that we can find α,β,γ to satisfy the boundary conditions u ( ) = u¯ ,u ( ) = u¯ as well as 1 1 2 2 −∞ ∞ the flux condition J = k (u¯ u¯ ). We obtain (4) . s 1 2 2 − 4 2.3 Homogenization of the sharp interface formulation Consider the situation where the domain Ω is made of a periodic microstructure at a scale ε << 1. Specifically,letY betheunitcubeconsistingoftwosubdomainsY andY separatedbyaninterface 1 2 Σ; Y = Y Y Σ. We assume that Ωε = ε(a +Y ), Ωε = ε(a +Y ), and ε = ε(a +Σ). 1 2 1 i i 1 2 i i 2 i i ∪ ∪ ∪ ∪ S ∪ We assume that the equations (4) hold in this domain with the reaction coefficient of order ε: i.e., k = εK for some K > 0 independent of ε. Peter and Böhm [27] (also see Auriault and Ene [6]) s s s show that this periodic system can be homogenized, and the homogenized equations are given by (1) where k ,k are given by the usual unit cell problem of diffusion and 1 2 f = f = K Area(Σ)(u¯ u¯ ) 1 2 s 1 2 − − where u¯ is the solution to the unit cell problem. Therefore, k ,f depend not only on the volume i i i fraction but also other aspects of the microstructure. However, we may view (2) and (3) as simple models for these. 2.4 Variational formulation The following theorem provides a variational formulation of the problem (1) above. Theorem 2.1. Let Ω Rn be a bounded, connected open set with Lipschitz boundary, ⊂ χ = χ L2(Ω;[0,1]) ∈ X { ∈ } be a given design and λ R. The problem ∈    (cid:90) 1 (cid:88) 1  inf L(u,χ) = k u 2+ χ(1 χ)u Au λχdx : u i i 2 |∇ | 2 − · − ∈ V  Ω  i=1,2 where = v H1(Ω;R2) : v = u∗ on ∂ Ω, i = 1,2 i i i V { ∈ } attains its minimum. Further, the minimum is unique and satisfies the Euler-Lagrange equation   (cid:90) (cid:88)  ki ui ϕi+χ(1 χ)ϕ Audx = 0 (7) ∇ ·∇ − · Ω i=1,2 for all ϕ = v H1(Ω;R2) : v = 0 on ∂ Ω, i = 1,2 . 0 i i ∈ V { ∈ } Proof. Set inf L(u,χ) : u = m { ∈ V} and observe that because our integrand is finite and satisfies the growth conditions λχ(x) f(x,v,ξ) c(1+ v 2+ ξ 2), − ≤ ≤ | | | | we have that < m < + . Let uν be a minimizing sequence, i.e. L(uν,χ) m as ν . For −∞ ∞ → → ∞ ν sufficiently large, (cid:90) m+1 L(uν,χ) γ uν 2 +γ uν 2 γ (x) dx γ uν 2 γ ≥ ≥ 1(cid:107)∇ (cid:107)L2 2(cid:107) (cid:107)L2 − | 3 | ≥ 1(cid:107)∇ (cid:107)L2 − 4 Ω 5 with γ > 0 independent of ν since Ω is bounded. It follows that k uν γ . (cid:107) (cid:107)W1,2 ≤ 5 appealing to our version of Poincaré’s inequality (Lemma 2.2 below). We deduce that there exists a u¯ and a subsequence (still denoted uν) that converges weakly in W1,2: uν (cid:42) u¯ in W1,2 as ∈ V ν . It follows from the convexity of the integrand (since k ,k ,k > 0) that the functional is 1 2 s → ∞ sequentially weakly lower semicontinuous. Therefore, liminfL(uν,χ) L(u¯,χ) ν→∞ ≥ and hence u¯ is a minimizer of (P). A simple calculation shows that any minimizer satisfies the Euler-Lagrange equation (7). We prove the uniqueness of the minimum by contradiction. Suppose L(u,χ) = L(v,χ) = m. Then,   (cid:90) 1 (cid:88) 1  ki( ui 2 vi 2)+ χ(1 χ)(u Au v Av)dx = 0. 2 |∇ | −|∇ | 2 − · − · Ω i=1,2 Further, since u,v , u v . Therefore, from the Euler-Lagrange equation (7) for v, we 0 ∈ V − ∈ V conclude   (cid:90) (cid:88) (cid:88)  ki vi 2+χ(1 χ)v Av ki ui vi χ(1 χ)v Audx = 0 |∇ | − · − ∇ ·∇ − − · Ω i=1,2 i=1,2 Adding these two equations,   (cid:90) 1 (cid:88)  ki ui vi 2+χ(1 χ)(u v) A(u v)dx = 0. 2 |∇ −∇ | − − · − Ω i=1,2 It follows that u = v a.e. and u v = ψ(x) 1,1 . Together, we conclude that ψ is constant i i ∇ ∇ − { } and from the boundary condition that ψ = 0. Thus u = v, giving us a contradiction. We have used the following lemma. Lemma 2.2. (Poincaré’s inequality, adapted from [18]) Let Ω and be as in the theorem above. V There exists a constant c, depending only on n and Ω such that u c u (cid:107) (cid:107)L2 ≤ (cid:107)∇ (cid:107)L2 for each function u . ∈ V Proof. We argue by contradiction. Were the stated estimate false, there would exist for each positive integer k a function uk satisfying ∈ V (cid:13) (cid:13) (cid:13) (cid:13) (cid:13)uk(cid:13) > k(cid:13) uk(cid:13) . (cid:13) (cid:13) (cid:13) (cid:13) L2 ∇ L2 We renormalize by defining uk u∗,uk u∗ vk := (cid:13) { 1 − 1 2 − 2}(cid:13) (cid:13) uk u∗,uk u∗ (cid:13) { 1 − 1 2 − 2} L2 6 and note that vk and vk = 1. It follows that 0 L2 ∈ V || || (cid:13) (cid:13) 1 (cid:13) vk(cid:13) < . (cid:13) (cid:13) ∇ L2 k In particular the functions vk ∞ are bounded in H1. It follows (e.g., [16], Thm. 12.11) that { }k=1 there exists a subsequence vkj ∞ vk ∞ and a function v L2(Ω) such that vkj v in L2. { }k=1 ⊂ { }k=1 ∈ → Further, the strong convergence implies that v and v = 1. On the other hand, the bound ∈ V0 (cid:107) (cid:107)L2 on vk from above implies that v = 0 a.e., and that v is constant since Ω is connected. Since ∇ ∇ v , v = 0 on Ω contradicting the conclusion v = 1. ∈ V0 (cid:107) (cid:107)L2 3 Optimal design problem Weseektofindthearrangementofthetwophaseswithagivenvolumeofphase1,v,thatmaximizes the normalized flux through the material: (cid:26)(cid:90) (cid:90) (cid:90) (cid:27) O := sup u∗k u nˆdA u∗k u nˆdA : χ , χdx = v . (8) 1 1 1 2 2 2 ∇ · − ∇ · ∈ X ∂1Ω ∂2Ω Ω Note that k u nˆ gives the inward flux per unit area of species i into Ω. We normalize each flux i i ∇ · by the prescribed concentration. Integrating by parts, using the variational characterization of the governing equations, and introducing a Lagrange multiplier to enforce the constraint on the given volume of phase 1, yields    (cid:90) 1 (cid:88) 1  O = sup inf L(u,χ) = k u 2+ χ(1 χ)u Au λχdx . (P) i i χ∈Xu∈V Ω 2 |∇ | 2 − · −  i=1,2 4 Characterization of the optimal design problem 4.1 Saddle point theorem We begin by showing that we can exchange the order of finding the supremum and infimum in the saddle point problem (P). Theorem 4.1. There exists v¯ , χ¯ such that ∈ V ∈ X L(v¯,χ¯) = sup inf L(v,χ) = inf supL(v,χ). χ∈Xv∈V v∈Vχ∈X for the saddle point problem (P). The proof of this draws from the following theorem adapted from Ekeland and Témam [17]. Theorem 4.2 (Proposition 2.4 of [17]). Suppose two reflexive Banach spaces V and Z satisfy (i) V is convex, closed and non-empty, A ⊂ (ii) Z is convex, closed and non-empty. B ⊂ Further let the function L : R satisfy A×B (cid:55)→ 7 (iii) u , p L(u,p) is concave and upper semicontinuous, ∀ ∈ A → (iv) p , u L(u,p) is convex and lower semicontinuous, ∀ ∈ B → (v) there exists p for bounded such that 0 ∈ B B lim L(u,p ) = + . 0 u∈A ∞ (cid:107)u(cid:107)→∞ Then L possesses at least one saddle point on . A×B We apply this theorem with L as in problem (P), V = W1,2(Ω;R2), Z = L2(Ω;[0,1]), = A V and = . Clearly, V and Z are reflexive Banach spaces as required by the theorem above. We B X now show that these satisfy the rest of hypothesis (H). Proposition 4.3. Both and are convex, closed, and non-empty. X V Proof. The point of concern is showing that our space is in fact closed. So consider a sequence X χ such that χ (cid:42)χ in L2. We seek to show that the limit function χ . From the definition i i of ∈,Xχ 1. Thus, we can pick a subsequence χ of χ such that χ∈(cid:42)X∗ χ¯, in L∞ as k . X (cid:107) i(cid:107)L∞ ≤ ik i ik → ∞ It follows χ (cid:42)χ¯, in L2 as k . Therefore, χ = χ¯ and esssupχ 1. Similarly we can show i k → ∞ ≤ essinfχ 0. Thus, χ and is closed. ≥ ∈ X X Proposition 4.4. For each χ , v L(v,χ) is convex and lower semicontinuous. ∈ X (cid:55)→ Proof. This follows trivially from the fact that the integrand in L is a sum of a positive definite quadratic term in u and a positive semidefinite quadratic form in u. i ∇ Proposition 4.5. For each v , χ L(v,χ) is concave and upper semicontinuous. ∈ V (cid:55)→ Proof. This follows trivially from the fact that the integrand in L is a sum of a positive definite quadratic term and two linear terms in χ. Proof of Theorem 4.1. From the aforementioned propositions, we have satisfied requirements (i) − (iv) of the theorem. To show (v), set χ(x) = 1/2. We have   (cid:90) 1 (cid:88) 1 1 L(u,0.5) =  k¯i ui 2+ u Au λ dx 2 |∇ | 8 · − 2 Ω i=1,2   (cid:90) 1 (cid:88) 1  k¯i ui 2 λ dx ≥ 2 |∇ | − 2 Ω i=1,2 c u c c u c 1 L2 2 3 H1 2 ≥ ||∇ || − ≥ || || − for suitable positive constants c , where we use the derived form of Poincaré’s inequality in the final i step. The requirement (v) follows. 8 4.2 Explicit characterization We are now ready to obtain the explicit characterization of the optimal design problem (P). Theorem 4.6. We have (cid:90) O = min W(v, v)dx. v∈V Ω ∇ where 1 (cid:0) ξ 2k + ξ 2k (cid:1) (v,ξ) , 2 | 1| 21 | 2| 22 ∈ R0     (cid:0)(cid:80) ∆k ξ 2(cid:1)2+2(cid:80) ξ 2(k (k +k ) 2λ∆k )+(k 2λ)2 W(v,ξ) = i i| i| i| i| v i1 i2 − i v − (v,ξ) ,  8kv ∈ R      1 (cid:0) ξ 2k + ξ 2k 2λ(cid:1) (v,ξ) 2 | 1| 11 | 2| 12− ∈ R1 with (cid:88) = (v,ξ) : ∆k ξ 2 2λ k , 0 i i v R { | | − ≤ − } i (cid:88) = (v,ξ) : k < ∆k ξ 2 2λ < k , v i i v R { − | | − } i (cid:88) = (v,ξ) : ∆k ξ 2 2λ k , 1 i i v R { | | − ≥ } i and ∆k = k k , k = k (v v )2. i i1 i2 v s 1 2 − − The function W is shown in Figure 4.2 as a function of ξ for a fixed v with various parameters. The shaded regions indicate the gradients for which mixed phases (χ (0,1)) occur. Note that ∈ mixed phases occur where the gradients of both species are comparable in magnitude, and pure phases occur otherwise. Proof. For v R2,ξ R2×2,χ R, set ∈ ∈ ∈ W(v,ξ,χ) = 1 (cid:88)(χk +(1 χ)k ) ξ 2+ ksχ(1 χ)v Av λχ i1 i2 i 2 − | | 2 − · − i=1,2 and W(v,ξ) = max W(v,ξ,χ). (9) χ∈[0,1] In light of the saddle point theorem, (cid:90) (cid:90) O = inf sup W(v, v,χ)dx = inf W(v, v)dx. v∈Vχ∈X Ω ∇ v∈V Ω ∇ It remains to compute W. To that end, note that for a fixed v,ξ, W(v,ξ,χ) is quadratic in χ and ∂W (v,ξ,χ) = 0 ∂χ 9 Figure 2: Contour plot of W for fixed v,λ. The shaded regions indicate where mixed phase (χ (0,1)) occurs. (a) k = k = k = 1,k = k = 0.1,(v v )2 = 1,λ = 0. (b) Parameters 11 22 s 12 21 1 2 ∈ − as in (a) except k = 5, (c) Parameters as in (a) except λ = 1, (d) Parameters as in (a) except 11 (v v )2 = 10. 1 2 − 10