Computational design of flow fields for vanadium redox flow batteries via topology optimization

Journal of Energy Storage 26 (2019) 100990

Contents lists available at ScienceDirect

Journal of Energy Storage journal homepage: www.elsevier.com/locate/est

Computational design of flow fields for vanadium redox flow batteries via topology optimization

T

Chih-Hsiang Chen, Kentaro Yaji⁎, Shintaro Yamasaki, Shohji Tsushima, Kikuo Fujita Department of Mechanical Engineering, Graduate School of Engineering, Osaka University, 2-1, Yamadaoka, Suita, Osaka 565-0871, Japan

ARTICLE INFO

ABSTRACT

Keywords: Redox flow battery Flow field design Topology optimization Mass transfer effect

Vanadium redox flow battery (VRFB) is a rechargeable battery, which has attracted attention as a next-generation electrochemical energy storage system. It is widely known that achieving high performance in terms of power density is critical for the commercialization of VRFBs. The aim of this paper is to propose a computational design approach for automatically generating an optimized flow field design of a VRFB to achieve high performance without relying on the designer’s intuition. To realize this, we focus on generating a freeform configuration of the flow field in a VRFB via topology optimization, which is known as a powerful design tool that is based on numerical optimization. In this study, to improve the mass transfer effect in a VRFB, we formulate the topology optimization problem as a maximization problem of the electrode surface concentration in the negative electrode during the charging process. We demonstrate through numerical investigation that a topology-optimized flow field can be obtained. As a result, it is revealed that the interdigitated flow field is an optimal flow field of the VRFB under the investigated operating conditions.

1. Introduction Over the past years, there has been growing interest in renewable energy since the increase of carbon emission poses a great threat to the environment [1]. Despite the significance of renewable energy, the intermittent characteristics of sources, such as solar, wind or water, is a fatal drawback for renewable energy, which leads to the increased uncertainty in the supply of electricity [2]. To address this issue, one promising solution is to regulate the power delivery via energy storage technology. Among the energy storage systems, vanadium redox flow batteries (VRFBs) attract a lot of attention due to the advantageous features: scalability, low cost and long cycle life [3–5]. However, achieving high performance in terms of power density is a critical issue for cost-effectiveness of VRFBs [6]. For a given set of electrochemical conditions, it is known that the performance of battery devices depends on various factors such as material properties of the electrode [7–9] and the type of catalyst [10–13]. In addition, the performance of VRFBs depends on mass transfer losses, where the mass transfer effect can be ameliorated by different flow fields [14]. Aaron et al. [15] first reported the effectiveness of using a serpentine flow field design in a VRFB. The results revealed that the limiting current density and peak power density are improved in comparison with a flow-through type VRFB. Xu et al. [16]

⁎

numerically studied the performance of VRFBs with several types of flow fields, namely, flow through, serpentine, and parallel flow fields. Messaggi et al. [17] empirically and numerically compared the serpentine and interdigitated flow field designs, and revealed that the serpentine flow field allows better performances under all the investigated operating conditions. On the other hand, several research groups have demonstrated that interdigitated flow fields can further improve the performance of VRFBs owing to the enhanced mass transfer effect [9,18,19]. From these results, it is clear that the appropriate type of flow field in VRFBs varies depending on the operating condition. Hence, it is required that the appropriate flow field is designed for each operating condition. To date, although various types of flow fields have been proposed to improve the performances of VRFBs, the design of the flow field remains laborious owing to the complicated physical and chemical mechanisms. To overcome this issue, numerical optimization can be used as a promising strategy for automatically generating optimal designs of VRFBs. Li et al. [20] studied the optimal charging conditions of VRFBs with respect to input current and electrolyte flow rate. Yue et al. [21] proposed a novel trapezoid flow battery and identified the optimal parameter settings of its geometrical shape for improving the mass transfer effect. You et al. [22] investigated the optimal number and size of channels in an interdigitated flow field. The degree of design freedom

Corresponding author E-mail address: [email protected] (K. Yaji).

https://doi.org/10.1016/j.est.2019.100990 Received 25 July 2019; Received in revised form 20 September 2019; Accepted 28 September 2019 2352-152X/ © 2019 Elsevier Ltd. All rights reserved.

Journal of Energy Storage 26 (2019) 100990

C.-H. Chen, et al.

Nomenclature A a c D

df F F I i0 J j K k KCK n p q R s T U u U0 x

anodic transfer coefficient cathodic transfer coefficient design-dependent inverse permeability (m 2 ) porosity overpotential (V) conductivity of the electrolyte (S m 1) dynamic viscosity (Pa s) material domain (m3) domain of the porous electrode (m3) electrolyte potential (V) electrode potential (V) design field conductivity of the electrode (S m 1) characteristic function

a c

