Contribution to the validation of best estimate plus uncertainties coupled codes for the analysis of NK-TH nuclear transients

ADVERTIMENT . La consulta d’aquesta tesi queda condicionada a l’acceptació de les següents condicions d'ús: La difusió d’aquesta tesi per mitjà del servei TDX (www.tesisenxarxa.net) ha estat autoritzada pels titulars dels drets de propietat intel·lectual únicament per a usos privats emmarcats en activitats d’investigació i docència. No s’autoritza la seva reproducció amb finalitats de lucre ni la seva difusió i posada a disposició des d’un lloc aliè al servei TDX. No s’autoritza la presentació del seu contingut en una finestra o marc aliè a TDX (framing). Aquesta reserva de drets afecta tant al resum de presentació de la tesi com als seus continguts. En la utilització o cita de parts de la tesi és obligat indicar el nom de la persona autora. ADVERTENCIA. La consulta de esta tesis queda condicionada a la aceptación de las siguientes condiciones de uso: La difusión de esta tesis por medio del servicio TDR (www.tesisenred.net) ha sido autorizada por los titulares de los derechos de propiedad intelectual únicamente para usos privados enmarcados en actividades de investigación y docencia. No se autoriza su reproducción con finalidades de lucro ni su difusión y puesta a disposición desde un sitio ajeno al servicio TDR. No se autoriza la presentación de su contenido en una ventana o marco ajeno a TDR (framing). Esta reserva de derechos afecta tanto al resumen de presentación de la tesis como a sus contenidos. En la utilización o cita de partes de la tesis es obligado indicar el nombre de la persona autora. WARNING. On having consulted this thesis you’re accepting the following use conditions: Spreading this thesis by the TDX (www.tesisenxarxa.net) service has been authorized by the titular of the intellectual property rights only for private uses placed in investigation and teaching activities. Reproduction with lucrative aims is not authorized neither its spreading and availability from a site foreign to the TDX service. Introducing its content in a window or frame foreign to the TDX service is not authorized (framing). This rights affect to the presentation summary of the thesis as well as to its contents. In the using or citation of parts of the thesis it’s obliged to indicate the name of the author TECHNICAL UNIVERSITY of CATALONIA Contribution to the validation of best estimate plus uncertainties coupled codes for the analysis of NK-TH nuclear transients by Raimon Pericas A thesis submitted in partial fulfillment for the degree of Doctor of Philosophy in the ETSEIB Physics and Nuclear Engineering Deparment March 2015 Declaration of Authorship I, RAIMON PERICAS, declare that this thesis titled, ‘Contribution to the validation of best estimate plus uncertainties coupled codes for the analysis of NK-TH nuclear transients’ and the work presented in it are my own. I confirm that: This work was done wholly or mainly while in candidature for a research degree at this University. Where any part of this thesis has previously been submitted for a degree or any other qualification at this University or any other institution, this has been clearly stated. Where I have consulted the published work of others, this is always clearly attributed. Where I have quoted from the work of others, the source is always given. With the exception of such quotations, this thesis is entirely my own work. I have acknowledged all main sources of help. Where the thesis is based on work done by myself jointly with others, I have made clear exactly what was done by others and what I have contributed myself. Signed: Date: iii “The block of granite which was an obstacle in the pathway of the weak, became a stepping-stone in the pathway of the strong.” Thomas Carlyle TECHNICAL UNIVERSITY of CATALONIA Abstract ETSEIB Physics and Nuclear Engineering Deparment Doctor of Philosophy by Raimon Pericas The typical conservative nuclear safety margins limit the actual industrial needs to increase the Nuclear Power Plant (NPP) power production. Best Estimate Plus Uncertainty (BEPU) calculations are the most advanced tool in nuclear system codes analysis. This technique is superior when compared to the old conservative methodology, where the safety margins were established by experts, operation hypothesis and conservative assumptions. The BEPU methodology is capable of providing a solution in terms of increasing the nuclear power production without compromising the safety margins. This study presents a comparison between the BEPU methodology and the Conservative Bounding methodology. Within the framework of safety analysis with nuclear system codes, neutron kinetics and thermal hydraulics (NK-TH) calculations are also the most advanced tool, and they are specially indicated for those transients which involves asymmetrical core conditions and return to critically scenarios. Main Steam Line Break (MSLB) in Asc´o (NPP) fits these pre-conditioners and thus is the selected transient for the present report. Some code improvements were needed when validating the used models, those improvements are presented in this study also. Finally, moreover the BEPU analysis with NK-TH coupled codes, the present study also shows a methodology of XS library creation valid for any point of the cycle life of the studied reactor. Acknowledgements I am using this opportunity to express my gratitude to everyone who supported me throughout the course of this Ph.D project. I am thankful for their aspiring guidance, invaluably constructive criticism and friendly advice during the project work. I am sincerely grateful to them for sharing their truthful and illuminating views on a number of issues related to the project. Foremost, I would like to express my sincere gratitude to my advisors Prof. Francesc Revent´os and Prof. Llu´ıs Batet for the continuous support of my Ph.D study and research, for their patience, motivation, enthusiasm, and immense knowledge. Their guidance helped me in all the time of research and writing of this thesis report. I could not have imagined having better advisors and mentors for my Ph.D study. Besides my advisors, I would like to thank the rest of my thesis committee: Prof. Kostadin N. Ivanov, Prof. Maria del Carmen Pretel S´anchez, Prof. Rafael Mir´o Herrero, Prof. Maria N. Avramova and Dra. Patricia Pla Freixa for their encouragement, insightful comments, and hard questions. My sincere thanks also goes to Prof. Kostadin Ivanov for his hospitality, for offering me the summer internship opportunities in his group and leading me working on the neutron kinetics part of the project. I thank my fellow lab mates in Technical University of Catalonia for the stimulating discussions, coffee talks and for all the fun we have had in the last years. Also I thank my friends in Universitat de Vic, Miquel Caballeria, Carme Vern´ıs, Manel Vilar and Josep Ayats for encouraging me through all the project duration. I thank Rafael Mendizabal from Consejo de Seguridad Nuclear for their funding support and Asociaci´on Nuclear Asc´o - Vandell´os II for providing data necessary for this project. I also thank, Dr. Rafael Mir´o from Universitat Polit`ecnica de Val`encia for his help in coding part of the project. Special thanks to Dr. Chris Allison and all Innovative Systems Software community (Zheng, Brian, Judy, Jenny, Peggy, Larry and Dick) for trusting me, for their help, advice, hospitality and support. Last but not the least, I would like to thank all my family for supporting me spiritually throughout this long journey. To all thanks for helping and supporting me through all the project. Raimon. vi Abbreviations ACR Advanced CANDU Reactor AFW Auxiliary Feed Water ATHLET Analysis of Thermal-hydraulics of Leaks and Transients ATWS Anticipated Transient Without SCRAM ANM Analytic Nodal Method BE Best Estimated BEBCC Best Estimated Base Case Calculation BEPU Best Estimated Plus Uncertainty BiCGSTAB Bi-Conjugate Gradient STABilized algorithm BILU3D Blockwise Incomplete LU preconditioner BOC Begin OfCycle BWR Boiling Water Reactor CCFL Counter Current Flow Limitation CHF Critical Heat Flux CIAU Code with capability of Internal Assessment of Uncertainty CMFD Coarse Mesh Finite Difference CNCC Corrective Nodal Coupling Coefficient CRISSUE Critical Issues in Nuclear Reactor Technology CSN Consejo de Seguridad Nuclear DAKOTA Design Analysis Kit for Optimization and Terascale Application DBA Design from Basis Accident DNBR Departure from Nucleate Boiling Ratio DOE Department OfEnergy ECCS Emergency Core Cooling System EFPD Equivalent Full Power Days xv Abbreviations xvi EOC End OfCycle ESBWR Economic Simplified Boiling Water Reactor FMFD Fine Mesh Finite Difference FoM Figure ofMerit FW Feed Water GenPMAXS Generation of the Purdue Macroscopic XS sets GET Grup d’Estudis Termohidr`aulics GMRES Generalized Minimum Residual method HPIS High Pressure Injection System HS Heat Structure IDN Informe de Dis˜neo Nuclear ITF Integral Test Facility LBLOCA Large Break Loss OfCoolant Accident LOCA Loss OfCoolant Accident LPIS Low Pressure Injection System LWR Light Water Reactor MCP Main Coolant Pumps MFW Main Feed Water MOC Mid OfCycle MSIV Main Steam Isolation Valve MSLB Main Steam Line Break NEA Nuclear Energy Agency NEM Nodal Expansion Method NESTLE Nodal Eigenvalue, Steady-state, Transient, Le core Evaluator NK-TH Neutron Kinetics and Thermal Hydraulics NPP Nuclear Power Plant NRC Nuclear Regulatory Commission OECD Organization for Economic Co-operation and Development PARCS Purdue Advanced Reactor Core Simulator PCT Peak Cladding Temperature PDF Probability Density Function PIRT Phenomenon Identification and Ranking Tables PSU Pennsylvania State University Abbreviations xvii PWR Pressurizer Water Reactor RDFMG Reactor Dynamics and Fuel Management Research Group RELAP Reactor Excursion and Leak Analysis Program SETF Separate Effects Test Facility SBLOCA Small Break Loss OfCoolant Accident SCRAM Safety Control Rod Axe Man SNAP Symbolic Nuclear and Analysis Program SP Special Processes TH Thermal Hydraulics TMI Thee Mile Island TPEN Triangle-based Polynomial Expansion Nodal TRAC Transient Reactor Analysis Code TRACE TRACE/RELAP Advanced and Computational Engine VVER Vodo-Vodyanoi Energetichesky Reaktor UAM Uncertainty Analysis in Modelling UPC Universitat Polit`ecnica de Catalunya XS Cross Sections Symbols symbol name unit adistance (m) asubscript: absorption cross section (-) Cm kmnode precursor density (mol/m3) Dgdiffusion coefficient (cm) Di(r, t) concentration of the decay heat precursors in decay heat group i(J/cm3) dsubscript: delayed neutrons (-) Eenergy (J) eff superscript: effective (-) Ffission matrix (-) fsubscript: fission cross section (-) Gneutron energy group Giheavy metal loading in ith region (kg) Gcheavy metal loading in the core (kg) gsubscript: fast/thermal neutron group (-) g→g′subscript: scattering group 1 to 2 cross section (-) Hjenergy of the decay heat precursor concentration in group j(J) henthalpy (J)  Jneutron current (neutron/m2s) Jm+ − gw mnode surface averaged net current (neutron/m2s) Jnumber of decay heat groups (-) keff effective neutron multiplication constant (-) Ksubscript: delayed neutron group (-) xix Symbols xx Lggroups leakage term (neutron/m2s) lsubscript: liquid (-) Mnon-fission matrix (-) msuperscript: mnode (-) Nl i(t) nuclei number density of isotope i(atom/m3) Ppower W (Js−1) ppressure (Pa) psubscript: prompt neutrons (-) q′heat transfer rate per unit volume (W/m2) q heat flux (W/m2) Sgdelayed neutron source term (neutron/group) ssubscript: scattering cross section (-) tsubscript: transport cross section (-) v velocity (m/s) −→ Wfluid velocity (m/s) wsubscript: wall (-) αcontrol rod fraction (-) βconfidence level (%) βeff effective fraction of the delayed neutrons (-) ∆Bicore average burnup increment in one step (MWd/kg) Γ volumetric mass exchange rate (kg/m3s) γprobability (%) γl ieffective yield of isotope i(atoms/fission) λH jdecay constant for decay-heat group j(sec−1) λl idecay constant of the isotope i(sec−1) νaverage number of neutrons produced per fission (neutron/fission) Σ macroscopic cross section (cm−1) Σsgg′groups-to-group scattering cross section (cm−1) σmicroscopic cross section (cm−1) ρdensity (kg/m3) Symbols xxi ωangular frequency (rads−1) ϕm gmnode averaged neutron flux (cm−2sec−1) Ψmmnode fission source term (neutron/fission) χgaverage fission spectrum (cm−1) Ω unit vector in direction of motion (solid angle) ζk(t) decay constant of the decay heat group i(sec−1) ’superscript: incident neutron energy and direction (-) ±flux direction (-) A la meva fam´ılia xxiii Chapter 2. Background 7 –A general description of the fuel assemblies is given in this section, all the needed values are reported here. These values cover either geometry considerations and neutron modeling, which means the number of prompt and delayed groups, the decay heat constants and other classical considerations needed when modeling the neutron kinetics part of the reactor. Composition maps for the 2D and 3D assembly types is given here. Finally there is a cross section library is facilitated to the participants. •Thermal-hydraulic data –Some geometrical specifications over the different thermal hydraulic components are given in this section. Reactor vessel, vent valves, steam generators, steam lines, feed water system, reactor coolant pumps...are some of the specified components in this section. Break modeling structure is also released here. Finally a set of boundary conditions in terms of Temperature, Pressure and Mass Flow is given in here. •NK-TH coupling guidelines –Some guidelines over the mapping composition for the coupled calculation are given here. Some examples of mapping identification are also released to the participants. •MSLB scenario –Finally a detailed description of the benchmark scenario is given at the end of this first volume. This scenario description includes: Initial steady state conditions; Point kinetics model input and Transient calculations Since a benchmark consists in a comparison between different techniques, user and codes in this particular case, that is why at the end of this first volume the output requested values were listed. The basis of the future comparisons are being set up in this first stage of the benchmark The second report consists in a point kinetics plant simulation. Such simulation models the primary and secondary systems. The aim of the second exercise of the OECD PWR MSLB benchmark [16–19] was to test the thermal-hydraulic system response. Each participant was provided with compatible point kinetics model inputs that preserve axial and radial power distribution, and scram reactivity obtained using a 3-D core neutronics model and a complete system description. First exercise was selected because traditionally the PWR MSLB transient was being modeled with the point kinetics approximation. By choosing the point kinetics approximation, several extremely conservative assumptions has to be taken. These assumptions are generally taken in order to account for Chapter 2. Background 8 the asymmetry in the core region that takes place during the transient. By considering these conservative assumptions the analysis becomes very limited in terms of the total power upgrades and extension of fuel cycles analysis. Also by considering a point kinetics approach, the spatial changes of the power density could not be capt by the nature of the approximation itself. This point kinetics plant simulation exercise intends to provide a detailed description of the simulated main steam line break transient specified for OECD PWR MSLB benchmark problem [16–19]. To overcome some of the limitations of the point kinetics approach, the reactivity feedback components were spatially weighted in radial and axial directions. Nevertheless some parameters should be preserved when running a point kinetics approach is made, in order to make the obtained results comparable with the 3D NK.TH approach. These parameters include: Tripped rod worth; Radial power distribution; Axial power distribution; moderator temperature coefficient; Doppler coefficient and some Kinetics parameters. In the same way some initial boundary conditions need to be identical, these conditions are: Power level; Boron concentration; Axial power shape of the rods; Xe distribution and moderator temperature. A list of neutronic parameters and transient assumptions was distributed to the participants. Finally a standard techniques for comparison data was established in order to compare different calculations. This standard methodology for comparing date consists in 4 steps: •Step 1: Isolate points of interest •Step 2: Calculate mean values and standard deviations •Step 3: Identify outliers and recalculate mean, if necessary •Step 4: Determine and report the deviation and figure of merit for each participants value At the end of the second report, a multiple comparison in between the different participants calculations was made. The key analyzed points where: Break mass flow rate, Reactor power; Pressure; temperatures; Reactivity and steam generator mass. First evaluation of the problem was achieved in this stage of the OECD PWR MSLB benchmark [16–19]. Third report consists in a coupled 3-D neutronics thermal-hydraulics evaluation of core response calculation. The aim of this third exercise of the OECD PWR MSLB benchmark [16–19] was to test the neutronics response to imposed thermal-hydraulic conditions. Each participant was provided with transient boundary conditions (radial distribution of mass flow rates and liquid temperatures at the core inlet, and radial averaged pressure versus time at both the core inlet and outlet), the initial axial liquid velocities, Chapter 2. Background 9 the initial axial distribution of liquid temperatures and a complete core description. When using a 3D approach analysis, all the above mentioned extremely conservative assumptions taken into account when performing a point kinetics calculation, are not necessary. Thus the new 3D NK approach may provide a margin of return to critically status compared over the point kinetic approach. This margin may contribute to the improvement of the operational flexibility and nuclear power plant performance. In the same line as the previous exercise, there were some general specifications released to the participants of the OECD PWR MSLB benchmark [16–19] which need to be taken into account when performing the 3D NK-TH analysis in this case. These specifications cover now: Core neutronics model; Cross section library; NK-TH coupling; initial steady state conditions and transient calculations. Since the aim of this exercise is to test the neutron kinetics model over a fixed boundary conditions, the core TH boundary conditions model was made by defining an inlet condition at the vessel bottom and an outlet condition at the vessel top. The vessel in this case represents an isolated core with boundary conditions at its bottom and top. These boundary conditions are Inlet mass flow rates, Temperatures and inlet/outlet pressures. Those conditions where taken from a TRAC-P/NEM coupled calculation. After all these specifications another multiple comparison was made in the same way as the previous exercise plus a statistical analysis of normalized parameters. Some conclusions were taken at the end of this exercise. Detail of the core modeling and the coupling scheme turns out to be some significant parameters which may cause some deviation in the results. Also the spatial decay heat modeling plus the density and doppler temperatures correlations used by the thermal hydraulics code had a noticeable effect on the different calculation deviations. Fourth part of the project is a best-estimate 3D NK-TH coupled core-plant transient model. Last exercise of the OECD PWR MSLB benchmark [16–19] was to simulate the entire transient and combine the first two exercises, fully testing the thermalhydraulic/neutronic coupling. The different coupled codes predictions are compared and evaluated in regard to: time and value of the power peak before reactor trip; time and value of a power peak after reactor trip; Whether the system remains critical after the momentary return to power (if it occurs) for the transient duration. AS in the previous exercises several boundary conditions for the steady state and transient calculations where released, the difference in this exercise was the completely NK-TH feedback of the proposed exercise. As it was made in previous exercises, some multiple comparisons were made. The comparisons were basically made in two ways: Standard techniques for comparison of results and Statistical analysis of normalized parameters within these techniques, different key parameters were evaluated and compared. As general conclusions for the OECD PWR MSLB benchmark [16–19], a specific list of relevant parameters for comparison in each exercise has been finally determined. More Chapter 2. Background 10 in detail it was determined that for the system behavior prediction in OECD PWR MSLB benchmark [16–19], key parameters were: The SG masses; The break flow rates; The coolant and fuel temperatures and the power. Also there is a big dependency on the TH core modeling. For the MSLB transients there is less dependency on the radial refinements of the neutronic model, this could be different for different scenarios. As conclusions applicable to the present study, OECD PWR MSLB benchmark [16–19] was used as a guideline of how to model and behave over the different stages of the present study performed calculations. From the point kinetics input to the 3D coupled input going through the lattice physics code and cross section generation, all the knowledge learned from OECD PWR MSLB benchmark [16–19], was fully applicable to the present study. 2.2.1 PWR MSLB transient The MSLB is the transient chosen for analysis in the present study. Within the GET group in Technical University of Catalonia there was some gained experience in MSLB scenarios on PWR’s due the participation to the OECD PWR MSLB benchmark [16–19] and OECD NEA PKL-2 [50, 51] project but also due some published articles [58] which give some consistent background on the study of the MSLB scenario in Asc´o NPP. Such knowledge was used in order to reproduce the MSLB scenario with the new developed models. A brief description and some results from the BEBCC (Best Estimate Base Case Calculation) selected is given in this section. This BEBCC is going to be used as a base case for the further comparison with the different methodologies and also as a base line for building the BE plus conservative assumptions case and BEPU case. Posterior analysis with conservative methodology and BEPU methodology are based in the transient presented a later sections of the present study. With all the related knowledge explained above the author consider that it is perfectly suitable to explain the selected transient and to show the BEBCC calculations in this chapter of the present study. A double ended MSLB in loop 2 is the initiating event. Immediately after the break, the high differential pressure between the steam lines causes activation of the high pressure injection systems. At the same time the turbine and the reactor are shut down. For this calculation we have postulated a control rod stuck in the withdrawn position during SCRAM. The high differential heat transfer ratio between the broken loop and the intact loop causes temperature and coolant density asymmetries in primary system, which is propagated into the core. There is some mixing effect between the three loops flows into the lower plenum, but, the cooler water mainly enters into the core region where the control rod is stuck in its fully withdrawn position. There is an increase of the total reactivity, mainly due to the density changes of the coolant (moderator). Table 2.1 Chapter 2. Background 11 shows the sequence of the main events for the calculated main steam line break scenario. Steady state values agreement can be found in table 4.4. Table 2.1: Sequence of events in MSLB scenario. Event Time (s) Double-ended loop 2 Main Steam Line Break opens 15.00 High differential pressure between steam lines signal. 15.05 Safety Injection Signal. SCRAM signal Steam isolation 15.50 Steam generator 2 empties 151.00 Manual AFW turbopump trip 195.00 Manual regulation of AFW valves (15%) 285.00 End of simulation 300.00 The power remains low and decreases quickly as is shown in figure 2.1. The boron injection from safety injection systems and the reactor scram reduce the power during the transient. Nevertheless, in the paper we are taking a closer look at the reactivity increase even if it remains at negative values. The reasons for the increase in reactivity are the local moderator density and the fuel temperature changes as well as the stuck control rod in its withdrawn position. The total reactivity evolution as function of time plus reactivity components are shown in figure 2.2. The three phenomena occur in the core region where a control rod is stuck in the withdrawn position and also where the main part of the coolant flow coming from the broken loops passes through. Figure 2.3 shows the 3D power distribution at steady state condition. The rest of the 3D graphics (figures 2.4, 2.5, 2.6, 2.7 and 2.8) show the evolution of the power during the SCRAM time when the total power decreases from 100% to almost 10% in 25.0 seconds window. These figures are divided in steps of 5.0 seconds wide each step. Notice the Z axis is a relative power. Also the influence of the control rod banks can bee seen in these plots, thus some local depression of the power is observed in the places where the control rod banks are inserted. Chapter 2. Background 12 0 50 100 150 200 250 300 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0 50 100 150 200 250 300 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Power Tim e (sec.) Total power Figure 2.1: Total Power Best Estimate Base Case Calculation. 0 50 100 150 200 250 300 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 0 50 100 150 200 250 300 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 reactivivty ($) Tim e (sec.) Total reactiv ity Moderator density reac. Boron reac. Fuel temp reac. Figure 2.2: Total reactivity and its components in BEBCC. Chapter 2. Background 13 Figure 2.3: 3D power distribution at steady state. Figure 2.4: 3D power distribution during transient step 1. Chapter 2. Background 14 Figure 2.5: 3D power distribution during transient step 2. Figure 2.6: 3D power distribution during transient step 3. Chapter 2. Background 15 Figure 2.7: 3D power distribution during transient step 4. Figure 2.8: 3D power distribution during transient step 5. 2.3 PKL project As it has been mentioned before the MSLB scenario is also been studied by the author of the present study by the participation in the OECD NEA PKL-2 [50, 51] project. OECD Chapter 2. Background 16 NEA PKL-2 [50, 51] project is an extensive test programme which aims to investigate PWR’s design concepts and PWR’s safety issues. The OECD NEA PKL-2 [50, 51] project is mostly focused on boron precipitation processes and complex heat transfer mechanisms which may occur after some postulated scenarios. PKL facility consist in a AREVA’s owned facility located in Germany which performs a serial of tests which are benchmark by several organizations world wide. OECD NEA PKL-2 [50, 51] project is been carried out for several years in different phases: •G1: Systematic investigation of the heat transfer mechanisms in the SGs in presence of nitrogen, steam and water (2 tests, performed in July and August 2008) •G2: Cool-down procedures with SGs isolated and emptied on the secondary side (1 test, 3 runs performed in December 2008) •G3: Fast cool-down transients (main steam line break) (1 test, performed in July 2009) •G4: Accident situation under reflux condenser conditions for new PWR design concept (1 test with two runs, performed in December 2009) •G5: Boron precipitation following large break loss of coolant accidents. •G6: RCS cool-down with void formation in RPV upper head (1 test performed in April 2011) •G7: Counterpart Test with ROSA / LSTF on small break LOCA with Accident Management procedures (1 test performed in July 2011) From the previous list it can be observed that several phenomenon will be studied in the framework of the OECD NEA PKL-2 [50, 51] project. From SG’s heat transfer mechanisms to cool-down scenarios (procedures and fast cool-down transients) also Boron precipitation after LB-LOCA and studies over the RCS cool-down with the presence of void in the RPV upper head are covered in this project. Literature, simulations and analysis over this project are very extensive. Nevertheless in the framework of the present study, G3.1 test participation become relevant to the author of this thesis report because it has helped to achieve a better understanding of the phenomena carried out during the MSLB scenario specially with the heat transfer mechanisms. Participation to the PKL-2 G3.1 [50, 51] test was made by using a RELAP5 3.3 [41–49] model held by the GET in Technical University of Catalonia. The important phenomena which can be observed in this scenario are for the primary side: Heat transfer to secondary side; Cool down rate and temperature distribution in U-tubes; Natural Chapter 2. Background 23 Figure 2.13: PKL-2 test G3.1 Steam Generator pressure. Figure 2.14: PKL-2 test G3.1 Steam Generator outlet pressure. Chapter 2. Background 24 Figure 2.15: PKL-2 test G3.1 Hot Leg temperature. Figure 2.16: PKL-2 test G3.1 Delta temperature. Chapter 2. Background 25 Figure 2.17: PKL-2 test G3.1 Break mass flow rate. By the end of the participation of this test benchmark project, several knowledge was achieved by the author in order to be applied on the future parts of the present thesis report. Within this list we find: Brake modeling issues where applied when modeling the Asc´o NPP MSLB scenario; Key parameters here described where also checked and take it into account in the posterior calculations made for this thesis report; Heat transfer (Primary to secondary) mechanism was also well identified and specially studied due its impact to the return to critically behavior for the postulated BEPU MSLB transient; ECCS injection and FW behavior was also specially studied and take it into account in the next calculations. 2.4 CRISSUE project CRISSUE-S project [20–22] is and international effort made with the collaboration of several institutions and groups in order to establish some guidelines about the actual LWR NPP system modeling. The objectives of the CRISSUE-S project [20–22] can be summarized as follows: To establish a state-of-the-art report on the subject; To provide results of best-estimate analysis of complex transients in existing reactors; To provide recommendations to interested organizations; To identify areas of the NPP design for which the design/safety requirements can be relaxed. CRISSUE-S project [20–22] was divided in three parts: Chapter 2. Background 26 •CRISSUE-S WP1: Data requirements and databases needed for transient simulations and qualifications •CRISSUE-S WP2: Neutronics/Thermal-hydraulics Coupling in LWR technology: State of the art report •CRISSUE-S WP3: Achievements and recommendations report One of the final targets for CRISSUE-S [20–22] activity consists identify and propose an available list of coupled 3D NK-TH. Once such list is formed CRISSUE-S [20–22] activity will be used to defend “acceptability” (or required precision) thresholds for the results of these analysis. The obtained list of transients is specific to the different NPP types such as PWR, BWR and VVER. The acceptability thresholds for calculation precision are general in nature and are applicable to all LWR’s. Finally it is important to remark the creation of a database for the main results of the 3D NK-TH coupled calculations. Following list shows the transients to which under the CRISSUE-S project [20–22] are recommended to be studied with 3D NK-TH tools •PWR transients – MSLB Initiation event is a double ended guillotine break in a main steam line. SG depressurization is followed by cold water injection in the primary side, which leads cold water trough the core. As a results positive reactivity is noticed in a core region. Even there is partial mixing at the lower plenum, the cold water is causing positive reactivity in on part of the reactor, causing some asymmetries in terms of radial power distribution. – LOFW-ATWS Initiating event is in this case the suddenly blockage of the FW pumps. This event leads to a increase of temperature in the primary loop, such increase combined with the modification of moderator density and Doppler effect, contribute to a power decrease. – CR ejection Sudden CR ejection will lead a regional increase of reactivity, which will be a good scenario to be studied with 3D NK-TH coupling techniques due the asymmetries generated for such event. – LBLOCA-DBA Initiating event in this case is the double ended break in the cold leg in between the reactor coolant pump and the reactor pressure vessel. In this scenario, widely used for licensing, 3D NK-TH connection is justified when appears a need to quantify the conservatism introduced by the highly conservative peak factors (PF) for linear power that cause high values for peak cladding temperature (PCT). Chapter 2. Background 27 – Incorrect connection (start up) of an inactive (idle) loop In that scenario the idle loop is assumed to have de-borated water which might lead to asymmetries into the reactor pressure vessel. – MSLB-ATWS even this is a classical DBA transient, its recommendation to be analyzed with 3D NK-TH techniques derives from the bounding nature that this transient might have over the input reactivity of the core and the consideration that the core integrity is predicted in these scenarios. – SBLOCA-ATWS Typical transient (TMI-type) which is originated by small break in the primary loop. 3D NK-TH analysis is justified here due the injections of de-borated water which might flow through the core in specific regions if there is no sufficient mixture in the lower plenum. •BWR transients – TT without condenser bypass available Initiation event is a positive pressure wave which propagates from the turbine isolation valve to the reactor pressure vessel getting into the core from the top and bottom. This cause the void fraction to collapse and such collapse causes a positive reactivity effect. – LBLOCA-DBA Rupture of a recirculation line is studied here, same considerations as the ones made for PWR’s are valid in here in order to propose this transient to be analyzed with 3D NK-TH techniques. – CR ejection Sudden CR ejection will lead a regional increase of reactivity, which will be a good scenario to be studied with 3D NK-TH coupling techniques due the asymmetries generated for such event. – FW temperature decrease-ATWS Malfunction of the FW pre-heaters is supposed here, The cold FW is reaching the core causing some positive reactivity effect. – MCP flow rate increase A sudden increase of the reactor coolant flow due the malfunction of the main coolant pumps is supposed here. – MSIV closure-ATWS Also like TT without condenser bypass available transient, the closure of the MSIV might lead to a void fraction collapse again like FW temperature decrease-ATWS scenario a positive reactivity effect is caused by such collapse. – Stability analysis This transient is being widely investigated for the nuclear scientific community during several years. the application of the 3D NKTH techniques is also being widely proved. the recommended transient for CRISSUE-S program [20–22] can be originated at nominal power and include MCP trip that brings the plant into the exclusion region of the BWR flow map. Chapter 2. Background 28 – Stability analysis-ATWS Same considerations as the previous scenario but assuming a failure of the SCRAM system. •VVER transients – MSLB since the ratios between the SG water and primary circuit in VVER’s are large than in PWR’s, MSLB scenario is expected to less severe. Nevertheless this scenario is being selected to be studied under 3D NK-TH analysis. – LOFW-ATWS There is no special different between this scenario in PWR’s and VVER’s. Also, due the large amount of water, VVER’s scenario is less severe than PWR’s scenario. – CR ejection Sudden CR ejection will lead a regional increase of reactivity, which will be a good scenario to be studied with 3D NK-TH coupling techniques due the asymmetries generated for such event. Same considerations as the ones made for PWR’s. – LBLOCA-DBA Sudden CR ejection will lead a regional increase of reactivity, which will be a good scenario to be studied with 3D NK-TH coupling techniques due the asymmetries generated for such event. Same considerations as the ones made for PWR’s. – Incorrect connection (start up) of an inactive (idle) loop Since the VVER’s are equipped with main isolation valves in hot leg and cold leg, this might introduce few differences on the scenarios, compared with PWR’s – MSLB-ATWS Its recommendation to be analyzed with 3D NK-TH techniques derives from the bounding nature that this transient might have over the input reactivity of the core and the consideration that the core integrity is predicted in these scenarios. Same considerations as the ones made for PWR’s. – SBLOCA-ATWS Typical transient (TMI-type) which is originated by small break in the primary loop. 3D NK-TH analysis is justified here due the injections of de-borated water which might flow through the core in specific regions if there is no sufficient mixture in the lower plenum. Same considerations as the ones made for PWR’s. – Isolation of one loop (ATWS) since VVER’s are equipped with main isolation valves, this new transient is also considered to be analyzed with 3D NK-TH techniques. This scenario is also considered as complement of Incorrect connection (start up) of an inactive (idle) loop. when analyzing any nuclear reactor system transient, there is a list of key parameters to check if they are under the acceptance criteria based on the licensing, experts judgment Chapter 2. Background 29 and safety requirements. CRISSUE-S program [20–22] also produced a reasonable minimum number of quantities of interest when performing a 3D NK-TH transient analysis. This list of quantities might vary depending on the type of transient performed and the type of reactor analyzed. Nevertheless a minimum common list is showed below. Next list will give a vision of the general quantities and its associated errors. •Reactor pressure vessel peak peak of pressure. Acceptable threshold quantity error is 10% of the nominal pressure of the considered system. Acceptable threshold time error is 100% of the BE value. •Time of occurrence of the RPV pressure’s peak. Acceptable threshold quantity error is 2% of the nominal pressure of the considered system. Acceptable threshold time error is 100% of the BE value. •Peak total power if applicable. Acceptable threshold quantity error 100% of the nominal or 300% from the initial, if initial power is smaller than nominal power. Acceptable threshold time error is 100% of the BE value. •CHF or DNB occurrence time. Acceptable threshold quantity error is 20% of the nominal or 100% from the initial, if initial power is smaller than nominal power. Acceptable threshold time error is 20% of the BE value. •PCT occurrence time. Acceptable threshold quantity error is 150K. Acceptable threshold time error is 20% of the BE value. •Maximum fuel temperature and occurrence time. Acceptable threshold quantity error is 200K. Acceptable threshold time error is 20% of the BE value. •Total energy released to the fluid during the transient. Acceptable threshold quantity error is 10% of the energy released to the fluid or 100% of the energy released to the fluid if the initial power is smaller the nominal. Acceptable threshold time error is 20% of the BE value. •Maximum in % of the core in terms of heat transfer area where at any time rod surface area is bigger than 1000K. Acceptable threshold quantity error is 10% of the heat transfer area. Acceptable threshold time error is 20% of the BE value. •Maximum in % of the core in terms of the volume occupied by fuel pins where at any time the fuel temperature is bigger than 3000K. Acceptable threshold quantity error is 10% of the volume occupied. Acceptable threshold time error is 20% of the BE value. Chapter 2. Background 30 Fuel effect is being identified as one of the big contributors to the LWR’s 3D NK-TH analysis. There are multiple fuel factors which could effect on the general 3D NK-TH analysis. Among these factors we find: burnup, power distribution, materials; power peaks; history exposure; general thermo-physical properties of the fuel and fuel failure mechanisms. Concerning the thermo-physical properties of the fuel materials, the recommended library to be used is MATPRO. Fuel failure mechanisms are been also identified as a big contributors in 3D NK-TH analysis. Besides the fuel degradation takes place during the life of the reactor, it could also appears some fuel degradation during the transient time. CRISSUE-S program [20–22] identify and explain the effects of the following fuel failure mechanisms: Manufacturing defects, Primary hydriding; Pellet-clad interaction; Corrosion; Dry-out, Cladding collapse; Grid-rod fretting; debris fretting; Baffle jetting and Assembly damage. The best estimated approach was also made in CRISSUE-S program [20–22]. First a BE versus conservative approach was made in order to evaluate the uncertainties. Once the uncertainties were identified, those were classified in fuel-related uncertainties. These ones where identified as radiolysis in fast reactivity transients; Dynamic sub-cooled boiling; Dynamic CHF; Volume void weighting on heat transfer surface for two fluids; Spacers with mixing vanes. The other source of uncertainties were identified as the uncertainties related to other phenomena or to components. The list for that type of uncertainties is: Valve characteristics; Frictional and discrete pressure losses; Phase separation at tee’s; high transient thermal flux and positive pressure pulse propagation. Last type of uncertainties are the ones related to models and codes. The way of the heat transfer inside the pin; the fuel modeling and the use of the neutron diffusion equation and the associated uncertainties due the methodologies used to solve it and the assumptions taken in order to simplify the problem determine this uncertainties type classification. CIAU method was used to determine the uncertainties, besides that in the framework of CRISSUE-S program [20–22] the CIAU method was extended a number of neutron kinetics parameters which contain uncertainties on the basis of a NK-TH calculation. These parameters are: Rod worth ±10% or ±15% (depending on the reference; Fraction of delayed neutrons β±5%; Doppler coefficient ±20%; Moderator coefficient ±30%; fuel heat capacity ±10% (this is relevant to the TH parameter); Change in the reactivity unit per change in the fuel and moderator temperature when fuel an moderator are ate the same temperature ±3.6x10−4∆ρ/◦C; Critical boron concentration at 100% of the core power ±50ppm; power distribution (at intermediate level and at 100% power) ±0.1∗relative power density for each measured fuel assembly. CRISSUE-S program also identify a list of tools capable to perform the 3D NK-TH analysis. Starting with the thermal hydraulic codes: ATHLET; RELAP5 (NRC version)[41– 49]; RELAP5-3D [31–35] (DOE version); CATHARE-2; TRAC-P, TRAC-M; TRAC-B Chapter 2. Background 31 and POLCA-T. Available neutron kinetics codes: DYN3D; NEM; NESTLE; PARCS [4, 5] and QUABOX. Cross sections generator codes: CASMO, HELIOS and SCALE. Finally some coupled systems like: TRACE-PARCS [1–3] and [4, 5]; RELAP5-PARCS [41–49], [4, 5] and , RELAP5-NESTLE [31–35] and SIMTRAN. Finally a database with the transient analysis results and qualifications requirements is made. CRISSUE-S program [20–22] sets up a list of requirements that should be taken into account when doing safety analysis with 3D NK-TH codes. These requirements are classified in four points. In the first point the CRISSUE-S program [20–22] gives some guidelines over the level of detail of the input deck. Second point gives guidelines and requirements for the thermal hydraulic nodalization with some acceptability criteria for the thermal hydraulic nodalization at steady state and transient steps. Third step gives the guidelines for the neutron kinetics input deck requirements and qualification. Last point is about the qualification requirements for the coupled input. Once all the qualification requirements are been exposed, CRISSUE-S document [20–22], exposes a list of transient related general acceptance criteria to be used in different transients evaluated during the stages of the program. These general acceptance criteria cover BWR stability transient, ATWS transients and rod ejection event transients. This issue is used as a linkage between the regulatory bodies (licensing works) and 3D NK-TH techniques. At the end the knowledge achieved from CRISSUE-S program [20–22] is being useful for industry, regulatory bodies and researchers. A base guideline is being setup from the point of view of 3D NK-TH techniques. The present study was based since the beginning over the CRISSUE-S project recommendations [20–22]. Starting with the model development (Thermal hydraulic model, Neutron kinetic model and coupled model) and continuing with the BEPU analysis considerations. Key parameters to be studied in the present report where also selected within the framework of the CRISSUE-S program [20–22]. Finally base case transient used in the present study was also selected according the CRISSUE-S program [20–22] recommendations. 2.5 PIRT’s studies Phenomenon Identification and Ranking Tables (PIRT’s) [10–15] technique is a structured process to identify safety-relevant/safety-significant phenomena and assess the importance and knowledge base by ranking the phenomena in order to meet some decisionmaking objective. PIRT has been applied to many nuclear technology issues including nuclear analysis in order to help guide research or develop regulatory requirements. PIRT methodology was developed in the 1980’s and it has been widely used ever since. Chapter 2. Background 32 Some decisions where taken during the development of the present work specially when selecting the relevant parameters to be uncertainty propagated in the BEPU calculations performed at the latest sections of this study. These decisions where taken over the extensive literature used to identify the relevant phenomenon involved in the studies scenario, but also where taken over the decision of the advisors (i.e. Ph. D advisors in NK and TH parts respectively) and over the suggestions of the current author of the present thesis report. In that way Literature recommendations, expert advice and user experienced tests where the contributors of the PIRT selection presented in the BEPU calculations chapter. In this section a demonstration of a PIRT’s methodology is given in order to make clear the selection made when performing a BEPU analysis at later phases of this report. The PIRT’s process starts by identifying different phenomena and ranking them by using some criteria which will generate a table from where will be easy to identify the most important phenomenon related to a specific issue. During this phenomena identification, uncertainties associated to each phenomenon needs also to be identified. The selected phenomenon are conditions of a particular reactor, system, component, a physical or engineering approximation, a reactor parameter, or anything else that might influence in the selected analyzed scenario. Each phenomenon is characterized by two three-leveled scales. First three-leveled scale is called Importance, which determine the relevance from each phenomenon over the figure of merit (Figure of Merit, FoM: represents the most relevance time trends in the studied scenario). The levels for this scale are High/Medium/Low.High implies that the ranked parameter has control (i.e. big impact) over the FoM, thus its accuracy should be very high not to introduce big perturbations. Medium implies that the ranked phenomenon has a moderate impact over the FoM, and its accuracy is not as critical as the previous group. Finally Low tells that the phenomenon has no impact or minimal over the FoM. Second three-leveled scale is called Knowledge, which determine the knowledge over the phenomenon. The levels for this scale are Known/Partially Known/Unknown. These levels of knowledge are well quantified and Known means fully or almost fully known (more than 75% of what we could expect to know). Partially Known means knowledge base is moderate (25% to 75% of the knowledge base is established). Finally Unknown means knowledge base is low (less than 25% of the knowledge base is established). As a conclusion from the last scale, if there is any phenomenon tag as Known there is no suggested research to this phenomenon, on the other way if there is any Unknown phenomenon a research over this phenomenon is a priority excluding the case where this phenomenon is being ranked as Low importance in the previous scale. Finally a Partially Known phenomenon implies the a research over this phenomenon is suggested if in the previous scale the same phenomenon was tag as High importance phenomenon. There are several existing PIRT’s applications in thermal-hydraulics, severe accidents, fuels, materials degradation, and nuclear analysis. For each case there is a different objective and the Chapter 2. Background 39 steps will study the uncertainties propagation over the thermal-hydraulics BWR system performance interaction with the Peach Bottom Turbine Trip and BEMUSE-3 experimental data and Coupled neutronics kinetics thermal-hydraulic core/thermal-hydraulic system BWR performance interaction with the Peach Bottom Turbine Trip experimental data and Peach Bottom stability performance interaction with EOC2 and EOC3 experimental data. It is recommended to use as much experimental data as possible when performing each one of the previous mentioned nine steps. By following the proposed methodology steps at the end of the benchmark it is intend to held a mixture of information from ITF, NPP and analytical data which will be compared with the current uncertainty methods and as a result will produce some benefits in the different approaches to arrive at some recommendations and guidelines. As can be seen for the structure of the above mentioned nine steps, the project is quite ambitious and it intend to cover the full scope of uncertainties sources. To the above mentioned tasks in an effective way the OECD UAM project [23] is structured in the three phases with three exercises within each phase: Phase I (Neutronics Phase) •Exercise I-1: Derivation of the multi-group microscopic cross-section libraries (nuclear data and covariance data, selection of multi-group structure, etc.). •Exercise I-2: Derivation of the few-group macroscopic cross-section libraries (energy collapsing, spatial homogenization of cross-sections and covariance data, etc.). •Exercise I-3: Criticality (steady state) stand-alone neutronics calculations with confidence bounds (keff calculations, diffusion approximation, etc.) Phase II (Core Phase) •Exercise II-1: Fuel thermal properties relevant for transient performance. •Exercise II-2: Neutron kinetics stand-alone performance (kinetics data, space-time dependence treatment, etc.). •Exercise II-3: Thermal-hydraulic fuel bundle performance. Phase III (System Phase) •Exercise III-1: Coupled neutronics/thermal-hydraulics core performance (coupled steady state, coupled depletion, and coupled core transient with boundary conditions) Chapter 2. Background 40 •Exercise III-2: Thermal-hydraulics system performance •Exercise III-3: Coupled neutronics kinetics thermal-hydraulic core/thermal-hydraulic system performance By looking at the previous list it is easy to identify that the present study should be situated at the phase III (System Phase). In this thesis report, the BEPU methodology has being also applied into the system phase and it was omitted when creating the cross section library, due a separate work is held in the GET group in that knowledge area. Since the project is still going on, there is no general conclusions which could give as some ideas to apply when performing a full scope BEPU analysis. Besides that there is an expected impact and benefits which may come from due the OECD LWR UAM benchmark [23] activity and will contribute to the LWR’s safety and licensing. The expected points are: Systematic identification of uncertainty sources; Systematic consideration of uncertainty and sensitivity methods in all steps. This approach will generate a new level of accuracy and will improve transparency of complex dependencies; All results will be represented by reference results and variances and suitable tolerance limits; The dominant parameters will be identified for all physical processes; Support of the quantification of safety margins; The experiences of validation will be explicitly and quantitatively documented; Recommendations and guidelines for the application of the new methodologies will be established. At the conclusion of the present study, the models will be ready to perform a full scope BEPU analysis, for future work it will be left to implement the actual models over the OECD LWR UAM benchmark [23] derived methodologies. Chapter 3 Codes and Models Different qualified tools are being used when performing the present study. Some description from all of the nuclear codes used is given in the first part of this chapter. Basic field equations from the lattice physics, cross section generation and treatment, core simulator, thermal-hydraulic and uncertainty propagation codes are given here. The reader can get and idea of how complex are the problems to be solved and which are the assumptions taken by the codes in order to obtain elegant and satisfactory solutions to each phenomena which needs to be simulated. In the second part of this chapter the author’s developed models are presented. These models where built from scratch by using the expertise gained in different contributions made in the chapter’s two mentioned reference benchmarks, such: OECD LWR UAM benchmark [23]; OECD PWR MSLB Benchmark [16–19]; OECD NEA PKL-2 [50, 51] project and CRISSUE-S project [20–22]. Also Ph.D advisors guidance was an important asset here when building the different models. Finally some trial and error method was also used in order to obtain the best optimized model for each particular case. Some deficiencies where detected and some ways of improving the different models are given at the conclusions chapter. Nevertheless due the complexity of the BEPU calculations, the used models have being proved as the most effective ones with used computing machines. 3.1 Brief description of the codes This section presents a review of the state of the art tools used for this Ph.D study. In that sense a brief description of the used computer codes is given here. The present study involves several codes. Since the final computation will be a Best Estimate Plus Uncertainties calculation in a coupled 3D NK-TH model, a thermal hydraulic system plant code is required to model full plant specifications. TRACE v5.0 patch 2 [1–3] was 41 Chapter 3. Codes and Models 42 the chosen code for the present study. A core simulator code is also required in order to simulate the code behavior under the 3D kinetics perspective. PARCS v3.0 [4, 5] which is internally coupled and compiled as one executable file, with TRACE v5.0 patch 2 [1–3] is the neutron kinetic code used. The cross section library will be the source of power used for the core simulator code. The cross section library contains all the core specifications along the life cycle of the core. An external code is needed to perform such calculations. HELIOS-1.9 [27–29] is the lattice physics code used for such purpose. Once the cross section library is created, specific code call GenPMAXS v5.0 [30] is used in order to convert data from the HELIOS-1.9 [27–29] output file to the PARCS v3.0 [4, 5] input file. No Uncertainty specifications have been considered at this point. The DAKOTA [6–9] code was chosen to perform the perturbation of several parameters from the thermal-hydraulic code and from the neutron kinetics code. Finally, all those steps have been performed, within the framework of the SNAP v.2.2.1 [60] platform. This is a visual interface that allows the user to build models; change specifications; launch calculations and essentially work with all the above mentioned codes, (except HELIOS-1.9 [27–29]) under same Windows mask. 3.1.1 TRACE TRACE [1–3] TRAC/RELAP [41–49] Advanced Computational Engine is the selected thermal-hydraulic code for the present study, this section is giving a brief description of the code operation procedures. The version used in the current study was the TRACE v5.0 patch 2 [1–3]. TRACE [1–3] is a Best Estimate code designed to perform analysis over the different scenarios for the different types of the LWR’s. As a thermal hydraulic system code TRACE [1–3] is capable to simulate all different thermal-hydraulic phenomena that may occur in test facilities and full scale reactors. The code characteristics which are used to predict different phenomena like multidimensional two phase flow, nonequilibrium thermodynamics, generalized heat transfer, reflood, level tracking, reactor kinetics, comprehensive heat transfer capability TRACE [1–3] system code is organized in cards. Within the cards the user can model the different features of the components. There are several thermal-hydraulic components available when modeling with TRACE [1–3] code: PIPE, VALVE, PUMP, PLENUM, PRIZER, CHAN, TEE, TURB, VESSEL, CONTAN and SEPD by using a combination of these components the user can build the thermal-hydraulic part from a full plant nuclear reactor system. Heat conduction properties are modeled by using: HTSTR and REPEAT-HTSTR. These elements are what we call passive heat structures elements. To produce/release heat to the fluid POWER component coupled to a HTSTR is used. There are also a FLPOWER component which is able to deliver power directly to the fluid as it can happen into waste Chapter 3. Codes and Models 43 transmutation facilities. RADENC component are used to simulate radiation enclosures between multiple arbitrary surface. Finally FILL and BREAK components are used to create boundary conditions to the system such mass flow rates or pressures. Besides the previous list of components the code has a bigger list of CONTROL SYSTEM components which are used to simulate the plant control systems. Control system is organized in Variables,Control Blocks and Trip System. Starting with variables there are different types like: Controlled variables used for example to know the pressure in a tank; Manipulated variables (i.e. modify some conditions to achieve a desired value, valve area for example). Next part is the control block.Control blocks are functions which operate over the signals to generate a output signal according the selected function. Variables information can be modified from one block to the other with a mathematical function here is where the transfer functions take and important role. Transfer functions are functions which define gain, delays, arithmetic relationships. . . from one control block to the other one. Usually there is some feedback effect in the flow path where the output signal of one control block is used to feed the input path of itself in order to reduce the produced error. Finally Trips are ON/OFF switches that can be used to generate a hardware action (i.e. open/close a valve), to define a signal’s status or to define a blocking or coincidence trip. All the previous control system are organized in a system called control block diagram which is giving the relationship between the above signals and the flow path of the information in between the control systems. Thermal-hydraulic system in a LWR is very complex and thermal-hydraulic codes need to reproduce such system with enough accuracy to perform validated analysis over the different plants and scenarios. Several phenomena are involved within a thermalhydraulic LWR system. Following list gives an idea of all the physical phenomena which are considered in TRACE [1–3] code analysis: ECC downcomer penetration and bypass, including the effects of countercurrent flow and hot walls; lower-plenum refill with entrainment and phase-separation effects; bottom-reflood and falling-film quench fronts; multidimensional flow patterns in the reactor-core and plenum regions; pool formation and countercurrent flow at the upper-core support-plate (UCSP) region; pool formation in the upper plenum; steam binding; water level tracking; average-rod and hot-rod cladding-temperature histories; direct injection of sub-cooled ECC water, without artificial mixing zones; critical flow (choking); liquid carryover during reflood; metal-water reaction; water-hammer pack and stretch effects; wall friction losses; horizontally stratified flow, including reflux cooling; gas or liquid separator modeling; non-condensablegas effects on evaporation and condensation; dissolved-solute tracking in liquid flow; reactivity-feedback effects on reactor-core power kinetics; two-phase bottom, side, and top offtake flow of a tee side channel; reversible and irreversible form-loss flow effects on the pressure distribution. There are some limitations of use when working with TRACE Chapter 3. Codes and Models 44 [1–3] code. Typically these type of system codes are only applicable in their assessment range of values. Notice that TRACE [1–3] is been qualified to analyze ESBWR design, conventional PWR and BWR Large and Small break LOCA. At this point TRACE [1–3] code is not being validated against BWR stability analysis or other operational transients. What is needed to model and obtain a realistic solution from thermal hydraulic system: •Simplified Vapor/Liquid balance equations (energy, mass and momentum). •State relationships –Relationships between thermal hydraulic variables for and specific fluid, for example water or heavy water. –Library with all the thermophysical and thermodynamical properties (β, k, Cp, ... ). •Jump conditions –Link to de decoupled balance equations. –to express continuity of mass, momentum and energy. –Γf=−Γg •Closure equations –List of correlations that computes independently interphase interactions such mass/heat exchange and dragging as well as wall-to-fluid heat exchanges and frictions. –Correlations set empirically through separate effects tests facilities SETS. –Validations through SETS’s and integral test facilities ITF’s. Best Estimate codes balance equations result from a simplification of Eulerian equations. Navier Stokes equations (no viscosity) + incompressibility: ∂(ρkΨk) ∂t =−∇·(ρkΨk Wk)−∇· JΨ,k +ρkSΨ,k (3.1) First term, left handed, is variation in time, first on the right hand is convection due the fluid motion, second term on the right hand is the diffusion term last term is the volumetric source term. With the following assumptions: •Mass balance equation Ψk= 1,  JΨ,k = 0 and SΨ,k = 0 Chapter 3. Codes and Models 45 •Momentum balance equation Ψk= Wk, JΨ,k =pk I−τkand SΨ,k =g •Energy balance equation Ψk=uk+ W2 k 2, JΨ,k =q′′ + (pk I−τk)· Wk and SΨ,k =q′′′ ρk+g  Wk Simplifications applied to the Eulerian equations are: •Space averaging over the control volume which neglects the turbulent fluctuations. •1D motion fluid which implies that local gradients and fluxes are not considered. Besides that, TRACE [1–3] code has a special solution with (3D) equations that can be applied in some components like vessel. •There are independent de coupled fluid phases which implies that there is no liquid and vapor interactions. •Hyperbolic solution. •Time. •Added artificial viscosity terms on the computations. the fluid field equations required to be solved in this type of systems are mentioned next. Mass equations 3.2, Energy equations 3.3, Momentum equations 3.4. ∂ρ ∂t +∇·(ρv) = 0 (3.2) ∂ρv ∂t +∇·(ρv2) + ∇(p)−∇·(T) =  Fext (3.3) ∂[ρ(u+1 2v2] ∂t +∇·[ρv(h+1 2v2]−∇·(T·v) = Qint +Qext + F·v (3.4) Unknowns from the equations 3.2, 3.3 and 3.4 are henthalpy, ppressure and v velocity. Last equations are averaged in time and volume for single phase gas, single phase liquid and combined with interface jump conditions. The fluid at each node is considered with single velocity, single energy and single pressure. Equations must be solved for liquid and gas phase, so the problem ends up with six field equations, three for liquid phase and three for gas phase. Non-condensable gasses and solute liquid are also considered in TRACE [1–3] code. TRACE [1–3] is capable to model on non-condensable gas as a regular option, but it can Chapter 3. Codes and Models 46 support multiple gas species if required. Equation 3.5 is the non-condensable mixture gas equation and by using it the mechanical equilibrium is assumed this is to assume that the non-condensable gas mixture is in thermal equilibrium with present steam and to move at same velocity. ∂(αρa) ∂t +∇·[αρava] = 0 (3.5) ∂[(1 −α)mρl] ∂t +∇·[(1 −α)mρlvl] = 0 (3.6) TRACE [1–3] code also includes a mass-continuity equation 3.6 for a solute moving within the liquid field, where mis the solute concentration (mass of solute/unit mass of liquid water). The solute concentration is not affecting hydrodynamics directly, but its effects over the reactivity feedback could have some effect over the hydrodynamics. More physics phenomena considered in the code are the drag models. The liquid field and gas field momentum equations, include terms of the interfacial shear force and wall drag force. Drag coefficients Ciinterfacial; Cwl wall liquid; Cwl wall gas are required to solve the closure equations. Different values of the coefficient are applied depending on the flow regime, there are different considered flow regimes depending on vertical flow, horizontal flow. Figure 3.1 shows a schematic representation of the different flow regimes. Both the models for interfacial drag coefficient and wall drag coefficient are dependent upon the flow regime. Figure 3.1: Schematic representation of the vertical flow regimes. Chapter 3. Codes and Models 47 The general equation which describes the heat conduction process in an arbitrary geometry is the following: ∂(ρCpT) ∂t +∇·q =q′′′ (3.7) assuming constant ρCpTand expressing the heat flux as temperature gradient by the Fourier’s law, 3.8: q =−k∇T(3.8) Equation 3.7 becomes: ρCp ∂(T) ∂t =∇·(k∇T) + q′′′ (3.9) Equation 3.9 its the reference equation for the heat conduction treatment within the TRACE [1–3] code. Besides that, there are different geometries which can contribute to release heat flux to the fluid. Therefore the 3.9 equation, must be applied to those different geometries to find the correct expression in order to determine the coupling of heat conduction field equation to the thermal hydraulic field equation. The geometries included in the TRACE [1–3] code are: Cylindrical walls, slabs and cire fuel rods. TRACE [1–3] code also account for the interfacial heat transfer and for the wall heat transfer. Different models and approximations are used for each geometry. Starting with the interfacial heat transfer, these models are needed for the mass and energy closure equations. In TRACE [1–3] code the interfacial mass transfer rate per unit volume, Γ, is expressed as the following equation: Γ = Γi+ Γsub (3.10) Which means the sum of mass transfer rates from the interfacial heat transfer and from the sub-cooled boiling. There are different considerations depending on the type of flow regime and heat transfer stage where the fluid and the wall are encountered. Pre-CHF interfacial heat transfer models which describe the interfacial heat transfer before the CHF occurs. Stratified flow interfacial heat transfer models, which is used for horizontal and inclined pipes where the flow becomes stratifies at low velocity conditions, as gravity and might cause the phases to separate. Post-CHF interfacial heat transfer which describes the interfacial heat transfer for the inverted flow. These situations may happen when the surface temperature is too hot for the liquid phase to contact the wall. Chapter 3. Codes and Models 48 also non-condensable gases effects are considered because its presence effects the mass transfer processes of condensation and evaporation. Also wall heat transfer considered in TRACE [1–3] with different models depending on the heat transfer situation. Wall heat transfer models are required for the mass and energy closure equations. Following equations 3.11 represents the heat transfer rate per unit volume from wall to the liquid and from wall to the gas-vapor mixture. q′′′ wl =hwl ·(Tw−Tl)·A′′′ w q′′′ wsat =hwsat ·(Tw−Tsat)·A′′′ w q′′′ wg =hwg ·(Tw−Tg)·A′′′ w (3.11) Where hwl wall to liquid heat transfer coefficient; hwl is the heat transfer coefficient for the direct boiling of the liquid; hwg wall to gas heat transfer coefficient and the wall heat transfer surface area per unit volume is A′′′ w= 4/Dh. TRACE [1–3] code holds a library of heat transfer correlations plus a selection algorithm which are used to calculate these heat transfer coefficients. By joining all these library correlations and algorithms the code is producing a continuous boiling curve where the more realistic heat transfer coefficient is selected at each step. The models can be divided in the following parts: Pre-CHF transfer (models for wall -liquid convection, nucleate boiling and subcooled boiling); CHF transfer (models for peak heat flux in nucleate boiling heat transfer regime and the wall temperature at which occurs); Minimum film boiling temperature (the temperature above which wall liquid contact does not occur); Post-CHF transfer (models for transition and film boiling heat transfer) and Condensation heat transfer (models for film boiling condensation and the non-condensable gas effect). Finally, after all these balance and closure equations there are some additional special processes correlations are added in the system to simulate and compute some local phenomena that might have a significant impact on the global thermal hydraulic behavior of the system. Those special processes correlations are representing phenomena that will not be capt in the solution because of the simplifications made on the balance and closure equations. These processes are used depending on the particular fluid conditions at each time step, the user needs to activate the computation of each special process where he (i.e. volume, junction, pump. . . ) thinks it will be required. In that sense experience and knowledge of the problem will tell where and when activate or deactivate such capabilities of the code. Basically these special processes SP, are recomputing some particular thermal hydraulic parameters once the system equations are been solved. Following list describes some of the special processes in TRACE [1–3] code: Chapter 3. Codes and Models 55 expensive nodal calculations during some time steps in the transient. There is also a transient fixed source problem at each time point in the transient. To solve the spatial discretization, there are different solution kernels available like, ANM and NEM those include the most popular LWR two group nodal methods. In order to minimize the computational time, there are also some well known computational methods for example, the solution of the CMFD linear system is obtained using a Krylov subspace method. The eigenvalue calculation to establish the initial steady-state is performed using the Wielandt eigenvalue shift method. PARCS [4, 5] code is written in FORTRAN90 language. Its portability has been tested on various platforms and operating systems, to include SUN Solaris Unix, DEC Alpha Unix, HP Unix, LINUX, and various Windows OS. The following list shows more detailed explanation of the PARCS features, as mentioned above the aim of this work is not in the field of Computer Science, nevertheless a brief description of the capabilities of the code is given. 3.1.2.1 PARCS calculation features •Eigenvalue problem –PARCS [4, 5] code is able to solve two kinds of neutronic problems. Eigenvalue problem and fixed source problem. Neutron flux solver needs to solve the nodal balance equation under the cartesian geometry approximation. 1 Vm g dϕm g dt =1 keff χpg G ∑ g=1 νpgΣm fgϕm g+χdg K ∑ k=1 λkCm k+ + G ∑ g′=1 Σm gg′ϕm g′−∑ u=x,y,z 1 hm u (Jm+ gu −Jm− gu )− m ∑ tg ϕm tg (3.15) and dCm k dt =1 keff G ∑ g=1 νdpkΣm fgϕm g+λkCm k(3.16) where Cm kis the precursor density, ϕm gis the node averaged flux and Jm+ − gw is the surface averaged net current. The index Gstands for the neutron energy group and Kindex is representative of the delayed neutron group. Finally plus and minus signs determine the flux direction pand dsubscripts stand for prompt and delayed neutrons. The eigenvalue keff it is determined during the steady state process and it is kept it constant during the transient calculation. There is no difference between delayed and prompt neutrons when solving Chapter 3. Codes and Models 56 steady state problem, also all of the derivative parts from equations 3.15 and 3.16 become zero. Thus the new equation can be written as follows: Mϕ =λFϕ ≡1 keff Fϕ (3.17) Where Fis the fission matrix and Mis the non-fission matrix. Equation 3.17 is solved in the code by using the fission source iteration method. Such solution is achieved via Wielandt eigenvalue shift method rather then common Chebyshev polynomial method. This is done this way because the Wielandt solution will be helpful when solving the transient fixed source problem on the next steps. •Transient fixed source calculation –Fixed source calculation can be used in different occasions, besides that in PARCS [4, 5] code is commonly used for the spatial kinetics problem. Time dependent solution for the equation 3.15 is solved on this step of the calculation. Nevertheless, as mentioned several times, PARCS [4, 5] code is very LWR orientated and there are several approximations which are taken according to that orientation. One of these orientated approximations taken can be seen in transient fixed source calculation where several approximations are made in order to make things more suitable to the problem and to include the use for the two energy groups. The simplifications are the following: ∗χp1=χd1= 1.0 and χp2=χd2= 1.0 ∗There is no dependence on the delayed neutron precursor yields on neutron energy ∗There is no up scattering Σm 21 = 0 ∗νdgkΣm fg =βm kνΣm fg and νpgΣm fg = (1 −βm k)νΣm fg where βm≡∑K k=1 βm k With the above approximations the two group kinetics equation takes the following form: 1 Vm g dϕm g dt =Rm g=   (1 −βm)Ψm+Sm d−Lm 1−Σm r1ϕm 1g= 1 Σm 12ϕm 1−Lm 2−Σm r2ϕm 2g= 2 (3.18) and dCm k dt =βm KΨm−λKCm k(3.19) Chapter 3. Codes and Models 57 where from the above equations we have the following definitions; Sgdelayed neutron source ; Lggroups leakage; Ψ total fission source term; These definitions take the following forms: Ψm≡1 keff 2 ∑ g=1 νΣm fgϕm g Sd m≡ K ∑ k=1 λkCm k Lg m∑ u=x,y,z Lm gu Lm gu ≡1 hm u (Jm+ gu −Jm− gu ) (3.20) The transient fixed source problem is formulated from the two groups kinetic equations by applying temporal discretization methods, in here CMFD formulation and two node nodal method are used to solve the flux equations. At the end the resulting transient fixed source problem will contain only node average fluxes, as the unknowns. •Numerical solution methods –In PARCS [4, 5] the primary solution algorithm is based on the nonlinear Coarse Mesh Finite Difference (CMFD) formulation. In the CMFD method the core is discretized into coarse mesh, typically the size of a fuel assembly. Finite difference discretization is applied between mesh. Each balance equation from each node is coupled to each balance equation of each one of the neighboring nodes using the leakage term. This nodal coupling is solved using a non linear nodal method, where the interface current between two nodes is represented by the average fluxes of the two facing nodes. At the end the PARCS [4, 5] code is solving a matrix system like equation 3.21 where the unknowns are the node average fluxes. The system is called Coarse Mesh Finite Difference Method for a transient fixed source problem. The system can be solved by using any iterative linear system solution method, Kyrlov subspace method, BiCGSTAB,GMRES and BILU3D pre-conditioners are used to solve the CMFD matrixes. Anϕn=Sn(3.21) These Kyrlov subspace methods are one of the most effective ways on solving linear systems. A consecution of the different above mentioned algorithms constitutes the solution methodology for the CFMD numerical method. The Chapter 3. Codes and Models 58 algorithm methodology goes beyond the scope of the present study and it has only being mentioned here in order to satisfy reader’s curiosity. More information can be found in the following references [62, 63] •Nodal diffusion methods –Nodal methods are the primary means used in PARCS [4, 5] to obtain higher order solutions to the neutron diffusion equation. Within the framework of the CMFD formulation the nodal method is used to solve the two node problem and to update the nodal coupling coefficient. The ANM (Analytic Nodal Method) in PARCS [4, 5] has been used most frequently within the Light Water Reactor (LWR) industry to solve the two group diffusion equation. For the previous computational solutions the CNCC corrective nodal coupling coefficient was assumed to be known. Besides that CNCC should be computed during the nonlinear iteration involved in the process of the CMFD and the two-nodes calculation. CNCC is determined when the interface current obtained by the CMFD is the same as the nodal interface current obtained from the two node calculation. The two node problem is a (1D) problem, for which the analytic solution is readily obtainable. In PARCS [4, 5] code the (1D) diffusion equation is obtained through the integration of the (3D) steady state neutron diffusion equation over the transverse plane also there is a common approximation used in the all transverse-integrate nodal methods which consists to assume a quadratic spatial variation of the transverse leakage. The NEM Nodal Expansion Method is used to find the solution for the steady state multi-group diffusion equation in cartesian geometry. The principal features of the polynomial nodal method are the quadratic expansions for the (1D) transverse integrated flux and for leakage model for the transversal leakage. The general multi-group neutron diffusion equation is written as:  ∇·Dg ∇ϕg+ Σtgϕg= G ∑ g=1 Σsgg′ϕg′+χg k G ∑ g=1 vg′Σfg′ϕg′(3.22) where Dgis the diffusion coefficient in (cm); ϕgis the neutron flux in (cm−2sec−1); Σtg is the total macroscopic cross section (cm−1); Σsgg′is the groups-to-group scattering cross section (cm−1); χgis the fission neutron yield; kis the multiplication factor (i.e. critical eigenvalue); vgis the average number of neutrons created per fission and Σfg is the macroscopic fission cross section (cm−1). •Xenon/Samarium calculation Chapter 3. Codes and Models 59 –Xenon and Samarium concentrations are also taken into account into PARCS [4, 5] code. Such concentrations will have an important effect in those transients where the power is shifted from upper to lower levels and vice versa. Reactivity variations will occur also due the presence of high concentrations of Xenon and Samarium isotopes. Xenon and Samarium precursors are Iodine and Promethium respectively. Time depletion equations of the fission products of the both chains looks like: dNl I(t) dt =γl I G ∑ g=1 Σl fg(t)ϕl g(t)−λl INl I(t) (3.23) dNl Xe(t) dt =λl INl I(t) + γl Xe G ∑ g=1 Σl fg(t)ϕl g(t)−λl XeNl Xe(t)− − G ∑ g=1 σl Xe,ag(t)ϕl g(t)Nl Xe(t) (3.24) Equations 3.23 and 3.24 are for the Xe135 and I135 decay chain. Equations 3.25 and 3.26 are for the Sm149 and Pm149 decay chain. dNl Pm(t) dt =γl Pm G ∑ g=1 Σl fg(t)ϕl g(t)−λl PmNl Pm(t) (3.25) dNl Sm(t) dt =λl PmNl Pm(t)− G ∑ g=1 σl Sm,ag(t)ϕl g(t)Nl Sm(t) (3.26) where Nl i(t) is the nuclei number density of isotope i,σl i,ag(t) is the group-wise microscopic cross section of the isotope i;γl iis the effective yield (atoms/fission) of isotope ifinally λl iis the decay constant of the isotope i. Previous equations are used to calculate the time dependent densities of the Xe and Sm isotopes and its precursors. To obtain the steady state number densities it is necessary to integrate the previous equations which lead to the following ones: Nl I,∞=γl I∑G g=1 Σl fgϕl g λl I (3.27) Nl Xe,∞=λl INl Xe,∞+γl Xe ∑G g=1 Σl fgϕl g λl Xe +∑G g=1 σl Xe,agϕl g (3.28) Equations 3.27 and 3.28 are for the Xe135 and I135 decay chain. Equations 3.29 and 3.30 are for the Sm149 and Pm149 decay chain. Chapter 3. Codes and Models 60 Nl Pm,∞=γl Pm ∑G g=1 Σl fgϕl g λl Pm (3.29) Nl Sm,∞=λl PmNl Pm,∞ ∑G g=1 σl Sm,agϕl g (3.30) To obtain transient concentrations (t+∆t) is added in the last set of equations. After all these calculations are finished, the resulting number densities are used to update the macroscopic cross section as it is shown in the equation 3.31: σl ag =σl ag + ∆σl Xe,ag + ∆σl Sm,ag (3.31) where: ∆σl Xe,ag = Σl Xe,agNl Xe and ∆Σl Sm,ag =σl Sm,agNl Sm (3.32) If the Xe and Sm contributions are not been calculated by the user, when running the lattice physics code, PARCS [4, 5] code, holds a list of default values for Xe and Sm contribution, either number densities either corrections to be added at each cross section. Also the code has the capability to run different calculations with or without such contribution (i.e. Xe in equilibrium or non-equilibrium, Sm present or not present). In the present study Xe and Sm contribution was included when building the cross section library. •Neutron Transport Methods –Although PARCS [4, 5] code is basically a neutron diffusion code, diffusion equation is not accurate enough for some of the computations required. PARCS [4, 5] code holds a multi-group transport calculation capability. Spherical harmonics method PNis the most common method used to solve the multi-group transport equation. Equation 3.33 is the steady state Boltzmann transport equation without an external source: Ω·∇Ψ(r, Ω, E)+Σt(r, E)Ψ(r, Ω, E) = =∫dΩ′∫dE′Σs(r, Ω′→Ω, E′→E)Ψ(r, Ω′, E′) + 1 4πSf(r, E)(3.33) where Sf(r, E) = χ(E)∫dE′νΣf(r, E′)ϕ(r, E′) (3.34) Chapter 3. Codes and Models 61 and ϕ(r, E) = ∫dΩ′Ψ(r, Ω′E) (3.35) Since PNis referred to one dimensional problem, the generalization of this method to a (3D) problem is known as SPNapproximation. To do that several changes need to be done in the original PNequations. In that sense, the three dimensional P1can be built from the P1one-dimensional equations by replacing ∂/∂x operator in one dimensional n= 0 with the divergence operator ∇; replacing ∂/∂x operator in one dimensional n= 1 with the gradient operator ∇; considering zeroth-order Legendre moment of the angular flux ϕ0as scalar and considering first-order Legendre moment of the angular flux ϕ0as vector. For SPNthe relations between the geometries are extrapolated keeping in mind to replace ∂/∂x operator in one dimensional for even nwith the divergence operator ∇; replacing ∂/∂x operator in one dimensional for odd nwith the gradient operator ∇; considering even-order Legendre moment of the angular flux as scalar and considering odd-order Legendre moment of the angular flux as vector. This methodology is well known process used to solve multi-group transport equation. PARCS [4, 5] code has its own SPN development and in the actual versions that methodology is truncated at SP3 for N > 3. This ends up with the following SP3equations: ∇·ϕ1g+ Σrgϕ0g=S0g 2 3∇·ϕ2g+1 3∇·ϕ0g+ Σtrgϕ1g= 0 3 5∇·ϕ3g+2 5∇·ϕ1g+ Σtgϕ2g= 0 3 7∇·ϕ2g+ Σtgϕ3g= 0 (3.36) where S0g=∑ g′ Σsg′gϕ0g′+χg keff ∑ g′ νΣfg′ϕ0g′(3.37) Equations 3.36 are the time independent equations resulting for the SP3 method, for the time dependent solution the equations 3.38 show the addition of the time derivative terms necessary for the time dependent solution. Chapter 3. Codes and Models 62 1 υ ∂ϕ0g ∂t +∇·ϕ1g+ Σrgϕ0g=S0g 1 υ ∂ϕ1g ∂t +2 3∇·ϕ2g+1 3∇·ϕ0g+ Σtrgϕ1g= 0 1 υ ∂ϕ2g ∂t +3 5∇·ϕ3g+2 5∇·ϕ1g+ Σtgϕ2g= 0 1 υ ∂ϕ3g ∂t +3 7∇·ϕ2g+ Σtgϕ3g= 0 dCk dt =−λkCk+β keff ∑ g′ νΣfg′ϕ0g′ (3.38) where S0g=∑ g′ Σsg′gϕ0g′+χg keff (1 −β)∑ g′ νΣfg′ϕ0g′+χdg ∑ k′ λkCk dCk dt =−λkCk+βg keff ∑ g′ νΣfg′ϕ0g′ (3.39) Fine Mesh Finite Difference Method, FMFDM is introduced in PARCS [4, 5] for solving the SP3equations. This becomes convenient when computing FA with heterogenous conditions. Some inefficiencies in terms of the accuracy of the solution and computational time where encountered when using the FMFDM and that is why PARCS [4, 5] code has an advanced nodal solver method called NEM Nodal Expansion Method. •Hexagonal modal methods –All the analyzed systems in the present study, were based on cartesian geometry solutions, nevertheless PARCS [4, 5] code, holds the capability of solving the neutron diffusion equation for the Hexagonal geometries (i.e. VVER reactors). The hexagonal nodal method needed to solve the neutron diffusion equation is made with TPEN Triangle-based Polynomial Expansion Nodal. TPEN solves two transverse-integrated neutron diffusion equations for a hexoctahedron node. Since this capability goes beyond the aim of the present study it is just mentioned here without getting into much detail than adding it into the list of the capabilities of the code. •Fuel depletion analysis –Depletion capability was added into the code in order to make it able to perform fuel cycle analysis. Burnup history and power needs to be entered to the code to perform such analysis. These information is entered via GenPMAXS [30] and needs to be computed when generating the cross section library with Chapter 3. Codes and Models 63 the lattice physics code. This information was added on the computation in the present study. The burnup capability is structured in steps. Each step is taking the power from each node to compute the advance burnup from each node. With this capability the user is able to run the system through all the stages of the life cycle of the reactor. The user can specify different burnup distributions for different FA and PARCS [4, 5] will compute the advancement step for the cycle of the reactor. This calculation will give as a result a new distribution in terms of burnup and cross section information that the previous one. The user can pick any point from the pre-desired steps to start point of a transient analysis. The burnup distribution is computed using the fluxes provided by PARCS with the following equation: ∆Bi= ∆Bc Pi Gi Pc Gc (3.40) where ∆Biis the core average burnup increment in one step, specified in the depletor input, this is an user decision. Giis the heavy metal loading in ith region. Gcis the total heavy metal loading in the core. Piis the power in the ith region and finally Pcis the total power in the core. History variables needs to be balanced in this section as well, the pass history of each fuel assembly will have an effect on the future burnup distribution. The considered history contributions are: Control rod history (HCR); Moderator density history (HMD); Soluble boron history (HSB); Fuel temperature history (HTF); Moderator temperature history (HTM). All these history variables are defined as weighted quantities as follow: HCR(Bp+ ∆B) = ∫Bp+∆B 0α(B)dB Bp+ ∆B=HCR(Bp)×Bp+α∆B Bp+ ∆B(3.41) where HCR is the control rod history for this example but it could be any one of the above mentioned history variables. Bpis the burnup at the beginning of this step and ∆Bis the burnup increment. αis the rodded fraction during this step, this procedure is used for the rest of history variables in order to obtain the burnup and power distribution for the next burnup step. If this feature is used the cross section will have two main contributions: ∗Contribution coming from the instantaneous variables, i.e. control rod insertion, moderator density, moderator and fuel temperature and soluble boron concentration. Chapter 3. Codes and Models 64 ∗Contribution coming from the history variables, i.e. control rod history, moderator density history, soluble boron concentration history and finally fuel and moderator temperature. The code will produce the resulting cross sections sets in two steps. First it will consider the all the instantaneous values of all dependent variables for a specified history state of a particular region to produce a cross section set depending on the history variables, this job is made by the DEPLETION module. In the Second stage the code will produce the cross section set by taking the cross section generated by the DEPLETION module and correcting it with the current instantaneous variables. •Decay heat calculation –PARCS [4, 5] is also capable of performing core depletion analysis. Burnup dependent macroscopic cross sections are read from the PMAXS file prepared by the code GENPMAXS [30] and the PARCS [4, 5] node-wise power is used to calculate the region-wise burnup increment for time advancing the macroscopic cross sections. Details of the PMAXS file and the GENPMAXS [30] code are provided in the GENPMAXS [30] description section. The amount of computed heat that will be released during SCRAM situation will depend on the burnup history of the FA which conforms the core, different enrichments and FA positions will have its importance in this calculation. Such option is being activated on the present study. Following equation gives the volumetric heat density with the decay heat contributions considered: qt(r, t) = (1 −αt) G ∑ g=1 kgΣfg(r, t)ϕg(r, t) + I ∑ i=1 ζiDi(r, t) (3.42) where Di(r, t) is the concentration of the decay heat precursors in decay heat group i(J/cm3); ζiis the decay constant of the decay heat group i(sec−1); αt=∑I i=1 αiis the total fraction of the fission energy appearing as decay heat where Iis the total number of decay heat groups; Finally αiis the fraction of the total fission energy appearing as decay heat for decay heat group i. After some modifications and simplifications, the concentration of the decay heat precursors D(r, t) becomes: Di(r, tn+1) = Di(r, tn)eζi∆t+αi ζi [1 −eζi∆t] G ∑ g=1 κgΣfg(r, tn)ϕg(r, tn) (3.43) •Pin power Calculation Chapter 3. Codes and Models 71 Figure 3.3: General information exchange flow diagram for a GenPMAXS code. to find a intermediate value from a two or more computed points is needed ti interpolate in between the cross section library values. This variation of the cross section typically is treated with the equation 3.49 where the resulting cross section has a contribution of five state parameters, which are: αControl rod insertion; Tf Fuel temperature; Tm Moderator temperature; Dm Moderator density and Sb Soluble boron concentration. This treatment is known as Partial Derivatives Model. Beside the Partial Derivatives Model there are other cross sections formalisms. Multi-dimensional Table: Piece-Wise lineal interpolation is another method based in the following equation 3.57: Σ(Tf, T m, Dm, Sb, α) = (1 −α)Σunrod(Tf, Tm, Dm, Sb, α)+ +αΣrod(Tf, Tm, Dm, Sb) (3.57) Where the indexes rod and unrod are referring to the computed values with and without control rods. If the desired XS value subindex iin the equation 3.58 is a non-existing pre-calculated value, its value will be obtained by the interpolation between Tfaand Tfbpre-calculated points in Multi-dimensional Table: Piece-Wise lineal interpolation method. Σi(Tf, T m, Dm, Sb) = T f −T fb Tfa−TfbΣi(Tfa, Tm, Dm, Sb)+ +Tfa−Tf Tfa−TfbΣi(Tfb, Tm, Dm, Sb) (3.58) With this method first a generation of the base XS with rodded and unrodded features is needed. Then different XS are calculated for different parameters as it can be seen in figure 3.4. When the required XS is in-between the calculated ones, a lineal interpolation, equation 3.58 is made. Chapter 3. Codes and Models 72 Figure 3.4: Multi-dimensional table XS treatment scheme. Different codes have different approximations and different ways of treating the cross sections, when the desired value is not a pre-computed one. The accuracy of the method used will depend on the number of contributors to the final cross section value and to the type of interpolation used when finding the desired value. other codes might use another methodology named Multiple tables method. Where the resulting cross section is a combination of a different parameters including some history effects of the fuel. GenPMAXS [30] method is a combination of the previous mentioned methods in order to be enough accurate without loosing a lot of time in terms of computational time costs. In GenPMAXS [30] method the cross section will have a contribution of the three main factors such: State variables; History variables; Neighboring contribution all of these factors are considered (with and without control rod). Depletion is also considered in GenPMAXS [30] method. Depletion capacity led burn the fuel assemblies and move along the core reactor life cycle. To burn the different fuel assemblies, historical data from each region is required, once the burnup step is achieved, the historical data will be used to recalculate the new power distribution. Throughout this process, TH feedback is constantly given at each time step. PMAXS [30] is structured in a macroscopic cross section format to be read by the PARCS [4, 5] depletion routine. The code is structured in state variables, such as: •The control rod poison (CR) •Density of Coolant (DC) Chapter 3. Codes and Models 73 •Soluble Poison concentration in Coolant (PC) •Temperature of Fuel (TF) •Temperature of Coolant (TC) •Impurity of Coolant (IC) •Moderator parameters such density, soluble poison concentration, temperature and impurity (DM, TM, PM, IM) The cross sections are also functions of burnup (B) and other history variables. The 5 history variables are: •Control rod history (HCR) •Coolant density history (HDC) •Coolant soluble poison history (HPC) •Fuel temperature history (HTF) •Coolant temperature history (HTC) The history variables together with the fuel burnup determines the history state, H= [h1, ..., hnh, B], where nh is number of history variables used in PMAXS [30]. Cross sections can also vary depending on the conditions of the neighboring assemblies. Because the absorption cross sections of Xenon and Samarium are considerably larger than for other isotopes and are strongly dependent on the flux level of each node, the absorption cross sections for Xenon and Samarium are represented by their microscopic cross sections and number densities. The representation of the macroscopic cross section at a certain state is given by: Σl(C, S, N, H) = ΣE,l(C, S, N, H)+Nl Xeσl Xe(C, S, N, H)+Nl Smσl Sm(C, S, N, H) (3.59) The above macroscopic cross section expression can be used for absorption, fission, transport and scattering respectively. But does not include Xenon and Samarium and that is why they are added in the two extra terms on the right side of the equation. Superscript land the various subscripts denote the node index and isotope name respectively. C, S, N and H are 4 sets of state variables. Cconcerns about the insertion fraction for each node and is provided for each control rod composition, C= [c1, ..., cNc]. If the Chapter 3. Codes and Models 74 ciis flux weighted then the control rod effect is non-linear. If flux weighting is not available or not used, then the control rod effect is linear. PARCS [4, 5] provides for flux weighting by solving a“3-node” problem with very fine mesh for the node with the partially inserted control rod. This is the standard method for treating the “control rod cusping” effect. Alternatively, for special cases where the standard method is not possible, flux weighting is provided by using several branches for the control rod compositions in which different rod fractions are used to represent the non-linear effect of a partially inserted control rod in a node. In this treatment, PMAXS [30] would contain branches for 0 < c ≤1 for the first control rod composition and for other control rod compositions if necessary. Srepresents the state variables of the current nodes except for the control rod state. Those state variables are: Density of Coolant (DC); Soluble Poison concentration in Coolant (PC); Temperature of fuel (TF); Temperature of Coolant (TC); Impurity of Coolant (IC); Moderator parameters such density, soluble poison concentration, temperature and impurity (DM, TM, PM, IM). The Nindex contains information for 4 pairs of neighboring assemblies in the plane. Hrepresents the history state which contains history information, these are: Control rod history (HCR); Coolant density history (HDC); Coolant soluble poison history (HPC); Fuel temperature history (HTF) and Coolant temperature history (HTC). The macroscopic cross sections in PARCS [4, 5] are constructed with the assumption of a linear superposition of the partial cross sections on a base reference state. Such structure is taken into account when computing the cross section library with HELIOS-1.9 [27–29]. Next equation is showing how the code is representing the macroscopic cross-sections: Σl(C, S, N, H) = c01Σ( C1 c01 , S, N, H) + Nc ∑ i=2 ciΣ(ci, S, N, H) (3.60) c01 represents the sum of the unrodded fraction in one node and the first composition fraction in the same node. c01 = 1 − Nc ∑ i=2 ci(3.61) Σ(Cr, S, N, H) = Σr(H) + Cr ∂Σ ∂Cr(Cr/2,H) + NS ∑ j=2 ∆Sj ∂Σ ∂Sj(Cr,Sm j,Nr,H) + + 4 ∑ j=1 (Nn ∑ k=1 nj,k ∂Σ ∂nk(Cr,S,Nm j,k,H))(3.62) Chapter 3. Codes and Models 75 Here the reference cross section will receive contributions from the control rod insertion in different nodal locations; This is the second term on the right side, from the previous equation 3.62. The next term on the right side, considers the contributions of the other independent variables such density of the coolant, soluble poison concentration in coolant, temperature of fuel, and the temperature and impurity of coolant. The last term considers the four neighboring nodes to the computed one. This is where the state variables for the neighboring nodes contribute to the modification each cross-section. As can be seen from the previous equation, each computed cross-section at each node, will be a contribution from a reference state plus control rod insertion; plus independent state variables and finally plus the neighboring nodes independent state variables. Historical variables are taken into account for each of the terms; note the H index at each term from the previous equation. The partial derivatives of cross sections are calculated at the midpoint in between the reference state and the actual state for the current node. These partial are obtained by a piecewise interpolation of the pre-tabulated data using a “tree structure”. The variables which can modify the reference cross sections are distributed in three large groups: a) The control rod fractions, b) Variables of the current node and c) Variables of the neighbors. Each group is treated in the following order: a) First a check for the control rod fraction, b) Second to account for the independent variables, c) Third to account for the neighboring independent variables, and d) Finally consideration of the historical contribution. Such order is followed as a response to a study, performed by the code authors in order to check for the impact relevance to the reference cross-sections by the previous mentioned groups. Once the methodology of cross sections correction contributions and computing is described, next is to describe “tree structure”. Tree structure is the methodology used in order to store multiple cross sections sets in a proper manner and to be consistent with the interpolation methodology used by the code to finally obtain the desired cross section. At this level PARCS v3.0 [4, 5] reads the branch information provided by GenPMAXS [30] according to the branches computed by HELIOS-1.9 [27–29] and constructs a tree structure, where all branches are present and the partial derivatives are computed, in between each reference pre-computed state. In order to explain the tree structure methodology, first logical step is to describe the branch structure. Branches are used to compute and store information from every state. Essentially at each branch case, the same cross section as the reference state is used, but with at least one parameter modification. This modification gives a different value to the cross section and constitutes a branch. From all the proposed modification ranges for the state parameters, there will be one base branch and then the subsequent branches which compute the cross sections over all ranges for the other variable state parameters. For every branch there is a single modification to each of the previous parameters. The difference between the reference branch and the modified branch is computed and stored Chapter 3. Codes and Models 76 through a partial derivative. Those partial derivatives are the midpoint between branch state and its state. The way to compute them is through following equation: ∂Σ ∂xkXm =Σ(Xi)−Σ(XB(i)) xi k−xr k (3.63) where: •Xi= (xi 1, xi 2, xi 3, ...xi n) represents the state variables for one state. •XB(i)represents the variables for the base state from the same previous state. •Xmrepresents the variables which are located at the midpoint between the base state and each state constituting the branch. •xkare the branch variables for each state. Figure 3.5 show a common scheme of the GenPMAXS [30] tree structure organization. The different dependence over the historical variables, instantaneous variables and burnup is structured in three levels in this example. Two main branches are clearly identified after the historical dependence level. These two branches are describing with and without control rod states. After this level, several modifications of the other state variables are made. In this particular case density of the moderator, boron concentration and fuel temperature. There is two options for each before mentioned variable state. Finally each branch is computed with a collection of different burnup points. This is a very simplified scheme but useful to illustrate the tree structure scheme. once the scheme is clear is time to illustrate the cross section computation methodology, let’s take an example with six states, that means there are six cross sections provided by the lattice physics code, HELIOS-1.9 [27–29] used in the present study. Those states can be “Ref”, for the unrodded reference state Cr1, which represents the rodded version of the reference state. Then there are two computations for the unrodded side without rods, where the coolant temperatures, TC1and TC2, are modified. Finally two additional temperature modifications, TC3and TC4, in the rodded branch, are made. Figure 3.6 is the representation of the above mentioned structure. In this case Reference cross sections are stored in “Ref” state. The Cr1branch stores the difference between the reference state and the rodded state, due the control rod insertion. The cases with modified coolant temperatures are placed in each branch consequently to their modifications and taking into account the control rod insertion. Finally the vertical “T C1T,TC2T, TC3Tand TC4Trepresented by the partial derivative computation are placed in the Chapter 3. Codes and Models 77 Figure 3.5: Example of a tree structure scheme. Figure 3.6: Branch structure example scheme. Chapter 3. Codes and Models 78 midpoint between the reference state and each modified state. Those partial derivatives are computed using the previous equation 3.63. Once the branch structure is build, the branches are listed sequentially in the input file. Once PARCS v3.0 [4, 5] has read all the cases it will automatically build the tree structure for the partial derivatives. It is remarkable that there is no need for symmetry in terms of the points at each branch. The misaligned points in figure 3.6, are drawn on purpose to show this fact remarkable. That fact makes the whole thing change the name from branch structure to tree structure, since not all the branches are necessarily forming a regular grid. The example in figure 3.7, explains the steps performed by the code in order to compute the desired cross section. In this case the point 1 is the place (in terms of variable states) where the cross section is required. Note this point is a partially rodded cross section because it is in between the two lines. Σ(c, TC) = Σr+c∂Σ ∂Cr(Cr1/2) + (TC −TCr)∂Σ ∂TC (c,TCT) (3.64) The above equation is the one used to compute the cross section in the desired state 1. The left side term is the desired cross section. First term on the right side from the previous equation is the cross section at the reference state. Second term on the right side, is the control rod contribution. (Note sub index 2 which indicates the position in terms of the control rod insertion where the cross section must be computed.) The last term on the right side is the contribution due from the coolant temperature term. Since first term is has already been computed, the second and third terms still need to be computed. Second term is computed as difference between reference state and control rod state. The results of this computation are point number 2. The key point of the process is to compute point number 3. Essentially, the partial derivative is obtained by a linear interpolation between the four surrounding partial derivatives with respect the coolant .temperatures from the two branches, TC1T,TC2T,TC3Tand TC4T, as shown in equation 3.65. ∂Σ ∂TC (c,TCT) =w1 ∂Σ ∂TC (0,TC1T) +w2 ∂Σ ∂TC (0,TC2T) + +w3 ∂Σ ∂TC (1,TC3T) +w4 ∂Σ ∂TC (1,T C4T)(3.65) Weights for the four points are determined by linear interpolation using following equations: Chapter 3. Codes and Models 79 w1= (1 −c)TC −TC2 TC1−TC2 w2= (1 −c)(1−TC −TC2 TC1−TC2) w3=cTC −TC4 TC3−TC4 w4=c(1 −TC −TC4 TC3−TC4 ) (3.66) Figure 3.7: Cross section computation example with tree structure. Within GenPMAXS [30] environment, all the history variables, except burnup, are treated with partial derivatives with respect to those history variables. The same type of computations as explained for a regular state variables are used for history variables. In a picture representation of this case a multiple layers represent the different tree structures for the different history variables considered. The code will perform the linear interpolations described above, to obtain the cross section in at a determined point. Burnup dependence of the cross section is treated with a piece wise linear interpolation. The following equation is shows the form of this piece wise linear interpolation. Σi(Hj, B) = Bi,j k−B Bi,j k−Bi,j k−1 Σi(Hj, Bi,j k−1) + B−Bi,j k−1 Bi,j k−Bi,j k−1 Σi(Hj, Bi,j k) (3.67) where: •Σirepresents the cross section data in ith branch •Hjis history state of jth history case •Bi,j kis first burnup point in ith branch of jth history case which be greater than B Using these features, the code is able to interpolate and generate a cross section between the pre-computed cross sections which constitute the initial tree structure. Special requirements are needed to make the HELIOS-1.9 [27–29] output file readable by Chapter 3. Codes and Models 80 GenPMAXS v5.0. [30]. Some specifications in the ZENITH-1.9 [27–29] output process are required, other ways the computed cross sections are useless. GenPMAXS [30] reads all the reference cross sections computed by HELIOS-1.9 [27–29] and generates the partial derivatives between the reference states. Generation of the partial derivatives is the first step to needed to generate the GenPMAXS [30] output file. Essentially, the HELIOS-1.9 [27–29] output file has to written in a specific format so that GenPMAXS [30] can read all the characters on the ASCII output file. This step is quite tedious and requires a lot of trial and error methodology since the manuals are not very clear in this section. There is list of required keywords that need appear on the ZENITH-1.9 [27–29] output file so GenPMAXS [30] can identify the following data and process it. Such list of keywords is showed in following tables 3.1 and 3.2. Chapter 3. Codes and Models 87 3.1.5 SNAP SNAP [60] (the Symbolic Nuclear Analysis Package) its an interface created by NRC. Small definition of its capabilities will be given in this chapter. Since the aim of thesis is to use this software as a working tool, this chapter will not go in deeply detail on the features of the software. In here is intend to show the tool and to illustrate its basic features. More information can be found in [60]. SNAP [60] interface intends to facilitate the task of performing Nuclear analysis with system codes. Those analysis can go from the most common thermal hydraulic analysis (using TRACE [1–3] or RELAP5 [41–49] codes) to the more complex analysis with BEPU metrologies involved. Almost all the nuclear codes from NRC are suitable to be use by SNAP [60] platform. SNAP v 2.2.1 [60] was the version used in the present study. TRACE [1–3], RELAP5 [41– 49], FRAPCON, FRAPTRAN, PARCS [4, 5], DAKOTA [6–9], SCALE, MELCOR and CONTAIN are the supported codes for this version. In one side SNAP [60] platform is able to read existing input decks form each one of the above mentioned list of codes. On the other side SNAP [60] is able to create input files from scratch from each one of the above mentioned codes. By using SNAP [60] the user takes the advantage of moving from a ASCII input file to a more comfortable SNAP [60] template. Such template is trying to represent with shapes, colors and figures, whatever is intend to be simulated inside the ASCII file. So if we take as and example TRACE [1–3] or RELAP5 [41–49] which are thermal hydraulic system codes, when using SNAP [60], what is shown to the user is a scheme of pipes, pumps, valves, etc. . . that represent the system modeled. Figure 3.9 show the typical view of TRACE [1–3] input on SNAP [60] platform. On this common view from figure 3.9, the screen is divided in four windows. As usual on the top there is a menu where common functions from every computer program are included. Starting on the upper left there is a window that is giving the options from the different parts of the input deck (TRACE [1–3] in this case) Inside there are several menus that are grouping the different features form an input file. In here (keep in talking on TH system code area) the user will find general specifications on time steps, TH components, Control system components, Heat structures, connections between different elements. . . Second window, left bottom, is showing the inner menu on the above selected option. Lets say the user is looking at one pipe, selected on the above menu, the specifications like dimensions, flow areas, orientations. . . will be placed here. As a good capability at the end of each parameter to be defined there is a question mark. If the user clicks on this question mark, information from the user’s manual (linked previously to SNAP [60]) appears on the screen. This information is saving a lot of time when building a new input file, specially for the beginners. On the right side upper window a representation of the input deck is shown, typically this representation Chapter 3. Codes and Models 88 Figure 3.9: Example of SNAP window appearance. is trying to preserve dimensions in the picture it self, also is showing the connection between the elements. This window can hold multiple tabs, for instance in the case of TH system codes, one of the tabs may contain all the TH elements, the next one all the Control system, finally there will be another one which may contain the Heat Structures system. Last window, right en below, its giving some information about the job performed by the user, any errors, misleads, malfunctions etc. . . are shown here. SNAP [60] is using a modular plug-in design. This capability is structured in different plug-in connectors between SNAP [60] and the list of available codes for SNAP [60] platform. Once the user have the code and the SNAP [60] platform, the correct plug-in connector will allow the SNAP [60] to read, write, import and export files form a specific code. Beyond those capabilities of working with inputs, SNAP [60] platform is also able to perform restart input files. SNAP [60] allows to launch a calculation by using its own tool called calculation server. This tool is a linkage between SNAP created input and the executable file from each code. This calculation server will take the SNAP [60] developed input file, convert it into ASCII file with all the selected code specifications and it will launch it against the exactable file at the specific folder. This capability becomes very comfortable specially when running BEPU analysis, since a minimum of 59 cases are required. In that case SNAP [60] is taking care of preparing the 59 inputs and executing them in separate folders. Chapter 3. Codes and Models 89 SNAP [60] tool is continuously in development that is why newer versions are coming every couple of months. Meanwhile the supported codes are improving and some bugs from SNAP [60] old versions are reported, new versions with new features are released. In that sense the user needs to be aware on the selected version, and it’s clearly impossible to catch up with the newer versions specifications and modifications. That is why in the present study, one version was selected and all the work was performed with that “user frozen” version. Nevertheless a minimum effort will be required at he present to update all the job done to the latest version. There are good and bad consequences of using SNAP [60] tool, as mentioned above good consequences are enormous, in the sense that the user by using SNAP [60], is getting the whole picture of the problem very easy and can self-learn a lot about the used code, just by building its own input file. As a bad consequence the user is loosing a little bit track of the input deck and a lot of selections come by default (typically the user is not paying much attention on those) so in that sense the ASCII user was getting more knowledge since the beginning. Nevertheless SNAP [60] is a wonderful tool that is been used world wide and it has contrasted reliability. 3.1.6 DAKOTA The DAKOTA [6–9] (Design Analysis Kit for Optimization and Terascale Applications) is the selected code for the Uncertainty propagation with the coupled 3D NK-TH calculations required in the present study. DAKOTA [6–9] an internal research and development activity at Sandia National Laboratories in Albuquerque, New Mexico. A primary goal for DAKOTA [6–9](Design Analysis Kit for Optimization and Terascale Applications) development is to provide a systematic and rapid means to obtain improved or optimal designs or understand sensitivity or uncertainty using simulation-based models. These capabilities generally lead to improved designs and system performance in earlier design stages, alleviating dependence on physical prototypes and testing, shortening design cycles, and reducing product development costs. DAKOTA [6–9] code is a toolkit which provides a flexible and extensible interface between simulation codes and iterative analysis methods. DAKOTA [6–9] contains: Algorithms for optimization with gradient and non-gradient-based methods; Uncertainty quantification with Sampling, Reliability, and Stochastic expansion methods; Parameter estimation with nonlinear least squares methods; and Sensitivity variance analysis with design of experiments and parameter study methods. These capabilities may be used on their own or as components within advanced strategies such as surrogate-based optimization, mixed integer nonlinear programming, or optimization under uncertainty. By employing object-oriented design to implement abstractions of the key components required for iterative systems analysis, Chapter 3. Codes and Models 90 the DAKOTA [6–9] toolkit provides a flexible and extensible problem-solving environment for design and performance analysis of computational models on high performance computers. This chapter about DAKOTA [6–9] It is not intended to be as a comprehensive theoretical treatment. Rather, this section is intended to summarize a set of DAKOTA-related [6–9] capabilities and functions over the uncertainty quantification and optimization. General flow diagram from DAKOTA [6–9] is shown in figure 3.10 . Figure 3.10: DAKOTA flow information chart. DAKOTA [6–9] it is constituted with a big variety of iterative methods and strategies. It also has a lot of flexibility in order to interface with almost any simulation code. The following list explains about the variety of the DAKOTA [6–9] algorithms which compose the code: •Parametric Studies. Parameter studies employ deterministic designs to explore the effect of parametric changes within simulation models, yielding one form of sensitivity analysis. •Design of Experiments. Design and analysis of computer experiments techniques are often used to explore the parameter space of an engineering design problem, for example to perform global sensitivity analysis. •Uncertainty Quantification. Uncertainty quantification methods (also referred to as nondeterministic analysis methods) compute probabilistic information about response functions based on simulations performed according to specified input parameter probability distributions.This feature is the one used from AKOTA [6– 9] in the present study. Chapter 3. Codes and Models 91 •Optimization. Optimization solvers built in order to minimize cost or maximize system performance, as predicted by the simulation model, subject to constraints on input variables or secondary simulation responses. •Calibration. Calibration algorithms are orientated in the way to maximize agreement between simulation outputs and experimental data. They are used solve inverse problems. As it has been mentioned, all the above features the selected one, which fits the needs for the present study, is the one about the Uncertainty Quantification. At a high level, uncertainty quantification or also known as nondeterministic analysis is the process of characterizing input uncertainties, forward propagating these uncertainties through a computational model, and performing statistical or interval assessments on the resulting responses. This process determines the effect of uncertainties and assumptions on model outputs or results. In DAKOTA [6–9], uncertainty quantification methods specifically focus on the forward propagation part of the process, where probabilistic or interval information on parametric inputs are mapped through the computational model to assess statistics or intervals on outputs. The aleatory Uncertainty Quantification methods in DAKOTA [6–9] include various sampling-based approaches, following list enumerates all the supported sampling-based approaches. •Latin Hypercube Sampling. In here Monte Carlo (random) sampling and Latin Hypercube sampling methods are supported. •Reliability Methods. This algorithm includes both global and local reliability methods. Global reliability methods are designed to handle non-smooth and multimodal failure surfaces, by creating global approximations based on Gaussian process models. Local methods include 1st and 2nd order of the Mean value method and most probable point method. Also include first and second order of advanced mean value method. •Stochastic Expansion Methods. Rather than estimating point probabilities, stochastic expansion methods form an approximation to the functional relationship between response functions and their random inputs. •Importance Sampling. This method method allows the user to estimate statistical quantities such as failure probabilities in a way that is more efficient than Monte Carlo sampling. •Adaptive Sampling. The idea of the adaptive sampling is to construct a surrogate model that can be used as an accurate predictor of an expensive simulation. Chapter 3. Codes and Models 92 •Interval Analysis. Interval analysis is often used to model epistemic uncertainty. In interval analysis, one assumes that nothing is known about an epistemic uncertain variable except that its value lies somewhere within an interval. •Dempster-Shafer Theory of Evidence. The objective of Evidence theory is to model the effects of epistemic uncertainties. Epistemic uncertainty refers to the situation where one does not know enough to specify a probability distribution on a variable. •Bayesian Calibration. In Bayesian calibration, uncertain input parameters are described by a prior distribution. The priors are updated with experimental data, in a Bayesian framework that involves the experimental data and a likelihood function which describes how well each parameter value is supported by the data. Several actions need to be taken into account when performing a uncertainty quantification analysis. Within DAKOTA [6–9] framework several options are available. Typically the choice of uncertainty quantification method depends on how the input uncertainty is characterized, the computational budget, and the desired output accuracy. Some user guidelines within DAKOTA [6–9] capabilities are shown in figure 3.11 . Figure 3.11: Guidelines for Uncertainty Qualification method selection. DAKOTA [6–9] is also coupled to SNAP [60] platform. This coupling is made due a plug-in communicator. This way make the things easier for the user in terms of sampling and also in terms of input construction. The user needs to selects the desired quantities which are going to be perturbated with uncertainties and apply over them the PDF’s, the mean value the standard deviation and the maximum and minimum in case there is any. Also the desired level of confidence and probability is required. By selecting these quantities it will determine the number of cases to be executed. Once this is done Chapter 3. Codes and Models 93 DAKOTA [6–9] is going to apply the input values over the selected sampling algorithm and it will bring as a output a list of a modified inputs to be executed by the code or codes if there is a coupled calculation. Also some statistics information about the parameters and their deviations is printed at the end of the runs. Since the framework of the present study is quite global within the safety analysis, no more details about DAKOTA [6–9] tool will be given in this section. DAKOTA [6–9] has been used as a tool in order to propagate the uncertainties over the 3D NK-TH coupled analysis, its improvement or detailed behavior and algorithms goes beyond the aim of the present study. 3.2 Model description Nuclear system coupled calculations, involve a minimum of two codes: a Neutron Kinetics (NK) code, for simulating the core behavior, and a Thermal Hydraulics (TH) code for the coolant system modeling. As it has been mentioned, PARCS [4, 5] and TRACE [1–3] are respectively the codes chosen in the present study for representing each model. This section describes the NK and TH models developed for this study plus the coupling assumptions made in the coupled model. This section also describes the methodology learned and used and tagged as a “Know How” building a cross-section library. All the required steps are described in deeply details and it leads the reader of how to perform a collection of lattice physics calculations which will lead to development of a “whole cycle” cross section library. This paricular sub-section constitutes one of the big achievements of the present report. When building these models different assumptions were taken according to the experience gained on the participation to the Chapter 2 mentioned Benchmarks but also some bibliography of similar works performed previously was consulted see reference [64]. 3.2.1 Thermal hydraulic model Asc´o NPP is a 3 loops PWR with 2900 MW at full power. The TRACE [1–3] model completely reproduces the whole NPP system. TRACE V5 patch2 [1–3] is the version of the code used in the present study. The model is been validated against a 50% loss of load transient, typically used at UPC to validate full plant models [36–40], since there is existing plant data from such transient. Also the mentioned models are been used in several fields of thermal-hydraulic research area of study for the GET group such the work performed in scaling field, see: [65, 66]. In a coupled 3D NK-TH code calculation, the most relevant part of the thermal hydraulic model is the vessel. A 3D vessel component model in TRACE [1–3] has been implemented for the present study. Chapter 3. Codes and Models 94 There are different types of vessel models that could be used for the present study. Figure 3.12 shows several different approaches that can be used to model the vessel. On the left side a representation of a parallel channel vessel is shown, in the center a regular single pipe model commonly used for thermal hydraulic analysis is shown and finally on the left side the model is a combination between parallel channels and 3D volumes. The variety of used vessel model used depends on the target of the study. The variety of models can go from 1D vessel, with only 1D volumes, going over to a pseudo 3D model, which is a combination of a parallel channels plus surrounding volumes to active core, to finally real 3D vessel component, which was the one selected for the present study. For a 3D NK-TH coupled calculation with MSLB scenario with high asymmetry in core parameters during the transient a full 3D vessel component was selected as the best option to reproduce with high accuracy and quality the NK-TH feedback and the return to critically event during the late phase of the selected scenario. Figure 3.12: Different vessel model types. The chosen vessel model has 15 axial layers, 6 azimuthal sectors and 5 radial rings. The three lower axial nodes represent the lower plenum. The next six axial nodes represent the active core, the center region, the down comer and the bypass for the external regions. The top layers describe the upper head and the upper plenum of the vessel. Figure 3.13 shows an axial cut of the vessel representation. The axial core region (lighter area) is subdivided radially for each layer in 18 TH cells formed by overlapping, three rings and six sectors. As a result there are eighteen TH cells for each axial layer in the active core. The outer rings represent the down comer and the bypass along the active core height. Below and above the active core region, the thirty TH cells formed by overlapping the azimuthal sectors and radial rings have a different meaning as mentioned above. The height of each axial node in the active core is 0.609 m. The total active core axial height is 3.654m. In terms of the thermal hydraulic model, the core region consists of 6 axial nodes and 18 radial cells (nodes) at each axial layer (node). It is important also to note that the real core has Cartesian geometry, due the fuel assemblies, but the used 3D Chapter 3. Codes and Models 95 vessel component has cylindrical geometry, that is why some assumptions were taken when the mapping input decks where developed. Figure 3.13: Used Vessel component scheme. The rest of the 1D plant model remains the same as for a non-coupled system calculation. The three loops plus the pressurizer are included in the primary circuit representation. Three main steam lines are modeled in the secondary circuit representation. The Main Feed Water and Auxiliary Feed Water systems are also modeled for each loop. In terms of the safety injection systems, there are three accumulators, three LPIS and three HPIS systems, with later six modeled with FILL components. Finally a huge control block system (more than 1400 components) is included based on the developed UPC RELAP5 [41–49] Asc´o NPP model, which has been validated and used for more than 20 years for Asc´o NPP calculations [36–40]. The aim of the control block system is to reproduce with accuracy the plant response to different transients. The control block system has been increased with the addition of the control rod position control block. Since the model has the 3D capability and the validation of the model has been performed with a loss of load transient, where the control rod position is setup as a response of the thermal hydraulic parameters which pass the information and the position to the core simulator code which places the control rod at the proper position, according the signal coming from the thermal hydraulic code. This feature was not available in the releases NRC version thus some code modifications where need in order to achieve such capability. Code modifications made are presented in next Chapter of the present study in the model validation sub-section. Table 3.3 shows the TH model specifications in terms of the quantities of the used components. Chapter 3. Codes and Models 96 Table 3.3: TH model specifications Component Quantity Fills 7 Breaks 16 Pipes 94 Pumps 3 Separators 3 Single junctions 15 Valves 19 Vessels 1 Control systems 1455 Heat structures 175 Power components 162 3.2.2 Neutron kinetics model As mentioned in the previous section, the Asc´o NPP is a three loop PWR with 2900MW of thermal power. The core is modeled neutronically with PARCS v3.0 code. There are a total of 157 fuel assemblies in the core with a 17x17 pin array for each fuel assembly. The detail of modeling is one node per assembly in radial plane which results in 157 radial fuel nodes, plus 64 radial reflector nodes, which gives a total of 221 radial nodes for each axial level. Axially the FA is divided in 24 + 2 nodes, 24 for the core active region and 2 for the bottom and top reflectors. The height of the neutronic nodes in the FA’s is varying with smaller nodes in the lower and upper regions and larger nodes in the central region. This modeling reproduces with greater accuracy the material and thus crosssection variation along the axial height. in that sense smaller nodes are introduced in the areas where the cross-section variation is larger, while larger nodes are introduced where cross section variation is smaller. There are 6 control rod banks. In terms of the crosssection there are 648 + 2 different compositions, which means 650 nodes where the cross section is evaluated, and they might give different feedback contribution to the thermalhydraulic nodes. The cross-section library has been generated with the lattice physics code HELIOS-1.9 [27–29] using IDN-Asc´o cycle 11, 12 and 13 [24–26] specifications, which are the technical reports coming from the NPP different cycles. Table 3.4 shows the general NK model specifications in terms of the quantities of the used components. Figure 3.14 represents a radial core assembly layout, in here a 27 different types of fuel assemblies can be observed, also the reflector position. Finally table 3.5 shows the core Chapter 3. Codes and Models 103 conform to the fuel assembly. With the HELIOS-1.9 [27–29] code each fuel assembly is a matrix of different cells. Figure 3.17 shows the typical matrix used for modeling one regular fuel assembly, in between comas, the different nomenclature of each cell can be observed. Figure 3.17: 17x17 HELIOS matrix used to model a regular Asc´o NPP fuel assembly. Different descriptions from the different types of fuel assemblies will be needed in order to completely model the core. After obtaining all necessary information from the nuclear power plant technical report, our core will contain four different types of fuel assemblies. Those are: •Norm FA: Regular fuel assemblies, no control rod in it, no gadolinium fuel pins, no instrumentation. Only different enrichment grades can be considered in here. •Rodded FA: Fuel assemblies which contain control rods in the case of the control rod insertion. •Gd FA: Fuel assemblies which contain a poison material such as Gd2O3. Such fuel poisoned pins can be placed in different positions across the fuel assembly matrix. This fact leads to a large variety of this type of fuel assembly. Chapter 3. Codes and Models 104 •Reflector: Reflector type contains the reflector elements placed at the periphery of the core. There is no difference between upper and lower reflectors, this means they are modeled with same features. Once all the different fuel assembly types are represented, the next step is to determine the computing ranges for the different core parameters. This includes determining the fuel temperature, moderator temperature, moderator density and boron concentration ranges. A reference state is declared at this point. Experience, user skills and scenario knowledge are helpful to determine the ranges above and below such reference state. Two important issues are required for such selection. The first is to ensure that the selected range covers at least all the situations to be reproduced. In our model we like to cover all situations from the lower temperatures and conditions found in MSLB scenarios to higher temperatures and conditions found in ATWS scenarios. Nominal conditions are placed in between such range. See figure 3.18. Figure 3.18: Purposed range to be covered with the cross sections library values. The next important issue, when fixing the computational ranges, is to consider when defining the parameter range is the spread the computation points. They should be far enough apart to minimize the computation points and close enough so that large Chapter 3. Codes and Models 105 interpolations between two computed nodes are not needed. Figure 3.7 shows the selected points used in order to compute the cross section library. As mentioned above the selection is made with the users experience and it can be useful for future cross sections libraries. The selection of the computational nodes is one of the most challenging and important parts from the creation of the cross section library. Some trial and error was needed before starting the real calculations in order to have a good balance between computational points (detail description) and computational time (time needed to obtain the cross section library). Once this selection is made, next consideration concerns about the power. The power density at each fuel assembly needs to be taken into account in order to compute the cross sections. 40.106 W/grU is the input power density for the fresh fuel elements, different power densities are considered for old fuel assemblies that came from other cycles. Thus a power increase was applied to the modeled core in previous cycles. Xenon treatment is accomplished by defining three states of Xenon, these states are the non-equilibrium state, equilibrium state and quasi-equilibrium state. Finally burnup steps are defined. To define the burnup steps, several considerations need to be taken. First a small burnup step is required to account for the xenon equilibrium time, in the case reported in this paper is up to 150 MWd/t, The following steps are differently spaced along the fuel assembly life. Again, user experience is the determinant when selecting where the stepwise and again considerations about the interpolation and the maximum allowed burnup are taken in each selection. Same kind of the above mentioned trial and error tests was used here before determining the burnup steps. The selected burnup steps are (0, 150, 500, 1000, 2000, 4000, 10000, 40000, 47000 and 57000) MWd/t. This chain was the one used for regular fresh fuel assemblies. Different burnup steps were taken on the wasted fuel assemblies coming from the previous cycles, due the power increase suffered in the Asc´o nuclear power plant. Chapter 3. Codes and Models 106 Figure 3.19: TREE structure scheme. Once reference state and the other states (burnup steps, temperature ranges, density ranges, boron concentrations, Xenon, power) are defined, the next step is to build a TREE operator. The TREE operator is a list of all the states to be computed. Figure 3.19 shows the initial part of the TREE structure. Each reference state is considered at every burnup step, and for each burnup step, there is a list of combinations concerning the parameter variations. This structure is prolonged for each reference state along all the burnup steps. Such operation has to ensure that all the possible combinations from all the above listings is taken into account. This ends up with a long list with approximately more than 1000 possible combinations of the previous states. Moreover, three main moderator reference temperatures are selected as a base case for this tree structure, Low, Medium, and High according to that three calculations are performed for each fuel assembly, with more than 1000 states computed inside each case. Every computed state is calculated over all burnup steps. More than one hundred days were need to compute all the cross section sets. Chapter 3. Codes and Models 107 Table 3.7: List of the computed points for each core parameter Fuel temperature (K) Nominal 898.15 TFu1 400.00 TFu2 1400.00 TFu3 1800.00 TFu4 2200.00 TFu5 2800.00 Moderator temperature (K) Nominal 579.65 Low 565.45 High 323.00 Cold 600.13 Tmod1 330.00 Tmod2 450.00 Tmod3 525.00 Tmod4 600.00 Moderator density (gr/cm3) Nominal 0.70137 Low 0.74264 High 0.66212 Cold 0.99804 Dmod1 0.01 Dmod2 0.3 Dmod3 0.55 Dmod4 0.65 Dmod5 0.75 Dmod6 1.00 Boron PPM (mg/kg) PPMN 0.001728 PPML 0.000475 PPMH 0.000533 PPM1 0.00001 PPM2 0.0022 Chapter 3. Codes and Models 108 AUROROA-1.9 [27–29] is the input processor subprogram, which inputs all specifications inside HELIOS-1.9 [27–29]. The HELIOS-1.9 [27–29] computing process generates a huge output file with more information than required for the creation of the cross section library. This information is saved in a .hrf file, which contains ASCII codifications and it is unreadable by regular text file reader programs. The user must build another input deck using ZENITH-1.9 [27–29], another HELIOS-1.9 [27–29] subprogram, in order to call all the requested values that will be used to conform the cross section library. Also the ZENITH-1.9 input deck is where all the homogenized cross sections are collapsed in two groups and over the material homogenization. As it has been mentioned in the present study, the common limit of 0.625 eV is being used in order to separate the fast and thermal groups. Also in the ZENITH-1.9 [27–29] processor a list of the output parameters that will appear in the output file must be declared. Figure 3.8 shows such list. Table 3.8: ZENITH-1.9 output parameters list Keyword Meaning D1 Diffusion coefficient group 1 D2 Diffusion coefficient group 2 SIGA1 Group 1 absorption macroscopic cross section SIGA2 Group 2 absorption macroscopic cross section SIGS Scattering 1 to 2 macroscopic cross section SIGF1 Group 1 fission macroscopic cross section SIGF2 Group 2 fission macroscopic cross section SIGNF1 Neutron produced by fission in group 1 SIGNF2 Neutron produced by fission in group 2 Flux1 Group 1 neutron flux Flux2 Group 2 neutron flux ADF1 Group 1 assembly discontinuity factor ADF2 Group 2 assembly discontinuity factor VELOC1 Group 1 absorption velocity VELOC2 Group 2 absorption velocity After extracting all necessary information from the ZENITH-1.9 [27–29] output file, all the files are saved in a proper manner. In the present study cycle 13 of the Asc´o NPP is being reproduced in a core configuration, where up to 27 different fuel assemblies were modeled. Some of the fuel assemblies were fresh fuel, some came from previous cycles, and one fuel assembly came from the first cycle of the reactor core. For each different fuel assembly, there will be up to a minimum of three calculations for each of Chapter 3. Codes and Models 109 the High, Medium and Low reference states. If the Fuel assembly contains control rods, this will add three additional calculations, for the High, Medium and Low reference state assembly configurations. There are 8 different kinds of fresh fuel assemblies, from 15A to 15H, with different enrichments considered in this study. As previously mentioned one fuel assembly came from the first cycle, one fuel assembly came from the 11th cycle, (number 11), one fuel assembly came from the 10th cycle (number 12) and 8 types of fuel assemblies came from the 12th cycle, (number 14). Tables 3.9 and 3.10 show all the calculations performed. Chapter 3. Codes and Models 110 Table 3.9: Fuel assemblies calculations for Asc´o NPP cycle 13 CYCLE Fuel assembly type Reference state Without CR With CR Calculations 13 15A(MAEF+IFM-ZR) Low x Medium x High x 15B(MAEF+IFM-ZR) Low x Medium x High x 15C(MAEF+IFM-ZR) Low x x Medium x x High x x 15D(MAEF+IFM-ZR) Low x Medium x High x 15E(MAEF+IFM-ZR) Low x x Medium x x High x x 15F(MAEF+IFM-ZR) Low x Medium x High x 15G(MAEF+IFM-ZR) Low x Medium x High x 15H(MAEF+IFM-ZR) Low x x Medium x x High x x 1 1(STD) Low x Medium x High x As it can be seen from tables 3.9 and 3.10, there are 25 different calculations, each with three reference states, plus reflector input decks which make close to 80 calculations to obtain all the necessary information required to build the cross section library for Cycle 13 of the Asc´o NPP using HELIOS-1.9 [27–29] and the following self-developed methodology for building a cross section library. Almost one day was needed for each computation with the used LINUX servers. The described methodology was self-developed for the present study. Once all the computations are finished, the output files will work as input files of GenPMAXS-v5.0 [30]. The idea for making all computations go through GenPMAXS-v5.0 [30] is to be able to provide all needed information to core simulator code PARCS-v3.0 [4, 5]. In this step, all information is collapsed into 27 different files representing the 27 different fuel assemblies in the core, plus 1 file which represents the reflector. First GenPMAXS-v5.0 [30] reads all the information from every ZENITH-1.9 Chapter 3. Codes and Models 111 Table 3.10: Fuel assemblies calculations for Asc´o NPP cycle 13 (Continuation) CYCLE Fuel assembly type Reference state Without CR With CR Calculations 10 12B(AEF+IFM) Low x Medium x High x 11 13(AEF+IFM) Low x Medium x High x 12 14A(AEF+IFM-ZR) Low x x Medium x x High x x 14B(AEF+IFM-ZR) Low x Medium x High x 14C(AEF+IFM-ZR) Low x x Medium x x High x x 14D(AEF+IFM-ZR) Low x x Medium x x High x x 14E(AEF+IFM-ZR) Low x Medium x High x 14F(AEF+IFM-ZR) Low x Medium x High x 14G(AEF+IFM-ZR) Low x Medium x High x 14H(AEF+IFM-ZR) Low x Medium x High x output file and converts it into a GenPMAXS-v5.0 [30] file. Second the GenPMAXSv5.0 [30] files can be merged. First the three reference state configurations are merged and then these files are merged with the rodded states. At the end there are 27 files, where some contain control rod calculations and others are without control rods calculations. The PARCS [4, 5] code structure will select a file to represent each different fuel assembly. In case of control rod insertion, nothing needs to be done, since all the fuel assemblies which might contain the control rods already hold the information to correct the cross sections according the control rod insertion. In that sense if there is some node which has a control rod inserted during a transient, the code will pick from the GenPMAXS [30] file the required correction to the control rod insertion and this Chapter 3. Codes and Models 112 effect will lead to a increase of absorption cross section and neutron flux reduction. Uncertainty deviations have not been considered when building the present cross section library. It is important to have in mind that for a complete Best Estimate Plus Uncertainty analysis in a coupled 3D NK-TH model, uncertainty considerations have to be taken into account since the creation of the cross section library, through over all the steps which compose a coupled calculation. In the present case of study Uncertainties have only been considered over the thermal hydraulic model parameters and neutronic core code parameters. OECD UAM project [23] project intends to carry out the uncertainty propagation over all the phases from lattice physics phase to coupled 3D TH-NK phase. Early phases of the project are still going on. General results (i.e. global Uncertainty propagation) will be discussed in coming years. UPC GET group is actively participating in this international project. 3.2.4 Coupled model The thermal hydraulic model has been modified in order to meet the neutronics model requirements. In the present model, there are 157 Heat Structures (HS) in the thermalhydraulic core region. The equivalence is one HS to one FA. There are 18 radial thermalhydraulic cells in each axial thermal-hydraulic active core layer. Figure 3.20 shows the assignment (mapping) in each axial layer of the active core and reflector TH cells to neutronics nodes. Every color area represents a thermal-hydraulic cell. In terms of the axial nodalization, there are a different number of nodes for each model. Due the axial cross-section variation, there are 24 non-equidistant axial nodes for the HS and for the neutronic models however there are 6 equidistant nodes for the hydraulic model. Consequently the axial mapping between the HS and each of the neutronics nodes is not one to one. Some sensitivity studies about this were made, when validating the model. Essentially a one to one axial distribution model was built and tested. These results are discussed in next Chapter. Notice, in the regular model, the neutronic nodalization is finer than the coarser thermal-hydraulic nodalization. In that sense, it should be taken into account that several neutronic nodes are receiving the same thermal-hydraulic information and vice versa. That is one thermal hydraulic node is receiving and averaging power information coming from different neutronic nodes. All information is contained in a mapping file. This file is responsible for the good agreement and exchange of information between the thermal-hydraulic code TRACE [1–3] and the reactor physics code PARCS [4, 5]. Chapter 4. Model Validation and Improvement 119 Table 4.2: Example of Asc´o NPP power vs. control rod position in steps Power Bank C Bank D (%) (withdrawn steps) (withdrawn steps) 0.0 113 0 5.0 114 0 10.0 125 0 13.0 131 0 15.0 136 4 20.0 147 16 25.0 158 27 30.0 168 38 35.0 179 49 40.0 190 60 45.0 201 71 50.0 212 83 55.0 223 94 55.7 225 95 60.0 225 105 65.0 225 116 70.0 225 127 75.0 225 138 80.0 225 149 85.0 225 161 90.0 225 172 95.0 225 183 100.0 225 194 Chapter 4. Model Validation and Improvement 120 Figure 4.1: Example of Asc´o NPP power vs. control rod bank position in steps. In the previous examples, figure 4.2 and figure 4.1 an actual representation from Asc´o NPP of the control rod bank overlapping movement which begins with the insertion of Bank D immediately after the power decreases below 1.0 of the nominal power. Note that Bank C remains withdrawn until the power is slightly below 0.6 of the nominal power. Banks A, B, SA and SB start to be inserted consecutively after banks D and C which are the first ones to be inserted. The same behavior can be seen in Table 6 where the total power is represented in terms of the total withdrawn steps. This feature cannot be captured with the card system SCRAM implemented in the PARCS [4, 5] code. Once the SCRAM card is activated all control rod banks get inserted into the core at the same time, just the previously marked as stuck control rod banks, remain completely withdrawn in the transient. As it has been mentioned before, within the framework of the present study, it is important to have a completely validated model in terms of Thermal hydraulic and Neutronic behavior. To achieve that purpose it has been considered a must, to solve the problem of the dynamic control rod movement between TRACE v5.0 patch2 [1–3] and PARCS 3.0. [4, 5]. What is intended to solve in this endeavor, it is a way to compute or to assign at every single control rod bank a step position against time. Such position must be computed in TRACE [1–3] code as an answer to certain thermal-hydraulic conditions such: •Manual stop I.S. manual Chapter 4. Model Validation and Improvement 121 •High flux. High set point RI •High flux. Low set point RP •High flux. High set point RP •High neutronic flux oscillation •C-3, OTDT •C-4, OPDT •Low mass flow one loop •Low mass flow 2/3 loops •High pressure in pressurizer •Low pressure in pressurizer •High level in pressurizer •Very low level in one steam generator •Automatic safety injection •Turbine trip The above list is the actual list from Asc´o NPP model which has been used in successful transient analysis and model validations over the last twenty years in Technical University of Catalonia. Some of the previous signals lead the plant to SCRAM status. After considering those signals, the IDN from Asc´o NPP [24–26] also needs to be checked, so the position of each control rod bank can be determined in terms of each thermal hydraulic condition above mentioned. With this new feature the coupled calculation there is full feedback, and every transient might be able to be reproduced without considering the control rod position in advance. Notice that a huge logic and control system was required in order to compare different signals and give as a result a final control rod position for each control rod bank. Both codes are coupled and compiled under the same executable file. Nevertheless an external MAPTAB, mapping file, is needed to assign the matrices between thermal hydraulic nodes and neutronic nodes. Also the weighting factors between both side nodes are fixed in this file. some special features are described in the heading part of the file. Such file is structured in several parts, Heading part; Cards part and Assignment part. Chapter 4. Model Validation and Improvement 122 ** PARCS Mapping for V3DcylHS and Cartesian - x[17], y[17], z[26] * * *d: Generated from PARCS model ssAsco2010.inp * * Doppler Feedback %DOPL LINC 0.7 * * SCRAM Trip *%TRIP *34 * %CRSIG 9111 1 9112 2 9113 3 9114 4 9115 5 9116 6 * * Reflector Properties %REFLPROP *ctemp ftemp cden cvoid ppm 584.65 898.15 701.37 0.0 1728.0 * * Volume Number Table %TABLE1 100 20 3 1 1.0 100 20 3 2 1.0 100 20 3 3 1.0 . . . . 100 23 10 5745 1.0 100 23 10 5746 1.0 * * Heatstructure Number Table %TABLE2 933 1 1 1 0.0 933 1 1 2 0.0 933 1 1 3 0.0 933 1 1 4 0.0 933 1 1 5 0.0 Chapter 4. Model Validation and Improvement 123 933 1 1 6 0.0 933 1 1 7 0.0 933 1 1 8 0.0 1001 1 1 9 0.0 1002 1 1 10 0.0 . . . . 1155 1 24 5736 0.0 1156 1 24 5737 0.0 1157 1 24 5738 0.0 933 1 24 5739 0.0 933 1 24 5740 0.0 933 1 24 5741 0.0 933 1 24 5742 0.0 933 1 24 5743 0.0 933 1 24 5744 0.0 933 1 24 5745 0.0 933 1 24 5746 0.0 * end of data Above lines constitute an example of the mentioned MAPTAB input file. The heading part is related to the specifications and title of the file. The cards part is related to the specific calculations for the Doppler effect, SCRAM trip from the thermal hydraulic code which will begin the SCRAM process into PARCS [4, 5] (Notice, this feature is being disabled with the leading asterisk, this is due the new dynamic control rod feature, will take SCRAM situation into account also. Nevertheless with the modified code, both features can perfectly coexist without any controversy). Next lines give some reflector properties and then there is the card %CRSIG, which is the one introduced in order to model the dynamic control rod movement capability. In the final part of the MAPTAB file, there are the assignment cards. These cards are structured in two parts, %TABLE1 cards, where assignment and weighting factors between thermal hydraulic volumes (Vessel nodes, notice vessel component number is 100) and neutronic nodes are given. Second part is %TABLE2, where assignments and weighting factors between the heat structures (notice HS start with number 1001 and number 933 is saved for the reflector) and neutronic nodes are given. Notice the neutronic nodes number ascends up to 5746, this equals 221 radial nodes times 26 axial neutronic levels. Chapter 4. Model Validation and Improvement 124 As it has been mentioned this new feature is been introduced into the code by using the %CRSIG card, which essentially gives in pairs a number of the control signal in TRACE [1–3] and a control rod bank assigned with such signal. Those control rod banks hold numerical names thus in our case banks A, B, C, D, SA and SB have been renamed 1, 2, 3, 4, 5 and 6 respectively. The source codes from TRACE v5.0 patch 2 [1–3] and PARCS 3.0 [4, 5] have been modified. In the following lines a detailed explanation of the all modifications and checks to the source code made are described. Subscript “Rai!” refers to some modified lines made in the source code. Other routines have also some modifications, but in the below description only TdmrErrorCheckM.f90 are showed. Adjustments made: Source code file: Pdmr mapM.f90. 1. nfields(line) modification to meet the number of fields that will be in %CRSIG card. Source code file: Pdmr initM.f90. 1. initcrp(i) divide at the end by ncrbstep. Source code file: Pdmr timeM.f90. 1. To check newcrp(i) in line 68. 2. To check crbpos(sgvbank(i)) in line 81. Source code file: TdmrInitCalcM.f90. 1. To check initcrp(i) and to check r8bufn(i), line 110. Source code file: TdmrTimeCalcM.f90. 1. To check newcrp(ii) in line 79. 2. To check initcrp(ii) in line 79. 3. To check csSig(jj)%presVal in line 81. 4. To check r8bufth(1+ii) in line 83. Source code file: TdmrErrorCheckM.f90. Chapter 4. Model Validation and Improvement 125 1. To check the following loops. •Check to see if new control rod bank positions are outside range. •Line 131 to line 146. ( ... ) WRITE (mesg(1), 2109) isgv(ii) 2109 FORMAT (“Processed Signal Variable #”, i5) CALL TDMRStatusMesg( mesg, dim=1, status=TDMRSTAT) IF (IABS(csSig(jj)%icn1) .NE. sgvbank(ii)) THEN !Rai CALL error(1, ’Fatal* wrong cnt. rod group ID’) !Rai RETURN !Rai END IF !Rai IF (initcrp(ii) .LT. 0.0D+00 .OR. initcrp(ii) .GT. 1.0D+00) THEN crcntl = .FALSE. WRITE (mesg(1), 2110) initcrp(ii) 2110 FORMAT (“Initital control rod bank position was out of accceptable range: ”, 1pe20.12) CALL TDMRStatusMesg( mesg, dim=1, status=TDMRWARN) END IF GO TO 21 21 CONTINUE ( ... ) Source code file: TdmrCommM.f90. 1. To check the structure value r8bufn = pbuf%gi2th( nbuf+1: nbuf+dimbuf(6)) line 212. Source code file: TransDriveM.f90. 1. To check the loop. •determine current crbank . Source code file: Pdmr commM.f90. Chapter 4. Model Validation and Improvement 126 1. To check SUBROUTINE pdmr comm copybufto(). PARCS [4, 5] input considerations and actions: 1. To add the card CR AXINFO third field ncrbstep, Control rod full insertion position from the bottom of the problem geometry, cm. 2. To consider that in PARCS [4, 5]: •0 steps, means completely inserted. •### steps number of the withdrawn steps. TRACE [1–3] input considerations and actions: 1. To create a function like: •idcb,icbn,icb1,icb2 •Such function has to change the number of the steps and normalize to 1. To consider that in TRACE [1–3]: –1.0 means completely inserted (This is 0 steps in PARCS [4, 5]). –0.0 means completely withdrawn (Maximum number of steps in PARCS [4, 5]). •To create a control signal in TRACE [1–3] who reads the above created function. –idsv = number which will go into MAPTAB file. –isvn = 16. –ilcn = (negative) function created before, where the control rod movement is inside, it goes from 0.0 to 1.0. –icn1 = Bunk number which will be moved under the previous parameters. ∗Important: It is not possible to modify this parameter inside the SNAP [60] platform, ASCII modification is required. MAPTAB file considerations and actions: 1. To add %CRSIG card. 2. To add in the under line of %CRSIG card the number of the control signal followed by each controlled control rod bank. They come in pairs until the last bank to be controlled. Chapter 4. Model Validation and Improvement 127 Final check: 1. To check with the MOVE BANK card. 2. MOVE BANK 1 0.0 225.0 20.0 225.0 25.0 125.0 30.0 100.0 50.0 50.0 3. Double check to be sure that same results appear when: •Signal coming from TRACE [1–3]. •MOVE BANK card from PARCS [4, 5]. 4. Same results in both cases. With all of those checks and modifications the new compiled code which holds, the capability of the dynamic control rod movement, has been developed. Moreover, the %SCRAM card in the MAPTAB file can be disabled, since the %CRSIG gives the position of every control rod bank at each time step, there is no need for the SCRAM feature coming from PARCS [4, 5]. In the new compiled version of the code, the SCRAM signal can be setup in the thermal-hydraulic side of the coupled calculation, as an answer of to the typical signal which leads to SCRAM, (such it is presented at beginning of the actual section) in a non-coupled calculation, replicating the signals in the list mentioned previously. Before going further, it is necessary to ensure that the modified version of the code is not overlapping any of the previous features. That is why the final check from the above mentioned steps is very important. Once the modifications in the code are solid, and work with any test scenario, it will be necessary to implement a control system into the thermal-hydraulic code. Such system will check the different parameters from the plant model and will assign the control rod bank step position for every time step and for each different control rod bank, depending on the status of the plant. Such control system also has to contain the SCRAM capability, in case of an eventual SCRAM event. At the end what are only by passed to PARCS [4, 5] are the steps for every control rod bank. At the end only the steps for each control rod bank are passed to PARCS [4, 5]. This control block system has been adapted from the one previously developed by the author for RELAP5-3D [31–35] code model. In the following figure and scheme of the control block system build for that porpoise can be seen. As it can be seen in figure 4.2, the control rod logic is quite complicated, because it involves different parameters and different actuation over different system signals. As a basic definition the control logic is looking at the core temperatures, core power and turbine power, to adjust the control rod position of the first inserted control rod bank (Bank D in this particular case) over the pre determined values of position coming from the Asc´o NPP specifications. Chapter 4. Model Validation and Improvement 128 Some response dead bands are also incorporated into this model in order to avoid some oscillations which might lead to some instabilities. In this way if there is some over power in the steady state conditions, the logic will tell to insert some positions of the control rod bank until the new balance is reached. Vice versa, if there is under full power situation, in steady state conditions, some steps will be required to be withdrawn in order to gain some power production. Control rod movement logic is also linked to the SCRAM signal which might come from different situations (They have been listed in previously). Once the position of the first control rod bank is determined, the other control rod banks position will be determined over the position of the first one. their insertion priority is also determined on the Asc´o NPP specifications. Control rod bank position is generated from 0 to 1 value in TRACE [1–3] and it is required in steps for PARCS [4, 5], some conversion from one side to the other “nomenclature” was also needed. With this final step working the coupled TRACE/PARCS [1–3] and [4, 5] model is ready to be validated. Figure 4.2: TRACE control rod bank position control system scheme. Next figures 4.3, 4.4, 4.5, 4.6, 4.7 and 4.8, show the validation process, between the 50% loss of load of plant data and the calculations from the new compiled version of the code. Typically the transient is computed as a restart file from a 12000 seconds calculation which has been carried out in order to obtain steady state conditions. A general good agreement is shown in all the figures. Even though the agreement is not 100% precise in some figures, what is important to achieve with our model, is the move from one stable region to another. This is in our case, steady state at full power at the beginning off the transient, to 50% of the full power at the end of the transient. Also looking at the