Optimal design of a micro-tubular fuel cell
G. Delgado1 ∗.
March 26, 2018

1

CMAP, École Polytechnique, Palaiseau, France.

arXiv:1603.08260v1 [] 27 Mar 2016

Abstract
We discuss the problem of the optimal design of a micro-tubular fuel cell applying an inverse
homogenization technique. Fuel cells are extremely clean and efficient electrochemical power generation devices, made up of a cathode/electrolyte/anode structure, whose energetic potential has not
being fully exploited in propulsion systems in aeronautics due to their low power densities. Nevertheless, thanks to the recent development of additive layer manufacturing techniques (3D printing),
complex structures usually impossible to design with conventional manufacturing techniques can be
constructed with a low cost, allowing notably to build porous or foam-type structures for fuel cells.
We seek thus to come up with the micro-structure of an arrangement of micro-tubular cathodes which
maximizes the contact surface subject to a pressure drop and a permeability constraint. The optimal
periodic design (fluid/solid) emerges from the application of a shape gradient algorithm coupled to a
level-set method for the geometrical description of the corresponding cell problem.

Keywords: Shape and topology optimization; Solid oxide fuel cells; Level-set method; Homogenization.

1

Introduction

Fuel cells are energy conversion devices which can continuously convert chemical energy into electrical
energy and heat, without involving direct combustion. This feature offers many advantages over traditional power sources such as improved efficiency, greater fuel diversity, high scalability, no moving parts
(hence less noise and vibration) and lower impact on the environment [21].
A fuel cell is a fairly simple device, mainly composed of three elements: an anode, a cathode, and
an electrolyte between the two electrodes. The two electrodes are connected together by an electrical
circuit. On the surface of these electrodes, electrochemical reactants react through half-redox reactions,
producing (or consuming) ions, electrons, and in most cases, heat. Ions pass through the electrolyte
meanwhile electrons are “channeled” in the electric circuit and then routed to the second electrode to be
consumed.
The efficiency of current fuel cells, which ranges from 40% to 60%, is higher than thermal systems such
as gas turbines since their operation is not constrained by any theoretical thermodynamic limitation as the
maximum Carnot efficiency. Furthermore, coupled to a gas turbine at a high temperature (800−1000o C),
the spawned hybrid system can achieve really high efficiencies (near 85%), producing electricity from the
waste of heat. This feature makes high temperature operating fuel cells an interesting complement to
aircraft engines [38, 48, 40, 39, 50].
Among the different types of fuel cells (see Fig. 1), one can recognize the Solid Oxide Fuel Cells (SOFC)
as a particularly appealing model for high temperature applications in hybrid systems. SOFC possesses
various advantages w.r.t other fuel cells such as a solid electrolyte allowing different geometries and
shapes, good performance and durability, and high operating temperature (800 − 1000o C) for reforming.
This last attribute has the potential of hybridization with a gas turbine, as it is shown in Fig. 3.
The cathode of an SOFC is usually an alloy of lanthanum, strontium, and manganese oxide. On its
surface takes place the half reduction reaction of oxygen, producing an oxygen ion. The electrolyte is
e.g. made of YSZ (yttrium stabilized zirconia). This allows the ion transport to the anode (nickel mixed
with YSZ), in contact with which it reacts with hydrogen to produce water and electrons. Finally, these
∗ Current

address: IRT SystemX, Palaiseau, France (gabriel.delgado@irt-systemx.fr)

1

Operating temperature (°C)

Electrical efficiency [%]

90

AFC
60

SOFC

PEMFC
MCFC
PAFC

1000

SOFC
MCFC
500

PAFC

AFC
PEMFC

100

30

1

10

100

1000

1

10000

10

100

1000

10000

Typical stack size (kW)

Typical stack size (kW)

Figure 1: Different types of fuel cells: Polymeric Electrolyte Membrane (PEMFC), Alkaline (AFC), Phosphoric
Acid (PAFC), Molten Carbonate (MCFC) and Solid Oxide (SOFC). Reprinted with permission from
[34].

