MES = [1;2;3;4;5;6;7;8;9;10;11;12]; DIADELMES = [31;28;31;30;31;30;31;31;30;31;30;31]; DIADELANO = [16;47;75;106;136;167;197;228;259;289;320;350]; % Gdmo % DECLI: ANGULO DE DECLINACION R = (2*pi/365)*(DIADELANO-1); DECLI = (0.006918-0.399912*cos(R)+0.070257*sin(R)-0.006758*cos(2*R)+0.000907*sin(2*R)-0.002697*cos(3*R)+0.00148*sin(3*R)) ANGULO_DECLINACION = DECLI*180/pi % Gdm0 = radiacion media diaria mensual sobre superficie horizontal Kwh/m2 % Gdm0 =[2.4550 3.0960 4.5840 5.2440 6.7540 7.2070 7.0430 6.3540 5.0690 3.8380 2.4890 2.0780]; % L = latitud despues fi L = 33.22 fi = L*pi/180; % ANGULO DE SALIDA Y PUESTA DE SOL Ws = acos(-tan(DECLI)*tan(fi)) salidaypuesta = Ws*180/pi; keyboard % Z : DURACION DEL DIA Z = salidaypuesta*(2/15) % Radiación media diaria mensual extraterrestre: Gext0 % Gsc = cantidad de energía que llega al exterior de la atmósfea por unidad % de área(perpendicular a los rayos solares) y tiempo; W/m2 Gsc = 1367; Gext = Gsc*(1+0.033*cos(360*DIADELANO/365)); for p=1:1:length(DIADELANO); Gext0(p)=(24/pi)*Gsc*(1+0.033*cos(360*DIADELANO(p)/365))*((cos(fi)*cos(DECLI(p))*sin(Ws(p)))+(Ws(p)*sin(fi)*sin(DECLI(p)))); end % KT índice de nuubosidad mensual o de claridad mensual % radiación del dia promedio horizonttal medido kWh/m2 sobre la radiacion % del dia promedio extraterrestre KT = Gdm0./Gext0'; keyboard % Tal KT, uno para cada més, se usa para determinarD/Gdm0 = relacion entre la componente directa y difusa = dir_dif for n=1:1:length(KT) dir_dif(n) = 1.391-4.0273*KT(n)+5.541*KT(n)^2-3.108*KT(n)^3; end % los horarios de salida y entrada de sol antes cálculados son para una % superficie horizontal. Si la superficie esta inclinada un angulo beta % entonces ambos horarios se modifican y se habla de horario de salida % aparente y entrada aparente del sol: Wsa % incli es el ángulo de inclinacion posible supuesto variable entre ALAFAMIN = 15 y % ALFAMAX = 55º BETAMIN = 15; BETAMAX = 55; incli=BETAMIN:1:BETAMAX; incli = incli*pi/180; for m=1:1:length(MES) for angu=1:1:length(incli) Wsa(angu,m)=acos(-tan(DECLI(m))*tan(fi-(incli(angu)))); end end % se calcula la relacion entre la radiación directa sobre una superficie inclinada (para % distintos incli) y la radiación directa sobre una superficie horizontal (esto es % funcion de la transmitancia atmósferica pero se supone como que la atmósfera no existe for m=1:1:length(MES) for angu=1:1:length(incli) if Ws(m)< Wsa(angu,m) Wn = Ws(m); else Wn = Wsa(angu,m); end Num = (cos(fi-incli(angu))*cos(DECLI(m))*sin(Wn))+(Wn*sin(fi-incli(angu))*sin(DECLI(m))); Den = (cos(fi)*cos(DECLI(m))*sin(Ws(m)))+(Ws(m)*sin(fi)*sin(DECLI(m))); Rd(angu,m)=Num/Den; end end keyboard %buscar el angulo para el cual Rd es máximo y el valor correspondiente de %Rd for t=1:1:length(MES) Vec=Rd(:,t); [Y,I]=max(Vec); VALORMAX(t) = Y; BETAMAX(t) = incli(I); end keyboard % se pasa a calcular la radiacion total RTOTMAX como la suma de la RDIRE más % la RDIFU mas la RALBE, para cada mes y para el ángulo que ha hecho máximo % Rd; ro = indice de reflectancia del terreno ro = 0.2; for z=1:1:length(MES) RDIRE(z)=(1-dir_dif(z))*VALORMAX(z); RDIFU(z)=dir_dif(z)*((1+cos(BETAMAX(z)))/2); RALBE(z)=ro*((1-cos(BETAMAX(z)))/2); RTOTMAX(z)=RDIRE(z)+RDIFU(z)+RALBE(z); end keyboard % una aproximación para determinar cual seria el mejor angulo de % inclinacion fija sería con el valor medio de los angulos que dan el % máximo : el valor dedio de BETAMAX BETAFIJO = mean(BETAMAX)*180/pi; % y en tal caso la radiacion total para tal angulo sería RTOTMED %los valores de Rd seran loa correspondientea a tal ángulo BETA = round(BETAFIJO); FILA = BETA+1-BETAMIN; for z=1:1:length(MES) RDIREMED(z)=(1-dir_dif(z))*Rd(FILA,z); RDIFUMED(z)=dir_dif(z)*((1+cos(BETA*pi/180))/2); RALBEMED(z)=ro*((1-cos(BETA*pi/180))/2); RTOTMED(z)=RDIREMED(z)+RDIFUMED(z)+RALBEMED(z); end keyboard % para el caso de angulo fijo (sin seguidor) el calculo de la produccion de % generador, en este caso diario, implica multiplicar el RTOTMED por Gdm0 GdmFIJO = RTOTMED'.*Gdm0 % Rendimiento del seguidor solar para cada mes: eta = [26.02 21.32 28.13 29.74 33.91 41.00 44.00 38.61 27.88 26.72 23.69 21.58]; % con lo que los kWh/m2 por dia promedio para cada mes con seguidor sera: Gdmseg = GdmFIJO.*eta/100+GdmFIJO; keyboard % si se compara esto con Gdm0, tal incremento porcentual será: INCREMENTO = (Gdmseg-Gdm0)*100./Gdm0; keyboard % y el incremento porcentual promedio INCREPROMEDIO = mean(INCREMENTO); keyboard %graficamente se muestra Gdm0, GdmFIJO y Gdmseg a lo largo del año figure(1) plot(MES,Gdm0,'b',MES,GdmFIJO,'r',MES,Gdmseg,'g');grid xlabel('MES DEL AÑO'); ylabel('KwH/DIA M2'); %legend('HORIZONTAL','angulo fijo','seguidor'); % los valores mensuales de radiacion incidente (W/m2) para una superficie % inclinada y para una superficie con seguidor los kwh/m2 por dia llevados a mes y % dividido por la cantidad de horas del dia, o, lo que es lo mismo los % kWh/m2 dia dividido 24 Wm2fijo = GdmFIJO/24; Wm2seg = Gdmseg/24; keyboard % rendimientos de los módulos %un dato fundamental es la temperatura ambiente % se plantea que se conoce la temperatura ambiente media diurna Tamd para % cada dia representativo de cada mes Tamb = [12.7 14.5 17.8 19.1 22.4 26.6 28.8 29.0 25.6 22.0 16.6 13.6]; % se supone, en realidad es un dato dado por el fabricante d elos paneles, % que la variación de la tension de circuito abierto (dVco) con la temperatura de la célula (dTc) es % dVco_dTc en mv/ºC (son los voltios que varía la tal tension por cada % grado de incremento de la temp., en este caso variación negativa dVco_dTc = -0.129; % otro dato propio del panel es el TONC (Temperatura de Operacion Nominal de la Celula, a % 800 W/m2, distribucion espectral 1,5, 20ºC y velocidad del viento de 1 m/s) TONC = 45; % A partir de lo cual se puede calcular una constante C C = (TONC-20)/800; % con esto se puede calcular la temp de operación de celula (Tc), en este caso % para la temp media mensual diurna y de la intensidad de radiación solar % media mensual. A partir de esto se calcula la variación de la tension de % co DeltaVco % PARA EL CASO FIJO for z=1:1:length(MES) TcF(z) = C*Wm2fijo(z)+Tamb(z); % variacion de la temperatura respecto de la condicion estandard DeltaTF(z) = TcF(z)-25; % Variacion de la Vco por temperatura DeltaVcoF(z) = dVco_dTc*DeltaTF(z); end % PARA EL CASO CON SEGUIDOR for z=1:1:length(MES) TcS(z) = C*Wm2seg(z)+Tamb(z); % variacion de la temperatura respecto de la condicion estandard DeltaTS(z) = TcS(z)-25; % Variacion de la Vco por temperatura DeltaVcoS(z) = dVco_dTc*DeltaTS(z); end % se supone que en funcion de tales temperaturas se determina el % rendimiento de la célula (objeto de otro problema); aqui se toma los % siguientes valores de rendimiento para cada mes % rendi Fijo etaTEM_F=[1.03 1.02 1.00 0.99 0.98 0.96 0.95 0.95 0.96 0.98 1.01 1.03]; % rendi Seguidor etaTEM_S=[1.02 1.01 0.99 0.98 0.96 0.94 0.93 0.93 0.95 0.97 1.01 1.02]; % SE SIGUE CON LOS RENDIMIENTOS % % PERDIDAS POR PUNTO DE TRABAJO: SE CONSIDERA LA POSIBILIDAD DE DESVIACION % DE TRABAJA EN EL "PUNTO DE MAXIMA POTENCIA" O PREICISION DEL DISPOSITIVO % PARA TRABAJAR EN TAL PUNTO; SE CONSIDERA QUE LA MAXIMA PERDIDA POR ESTE % MOTIVOP PUEDE SER DEL 2% etaPT = 0.98; % PERDIDAS POR DISTORSION O CONEXIONADO: LOS PANELES TIENEN SUS PROPIAS % DESVIACIONES RESPECTO DEL VALOR NOMINLA DE POTENCIA, PUEDE QUE EL % FABRICANTE INDIQUE CUAL ES TAL DISPERSION, SE SUPONE UN 5% etaDIS = 0.95; % PERDIDAS POR SUCIEDAD: PUEDEN SER NULASDESPUES DE UNA LLUVIA O DE HASTA % UN 8% CON MUCHO POLVO, SE SOMA UN VALOR MEDIO DE 3% etaSUC = 0.93; % PERDIDAS POR REFLECTANCIA ANGULAR: PUEDE SER NULA PARA EL MEDIODIA Y UN % PAR DE HORAS ALREDEDOR DEL MEDIODIA O TAMBIEN PARA EL CASO DE SEGUIDORES % A DOS EJES. ES MAS IMPROTANTE A MEDIDA QUE CRECE LA LATITUD; POR LO TANTO % SON DISTINTOS SEGUN SEA LA INSTALACION FIJA O CON SEGUIDOR etaREF_F = 0.97; etaREF_S = 0.97; % PERDIDAS POR SOMBREADO: SUPONIENDO QUE SE EVITA LA DE LOS PROPIOS PANELES % QUEDARIA POR REVISAR LA DE OTROS OBSTACULOS O LA DE LOS SEGUIDORES; EN % ESTE ÇULTIMO CASO SE SUPONE INEVITABLE Y SE AFECTA LIGERAMENTE: etaSOM_F = 1; etaSOM_S = 0.97; % RENDIMIENTOS DE OTROS COMPONENTES % RENDIMIENTO DEL INVERSOR: EL FABRICANTE DA LOS RENDIMIENTOS Y % EVENTUALMENTE PARA DISTINTOS ESTADOS DE CARGA etaINV = 0.96; % RENDIMIENTO DEL CABLEADO: PERDIDAS OHMICAS TANTO EN CC (1.5%) COMO EN CA % (2%) etaCC = 0.985; etaCA = 0.98; etaCAB = etaCC*etaCA; % RENDIMIENTO DEL TRANSFORMADOR etaTRAFO = 0.97; % FACTOR DE RENDIMIENTO GLOBAL (PR): PRODUCTO DE TODOS LOS RENDIMIENTOS % PANELES FIJOS PR_F = etaPT*etaDIS*etaSUC*etaREF_F*etaSOM_F*etaINV*etaCC*etaCAB*etaTRAFO*etaTEM_F; % SEGUIMIENTO PR_S = etaPT*etaDIS*etaSUC*etaREF_F*etaSOM_F*etaINV*etaCC*etaCAB*etaTRAFO*etaTEM_S; figure(2) plot(MES,PR_F,'r',MES,PR_S,'b');grid xlabel('MES DEL AÑO') ylabel('RENDI GLOBAL') keyboard % el rendimiento energético total PRMedio_F = mean(PR_F)*100; PRMedio_S = mean(PR_S)*100; % se esta en condiciones de determnar la producion de energia del generador % fotovoltaico en funcion de la POTENCIA PICO TOTAL INSTALADA kWp_INS kWp_INS = 1000; %1MWpico instalado % Energía diaria EDIARIA_F = kWp_INS*GdmFIJO.*PR_F; EDIARIA_S = kWp_INS*Gdmseg.*PR_S; % energia mensual EMES_F = EDIARIA_F.*DIADELMES; EMES_S = EDIARIA_S.*DIADELMES; % ENERGIA ANUAL EANUAL_F = sum(EMES_F) EANUAL_S = sum(EMES_S) %RATIO_F = EANUAL_F/kWp_INS %RATIO_S = EANUAL_S/kWp_INS