Optimal design of a model energy conversion device

TL;DR

This paper models the optimal material arrangement in energy conversion devices, focusing on a solar-driven thermochemical process, revealing that fine phase mixtures naturally emerge as optimal through a phase-field approach.

Contribution

It formulates a coupled reaction-diffusion system for multi-material energy devices and introduces a phase-field method for optimal design, highlighting the emergence of fine material mixtures.

Findings

Fine phase mixtures are optimal in the material arrangement.

The phase-field model effectively captures the optimal design.

Numerical examples demonstrate the practical application of the approach.

Abstract

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 different materials can be critical in the performance of an energy conversion device. In this paper, we study a model problem motivated by a solar-driven thermochemical conversion device that splits water into hydrogen and oxygen. We formulate the problem as a system of coupled multi-material reaction-diffusion equations where each species diffuses selectively through 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 arise naturally. To explore this further, we introduce a phase-field formulation of the optimal design problem, and numerically study selected examples.

Tables

Table 1. Table 1: Contributions to the energy functional L ​ ( u , χ ) 𝐿 𝑢 𝜒 L(u,\chi) , the phase field regularization, the flux J i subscript 𝐽 𝑖 J_{i} of each species calculated at the boundary, and the reaction (right hand side of Eq. ( 1 ) 1 ) integrated on the domain for v = 0.5 𝑣 0.5 v=0.5 .

$k_{11}$	$k_{22}$	Transport Energy	Reaction Energy	Phase Field	$J_{1,\text{in}}$	$J_{2,\text{out}}$	Reaction
		$\int \frac{1}{2}\sum k_i{\left|\nabla u_i\right|}^2$	$\int \frac{1}{2}\chi \left(1-\chi \right)u\cdot Au$		$\int k_1\nabla u_1\cdot \hat{n}$	$-\int k_2\nabla u_2\cdot \hat{n}$
0.1	0.1	0.0451	0.0023	0.0203	0.0948	0.0948	0.0949
0.1	1	0.1706	0.0285	0.0748	0.3866	0.3977	0.3983
0.1	10	0.2952	0.2022	0.0779	0.9610	0.9947	0.9948
1	0.1	0.1706	0.0285	0.0550	0.3977	0.3864	0.3983
1	1	0.4276	0.0340	0.1070	0.9202	0.9201	0.9232
1	10	0.9044	0.3595	0.1502	2.5015	2.5256	2.5278
10	0.1	0.2953	0.2021	0.0574	0.9946	0.9602	0.9947
10	1	0.9044	0.3596	0.1150	2.5257	2.5011	2.5278
10	10	2.2730	0.9990	0.1699	6.5257	6.5254	6.5440

Equations