electrons are conducted and used by an electric device (Figure 2). Both the anode and the cathode must
be very porous so as to allow the transport of the fuel and the oxygen, respectively.
SOFC can be designed following many geometric configurations. However, the most common designs
are the planar and the tubular ones. Each design offers advantages and drawbacks depending of the
application requirement. Planar design configuration has a low physical component volume profile and
short current path between single cells, allowing higher power densities, meanwhile tubular designs have
e.g. high thermo-mechanical properties, simple sealing requirements and good thermal shock resistance
[19].
Despite of the impressive energy efficiencies achieved by the SOFC, the foregoing geometries lack of
an essential property required by any aircraft system: high gravimetric and volumetric energy densities.
Indeed, fuel cells are still too heavy to propel any large aircraft since they have a lower power density
when compared with conventional turbines [25]. For instance, the gravimetric power density (kW/kg) of
a SOFC compared to the turbojet CFM-56 (one of the most popular turbojets in the world and often
used as a reference for studies of this type) is at least five times smaller.
In order to improve the fuel cell efficiency, various alternatives have emerged, such as miniaturization
of the structure through a micro-tubular configuration [19] and the use of porous materials (foam) in
the design of some components (e.g. bipolar/end plates) [33]. In general these techniques take advantage of recent developments in additive layer manufacturing or 3D printing procedures, which allow the
construction of extremely complicated 3D structures from a CAD model.
Thus, the design of fuel cells with very small features, leading notably to periodic patterns, becomes
a foremost challenge in the construction of future technologies (see Figure 4). An adapted tool for this
purpose is inverse homogenization. The homogenization method for topology optimization, which has
been successfully implemented in structural optimal design (see e.g. [1, 7, 15, 5, 9]), consists in admitting
composites regions where the void appears at a micro-scale. The theoretical foundations are given in
[32, 51]. Meanwhile in homogenization the effective properties of a material are found from the microstructure, in inverse homogenization the micro-structure does not exist initially but we seek to come
up with a micro-structure with prescribed or extreme homogenized properties. The design of materials
with extreme or prescribed properties using inverse homogenization in elasticity, fluid mechanics, wave
propagation, etc., is well known in the literature [14, 46, 45, 23, 47, 53, 22, 8].
In the process of finding the optimal micro-structure, a topology optimization problem must be solved
within the so-called cell problem. The cell problem defines the link between the physical macroscopic
(effective) properties and the micro-structure characterized by a certain geometry. Among the most
popular methods in topology optimization for fluid flow problems, the level-set method for shape optimization [4, 13, 44, 52] arises as a viable, robust and efficient alternative to more standard density-based
approaches [1, 3, 7]. First introduced in [37], the level-set method has the advantage of tracking the interfaces on a fixed mesh, easily managing topological changes without any need of re-meshing. Allied to the
Hadamard method of shape differentiation, the level-set approach gives a better description and control
2

e-

e-

Fuel

Oxydant

Byproduct

Byproduct

Figure 2: Triple structure of a solid oxide fuel cell. The O2 from the air is consumed in the cathode, meanwhile
the O2− ions liberated from the reduction reaction travel through the electrolyte to the anode, where
the oxidation reaction of the fuel (H2 ) takes place, liberating heat, water and electricity. Reprinted
with permission from [34].

of the geometrical properties of the interface without need of any intermediate density, avoiding typical
drawbacks such as intermediate density penalization and possible spurious physical behavior during the
optimization process.
In the present article we propose an optimal configuration of a micro-tubular fuel cell, whose cathode
constitutes a tubular periodic structure designed by inverse homogenization. This particular structure
leads to a maximal surface micro-structure subjected to a permeability constraint for the air phase and a
pressure drop constraint for the fuel phase. In Section 2, the main features of the problem are described,
namely the physical modeling of a porous fuel cell, the homogenization of the equations governing the
system and the optimization problem we contemplate to solve. Then the shape gradients of the functionals
involved in the optimization problem are recalled in Section 3. Finally an example of an optimal periodical
micro-structure (fluid/solid) is detailed in Section 4, emerging from the application of the level-set method
for topology optimization to the corresponding cell problem.

2
2.1

Problem setting
Physical modeling of a porous fuel cell

A generic SOFC, regardless of the geometrical configuration (planar, tubular, monolithic, etc), is always
composed of two porous electrodes (anode and cathode), a dense electrolyte, an anodic and cathodic
gas channel and two current collectors. However, for the sake of simplicity, we propose to study a
reduced model of a micro-tubular SOFC [19], as described in Figure 5, focusing on the three former
components. Thus, we consider a periodically perforated domain, where the air flows freely around a
periodic collection of cathode/electrolyte/anode tubes. The air cannot penetrate inside the tubes and
the fuel flows orthogonally through them. The air transports several species but we concentrate on the
O2 , which is consumed in the electrochemical reaction on the surface of the cathode.
Assumptions.

Now we list the main physical assumptions of our model, similarly to [11].

(1) Steady state velocity of the air (also called slow steady state).

3

Current

2
Combustion
Chamber

1

Fuel Cell
3

4

High Pressure
Compressor

Current Generator

High Pressure
Turbine

Low Pressure
Compressor

Low Pressure
Turbine

Air flow
Gas flow
Semi-Exhaust gas flow
Exhaust gas flow

Figure 3: Example of an hybrid system. Reprinted with permission from [34].

Power density

Foam- like
geometries

? ? ?
?
?

Monolithic
geometry
Planar
geometry
Tubular
geometry

•

Existing geometries

Active surface area per volume

Figure 4: Evolution of power density vs active surface of SOFC. Reprinted with permission from [34].

