ML20235F200

Final Response to FOIA Request for Documents Ref in NUREG-1150.Forwards App C Document.Document Also Available in Pdr.App D Document Totally Withheld (Ref FOIA Exemption 5)

ML20235F200
Person / Time
Issue date:	09/23/1987
From:	Grimsley D NRC OFFICE OF ADMINISTRATION & RESOURCES MANAGEMENT (ARM)
To:	Sholly S MHB TECHNICAL ASSOCIATES
References
FOIA-87-416, RTR-NUREG-1150 NUDOCS 8709280449
Download: ML20235F200 (4)

Text

{{#Wiki_filter:_ _ _ _ _ _ - _ _ - _ _ _ _ 9 FOIA-87-416 ,f, ea s s% \\ REseO~st ryet i RESPONSE TO FREEDOM OF XI*"^' I l " ' ^" INFORMAT:ON ACT (FOIA) REQUEST g, SEP 2 31987 eseee DOCKET kuM8ERt$\\ W epphesanl REQUESTER Mr. Steven C. Sholly PART 1.-RECORDS RELEASED OR NOT LOCATED (See checked bones) No agency records subpct to the sugeost have been located. X No additional agency records subpct to the request have been located, Agency records subpect to the request tibet are identified in Appendix are already availab6e for public lospection and copying in the NRC Public Document Room, 1717 H Street, N.W., Westungton, DC. j b are being made avdaNe for pudic inspection and copying in the NRC Public Document Agency records subM to the request that are identifed in Appendix X Rt.om,1717 H Street, N.W., Washington, DC, in a folder under the FOIA number and requester name. The nonproprietary version of the propostl(s) that you agreed to accept in a telephone conversation with a rnember of my staff is now being made avalable for public inspectu>n and coying at the NRC Pubhc Document Room,1717 H Street, N W, Washington, DC, in a folder under the FOIA number and requester name. Enclosed is informaton on how you may obtain access to and the charges for copying records placed in the NRC Public Document Room,1717 H Street, N.W., Washington, DC. Agency records subpct to the rowest are enclosed. Any applicable charge for copes of the records provided and payment procedures are noted in the comments section. Records subpect to the request have been refened to another Federal agency (ies) for review and direct response to you. In view of NRC's response to the roovest, no further action is being taken on appeal letter dated PART lt.A-INFORMATION WITHHELD FROM PUBLIC DISCLOSURE Certain information in the requested records a bang withheid from public disclosure pursuant to the FOIA exemptions described in and for the reasons stated in Part II, sec-y tions B, C, and D. Any ro6eemed portons of the documents for which only part of the record is being withheld are being made available for public mspectici and copying in the NRC PubFc Document Room,1717 H Street, N.W., Washington, DC, in a folder under ttus FOIA number and requester name. Commsats 8709280449 870923 PDR FOIA SHDLLYB7-416 PDR d*1 i SIGN R D88tECTOR, OtVtsi O ULIS m CORDS hW. 4 NRC FORM 464 (Pan in to sei

07*4lb DATE;SEP 2 3 1987 FREEDiM OF INFORMATION ACT MSPONSE FOIA NUMBER (SI: PAMT 11 Q-APPLICABLE FOIA EXEMTTIONS t Records subject to the request that are described in the enclosed Appendices O are being withheld in their eritrety or in part under FOIA l Exemptions and for the reasons set forth below pursuant to 5 U.S.C. 552(b) and 10 CFR 9.5(a) of NRC Regulations. l l

1. The withheld informaton is property classifed pursuant to Executree Order 12356 (EXEMPTION 1) l

2. The withheld information relates solely to the internal personnel ndes and procedures of NRC. IEXEMPTION 2) i

3. The wttm.M information is specifically exegted from public dehre by statute indicated: (EXEMPTION 3)

Sectk>n 141144 of the Atomic Energy Act which prohibits the declosure of Restricted Data'or Formerty Restreted Data (42 U.S.C. 2161-2165). Section 147 of the Atomic Energy Act which prohibits the ^=< hare of Unclassified Safeguards Information 142 U.S.C. 21671.

4. The withheld information is e trade secret or commercial or fnsncial information that is being withheld for the reason (s) indicated: (EXEhfFTION 4)

The informaton is considered to be confidential business (propnetary) informaton. The information is considered to be propnetary information pursuant to 10 CFR 2.790(d)(11. The informaton was submitted and received in confidence from a foreign source pursuant to 10 CFR 2.790tdH2).

5. The wrthheld information consists of interagency or intraagency records that are not available through discovery dwing liigation. Oscionse of predecisional information X

wouid tend io innibit tne op.n and fr.nx escnange of ide.s e.sermai to ine deliberative process. where records are thheid in their om,ety. tne facis.re inestricabiy intertwmed with the oredeciseonel informaton. There also are no reennnably segregable factual portons because the release of the facts noidd permit an indirect inquiry into the predecesonal process of the agency. IEXEMPTION 5)

6. The withheld information is exempted from public disclosure because its disclosure would result in a clearly unwarranted invasion of perunaf privacy. (EXEMPTION 6)

7. The withheld information conssts of investigatory records comp 4ed for law enforcement purposes and is bemg withheld for the reason (si sidkcated. (EXEPPTION 7)

Disclosure woatJ interfere with an enforcement proceedmg because it could reveal the scope, direction, and focus of enforcement eMcr1s. and thus could possibly allow them to take action to shield potential wrongcbrig o' a volation of NRC requirements from invesugators. IEXEMPTIOh 7uu) Disclosure would constitute an unwarranted invasion of personal privacy (EXEMPTION 7(CH The informaton consists of names of individuals and other ' formaton the disclosure of which would reveal identroes of confidentel saurces. (EXEMPTION 7(D)) b o i PART ll.C-DENYING OFFICIALS Pursuant to 10 CFR 9.9 and 'or 9.15 of the U.S. Nuclear Regulatory Commason regulations, it has been determined that the information withhoid is mampe from production or disclosure, and that its production or disclosur6 is contrary to the public i iterest. The persons responsible for the denial are those officials identifed below as deyng officials and the Director, Division of Rules and Records. Office of Administration, for any dermals fut may be appealed to the Executive Director for Operations tEDOI. OENylNG OFFICIAL TITLE / OFFICE RECORDS DENIED APPELLATE OFFICIAL ~ Dr. Thomas E. Murley Director, Office Of Nuc' ear Appendix D X Regulatory Regulation PART 11 D-APPEAL RIGHTS Th) denial by each denying official identified in Part it.C may be appealed to the Appellate Official identified in that section. Any such appeal must be in writing and must be made within 30 days of receipt of this respoinse. Appeals must be addressed as appropriate to fhe Executive Dinsetor for Operations or to th3 Secretary of the Commission, U.S. Nuclear Regulatory Commission Washington, DC 20566, and shou 6d clearly state on the erwelope and in the letter thit it is an " Appeal frorn an Initial FOIA Decision." NRC FOAM es4 (Part a U.S. NUCLEAR REGULATORY COMMISSION sesi FOlA RESPONSE CONTINUATION

Re: FOI A-87-416 APPENDIX C RECORDS MAINTAINED'IN THE PDR UNDER THE ABOVE RE0 VEST NUMBER ' NUMBER' DATE DESCRIPTION ITEM X' '1. Undated Dr' aft Preliminary Report' for' Comment, "MELPROG-PWR/ MODO: A Mechanistic Code for Analysis of Reactor Core Melt Progression and Vessel Attack Under Severe Accident Conditions" by W. J. Camp, et al. (179pages) e.

Re: FOIA-87-416 APPENDIX D DOCUMENTS TOTALLY WITHHELD EXEMPTION 5 ITEM C 1. 8/1/85 Memorandum from T. P. Speis to H. R. Denton,

Subject:

" Recommendations Resulting from Zion Risk Inquiry." (80 pages) l ) l. l 1 r b

TECHNICAL ASSOCIATES TECHNICAL CONSULTANTS ON ENERGY & THE ENVIRONMENT Dale G. Bradenbaugh 1723 Hamilton Avenue-Suite K Richard B. Hubtsard San Jose. California 95125 y Gregory c. Minor Phone:(408) 266 2716 fMEDOM OF INFORMATION 2 July 1987 N.M EQEST N /A / 7-W 4, Mr. Donnie H. Grimsley, Director /n 's-~~ pf//> Division of Rules and Records 7 / Office of Adminis4 ration and Resources Management U.S. Nuclear Regulatory Commission Washington, D.C. 20555 RE: Request for Certain Documents Referencedin NUREG-1150 and Related Contractor Reports

Dear Mr. Grimsley:

Pursuant to the Freedom of Information Act and 10 CFR Part 9, Subpart A, " Freedom of Information Act Regulations", please make available at the Commission's Washington, D.C., Public Document Room single copies of the following records: ff A. Sandia National Laboratories memorandum dated 2/25/86 from D.A. Powers and D. Williams to Source Term Peer Review Team, " Radionuclides Release and Aerosol Generation During Ejection of Core Debris from a Pressurized Reactor Vessel". [ Reference #20 in NUREG/CR-4551, Vol. 5, Evaluation of Severe Accident Risks and Potential for Risk Reduction: Zion Power Plant, Brookhaven National Laboratory, Draft for Comment, February 1987) B. Sandia National Laboratories, " Instructions for Programs EVNTREISS and LAHRS', February 1986. [ Reference #15 in and Po/CR-4551, Vol. 5, Evaluation of Severe Accident Risks NUREG tential for Risk Reduction: Zion Power Plant, Brookhaven National Laboratory, Draft for Comment, February 1987) C. NRC rnemorandum dated 8/1/85 from T.P. Speis to H.R. Denton,

Subject:

' Recommendations Resulting from Zion Risk Inquiry". [ Reference #5 in NUREG/CR-4550, Vol. 7, Analvsis 1/

As used here, " records" has the definition provided in 10 CFR 9.3(b) of he Commission's regulations. Furthermore, ' records" are considered to le those in the possession of the NRC, its contractors, its subcontractors, or c.ners as provided for in 10 CFR 9.4. g +. > O h s n

I' s 02-t: of Core Damaae Freauencv From Intemal Events: Zion Unit 1. Sandia National Laboratories, October 1986] l D. F.T. Harper, et al., Analvsis of Core Damaae Freauency From Intemal Events: Methodoloav Guidelines. Sandia National l Laboratones, NUREG/CR-4550, SAND 86 2084, Vol. 1. EReference #3 in NUREG/CR 4550, Vol. 4, Analysis of Core bamaae Freauencv From Internal Events: Peach Bottom. Unit ~ 2, Sanda National I. laboratories, October 1986] E. J. Minarick, "BWR Event V," Presentation at ASEP SCG/0956 - NRC Meeting, 5/22/85. Vol. 4, Analvsis of Core [ Reference #12 in NUREG/CR 4550, Damaae Freauencv From Internal Events: Peach Bottom. Unit 2. Sandia Nabonal Laboratories, October 1986] F. A.J. CaE, et al., "LaSalle County Station ProbabEistic Safety I Analyss," General Electric Company, NEDO-31085, Class I, November 1985. Reference #16 in NUREG/CR-4550, Vol. 4, Analvsis of Core [Damaae Freauency From Intemal Events: ~ Peach Bottom. Unit 2. Sandia National I.aboratones, October 1986]- G. Letter dated 6/18/85 from G.J. Boyd ( and Reliability Optimustion Sennees, Inc.) to F.T. Harper National Laboratories). [ Reference #21 in NUREG/ R-4550,~ Vol. 4, Analvsis of Core Damaae Freauencv From Intemal Events: Peach Bottom. Un/t 2. S~andia National Laboratones, October 1986] H. Letter dated 7/2/05 from F.T. Harwr and G.J. Kolb (Sandia Nabonal Laboratories) to "PRA Experts',

Subject:

" Subtle interactions Found in Past PRAs and PRA Related Studies'.

i Reference #22 in NUREG/CR-4550, Vol. 4, Analvs/s of Core bamaae Freauency From Intemal Events: Peach Bottom. Unit 2 Sandia National Laboratories, October 1986] l. A.D. Sarain, Accident Seouence Evaluation Proaram Human . ""-/ Analvsis Procedure. Sandia National Laboratories, NUREG/CR4772, SAND 86-1996. [ References #25 and #26 in. NUREG/CR 4550, Vol. 4, Analvsis of Core Damaae Freauencv From Intemal Eventn: Peach Bottom. UnN 2. Sandia National Laboratories, Octobet 1986] - J. H.R. Denton memorandum dated 10/ 84, "NRR Office Letter No.16, Revision 2-Regulatory sis Guidelines," U.S. NRC. [Section 3, Reference #6 in NUR -4551, Vol. 4, Evakation of Severe Accident Risks and Potential lbr Risk Reduction: Grand Gulf. Unit 1. Sandia National Laboratones, April 1987] l 1

+ -. K. W.J. Lukas, et al., A Human Reliabilltv Analvsis for the ATWS Acc/ dent Seouence with MSIV Closure at the Peach Bottom Atom /c Power Station. Brookhaven National Laboratory, May 1986. [ Reference #29 in NUREG/CR-4550, Vol. 4, Analvs/s of Core Damaae Frecuencv From Internal Events: Peach Bottom. Unit 2, Sandia National Laboratories, October 1986] 1 L R.J. Dallman, et al., Severe Accident Seouence Analvsis Proaram - Anticloated Transient Without Scram Simulaelons for Browns Ferry Nuclear Plant Unit 1. Idaho Nabonal Eng neering Laboratory (EG&G Idaho NUREG/CR 4165, EGG-2379,- February 1985. [ Reference),#34 in NUREG/CR-4550, Vol. 4, Analvsis of Core Damaae Freauencv From Internal Events: Peach Bottom. Unit 2. Sandia National Laboratories, October u 1986) M. Philadel hia: Electric Company comments on Draft NUREG CR-4550, Vol. 4.

4550, ol. 4, Analvsis of[ Sag, Aspendix C in NUREG/CR-Core Damaae Freauencv From Intemal Events: Peach Bottom. Unit 2. Sandia National

~ Laboratories, October 1986) N. M.T. Drouin, et al., Analvsis of Core Damaae Freauencv frorn intamal Events: Grand Gulf. Unit 1, ~ Sandia Nationa Laboratories, NUREG/CR-4550, Vol. 6, November 1986. Section 1, Reference #4, in NUREG/CR-4551, Vol. 4, ivaluation of Severe Accident Risks and the Potential for Risk Reduction: Grand Gulf. Unit 1. Sandia Natonal Laborator:es, Aprl1987] O. J.M. Griesmeyer, User's Guide for the EVNTRE Comouter Cada, Sandia National Laboratories. [ ion of Severe Accident Section 2, Reference #5 in NUREG/CR-4551, Vol. 4, Evaluat Risks and the Potential lbr Risk Re duction: Grand Gulf. Unit 1, Sandia Natonal Laboratones, Apri 1987) P. M.L Corradini " Additional Comments on the Peach Bottom and Grand Gulf Containment-Event Trees," letter repoit to Mark Cunningham, U.S. NRC, 10/13 Reference -#7 in NOREG/CR 4551, Vo/86.l. 4, Ev[aluafkrz of Sectm 2, Senere Accident Risks and the Potential ibr Risk Reducti.go; j Grand Gulf. Un/t 1. Sandia National Laboratones, April 1987) O. Brookhaven National Laboratories, Rev/sw of the Source-Term Code Package. Section 2, Reference #10 in NUREG/CR-4551, Vol. 4, Eva[luation of Severe Accident Risks and the i Potential for Risk Reduc 60n: Grand Gulf. Unit 1, Sandia Nabonal Laboratories, Apri 1987) [

R. Sandia National Laboratories, Probab/ list /c Risk Uncertiantv Evaluation Prooram (PREUP). Section 2, Reference #12 in NUREG/CR-4551, Vol. 4, Evalua(tion of Severe Accident Risks and the Potential for Risk Reduction: Grand Gulf. Unit 1. Sandia National Laboratories, April 1987] S. F.T. Harper, "ASEP Data Reevaluation,' memorandum to Senior Consultant Group, et al., Sandia National Laboratories, 3/15/85. [Section 2, Reference #18 in NUREG/CR-4551, Vol. 4, Evaluation of Severe Accident Risks and the Potential for Risk Reduction: Grand Gulf. Unit 1. Sandia National Laboratones, April 1987) T. A.S. Benjamin and F.T. Harper, Value Imomet Investicetion of Filtered-Vented Containment Svstems and Other Safetv Ootions for a BWR Mark I Containment. Sandia Natsonal Laboratories, NUREG/CR-4065 _ (Draft), Septebmer 1984. 'Section 3, Reference #5 in NUREG/CR-4551, Vol. 4, & valuation of Severe Accident Risks and the Potentialibr Risk Reduction: Grand Gulf. Unit 1. Sandia Natnonal Laboratories, April 1987) U. J.L Dooley,, et al., Mitloation Svstems for Mark Il Reactors. R&D Associates, RDA-TR-127303-001, May 1984. [luation of Appendix F, Reference #12, in NUREG/CR-4551, vol 4, Eva Severe Accident Risks and the Potentist for Risk Reduction: . Grand Gulf. Unit 1. Sandia National Laboratories, April 1987) V. Letter dated 6/4/86 from K.D. Bergeron (Sandia National to T.M. Lee (NRC). [Section 4, Reference #4 in laboratories)4551, Vol. 2, Evaluation of Savare Accident Risks NUREG/CR-at the Potential for Risk Reduction: Seouovah Power Station. Uhlt 1. Sandia National Laboratories, February 1987) W. latter dated 3/19/86 from G. Greene (Brookhaven National Laboratories) to A.S. Benjamin (Sandia Nabonal Laboratories). Appendix B, Reference #20 in NUREG/CR-4551, Vol. 2, kluation of Severe Accident Risks at the Potential for Risk Reduction: Seouovah Power Station. Unit 1 Sandia National Laboratones, February 1987) W.J. Camp, Progression and Vessel Attack Under Severeet al., A Mechanist X. Core Melt L Accident Condit7ons, Sandia National Laboratories, SAND 85-0237, draft. (Appendix J.2, Reference J.2.12 in NUREG-1150, Vol. 3, Reactor Risk Reference Document. February 1987) In eadi case above, the requested documents are references in other published NRC reports (references are provided). /f there is a more recent version of the document than Ested in the reports, please provide the most recent version rather than the one listed in he referenced reports.

,e i l If you or any members of your staff have any questions conceming this request, please contact the undersigned directly by telephone at (408) 266-2716. Partial l responses would be appreciated as documents become available. Your prompt attention to this request will be apprecated. I Sincerely, Steven C. Sholly Associate Consultant i i l 1

y q (- DRAFT PRELIMINARY REPORT FOR COMMENT l l l* FIN NO. A-1342 3 l4 l ~. <. MELPROG-PWR/ MODO: A Mechanistic Code for Analysis of Reactor Core Melt Progression and Vessel Attack under Severe Accident Conditions Prepared by Reactor Safety Theoretical Physics Division Sandia National Laboratories and Safety Code Development Group, Energy Division Los Alamos National Laboratory for U.S. Nuclear Regulatory Commission NOTICE THIS DRAFT PRELIMINARY REPORT IS ISSUED TO USNRC SEVERE FUEL DAMAGE PROGRAM PARTICIPANTS ONLY This report was prepared in contemplation of Commission action. It may not have received patent review and may contain information received in confidence. Therefore, the contents of this report should neither be disclosed to others nor reproduced, wholly or partially, unless written permission to do so has been obtained from the appropriate USNRC office. The recipient is requested to take the necessary action to ensure the protection of this report. p3p \\ .) j

DRAFT PRELIMINARY REPORT FOR COMMENT l FIN NO A-1342 MELPROG-PWR/ MODO: A Mechanistic Code for Analysis of Reactor Core Welt-Progression and Vessel Attack under Severe Accident Conditions Responsible Staff: W. J. Camp + Project Leader M. F. Young + Structural Mechanics J. L. Tomkins+* Structural and Radiation Heat Transfer J. E. Kelly + Debris Bed Analysis P. J. Maudlin + Fluid Dynamics R. J. Henninger** TRAC-MELPROG Interface and Fuel Rod Analysis ) l 1 l l Sandia National Laboratories, Albuquerque, NM 87185 +- NM Science Applications International Corp., Albuquerque, j Los Alamos National Laboratory, Los Alamos, NM 87545 I I

Abstract l The MELPROG computer code is being developed to _ provide detailed, best-estimate, coupled analyses of all ~ the major phenomena involved in the course of a severe accident progression. The initial version of this code, MELPROG-PWR/ MODO, is described in this report. I s The purpose of MELPROG is to provide at any time during an accident sequence a description of (1) the state of the reactor core and surrounding in-vessel environment, and (2) the disposition of core materials (in particular, fission products) contained within the reactor coolant system boundary. HELPROG is coupled to the TRAC-PF1 RCS thermal-hydraulics code to provide an integrated analysis of the behavior of core, vessel, and HELPROG reactor coolant system during severe accidents. treats core degradation and loss of geometry, debris formation, core melting, attack on supporting structures,

slumping, melt / water interactions and vessel' failure.

The key element in MELPROG is the use of detailed modeling for the entire demage progression and failure -sequence, Emphasis is also placed on the rates of

hydrogen, steam and fission product.formatien, and transport to containment during damage progression.

1 l 'l l ) l l iii,1v

o TABLE OF CONTENTS Page .1'. INTRODUCTION AND CODE DESIGN 1 1.1 Introduction 1 1.2 MELPROG Design and Structure 3 1.3 References 9 2. FLUID DYNAMICS 11 2.1 Introduction 11 2.2 Conservation Equations 11 2.3 -Constitutive Equations 14 2.4 Mass-Transfer Functions 17 2.5 Solution Technique 19 2.6 TRAC /MELPROG Link 20 2.7 References 20 3 FUEL AND CONTROL RODS 22 31 Introduction 22 32 Model Descriptions 22 3 2.1 Cladding Rupture and Ballooning Models 24 3 2.2 zircaloy Oxidation. Hydrogen Generation, 24 and Embrittlement 323 Control Rod Meltdown Models 29 3 2.4 Fuel' Cladding Liquefaction Phenemens 30 3 2.5 Fuel Rod Failure and Debris Formation 32 33 Fuel Pod Temperature Analysis Model 34 3.4 Materials Properties, Gap Conductance. 39 and Cladding Stress 35 References 41 4. STRUCTURES 43 4.1 Introduction 43 4.2 Structure Heat Transfer 43 4.3 Structure Types 43 4.3 1 Wall Type Structures 46 4.3 2 Plate Type Structures 46 4.3 3 Column Type Structures 48 434 General Type Structures 48 4.4 Mesh Arrangement 48 4.5 Material Properties 49 4.6 Conduction Heat Transfer Solution Technique 49 I v

4.6.1 Formulation'for Material Nodes 49 4.6.2 Surface Nodes 50 4.6 3 Matrix Solution. 51 4.7 Melting;and Ablation -52 L 4.8 Structural Mechanics'Modeling 54 4.9 Modeling Design Considerations 54 4.10' Mechanics Models 55 1. 4.10.1 Lower Core Grid Plate 57 4.10.2 Elliptic Plate 60 62 4.10.3 Cylinder 4.10.4 Columns 63 4.10.5 Example 63 64 4.11 References 65 5. BRIS MODULE 65 5.1 Introduction 66 52 overview of the Module 53 Debris Bed Formation and Release Model 68 '5.4 Post-Dryout Debris Bed Model 71 5.4.1 Debris Heat-up Model 71 5.4.2 Melt Relocation Model 74 5.4.3 Debris Phase Diagram 84 5.4.4 Zr oxidation Model 85 87 5.5 References 6. ' RADIATION HEAT TRANSFER iS 8 88 6.1 Introduction-88 6.2 Model Theory and Design 6.3 Rad.tation Heat Transfer Model Implementation 92 6.3 1 Geometry and Fluid Fields Modeled 92 6.3 2 Solution Procedure 93 6.3 3 View Factor Calculation 95 6.3.4 Material Properties and Mean Transmission Coefficients 97 6.3 5 Summary of Basic Modeling Assumptions 98 6.3 6 RADIATION Mocule Interaction with other Modules 98 98 6.4 References vi

7 SAMPLE CALCULATION 99 7.1 Introduction 99 7.2 Problem Setup 99 73 Results of S1D Test Problem 102 731 FLUIDS Module 102 7.3 2 PINS Modul0 108 7.3 3 DEBRIS Module 110 7.3.4 RADIATION Module 111 735 STRUCTURES Module 111 7 3.6 Summary 113 } 7.4 References 114 l l APPENDIX A " CODE INPUT DESCRIPTION 115 A.1 Introduction 115 A.2 General Description 115 A3 Input cards 115 l A.3 1 General Input 115 A.3 2 Fluid Dynamics Input 117 A.3 2.1 Fluids Scalar Input 117 A.3.2.2 F1 tr i d Dynamics Array Input 118 A.3 2.3 Fluids Initial Conditions 118 A.3 2 3 1 Field One - Steam 118 and Hydrogen A.3.2.3 2 Field Two - Water 119 A.3 2 3 3 Field Three - Corium 119 A.3 2.4 Fluid System Boundary Conditions 119 A.3 3 Rod Model Input 120 A.3 3.1 Rod Model Sesinr Input 120 A.3 3 2 Rod Model Array Input 121 A.3.4 Structures Input 122 A.3 5 Radiative Heat Transfer Input 125 A.3.6 Debris Input 126 A.4 General Rules for Problem Restarts 127 A.5 Input cards 127 't A.S.1 General Input 127 A.5.2 Fluid Dynamics Input 129 vii

A.5.2.1 Fluids Scalar Input 129 A.5.2.2 Flv!r Dynamics Array Input 130 A.5.2 3 Fluida Initial Conditions 130 Steam 130 l A.S.2 3 1 Field One and Hydrogen A.5.2 3 2 Field Two - Water 131 A.5.2 3 3 Field Three - Corium 131 A.S.2.4 Fluid System Bounddry Conditions 131 A.5.3 Rod Model Input 132 A.5.3.1 Rod Model Sc31ar Input 132 A.5.4 Structures Input 133 A.S.5 Radiative Heat Transfer Input 133 A.5.6 Debris Input 133 A.6 References 134 APPENDIX B - CODE OUTPUT DESCRIPTION 135 B.1 Introduction 135 B.2 FLUIDS Module Printout Block 136 B.3 PINS Module Printout Block 136 B.4 STRUCTURES Module Printout Block 136 B.5 DEBRIS Module Printout Block 137 B.6 RADIATION Module Printout Block 137 B.7 Sample Output 137 APPENDIX C d CODE STRUCTURE 141 C.1 Introduction 141 C.2 FLUIDS Module 141 C.3 PINS Module 145 C.4 STRUCTURES Module 148 C.4.1 Heat Transfer Routines 148 C.4.2 Mechanics Routines 150 C.4 3 Flow Areas and Hydraulie Diameters 153 C.5 DEBRIS Module 153 C.6 RADIATION Module 155 C.7 TRAC /MELPROG Link Computational Flow 156 C.8 MELPROG Data Input and Initialization l Module (INPUT) 156 C.9 Data Management Structure Used in MELPROG 161 viii t

LIST OF FIGURES Illung gagg 1.1. Modular Structure of MELPROG 4 3.1. Comparison of Burst Temperature / Pressure 25 Correlations at Low Temperature Ramp Rates 32 Comparison of Burst Temperature / Strain 25 Correlations at Low Temperature Ramp Rates 33 Cross Sect}on of Fuel Rod with Oxidized Cladding 26 3.4 Zr-UO Binary-Phase Diagram 31 2 3.5 Finite-difference Noding Scheme Used in 35 Fuel Rod Model 4.1 Typical Plate and Wall Type Structures 45 4.2 Structure Arrangement in a Cell 47 4.3 Procedure for Melting and Ablation 53 4.4 Westinghouse PWR Schematic 56 4.5 Cell of Westinghouse Lower Core Grid Plate, 58 Top View 5.1 Geometry of Infinite Reservoir Case 76 5.2 Slug Penetration into Bed 78 5.3 Simple Hydrodynamic Model Compared to PLUGM 81 5.4 Simple Thermal Model Compared to PLUGM 82 6.1 Typical Geometries Treated in RADIATION Module 94 7.1 S1 D Vessel Nodal.t zation 1 01 7.2 Vapor Volume Fraction Distribution 104 l l 73 Corium volume Fraction Distribution 106 l 7.4 Corium Temperature Distribution 107 1. 7.5 Cladding Wall Temperature in Ring 1 109 7.6 Core Barrel Temperature Distribution 112 ix

i; 142 C.1 Overall MELPROG Structure C.2 Functional Diagram of FLUIDS Module Routines 143 C.3 Functional Diagram of Intact Fuel Rod Model 147- ~ C.4 Functional Diagram'of Structure Heat 149 Transfer. Routines C.5 Functional' Diagram of Structure Mechanics 149 Routines. C.6 Functional Diagram of DEBRIS Modul,e Routines 154 C7 Functional Diagram of RADIATION 154 Module Subroutines C.8 Outer Iteration Flow Chart for TRAC /MELPROG 157' C.9 Functional Diagram of INPUT Module Subroutines-158 LIST OF TABLES Pagg I*bla 2.1 Materials Conserved in the Fluid 13 Dynamics Solution 31 Fuel and Control Rod Morphological Forms 23 28 32 Cladding Layer Growth Rate Constants 44 4.1 Examples of Structure Types 86 5.1 Summary of DEBRIS Phase Change Rules 6.1 View Factor Matrix for a Typical Type 1 Cell 96 7.1 Summary of Events in S1D Calculation 103 151 C.1 Ring Load Table Data Structure 160 C.2 Geometry Table Data Structure 160 C.3 Loading Table Data structure C.4 Transfer Processes Needed to Obtain FLUIDS 163 Data l 1 l i

1. INTRODUCTION AND CODE DESIGN 1.1. Introduction Since the Three Mile Island accident, increased emphasis has been placed within the reactor safety community on analysis of severe reactor accidents. In particular, it has been recognized that such accidents represent the dominant public risk from nuclear reactors.- The U.S.N.R.C. has made the development of mechanistic-models for severe accident progression a major priority.'[1.13 .The purpose of these models is to provide

detailed, best" estimate, coupled analyses of all the major phenomena involved in the course of the accident.