area (m2) specific surface area (m 1) concentration of vanadium species (mol m 3 ) diffusivity (m2 s 1) fixed design domain (m3) mean fiber diameter (m) Faraday constant (C mol 1) body force (N m 3 ) applied current (A) exchange current density (A m 2 ) objective function (mol m 3 ) transfer current density (A m 2 ) permeability (m2) reaction rate constant (m s 1) Carman–Kozeny constant unit normal vector pressure (Pa) tuning parameter for convexity of αρ gas constant (J mol 1 K 1) source term of vanadium species (mol m 3 s 1) temperature (K) open-circuit potential (V) velocity (m s 1) equilibrium potential (V) position (m)

αρ ϵ η e

μ Ω

p

e

ρ

s

s

χΩ Subscripts

V 2+ V 3+

V 2 + ion V 3+ ion

Superscripts eff fic in out s

effective properties fictitious properties inlet outlet surface of the electrode

Greek symbols α

inverse permeability (m 2 )

was limited in these previous works because they optimized a small number of design variables based on parametric optimization. Therefore, the optimal type of flow field under an operating condition could not be identified in general, as the type of flow field had to be determined in advance. Topology optimization is an optimization approach for realizing a freeform design without a promising initial guess [23]. It expresses an optimization problem as a material distribution problem in a given design domain and then derives a freeform configuration. In topology optimization, design variables individually correspond to the material existing at each position for realizing the optimization of not only the shape but also the topology [24]. Due to its high degree of design freedom, topology optimization has been applied to various problems, e.g., structural stiffness problems [23,25], thermal problems [26,27], and electromagnetic problems [28–30]; furthermore, applications to practical device designs have also been attracted attention in micro actuator design [31], fuel cell design [32–34] and so on. For fluid problems, Borrvall and Petersson [35] first proposed a topology optimization method to minimize power dissipation in Stokes flow, and this has been expanded to laminar Navier–Stokes flow problems [36–38] and turbulence problems [39,40]. The fluid topology optimization has been applied to multiphysics problems such as fluid-structure interaction problems [41,42], forced convection problems [43–46], natural convection problems [47–49] and turbulent heat transfer problems [50,51]. Recently, Yaji et al. [52] proposed a topology optimization method for the design of flow fields in VRFBs. In their approach, instead of the formula of practical electrochemical reactions, a simplified formula is introduced as a two-dimensional model. They indicated that topology optimization can be a promising approach for realizing the freeform design of the flow field in a VRFB, and clarified that the optimized configurations tend to be of the interdigitated flow field type under the

investigated operating conditions. However, this result may not be always correct owing to the use of the two-dimensional pseudo-analysis model for simplicity. As a more comprehensive study, this paper aims to construct topology optimization for flow fields in VRFBs based on a three-dimensional model incorporating with electrochemical kinetics. Referring to the models presented by several researchers [53–55], a three-dimensional analysis model of a negative electrode in a VRFB is introduced. We demonstrate through a numerical example that the proposed approach enables the automatic generation of a novel flow field configuration and discuss the efficacy of the proposed approach. 2. Mathematical model 2.1. Model assumptions A schematic diagram of a typical redox flow battery is shown in

Fig. 1. Schematic diagram of a vanadium redox flow battery. 2

Journal of Energy Storage 26 (2019) 100990

C.-H. Chen, et al.

Fig. 1. The positive and negative electrodes are separated by the ion exchange membrane, which only allows protons to penetrate. We suppose that the electrodes compose of the carbon fiber electrode and flow channel. Note that the use of flow channel enables the reduction of pressure loss in comparison with the case of only using the carbon fiber electrode [16]. When the electrolyte stored in tanks circulates through the positive and negative electrode separately by pumps, the electric energy is released or stored by the electrochemical reactions in the electrodes. The main reactions can be described as follows.

Positive electrode: VO2 + + H2 O

Negative electrode: V 3 + + e

e

Charge Discharge

Charge Discharge

V2+

VO+2 + 2H+

(1)

j = i0

(2)

=

(3)

)

2u

+ F = 0.

c V2 +

exp

aF RT

,

(9) (10)

e

(11)

U,

c 3+ RT ln V , F c V2 +

(12)

p = pin on the inlet,

(13)

p = pout on the outlet,

(14) (15)

u = 0 on the remaining boundaries,

is the given pressure value on the inlet, and is the given where pressure value on the outlet. The boundary conditions for species conservation are shown below:

pin

ci = ciin on the inlet, (6)

Dieff ci·n = 0 on the remaining boundaries,

where df is the fiber diameter, ϵ is the porosity of electrodes, KCK is the Carman–Kozeny constant described by the characteristic of the fibrous material. In addition, instead of the electrolyte flow described by Stokes equation and Darcy’s law respectively, Brinkman equation that combines Stokes equation with Darcy’s law can be used to describe the electrolyte flow in the mixture of different porous medium generally and is formulated as

