diff options
Diffstat (limited to 'Donjon/src/THMCDI.f')
| -rw-r--r-- | Donjon/src/THMCDI.f | 209 |
1 files changed, 209 insertions, 0 deletions
diff --git a/Donjon/src/THMCDI.f b/Donjon/src/THMCDI.f new file mode 100644 index 0000000..89f00af --- /dev/null +++ b/Donjon/src/THMCDI.f @@ -0,0 +1,209 @@ +*DECK THMCDI + FUNCTION THMCDI(T2K,T1K,BURN,POROS,FRACPU,ICONDF,NCONDF,KCONDF, + > UCONDF,IFRCDI) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Compute the thermal conductivity integral of UOX or MOX fuel. +* +*Copyright: +* Copyright (C) 2012 Ecole Polytechnique de Montreal. +* +*Author(s): +* A. Hebert, V. Salino +* +*Parameters: input +* T2K final temperature in Kelvin. +* T1K initial temperature in Kelvin. +* BURN fuel burnup in MWday/tonne. +* POROS fuel porosity. +* FRACPU plutonium mass fraction in fuel. +* ICONDF fuel conductivity flag (0=Stora-Chenebault or COMETHE/ +* 1=user-provided polynomial + inverse term). +* NCONDF degree of user-provided fuel conductivity polynomial. +* KCONDF polynomial coefficients for fuel conductivity in W/m/K^(k+1) +* (except for the two last coefficients which belongs to the +* inverse term). +* UCONDF required unit of temperature in polynomial for fuel +* conductivity (KELVIN or CELSIUS). +* IFRCDI flag indicating if average approximation is forced during +* fuel conductivity evaluation (0=default/1=average +* approximation forced). +* +*Parameters: output +* THMCDI thermal conductivity integral in Watt/m/K. +* +*Reference: +* A. Poncot, "Assimilation de donnees pour la dynamique du xenon dans +* les coeurs de centrale nucleaire", Ph.D Thesis, Universite de +* Toulouse, France, 2008. +* +*----------------------------------------------------------------------- +* + IMPLICIT NONE +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER ICONDF,NCONDF,IFRCDI + REAL T1K,T2K,BURN,POROS,FRACPU,KCONDF(NCONDF+3),THMCDI + CHARACTER UCONDF*12 +*---- +* LOCAL VARIABLES +* NPAS number of rectangles in the quadrature +* DT rectangle width +* T2T1 temperature difference +* DTMIN cutoff criterion for selecting the approximation +* FPI burnup correcting factor +* CIRRA burnup correction constant +* HV* coefficients of the Stora-Chenebault correlation +* HK* coefficients of the Comethe correlation +*---- + INTEGER NPAS,I,K + REAL T1,T2,DT,TM,DTMIN,T2T1,FPI,TT,TMK,TEMP,FTP,CINT,FP,TTK + REAL HV1, HV2, HV3 + REAL HK1, HK2, HK4, HK5 + REAL ZKELV,CIRRA +* + PARAMETER ( ZKELV=273.15 ) + PARAMETER ( HV1= 1.3324E-08 , HV2 = -4.3554E-05 , + & HV3 = 5.8915E-02 ) + PARAMETER ( HK1= 40.05 , HK2 = 129.4 , HK4 = 0.8 , + & HK5 = 0.6416E-12 ) + PARAMETER ( CIRRA= 0.124E-02 ) +* + REAL A + DATA NPAS /10/ + DATA DTMIN /10./ +* + IF(MIN(T1K,T2K).LE.0.0) THEN + CALL XABORT('@THMCDI: NEGATIVE TEMPERATURE.') + ENDIF + T1=T1K-ZKELV + T2=T2K-ZKELV +* + T2T1 = T2-T1 + DT = T2T1/NPAS + TM = (T1+T2)/2.0 + IF(ICONDF.EQ.1) THEN +* User-given conductivity, as a function of temperature + IF((ABS(T2T1).LT.DTMIN).OR.(IFRCDI.EQ.1)) THEN +* Use the average value approximation + THMCDI=0.0 + IF(UCONDF.EQ.'KELVIN') THEN + TMK = TM + ZKELV + DO K=1,NCONDF+1 + THMCDI=THMCDI + KCONDF(K)*TMK**(K-1) + ENDDO + THMCDI=THMCDI + KCONDF(NCONDF+2)/(TMK-KCONDF(NCONDF+3)) + ELSE + DO K=1,NCONDF+1 + THMCDI=THMCDI + KCONDF(K)*TM**(K-1) + ENDDO + THMCDI=THMCDI + KCONDF(NCONDF+2)/(TM-KCONDF(NCONDF+3)) + ENDIF + ELSE +* Use the rectangle quadrature approximation + TT=T1-DT*0.5 + CINT=0. + DO I=1,NPAS + TT=TT+DT + IF(UCONDF.EQ.'KELVIN') THEN + TTK = TT + ZKELV + DO K=1,NCONDF+1 + CINT=CINT + KCONDF(K)*TTK**(K-1) + ENDDO + CINT=CINT + KCONDF(NCONDF+2)/(TTK-KCONDF(NCONDF+3)) + ELSE + DO K=1,NCONDF+1 + CINT=CINT + KCONDF(K)*TT**(K-1) + ENDDO + CINT=CINT + KCONDF(NCONDF+2)/(TT-KCONDF(NCONDF+3)) + ENDIF + ENDDO + THMCDI=CINT/NPAS + ENDIF + ELSE IF(FRACPU.GT.0.) THEN +* Use the Comethe correlation for MOX fuel + FPI=CIRRA*BURN + IF((ABS(T2T1).LT.DTMIN).OR.(IFRCDI.EQ.1)) THEN +* Use the average value approximation + IF(TM.GT.1000.0) THEN + A=2.0 + ELSE + A=2.58-0.58E-03*TM + ENDIF + FP=(1.0-A*POROS)/(1.0-A*0.05) + TMK = TM + ZKELV + TEMP = HK2 + (1.0 + HK4*FRACPU*1.E-02) * TMK + FTP = FP * (HK1/TEMP + HK5*TMK*TMK*TMK) *100.0 + IF(TM.EQ.0.) THEN + THMCDI=FTP + ELSE + THMCDI=1.0/(1.0/FTP+FPI/TM) + ENDIF + ELSE +* Use the rectangle quadrature approximation + TT=T1-DT*0.5 + CINT=0. + DO I=1,NPAS + TT=TT+DT + IF(TT.GT.1000.0) THEN + A=2.0 + ELSE + A=2.58-0.58E-03*TT + ENDIF + FP=(1.0-A*POROS)/(1.0-A*0.05) + TTK = TT + ZKELV + TEMP = HK2 + (1.0 + HK4*FRACPU*1.E-02) * TTK + FTP = FP * (HK1/TEMP + HK5*TTK*TTK*TTK) *100.0 + IF(TT.EQ.0.0) THEN + CINT=CINT+FTP + ELSE + CINT=CINT+1.0/(1.0/FTP+FPI/TT) + ENDIF + ENDDO + THMCDI=CINT/NPAS + ENDIF + ELSE +* Use the Stora-Chenebault correlation for UOX fuel +* (also called the "HGAP Variable 88" correlation) + FPI=CIRRA*BURN + IF((ABS(T2T1).LT.DTMIN).OR.(IFRCDI.EQ.1)) THEN +* Use the average value approximation + IF(TM.GT.1000.) THEN + A=2.0 + ELSE + A=2.58-0.58E-03*TM + ENDIF + FP=(1.0-A*POROS)/(1.0-A*0.034) + FTP=FP*(HV3+HV2*TM+HV1*TM*TM)*100.0 + IF(TM.EQ.0.0) THEN + THMCDI=FP*HV3*100.0 + ELSE + THMCDI=1.0/(1.0/FTP+FPI/TM) + ENDIF + ELSE +* Use the rectangle quadrature approximation + TT=T1-DT*0.5 + CINT=0. + DO I=1,NPAS + TT=TT+DT + IF(TT.GT.1000.) THEN + A=2.0 + ELSE + A=2.58-0.58E-03*TT + ENDIF + FP=(1.0-A*POROS)/(1.0-A*0.034) + FTP=FP*(HV3+HV2*TT+HV1*TT*TT)*100.0 + IF(TT.EQ.0.0) THEN + CINT=CINT+FP*HV3*100.0 + ELSE + CINT=CINT+1.0/(1.0/FTP+FPI/TT) + ENDIF + ENDDO + THMCDI=CINT*DT/T2T1 + ENDIF + ENDIF + RETURN + END |