(2) Laminar and incompressible flow. This assumption derives from the low gas speed in the SOFC gas
channels, where the density variation of each specie is not related to compressions/expansions, but
rather caused either internally by heat release of chemical reactions or externally by wall heating
and by mass variation. Furthermore, considering that the gas speed in SOFC gas channels is always
very low (low Mach number < 0.3), it is a common practice to assume a laminar flow in the gas
channels, meaning that the non-linear term in the momentum conservation equation (inertia term)
is negligible with respect to the viscous term.
(3) Isothermal state. Since the temperature of the SOFC is quite high (800o − 1000o C), the local
variations can be underestimated. This allows in particular to avoid the dependence of the diffusion
tensor and the reaction ratio of O2 with respect to the temperature.
(4) The electrochemical reactions are confined to the electrode-electrolyte interface. The place where the
electrochemical reaction takes place, the so-called triple-phase-boundary, is the site where ions, electrons and gas coexist, thus enabling the redox reactions. Since this place represents a small portion
of the entire electrode domain and the high electronic conductivity of the electrodes compared to the
ionic conductivity, the redox reactions are likely to take place very close to the electrode-electrolyte
interface.
(5) The cathode/electrolyte/anode is a lumped structure so it can be treated as one interface.
Generally speaking, the physic of the SOFC is modeled as a complex system where chemical reactions,
electrical conduction, ionic conduction, gas phase mass transport, and heat transfer take place simultane4

Figure 5: Micro-tubular design with arbitrary shaped tubes proposed in [34] (reprinted with permission). See
also the similar micro-tubular SOFC design of AIST [36].

ously and are tightly coupled [28, 29, 17]. Nevertheless, we simplify the physical description by focusing
only on the gas mass transport and the chemical reaction on the cathode.

ΓN
Ω
ΓD

Γ

ΓN
Figure 6: Transversal section of the porous SOFC and boundary conditions. The air flows from the left to the
right. The fuel penetrates orthogonally through the figure inside the red tubes. Γ represents the
anode/electrolyte/cathode lumped structure.

Consider Fig. 6. Let Ω ⊂ Rd be the porous volume of the fuel cell where the air flows freely, delimited
by the boundaries ΓD , ΓN and Γ . The air is injected through ΓD , meanwhile ΓN is impermeable. A
periodic arrangement of anode/electrolyte/cathode tubes constitutes the interface Γ . The redox reaction
takes place on Γ and the adimensional parameter  corresponds to the ratio between the characteristic
size of each tube and a macroscopic characteristic length.
Conservation and constitutive laws. The fluid adimensionalized equations (mass conservation and
momentum conservation) plus the specie (O2 ) diffusion-convection-(surface)reaction adimensionalized
equation respectively read [6]

div(u ) = 0
x ∈ Ω ,



−2 µ∆u = ∇p
x ∈ Ω ,
(1)
∂X
+
u
·
∇X
=
λ∆X
t > 0, x ∈ Ω ,





 ∂t
X = Xinit
t = 0, x ∈ Ω ,
where X , Xinit represent the current and the initial concentrations of O2 in the air, u is the local velocity
of the air, p the local pressure, µ the viscosity of the air and λ the diffusion coefficient. The parameter
 as it was explained above, corresponds to the ratio between the characteristic size of each tube and a

5

macroscopic characteristic length. We remark that in this model, the diffusion and the convection of the
O2 in the transport equation are equilibrated at the macro-scale.

Boundary conditions

X = XD
t > 0, x ∈ ΓD ,
 u = uD
u · n = 0 D ∇X · n = 0
t > 0, x ∈ ΓN ,

1
u = 0
∇X · n = − 4e
R(X ) t > 0, x ∈ Γ ,

(2)

where (uD , XD ) are the velocity of the air and the concentration of O2 at the inflow boundary ΓD , e is
the electron charge and R(X ) is the Butler-Volmer reaction term [30], which can be simplified [42, 11]
to

zF 
η .
R(X ) = j0 (Xm , T ) exp − αc
RT
Here, j0 corresponds to the exchange current density and depends linearly upon Xm , m is the order of
reaction, η is the over-potential with respect the Nerst equilibrium potential, αc is the cathodic transfert
coefficient, T the temperature, z the number of electrons involved in the electrode reaction, R the universal
gas constant and F the Faraday constant. All the above parameters are assumed to be constants.
We remark that the coefficient m is positive and in order to prove the existence of a weak solution of
(1), (2), m should be smaller than a certain critical value m̄, which depends on the dimension d of the
space. By taking e.g. m ≤ 1, R(X ) ∈ L2 (Γ ) and the variational formulation of (1), (2) will be well
posed. Typical values for m are e.g. 0.4 in [31], for a composite cathode made of 60% volume of LSM
(lanthanum-strontium-manganite) and 40% YSZ, and 0.25 in [16, 11]. Thus, depending on the physical
characteristics of the cathode, the reaction term R is in general non-linear.

2.2

Homogenized system