For the in' vessel and reactor coolantasystem (RCS) aspects of severe accidents, the U.S.N.R.C. is developing the severe fuel damage analysis package. This is a two' fold model development effort. On the one hand, the Severe Core Damage Analysis Package (FCDAP[1.23) model is being developed at the Idaho National Engineering Laboratory by EG&G Idaho, Inc. This code provides an extremely detailed treatment of rod behavior and degradation of light water reactor (LWR) cores during the initiation and heat-up phase of LWR accidents. SCDAP is to utilize the RELAP-S[1 3] RCS thermal-hydraulics code to provide an.inte6 rated treatment of the reactor and coolant system during the initiation and heatdup phases of the accident. On the other

hand, the MELPROG (MELT Progression) mode 1[1.4.1.53 is being developed at Sandia and Los Alamos National Laboratories.

HELPROG is coupled to the TRAC-PF1 RCS thermal-hydraulics code [1.63 to provide an integrated analysis of the behavior of core, vessel, and RCS during severe accidents. MELPROG provides a detailed core model, so that it can provide accident analysis capability _ from initiation through vessel failure.

However, MELPROG is not intended to model the heatbup phase in as detailed a manner as SCDAP.

Rather, the major emphasis in the overall MELPROG effort is on the treatment of core degradation and loss of geometry, debris formation, core melting, attack on supporting structures, slumping, melt / water interactions and vessel failure. The key element in MELPROG is detailed modeling of the entire damage progression and failure sequence. Emphasis is also placed on the rates of hydrogen, steam and fission product formation, and transport to containment during damage progression. An extensive discussion of the goals and design concepts for MELPROG has already been presented in MODELING OF LIGHT-WATER REACTOR SYSTEMS DURING SEVERE ACCIDENTS: MELPROG PERSPECTIVE, (Camp and Rivard, eds.)[1.4]. Herein we first review the MELPROG l: 'g program design, and then d-in the chapters to follow '- discuss I the detailed design and structure of the first version of MELPROG, MELPROG-PWR/ MODO. 1 ) 1

The purpose of MELPROG is to provide at any time during the accident sequence (i.e., from the initial core-damaging stage to the stage at which the vessel fails and core materials are discharged to the containment) a description of (1) the state of the reactor core and surrounding in-vessel environment, and (2) the disposition of core materials (in particular, fission products) contained within the reactor coolant system boundary. ~ The principal processes which are to be modeled by MELPROG include: 1. Vessel thermal-hydraulics and core heat-up; 2. Liquefaction and melting of reactor core materials; l 3 Solidification and freezing of fuel, cladding, control, and structural materials; 4 Fragmentation and relocation of liquefied and solidified materials; 5. Stressing and ablation of maj or structural features including core supports and pressure vessel; 6. Energy-generation processes, including decay heating and chemical oxidation; 7. Energy-transfer and mass-transport processes, including conductive, convective, and radiative heat

transfer, steam explosions, and fluid flow; and 8.

Associated thermochemical processes which affect the re-lease, transport, and deposition of fission products. To determine what models are needed to fashion MELPROG and also what attributes are required of these

models, it is necessary to consider:

1. The differences in various features of the designs of pressurized water reactor (PWR) and boiling water reactor (BWR) systems that make it necessary to provide design-specific models of the in-vessel processes; 2. The PWR and BWR severe-accident sequences that initiate core-damaging processes and that prescribe the nature and course of development of these processes; and 3 The treatment to be accorded these processes in terms of both the level of coverage (i.e., the number of potential accident sequences and therefore of potential processes covered) and the level of treatment (i.e., the amount of detail elaborated in modeling these processes). j 1 i 4 l 2 i l j

MELPROG consists of several explicitly linked modules which, in turn, are comprised of models that represent and treat the physical processes that occur during a severe core-damaging accident. These modules only interact through well defined interfaces. Any module can be run-stand-alone if suitable interfaces are provided. Similarly, any module can be removed ~' and replaced by dummy interfaces and the code will still run. Although the general functions performed by the modules in MELPROG are nearly equivalent for PWR and BWR

systems, the specific models that are needed to form these modules are not.

Due to the substantial differences in.the design of the two raystems, distinctly different models are required to treat the reactor-specific systems and components because these can have significant impact on the timing and importance of various accident sequences. As a

result, the initial MELPROG code development effort has been focused on the development of models for PWR designs.

1.2. HELPROG Design and Structure The MELPROG-PWR/ MODO structure consists of five main modules coupled explicitly through a driver program that handles overall program flow (see Figure 1.1). The main modules are

FLUIDS, PINS, RADIATION, STRUCTURES, and DEBRIS.

HELP ROG-PWR / MOD 1. will include at least three additional modules: IFCI, FISSION PRODUCTS and EJECT. In addition to these modules, all versions of MELPROG contain an implicit numerical interface to TRAC-PF1 (TRAC-LINK) which enables the entire reactor coolant system to be modeled. The roles'or each module and the TRAC-LINK routine are: 1. FLUIDS: Treats multi-component, multi-field compressible flow and heat transfer within the reactor vessel. In i PWR/ MODO, FLUIDS uses a one-dimensional model and treats three fields: liquid coolant, vapors, and corium. (Note that is this report the terms "corium", " debris", and " rubble" used interchangeably.) In PWR/ MOD 1 and in all BWR versions a two-dimensional four-field treatment will be used and the transport of fission products will be calculated. The use of four fields allows the corium field to be divided into a liquid corium field and a solid corium field. 2. PINS: Treats energy. generation, heat

transfer, thermochemistry, liquefaction and melting of fuel and cladding, cladding mechanical response, and rod failure and breakup.

PINS provides a simplified treatment of the rod heat-up and failure processes considered in - g detail in the SCDAP model. In PWR/ MOD 1, better failure l models and detailed candling models will be included. l In all BWR versions, r. significantly differeat, approach l 3 l

TRAC UNK 1 t DRIVER l 1r 1r 1r 1r 1r RADIATION PINS FLUIDS STRUCTURES DEBRIS l 1r 1r VICTORIA IFCI 1r EJECT Figure 1.1 Modular Structure of MELPROG p 4 w-----_-___-____________

t-will be needed because of the existence of can walls and control rod blades. 3 ' RADIATION: Treats radiation heat transfer within the core and vessel. The model uses two-dimensional cells (R-Z) which are' diffuse emitters and' absorbers. Within each cell, rods, debris, and structures are treated by-three-dimensional view factors which change to acccamodate geometry changes. 4 STRUCTURES: Treats mechanical -and thermal behavior of all structures within the vessel, as well as the vessel itself. A number of generic structure types are used, such as perforated plates, columns, cylindrical shells, etc., to represent specific structures. The structure package currently accommodates PWR's (e.g., Babcock and

Wilcox, Westinghouse, and Cembustion Engineering).

Structural heat transfer is via one-dimensional con-duction, convection, an'd phase-change (ablation models). Radiation to and from structures is handled by RADIATION in a consistent fashion. Ablation mass transport from structures is treated in a consistent - f ashion by the STRUCTURES and: FLUIDS modules. In PWR/ MOD 1 and all subsequent

versions, two-dimensional heat transfer models and simultaneous crusting and ablation will be treated.-

Also, in BWR versions, the structure package will have some significant differences from that for any of the PWR designs. 5. DEBRIS: Treats the formation of rubble or debris beds throughout the vessel, using simplified models. Once a rubble bed is

formed, DEBRIS handles heat transfer within the bed.

However, FLUIDS continues to treat coolant

boiling, coolant dry-out, reflooding, and quenching.

For a dried-out

bed, DEBRIS calculates material liquefaction and relocation, refreezing and crust / blockage formation.

DEBRIS determines heat fluxes te underlying structure and heat transfer from debris to the FLUIDS coolant field as well. 6. VICTORIA: Treats release and transport of fission products within the vessel. The ultimate driving forces for release of fission products are the differences I between eq.111brium and existing partial vapor pressures of all relevant molecular species in the overlying vapor. A number of limitations on release are modeled including diffusion within the fuel / cladding

matrix, va;1rization kinetics and transport across boundary

,,I layers. Condensation onto fuel,

cladding, structures and pre-existing aerosols is treated.

The aerosol formation, agglomeration, and settling processes are also modeled[1.7]. This module is not included in 5

PWR/ MODO, but will be included in PWR/ MOD 1 and all subsequent versions. 7. IFCI: This integrated fuel coolant interaction model treats explosive and non-explosive melt / water inter-actions, including detailed modeling of the detonation process itself. IFCI useo correlational models for fuel / coolant inter-mixing, film

boiling, and fuel fragmentation.

These models are used in a finely noded FLUIDS treatment to calculate detonation wave propaga-

tion, vapor generation and post-explosion expansion.

This module is not included in PWR/ MODO, but will be in-cluded in PWR/ MOD 1 and all subsequent versions. 8. EJECT: Treats release of materials from the vessel associated with various vessel failure modes. FJECT handles the entire spectrum of release modes from gravity slumps to high pressure blowdowns. For the

latter, EJECT calculates melt fragmentation by effervescence, atomization and hydrodynamic fragmen-tation.

It then calculates the aerosol fraction and size distribution delivered to the containment atmosphere. This module is not included in PWR/ MODO, but will be included in PWR/ MOD 1 and all subsequent versions. 9. TRAC-LINK: Provides an implicit numerical interface to the TRAC-PF1 reactor coolant system thermal-hydraulics model by allowing MELPROG to be treated as a CORE or VESSEL module in TRAC-PF1. This link effectively makes MELPROG part of TRAC.

However, until T R A C-P F1 is generalized suitably, only liquid coolant and the vapor field can be transported through this link.

Corium materials passed through the interface to TRAC-PF1 are lost from the problem. Although aerosols and fission product vapors can, in princip1b, be transported within the vapor field, the important transport and deposition procesbes are ignored. Future versions of TRAC are expected to address this problem. As discussed above, MELPROG-PWR/ MODO, the version of MELPROG j discussed herein, incorporates only the link to TRAC plus the j

namely, the 1-D FLUIDS,

PINS, first five of these modules RADIATION, STRUCTURES, and DEBRIS modules.

Although the FLUIDS module used in PWR/ MODO is one-dimensional, the remaining modules use two-dimensional treatments in anticipation of incorporation of the two-dimensional FLUIDS treatments in subsequent versions. i It should be noted that the PWR/ MODO FLUIDS, PINS, and RADIAT79N modules have been taken over from the MIMAS code effort [1.82 at Los Alamos National Laboratory. However, except for the PINS module which has remained essentially unchanged, these modules have been extensively improved and extended for use in MELPROG. 6

STRUCTURES and DEBRIS [1.93 represent completely new development efforts as do the two-dimensional FLUIDS [1.10] package, VICTORIA, IFCI, and EJECT [1.11] which are being developed for subsequsnt MELPROG versions. In a typical accident simulation, the calculation begins with an initially intact core and vessel. A break in the coolant system or loss of heat removal capability is specified. In either case, coolant loss from the system is calculated by the FLUIDS module. In a complete system simulation, TRAC-LINK is employed to couple the vessel with a TRAC-PF1 model of the RCS. In some cases, a set of boundary conditions is used to simulate the TRAC calculation. This latter approach is most appropriate if a detailed TRAC of dELAP calculation of the coolant blowdown or boil-off sequence is available. As the coolant level drops below the top of the core, rapid core heat-up begins due initially to decay heat generated in the fuel rods. Above 1273 K, rapid oxidation of the Zircaloy cladding begins near the top of the core. This contributes even more heat to the core. In

addition, it generates large quantities of hydrogen gas.

In many cases, ballooning of the cladding will occur. In any case, a " burn front" will dev.elop in the upper regions of the core and move upstream to the' lower portions of the core as the accident progresses. The upper regions of the core are often sufficiently hot and damaged to fail even before water has completely boiled out of the core region. The heatup of the core, cladding oxidation, and rod failure are all calculated by the PINS module in concert with the FLUIDS module. While the core is heating up, heat is transported to surrounding structures by FLU 1DS (convection heat transfer) and especially by the RADIATION module. Both the core barrel and the upper core grid plate heat up rapidly. Indeed, the STRUCTURES module calculates the thermal and mechanical response to these loads. For some accidents, the upper grid plate is calculated by STRUCTURES to fail fairly early in the accident. During the heat-up process, the first failures in the core are typically the control rods. As these fail, they form debris which is passed by PINS to the FLUIDS module for transport and possible formation of debris beds. As rods begin to fail, PINS also adds this debris to the FLUIDS corium field. A simple debris bed formation model is employed to form stationary debris beds in the lower core region (as well as in the lower head). When debris beds are formed. FLUIDS passes control of the corium involved to the DEBRIS module, i i Ordinarily, the debric will hang up on lower pin stubs l*- either at grid spacers, or in narrowed or blocked sections where candled cladding has refrozen. If water remains in the region 7 l

where' debris beds form, FLUIDS maintains control of the heat transfer from the debris, and performs an average calculation for each fluid node.

However, if the region is initially dry or becomes dry later, DEBRIS would perform a finely-meshed analysis of debris heat-up,

melting, relocation,

crusting, and crust melting.

DEBRIS also calculates the thermal loads to underlyinE structures. If a debris bed begins to move through failure of its supporting structure, then DEBRIS releases the debris bed and returns control to the FLbIDS module. Eventually, the debris beds so formed become largely molten and the thermal loads on the grid plate become sufficient to cause local failures of the grid plate. This failure process as well as thermal processes in the structures are modeled on a cell-by-cell basis by STRUCTURES. After grid plate failure, the partially molten corium slumps into the lower head (under control of FLUIDS) and is quenched. During the quenching process copious quantities of steam and hydrogen are generated. These are transported by FLUIDS. Note that one should expect a steam explosion at the point of fuel slumping into the water pool. However, MELPROG-PWR/ MODO does not contain IFCI, the module under development for subsequent

versions, which will treat such explosions.

Finally, the quenched debris heats-up, boils the remaining coolant dry, and remelts under renewed control of DEBRIS. The molten corium melts through the vessel wall. This process is again calculated by STRUCTURES using the thermal fluxes generated by DEBRIS. At this stage, melt begins to exit the vessel. This process is not modeled in detail by MELPROG-PWR/ MODO. That is, the hydrodynamics of the exiting melt is not treated in the current version. The EJECT module, which currently exists in stand-alone

form, will provide that detailed treatment in subsequent versions.

A couple of points need to be noted concerning this " typical" scenario.

First, the code is not limited to the scenario described above.

Rather, that scenario is fairly typical of the scenarios calculated to date by PWR/ MODO, namely the S1 D (medium-break, loss of coolant accident) and the station blackout accident (TMLB').

As more scenarios are treated, and especially when the two-dimensional version of FLUIDS is employed for boil-off accidents, we may expect significant deviations from i the above scenario. Second, in subsequent versions of MELPROG, the VICTORIA fission product module will sciculate release of fission products from pins, debris, and molsen corium in concert with the DEBRIS, PINS, and IFCI modules, and transport of fission products in concert with the FLUIDS module. In the remaining sections of this document, we turn first in Sections 2 through 6 to detailed descriptions of the five PWR/ MODO modules. The TRAC-LINK interface is discussed in 8 h l L____-_-_____ l

Section 2,. Then in Section 7, we describe the tecting of MELP ROG-PWR /MO DO performed in the context of a small-break loss of coolant accident (LOCA) for a ' PWR. To be specific, a S1 D small-break LOCA accident in the Zion pressurized water reactor plant has been analyzed in some detail. That analysis and its implications for MELP ROG ' assessment are presented in Section 7.

Finally, code input and output descriptions are given. in Appendices A-and B,

respectively, while Appendix C presents a dota11ed layout of code structure. 1 1.3 Refegences j i 1.1. U.S.N.R.C., The Nuclear Power Plant Severe Accident Re- ~ UREG 6960 (January 1553). ~ N I search Plan, 1.2. G. A. Berna et al., SCpAPfMQpiflq1__&_Cgmgutgt_Cgde_[gt the Analysis _of__ LWR, Vessel _ Behavior _Duritg Severe _ Accident j Transients, IS-SAAM-64-002, EG&G Idaho, Inc., Idaho Falls, ID, (July 1984). _ ode Manual, Volumes 1 -) 1.3 V. H. Ransom, et. al., RELAP5fMOD1 C and 2 NUREG/CR-1326, EGG-2070 EG&G Idaho, Inc., Idaho Falls, ID, (1980).

1. 4. -

W. J. Camp and J. B. Rivard, Modeling of _ LWR _Sys t em s Dur-I ing Severe Accidents: The MELPROG Perspective, SANDF3-1205, Sandia National Laboratories. Albuquerque, NM, (to be published). z 1.5. M. F. Young et al., "HELPROG Code Development and Methods," Proceedings _ International _ Meeting _on_ LWR _ Severe 2 Accident Evaluation, Vol. 1, p.2.5-1, Cambridge, MA ---~~-~~- TTVETT. i ~1,6. T R A C-P F1 : An Advanced Best-Estimate Computer Program fog Pggssugized Water Reactog Analysis, LA-9944-MS, NUREG/ l CR-3537, Los Alamos National Laboratory, Los Alamos, NM, 1 } (Feb. 1984). ) 1.7. W. J. Camp, " Release of Fission Products from Fuel during i the In-Vessel Phases of Severe Nuclear Reactor Accidents," European Applied Research Report - Nuclear Science Technologr. Vol. 5.

p. 1605 (1984).

't 1 i h j I l 9 l I I L j

~ 1.8. P. J.. Maudlin et al, Light-Water Reactor Degraded-Core i Cooling Program Technical Note (A Damage Assessment of j TMI-2, LA-UR-83-737, Los Alamos National Laboratory, Los ' i Alamos, NM, (1983). ] 1.9. J. E. Kelly, " SPARTA: A Core Debris Meltdown Model," Internal Sandia Internal Communication. i 1.10. J. F.

Dearing,

"A Four-Field Model of PWR Degraded Cores," LA-UR-84-3587, submitted to the 3rd International . Conference on Reactor Thermal-Hydraulics (to be published). 1.11. M. Pilch, Sandia Internal ~ Communication. 4 0 10 s_

2. FLUID DYKAMICS 2.1. Introduction The fluid dynamic modeling and methods described in this section are based tpon those of the TR AC-PF1[2.13 computer code development effort. The form of the conservation equations and the solution techniques used in MELPROG are straightforward extensions of those in TRAC that have undergone many man years of evolution. The module which controls the fluid dynamics calculation is called FLUIDS. Although this section will briefly describe the fluid i'y na mi c s used in

MELPROG, the interested reader is referred to Reference 2.1 for more detail.

It should also be noted that the discussion herein is limited to the one-dimensional model actually employed in MELPR00-PWR/ MODO. 2.2. Conservation Eguations The fluid-dynamic conservation equations solved in the FLUIDS module are one-dimensional versions of the following general form B(ag p gj) 6IT (2.ta)

• V'( l'ij"i}

" k-1 ijk Bt g[ l' i ~~IVP t ik ~ "1 P ig (2.1b) f og p kei B(ape)g 4 V.(apev)g - + g ~P[ Bag + V.(av)g) + qik ijk jk h (2.1c) + k=1j=1 Eqs. 2.1a through 2.1c conserve

mass, momentum, and internal energy, respectively.

The symbols in these equations are defined as follows: fluid field index, 1 = caterial index, J = field index for fluid fields and stationary fields, k = volume fraction of field 1, o = g p g macroscopic density for field 1. = 11

g) macroscopic-density for material j in field 1, p = vg velocity vector for. field i, =

time, t

V. divergence-operator. = -V + gradient operator, gravitational acceleration, g = f ijk= mass transfer function for material j from field k to field 1 total. static' pressure., P = voluoetric friction force vector between fields i and k, f = ik eg-internal energy for field 1, = egj internal energy for material j in field 1, = gik volumetric heat transfer between i and k, or power = density (when i=k), material j-enthalpy, internal energy, or reaction heat for h = jk field k. Note that index k ranges from 1 to 6: 1= vapor field, 9= water field, 3=corium field, 4= stationary fuel rod field, 5= stationary structure field, and 6-stationary debris bed field; and index i ranges from 1 to 3: 1= vapor field, 2= water field, and 3=corium field. Table 2.1 itemizes the materials that are conserved in the fluid-dynamics solution; that is, a continuity equation of the form of Eq. 2.1a is solved for each of the nine materials (the index j ranges.from 1 to 9) listed in. Table 2.1. The third column of this table indicates the field of motion of each material and the fourth column indicates the energy field associated with each material. Currently, a momentum equation of the form of Eq. 2.1b is solved for each of the three fields: vapor, water, and corium. Also an energy equation of the form of Eq. 2.1c is solved for each of the same three fields. In other words, materials 1 and 2 ( j = 1-to 2) are combined to form the vapor

field, for which motion and energy conservation is performed.

For the same purpose, material 3 (j=3) forms the 'I water

field, and materials 4

through 9 (j=4 through 9) are combined to form the corium field. The dependent variables being solved for within the MELPROG fluid-dynamics ares three volume fractions (vapor,

water, corium),

three velocities (vapor,

water, corium),

three temperatures (vapor, water, corium), a total pressure (mechanical nonequilibrium between fields ignored), a hydrogen partial

pressure, and macroscpoic densities for six fuel rod / structure materials (UO Zr, Zr0 stainless steel, control rod absorber 2,

2, material, and eutectic). This is a total of seventeen unknowns that are solved for using fifteen conservatica equations, - ) appropriate equations of. state, phase change relationships, j -constitutive equations for friction and heat-transfer L coefficients, and the following constraints on volume fractions, 12 l l L__

Table 2.1 Materials Conserved.in the. Fluid Dynamics-Solution. Material. Velocity Energy Number Name Field Field 1 Steam Vapor Vapor 2 Hydrogen Vapor Vapor 3 Water Water . Water j 4 UO Corium Corium 2 5 Zr Corium Corium 6 Zr0 Corium Corium 2 7 Steel Corium Corium 8 Control Corium Corium 9 Eutectic Corium Corium t' 13 4 L

partial' pressures, field'. macroscopic . densities, and-field internal energies:- 1 - a3 + a2 *T"3 (2.2a). -P'- Pateam + P (2.2b) H. 2 (2.2c) p g pij (2.2d> p3 "3 a 3J - (2.2e) P p1j'1j (2.2f) e3 = j 1 p1 e2 '23' (2.2g) e3.-}(.3J'3J (2.2h) 3-93 Eq. 2.2a conserves total' fluid flow volume, Eq. 2.2b sums the partial. pressures to.' form the total pressure.. and Eqs. 2.2c through - 2.2h. define ' the field macroscopic densities ~ and internal energies 'in. terms of the material macroscopic densities and material internal energies.respectively. 23 Sans111u11y..g Egya11gns_ For the sojution of Eqs. 2.1, relationships for the friction . force vector (f ik) and the vo lu me t r i,c heat transfer term (qik) are needed. In particular, resistance and heat transfer Ocefficients are needed for the following six combinations: 1. vapor and water, 2. vapor and stationary fields (i.e., fuel rods, structure, and debris beds), 3 water and stationary fields, ~ l-4 water and corium, i 14

5 vapor and corium, 6. corium and stationary fields. The actual values of these coefficients are, of course, geometry and flow regime dependent. It should be noted at this point that the MELPROG-PWR/ MODO constitutive equations do not address the steam explosion regime associated with the above fourth combination. The modeling of this regime will be addressed in MELPROG-PWR/ MOD 1. \\ The interfacial friction force between two fields is defined in MELPROG as =F (ap)g C AV lAV ikl (2 3) f ik g fik ik where F is a geometry factor (for example, F is 1/(2DH) IOP hydhaulic diameter), internaf-flow wall friction, where D is a g (op)g is fluid a macroscopic density, C is a friction factor rik or drag coefficient, AV is a relative velocity between the ik fields, and 2 is a unit vector pointing in the axial direction. The main difficulty here is obtaining the resistance coeffi-for the above six combinations. cients, Cfik, The procedure and software for determining the resistance coefficients for the first three combinations have been extracted directly from the TRAC-PF1 code. The interested reader is referred to Ref erence 2.1 for the details. The resistance coefficient for the fourth and fifth combinations is determined using drag coefficient correlations for a sphere moving in an infinite medium [2.23; this assumes the corium field consists mainly of solid particles in those cases involving either direct corium-water or corium-vapor contact. The resistance coet-ficients for the sixth combination is assumed to be zero in PWR/ MODO, thus neglecting the friction interaction between stationary fields and the corium field, respectively.

However, the DEBRIS module does account for the interaction of the corium with structures in that debris beds can form on horizontal surfaces.

The interfacial heat transfer term between two fields is defined with the following familiar form: 's (2.4) gig -Ag h AT g 3g ig 15

where A is the interfacial heat transfer area between fields i ik and k per unit

volume, h

is the heat transfer coefficient ip between fields i and k, and AT is the temperature difference gp between fields i and K. When indices i and k are equal then q becomes an power generation term (e.g., decay heating orchemickk reaction energy). It should be noted that the heat transfer terms from the fuel rods to the coolant are calculated by the PINS module and simply transferred to the FLUIDS module as an explicit term.

Likewise, the heat transfer from structures to the coolant is calculated in the STRUCTURES module and transferred to FLUIDS.

This procedure is necessary to ensure energy conservation between modules. Consequently, FLUIDS only cciculates the interfacial heat transfer for the first, fourth and fifth combinations above. However, FLUIDS does calculate the heat transfer coefficients for all of the above six combinations. This process is similar to the above determination of resistance coefficients. The procedure and software for determining heat transfer coefficients for the first three combinations is again extracted directly from the TRAC-PF1 code (see References 2.1 and 23 for details). The heat transfer coefficient for the sixth com-bination is currently assumed to be zero, thus neglecting in PWR/ MODO any treatment by FLUIDS of the heat transfer between stationary structure and the corium field. However, the DESRIS and STRUCTURES modules do account for heat transfer from debris to structures. The heat transfer coefficients for the fourth and fifth combinations are determined in the following manner. If the corium field is dispersed (i.e. low volume fraction), then the coefficients are calculated as though the corium consisted of individual spherical particles. If the corium is not dispersed, then the coefficients are determined using a flow regime map analogous to that in TRAC-PF1. The model considers the following five heat transfer regimes: (1) Single Phase Liquid Regime, (2) Nucleate Boiling Regime, (3) Transition Boiling Regime, (4) Film Boiling Regime, and (5) Single Phase Vapor Regime. In the first

regime, the debris-to-liquid heat transfer coefficient is determined from the correlation of Whitaker[2.4]

I and the debris-to-vapor coefficient is set to zero. This regime I 1 16

covers' the situations where'there is.no vapor present (i.e. up to the point where boiling begins). L In the. second regime,- the debris-to-liquid coefficient is ' determined' using the ~ correlation found in the TRAC-PF1 model. 'Again... the. debris-to-vapor coefficient is set to zero. The boundaries for this regime are from the point where boiling beg ins, up., to the dryout point. The. dryout L point is determined-using the Lipinski mode 1[2.53. o .The third regime covers the transition from nucleate boiling to film-boiling. The coefficients in.this regime are calculated using the correlations. found in the TRAC-PF1 model with the

effective diameter and Reynolds number being defined in an appropriate manner.

The boundaries for this regime are from the dryout point.up to the point where the debris temperature exceeds the minimum stable film boiling temperature. l In the fourth regime. 'the de b ri s-t o-li q u i d heat transfer coefficient is set to zero. The debris-to vapor coef ficient - is determined' ;using.the film boiling correlations' found in the TRAC-PF1 heat-transfer model. .The boundaries of this regime are from the minimum stable film boiling point up to the point where 'all the liquid is vaporized. In the fif th regime, the debris-to-liquid coefficient is set-to zero, and the debris-to-vapor coefficient is calculated from the Whitaker model. This regime is used only when there is no liquid present. 'The regime: mapping is relatively simple, but in view of the lack of good experimental data it seems to be adequate. Certainly, the division into regimes seems appropriate. However. -the modification and application of rod bundle correlations to porous media flow.is. questionable and will be improved in future versions. 2.4. Mass: Transfer Functions Before the conservation Eqs. 2.1 can be

solved, the are needed to describe the loss mass-transfer functions (f ijk) and gain of mass that occurs during phase

change, chemical reactions, or material fragmentation.

For example, there is an inj ection ' of hydrogen into the vapor field during the process of cladding oxidation, and a corresponding loss of water vapor from a the steam component of the vapor field. This H gain and H0 J 2 p loss is accounted for in the conservation equations by the mass transfer functions f ijk* j ] ) 17

Currently there are a total of nine unique mass-transfer functions in the MELPROG-PWR/ MODO conservation equations. These ares 1. H O vaporization and. condensation, 2 2. H generation due to metal oxidation, 2 H O vapor loss due to metal oxidation, 3 2 4 UO2 melting, freezing, and fragmentation, S. Zr melting, freezing, oxidation, and fragmentation, 6. Zr0

melting, freezing, formation by oxidation, and 2

fragmentation, 7. stainless steel melting and freezing, 8. control rod melting and freezing, and 9. eutectic formation and freezing. The mathematical form for each of these mass transfer functions is briefly described next. Note that stainless steel oxidation, which becomes important at high temperatures, is not included in MELPROG-PWR/ MODO. The form of the mass transfer function that models HO 2 vaporization and condensation is discussed in detail in Reference 2.1. Tnat discussion will not be repeated here. The H 0, Zr, and Zr0 mass parameters which determine the H 2 2 2, transfer functions associated with the metal-oxidation reaction are discussed in Section 3 The remaining mass transfer functions which deal with fuel or control rod melting, fragmentation, or freezing have the form (dropping the indices for convenience): T - m s(t) (2.5) i o where s(t) is a smoothing function that determines the rate at which the mass transfer takes place, m is the total quantity of o mass in a mesh cell that is involved in the mass transfer, and t is the time elapsed since the mass transfer was initiated and is also cell-dependent. For example, when the control rod mass mo

0) in a particular numerical cell in the reactor melts (at t =

core, Eq. 2.5 represents the mass transfer of this material into i 18

-the corium field in FLUIDS at a rate determined by the smoothing function s(T). The function s(t) is chosen such 'that its integral from zero to infinity is equal to 1.0. This smoothing helps eliminate numerical instabilities caused by rapid phase change. 2.5. Solution Technigues ,o The solution of the conservation equations represented by Eqs. 2.1 is achieved by using the Stability Enhancing Two-Step (SETS) numerical method [2.6]. This method has some similarity to the historic predictor-corrector class of. numerical solution techniques, and is designed to propagate the information needed for stability with minimal implicit coupling between spatial cells. The-practical advantage of using the Two-Step method is its inherent ability to exceed the material Courant limit by orders of magnitude in many situations. The interested reader is directed to Reference 2.6 for more detail. Ten of the fifteen conservation equations discussed in Section 2.2 are solved using the Two-Step method; these are the continuity, momentum, and energy equations for the vapor, water, The and corium

fields, and a hydrogen continuity equation.

remaining continuity equations for UO Zr, ZrO stain ~1ess 2, 2, steel, control rod material, and eutectic are solved in either an uncoupled or an explicit fashion external to the Two-Step method. An additional constraint is imposed on the solution of the conservation equations to simulate the formation of debris beds. I Once fragmented or molten debris is generated and begins to j relocate to some level to form a debris bed, the corium void fraction is partially constrained not to exceed a maximum packing l fraction in any given numerical cell. This constraint is imposed by monitoring the corium volume fraction for each cell and when it exceeds a specified limit (typically about 445 if the corium is solid), then the appropriate corium velocity is forced to

zero, implying that the cell is full of debris.

This action insures that a cell will not overfill, but still allows the corium volume fraction for the cell to change as the debris cools, melts, or fragments further. The reader is referred to Section 5 for a detailed description of the debris bed treatment. The actual maximum packing fraction allowed for a cell is determined by assuming a debris size distribution and using a mechanistic debris bed packing model[2.73. This model follows the dynamics of non-random packing of fuel and control debris around and above surviving rod segments, and on top of a porous floor (such as a grid-spacer or lower support plate). The i irregular shape of the debris particles is also accounted for in this model. The interested reader is referred to Reference 2.7 for the details. It should be noted, that if the debris is ( 19

J f molten, then a maximum packing fraction of 0.9 is assumed for the cell. 2.6 TRAC-MELPROG Link The TRAC-LINK methodology was developed to allow MELPROG to ~ be coupled to the entire reactor coolant system which can be modeled with TRAC-PF1. This coupling is only between the fluid dynamics in MELPROG and TRAC. That is, MELPROG acts as a core or vessel component for TRAC. Since the coupling is fully implicit, iteration between TRAC and MELPROG on each time step is needed to obtain a solution. Currently, the TRAC-MELPROG link only communicates infor-mation about the vapor (composed of steam and hydrogen) and water fields between the two codes. These are the two fluid fields which are modeled by TRAC-PF1. Other fluid components, such as the fission p ro d c. c t s and core debris, cannot be followed when they leave the HELPROG domain due to limi t a tior.s in TRAC. Proposed TRAC 1rprovements will allow the transport of all the MELPROG fluid fields and components. While running in a linked mode, the computational flow is controlled by the TRAC code. This means that MELPROG is treated as a one-dimensional core or vessel component. As with any one-dimensional compcnent, all information necessary for the network solution is communicated via the TRAC-PF1 boundary array. Tne MELPROG code is written such that a stand-alone MELPROG calculation can be run with no changes other than including a main program to start the computation. If a combined TRAC-MELPROG calculation is desired, one merely loads the two codes togethar and TRAC assumes control of the computation. When loaded in this fashion, a stand-alone TRAC calculation is also possible. The details of the computational flow are given in Appendix C. 2.7. Referencea

2.1. TRAC-PF1

An Advanced Best-Estimate Computer Program for Pressurized Water Reactor _Ana1ysis, LA-994#-MS, NUREG/ CR-3567, Los Alamos National _ Laboratory, Las A3amos, NM, (Feb. 1934). 2.2. J. W. Daily and D. R. F. Harleman, Fluid Dynamics, Addison-Wesley Publishing Company (1966). 20 t

1- <2.3.- D. A..Mandell, The TRAC-PD2 Boiling Curve, LA-UR-80-3679, ) Los Alamos National Laboratory, Los Alamos, NM, (January 1981). 2.4 S. Whitaker, " Forced Convection Heat Transfer Correlations for Flow in Pipes,-Past Flat Plates, Single Cylinders, Single Spheres, and for Flow in Packed Beds and Tube Bundles", AIChE Journal, Vol 18, No. 2,

p. 361 (1972).

• o 2.5.

R. J. Lipinski, A Model for Boiling and Dryout in Particle Beds, NUREG/CR-2546, SANDB2-0765, Sandia National Labora-tories, Albuquerque, NM, (1982). 2.6. J. H. Mahaffy, A Stability Enhancing Two-Steg Method for One-Dimensional Two-Phase Flow. LA-7951-MS, Los Alamos t National Laboratory. Los Alamos, NM, (1979). 2.7. J. A. Moore, A Mechanistic Debris Bed Packing Model, Nuc1 Technol, Vol 67, p. 66 (1984), i 1 .o. I 21

1 3 FUEL AND CONTROL RODS r-l 3

• 1 '-

1Elngdugligg The. PINS module in MELPROG contains the fuel and control rod behavior models. These models calculate heat transfer and morphological changes for the fuel rods as well as the melting of the control rods, PINS interacts with the other MELPROG' modules, principally the FLUIDS, RADIATION, and DEBRIS modules, through various heat and mass transfer terms. The various fuel rod mor,phological f orms that should be considered in a fuel rod behavior module are listed and described in Table 3 1. This list is based on various experimental observations [3 1.3.23. Only the first five of these are included in the present PINS module, which was inherited from the MIMAS code [3 33 and has' remained virtually unchanged since then. It is worth noting that MIMAS was designed to analyze TMI-2 type accidents which involve fairly slow core heating and reflood. Consequently, some of the potential fuel rod failure mechanisms were treated in a very approximate manner. The models will be extended in MELPROG-PWR/ MOD 1 and.other future versions of MELPROG. As mentioned in Section 2, MELPROG divides the core region into a number of axial and radial numerical cells. Within each

cell, only one representative fuel rod as well as one representative control rod is modeled.

This implies that within-that cell 'no temperature distribution among the fuel rods is calculated; that is, all rods behave as the average rod. In analyzing fuel rod behavior, a radial heat conduction calculation is performed in which both axial and radial core power and burnup distributions are considered.

However, axial conduction heat transfer is ignored.

In the current version, there is no heat transfer calculation performed in the control rods. In the following sections, the physical modeling included in the PINS module is discussed. In addition to this, the numerical methods used to calculate-the heat transfer in the fuel rods are described. i 3.2. 3pdg1_Degsgig11ggg This section describes the principal models used to analyze fuel and control rod behavior. Only the models which are currently included in the code are presented. l 22

Table 3 1 Fuel and Control Rod Section Morphological Forms Form Description Intact Intact geometry with only thermal expansion considered Ballooned Breached cladding with gross cladding expansion and possible fuel pellet disruption l Contained' Liquified fuel ' contained within Zr02 shell Eutectic Embrittled 8-Ir contains sufficient oxygen such that embrittlement criteria are met { j Disintegrated Fuel rods transformed to debris upon . quenching or melting Atrophied Partial pellet stack which remains after liquified fuel and cladding drain away Not included in MELPROG-PWR/ MODO Perforated Perforated cladding surrounding atrophied pellet stack and say.be emitting liquitied fuel Not included in MELPROG-PWR/ MODO Candled Fuel rod coated with molten or frozen debris Not included in MELPR00-PWR/ MODO 1 23

a i I ) 3 2.1. Cladding Rupture and Ballooning Models l i The ballooning behavior of the fuel rod cladding is determined in MELPROG from the model in NUREG-0630[3.4], which appears to be the most complete of those models that are available at present. This model is based principally on data i obtained at Argonne National Laboratory and Oak Ridge Natjonal Laboratory in single-rod burst tests. The NUREG-0630 model has been compared against two other models to estimate the range of uncertainty of the data used in the models. These other models are the REBEKAC3 2] model, based on data obtained at KfK in Germany, and the MATPRO-11[ 3 1] model. The ballooning rupture model consists of two parts: a correlation for burst temperature versus burst

stress, and a correlation for burst hoop strain versus burst temperature.

Comparisons of the correlations for each of the three models are shown in Figures 3 1 and 3.2. In the PINS

module, application of the NUREG-0630 correlations involves a two-step procedure.

In the first step, the cladding stress is compared to the burst hoop stress from the first correlation. If the cladding stress exceeds the burst

stress, then the cladding is assumed to have ballooned and ruptured.

In the second step, the second correlation is then used to determine the strain at rupture. The strain at rupture is multiplied by 0.56 to account for interference effects between the adjacent fuel rods present in a real rod bundle that are not accounted for in the single-rod correlations [3.4.3 5.3.6]. The resulting multiple-rod strain is then used to determine the new channel flow area in the ballooned axial section. An additional restriction is also imposed which states that the total flow area cannot be less than 5% of the original area [3.4]. Although it is recognized that azimuthal temperature j variations can considerably reduce the burst strain (to 40% at AT of 100 K), this effect is not included in the current MELPROG model. 3.2.2. Zircaloy Oxidation, Hydrogen Generation, and Embrittlement Zircaloy oxidation and embrittlement is modeled in PINS by using a version of the COBILD routine from the MATPRO library [3 1 ]. Figure 3.3 shows the cross section of a fuel rod that has been oxidized in excess steam above 1200 K for a period of

time, and that his also experienced fuel-cladding c on t a c t [ 3.1 ].

At this temperature, the unaffected zirconium 1 f (shown in Figure 3.3 as the wide center layer) has transformed to } the relatively ductile 6-Zr phase. The outermost layer has oxidized to Zr0 Immediately adjacent to this layer, another 2 24

W l mnuevar un arc - run 9 00 ,= tup 0 - F900-4 f asaru I- %,s 900 'N,' ~ esatc. SMD 700* wtumo-e 600 20 M 40 50 60 70 40 90 900 110 O to StatST pettsuutt (tae) Figure 31 comparison or Burst Temperature / Pressure correlations at Low Temperature Ramp Rates. 4 CLAD AT - 0* K LOW HEAT RATES u0-00- /~h no- [ Katu ,,o. g, ,/ WTHto { '0 " g eo-

• ?0 -

a c. l: \\ / ou0 ,0 /. 40- / \\&/ \\ 50 \\ to- \\ to d 0 600 700 600 900 1000 floo l most tie.*c l l [ Figure 3.2 comparison or Burst Temperature / Strain correlations at t Low Temperature Ramp Rates. l 1 i i 25

e - ~

2r0,

.-reto) f e-xnou

• -3d*k

_. (U.Zr) i I l Figure 3 3 Cross Section of Fuel Rod with Oxidized Cladding [3 2] l 1 i I l 1 1 l i l 26 '1 i

layer exists in which a high oxygen concentration has caused the zirconium lattice-to assume the low-temperature alpha-configuration. Oxygen diffuses through both of these outer layers into the B-Zr region. .Since the 8-Zr layer is the most ductile of all those shown, this region is assumed to carry all of the cladding load. As the 8-Zr layer absorbs oxygen by diffusion, it becomes embrittled, eventually breaking up under modest thermal sheck loadings. The configuration of the cladding layers adjacent to the fuel is currently open to some question. The inner cladding layer configuration assumed in the current modeling (and shown in Figure 3.3) consists of a metallic uranium-zirconium layer sandwiched between two oxygen otabilized a-Zr layers [3 1]. The oxygen needed to form these layers is assumed to be supplied by the uranium dioxide fuel. Assuming availability of an excess steam concentration. the growth rates of the Zrop and oxygen-stabilized a-Zr layers follow parabolic kinetics, that is..ts16,68.16 =fexp (-B/RT) (3.1) i where L= layer thickness A = cor.stant B = activation energy R = gas constant T = absolute temperature i time. t = Values for A and B are given in Table 3.2. The values for temperatures below 1850 K are taken from Reference 31. For l temperatures above 1850 K, the correlation of Urbanic and Heidrick[3.7] is used. i The growth rate of the Zr0 layer is important because it 2 determines the hydrogen production rate according to the chemical t equations i 2H 0+ Zr = 2H2 + Zr02 Q (3.2) 1 + p where 4 Q = heat of reaction = 6.5 x 10' J/kg Zr. 27

l' Table 3 2 .Zr Oxidation Rate Constants Layer A B ~ Zr02 1270 K < T < 1850 K 1.126 x10-6,2/s

1. 5 02 x105 afmot T> 1850 K 2.07 4 x10-6,2/s 1 3 3 2x105 afmot Oxygen Uptake 1270 K < T < 1850 K 16.8 (kg/m )2/s

1. 6 6 8 x105 afmot 2

T > 1850 K 30.9 (kg/m2)2/s 1 381 x105 afmot m /s

2. 0 5 x105 afmoi 2

Outer a-Zr

7. 615 x10-5 Inner a-Zr 7.0x10~S m /s

1. 8 4 x105 afmot 2

i a-Zr Nex t to Fuel 3 2x10-5,2 /s 2.05x105 afmoi-l k i l l 28 2

l For Qituationo whore otoichiomotric or exceco steam is available, the heat of the steam-zirconium reaction can be much greater than the fission product decay heat of the f u e l,.

However, the reaction can be limited if the steam partial pressure is too low.

In this case, the reaction is said to be " steam-starved". In .MELPROG, if the steam mole fraction is less than 0.01, the oxidation reaction is assumed to be steam-starved and the reaction ceases. By combining Equations 3.1 and 32, the H production rate 2 from Zircaloy oxidation as a function of temperature is calculated in MELPROG. This H is added to the vapor field in j 2 1 FLUIDS by using a mass transfer function. Since the brittleness of B-Zr is correlated with its oxygen i

content, it is necessary to

' calculate the rate of oxygen diffusion into the 8-Zr cladding layer. In MELPROG, this is done' i by numerically solving the Fick's law diffusion equation for the j cladding geometry, j i D(d ) (3 3) + l I 1

Here, C

is the oxygen concentration, D is the diffusion l coefficient for oxygen in zirconium, t is time, and r is radial l distance into the B-Zr layer. The value of D, as well as the exact numerical scheme for solving Eq. 33 are taken from Reference 3 1. When fuel-cladding contact occurs, the boundary conditions used in solving Eq. 33 assume high oxygen concentrations on both sides of the B-Zr layer. Otherwise, oxygen is assumed to diffuse into the B-Zr only from the outer surface. 3.2 3 control Rod Meltdown Models j The PINS module currently uses a rather simple approach in modeling the meltdown behavior of the stainless steel

clad, Ag-In-Cd control rods which are contained in Zr guide tubes.

In the model, temperature is considered the controlling parameter in determining the meltdown behavior. As

such, the following thermal conditions are assumed in the modeling approach,

First, the temperature of the Ag-In-Cd absorber

material, stainless steel cladding, and Zr guide tube are considered uniform within i

the entire core cell and are taken as equal to the local steam temperature.

Second, the melting points of the various constituents (median value of the solidus-liquidus) are 1070 K

] l i I 29 i l l

F for the Ag-In-Cd alloy and 1700 K for the stainless steel. Finally, the beta-to-alpha and Zro oxidation products in the Zr p guide tubes are modeled in a manner similar to the fuel rod cladding oxidation, but the guide tube temperature is also assumed to be equal to the local steam temperature. For the assumed condition of uniform temperature, complete 1070 K) is melting of the Ag-In-Cd absorber material (T = predicted prior to melting of the stainless stee3 cladding (T, - 2 1700 K). However, until cladding rupture occurs via melting, the molten Ag-In-Cd is assumed to remain " bottled" within the intact cladding. This assumption is supported by studies of the Ag-Fe system that indicate that As and, Fe are immiscible in the liquid state and therefore exclude allocability (see References 3.8 and 3 9) and low temperature eutectic formation. Cladding rupture via melting can therefore be expected to occur at or near the stainless steel melt temperature. For the present, the cladding melt criterion is considered the only mode of cladding failure. For the core axial cell where the cladding melting point is first predicted, local cladding failure is

assumed, at which time release of the " bottled" molten Ag-In-Cd within that cell and above cells (which contain molten Ag-In-Cd) occurs.

This breach / release mode is one of catastrophic failure, where no resistance to molten Ag-In-Cd release is offered by either the cladding or guide tube. Upon control rod cladding failure, the mass of the molten control material, cladding, and guide tube are transferred to the debris field in FLUIDS by means of a mass transfer function. Although the molten Ag-In-Cd may interact with surrounding fuel r o d s [ 3.10, 3.113, no further interactions of the molten control material is considered in the present version of MELPROG-PWR/ MODO. 3 2.4 Fuel-Cladding Liquefaction Phenomena layer and the fuel pellet When the material between the Zr02 melts, the molten cladding can dissolve a portion of the UO 2 pellet and form " liquefied fuel". The subsequent behavior of the section of the fuel rod containing liquefied fuel is a function of the prior history of the rod, and is generally known only qualitatively. However, since the behavior of fuel rod sections containing liquefied fuel is important to the overall reactor response to a degraded core accident, models for liquefied fuel phenomena are included in MELPROG, with the understanding that these models will be updated with new data. In the PINS module, the models for liquefied fuel phenomena assume that the oxide layer existing on the cladding outer surface contains the molten cladding / liquefied fuel until the 2 by the available molten cladding is dissolution of the UO complete. The solubility of UO in the molten cladding is given p by the Zr-UO binary phase diagram shown in Figure 3.4, which is 2 30

l 1 l I 3400 h + We sapo-Q 3000-Li v L+b , anoo-s h D m-1 u. u,w s z4oo-g, ' cu=N i E E m-k asoo o 02 e4 os os Mole Fraction UO Figure 3.11 Zr-UO2 Binary Phase Diagram.

• d>

l l 31 l ) i J

l based on data by Politis[3 12]. In this diagram, the variation } in the melting point of zirconium with oxygen content is I represented by the solid and dashed phase boundaries between zero' l and 0.05 U0 mole fraction and 2100 to 2200 K. If the fuel and 2 cladding are in contact, the cladding is assumed to contain sufficient oxygen to elevate its melting point to 2200 K. However, the melting of the U-Zr metallic layer, which occurs far below 2200 K, is not modeled. The liquidus temperature increases with increasing U0 mole frection up to approximately 2670 K, 2 above which two liquid phases are formed. Above this is approached, so the temperature, the melting point of ZrO2 layer is no longer. valid. initial assumption of a containing Zr02 Heat-of-fusion effects are included for the melting of the non-oxidized cladding, but the heat of solution of UO in molten 2 Zr is not included at this time. Oiven a volume of molten Zr in a particular cladding element and the temperature of that

element, the data in Figure 3 4 are used to calculate the fraction of the pellet radius that would dissolve in the molten cladding at equilibrium.

It is then assumed that this equilibrium amount of fuel is dissolved into the molten Zr. This process is repeated on each time step so that the liquefied fuel volume in the cladding segment is calculated as a function of time. 3 2.5. Fuel Rod Failure and Debris Formation In 'the PINS

module, the failure of the fuel rods is considered to occur by one of four mechanisms.

The failure mechanisms are based on experimental results and can be classified as: 1. Thermal Shock Failure, 2. Fuel-Cladding Contact Failure, 3 Fuel-Cladding Eutectic Formation, and 4 Complete Melting of the Cladding. The meaning of each of these failure mechanism is discussed below. Based on the MATPRO correlations, two types of thermal shock failure are considered in the PINS module. The failure mode depends on the cooling rate. The two types of cooling rates are fast cooling and slow cooling. Fast cooling is defined as a quench of the cladding from a temperature above the a-S Zr direct phase transition temperature with cooling rates of 100 K/s or more. Slow cooling occurs when the cladding is cooled at about j f 32

5 K/s to a temperature below the phase transition temperature, ) l and then quenched. l By carefully correlating analytical results with observed cladding behavior for these cooling rates, a set of criteria has been constructed relating 6-Zr oxygen content and temperature history with breakup during coo?ing. It was found that under fast cooling conditions, the cladding broke up if I 1. the cladding temperature ever exceeded 1700 K, and 2. oxygen concentration in the 6-Zr exceeded 905 satura-tion, or if 3 oxygen concentration in the B-Zr exceeded 65%, by weight. The cladding was found to shatter under slow cooling conditions if at least half of the cladding thickness contained more than one weight percent of oxygen. These criteria and the calculated oxygen content of the 6-Zr cladding layer are used in PINS to determine the fuel rod breakup by thermal shock. The second failure mechanism results when there is fuel-cladding contact and the cladding is embrittled. In this case the fuel would have to swell against the embrittled cladding to cause failure. As implemented in the PINS

module, the cladding must be more than 80% oxidized when the contact occurs in order to have failure.

Otherwise, the cladding is assumed to plastica 11y deform. The remaining two failure mechanisms occur at high temperature and deal with melting phenomena. These mechanisms have been formulated by using the experimental observations of Hagen[3.13,3 143. From these observations the following scenarios have been proposed. cladding 1. With rapid fuel rod heatup rates, the Zr0 2 layer is relatively thin when the liquefied fuel forms. The volume expansion of the liquefied fuel apparently breaks up the ZrO

layer, followed by downward 2

relocation of the liquefied fuel and ZrO2 fragments. An atrophied fuel pellet stack may remain behind. cladding 2. With moderate fuel rod heating rates, the ZrO2 layer is somewhat thicker and stronger during the liquefied fuel formation. The liquefied fuel leaks out and downward through perforations in the oxide layer, leaving behind an intact, but badly damaged, fuel rod structure. 3 With low fuel rod heating rates, the oxide cladding I 33 1

l i layer becomes thick enough to contain the -liquefied )

fuel, so that the oxidized fuel rod remains intact throughout the transient.

The third failure mechanism attempts to simulate the first scenario. This mechanism involves fuel-cladding eutectic formation at a time when the oxide layer is thin. As implemented ~ in PINS, the oxide layer is considered to be thin if it is less than half the thickness of the original cladding. The fourth . failure mechanism accounts for the third scenario. In this case, the oxidized cladding is assumed to be strong enough to support the liquefied fuel expansion and failure As only occurs by the melting of this heavily oxidized cladding. . implemented in PINS, if the cladding is more than 50% oxidized, then the cladding will not fail until it melts or experiences a thermal shock. When failure occurs, the current version assumes complete disintegration of fuel rods in the given cell. This assumption precludes a candling type of failure which, however, will be included in MELPROG-PWR/ MOD 1. The mass of the failed rods is transferred to the debris field in FLUIDS. Also, the composition 18 of the debris mass, which consists of Zr metal, ZrO, and_UO2' 2 calculated and transferred to FLUIDS. It is assumed that the distribution of the debris particles from the disintegrated fuel rods can be represented by an average particle size. For the thermal shock failure mechanism, a particle size of 0.002 m is For the other failure mechanisms, a particle size of used. 0.005 m (approximately one-half the fuel rod diameter) is used. 33 Fuel Rod Temperature Analysis Model Subroutine TEMCALC, which is called from subroutine TMPREP, calculates the temperature field in intact and ballooned fuel rod elements based on a one-dimensional radial finite-difference approximation to the general conduction equation, (3.4) f[ (r k )) + q = po p The finite difference equations are solved in a conservative manner, sequencing from left to right, where the left boundary corresponds to the smallest radius. Nodal points are positioned l on material interfaces and material properties are evaluated l ~ between nodes (see Figure 3 5). l 34 l 1

l. ' l j som h + Uc l 3a00-Q 3000-Li i v L+4 , smoo-k Q Noo- \\ Le + (Wtrka 5 g 2400-(U2' N E " ~ $ 802-k 3a00-.-2stobt e-h9) + UQn seco-o em a4 os os Mole Fraction 00: Figure 3.4 2r-Uo Binary Phase Diagram. 2 k ) { i I s l 31 J

5 K/s to a temperature below the phace transition temperature, and then quenched. By carefully correlating analytical results with observed cladding behavior for these cooling rates, a set of criteria has been constructed relating B-Zr oxygen content. and temperature history with breakup during cooling. It was found that under fast cooling conditions, the cladding broke up if l 1. the cladding temperature ever exceeded 1700 K, and l 2. oxygen concentration in the B-Zr exceeded 905 satura-tion, or if 3 oxygen concentration in the 6-Zr exceeded 65%, by I l weight. I The cladding was found to shatter under slow cooling conditions j I if at least half of the cladding thickness contained more than one weight percent of oxygen. These criteria and the calculated oxygen content of the B-Zr cladding layer are used in PINS to determine the fuel rod breakup by thermal shock. l The second failure mechanism results when there is fuel-cledding contact-and the cladding is embrittled. In this case the fuel would have to swell against the embrittled cladding to cause failure. As implemented in the PINS

module, the cladding must be more than 80% oxidized when the contact occurs in order to have failure.

Otherwise, the cladding is assumed to plastica 11y deform. The remaining two failure mechanisms occur at high temperature and deal with melting phenomena. These mechanisms have been formulated by using the experimental observations of Hagen[3 13,3 143. From these observations the following scenarios have been proposed. 1. With rapid fuel rod heatup rates, the Zr0 cladding 2 layer is relatively thin when the liquefied fuel forms. The volume expansion of the liquefied fuel apparently breaks up the Zr0

layer, followed by downward 2

fragments. An relocation of the liquefied fuel and Zr02 atrophied fuel pellet e;.ock may remain behind. cladding 2. With moderate fuel rod neating rates, the Zrop layer is somewhat thicker and stronger during the liquefied fuel formation. The liquefied fuel leaks out and downward through perforations in the oxide layer, leaving behind an intact, but badly damaged, fuel rod structure. 3 With low fuel rod heating rates, the oxide cladding 33

l l ) l i l C.L r, ~t l l Art Ar. 4 n ei a i < >i ,i,< i ja V 1 / \\ l Node l 4 J \\ ) r-i Cap Figure 3 5 Finite-dif ference Noding Scheme f or MELPROG Fuel Rod Thermal Model e 35 l

L For the innermost fuel pellet node (1-1), the finite difference equation-is 2 ar (pe p) k k 3 g, + y (r,4r; + _g_ ) __gg__1 ) r, + 1 +_"_3 /2 [_P_3 / 2 3 1 n 7n+1 g, 2. ^((]((f[{21ri+qp*1] (r3 Ari (3 5) = -- with the boundary condition - rk II l =0 (3.6) dr g,3 at the pellet inner radius. surface. In these equations, material thermal conductivity in node 1 k = g = product of material density and specific heat in node 1, q'p)g (po internal heat generation in node 1, = i time step size.. at = n = time step level, i.e., n = present time, n+1 = present time + At, temperature of node 1, and T = g Arg = mesh spacing for node 1. For interior nodes (1 < 1 < NF) in the fuel pellet, the finite-difference equation is r k r ~+1, { r i-1/2_k k Pi-1/2 g.1 n g.j + __1_+ 1 / 2 g g arg.3 Arg.3 org 2 Of 1 i 1 777((rg Arg_3 -_g_-) (pop)g.g + e w 36

l f i-, (rg'Arg!+ ) (pop) g') } 'T[ - -5r + + i ~ (rg Arg.g )s { Ty+q[M] =- 4 (rg Arg + ).[ Ty+q[*I) (3.7) + where NF is the number of nodes in the fuel. The' gap that exists between the fuel and the. cladding in fuel rods is treated by explicit noding on fuel and cladding-surf aces and a heat-transfer coefficient between these nodes. Stored energy and internal heat generation.in the gap region are . neglected. The finite difference equation for the outermost fuel

pellet node (1-NF) is k NF-1 rNF-1/2 kyp.1 Tn+1

( rNF-1/2 + ryp hgap yp.3 NF-1 aryp.g AP NF-1) (PCp)NF-1) Tn+1 AP 1y (ryp Aryp.g 4 6t yp h Tyh3

• PNF gap ArNF-1) ( -(pop)NF-1 1

In y (rNF APNF-1 ~ 4 At NF qQ*l.g) (3.B) + and the boundary condition at the outer pellet surface is 37

- ~ k l =hgap (TNF TNF+1) (3.9) is the conductance of the pellet-cladding gap (or

Here, h *[r f a c e ).

s The method used for determining the magnitude contact of h is discussed in the next section. gap In,the cladding region,_the internal heat generation rate is redefined to include metal-water reaction 'and radiation heat transfer when appropriate. For the innermost cladding node (i=NF+1), the finite difference equation is rNF+ 3 /2 k yp.3 h Tn+1 [P NF+1 h ~ P NF+1 gap NF Eap NF+1) ~(DCp)NF+1 P 1 ) Tn+1 + y (rNF+1 APNF+1.* NF+1 ~ a ag k PNF+3/2 NF+1

• 1 7"NF+2

) 3P NF+1 + ~~ F+1) [(@Cp)hF+1 ArN 1 ~ y (ryp,j ArNF+1 Tn g 6t NF+1 +qypi] (3 10) Tne finite difference equation for interior cladding nodes (NF+1) < 1 <N is identical to Eq. 37. For the outermost cladding node (i=N), the finite difference equation is kN-1 PN-1/2 NN-1 I.g [rN-1/2 ~ ArN-1 N-1 ap l (P APN-1 g N ) + g))Ty+1 ~ + f PN (hg +h as 38

t i -j (r, ar,_t - '"h') I Ejh-1Ty+q"'g*j) = y_-T[1)) -(3.11) t

• I ),+ h (f

T + ry [hg (f TN~Tg g t i. with'the boundary. condition at the outer surface given by l-l l'l -k l = hg (Ty.- Tg) + h ~(TN-Tg) (3 12) g i=N -l where .j -hg --heat transfer coefficient to liquid, i h, = heat transfer coefficient-to gas, t-liquid temperature, I = f f,: = gas temperature, T = steady-state flag (1 for steady-state calculation, O otherwise),'and transient flag (1'for transient calculation, f = t O otherwise). 0), a. Note that for a steady-state. calculation.(f 1 and f t (=f = ss fully' implicit' form is used. Under transient conditions =0 as t - a semi-implicit form is obtained. l}' and f i 1 tridiagonal Ty*3 [the-f amiliar These equations can be written 3 15]. This system of form for.the nodal temperatures Oguations can be solved for the temperatures T"*I by a variety of tridiagonaf matrix inversion nethods. In MELPR00, a standard scheme is-used. l and Cladding _ Stress 3.4 Materials Progerties _ Gap _ Conductance 2 2 This section discusses a number of models whose results are used by various subroutines called from INTACT. Most of the material property models are for thermal properties, and therefore supply. values for the finite-difference equations in the TEMCALC routine. Other models for fission gas

release, cladding plenum

pressure, and differential fuel-cladding expansion are used to calculate f uel-cladding gap conductance and cladding stress.

l Models for the fuel and cladding material thermal properties .in. MELPROC are from Reference 31. The models for UO and 4 - Zircaloy cover a temperature range extending beyond thesaferial i 1 39

l' h l The thermal conductivity-for Zrop.as a function L celting' points. of. - tem pera tur e was also taken from Reference 34 the specific density of Zr0 in MELPROG were derived from; data heat.and Ref erence. 3 16.2 used As a first approximation, the ? at of given.'in fusion off Zircaloy is approximated in the TEMCALC analysis by assigning an artificially high specific heat to the cladding when the. melting point is reached. This artificial specific heat equals one tenth of the Zircaloy heat of fusion. When the molten i Zircaloy -(with the high specific heat) reaches a temperature 10.K above..the melting point, the actual molten Zircaloy specific heat value from Reference 34 is used. Best-estimate values for the l-thermal properties of liquefied f uel are assumed to equal the ~ cedian between the corresponding molten zirconium and molten UOp 'i values. To calculate the size of the fuel-cladding

gap, the differential fuel-cladding - expansion is calculated in MELPROG using the models from Reference 3.4 for fuel and cladding thermal expansion and fuel-fission product swelling as ~ a function of temperature.

Once the gap size is calculated, subroutine GAPCON (h as a function' cond uct an cka p ) calculates the gap heat-transfer coefficient fuel-cladding of three components: gap. gas interfacial contact, and fuel-cladding thermal radiation. The equation used is h =h +h (3.13) gap gas radiation where h =kgas/(rgap + 6) (3.14) gas oF (Tf - T$) (3 15) h r a d i a t i o n " "~~(~T[~- Tc) ~ and 1)] (3.16) F= 1/( 1/cf + (R /Re)(1/c f c Subscripts f and c refer to fuel and cladding, respectively, e is 40

I 1 10-gs the inner F to the-Stefan-Boltzman constant, and the quantity m is used for i outer cylinder view factor. A value of 4.4 = 6, which includes the mean surface roughness of the fuel and cladding. The values for k are calculated from the GTHCON l as codel in Reference 3.4 whjch includes the effect of conductivity I .] degradation for small gaps. The amount of fission product gas released from the. fuel is calculated with the FGASRL model of Reference 3 4. To calculate j l the total gas pressure in.the fuel rod, the MELPROG code adds j l this released fission gas to the initial rod fill gas, which is cssumed to _be 3 2 MPa at room temperature. The fuel rod gas pressure is calculated as a-function of time using the ideal gas i law. 3 l l The cladding stress is calculated with the thin cylinder j j l formula l Pd (3 17) a = y-1 where, o is the cladding stress, P is the difference between the rod internal gas pressure and the coolant pressure, d is the mean rod diameter, and t is the 6-Zr layer thickness. Cladding stress caused by fuel-cladding contact is not calculated explicitly. Whenever the fuel outer diameter is calculated to be greater than the cladding inner diameter, the cladding is assumed to be pushed out by the fuel, provided the cladding does not fail via one.of the failure mechanisms given above. 3 5. References 31, MATPRO-7ersion 11 (Revision 1)2 A Handbook of Materials 1 Properties for Use in the Analysis of Light Water Reactor j Fuel Rod Behavior, NUREG/CR-0947, TREE-12E6 (1956). j 1 ) 3 2. F. J. Erbacher, " LWR Fuel Cladding Deformation in a LOCA and Its Interaction With Emergency Core Cooling", Proc. / o[ ANS Tggical Mpeting on Safety Aspects g[ Fuel Behavlog, i Vol. 2, pp. 100-113, Sun Ve11ey, Idaho ( August 1981). j i t 33 P. J. Maudlin et al, Light-Water Reactog Degraded-Core i c, ogling Pgggtam Technical Note A Damage Assessment of l TMI-2, LA-UR-63-737, Los Alamos National Laboratory, Los Alamos, NM, (1983). i ( 41 I i

1 3.4 D. A. Powers and R. 0..Meyer, Cladding Swelling and i Rupture Models for LOCA Analysis..NUREG-0630 (1980). j l 3.5. A. W. Longest, Multired Butst Test Program Progress Report j 11{1, N U R E G/ C R-2 3 6 6 vol. 1, 10C faggarr _dgge t ORNL/TM-605I, Oak Ridge Nationti Laboratory, Oak Ridge, TN, (December 1981). 3.6. J. L. Crowley, Multirod Burst Test Program Progress Report f0C EMll ~_EfEEE51P 1961, NUREG/CR-2356 Vol. 2, ORNL/TH-6190, 03k Ridge National Laboratory, Oak Ridge, TN, (llarch 1982). 3 7. V. F. Urbanic and T. R. Heidrick, "High Temperature i Oxidation of Zircaloy-2 and Zircaloy-4 in Steam", J. of d ~~ ~~~ Nuc. Mat. 75 pp. 251-261 (1978). j l 3.8. M. Hanson, Constitution of Binary A11 ors, McGraw-Hill, New l York (1958). j i 3.9. R. P. Elliott, Constitution of Binatr A11ogs : 1st f Supplement, McGraw-Hill, New York (1965). j i 3 10. G. W. Parker, G. E. Creek, - and A. L. Sutton, " Aerosol 1 Release Transport", ORNL Quar _terly Report for Jult-Segt. l 1981, Oak Ridge National Laboratory, Oak Ridge, TN. 3 11. G. W.

Parker, G.

E. Creek, and A. L. Sutton, " Aerosol Release and Transport", ORNL Qgattet g Report for Oct-Dec. i 1981, Oak Ridge National Laboratory, Oak Ridge, TN. 1 1 3 12. C. Politi s, Unt ers uellungen in Dreishoffsgstem Uran-Zirkog-Sauerstoff, Kf K 2157 Kernforschungszentrum, Federal Republic of Germany, (1979). 3 13 S. Hagan and H. Malauschek, " Bundle Experiments on the Meltdown Behavior of P.W.R. Fuel Rods", Transactions of the American Nuclear Society 11, pp. 505-506, San I Fra-ncisco, CA (1979). 3 14 S. Hagen, "Out-of-Pile Experiments on the High Temperature Behavior of Zry-4 Clad Fuel Rods," Proceedings of the 6th International Conference on Zirconium in the Nuclear IndustrI, Vancouver, D. C., Canada (1982). ~~~~~~~ ) 3 15. B. Carnahan, H. A.

Luther, J.

O. Wilkes, Agglied Numet cal i Methods, John Wiley and Sons, New York (1969). 3 16. R. C. Weast and M. J. Astle, eds., CRC Handbook of ChemistrI and Physics, CRC Press, Boca Raton, FL (1979). l l i 42

1 .4 STRUCTURES I l 4.1. Introduction 1 The STRUCTURES module consists of two maj or parts: heat transfer. and mechanical models. A description of the heat i transfer models is given in Sections 4.2 through 4.7 and the mechanical models are discussed in Sections 4.8 through 4.10. A brier. description of the major subroutines in the STRUCTURES module is given in Appendix C. 4.2. Structure Heat Transfer The structure heat transfer module is intended to calculate all heat transfer processes for the reactor vessel and all in-vessel structures excluding the reactor core (i.e. fuel and control rods). The intent is to treat both PWR and BWR designs; however, MELPROG-PWR/ MODO is limited to PWR designs. The module has been designed to treat both one-dimensional and two-dimensional conduction heat

transfer, ablation and bulk

melting, crust formation and bulk freezing, and plateout and deposition of fission products and aerosols.

In the PWR/ MODO version, one-dimensional thermal conduction as well as ablation and bulk melting are treated. Crust formation and freezing of bulk fluid components onto structures are currently under development and will be incorporated in later releases of MELPROG. Two-dimensional conduction heat transfer will also be provided for certain types of structures in the future. The l framework and all necessary input are included in PWR/ MODO for both crust formation and two-dimensional heat transfer. Boundary i conditions to the one-dimensional fluids (steam and/or hydrogen, j water, and corium) module, to the radiation heat transfer module, j and to the detris module are also treated in the current version. l 1 4.3. Stru_gt_gg T.2ge s 1 The in-vessel structures are categorized as walls, plates, j and special structures as listed in Table 4.1. Walls in a PWR I are the core baffle, the core barrel, the thermal shield, and the vessel wall. PWR plates include the veseel top and bottom heads, the lower core support plate, the lower core grid plate, the upper core grid plate, the upper support

plate, and others.

Special structures currently include columns and a general structure. Each of these categories can be expanded to treat additional specific structures, for example BWR structures. . Typical plate and wall type structures are shown in Figure 4.1. 's ^ 43

Table 4.1 Examples of Structure Types Wall Plate Column General Core Baffle Vessel Bottom. Support Column Upper Plenum Structure Core' Barrel-Lower Core Guide Tubes Steam Separators Support Plate Thermal Shield Diffuser Plate Vessel Wall Lower-Grid Plate Upper Grid Plate Upper Core Support Plate Vessel Head i a 1 I l 1 l l 44 i i I t

i ur < Idi d -r Z itsep/ O F< h Y e a 1 v M rN< L ta k a da3 1 32 I N II E S E 7V13W N g N 3SVG 1 M e ~ },. ) 45

Structures are set-up on a two-dimensional cell basis. Currently a maximum of five structures can be input in a particular cell and a maxieum of six structures can be treated. The extra space is provided to treat movement of failed l structures from one cell to another. The limits (5 input and 6 total) can be easily changed and will be made part of the input in future versions of MELPROG. 4.31. Wall Type Structures All current wall structures are annu11. These may be placed at either the inner, the outer, or at both the inner and outer radial surface (s) of the computational cell as illustrated in Figure 4.2. Surface 1 of a wall always faces toward the cell center and always uses the fluid conditions (i.e. the temperature and heat transfer coefficient) in the cell to formulate its boundary conditions. Surface 2 always coincides with the cell boundary and always faces the adjoining cell or the problem boundary. If surface 2 faces the adjoining cell, then the fluid conditions in the adj oining cell are used to formulate its boundary condition. If surface 2 is at the problem boundary, three optional boundary conditions are available. First, it may use the fluid conditions in its cell to formulate its boundary condition. Second, it may use an adiabatic boundary condition. Third, it may use an input, constant sink temperature with an input, constant heat transfer coefficient. If the radiation heat transfer module is activated, then radiation heat transfer is calculated for both surfaces unless surface 2 is adiabatic. Currently walls do not have holes in them. 4.3 2. Plate Type Structures i All plates Ere assumed to be horizontal and flat for heat transfer calculations. Plates may or may not have holes depending on the input. In PWR/ MODO all holes in plates are assumed to be round. In th7 future, r e ct an gul ar holes will be treated explicitly. A plate can be input at the bottom, the top, or at both surface (s) of a computational cell. Surface 1 of a plate always faces into the cell and always uses the fluid conditions in the cell to formulate its boundary conditions. Sv.rface 2 always faces the neighboring cell or the problem bounac y and is located on the cell boundary. As in the case of walls, if surface 2 faces an adjoining cell, it uses the fluid conditions in the adjoining cell to formulate its boundary condition. Likewise, if surface 2 is at the problem boundary its boundary condition may either use its cell's fluid conditions, be adiabatic, or use an input, constant sink temperature with an input, constant heat transfer coefficient. The one-dimensional heat transfer model for plates with holes assumes that half of the hole surface area is part of the upper surface and the other 46

SURTACE 2 ,.........................y........................, h l l l l l l l l l d' d l SURTACE 2 ---+ c SURTACE 1'S 9---- SURFACE 2 l E 2 l l 6 U / ///// l .......................... h......................... w_.- W ESH CE LL SURTACE 2 R Figure 4.2 Structure Arrangement in a Cell ' s 47

' halt is part' of-the' lower surface. The ' thickness of the plate for heat transfer is the plate actual volume divided by half of its' total: surface area. If activated, radiation heat transfer is treated on both surfaces except when surface 2 is specified to be adiabatic:in which case only surface 1 radiates. If a debris bed rests on a plate, then ' surf ace 1 uses the debris conditions to formulate its boundary conditions. f4.3 3 Column Type Structures -A column is a vertical cylindrical structure that may or.may not be hollow .however, it must always have a finite inside radius. Its outer surface is surface 1 and it. sees the fluid conditions. in its computational cell. The inner surface is surface 2 and it may be adiabatic (for example a solid cylinder) or it may see the same fluid conditions -as surface 1. In the current. version, there are no radiation heat transfer boundary conditions on columns. 4.3.4. General Type Structures i A general type structure is a slab with no particular geometric orientation or location in a computational cell. The structure is defined by a heat transfer thickness and a total surface area. The' -surface area is divided into equal area surfaces 1 and 2. Surface 1 always sees the cells fluid conditions and surface 2 may either see the cells fluid conditions, be adiabatic, or see an

input, fixed boundary temperature through an input, constant heat transfer coefficient.

General structures do not have radiation heat transfer boundary conditions. General structures are useful for treating the heat capacitance of structures that have complicated geometries which are not explicitly modeled in MELPROG. They are also useful for treating structures for which only the mass and surface area are known. 4.4 Mesh Arrangemen Any structure may be made up of one or more materials. Each material corresponds to a " region". The number of regions may vary dynamically as melting and crust formation occur. Each region is divided into one or more heat transfer nodes. l Initially, all nodes in a region are of equal mass; however. l during melting or crust formation node masses may change. L Temperatures are calculated at the surface of each region and at the mass center of each node. 48

4.5. Malggial_ Properties Temperature dependent thermal conductivity and heat capacity are treated for each structural material. Melting point or melting range and heat of fusion data are also treated. Because

mass, instead of

volume, is treated in the heat transfer formulation, temperature dependent densities are not needed.

4.6. Conduction Heat Transfer Solution Technique _ 4.6.1. Formulation for Material Nodes A modified Crank-Nicolson finite difference scheme [4.13 is used for both the one and two-dimensional heat transfer formulations for material nodes. This formulation is: I

• l N

- E_) y ( T"1_ -Ty) (1 - 8) 1 [ knIJ [j()nIJ (T"J - 7"l*l))

• l (sC

= ax at j-1 N B[ [kyj( )[3(T]-Ty)) + j=1 +myQy (4.1) where I temperature in node i at time n+1 (K), l Ty*l = l temperature in node i at time n (K), Ty = the specific heat of node i at time n (J/Kg), Cy = p the mass of node i at time n (Kg), my = the time step (s), at = the forward weighting factor B = 05 8 5 0.49999, the total number of temperature nodes, l N = i 49 i

kyj' the thermal conductivity between = nodes 1.and j at time n (W/m K), A (OX) the heat transfer geometry factor = between nodes i and j at time n (m), and Qy the power per unit mass in node i = at time n (W/kg). The forward weighting factor (B),- 'obtained from a linear stability analysis, is calculated as follows: ) 8 = minimum (8 82* 1, '8 .B N 3, = g--- p ) y / A t (mC where S i [ kyj (7 )[j ja1 4.6.2. Surface Nodes Surface nodes or temperatures are treated somewhat differ-ontly because they have no heat capacitance. The formulation for these nodes is: N ^ I [ k?3(7-)?3(T " -Tf*'))+Ao? (4.2) o= j 1 j=1 where the terms are the same except for Q[=thesurfacehgatflux,usuallyradiationheat transfer (W/m ), and 2 A g = the surf ace area of node 1 (m ), At the boundaries, only the surface node equation is involved. Radiation heat transfer is treated as a surface source. When heat transfer to the fluid components (steam and/or hydrogen, water, and corium in MODO) is involved the surface node equation becomes: Ay(Tyd -Tyd) =Hy'l Ay(Tj'I -Ty*I) + H 0 g is e 50

Ay(Ty'l -Ty*l). ((ky)( )y)(T)*1 -Ty*l)) + + Hcs j=1 (4.3) l AQy + 3 where the' terms are ara defined above except as follows e the heat transfer coefficient from steam H "'*l 2 = E and/or hydrogen to structure at time n*t (W/m K), the temperature of gas at time n+1 (K), Tj*l = form water to the heat transfer coefficient H n+1 = 8 structure at time n+1 (W/m K). the temperature of water at time n+1 (K), Ty*l = from corium to H "*l the heat transfer coefficient structure at time n+1 (W/m K), and the temperature of corium at time n+1 (K). Tn+1 = Currently heat transfer between corium and structures is not included except f or debris resting on structures and by radiation - from fluidized corium or radiation from a debris bed. The heat transfer coefficients described above are supplied by the FLUIDS module. For debris bed to structure heat transfer, the heat transfer coefficient is supplied by the DEBRIS module. 4.6.3 Matrix Solution For one-dimensional heat transfer the equations can be arranged-to form a tridiagonal matrix which is easily solved. A simple one material node example of the equations in tridiagonal i form ist 12j'1 7y'l T --R C C g 33 233*' a (* a) c,,Ty+1 c,,Tg+' 7 c 2 33g+1 c,Tg+1 R T c 3 3 + His + H,) 12( )12 *A1(H where C =k e gs gg 1 1 51 L_-_- _ _

l' 1 12(IV)12 ~ C --k 12 isy*l csy*) gsj'I A Qj T + T +H T +H A3(H R 3 g (I-8)C12 C 21 (mCp)2 C -C C22 " ~~{g pg 23 23(Zi)23 C23 - (1-6) k -Ty)] R -6 [-C12(Ty-T2) -C23(T3 2 (mCp)2 Tn + mpQ2 * - at 2 C -C dl~8) 32 23 +His + He3) C33 - C32

• A3(Hgs esy*l)

Ty*l gsy'l +AQ' T +H T +H and R3 - A (H 33 y3 3 Extension af this equation to many materials and multi-node cases is straightforward. 4.7. Melting and Ablation Melting and ablation are handled by an iteration on-the effective specific heat. The heat of fusion is assumed to act like a large step change in the specific heat over a temperature range of one or more degrees Kelvin. The iteration procedure consists of calcul'ating the temperature distribution; checking for temperatures that are in or are passing through the melting range; correcting these temperatures to account for the heat of fusion; modifying the specific heat for these nodes to account for the heat of fusion; and then checking for convergence of the these temperatures. If the temperatures are not converged and the iteration limit has not been reached the procedure is continued by recalculating the temperatures with the modified specific heat. This procedure is shown in Figure 4.3. 52

CALCULATE TEMPERATURES u 4 MELTING OR FREEZING yes u ADJUST TEMPERATURES ADJUST SPECIFIC HEAT no u CHECK CHECK ITERATION E TEMPERATURE UMIT CONVERGENCE yes yes U PRINT MESSAGE Av Tigure 4.3 Procedure for melting and ablation 53 L_______

4.B. Stguetugal Mechanics Modeling The objective of the structural mechanics models is to predict failure of major in-vessel structures due to sechanical loading. The specific structures to be considered fall generally into three categories: cylindrical walls, flat perforated plates supported at regular intervals, or elliptic plates supported at the edges, with or without perforations. There are also several special structures that need to be considered, such as support columns (Westinghouse PWR) or support girders (Babcox & Wilcox l PWR), hesdbolts, etc., whose presence depends on the particular l reactor design. The current structures mechanics module contains simplified models describing the three general categories mentioned above, plus a support column model, which is specific to the Westinghouse PWR design. 4.9. M_o_g l i n_g D e s i gn, C o n s i d e_r a t i o n s_ There are several factors present in reactor designs that I guided the selection of the mechanical failure models. The most l important factor is that in-vessel structures, while supporting great weight, are actually under very little stress. Thi s i s d ue to the fact that reactors are designed to minimize deformation of reactor support structures under normal operating conditions. The exception to being under low stress is the vessel itself, j which is under high internal pressure. The low stress means that j structures are at high temperatures at failure under most i conditions, indicating that a plastic or creep failure criterion l is applicable. The exception to failure at high t empe rat ur e is l I the possibility of transitr_t high pressure loading from a steam f spike or vapor explosion. This could cause a ductile-fracture type of failure, the most probable location for such a failure l being either in the lower plenum or at the headbolts[4.23 I l Another factor in the choice of mechanical models was the desire to avoid detailed calculation of stress-strain distri-butions in the structures. Calculation of inelastic mechanical behavior can be time-consuming in terms of computer resources, usually requiring a finite element or finite difference solution of the mechanics equations. The route taken for MELPROG-PWA/ MODO was to use stresses calculated at the known peak stress locations i by a n.a l y t i c formulae. This can be done by noting that very l little deformation will take place under the given low loading conditions until the structure is fully plastic at the peak stress locations. These locations can then be analyzed using the simple stress distributions present under fully plastic conditions, rather than calculating the detailed elastic plastic distribution. Before the stress level reaches the yield stress, the stress distribution is known from the elastic s ol ut i o n ; however, under elastic conditions, the structure will not deform significantly. 54 l i

l The models chosen predict only the stress and not the I strain; this requires choosing a failure criterion that does not depend on strain, but that considers both plastic behavior and high-temperature creep. A well-known criterion satisfying this requirement is the Larson-Miller Parameter (LMP) corre-lation[4 33, which relates time-to-failure at a given stress and temperature. The LMP correlation also has the advantage of allowing easy extrapolation of rupture times at a given stress and temperature to other temperatures, which was felt to be desirable due to the limited database available. The LMP correlation was chosen as the main failure criterion for high temperature, low load failure in MELPROG/ MODO. The life fraction rule is used in conj unc tion with the LMP correlation to account for changing rupture times as load and temperature conditions' change with time. The LMP correlation in SI units and associated life fraction rule are LMP = T(log go(tR) - 16.44), (4.5a) e [t dt (4.5b) - 1 J tR O where LMP = the Larson-Miller Parameter, T' = the temperature (K), t = the total time to failure (sec), and t = the time to rupture (sec). g It is intended that MELPROG/ MOD 1 will include other failure criteria as well as additional mechanics models. The primary additional failure criterion will be a maximum strain criterion this will also require upgrading the mechanics models to allow reasonable strain calculations to be performed. A ductile-fracture failure criterion will also be included, principally for use by the headbolt mode 1[4.33. l l l 4.10. Mechanics Models The following section describes the models used for the ) different types of structures in the Westinghouse PWR, including some examples of typical stresses. First, it is helpful to show a realistic example of the structures to be modeled and their j relationship to each other. A stylized representation of the major Westinghouse in-vessel structures is shown in Figure 4.4 55

i i UPPER VESSEL HEAD i l ,,n . U??EP. CRID. PLATE e VESSEL WALL CORE CORE BARREL 51 l l l l I 111111 l l I l 11 l l l I l l l l l l l l I l I LOWER GRID PLATE CORE SUPPORT PLATE l I LOWER VESSEL HEAD Figure 4.4 Westinghouse PWR Schematic. ~ l 56

l cora barrel is cttached at its upper edge to the head flange. The vessel bottom is supported at its odge by the vessel wall. The wessel wall itself is supported at the head flange by external support structures, which do not concern us here. The top core grid plate hangs on support columns from an upper core support

plate, which.is supported by its rim from the head

' flange. 4.10.1. Lower Core Grid Plate The grid plate is supported at regular intervals by support columns in a square pattern. It is perferated by holes, also in c square pattern (Figure 4.5). The presence of the holes can be accounted for by modifying the mechanical properties of the plate and then using analysis methods for a solid plate with these "oquivalent" mechanical properties [4.43. In the context of plastic analysis, the major effect of the holes is to reduce the load-bearing cross-sections. The plate can be treated by analyzing one cell, consisting or f our columns and the associated oquare section of plate [4.53. The procedure used is to find the soments at the center and corners of this cell, using elastic solutions for a flat plate supported at the corners by columns with the mechanical properties a dj us ted for the presence of the holes. The moments are compared to the maximum plastic moment M t when a moment is less than M the average stress is taken as o p, ohe half the peak stress: when a moment equals or exceeds M the overage stress is taken as equal to the yield stress at the average plate temperature. This average stress is input to the 1 LMP criterion. Since the average temperature is used in the l failure correlation, temperature distributions through the plate are not accounted for in the present model. This means that a plate may have a fairly low average temperature, but a temperature distribution ranging from molten on one side to near water saturation temperature on the other side. However, large deformations will not occur until plastic or high creep rate conditions exist in the entire plate thickness. Therefore, this lack of information is not considered a seriout drawback at the present level of detail in the mechanics models. The actual oquations used are as follows: 2 M = Sqb /c (4.6a) e M,- -q in( ) - a)(1+v,)/c (4.6b) G 4 W t 57 L____-______

] n i ) i e a i ) 1 ,/ ,f 1 i

SUPPORT

,! COLUMN ~........

• 1 I

l HOLE GRID PLATE a SECTION i ( I i 1 ~. ,/ / 's. e l l ~. ' Figure 4.5 Cell of Westinghouse Lower Core Grid Plate Top View. 4 ' 'l I 58

'6M (4.7a) o a --2 h d-a (4.7b) c = where the elastic moment at center of plate (N), M = the. elastic moment at edge of plate (N), c M,' - the. length of plate edge (column pitch) (m), b

== the diameter of support column (m), c the load per unit area of e'. ate (Pa), q = = the effective-Poissons rat. v, a = a constant = 0.811 for squar., late cell, B' = a constant = 0.0331 for square plate cell, the ligament efficiency, c = d = the hole pitch (m), M .= either soment (N), c = the peak stress at either location (Pa), j the plate thickness (m), and h = the hole diameter (m). a = The load q in the above equations is related to the actual load Q by the formula i 2 d (4.8) 4 Q d2 2 3 L where Q is the load per unit cell area. The maximum plastic moment is determined by setting the maximum stress in the usual clastic equation equal to the yield stress and multiplying by the ratio between the maximum plastic moment and the moment at incipient yield (1.5): 2 1.5cy h (4.9) M = 7 p where M = the plastic moment of the plate, and p oy = the yield stress. 59

( 4.10.2. Elliptic Plate The elliptic plate model describes an elliptic perforated shell, with allowance made for the possible presence of a central L hole, as in the Westinghouse lower core support plate. This model is be used for structures such as the vessel top,. vessel bottom, or core support plate. The equations used are those for an elliptic shel1[4.3] with the presence of perforations taken into account by modifying the mechanical properties. The central hole, if present, is modeled by applying a stress concentration factor. The model is used-to evaluate the longitudinal and azimuthal stresses at the outside of eech radial ring containing the plate. The shear stress at the crown is also evaluated if there is no central hole. The maximum of the edge and crown shear stresses and the local yield stress is used as the input stress to the LMP criterion. The applied load is the average plate load corrected for the holes, as in Equation 4.9 above. If no holes are present, then the pressure differential across the two surfaces is also applied as a load. The elliptic shell equations are P 2 og = Q'g, (4.10a) o = Q' h (1 ), (4.10b) g 2rg Q' =P+ q, (4.10c) r2 = ((A2-R2)[ )2, p2]), (4,33,) 2 r3 =r (4.11b) g. l I O E o i 60 1

j. I l I where I o g. the longitudinal stress at outer edge of radial ring, I = the azimuthal stress at same location. e6 = A = the radius of plate major axis (radial) - l B = the radius of plate minor axis (vertical), R = the outer radius of radial ring, o P = the pressure, and the total load. Q' = ~ The shear stress at R is then given by o og.- 0 0 (

• 1

} to= 2 where t is the shear stress at the edge of the radial ring. o The equations for the crown (with no central hole) are i =f. (4.13a) o p og - o p (4.13b) t = o 2 i where o = the average radial stress, p t = the shear stress at crown of shell. o If a central hole is present, then there is no crown and the edge shear stress for the edge of the innermost radial ring is multiplied by a stress concentration factor given as (.g-)2 C K 1 + (4.14) o where o K = the stress concentration factor, C = the radius of central hole. 61

The above equations can also be used for a flat plate supported at the edges by letting the minor axis B approach zero. Use of the shell equations for thick plates such as the core support plate is somewhat dubious, since it neglects bending moments. However, since its application is at high-temperature, plastie q conditions, this may not be too serious a limitation. The level l of detail is commensurate with the level of the other mechanic _ i models. l 4.10.3 cylinder The thick wall equations for a cylinder are used to evaluate the stresses on the inner surface (the peak stress location for the hoop stress)[4.33. The average axial stress is also determined. These stresses are then used to form the Tresca shear stress, which is input to the LMP criterion. The thick wall equations are 2 = qRf (1 - (R ) ), (4.15a) o p 2 (R ) ), (4.15b) e = qR (1 + o - Rf, (4.15e) 2 2 h =R where the radial stress at inner wall, o = p = the azimuthal stress at inner wall, 0 0 q = the radial load, R = the radius of outer wall, o R = the radius of inner wall. g The Tresca shear stress t is formed from the maximum and minimum of the principal stresses: "I", (4.16a) t = z), (4.16b) (o 0 O o = max p, 0' max a 62

z). (4.16c) min "'"I"'(O r' 'O' O o l I 4.10.4. Columns ~ The~ column-model is currently set up for cylindrical columns ~ loaded in compression. specifically the core support columns in a ~ Westinghouse PWR. An extended Euler buckling model.is used[4.6): 2 5 EI g '(#*17*) i P = ep 2 b s(ch-cf) (4.17b) 'I = 4 l where P = the critical. buckling load, er E = the tangent modulus, f g L. = the length of column, .? I - the moment of inertia of the column cross-section, = the outer column radius, j c o = the inner. column radius. l eg 4.10.5. Example Using data from the Zion FSAR[4.7], a cell for the lower grid plate is loaded by four fuel assemblies, resulting in a weight on the cell of 2600 kg. The approximate dimensions of the cell are as follows: h = 0.051 m, a = 0.074 m, b = 0 304 m, 0.086 m, c = 0.152 m. d = These numbers result in a ligament efficiency c of 0.4, a load Q of 0.275 MPa, and an area ratio (cell area / actual area) of 1.23 The actual load per unit area q is then equal to 1.230, or 1 0.34 MPa. The effective Poissons ratio for this case is 0.33, J roughly the same as that of the actual material, 304ss, at low j 63

temperatures. At high temperatures near melting, it is 0.44. l [ Using this high temperature value, we find M

• 2000 N*

c M, = -4100 N, M = 0.0006cy. p To put the above results in perspective, the yield stress ey of PWg) N. is 130 MPa, 304SS at 550 K (operating temperature for a resulting in a maximum plastic moment of 7.6x10 The grid plate is obviously under very little stress and will fail only at ' temperatures very close to the melting point. The yield stress at 1273 K has dropped to 20 MPa.

However, it must drop to 4.8 MPa before the center of'the plate cell is totally plastic.

Based on a linear extrapolation from 20 MPa at 1273 K to O MPa at 1700 K, this stress corresponds to an average temperature about 100 K less than the melting point. Therefore, grid plate failure is most likely to occur only as its melting point is approached. 4.11. References 4.1. B. Carnahan, H. A.

Luther, J.

O. Wilkes, Applied Numerical Methods,, John Wiley and Sons, New York (1969). 4.2. D. V. Swenson and M. L. Corradini, Monte Carlo Analysis of ~~ LWR Steam Explosions, NUREG/CR-2307, SAND 8I-1092, Sandia National Laboratories. Albuquerque, NM (October 1981). ~ 4.3 J. F. Harvey, Pressure Component Construction, Van Nostrand Reinhold (1980). 4.4 W. J. O'Donnell and B. F. Langer, " Design of Perforated Plates", Trans. of_ the ASME,.T. of Eng. for Industry, p307 (August 1962). 4.5. S. Timoshenko and S. Weinowsky-Krieger, Theory ~of Platos ~ and Shells, McGraw-Hill (1959). 4.6. T. H. Lin, Theory of Inelastic Structures, John Wiley and Sons, Inc., New York (1968). 4.7. Zion Nuclear Station, Finel Sa f e t'; Analysis Report, Commonwealth Edison Co. (1971). O I I l 64

l l 1 l-5 DEBRIS Module 5.1. Introduction This section describes the physical models and numerical methods in the core debris analysis module. DEBRIS. This module is designed to analyze the formation of debris beds and the subsequent behavior of the debris. Since the debris has a strong influence on the integrity of core structural components and on the behavior of fission products, it is essential that the behavior and state of the debris be analyzed. The term " core debris" is meant to represent all forms of failed fuel

rods, control

rods, and structural components.'

During the course of an accident, the debris can take on many forms. For example, the debris can be either a particulate or cohesive and it may be stationary or moving. Each of these different states may require special analysis. In order to facilitate analysis of the debris behavior, three regimes for the state of the debris are considered in the modeling. The first is the stationary, pre-dryout regime in which the debris is not moving and its temperature is less than the saturation temperature. This case occurs, for

example, when debris accumulates on grid spacers or freezes and jams between intact fuel rods and liould water is still present.

For _ this regime, the debris can take the form of a particulate debris bed or a cohesive agglomeration of debris. Either type of debris form is referred to as a debris bed in this report. The second regime is the stationary, post-dryout regime. This regime is similar to the first regime in that the debris is stationary.

However, in this case there is no liquid water present to cool the debris and it is assumed that the debris has dried-out.

Consequently, the debris can overheat and melt. The melting can, in turn, lead to compaction of the debris and a rapid attack on the supporting structures. The third regime is the non-stationary regime which covers all cases in which the debris is moving. This case occurs, for l

example, when eithe,r molten or solid debris is moving down j

through the core or'when particulate debris is swept out of the core. In this regime, the debris is not considered to form a debris bed. The DEBRIS module is comprised of two main parts. Each part considers a different aspect of the debris behavior. The two parts are 1. the Debris Bed Formation and Release Model, 65

2. the Post-Dryout Debris Bed Model. The debris bed formation and release model predicts how debris beds will form and whether or not the debris is stationary. The post-dryout debris bed model analyzes the stationary, post-dryout regime. This model is specifically conceried with calculating debris heat up, melting, relocation, and crusting, as well as crust remelt, within the debris. If tne debris can be classified as being in either the non-stationary or pre-dryout

regimes, then, as discussed in Section 2, the FLUIDS module calculates the behavior of the debris.

That is, FLUIDS analyzes these regimes by solving mass, energy, and momentum equations for the debris. This analysis is not detailed in that FLUIDS uses relatively large control volumes i in its analysis, l The bulk of the DEBRIS module deals with the post-dryout modeling. There are three phases of the post-dryout behavior of debris beds that the model is designed to calculate. First, the model calculates the transient heat-up of dry debris beds. This analysis considers the heat sources in the debris due to decay heat and oxidation as well as conduction, convection, and I radiation heat transfer in the debris bed. From this analysis the temperature distribution as a function of time can be calculated for the debris bed. Second, the model calculates the melting of debris and relocation of the molten debris within the debris bed. This analysis uses a simplified ternary phase diagram to predict the formation of molten debris. The relocation step uses simplified laws for porous flow to predict the behavior of the molten material. From this analysis, the mass and compositional distribution of the debris is calculated as a function of time.

Third, the model provides boundary conditions for the STRUCTURES heat transfer model and input to the RADIATION module.

The details of the DEBRIS stodule are discussed in this section. The discussion describes the current modeling and indicates the limitations of the modeling. 5.2. Overview of the Module The DEBRIS module has as its responsibility the calculation l of the behavior of the debris. The debris can be modeled in detail or an average analysis (performed in FLUIDS) can be used. The decision as to which type of modeling to use is based on the debris regime as discussed above. In order to better understand the treatment of debris and debris beds, it is useful to trace through the calculational sequence in the DEBRIS module. l 66

If there is no debris present in the problem, then this -Dodule is not active. However, once debris is f ormed, the. DEBRIS module is activated and the debris bed-f ormation and release model is. entered. This model is responsible'for 1. determining'the location of debris beds, l'- 2. defining the state of the debris beds, and 3 determining if the debris beds should release. The model examines.the FLUIDS data arrays to determine. the location of debris throughout the vessel. This information is used.to determine whether or not a detailed calculation should be q performed. The debris must be dry and-stationary before a l detailed calculation is initiated. If it is determined.that a q dry debris bed has. formed, then this model sets a flag to begin a detailed calculation and determines the initial and boundary conditions for the bed. Otherwise, no flag is set and the -analysis of the debris its perf orred by FLUIDS, If the formation model determines that a detailed calculation should not be. perf ormed, then the analysis of the debris behavior is only performed in the FLUIDS module. The analysis in FLUIDS is not detailed in the sense that the geometrical scale is large. The justification for using,this approach is two-fold. First, if the debris is not dry, then its temperature is relatively low and its behavior i s. benign. Second, if the debris is moving, it is probably moving fast and will eventually come to rest. The behavior during this short time scale does not need to be analyzed in detail. Hence, the calculation performed in FLUIDS is sufficient. If a dry debris bed is formed, then the post-dryout debris bed model calculates the debris behavior. This model combines a mixture energy equation for the debris bed with equations for seit relocation. One-dimensional geometry (axial) is used in the model since scoping calculations indicated that radial heat flow in dry debris would not be significant. The energy equation is time dependent, and includes internal heat generation (both chemical oxidation and decay heat), conduction, convection and radiation heat transfer. In the case that a detailed calculation is to be performed, the following steps would be taken. First, the bed's geometry., boundary conditions and initial conditions are determined and then the actual calculation begins. The basic equation solved at each time step is a mixture energy equation. This equation is written in terms of an equilibrium temperature and, by solving this equation, the temperature-distribution within the bed can be calculated. In this model, it is assumed that the gas and solid are in equilibrium. This assumption is reasonable in view of the large surface area of the debris and tightly packed nature of the debris bed which result in high flow resistance and concomitant 67

l i 1 low velocities. It is also implicitly assumed that the radiation heat transfer within the bed can be accounted for by using an effective conductivity model. Based on several in pile celt progression experiments, this assumption is quite rea-sonable[5.13. After the energy equation is solved, the melt formation and relocation part of the model is entered. Using simple phase diagrams to describe the various phases which may be present, the state of the debris (i.e. its composition and temperature) and the amount of liquid (if any) within a bed can be calculated. If liquid exists, it is relocated using simplified porous media equations. Following relocation, the melt within the bed is equilibrated with the solid which may cause freezing of the melt. The porosity distribution along with the composition distribution are also updated by again referencing the phase diagram. In the current

version, debris components which are accounted for Zr(0)-UO stainless

steel, and the include UOp,-

Zr,

Zr02, 2,

Ag-In-Cd alloy control rod material. After the relocation

step, material properties for the debris are recalculated.

This marks the end of a time step. At this point new temperatures and porosities have been calculated for the bed. This final result is then used as the initial condition for the next time step. During each tize step of the detailed calculation, the debris bed formation and release model performs three additional functions. First, it determines whether debris is being added to the top of the bed or added internally (e.g. from rod break-up) and updates the geometry accordingly.

Second, the boundary conditions for the bed are updated based on FLUIDS RADIATION and STRUCTURES information.

Finally, the debris release part of the model determines if the bed supporting structure has failed. If this occurs, the detailed calculation ceases and the bed is allowed to relocate. The above discussion has simply traced through the calculational sequence. The specific details of the models are described in the following sections. 5.3 Debris Bed Formation and Release Model The debris bed formation and release model is concerned with identifying the locations of debris

beds, determining their initial and boundary conditions, and determining when they.should release.

In the model, the area of greatest uncertainty lies in l determining how the debris progresses from its initial formation J stage to the point where it can be called a debris bed. Certainly, this is one area where current and future experimental 68

activities are crucially needed to i mpr ove the current state of knowledge. To understand the problem that is being addressed, it is necessary to review how debris is forr.ed. There are three potential sources for debris: fuel rod

failure, control rod

failure, and structural component failure.

However, in MELPR00-PWR/ MODO the debris formed by structural component failure is not accounted for.

The fuel rod and control red failure mechanisms are the same as in the original MIMAS code [5.23. The control rods are assumed to fail when the gas temperature reaches 1700 K which represents the point at which their steel cladding will fail and release the molten internal absorber material. For the fuel rods, there are four criteria for failure. These criteria are 1. shattering upon quench, 2. fuel-cladding contact if the cladding is brittle, 3 partial cladding melting if oxide layer is thin, and 4 complete melting of the cladding. Further detalla of these limits are given in Section 3 Once failure occurs, the mass and energy of the fuel rods or control rods an well as their effectise diameter (or surface-to-volume ratio) are supplied to the FLUIDS module. FLUIDS then begins to transport the debris and to calculate the heat transfer between the debris and the coolant. In general, the debris will cool and move downward, usually through intact rod regions. At each time step, the debris formation model examines the FLUIDS data arrays to determine if a debris bed has formed. There are three criteria which can lead to a bed being formed. The first criterion is the most obvious and states that if the debris accumulates on an impermeable plate, then a debris bed will begin to form. The second criterion is similar and states that if the cell below a given cell is full of debris, then the i upper cell can also be considered to be a debris bed cell. The I third criterion is based on the local debris fraction exceeding a specified limit. Currently, debris beds are allowed to form in any cell in which the local debris fraction exceeds limits as defined by Moore [5.3]. These limits are based on experimental and numerical studies of debris packing fractions for THI-type geometry. Typically, the packing fraction is approximately 0.44 However, if the debris is entirely molten, a higher limit of 0.9 is used. The application of the first and second criteria are straightforward, but the application of the third criterion requires the calculation of the local debris fraction in each cell (i.e. the two dimensional distribution). The local debris fractions are calculated in the DEBRIS module by solving l 69 l I l

l auxiliary continuity equations. In solving these equations, it is assumed that the debris velocity is radially uniform and that there is no radial motion of the debris. The formation model uses these three criteria to identify 1 the locations of debris beds in the vessel. If a cell contains a bed, then the debris velocity at the bottom of that cell is set to Zero. Once a bed has formed in a bell, the model examines the conditions of the debris to determine whether or not it is dry and static. For pre-dryout or non-stationary conditions, the heat transfer and momentum exchange between the debris and the liquid and vapor are calculated exclusively by FLUIDS. If the debris is dry and static, then the bed formation model obtains the initial temperature, porosity, and composition of the bed I l based on values obtained from FLUIDS. After these conditions have been obtained, the model nodalizes the bed according to user specification. The boundary conditions for the bed are obtained by searching STRUCTURES and RADIATION data arrays to determine j either an appropriate boundary temperature or surface heat flux. l The vapor flow rate through the bed is obtained from FLUIDS. The initial conditions for a bed are only obtained 'once. However, on subsequent time steps, the boundary conditions and vapor flow rate through the bed may vary. Consequently, these quantities are continually updated. Changes in the geometry of the debris bed are also accounted for. Specifically, if debris is added to the top of the bed, the ) model allows for the bed height to grow with time by adding nodes to the bed as necessary. Also, if beds form among intact rods and the rods subsequently fail, the debris is added internally to the bed by adjusting the composition, flow area, and porosity of the bed. 1 In addition to forming beds and updating their geometry, the l model also determines when beds should release. The word " release" means that the beds are no longer supported and, hence, i can begin to move. Obviously, if the structure which is supporting a debris bed f ails, then it is possible for the bed to begin to move. This represents a transition from the stationary l regime to the non-stationary regime so that detailed calculations 0f course, when the debris again cease when the bed is released. 2 comes to rest, it may form a new bed.

However, the detailed l

information about the bed behavior is not retained during this transition. The release of a bed is controlled by examining the STRUCTURES and PINS data arrays on-each time step to determine if the bed supporting structur es have failed. STRUCTURES provides 70

i information concerning the-failure of support plates which would be one mechanism f or-bed release.. If a bed is resting on.a Jetructure which fails, then the bed will release. Similarly, PINS provides.the. rod failure information. 'If a bed is formed by i freezing between rods and the rods subsequently fail, then the ] bed will release. .j l ~5.4. Post-Dryout-Debris Bed Model 4 The post-dryout debris bed model ' analyzes the-behavior of debris beds after they become dry. This model. is specifically concerned with calculating debris heat up, melting, relocation, Ond' crusting, as well-as crust remelt, within the debris'.and also the heat transfer to supporting structures. To : accomplish this task, the model divides the analysis into two sections. First. the model calculates the transient heat-up of. the dry debris. This - analysis considers the heat sources in the debris due to decay heat and oxidation as well as conduction, convection and radiation heat transfer in the debris. From this analysis, the a -temperature distribution as a function of time can be calculated j for the debris. Second,tne model calculates the melting of the ] dobris and relocation of the molten debris within the debris bed. l This analysis'uses phase' diagrams for the various components to predict the formation of-molten debris. The relocation step involves using simplified laws for porous flow to predict the

behavior of the molten material.-

From this analysis, the mass and compositional distribution of. the debris is calculated as a function of time. By knowing the state of the debris..the heat transfer.to supporting structures can be calculated. This section discusses the details of the post-dryout debris bod model. Both.the heat-up part as well as the seit relocation part are discussed. Additional modeling aspects of the problem j cro also described. J 5.4.1.- Debris Heat-up Model As discussed above a mixture energy equation is solved in tho model. to determine the temperature distribution in the debris. This equation is obtained by ' combining the energy equation f or the debris with that of the gas and is given as (((1 c)pCp), + (cpCp)g) BT/Bt - k,ff 92 T - cGC BT/Bz + Q (5.1) l p l: where (pC3)s is the effective volumetric specific heat of the colid (and/or liquid) and is given by j l .) 71 L________-_____

(5.2) (pCp)s - 1(opCp)g where the summation is over all components of the bed. In this ~ equation (and in this section), the symbols used have the 3 following' meanings: c is porosity, p is density (kg/m ), C is S (s),2),ff specific heat (J/kg-K), T is temperature (K), t is time is the effective conductivity (W/m-K) G is mass flux (k s-m 2 is the axial dimension (m), Q is volumetric power (W/m ), and a is the component volume fraction. The term on the left hand side represents the rate of change of the mixture energy (both solid and gas.). The first term on the right hand side represents the energy gained by conduction. An effective conductivity is used which allows for both conduction and radiation heat transfer effects to be included. The next term on the right hand side represents the energy lost by a convective flow of gas. The last term on the right hand side represents the energy gained by internal heat generation. This equation implicitly assumes that the gas and solid are in equilibrium and that an effective conductivity can be used to represent both conduction and radiation heat transfer within a debris bed. However, since porosity variations are considered in the

model, this model would apply equally well to either particulate or cohesive debris.

The model also assumes that the vapor flow rate is uniform through the bed. A standard finite difference approximation for the differential equation is used. This approximation is obtained by first integrating the equation over the volume of a node. The temperature of a node is defined at the node center and represents the space averaged temperature of the node. In the following equations the center of a node has axial index j. Integrating over the volume of a node and rearranging one l obtains AZ ( (((1-c)pC )s + (cpCp)g)BT/Bt - Q) p 3 [T]lJ+i BT / Bz)l J +i =0 (5.3) A [cGC - A [keff + p l 2 where A is the area (m ) and 2 is given by 72

(5.4) Z= zj.g - zj.g This equation represents the volume average energy equation for a node in the problem. S t an d a r.d differencing approximations are then used to arrive at the finite difference form of the energy equation. The finite difference equation solved in the code is given by (cpCp)g)(T-TU)/At - Q] AZ[(((1-c)pCp), + - A[ (k/az)).g(T) 3-T)) (k/Az)).g(T)-T.g)) j 0 (5.5) + AcGCp(T).g - T).g) = where k)k) 3(Az) + Azj.3) (5.6) (k/Az)) g = j z) ) + k) 3 Az) A k (5.7) Az) = zj.g - zj.g and (5.8) Azj.g = zj.3 - zj The temperature, T. is treated in a fully implicit sense in that all temperatures except Tn are evaluated at the new time step (n+1). This equation closely resembles the transient heat conduction equation. The main difference is the addition of the term representing forced convection of the gas. The set of equations obtained by writing an equation for each node can be formulated into a tridiagonal set of equations which can be solved using a straightforward Caussian elimination procedure. \\ l. 73 L-_____-______-__

