?? fm_sofc.m
字號:
function fm_sofc(flag)% FM_SOFC define Solyd Oxyde Fuel Cells%% FM_SOFC(FLAG)% FLAG = 0 -> initialization% FLAG = 1 -> algebraic equations% FLAG = 2 -> algebraic Jacobians% FLAG = 3 -> differential equations% FLAG = 4 -> state Jacobians% FLAG = 5 -> non-windup limiters%%Author: Federico Milano%Date: 25-Nov-2003%Version: 2.0.0%%E-mail: fmilano@thunderbox.uwaterloo.ca%Web-site: http://thunderbox.uwaterloo.ca/~fmilano%% Copyright (C) 2002-2005 Federico Milano%% This toolbox is free software; you can redistribute it and/or modify% it under the terms of the GNU General Public License as published by% the Free Software Foundation; either version 2.0 of the License, or% (at your option) any later version.%% This toolbox is distributed in the hope that it will be useful, but% WITHOUT ANY WARRANTY; without even the implied warranty of% MERCHANDABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU% General Public License for more details.%% You should have received a copy of the GNU General Public License% along with this toolbox; if not, write to the Free Software% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307,% USA.global Sofc DAE Bus PV SWIk = DAE.x(Sofc.Ik);Vk = DAE.x(Sofc.Vk);Vb = DAE.V(Sofc.bus);pH2 = DAE.x(Sofc.pH2);pH20 = DAE.x(Sofc.pH20);pO2 = DAE.x(Sofc.pO2);qH2 = DAE.x(Sofc.qH2);m = DAE.x(Sofc.m);Sn = Sofc.con(:,2);Vn = Sofc.con(:,3);Te = Sofc.con(:,4);TH2 = Sofc.con(:,5);KH2 = Sofc.con(:,6);Kr = Sofc.con(:,7);TH20 = Sofc.con(:,8);KH20 = Sofc.con(:,9);TO2 = Sofc.con(:,10);KO2 = Sofc.con(:,11);rHO = Sofc.con(:,12);Tf = Sofc.con(:,13);Uopt = Sofc.con(:,14);r = Sofc.con(:,17);N0 = Sofc.con(:,18);E0 = Sofc.con(:,19);RTon2F = 0.5*8.314*Sofc.con(:,20)/96487;Pref = Sofc.con(:,21);Vref = Sofc.con(:,22);xt = Sofc.con(:,26);Km = Sofc.con(:,27);Tm = Sofc.con(:,28);m_max = Sofc.con(:,29);m_min = Sofc.con(:,30);switch flag case 0 % initialization % Check time constants idx = find(Te == 0); if idx fcwarn(idx, ['Time constant Te cannot be zero. ', ... 'Te = 0.001 s will be used.']) Sofc.con(idx,5) = 0.001; end idx = find(TH2 == 0); if idx fcwarn(idx, ['Time constant TH2 cannot be zero. ', ... 'TH2 = 0.001 s will be used.']) Sofc.con(idx,7) = 0.001; end idx = find(TH20 == 0); if idx fcwarn(idx, ['Time constant TH20 cannot be zero. ', ... 'TH20 = 0.001 s will be used.']) Sofc.con(idx,10) = 0.001; end idx = find(TO2 == 0); if idx fcwarn(idx, ['Time constant TO2 cannot be zero. ', ... 'TO2 = 0.001 s will be used.']) Sofc.con(idx,12) = 0.001; end idx = find(Tf == 0); if idx fcwarn(idx, ['Time constant Tf cannot be zero. ', ... 'Tf = 0.001 s will be used.']) Sofc.con(idx,15) = 0.001; end idx = find(Tm == 0); if idx fcwarn(idx, ['Time constant Tm cannot be zero. ', ... 'Tm = 0.001 s will be used.']) Sofc.con(idx,28) = 0.001; end % Find and delete PV generators check = 1; for j = 1:Sofc.n idx_pv = []; idx_sw = []; if PV.n idx_pv = find(PV.bus == Sofc.bus(j)); end if SW.n idx_sw = find(SW.bus == Sofc.bus(j)); end if isempty(idx_pv) & isempty(idx_sw) fm_disp(['Error: No generator found at bus #', ... int2str(Sofc.bus(j))]) check = 0; elseif ~isempty(idx_pv) PV.con(idx_pv,:) = []; PV.bus(idx_pv) = []; PV.n = PV.n - 1; else fm_disp(['Warning: fuel cells cannot ', ... 'substitute slack buses.']) fm_disp(['Fuel cells will be initialized but ', ... 'eigenvalue analyses and time domain']) fm_disp(['simulations will produce errors ', ... '(singularity in the Jacobian matrix).']) end end if ~check fm_disp('Fuel cells cannot be properly initialized.') return end % Get generator powers Pfc = Bus.Pg(Sofc.bus); Qfc = Bus.Qg(Sofc.bus); % Define reference active power Sofc.con(:,21) = Pfc; % Fuel Cell stae variable initialization DAE.x(Sofc.Ik) = Sn.*Pfc./Vb./Vn; Ik = DAE.x(Sofc.Ik); DAE.x(Sofc.qH2) = 2*Kr.*Ik./Uopt; qH2 = DAE.x(Sofc.qH2); DAE.x(Sofc.pO2) = (qH2./rHO-Kr.*Ik)./KO2; pO2 = DAE.x(Sofc.pO2); DAE.x(Sofc.pH2) = (qH2 - 2*Kr.*Ik)./KH2; pH2 = DAE.x(Sofc.pH2); DAE.x(Sofc.pH20) = 2*Kr.*Ik./KH20; pH20 = DAE.x(Sofc.pH20); DAE.x(Sofc.Vk) = -r.*Ik./Vn + ... N0.*(E0+RTon2F.*log(pH2.*sqrt(pO2)./pH20))./Vn; Vk = DAE.x(Sofc.Vk); % Base power and voltage Sofc.con(:,23) = (Ik.*Vk.*Vn./Sn)./Pfc; Sofc.con(:,24) = Vb./Vk; % Adjust value of control type Sofc.con(:,25) = ~Sofc.con(:,25); idx = find(Sofc.con(:,25)); if idx, Sofc.con(idx,25) = Vk(idx); end % Initialize tap ratio DAE.x(Sofc.m) = sqrt(((xt./DAE.V(Sofc.bus)./Vk./Sofc.con(:,24)).^2) ... .*(Pfc.^2+(Qfc+(DAE.V(Sofc.bus).^2)./xt).^2)); % Define reference AC voltage Sofc.con(:,22) = DAE.V(Sofc.bus)+DAE.x(Sofc.m)./Km; fm_disp('Initialization of Solid Oxyde Fuel Cells completed.') case 1 % algebraic equations DAE.gp = DAE.gp - sparse(Sofc.bus,1,Ik.*Vk.*Vn./Sn.*Sofc.con(:,24),Bus.n,1); Vt = m.*Vk.*Sofc.con(:,24); Vs = DAE.V(Sofc.bus); Q = -Vs.*(Vs - Vt.*sqrt(1-(xt.*Ik./Vs.*Vn./Sn./m).^2))./xt; DAE.gq = DAE.gq - sparse(Sofc.bus,1,Q,Bus.n,1); case 2 Vt = m.*Vk.*Sofc.con(:,24); Vs = DAE.V(Sofc.bus); sq = sqrt(1-(xt.*Ik./Vs.*Vn./Sn./m).^2); dQdv = -2*Vs./xt+Vt./xt.*sq+0.5*Vs.*Vt./xt.*(2*((xt.*Ik.*Vn./Sn./m).^2)./(Vs.^3))./sq; DAE.J22 = DAE.J22 - sparse(Sofc.bus,Sofc.bus,dQdv,Bus.n,Bus.n); case 3 % differential equations Pref = Sofc.con(:,21).*Sofc.con(:,23); Umax = 0.5*Sofc.con(:,15).*qH2./Kr; Umin = 0.5*Sofc.con(:,16).*qH2./Kr; Input = Sn.*Pref./Vk./Vn; idx = find(Sofc.con(:,25)); if idx, Input(idx) = Sn(idx).*Pref(idx)./Sofc.con(idx,25)./Vn(idx); end idx = find(Input > Umax); if idx, DAE.f(Sofc.Ik(idx)) = (Umax(idx) - Ik(idx))./Te(idx); end idx = find(Input < Umin); if idx, DAE.f(Sofc.Ik(idx)) = (Umin(idx) - Ik(idx))./Te(idx); end idx = find(Input <= Umax & Input >= Umin); if idx, DAE.f(Sofc.Ik(idx)) = (Input(idx) - Ik(idx))./Te(idx); end DAE.f(Sofc.Vk) = (-Vk-r.*Ik./Vn+N0.*(E0+RTon2F.*log(pH2.*sqrt(pO2)./pH20))./Vn)./1e-3; DAE.f(Sofc.pH2) = ((qH2-2.*Kr.*Ik)./KH2-pH2)./TH2; DAE.f(Sofc.pH20) = (2.*Kr.*Ik./KH20-pH20)./TH20; DAE.f(Sofc.pO2) = ((qH2./rHO-Kr.*Ik)./KO2-pO2)./TO2; DAE.f(Sofc.qH2) = (2.*Kr.*Ik./Uopt-qH2)./Tf; DAE.f(Sofc.m) = (Km.*(Vref-DAE.V(Sofc.bus))-m)./Tm; idx = find(m >= m_max & DAE.f(Sofc.m) > 0); if idx, DAE.f(Sofc.m(idx)) = 0; end DAE.x(Sofc.m) = min(DAE.x(Sofc.m),m_max); idx = find(m <= m_min & DAE.f(Sofc.m) < 0); if idx, DAE.f(Sofc.m(idx)) = 0; end DAE.x(Sofc.m) = max(DAE.x(Sofc.m),m_min); case 4 % state variable Jacobians Pref = Sofc.con(:,21).*Sofc.con(:,23); % DAE.Fx Umax = 0.5*Sofc.con(:,15).*qH2./Kr; Umin = 0.5*Sofc.con(:,16).*qH2./Kr; Input = Sn.*Pref./Vk./Vn; idxc = find(Sofc.con(:,25)); if idxc Input(idxc) = Sn(idxc).*Pref(idxc)./Sofc.con(idxc,25)./Vn(idxc); end idx = find(Input > Umax); if idx DAE.Fx = DAE.Fx + sparse(Sofc.Ik(idx),Sofc.qH2(idx), ... 0.5*Sofc.con(idx,16)./Kr(idx),DAE.n,DAE.n); end idx = find(Input < Umin); if idx DAE.Fx = DAE.Fx + sparse(Sofc.Ik(idx),Sofc.qH2(idx), ... 0.5*Sofc.con(idx,16)./Kr(idx),DAE.n,DAE.n); end idx = find(Input <= Umax & Input >= Umin); idxc = find(~Sofc.con(idx,25)); if idxc, idx = idx(idxc); end if idx DAE.Fx = DAE.Fx + sparse(Sofc.Ik(idx),Sofc.Vk(idx), ... -Sn(idx).*Pref(idx)./Vk(idx).^2 ... ./Vn(idx)./Te(idx),DAE.n,DAE.n); end DAE.Fx = DAE.Fx + sparse(Sofc.Ik,Sofc.Ik,-1./Te,DAE.n,DAE.n); DAE.Fx = DAE.Fx + sparse(Sofc.Vk,Sofc.Vk,-1e3,DAE.n,DAE.n); DAE.Fx = DAE.Fx + sparse(Sofc.Vk,Sofc.Ik,-r./Vn./1e-3,DAE.n,DAE.n); DAE.Fx = DAE.Fx + sparse(Sofc.Vk,Sofc.pH2,N0.*RTon2F./pH2./Vn./1e-3,DAE.n,DAE.n); DAE.Fx = DAE.Fx + sparse(Sofc.Vk,Sofc.pH20,-N0.*RTon2F./pH20./Vn./1e-3,DAE.n,DAE.n); DAE.Fx = DAE.Fx + sparse(Sofc.Vk,Sofc.pO2,0.5*N0.*RTon2F./pO2./Vn./1e-3,DAE.n,DAE.n); DAE.Fx = DAE.Fx + sparse(Sofc.pH2,Sofc.Ik,(-2.*Kr./KH2)./TH2,DAE.n,DAE.n); DAE.Fx = DAE.Fx + sparse(Sofc.pH2,Sofc.pH2,(-1)./TH2,DAE.n,DAE.n); DAE.Fx = DAE.Fx + sparse(Sofc.pH2,Sofc.qH2,(1./KH2)./TH2,DAE.n,DAE.n); DAE.Fx = DAE.Fx + sparse(Sofc.pH20,Sofc.Ik,(2.*Kr./KH20)./TH20,DAE.n,DAE.n); DAE.Fx = DAE.Fx + sparse(Sofc.pH20,Sofc.pH20,(-1)./TH20,DAE.n,DAE.n); DAE.Fx = DAE.Fx + sparse(Sofc.pO2,Sofc.Ik,(-Kr./KO2)./TO2,DAE.n,DAE.n); DAE.Fx = DAE.Fx + sparse(Sofc.pO2,Sofc.pO2,(-1)./TO2,DAE.n,DAE.n); DAE.Fx = DAE.Fx + sparse(Sofc.pO2,Sofc.qH2,(1./rHO./KO2)./TO2,DAE.n,DAE.n); DAE.Fx = DAE.Fx + sparse(Sofc.qH2,Sofc.Ik,(2.*Kr./Uopt)./Tf,DAE.n,DAE.n); DAE.Fx = DAE.Fx + sparse(Sofc.qH2,Sofc.qH2,(-1)./Tf,DAE.n,DAE.n); Vt = m.*Vk.*Sofc.con(:,24); Vs = DAE.V(Sofc.bus); sq = sqrt(1-(xt.*Ik./Vs.*Vn./Sn./m).^2); dQdi = 0.5*Vs.*Vt./xt.*(-2*((xt.*Vn./Sn./m./Vs).^2).*Ik)./sq; dQdv = m.*Sofc.con(:,24).*Vs./xt.*sq; DAE.Gx = DAE.Gx - sparse(Sofc.bus+Bus.n,Sofc.Vk,dQdv,2*Bus.n,DAE.n); DAE.Gx = DAE.Gx - sparse(Sofc.bus+Bus.n,Sofc.Ik,dQdi,2*Bus.n,DAE.n); idx = find(m <= m_max & m >= m_min & DAE.f(Sofc.m) ~= 0); if ~isempty(idx) DAE.Fx = DAE.Fx + sparse(Sofc.m,Sofc.m,-1./Tm,DAE.n,DAE.n); DAE.Fy = DAE.Fy + sparse(Sofc.m,Sofc.bus+Bus.n,-Km./Tm,DAE.n,2*Bus.n); dQdm = Vk.*Sofc.con(:,24).*Vs./xt.*sq + ... 0.5*Vs.*Vt./xt.*(2*((xt.*Ik.*Vn./Sn./Vs).^2)./(m.^3))./sq; DAE.Gx = DAE.Gx - sparse(Sofc.bus+Bus.n,Sofc.m,dQdm,2*Bus.n,DAE.n); end % DAE.Gx DAE.Gx = DAE.Gx - sparse(Sofc.bus,Sofc.Ik,Vk.*Vn./Sn,2*Bus.n,DAE.n); DAE.Gx = DAE.Gx - sparse(Sofc.bus,Sofc.Vk,Ik.*Vn./Sn.*Sofc.con(:,24),2*Bus.n,DAE.n); case 5 % non-windup limiters idx = find((m >= m_max | m <= m_min) & DAE.f(Sofc.m) == 0); if ~isempty(idx) global Settings k = Sofc.m(idx); DAE.tn(k) = 0; DAE.Ac(:,k) = 0; DAE.Ac(k,:) = 0; if Settings.octave DAE.Ac(k,k) = -eye(length(idx)); else DAE.Ac(k,k) = -speye(length(idx)); end endend% -------------------------------------------------------------------% function for creating warning messagesfunction fcwarn(idx, msg)fm_disp(strcat('Warning: FC #',int2str(idx),msg))?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -