Discussion Forum

Hi all:

I am trying to model an electrochemical surface (i.e. disk electrode) reaction in order to predict i-V curves for a deposition process and I do not have the Chemical Engineering Module. So, I am trying to formulate the problem using coefficient form PDE modes (to represent the mass-transport/potential field Nernst-Planck Equations for each species) and the Conductive-dc mode for the potential.

The Nernst-Planck equations involve the gradient of the potential field. The way I have formulated the Conductive-dc mode for the potential, it is a scalar multiple of the mass transport equations (i.e. [mols/area*time]*Faraday constant = current density). Therefore; for n species, I have n+1 unknowns when counting the potential field. This is inherently un-solvable.

However, I have an electroneutrality constraint that says the sum of the species at all points times each species' respective valence equals zero. This problem now is do-able, but I cannot implement this constraint for the life of me!!

I have tried using an Integration Coupling Variable to make the aforementioned sum of species*valence terms globally available. I then tried minimizing this variable to zero by altering one of the non-electrochemically active species as a variable using the Global Equations feature. When I include PDE modes for all n species, I have an error involving a duplicate variable name (i.e. the non-active species concentration is the dependent variable of the PDE mode and the variable in the Global Equations feature).

When I try the same using PDE modes for n-1 species (only electrochemically-active species) and solving for the last species with electroneutrality, I get an error stating that the gradient of the implied species (that calculated with by the constraint) cannot be calculated (i.e. c_SO4zz (axial gradient of sulfate concentration) is undefined). These gradients are terms in the potential field governing equation.

It seems that this is a fully-determined problem and that I should be able to solve it, but I cannot figure out the ways of COMSOL to do so.

I've attached my model in case anyone with some insight would be able to help me. It is an axisymmetric geometry with an electrode disk core in the center of an insulating shell. The boundary conditions are fixed concentrations at the top of the simulation rectangle. At the electrode, the Butler-Volmer equation modified to express in units of mass flux based on stoichiometry of the reaction. The potential BC's are set as fixed potential far from the surface and set as current flux (determined by the Butler-Volmer equation) at the disk.

Would really appreciate some feedback!