# Modeling and Behavior Analysis of a Membraneless Fuel Cell

Authors: Ciro-Filemon Flores-Rivera
Journal: ISRN Applied Mathematics
 (2012)

---

## Abstract

Membraneless fuel cells are examples of microelectromechanical systems (MEMSs) that can be considered as alternate energy sources. Applications include microfluidic-based devices like miniaturized laboratories, sensors, or actuators to be used in medicine or agronomy. This paper presents a mathematical model for this type of cells based on the governing physical laws. It includes fluid dynamics, electric charge distribution and electrostatics modeled by the Navier-Stokes, Nernst-Planck, and Poisson equations, respectively. A robust numerical algorithm is proposed to solve the model. Two cases are discussed: allowing electrochemical reactions on one of the electrodes and the simpler situation of null exchange current density. An initial characterization for the behavior of membraneless fuel cells is achieved concerning to  prevalence of velocity and electric field, use of non-Newtonian fluids, relationship to initial conditions for some variables, general profile for conductivity and electric density, and linear dependence on current density under specific conditions.


---

## Body

## 1. Introduction

In the past two decades, there has been an intense effort in the research of microelectromechanical systems (MEMS) [1], due to the broad spectrum of their applications in engineering or life sciences, including miniaturized laboratories (labs-on-a-chip). Microfuel cells are good examples of MEMS in the search of alternate energy sources and commonly are essential components of industrial labs-on-a-chip. A fuel cell is a device that converts chemical energy from a fuel and an oxidant into electric energy, by means of oxidoreduction reactions. Fundamental physical and electrochemical processes involved are far from being well understood. Investigation on the behavior of these devices is at an early stage; nevertheless, actual fuel cells have been built showing their feasibility, even though their efficiency is limited [2]. This motivates research like present work toward a better understanding of fuel cell  operation.The basic design of a fuel cell considers three main components which are sandwiched together: anode, electrolyte, and cathode [3]. An external load is connected between both electrodes. Usually the electrolyte is made up of two fluids separated by a polymer membrane, which allows ions transport (H+) from the anode to the cathode [4]. Assembly of the membrane increases the cell cost, and its deterioration implies the necessity of frequent replenishments. An attractive alternate design is offered by membraneless fuel cells with obvious cost reduction. In these devices separation between reactant fluids is achieved only by stable laminar flows. Therefore, it is important to determine conditions to maintain stability on the flow for the proper functioning of the fuel cell [5].Figure1 shows the simplified model of a membraneless fuel cell assumed in present work, according to the basic geometry proposed by [2, 6]. Reactants are fed from left input and flow in a laminar way through a channel with an anode on the fuel side and a cathode on the oxidant side. Reaction products flow out from the right output while electricity is produced electrochemically in the electrolyte. Many combinations of fuels and oxidants are possible. For instance, a hydrogen fuel cell uses hydrogen as fuel and oxygen as oxidant. Other fuels include hydrocarbons and alcohol, whereas other oxidants could be chlorine or chlorine dioxide [7]. Electricity is used by an external load (or resistance) connected between the electrodes, usually made of platinum. In this way the fuel cell becomes a real source of energy.Figure 1
Simplified model of a membraneless fuel cell. Fuel and Oxidant flow from left to right and are separated by laminar flow properties. Oxidoreduction reactions generate ions (H+) which are allowed to traverse from anode to cathode producing an electric current that is used by an external load.Present paper is a continuation of [8] and is organized as follows. Section 2 establishes the physical and mathematical model in terms of a system of coupled, nonlinear, time-dependent, partial differential equations (PDEs), taking into account the physical laws involved, which come from fluid dynamics, electric charge distribution, and electrostatics. The corresponding boundary conditions and the numerical strategy followed to solve such system of equations are given as well. In Section 3, results of numerical experiments are presented for two cases: first when a null exchange current density is assumed at both electrodes  and second when an Arrhenius-like electrical reaction is allowed at the anode. This represents the main contribution of this work: a complex mathematical model to describe the behavior of a membraneless fuel cell, under different circumstances that will be addressed. Furthermore, several initial conditions are considered for both cases, and some technical considerations about the finite-element method (FEM) used for the numerical solution are presented. Section 4 discusses main results and presents the findings achieved by this research. Finally, some concluding remarks are given in Section 5.


## 2. Physical and Mathematical Model

To develop a suitable mathematical model for a membraneless fuel cell as that shown in Figure1, some simplifications will be assumed. Since electrochemical reactions among charge carriers in the reactants and electrodes may be very complex, depending on the specific fluids considered, the model incorporates the fundamental phenomena. It will be considered that the anionic and cationic charge carriers in fuel and oxidant are the same [8], meaning that both reactants carry the same charge excess (z±), then(2.1)1=z+(fuel)=z+(oxidant)=-z-(fuel)=-z-(oxidant).Furthermore, diffusivity (D) and mobility (w) coefficients for fuel and oxidant will be considered as the same. Both fluids are supposed to be incompressible and Newtonian (even though this last assumption will be relaxed later), having the same dynamic viscosity η and the same mass density ρ. Table 1 summarizes the specific values for the parameters appearing in the model.Table 1
Parameters used in the simulations of the membraneless fuel cell: first without reaction at electrodes and second considering an Arrhenius-like electric reaction on the anode. Most of them have been chosen according to experimental designs reported in [17, 18].
ParameterValueDescriptionρ1000 kg/m3Density of the fluidη0.001 Pa·sDynamic viscosity of the fluidU00.1 mm/sMean inflow velocityεr78.3Relative electric permittivity of the fluidζ−0.07 VZeta-potential on the wall fluid boundaryw8.2 × 10−13 mol·s/kgMobility coefficientD2 × 10−9 m2/sDiffusion coefficientc00.5 × 10−10 mol/m3Initial concentrationϕ02 × 103 VInitial electric potentialT293.15°KAbsolute temperatureκ4.459 × 109Constant for electric displacementi00.05 A/m2Exchange current density at anodeαa,αc0.5Anodic, cathodic charge transfer coefficientsThis section starts defining the mathematical domain for the membraneless fuel cell. Then, it is established the set of governing equations with appropriate boundary conditions. Finally, the numerical strategy followed to solve the mathematical model is presented.


### 2.1. System of Governing Equations

DomainΩ representing the simplified model of a membraneless fuel cell as shown in Figure 1, consists of a rectangle: 10 mm large and 2 mm height, as it is shown in Figure 2. This rectangle is placed in a Cartesian coordinate system xy, the lower left corner being the point (0,-0.001). There are four boundaries identified as follows: vertical left and right boundaries Γ1 and Γ2 are the inlet and the outlet, respectively, whereas horizontal boundaries Γ3 and Γ4 are the electrodes. The lower electrode Γ3 is the cathode and the upper one Γ4 is the anode. Two different fluids are entering the right by the inlet: the upper half corresponds to the fuel, while the lower half to the oxidant. Redox reactions between these two reactants start producing ions, and an electric current begins to flow in the vertical direction because of the electrodes. Three physical laws are involved and serve as basis to depict the phenomena and to build up the proper mathematical model. They are described next.Figure 2
DomainΩ and boundaries. Γ1 is the inlet, Γ2 is the outlet, Γ3 and Γ4 are the electrodes. Fluids flow from left to right. The upper fluid in the anode side is the fuel, while the cathode side one is the oxidant. Electric current flows in the vertical direction.


#### 2.1.1. Fluid Dynamics (Navier-Stokes’ Equations)

Starting with momentum and mass balance in terms of stresses—in the micrometer scale of the membraneless fuel cell, dynamics of fluids is modeled by the generalized Navier-Stokes equations (2.2)-(2.3), in terms of transport properties and velocity gradients [9]. The term within brackets is the total stress tensor and will be used later in the establishment of proper boundary conditions. Here η is the dynamic viscosity, ρ is the fluid density, u the velocity field, p the pressure, and f the external force per unit volume. Equation (2.2) allows the possibility of considering variable viscosity of the fluid, that is, non-Newtonian fluids. This analysis will be addressed in Section 4. Note that velocity u is a time-dependent variable:(2.2)ρ∂u∂t+ρ(u⋅∇)u=∇⋅[-pI+η(∇u+(∇u)t)]+f,(2.3)∇⋅u=0.


#### 2.1.2. Charge Distribution (Nernst-Planck’s Equations)

According to [10], the distribution of electric charge in a channel filled with two moving dilute electrolyte solutions is modeled by the Nernst-Planck equations (2.4). D1,2 are the diffusion coefficients, w+,- are the mobility coefficients, z+,- are the valence numbers as stated in   (2.1), F is the Faraday’s constant, ϕ is the electric potential, and c1,2 are the charge molar concentrations of reactants (positive and negative, resp.). Note that these concentrations are time-dependent variables
(2.4)∂c1∂t+u⋅∇c1=∇⋅(D1∇c1+w+z+Fc1∇ϕ),∂c2∂t+u⋅∇c2=∇⋅(D2∇c2+w-z-Fc2∇ϕ).


#### 2.1.3. Electrostatics’s (Poisson Equation)

