C*************************************************************** SUBROUTINE TYPE62(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C Version from: 11/16/89 BY J. ECKSTEIN C SAMPLE INPUT AND OUTPUT ADDED BY N. BLAIR JULY 1992 C--------------------------------------------------------------- C This subroutine represents a four parameter model of C a Photovoltaic array. It is capable to predict the C complete current-voltage characteristic over the entire C operating voltage range of a flat-plate, non-sunconcen- C trated collector. While a series resistance is taken C into account, a shunt resistance is assumed to be infinite C and thus neglected in the model. C A routine is implemented which determines the voltage and C current at maximum power point. C The operating current is found for a given voltage, C irradiance and ambient temperature as input. C An option is provided to evaluate the series resistance with C the bisection method, if not given as input. C To overcome convergence problems which appear when simul- C ating direct coupled systems, a bisection method type of C convergence promotion is included. It is turned off, if C input FLAG is equal zero. C************************************************************** C************************************************************** C**** Definition of the variables: C**** TRNSYS specific variables: C XIN == input array C OUT == output array C PAR == parameters C TIME == simulation time C T,DTDT == not used in this component C S == storage array C NSTORE == dimension of S C IAV == pointer within S C ISTORE == index C INFO == array to use TRNSYS internal information C**** component specific variables: C A == completion factor C AREA == collector area [m^2] C CURRENT == function called current C DUMMY == auxiliary variable for convergence promotion C EFFREF == reference max. power efficiency C EG == bandgap energy [eV] C FF == fill factor C FIRST == auxiliary variable; allows that at first timestep C with solar radiation, V = VMP, otherwise the first C guess for V is V from the previous timestep C FLAG == switch: if FLAG=1 then convergence promotion is on C if FLAG=0 then convergence promotion is off C GAM == curve fit factor C I == current [amps] C IL == light current [A] C ILR == reference light current [A] C IMP == current at max. power [A] C IMR == reference current at max. power [A] C IO == reverse saturation current [A] C IOR == dito at reference [A] C ISC == short circuit current [A] C ISCR == dito at reference [A] C MEMO == memorizer if convergence promotion is on or off C MISC == temperature coefficient: short circuit current [A/K] C MVOC == temperature coefficient: open circuit voltage [V/K] C N == pointer to a relative adress in program C NP == number of modules in parallel C NS == number of modules in series C P == power [W] C PMAX == power [W] C Q_BZ == electron charge/Boltzmann constant [C*K/J] C TA == ambient temperature [Kelvin (K)] C TANOCT == ambient temperature at NOCT [K] C TAU_AL == transmittance-absorptance product C TC == cell temperature [K] C TCR == reference cell temperature [K] C TCNOCT == cell temperature at NOCT [K] C NCS == NUMBER OF CELLS IN SERIES WITHIN THE SYSTEM C SUN == irradiance [W/m^2] C SUNR == reference irradiance [W/m^2] C SUNNOCT == irradiance at NOCT [K] C V == voltage [volts] C VMP == max. power point voltage [V] C VMR == reference max. power point voltage [V] C VOC == open circuit voltage [V] C UTIL == utilization of the array: ratio of actual power C to max. power C**** variables used to determine the series resistance C ANEW,ALOW,AUP == A is the completion factor, the indexes C stands for the limits: lower and upper, C and for the current value: new C FNEW,FLOW,FUP == objective functions: at the interval C limits and at the current value C GAMNEW,GAMLOW,GAMUP == curve fit factor: at the interval C limits and at the current value C IONEW,IOLOW,IOUP == saturation current: at the interval C limits and at the current value C RS,RSNEW,RSLOW,RSUP == series resistance: at the interval C limits and at the current value C**** variables used in operating current calculation C F == objective function for Newton's method C FPRIME == first derivative of F C IOLD,INEW == iteration variables C**** variables used in maximum power evaluation C F1 == objective function C F1P == first derivative of F1 C IMXO,IMXN == iteration variables C C****************************************************** C C SAMPLE INPUT AND OUTPUT C C INPUTS SAMPLE VALUES C XIN(1) SUN (*3.6) 1400.0 C XIN(2) TA 24.0 C XIN(3) V 4.6 C XIN(4) FLAG 1.0 C C PARAMETERS C PAR(1) ISCR 1.82 C PAR(2) VOC 21.3 C PAR(3) TCR 298.0 C PAR(4) SUNR 1000.0 C PAR(5) VMR 17.8 C PAR(6) IMR 1.68 C PAR(7) MISC .0015 C PAR(8) MVOC -.073 C PAR(9) NCS 36.0 C PAR(10) NS 3.0 C PAR(11) NP 4.0 C PAR(12) TCNOCT 318.0 C PAR(13) TANOCT 293.0 C PAR(14) SUNNOCT 1000.0 C PAR(15) AREA .259 C PAR(16) TAU_AL .9 C PAR(17) EG .1 C PAR(18) RS 30.0 C C OUTPUTS C OUT(1) V -8.255 C OUT(2) I -0.4058 C OUT(3) P 3.35 C OUT(4) PMAX 3.35 C OUT(5) UTIL 1.0 C OUT(6) VMP -8.255 C OUT(7) IMP -0.4058 C OUT(8) VOC -16.44 C OUT(9) ISC 2.212 C OUT(10) FF 0.0 C C***************************************************************** C**** declaration of the variables: IMPLICIT NONE REAL PAR,TIME,T,DTDT,ICNTRL,S DOUBLE PRECISION XIN, OUT REAL I,CURRENT,V,IO,IOR,TA, & A,Q_BZ,GAM,IL,ISC,TC,VOC,NS,NP, & ISCR,VOCR,TCR,SUNR,VMR,IMR,MISC,EG, & ILR,SUN,SERCELL,TCNOCT,TANOCT,SUNNOCT, & AREA,TAU_AL,EFFREF,IMP,VMP,P,PMAX,MVOC, & UTIL,VMI,VMA,VOLD,FF REAL ANEW,ALOW,AUP,FNEW,FLOW,FUP,GAMNEW,GAMLOW,GAMUP, & IONEW,IOLOW,IOUP,RS,RSNEW,RSLOW,RSUP REAL IOLD,INEW,F,FPRIME REAL IMXO,IMXN,F1,F1P INTEGER INFO,FLAG,DUMMY,ISTORE,NSTORE,IAV,MEMO,N,FIRST DIMENSION XIN(4), OUT(10), PAR(18), INFO(15) C****************************************************************** COMMON /STORE/ NSTORE,IAV,S(5000) C**** store is used to store values from previous timestep COMMON /ARRAY/ GAM,TC,Q_BZ,IL,IO,RS,IMP C**** array is used to transfer data to function 'current' COMMON /PARAM/ SERCELL,TCR,IMR,VMR,ISCR,VOCR,MVOC,MISC,EG C**** param is used to transfer data to subroutine 'series' C***************************************************************** C----------------------------------------------------------------- C**** Set inputs SUN=XIN(1)/3.6 TA=XIN(2) V=XIN(3) FLAG=XIN(4) C----------------------------------------------------------------- C**** in this section a couple of checks on the inputs are C**** performed. This is done at the second and following C**** calls in timestep IF(INFO(7).GT.0)THEN DUMMY=INT(S(ISTORE+3)) MEMO=INT(S(ISTORE+4)) N=INT(S(ISTORE+5)) C**** all calculations are being scipped during time periods C**** with no insolation IF(SUN.EQ.0.)THEN V=0. I=0. C**** jump to the output section of the program GOTO 1000 ELSE C**** check on operation mode IF(FLAG.EQ.0)THEN C**** convergence promotion is off C**** if voltage is negative in system with battery storage, C**** cell is disconnected, i.e. output is zero current IF(V.LT.0.)THEN V=0. I=0. ELSE C**** normal operation C**** check if voltage greater than open circuit voltage VOC=OUT(8) IF(V.GT.VOC)THEN V=VOC I=0. ELSE I=CURRENT(V) ENDIF ENDIF C**** jump to the output section of the program GOTO 1000 ELSE C**** this part is just invoked if a system with battery C**** is operated and the battery has to be disconnencted C**** temporarily to protect the battery from overcharging C**** or deep discharge. That means the system is changed C**** to a direct coupled system. Since this can happen C**** after a few iterations running the battery system, C**** the convergence promotor, i.e. the number of C**** iteration involved in the algorithm has to be adressed C**** relatively using N IF(MEMO.EQ.0)THEN N=INFO(7) MEMO=1 V=VMP I=IMP GOTO 1000 ENDIF ENDIF ENDIF ENDIF C---------------------------------------------------------------- C**** initial call in simulation IF (INFO(7).EQ.-1) THEN INFO(6)=10 INFO(9)=0 C**** storage allocation: values from previous iteration C**** are stored in the S-array C**** S(ISTORE) -- VOLD is voltage from previous iteration C**** S(ISTORE+1) -- VMI is lower voltage limit of converg- C**** ence integral C**** S(ISTORE+2) -- VMA is upper limit on interval C**** S(ISTORE+3) -- DUMMY C**** S(ISTORE+4) -- MEMO C**** S(ISTORE+5) -- N INFO(10)=6 CALL TYPECK(1,INFO,4,18,0) ISTORE=INFO(10) S(ISTORE)=0. S(ISTORE+1)=0. S(ISTORE+2)=0. S(ISTORE+3)=0. S(ISTORE+4)=0. S(ISTORE+5)=0. Q_BZ=11604.45 C**** set parameters *************************** ISCR=PAR(1) VOCR=PAR(2) TCR=PAR(3) SUNR=PAR(4) VMR=PAR(5) IMR=PAR(6) MISC=PAR(7) MVOC=PAR(8) SERCELL=PAR(9) NS=PAR(10) NP=PAR(11) TCNOCT=PAR(12) TANOCT=PAR(13) SUNNOCT=PAR(14) AREA=PAR(15) TAU_AL=PAR(16) EG=PAR(17) RS=PAR(18) IF (RS.LT.0) THEN C**** in this case RS is not provided as a parameter and C**** has to be evaluated. This is done in subroutine "series" CALL SERIES(RS,Q_BZ) ENDIF C**** evaluation of the 3 unknowns at reference condition GAM=Q_BZ*(VMR-VOCR+IMR*RS)/(TCR*LOG(1.-IMR/ISCR)) ILR=ISCR IOR=ILR/EXP(Q_BZ*VOCR/(GAM*TCR)) C**** set up parameters for the entire array ILR=NP*ILR IOR=NP*IOR GAM=NS*GAM RS=(NS/NP)*RS A=GAM/(NS*SERCELL) MEMO=0 FIRST=0 GOTO 1000 C----------------------------------------------------------------- C**** first call in timestep ELSE IF (INFO(7).EQ.0) THEN VOLD=S(ISTORE) C**** evaluation of cell temperature from NOCT conditions EFFREF=IMR*VMR/(SUNR*AREA) TC=TA +(SUN*(TCNOCT-TANOCT)/SUNNOCT*(1.-EFFREF/TAU_AL)) C**** this part calculates how IL and IO vary with temp. and SUN IL=(SUN/SUNR)*(ILR+MISC*NP*(TC-TCR)) IF(IL.LT.0.0) IL=0.0 IO=IOR*((TC/TCR)**3)*EXP((Q_BZ*EG/(A))*((1./TCR)-(1./TC))) C**** Open circuit voltage VOC=GAM*TC/Q_BZ*LOG(IL/IO+1.) C**** Short circuit current ISC=IL C**** this part calculates the max-power current & voltage C**** using NEWTON's method to solve for the max-power voltage C**** and max-power current IMXO=0.0 C**** first guess for max. power point current IMXN=SUN/SUNR*NP*(IMR+MISC*(TC-TCR)) DO WHILE (ABS(IMXN-IMXO).GT.0.001) IMXO=IMXN F1=IMXO+(IMXO-IL-IO)*(LOG((IL-IMXO+IO)/IO)-IMXO*Q_BZ*RS/ > (GAM*TC))/(1.+(IL-IMXO+IO)*(Q_BZ*RS/(GAM*TC))) F1P=2.+(LOG((IL-IMXO+IO)/IO)-(Q_BZ*RS*IMXO/(GAM*TC)))/ > ((1.+(IL-IMXO+IO)*(Q_BZ*RS/(GAM*TC)))**2.) IMXN=IMXO-(F1/F1P) END DO IMP=IMXN VMP=LOG(1.+(IL-IMP)/IO)*(TC*GAM/Q_BZ)-IMP*RS IF ((VOC.GT.0.).AND.(ISC.GT.0)) THEN C**** Fill factor FF=VMP*IMP/VOC/ISC ELSE FF=0. ENDIF C**** all calculations are being scipped during time periods C**** with no insolation IF(SUN.EQ.0.)THEN V=0. I=0. ELSE C**** check on operation mode IF(FLAG.EQ.0)THEN C**** convergence promotion is off C**** if voltage is negative in system with battery storage, C**** cell is disconnected, i.e. output is zero current IF(V.LT.0.)THEN V=0. I=0. ELSE C**** normal operation C**** check if voltage greater than open circuit voltage IF(V.GT.VOC)THEN V=VOC I=0. ELSE I=CURRENT(V) ENDIF ENDIF ELSE C**** convergence promotion is on, i.e. direct coupled mode MEMO=1 N=INFO(7) IF(FIRST.EQ.0)THEN C**** at the first timestep of the simulation at which C**** solar radiation occurs, the first guess on V is C**** V at max. power point V=VMP I=IMP FIRST=1 ELSE C**** the first guess foe each timestep on V is the C**** voltage at operating point from the previous C**** timestep V=VOLD I=CURRENT(V) ENDIF ENDIF ENDIF PMAX=IMP*VMP DUMMY=0 GOTO 1000 C------------------------------------------------------------------ C**** second call in timestep, first iterative call ELSE IF (INFO(7).EQ.(N+1)) THEN C**** here begins the convergence promotion algorithm C**** the first guess was the voltage from the previous C**** timestep. The voltages for the next two iterations C**** are tracked before the promoter intervenes. Until C**** then a protection prevents the voltage to jump C**** out of bounds, i.e. if the input voltage is C**** greater than open circuit voltage it is set equal C**** to it. The convergence promotion algorithm is C**** basically the same as the bisection method, it C**** cuts an interval in half until a sufficient tolerance C**** is obtained. VOLD=S(ISTORE) VOC=OUT(8) IF(V.GT.VOLD)THEN IF(V.GT.VOC)THEN DUMMY=4 VMI=VOLD VMA=VOC V=(VMI+VMA)/2. ELSE DUMMY=1 VMI=VOLD ENDIF ELSE DUMMY=2 VMA=VOLD ENDIF I=CURRENT(V) GOTO 1000 C----------------------------------------------------------------- C**** third call in timestep ELSE IF(INFO(7).EQ.(N+2))THEN VOLD=S(ISTORE) VOC=OUT(8) IF(DUMMY.EQ.1)THEN IF(V.GT.VOLD)THEN DUMMY=3 ELSE VMA=VOLD IF(V.GT.VMI)THEN VMI=V ENDIF V=(VMI+VMA)/2. ENDIF ELSEIF(DUMMY.EQ.4)THEN IF(V.GT.VOLD)THEN VMI=VOLD IF(V.GT.VOC)THEN VMA=VOC ELSE VMA=V ENDIF V=(VMI+VMA)/2. ELSE VMA=VOLD IF(V.GT.VMI)THEN VMI=V ENDIF V=(VMI+VMA)/2. ENDIF ELSE IF(V.GT.VOLD)THEN VMI=VOLD IF(V.LT.VMA)THEN VMA=V ENDIF V=(VMI+VMA)/2. ELSE DUMMY=3 ENDIF ENDIF I=CURRENT(V) GOTO 1000 C---------------------------------------------------------------- C**** fourth and following calls in timestep ELSE IF(DUMMY.NE.3)THEN VMI=S(ISTORE+1) VMA=S(ISTORE+2) VOLD=S(ISTORE) VOC=OUT(8) IF(V.GT.VOLD)THEN VMI=VOLD IF(V.LT.VMA)THEN VMA=V ENDIF ELSE VMA=VOLD IF(V.GT.VMI)THEN VMI=V ENDIF ENDIF V=(VMI+VMA)/2. ENDIF C**** in this case iteration is done with TRNSYS ENDIF C--------------------------------------------------------------- I=CURRENT(V) 1000 CONTINUE P=I*V C**** updating and storing voltage values for the C**** convergence promotor C**** voltage is only stored at timesteps with incident radiation IF(SUN.NE.0.)THEN VOLD=V ENDIF S(ISTORE)=VOLD S(ISTORE+1)=VMI S(ISTORE+2)=VMA S(ISTORE+3)=DUMMY S(ISTORE+4)=MEMO S(ISTORE+5)=N C**** SET OUTPUTS ************************************ OUT(1)=V OUT(2)=I OUT(3)=P OUT(4)=PMAX IF(PMAX.NE.0.)THEN UTIL=P/PMAX ELSE UTIL=0. ENDIF OUT(5)=UTIL OUT(6)=VMP OUT(7)=IMP OUT(8)=VOC OUT(9)=ISC OUT(10)=FF RETURN 1 END C**************************************************************** FUNCTION CURRENT(V) C**** function computes the current for a given voltage C**** the implicit equation is solved using Newton's C**** method IMPLICIT NONE COMMON /ARRAY/ GAM,TC,Q_BZ,IL,IO,RS,IMP REAL IOLD,INEW,F,FPRIME,IMP,I,CURRENT,V REAL GAM,TC,Q_BZ,IL,IO,RS IOLD=0.0 C**** first guess is max.power point current INEW=IMP DO WHILE (ABS(INEW-IOLD).GT.1.E-3) IOLD=INEW F=IL-IOLD-IO*(EXP(Q_BZ*(V+IOLD*RS)/(GAM*TC))-1.) FPRIME=-1.-IO*RS*Q_BZ/GAM/TC*EXP(Q_BZ*(V+IOLD*RS) > /GAM/TC) INEW=IOLD-(F/FPRIME) END DO I=INEW CURRENT=I END C**************************************************************** SUBROUTINE SERIES(RS,Q_BZ) C**** determination of series resistance using bisection method IMPLICIT NONE COMMON /PARAM/ SERCELL,TCR,IMR,VMR,ISCR,VOCR,MVOC,MISC,EG REAL SERCELL,TCR,IMR,VMR,ISCR,VOCR,MVOC,MISC,EG,Q_BZ REAL ANEW,ALOW,AUP,FNEW,FLOW,FUP,GAMNEW,GAMLOW,GAMUP REAL IONEW,IOLOW,IOUP,RS,RSNEW,RSLOW,RSUP C**** parameters at the upper limit of the convergence interval RSUP=((SERCELL*TCR*LOG(1.-IMR/ISCR)/Q_BZ)+VOCR-VMR)/IMR AUP=1. GAMUP=SERCELL IOUP=ISCR*EXP(-Q_BZ*VOCR/(GAMUP*TCR)) C**** parameters at the lower limit of the convergence interval RSLOW=0.0 GAMLOW=Q_BZ*(VMR-VOCR)/(TCR*LOG(1.-IMR/ISCR)) ALOW=GAMLOW/SERCELL IOLOW=ISCR*EXP(-Q_BZ*VOCR/(GAMLOW*TCR)) DO WHILE ((ABS(RSUP-RSLOW)).GT.0.0005) RSNEW=0.5*(RSUP+RSLOW) GAMNEW=Q_BZ*(VMR-VOCR+IMR*RSNEW)/(TCR*LOG(1.-IMR/ISCR)) ANEW=GAMNEW/SERCELL IONEW=ISCR*EXP(-Q_BZ*VOCR/(GAMNEW*TCR)) FUP=-MVOC+(GAMUP/Q_BZ)*(LOG(1.+ISCR/IOUP)+(TCR/(ISCR+ > IOUP))*(MISC-ISCR*((Q_BZ*EG/(AUP*TCR*TCR))+3./TCR))) FLOW=-MVOC+(GAMLOW/Q_BZ)*(LOG(1.+ISCR/IOLOW)+(TCR/(ISCR+ > IOLOW))*(MISC-ISCR*((Q_BZ*EG/(ALOW*TCR*TCR))+3./TCR))) FNEW=-MVOC+(GAMNEW/Q_BZ)*(LOG(1.+ISCR/IONEW)+(TCR/(ISCR+ > IONEW))*(MISC-ISCR*((Q_BZ*EG/(ANEW*TCR*TCR))+3./TCR))) IF((FLOW*FNEW).LT.0.0) RSUP=RSNEW IF((FLOW*FNEW).GT.0.0) RSLOW=RSNEW GAMUP=Q_BZ*(VMR-VOCR+IMR*RSUP)/(TCR*LOG(1.-IMR/ISCR)) AUP=GAMUP/SERCELL IOUP=ISCR*EXP(-Q_BZ*VOCR/(GAMUP*TCR)) GAMLOW=Q_BZ*(VMR-VOCR+IMR*RSLOW)/(TCR*LOG(1.-IMR/ISCR)) ALOW=GAMLOW/SERCELL IOLOW=ISCR*EXP(-Q_BZ*VOCR/(GAMLOW*TCR)) END DO RS=RSNEW END