p+µ

RT

Fig. 2 shows the schematic diagram of the analysis domain and boundary settings. The pressure conditions are imposed on the inlet and the outlet, and the no-slip condition is applied on the remaining outer boundaries. The expressions are shown below:

(5)

, 2

c Vs2 +

2.3. Boundary conditions

where K is the permeability coefficient, which can be described by the Carman–Kozeny equation [56] as follows:

d f2 3 16K CK (1

cF

where U0 is the equilibrium potential.

where μ is the viscosity of electrolyte flow, p is the pressure, and u is the velocity, respectively. In addition to the flow channel, the porous electrode is also permeated with the electrolyte flow and the velocity in the electrode can be expressed by Darcy’s law:

K=

s

U = U0 +

(4)

p,

c V3 +

exp

where U is the open-circuit potential in the negative electrode, which can be estimated by Nernst equation as follows:

Based on the assumptions presented in the previous section, the governing equations incorporated in the numerical model are briefly introduced here. The detailed equations can be seen in the previous work [55]. The electrolyte flow passing through the flow channel in VRFBs is governed by Stokes equation and continuity equation as follows:

µ u= K

c Vs3 +

V 3 +},

Here, i0 is the exchange current density, η is the overpotential, c and a are the cathodic and anodic transfer coefficients, respectively, k is the reaction rate constant, R is the gas constant, T is the temperature, and cis is the species concentration at the surface of the negative electrode. The overpotential in (9) is the difference between the electrode potential s and the electrolyte potential e , which can be expressed as follows:

2.2. Governing equations

· u = 0,

(8)

= si,

i 0 = Fk (c V2 +) c (c V3 +) a .

process.

= 0,

i

where ci is the concentration of vanadium species i and si is the source term of species i due to electrochemical reactions, s V2 + = aj/ F and s V3 + = aj/ F . Here, a is the specific area of the electrode, j is the transfer current density, and F is the Faraday constant. In addition, Dieff = 1.5Di , where Di is a diffusion coefficient of the species i, is an effective diffusion coefficient described by the Bruggemann correction. The transfer current density j, which originates from electrochemical reactions, can be described using Butler–Volmer equation, as follows:

• The electrolyte flow is treated as stationary, incompressible and isothermal Stokes flow. • The dilute-solution approximation is used. • The side reactions are neglected in electrochemical reactions. • The migration phenomenon is ignored in the species transport

2u

2c

