Initial commit from Polytechnique Montreal

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