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
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
|
*DECK THMROD
SUBROUTINE THMROD(IMPX,NFD,NDTOT,MAXIT1,MAXIT2,ERMAXT,DTINV,
1 RAD,XX0,XX1,QFUEL,FRO,TSURF,POWLIN,BURN,POROS,FRACPU,ICONDF,
2 NCONDF,KCONDF,UCONDF,ICONDC,NCONDC,KCONDC,UCONDC,IHGAP,KHGAP,
3 IFRCDI,TC1,XX2,XX3,ZF)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Solution of the discretized thermal conduction equations in a single
* fuel rod in an axial slice of the fuel channel.
*
*Copyright:
* Copyright (C) 2018 Ecole Polytechnique de Montreal.
*
*Author(s):
* A. Hebert
*
*Parameters: input
* IMPX print parameter.
* NFD number of fuel discretized points in the cladded fuel rod.
* The last point of the discretization (i=NFD) is taken at the
* surface of the fuel pellet.
* NDTOT total number of discretized points in the cladded fuel rod
* with the radial zones in the cladding. The points which are
* located at i=NFD+1 and i=NDTOT are respectively taken at the
* inner surface of clad and in the center of external clad ring.
* MAXIT1 maximum number of conduction iterations.
* MAXIT2 maximum number of center-pellet iterations.
* ERMAXT convergence criterion.
* DTINV inverse time step. Equal to 1/DT in transient cases. Equal to
* 0 in steady-state cases.
* RAD fuel and clad radii (m).
* XX0 temperatures at time n-1 (K).
* XX1 estimate of the temperatures at time n (K).
* QFUEL volumic power in fuel at time n (W/m^3).
* FRO radial power form factors. All components are set to 1.0 for
* a constant power source in fuel.
* TSURF estimate of the external clad surface temperature at
* time n (K).
* POWLIN estimate of the linear power at time n (W/m).
* BURN fuel burnup in MWday/tonne.
* POROS fuel porosity.
* FRACPU plutonium percent fraction.
* 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).
* ICONDC clad conductivity flag (0=default/1=user-provided
* polynomial).
* NCONDC degree of user-provided clad conductivity polynomial.
* KCONDC polynomial coefficients for clad conductivity in W/m/K^(k+1).
* UCONDC required unit of temperature in polynomial for clad
* conductivity (KELVIN or CELSIUS).
* IHGAP flag indicating HGAP chosen (0=default/1=user-provided).
* KHGAP fixed user-provided HGAP value in W/m^2/K.
* IFRCDI flag indicating if average approximation is forced during
* fuel conductivity evaluation (0=default/1=average
* approximation forced).
*
*Parameters: output
* TC1 estimate of center-pellet temperature at time n (K).
* XX1 estimate of the temperatures at time n (K).
* XX2 first component of temperatures at time n (K).
* XX3 second component of temperatures at time n. The actual
* temperatures are given as XX2(:)+TSURF*XX3(:) where TSURF
* is the temperature of the external clad surface.
* ZF components of the linear power transmitted from clad to fluid.
* The linear power (W/m) is given as 2*PI*(ZF(1)-TSURF*ZF(2)).
*
*-----------------------------------------------------------------------
*
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER IMPX,NFD,NDTOT,MAXIT1,MAXIT2,ICONDF,NCONDF,ICONDC,NCONDC,
1 IHGAP,IFRCDI
REAL ERMAXT,DTINV,RAD(NDTOT),XX0(NDTOT),XX1(NDTOT),QFUEL,
1 FRO(NFD-1),TSURF,POWLIN,BURN,POROS,FRACPU,KCONDF(NCONDF+3),
2 KCONDC(NCONDC+1),KHGAP,TC1,XX2(NDTOT),XX3(NDTOT),ZF(2)
CHARACTER UCONDF*12,UCONDC*12
*----
* LOCAL VARIABLES
*----
REAL COEF(3)
*----
* ALLOCATABLE ARRAYS
*----
REAL, ALLOCATABLE, DIMENSION(:) :: DAR,ZK,CONDXA
REAL, ALLOCATABLE, DIMENSION(:,:) :: TRID
*----
* SCRATCH STORAGE ALLOCATION
*----
ALLOCATE(DAR(NDTOT),ZK(NDTOT),CONDXA(NDTOT),TRID(NDTOT,NDTOT+2))
*----
* COMPUTE ARs AND VOLUMES
*----
DAR(:NDTOT)=0.0
ARF=0.5*RAD(NFD)**2 ! at fuel radius
ARCI=0.5*RAD(NFD+1)**2 ! at internal clad radius
ARCE=0.5*RAD(NDTOT)**2 ! at external clad radius
DO I=1,NFD-1
DAR(I)=0.5*(RAD(I+1)**2-RAD(I)**2)
ENDDO
DO I=NFD+2,NDTOT
DAR(I)=0.5*(RAD(I)**2-RAD(I-1)**2)
ENDDO
*----
* COMPUTE THE THERMAL CONDUCTIVITY INTEGRALS AT TIME n-1
*----
ZK(:NDTOT)=0.0
DO I=1,NFD-1
ZK(I)=THMCDI(XX0(I),XX0(I+1),BURN,POROS,FRACPU,ICONDF,NCONDF,
> KCONDF,UCONDF,IFRCDI)
ENDDO
DO I=NFD+1,NDTOT-1
ZK(I)=THMGDI(XX0(I),XX0(I+1),ICONDC,NCONDC,KCONDC,UCONDC)
ENDDO
ZK(NDTOT)=THMGDI(XX0(NDTOT),TSURF,ICONDC,NCONDC,KCONDC,UCONDC)
*----
* COMPUTE CONDXA
*----
CONDXA(:NDTOT)=0.0
COEF(1)=0.0
COEF(2)=0.0
COEF(3)=0.0
DO I=1,NDTOT
IF(I.LE.NFD-1) THEN
CONDXA(I)=THMCCD(XX0(I),POROS,FRACPU)*DTINV*XX0(I)*DAR(I)
ELSE IF(I.GE.NFD+2) THEN
CONDXA(I)=THMGCD(XX0(I))*DTINV*XX0(I)*DAR(I)
ENDIF
ENDDO
*----
* ITERATIVE PROCEDURE
*----
ITERT=0
10 ITERT=ITERT+1
IF(ITERT.GT.MAXIT1) CALL XABORT('THMROD: CONVERGENCE FAILURE(1).')
*----
* COMPUTE THE THERMAL CONDUCTIVITY INTEGRALS AT TIME n
*----
ZK(:NDTOT)=0.0
DO I=1,NFD-1
ZK(I)=THMCDI(XX1(I),XX1(I+1),BURN,POROS,FRACPU,ICONDF,NCONDF,
> KCONDF,UCONDF,IFRCDI)
ENDDO
DO I=NFD+1,NDTOT-1
ZK(I)=THMGDI(XX1(I),XX1(I+1),ICONDC,NCONDC,KCONDC,UCONDC)
ENDDO
ZK(NDTOT)=THMGDI(XX1(NDTOT),TSURF,ICONDC,NCONDC,KCONDC,UCONDC)
IF(IHGAP.EQ.0) THEN
CALL THMGAP(POWLIN,BURN,HGAP)
ELSE IF(IHGAP.EQ.1) THEN
HGAP=KHGAP
ENDIF
*----
* BUILD THE TRIDIAGONAL SYSTEM
*----
TRID(:NDTOT,:NDTOT+2)=0.0
COEF(1)=0.0
COEF(2)=0.0
COEF(3)=0.0
DO I=1,NDTOT
TRID(I,NDTOT+1)=CONDXA(I)
IF(I.LE.NFD-2) THEN
ARI=0.5*RAD(I+1)**2
COEF(3)=4.0*ARI*ZK(I)/(DAR(I)+DAR(I+1))
TRID(I,NDTOT+1)=TRID(I,NDTOT+1)+QFUEL*FRO(I)*DAR(I)
ELSE IF(I.EQ.NFD-1) THEN
ARI=0.5*RAD(I+1)**2
COEF(3)=4.0*ARI*ZK(I)/DAR(I)
TRID(I,NDTOT+1)=TRID(I,NDTOT+1)+QFUEL*FRO(I)*DAR(I)
ELSE IF(I.EQ.NFD) THEN
RAVG=2.0*RAD(NFD)*RAD(NFD+1)/(RAD(NFD)+RAD(NFD+1))
COEF(3)=RAVG*HGAP
ELSE IF(I.EQ.NFD+1) THEN
COEF(3)=4.0*ARCI*ZK(I)/DAR(I+1)
ELSE IF(I.LE.NDTOT-1) THEN
ARI=0.5*RAD(I)**2
COEF(3)=4.0*ARI*ZK(I)/(DAR(I)+DAR(I+1))
ELSE IF(I.EQ.NDTOT) THEN
COEF(3)=4.0*ARCE*ZK(I)/DAR(I)
TRID(I,NDTOT+2)=TRID(I,NDTOT+2)+COEF(3)
ENDIF
COEF(2)=COEF(1)+COEF(3)
IF(I.GT.1) TRID(I,I-1)=-COEF(1)
IF(I.LE.NFD-1) THEN
TRID(I,I)=THMCCD(XX1(I),POROS,FRACPU)*DTINV*DAR(I)
ELSE IF(I.GE.NFD+2) THEN
TRID(I,I)=THMGCD(XX1(I))*DTINV*DAR(I)
ENDIF
TRID(I,I)=TRID(I,I)+COEF(2)
IF(I.LT.NDTOT) THEN
TRID(I,I+1)=-COEF(3)
COEF(1)=COEF(3)
ENDIF
ENDDO
ZWORK=COEF(3)
*----
* SOLVE LINEAR SYSTEM
*----
CALL ALSB(NDTOT,2,TRID,IER,NDTOT)
IF(IER.NE.0) CALL XABORT('THMROD: SINGULAR MATRIX')
*----
* SET TEMPERATURE AT TIME n
*----
ERR=0.0
IMAX=0
DO I=1,NDTOT
TNEW=TRID(I,NDTOT+1)+TSURF*TRID(I,NDTOT+2)
IF(ABS(XX1(I)-TNEW).GT.ERR) THEN
ERR=ABS(XX1(I)-TNEW)
IMAX=I
ENDIF
IF(ITERT.LE.20) THEN
XX1(I)=TNEW
ELSE
* perform under-relaxation
XX1(I)=0.5*(TNEW+XX1(I))
ENDIF
ENDDO
ZF(1)=ZWORK*TRID(NDTOT,NDTOT+1)
ZF(2)=ZWORK*(1.0-TRID(NDTOT,NDTOT+2))
IF(IMPX.GT.4) WRITE(6,100) ITERT,ERR,ERMAXT,IMAX
IF((ERR.LT.ERMAXT).AND.(ITERT.NE.1)) GO TO 20
GO TO 10
*----
* SCRATCH STORAGE DEALLOCATION
*----
20 DO I=1,NDTOT
XX2(I)=TRID(I,NDTOT+1)
XX3(I)=TRID(I,NDTOT+2)
ENDDO
DEALLOCATE(TRID,CONDXA,ZK,DAR)
*----
* COMPUTE THE CENTER-PELLET TEMPERATURE.
*----
TC=0.5*(XX1(1)+XX1(2))
ITERC=0
30 ITERC=ITERC+1
IF(ITERC.GT.MAXIT2) CALL XABORT('THMROD: CONVERGENCE FAILURE(2).')
TCOLD=TC
CC1=THMCDI(XX1(1),TC,BURN,POROS,FRACPU,ICONDF,NCONDF,KCONDF,
> UCONDF,IFRCDI)
CC2=THMCDI(TC,XX1(2),BURN,POROS,FRACPU,ICONDF,NCONDF,KCONDF,
> UCONDF,IFRCDI)
TC=(CC1*XX1(1)+CC2*XX1(2))/(CC1+CC2)
IF(ITERT.GT.20) TC=0.5*(TC+TCOLD)
DELTAA=ABS(TC-TCOLD)
IF(IMPX.GT.4) WRITE(6,110) ITERC,DELTAA,ERMAXT
IF((DELTAA.LT.ERMAXT).AND.(ITERC.NE.1)) GO TO 40
GO TO 30
40 TC1=2.0*XX1(1)-TC
RETURN
100 FORMAT(/15H THMROD: ITERT=,I5,1P,7H ERROR=,E12.4,5H EPS=,E12.4,
> 5H POS=,I5)
110 FORMAT(/15H THMROD: ITERC=,I5,1P,7H ERROR=,E12.4,5H EPS=,E12.4)
END
|
