Modeling of Interdigitated Electrodes and Supercapacitors with Porous Interdigitated Electrodes

As industrial and commercial demand for thin, small-footprint energy storage devices increases (e.g., for wearable electronics, wireless sensor networks, etc.), it is essential to develop theoretical descriptions for those devices in order to make effective design choices. Of particular interest are porous interdigitated electrodes, which can function as supercapacitors but which are nontrivial to analyze due to their unique geometry. In this paper, we develop a purely mathematical model to determine the dependence of a cell’s capacitance and/or resistance on electrode height, width, and spacing, and electrolyte layer height. We then extend this model to incorporate effects from the porous nature of the electrodes, and show that in many cases of interest, the system can be described using a simple equivalent circuit. © The Author(s) 2017

The advancement of miniaturized, wearable electronics has increased the demand for thin, flexible energy storage. Two-dimensional or in-plane electrode geometries can be fabricated in a single process and are more robust against bending when compared to stacked electrode geometries. 1,2 Supercapacitors are especially relevant to modern electronics, as they can readily meet high power demands associated with wireless communication. Here we develop an analytical model of interdigitated supercapacitor electrodes to understand the effect of device design on the performance of the planar supercapacitor.
Modeling supercapacitors with interdigitated geometries requires consideration of the nontrivial geometry itself as well as how it couples to the behavior of porous electrodes. While theories for these two components have been developed independently, we are unaware of a study which incorporates both. Wei first used the technique of conformal mapping to derive in closed form the capacitive properties of a periodic array of electrodes. 3 A subsequent extension of the theory by Igreja and Dias accounts for a finite dielectric layer height and substrate effects, 4 but not for a finite electrode height, which we explore in this paper.
The characteristics of the cell with interdigitated electrodes, described by El-Kady and Kaner, depend on the electrodes being porous. 1 There has been much analysis of porous electrodes since the 1960s, for example the works of Euler and Nonnenmacher, or Newman and Tobias. 5,6 A recent example is the modeling of porous electrodes of lithium batteries by Thomas et al. 7 In general the equations for such systems become difficult to solve in greater than one dimension, where the electrical potential gradient is nonzero in multiple directions. Here, we make use of the numerical solver COMSOL, and combine results from finite element analysis (FEA) with those from conformal mapping to provide insights on dependence of an interdigitated porous electrode cell's electrical behavior on its geometric characteristics.

