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.

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].

Figure 1 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.

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 Figure 1, 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 (๐‘งยฑ), then1=๐‘ง+(fuel)=๐‘ง+(oxidant)=โˆ’๐‘งโˆ’(fuel)=โˆ’๐‘งโˆ’(oxidant).(2.1)

Furthermore, diffusivity (๐ท) and mobility (๐‘ค) 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.

This 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 ๐‘ฅ๐‘ฆ, 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.

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, ๐ฎ the velocity field, ๐‘ the pressure, and ๐Ÿ 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 ๐ฎ is a time-dependent variable:๐œŒ๐œ•๐ฎ๎€บ๎€ท๐œ•๐‘ก+๐œŒ(๐ฎโ‹…โˆ‡)๐ฎ=โˆ‡โ‹…โˆ’๐‘๐ˆ+๐œ‚โˆ‡๐ฎ+(โˆ‡๐ฎ)๐‘ก๎€ธ๎€ป+๐Ÿ,(2.2)โˆ‡โ‹…๐ฎ=0.(2.3)

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). ๐ท1,2 are the diffusion coefficients, ๐‘ค+,โˆ’ are the mobility coefficients, ๐‘ง+,โˆ’ are the valence numbers as stated in (2.1), ๐น is the Faradayโ€™s constant, ๐œ™ is the electric potential, and ๐‘1,2 are the charge molar concentrations of reactants (positive and negative, resp.). Note that these concentrations are time-dependent variables ๐œ•๐‘1๐œ•๐‘ก+๐ฎโ‹…โˆ‡๐‘1๎€ท๐ท=โˆ‡โ‹…1โˆ‡๐‘1+๐‘ค+๐‘ง+๐น๐‘1๎€ธ,โˆ‡๐œ™๐œ•๐‘2๐œ•๐‘ก+๐ฎโ‹…โˆ‡๐‘2๎€ท๐ท=โˆ‡โ‹…2โˆ‡๐‘2+๐‘คโˆ’๐‘งโˆ’๐น๐‘2๎€ธ.โˆ‡๐œ™(2.4)

2.1.3. Electrostaticsโ€™s (Poisson Equation)

Based on electromagnetic theory [11] and assuming there is no a magnetic field and the electric field ๐„ is quasi-static, Maxwellโ€™s equations are โˆ‡ร—๐„=0; hence, ๐„=โˆ’โˆ‡๐œ™. Then, the Gauss law โˆ‡โ‹…๐„=๐œŒ๐‘ž/๐œ€0๐œ€๐‘Ÿ, where ๐œŒ๐‘ž is charge density and ๐œ€๐‘Ÿ is the relative electric permittivity of the electrolyte, implies the electric potential ๐œ™ must satisfy the Poisson equation (2.5):โˆ’๐œ€0๐œ€๐‘Ÿ=โˆ‡2๐œ™=๐œŒ๐‘ž.(2.5)

On the other hand, there are two important physical quantities for a fuel cell: electric conductivity ๐œŽ and charge density ๐œŒ๐‘ž, 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๐‘ค๎€ท๐‘1+๐‘2๎€ธ,๐œŒ(2.6)๐‘ž๎€ท๐‘=๐น1โˆ’๐‘2๎€ธ.(2.7)

Furthermore, the external force ๐Ÿ in (2.2) is the electromagnetic force given by๐Ÿ=๐œŒ๐‘ž๎€ท๐‘๐„=โˆ’๐น1โˆ’๐‘2๎€ธโˆ‡๐œ™.(2.8)