Based on electromagnetic theory [11] and assuming there is no a magnetic field and the electric field E is quasi-static, Maxwell’s equations are ∇×E=0; hence, E=-∇ϕ. Then, the Gauss law ∇·E=ρq/ε0εr, where ρq is charge density and εr is the relative electric permittivity of the electrolyte, implies the electric potential ϕ must satisfy the Poisson equation (2.5):(2.5)-ε0εr=∇2ϕ=ρq.On the other hand, there are two important physical quantities for a fuel cell: electric conductivityσ and charge density ρq, defined by (2.6) and (2.7), respectively, in terms of the charge molar concentrations. These quantities will be used to report results rather than concentrations because of their important physical interpretation:
(2.6)σ=F2w(c1+c2),(2.7)ρq=F(c1-c2).Furthermore, the external forcef in (2.2) is the electromagnetic force given by(2.8)f=ρqE=-F(c1-c2)∇ϕ.Therefore, recallingD1=D2=D and w+=w-=w and taking into account previous discussion, the complete set of coupled, nonlinear, time-dependent, partial differential equations (PDEs) governing the model for the membraneless fuel cell is given by (2.9)–(2.13) and must be satisfied in domain Ω(2.9)ρ∂u∂t+ρ(u⋅∇)u=∇⋅[-pI+η(∇u+(∇u)t)]-F(c1-c2)∇ϕ,(2.10)∇⋅u=0,(2.11)∂c1∂t+u⋅∇c1=∇⋅(D∇c1+wFc1∇ϕ),(2.12)∂c2∂t+u⋅∇c2=∇⋅(D∇c2-wFc2∇ϕ),(2.13)-ε0εr∇2ϕ=F(c1-c2).Although the system is complex, boundary conditions for the Navier-Stokes, Nernst-Planck, and Poisson equations are decoupled, as will be shown below. This allows the application of a convenient numerical strategy to solve it.


### 2.2. Boundary Conditions

The establishment of boundary conditions that properly describe the functioning of the membraneless fuel cell is a very important part in the modeling process. Based on the physical laws invoked above, they are set as follows. As usually,n represents the unitary outward normal vector.


#### 2.2.1. Fluid Dynamics

Laminar inflow is assumed for the fluids, and it is expected a close to zero pressure by the end of the channel. Since fluids flow in the horizontal direction and the channel is very thin, a nonslip condition is imposed on the electrodes. These conditions are in accordance with fluid dynamics theory at the scale of micrometers [12]. Details are given below.Inlet (Laminar Inflow)
Boundary conditions on the inlet are given by   (2.14). They simulate the effect of a fictious channel outside of the true geometry as shown in Figure 3. The subscript t emphasizes that such equations are solved for tangential components of the boundary. Lentr is large enough (≈5 × 10−3m) so that the flow can develop its laminar form, while pentr is solved so that the average velocity of the fluid through the boundary Γ1 meets U0. This implies that in the microscale, with the low Reynolds number used (≈0.1), the flow has a fully developed laminar profile on the inlet:
(2.14)Lentr∇t⋅[pI-η(∇tu+(∇tu)t)]=-npentr,∇t⋅u=0.Figure 3
Sketch of the physical situation simulated by a laminar inflow. Solid line is the actual computational domain while dashed line is a fictitious domain. The fluid, flowing from left to right, achieves a fully developed laminar profile on the inlet.Outlet (Normal Stress)
The boundary condition to be satisfied on the outlet is given by (2.15). This means that total stress in the tangential direction is zero. Under these circumstances pressure in the outlet is approximately zero:
(2.15)[-pI+η(∇u+(∇u)t)]n=0.Electrodes (No Slip)
Due to the microscale of the channel, fluid’s velocity is set to zero at these two walls. Equation (2.16) establishes this condition:
(2.16)u=0.


#### 2.2.2. Charge Distribution

According to [13], it is assumed that there is no normal flux of anions at the electrodes,   at the inlet, nor at the outlet. It corresponds to the basic model which will be referred in this work as Case A: null exchange current density. Later on a relaxation on the condition for the anode will be allowed, and cations will be assumed to undergo reactions at a rate of an Arrhenius-like form. This will be called in this paper Case B: nonzero exchange current density. Section 3 addresses details for this modified case. Then the boundary conditions for Case A are stated as follows.Inlet, Outlet, and Electrodes (Convective Flux)
This condition states that the transport of ions due to diffusive or electric migration effects, perpendicular to the four boundaries, is negligible. Thus, (2.17)-(2.18) model this condition:
(2.17)n⋅(-D∇c1-wFc1∇ϕ)=0,(2.18)n⋅(-D∇c2+wFc2∇ϕ)=0.


#### 2.2.3. Electrostatics

Inlet and outlet are considered to have zero electric charge because of the absence of electrodes. A very different situation happens on the horizontal frontiers. According to [14, 15], the existence of electrodes, the applied electric field on the vertical direction, and the flow of the electrolyte in the channel  lead—in the fuel cell microscale—to the well-known [12] phenomenon of the electrical double layer (EDL). It concerns two parallel layers of charge surrounding the electrodes. The first layer, called the Stern or compact layer, comprises ions adsorbed directly onto the electrode (wall) due to the host of chemical interactions. Its thickness is the order of one molecule diameter. The second layer, called the diffuse layer, is composed of ions attracted to the surface charge via the Coulomb force, electrically screening the first layer. The second layer is loosely associated to the wall, because it is made of free ions which move in the fluid (electrolyte) under the influence of electric attraction rather than being firmly anchored [12]. Its thickness is the order of 10−8 to 10−7m. The remaining part of the cell, which corresponds to almost all the channel height, is called the bulk. Figure 4 shows these components schematically without scales. The Gouy-Chapman-Stern theory [13, 15] states that ∇ϕ is a constant in the Stern layer and approaches zero in the diffuse layer. Then a linear drop voltage, that is, an electric displacement, can be assumed as appropriate boundary condition on electrodes. Details are presented next.Figure 4
Electrical double-layer (EDL) effect close to an electrode. The Stern layer is compact and comprises ions adsorbed directly onto the electrode. Its thickness is the order of one molecule diameter. The diffuse layer is composed by ions attracted to the surface via Coulomb-forces. Its thickness is the order of 10−8 to 10−7m. The bulk is the remaining part of the channel. Picture is drawn without scales.Inlet and Outlet (Zero Charge)
This condition assumes that the normal component of the electric field is zero and is given by(2.19)∇ϕ⋅n=0.Upper Electrode (Electric Displacement)
The normal component of the electric field on the anode is specified as stated in(2.20)∇ϕ⋅n=κ(ϕ0-ϕ),
where
(2.21)κ=1δγ,δ=ΔCrefCref,γ=λDh,Cref is the reference molar concentration, λD is the Debye length (thickness of the diffuse layer), ϕ0 the initial electric potential, and h is a half of the channel height. According to [16–18],  appropriate values for these parameters taken from experimental designs, are
(2.22)δ=1.121212×10-5,γ=2×10-5⟹κ=4.459×109.Lower Electrode (Electric Displacement)
In the case of the cathode, the normal component of the electric field on this boundary is specified as shown in(2.23)∇ϕ⋅n=-κ(ϕ0+ϕ).Therefore, the complete mathematical model for the membraneless fuel cell is given by the set of (2.9)–(2.13) to be satisfied in domain Ω, together with boundary conditions stated by (2.14)–(2.23) to be fulfilled on the corresponding boundaries.As extension of [8], present work is not limited to the case of electroneutrality (ρq=0) or electrochemical cells (u=0)  but, following proposal in [16], considers and gives solution to the complete complex system just described which models a more realistic situation. Furthermore, nonzero electrical reactions on the anode, that is, variations on   (2.11) and boundary condition (2.17), will be addressed as well.


### 2.3. Numerical Strategy for Solution

Since system (2.9)–(2.13) is highly nonlinear and time dependent, an effective numerical strategy must be implemented to integrate the included equations. Actually the system contains the following five unknowns: the bidimensional velocity field u, the pressure p, the charge concentrations c1 and c2, and the electric potential ϕ. The proposed numerical strategy is a robust combination of time integration by operator-splitting techniques [19, 20] and a finite element method for space integration [21].Then, ifΔt is the time discretization step and using the following notation(2.24)un≡u(x,nΔt),pn≡p(x,nΔt),c1n≡c1(x,nΔt),c2n≡c2(x,nΔt),ϕn≡ϕ(x,nΔt),
the numerical strategy can be summarized as follows.Given the initial valuesu0,p0,c10,c20, and ϕ0, find the solution at time n+1 solving (2.25)–(2.29), for n=0,1,2,…,(2.25)ρun+1-unΔt+ρ(un+1⋅∇)un=∇⋅[-pn+1I+η(∇un+(∇un)t)]-F(c1n-c2n)∇ϕn,(2.26)∇⋅un+1=0,(2.27)c1n+1-c1nΔt+un+1⋅∇c1n+1=∇⋅(D∇c1n+wFc1n∇ϕn),(2.28)c2n+1-c2nΔt+un+1⋅∇c2n+1=∇⋅(D∇c2n-wFc2n∇ϕn),(2.29)-ε0εr∇2ϕn+1=F(c1n+1-c2n+1).This algorithm quickly converges to a steady-state solution. According to [22], in the microscale of the fuel cell, it is enough to simulate at least 0.013 seconds in order to capture the dynamics of the phenomenon. In this work all of the simulations were run up to 1 second to assure a steady-state solution. Results are reported in the following section. Details on the FEM used for the numerical solutions are also included.


## 2.1. System of Governing Equations

DomainΩ representing the simplified model of a membraneless fuel cell as shown in Figure 1, consists of a rectangle: 10 mm large and 2 mm height, as it is shown in Figure 2. This rectangle is placed in a Cartesian coordinate system xy, the lower left corner being the point (0,-0.001). There are four boundaries identified as follows: vertical left and right boundaries Γ1 and Γ2 are the inlet and the outlet, respectively, whereas horizontal boundaries Γ3 and Γ4 are the electrodes. The lower electrode Γ3 is the cathode and the upper one Γ4 is the anode. Two different fluids are entering the right by the inlet: the upper half corresponds to the fuel, while the lower half to the oxidant. Redox reactions between these two reactants start producing ions, and an electric current begins to flow in the vertical direction because of the electrodes. Three physical laws are involved and serve as basis to depict the phenomena and to build up the proper mathematical model. They are described next.Figure 2
DomainΩ and boundaries. Γ1 is the inlet, Γ2 is the outlet, Γ3 and Γ4 are the electrodes. Fluids flow from left to right. The upper fluid in the anode side is the fuel, while the cathode side one is the oxidant. Electric current flows in the vertical direction.