{V 2 +,

In this study, only the negative electrode is considered for reducing the computation cost, and some basic assumptions are made for simplification, as follows [54]:

p+µ

Dieff

u· ci

pout

(16) (17)

(7)

u is the body force, where α is the inverse-permeability, Here, F = which is defined as = µ /K in the porous medium, whereas = 0 in the pure fluid domain. With the electrolyte flow containing vanadium species, the species transport needs to be considered in the numerical model, which can be expressed as follows:

Fig. 2. Schematic diagram of the analysis domain and boundary settings. 3

Journal of Energy Storage 26 (2019) 100990

C.-H. Chen, et al.

where n is the unit normal vector. For charge conservation, the VRFB is assumed to be operated in galvanostatic situation. During charging process, the flux conditions in the negative electrode can be described as follows:

s

eff s

s ·n

eff e

e ·n =

=

I on the wall 1, A

electrode surface. During charging process, the optimization problem can be defined as a maximization problem of average concentration of oxidized reactants at the electrode surface in a negative electrode. The expressions of this optimization problem are shown as follows:

maximize J =

(18)

I on the wall 2, A

subject to 0

(20)

Here, is the effective conductivity of the electrode, is the effective conductivity of the electrolyte, I is the applied current, and A is the electrode surface area. Note that the remaining boundary conditions for s and e are the no-flux conditions. eff s

eff e

3. Topology optimization for flow fields 3.1. Concepts of topology optimization for fluid problems

1 if x 0 if x

(21)

where x is the position in , and Ω expresses the material domain in . In this expression, (x) = 1 and (x) = 0 represent the material and void at x, respectively. Since the characteristic function χΩ is a discontinuous function, topology optimization problems typically require relaxation treatments. As the popular and simple way for relaxing topology optimization problems, the density approach [24] used in this study replaces the characteristic function with the continuous function, 0 (x) 1. To determine which points in should be fluid or solid, based on the previous research dealing with fluid topology optimization [35], the body force F is redefined using the fictitious body force, Ffic, as follows:

F fic =

fic

u

with

fic

=

q (1

) +q

fic,

1 for

d , x

(23)

.

The governing equations are solved by using the package COMSOL Multiphysics ® (version 5.4), which is based on the finite element method. The CFD module of COMSOL is used to calculate the fluid velocity and pressure in Eqs. (3)–(7). The concentration of vanadium species in Eq. (8) is calculated by using the Transport of Diluted Species, and the potential fields of the electrode and electrolyte are calculated by using the Weak Form PDE. The relative tolerance is set to 1.0 × 10 4 for all calculations in this study. The analysis model was validated by comparing with the result of the previous work [55]. The optimization algorithm is constructed on the basis of mathematical programming and is structured as follows:

, ,

(x)

p

3.3. Numerical implementation

Topology optimization aims to obtain the improved structural design in a specific domain with a given objective function and constraints. The main idea is to formulate a topology optimization problem as a material distribution problem, where the expression for material distribution in a fixed design domain is defined as follows:

(x) =

c Vs3 + d

Here, p is the domain of the porous electrode. The detailed expression for c Vs3 + is found in the previous work [55]. It should be noted that the overpotential η or the transfer current density j cannot be directly used as the objective function J for achieving sufficient performance under the investigated operating conditions in our study because of the numerical instability. This is because, in general, topology optimization problems cannot avoid multimodality, which means that there are many local minima in the solution space, as in nonlinear optimization problems. Note that the appropriate determination of the problem formulation should be performed in advance [57].

(19)

= 0 on the interface.

p

Step 1 The design variables are initialized, and all parameters shown in Table 1 are set. Step 2 The objective function J in (23) is evaluated by solving the governing equation via the finite element method. Step 3 If the objective function is converged, the iteration will terminate. Otherwise, the sensitivities—gradient of objective function Table 1 Parameter settings for electrode, electrolyte properties.

(22)

where αfic is the fictitious inverse-permeability used for expressing the as with the previous work [35], and q is a tuning solid domain parameter for fic . In addition, = 1 expresses the fluid domain with fic = 0, and = 0 expresses the solid domain with fic = fic . In the solid domain, since the fictitious body force is large enough compared to the fluid domain, it is difficult for fluid to pass through the solid domain. Therefore, according to the sensitivity information of the objective function, the fictitious body force is determined at each point in so that the structural design of flow channel can be obtained. Note that the fictitious inverse permeability, αfic, is different from α. That is, the former is used for expressing the solid domain in the fixed design domain , whereas the latter is used for expressing the porous electrode that is the non-design domain. In this study, q and fic are defined as 0.01 and 5α, respectively.

Parameter

Symbol

Value

Unit

Ref.

Porosity Specific surface area Carbon fiber diameter

ϵ a df

0.929 1.62 × 104

– m m

[58] [58] [58]

Electronic conductivity

Carman–Kozeny constant Electrode length Electrode width Electrode thickness

4.28 0.1 0.1

3 3

c in2 +

4.9 × 10 750

Inlet V 3 + concentration

c in3 +

750

V 2 + diffusion coefficient

D V2 +

2.4 × 10

4

D V3 +

2.4 × 10 7.8

4

Inlet V 2 + concentration

Cathodic transfer coefficient Anodic transfer coefficient Equilibrium potential Temperature Inlet pressure Outlet pressure

Applied current

4

1.76 × 10 5 1.0 × 103

3.0 × 10

μ

Reaction rate

In VRFBs, the mobility of ions in electrolyte is poor compared to the electrons, which means it is more difficult to reach the electrode surface for ions. With the flow channel embedded in a VRFB, the mass transfer effect is improved so that the concentration of reactants at the electrode surface will also increase. Therefore, whether the mass transfer effect is improved can be estimated by the concentration of the reactants at the

K CK L W te

Dynamic viscosity

V 3 + diffusion coefficient Ionic conductivity

3.2. Problem formulation

s

V

V

S m – m m m

[54] [16] [16] [16]

Pa s

[54]

mol m

3

[54]

mol m

3

[54]

m2 s

1

m2 s

1

S m

1

m s – – V K Pa

1

pin

1.7 × 10 7 0.5 0.5 −0.255 298 1.0 × 103

pout I

4.0

A

e

k

c a

U0 T

0

[54]

1

Pa

[55] [55]

Estimated [3]

Assumed Assumed [59] [54] Assumed Assumed [16]

Journal of Energy Storage 26 (2019) 100990

C.-H. Chen, et al.

Fig. 3. Iteration history of topology-optimized flow field expressed as ρ ≥ 0.5.

Fig. 4. Flow field configurations.

compression or injection molding for mass production. On the other hand, 5-axis milling machines or 3D printers can be used for limited production. From the result, it can be confirmed that the proposed approach enables the automatic generation of an optimized flow field design, and therefore it can be said that the interdigitated flow field is an optimal flow field of the VRFB, at least under the operating conditions investigated in this study. To investigate the performance of the optimized flow field design, we compare it with the reference models—parallel, and interdigitated flow fields—as shown in Fig. 4. The geometric parameters of the electrode are presented in Table 1. The width and thickness of the flow channel in the reference flow fields are 3 mm and the interval between the branches of the flow field is 9 mm. Table 2 shows the detailed mesh information for discretizing the flow fields in Fig. 4. Here, the optimized flow field design is recalculated using a body-fitted unstructured mesh built with the User-controlled mesh in the COMSOL settings for comparing with the reference models under similar mesh properties. Note that although the analysis domain is discretized using the uniform hexahedral elements during the optimization process shown in Fig. 3, we confirmed that the discrepancy of the objective function value between the uniform mesh and body-fitted mesh was less than 1%. The comparison results for different flow fields are shown in Figs. 5 and 6, which display the distributions of |u|, p, c V2 +, and c V3 + on the middle plane of the electrode and on the sectional planes, respectively. Note that although the minimum and maximum values of p, c V2 +, and c V3 + in Figs. 5 and 6 are respectively equal with respect to the color bars, the maximum value of |u| is different due to a large difference in the order between the flow channel and porous electrode. As shown in Figs. 5 and 6, the velocity of the electrolyte flow in the parallel flow field is small in the porous electrode because of the difficulty in distributing the electrolyte effectively. The under-rib convection of the parallel flow field is weak, especially on both sides of the electrode. In contrast, the under-rib convection is strong in the interdigitated flow field and topology-optimized flow field. Especially, in the optimized flow field, it can be confirmed that the velocity magnitude of the inner

Table 2 Mesh properties for discretizing the analysis domains of the flow fields in Fig. 4. Parameter

Value

Maximum element size

2.25 × 10

3 4

Maximum element growth rate

6.93 × 10 1.25

Minimum element size

with respect to the design variables—are calculated. Step 4 The design variables are redistributed in the fixed design domain using sequential linear programming, and the iteration will return to the second step. The detailed algorithm can be referred from the previous research dealing with two-dimensional topology optimization for VRFBs [52], as the above optimization procedure is same as that reported in the research. 4. Results and discussion We now investigate the efficacy of the proposed approach through a numerical example. The main objective of this study is to investigate whether the proposed approach can realize the automatic generation of the flow field design of a VRFB. Moreover, we compare the performance of the optimized result with those of the reference flow fields—parallel and interdigitated flow fields. Fig. 3 shows the iteration history of topology-optimized flow field, in which the electrolyte domain is expressed as ρ ≥ 0.5. The computational domain is discretized using 2.4 × 105 uniform hexahedral elements for all variables in this study. In the topology-optimized design shown in Fig. 3, some of the flow channels are not connected in the topology-optimized flow field, which can enhance the mass transfer effect and further reduce the mass transfer loss. We believe that the current collector with the optimized design can be machined via 5

Journal of Energy Storage 26 (2019) 100990

C.-H. Chen, et al.

Fig. 5. Comparison results of different flow fields on the middle plane of the electrode: (a) Velocity; (b) Pressure; (c) Concentration of V 2 + ; (d) Concentration of V 3 + .

electrode is larger than that of the reference flow fields, whereas the velocity magnitude of the inner flow channel is smaller. From this result, it is estimated that the optimized flow field aims to gain large under-rib convection by actively driving the electrolyte inside the electrode.

Table 3 shows the overpotential and maximum current density of the reference and optimized flow fields. Although the performance of the optimized flow field is the highest, its discrepancy from the reference flow fields is only approximately 1%, which cannot be regarded as an obvious difference. This is because the analysis model of the 6

Journal of Energy Storage 26 (2019) 100990

C.-H. Chen, et al.

Fig. 6. Comparison results of different flow fields on the sectional planes: (a) Velocity; (b) Pressure; (c) Concentration of V 2 + ; (d) Concentration of V 3 + .

numerical example is of a small scale. It should be noted that there is high possibility of observing a distinct difference between the optimized and non-optimized flow fields if a large-scale analysis model considering pilot-scale VRFBs is used [60]. To realize this enhanced

performance, it is clear that we need to overcome massive computational burden by utilizing large-scale parallel computers, which is beyond the scope of this work.

7

Journal of Energy Storage 26 (2019) 100990

C.-H. Chen, et al.

Table 3 Overpotential and maximum current density of the flow fields in Fig. 4. Flow field Overpotential (V)

Maximum current density (A m

2)

Parallel

Interdigitated

Optimized

−0.0204 20.8

−0.0201 20.8

−0.0198 20.7

[8] B. Li, M. Gu, Z. Nie, Y. Shao, Q. Luo, X. Wei, X. Li, J. Xiao, C. Wang, V. Sprenkle, W. Wang, Bismuth nanoparticle decorating graphite felt as a high-performance electrode for an all-vanadium redox flow battery, Nano Lett. 13 (3) (2013) 1330–1335. [9] J. Houser, J. Clement, A. Pezeshki, M.M. Mench, Influence of architecture and material properties on vanadium redox flow battery performance, J. Power Sources 302 (2016) 369–377. [10] B. Li, M. Gu, Z. Nie, X. Wei, C. Wang, V. Sprenkle, W. Wang, Nanorod niobium oxide as powerful catalysts for an all vanadium redox flow battery, Nano Lett. 14 (1) (2013) 158–165. [11] H. Wang, R. Liu, Y. Li, X. Lü, Q. Wang, S. Zhao, K. Yuan, Z. Cui, X. Li, S. Xin, R. Zhang, M. Lei, Z. Lin, Durable and efficient hollow porous oxide spinel microspheres for oxygen reduction, Joule 2 (2) (2018) 337–348. [12] X. He, S. Luan, L. Wang, R. Wang, P. Du, Y. Xu, H. Yang, Y. Wang, K. Huang, M. Lei, Facile loading mesoporous Co3O4 on nitrogen doped carbon matrix as an enhanced oxygen electrode catalyst, Mater. Lett. 244 (2019) 78–82. [13] K. Huang, J. Liu, L. Wang, G. Chang, R. Wang, M. Lei, Y. Wang, Y. He, Mixed valence CoCuMnOx spinel nanoparticles by sacrificial template method with enhanced ORR performance, Appl. Surf. Sci. 487 (2019) 1145–1151. [14] X. Zhou, T. Zhao, L. An, Y. Zeng, L. Wei, Critical transport issues for improving the performance of aqueous redox flow batteries, J. Power Sources 339 (2017) 1–12. [15] D. Aaron, Q. Liu, Z. Tang, G. Grim, A. Papandrew, A. Turhan, T. Zawodzinski, M. Mench, Dramatic performance gains in vanadium redox flow batteries through modified cell architecture, J. Power Sources 206 (2012) 450–453. [16] Q. Xu, T. Zhao, P. Leung, Numerical investigations of flow field designs for vanadium redox flow batteries, Appl. Energy 105 (2013) 47–56. [17] M. Messaggi, P. Canzi, R. Mereu, A. Baricci, F. Inzoli, A. Casalegno, M. Zago, Analysis of flow field design on vanadium redox flow battery performance: development of 3D computational fluid dynamic model and experimental validation, Appl. Energy 228 (2018) 1057–1070. [18] S. Tsushima, F. Kondo, S. Sasaki, S. Hirai, Efficient utilization of the electrodes in a redox flow battery by modifying flow field and electrode morphology, Proceedings of 15th International Heat Transfer Conference, (2014), https://doi.org/10.1615/ IHTC15.ecs.009326. [19] R.M. Darling, M.L. Perry, The influence of electrode and channel configurations on flow battery performance, J. Electrochem. Soc. 161 (9) (2014) A1381–A1387. [20] Y. Li, X. Zhang, J. Bao, M. Skyllas-Kazacos, Studies on optimal charging conditions for vanadium redox flow batteries, J. Energy Storage 11 (2017) 191–199. [21] M. Yue, Q. Zheng, F. Xing, H. Zhang, X. Li, X. Ma, Flow field design and optimization of high power density vanadium flow batteries: a novel trapezoid flow battery, AlChE J. 64 (2) (2018) 782–795. [22] X. You, Q. Ye, P. Cheng, Scale-up of high power density redox flow batteries by introducing interdigitated flow fields, Int. Commun. Heat Mass Transf. 75 (2016) 7–12. [23] M.P. Bendsøe, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Comput. Methods Appl. Mech. Eng. 71 (2) (1988) 197–224. [24] M.P. Bendsøe, O. Sigmund, Topology Optimization: Theory, Methods, and Applications, Springer, 2003. [25] K. Yaji, M. Otomori, T. Yamada, K. Izui, S. Nishiwaki, O. Pironneau, Shape and topology optimization based on the convected level set method, Struct. Multidiscip. Optim. 54 (3) (2016) 659–672. [26] Q. Li, G.P. Steven, O.M. Querin, Y. Xie, Shape and topology design for heat conduction by evolutionary structural optimization, Int. J. Heat Mass Transf. 42 (17) (1999) 3361–3371. [27] A. Iga, S. Nishiwaki, K. Izui, M. Yoshimura, Topology optimization for thermal conductors considering design-dependent effects, including heat conduction and convection, Int. J. Heat Mass Transf. 52 (11–12) (2009) 2721–2732. [28] T. Nomura, K. Sato, K. Taguchi, T. Kashiwa, S. Nishiwaki, Structural topology optimization for the design of broadband dielectric resonator antennas using the finite difference time domain technique, Int. J. Numer. Methods Eng. 71 (11) (2007) 1261–1296. [29] S. Yamasaki, T. Nomura, A. Kawamoto, K. Sato, S. Nishiwaki, A level set-based topology optimization method targeting metallic waveguide design problems, Int. J. Numer. Methods Eng. 87 (9) (2011) 844–868. [30] K. Nomura, S. Yamasaki, K. Yaji, H. Bo, A. Takahashi, T. Kojima, K. Fujita, Topology optimization of conductors in electrical circuit, Struct. Multidiscip. Optim. 59 (6) (2019) 2205–2225. [31] O. Sigmund, Design of multiphysics actuators using topology optimization—Part I: one-material structures, Comput. Methods Appl. Mech. Eng. 190 (49–50) (2001) 6577–6604. [32] H. Iwai, A. Kuroyanagi, M. Saito, A. Konno, H. Yoshida, T. Yamada, S. Nishiwaki, Power generation enhancement of solid oxide fuel cell by cathode–electrolyte interface modification in mesoscale assisted by level set-based optimization calculation, J. Power Sources 196 (7) (2011) 3485–3495. [33] R. Behrou, A. Pizzolato, A. Forner-Cuenca, Topology optimization as a powerful tool to design advanced PEMFCs flow fields, Int. J. Heat Mass Transf. 135 (2019) 72–92. [34] J. Onishi, Y. Kametani, Y. Hasegawa, N. Shikazono, Topology optimization of electrolyte-electrode interfaces of solid oxide fuel cells based on the adjoint method, J. Electrochem. Soc. 166 (13) (2019) F876–F888. [35] T. Borrvall, J. Petersson, Topology optimization of fluids in Stokes flow, Int. J. Numer. Methods Fluids 41 (1) (2003) 77–107. [36] A. Gersborg-Hansen, O. Sigmund, R.B. Haber, Topology optimization of channel flow problems, Struct. Multidiscip. Optim. 30 (3) (2005) 181–192. [37] L.H. Olesen, F. Okkels, H. Bruus, A high-level programming-language

5. Conclusions In this paper, we proposed a three-dimensional topology optimization method for the flow field in VRFBs. We demonstrated that the optimization problem can be formulated as a maximization problem of the electrode surface concentration of oxidized reactants during the charging process. Based on the proposed formulations, we investigated the efficacy of the proposed approach through a numerical example. As a result, we confirmed that topology optimization could be applied to the optimization problem of a VRFB by considering the electrochemical reaction kinetics. In addition, we compared the optimized flow field with the reference flow fields—parallel and interdigitated flow fields—in terms of the overpotential and maximum current density. Although the result indicated that the optimized flow field has slightly better performance, the investigated flow fields showed almost the same performance owing to the small-scale analysis models. The main contribution of this paper is it reveals that topology optimization can realize the automatic generation of a three-dimensional optimized flow field of a VRFB. As a result, we revealed that the interdigitated-type is an optimal flow field of the VRFB, at least under the operating conditions investigated in this study. In our future work, we plan to deal with a large-scale analysis model considering pilot-scale VRFBs to clarify the distinct difference between the optimized and non-optimized flow fields. It also remains a challenge for future research to perform necessary experiments to verify the actual performance of the topology-optimized flow field. Further insight into this aspect will be presented in our future work. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments This work is partially supported by a research grant from The Mazda Foundation. Supplementary material Supplementary material associated with this article can be found, in the online version, at 10.1016/j.est.2019.100990. References [1] S. Chu, A. Majumdar, Opportunities and challenges for a sustainable energy future, Nature 488 (7411) (2012) 294. [2] Z. Yang, J. Zhang, M.C. Kintner-Meyer, X. Lu, D. Choi, J.P. Lemmon, J. Liu, Electrochemical energy storage for green grid, Chem. Rev. 111 (5) (2011) 3577–3613. [3] E. Sum, M. Skyllas-Kazacos, A study of the V(II)/V(III) redox couple for redox flow cell applications, J. Power Sources 15 (2–3) (1985) 179–190. [4] E. Sum, M. Rychcik, M. Skyllas-Kazacos, Investigation of the V(V)/V(IV) system for use in the positive half-cell of a redox battery, J. Power Sources 16 (1985) 85–95. [5] P. Alotto, M. Guarnieri, F. Moro, Redox flow batteries for the storage of renewable energy: a review, Renew. Sustain. Energy Rev. 29 (2014) 325–335. [6] B. Dunn, H. Kamath, J.-M. Tarascon, Electrical energy storage for the grid: a battery of choices, Science 334 (6058) (2011) 928–935. [7] W. Wang, X. Wang, Investigation of IR-modified carbon felt as the positive electrode of an all-vanadium redox flow battery, Electrochim. Acta 52 (24) (2007) 6755–6762.

8

Journal of Energy Storage 26 (2019) 100990

C.-H. Chen, et al.

[38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48]

implementation of topology optimization applied to steady-state Navier–Stokes flow, Int. J. Numer. Methods Eng. 65 (7) (2006) 975–1001. S. Kubo, K. Yaji, T. Yamada, K. Izui, S. Nishiwaki, A level set-based topology optimization method for optimal manifold designs with flow uniformity in plate-type microchannel reactors, Struct. Multidiscip. Optim. 55 (4) (2017) 1311–1327. G.H. Yoon, Topology optimization for turbulent flow with Spalart–Allmaras model, Comput. Methods Appl. Mech. Eng. 303 (2016) 288–311. C.B. Dilgen, S.B. Dilgen, D.R. Fuhrman, O. Sigmund, B.S. Lazarov, Topology optimization of turbulent flows, Comput. Methods Appl. Mech. Eng. 331 (2018) 363–393. G.H. Yoon, Topology optimization for stationary fluid–structure interaction problems using a new monolithic formulation, Int. J. Numer. Methods Eng. 82 (5) (2010) 591–616. N. Jenkins, K. Maute, Level set topology optimization of stationary fluid-structure interaction problems, Struct. Multidiscip. Optim. 52 (1) (2015) 179–195. T. Matsumori, T. Kondoh, A. Kawamoto, T. Nomura, Topology optimization for fluid–thermal interaction problems under constant input power, Struct. Multidiscip. Optim. 47 (4) (2013) 571–581. K. Yaji, T. Yamada, S. Kubo, K. Izui, S. Nishiwaki, A topology optimization method for a coupled thermal–fluid problem using level set boundary expressions, Int. J. Heat Mass Transf. 81 (2015) 878–888. K. Yaji, M. Ogino, C. Chen, K. Fujita, Large-scale topology optimization incorporating local-in-time adjoint-based method for unsteady thermal-fluid problem, Struct. Multidiscip. Optim. 58 (2) (2018) 817–822. H. Kobayashi, K. Yaji, S. Yamasaki, K. Fujita, Freeform winglet design of fin-andtube heat exchangers guided by topology optimization, Appl. Therm. Eng. (2019) 114020. J. Alexandersen, N. Aage, C.S. Andreasen, O. Sigmund, Topology optimisation for natural convection problems, Int. J. Numer. Methods Fluids 76 (10) (2014) 699–721. P. Coffin, K. Maute, A level-set method for steady-state and transient natural convection problems, Struct. Multidiscip. Optim. 53 (5) (2016) 1047–1067.

[49] J. Alexandersen, O. Sigmund, N. Aage, Large scale three-dimensional topology optimisation of heat sinks cooled by natural convection, Int. J. Heat Mass Transf. 100 (2016) 876–891. [50] E. Kontoleontos, E. Papoutsis-Kiachagias, A. Zymaris, D. Papadimitriou, K. Giannakoglou, Adjoint-based constrained topology optimization for viscous flows, including heat transfer, Eng. Optim. 45 (8) (2013) 941–961. [51] S.B. Dilgen, C.B. Dilgen, D.R. Fuhrman, O. Sigmund, B.S. Lazarov, Density based topology optimization of turbulent flow heat transfer systems, Struct. Multidiscip. Optim. 57 (5) (2018) 1905–1918. [52] K. Yaji, S. Yamasaki, S. Tsushima, T. Suzuki, K. Fujita, Topology optimization for the design of flow fields in a redox flow battery, Struct. Multidiscip. Optim. 57 (2) (2018) 535–546. [53] A. Shah, M. Watt-Smith, F. Walsh, A dynamic performance model for redox-flow batteries involving soluble species, Electrochim. Acta 53 (27) (2008) 8087–8100. [54] D. You, H. Zhang, J. Chen, A simple model for the vanadium redox battery, Electrochim. Acta 54 (27) (2009) 6827–6836. [55] X. Ma, H. Zhang, F. Xing, A three-dimensional model for negative half cell of the vanadium redox flow battery, Electrochim. Acta 58 (2011) 238–246. [56] M.M. Tomadakis, T.J. Robertson, Viscous permeability of random fiber structures: comparison of electrical and diffusional estimates with experimental and analytical results, J. Compos. Mater. 39 (2) (2005) 163–188. [57] S. Yamasaki, K. Yaji, K. Fujita, Knowledge discovery in databases for determining formulation in topology optimization, Struct. Multidiscip. Optim. 59 (2) (2019) 595–611. [58] H. Zhou, H. Zhang, P. Zhao, B. Yi, A comparative study of carbon felt and activated carbon based electrodes for sodium polysulfide/bromine redox flow battery, Electrochim. Acta 51 (28) (2006) 6304–6312. [59] M. Pourbaix, Atlas of Electrochemical Equilibria in Aqueous Solutions, second ed., NACE International, Houston, 1974. [60] L. Arenas, C.P. De León, F. Walsh, Engineering aspects of the design, construction and performance of modular redox flow batteries for energy storage, J. Energy Storage 11 (2017) 119–153.

9