# Problem with Numerical DAE Implementation of Nernst-Planck Equations

I write to you all today in order to clarify if I implemented my model correctly into MATLAB.

I will outline the following: analytical governing equations, my discretization for the equations and boundary conditions, and outline the implementation into ode15s. Ultimately, my problem with implementation is that steady-state is not achieved as I reach an infinite accumulation with and without a volumetric rate term. I want to make sure that my implementation is correct before reevaluating the governing equations.

System description - arbitrary cation and anion react with gaseous species in porous electrode. Cell potential is specified. BV kinetics is utilized but not included here. Concentrations are mol/m^3. Fluxes are mol/m^2/s. Volumetric charge is A/m^3.

The gas balances work and are not included.

Governing equations: Discrete form of the equations: How to get boundary condition for anionic species: Numerical solving procedure - method of lines with ode15s. Finite differences with second order central and first order central. Electroneutrality equation applies but is not used in DAE solver. Total volumetric charge conservation is used as algebraic constraint. To deal with the imaginary points , I solve for the cation, anion, and potential value via the boundary values.

5 equations (3 listed + 2 not-listed) / node. n refers to number of nodes. Length of cell is on the range of a 100 microns.

Three big questions:

1. The boundary condition for the anionic species (i.e. a) is derived utilizing the electroneutrality as applied to the cationic boundary condition. Is this sufficient to specify the boundary value problem?

2. I have a 5xNode# by 5xNode# mass matrix with 1 for equations 1-4 on the diagonals and 0 for equation 5. Is this mass matrix constructed correctly?

3. The arbitrary reference potential is set as 0 V. I specify in my mass matrix a value of 1 for the charge balance as I set the potential. I set my yt for the charge balance at node n = 0 as the initial value is set as the boundary condition. Is this causing an instability problem or should the mass matrix be set to 0 for this equation?

Thanks for any help!