### 2.1.1. Fluid Dynamics (Navier-Stokes’ Equations)

Starting with momentum and mass balance in terms of stresses—in the micrometer scale of the membraneless fuel cell, dynamics of fluids is modeled by the generalized Navier-Stokes equations (2.2)-(2.3), in terms of transport properties and velocity gradients [9]. The term within brackets is the total stress tensor and will be used later in the establishment of proper boundary conditions. Here η is the dynamic viscosity, ρ is the fluid density, u the velocity field, p the pressure, and f the external force per unit volume. Equation (2.2) allows the possibility of considering variable viscosity of the fluid, that is, non-Newtonian fluids. This analysis will be addressed in Section 4. Note that velocity u is a time-dependent variable:(2.2)ρ∂u∂t+ρ(u⋅∇)u=∇⋅[-pI+η(∇u+(∇u)t)]+f,(2.3)∇⋅u=0.


### 2.1.2. Charge Distribution (Nernst-Planck’s Equations)

According to [10], the distribution of electric charge in a channel filled with two moving dilute electrolyte solutions is modeled by the Nernst-Planck equations (2.4). D1,2 are the diffusion coefficients, w+,- are the mobility coefficients, z+,- are the valence numbers as stated in   (2.1), F is the Faraday’s constant, ϕ is the electric potential, and c1,2 are the charge molar concentrations of reactants (positive and negative, resp.). Note that these concentrations are time-dependent variables
(2.4)∂c1∂t+u⋅∇c1=∇⋅(D1∇c1+w+z+Fc1∇ϕ),∂c2∂t+u⋅∇c2=∇⋅(D2∇c2+w-z-Fc2∇ϕ).


### 2.1.3. Electrostatics’s (Poisson Equation)

Based on electromagnetic theory [11] and assuming there is no a magnetic field and the electric field E is quasi-static, Maxwell’s equations are ∇×E=0; hence, E=-∇ϕ. Then, the Gauss law ∇·E=ρq/ε0εr, where ρq is charge density and εr is the relative electric permittivity of the electrolyte, implies the electric potential ϕ must satisfy the Poisson equation (2.5):(2.5)-ε0εr=∇2ϕ=ρq.On the other hand, there are two important physical quantities for a fuel cell: electric conductivityσ and charge density ρq, defined by (2.6) and (2.7), respectively, in terms of the charge molar concentrations. These quantities will be used to report results rather than concentrations because of their important physical interpretation:
(2.6)σ=F2w(c1+c2),(2.7)ρq=F(c1-c2).Furthermore, the external forcef in (2.2) is the electromagnetic force given by(2.8)f=ρqE=-F(c1-c2)∇ϕ.Therefore, recallingD1=D2=D and w+=w-=w and taking into account previous discussion, the complete set of coupled, nonlinear, time-dependent, partial differential equations (PDEs) governing the model for the membraneless fuel cell is given by (2.9)–(2.13) and must be satisfied in domain Ω(2.9)ρ∂u∂t+ρ(u⋅∇)u=∇⋅[-pI+η(∇u+(∇u)t)]-F(c1-c2)∇ϕ,(2.10)∇⋅u=0,(2.11)∂c1∂t+u⋅∇c1=∇⋅(D∇c1+wFc1∇ϕ),(2.12)∂c2∂t+u⋅∇c2=∇⋅(D∇c2-wFc2∇ϕ),(2.13)-ε0εr∇2ϕ=F(c1-c2).Although the system is complex, boundary conditions for the Navier-Stokes, Nernst-Planck, and Poisson equations are decoupled, as will be shown below. This allows the application of a convenient numerical strategy to solve it.


## 2.1.1. Fluid Dynamics (Navier-Stokes’ Equations)

Starting with momentum and mass balance in terms of stresses—in the micrometer scale of the membraneless fuel cell, dynamics of fluids is modeled by the generalized Navier-Stokes equations (2.2)-(2.3), in terms of transport properties and velocity gradients [9]. The term within brackets is the total stress tensor and will be used later in the establishment of proper boundary conditions. Here η is the dynamic viscosity, ρ is the fluid density, u the velocity field, p the pressure, and f the external force per unit volume. Equation (2.2) allows the possibility of considering variable viscosity of the fluid, that is, non-Newtonian fluids. This analysis will be addressed in Section 4. Note that velocity u is a time-dependent variable:(2.2)ρ∂u∂t+ρ(u⋅∇)u=∇⋅[-pI+η(∇u+(∇u)t)]+f,(2.3)∇⋅u=0.


## 2.1.2. Charge Distribution (Nernst-Planck’s Equations)

According to [10], the distribution of electric charge in a channel filled with two moving dilute electrolyte solutions is modeled by the Nernst-Planck equations (2.4). D1,2 are the diffusion coefficients, w+,- are the mobility coefficients, z+,- are the valence numbers as stated in   (2.1), F is the Faraday’s constant, ϕ is the electric potential, and c1,2 are the charge molar concentrations of reactants (positive and negative, resp.). Note that these concentrations are time-dependent variables
(2.4)∂c1∂t+u⋅∇c1=∇⋅(D1∇c1+w+z+Fc1∇ϕ),∂c2∂t+u⋅∇c2=∇⋅(D2∇c2+w-z-Fc2∇ϕ).


## 2.1.3. Electrostatics’s (Poisson Equation)

Based on electromagnetic theory [11] and assuming there is no a magnetic field and the electric field E is quasi-static, Maxwell’s equations are ∇×E=0; hence, E=-∇ϕ. Then, the Gauss law ∇·E=ρq/ε0εr, where ρq is charge density and εr is the relative electric permittivity of the electrolyte, implies the electric potential ϕ must satisfy the Poisson equation (2.5):(2.5)-ε0εr=∇2ϕ=ρq.On the other hand, there are two important physical quantities for a fuel cell: electric conductivityσ and charge density ρq, defined by (2.6) and (2.7), respectively, in terms of the charge molar concentrations. These quantities will be used to report results rather than concentrations because of their important physical interpretation:
(2.6)σ=F2w(c1+c2),(2.7)ρq=F(c1-c2).Furthermore, the external forcef in (2.2) is the electromagnetic force given by(2.8)f=ρqE=-F(c1-c2)∇ϕ.Therefore, recallingD1=D2=D and w+=w-=w and taking into account previous discussion, the complete set of coupled, nonlinear, time-dependent, partial differential equations (PDEs) governing the model for the membraneless fuel cell is given by (2.9)–(2.13) and must be satisfied in domain Ω(2.9)ρ∂u∂t+ρ(u⋅∇)u=∇⋅[-pI+η(∇u+(∇u)t)]-F(c1-c2)∇ϕ,(2.10)∇⋅u=0,(2.11)∂c1∂t+u⋅∇c1=∇⋅(D∇c1+wFc1∇ϕ),(2.12)∂c2∂t+u⋅∇c2=∇⋅(D∇c2-wFc2∇ϕ),(2.13)-ε0εr∇2ϕ=F(c1-c2).Although the system is complex, boundary conditions for the Navier-Stokes, Nernst-Planck, and Poisson equations are decoupled, as will be shown below. This allows the application of a convenient numerical strategy to solve it.


## 2.2. Boundary Conditions

The establishment of boundary conditions that properly describe the functioning of the membraneless fuel cell is a very important part in the modeling process. Based on the physical laws invoked above, they are set as follows. As usually,n represents the unitary outward normal vector.


### 2.2.1. Fluid Dynamics

Laminar inflow is assumed for the fluids, and it is expected a close to zero pressure by the end of the channel. Since fluids flow in the horizontal direction and the channel is very thin, a nonslip condition is imposed on the electrodes. These conditions are in accordance with fluid dynamics theory at the scale of micrometers [12]. Details are given below.Inlet (Laminar Inflow)
Boundary conditions on the inlet are given by   (2.14). They simulate the effect of a fictious channel outside of the true geometry as shown in Figure 3. The subscript t emphasizes that such equations are solved for tangential components of the boundary. Lentr is large enough (≈5 × 10−3m) so that the flow can develop its laminar form, while pentr is solved so that the average velocity of the fluid through the boundary Γ1 meets U0. This implies that in the microscale, with the low Reynolds number used (≈0.1), the flow has a fully developed laminar profile on the inlet:
(2.14)Lentr∇t⋅[pI-η(∇tu+(∇tu)t)]=-npentr,∇t⋅u=0.Figure 3
Sketch of the physical situation simulated by a laminar inflow. Solid line is the actual computational domain while dashed line is a fictitious domain. The fluid, flowing from left to right, achieves a fully developed laminar profile on the inlet.Outlet (Normal Stress)
The boundary condition to be satisfied on the outlet is given by (2.15). This means that total stress in the tangential direction is zero. Under these circumstances pressure in the outlet is approximately zero:
(2.15)[-pI+η(∇u+(∇u)t)]n=0.Electrodes (No Slip)
Due to the microscale of the channel, fluid’s velocity is set to zero at these two walls. Equation (2.16) establishes this condition:
(2.16)u=0.


### 2.2.2. Charge Distribution

According to [13], it is assumed that there is no normal flux of anions at the electrodes,   at the inlet, nor at the outlet. It corresponds to the basic model which will be referred in this work as Case A: null exchange current density. Later on a relaxation on the condition for the anode will be allowed, and cations will be assumed to undergo reactions at a rate of an Arrhenius-like form. This will be called in this paper Case B: nonzero exchange current density. Section 3 addresses details for this modified case. Then the boundary conditions for Case A are stated as follows.Inlet, Outlet, and Electrodes (Convective Flux)
This condition states that the transport of ions due to diffusive or electric migration effects, perpendicular to the four boundaries, is negligible. Thus, (2.17)-(2.18) model this condition:
(2.17)n⋅(-D∇c1-wFc1∇ϕ)=0,(2.18)n⋅(-D∇c2+wFc2∇ϕ)=0.