Therefore, recalling ๐ท1=๐ท2=๐ท and ๐‘ค+=๐‘คโˆ’=๐‘ค 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 ฮฉ๐œŒ๐œ•๐ฎ๎€บ๎€ท๐œ•๐‘ก+๐œŒ(๐ฎโ‹…โˆ‡)๐ฎ=โˆ‡โ‹…โˆ’๐‘๐ˆ+๐œ‚โˆ‡๐ฎ+(โˆ‡๐ฎ)๐‘ก๎€ท๐‘๎€ธ๎€ปโˆ’๐น1โˆ’๐‘2๎€ธโˆ‡๐œ™,(2.9)โˆ‡โ‹…๐ฎ=0,(2.10)๐œ•๐‘1๐œ•๐‘ก+๐ฎโ‹…โˆ‡๐‘1๎€ท=โˆ‡โ‹…๐ทโˆ‡๐‘1+๐‘ค๐น๐‘1๎€ธโˆ‡๐œ™,(2.11)๐œ•๐‘2๐œ•๐‘ก+๐ฎโ‹…โˆ‡๐‘2๎€ท=โˆ‡โ‹…๐ทโˆ‡๐‘2โˆ’๐‘ค๐น๐‘2๎€ธโˆ‡๐œ™,(2.12)โˆ’๐œ€0๐œ€๐‘Ÿโˆ‡2๎€ท๐‘๐œ™=๐น1โˆ’๐‘2๎€ธ.(2.13)

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, ๐ง 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 ๐‘ก emphasizes that such equations are solved for tangential components of the boundary. ๐ฟentr is large enough (โ‰ˆ5 ร— 10โˆ’3โ€‰m) so that the flow can develop its laminar form, while ๐‘entr 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: ๐ฟentrโˆ‡๐‘กโ‹…๎‚ƒ๎‚€โˆ‡๐‘๐ˆโˆ’๐œ‚๐‘ก๎€ทโˆ‡๐ฎ+๐‘ก๐ฎ๎€ธ๐‘ก๎‚๎‚„=โˆ’๐ง๐‘entr,โˆ‡๐‘กโ‹…๐ฎ=0.(2.14)

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: ๎€บ๎€ทโˆ’๐‘๐ˆ+๐œ‚โˆ‡๐ฎ+(โˆ‡๐ฎ)๐‘ก๎€ธ๎€ป๐ง=0.(2.15)

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: ๐ฎ=0.(2.16)

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: ๎€ท๐งโ‹…โˆ’๐ทโˆ‡๐‘1โˆ’๐‘ค๐น๐‘1๎€ธ๎€ทโˆ‡๐œ™=0,(2.17)๐งโ‹…โˆ’๐ทโˆ‡๐‘2+๐‘ค๐น๐‘2๎€ธโˆ‡๐œ™=0.(2.18)

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โˆ’7โ€‰m. 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.

Inlet and Outlet (Zero Charge)
This condition assumes that the normal component of the electric field is zero and is given by โˆ‡๐œ™โ‹…๐ง=0.(2.19)

