% % Unit fuel cell code, using the formulation developed in the UPC % workshop lectures % % This code is exercise #4 in those lectures % %------------------------ % While in code development, just in case clear param sol % parameters % physical parameters param.R = 8.314; param.F = 96485; % Electrochemistry fitted constants param.alpha = 1; param.E = 1.28; param.Cref = 40.9; param.i0 = 1e-7; param.delta = 8; param.Rm = 0.1; % Operating conditions param.T = 350; param.iave = 1; param.s = 1.1; param.Ctot = 100; % Numerical parameters N= 100; param.N = N; param.h = 1/N; Ntheta = 5; %------------------------- % Find theta=0 solution % First oxygen channel concentration Ctot = param.Ctot; C0 = 0.21 * Ctot; param.C0 = C0; sol.C = zeros(N+1,1) + C0; % Now cell voltage F = param.F; delta = param.delta; E = param.E; R = param.R; T = param.T; alpha = param.alpha; i0 = param.i0; Cref = param.Cref; Rm = param.Rm; iave = param.iave; I = iave; U = E - R*T/alpha/F*(log(I/i0)- log((C0-delta*I)/Cref))-Rm*I; sol.U = U; % Now constant nitrogen flux s = param.s; Qn = 0.79/0.21*s*iave/4/F; param.Qn = Qn; %-------------------------- % Now do continuation to find the solution % A more sophisticated code would use variable steps for the % continuation process vtheta = linspace(0,1,Ntheta); for icont = 2:Ntheta theta = vtheta(icont); final =0; if icont == Ntheta final = 1; end [fail sol] = solve_unit(sol, param, theta, final); if fail == 1 disp(sprintf('Solution Failure: trying increasing Ntheta')) disp(sprintf('Paused - hit control c')) pause end end h = param.h; xc = ((1:(N+1))-1)*h; xi = ((1:N)-0.5)*h; figure(1) plot(xc, sol.C, 'k-') axis([0 1 0 C0+1]) xlabel('x') ylabel('Channel O2 concentration') figure(2) plot(xi, sol.i, 'k-') imax = C0/delta; axis([0 1 0 imax]) xlabel('x') ylabel('Current density') disp(sprintf('Unit cell voltage: %d',sol.U));