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
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
|
*****************************************************************
* *
* Procedure : Steady.c2m *
* Purpose : reactor steady state computation including *
* neutron physics, thermalhydraulics, local xenon *
* effects and critical boron search feature *
* Author : V. Salino *
* *
* CALL : *
* Power Fmap Matex Thermo MicroF Burn := Steady *
* Fmap Matex Thermo MicroF Burn ReflXS Track :: *
* <<CritCB>> <<CBinit>> *
* >>CB<< >>keff<< ; *
* *
* Parameters: input *
* CritCB set to True to search critical boron concentration *
* computation, False otherwise (logical variable). *
* CBinit initialization for boron concentration (ppm) when *
* CritCB.EQ.True, otherwise boron concentration to *
* be used. *
* Prel relative power=requested power/nominal power. *
* InletTmod inlet coolant temperature (Celsius). *
* *
* Parameters: output *
* CB critical boron concentration if CritCB.EQ.True, *
* can be negative in which case CB is extrapolated *
* and should not considered physical but fictive. *
* keff k-effective, differs from 1 only if CB is *
* negative or if CritCB.EQ.False. *
* *
*****************************************************************
PARAMETER Power Flux Thermo MicroF Burn Fmap Matex SapUOX SapMOX
MacroRefl Track ::
::: LINKED_LIST Power Flux Thermo MicroF Burn Fmap Matex SapUOX
SapMOX MacroRefl Track ; ;
MODULE SCR: NCR: MACINI: RESINI: TRIVAA: FLUD: FLPOW: THM: EVO: SIM:
GREP: DELETE: ABORT: END: ;
*--
* Local objects and procedure
*--
LINKED_LIST FluxAv MacroTot MacroF System PowerIter
OldPower
;
PROCEDURE ComparePdist ;
*--
* Procedure arguments
*--
LOGICAL CritCB evo pred ;
INTEGER BUindex ;
REAL CB ;
:: >>CritCB<< >>CB<< >>evo<< >>pred<< ;
*--
* Local variables
*--
REAL CBmax := 2000.0 ;
REAL CBmin := 0.0 ;
REAL ErrorRho ErrorCB ErrorTFuel ErrorDCool ErrorTCool ErrorPdis ;
REAL EpsRho EpsCB EpsTFuel EpsDCool EpsPdis :=
0.5 0.05 0.5 5.0E-5 5.0E-4 ;
INTEGER Iter := 1 ;
REAL CBinterp CBp1 CBm1 keff DeltaRho Rho Rhom1 ;
LOGICAL CONV := $False_L ;
REAL DiffBorWorth ;
*--
* Thermalhydraulics parameters
*--
REAL dx := 21.5 ;
REAL Tot_tub := 6.6E-03 2.0 ** $Pi_R * 25.0 * ;
REAL Tot_pin := 4.7E-03 2.0 ** $Pi_R * 264.0 * ;
REAL asssect := dx dx * 1.E-04 * Tot_tub - Tot_pin - ;
REAL coresect := 157.0 asssect * ;
REAL sass := dx dx * 1.E-04 * ;
*--
* Initialized with a THM computation using a flat power distribution
*--
REAL Ptot := 2750.0 ;
REAL BundPow := Ptot 157.0 / 29.0 / 1.0E+3 * (* MW to kW *) ;
Fmap := RESINI: Fmap ::
BUNDLE-POW SAME <<BundPow>> ;
*--
* Iteration loop with feedback, flux, xenon and simplified
* thermalhydraulics
*--
IF pred THEN
FluxAv := Flux ;
ENDIF ;
*--
REPEAT
*--
* Determine CB to use for interpolation
*--
IF CB CBmin < THEN
EVALUATE CBinterp := CBmin ;
ELSEIF CB CBmax > THEN
EVALUATE CBinterp := CBmax ;
ELSE
EVALUATE CBinterp := CB ;
ENDIF ;
*--
* Thermalhydraulics computation, xenon saturation
*--
IF Iter 1 = THEN
Fmap := SIM: Fmap ::
EDIT 1
SET-PARAM 'C-BORE' <<CBinterp>> ;
ELSE
System := DELETE: System ;
Thermo Fmap := THM: Thermo Fmap :: EDIT 0
SET-PARAM 'C-BORE' <<CBinterp>> ;
! Evolution jusqu'a l'equilibre seulement si demande
! Sinon on va rechercher la CB crit sans etre a l'equilibre Xe
IF evo THEN
ECHO "bugg1" ;
Burn MicroF := EVO: Burn MicroF PowerIter ::
EDIT 0 FLUX_POW DEPL 100. DAY KEEP ;
ECHO "bugg2" ;
Burn MicroF := EVO: Burn MicroF PowerIter ::
EDIT 0 FLUX_POW SAVE 100. DAY KEEP ;
ENDIF ;
ENDIF ;
*--
* Saphyb Interpolation
* Historic effects are neglected for simplicity's sake
*--
MicroF := SCR: MicroF SapUOX SapMOX Fmap ::
EDIT 0 RES
MICRO LINEAR
TABLE SapUOX 'BURN'
MIX 4
INST-BURN
SET TF 526.85 ! Celsius
SET TCA 326.85 ! Celsius
SET DCA 0.659 ! g/cm3
SET ppmB 500.0 ! ppm
ADD 'TF' 526.85 MAP
REF 'BURN' SAMEASREF ENDREF
ADD 'TCA' 326.85 MAP
REF 'BURN' SAMEASREF ENDREF
ADD 'DCA' 0.659 MAP
REF 'BURN' SAMEASREF ENDREF
ADD 'ppmB' 500.0 MAP
REF 'BURN' SAMEASREF ENDREF
MICRO ALL
* Recherche poids Xenon : commenter/decommenter
! 'Xe135' *
! 'I135' 0.
'B10' *
'B11' *
* Temporairement, Sm/Pm sont remis a leurs valeurs de Saphyb
* (equivalent a ne pas les particulariser). Apres premier
* appel a Steady, on pourrait enlever cela pour les laisser
* evoluer librement.
! 'Nd147' 0.
! 'Pm147' 0.
! 'Pm148' 0.
! 'Pm148m' 0.
! 'Pm149' 0.
! 'Sm149' 0.
'*MAC*RES' 1.
ENDMIX
TABLE SapMOX 'BURN'
MIX 5
INST-BURN
SET TF 526.85 ! Celsius
SET TCA 326.85 ! Celsius
SET DCA 0.659 ! g/cm3
SET ppmB 500.0 ! ppm
ADD 'TF' 526.85 MAP
REF 'BURN' SAMEASREF ENDREF
ADD 'TCA' 326.85 MAP
REF 'BURN' SAMEASREF ENDREF
ADD 'DCA' 0.659 MAP
REF 'BURN' SAMEASREF ENDREF
ADD 'ppmB' 500.0 MAP
REF 'BURN' SAMEASREF ENDREF
MICRO ALL
* Recherche poids Xenon : commenter/decommenter
! 'Xe135' *
! 'I135' 0.
'B10' *
'B11' *
* Temporairement, Sm/Pm sont remis a leurs valeurs de Saphyb
* (equivalent a ne pas les particulariser). Apres premier
* appel a Steady, on pourrait enlever cela pour les laisser
* evoluer librement.
! 'Nd147' 0.
! 'Pm147' 0.
! 'Pm148' 0.
! 'Pm148m' 0.
! 'Pm149' 0.
! 'Sm149' 0.
'*MAC*RES' 1.
ENDMIX
CHAIN
B10 STABLE
B11 STABLE
I135 NG 0.
Xe135 NG 0. FROM DECAY 1.0E+00 I135
Nd147 STABLE
Pm147 STABLE
Pm148 STABLE
Pm148m STABLE
Pm149 STABLE
Sm149 STABLE
*--------- With Samarium -----------------------------------
Nd147 NG 0.
Pm147 NG 0. FROM DECAY 1.0E+00 Nd147
Pm148 NG 0. FROM NG 5.3E-01 Pm147
Pm148m NG 0. FROM NG 4.7E-01 Pm147
Pm149 NG 0. FROM NG 1.0E+00 Pm148
NG 1.0E+00 Pm148m
Sm149 NG 0. FROM DECAY 1.0E+00 Pm149
MACT NFTOT 0.
ENDCHAIN
;
MacroF := MicroF :: STEP UP "MACROLIB" ;
*--
* Flux and power computation
*--
MacroTot Matex := MACINI: Matex MacroRefl MacroF ;
MacroF := DELETE: MacroF ;
System := TRIVAA: MacroTot Track :: EDIT 0 ;
MacroTot := DELETE: MacroTot ;
IF Iter 1 = THEN
Flux := FLUD: Flux System Track :: EDIT 0
ADI 4 EXTE 1.0E-5 ACCE 5 3 ;
ELSE
Flux := FLUD: Flux System Track :: EDIT 0
RELAX 0.15 ;
ENDIF ;
GREP: Flux :: GETVAL 'K-EFFECTIVE ' 1 >>keff<< ;
EVALUATE Rho := 1. 1. keff / - 1.0E5 * ;
IF Iter 1 > THEN
PowerIter := DELETE: PowerIter ;
ENDIF ;
PowerIter Fmap := FLPOW: Fmap Flux Track Matex ::
EDIT 0 PTOT <<Ptot>> PRINT DISTR POWER ;
*--
* Compute convergence errors (absolute values) on :
* - reactivity (pcm)
* - boron concentration (ppm)
* - 3D power distribution (relative discrepancy, maximum)
* - fuel temperature distribution (K, maximum)
* - moderator density distribution (g/cm3, maximum)
*--
IF Iter 1 > THEN
EVALUATE ErrorRho := Rho Rhom1 - ABS ;
EVALUATE ErrorCB := CB CBm1 - ABS ;
Fmap OldPower PowerIter := ComparePdist :: >>ErrorPdis<< ;
OldPower := DELETE: OldPower ;
GREP: Thermo :: GETVAL 'ERROR-T-FUEL' 1 >>ErrorTFuel<< ;
GREP: Thermo :: GETVAL 'ERROR-D-COOL' 1 >>ErrorDCool<< ;
ENDIF ;
OldPower := PowerIter ;
*--
* Reestimate critical CB. 3rd iteration is waited, in order to have a
* not-too-false k-effective, in particular regarding thermalhydraulics
* and neutron physics convergence. Otherwise predicted CB would be too
* far from the truth.
*--
IF CritCB Iter 2 > * THEN
* To evaluate DiffBorWorth, we have 3 other possibilities :
* - Evaluate properly DiffBorWorth (internship ?), may be costly
* - Use the 2 previous computations to find (1 line with 2 points)
* (may be instable during the first iterations but effective
* near the end)
* - Use module FIND0 (may be long)
* Plot reactivity=f(CB), see which is going to be good
EVALUATE DiffBorWorth := -7.0 ; ! pcm/ppm
EVALUATE CBp1 := CBinterp Rho DiffBorWorth / - ;
ENDIF ;
*--
* Prints
*--
ECHO "#---" ;
ECHO "Iter" Iter ;
ECHO "keff=" keff "CB=" CB "Rho=" Rho ;
ECHO "CBinterp=" CBinterp ;
IF CritCB Iter 2 > * THEN
ECHO "DiffBorWorth=" DiffBorWorth "CBp1" CBp1 ;
ENDIF ;
ECHO "---" ;
*--
* Check convergence
*--
IF Iter 3 > THEN
ECHO "Discrepancies between the previous iteration and current"
"iteration" ;
ECHO " Rho (pcm) CB (ppm) TFuel (K) DCool (K) "
"Pdis (relat)" ;
ECHO ErrorRho ErrorCB ErrorTFuel ErrorDCool ErrorPdis ;
ECHO "" ;
ECHO "Convergence criterion (epsilon)" ;
ECHO " Rho (pcm) CB (ppm) TFuel (K) DCool (K) "
"Pdis (relat)" ;
ECHO EpsRho EpsCB EpsTFuel EpsDCool EpsPdis ;
ECHO "" ;
ECHO "Convergence reached ?" ;
ECHO " Rho (pcm) CB (ppm) TFuel (K) DCool (K) "
"Pdis (relat)" ;
ECHO "" ErrorRho EpsRho < " " ErrorCB EpsCB < " "
ErrorTFuel EpsTFuel < " " ErrorDCool EpsDCool <
" " ErrorPdis EpsPdis < ;
EVALUATE CONV := ErrorRho EpsRho <
ErrorCB EpsCB < *
ErrorTFuel EpsTFuel < *
ErrorDCool EpsDCool < *
ErrorPdis EpsPdis < * ;
ENDIF ;
*--
* Prepare next iteration
*--
IF CONV NOT THEN
EVALUATE CBm1 Rhom1 := CB Rho ;
IF CritCB Iter 2 > * THEN
EVALUATE CB := CBp1 ;
ENDIF ;
EVALUATE Iter := Iter 1 + ;
ENDIF ;
*--
IF Iter 100 > THEN
ECHO "Steady.c2m: maximum iteration reached (50)." ;
ABORT: ;
ENDIF ;
UNTIL Iter 4 = ;
*--
* Predicteur-correcteur
*--
IF pred THEN
PowerIter := DELETE: PowerIter ;
FluxAv := FLUD: FluxAv System Track :: EDIT 0
RELAX 0.5 ;
PowerIter Fmap := FLPOW: Fmap FluxAv Track Matex ::
EDIT 0 PTOT <<Ptot>> PRINT DISTR POWER ;
FluxAv := DELETE: FluxAv ;
ENDIF ;
*--
* Print Xe and I average fission yields. This average depends
* on burnup, thermalhydraulics and xenon conditions.
*--
REAL YieldI YieldXe ;
GREP: MicroF :: STEP UP 'DEPL-CHAIN'
GETVAL 'FISSIONYIELD' 1 >>YieldI<<
GETVAL 'FISSIONYIELD' 2 >>YieldXe<< ;
ECHO "Mean fission yield for I =" YieldI ;
ECHO "Mean fission yield for Xe=" YieldXe ;
Power := PowerIter ;
PowerIter := DELETE: PowerIter ;
System := DELETE: System ;
:: <<CB>> ;
:: <<keff>> ;
END: ;
|