### 2.2.3. Electrostatics

Inlet and outlet are considered to have zero electric charge because of the absence of electrodes. A very different situation happens on the horizontal frontiers. According to [14, 15], the existence of electrodes, the applied electric field on the vertical direction, and the flow of the electrolyte in the channel  lead—in the fuel cell microscale—to the well-known [12] phenomenon of the electrical double layer (EDL). It concerns two parallel layers of charge surrounding the electrodes. The first layer, called the Stern or compact layer, comprises ions adsorbed directly onto the electrode (wall) due to the host of chemical interactions. Its thickness is the order of one molecule diameter. The second layer, called the diffuse layer, is composed of ions attracted to the surface charge via the Coulomb force, electrically screening the first layer. The second layer is loosely associated to the wall, because it is made of free ions which move in the fluid (electrolyte) under the influence of electric attraction rather than being firmly anchored [12]. Its thickness is the order of 10−8 to 10−7m. The remaining part of the cell, which corresponds to almost all the channel height, is called the bulk. Figure 4 shows these components schematically without scales. The Gouy-Chapman-Stern theory [13, 15] states that ∇ϕ is a constant in the Stern layer and approaches zero in the diffuse layer. Then a linear drop voltage, that is, an electric displacement, can be assumed as appropriate boundary condition on electrodes. Details are presented next.Figure 4
Electrical double-layer (EDL) effect close to an electrode. The Stern layer is compact and comprises ions adsorbed directly onto the electrode. Its thickness is the order of one molecule diameter. The diffuse layer is composed by ions attracted to the surface via Coulomb-forces. Its thickness is the order of 10−8 to 10−7m. The bulk is the remaining part of the channel. Picture is drawn without scales.Inlet and Outlet (Zero Charge)
This condition assumes that the normal component of the electric field is zero and is given by(2.19)∇ϕ⋅n=0.Upper Electrode (Electric Displacement)
The normal component of the electric field on the anode is specified as stated in(2.20)∇ϕ⋅n=κ(ϕ0-ϕ),
where
(2.21)κ=1δγ,δ=ΔCrefCref,γ=λDh,Cref is the reference molar concentration, λD is the Debye length (thickness of the diffuse layer), ϕ0 the initial electric potential, and h is a half of the channel height. According to [16–18],  appropriate values for these parameters taken from experimental designs, are
(2.22)δ=1.121212×10-5,γ=2×10-5⟹κ=4.459×109.Lower Electrode (Electric Displacement)
In the case of the cathode, the normal component of the electric field on this boundary is specified as shown in(2.23)∇ϕ⋅n=-κ(ϕ0+ϕ).Therefore, the complete mathematical model for the membraneless fuel cell is given by the set of (2.9)–(2.13) to be satisfied in domain Ω, together with boundary conditions stated by (2.14)–(2.23) to be fulfilled on the corresponding boundaries.As extension of [8], present work is not limited to the case of electroneutrality (ρq=0) or electrochemical cells (u=0)  but, following proposal in [16], considers and gives solution to the complete complex system just described which models a more realistic situation. Furthermore, nonzero electrical reactions on the anode, that is, variations on   (2.11) and boundary condition (2.17), will be addressed as well.


## 2.2.1. Fluid Dynamics

Laminar inflow is assumed for the fluids, and it is expected a close to zero pressure by the end of the channel. Since fluids flow in the horizontal direction and the channel is very thin, a nonslip condition is imposed on the electrodes. These conditions are in accordance with fluid dynamics theory at the scale of micrometers [12]. Details are given below.Inlet (Laminar Inflow)
Boundary conditions on the inlet are given by   (2.14). They simulate the effect of a fictious channel outside of the true geometry as shown in Figure 3. The subscript t emphasizes that such equations are solved for tangential components of the boundary. Lentr is large enough (≈5 × 10−3m) so that the flow can develop its laminar form, while pentr is solved so that the average velocity of the fluid through the boundary Γ1 meets U0. This implies that in the microscale, with the low Reynolds number used (≈0.1), the flow has a fully developed laminar profile on the inlet:
(2.14)Lentr∇t⋅[pI-η(∇tu+(∇tu)t)]=-npentr,∇t⋅u=0.Figure 3
Sketch of the physical situation simulated by a laminar inflow. Solid line is the actual computational domain while dashed line is a fictitious domain. The fluid, flowing from left to right, achieves a fully developed laminar profile on the inlet.Outlet (Normal Stress)
The boundary condition to be satisfied on the outlet is given by (2.15). This means that total stress in the tangential direction is zero. Under these circumstances pressure in the outlet is approximately zero:
(2.15)[-pI+η(∇u+(∇u)t)]n=0.Electrodes (No Slip)
Due to the microscale of the channel, fluid’s velocity is set to zero at these two walls. Equation (2.16) establishes this condition:
(2.16)u=0.


## 2.2.2. Charge Distribution

According to [13], it is assumed that there is no normal flux of anions at the electrodes,   at the inlet, nor at the outlet. It corresponds to the basic model which will be referred in this work as Case A: null exchange current density. Later on a relaxation on the condition for the anode will be allowed, and cations will be assumed to undergo reactions at a rate of an Arrhenius-like form. This will be called in this paper Case B: nonzero exchange current density. Section 3 addresses details for this modified case. Then the boundary conditions for Case A are stated as follows.Inlet, Outlet, and Electrodes (Convective Flux)
This condition states that the transport of ions due to diffusive or electric migration effects, perpendicular to the four boundaries, is negligible. Thus, (2.17)-(2.18) model this condition:
(2.17)n⋅(-D∇c1-wFc1∇ϕ)=0,(2.18)n⋅(-D∇c2+wFc2∇ϕ)=0.


## 2.2.3. Electrostatics

Inlet and outlet are considered to have zero electric charge because of the absence of electrodes. A very different situation happens on the horizontal frontiers. According to [14, 15], the existence of electrodes, the applied electric field on the vertical direction, and the flow of the electrolyte in the channel  lead—in the fuel cell microscale—to the well-known [12] phenomenon of the electrical double layer (EDL). It concerns two parallel layers of charge surrounding the electrodes. The first layer, called the Stern or compact layer, comprises ions adsorbed directly onto the electrode (wall) due to the host of chemical interactions. Its thickness is the order of one molecule diameter. The second layer, called the diffuse layer, is composed of ions attracted to the surface charge via the Coulomb force, electrically screening the first layer. The second layer is loosely associated to the wall, because it is made of free ions which move in the fluid (electrolyte) under the influence of electric attraction rather than being firmly anchored [12]. Its thickness is the order of 10−8 to 10−7m. The remaining part of the cell, which corresponds to almost all the channel height, is called the bulk. Figure 4 shows these components schematically without scales. The Gouy-Chapman-Stern theory [13, 15] states that ∇ϕ is a constant in the Stern layer and approaches zero in the diffuse layer. Then a linear drop voltage, that is, an electric displacement, can be assumed as appropriate boundary condition on electrodes. Details are presented next.Figure 4
Electrical double-layer (EDL) effect close to an electrode. The Stern layer is compact and comprises ions adsorbed directly onto the electrode. Its thickness is the order of one molecule diameter. The diffuse layer is composed by ions attracted to the surface via Coulomb-forces. Its thickness is the order of 10−8 to 10−7m. The bulk is the remaining part of the channel. Picture is drawn without scales.Inlet and Outlet (Zero Charge)
This condition assumes that the normal component of the electric field is zero and is given by(2.19)∇ϕ⋅n=0.Upper Electrode (Electric Displacement)
The normal component of the electric field on the anode is specified as stated in(2.20)∇ϕ⋅n=κ(ϕ0-ϕ),
where
(2.21)κ=1δγ,δ=ΔCrefCref,γ=λDh,Cref is the reference molar concentration, λD is the Debye length (thickness of the diffuse layer), ϕ0 the initial electric potential, and h is a half of the channel height. According to [16–18],  appropriate values for these parameters taken from experimental designs, are
(2.22)δ=1.121212×10-5,γ=2×10-5⟹κ=4.459×109.Lower Electrode (Electric Displacement)
In the case of the cathode, the normal component of the electric field on this boundary is specified as shown in(2.23)∇ϕ⋅n=-κ(ϕ0+ϕ).Therefore, the complete mathematical model for the membraneless fuel cell is given by the set of (2.9)–(2.13) to be satisfied in domain Ω, together with boundary conditions stated by (2.14)–(2.23) to be fulfilled on the corresponding boundaries.As extension of [8], present work is not limited to the case of electroneutrality (ρq=0) or electrochemical cells (u=0)  but, following proposal in [16], considers and gives solution to the complete complex system just described which models a more realistic situation. Furthermore, nonzero electrical reactions on the anode, that is, variations on   (2.11) and boundary condition (2.17), will be addressed as well.


## 2.3. Numerical Strategy for Solution