Upper Electrode (Electric Displacement)
The normal component of the electric field on the anode is specified as stated in ๎€ท๐œ™โˆ‡๐œ™โ‹…๐ง=๐œ…0๎€ธโˆ’๐œ™,(2.20) where 1๐œ…=๐›ฟ๐›พ,๐›ฟ=ฮ”๐ถref๐ถref๐œ†,๐›พ=๐ทโ„Ž,(2.21)๐ถref is the reference molar concentration, ๐œ†๐ท is the Debye length (thickness of the diffuse layer), ๐œ™0 the initial electric potential, and โ„Ž is a half of the channel height. According to [16โ€“18], appropriate values for these parameters taken from experimental designs, are ๐›ฟ=1.121212ร—10โˆ’5,๐›พ=2ร—10โˆ’5โŸน๐œ…=4.459ร—109.(2.22)

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 ๎€ท๐œ™โˆ‡๐œ™โ‹…๐ง=โˆ’๐œ…0๎€ธ+๐œ™.(2.23)

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 (๐œŒ๐‘ž=0) or electrochemical cells (๐ฎ=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 ๐ฎ, the pressure ๐‘, the charge concentrations ๐‘1 and ๐‘2, 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 ฮ”๐‘ก is the time discretization step and using the following notation๐ฎ๐‘›โ‰ก๐ฎ(๐ฑ,๐‘›ฮ”๐‘ก),๐‘๐‘›๐‘โ‰ก๐‘(๐ฑ,๐‘›ฮ”๐‘ก),๐‘›1โ‰ก๐‘1(๐ฑ,๐‘›ฮ”๐‘ก),๐‘๐‘›2โ‰ก๐‘2๐œ™(๐ฑ,๐‘›ฮ”๐‘ก),๐‘›โ‰ก๐œ™(๐ฑ,๐‘›ฮ”๐‘ก),(2.24) the numerical strategy can be summarized as follows.

Given the initial values ๐ฎ0,๐‘0,๐‘01,๐‘02, and ๐œ™0, find the solution at time ๐‘›+1 solving (2.25)โ€“(2.29), for ๐‘›=0,1,2,โ€ฆ,๐œŒ๐ฎ๐‘›+1โˆ’๐ฎ๐‘›๎€ท๐ฎฮ”๐‘ก+๐œŒ๐‘›+1๎€ธ๐ฎโ‹…โˆ‡๐‘›๎€บ=โˆ‡โ‹…โˆ’๐‘๐‘›+1๎€ท๐ˆ+๐œ‚โˆ‡๐ฎ๐‘›+(โˆ‡๐ฎ๐‘›)๐‘ก๎€ท๐‘๎€ธ๎€ปโˆ’๐น๐‘›1โˆ’๐‘๐‘›2๎€ธโˆ‡๐œ™๐‘›,(2.25)โˆ‡โ‹…๐ฎ๐‘›+1๐‘=0,(2.26)1๐‘›+1โˆ’๐‘๐‘›1ฮ”๐‘ก+๐ฎ๐‘›+1โ‹…โˆ‡๐‘1๐‘›+1๎€ท=โˆ‡โ‹…๐ทโˆ‡๐‘๐‘›1+๐‘ค๐น๐‘๐‘›1โˆ‡๐œ™๐‘›๎€ธ๐‘,(2.27)2๐‘›+1โˆ’๐‘๐‘›2ฮ”๐‘ก+๐ฎ๐‘›+1โ‹…โˆ‡๐‘2๐‘›+1๎€ท=โˆ‡โ‹…๐ทโˆ‡๐‘๐‘›2โˆ’๐‘ค๐น๐‘๐‘›2โˆ‡๐œ™๐‘›๎€ธ,(2.28)โˆ’๐œ€0๐œ€๐‘Ÿโˆ‡2๐œ™๐‘›+1๎€ท๐‘=๐น1๐‘›+1โˆ’๐‘2๐‘›+1๎€ธ.(2.29)

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 though ๐‘ 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 ๐‘ก=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 ๐‘ฆ 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 ๐‘ฅ and ๐‘ฆ coordinates. Results are next.

3.1.1. Case A.1: ๐œ™(0)=๐œ™0๐‘ฆ

Figure 5 shows the three-dimensional image of the velocity field u. Surface and height correspond to |๐ฎ|. 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.

In Figure 6 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.

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 ๐‘ฅ=5โ€‰mm constant and varying ๐‘ฆ, 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.

Electric charge density ๐œŒ๐‘ž 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 ๐œŒ๐‘ž qualitatively agree with results reported in the literature [8, 16, 23]. Units are C/m3.

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)=๐œ™0๐‘ฅ๐‘ฆ

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)=๐œ™0๐‘ฅ๐‘ฆ. 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 for u 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. Figure 11 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 ๐‘ฅ=0. The cross-section for ๐œŽ at ๐‘ฅ=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.

Similar behavior appears in the plot of electric charge density, as shown in Figure 12. There is negligible density in the electrodesโ€™ half near the axis ๐‘ฅ=0. The cross-section of ๐œŒ๐‘ž 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 ๐œŒ๐‘ž in the basis scenario (c.f. Figure 10). Again, electroneutrality is achieved in the bulk.

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 ๐ด๐‘… stated by the charge transfer control theory [3] as shown in ๐ด๐‘…=๐‘–0๎€บ๐‘’(๐›ผ๐‘Ž๐น/๐‘…๐‘‡)(๐œ™โˆ’๐œ™0)โˆ’๐‘’โˆ’(๐›ผ๐‘๐น/๐‘…๐‘‡)(๐œ™โˆ’๐œ™0)๎€ป,(3.1) where ๐‘–0 is the anodic exchange current density, measured in A/m2, ๐›ผ๐‘Ž and ๐›ผ๐‘ are the charge transfer coefficients, anodic and cathodic, respectively (and depend on the specific couple fuel-oxidant used in the experiment), ๐น is the Faradayโ€™s constant, ๐‘… is the universal gas constant, and ๐‘‡ is the absolute temperature.