Before stating the homogenization of the system (1),(2), let us describe more precisely the assumptions on
the porous domain Ω [1]. As usual in periodic homogenization, a periodic structure is defined by a domain
Ω and an associated micro-structure, or periodic cell Y = (0, 1)d , which is made of two complementary
parts : the solid part ω and the fluid part Y \ω (see Figure 7). We assume that Y \ω is a smooth and
connected open subset of Y , identified with the unit torus (i.e. Y \ω, repeated by Y −periodicity in Rd ,
is a smooth and connected open set of Rd ). The domain Ω is covered by a regular mesh of size : each
cell Yi is of the type (0, )d , and is divided in solid part ωi and fluid part Yi \ωi , i.e. is similar to the
unit cell Y rescaled to size . The fluid part Ω of a porous medium is defined by
N ()

[

Ω = Ω\

ωi ,

i=1

where the number of cells is N () = |Ω|−d (1 + o(1)).

Y

Y \ω

ω

Figure 7: Unit cell of a porous medium.

6

Now according to [26, 27] for general heterogeneous catalysis and [42] for PEM fuel cells, the homogenization of the initial boundary value problem (1),(2) corresponds to the following system

div(u∗ ) = 0
x∈Ω


 u∗ = − K ∇p∗
x
∈Ω
µ
(3)
∗
∂X
∗
∗
∗
∗
 ∂t + u · ∇X = λdiv(D∇X ) + |∂ω|R(X ) x ∈ Ω

 ∗
X = Xinit
t = 0, x ∈ Ω,
and the boundary conditions


u∗ = uD
X ∗ = XD
t > 0, x ∈ ΓD ,
u∗ · n = 0 ∇X ∗ · n = 0 t > 0, x ∈ ΓN .

(4)

The functions u∗ , p∗ , X ∗ correspond to the homogenized velocity, pressure and O2 concentration, respectively. D is the effective porous media diffusion tensor, |∂ω| is the perimeter of the micro-fuel-tube
ω scaled to the unit cell and K is the permeability tensor.
The definition of the tensors K and D stems from the so-called cell problems [26], namely
Z
Z

 

Kij =
∇ui (y) : ∇uj (y)dy, Dik =
ei + ∇πi · ek + ∇πk dy
Y \ω

(5)

Y \ω

where the tensorial product “ :” represents the twice-contracted tensorial product, i.e. for A, B ∈ Rd×d
X
A:B=
Ak` Bk` ,
k,`

and the respective cell problems read

∇pi − ∆ui = ei , x ∈ Y \ω,



div(ui ) = 0, x ∈ Y \ω,
ui = 0, x ∈ ∂ω,



y → pi (y), ui (y) Y-periodic,





−div(∇πj + ej ) = 0, x ∈ Y \ω,
∇πj · n = −ej · n, x ∈ ∂ω,
y → πj (y) Y-periodic,

(6)

with (ei )1≤i≤d being the canonical basis of Rd .

2.3

Optimization problem

The topology optimization problem consists in finding the optimal lay-out ω ∗ ⊂ Y , representing the
scaled shape of the fuel tubes inside the unit-cell, which solves

max |∂ω|


ω⊂Y


s.t.
(7)

|ω| ≥ Cf |∂ω|



tr(K)
≥ kmin ,
d
where the tensor K is defined according to (5) and tr() is the trace operator. The first constraint
represents a lower bound Cf ≥ 0 of the fuel tube hydraulic diameter, which prevents
a drastic pressure
p
drop and a “too oscillating” boundary. The constant Cf must satisfy Cf ≤ 1/4π, according to the
isoperimetric inequality. The second constraint involves a lower bound kmin of the trace of K (which is
intended to give a measure of the permeability [22]) avoiding a high pressure drop in the cathode. In
other words, we try to find the shape ω with the largest perimeter (so it maximizes the factor on the
electro-chemical reaction term R of the homogenized transport equation in (3)) such that two constraints
of pressure drop inside and outside the tubes are fulfilled.
Let ω0 ⊂ Y be a regular fixed open set. Then the following result yields the existence of an optimal
solution of (7)

7

Proposition 2.1. Denote as Uad the collection of open subsets of Y close to ω0 in the sense of the
pseudo-distance


kT − IdkW 2,∞ + T −1 − Id W 2,∞ ,
dD (ω0 , ω) =
inf
T ∈D|T (ω0 )=ω

where D is the set of diffeomorphism

D = T such that (T − Id) ∈ W 2,∞ (Rd ; Rd ), (T −1 − Id) ∈ W 2,∞ (Rd ; Rd ) .
Furthermore, chose in particular the family of applications T = Id + θ, where θ ∈ W 2,∞ (Rd ; Rd ) and
kθkW 2,∞ < 1. Then problem (7) admits at least one optimal solution in Uad .
Proof. It is a classical result that the set Uad is compact w.r.t. the above topology, introduced by
F.Murat and J.Simon [35]. Furthermore since kθkW 2,∞ < 1, any minimizing sequence of (7) has a
convergent subsequence in W 1,∞ (Rd ; Rd ), due to the Sobolev embedding theorem. Finally the existence
of a solution of (7) stems from the continuity of the perimeter, the area and the permeability tensor K
w.r.t. kθkW 1,∞ , thanks to their respective shape differentiability. The concept of the shape derivative
will be recalled in the next section.

3

Shape sensitivity analysis

In this section we briefly recall the definition of the shape derivative and the gradients of the functionals
involved in (7) (the proofs can be found in [41]).
Shape differentiation is a classical topic that goes back to Hadamard [2, 24, 49]. Let the overall
domain Ω ⊂ Rd be fixed and bounded. Let ω ⊂ Ω be a smooth open subset which is variable. Indeed,
we consider variations of the type

Id + θ (ω) := {x + θ(x) for x ∈ ω} ,
with θ ∈ W 1,∞ (Ω; Rd ) such that kθkW 1,∞ (Ω;R2 ) < 1 and tangential on ∂Ω (i.e., θ · n = 0 on ∂Ω ; this
last condition ensures that Ω = (Id + θ)Ω). It is well known that, for sufficiently small θ, (Id + θ) is a
diffeomorphism in Ω.
Definition 3.1. The shape derivative of a function J(ω) is defined as the Fréchet derivative in W 1,∞ (Ω; Rd )
at 0 of the application θ → J Id + θ ω , i.e.
J

 
Id + θ ω = J(ω) + J 0 (ω)(θ) + o(θ)

with

lim

θ→0

|o(θ)|
=0,
kθkW 1,∞

where J 0 (ω) is a continuous linear form on W 1,∞ (Ω; Rd ).
Now we apply the above definition in the context of our cell problem (6).
Lemma 3.2. Let ω be a smooth subset of Y . Define
Z
Jsurf (ω) = |∂ω| =
ds and
∂ω

Jvol (ω) = |ω| =

Then the respective shape gradients at 0 in the direction θ read
Z
Z
0
0
Jsurf
(ω)(θ) =
θ · nH ds, Jvol
(ω)(θ) =
∂ω

∂ω

Z
dx.
ω

θ · n ds,

where H is the mean curvature of ∂ω defined by H = div(n).
Proposition 3.3. Define the spaces

1
H0,#
(Y \ω)d = v ∈ H 1 (Y \ω)d : v|∂ω = 0, y → v(y) is Y −periodic ,
(
)
Z
L20,# (Y \ω) =

q ∈ L2 (Y \ω) :

Y \ω

8

qdx = 0, y → q(y) is Y −periodic .

1
Let ui ∈ H0,#
(Y \ω)d and pi ∈ L20,# (Y \ω), i = 1, ..., d be the collection of solutions of the variational
formulation of the first cell problem in (6)
Z
Z


1
∇ui : ∇vi − div(vi )pi − div(ui )qi dx =
ei · vi dx, ∀(vi , qi ) ∈ (H0,#
(Y \ω)d , L20,# (Y \ω)). (8)
Y \ω

Y \ω

Define the cost function (which does not depend on the pressure)
Z
J(ω) =
j(x, u, ∇u)dx, where u = (ui )i=1...d ,
Y \ω

1
and the collection of adjoint states (Ui , Pi ) ∈ (H0,#
(Y \ω)d , L20,# (Y \ω)), such that
Z
Z



 ∂j
∂j
(x, ui , ∇ui ) · vi +
(x, ui , ∇ui ) : ∇vi dx,
∇Ui : ∇vi − div(vi )Pi − div(Ui )qi dx = −
∂∇ui
Y \ω
Y \ω ∂ui
(9)
1
for all (vi , qi ) ∈ (H0,#
(Y \ω)2 , L20,# (Y \ω)). Then J(ω) is shape differentiable at 0 in the direction θ and
the shape derivative reads

J 0 (ω)(θ) =

d 
n
X
∂j
∂ui o
∂Ui ∂ui
θ · nds.
j(x, u, ∇u) −
·
+
·n·
∂n ∂n
∂∇ui
∂n
∂ω
i=1

Z

Corollary 3.4. According to (5)
d

tr(K)
1X
1
=
Kii =
d
d i=1
d

Z

d
X

Y \ω i=1

|∇ui |2 dx,

so the shape derivative reads
d

1
1X
tr(K)0 (ω)(θ) = −
d
d i=1

Z

 ∂u 
∂ui ∂Ui 
i 2
+
·
θ · nds,
∂n
∂n ∂n
∂ω

with Ui solution of the adjoint problem
Z


∇Ui : ∇vi − div(vi )Pi − div(Ui )qi dx = −2

Y \ω

4

Z

Y \ω

1
∇ui : ∇vi , ∀(vi , qi ) ∈ (H0,#
(Y \ω)2 , L20,# (Y \ω))

Numerical results

Now we detail the numerical solution of problem (7) via the level-set method of Osher and Sethian [37]
for topology optimization. At each iteration of the algorithm, the shape of the tube ω in the unit cell Y
is parametrized via a level-set function

 ψ(x) = 0 for x ∈ ∂ω,
ψ(x) < 0 for x ∈ ω,

ψ(x) > 0 for x ∈ Y \ω,
defined on a fixed mesh in an Eulerian framework.
During the optimization process, the shape ω(t) is going to evolve according to a fictitious time t ∈ R+ ,
which corresponds to the descent stepping. As it is well known, the evolution of the level set function ψ
is governed by a Hamilton-Jacobi equation
∂ψ
+ V|∇ψ| = 0,
∂t

(10)

where V is the normal velocity. Equation (10) is posed in the whole domain Y , and not only on the
interface ∂ω. The main point in using a level set method is that it replaces the Lagrangian evolution of
the boundary ∂ω by the Eulerian solution of a transport equation in the whole fixed domain Y , easily
allowing topology changes. In practice the Hamilton-Jacobi equation (10) is solved by an explicit second
9

order upwind scheme (see e.g. [43]) on a fixed Cartesian grid. The boundary conditions for ψ are of
Neumann type. Since this scheme is explicit in time, its time stepping must satisfy a CFL condition.
Moreover, in order to regularize the level set function (which may become too flat or too steep), we
reinitialize it periodically by solving another Hamilton-Jacobi equation which admits as a stationary
solution the signed distance to the initial interface [43].
The choice of the normal velocity V in (10) is based on the shape derivatives computed in the previous
section accordingly to [4]. Generally speaking, the shape derivative of a functional J(ω) in the direction
θ ∈ W 1,∞ (Y ; Rd ) reads
Z
J 0 (ω)(θ) =

∂ω

T θ · nds,

(11)

where the integrand T (x) depends on the solution of the cell problem and possibly of some adjoint
equation. Since only the normal component of θ plays a role in (11), a descent direction for J is just a
vector field θ satisfying
Z
θ = Vn

and

J 0 (ω)(θ) =

∂ω

T Vds ≤ 0.

(12)

To ensure the decrease of J, the simplest choice is V = −T . However, T is a priori defined only on the
interfaces ∂ω while V must be defined in the entire domain Y . Therefore, the choice V = −T is implicitly
depending on some extension process of T . If such an extension is not obvious or if we want to regularize
the velocity fields, there is an alternate choice based on a different underlying scalar product (see e.g. [4],
[18]).
Since the interface ∂ω does not necessarily match with the mesh, in order to simulate an impenetrable
interface in (6), an artificial permeability (or porosity) term is added on the left hand sides of equations
(8) and (9), respectively
Z
Z
5
5
(u
·
v
)dx,
(Ui · vi )dx,
(13)
i
i
2
2
Y 2ρ
Y 2ρ
where 2ρ2 /5 represents the permeability of the medium [12, 20]. For topology optimization purposes,
ρ is a regular approximation of H(ψ), where H is the Heaviside function and min ρ(x) = δ > 0 [7, 10].
x∈Y

Typically this value is set up to δ = 0.001. We remark that thanks to the introduction of (13), all
integrals are thus calculated in Y .

Figure 8: Example of the 2D stokes solutions u1 , u2 (left and right respectively) of the cell problem for a circular
domain in (8). The black zone represents the air phase (transporting the oxygen) meanwhile the white
one the fuel phase. A solid interface lays between them. We remark that inside the fuel phase the
velocity of the air is almost zero, due to (13).

For the optimization algorithm we use an augmented Lagrangian method



2
tr(K)
µ1
µ2
tr(K)
2
L(ω, `, µ) = |∂ω| − `1 (Cf |∂ω| − |ω|) − `2 kmin −
+
(Cf |∂ω| − |ω|) +
kmin −
,
d
2
2
d
where ` = (`i )i=1,2 and µ = (µi )i=1,2 are lagrange multipliers and penalty parameters for the constraints.
The Lagrange multipliers are updated at each iteration n according to the optimality condition


