% % Stack fuel cell code, using the formulation developed in the UPC % workshop lectures % % This code is exercise #12 in those lectures % %------------------------ % While in code development, just in case clear param sol sols % 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 = 2.0; param.Ctot = 100; % Stack parameters M = 7; % There are 2M-1 cells in the stack param.M = M; param.sanom = 1.2; % anomalous stoich of centre cell #1 param.lambda = 100; % bipolar plate resistivity % Numerical parameters N= 100; param.N = N; param.h = 1/N; Ntheta = 5; %------------------------- % Find theta=0 solution for the background unit cell % 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; % Convenient place to here, although will only be used in the stack % computations later sanom = param.sanom; Qna = 0.79/0.21*sanom*iave/4/F; param.Qna = Qna; %-------------------------- % Now do continuation to find the base unit cell solution vtheta = linspace(0,1,Ntheta); disp(sprintf('Begin unit cell computation at base conditions')) 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('Unit Cell 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') fprintf('Unit cell voltage: %d \n \n',sol.U); %------------------------- % Now initialize the stack computation Cs = zeros((N+1),M); for m=1:M Cs(:,m) = sol.C; end Us = zeros(N,M)+ sol.U; sols.Cs = Cs; sols.Us = Us; % Use continuation to solve the stack problem disp(sprintf('Begin stack computation')) for icont = 2:Ntheta theta = vtheta(icont); final =0; if icont == Ntheta final = 1; end [fail sols] = solve_stack(sols, param, theta, final); if fail == 1 disp(sprintf('Stack Solution Failure: trying increasing Ntheta')) disp(sprintf('Paused - hit control c')) pause end end Us = sols.Us; is = sols.is; figure(3) plot(xi,Us(:,1),'k-',xi,Us(:,2),'k-',xi,Us(:,3),'k-',xi,Us(:,M),'k-') figure(4) plot(xi,is(:,1),'k-',xi,is(:,2),'k-',xi,is(:,3),'k-',xi,is(:,M),'k-')