Host of the material prope? ties are allowed to be temperature dependent and must be continually updated during a calculation. Additionally, the conductivity of the debris will also depend on the porosity and volume fraction of the constituent materials. Currently, two debris conductivity models are included in the code. Both of these have been described in Reference 5.1. The first is the Imura-Takegoshi/Vortmeyer mode 1[5.4,5.53 which has been demonstrated to represent an upper bound to the conductivity. The second model is the W111hite-Kun11-Smith /Luikov mode 1[5.6,5.73. This model yields a lower bound to the conductivity. Both models include the effects of conduction and radiation within the bed. The multi-component nature of the debris is accounted for through the use of the Maxwell mode 1[5.73. 5.4.2. Helt Relocation Model The melt relocation model is designed to calculate the redistribution of material wi thin' the debris bed after melting begins. The model combines simple porous media flow equations with mass and energy conservation equations. The treatment of i the relocation is relatively simple, but is based on physical j arguments and experimental observation. This section describea l the details of the modeling. l The physical picture which the model is attempting to simulate can be described as follows. As the debris heats, it may eventually reach a temperature ht which liquid begins to form. As the molten debris

forms, it begins to experience various forces within the bed which will govern its motion.

Gravity tends to move the melt downward, while viscous forces tend to prevent motion. If the melt Wets the remaining solid and its surface tension is high, then capillary forces tend to draw the melt into the bed in a spherically nymmetric manner. This behavior is analogous to that which one observes by touching a sponge to a pool of water. If wetting does not

occur, agglomeration of the melt may occur.

Currently, the model only considers the wetting case. In addition to these forces there are at least two other mechanisms which will limit the motion of the melt. First, it is possible for the melt to freeze in the colder regions of the bed. Due to the large surface area in the bed, freezing can be a rapid mechanism. Second, as the melt is drawn into the b e 'd it may plate out as a liquid film on the particles. Surface tension effects will allow the film to remain on the particles. The maximum amount of film deposition can be approximated from residual saturation data f or. porous media. Values of the film thfekness as a fraction ofl the bed porosity range f rom 0.75 f or a bed of 105 porosity to 0.12 for a bed of 60% porosity [5.83. I 74 l

Hence, film deposition can reduce the distance the melt can relocate.

In order to calculate the motion of the melt exactly, one would have to solve mass, energy, and momentum equations for the

liquid, as well as mass and energy equations for the solid.

These equations would need to be solved simultaneously and resulting set of equations would be quite complicated. In view of the uncertainties in the material properties and the phenomenological behavior of the melt, a highly sophisticated model is probably not warranted. Instead, simplifications to the governing equations have been investigated in order to gain insight into the physics of the melt behavior. In particular, the hydrodynamic and heat transfer effects have been studied separately. Considering only the hydrodynamic forces, the equation of motion for molten material entering a debris bed is given by obi (5.9) -pV/c - pg + pBV/Bt = 2 gsvelocity (m/s), u is viscosity (N-s/m ), oc In this equation, V is perm (m ), g is the gravitational acceleration (9.8 m/s'gability ), and o is surface tension (N/m). The term on the left hand side represents the acceleration of the melt. The first term on the right hand side is the viscous force term, the second term is the gravity force

term, and the thir-d term is the capillary force term.

This equation has been solved for the case of an infinite reservoir of UO liquid moving downward into a 2 bed (see Figure 5.1). The results show that the melt porous UO2 is rapidly accelerated due to the capillary force and then decelerated due to the viscous force. Within a few seconds the motion essentially ceases and the penetration distance is limited to a few centimeters. Gravity continues to move the melt downwcrd, but at a very slow rate. In a real case the " reservoir" of liquid is not infinite, but is relatively small. In fact, the liquid will most likely be drawn away from the melting zone as it forms. With a srall reservoir that can be depleted, the effect of the capillary force is eliminated once the melt is completely drawn inte the bed. ) This is because the driving force for capillary force is the difference between the permeability (or porosity) in the reservoir and that in the bed and if there is no reservoir there is no force (assuming uniform porosity in the bed). The resulting effect on the motion is that while the time scale for 75

l-l-l UQUID U02 RESERVOIR O "Z f U0 BED 2 i Figure 5.1. Geometry of Infinite Reservoir Case I 76 i

the motion remains

short, the penetration distance is much I

shorter relative to the infinite reservoir case. In view of the short time scale of liquid motion, the details of the melt motion are probably not important to calculate. Rather, the key parameter to calculate is the melt penetration distance. This distance can be estimated in a variety of ways. One way is to assume that the motion of the liquid ceases when it is entirely in the bed. This assumption is valid for upward flow or horizontal

flow, but is only approximately correct for downward flow.

Figure 5.2 illustrates l the case of a liquid slug moving into a porous bed. The initial 1ength of the slug is X the porosity of the reservoir is c and the porosity of the bed is cb. At any given time

tbe, o,

position of the leading edge of the slug, Xg, and the position of are given as the trailing edge of the slug, XT, t ,t (5.10) V gdt +X Xg= o ,t (5.11) Y dt XT* T When the position of the trailing edge equals the initial slug length, the motion is assumed to cease or in other words, ,tf (5.12) V dt X = o T is the final time. where tf If one further assumes that the flow rate at the leading edge equals the flow rate at the trailing edge, then one finds i (5.13) gb*YATp A V where A is the flow area and is proportional to c and the subscripts b and r refer to the bed and reservoir, respectively. This equation can be rearranged to find e 77

I 1 l o) Time = 0 l O ~ a MELT X. (c,) ~ BED (c.) b) Time = t (0 < f < t,) 0 6 Xr

• VT XL MELT U

.................}...... vt BED c) Time = t, 0 n a X3 X y t b X., BED Ir U Figure 5.2. Slug Penetration into Bed 78