tr(K)
n+1
n
n
`n+1
=
`
−
µ
(C
|∂ω
|
−
|ω
|)
and
`
=
`
−
µ
k
−
.
1
f
n
n
2
min
1
2
1
2
d
10

The penalty parameters are augmented every 5 iterations. With such an algorithm the constraints are
enforced only at convergence.
Fig. 9 shows different optimal periodic layouts starting from an intuitive circular arrangement. The
radius of the circle in Y was set up to r = 0.3, meanwhile the coefficients Cf = 0.15 and kmin = 0.011
were calculated so as both constraints in (7) were active within this layout. The first optimal design
(second row Fig. 9) corresponds to the optimal solution of the algorithm for the foregoing parameters.
Then the second design (third row Fig. 9) was established by reducing Cf and kmin to the half. The
perimeter gain w.r.t. the circular layout is 7% for the former design and 60% for the latter one.

Figure 9: Two micro-tubular optimal designs for different values of kmin and Cf starting from a circular layout.
The base cell results are on the left and the corresponding periodic structures on the right.

5

Summary and Outlook

The above results demonstrate the capabilities of topology optimization via a level-set method to
enhance the design of micro-tubular fuel cells for the aeronautic industry. Thanks the application of an
inverse homogenization technique and the level-set method, periodic optimal micro-tubular fuel cells with
a sharp contour can easily be designed and then manufactured by 3D printing. The foregoing study thus
suggests a promising use of these technologies in the future computer aided design of fuel cells.
Acknowledgments: The author would like to thank Airbus Group for its financial support in the
framework of his PhD thesis, Ch. Nespoulous and E. Moullet of Airbus Innovations for their important
contribution through the joint work developed during E. Mullet’s internship [34], and also G. Allaire of