Since system (2.9)–(2.13) is highly nonlinear and time dependent, an effective numerical strategy must be implemented to integrate the included equations. Actually the system contains the following five unknowns: the bidimensional velocity field u, the pressure p, the charge concentrations c1 and c2, and the electric potential ϕ. The proposed numerical strategy is a robust combination of time integration by operator-splitting techniques [19, 20] and a finite element method for space integration [21].Then, ifΔt is the time discretization step and using the following notation(2.24)un≡u(x,nΔt),pn≡p(x,nΔt),c1n≡c1(x,nΔt),c2n≡c2(x,nΔt),ϕn≡ϕ(x,nΔt),
the numerical strategy can be summarized as follows.Given the initial valuesu0,p0,c10,c20, and ϕ0, find the solution at time n+1 solving (2.25)–(2.29), for n=0,1,2,…,(2.25)ρun+1-unΔt+ρ(un+1⋅∇)un=∇⋅[-pn+1I+η(∇un+(∇un)t)]-F(c1n-c2n)∇ϕn,(2.26)∇⋅un+1=0,(2.27)c1n+1-c1nΔt+un+1⋅∇c1n+1=∇⋅(D∇c1n+wFc1n∇ϕn),(2.28)c2n+1-c2nΔt+un+1⋅∇c2n+1=∇⋅(D∇c2n-wFc2n∇ϕn),(2.29)-ε0εr∇2ϕn+1=F(c1n+1-c2n+1).This algorithm quickly converges to a steady-state solution. According to [22], in the microscale of the fuel cell, it is enough to simulate at least 0.013 seconds in order to capture the dynamics of the phenomenon. In this work all of the simulations were run up to 1 second to assure a steady-state solution. Results are reported in the following section. Details on the FEM used for the numerical solutions are also included.


## 3. Results

In this section results of the numerical simulations for two cases of the modeled membraneless fuel cell, with and without electric reaction at the anode, are shown. Then some technical considerations about the FEM used in the numerical solution of mathematical models are summarized. The particular values for parameters used in both models, following experimental designs reported in [17, 18], are listed in Table 1.Since fluid dynamics is not pressure driven, even thoughp was calculated, it is not shown in results because it is irrelevant in present study.Results after applying algorithm (2.25)–(2.29) to solve problem (2.9)–(2.13) in domain Ω with corresponding boundary conditions (2.14)–(2.23) are shown next.


### 3.1. Case A: Null Exchange Current Density

The first case considered in this paper is the simulation of a membraneless fuel cell without electric reactions on electrodes. The mathematical model corresponds to PDE’s (2.9)–(2.13) satisfied in domain Ω with appropriate boundary conditions (2.14)–(2.23). Solving this system using the algorithm (2.25)–(2.29), the following results are obtained. All of the results are reported at time t=1 second, assuring a steady state solution.In order to depict the fuel cell behavior under different operating conditions, two different situations are considered within Case A. They differ in the initial conditions imposed to the electric potentialϕ. In the first situation, named as Case A.1, the initial electric potential is a function only on the y coordinate; that is, it varies on the direction perpendicular to the electrodes. This direction is parallel to the applied electric potential. Meanwhile in the second situation, named as Case A.2, the initial electric potential depends on both x and y coordinates. Results are next.


#### 3.1.1. Case A.1:ϕ(0)=ϕ0y

Figure5 shows the three-dimensional image of the velocity field u. Surface and height correspond to |u|. Since a fully developed laminar profile was imposed as boundary condition on the inlet and due to the low Reynolds number (0.1) used in the model, the flow maintains a laminar behavior all the way in the channel. It also reflects the absence of pressure at the outlet. Units are mm/s.Figure 5
Velocity fieldu measured in mm/s. Surface and height correspond to |u|. Boundary conditions and the low Reynolds number used in the model assure the laminar behavior observed.In Figure6 the electric potential ϕ is depicted. Since a symmetric voltage is applied on electrodes, same but positive and negative values are achieved on these boundaries. A linear behavior is observed consistent with the Gouy-Chapman-Stern theory [13, 15]. Moreover, it is zero in the central part of the cell, called the bulk, as that theory predicts. Units are volts.Figure 6
Electric potential fieldϕ, measured in volts. A linear drop voltage is observed vanishing in the middle (the bulk) of the cell, consistent with theoretical predictions.Electric conductivityσ defined by (2.6) is shown in Figure 7. It can be observed the effect of the EDL formed near to both electrodes. In order to get a better picture of what happens with electric conductivity behavior, a cross section of σ computed at the middle of the fuel cell, that is, maintaining the coordinate x=5 mm constant and varying y, is shown in Figure 8. The double-layer effect is appreciated better. This picture is also qualitatively consistent with results reported in [23] where a one-dimensional case was studied. Note that electric conductivity is constant in the bulk and nonzero, again in accordance with theoretical predictions [13, 15]. This result is also in agreement with the two-dimensional case addressed in [8], and the order of magnitude of σ is commensurable with that shown in [16]. Units are S/m.Figure 7
Electric conductivityσ corresponding to Case A.1 where the initial electric potential depends only on the y coordinate which is the direction of the applied voltage: ϕ(0)=ϕ0y.Figure 8
Cross section of the electric conductivity in Figure7; at x=5 mm, σ is constant but not zero in the bulk: ϕ(0)=ϕ0y.Electric charge densityρq defined by (2.7) is depicted in Figure 9. It can be appreciated the effect of the symmetric applied electric potential. Charge density is zero in the bulk; that is, electroneutrality is achieved in that part of the fuel cell, as it is established by theoretical predictions [13, 15]. A better appreciation of electric density behavior is gotten by the cross-section shown in Figure 10, computed at the same place as it was for σ. As it happened with conductivity, results for ρq qualitatively agree with results reported in the literature [8, 16, 23]. Units are C/m3.Figure 9
Electric charge densityρq corresponding to Case A.1 where the initial electric potential depends only on the y coordinate which is the direction of the applied voltage: ϕ(0)=ϕ0y.Figure 10
Cross-section of the electric charge density in Figure9;  at x=5 mm, ρq is zero in the bulk meaning electroneutrality is achieved in that zone: ϕ(0)=ϕ0y.Therefore, Case A.1 can be considered as the basis scenario. The membraneless fuel cell modeled under these conditions behaves according to theoretical postulates, and it is a good extension of previous more specific cases studied by some authors [8, 16, 23]. In the present work, it also serves as validation for both cases: the proposed model and the method used to get the numerical solution. Thus, simulations following this scheme can be considered confident and reliable. Cases A.2 and B, presented below, introduce modifications on the basis scenario and are solved subject to the same treatment. They are useful to start a mental map showing the behavior of a membraneless fuel cell.


#### 3.1.2. Case A.2:ϕ(0)=ϕ0xy

An interesting case from the mathematical and physical standpoint is a slight variation in the initial condition for the electric potential. Now it is assumed that the initial electric potential depends on both directions of the plane:ϕ(0)=ϕ0xy. Even though the introduced change seems to be small, the main effect is on the electric conductivity and on the charge density, as it will be shown next.Doing the same procedure as in Case A.1, the numerical solution of the modified model shows there are no major variations in fluid dynamics and in electric potential behavior, since plots foru and ϕ are practically the same as those shown in Figures 5 and 6, respectively. Then, it is not worth repeating figures.Nevertheless, a significant change occurs in the shape of electric conductivity. Figure11 shows σ for Case A.2. The effect of the EDL is again appreciated, and σ is negligible in almost a half of each electrode, near the side x=0. The cross-section for σ at x=5 mm is practically the same as that shown in Figure 8, then it is not worth to present it. Electric conductivity is constant but nonzero in the bulk. Difference between maxima and minima values is as small as 0.006 × 10−13.Figure 11
Electric conductivityσ corresponding to Case A.2 where the initial electric potential depends on both directions of the plane: ϕ(0)=ϕ0xy.Similar behavior appears in the plot of electric charge density, as shown in Figure12. There is negligible density in the electrodes’ half near the axis x=0. The cross-section of ρq shown in Figure 13 has the same shape as in the previous Case A.1, but its maximum intensity is two orders of magnitude less than ρq in the basis scenario (c.f. Figure 10). Again, electroneutrality is achieved in the bulk.Figure 12
Electric charge densityρq corresponding to Case A.2 where the initial electric potential depends on both directions of the plane: ϕ(0)=ϕ0xy.Figure 13
Cross-section of the electric charge density in Figure12; at x=5 mm, ρq is zero in the bulk meaning electroneutrality is achieved in that zone: ϕ(0)=ϕ0xy.As far as the author’s knowledge, there are no results like these for electric conductivity and charge density reported in the literature.


### 3.2. Case B: Nonzero Exchange Current Density

