summaryrefslogtreecommitdiff
path: root/Donjon/src/THMCDI.f
blob: 89f00af6dc60cbcafdd7b61820b2c53121c80e54 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
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