11

Centre de Mathématiques Appliquées of École Polytechnique for many helpful suggestions after kindly
reading this manuscript.

12

References
[1] G. Allaire. Shape optimization by the homogenization method, volume 146. Springer, 2002.
[2] G. Allaire. Conception optimale de structures, volume 58. Springer Verlag, 2006.
[3] G. Allaire, E. Bonnetier, G. Francfort, and F. Jouve. Shape optimization by the homogenization
method. Numerische Mathematik, 76(1):27–68, 1997.
[4] G. Allaire, F. Jouve, and A.-M. Toader. Structural optimization using sensitivity analysis and a
level-set method. Journal of computational physics, 194(1):363–393, 2004.
[5] G. Allaire and R. V. Kohn. Optimal design for minimum weight and compliance in plane stress
using extremal microstructures. European journal of mechanics. A. Solids, 12(6):839–878, 1993.
[6] J.-L. Auriault, C. Boutin, and C. Geindreau. Homogenization of coupled phenomena in heterogenous
media, volume 149. John Wiley & Sons, 2010.
[7] M. Bendsøe. Methods for optimization of structural topology, shape and material. Springer Verlag,
New York, 1995.
[8] M. Bendsøe and O. Sigmund. Topology optimization: theory, methods and applications. Springer,
2004.
[9] M. P. Bendsøe and N. Kikuchi. Generating optimal topologies in structural design using a homogenization method. Computer methods in applied mechanics and engineering, 71(2):197–224, 1988.
[10] T. Borrvall and J. Petersson. Topology optimization of fluids in stokes flow. International Journal
for Numerical Methods in Fluids, 41(1):77–107, 2003.
[11] R. Bove and S. Ubertini. Modeling solid oxide fuel cells: methods, procedures and techniques, volume 1. Springer Science+ Business Media, 2008.
[12] H. Brinkman. A calculation of the viscosity and the sedimentation constant for solutions of large
chain molecules taking into account the hampered flow of the solvent through these molecules.
Physica, 13(8):447–448, 1947.
[13] M. Burger. A framework for the construction of level set methods for shape optimization and
reconstruction. Interfaces and Free boundaries, 5(3):301–330, 2003.
[14] J. E. Cadman, S. Zhou, Y. Chen, and Q. Li. On design of multi-functional microstructural materials.
Journal of Materials Science, 48(1):51–66, 2013.
[15] A. Cherkaev. Variational methods for structural optimization, volume 140. Springer, 2000.
[16] P. Costamagna, A. Selimovic, M. Del Borghi, and G. Agnew. Electrochemical model of the integrated
planar solid oxide fuel cell (ip-sofc). Chemical Engineering Journal, 102(1):61–69, 2004.
[17] V. A. Danilov and M. O. Tade. A cfd-based model of a planar sofc for anode flow field design.
International Journal of Hydrogen Energy, 34(21):8998–9006, 2009.
[18] F. de Gournay. Velocity extension for the level-set method and multiple eigenvalues in shape optimization. SIAM J. Control Optim., 45(1):343–367, 2006.
[19] R. De la Torre Garcı́a. Production of Micro-Tubular Solid Oxide Fuel Cells. PhD thesis, University of
Trento, http://eprints-phd.biblio.unitn.it/541/1/PRODUCTION OF MICRO-TUBULAR SOLID
OXIDE FUEL CELLS.pdf, 2011.
[20] L. Durlofsky and J. Brady. Analysis of the brinkman equation as a model for flow in porous media.
Physics of Fluids, 30(11):3329–3341, 1987.
[21] B. Gou, W. K. Na, and B. Diong. Fuel cells: modeling, control, and applications. CRC press, 2009.
[22] J. K. Guest and J. H. Prévost. Design of maximum permeability material structures. Computer
Methods in Applied Mechanics and Engineering, 196(4):1006–1017, 2007.