Another case considered in the paper  is the simulation of a membraneless fuel cell allowing electric reactions at one of the electrodes. Some information about these phenomena has appeared in the  specialized literature [24]. Case B concerns the situation when cations are subject to electric reactions at a rate AR stated by the charge transfer control theory [3] as shown in (3.1)AR=i0[e(αaF/RT)(ϕ-ϕ0)-e-(αcF/RT)(ϕ-ϕ0)],
where i0 is the anodic  exchange current density, measured in A/m2, αa and αc are the charge transfer coefficients, anodic and cathodic, respectively (and depend on the specific couple fuel-oxidant used in the experiment), F is the Faraday’s constant, R is the universal gas constant, and T is the absolute temperature.Assumingαa=αc=α that is, one step reaction  in which only one electron is transferred, it becomes into the hyperbolic reaction:(3.2)HAR=2i0sinh[αFRT(ϕ-ϕ0)].Equation (3.1) is known as the Butler-Volmer equation [14, 15] and represents an Arrhenius-like electrochemical reaction.In order to adjust the mathematical model to this situation, some changes must be introduced in the set of governing PDEs (2.9)–(2.13), in the corresponding boundary conditions (2.14)–(2.23) and in the numerical algorithm (2.25)–(2.29) as well. Specific arrangements are given next.Modified Equation for Charge Distribution
Equation (2.11) is replaced by (2.12′) where the hyperbolic Arrhenius reaction (oxidation) given by (3.2) is allowed in domain Ω:
2.12′∂c1∂t+∇⋅(-D∇c1-wFc1∇ϕ)=HAR-u⋅∇c1.Modified Boundary Condition on Anode
Also the convective flux boundary condition (2.17) on the anode is replaced by the flux condition (2.19′), where the flux is set to be as the hyperbolic Arrhenius reaction stated in   (3.2). The no-slip boundary condition (2.16) neglects the last term within parenthesis; nevertheless, it is a general formulation:
2.19′-n⋅(-D∇c1-wFc1∇ϕ+c1u)=HAR.
All of the remaining equations and boundary conditions remain the same as in the basis scenario. In order to maintain consistent initial values of the electric potentialϕ so that (2.12′) and (2.19′) are still valid for c1 since the right-hand side is nonzero (and ϕ dependent), the initial condition for the electric potential ϕ is now set to a constant ϕ0. Its value is given in Table 1. Thus, ∇ϕ vanishes letting HAR to drive the flux.Modified Step in Numerical Algorithm
Finally, the corresponding step in the numerical algorithm given by   (2.27) is replaced by  (2.29′), taking into consideration the following notation:
2.29′HARn≡HAR(x,nΔt)=2i0sinh⁡[αFRT(ϕn-ϕ0)],c1n+1-c1nΔt+un+1⋅∇c1n+1=∇⋅(D∇c1n+wFc1n∇ϕn)+HARn.
Therefore, introducing these three changes (2.12′), (2.19′), and (2.29′) in the set of governing PDEs, boundary conditions and numerical algorithm, resp.), the following results are obtained.
As it happened in Case A.2, there are no major variations in fluid dynamics and in electric potential, since plots foru and ϕ are practically the same as those shown in Figures 5 and 6, respectively. There is no worth presenting such graphs here.
Figure14 shows the electric conductivity for Case B. It can be appreciated the effect of the major electric activity on the anode. The order of magnitude of σ is the same as in the basis scenario. Now conductivity is not constant in the bulk and approaches asymptotically zero near the cathode. Conductivity behavior is dominated by the Arrhenius-like reaction on the anode.
The electric charge densityρq for Case B is shown in Figure 15. Again the influence of the Arrhenius-like reaction allowed on the anode is appreciated. It dominates density behavior. The order of magnitude for density is the same as in the basis scenario. Charge density is not zero in the bulk but approaches this value near the cathode.
As far as the author’s knowledge, this is the first time results like these for electric conductivity and charge density are reported.
In order to deepen  the behavior analysis of this membraneless fuel cell, variations in the anodic exchange current densityi0 were introduced to study their effect on σ and ρq. Letting i0 vary from 0.1 to 10,000 A/m2 and solving the model to get the corresponding solutions for electric conductivity and charge density, it was observed that both shapes remained basically the same; that is, they have the same qualitative behavior. The main effect occurred in the magnitude of those quantities and in the range from the minima to the maxima values achieved in the domain Ω. Figures 16 and 17 show the dependence of σ and ρq, respectively, when doing this sensitivity analysis on the anodic  exchange current density i0. It can be appreciated that the greater the current density, the greater the maximum achieved value and the broader the difference between the minima and the maxima values. Also a linear growth is observed in the set of maxima values. Minima values are almost constant at zero. This is true for both electric conductivity and charge density.
A slope of 2.7 × 10−9 is observed for the electric conductivity maxima values’ growth. The corresponding value for electric charge density is 3.5 × 10−2. These are also novelty results.Figure 14
Electric conductivityσ corresponding to Case B. An anionic, hyperbolic Arrhenius’ reaction is allowed in the domain Ω and on the anode. The initial condition for the electric potential is a constant ϕ(0)=ϕ0.Figure 15
Charge-densityρq corresponding to Case B. An anionic, hyperbolic Arrhenius’ reaction is allowed in the domain Ω and on the anode. The initial condition for the electric potential is a constant ϕ(0)=ϕ0.Figure 16
Dependence of electric conductivityσ on the anodic exchange current density i0. There is a linear growth of maxima values achieved in the domain Ω.Figure 17
Dependence of charge densityρq on the anodic  exchange current density i0. There is a linear growth of maxima values achieved in the domain Ω.


### 3.3. Considerations on Finite-Element Method

Some technical considerations about the FEM used to get the numerical solution of mathematical models for membraneless fuel cells depicted in Cases A.1, A.2, and B are summarized in Table2. Models were solved using a Latitude D630 laptop having 2 GB RAM and a 2 GHz processor. Quadratic elements were used to approximate each scalar component of velocity u, both charge molar concentrations c1 and c2, and the electric potential ϕ as well. Linear elements were used to approximate pressure p.Table 2
Technical considerations on FEM used to find out numerical solutions of mathematical models for membraneless fuel cells depicted in Cases A.1, A.2, and B. A refined mesh was needed in the last Case to numerically capture the Arrhenius-like electric reactions on the anode.
CaseNo. of mesh elementsNo. of equations solvedCPU time (sec)A.11,02217,04569.874A.21,02217,04573.048B1,70032,011105.39Time domain interval for simulation was set to[0,1] seconds taking a time discretization step of Δt=0.001. In order to numerically capture EDL effects, a graduated structured mesh was used with smaller element size in the order of 1 × 10−6 m next to both electrodes. But in Case B, it was needed a mesh refinement in the horizontal direction to make visible anodic Arrhenius’ reactions in the given space-time domain.To verify accuracy of the method to maintain mass conservation as a true condition, numerical divergence of the solution for the velocity fieldu was computed. In the three reported cases, it was in the order of 1 × 10−4. Given the complexity of the model, it can be considered as an acceptable value.


## 3.1. Case A: Null Exchange Current Density

The first case considered in this paper is the simulation of a membraneless fuel cell without electric reactions on electrodes. The mathematical model corresponds to PDE’s (2.9)–(2.13) satisfied in domain Ω with appropriate boundary conditions (2.14)–(2.23). Solving this system using the algorithm (2.25)–(2.29), the following results are obtained. All of the results are reported at time t=1 second, assuring a steady state solution.In order to depict the fuel cell behavior under different operating conditions, two different situations are considered within Case A. They differ in the initial conditions imposed to the electric potentialϕ. In the first situation, named as Case A.1, the initial electric potential is a function only on the y coordinate; that is, it varies on the direction perpendicular to the electrodes. This direction is parallel to the applied electric potential. Meanwhile in the second situation, named as Case A.2, the initial electric potential depends on both x and y coordinates. Results are next.


### 3.1.1. Case A.1:ϕ(0)=ϕ0y

Figure5 shows the three-dimensional image of the velocity field u. Surface and height correspond to |u|. Since a fully developed laminar profile was imposed as boundary condition on the inlet and due to the low Reynolds number (0.1) used in the model, the flow maintains a laminar behavior all the way in the channel. It also reflects the absence of pressure at the outlet. Units are mm/s.Figure 5
Velocity fieldu measured in mm/s. Surface and height correspond to |u|. Boundary conditions and the low Reynolds number used in the model assure the laminar behavior observed.In Figure6 the electric potential ϕ is depicted. Since a symmetric voltage is applied on electrodes, same but positive and negative values are achieved on these boundaries. A linear behavior is observed consistent with the Gouy-Chapman-Stern theory [13, 15]. Moreover, it is zero in the central part of the cell, called the bulk, as that theory predicts. Units are volts.Figure 6
Electric potential fieldϕ, measured in volts. A linear drop voltage is observed vanishing in the middle (the bulk) of the cell, consistent with theoretical predictions.Electric conductivityσ defined by (2.6) is shown in Figure 7. It can be observed the effect of the EDL formed near to both electrodes. In order to get a better picture of what happens with electric conductivity behavior, a cross section of σ computed at the middle of the fuel cell, that is, maintaining the coordinate x=5 mm constant and varying y, is shown in Figure 8. The double-layer effect is appreciated better. This picture is also qualitatively consistent with results reported in [23] where a one-dimensional case was studied. Note that electric conductivity is constant in the bulk and nonzero, again in accordance with theoretical predictions [13, 15]. This result is also in agreement with the two-dimensional case addressed in [8], and the order of magnitude of σ is commensurable with that shown in [16]. Units are S/m.Figure 7
Electric conductivityσ corresponding to Case A.1 where the initial electric potential depends only on the y coordinate which is the direction of the applied voltage: ϕ(0)=ϕ0y.Figure 8
Cross section of the electric conductivity in Figure7; at x=5 mm, σ is constant but not zero in the bulk: ϕ(0)=ϕ0y.Electric charge densityρq defined by (2.7) is depicted in Figure 9. It can be appreciated the effect of the symmetric applied electric potential. Charge density is zero in the bulk; that is, electroneutrality is achieved in that part of the fuel cell, as it is established by theoretical predictions [13, 15]. A better appreciation of electric density behavior is gotten by the cross-section shown in Figure 10, computed at the same place as it was for σ. As it happened with conductivity, results for ρq qualitatively agree with results reported in the literature [8, 16, 23]. Units are C/m3.Figure 9
Electric charge densityρq corresponding to Case A.1 where the initial electric potential depends only on the y coordinate which is the direction of the applied voltage: ϕ(0)=ϕ0y.Figure 10
Cross-section of the electric charge density in Figure9;  at x=5 mm, ρq is zero in the bulk meaning electroneutrality is achieved in that zone: ϕ(0)=ϕ0y.Therefore, Case A.1 can be considered as the basis scenario. The membraneless fuel cell modeled under these conditions behaves according to theoretical postulates, and it is a good extension of previous more specific cases studied by some authors [8, 16, 23]. In the present work, it also serves as validation for both cases: the proposed model and the method used to get the numerical solution. Thus, simulations following this scheme can be considered confident and reliable. Cases A.2 and B, presented below, introduce modifications on the basis scenario and are solved subject to the same treatment. They are useful to start a mental map showing the behavior of a membraneless fuel cell.


### 3.1.2. Case A.2:ϕ(0)=ϕ0xy

