clear clc /////////////////////////////////////////////////////////////////////////////// // PME3380 - Modelagem de Sistemas Dinâmicos //Grupo 05 - Modelagem de Turbina Eólica // Felipe Todaro Fleischmann.........9795293 // Mariana Claudino Pin..............9348664 // Paulo Montijo Bandeira............9348449 //Vitor Olavo Tonaco Alexandre.......9836176 /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// //A modelagem é dividida em 3 módulos: 1. Aerodinâmico, 2. Mecânico e 3. Elétrico //1. Módulo aerodinâmico //Consiste em modelar como o vento entrega energia e torque para as pá da turbina //eólica. //2. Módulo Mecânico //Consiste em modelar como a energia proveniente do vento é convertida em um //torque mecânico, que por sua vez é transmitido ao gerador por meio de uma caixa // de engrenagens funcionando como caixa redutora //3. Módulo Elétrico //Consiste em modelar o gerador, que faz a conversão da energia mecânica em elé- //trica para uso na rede elétrica // //Objetivo: estudar a geração de energia elétrica como base em alterações do //vento para posterior controle da potência gerada pela turbina, a ser mantida o //mais próximo possível da potência nominal sem acrretar danos às partes mecânica //e elétrica do sistema /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// //1. Módulo Aerodinâmico //Parâmetros do modelo aerodinâmico rho = 1.205; //densidade do ar [kg/m3] raio=38; //raio do rotor [m] area=raio^2*%pi; //Conversão da potência aerodinâmica em potência mecânica //Parâmetros para cálculo do coeficiente de potência c1 = 0.5; c2 = 116; c3 = 0.4; c4 = 0; c5=5; c6 = 21; betha = 12.1; // ângulo de pitch /////////////////////////////////////////////////////////////////////////////// //2. Módulo Mecânico //Parâmetros do modelo mecânico Jt =3*10^6; //Momento de inércia da turbina [kgm^2] Kt = 30*10^7; //Elasticidade do eixo da turbina [N/rad] Bt = 13*10^6; //Amortecimento do eixo da turbina [Ns/m] n = 67; //Fator de transmissão /////////////////////////////////////////////////////////////////////////////// //3. Módulo Elétrico //Parâmetros do modelo elétrico Jg = 130; //Momento de inércia do gerador [kgm^2] La = 0.0259; //indutância da armadura [H] Ra = 3.2919; //Resistência da armadura [Ohms] Km = 0.0928*100; //Constante eletromotriz [Vs/rad] Va=600; /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// // Parâmetros iniciais e de simulação v_vento_operacao = 12.1; Va = 600; u_operacao=[v_vento_operacao; Va] x_operacao=[0;1.27;0;1.27*n;u_operacao(2)/Ra] x0=x_operacao; // Tempo ti=0; tf=30; passo=0.01; t = [ti:passo:tf]; /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// //4. Modelo não-linear function v_w=vel(t) //T_g0=3; //T_g=2; //u=20; v_w=12.1; // V_BASE //if t>T_g0 then v_w=12.1+0.5*u*(1-cos(3.14*2*(t-T_g0)/T_g)) //end //if t>T_g0+T_g then v_w=12.1 //end //if t> T_g0 then v_w = 12.1+(0.5*cos(%pi/3*t)+0.3*cos(%pi/2*t+%pi/10)) // end //if t>20 then v_w=12.1+(t-20)/4 //end //if t>30 then v_w=14.6 //end endfunction function [y_ponto] = turbina(t,y) lambda=y(2)*raio/vel(t); i= (1/(lambda+0.08*betha)-0.035/(betha^3+1))^-1; Cp = c1*(c2/i - c3*betha-c5)*exp(-c6/i) y_ponto1=y(2); y_ponto2=(%pi*0.5*rho*raio^2*vel(t)^3*Cp/y(2)-Kt*(y(1)-y(3)/n)-Bt*(y(2)-y(4)/n))/Jt; y_ponto3=y(4); y_ponto4=(-y(5)*Km+(Kt*(y(1)-y(3)/n)+Bt*(y(2)-y(4)/n))/n)/Jg; y_ponto5=(y(4)*Km-Ra*y(5)-Va)/La; y_ponto=[y_ponto1;y_ponto2;y_ponto3;y_ponto4;y_ponto5]; endfunction /// Simulação sol = ode(x0, t(1), t, turbina); // Parâmetros aerodinâmicos for j=1:1:length(t) velw(j)=vel(t(j)); lambda(j)=sol(2,j)*raio/velw(j); i(j)= (1/(lambda(j)+0.08*betha)-0.035/(betha^3+1))^-1; Cp(j) = c1*(c2/i(j) - c3*betha-c5)*exp(-c6/i(j)); P_w(j)= %pi*0.5*rho*raio^2*(velw(j))^3*Cp(j);// Potência aerodinâmica P_ar(j)= %pi*0.5*rho*raio^2*(velw(j))^3;// Potência total T_w(j)=P_w(j)/sol(2,j); // Torque aerodinâmico P_e(j)=Km*sol(4,j)*sol(5,j); // Potência gerada end // Parâmetros aerodinâmicos de referência lambda_i = 2; lambda_f = 16; h_l = 0.001; lambda_r = [lambda_i:h_l:lambda_f]; betha_r = [0;5;10;15;20]; i_r=zeros(length(lambda_r),length(betha_r)) for p=1:1:length(lambda_r) for q=1:1:length(betha_r) i_r(p,q)= (1/(lambda_r(p)+0.08*betha_r(q))-0.035/(betha_r(q)^3+1))^-1; Cp_r(p,q) = c1*(c2/i_r(p,q) - c3*betha_r(q)-c5)*exp(-c6/i_r(p,q)); end end /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// // 5. Modelo Linear derivative1=-(area*c2*rho*raio*u_operacao(1)^2*%e^(-c6*(1/((raio*x_operacao(2))/u_operacao(1) + 1.088) - 1.74464*10^(-5)))*c1*(c2*(1/((raio*x_operacao(2))/u_operacao(1) + 1.088) - 1.74464*10^(-6)) - c5))/(2*Jt*x_operacao(2)*((raio*x_operacao(2))/u_operacao(1) + 1.088)^2) - (area*rho*u_operacao(1)^3*%e^(-c6*(1/((raio*x_operacao(2))/u_operacao(1) + 1.088) - 1.74464*10^(-6)))*c1*(c2*(1/((raio*x_operacao(2))/u_operacao(1) + 1.088) - 1.74464*10^(-6)) -c5))/(2*Jt*x_operacao(2)^2) + (area*c6*rho*raio*u_operacao(1)^2*%e^(-c6*(1/((raio*x_operacao(2))/u_operacao(1) + 1.088) - 1.74464*10^(-5)))*c1*(c2*(1/((raio*x_operacao(2))/u_operacao(1) + 1.088) - 1.74464*10^(-6)) - c5))/(2*Jt*x_operacao(2)*((raio*x_operacao(2))/u_operacao(1) + 1.088)^2) A=[0 1 0 0 0; -Kt/Jt derivative1-Bt/Jt Kt/(n*Jt) Bt/(n*Jt) 0; 0 0 0 1 0; Kt/(n*Jg) Bt/(n*Jg) -Kt/(n^2*Jg) -Bt/(n^2*Jg) -Km/Jg; 0 0 0 Km/La -Ra/La]; derivative2=(area*c2*rho*raio*u_operacao(1)*%e^(-c6*(1/((raio*x_operacao(2))/u_operacao(1)+1.088)-1.74464*10^(-5)))*c1*(c2*(1/((raio*x_operacao(2))/u_operacao(1)+1.088)-1.74464*10^(-6))-c5))/(2*Jt*((raio*x_operacao(2)/u_operacao(1)+1.088)^2))+(3*area*rho*u_operacao(1)^2*%e^(-c6*(1/((raio*x_operacao(2))/u_operacao(1)+1.088)-1.74464*10^(-5)))*c1*(c2*(1/((raio*x_operacao(2))/u_operacao(1)+1.088)-1.74464*10^(-6))-c5))/(2*Jt*x_operacao(2))-(area*c6*rho*raio*u_operacao(1)*%e^(-c6*(1/((raio*x_operacao(2))/u_operacao(1)+1.088)-1.74464*10^(-5)))*c1*(c2*(1/((raio*x_operacao(2))/u_operacao(1)+1.088)-1.74464*10^(-6))-c5))/(2*Jt*((raio*x_operacao(2))/u_operacao(1)+1.088)^2) B = [0 0; derivative2 0; 0 0; 0 0; 0 -1/La] C=[0 1 0 0 0; 0 0 0 1 0; 0 0 0 0 1; 0 0 0 Km*x_operacao(5) Km*x_operacao(4)]; turbina_linear = syslin('c',A,B,C) /////////////////////////////////////////////////////////////////////////////// // 6. Simulação // Entradas B1 = B(:,1);// Velocidade do vento B2=B(:,1);// Tensão // Saídas C1 = C(1,:); // Velocidade angular do rotor C2 = C(2,:); // Velocidade angular do gerador C3 = C(4,:);// Potência gerada C4 = C(3,:);// Corrente da armadura D = 0; turbina_vento_w_t = syslin('c',A,B1,C1,D);// turbina_vento_w_g = syslin('c',A,B1,C2,D); turbina_vento_pot= syslin('c',A,B1,C3,D); turbina_vento_corr = syslin('c',A,B1,C4,D); [delta_w_t,delta_x_w_t] = (csim(velw',t,turbina_vento_w_t)); [delta_w_g,delta_x_w_g] = (csim(velw',t,turbina_vento_w_g)); [delta_pot,delta_x_pot] = (csim(velw',t,turbina_vento_pot)); [delta_corr,delta_x_corr] = (csim(velw',t,turbina_vento_corr)); y_w_t=delta_w_t+x0(2); y_w_g=delta_w_g+x0(4); y_pot=delta_pot+x0(4)*x0(5)*Km; y_corr=delta_corr+x0(5); for j=1:1:length(t) betha= 0.2105*180/%pi; lambda_l(j)=y_w_t(j)*raio/velw(j); i(j)= (1/(lambda_l(j)+0.08*betha)-0.035/(betha^3+1))^-1; Cp_l(j) = c1*(c2/i(j) - c3*betha-c5)*exp(-c6/i(j)); end //Diferença entre modelo linear e não linear diff_wrot=sol(2,:)-y_w_t; diff_wg=sol(4,:)-y_w_g; diff_pot=P_e'-y_pot; /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// // 6. Função de transferência FT1=ss2tf(turbina_vento_w_t);// Função de transferência entrada vento e saída velocidade do rotor FT2=ss2tf(turbina_vento_pot);// Função de transferência entrada vento e saída Potência gerada /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////// // 7. Gráficos set(gda(), "grid_style",[7 7], "grid", [1 1]*color("grey70")); title(gda(), "fontsize", 3, "fontname", "helvetica bold"); xlabel(gda(), "fontsize", 2,"color", "black", "fontname", "helvetica"); ylabel(gda(), "fontsize", 2, "fontname", "helvetica"); xset('window', 1) plot(t,y_w_t,'k','LineWidth',2);//linear plot(t,sol(2,:),'r','LineWidth',2);// não linear //plot([15 15],[1 6],'b-.') // Traço para visuzalizar //plot([19 19],[1 6],'b-.')// Traço para visuzalizar hl=legend(['Linear';'Não linear'],4); xtitle('Velocidade angular do rotor no tempo','tempo(s)','Velocidade angular do rotor (rad/s)'); xset('window', 2) plot(t,y_w_g, 'k','LineWidth',2);//linear plot(t,sol(4,:),'LineWidth',2);//não linear hl=legend(['Linear';'Não linear'],4); xtitle('Velocidade angular do gerador no tempo','tempo(s)','Velocidade angular do gerador (rad/s)'); xset('window', 3) plot(t,y_pot,'r','LineWidth',2);//linear plot(t,P_e,'k','LineWidth',2);// não linear plot(t,P_ar,'b','LineWidth',2);// não linear //plot(t,P_w,'g','LineWidth',2);// não linear hl=legend(['Potência linear gerada';'Potência não linear gerada';'Potência aerodinamica'],2); xtitle('Potência','tempo(s)','Potência (W)'); xset('window', 4) plot(t,y_corr,'r','LineWidth',2);//linear hl=legend(['Linear';'Não linear'],4); xtitle('Corrente da armadura no tempo','tempo(s)','Corrente (A)'); xset('window', 5) plot(lambda_r,Cp_r(:,:),'LineWidth',2);// Cp padrão hl=legend(['Tetha=0º';'Tetha=5º';'Tetha=10º';'Tetha=15º';'Tetha=20º'],1); a=get("current_axes") a.data_bounds=[2,-0.1;16,0.5]; xtitle('Coeficiente de Potência','Lambda','Cp'); xset('window', 6) plot(t,velw,'k','LineWidth',2);//velocidade do vento simulada //a=get("current_axes") //a.data_bounds=[2,-0.1;6,50]; xtitle('Velocidade do vento','tempo(s)','Velocidade (m/s)'); xset('window', 7) plot(lambda,Cp,'b',lambda(1),Cp(1),'r-*','LineWidth',2)// Cp simulado xtitle('Coeficiente de Potência','Lambda','Cp')