13

[23] J. Haslinger and J. Dvorak. Optimum composite material design. RAIRO-M2AN Modelisation Math
et Analyse Numerique-Mathem Modell Numerical Analysis, 29(6):657–686, 1995.
[24] A. Henrot and M. Pierre. Variation et Optimisation des formes, Une analyse géométrique, volume 48
of Mathématiques et Applications. Springer, 2005.
[25] M. F. Hordeski. Hydrogen and Fuel Cells: Advances in Transportation and Power. The Fairmont
Press, Inc., 2009.
[26] U. Hornung. Homogenization and porous media, volume 6. Springer, 1997.
[27] U. Hornung and W. Jäger. Diffusion, convection, adsorption, and reaction of chemicals in porous
media. Journal of differential equations, 92(2):199–225, 1991.
[28] S. Hosseini, K. Ahmed, and M. O. Tadé. Cfd model of a methane fuelled single cell sofc stack for
analysing the combined effects of macro/micro structural parameters. Journal of Power Sources,
2013.
[29] M. Hussain, X. Li, and I. Dincer. Mathematical modeling of planar solid oxide fuel cells. Journal of
Power Sources, 161(2):1012–1022, 2006.
[30] R. J. Kee, H. Zhu, R. J. Braun, and T. L. Vincent. Chapter 6 - modeling the steady-state and
dynamic characteristics of solid-oxide fuel cells. In K. Sundmacher, editor, Fuel Cell Engineering,
volume 41 of Advances in Chemical Engineering, pages 331 – 381. Academic Press, 2012.
[31] B. Kenney and K. Karan. Mathematical micro-model of a solid oxide fuel cell composite cathode.
Proceedings hydrogen and fuel cells, pages 1–11, 2004.
[32] R. V. Kohn and G. Strang. Optimal design and relaxation of variational problems; Parts I, II, III.
Communications on Pure and Applied Mathematics, 39(1):113–137, 1986.
[33] A. Kumar and R. Reddy. Modeling of polymer electrolyte membrane fuel cell with metal foam in
the flow-field of the bipolar/end plates. Journal of power sources, 114(1):54–62, 2003.
[34] E. Moullet. Optimisation du design d’une pile à combustible. Technical report, École de Mines,
EADS, 2013.
[35] F. Murat and J. Simon. Quelques résultats sur le contrôle par un domaine géométrique. VI Laboratoire d’Analyse Numérique, 1974.
[36] National Institute of Advanced Industrial Science and Technology. https://www.aist:go:jp/index,
2009.
[37] S. Osher and J. Sethian. Fronts propagating with curvature-dependent speed: algorithms based on
hamilton-jacobi formulations. Journal of computational physics, 79(1):12–49, 1988.
[38] J. Palsson, A. Selimovic, and L. Sjunnesson. Combined solid oxide fuel cell and gas turbine systems
for efficient power and heat generation. Journal of Power Sources, 86(1):442–448, 2000.
[39] B. Roth and R. Giffin III. Fuel cell hybrid propulsion challenges and opportunities for commercial
aviation. AIAA Paper, 6537:2010, 2010.
[40] S. Samuelsen. Fuel cell/gas turbine hybrid systems. ASME International Gas Turbine Institute,
2004.
[41] S. Schmidt and V. Schulz. Shape derivatives for general objective functions and the incompressible
navier-stokes equations. Control and Cybernetics, 39(3):677–713, 2010.
[42] M. Schmuck and P. Berg. Homogenization of a catalyst layer model for periodically distributed pore
geometries in pem fuel cells. Applied Mathematics Research eXpress, 2013(1):57–78, 2013.
[43] J. Sethian. Level set methods and fast marching methods: evolving interfaces in computational
geometry, fluid mechanics, computer vision, and materials science, volume 3. Cambridge university
press, 1999.