An interesting case from the mathematical and physical standpoint is a slight variation in the initial condition for the electric potential. Now it is assumed that the initial electric potential depends on both directions of the plane:ϕ(0)=ϕ0xy. Even though the introduced change seems to be small, the main effect is on the electric conductivity and on the charge density, as it will be shown next.Doing the same procedure as in Case A.1, the numerical solution of the modified model shows there are no major variations in fluid dynamics and in electric potential behavior, since plots foru and ϕ are practically the same as those shown in Figures 5 and 6, respectively. Then, it is not worth repeating figures.Nevertheless, a significant change occurs in the shape of electric conductivity. Figure11 shows σ for Case A.2. The effect of the EDL is again appreciated, and σ is negligible in almost a half of each electrode, near the side x=0. The cross-section for σ at x=5 mm is practically the same as that shown in Figure 8, then it is not worth to present it. Electric conductivity is constant but nonzero in the bulk. Difference between maxima and minima values is as small as 0.006 × 10−13.Figure 11
Electric conductivityσ corresponding to Case A.2 where the initial electric potential depends on both directions of the plane: ϕ(0)=ϕ0xy.Similar behavior appears in the plot of electric charge density, as shown in Figure12. There is negligible density in the electrodes’ half near the axis x=0. The cross-section of ρq shown in Figure 13 has the same shape as in the previous Case A.1, but its maximum intensity is two orders of magnitude less than ρq in the basis scenario (c.f. Figure 10). Again, electroneutrality is achieved in the bulk.Figure 12
Electric charge densityρq corresponding to Case A.2 where the initial electric potential depends on both directions of the plane: ϕ(0)=ϕ0xy.Figure 13
Cross-section of the electric charge density in Figure12; at x=5 mm, ρq is zero in the bulk meaning electroneutrality is achieved in that zone: ϕ(0)=ϕ0xy.As far as the author’s knowledge, there are no results like these for electric conductivity and charge density reported in the literature.


## 3.1.1. Case A.1:ϕ(0)=ϕ0y

Figure5 shows the three-dimensional image of the velocity field u. Surface and height correspond to |u|. Since a fully developed laminar profile was imposed as boundary condition on the inlet and due to the low Reynolds number (0.1) used in the model, the flow maintains a laminar behavior all the way in the channel. It also reflects the absence of pressure at the outlet. Units are mm/s.Figure 5
Velocity fieldu measured in mm/s. Surface and height correspond to |u|. Boundary conditions and the low Reynolds number used in the model assure the laminar behavior observed.In Figure6 the electric potential ϕ is depicted. Since a symmetric voltage is applied on electrodes, same but positive and negative values are achieved on these boundaries. A linear behavior is observed consistent with the Gouy-Chapman-Stern theory [13, 15]. Moreover, it is zero in the central part of the cell, called the bulk, as that theory predicts. Units are volts.Figure 6
Electric potential fieldϕ, measured in volts. A linear drop voltage is observed vanishing in the middle (the bulk) of the cell, consistent with theoretical predictions.Electric conductivityσ defined by (2.6) is shown in Figure 7. It can be observed the effect of the EDL formed near to both electrodes. In order to get a better picture of what happens with electric conductivity behavior, a cross section of σ computed at the middle of the fuel cell, that is, maintaining the coordinate x=5 mm constant and varying y, is shown in Figure 8. The double-layer effect is appreciated better. This picture is also qualitatively consistent with results reported in [23] where a one-dimensional case was studied. Note that electric conductivity is constant in the bulk and nonzero, again in accordance with theoretical predictions [13, 15]. This result is also in agreement with the two-dimensional case addressed in [8], and the order of magnitude of σ is commensurable with that shown in [16]. Units are S/m.Figure 7
Electric conductivityσ corresponding to Case A.1 where the initial electric potential depends only on the y coordinate which is the direction of the applied voltage: ϕ(0)=ϕ0y.Figure 8
Cross section of the electric conductivity in Figure7; at x=5 mm, σ is constant but not zero in the bulk: ϕ(0)=ϕ0y.Electric charge densityρq defined by (2.7) is depicted in Figure 9. It can be appreciated the effect of the symmetric applied electric potential. Charge density is zero in the bulk; that is, electroneutrality is achieved in that part of the fuel cell, as it is established by theoretical predictions [13, 15]. A better appreciation of electric density behavior is gotten by the cross-section shown in Figure 10, computed at the same place as it was for σ. As it happened with conductivity, results for ρq qualitatively agree with results reported in the literature [8, 16, 23]. Units are C/m3.Figure 9
Electric charge densityρq corresponding to Case A.1 where the initial electric potential depends only on the y coordinate which is the direction of the applied voltage: ϕ(0)=ϕ0y.Figure 10
Cross-section of the electric charge density in Figure9;  at x=5 mm, ρq is zero in the bulk meaning electroneutrality is achieved in that zone: ϕ(0)=ϕ0y.Therefore, Case A.1 can be considered as the basis scenario. The membraneless fuel cell modeled under these conditions behaves according to theoretical postulates, and it is a good extension of previous more specific cases studied by some authors [8, 16, 23]. In the present work, it also serves as validation for both cases: the proposed model and the method used to get the numerical solution. Thus, simulations following this scheme can be considered confident and reliable. Cases A.2 and B, presented below, introduce modifications on the basis scenario and are solved subject to the same treatment. They are useful to start a mental map showing the behavior of a membraneless fuel cell.


## 3.1.2. Case A.2:ϕ(0)=ϕ0xy

An interesting case from the mathematical and physical standpoint is a slight variation in the initial condition for the electric potential. Now it is assumed that the initial electric potential depends on both directions of the plane:ϕ(0)=ϕ0xy. Even though the introduced change seems to be small, the main effect is on the electric conductivity and on the charge density, as it will be shown next.Doing the same procedure as in Case A.1, the numerical solution of the modified model shows there are no major variations in fluid dynamics and in electric potential behavior, since plots foru and ϕ are practically the same as those shown in Figures 5 and 6, respectively. Then, it is not worth repeating figures.Nevertheless, a significant change occurs in the shape of electric conductivity. Figure11 shows σ for Case A.2. The effect of the EDL is again appreciated, and σ is negligible in almost a half of each electrode, near the side x=0. The cross-section for σ at x=5 mm is practically the same as that shown in Figure 8, then it is not worth to present it. Electric conductivity is constant but nonzero in the bulk. Difference between maxima and minima values is as small as 0.006 × 10−13.Figure 11
Electric conductivityσ corresponding to Case A.2 where the initial electric potential depends on both directions of the plane: ϕ(0)=ϕ0xy.Similar behavior appears in the plot of electric charge density, as shown in Figure12. There is negligible density in the electrodes’ half near the axis x=0. The cross-section of ρq shown in Figure 13 has the same shape as in the previous Case A.1, but its maximum intensity is two orders of magnitude less than ρq in the basis scenario (c.f. Figure 10). Again, electroneutrality is achieved in the bulk.Figure 12
Electric charge densityρq corresponding to Case A.2 where the initial electric potential depends on both directions of the plane: ϕ(0)=ϕ0xy.Figure 13
Cross-section of the electric charge density in Figure12; at x=5 mm, ρq is zero in the bulk meaning electroneutrality is achieved in that zone: ϕ(0)=ϕ0xy.As far as the author’s knowledge, there are no results like these for electric conductivity and charge density reported in the literature.


## 3.2. Case B: Nonzero Exchange Current Density

Another case considered in the paper  is the simulation of a membraneless fuel cell allowing electric reactions at one of the electrodes. Some information about these phenomena has appeared in the  specialized literature [24]. Case B concerns the situation when cations are subject to electric reactions at a rate AR stated by the charge transfer control theory [3] as shown in (3.1)AR=i0[e(αaF/RT)(ϕ-ϕ0)-e-(αcF/RT)(ϕ-ϕ0)],
where i0 is the anodic  exchange current density, measured in A/m2, αa and αc are the charge transfer coefficients, anodic and cathodic, respectively (and depend on the specific couple fuel-oxidant used in the experiment), F is the Faraday’s constant, R is the universal gas constant, and T is the absolute temperature.Assumingαa=αc=α that is, one step reaction  in which only one electron is transferred, it becomes into the hyperbolic reaction:(3.2)HAR=2i0sinh[αFRT(ϕ-ϕ0)].Equation (3.1) is known as the Butler-Volmer equation [14, 15] and represents an Arrhenius-like electrochemical reaction.In order to adjust the mathematical model to this situation, some changes must be introduced in the set of governing PDEs (2.9)–(2.13), in the corresponding boundary conditions (2.14)–(2.23) and in the numerical algorithm (2.25)–(2.29) as well. Specific arrangements are given next.Modified Equation for Charge Distribution
Equation (2.11) is replaced by (2.12′) where the hyperbolic Arrhenius reaction (oxidation) given by (3.2) is allowed in domain Ω:
2.12′∂c1∂t+∇⋅(-D∇c1-wFc1∇ϕ)=HAR-u⋅∇c1.Modified Boundary Condition on Anode
Also the convective flux boundary condition (2.17) on the anode is replaced by the flux condition (2.19′), where the flux is set to be as the hyperbolic Arrhenius reaction stated in   (3.2). The no-slip boundary condition (2.16) neglects the last term within parenthesis; nevertheless, it is a general formulation:
2.19′-n⋅(-D∇c1-wFc1∇ϕ+c1u)=HAR.
All of the remaining equations and boundary conditions remain the same as in the basis scenario. In order to maintain consistent initial values of the electric potentialϕ so that (2.12′) and (2.19′) are still valid for c1 since the right-hand side is nonzero (and ϕ dependent), the initial condition for the electric potential ϕ is now set to a constant ϕ0. Its value is given in Table 1. Thus, ∇ϕ vanishes letting HAR to drive the flux.Modified Step in Numerical Algorithm
Finally, the corresponding step in the numerical algorithm given by   (2.27) is replaced by  (2.29′), taking into consideration the following notation:
2.29′HARn≡HAR(x,nΔt)=2i0sinh⁡[αFRT(ϕn-ϕ0)],c1n+1-c1nΔt+un+1⋅∇c1n+1=∇⋅(D∇c1n+wFc1n∇ϕn)+HARn.
Therefore, introducing these three changes (2.12′), (2.19′), and (2.29′) in the set of governing PDEs, boundary conditions and numerical algorithm, resp.), the following results are obtained.
As it happened in Case A.2, there are no major variations in fluid dynamics and in electric potential, since plots foru and ϕ are practically the same as those shown in Figures 5 and 6, respectively. There is no worth presenting such graphs here.
Figure14 shows the electric conductivity for Case B. It can be appreciated the effect of the major electric activity on the anode. The order of magnitude of σ is the same as in the basis scenario. Now conductivity is not constant in the bulk and approaches asymptotically zero near the cathode. Conductivity behavior is dominated by the Arrhenius-like reaction on the anode.
The electric charge densityρq for Case B is shown in Figure 15. Again the influence of the Arrhenius-like reaction allowed on the anode is appreciated. It dominates density behavior. The order of magnitude for density is the same as in the basis scenario. Charge density is not zero in the bulk but approaches this value near the cathode.
As far as the author’s knowledge, this is the first time results like these for electric conductivity and charge density are reported.
In order to deepen  the behavior analysis of this membraneless fuel cell, variations in the anodic exchange current densityi0 were introduced to study their effect on σ and ρq. Letting i0 vary from 0.1 to 10,000 A/m2 and solving the model to get the corresponding solutions for electric conductivity and charge density, it was observed that both shapes remained basically the same; that is, they have the same qualitative behavior. The main effect occurred in the magnitude of those quantities and in the range from the minima to the maxima values achieved in the domain Ω. Figures 16 and 17 show the dependence of σ and ρq, respectively, when doing this sensitivity analysis on the anodic  exchange current density i0. It can be appreciated that the greater the current density, the greater the maximum achieved value and the broader the difference between the minima and the maxima values. Also a linear growth is observed in the set of maxima values. Minima values are almost constant at zero. This is true for both electric conductivity and charge density.
A slope of 2.7 × 10−9 is observed for the electric conductivity maxima values’ growth. The corresponding value for electric charge density is 3.5 × 10−2. These are also novelty results.Figure 14
Electric conductivityσ corresponding to Case B. An anionic, hyperbolic Arrhenius’ reaction is allowed in the domain Ω and on the anode. The initial condition for the electric potential is a constant ϕ(0)=ϕ0.Figure 15
Charge-densityρq corresponding to Case B. An anionic, hyperbolic Arrhenius’ reaction is allowed in the domain Ω and on the anode. The initial condition for the electric potential is a constant ϕ(0)=ϕ0.Figure 16
Dependence of electric conductivityσ on the anodic exchange current density i0. There is a linear growth of maxima values achieved in the domain Ω.Figure 17
Dependence of charge densityρq on the anodic  exchange current density i0. There is a linear growth of maxima values achieved in the domain Ω.


## 3.3. Considerations on Finite-Element Method

Some technical considerations about the FEM used to get the numerical solution of mathematical models for membraneless fuel cells depicted in Cases A.1, A.2, and B are summarized in Table2. Models were solved using a Latitude D630 laptop having 2 GB RAM and a 2 GHz processor. Quadratic elements were used to approximate each scalar component of velocity u, both charge molar concentrations c1 and c2, and the electric potential ϕ as well. Linear elements were used to approximate pressure p.Table 2
Technical considerations on FEM used to find out numerical solutions of mathematical models for membraneless fuel cells depicted in Cases A.1, A.2, and B. A refined mesh was needed in the last Case to numerically capture the Arrhenius-like electric reactions on the anode.
CaseNo. of mesh elementsNo. of equations solvedCPU time (sec)A.11,02217,04569.874A.21,02217,04573.048B1,70032,011105.39Time domain interval for simulation was set to[0,1] seconds taking a time discretization step of Δt=0.001. In order to numerically capture EDL effects, a graduated structured mesh was used with smaller element size in the order of 1 × 10−6 m next to both electrodes. But in Case B, it was needed a mesh refinement in the horizontal direction to make visible anodic Arrhenius’ reactions in the given space-time domain.To verify accuracy of the method to maintain mass conservation as a true condition, numerical divergence of the solution for the velocity fieldu was computed. In the three reported cases, it was in the order of 1 × 10−4. Given the complexity of the model, it can be considered as an acceptable value.


## 4. Discussion

Analyzed Cases A.1, A.2, and B can serve to start building up a mental map of the behavior of membraneless fuel cells. Having numerical scheme (2.25)–(2.29) and proper variations available, it becomes a useful tool to explore certain behaviors of this type of cells under different conditions. Several of these conditions were studied in present work and main findings of this research are presented in this section, together with some remarkable facts observed in the mentioned cases. Even though generalizations from particular cases are invalid, these observations/speculations could serve as guidelines to orient/confirm theoretical future research in order to reach a better understanding of this phenomenon. Similar boundary conditions must be remained.Velocity fieldu and electric potential ϕ are practically the same in the three cases. It can be speculated,  fluid dynamics and electric displacements are dominant; thus, u and ϕ will behave as in the basis scenario. Therefore, it is expected a laminar flow all the way in the channel and a linear voltage drop between symmetric electric displacements.Since (2.2) has been formulated in terms of the total stress tensor, it allows the possibility of considering fluids with variable viscosity η, that is, non-Newtonian fluids. Influence of electric charge molar concentrations c1 and c2 on viscosity in the form stated by   (4.1) is usually observed in solutions of larger molecules [25].(4.1)η=η0[1+β(c1+c2)2],
where β is a parameter used for sensitivity analysis and η0 can be set to 0.001 Pa·s, as in the basis scenario. Note that the term within parenthesis is in the order of 1 × 10−20. Letting β vary from 0.5 × 1020 to 1,000 × 1020, it was observed that solutions for u, ϕ, σ, and ρq were the same all the range in the three cases. Then it can be speculated that viscosity of fluids does not influence the membraneless fuel cell behavior.The effect of different initial conditions of the electric charge molar concentrations on the solution of the model was also analyzed. There were considered three different initial values forc1 and c2, as follows: c0,c0x, and c0xy. Value for c0 is shown in Table 1. Again, it was observed the same solution for u, ϕ, σ, and ρq in the three cases. Then, it is expected that the initial condition of charge concentrations does not influence the behavior of the membraneless fuel cell.After reviewing Cases A.1 and A.2, the following can be said. In absence of electrochemical reactions on the electrodes, it is expected that the electric conductivity will have a “W” shape-cross section like that shown in Figure8: flat on the bulk and with high positive gradients near electrodes. Also, charge density cross-section is expected to behave as Figure 10 shows: zero on the bulk and with high opposite sign gradients near electrodes. The high gradients are the effect of the EDL next to anode and cathode.If electrochemical reactions are allowed on at least one of the electrodes, special care must be placed in the establishment of initial conditions for the electric potentialϕ. Otherwise, (2.12′) or (2.19′) might not be satisfied. For instance, if ϕ(0) is set to ϕ0x or ϕ0xy in Case B, inconsistent initial values are found and the mathematical system cannot be solved. Thus, simple initial conditions for ϕ, as constants, are suitable in membraneless fuel cells having electrochemical activity on the electrodes.Finally, the analysis of the last two figures in Section3 allows speculating the following. There is a linear dependence between maxima values of electric conductivity σ and charge density ρq on the exchange current density i0, while minima values are almost constants about zero in membraneless fuel cells with electrochemical activity on one electrode.All mentioned features for membraneless fuel cells are plausible and deserve attention when planning design and manufacturing. Some theoretical statements as the Gouy-Chapman-Stern [15] or charge transfer control theory [3] have been confirmed and/or somehow extended.


## 5. Conclusions

Membraneless fuel cells are interesting and useful examples of MEMS because of their applications in engineering and life sciences, and are usually essential components of labs-on-a-chip. There is intense research to explain fundamental physical and electrochemical processes involved in the operation of this type of cells. This work is part of this effort.Main contribution of this paper lies in three subjects. First it is the design of a physical and mathematical model by means of a set of PDEs with appropriate boundary conditions, to depict the behavior of a membraneless fuel cell. Second is the implementation of a robust numerical algorithm to solve that model; it could be used in subsequent studies. Third is the characterization of the behavior of a membraneless fuel cell in terms of several facts like  conditions to get dominant velocityu and electric potential ϕ,  independence of u, ϕ, σ, and ρq with respect to variable viscosity η and initial conditions for charge molar concentrations c1 and c2, general shape for cross-sections of σ and ρq, awareness in the establishment of consistent initial conditions for ϕ with the Arrhenius-like reactions on anode, and linear dependence of σ and ρq maxima values on the exchange current density when electrochemical reactions are allowed on anode. This work is consistent with Gouy-Chapman-Stern’s theory [15] and charge transfer control theory [3]. In consequence a better understanding of the phenomenon has been achieved.Other contributions of this work lie in the completion and extension of previous results reported in [8, 16], where more specific problems were addressed.Applications of present study are devices based on microfluidics. Examples can be found in labs-on-a-chip used in medicine  or sensors and actuators for agronomy purposes. Contributions of this study can also be useful to optimize design and manufacturing of fuel cells.The characterization of the behavior of membraneless fuel cells is far to be complete, therefore complementary research is necessary. Future work can be oriented in several directions like the following. Allow a general Arrhenius-like reaction on the anode, not necessarily a hyperbolic one. Consider reactions on both electrodes and/or allow other types of reactions. Include the study of other geometries that favor the separation of reactants along the cell. Extend used geometry to the three-dimensional case. Include some other fuel/oxidant combinations.A missing consideration in present study is the creation of products resulting from oxidoreduction reactions. The inclusion of this phenomenon could lead to more realistic simulations. Nonetheless, it would imply the consideration of proper chemical laws, becoming the mathematical model more complex and difficult to handle.



---

*Source: 695167-2012-02-09.xml*