VA

V C

Tr (5.14) Tp_ Vg-A C = b b-The position of the leading edge at the time when the trailing edge enters the bed is given by ,t f yCor+X, (5.15) Y dt =- XL, max ", L Hence, the maximum penetration distance is given as XCor M.16) max " X, max - X = X L n c Therefore, by simply knowing the porosities of the melting and non-melting regions, the maximum melt penetration distance can be estimated. For downward flow, gravity will continue to move the melt even after the capillary force is lost. However, as discussed ~ above, a liquid film will be deposited on the particles as the melt moves through the bed. If the bed is large enough, the melt may'be totally depleted as a film. A second criterfon for the melt penetration distance can be derived it.it is assumed that the melt is completely deposited as a film. In this case, a mass balance of the melt yields the following expression for the penetration distance XCor (5.17) X max g y I where Y is the fraction of the flow area which contains the film. To assess the validity of these two models, comparisons have been. made between the models and calculations made with PLUGM[5.9]. This computer code solves the coupled mass, energy, and momentum equations for one-dimensional flow of liquid in a i porous-medium. In this code freezing of liquid and film deposition are treated. However, axial conduction in the debris is not calculated. In these comparisons, Eq. 5.16 has been used t., l } 79

fer.the upward flow cases and Eq. 5.17 with Y = 0.2 has been used f or the downward flow cases. As seen in Figure 5 3, the agreement is quite satisfactory for upward flow. For downwarc flow'the agreement is not as good, but is still quite acceptable.

Hence, Eqs. 5.16 and 5.17 can be used to estimate the 1

hydrodynamic penetration distance of melt into a porous bed. 'l The above discussion has only considered hydrodynamic effects. In an actual bed, there will be temperature gradients ttat will allow the melt to freeze and possibly form blockarvs in the bed. The hydrodynamic limits may be viewed then to be absolute maximum limits since thermal effects will shorten the penetration distance. In view of these potential thermal

effects, a third criterion has been derived to estimate the penetration distance.

A simplified steady-state heat balance for molten debris that is freezing is given as (5.18) (1-c)(pC )3 ( T, - T g) cp h ts = g p where T, is the melt temperature, T is the local be d 3 is the heat of fusion. This equation can be temperature, and his rearranged to find the maximum bed temperature, T that will

max, freeze all the adjacent liquid cp h is g

=T,- (1-c)( Cp)s max If the bed temperature is less than T then complete freezing

max, is possible.

If the heat transfer rate is very large, then this temperature limit can be used to estimate the position where complete freezing and, hence, maximum penetration will occur. Due to the large surface area in the bed, the heat transfer is expected to be large. Hence, with the calculated temperature distribution and this temperature limit, the maximum penetration distance based on thermal effects can be calculated. i To assess this model, PLUGM predictions were compared to the model predictions for downward flow cases. A tem pe r atur e } gradient of 20 K/mm was 'as e d, and the freezing penetration distance as a function of porosity was calculated. For PLUGM, three different reservoir neights were used to assess the influence of this parameter. As seen in Figure 5.4, the simple 80

25 i i i i i i 8 LEGEND T ELT SIMPLE MODEL (UPFLOW) X, SIMPLE MODEL (DOWNFLOW) 1 o PLUOM (UPFLOW AND X, =1 mm) Of.:BRIS 20 O PLUGM (DOWNFLOW AND X, = 1 mm) WED (tm) 1 15 o I M D D 10 6 \\ \\ \\ \\ \\ \\b 5 \\\\ N~%ws i i i i I O O 10 20 30 40 50 60 70 (B N Figure 5.3 Simple Hydrodynamic Model Compared to PLUGM 81

40 LEGEND 30 - SIMPLE THERMAL MODEL o PLUGM X, = 1 mm 3 0 PLUGM X, = 5 mm A PLUGM X, = 50 mm E S 20 - M i. A 0 10 e 6 0 o o o k l l 1 1 0 0 10 20 30 40 50 60 POROSITY (%) Figure 5.21 Simple Thermal Model Compared to PLUGM 82

thermal model yields an upper bound for the penetration distance The reason j f or all porosities, even for large reservoir heights. for this is that partial freezing of the liquid combined with film deposition leads to a faster depletion of the melt and limits the penetration. In the X -1 mm case, hydrodynamic o l effects tend to dominate the penetration distance for porosities greater than 305. In this case, the hydrodynamic limit yields a better limit. The criteria developed for predicting the penetration distance have been shown to provide an upper limit to the penetration distance calculated using the fully coupled mass,

energy, and momentum equations found in PLUGM.

While not aneurately predicting the correct penetration distance in all these criteria are probably sufficient to estimate it. In

casas, the relocation model, each of these criteria is evaluated and the minimum distance is used in the calculation of the melt progression.

With these simple equations for the melt penetration distance, it is now possible to describe how the melt relocation model works. The calculational sequence in this model begins after the temperature distributions have been calculated. With the calculated temperature distribution, the amount of liquid (if any) is calculated in the phase diagram routine. If liquid is present, the melt relocation routine is entered. There are two possible paths for calculating the relocation depending on whether or not a stable crust has formed. If no crust has formed, then the size of the melting zone is determined and the mass and energy of the molten material is calculated. The upward penetration distance from the top of molten zone is then calculated and an attempt is made to move 50% of the melt this distance. This number is based on the assumption that the capillary force is symmetric. The porosity distribution is updated and if a blockage

occurs, then the excess melt is reserved for downward flow.

q f Following the upward relocation, the downward penetration I distance from the bottom of the molten zone is calculated. An attempt is then made to move the remaining melt this distance. Again porosities are updated and if a blockage occurs, the melt is backfilled on top of the blockage. The second option for relocation occurs after a stable crust has formed. In this case, the model does not allow melt to move through the crust. The model first determines the amount of melt between the crust and the top of the bed. The upward penetration distance is calculated and relt is moved upward this distance. } The melting zone below the crust is then calculated and the downward penetration distance is determined. The melt is moved .+ 83

downward this distance and the porosities are updated. If a j blockage occurs the melt is backfilled on top of the blockage. f Following the relocation

step, the temperatures of the liquid and solid in each node are equilibrated.

The equilibration is assumed to' occur instantaneously which is reasonable in view of the large surface area in the bed. Once equilibrated, the phase diagram routine is used to correctly j define the' state of the bed. ] 1 The relocation model also has a few additional features which are worth noting. First, there is an upper limit on the porosity of a node. If the porosity exceeds this limit, then the node is assumed to collapse down into the node below.

Second, mass conservation equations are solved for all components of the bed.

This means that the distribution of UO Zr0 ZP' i 2, 2' Zr(0)-UO steel and control rod material within the bed are 2, calculated. Since the melting occurs over a relatively long period of time, many stepwise movements of the molten material will occur. In an integral sense, the sum of all these steps should yield a good approximation of the actual melt motion. 5.4.3 Debris Phase Diagram The phase diagram routine in the DEBRIS module uses the temperature and composition of the debris to predict the state of the debris. The model is relatively

simple, relying on simplified ternary and binary phase diagrams.

The U-Zr-O system is f airly complicated, but the limited knowledge in this area has been combined to produce a simple but functional method of determining the state of the

debris, and in particular the melting and freezing of the debris constituents.

In the following discussion the term " melting" is used to mean liquid formation and the term " freezing" is used to mean solid formation. The current model considers melting and/or freezing of the following constituents: 1. Ag-In-Cd control rod alloy, 2. stainless steel, 3 Zr metal, 4 UO 2' 5. Zr0 2' and 6. UO -Zr02, 2 7. Zr(0)-UO p. 84 e

The first three constituents are treated as simple single phase melting and freezing processes. Although the steel and control rod material are alloyt, they are assumed to form ideal liquid system is treated as a simple and solid solutions. The UO -Zrop 2 binary eutectic which forms a liquid with a 50-50 composition. This liquid will form if UO2 and Zr0 are the only solids present 2 in a given node and the temperature is high enough. The system is treated as a binary system in which the Zr(0)-UOp composition of the liquid is temperature dependent (see Figure 3 4). Liquid is formed by first melting the Zr metal and into it. The amount of UO dissolved 2 then dissolving the UO2 will depend on the local temperature. Table 5.1 summarizes the phase rules used in the model. If, the composition and temperature of a node is given, then the amount of solid and liquid in the node can be calculated. 5.4.4 Zr oxidation Model As a debris bed forms, the amount of Zr metal in the bed can be significant. If there is a supply of steam to the bed, the oxidation of the Zr can provide a significant heat source. In fact, if the cladding is rubblized, it will have a high surface to volume ratio and this will lesd to very large oxidation heat source. The oxidation model examines the composition of the bed, the amount of steam and hydrogen flowing through the bed, and the temperature of the bed to determine the oxidation rate. The Zr particulate is assumed to be represented by an effective spherical diameter, Once the rate of oxidation is determined, it is multiplied by the surface to volume ratio and heat of formation to determine the power source. The kinetics of the oxidation are determined by using the Cathcart model for temperatures up to 1800 K and the Urbanic model for higber temperatures [5.23. This kinetics modeling is the same as used in the P'NS module. These kinetic models simply determine the rate of oxide production given an unlimited steam supply. The oxidation model uses this rate as an upper limit. If the steam supply or amount of 2r metal ja insufficient to sustain this rate, then the rate is modified accordingly. The possi.. 'ity of hydrogen blanketing is also considered as a limiting mecnanism. Also, blockages in the bed are assumed to steam starve regions of the bed above their location. 85

I Table 5.1 l Summary of DEBRIS Phase Change Rules Component Melting Temperature Type of Melting l Ag-In-Cd 1070 K Single Component Steel 1700 K Single Componer3t Zr Metal 2200 K Single Component 2200-2650 K-UO Dissolves into Zr(O)-UOp 2 Molten Zr 2900 K Ideal Binary Solution Ir0 -UO2 2 3 86

l-5.5 References i 51. J. E. Kelly, J. T. Hitchcock, and M. L '. Schwarz, " Heat Transfer Characteristics of-Dry Porous Particulate Beds with Inter'nal' Heat Generation", Proceedings _of the - j ASME/JSME' Thermal Engineering Joint Conference, Vol. 4, ~~ U~I5~(I955N~~ 5.2. P. J. Maudlin et al, Light-Water Reactor Degraded-Core 3 Cooling _ Program Technical Note A Damage Assessment of TMI-2, LA-UR-53-737, Los Alamos National Laboratory, Los Alamos, NM,.(1983). 5.3 J. A. Moore, A Mechanistic Debris Bed Packing Model, LA -9 619 -MS (1953), Los Alamos National Laboratory, Los Alamos, NM. 5.4. S. Imura and.E. Takegoshi, "Effect of Gas Pressure on the Effective Thermal Conductivity of Packed Beds," Heat Transfer Japanese _Research, Vol. 3, No. 4, p. 13 (1974). 5.5. D. Vortmeyer, " Radiation in Packed Solids," 6th Int. Heat Trans. Conf., Toronto (1978). 5.6. G. P. W111hite, D. .Kun11, and J. M. Smith, " Heat Transfer in Beds of Fine Particles (Heat Transfer Perpendicular to Flow)," A.I.Ch.E. Journal, Vol. 8, No. 3, p.. 340 (1962). 5.7. A. V,

Luikov, A.

G.

Shashkov, L.

L. Vasiliev, and Yu. E. Fraiman, " Thermal Conductivity of Porous Systems," Int. J. Heat Mass Transfer, Vol. 11, p. 117 (1968). O erations, J. Wiley and Sons, 5.8. Brown and Associates. Unit E New York, NY, p. 223 (1950). 5.9. M. Pilch and P. K. Mast, PLUOM: A_ Coupled Thermal-Hydraulic Computer Model f or Fgeezing_ Melt Flow in a - Channel, NURE G/ CR-3190 SANDE2-15F0, Sandia National Laboratories, Albuquerque, NM, (1983). i me 87

6. RADIATION HEAT TRANSFER 6.1. Introduction A thermal-radiation heat transfer model for analyzing 1 . postulated LWR degraded core accidents must be able to treat a wide variation of core conditions. Among the specific core l conditions that must be modeled are intact

geometry, partial disruption, rubble

bed, slumped fuel

rods, and various combinations of these.

Also, the model must account for absorption and emission of thermal radiation in the water, vapor, and core debris that may exist in the core.

To meet these modeling objectives, a thermal-radiation heat transfer model based on a combination of a net radiation i enclosure mode 1[6.13 and a flux model was developed. Each numerical cell is treated as a net radiation enclosure with internal structures and containing a homogeneous, non-scattering, absorbing, and emitting medium. The internal structures are modeled as surfaces that are diffuse emitters, absorbers, and reflectors. The components in the medium are treated as diffuse emitters and absorbers. The numerical cell boundaries are modeled as " black" surfaces that absorb all incident thermal radiation. The thermal radiation absorbed by'the numerical cell boundary is treated as a directed flux that enters from the adj oining cell, The cell boundary fluxes are the means by which cell-to-cell thermal radiation heat transfer is accomplished. The model treats each numerical cell as three-dimensional, that is three-dimensional view factors are used. The coupling to adjacent ~ cells may be one, two, or three-dimensional.

However, as implemented in MELP ROG-FWR/ MODO the coupling is two-dimensional (r,z).

6.2. Model_ Theory and Design l For an enclosure, the net radiative flux leaving surfaw k in the wavelength interval between 1 and 1 +dA is (6.1) dql,k = dqAo,k dq11,k - total rbdiative flux leaving surface k where dqAo,k in the wavelength interval between 1 j and A +dA and total radiative flux incident on surface k dq11,k = in the wavelength interval between 1 and 1 + dA. 88 i

Also, the total radiative flux leaving surface k

in the wavelength interval between A and A + dA is (6.2) dA

• P A,k dqAi,k
• E A,k 'Ab,k dqAo,k emissivity of surface k at wavelength 1,

= where C A,k black-body emission term for surface k in elb,k the wavelength interval between A and A + dA, and 4 reflectivity of surface k at wavelength A. = pA,k Eqs. (6.1) and (6.2) can be combined to eliminate the dqAi,k term. The resulting equation is dql,k ( 'Ab,k dA

---

(1 pA,k) dqlo,k )

(6.3) PA,k A,k If the term dqAi,k in Eq. 6.1 is expressed as the sum of the radiation leaving all other surfaces that is incident on surface k plus the radiation from the vapor that is incident on surface k, then Eq. 6.1 can be written as a function of dq o,p: g N (dq o,j I,j-k dql,k dqAo,k -[ T -j A k g a j=1 F -j al,j-k ) (6.4) + eAb g dl k F -j view factor from surface k to surface j, = where k T,j-k = mean transmission coefficient for medium A between surfaces j and k at wavelength 1, eAb s - blsek-body emission for medium between surfaces j and k in the wavelength interval dl, between 1 and A + E,j-k = mean absorption coefficient for medium A between surfaces j and k at wavelength A, and N = number of surfaces in the enclosure. 69 89 1 ~^-_-,-__m_...__

Combining Eqs. (6.3) and (6.4) to eliminate the dql,k terms l results in a set of equetions in which the dqin terms are the dependent variables. This set of equations is given by N I,j-k (dqxo,) Fj A dA

• P A,k K

A,k 'Ab,k dqAo,k " C j=1 I,j-k ) (6.5) F -j dA + elb,g A R where k varies from 1 to N. The set of equations defined by i Eq. 6.5 is linear and complete if the source terms elb k and are known and if the view

factors, mean abso,rption elbNicients, coe mean transmission coefficients, emissivitles, and reflectivities are constant.

For these assumptions to be valid, the temperatures of the various surfaces within the enclosure must not change significantly during a time step. In this model, all surfaces are assumed to be diffuse emitters, absorbers, and reflectors and to behave according to Kirchhoff's law. Therefore, the following relationships are valid: (6.6) F -j - A F -k A k j j k (6.7) pA,k = 1.0 - cA,k (6.8) and I,j-k A,k-j A (6.9) I,j-k 1.0 - T,j-k j = A A f depends on the T,djkk E,bfacesand The calculation of A x y an and on the view su j intervening medium between factor between surfaces j and k. Therefore, for a medium with a macroscopic abscrption cross section, ag, for thermal radiation j at wavelength 1, a transmission coefficient can be defined as A (6.10) 3(s) -e t where s is the path length. The mean transmission coefficient is l I then defined to be 90

I l A(s) cosB C088 I 1 T R j dA) dA (6.11) ) t K l j-k 2 A F -k. A <A) Jj v s l-g Whare t i k cosSj cosS g dA dA (6.12) F -k " T ' A ' A j 2 j k J w a k j 10 the view factor between surfaces j and k and A3 and Ak are the oroas of surfaces j and k, respectively. The view factors, mean transmission coefficients, and mean absorption coefficients are calculated as a function of time and space. The numerical cell boundaries are modeled as ' black"

Ourfaces, that is they have emissivities of 1.0 and i

roflectivities of 0.0. Therefore, for the cell boundaries, Eq. 6.5 can be simplified to I (6.13) dq10,k = dqAb,k Whare dqAb.k is the flux entering the enclosure through the cell boundary f rom an adj oining cell. Application of Eq. 6.13 results in a reduction in the order of the matrix that must be solved for occh cell by the number of cell boundaries. However, because is an unknown, an iterative solution procedure sust be dqlb'hortheboundaryfluxes. uned ~ Once the values of dqx have been determined, the net transfer from tSe structure surfaces and from the radiation heat codium can be determined. For surface k the net radiation haat transfer leaving it is given by Q

• A (dq10,k - dq33,g) dl (6.14) k k

l area of surface k, and where A = k 91 l - ~ _ _ - _ _ _ _ _ _ _ _ _ _

1 l i j N l I,j-k (d4 o,j dq11,k;" 1 F -j A K A j-1 E,j-k) (6.15) F -j dA + eg,g A k For. the medium, the net radiation heat transfer leaving it is determined from the net radiation heat transfer leaving-the surfaces. The net radiation heat transfer from the medium is given by N (6.16) l EO Q "~ k g k =1 6.3 Radiation Heat Transfer Model Implementation 6.31. Geometry and Fluid Fields Modeled Two-dimensional numerical cells in an R-Z cylindrical mesh are treated by the radiation heat transfer model. Within a numerical cell, the PWR/ MODO model can treat any combination of the following structures: 1. a wall at the cell inner radius, 2.

a. wall at the cell outer radius, 3

a plate at the cell bottom, 4. fuel rods, and 5. no structures. Additionally, the top surface of a debris bed is treated as a structure surface if it is present in the cell. Water, steam, and corium are treated as a homogeneous mixture in cells without a debris bed, while only water and steam are ' treated as a homogeneous mixture in cells with a debris bed. The above list of structures does not include columns or " general" structures. Columns will be treated in future versions of MELPROG

however, because detailed geometry is not provided fcr

" general" structures they will not be treated in the radiation model. Other structures that might be added to MELPROG in the future will be considered for inclusion in the radiation heat transfer calculation at that time. Figure 6.1 shows typical geometries that are currently treated by the RADIATION module. 92

l l l ( ( 6.3 2 Solution Procedure The radiation heat transfer solution is partitioned into several steps. The first step is to calculate black body I emission terms for each structure surface and a mixture term for ) the fluid fields. Currently only one energy group is used. Therefore,y these terms are cT for each surface, while for the fluids, c7 is calculated for each field and then averaged using a weighting factor based on its macroscopic absorption cross section.

Also, as part of the first

step, the individual macroscopic absorption cross sections and a mixture absorption cross section are calculated.

In the second part of the solution procedure, view facters q are calculated (if there have been geometric changes in a cell) and mean transmission coefficients are calculated. Details of ) the view factor calculation are presented in Section 6.3 3 and l i details of the mean transmission coefficient calculation are ] presented in Section 6.3.4 Part three of the solution assembles the coefficient matrix and calculates a portion of the resultant vector for each cell. The resultant vector calculation is divided into parts because some of the terms that make it up change during the iteration performed in part five. In part four, the coefficient matrix assembled for each cell in part three is inverted. 1 fis s ; an iteration on the cell boundary fluxes is In partThe ;.eration is necessary because the radiative flux performed. entering a cel from an adjoining cell depends on the solution for the fluxes J. that adjoining cell just as the radiative flux calculated to i r. leaving a particular cell affects the solution in the cell that the flux is entering. The first step in the iteration is to calculate the changing part (the part affected by cell boundary fluxes) of the resultant vector for each cell. The solution vector for a cell is then calculated by multiplying the inverted matrix by $ resultant vector. The new cell boundary fluxes (i.e., fluxe

hat are leaving the cell) are calculated and a check on conve. ence is made between iterations.

( For the first iteration the fluxes from the previous time step are used for the convergence check.) The iteration continues until all boundary fluxes have met a user input convergence criteria (usually about 1.0E-5) or a user input iteration limit has been reached (usually about 50). Convergence is normally achieved after four or five iterations, however, when there have been significant geometry changes substantially more iterations may be required. 93

1. r:

U\\ m" ~# O H H l H lIl j [ \\ l y \\ /5 e[l l HW m"! + e l i e n E u T d a "! S A o e R i e L M M nl l l P l O rD [ T T N . d T O k O O I l< ~fQP e P B TA f l L I U E D r S S A S R R E E V W n i O L de LEL t SL a S A er EW v T s e N i h r t e O L m Yl$]/_. E E o R R e O R G I C A B lac i D p A y E I T T m .6 A L 1 P 1 R $ EP e P r hT i F U ug i $f C S A D i i TO i NR I Y e* lll

The final part o? the radiation neat transfer solution is l the calculation of net radiative fluxes to each structure surface and to each fluid field. These calculations are performed for each two-dimensional numerical cell. 6.3 3 View Factor Calculation As stated in Section 6 3 1 not all possible structures are treated by the RADIATION module in PWR/ MODO. The viek factor calculation can be extremely complicated even for relatively simple geometries. To timplify this problem somewhat, certain structures have been neglected (currently columns and " general" structures) and some assumptions about the location of other structures in the two-dimensional mesh have been made. Walls are always assumed to coincide with either the inner or outer radius of a cell. Plates are always assumed to coincide with either the top or bottom boundary of a cell. In cells with fuel rods, it is assumed that there are a sufficient number of fuel Pods so that there is no line of sight path through the rods in the radial direction. In PWR/ MODO the view factor calculation has been organized to treat two types of cells: cells without fuel rods (type 1) and cells with fuel rods (type 2). Within the assumptions made about structure locations and the treatment of radiation heat transfer across cell boundaries, all the view factors for type 1 cells can be calculated using analytical solutions, however, for type 2 cells some approximations must be made. There are three analytical view factor functions used in PWR/ MODO. These include a cylinder's inner surface to an internal cylinder's outer surface, a cylinder's inner surface to itself shadowed by an internal

cylinder, and a

disk to a parallel disk of equal diameter. View factors between surface 2 (see Figure 4.2) and its adjacent cell boundary are defined to be 1.0. Therefore, surface 2 of a structure always radiates to an adjoining cell except when it is on a problem boundary. Walls and plates l (including plates with holes) are assumed to be opaque to all l thermal radiation. Plates with holes will be properly treated in future code versions. The thickness of walls, plates, and debris beds are accounted for in the calculation of view factors, for example, as a debris bed increases in depth in a cell the view factors in that cell will change. For a type 1 cell, only two view factors must be explicitly calculated. These are a cylinder to an internal cylinder and a cylinder to itself shadowed by an internal cylinder. All of the other view factors that are needed for a type 1 cell are determined by applying Kirchhoff's law (Eq. 6.6), by requiring that each view factor be equal to or greater than zero, and by requiring that the sum of all view factors from any given surface be equal to 1.0. Table 6.1 shows how the view factor matrix 95 l

M Le O O. O. O. O. O. O. O. ew aD O e O O O O O O 8e I an a cm 4 K E K Le M. O. M. O. O. M. O. O. ew O O G O C O CD g g Ma = w w h. to H H U H. a O O O O 80 O I d H E E H A A A ea M l l a W L U E o o O. O. O. O. O. O. O. O. m eg a. a O O O e O O O O D U 5 e4 ~ N y C O. O. O. O. O. O. O. O. e 8 U W L O O O O O O O L g 3 3 o oe M a w ae o an 4 w O e4

== A 4 g 30 A O h A O. A. O. O. AD 3 L g N O 3.. g 9 O O O .O O 5 E 4 W trj In. A In. A >m cc to.

• c 3 eo b

ec 5% !O" x E E *e b U ea O. O. O. O. O. O. O. O. L ec w L g e 4 c i. E e = g O O O O O O O h $g e36 su u O ha eMaa e NgC a aM MMM y SEREa5 U gy e E. O. E. O.. O. E. O. e m = O e===c R execce O O O z h 3 n. M G M k k k lm. M& NN N N S tEEE y g e g - - so en Ms3AHK w N e N e3 h3 E d d d d a L L s. i o a c e, av r 3 5 2 3 O O ce 3e O M M M t/) 8D H MM QR E E Es e bM s a. a aM A 1 I 1 1 i I 96

would look for a cell with an outer wall and a bottom plate (for example a cell with the core barrel and diffuser plate). If this sample cell had a debris bed in it in addition to the other s rfaces, the view factor calculation would be performed in the same way except that all of the view factors to surface 1 of the bottom plate would be to the debris bed surface. The view factor for surface 1 of the bottom plate to itself would be defined to be 1.0. In a type 2 cell (a cell with fuel rods), there are some l approximations that must be made to utilize the three analytical I vtew factor solutions discussed above. First, the fuel rods are { assumed to approximate a annulus whose outer diameter is one-half tme fuel rod pitch smaller than the outer cell radius or the imner radius of an outer wall and whose inner diameter is ome-half the fuel rod pitch larger than the cell inner radius or the outer radius of an inner wall. Second, the fuel rods in a cell are assumed to be of sufficient number that there is no line of sight throtagh them in the radial direction. Third, the view factor from the cell bottom to the cell top through the fuel rods is approximated by the view factor from a disk, having an area equal to the total cell bottom surf ace area divided by the number of fuel rods in the cell minus the area of one fuel rod, to a disk having the same area at the cell top. By using these three approximations, together with Kirchhoff's law (Eq. 6.6), the requirement that each view factor be equal to or greater than zero, and the requirement that the sum of all view factors for a surface be 1.0, a consistent set of view factors for type 2 cells are calculated. This procedure is valid for PWR whole core calculations; however, for BWRs or small bundle experiments the approximations are invalid. A more detailed numerical view factor solution is required to adequately treat BWRs and small fuel rod bundles. The view f actor calculation in PWR/ MODO is dynamic, that is, whenever there are geometry changes in a cell the view factors rer that cell are recalculated. Examples of situations when view factors are recalculated are when a structure melts, when fuel rods balloon or disintegrate, and when debris beds form or change in thickness. 634 Material Properties and Mean Transmission Coefficients The RADIATION module uses constant emissivities for each structure surface and for debris bed surfaces. The emissivities are user input. The water and corium mass absorption cross sections, which are user input constants, are also required. For steam, a temperature and pressure dependent function is used to calculate a Planck mean absorption cross section. Hydrogen is treated as transparent to thermal radiation. In future versions of MELPROG, temperature dependent emissivities will be provided. 97

l 1 i The' current version treats only~one energy groups however, the model is designed to treat multiple energy groups and could do so with additional emissivity and absorption cross section data and this option may be included in future versions. The calculation of mean transmission coefficients using l Equation 6.11 is extremely difficult and can not be done analytically for the geometries considered in the RADIATION module. Therefore, Equation 6.10 is used with an average path length (s). The path length calculation is dependent on the cell -geometry and varies as the cell geometry changes during an accident calculation. 6.3 5 Summary of Basic Modeling Assumptions The RADIATION module assumes that all surfaces are diffuse emitters, absorbers, and reflectors; that the fluid fields are homogeneously mixed: and that radiation across numerical cell boundaries is diffuse. Furthermore, it assumed that it is not necessary to treat every surface,. that fuel rods are always vertical and in sufficient numbers to prevent line of sight viewing through them, that there are no holes in plates or walls, that all plates are flat, and that all walls are cylindrical. Also, it is assumed that all temperatures are constant within a time step and, therefore, do not change substantially between time steps. 6.3.6 RADIATION Module Interaction with other Modules The RADIATION module interacts with each of the other calculation modules in MELPROG. From the' FLUIDS module, it gets densities and temperatures and it returns volumetric energy 'i transfer rates. From the PINS

module, it gets fuel rod temperatures and dimensions and it returns radiative fluxes.

From the STRUCTURES module, it gets geometry and temperatures and it returns radiative fluxes. Finally, from the DEBRIS module, it gets debris bed geometry and temperatures and it returns radiative fluxes. 6.4. References 6.1. R. Siegel and J. R. Howell, Thermal Radiation Heat Transfer, McGraw-Hill Book Company, pp.550-558 (1972). e 98 i-

w k 7.: SAMPLE CALCULAT20N- .7.1. Introduction As a demonstration' of the. capabilities of MELPROG,.a sample ' calculation.has been perf ormed. The -intent of this-effort was -two-fold. First, ' no' f ully-coupled meltdown calculation had ever been performed. Consequently, this. case served.to demonstrate thati MELPROG could Lperf orm. such a calculation and also to show the' importance of performing a coupled calculation.

Second, while the-various modules in MELPROG had.been tested separately, the ' combined code and the' coupling between modules had not been fully tested.

Hence, this problem served to' identify and resolve inconsistencies between modules. 'The test problem which was chosen was a simulation of an S1D accident sequence. This sequence represents a medium-size loss of coolant accident which results in rapid depressurization - and It is assumed that emergency core cooling systems core uncovery. '(ECCS) are disabled so the core eventually melts down. Since the purpose of the -calculation was mainly _ to demonstrate the code capabilities, simplifications to both the geometrical .model and accident sequence were made. In .particular,. since this problem involves a rapid loss of the primary coolant, it was felt that after the initial transient, feedback from the reactor coolant system (RCS) to the vessel. would be minimal, and that the reactor vessel could be decoupled from the remainder of the RCS without substantial error being introduced. Furthermore, a f airly detailed RELAP-5 aisulation of the loss of coolant phase was available.to compare the transient with[7.1h Other simplifications were made to the geometrical model of the reactor vessel and to the boundary conditions for the problem. In-the following sections, a description of the calculation is given. This includes a discussion of the geometrical model which was used and a summary of the calculated results. 7.2. Pro m e Setup, As noted above, the calculation was simplified such that the modules in MELPROG would be exercised without unnecessary cost. The simplifications involved using both a simple geometrical model as well as using approximations for the boundary conditions. The simplifications to the geometrical r.odel l *~ involved several aspects. First, the model only included those portions of the reactor vessel from the vessel lower head to the upper grid support plate. This means that the upper plenum was not modeled. The downcomer was modeled, but due to the 1-D 99

nature of the FLUIDS module, the flow of liquid and vapor was always in the same direction'as in the core region of the vessel. j Hence, the heat transfer in the downcomer could not be accurately ) modeled. (This problem could be circumvented to a certain extent ] by noding the downcomer separately, tying its thermal boundary i conditions to the core barrel and vessel on its inner and outer {

surfaces, respectively, and reversing gravity.

However, this solution would still not be totally adequate.) The problem was set up with four radial rings (three in the core region and for the downcomer) and seven axial nodes (five in the core region and two in the lower plenum).

The vessel nodalization model for this { test problem is illustrated in Figure 7.1. The size of the volumes is large, but sufficient detail is present to evaluate the performance of the code. Furthermore, the heat transler in the DEBRIS and STRUCTURES modules use a finer mesh for those processes requiring it. Structures modeled include the vessel t wall and vessel bottom head, the core barrel, the core support J plate, support columns (not shown) -and upper and lower core grid plates. Also, instrument penetrations of the lower head were not treated in this calculation, although these may be included at the option of the code user. In performing this simulation, the link to TRAC-PF1 was not utilized in order to avoid unnecessary complications.

Instead, the Reactor Coolant System (RCS) was represented by a s'e t of q

upper-plenum pressure, liquid fraction, and temperature boundary conditions obtained from a RELAP5 calculation of this i a c c i d e n t [ 7.1 ]. A simulated flew coastdown was used for the lower boundary condition. These simplification have repercussions on the details of the calculated sequence, since no active feedback to the RCS hydraulics from in-vessel processes is included. In particular, any material crossing the upper boundary was lost j from the calculation.

However, the use of these boundary conditions presents no limitations on the exercise of individual modules, nor on the overall code itself.

The SID sequence involves a high pressure blowdown through a 66-102 mm diameter break with scram and pump trip. The upper plenum pressure behavior can be characterized as a sharp drop to saturation pressure over about a 50 second time span followed by continuous steam generation and slow pressure drop as water flashes out of the reactor. Since the upper plenum was treated as a boundary condition, it was necessary to simulate the voiding of the upper plenum. This was accomplished by linearly decreasing the liquid fraction at the upper plenum boundary to l zero in about 80 seconds. This rate of voiding was arbitrarily / chosen and will certainly have some influence on the timing of { events during the accident sequence. 100 l L__ _ i

i l l y~ m y yg ,- 3 1, .c,- 7 l l l 1 ....g.. i.. +.. < .. p.,..g... ...i..., ... i..< I e l i j i l l ( . 4., . 4... i ...t.....;...... 8 l a l El E! E E H.Eit i = = = = = =. ..9.....g... ..j... 9.....< 4 l

• YESD. IAU.

...j ....J,.....:. ( 1 .. 1.. ....Y.... i 3 my ga pun 7,7 w m rirss,w n usu usrr

wvivi, i

l ) 2 l YH/M, / M/H/HHHH//MHMH/A =.. ante suPEwt 9 LATE I I i I + l 1 e g w 'v a a l e/ w i r a / u u ea u / a s/ w /r/a u ro e.--VESSEL 30Tf0M j 4 3 2 1 1 2 3 4' j It!NGS 4 l Figure 7.1 SID Vessel Nodalization ] 101

7.3 Results of S10 Test Problem The results of this test problem are summarized in the following sections. The results presented are intended simply to illustrate the calculational potential of MELPROG. In order to clarify the results, the details of the calculation are discussed on a module by module basis. Table 7.1 hi ghl i gh t s the calculation for each of the in-vessel modules (except for the RADIATION module). 'This summary gives a good indication as to where and when various phenomena are being calculated. The main events in each module t re discusred below. 7.3.1. FLUIDS Module As discussed in Section 2, the FLUIDS module predicts the fluid dynamic behavior in the problem. In particular, it calculates the volume fraction, temperature, and velocity of each field as a function of time and space. Since the modeling is one-dimensional, only radially averaged quantities are calculated. The main events calculated by FLUIDS are listed in Table 7.1. In this table, the maximum temperatures whien are listed correspond to the tiaximum corium temperature in any node. Of particular inter (st during thfe test sequence is the production and transportt of vapor. Since the hydrodynamic solution is directly derived from the TRAC code development effort, it is felt that the production and transport of vapor are well predicted by MELPROG. The prediction of these quantities for this calculation are best illustrated by the vapor volume fraction distribution shown in Figure.7.2. (Note that the sum of the liquid, vapor, and corium volume fractions is always equal to 1.0 and the point Z=0 is at the bottom of the vessel.) Flashing of water to steam is predicted to begin relatively early in the sequence when the system pressure drops to the saturation pressure. The lack of a model for the upper plenum accelerates the events in the sequence by at least 400 s. Core uncovering begins shortly after flashing and at 560 s the core is 80% uncovered and the maximum steam temperature is 1100 K. The steam continues to heat and fuel rod oxidation begins shortly after this point. The water in the vessel continues to vaporize due to the depressurization and heat transfer and at 900 s the core is completely uncovered. From this point until later in the accident, the vapor behavior is benign as the water in the lower plenum slowly vapor 12es.

However, at 3375 s, grid plate failure leads to a

1 debris-coolant interaction resulting in rapid vaporization and a ~ j steam spike. As discussed previously. HELP ROG-PWR / MODO does not include a mechanistic fuel-coolant interaction model.

Rather, heat transfer between the debris and coolant is calculated, but no rapid fragmentation is accounted for.

102 l

er R R de u l 0 0 e e e i 6 0 l N l a 4 8 l S l F 1 1 a m ER a = F K o U R F l e e R K K t T K e t t e 0 t C 0 e r a a t 0 0 0 0 o 0 t r l l a 0 2 5 4 B U m 1 a a p p l 5 5 6 1 R 1 l B P = = = = l T P d d d d d d e =, =, S e l i d e a e a s hhs , p r e r i e e e , o o g g r h h T 7 T C T T G T T T T V n s o e n t n i t n i g e o g e o s n t e n n mB p m e t B e e m r t s o t t n g o o l e B g l l mi o n C F o h n o o M c n i M M it 5 t d l n I t n o 4 e l 4 Ml a 1 e a o 0 u e e 3 3 l R M P 7 Q m M 6 7 9 u E g r c D h e n n s e o e e s l i e i i P i i i l e e t t r r r r r l l b b d b b b a Mbe e o e e e e e e C D M M D D B D D D D 1 l S e 7 n n i e i ge l s B b s r s a t u n e c i r e d d d d d d T n e e e e e e c g e h e l i l l l l 0 e 1 v t l i i i B d 1 Mle a e a a a E g e e F F F F F F n n l F i o t f n t t d g e s s e s e t o o t i o n n n n n n n S o a r m i i uiiii P P P P P F N l s b y I l t n l Mf r P e s a o f f f f f o o o o o o a B O r e t m n n d n &0 5 6 2 0 2 m i o o 2 4 6 7 8 9 u P P m C S sn y i r g d D e d e s B e r r y r e u l e g e e c e n n e o e t i i o c O e g r c n l e e n E U K K R R e p K K K K R S e U k m t 0 0 0 C t o 5 0 0 0 0 S o Mt D g c 2 0 0 0 p C 8 0 5 2 9 e 2 4 4 4 8 2 5 4 7 7 I n n 0 Ul 2 2 2 2 l 2 2 2 2 2 U i U 9 =, =, =, =, =, =, =, =, =, =, n e L h e e e F s e e ,, e e T. T, 7, T, T, ,r ,, t w a r r o l o o F C C T C T T T T S T E) 0 5 0 0 0 5 5 0 0 0 0 0 0 5 0 0 0 0 0 0 Mo 5 3 0 6 0 1 2 0 0 0 0 0 0 7 9 0 0 3 0 0 I( 1 4 8 8 8 9 2 9 3 7 2 3 3 4 7 4 4 9 T 1 1 2 2 3 3 3 3 3 4 6 6 pow

i J l ] ~ I L8" >I H 2e. n W w~yg 00-l 3 510 Time (5) M E'4 l f I l Figur>e 7.2 Vapor Volume Fraction Distribution l l e n e 104

At the time of this

event, the lower plenum contained approximately 12500 kg of saturated water at a

preocure of 1.65 MPa. Approximately 36000 kg of molten debris with a saximum i temperature of 2500 K interacted with this water. To calculate this event, extremely small time steps were required, as is typical for hydrodynamic codes under conditions of rapid phase I change and shock generation. The quenching of the debris and subsequent boiling of the water led to complete voiding of the vessel at 3380 s. It should be noted that one may realistically expect a vapor explosion to occur under these conditions. The lack of fuel-coolant interaction (FCI) models in PWR/ MODO prevents the code fro; calculating the effects of such an explosion. This shortcoming will be removed from subsequent code versions which will incorporate detailed FCI mo tle l s. In addition to calculating the vapor fraction distribution, the FLUIDS module also calculates the debris fraction distribution. As seen in Figure 7 3, the debris begins to form at 800 s. As the debris forms, it moves down and freezes on top of the grid plate. It continues to build up on the grid plate. After grid plate failure, the debris moves into the lower plenum where it again builds up. The debris temperatures during this sequence are illustrated in Figure 7.4 It is seen that the debris approaches an equilibrium condition before grid plate failure. The coupling of the debris heat transfer to the structures allow for this equilibrium condition to occur. j Overall, the FLUIDS module performed quite satisfactorily for this test calculation. No major numerical pr obl em s were encountered during the calculation. In fact, only two events caused the pressure solution method to fail. The first was the rapid vaporization caused by the debris-water interaction. The i second was the rapid H generation caused by Zr oxidation. In l p both cases, by using a small time step, convergence could be obtained. Both of these problems indicate that some portions of ~ the pressure iteration are explicit in nature and small enough time steps must be used to prevent divergence. l i The main difficulty with FLUIDS is its one-dimensional restriction. This limitation makes it impossible to model the important two-dimensional effects such as natural circulation in the vessel. Also, it is difficult to model the entire vessel consistently. For example, the downcomer cannot be modeled in a convenient manner that would allow proper treatment of flow anc l heat transfer. Another problem is tracking the debris on a two-l dimensional basis.

Clearly, the debris is not formed or l

transported uniformly across the entire core. Rather, there is a strong radial dependence to its formation and transport. The current FLUIDS module cannot accurately determine this dependence.

Hence, there is a

strong need to have a two-dimensional FLUIDS module. e 105

1 I a-b / \\ i $e >/f^N j S S" ) l I

1. m, (s) fi ss'

<a 2(g)4 Figure 7.3 corium Volume Fraction Distribution I l 106 4 f

g10* I l W 0 p G y9~ J M af8 P* ff l Yg(5) SO' W q 2(kf* Figure 7.4 Corium Temperature Distribution l.. 107 L----_-__-_-

7.3 2. PINS Module As discussed in Section 3, the PINS module is responsible for predicting the fuel and control rod behavior. Specifically, in each core region ring, a representative fuel and control rod are modeled. A temperature distribution is calculated in the fuel rod and the control rod is assumed to be at the coolant temperature. When the representative fuel or control rod in a given node is predicted to fail, all the corresponding rods in that node fail and their associated mass and energy are transferred into the debris field in FLUIDS using the mass transfer functions discussed in Section 2. In the test problem, the control rods begin to fail at 81 5 s when the vapor temperature reaches 1700 K. The mass and energy of the. control rods (i.e., the steel cladding and absorber material) is transferred to the FLUIDS corium field via a mass transfer function. As the sequence continues, the control rods continue to fail throughout the vessel. The fuel rods begin to overheat as soon as core uncovering is initiated. Figure 7.5 illustrates the cladding surface temperature distribution for the rod modeled in ring 1. ( A value of 0K for the temperature indicates that there are no rods at that location.) This figure shows the rapid heating of the rod as cooling is lost and Zr oxidation begins. Ballooning in the top of the core occurs at 400 s. Oxidation and subsequent rapid heating of the rods begin to occur at 560 s. The rods become embrittled at 800 s. At this point, a quench would shatter the rods. Melting of the Zircaloy cladding begins at 815 s and the first rod failure occurs at 825 s by complete melting of the clad. Following this point, rod failure occurs in other cells as the sequence continues. Rod failure can be seen by a sudden decrease in temperature to zero. At 1900 s, 60% of the rods have failed. The model continues to calculate the behavior of the remaining intact rod stubs until complete failure of all the rods occurs. The PINS model seems to perform satisfactorily up to the point of rod failure. The modeling is not overly detailed, but most of the key phenomena are accounted for. Since the heat transfer is calculated using models from the TRAC code, it is felt that this phase of the accident is predicted reasonably well. The difficulty with the model is that it assumes that when cladding failure occurs the fuel rods are fragmented and the debris which is formed is of a particulate nature. Recent experiments tend to indicate that the fuel rods may not fragment upon cladding failure [7.2]. This fragmentation assumption is a carryover from MIMAS where the principal application was the analysis of TMI-2. Hence, the major weakness in the rod model is that it does not properly treat the candling phenomena.

However, i

108

W~ f j ~ U a-U g p s00 y-g as8 trne (s) (% i Figure 7 5 cladding Wall Temperature in Ring 1 109

l 1 a candling model is being developed and will be added in future code versions. 7.3 3 DEBRIS Module As discussed in Section 5. the DEB RIS module is responsible for predicting the debris behavior after the debris is created. While FLUIDS calculates the average

mass, temperature, and velocity of the debris on a

one-dimensional

basis, DEBRIS calculates detailed temperature and compositional distributions in the debris on a two-dimensional basis.

This information is important for providing the correct coupling of heat flow to supporting structures such as the grid plate and vessel bottom head. In the current problem, the debris formed by rod failure is predicted to move downward and freeze in the lower sections of the debris does'not move through the grid plate the core.

Hence, until grid plate failure begins.

The debris accumulates on the grid plate and the first debris bed is formed at 825 s. At this point, the DEBRIS module begins to calculate. Core uncovering leads to bed dryout and the debris begins to overheat. Since there is steam production below the bed, Zr oxidation occurs in the

bed, By 900 s, debris melting begins leading to some compaction of the bed.

As rod failure continues, the resulting debris is added to the bed and the radial extent of the bed increases. The melting, relocation, and freezing of the debris in the bed lead to blockages in the bottom of the bed which effectively insulate the molten regions from the grid plate. At 2300 s, a molten pool is predicted to exist in the bed. The downward heat flux from the bed to the grid ultimately results in grid plate failure at 3375 s. At this time., there is approximately 130000 kg of debris on the grid plate and of this amount, 755 is molten. The composition of the debris, based on the DEBRIS module phase diagram, is 12500 kg Zr, 3500 kg Zr02 and 114000 kg UO where all of the UO is in the form of the 2 2 eutectic. As stated above, approximately 36000 kg of Zr(0)-UO2 the debris interacts with the water in the lower head to produce a steam spike and to quench some of the debris. However, most of the debris remains hot and the remaining water quickly boils away.

Hence, the debris is left uncooled and melting and relocation are predicted to occur again.

At the time the vessel bottom head failure (6900 s), nearly 735 of the debris is again molten. Overall, the DEBRIS module performed satisfactorily. The heatup and meltdown of the debris seemed to be consistent. The main problem with the module lies in the debris bed formation and rod release model.

Again, the lack of a rod candling model

^ forces the model to treat the debris as if it were in the form of 110 L ____ __ _-_ _ _

a particulate. This only has a major effect on the formation model since the rules for bed formation will be different for the two situations (i.e. particulate debris versus candled debris). However, as the debris melts, the model does account for the cohesive nature of the debris. Consequently, the formation model to be expanded to handle the candling phenomena. needs 7.3.4 RADIATION Module As discussed' in Section 6, the RADIATION module is responsible for predicting the radiative heat transport in the vessel. As such, it calculates the radiative energy exchanged between the various modules. Although maj or events cannot be associated with this

module, radiation heat transfer significantly affects the heat flow distribution and timing of in the vessel.

In fact, in this calculation it is the events For primary medule which couples all the other modules together.

example, the failure of the outer ring of pins is delayed lost by these outer relative to the inner rings due to the energy pins to the core barrel.

Furthermore, for stagnant flow conditions the primary heat transfer mechanism to the core barrel and top grid plate is radiation. (Recall that with a 1-D fluid dynamics solver, natural convection is ignored.) Radiation heat transfer also acts to cool the top surface of the debris bed and strongly impacts the heating of the steam. Another important feature of this module is its ability to modify the view factors structures and pins fail. That f eature worked well throughout asthe test calculation. 7.3 5. STRUCTURES Module As discussed in Section 4 the STRUCTURES module is responsible for calculating the thermomechanical response of the structures in the problem. In each structure modeled, a detail j heat transfer calculation is performed and, based on this, the j l mechanical response of the structure is determined. Since there is little or no heat source in any_ structure, the the temperature increase of a given structure is governed by heat transfer from the vapor to the structure and by radiation heat transfer. This leads to a lag in the heating of the l As seen in Table 7.1, structures relative to the vapor and pins. the maximum temperature of any structure is 1000 K at 1200 s. Initially, the core barrel is the hottest structure. The temperature distribution in the core barrel is illustrated in Figure 7.6.

However, the core barrel does not begin to fail until 2300 s.

Rather, it is the top plate which mechanically o

f ails first at 1900 s. I 1 l i 1 1 111 i a

i I l 1 j-1 ggDD" \\ Eg># yson ~P 4 I so88 as s C in (8) e 4 (% Figure 7.6 Core Barrel Temperature Distribution I 112 L ---- ------

The grid plate heat up is controlled by heating from the debris bed on top of it and from cooling by steam flow from below. Average temperatures of the grid plate-at various times are given in Table 7.1. Failure occurs at 3375's due to melt ) through. At this point.in the sequence, the lower vessel head is still at the saturation temperature and, hence, is relatively

cold, The. increase in lower vessel head temperature (Thead) is given in Table 7.1.

Failure of the vessel head occurs by melting at 6900 s.- At this point the problem was terminated. The STRUCTURES module also performed satisfactorily for this problem. However, two possible inconsistencies were identified.

First, structures only failed when their average temperature reached the melting temperature.

This may or may not be the actual failure mode. Consequently, further study of this situation is warranted. Second, the melted structures were not transferred into the FLUIDS debris field. For consistency, this transfer should occur and will be incorporated into future code versions. 7.3.6. Summary On a integral basis, at the time of vessel bottom head Nd been evolved which represents oxidation fa13ure, 486 kg of H2 of approximately 50% of the Zircaloy in the core. Approximately 106000 kg of molten debris would have been ejected which represents 735 of the core inventory. While these numbers and the timing of specific events are certainly not a best-estimate due to the simplifications

involved, the results of this calculation demonstrate the potential capabilities of MELPROG.

Overall, the performance of the modules in MELPROG for this S1D test problem calculation has demonstrated that the current code contains most of the features necessary to calculate a meltdown progression. The code is able to analyze in an integrated manner the important in-vessel phenomena anticipated for this accident sequence with the exception of fission product release and transport, which is currently being added. The code for this test problem. The timing performed largely as designed of events as well as the interaction of the various modules seemed to be consistent. While the simplifications in the modeling of the reactor core and vessel mean that this particular calculation is not a best-estimate, it should be clear, that the major components for a best-estimate calculation currently exist in the code. One item that should be noted is the neglect of the accumulator water, which was not modeled in this test problem. For this accident sequence, the accumulator would have discharged at approximately 1100-1200 s into the sequence, when the vessel + At this point in time, the core is pressure dropped to 4.5 MPa. 113

uncovered and has i arted to oxidize. The main effect of the accumulator water sould probably be to enhance solid-state breakup of the fuel rods, as in TMI-2, as well as delaying the subsequent events in the test problem. 7+4+ EfffCggggs 7.1. S. M. Modro et al., "RELAP5 Analysis of LOFT and Zion Nuclear Power Plant Small ' Break LOCA's," Proceedings of the International Meeting on Thermal Nuclear Reactor Sa f et_g, Chi ca~go, IL, NUREG/CP-0027, Vol.2, p.839 (1982). 7.2. J. B. Rivard et al., DF-i l Results and Analys i s,. S AND85 -07 82, Sandia National Laboratories, Albuquerque, NM, (to be published). v S l 114

Appendix A. Input Instructions A.1. Introduction In this appendix, the input instructions for MELPROG are given. Both the format of the input as well as descriptions of the' input variables are discussed. If there is a recommended value for a particular variable, this is noted by enclosing the value in parentheses following the variable description. A.2. General Description The following describes the input as of December 1, 1984. The input to MELPROG is in free format and uses the FUNRD[A.1] input processor. Numbers are only read as reals. A " ' " in column 1 of the input line indicates a literal. A" $" in column 1 indicates a comment card, which may be inserted anywhere in the input deck to add readability. n8v", where a Repetition of values may be done with the form n is the number of times the value v is to be read. Separators between values can be blanks or commas. Repetition of groups of

• • • )'

where n is values is done with the form n'(mj'vj, m2'V2' the number of times the group in () is repeated and the form n 'v means repetition of v n times. x g x g A restriction on the free form input is that the input buffer in MELPROG can take only 40 values per lines therefore something like "60'0.01" should be broken up into 40*O.01 on one line and 20'O.01 on the next line. ~ A.3 Input cards i A. 3 1. General Input 1. ITITLE (20A4) A. ITITLE = problem title (up to 80 characters). i 2. IRESTRT A. IRESTRT = the restart switch 0 means no restart. a. IRESTRT = b. IRESTRT = N means read restart dump N from file (MIMDMP). c. IRESThT = -1 means read final restart dump from file (MIMDMP). l** l 3 DMPINT A. DMPINT = restart dump interval (sec). f,, l l 1 115 l t

NOTES DMPINT is not input when MELPROG is linked to-TRAC. 'The TRAC restart dump interval is used when linked to TRAC. 4. IPRNT . number of time' intervals between A. IPRNT = printouts. 5. IPRTF, IPRTP, IPRTR, IPRTS..IPRTD A. IPRTF = fluids print. flag (0=off, 1=on). B. IPRTP = pina print' flag (same). C. IPRTR. = radiation print flag (same). D. IPRTS. = structuree print flag (same). E. 'IPRTD = debris print flag-(same). 6. KMAX A. KMAX = number of axial mesh cells 7 NRING A. .NRING = number of radial' rings 8. NCORE1 A.. NCORE1 = fluids cell number for core bottom cell 9. NCORE2 A. NCORE2 = fluida cell number for core top cell

10. LRAD A.

LRAD = radiation switch (0=off, 1=on)

11. TSTART A..

TSTART = the starting time (sec).

12. ENDTIM A.

ENDTIM = problem end time (sec)

13. DELTAT A.

DELTAT = the initial time step (sec).

14. NTIM number of time step pairs in time step A.

NTIM = table (22).

15. (TIM (I), DTIM(I). I-1,NTIM)

A. TIM (I) = Ith time at which the following time-step [ begins (sec). l Ith time step which begins at time TIM (I) B. DTIM(I) = (sec).

16. SS pin model steady state switch (0=0n, A.

SS = 1=0ff). If-this is set on, MELPROG will attempt to y reach an equilibrium state for temperature and fluid flow. NOTE: If the equilibrium temperature for the particular set of initial conditions is high enough to cause significant transient effects, e.g., clad oxidation, no equilibrium will be reached the problem should be run as a transient starting with the input as initial conditions.

17. QLIN A.

QLIN = total reactor power (w).

18. N0Q l

A. NOQ = number of pairs in the power table.

19. ((TIM (K), POW (K)), K-1,N0Q) 116

A. TIM (K) = time point K in power history table.(sec). multiplies B. POW (K) = relative power at time K total reactor power. A.3 2. Fluid Dynamics Input A. 3. 2.1 Fluids Scalar Input

20. ERROR 1 Error cri ter f.on 1

for Newton-Raphson A. ERROR 1 i = iteration. This is the relative convergence criterion on the pressure and void fractions.(0.001-0.1) l

21. ERROR 2 Error criterion 2

for Newton-Raphson A. ERROR 2 = iteration. This is the absolute convergence criterion on temperature in degrees K. (0.1-10.K)

22. ERROR 4 time step expansion and contraction A.

ERROR 4 = parameter. This is actually the percentage of the current time step by which the time step can be reduced or increased in one step.(0.0"99.0). ]

23. EPSA maximum fractional change in alpha over A.

EPSA = time step (0.05-0.1). 24 EPST A. EPST = maximum fractional change in temperature over time step (0.01-0.1).

25. IFLOOD A.

IFLOOD = reflood switch a. IFLOOD = 0 means no reflood b. IFLOOD = 1 means reflood ~

26. ITERMIM A.

ITERMIM = minimum number of Jacobi iterations (5).

27. ITERMAX A.

ITERMAX = maximum number of Jacobi iterations (20).

28. AL10 low alpha limit for vapor (1E 1E-3).

A. AL10 = If the volume fraction of vapor is below this value the field is turned off. l

29. AL20 low alpha limit for water (1E 1E-3).

A. AL20 = If the volume fraction of water is below this value the field is turned off.

30. AL30 low alpha limit for corium (1E-6 " 1E-3).

A. AL30 = l If the volume fraction of corium is below this value the field-is turned off. l

• =

31. ETA l

A. ETA = not currently used but still read ~~

32. BETA 117 l

l

A. BETA = Courant time step factor (2.0-10.0). .This is the factor by which MELPROG is allowed to exceed the Courant condition. 33 IFILL inlet boundary condition flag. A. IFILL = a. IFILL = 0 means loop boundary condition is to be specified. b. IFILL = 1 means velocity boundary condition is ~ to be specified. A.3 2.2 Fluid. Dynamic Array Input

34. (DELTAZ(I),I-1,KMAX+2)

A. DELTAZ(I)= mesh cell length for cell I (a).

35. ((FRCTN(I,J),I=1,KMAX+2).J-1,3)

A. FRCTN(I,J)- additive friction factor for cell I'and field J.

36. (ABYP(I),I=1,KMAX+2)

A. ABYP(I) = bypass flow area for cell I (m*

• 2).

This must be greater than zero in cells having no pins. Note that setting the bypass flow to zero does not shut off heat transfer from the radial wall to the bypass channel.

37. (NPINS(J),J-1,NRING)

A. NPINS(J)= number of fuel rods in radial ring J.

38. (NCNROD(J) J-1,NRING)

A. NCNROD(J)= number of control rods in radial ring J.

39. (NPOISN(J),J-1,NRING)

A. NPOISN(J)= number of poison rods in radial ring J.

40. (PACDEN(I),I-1,KMAX+2) (KMAX+2 values)

A. PACDEN(I)= maximum packing fraction for corium debris in cell I (0.44). i

41. (RI(J),J-1,NRING), (RO(J) J-1,NRING) (2'NRING values)

I A. RI(J) = inner radius of ring J (m). B. R0(J) - outer radius of ring J (m). A.3 2 3 F1' aids Initial Conditicus A.3 2 3 1. Field One ' Steam and Mydrogen i I

42. (VG(I) I=1,KMAX+1) (KMAX+1 values)

) initial vapor velocity for cell I starting ] A. VG(I) = at cell I-1 (m/s). i

43. (TG(I),I 2,KMAX+1) (KMAX values)

) initial vapor. temperature for cell I A. TG(I) = starting at I-2 (K). I 44 ( ALFG(I),I=2,KM AX+1) (KMAX values) initial vapor void fraction for cell I A. ALFG(I) = starting at I-2. l 118 j I L

1 l A.3 2 3 2. Field Two - Wat3r

45. (VL(I),I-1,KMAX+1) (7 MAX +1 values)

A. VL(I) initial water velocity for cell I starting = at cell 1 (m/s).

46. (TL(I),I=2,KMAX+1) (KMAX values)

A. TL(I) initial water temperature for cell I = starting at cell 2 (K). J

47. (ALFL(I),I-2,KMAX+1) (KMAX values) 1 A.

ALFL(I) initial water void fraction for cell I = starr.ing at cell 2. A.3 2.3 3. Field Three - Corium I

48. (VC(I),I-1,KMAX+1) (KMAX+1 values)

A. VC(I) initial corium velocity in cell I starting = in cell 1 (m/s).

49. (TC(I),I-2,KMAX+1) (KMAX values) initial 'corium temperature in cell I

A. TC(I) = starting at cell 2 (K).

50. (ALFC(I),I=2,KMAX+1) (KMAX values)

A. ALFC(I) initial corium void fraction in cell I = starting at cell 2.

51. (P(I).I=2,KMAX+1) (KMAX values) i initial pressures (Pa).

A. P(I) = NOTE: it is best to start the problem with axially uniform initial velocities if generating a steady-state solution. A.3 2.4 Fluid System Boundary Conditions ~

52. NVTAB number of inlet velocity boundary A.

NVTAB = condition table pairs. 53 ((TIM (K,J), VB(K,J)), K=1,NYTAB), J-1,NFIELD) A. TIM (K,J)- time point K for field J in velocity table (sec). B. VB(K,J) inlet velocity at time point K for field J = (m/s). NOTE: To avoid evaporation / condensation problems when either the steam or water void fractions are very small, set the inlet velocity of the field with the small void fraction to a small value (1ea10).

54. NTTAB A.

NTTAB = number of inlet temperature table pairs. a

55. ((TIM (K,J), TB(K.J)), K-1,NTTAB), J-1,NFIELD)

A. TIM (K,J)- time point K for field J in temperature l table (sec). inlet temperature at time K for field J B. TB(K,J) = (K). l l 119

56. NATAB A.

NATAB' = number of inlet alpha table pairs.

57. ((TIM (K,J), ALF(K,J)), K-1,NATAB), J-1 HFIELD)

A. TIM (K,J)= time point K for field J in alpha table (sec). B. ALF(K,J)= inlet void fraction at time K for field J.

58. NPTAB A.

NPTAB = number of pressure drop table pairs. 59.-((TIM (K), PH(K)), K-1,NPTAB) time point K for pressure drop table A. TIM (K)- = (sec). B. PH(K) = pressure drop at time point K (Pa). This pressure drop is applied between the inlet and outlet boundaries of the problem if the inlet loop ' boundary condition is set (IFILL=0).

60. NPBTAB number of outlet pressure boundary A.

NPBTAB = condition pairs.

61. (TIM (K), PB(K)), K-1,NPBTAB) time point K for pressure boundary table A.

TIM (K) = (sec). outlet pressure at time K on boundary B. PB(K) = (Pa). A.3 3 Rod Model Input A.3 3 1 aod Model scalar Input

62. NOFUEL number of fuel radial nodes.

A. NOFUEL =

63. NOCLAD number of clad radial nodes.

A. NOCLAD = i

64. TLM Temperature error criterion (10.0K).

{ A. TLM = {

65. CLADO outer cl.ad diameter (s).

j A. CLADO = ]

66. CLADT clad thickness (a).

A. CLADT =

67. FUELO fuel pellet outer diameter (m).

A. FUELO =

68. FUELI fuel pellet center hole diameter (m).

j A. FUELI = i

69. PICHI fuel pin unit cell PICHI (m).

This unit A. PICHI = cell dimension is used by the fuel pin model to determine effective flow area per pin, rather than the flow areas entered in the mesh cell dimension section below. f

70. CRNDD0 control rod diameter (m).

) A. CRNDD0 = )

71. POISD0 i

120

poison rod diameter (m). A. POISD0 =

72. BUAVG A.

BUAVG = average Burnup (MWD /HT).

73. THEL 1 l

A. THEL 1 = Zr-rich eutectic melting point 3200K). ) 74 THEL 2 A. THEL 2 = Eutectic melting point (2200K).

75. THEL 3 A.

THEL 3 = Zr-lean eutectic melting point (2200K).

76. TMEL4 A.

TMEL4 = Zr0 281 ting point (2650K). 2

77. BETA 1 disintegration rate time constant A.

BETA 1 = (0.1 sec'l).

78. FFAIL A.

FFAIL = fraction of rods that fail (1.0).

79. VFMOLT A.

VFMOLT = fraction of fuel that is molten.

80. RATOM A.

RATOM = oxygen to metal ratio.

81. TFMELT A.

TFMELT = fuel melting point (3000 K).

82. DELTSL solidus to liquidus AT.

A. DELTSL =

83. PUCONC Pu concentration.

A. PbCONC =

84. FRACTD fraction theoretical density of A.

FRACTD = fuel (0.95). 85, CPMELT (500 J/kg/K). specific heat of molten UO2 A. CPMELT =

86. RUO20 A.

RUO20 = room temp density of UO2 (10500 kg/m**3). ~ 1 A.3 3 2. Rod Aodel Array Input + 1 NOTE: NOAXAL = NCORE2 - NCORE1

87. (TSS(I), I=1,NOAXAL) (NOAXAL values)

A. TSS(I) = mverage steady state axial temperature (1500 K).

88. ((PSHP(I,J), I-1,NOAXAL), J-1,NRING)

A. PSHP(I,J)= reactor power shape for cell I and ring J.

89. FFLUSC effective fast fluence for strength A.

FFLUSC = coefficient (1.5e22)

90. FFLUSH effective fast fluence for strain

) A. FFLUSH = hardening (1.5e22) I I

91. (CWSC(I), I-1,NOAXAL) (NOAXAL values)

~ 121 i

A. CWSC(I) = _ effective-cold work for strength coefficient in cell I.

92. (CWSH(I), I=1,NOAXAL) (NOAXAL: values)

A. CWSH(I) = effective cold work for strain hardening. '93. (((TFUI(K.J), K=1,NOFUEL), (TCLI(K,J), K-1,NOCLAD)), J-1,NOAXAL) A. -TFUI(K,J)= initial fuel temperatures (K). B. TCLI(X,J)- initial clad temperatures (K).

94. NBALL' A.

NBALL = number of axial' elements around burst element that balloon (1).

95. (ZRCON(J), J-1,NRING)

A. ZRCON(J)= mass of zire associated with control rods in ring J (kg). 96.~(ABCON(J), Ja1,NRING) A. ABCON(J)- mass of absorber material associated with ^ control rods in ring J (kg). 97 (SSCON(J),-J-1,NRING) A. $$ CON (J)= mass of stainless steel associated with control rods _in ring J-(kg).

98. (ZRPOIS(J), J-1,NRING)

A. ZRPOIS = mass of zire associated with poison rods-inLring J (kg).

99. (POPOIS(J), J-1,NRING)

= mass of poison material associated with i A. POPOIS I poison rods in ring J (kg). l A.3.4 Structures Input 100. MAXMOD,STCNVG SFCNVG.ITCRST,ITMELT.IFAILS,NMATMX A. MAXMOD = maximum number of unique structure models to be input. freezing and melting convergence B. STCHVG = ceiterion. C. STCNVG = crust freezing and melting convergence 4 criterion. D. ITCRST' = maximum iterations on crust. E. ITMELT = maximum iterations on structure melting. structure mechanics calculation switch F. IFAILS = (0=off, 1=on). G. NMATMX = maximum number of materials for which mechanical properties are input on file MIMMCH. NMATMX must be at least 1. Repeat the following card for each cell 101. (IDSTRC(I) I-1,5) A. IDSTRC(I)= structure identification number for Ith structure in cell. Enter 0 (zero) for no structure 7 I Structure _ model input must be present for each unique structure entered on the cell card. The model data is t ; l input in the order given on the cell card. Model data l 122 j l

f l t for a structure previously entered does not need to be repeated. 102. IDSTRC, ISTYPE, IHTTYP, MCHMAT, MCHMOD, IBNDOP, IPOWOP, NREGNS, MAXNOD, IDCEL, IDSUP, STTEMP A. IDSTRC = structure model identification number as entered on cell card. B. ISTYPE = structure type identification number. a. ISTYPE = 1-9 means core structures (pins, etc). 10-29 mean's wall type structures (core b. ISTYPE = l-barrel, vessel wall). c. ISTYPE = 30-49 means plate type structures. d. ISTYPE = 50669 means special structures (head I bolts, columns, etc). i C. IHTTYP = heat transfer model switch a. IHTTYP = 0 means no heat transfer calculation. 1 means 1D finite difference. b. IHTTYP = c. IHTTYP = 2 means 2D finite difference. mechanics material property identification D. MCHMAT = number. a. MCHMAT = 0 will default to 1. b. MCHMAT > 0 will use the mechanical properties for material MCHMAT. [ E. MCHMOD = mechanical model selection switch. a. MCHMOD = 0 means no mechanical calculation done for this structure. b. MCHMOD = 1 means grid plate cell model. 2 means elliptic perforated plate c. MCHMOD = model. d. MCHMOD = 3 means unpressurized cylinder model. e. MCHMOD = 4 means pressurized cylinder model. f. MCHMOD = 5 means eggerate plate cell model. l + F. IBNDOP = heat transfer boundary. condition flag. a, IBNDOP = 0 means all surfaces see cell fluid. 1 means inner surface sees cell fluid b. IBNDOP = and outer surface is adiabatic. j c. IBNDOP = 2 means inner surface sees cell fluid i and outer surface sees fixed temperature through l constant heat transfer coefficient. G. IPOWOP = power option flag. a. IPOWOP = 0 means no internal power generation l b. IPOWOP = 1 means uniform power generation per unit mass. l H. NREGNS = number of structure material regions to be input. ( I. MAXNOD = maximum number of temperature nodes to be divided among all material regions. uni qu e identification number for J. IDCEL = structure. This number is used on output and internally to identify a structure for mechanical l calculations, e g., IDCEL = 1 may refer to the bottom 123 ) l 1 \\

head. These numbers should start with 1 and increase sequentially. identification number for the structure K. IDSUP = supporting structure IDCEL. An identification number >= 900 indicates an external support outside the problem domain: a. 900 means support at outer top corner of cell. b. 901 means support at outer bottom corner. c. 902 means support at inner top corner. d. 903 means support at inner bottom corner. initial structure temperature (K). L. STTEMP = .The following card is input only if IBNDOP = 2. 103. STBND, STHTCF A. STBND = boundary sink temperature (K). boundary heat transfer coefficient B. STHTCF = (W/m**2/K).. The following card is entered for each material region (NREONS times). "XXXX" stands for input that varies depending on the type of structure. 104. MATID, NNODES, XXXX region material thermal property A. MATID = identification number. number of finite difference nodes in B. NNODES = region. C. XXXX = geometry input varies depending on type of structure, s. For ISTYPE = 10-29, wall data. RADIN, RADOUT RADIN = inside radius of region (m). RADOUT = outside radius of region (m). l b. For ISTYPE = 30-49, plate data. The following is for MCHMOD = 1

RADIN, RADOUT,

THICK, HOLED,

HPTCH, COLMOD, COLPTCH RADIN = inside radius of plate region in cell (m).

RADOUT = outside radius of plate region in cell (m). THICK = plate thickness (m). HOLED = hole diameter for perforated plate (m). HPTCH = hole pitch (m). COLMOD = column outside diameter for plates supported on columns (m). COLPTCH = column pitch (m). The following is for MCHMOD = 2 RADIN, RADOUT, THICK, HOLED, HPTCH, RMINOR inside radius of plate region in cell RADIN = (m). l i f 124

-RADOUT = outside radius of plate region in cell (m). THICK = plate thickness (m). HOLED = hole diameter for perforated plate (m). HPTCH = hole pitch (m). RMINOR = minor axis length for entire plate (m). c. For special structures. ISTYPE = 50 column data COLMOR. COLMIR, XLNGTH, COLMP 'COLMOR = column outer radius (m). COLMIR = column inner radius (m). XLNGTH = actual length of column-(m). COLMP = column pitch (m). d. ISTYPE = 51 wall-plate joint model RADIN, RADOUT, THICK RADIN = inside wall radiue (m). i RADOUT = outside wall radius (a). THICK = plate thickness at j oint (m). e. ISTYPE = 52 Bolt model BOLTOD, XLNGTH. THRDTP,.CRCLOD, SPACE BOLTOD = bolt shank diameter (m). XLNGTH = bolt shank length (m). THRDTP = thread type. diameter of circle on which bolts are CRCLOD = placed (m). SP ACE = bolt spacing along circle perimeter (m). f. ISTYPE = 53 general structure (heat transfer only) ATOT, THICK 2 ATOT = total surface area (m ) THICK = thickness (m) A.3 5 Radiative Heat Transfer Input 105. NGROUP A. NGROUP = number of radiation groups (1). 106. ITRMAX A. ITRMAX = max number of iterations (50). 107. RCONV A. RCONV = radiation convergence criterion (1.0E-4). 108. (EMISS(K), K-1,6) A. EMISS(K)= emissivity of structure K. The structures ares pins, outer wall, inner wall, bottom plate, top f plate, and debris. 109. (ARH0(K), K-1,3) l A. ARH0(K) = mass absorption coefficient for field K. r. L< \\ I l 125 1 1

i. l A.3 6 -Debris Input 1 110. LDOFF A. LDOFF = debris bed on/off (0/1) switch. 111. NZMAX, NBEDM, ICOND A. NZMAX = maximum number of nodes in a bed. B. NBEDM = maximum number of beds. C. ICOND = effective conductivity model indicator. a.- ICOND=1 means use. Imura'Takegoshi/Vortseyer model b. ICOND=2 seans use W111hite-Kunii-Smith /Luikov model 112. PORMAX, DZMIN, DEFF, SO A. PORMAX = maximum porosity in a bed (0.55). B. DZMIN = minimum mesh size (.05 m ). C. DEFF = effective particle diameter (0.005 m). D. SO = initial specific power (w/kg). l l i e ~ f 126 )

A.4 General Rules for Problem Restarts ) i Two data files are required to restart a MELPROG calculation, MIMDMP and MIMIN. The MIMDMP file is a binary i restart file written by MELPROG. When MELPROG is run in the I restart mode, file MIMDMP is used as the initial restart data l source and as the subsequent restart output file. Therefore, a copy of each restart file should be made to prevent destruction of the previous restart file. The user input data file for a restart is MIMIN (the same j file name as the problem initial input file) and the input j required is similar to that required for the initial problem output. Major differences are that very little PINS module data is input and no STRUCTURES module data is input. All of the ) general input, FLUIDS module input, RADIATION nodule

input, and DEBRIS module input are required.

Restart input data is of three types data that must be ) entered, data that at the user's option may be entered or entered i as 1.1E37, and data that should always be entered as 1.1E37. Entering a value of 1.1E37 tells MELPROG to use for this entry the value from the restart file, MIMDMP. The type of data,needed for each entry in the restart input file is indicated in the input instructions as " enter", as "either", or as "1.1E37" respectively for each type of data. A.S. Restart Input Cards A.5.1. General Input 1. ITITLE (20A4), (enter). A. ITITLE = problem title (up to 80 characters). 2. IRESTRT (enter). A. IRESTRT = the restart switch .a. IRESTRT = 0 means no restart. b. IRESTRT = N means read restart dump N from file (MIMDMP), c. IRESTRT = *1 neans read final restart dump from file (MIMDMP). 3 DMPINT (enter). A. DMPINT = restart dump interval (sec). NOTE: DMPINT is not input when HELPROG is linked to TRAC. The TRAC restart dump interval is used when linked to TRAC. 4 IPRNT, (enter). number of time intervals between A. IPRNT = printouts. 127

5. IPRTF, IPRTP, IPRTR, IPRTS, IPRTD, (enter). A. IPR 17 = fluids print flag (0=off, 1=on). B. IPRTP = pins print flag (same). C. IPRTH = radiation print flag (same). D. IPRTS = structures print flag (same). E. IPRTD = debris print flag ( sase ). 6. KMAX (enter). A. KMAX = number of axial mesh cells 7. NRING, (enter). A. NRING = number of radial rings 8. NCORE1 (enter). A. NCORE1 = fluids cell number for core bottom cell 9. NCORE2 (enter). A. NCORE2 = fluida cell number for core top cell

10. LRAD (enter).

A. LEAD = radiation stfitch (0=off, 1=on)

11. TSTART (either).

A. TSTART = the starting time (sec).

12. ENDTIM (enter).

A ENDTIM = problem end time (sec)

13. DELTAT (1.1E37).

A. DELTAT = the initial time step (sec). 14 NTIM (enter). A. NTIM number of time step pairs in time step = table (22).

15. (TIM (I), DTIM(I), I-1,NTIM) (enter).

A. TIM (I) = Ith time at which the following time step begins (sec). B. DTIM(I) Ith time step which begins at t.ime TIM (I) = (sec).

16. SS (enter).

A. SS pin. model steady state switch (0=0n, = 1=0ff). If this is-set on, MELPROG will attempt to reach an equilibrium stcte for temperature and fluid flow. NOTE: If the equilibrium t<asperature for the particular set of initial conditions is high enough to cause significant transient effects, e.g., clad oxidation, no equilibrium will be reached; the problem should be run as a. transient starting with the input as initial conditions.

17. QLIN (enter).

A. QLIN = total reactor power (w).

18. N00 (enter).

A. N0Q = number of pairs'in the power table.

19. ((TIM (K), POW (K)), Ka1,N0Q) (enter).

A. TIM (K) = time point K in power history table (sec). B. POW (K) = relative power at time K multiplies total reactor power. 128

i A.5.2. Flu 1> !nnut A. 5. 2.1. Fluids Scalar Input

20. ERROR 1 (either).

Error criterion 1 for Newton-Raphson 3 A. ERROR 1 = iteration. This is the relative convergence criterion on the pressure and void fractions.(0.001-0.1)

21. ERROR 2 (either).

Error criterion 2 for Newton-Raphson A. ERROR 2 = iteration. This is the absolute convergence criterion on temperature in degrees K. (0.1-10.K)

22. ERROR 4 (either).

time step exp2nsion and contraction A. ERROH4 = parameter. This is actually the percentage of the current time step by which the time step can be reduced or increased in one step.(0.0-99.0).

23. EPSA (either). maximum fractional c h a r.g e in alpha over A.

EPSA = time step (0.05-0.1).

24. EPST (either).

A. EPST = maximum fractional change in temperature over tire step (0.01'O.1). 25 IFLOOD (enter). A. IFLOOD = reflood switch a. IFLOOD = 0 means no reflood b. IFLOOD = 1 means reflood 26..ITERMIM (enter). A. ITERMIM = minimum number of Jacobi iterations (5).

27. ITERMAX (enter).

A. ITERMAX = maximum number of Jacobi iterations (20).

28. AL10 (either).

~ low alphs limit for vapor (IE 1E-3). A. AL10 = If the volume fraction of vapor is below this value the field is turned off.

29. AL20 (either).

low alpha limit for water (IE 1E-3). A. AL20 = If the volume fraction of water is below this value the field is turned off.

30. AL30 (either).

low alpha limit for corium (1E'6 6 1E-3). A. AL30 = If the volume fraction of corium is below this value the field is turned off.

31. ETA (either).

A. ETA = not currently used but still read

32. BETA (either). Courant time step factor (2.0-10.0).

This A. BETA = is the factor by which MELPROG is allowed to exceed the Courant condition. L

33. IFILL (enter).

inlet boundary condition flag. I A. IFILL = 129

O I a. IFILL'- O means loop boundary condition is to be specified. I f b. IFILL = 1 means velccity boundary condition is .to be specified. A.5.2.2. Fluid Dynamic Array Input 34, (DELTAZ(I) I-1,KMAX+2) (1.1E37) A. DELTAZ(I)= mesh cell length for cell I (m).

35. ((FRCTN(I,J),I-1,KMAX+2),J-1,3) (either).

A. FRCTN(I,J)= additive friction factor for cell I and field J. j

36. (ABYP(I),I-1,KMAX$2) (1.1E37).

j A. ABYP(I) = bypsse flow aist. far/ cell I (m'#2). This must be greateiCthAn.tero in cella having no pins. Note that settirg.the bypaAa' flow to zero does not shut off heat tranhfer from the' radial wall to the bypass channel. 37 (NPII:S(J),J-1,NRMG). 01.1J 3 D - A. NPINS(J)- number. or f'act toits jn c60tal ring J.. 3

38. (NCNROD(J),J-1,NRINC) (1,1E3'7).

J A. NCNROD(J)= number' or, con'.rof rod & in radial ring J. i

39. ( N P OIS N ( J ), J-1, N RING ) '( 1. 0231 )..

A. NPOISN(J). number of' 3613en rods in radial ring J. i j

40. ( P A CDE N ( I), I = 1, KM Ay (2),( 1.1 EJ 7 )'w A.

PACDEN(I)- max'imum packinr,, fraction for corium 1 debris in cell I !4.44)w q

41. (RI(J) J-1,NRING),

( 3 0 ( J ), Ja 1, N T.1N C ), (2*NRING values) 1 (1.1E37) inner radius of ring J (m). .A. RI(J) = outer radius of rinccJ (t). B. R0(J) = A.5.2 3 Fluids Initial Conditions i I A.S.2 3 1. Field One - Steam and Mydr' ogen

42. (VG(I),I=1,KMAX+1), (KM AX+1 values ) (1.1E37).

A. VG(I) initial vapor' velocity for cell I starting = at cell I-1 (m/s). I

43. (TG(I).I-2,KMAX+1), (Khl+ X values ) (1.1E37).

initial vapor. temperature for cell I A. TG(I) = starting at I-2 (K).

44. (ALFG(I),1=2,KMAX+1), (KMAX values) (1.1E37).

initial vapor void fraction for cell I A. ALFG(I) = starting at I-2. ] i 8 a 130 l 1 l

h A.S.2 3 2. Field Two'" Water l

45. (VL(I),I-1,KMAX+1), (KMAX+1 values), (1.1E37).

A. VL(I) = initial water velocity for cell I starting at cell 1 (m/s).

46. (TL(I),I-2,KMAX+1), (KMAX values) (1.1E37).

l A. TL(I) initial water temperature for cell I = l starting at cell 2 (K). ~ l

47. (ALFL(I),I=2,KMAX+1), (KMAX values) (1.1E37).

l A. ALFL(I) initial water void fraction for cell I = starting at cell 2. A.L.2.3 3. Field Three - Corium

48. (VC(I),I=1,KMAX+1), (KM AX+1 values ) (1.1E37).

A. VC(I) = initial corium velocity in cell I starting in cell 1 (m/s).

49. (TC(I).I=2,KMAX+1), (KMAX values) (1.1E37).

A. TC(I) initial 'corium temperature in cell I = starting at cell 2 (K).

50. (ALFC(I) I=2,KMAX+1), (KMAX values) (1.1E37).

A. ALFC(I) initial corium void fraction in cell I = starting at cell 2.

51. (P(I) I=2,KMAX+1), (KMAX values) (1.1E37).

initial pressures ( Pa ). A. P(I) = NOTE: it is best to start the problem with axially uniform initial velocities if generating a steady-state s ol uti on. A.5.2.4. Fluid System Boundary Conditions

52. NVTAB (enter).

number of inlet velocity boundary I A. NVTAB = condition table pairs.

53. ((TIM (K.J), VB(K.J)), K-1,NVTAB), J-1,NFIELD) (enter).

A. TIM (K,J)= time point K for field J in velocity table (sec). B. VD(K,J) inlet velocity at time point K for field J = (m/s). NOTE: To avoid evaporation / condensation problems when either the steam or water void fractions are very small, velocity of t(he field with the small void set the inlet small value 1e-10). action to a 54, crYTTAB (enter). A. NTTAB = number of inlet temperature table pairs.

55. ((TIM (K,J). TB(K.J)), K-1,NTTAB), J-1,NFIELD) (enter).

'o A. TIM (K.J)= time point K for field J in temperature table (sec). l inlet temperatvee at time K for field J l B. TB(K,J) = (K). l 131

L. 56.-NATAB (enter). A. NATAB = number of inlet alpha table pairs. 57..((TIM (K,J), ALF(K,J)), K-1,NATAB), J-1,NFIELD) (enter). A. TIM (K.J). time point K for field J in alpha table (sec). B. ALF(K,J)= inlet void fraction at time K for field J.

58. NPTAB (enter).

A. NPTAB = number of pressure drop table pairs.

59. ((TIM (K), PH(K)), K-1,NPTAB) (enter).

A. TIM (K) = time point K for pressure drop table .(sec). pressure drop at time point K (Pa). This B. PH(K) = pressure drop is applied between the inlet and outlet' boundaries of the problem if the inlet loop boundary condition is set (IFILL=0).

60. NPBTAB (enter).

number of outlet pressure boundary A. NPBTAB = condition pairs.

61. (TIM (K). PB(K)), K-1,NPBTAB) (enter).

time point K for pressure boundary table A. TIM (K) = (sec). B. PB(K) = outlet pressure at time K on boundary (Pa). A.5.3 Rod Model Input A.5 3 1. Rod Model Scalar Input

62. NOFUEL (enter).

~ Note: NOFUEL must be the same as in the original input. A. NOFUEL = number of fuel radial nodes.

63. NOCLAD (enter).

J Note: NOCLAD must be the same as the original input. j A. NOCLAD = number of clad radial nodes.

64. TLM (either).

A. TLM = Temperature error criterion (10.0K).

65. CLADO (1.1E37).

A. CLADO = outer clad diameter (m).

66. CLADT (1.1E37).

A. CLADT = clad thickness (m).

67. FUELO (1.1 E37).

A. FUELO = fuel pellet outer diameter (m).

68. FUELI (1.1 E37).

A. FUELI = fuel pellet center hole diameter (m).

69. PICHI (1.1 E37).

fuel pin unit cell PICHI (m). This unit A. PICHI = cell dimension is used by the fuel pin model to I 132 l l

determine effective flow area per pin, rather than the flow areas entered in the mesh cell dimension section below.

70. CRNDDO (1.1E37).

A. CRNDD0 = control rod diameter (m).

71. POISD0 (1.1E37).

A. .POISDO = poison rod diameter (m).

72. BUAVG (1.1E37).

A. El' AV O = average Burnup (MWD /MT).

73. THEL 1 (either).

A. THEL 1 = Zr' rich eutectic melting point (2200K).

74. THEL 2 (either).

A. TMEL2 = Eutectic melting point (2200K).

75. THEL 3 (either).

A. THEL 3 = Zrk1'ean eutectic melting point (2200K).

76. THEL 4 (either).

A. 'TMEL4 = Zr02 melting point (2650K).

77. BETA 1 (1.1E37).

disintegration rate time constant A. BETA 1 =*l). (0.1 sec

78. FFAIL (1.1E37).

A. FFAIL = fraction of rods that fail (1.0). A.5.4. Structures Input

• No Structures Input on Restarts A.5.5.

Radiative Heat Transfer Input

79. NGROUP (enter).

A. NGROUP = number of radiation groups (1). ~

80. ITRMAX (enter).

A. ITRMAX = max number of iterations (50).

81. RCONV (either).

A. RCONV = radiation convergence criterion (1.0E*4).

82. (EMISS(K), K-1,6) (either).

A. EMISS(K)= emissivity of structure K. The structures ares pins, outer wall, inner wall, bottom plate, top plate, and debris.

83. (ARH0(K), K-1,3) (either).

A. ARH0(K) = mass absorption coef ficient for field K. A.5.6. Debris Input

84. LDOFF (enter).

A. LDOFF = debris bed on/off (0/1) switch.

85. NZMAX, NBEDM, ICOND (enter).

A. FZMAX = maximum number of nodes in a bed. j B. NBEDM = maximum number of beds. C. ICOND = effective conductivity model indicator. 133 i U_

L a. .ICOND=1 means use Imura"Takegoshi/Vortmeyer model b. ICOND=2 means use Willhite8Kun114 Smith /Luikov model

86. PORMAX,.DZMIN, DEFF, SO (enter).

A. PORMAX = maximum porosity in a bed (0.55). B. DZMIN = minimum mesh. size (.05 m). C. DEFF = effective particle diameter-(0.005 m). D. SO = initial specific power (w/kg). A.6. References A.1. F. A. Koontz, FUNRD A A Free Form Input Subroutine, Report No'. 087260, Tenn. Valley Auth., Knoxville, TW (1975). e 134

r E Appendix B. Output Description B.1. ' Introduction .In this a ppe nd i x', a description of the MELPROG output is ~given. .The information given herein briefly describes the MELPROG formatted output file (called MIMOUT). which contains.a The printed history of the more'important calculational results. user should be cautioned that the. output file undergoes constant evolution and will change with time. The ~ sample output given at the end of this appendix will be helpful to-the user in understanding this description. e 'As mentioned in the' input description, the input. processor uses a free-form input. routine. The input comments'are echoed to the output file, and thus appear as the first items in the MIMOUT file. After these' echoed comments the results are written out in a sequence of blocks in the following order: 1. FLUIDS printout block, 2. PINS printout block, 3 STRUCTURES printout block, 4. DEBRIS printout block (if a debris bed exists), and 5. RADIATION printout block. All the results printed out in these blocks are in SI units unless otherwise explicitly indicated. This sequence of blocks is only written to the output file every "IPRNT" time steps, where IPRNT is an input parameter. Between these output sequences, convergence information is written out after each time step. Contained in this single line of print is the time step number, the current time in seconds, the time step, the percent pressure error, the number of Newton iterations needed for convergence, the maximum percentage change in volume fraction over the step, and the maximum percentage change in temperature over the step. This line of information is also written out to the user's terminal if MELPROG is being executed interactively. j 1 135 I

B.2. FLUIDS Module Printout Block As indicated above the first block of information written to MIMOUT is the FLUIDS module print. This block contains temperature,

pressure, volume fraction - velocity, macroscopic density, macroscopic energy density, and microscopic density for each of the fields (1-vapor field, 2= water field, 3-corium field, 4-fuel rods, 5-strue'ure).

Also included in this field print is the hydrogen part'a1

pressure, hydrogen

density, saturation temperature, steam' water mass transfer rate (in the units kg/s-m'),

Zr macroscopic density in the corium field, and the corium mass transfer rate (in the units kg/s-m'). The FLUIDS block continues with friction coefficients (in the units kg/m") between the various fields and structures, heat transfer coefficients (in the units W/K-m 8) between the various fields and structures, and heat transfer rates (in the units W/m') between the various fields and structures. The last item in the FLUIDS block is the geometry information interfacial areas between the various fields and structures, the average diameter of the corium particles, and the fluid flow areas. All of the above FLUID block results are printed out by axial level where level 1 is a boundary condition, and levels 2 through KMAX+1 represent real fluid axial cells. B.3 PINS Module Printout Block The PINS module output (i.e. fuel rod output) appears next in MIMOUT and' is divided into radial-ring sections (NRING sections). The information printed in each of these includes the f uel rod temperature distribution, the fluids flow area through a given ring, diametral gaps for the fuel rod, gap conductances (in the u n i t.s W/m*K), fission gas release (in the units kg-mole), cladding circumferential stress, cladding oxide layer thickness, cladding oxygen fraction in the S'Zr region, 88Zr thickness in the

cladding, linear power due to Zr oxidation, hydrogen generation (in the units kg/m-s), radially pin-averaged fuel temperatures, total hydrogen mass generated by all fuel rods in the core, and the rod power density distribution (in the units W/m').

All of the above fuel rod information is printed out by core axial level except for the hydrogen mass. B.4. STRUCTURES Printout Block The STRUCTURES module output block appears next in MIMOUT. This block contains geometry and temperature distribution information about each structure type that exists in each 2D mesh l cell. Heat

transfer, boundary condition, and material composition information is also printed for each structure type.

The reader is referred to the input description for an explanation of the available options associated with these items. ) 136

f ' At the end of the structures block, a summary of the structures is printed containing location and description of the different structure types, failure indicators, average temperatures, loads, and stresses. B.S. DEBRIS Module Printout Block The DEBRIS module block appears next in MIMOUT provided-at least one debris bed exists during the current time step. The output consists of a section of print for each existing bed. If no beds exist during the current time step, then no debris bed print will appear in the output. The printout for a particular bed contains typical debris bed parameters, bed location and

size, temperature distribution, volume and solid f ra ct i on' distributions of various materials, steam fraction in the bed, and volumetric heat generation due to decay heat and Zr oxidation.

B.6. RADIATION Module Printout Block The RADIATION module block is the last in the sequence, and is divided into two sections. The first section is the radiation volumetric heat source deposition (in the units W/m') into each of the three fluid fields. These heat sources are listed by axial level for each radial ring (i.e., levels first and then rings). The second section lists the radiation heat fluxes (in the units W/m ) to the various structure surfaces and fuel rods. 8 B.7. Sample Output ~ l i 137

Sample Convergence Information e s,t. u.E m eat n.ece sine ep.0

aa e 4.006 t.03 1. 3.f f t.Os 4 3.9 0

.t.09 30390 3.M 68 03 S.39 64 02 4.643.S.t.0S t.03 1.3 2H.04 303SS 3. M.6t *03 S.29a 3.4 *038 03 1. 64 30399 3.M 3 g.03 1,30 4g.0a it. .64 t.0g e a.C.$44.e3 1.3 6.et Os 30393 3.M.88.03 S.39403 S,39391 02 1.10394 09 e..O 34 0a 4 4.0 0.t.03 T.30 e 4.066 303SS 3.M 48 03 S.39344 02. 61 6t.0S 20364 3.M t.03 T.309%E.Os i 9.66628 09 30794 3.M.9t*03 S t93 t.03 l e e..M 4 03 1.20068 0a 4 4.t.31 91 03 S.29328 03 t.41814 05 t.03 1.301 4 04 30391 3.M20t *03 S.39298 02 1.64 74 09 t.03 T.301 4 04 303S9 3.M3 9t*03 S.292St.02 9.4.1994 09t.61944.09 e 4.18.34 83 T.300St.04 e 4.tSt.9t.e3 ?.306M 04 30368 3.M304 03 S.392?t.03 f.6 02t. 4 4.30 30MO 3.M2.t.03 S.392N.03 3t.05 e e. 38 ts.t.es 1. 30a t t.Os 09 20M t 3.Ma tt.03 S.3920t *02 9..e79.9t 09 e e.350.t.es 1.303.f.0s 78 03 41 t.93 ?.309 4 30M3 3.M228 03 $,39.St.02 1.61078 05 e 4 213 04 30363 3.M2J4 03 S.39.H.02 30364 3.M2M *03 S.39 0.818 t*06 e e.796st.83 f.302M.Os Sample FLUIDS Printout Block ..................e n m..... .n. u..... .....n... .....u ...n......u... .u...... m. o.o... .....n. n .. n. .u...... ..n.... mu. .u. o m ....c on ..m.,.. . u... n. .u ,u.. .m.. . n.. -.. u.. ...........u.... .m....n ........................m.... ...o..u..... u.m... ...c..... ......... m ...u... ...u.. .m

c. u.

......m. .nu... m. . n.m. u. c. m. ..m..m... ..m...n . u.u..c .n ...u... .n m..... .. n..u... .u .m.. m,.. m.. m.n.. ..............u.n .....u.... .a..a......u...... ..mn c -.. m. ...u ...u... n . m m.c ...c u.. ...u.... .c

n...

.u. .m .um ... u n.m. ..m,..u.. u.u...c...m...n..

m..u.....

u m, c.. ..n...u.... .mn...c.. .m ..u...... .c. u .........m................. .m...........m. ....m. .m.......

m...

...u..c..... j . n.. n. i i ..m... i u.. un.. .. m u.. 1 i ...m.o.n ....m... n ...n.u.... ..u.u. .. m.n...... .u .u........ ...............................u........................n................

a..a..

.. n . co u. .. mu . m.. n. 1 .-......... m.n. ...m.., ..un.. .u.... ..m n . mu ..............o... . a.. ..n...... ..m .m. ..n .m .. m. .m ..um n. .... mn. n. . u.. n. . u... l ..n.. I ..u ...m. l l 138 l I"'."^-- -______m_

Samplo PINS Printout Block W18V1 80s GleG 3 E eget (La*Sto teawlltet flat

• 334 t.6e SE C

-= to salat ase ea01st tieplea?Uet 0151eleufl0N Ots. a elet.8 tiet.3 tie 6 9 4146.5 8184.0 t186.3 side.e 41a3.3 0143,3 v143.3 LDw sets ete 400

• o*.2 0

660188808 3 .309528*01 3 .et946t*00 4 .te264t*0t S.94tlet*01 sent taat see waLWil f ee taCm eal6L PIN Sitt3DN

• a 4

.040938 Os 48 CIndDUCleasCt Watutt fpe E ACN Aalek Sf C180N*W/5 06.s 4 .ptlett*08 Stet fl1$30N Dal 98Lt4580 Face ta(M aalaL $tC1l0N-s@g 9 401998=03 tap 0les CleCWM. $10885 foe LaCM salat SICilew. pe t O. SIM 44Vte 1MICENttl F0e f aCN ansat 5f Ctless

• u e

932848 03 % teactiew Of $N 3 2e toe taCM PIN antal SIC?t0N teet t O. 3 *t ect et Satyeafl0N IN Sita PMalt 888 ealet lfC1* PCT 40. l tie.Fe THICENill 80s teCM eal6L SICil0N

• m 9

.St?2M 03 Saute # rowte Out 10 le Calcat3DN Pts anf at SECT w/m t .52130t*02 USSOGEN e4Assp FUNCilON 908 (ACH 4:14L StCTION t.844415 02 Wer f essC180NS Fee 3.C5.se.st. 6ND TE

• st/m**3/ SIC S

0. O. O. O. O. 401a449' avreact0 tutt TEmptsatystl DtG. a 9.ttelet*04 19tet Of 4628*02 es m2 met ettN ettralt0 Out 10 78 cascat 0N slet ese ma0l&L Ponte Otaal Olitelpu110N

• w/m**3

.322ett*01 .32241t*01 g222a?t*01.222878*01.22341t*01.2224?t*01 C. .31289t*01 *.94249t*0t Sample STRUCTURES Printout Block Steuttuettl) IN W!tse d 484 6alat af get ? j SteuttWet tiet to tef at te*Nltre DetleN g 30uese? CosettitN 08159N 3 MCION samatta lhett$ lek 80tN1lf 3 Cat l0N t&AISta t S es00t 'ItMPt e a twat Souncast a

3. 49tM*00 3.994 tt *02

• 1.tetit*02 3.3stSt*00 8.0724t*02 3.sel2t*00 S.1049t*02

3. 320M *00 5.9344t*03 a

3.St0tt*00 9.3953t*02 3.a t tat *00 9.3 4138 *03 falLWet SUTPUT 3 d tot es00th samt tatt f La0 LIFT f eactieN svit fiesP 4040 Stells 9 4 1 e Sett Ple?4 9 9. 8.0088*02

• 3.193t*04 3.0884*00 e 3 1

e pett PLatt 0 0. S.00st*03 f.22ft*03 f.042t*0e 139

Samplo DEBRIS Printout Block 3 mm =004 se. e e t =c.amt e. tov sect =0, avt. eM0 Sat . Deva tavt . un. M D mWete. e M.ets s. 3.t t,e.a e attwt. .e 3 *= Da .0003 s ,0 tat t o,tetc uC3 90,ulee,... in.3nt *0. i sera 6 0.ct ..netit*0s, 3.M t. no.3n. t*03 se nOS 10546e.tn tw*03 se , t.5, a. 46 . 3.s 0 Se 50 '3e nell. 3.et??at*03 se 2003 mall. 6.3070H +03 se WO3 mall . 3.46491t*04 EG tuttCitt mall. 4.?geest*03 sg 90tuut feaC189NS SOLIO Feattlewt EWittitt feaCilowl east 0 Star IPS WO3 le 2001 tut. d3 fe fact tut. Se feet WO3 M30v f tem 43 .e tt 3304 .St. 0.000 0,000 0.000 0.000 9.800 0.000 0.000 4.000 .ett c 000 . tot .860 .seet et . Des 3018 .038 998 0.000 0.000 .003 9.800 0.000 0.000 0.000 .ett 0.000 409 .s60 . east to . set Slot. 0.000 .996 ' O.000 .000 0.o00 9.e00,0,000 9.000 0.000 .ett 0.000 406 .860 0.e000 g .eae 3404. 0.000 .996 0,000 004 0.000 9.800 0.000 0.000 0.000 .ett 0.000 909 .040 0.cooO O .ges 3900. .330 .He 0.000 .033 .000 .013 0.000 1903 0.000 .ett .000 .tol .060 . eses T .esa 3900. 0.000 014 0.000 .033 503 .918 0.000 . Set 0.000 .ett 0.000 .tM .860 0.eo00 .ese 3900. 0.000 .910 0,000 .030 0.000 434 0.000 .340 0.000 .999 .000 40% .660

0. coop a

g ,gg g,00, 0.000 .ete 0.000 00s .381 .998 0.000 .918 0.000 . set .000 004 .860 0.eupo a . set 3900. 0.000 004 0.000 .ett .090 . ate 0.000 .ete 0.000 ,etc .030 . toe .e40 0.e000 3 . 3m 3993. 0.000 .114 0.000 .030 .See S.000 0.000 0.000 0.000 .003 .0tt .404 .860 0.e000 3 430 3408 0.000 .399 0.000 .043 .ses e.000 0.000 1.000 0.000 .ett .086 .tol .860 0.0000 t .e=0 in3. 0.003 .330 .tu .00* .sce 9.e00 1.000

4. coo t.000

.en .Dit .i0s .no 0.eooO Sample RADIATION Printout Block f tute 840tallem estat lovett, s/e**3 +3.346est*00 +3.elltat*09 6. 3.83003t*03 9.39458t*04 c. 4.193ect*03 6. O.

• t.310938*06 0.

O. et.34304t*06 c. 9. 4.37313t*03 9. O. 4.943 Sit *04 0. O. S.87947t=0s f.StHet*00 C. l.91834t*03 3.06099t*04 4. i S.199334*04 0. O. J

• 3.t3644t*06 0.

9. l 30 0. O.

• 0. 0. t. 4.1 t
• 06 t.

0 C. O. 3.04337t*04 9. O. l 3.313a0t*03 3.08001E*03 0. J f.6896tt*03 3.St9Mt *04 c.

• T.e31 tot *04 9.

O.

• 3.35ttet*06 6.

O. 3.319 6 9t eM 0. 9.

• 3.019tet*03 9.
• 3.73 8 t M.05 6.86e?3t*04 9.

O. O. 9. O. O. O. C. O. O. O. i O. O. O, 0. O. O.* O. O. 6. O. O. O. Baeta130p setet FtW1tl TO SteWCfuMS. w/m**3 l

• t.34448t*06 *t.40139t*03 1
• S.4338tE*01
• t. M3eet *0 4 0,

as.93te28 0e t.30694t*01 d.Se390E*06 5.03 e34t *M

• t 0 tttet *M t.3933tt*06 *t.99900E*06
• t.2et60t*00
• t.60137E*03

+3.43333t*06

• t.M 36et*04
• e,03 709t *03 9.
• t.3 tHtt *06 0.6533et*06 S.tS t39t *M 4.193368*06
• S.501e3t*06

%.366tli*06 *9.990335*0%

• t.Settet*0t
• t.40?21E*03
• T. 3 9 3 3 el
• 00 4.e460e8*03 3.83333t*04 *e.3636t4*00 l
• 3.93324t*03 6.8138 t8 04 1.esteat*04 9.

9.d Mott*04 3.68939t*05 4.07334t*09 *0.83eett*03 S.See90t*06

• 4.S t3 338 +M
• e.0999tt*03 8.18446t*06

+3.093 tit *0e 8.38413t*06 at.0910tE*05 e 0.9t?ltt*03 *t.38393t*03 f.tt934t*03 *e.s07ett*03 S.00106t*03 at.3ttett*03 3.19391t*04 *t.33M64*03 9.90663t*04

• l.397944*03 3.e3646t*06 *t.30306t*03 3.446474*06 ol.3?t438*03 8.930138*04
• t.3tM1 fen 3 i

140 L___-_____

l Appendix C. Code Structure C.1.. Introduction E This appendix describes the functional layout and operation of each.of the. main MELPROG modules. For each module, the purpose and operation of the main subroutines are described. In addition to describing the modules, the TRAC-MELPROG link, the initialization section and the data storage system are also discussed. Figure C.1 illustrates the overall structure of MELPROG. The main MELPROG driver calls each of the module controlling subroutines. At the beginning of a problem, the initialization routine MPRGIN is called. Then, on each time step the other modules are called'in the following order: FLUIDS (FLDYN), PINS (PINZ), STRUCTURES (STRUCT). DEBRIS (DEBRIS), and RADIATION (RADMPG). The driver also checks for end of problem

time, performs graphical output by calling GRF every IPRNT time stops, and writes restart dumps at DMPINT intervals.

The input for a problem should be stored on a file named MIMIN. The output from a problem will appear on a file MIMOUT. If a restart is performed, the dump file should be named MIMDMP. C.2. FLUIDS Module This section describes the functional layout, purpose and operation of the FLUIDS module. The main FLUIDS module driver is FLDYN. This driver calls three o*her principal routines: POINT,

FLUIDS, and TSTEP.

POINT is the routine wherein the FLUIDS module pointers are set. TSTEP is responsible for time step Control. The FLUIDS subroutine is the run-time driver which calls the { operating routines and its functional arrangement is shown in Figure C.2. The first routine, FA, calculates the fluid flow area for each mesh cell. CDIAM then sets the corium j characteristic length, based on local hydrodynamic conditions and 1 previous history of the corium (e.g., if it were molten). CDIAM j also sets the maximum corium packing fraction for each cell. The next pair of routines, BC and

SETBCS, set boundary conditions for the FLUIDS module.

BC extracts the boundary conditions for the current time step from tables of the boundary velocities, void fractions, and temperatures versus time. SETBCS sets up the actual boundary conditions for the appropriate fluids variables. l I 141

4 TRAC LINK DRIVER 1r 1r 1r IP If RADIATION. PINS FLUIDS STRUCTURES DEBRIS 1r if VICTORIA IFCI 1r EJECT t l Figure C.1 Overall MELPROG Structure 1 l 142 l

R.UIDS rTA

ODIAM

--*BC

SETBCS

-* GAMMA ~ -* KNOB -*SETPROP -*PERIM llc -*MOPRED -*MOSTAB -*MOCOEF -*SETBCS -*SETPROP

EAT

-* GAMMA 3 --*JACCO --*NEWRHS - *JACSET -*MATSOL ~ - *CONVRG -* ADJUST -*SETBCS l -*SETPROP -* GAMMA

0

-*MSTAB

CSTAB

-*ZRHEAT - *CNTNUTY - *SOLNPRT 200NPRT Figure C.2 Functional Diagram of FLUIDS Module Routines 143

i OAMMA-is responsible for.. setting up.the various mass transf er ' terms, the principle ones being the' wa ter-steam phase change and the fuel / control roduto-corium mass transfer. KNOB is a control-routine for the fluid dynamics e qua t io n solver. It selects the~ fields in each mesh cell that are not present (void fraction below a small value, e.g...AL10) and also do not have a source term (mass transfer, convective transport) l that would increase the void fraction. These fields-have their l-corresponding mass, momentum and energy equations turned off within the affected cell. (LKNOBS is the mass and energy equation pointer, LVON is the pointer for the. momentum equation). SETPROP is then called to pick up the old properties and' I their derivatives with respect to pressure, void fraction, and temperature (IOLD is the pointer to the'old values). The next two routines, PERIM and FRIC, calculate the effective interfield. and ' field-to-structure perimeters and friction

factors, respectively.

The next

routine, MOPRED, is the first step for the Stability Enhancing Two-Step (SETS) solution t e chni qu e.

MOPRED uses old values in the convective terms of the momentum equations to decouple them spatially. These equations are then solved simultaneously on a celldby* cell basis, treating the interfield momentum. transfer-terms implicitly. This produces the " predictor" velocities. MOSTAB then uses the predicted velocities in the interfield terms' to decouple the momentum equations between fields. The equations are then solved on-a field-by= field basis, this time treating' the convective terms between cells implicitly. This produces the " stabilizer" velocities. MOCOEF.is used to solve the momentum equations for velocity as a function of pressure gradient explicitly. The stabilizer velocities are used in the interfield terms to decouple the fields. This step results in the velocities being expressed in the form vfl/2-f(Py*l,P[)) These expressions are then substituted.directly into the mass and energy equations for the pressure iteration. The pressure iteration begins by calling SETPROP, this time using new properties (pointer is INEW). HEAT calculates the heat transfer coefficients, GAMMA updates the mass transfers, and Q ~ calculates the heat transfer rates, based on the coefficients from HEAT. The actual Newton-Raphson pressure solution is performed by the routines JACCO,

NEWRHS, JACSET, and MATSOL.

JACCO calculates the Jacobian matrix coefficients, which are 144

L derivatives of the field mass and energy equations with respect to the independent variables (pressure, void fractions, and temperatures). NEWRHS sets up the right* hand-side of the matrix equation, which consists of the field mass and energy equation 1 residuals evaluated with the values of the independent variables at the current iteration. JACSET sets up the matrix equation for the cell-wise decomposition step for each cell. MATSOL then partially solves the matrix equations for each cell. The resulta j l.c' of this step are used to set up a tridiagonal matrix involvitg I only the pressures as independent variables, which is solved for the pressures. The cell pressures are then substituted back into the cell matrix equations to find the other independent i variables, and into the velocity equation from MOCOEF to find the new velocities. CONVRG checks for convergence of the pressure iteration by comparing the changes in the independent variables

• l to the error bounds (input).

The routine returns a value in j i ICONV of either 1 (continue) or 2 (converged). The old solution guess is then updated in ADJUST. ADJUST also checks for various pathological results, such as negative void fractions or zero temperature. If one of these is detected, ICONV is set to 3.- l At .this

point, FLUIDS either continues the pressure iteration by looping back to JACCO (ICONV = 1), exits back to the main driver program to have the time step reduced i n' TSTEP (ICONV - 3), or exits the' pressure iteration loop (ICONV = 2 ).

In the latter case, SETBCS, SETPROP, GAMMA, and Q are then called to update boundary conditions, fluids properties and transfer terms to the new time level, using the new fluids values. MSTAB and ESTAB are then called to perform the final stabilizing solution of the mass and energy equations, respectively. This ends the SETS solution of the fluids equations. ZRHEAT is used to calculate the oxidation of the Zircaloy component of the corium field. CNTNUTY solves the mass transport ^ equations for the various components of the fluid fields, such as Zr, 2r0 steel, etc., using the new fluid velocities and mixture r 2' densities. l The final two routines, SOLNPRT and CONPRT, output FLUIDS l module quantities on the print file whenever the print flag is set (every IPRNT time steps, or the end of the problem time). C.3 PINS Module This section describes the functional layout, purpose, and operation of the PINS module. The driver for this module is PINZ. This routine calls the main fuel and control rod ) subroutines.- The operational routines which are called by PINZ are DEBSRC, INTACT, and CONROD. 1 l s q 1 145

If there are disintegrated fuelt rod sections present (NODEBR or NOMELT ) 0), then DEBSRC is called to ' calculate the-mass transfer. rate to the corium~ field. If intact fuel rod sections are present' (NOINTC > 0), -then INTACT is called to perform the intact fuel rod. calculation. The final routine, CONROD, checks the control rod failure criterion to see if any' control rods _have

failed, If failure occurs, then CONROD ' calculates the mass transfer rate of control rod material to the corium field.

The most~important subroutine in the PINS module'is INTACT. The functional layout of this routine, which consists of the intact fuel rod model, is shown in Figure C.3.- First, radially averaged fuel and clad temperatures, used by several of-the models, are calculated in THINTG. Thermal saterial properties for fuel, clad, and eutectic are found in routines FUELPRP, CLADPHP, and UTKPRP. EXPAN is then called to determine fuel and l clad thermal expansion,- and the gap width considering' only thermal expansion and fission product swelling. COBILD calculates the cladding oxidation and oxygen l l absorption by'the clad 8-phase. The associated routines OXSAVE and OXSET.are usec to set up the parameters used in the COBILD model. l FOASRL. finds the fission gas release from the fuel matrix during ' the problem (not called if no burnup). A check is then made in INTACT to see if the cladding has failed. If so, the pin i internal pressure is set to the system pressure. PRSSTR uses the internal and system pressures to adjust the clad elastic strain. l YNGMOD returns the Young's Modulus for the clad as a function of temperature, fluence, and oxygen concentration. The strain is then used to update the fuel-clad gap and channel flow area. In the case of fuel-clad interference (i.e., the calculated gap is negative), CONMAIN makes a correction to the fuel and clad strains so that the gap is zero. GAPCON is then called to get the effective gap thermal conductivity in each axial section and GASCON is called to provide the fill gas thermal conductivity. QCALC determines the volumetric heat generation rate due to fission product decay. TMPREP then prepares the fuel rod variables to pass to TMCALC, which performs the solution of the radial heat conduction equation in the fuel rods for each axial 1 mesh. At this point, a comparison is done in COMPARE between l extrapolated temperatures (from THINTG) and new temperatures from TMPREP. If the difference is too great or if a clad radial node has just melted in the current time step (checked by CMELT, UTKTK and MELCHK), a loop is made back to THINTG; otherwise, the calculation is continued by calling WIPE. 1 146 lL_ _ _ __ _

INTACT

THINTG l

--+FUELPRP

OLADPRP

--*UTKPRP --*EXPAN ) --*OXSAVE OOBILD 1 - *OXSET --*fGASRL --+PRSSTR

YNGMOD l

--*CONMAIN f i --+GAPCON -GASCON --*OCALC -*TMPREP

TWCALC l

--* COMPARE OKLT --*UfKTK -{ --+MELCHK I f 0URST--*CKMN --+ WIPE M:BRTTL --*H2 GEN

!" POUT

-*CONRITE -*DISINT

DISSET Figure C.3 Functional Diagram of Intact Fuel Rod Model I

147 )

WIPE sets some flags if cladding has failed due to ballooning or one of the fuel rod failure criteria during the previous time step. It also calls BURST, which calculates clad ballooning and rupture, and CBRTTL, which determines the onset of clad embrittlement by oxygen uptake. CKMN provides parameters f or the clad e quation-of dstate to BURST. H2 GEN calculates the hydrogen mass transfer rate to the fluids vapor field and integrates the rate to get the total hydrogen production. FPOUT finds the fission product release rates to the vapor field. If this is a print step, CONRITE is called to output fuel rod variables. 1 DISINT contains the various fuel rod disintegration criteria. A flag is set for each axial section that satisfies one or more of these criteria. DISSET shuffles the data arrays used by INTACT. { C.4. STRUCTURES Module This section discusses the STRUCTURES module which controls the thermal and mechanical calculations for structures. The routines in this modale are divided into two functional sets: heat transfer routines and mechanics routines. Each of these sets is discussed separately in the following sections. I Subroutine STRUCT is the driver routine for the STRUCTURES module. STRUCT is divided into five parts. Part 1 performs data transfer from other modules to the STRUCTURES module. Part 2 consists of calling subroutine STRCHT to perform structure heat transfer calculations. Part 3 performs the structure mechanical calculations. Part 4 calls subroutines STGEOM and STRSET to calculate changes in structure geometry and fluid flow areas and hydraulic diameters, respectively. Part 5 calls subroutine OUTSTR to print structures information. C.4.1. Heat Transfer Routines This section describes the routines which perform the structure thermal analysis. The calling sequence for structure heat transfer calculations is shown in Figure C.4. Subroutine STRCHT performs the heat transfer solution for each structure. This routine calls subroutine STRPRP to determine thermal material properties and calls subroutine BCSSET for boundary conditions on each structure surface. STRCHT sets up the heat transfer matrix which 18 then solved by subroutine T R Y S LV. After solving for the temperature distribution, it l 148 4

) .I I i 1 I i STRUCT

STRCHT

STRPRP--*MSTRCT I

-*BCSSET l -+TRYSLV I - *REMESH j -*STGEOM

STRPRP--*MSTRCT

) !l i 1 l Figure C.4 Functional Diagram of St,ucture Heat Transfer Routines i l i STRUCT

SMOVE

-*DEBID 1 1-* LOAD l -

• STRESS MEC::

. MCHPRP MPFAIL -+GPLATE --*ELUPS --*CYUN --+DYNAM NEWSUP -+ COLUMN M:NGLOD -+LM P --+STRSET Figure C.5 Functional Diagram of Structure Mechanics Routines j I i 149 i

checks for melting and performs an iteration on the specific heat if melting is occurring. If a structure is melting, subroutine REMESH is called to calculate changes in the structure geometry and to determine if failure by complete melting has occurred.

Finally, STRCHT calculates the net heat transfer from the i

structures in a cell to each of the fluid components in that cell. Subroutine STRPRP uses subroutine MSTRCT to calculate thermal properties for use in the structure heat transfer calculation. STRPRP is called by subroutines STRCHT and STGEOM. Subroutine BCSSET calculates thermal boundary conditions for heat transfer between structures and fluid components. BCSSET treats one boundary condition on surface 1 and three boundary conditions on surface 2 of each structure. BCSSET is called from i subroutine STRCHT. 1 Subroutine TRYSLV is a matrix solution routine for non-symmetric tridiagonal matrices. TRYSLV does not destroy the original matrix and it is called from subroutine STRCHT. i Subroutine REMESH is used to calculate changes in a structures geometry that result from melting. REMESH determines mass transfer for structural material to the fluids. It recalculates the nodal heat transfer mesh in a structure that is melting and it determines structural failure by complete melting. REMESH is called from subroutine STRCHT. Subroutine STGEOM calculates the initial nodal heat transfer mesh for each structure. In

addition, it calculates heat transfer geometric factors and nodal masses for the

~ initialization and whenever the nodal heat transfer mesh is

changed, for example by melting.

STGEOM is called from subroutine MPRGIN during initialization and from subroutine STRUCT at eacn time step. STGEOM calls subroutine STRPRP to get material densities. C.4.2. Mechanics Routines This section desv.ibes the functional layout, purpose and operation of the mechanics routines in STRUCT. The routines are set up functionally as shown in Figure C.S. The first routine, SMOVE, controls movement of structures that failed during the previous time step in the calculation. The information for moving the failed structures, which consists of the destination and original cell numbers and the structure identification (ID) number, is actually calculated during the previous time step, so that the function of SMOVE is to complete the move by switching data pointers. At the completion of SMOVE, 150

MELPROG -thinks that the failed structure is in the destination

cell and'that'it is no longer'in the original cell.

SMOVE also, l-changes the f ailure. flag ISFAIL so that the-structure appears l . intact to the heat transfer routines.

Thus, heat transfer l:

' calculations will continue in the new cell. (l ~ DEBID is the first step in determining the loads on the l-structures. It is responsible for determining where debris is l.. located and which structure is supporting the debris, if any. q The ID numbers of the supporting structures are stored in an array,. LSDBAS. The routine also assigns uni qu e ' identifying numbers to each separate debris bed, taking into account possible dividing walls anu plates. The next

routine, LOAD, sets up the loads for the:

f structures. This is done in three passes. The first pass forms i ~ cell loads on each structure from hydraulic forces, the weight of the part of the structure in that cell, fuel rods being directly supported, or debris, using the support ID numbers calculated in DEBID. Loading forces are determined for both the ' radial and -axial directions. The structures for which loads are found include any failed and relocated pieces of structure, which are termed " lumps" and identified by unique ID numbers greater than NST+1. Pass 2 uses the cel'1 loads found in pass 1 to initialize the ring table. This. table contains a row for each intact or failed i structure and has NSR rows. The row contains the total forces in the radial and axial directions as well as the cell specific forces for each cell that the structure is in. These forces will include the contributions f rom supported structures and all parts of the supporting structure at the end of pass 3: pass 2, however, initialized the ring table to only the self-load for each structure. The data storage arrangement is shown in Table C.1. Each row (structure) has entries for the total radial and exial forcas in the first two positions. The following entries are either total forces in an exial node for axial structures (e.g., walls or columns) or total forces in a radial ring (plates). Table C.1. Ring Load Table Data Structure IDCEL: F tot F tot F 1 F 1 ... F NF N p g p 2 7 2 1: 2: t e q I NSR 151

Pass 3 completes the ring table setup by adding in the forces from the supported structures, using the load table to determine which structures to add. Information about the type of structure is also used to determine the appropriate way to apply the loads. The STRESS routine uses the ring table to get the loads in a cell.for each structure. These loads are passed, along with the structure average temperature., to

MECH, which is a

control routine for the mechanics models. MECH calls MCHPRP to get the mechanical properties at the average structure temperature and i then calls the appropriate mechanical model to evaluate the ) stresses. MECH then calls

LMP, The LarsondMiller Parameter routine, to find the predicted rupture time and update the life fraction.

If the life fraction exceeds 1.0, indicating that the rupture lifetime has been used up, the structure failure flag is set and passed back to STRESS. This is done for each cell containing the structure. STRESS calls UPFAIL after MECH if the failure flag has been set. UPFAIL simply makes sure that the failure flags are also j set for the structure in connected cells, if those parts of the i structure were supported by the failed section of structure. For l instance, if a wall supported from above fails in a node, UPFAIL sets the failure flags in the sections of wall below the failed node. 1 The next main

routine, DYNAM, is responsible for the I

dynamics of moving failed pieces of structure from cell to cell. Currently, DYNAM is a very simple model that assumes all structure pieces eventually end up falling under gravity. The model checks the cells below a failed piece of structure until it encounters an intact structure capable of supporting the failed ~ pieces a support structure is defined to be a plate or column. This structure will receive the load from the failed piece in subsequent time steps. If there is no intervening obstruction, such as a previously failed structure (now called a lump), debris, or fuel rods, then the failed piece is set up to be moved to the cell of the new support structure (or to the cell above it, in the case of the support being a top plate) on the next time step by SHOVE. All of this phase of moving the new lump is done in NEWSUP, called from DYNAM. In the event that no new support is found, the lump will be moved to the bottom cell in the problem. If the support is in the same cell as the lump, no move is performed, but the load is transferred. The move from the old cell to the new is set up by storing pointers to the faileddpiece data and to the new cell in the LSLUMP array for use by SMOVE. The mechanism for changing the loading pattern is to merge ~ l the load chain associated with the structure from which the failed piece came with that of the new support. The lump is l 152

assigned a new structure ID number greater than NST+1 and placed in the support structure's chain directly before the merged . chain.. This means that any structure supported by the failed piece will now load the new support. The chain merging is performed in CNGLOD. .L C.4.3 Flow Areas and Hydraulic Diameters . =. Subroutine STRSET dynamically calculates fluid flow areas and_ hydraulic diameters for use throughout MELPROG. These calculations are performed here rather than in the FLUIDS module because the STRUCTURES module contains the detailed geometric information needed to perform the calculations. Subroutine FA in-the FLUIDS module transfers the flow area and hydraulic diameter into a form used by the FLUIDS module. C.S. DEBRIS Module This section-discusses the main subroutines in the DEBRIS module. Subroutine DEBRIS controls the DEBRIS module ' calculation. -This routine calls the other routines which perform the initialization and analysis of the debria beds. Figure C.6 -l illustrates the calling sequence for the major subroutines in i this module. The first routine, ALFLOC, is used to determine the corium l volume fraction distribution on a two-dimensional basis. This I routine obtains most of its information from FLUIDS arrays. With this information and the approximations discussed in Section 5.3, the local corium volume fractions can be calculated. The second routine, DBFORM, decides whether or not a debris bed is forming based on an empirical blocking correlation. This two dimensional corium volume fraction a routine uses the distribution calculated in subroutine ALFLOC. 1 l If a debris bed forms, then the third routine,

SEARCH, determines whether dryout has occurred.

If the bed is dry, then l subroutine FINDIC is called to determine the initial conditions. -l After a bed has formed subroutine BEDBC is called to determine the boundary conditions. When debris beds exist, subroutine DBSRC prepares the arrays which transfer mass and energy to FLUIDS. The first call to this routine initialized the arrays and the second call actually sets the values. Subroutine DBCALC controls the debris bed calculation. This routine sweeps through all debris beds calling subroutine T RA NS 153

, :1 - y.. DEBRIS

ALR.OC a

---*0BFORM --* SEARCH --+flNDIC --*BEDBC --*DBSRC --*0BCALC

TRANS i

i --+0BSRC 1 --*DBINTR 1 Figure C.6 Functional Diagram of DEBRIS Module Routines i 1 1 l RADWPC MADPNT 1 --*fKPM --*VIEWFC --*SGECO --*SGESL --+WRTLCM ~ 1 Figure C.7: Functional Diagram of RADIATION Module Subroutines 154

l for 'each one. TRANS' solves the detailed heat transfer and waterial relocation equations. After solving the equations at a given time step, subroutine DBINTR is called to obtain information to be transferred to other modules. i - C.6. RADIATION Module l The RADIATION module consists of two main subroutines, l RADMPG and VIEWFC. There ar6 four additional subroutines: an l input and initialization subroutine called RADIN; two pointer I calculation subroutines called RADPNT and RDPNTT and a subroutine to calculate the Planck mean absorption cross section for steam called FKPM. Subroutines RADIN and RDPNTT are called from the INFUT module. RADMPG is the driver subroutine for the ~ RADIATION module. Figure C.7 shows the calling sequence for the RADIATION hodule. Subroutine RADMPG is divided into six sections. There are I four parts to the calculation sequence, a data transfer section and an output section. The data transfer section transfers needed information from each of the other calculational modules to the RADIATION module. The first part of the calculation sequence includes calculation of absorption coefficients, l emissivities, black body emission spectra. Subroutine VIEWFC is j called to calculate view factors and mean transmission i coefficients. In the second part, a coefficient matrix and part of'the resultant vector are assembled for each cell. Also, the coefficient matrices are inverted. In the third

part, an iteration on the cell boundary fluxes is performed.

The changing parts of the resultant vector are calculated and 'the total flux leaving each surface is calculated. The fourth part calculates the net radiation fluxes to each surface and to each fluid component. The final section prints radiation heat fluxes and stores permanent radiation data. RADMPG calls subroutines FKPM, VIEWFC, SGECO, and SOESL. Subroutines SGECO and SGESL are system utility routines used in the matrix solution. Subroutine VIEWFC calculates view factors and mean transmission coefficients. View factors are calculated the first time that VIEWFC is called and whenever the geometry of a cell ( changes. Mean transmission coefficients are calculated each time VIEWFC is called. VIEWFC is called from RADMPG once per time step. VIEWFC calls no subroutines. l t 155 I'

l i C.T. TRAC /MELPROG Link Nmputational Flow Solution of the fluid dynamics equations by means of the SETS method requires at least six passes through each one-dimensional TRAC component for every time step. The first two passes are made in the prepass step to solve the stabilizer e motion equations. This step first obtains the velocities interior to a component as a function of the junction velocities. When the core is being handled by MELPROG, this function is accomplished by the subroutine DCORE1 that in turn uses the section of FLUIDS that calla MOPRED and HOSTAB. Next, the motion equations are applied and a closed set of equations is obtained for the junction velocities. This system is then solved for the junction velocities. The interior velocities are then obtained by a backward substitution using the known junction velocities. For a MELPROG-treated core, this is handled by a call to DCORE1 that in turn calls DCORE10. The next two or more passes are to solve the basic hydrodynamics equations. The flow for the so' called outer iteration is shown in Figure C.8. A single complete pass through this procedure provides the solution for the linearized finite-difference equations. Subsequent passes through the procedure for the same time step produce a Newton iteration on the non linear dif ference equations. a The final two passes made in the postpass are to solve the stabilizing mass and energy equations. The solution of these equations is quite similar to that of the stabilizing motion equations and is handled for a MELPROG core simulation by the section of FLUIDS that calls MSTAB and ESTAB on the first pass. The back substitution pass to complete the solution is handled by DCORE3D. Before returning to TRAC, the final MELPROG fluids pass is followed by calls to the MELPROG routines that perform the rod, radiation, debris bed, and structures calculations. C.B. HELPROG Data Input and Initialization Module (INPUT) MELPROG is a modular code in which the data input and initialization are also modular. Data input and initialization for all modules is controlled by subcoutine MPRGIN:

however, specific input for each calculational module is through its own input subroutine.

Similarly, i ni tial i za t io n for each calculational module is through one or more subroutines written specifically as part of that module. The calling sequence for the major subroutines that are part of the INPUT module is shown in Figure C.9. A brief description of subroutine MPRGIN and the other input routines for each module are given in this section. In addition to acting as the driver routine for the INPUT module, subroutine MPRGIN serves as the input routine for the " general input". MPRGIN calls several subroutines to read and initialize the input for the various calculational modules. As 156

4 I PREPASS ICOMP=1 ICOWPzt YES YES ir g WEL8 ROC WELPROG NO I' BRANCN BRANCN DCORE2 DCORE2 BACK EISN" FORWARD CO NT D RAC COWPONENT DCME2C-l DCORE29 E$DS SUB Tl ON R FORWARD 3, UPDATE EUWlNAllON BOUNDARY UPDATE o ARRAYS SOUNDARY DCORE2D 1r I ARRAY [NTRYIN Rut 05-UPDATE DEPENDENT ICOMPs VARIABLES ICOMP+1 ICOMP-ICOWP+1 LAS NO LAS NO l COMPONENT - YtS f YES CONVDtGED YES POSTPASS 1 j , Figure C.B. Outer iteration flow chart for TRAC /MELPROG. i I t 157 L_____-_________

k. l. I l l MPRGIN 20 SETUP -*VLDIN --*RODIN r'NITIAL ' ~ - +1NITZ O NITGS LDTINT '. --*STRCIN --*STGEOM --*STINIT --+STRSET --+MECHIN rUNK --*LMPINT --+RADIN --+0BINP --+fA --*0C --+SETBCS --+SOLNPRT --*SETPROP

PROPRNT

-+ KNOB Figure C.9: Functional Diagram of INPUT Module Subroutines-O 158

r part of the input processing, MPRGIN checks for sufficient storage space in the X array for each module. (Data storage in the X array is discussed in section C.9.) If any module is found to require more space than is available with the existing array dimensions, an error message is printed which indicates the amount of storage space available and the amount required. If it is found that insufficient space exist for any module, the calculation is stopped af ter completion of the input processing. t Subroutine FLDIN is the data input routine for the FLUIDS module input. FLDIN reads FLUIDS scalar and array data and also reads FLUIDS boundary conditions. FLUIDS initialization is performed after all input has been read for each module. FLUIDS module subroutines FA (flow area and hydraulic diameter calculation), BC and SETBCS (boundary conditions). SETPROP (fluid properties), and KNOB (field equation control) are es11ed from subroutine MPRGIN to initialize the FLUIDS module. Subroutine RODIN is called to read the PINS module input. After the input has been read, the PINS module initialization is performed. Subroutine INITIAL is called to set up the initial node spacing and temperatures in the fuel rod model. Subroutine INITZO is called to initialize the variables used in the clad oxidation and oxygen diffusion routine (COBILD). Fission product concentrations in each axial node are initialized in FPPREP. Parameters used by the fission gas release and fuel rod plenum { pressure models are set up in subroutine INITGS. This routine is skipped if fresh fuel rods are specified by setting the average i burnup to zero (BUAVG 0). Subroutine OLDTINT then calculates clad temperature for each axial cell. an average Subroutine STRCIN is the main input routine for STRUCTURES module

input, teat transfer and mechanical data.

Mechanical property data are read by subroutine MECHIN. In addition to reading data STRCIN calculates array space requirements for the structure heat transfer and mechanical data arrays. The array j lengths are used in the calculation of pointers to locations in i the X array (see section C.9). Initialization of data for the structure heat transfer is performed in subroutine STGEOM (see section C.8). Subroutine STRSET is called as part of STRUCTURES module initialization to calculate flow areas and hydraulic diameters (see section C.4) for use in the FLUIDS module initialization. Several subroutines are used to initialize the structure mechanics data. The first subroutine in this initialization is STINIT. STINIT processes the input problem description to set up a table which specifies the geometry for all structures. The geometry table (Table C.2) contains a row for each distinct o structure. The rous are accessed by the structure identification numbers IDCEL. There are NST of these identification numbers, o r.e for each structure, which are assigned in the input for the 159

Table C.2. Geometry Table Data Structure l IDCEL:. Type Model Material Geometry Data 1: 22 NST: problem. Each row of the geometry table contains the mechanical type number., the mechanical model

number, and geometric information describing the structure for that row.

The next initialization routine, STSET, clears the loading table and assigns some temporary working space for use by the LINK routine. It then calls LINK, which forms a table.that specifies the' loading for each structure. The loading table is constructed from the information given in the input via IDCEL (the structure ID number) and IDSUP (the support structure ID number). The loading table (Table C.3) contains a row for each structure, as in the geometry table. Each row describes a load chain, that is, which structures load which other structures. LINK. forms these load chains using a recursive procedure, given the numbers IDCEL and IDSUP for each cell in the problem. At the end-of

LINK, the load chains describe the complete loading pattern in the problem.

For example, the' core barrel supports the core support plate, which in turn holds up the support

columns, which support the lower core grid plate. - etc.

An identification number greater than or equal to 900 is used in the input to indicate an external support for a structure, as for instance the support for the reactor vessel. MECHIN reads in the mechanical property tables for NMATMX materials from an alternate input file (MIMMCH). These properties are the elastic modulus, yield strength, ultimate Table C.3. Loading Table Data Structure IDCEL: Load ID 1 Load ID 2 Load ID N 1: 2: NST: 160

ten'sile

strength, Poisson's

ratio, and mean coefficient of thermal expansion as functions of temperature.

Following these tables, the rupture time table is read in as a function of stress and temperature. It is necessary to have at least one rupture time for each stress point entered. The rest of the. times for the given stress (at the other temperature points) can be zero. The unknown times are then filled in by the LMPINT routine using interpolation or extrapolation as appropriate. This completes initialization of the mechanics routines. Subroutine RADIN reads the RADIATION module input data. It also calculates array spaco requirements and a special pointer array for location of two-dimensional cell radiation boundary flux 3s. RADIN: is called by MPRGIN. Subroutine DBINP is used to read the DEBRIS module input data. DBINP also initia11zes a number of constants and flags used throughout the DEBRIS module. C.9. Data Management Structure Used in MELPROG Data storage in MELPROG is divided into two categories, permanent and temporary. Permanent data are data that must be kept from time step to time step or that need to be passed from one module to another. Temporary data are data that are transferred into a module from another module and data needed during a time step but not from one time step to another (working storage). All permanent data arrays, referenced by pointers, are stored in a single large array called XLCM. Access to data stored in XLCM is by blocks that correspond to each of the calculation modules. The permanent data are transferred to the working array X for use in a given module. All calculations using array storage, with a few exceptions, are performed in the X array. No direct access to the XLCM array is allowed. On a CDC 7600 computer, the XLCM array is stored in Large Core Memory (LCM) while on a CRAY computer, or any single memory machine, it is stored in main memory. The X array is stored in' Small Core Memory (SCM) on a CDC 7600 and in main memory on CRAY or other single memory computers. Subroutine RWLCM is used to perform the transfer from the XLCM array to the X array and from the X array to the XLCM array. Three entry points (RWLCM, REDLCM, and WRTLCM) are used. RWLCM is used for initialization of the storage in the XLCM array. REDLCM is used for reading data from the XLCM array to the X array. WRTLCM is used to write data f rom the X array to the XLCM array. All data for a particular module is transferred to or l from the XLCM array in a single block except for PINS module ~ data. PINS module data is transferred by blocks that correspond to each radial ring. Currently, the data transfer is necessary for both single memory and multi' level memory computers; however, 161 i

in future versions of MELPROG the data management will be changed for single memory computers. The actual cost of the data transfer on a CRAY is approximately 15 of total computer cpu time. All of the entry points (RWLCM, REDLCM, and WRTLCM) use the same. argument list as. demonstrated here by a call to RWLCM: s CALL RWLCM(NBLOCK, IRING IXADDR, LENGTH) where NBLOCK = internal I/O unit number. IOFLDS for the FLUIDS

module, IOPINS for the PINS module, IOSTRC for the STRUCTURES module, IORAD for the RADIATION

module, and IODBED for the DEBRIS module.

IRING = 1 for all modules, except PINS module. the ring For the PINS module, IRING = number (1 to NRING). IXADDR = beginning address in the X array to transfer data to or from. LENGTH = number of words to be transferred to or from the X array. For data transferred to the X array, LENGTH is returned, while for transfers from the. X array to the XLCM array, LENGTH'must be supplied. lithen data is transferred to the XLCM array from the X array, the ] value of LENGTH must not exceed the value used to initialize the j ILCM array. If a call is made to entry point WRTLCM with a value of LENGTH greater than the value used to initialize the storage, the code will stop on an error. The number ol' modules for which data must be transferred is dependent on each module. The-FLUIDS and RADIATION modules interact with all. the other modules and, therefore, must get information from each of them. The PINS module interacts with the FLUIDS and RADIATION modules. The STRUCTURES module interacts with the FLUIDS, RADIATION,' and DEBRIS modules, while the DEBRIS module interacts with the FLUIDS, PINS, RADIATION, and STRUCTURES modules. Data storage in the X array is divided into three groups. q The groups are permanent data, transfer data, and working space. I I P. the FLUIDS, PINS, and RADIATION modules, the permanent data is c stored first. The transfer data is stored next and the working space storage is at the end. In the STRUCTURES and DEBRIS nodules, the order of storage'is transfer data, permanent data, i 162 i j j

Institute Nuclear Energy Research Lungton Taiwan Republic of China Attn: S. I.. Chang Reactor Centrum Nederland P.O. Box 1 1755 ZG Petten The Netherlands Attn K. J. Brinkman Department. LWR Fuel Belgonucleaire Rue de Champde Mars. 25 B-1050 Brussels, Belgium - Attn: H. Bairiot, Chief Direzione Centrale della Sicurezza Nucleare e della Protezione Sanitaria' (DISP) Entre'Nazionale Energie Alternative (ENEA) Viale Regina Margherita, 125 Casella'Postale N. 2358 1-00100~ Roma A.D., Italy Attn: G. Petrangeli I b 4 s'

...oes...... .,o.. ....o,..~....--.......... g.o E/,"a'# CIBUOGRAPHK: DATA SHEET . n.. o,....... . in 6..=o.a. iiv 6. MELPROG-PWR/MD00: A Mechanistic Code 3 6 8. va *.-- j for Analysis of Reactor Core Melt Progression l and vessel Attack Under Severe Ar-ident Conditions anuan 85 "*' W. J. Camp, J. E. Kelly, P. J. Maudlin, J. L. Tomkins,* M. F. Young, R. J. Henninger** + May 1985 - ...ca ci s . x s -i... .g.g....o ......oo u,,.,. c Reactor Safety Theoretical Physics Sandia National Laboratories NRC FIN No. P.O. Box 5800 A-1342 Albuquerque, NM 87185 ..oo.o..a.i.o.. ..o.......oo uu..,.c Fuel Systems Research Branch I U.S. Nuclear Regulatory Commission Washington, D.C. 20555 ..a *oc co 'au "--- -- t) Sv $...W1.. 0.of.b

• S.A.I.,

Inc.

• • Los Alamos National Laboratory

n. c, au _,

, The MELPROG computer code is being developed to provide detailed, best-estimate, coupled analyses of all the major phenomena involved in the course of a severe accident progression. The initial version of this code, l MELPROG-PWR/ MODO, is described in this report. l The purpose of MELPROG is to provide at any time during an accident se-quence a description of (1) the state of the reactor core and surrounding in-vessel environment, and (2) the disposition of core materials (in l particular, fission products) contained within the reactor coolant system i boundary. MELPROG is coupled to the TRAC-PF1 RCS thermal-hydraulics code to provide an integrated analysis of the behavior of core, vessel, and reactor coolant system during severe accidents. MELPROG treats core de-gradation and loss of geometry, debris formation, core melting, attack on supporting structures, slumping, melt / water interactions and vessel fail-ure. The key element in MELPROG is the use of detailed modeling for the entire damage progression and failure sequence. Emphasis is also placed on the, rates of hydrogen, steam and fission product formation, and trans-port to containment during damage progression. l .. ooco ..o............... u.c. - o., .. p,.,g.;. U ,. u e.. c..n. c..c-(7as.ye . e. ws. wo.....o. o v.. U = (Tn. w. es U ... o... o i. ig Sm ([ e _ - _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ___}}