$$\{\begin{array}{ll}\nabla \cdot k_i\nabla u_i=f_i,& \text{in}\Omega ,\\ k_i\nabla u_i\cdot \hat{n}=0& \text{on}\partial \Omega \setminus {\partial }_i\Omega \\ u_i=u_i^{\ast }& \text{on}{\partial }_i\Omega ,\end{array}$$

$$k_1(x)=k_{11}\chi (x)+k_{12}(1-\chi (x)),k_2(x)=k_{21}\chi (x)+k_{22}(1-\chi (x));$$

$$f_1=-f_2=\chi (1-\chi )k_s(u_1-u_2)$$

$$A=k_{s}\left(\begin{array}[]{cc}1&-1\\ -1&1\end{array}\right).$$

$$\{\begin{array}{ll}\nabla \cdot K_i\nabla u_i=0& \text{in}{\Omega }_i\\ -K_i\nabla u_i\cdot \hat{m}=k_s(u_1-u_2)& \text{on}\mathcal{S}\\ u_i=u_i^{\ast }& \text{on}{\partial }_i\Omega \\ K_i\nabla u_i\cdot \hat{n}=0& \text{on}\partial \Omega \setminus {\partial }_i\Omega \end{array}$$

$$\nabla \cdot k_i^{\eta }\nabla u_i^{\eta }={\eta }^2{f}_i^{\eta },i=1,2\text{in}\Omega .$$

$$\nabla \cdot k_i\nabla u_i=f_i,i=1,2$$

$$u_i=\alpha U_i+\beta V_i+\gamma$$

$$u_1\to \alpha +\beta +\gamma \text{mbox}as{x}_1\to -\infty ,u_2\to \alpha -\beta +\gamma \text{mbox}as{x}_1\to \infty .$$

$$J=[[-k_i\nabla u_1\cdot e_1]{]}_{-\infty }^{\infty }=K_1{u}_1^{\prime }(-\infty )=[[k_i\nabla u_2\cdot e_1]{]}_{-\infty }^{\infty }=-K_2{u}_2^{\prime }(+\infty )=\alpha J_U+\beta J_V$$

$$f_1=-f_2=K_s\text{mbox}Area(\Sigma )({\bar{u}}_1-{\bar{u}}_2)$$

$$\chi \in \mathcal{X}=\{\chi \in L^2(\Omega ;[0,1])\}$$

$$\inf \left\{L(u,\chi )={\int }_{\Omega }\frac{1}{2}\sum _{i=1,2}{k}_i|\nabla u_i{|}^2+\frac{1}{2}\chi (1-\chi )u\cdot Au-\lambda \chi dx:u\in \mathcal{V}\right\}$$

$$\mathcal{V}=\{v\in H^1(\Omega ;{\mathbb{R}}^2):v_i=u_i^{\ast }\text{ on }{\partial }_i\Omega ,i=1,2\}$$

$${\int }_{\Omega }\left(\sum _{i=1,2}{k}_i\nabla u_i\cdot \nabla {\phi }_i+\chi (1-\chi )\phi \cdot Au\right)dx=0$$

$$\inf \left\{L(u,\chi ):u\in \mathcal{V}\right\}=m$$

$$-\lambda \chi (x)\le f(x,v,\xi )\le c(1+|v{|}^2+|\xi {|}^2),$$

$$m+1\ge L(u^{\nu },\chi )\ge {\gamma }_1\text{norm}{\nabla u^{\nu }}_{L^2}^2+{\gamma }_2\text{norm}{u^{\nu }}_{L^2}^2-{\int }_{\Omega }|{\gamma }_3(x)|dx\ge {\gamma }_1\text{norm}{\nabla u^{\nu }}_{L^2}^2-{\gamma }_4$$

$$\text{norm}{u^{\nu }}_{W^{1,2}}\le {\gamma }_5.$$

$$\underset{\nu \to \infty }{\mathrm{liminf}}L(u^{\nu },\chi )\ge L(\bar{u},\chi )$$

$${\int }_{\Omega }\left(\frac{1}{2}\sum _{i=1,2}{k}_i(|\nabla u_i{|}^2-|\nabla v_i{|}^2)+\frac{1}{2}\chi (1-\chi )(u\cdot Au-v\cdot Av)\right)dx=0.$$

$${\int }_{\Omega }\left(\sum _{i=1,2}{k}_i|\nabla v_i{|}^2+\chi (1-\chi )v\cdot Av-\sum _{i=1,2}{k}_i\nabla u_i\cdot \nabla v_i-\chi (1-\chi )v\cdot Au\right)dx=0$$

$$\frac{1}{2}{\int }_{\Omega }\left(\sum _{i=1,2}{k}_i|\nabla u_i-\nabla v_i{|}^2+\chi (1-\chi )(u-v)\cdot A(u-v)\right)dx=0.$$

$$\text{norm}{u}_{L^2}\le c\text{norm}{\nabla u}_{L^2}$$

$$\text{norm}{u^k}_{L^2}>k\text{norm}{\nabla u^k}_{L^2}.$$

$$v^k:=\frac{\{u_1^k-u_1^{\ast },u_2^k-u_2^{\ast }\}}{\text{norm}{\{u_1^k-u_1^{\ast },u_2^k-u_2^{\ast }\}}_{L^2}}$$

$$\text{norm}{\nabla v^k}_{L^2}<\frac{1}{k}.$$

$$O:=\sup \left\{{\int }_{{\partial }_1\Omega }{u}_1^{\ast }{k}_1\nabla u_1\cdot \hat{n}dA-{\int }_{{\partial }_2\Omega }{u}_2^{\ast }{k}_2\nabla u_2\cdot \hat{n}dA:\chi \in \mathcal{X},{\int }_{\Omega }\chi dx=v\right\}.$$

$$O=\underset{\chi \in \mathcal{X}}{\sup }\underset{u\in \mathcal{V}}{\inf }\left\{L(u,\chi )={\int }_{\Omega }\frac{1}{2}\sum _{i=1,2}{k}_i|\nabla u_i{|}^2+\frac{1}{2}\chi (1-\chi )u\cdot Au-\lambda \chi dx\right\}.$$

$$L(\bar{v},\bar{\chi })=\underset{\chi \in \mathcal{X}}{\sup }\underset{v\in \mathcal{V}}{\inf }L(v,\chi )=\underset{v\in \mathcal{V}}{\inf }\underset{\chi \in \mathcal{X}}{\sup }L(v,\chi ).$$

Full text

1 Introduction

The efficiency of fuel cells, batteries and thermochemical energy conversion devices depends on inherent material characteristics that govern the complex chemistry and transport of multiple species as well as the spatial arrangement of the various materials. Therefore, optimization of the spatial arrangement is a recurrent theme in energy conversion devices. Traditional methods of synthesis offer limited control of the microstructure and there has been much work in advanced imaging for these uncontrolled microstructures (e.g., [35]) and optimizing gross features. However, the growing ability for directed synthesis [31, 23, 33, 26] allows us to ask the question of what microgeometries are optimal for particular applications. In this sense we direct the problem to one of optimal design where we are not limited by the imagination in determining new microstructures but instead allow for the underlying physical behavior and optimization techniques to direct architecture 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 solid oxide fuel cells, the efficiency is largely determined by reactions at triple phase boundaries and 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 balance between surface reaction and transport [34, 28]. Similarly, the importance of microstructure on anode performance [5, 30, 15, 13] has also been established. The morphology of materials used in 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 exposed to inert gas at high temperature, generating oxygen vacancies in the structure, and reactant gas ($\text{H}_{2}\text{O},\text{ CO}_{2}$), at moderate temperature, releasing fuel upon reoxidation the oxide [31]. The lack 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.

There is a large literature in the study of optimal design problems, especially seeking to minimize 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 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_{1}$ and $u_{2}$ is governed by the following reaction diffusion equations for $i=1,2$:

[TABLE]

where the isotropic conductivities are

[TABLE]

with $k_{11},k_{22}>>k_{12},k_{21}>0$, and the sources are

[TABLE]

for $k_{s}>0$ and $\chi:\Omega\to[0,1]$. Briefly, we have a two-material system and $\chi$ describes the volume fraction of material 1 (say solid phase). Chemical species 1 (say bound oxygen) diffuses preferentially in material 1 ($\chi=1$) while species 2 (say oxygen gas) diffuses preferentially in material 2 (say pore, $\chi=0$). The species react and therefore there is a source at the interface $\chi\neq 0,1$. For future use, we write the source as $f=\chi(1-\chi)Au$ where $f=\{f_{1},f_{2}\}$, $u=\{u_{1},u_{2}\}$ and

[TABLE]

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 all possible arrangements $\chi$ in Section 3. We show that this gives rise to a saddle point problem. 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 $\Omega\subset{\mathbb{R}}^{n}$ with Lipschitz boundary separated into two regions $\Omega_{1}$ and $\Omega_{2}$ by an interface $\mathcal{S}$ show in Figure 1. We consider the diffusion of one species with concentration $u_{1}$ in region $\Omega_{1}$ with isotropic diffusivity $K_{1}>0$, and a second species with concentration $u_{2}$ in region $\Omega_{2}$ with isotropic diffusivity $K_{2}>0$. The two species meet at the interface and react with reaction rate $k_{s}>0$. The boundary of $\Omega$ is divided into three regions $\partial\Omega=\partial_{1}\Omega\cup\partial_{2}\Omega\cup\partial_{0}\Omega$ where $\partial_{i}\Omega\subset\partial\Omega_{i}$. The concentration of species $i$ is held at a prescribed value $u_{i}^{*}$ on $\partial_{i}\Omega$ while $\partial_{0}\Omega$ is insulating. This is described by the following system of equations:

[TABLE]

for $i=1,2$ where $\hat{m}$ represents the normal to $\mathcal{S}$ pointing from $\Omega_{1}$ pointing to $\Omega_{2}$, and $\hat{n}$ represents the outward normal to $\partial\Omega$.

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 $\chi$ be the characteristic function of $\Omega_{1}$ as defined in Section 2.1. Let $\chi^{\eta}$ be the mollification of $\chi$ with a standard mollifier at length-scale $\eta$: $\chi^{\eta}=\varphi^{\eta}*\chi$ where $\varphi^{\eta}(x)=\eta^{-n}\varphi(x/\eta)$. Let $k_{1}^{\eta},k_{2}^{\eta}$ be as in (2) with $\chi=\chi^{\eta}$ and $k_{12}^{\eta}=(1-\exp(-\eta))k_{12}$, $k_{21}^{\eta}=(1-\exp(-\eta))k_{21}$, and $f_{i}^{\eta}(x)=f_{i}(x/\eta)$. Let $u_{i}^{\eta}$ solve

[TABLE]

First consider the outer expansion $\eta\to 0$, and note that (5) formally gives (4)1,3,4 in $\Omega_{1,2}$. Further, note that $u_{1}$ (respectively $u_{2}$) is indeterminate on $\Omega_{2}$ (respectively $\Omega_{1}$). However, this outer expansion does not give any condition on the interface $\mathcal{S}$. To obtain this condition, denote the limiting values on the interface to be $\bar{u}_{1},\bar{u}_{2}$. We seek to relate these to the flux as in (4)2.

Now consider the inner expansion. Pick a point $x_{0}\in{\mathcal{S}}$ and change variables $x\mapsto(x-x_{0})/\eta$. We obtain

[TABLE]

where $k_{i}=k_{i}^{1}$. Further, as $\eta\to 0$, $\chi$ and hence the solution depend only on one dimension that is normal to the interface. We take this direction to be $x_{1}$ by changing variables if necessary. Let $U_{i}$ solve (6) for the boundary conditions $(u_{1},u_{2})(x_{1})\to(1,0)$ as $x_{1}\to-\infty$ and $(u_{1},u_{2})(x_{1})\to(0,1)$ as $x_{1}\to\infty$, and $V_{i}$ solve (6) for the boundary conditions $(u_{1},u_{2})(x_{1})\to(1,0)$ as $x_{1}\to-\infty$ and $(u_{1},u_{2})(x_{1})\to(0,-1)$ as $x_{1}\to\infty$. Note that

[TABLE]

also solves (6) for any arbitrary scalars $\alpha,\beta,\gamma$, and satisfies the boundary conditions

[TABLE]

Further, by integrating (6), we find that the flux

[TABLE]

where $J_{U}$, $J_{V}$ are the fluxes associated with the solutions $U$ and $V$ respectively. It is easy to verify that we can find $\alpha,\beta,\gamma$ to satisfy the boundary conditions $u_{1}(-\infty)=\bar{u}_{1},u_{2}(\infty)=\bar{u}_{2}$ as well as the flux condition $J=k_{s}(\bar{u}_{1}-\bar{u}_{2})$. We obtain (4)2.

2.3 Homogenization of the sharp interface formulation

Consider the situation where the domain $\Omega$ is made of a periodic microstructure at a scale $\varepsilon<<1$. Specifically, let $Y$ be the unit cube consisting of two subdomains $Y_{1}$ and $Y_{2}$ separated by an interface $\Sigma$; $Y=Y_{1}\cup Y_{2}\cup\Sigma$. We assume that $\Omega^{\varepsilon}_{1}=\cup_{i}\varepsilon(a_{i}+Y_{1})$, $\Omega^{\varepsilon}_{2}=\cup_{i}\varepsilon(a_{i}+Y_{2})$, and ${\mathcal{S}}^{\varepsilon}=\cup_{i}\varepsilon(a_{i}+\Sigma)$. We assume that the equations (4) hold in this domain with the reaction coefficient of order $\varepsilon$: i.e., $k_{s}=\varepsilon K_{s}$ for some $K_{s}>0$ independent of $\varepsilon$. Peter and Böhm [27] (also see Auriault and Ene [6]) show that this periodic system can be homogenized, and the homogenized equations are given by (1) where $k_{1},k_{2}$ are given by the usual unit cell problem of diffusion and

[TABLE]

where $\bar{u}_{i}$ is the solution to the unit cell problem. Therefore, $k_{i},f_{i}$ depend not only on the volume 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 $\Omega\subset\mathbb{R}^{n}$ be a bounded, connected open set with Lipschitz boundary,

[TABLE]

be a given design and $\lambda\in{\mathbb{R}}$. The problem

[TABLE]

where

[TABLE]

attains its minimum. Further, the minimum is unique and satisfies the Euler-Lagrange equation

[TABLE]

for all $\varphi\in{\mathcal{V}}_{0}=\{v\in H^{1}(\Omega;\mathbb{R}^{2}):v_{i}=0\text{ on }\;\partial_{i}\Omega,\,i=1,2\}$.

Proof.

Set

[TABLE]

and observe that because our integrand is finite and satisfies the growth conditions

[TABLE]

we have that $-\infty<m<+\infty$. Let $u^{\nu}$ be a minimizing sequence, i.e. $L(u^{\nu},\chi)\rightarrow m$ as $\nu\to\infty$. For $\nu$ sufficiently large,

[TABLE]

with $\gamma_{k}>0$ independent of $\nu$ since $\Omega$ is bounded. It follows that

[TABLE]

appealing to our version of Poincaré’s inequality (Lemma 2.2 below). We deduce that there exists a $\bar{u}\in\mathcal{V}$ and a subsequence (still denoted $u^{\nu}$) that converges weakly in $W^{1,2}$: $u^{\nu}\rightharpoonup\bar{u}$ in $W^{1,2}$ as $\nu\rightarrow\infty$. It follows from the convexity of the integrand (since $k_{1},k_{2},k_{s}>0$) that the functional is sequentially weakly lower semicontinuous. Therefore,

[TABLE]

and hence $\bar{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,\chi)=L(v,\chi)=m$. Then,

[TABLE]

Further, since $u,v\in{\mathcal{V}}$, $u-v\in{\mathcal{V}}_{0}$. Therefore, from the Euler-Lagrange equation (7) for $v$, we conclude

[TABLE]

Adding these two equations,

[TABLE]

It follows that $\nabla u_{i}=\nabla v_{i}$ a.e. and $u-v=\psi(x)\{1,1\}$. Together, we conclude that $\psi$ is constant and from the boundary condition that $\psi=0$. Thus $u=v$, giving us a contradiction. ∎

We have used the following lemma.

Lemma 2.2**.**

*(Poincaré’s inequality, adapted from [18])

Let $\Omega$ and ${\mathcal{V}}$ be as in the theorem above. There exists a constant $c$, depending only on $n$ and $\Omega$ such that*

[TABLE]

for each function $u\in\mathcal{V}$.

Proof.

We argue by contradiction. Were the stated estimate false, there would exist for each positive integer $k$ a function $u^{k}\in\mathcal{V}$ satisfying

[TABLE]

We renormalize by defining

[TABLE]

and note that $v^{k}\in{\mathcal{V}}_{0}$ and $||v^{k}||_{L^{2}}=1$. It follows that

[TABLE]

In particular the functions $\{v^{k}\}_{k=1}^{\infty}$ are bounded in $H^{1}$. It follows (e.g., [16], Thm. 12.11) that there exists a subsequence $\{v^{k_{j}}\}_{k=1}^{\infty}\subset\{v^{k}\}_{k=1}^{\infty}$ and a function $v\in L^{2}(\Omega)$ such that $v^{k_{j}}\rightarrow v\text{ in }L^{2}.$ Further, the strong convergence implies that $v\in{\mathcal{V}}_{0}$ and $\norm{v}_{L^{2}}=1$. On the other hand, the bound on $\nabla v^{k}$ from above implies that $\nabla v=0$ a.e., and that $v$ is constant since $\Omega$ is connected. Since $v\in{\mathcal{V}}_{0}$, $v=0$ on $\Omega$ contradicting the conclusion $\norm{v}_{L^{2}}=1$.

∎

3 Optimal design problem

We seek to find the arrangement of the two phases with a given volume of phase 1, $v$, that maximizes the normalized flux through the material:

[TABLE]

Note that $k_{i}\nabla u_{i}\cdot\hat{n}$ gives the inward flux per unit area of species $i$ into $\Omega$. We normalize each flux 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

[TABLE]

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 $\bar{v}\in{\mathcal{V}}$, $\bar{\chi}\in{\mathcal{X}}$ such that

[TABLE]

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

1. (i)

$\mathcal{A}\subset V$* is convex, closed and non-empty,* 2. (ii)

$\mathcal{B}\subset Z$* is convex, closed and non-empty.*

Further let the function $L:\mathcal{A}\times\mathcal{B}\mapsto\mathbb{R}$ satisfy

1. (iii)

$\forall u\in\mathcal{A},\,p\rightarrow L(u,p)$* is concave and upper semicontinuous,* 2. (iv)

$\forall p\in\mathcal{B},\,u\rightarrow L(u,p)$* is convex and lower semicontinuous,* 3. (v)

there exists $p_{0}\in\mathcal{B}$ for $\mathcal{B}$ bounded such that

[TABLE]

Then $L$ possesses at least one saddle point on $\mathcal{A}\times\mathcal{B}$.

We apply this theorem with $L$ as in problem (P), $V=W^{1,2}(\Omega;\mathbb{R}^{2})$, $Z=L^{2}(\Omega;[0,1])$, ${\mathcal{A}}={\mathcal{V}}$ and $\mathcal{B}=\mathcal{X}$. Clearly, $V$ and $Z$ are reflexive Banach spaces as required by the theorem above. We now show that these satisfy the rest of hypothesis (H).

Proposition 4.3**.**

Both $\mathcal{X}$ and $\mathcal{V}$ are convex, closed, and non-empty.

Proof.

The point of concern is showing that our space $\mathcal{X}$ is in fact closed. So consider a sequence $\chi_{i}\in\mathcal{X}$ such that $\chi_{i}{\rightharpoonup}\chi$ in $L^{2}$. We seek to show that the limit function $\chi\in\mathcal{X}$. From the definition of $\mathcal{X}$, $\norm{\chi_{i}}_{L^{\infty}}\leq 1.$ Thus, we can pick a subsequence $\chi_{i_{k}}$ of $\chi_{i}$ such that $\chi_{i_{k}}\overset{\ast}{\rightharpoonup}\bar{\chi},$ in $L^{\infty}$ as $k\to\infty$. It follows $\chi_{i_{k}}{\rightharpoonup}\bar{\chi},$ in $L^{2}$ as $k\to\infty$. Therefore, $\chi=\bar{\chi}$ and $\operatorname*{ess\,sup}\chi\leq 1$. Similarly we can show $\operatorname*{ess\,inf}\chi\geq 0$. Thus, $\chi\in{\mathcal{X}}$ and $\mathcal{X}$ is closed. ∎

Proposition 4.4**.**

For each $\chi\in\mathcal{X},$ $v\mapsto L(v,\chi)$ is convex and lower semicontinuous.

Proof.

This follows trivially from the fact that the integrand in $L$ is a sum of a positive definite quadratic term in $\nabla u_{i}$ and a positive semidefinite quadratic form in $u$.

∎

Proposition 4.5**.**

For each $v\in\mathcal{V},$ $\chi\mapsto L(v,\chi)$ is concave and upper semicontinuous.

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 $\chi$. ∎

Proof of Theorem 4.1.

From the aforementioned propositions, we have satisfied requirements $(i)-(iv)$ of the theorem. To show $(v)$, set $\chi(x)=1/2$. We have

[TABLE]

for suitable positive constants $c_{i}$, where we use the derived form of Poincaré’s inequality in the final step. The requirement $(v)$ follows.

∎

4.2 Explicit characterization

We are now ready to obtain the explicit characterization of the optimal design problem (P).

Theorem 4.6**.**

We have

[TABLE]

where

[TABLE]

with

[TABLE]

and

[TABLE]

The function $\overline{W}$ is shown in Figure 2 as a function of $\xi$ for a fixed $v$ with various parameters. The shaded regions indicate the gradients for which mixed phases ($\chi\in(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\in{\mathbb{R}}^{2},\xi\in{\mathbb{R}}^{2\times 2},\chi\in{\mathbb{R}}$, set

[TABLE]

and

[TABLE]

In light of the saddle point theorem,

[TABLE]

It remains to compute $\overline{W}$. To that end, note that for a fixed $v,\xi$, $W(v,\xi,\chi)$ is quadratic in $\chi$ and

[TABLE]

has an unique solution $\chi=\chi^{*}$. A simple calculation reveals

[TABLE]

Notice that

[TABLE]

A longer, but straightforward, calculation reveals that

[TABLE]

Similarly,

[TABLE]

Now, we can verify by explicit calculation that

[TABLE]

We obtain the desired result by recalling (9), rewriting

[TABLE]

where

[TABLE]

and using (4.2).

∎

5 Phase-field formulation of the optimal design problem

The min-max problem based on the functional $L$ is difficult to solve numerically due to the fact that $\chi$ is only in $L^{2}$ and because of the constraint $\chi\in[0,1]$. The relaxed functional is also difficult to solve numerically since $\overline{W}$ is not strictly convex. Therefore, we now pursue an alternative approach to the optimal design problem that is amenable to numerical treatment. We regularize the functional $L(u,\chi)$ by adding the $L^{2}$ norm of $\nabla\chi$ and requiring $\chi\in H^{1}$. We also replace the constraint $\chi\in[0,1]$ with a penalty. Finally, from a practical point of view, it would also be beneficial to have solutions that prefer the pure phases $\chi\in\{0,1\}$. Therefore, we add a term to the energy that penalizes any deviation from this set.

We consider the functional

[TABLE]

where

[TABLE]

has two wells at $\chi\in\{0,1\}$.

The additional terms in parenthesis form the integrand of the Allen-Cahn functional [2]. Minimizers of this functional partition the domain into regions where $\chi\approx 0$ and $\chi\approx 1$ separated by transition layers with thickness $\sim\sqrt{\beta/\alpha}.$ In our setting, we expect this to be modified by the transport energy.

We seek to find the saddle point by considering a gradient flow:

[TABLE]

for every $\varphi,\psi\in H^{1}(\Omega;\mathbb{R}^{N})$ subject to the appropriate boundary conditions where $\langle\cdot\rangle$ denotes the $L^{2}$ inner product and $d_{\chi},d_{u}>0$ are the inverse mobilities. We obtain the following system of equations:

[TABLE]

6 Numerical study of the optimal design problem

We have implemented the phase field formulation of the optimal design problem (11, 12) using the commercial software COMSOL [14]. All our simulations are in two dimensions ($n=2$). We work with non-dimensional units where the size of the domain, the concentration at a boundary and the (diagonal components of the) diffusion coefficient are $\mathcal{O}(1)$. We discretize the problem spatially using linear finite elements generated by Delaunay triangulation, and integrate the resulting ordinary differential equation in time by using the backward differentiation formula. We impose the volume constraint as a global constraint that is built into COMSOL. Additionally, we impose a point-wise constraint restricting $\chi\in[0,1]$. We typically begin with an initial guess of uniform $\chi$, and run the simulations until an apparent steady state is reached (i.e., when the right hand sides of (11, 12) become small compared to a given tolerance). The simulations can get stuck in local optima, but we try to avoid this by doing parameter sweeps and studying additional initial conditions.

6.1 Square reactor

We begin with a square domain, $\Omega=(0,1)^{2}$, shown in Figure 3(a). We prescribe $u_{1}=1$ on the left face $\partial_{1}\Omega=\{0\}\times(0,1)$ corresponding to a source of species 1, $u_{2}=0$ on the right face $\partial_{2}\Omega=\{1\}\times(0,1)$ corresponding to a sink of species 2, and zero-flux boundary conditions otherwise. We also impose a zero flux boundary condition on our phase-field variable $\chi$.

The resulting optimal design $\chi$ is shown in Figure 3(b) for the parameters

[TABLE]

This simulation had a mesh with 67068 elements, took $50$ non-dimensional units of time over 845 time steps and the $L^{2}$ norm of the time derivative of $\chi$ is $3.945\times 10^{-4}$ at the end of the simulation. We have verified that the design does not change by refining the mesh and driving the $L^{2}$ norm of the time derivative of $\chi$ to $10^{-12}$.

The resulting design has a clear intuitive explanation. Given the boundary conditions, the design seeks to draw in species 1 from the left, react it in the center to convert species 1 to species 2, and expel species 2 at the right. Therefore, the design puts material 1, which has a high diffusivity of species 1, on the left so that it can easily transport species 1 from the source to the interface where the reaction consumes it. Material 2, which has a high diffusivity of species 2, is placed on the right so that it can easily transport species 2 from the interface, where the reaction generates it, to the sink. The design maximizes the reaction by creating a zig-zag interface between the two materials.

We now begin a parameter study for the same problem. Figure 4 shows the resulting designs for a volume fraction $v=0.5$ for various diffusivities $k_{11},k_{22}$. Figures 5 and 6 show the corresponding concentration fields $u_{1}$ and $u_{2}$ respectively while Figure 7 shows the corresponding reactions.

We begin at the center for the case $k_{11}=k_{22}=1$, which is what we described earlier. Decreasing both diffusivities by moving up on the diagonal to $k_{11}=k_{22}=0.1$ leads to a similar segregation of the material but the interface is sharper and straight. On the other hand, increasing both diffusivities by moving down the diagonal to $k_{11}=k_{22}=10$ still segregates the material, but in a very diffuse manner with an almost constant gradient. Note that the interface width changes despite the fact that length-scale, $\sqrt{\beta/\alpha}$, predicted by the phase-field alone is held fixed. This is because of the relative importance of the diffusion and the reaction. When the diffusivities are both small, $k_{11}=k_{22}=0.1$ as in the upper-left, the reaction is relatively easy and diffusion difficult. Thus one only needs a narrow region for the reaction, saving much of the pure material for optimal transport. Conversely, when the diffusivities are both large, $k_{11}=k_{22}=10$ as in the bottom-right, the reaction is relatively difficult and diffusion easy. Thus, one creates a very diffuse interface to optimize the reaction.

We now turn to the situation when the diffusivities are different. Consider the case when $k_{11}=1,k_{22}=0.1$ as shown on the top-center. The diffusion of species 1 is considerably easier than that of species 2. Therefore, it is advantageous to have the reaction close to the sink. Species 1 is transported by the long arms of material 1 (red) which protrude from the left to the right where it reacts very close to the sink, thereby reducing the distance that species 2 has to be transported. The excess material 2 (blue) is ‘hidden’ on the left in arms that do not participate in the transport. The case $k_{11}=10,k_{22}=0.1$ shown on the top-right is similar with a slightly wider interface since reaction is more difficult compared to the transport. The case $k_{11}=10,k_{22}=1$ shown on the right-middle is also similar except the interfacial region is even wider. The cases $k_{11}=0.1,k_{22}=1$; $k_{11}=0.1,k_{22}=10$ and $k_{11}=1,k_{22}=10$ are the analogous, with the roles of material 1 and 2 reversed.

The phase-field functional, the domain, and the boundary conditions have a symmetry, and we examine if the resulting designs reflect this symmetry. Specifically, note that if $\{u_{1},u_{2},\chi\}$ is a solution for a problem with $k_{1},k_{2}$ on the square domain, then $\{1-u_{2},1-u_{1},1-\chi\}$ is a solution for a problem with $k_{2},k_{1}$ on the square domain obtained by changing $x$ to $1-x$. We see that our designs reflect this symmetry. Specifically, compare the case $k_{11}=0.1,k_{22}=1$ and the resulting design $\chi_{1}$ shown in middle-left of Figure 4 and the case $k_{11}=1,k_{22}=0.1$ and the resulting design $\chi_{2}$ shown in top-center of Figure 4. We see that $\chi_{1}(x,y)\approx 1-\chi_{2}(1-x,1-y)$.

Table 1 shows how the different contributions to the energy change for the various cases. It also shows how the flux varies. Further, it shows the the flux at the source, sink, and reaction zone all agree.

Figure 8 shows the designs for the same parameters, but for a volume fraction $v=0.3$. The designs are similar, except the interface is more to the left.

6.2 Cylindrical reactor

Many reactors designed for thermochemical conversion devices implement a cylindrical ceramic structure that allow for even heating and easy transport of recant gas. Thus, for the second example we look at an annular structure where the inner edge with $r=0.2$ is held as the source of the first chemical species ($\partial_{1}\Omega$ where $u_{1}=1$) and the outer at $r=1$ is set as a sink for the second ($\partial_{2}\Omega$ where $u_{2}=0$). We consider the same parameters as (13). The resulting design is shown in Figure 9(b). The first species enters from the inside, reacts and converts to the second species which exits from the outside. Thus, we see much of the first material on the inside and the second on the outside. Further, to enable sufficient reaction, the interface region is graded. If we decrease the phase field coefficients by an order of magnitude, we obtain the design in Figure 9(a) where the mixed region increases as the penalty for deviating from the pure materials is reduced. On the other hand, increasing the phase field coefficients by an order of magnitude yields the design in Figure 9(c). Indeed, here, the penalty for deviation from the pure phases increases and therefore the interface becomes corrugated allowing sufficient reaction.

The second row of Figure 9 show the analogous result when the diffusivity is reduced by an order of magnitude. Transport is now harder compared to the reaction, and therefore nearly pure phases dominate to ensure transport and complex interfaces are avoid due to the phase field. Again, increasing the phase field parameters promotes pure phases. The final row of Figure 9 show the results for unequal conductivity. Since the transport of first species is easier, material 1 forms long arms to transport the first species to close to the outlet where the reaction takes. Further, increasing the phase field parameters promotes pure phases and leads to fewer arms.

6.3 Periodic cellular reactor

It is common to construct reactors as a periodic tubular array where the first species enter the reactor through one set of tubes while the second species is extracted from the reactor with a different set of tubes. Looking at a cross-section, one sees a square array of inlets and a square array of outlets. This motivates our next example where the reactor is taken to be periodic with the unit cell shown in Figure 10(a). The source is at the corners of the cell while the outlet is at the center. We look for a periodic design to optimize the transport as before. The resulting unit design for the parameters shown in (13) is shown in Figure 10(b). It is repeated periodically in Figure 10(c).

Acknowledgement

It is a pleasure to acknowledge many interesting discussions with Sossina M. Haile and Robert V. Kohn. We gratefully acknowledge the financial support of the National Science Foundation through the PIRE grant OISE-0967140.