Theory and Results
Interdigitated geometry.-For an electrical cell comprised of two conducting terminals separated by some dielectric or conductive medium of permittivity or conductivity σ, it is possible to find the capacitance (C) or resistance (R) of the cell, respectively, via the following relations: Here, V is the voltage difference between the two terminals, E is the electric field as a function of position, and S is the surface of one of the terminals. For an ideal parallel plate capacitor (with no fringe fields) of plate area A and plate separation d, these integrals are easily evaluated when and σ are constants in space: Therefore, if the (possibly nontrivial) cell geometry of interest can be transformed into a geometry resembling a parallel plate configuration in such a way that preserves the orthogonality of electric field lines and equipotential lines (as demanded by the solution to Laplace's equation which specifies the voltage at every point in the cell), then that transformation will easily provide the desired values of resistance and capacitance. Conformal mappings, or transformations of the 2D plane which preserve right angles, can be used to this end. Specifically, the Schwarz-Christoffel equation maps the upper half of the plane to any arbitrary polygon with specified interior angles. 8 Olthius et al. used this method to derive the dependence of capacitance on the ratio of electrode width to electrode spacing for interdigitated, twodimensional electrodes of specified number and length, assuming an infinitely thick medium above the electrodes. 9 Here, we investigate in addition the effects of finite electrode and electrolyte layer heights. Following the assumptions made in that paper, we take the electrodes to be of a sufficient length, L, that fringe effects at the top and bottom buses in Fig. 1a due to the gap δ can be ignored. We also assume the number of electrode couples, N , is large enough that fringe effects at either end of the device may also be neglected. This establishes a single symmetry unit as the basis of analysis, as shown in Fig. 1b. Comments on the magnitude of error introduced by these assumptions are given at the end of the section.
If we denote R unit and C unit as the resistance or capacitance (respectively) per unit length derived from the single cell, then the resistance or capacitance of the overall device will be: Mappings from the upper half-plane to both the desired geometry and the equivalent parallel plate geometry are shown qualitatively in  If we specify the coordinates of the original upper half-plane as z, those of the desired cell geometry as z , and those of the equivalent parallel plate geometry as z , then the mathematical transformations are given as: The integral in Eq. 8 can be transformed into an elliptic integral via a change in variables. 9 Unfortunately the integral in Eq. 7 does not simplify to a well-known form, but it can similarly be transformed into an expression which allows integration between points of discontinuity (See Appendix). Since the integrals in Eqs. 7 and 8 are not analytically solvable, it is not possible to start with a desired cell geometry and from there derive its capacitance and resistance parameters. Instead, both integrals were evaluated numerically using the 'integral' function in MATLAB for a wide range of combinations of values for (A, B, C, D), giving corresponding values (A , B , C , D ) and (A , C ) for the two different mappings. This produced a "database" of geometries which could be searched for a given cell geometry of interest, with corresponding R unit , C unit given by: It should be noted that although the mapping given by Eq. 7 produces a unique result ( A , B , C , D ) for any single (A, B, C, D), that result may be scaled by any constant factor without affecting its C unit .
The effects of each parameter in the four-dimensional parametric space were visualized by selecting several different baseline geometries in the z space, as shown in Fig. 3. For each of the geometries, the effect of varying electrode spacing, electrode height, electrode width, and electrolyte layer height (or A , B , C , and D , respectively) in isolation is shown in the graphs in Fig. 4 (The quantities A and B correspond to S 2 and W 2 given in Fig. 1). Each color represents the corresponding baseline geometry in Fig. 3. Normalized capacitance is given by C A , and normalized dimensions are given as a fraction of the sum of all dimensions; for example, the normalized A' dimension, A norm , is  Table I. As would be expected, capacitance increases with increasing electrode height (B') and width (C'), although the latter dependence is weak, and decreases with electrode spacing (A'). The topmost geometry in Fig. 3 is considered in greater detail for the remainder of the paper.
Results were verified using FEA in COMSOL, to excellent agreement. A comparison of electric field maps between the two methods is shown in Fig. 5. For the particular geometry shown, the unit capacitance of the cell with a vacuum dielectric was found to be 3.3619 pF/m and 3.3585 pF/m via the conformal mapping and FEA methods, respectively.
One advantage of conformal mapping over numerical analysis is the degree of precision to which cell capacitance/resistance values can be calculated. As shown in Fig. 6, decreasing the mesh size of the numerical simulation results in calculated capacitance values which converge to that found via numeric integration in MATLAB. The numerically integrated values were constant to full available precision for a wide range of specified tolerances, which means the functions of interest are smooth enough to be calculated with minimal error. Furthermore, the time required to perform the transformation integrals for a single cell geometry on MATLAB is on the order of 50 ms, whereas the equivalent calculation takes roughly ten times that on COMSOL (both performed on a PC with an Intel(R) Core(TM) i7-6800K 3.40 GHz CPU and 32.0 GB memory). Therefore, in a situation where it is desirable to systematically analyze the general dependence of a cell's electrical parameters on its dimensions (for example as shown in Fig. 4), this method is far more efficient.   The mapping given by Eq. 6 assumes the electric field to be completely confined to the electrolyte layer, which becomes true in the limit where electrolyte / air → ∞ A simulation was carried out in COMSOL to assess the magnitude of the error resulting from this approximation, and it was found that the true capacitance (accounting for imperfect confinement of the electric field to the electrolyte layer) decreased by 14% with a electrolyte / air ratio of 10 for the geometry shown in Fig. 5. However, the calculation for cell resistance remained quite accurate when the surrounding medium was nonconductive (as would be the case in many applications).
Additional sources of error in the model described above comes from the neglect of fringe capacitance. This occurs at the two electrodes on the far ends of the whole device, as well as the end of each leg (shown by the gap δ in Fig. 1a). Based on the weak dependence of normalized capacitance on C norm observed in Fig. 4, we expected the  former effect to be insignificant. This was confirmed by simulating a unit cell in COMSOL with one electrode region being twice the length of the other (as would be the case at device ends). The change in capacitance was less than 0.1%. This result suggests that Eqs. 5 and 6 might be more accurately expressed by replacing N with (N −1), with the resulting fractional error depending on N . It also suggests that the fractional error due to the approximation of infinitely long electrodes can be reasonably estimated in the case that δ ≈ S as ≈ W W +L . For the devices in Ref. 1 we estimate that this will be roughly on the order of 10%.