14

[44] J. A. Sethian and A. Wiegmann. Structural boundary design via level set and immersed interface
methods. Journal of computational physics, 163(2):489–528, 2000.
[45] O. Sigmund. Design of material structures using topology optimization. PhD thesis, Technical
University of Denmark Denmark, 1994.
[46] O. Sigmund. Materials with prescribed constitutive parameters: an inverse homogenization problem.
International Journal of Solids and Structures, 31(17):2313–2329, 1994.
[47] O. Sigmund and S. Torquato. Design of materials with extreme thermal expansion using a three-phase
topology optimization method. Journal of the Mechanics and Physics of Solids, 45(6):1037–1067,
1997.
[48] S. Singhal. Advances in solid oxide fuel cell technology. Solid state ionics, 135(1):305–313, 2000.
[49] J. Sokolowski and J.-P. Zolésio. Introduction to shape optimization. Springer, 1992.
[50] C. J. Steffen Jr, J. E. Freeh, and L. M. Larosiliere. Solid oxide fuel cell/gas turbine hybrid cycle
technology for auxiliary aerospace power. Proceedings of the ASME Turbo Expo, Reno-Tahoe, United
States, 2005.
[51] L. Tartar. The general theory of homogenization: a personalized introduction, volume 7. Springer,
2009.
[52] M. Y. Wang, X. Wang, and D. Guo. A level set method for structural topology optimization.
Computer methods in applied mechanics and engineering, 192(1):227–246, 2003.
[53] S. Zhou and Q. Li. Computational design of multi-phase microstructural materials for extremal
conductivity. Computational Materials Science, 43(3):549–564, 2008.

15