Assuming ๐›ผ๐‘Ž=๐›ผ๐‘=๐›ผ that is, one step reaction in which only one electron is transferred, it becomes into the hyperbolic reaction:HAR=2๐‘–0sinh๎‚ƒ๐›ผ๐น๎€ท๐‘…๐‘‡๐œ™โˆ’๐œ™0๎€ธ๎‚„.(3.2)

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 ฮฉ: ๐œ•๐‘1๎€ท๐œ•๐‘ก+โˆ‡โ‹…โˆ’๐ทโˆ‡๐‘1โˆ’๐‘ค๐น๐‘1๎€ธโˆ‡๐œ™=HARโˆ’๐ฎโ‹…โˆ‡๐‘1.(2.12')

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: ๎€ทโˆ’๐งโ‹…โˆ’๐ทโˆ‡๐‘1โˆ’๐‘ค๐น๐‘1โˆ‡๐œ™+๐‘1๐ฎ๎€ธ=HAR.(2.19')
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 ๐‘1 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: HAR๐‘›โ‰กHAR(๐ฑ,๐‘›ฮ”๐‘ก)=2๐‘–0๎‚ƒsinh๐›ผ๐น๎€ท๐œ™๐‘…๐‘‡๐‘›โˆ’๐œ™0๎€ธ๎‚„,๐‘1๐‘›+1โˆ’๐‘๐‘›1ฮ”๐‘ก+๐ฎ๐‘›+1โ‹…โˆ‡๐‘1๐‘›+1๎€ท=โˆ‡โ‹…๐ทโˆ‡๐‘๐‘›1+๐‘ค๐น๐‘๐‘›1โˆ‡๐œ™๐‘›๎€ธ+HAR๐‘›.(2.29')
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 for u and ๐œ™ are practically the same as those shown in Figures 5 and 6, respectively. There is no worth presenting such graphs here.
Figure 14 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 ๐œŒ๐‘ž 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 density ๐‘–0 were introduced to study their effect on ๐œŽ and ๐œŒ๐‘ž. Letting ๐‘–0 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 ๐œŒ๐‘ž, respectively, when doing this sensitivity analysis on the anodic exchange current density ๐‘–0. 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.

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 Table 2. 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 ๐‘1 and ๐‘2, and the electric potential ๐œ™ as well. Linear elements were used to approximate pressure ๐‘.

Time domain interval for simulation was set to [0,1] seconds taking a time discretization step of ฮ”๐‘ก=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 field u 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 field u 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 ๐‘1 and ๐‘2 on viscosity in the form stated by (4.1) is usually observed in solutions of larger molecules [25].๐œ‚=๐œ‚0๎‚ƒ๎€ท๐‘1+๐›ฝ1+๐‘2๎€ธ2๎‚„,(4.1) 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 ๐œŒ๐‘ž 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 for ๐‘1 and ๐‘2, as follows: ๐‘0,๐‘0๐‘ฅ, and ๐‘0๐‘ฅ๐‘ฆ. Value for ๐‘0 is shown in Table 1. Again, it was observed the same solution for u, ๐œ™, ๐œŽ, and ๐œŒ๐‘ž 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 Figure 8: 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 ๐œ™0๐‘ฅ or ๐œ™0๐‘ฅ๐‘ฆ 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 Section 3 allows speculating the following. There is a linear dependence between maxima values of electric conductivity ๐œŽ and charge density ๐œŒ๐‘ž on the exchange current density ๐‘–0, 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 velocity u and electric potential ๐œ™, independence of u, ๐œ™, ๐œŽ, and ๐œŒ๐‘ž with respect to variable viscosity ๐œ‚ and initial conditions for charge molar concentrations ๐‘1 and ๐‘2, general shape for cross-sections of ๐œŽ and ๐œŒ๐‘ž, awareness in the establishment of consistent initial conditions for ๐œ™ with the Arrhenius-like reactions on anode, and linear dependence of ๐œŽ and ๐œŒ๐‘ž 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.

Acknowledgments

This work was partially sponsored by the University of California at Berkeley and was completed during a research visit of the author. He gratefully thanks the Berkeley Sensor and Actuator Center and the Summer Sessions Department.