Porous electrodes.-
In what follows we make the usual assumption that the matrix phase of the electrode and the electrolyte therein can be treated as two superimposed continua with macroscopically uniform properties. 6 We further assume that conduction in the matrix obeys Ohm's law with some effective conductivity. Finally, we assume that conduction in the electrolyte will also obey Ohm's law. This last assumption is made because we are primarily interested in electrolytes that are ionic liquids consisting of one cation, one anion and no solvent; consequently, electroneutrality mandates there can be no concentration gradients, outside double layers, and therefore conduction in the electrolyte is only by migration. The assumption is also valid for electrolytes, with a solvent, that are well supported electrolytes. The assumption would be an approximation in other cases.
We are considering a simple, 2-dimensional cross section of an interdigitated electrode as shown below in Fig. 7. We continued to study the geometry shown in Fig. 5 since it approximates that of a real device of interest, for example as fabricated in Ref. 1.
Regions I and III are treated as a two-phase system consisting of an ionic liquid electrolyte and a "perfectly conducting" matrix of porous material. Both are assumed to be continuous throughout the regions so that current may flow at any point either in the x or y directions in the electrolyte phase, or into the electrode phase. Region II is the bulk electrolyte. In that region we assume no accumulation of charge, and therefore the potential obeys Laplace's equation: For the electrodes, we consider the matrix/electrolyte interface to be purely capacitive (no faradaic or leakage current), and therefore the transfer of current from electrolyte to matrix at any point in the region is proportional to the time derivative in voltage difference between those two phases. Specifically, the current density in the electrolyte is: And therefore, Where C spec is the capacitance per unit volume of the matrix material and σ is the effective conductivity of the electrolyte within the matrix. We make the simplifying assumption that the matrix material is highly conductive, such that its potential is independent of position. Furthermore we suppose that the potential applied to the elctrodes  remains constant over time (for t > 0), and thus Eq. 11 simplifies to a PDE in one variable: The constant K is now equal to Cspec σ . For porous carbon-based materials, values cited in the literature for specific capacitance vary from 35 F/cm 3 to 112 F/cm 3 but are generally on the order of 100 F/cm 3 . 10 The conductivities of the ionic liquid electrolytes of interest in our investigation also vary, but are generally on the order of 1 S/m for salts with imidazolium-based cations, which are commonly used in devices of interest for this study. 11 These values were used to obtain the results below.
Spatial boundary conditions specify no component of the electric field perpendicular to outer boundaries, which are either electrically insulated (top and bottom) or symmetry planes (sides): The initial voltage distribution at t = 0 was obtained by setting e to 0 V and 1 V in the two electrodes, and solving for the stationary voltage distribution in the bulk electrolyte. To examine this, consider the voltage of the matrix to be switched on instantaneously so that the expression ∂( e − m ) ∂t equals δ(t = 0). At any point in the electrode regions, the charge can either move in the cross-sectional plane (equivalent to an infinite number of differential resistors in parallel), or to the boundary with the matrix phase (equivalent to a differential capacitor in parallel with the resistors). Because the "frequency" of the delta function is infinite, the capacitor will "look" to the charge like a short, and all current in the region will go to charging the capacitor in that instant. Therefore there will be no voltage drop in the crosssectional plane (since the capacitor is initially uncharged), and we can set e = m at t = 0 within the electrodes. The initial distribution is shown in Fig. 8.
The electric field strength (directly proportional to current density norm by Ohm's Law) is shown as a function of time over the discharging of the capacitor in Fig. 9. It is clear that most of the current is concentrated at the vertical edges of the electrodes, which corroborates the steep increase in capacitance versus electrode height observed in Fig. 4.

Equivalent circuit model for porous interdigitated electrodes.-
Of interest in this study is whether, and to what extent, a simplified equivalent circuit model can describe the behavior of the porous interdigitated electrode system. In order to investigate this, the parameters for electrolyte conductivity and electrode specific capacitance were varied over a range encompassing the baseline values of 1 S/m and 100 F/cm 3 (respectively), and exponential curves were fitted to the FEA output for current through the cell versus time. For the equivalent circuit shown in Fig. 10, which includes two capacitive and resistive elements (roughly representing the porous electrode matrices and the electrolyte in the electrode, respectively) in series with a third resistive element (roughly representing the resistance of the bulk electrolyte in the cell), the time constant would be given as: because the capacitors add reciprocally in series. As a first estimate, C was taken as the specific capacitance of the electrode material times the area of the electrode region. R 2 was found using the conformal mapping methods described in Interdigitated geometry section (for non-porous electrodes at a fixed potential difference of 2 V). This left R 1 as an unknown parameter to be determined by the simulation results. Two cases were examined: (1) In which the conductivity of the electrolyte was taken to be equal in both the bulk region (II in Fig.  7) and in the wetted electrode (I and III in Fig. 7); and (2) In which the resistivity of the electrolyte in the wetted electrode was specified as two orders of magnitude greater than that in the bulk region, as found in recent experiments. 12 Fig. 11 shows results from the first case, which confirm simple RC-circuit behavior within the range of electrolyte resistivities and electrode specific conductances examined, viz. linear dependence of the decay time constant on those parameters. Regression lines were constrained to pass through the origin in order to give a semiqualitative description of how well the simulation data conformed to Eq. 17. The values of R 1 and R 2 , calculated as described above, are both found to equal 3.3 · m, which indicates that resistance of the bulk electrolyte and electrolyte within the electrode contribute in equal measure to the cell's equivalent series resistance. Fig. 12 shows results from the second case. Although the decay time constant still varies nearly linearly with electrolyte resistivity and electrode specific conductance, the coefficient of determination for the linear regression is notably lower, which suggests deviation from the simplified equivalent circuit model. However, the calculated value of R 1 is 5.6 · m, indicating that even the drastic increase in resistance of the electrolyte in the electrode region does not increase the overall estimated series resistance at a comparable scale.
The difference observed between the two cases may be due in part to non-uniform charging/discharging throughout the electrode when the conductivity of the electrolyte in that region is severely limited. Fig. 13 shows surface maps of charge/discharge current    density throughout the cell generated by FEA simulations. These maps confirm that this phenomenon indeed occurs, even for a relatively small electrode height.
We lastly investigated the effects of varying only the electrode height (given by B') on the charge/discharge time constant of the cell. This dimension was swept along a range of values from 5 to 150 μm, from an original value of 10 μm. This simulation was again performed for the same geometry, in the two cases outlined above for conductivity of the electrolyte in the electrode. Results are shown in Figures 14 and 15. Close to the original value of B', the time constant varies approximately linearly with B', but at larger electrode heights the relationship appears asymptotic and even shows a negative correlation for the case of decreased conductivity of the electrolyte in the electrode. This behavior can be rationalized by the fact that the height of the electrode in the geometry of interest is quite small relative to the other dimensions of the cell, and therefore its normalized value B norm , given by B A +B +C +D , does not vary significantly (in absolute terms) with B'. For example, for the linear regime displayed in the inset of Fig. 13, B norm varies from roughly 1% to 3%. Upon examination of Fig. 4, it is clear that the cell constant does not vary significantly in this range. The main effect on the time constant, then, comes from the change in cross-sectional area of the porous electrode, which affects C and R 1 only. However, as B norm increases past this linear regime, the effects of decreased R 2 in addition to increased R 1 and C result in a leveling-off of the time constant dependence on electrode height. We hypothesize that the negative slope beyond the linear regime in Fig. 15 is due to a restriction of capacitative behavior to the edges of the electrode (Fig. 13), which means that C does not scale as rapidly as the area of the electrode with increasing electrode height.

Conclusions
The Schwarz-Christoffel conformal mapping can be used to analyze the relationship between a cell's capacitance/resistance and its geometric parameters (electrode height, width, and spacing, and electrolyte layer height) for an interdigitated electrode configuration. Results from this mathematical technique are confirmed by -and even show signs of improvement in speed and accuracy over -FEA. Therefore, given specific practical constraints (for example on cell volume, amount of electrode material to be used, etc.), one can optimize the cell's electrical characteristics.
By combining the results of the above model with concepts from porous electrode theory, we shwed that in many circumstances, a supercapacitor made of porous interdigitated electrodes can be effectively modeled using a simple equivalent circuit, even in the realistic scenario of greatly reduced electrolyte conductivity within the porous matrix of the electrode. Furthermore, FEA simulations in which electrode height is varied in isolation provide results which agree with qualitative predictions from conformal mapping, and also suggest a point of diminishing return with regard to maximizing the cell's effective RC time constant. Overall, our models provide quantitative metrics which can be used to inform the design and fabrication of porous interdigitated electrode supercapacitors in order to achieve optimal performance.

Appendix: Integral Transformations
Although the integral in Eq. 7 is divergent at the desired boundaries of integration, it diverges at those points more rapidly than z −1 , which means the integrals converge to a finite value. The following transformations were performed in order to calculate the definite integrals numerically: The mapping between z and z" coordinates (Eq. 8) is given in Ref. 9.