




















 
















 





    




  

















  

















































࿖㓙න૏♽㧔5+㧕

⴫㧝㧚5+ၮᧄන૏
5+ၮᧄන૏
ၮᧄ㊂ฬ⒓⸥ภ㐳ߐࡔ ࡯ ࠻ ࡞ O
⾰㊂ ࠠࡠࠣ࡜ࡓ MI
ᤨ
㑆
⑽
U
㔚
ᵹࠕ ࡦ ࡍ ࠕ #ᾲജቇ᷷ᐲ ࠤ ࡞ ࡆ ࡦ ‛ ⾰ ㊂ࡕ࡞ OQNశᐲ ࠞ ࡦ ࠺ ࡜ EF⴫㧞㧚ၮᧄන૏ࠍ↪޿ߡ⴫ߐࠇࠆ5+⚵┙න૏ߩ଀
5+ၮᧄන૏
⚵┙㊂ฬ⒓⸥ภ㕙Ⓧᐔ ᣇ ࡔ ࡯ ࠻ ࡞O૕Ⓧ┙ ᴺ ࡔ ࡯ ࠻ ࡞Oㅦ ߐ 㧘 ㅦ ᐲࡔ ࡯ ࠻ ࡞ Ფ ⑽OUടㅦᐲࡔ ࡯ ࠻ ࡞ Ფ ⑽ Ფ ⑽
OU
ᵄ
ᢙᲤ
ࡔ
࡯
࠻
࡞
O
ኒ ᐲ 㧔 ⾰ ㊂ ኒ ᐲ 㧕 ࠠࡠࠣ࡜ࡓᲤ┙ᴺࡔ࡯࠻࡞ MIO
⾰㊂૕Ⓧ㧔Ყ૕Ⓧ㧕 ┙ᴺࡔ࡯࠻࡞Ფࠠࡠࠣ࡜ࡓ OMI㔚ᵹኒᐲ ࠕࡦࡍࠕᲤᐔᣇࡔ࡯࠻࡞
#O⏛ ⇇ ߩ ᒝ ߐࠕ ࡦ ࡍ ࠕ Ფ ࡔ ࡯ ࠻ ࡞#O㧔 ‛ ⾰ ㊂ ߩ 㧕 Ớ ᐲ ࡕ ࡞ Ფ ┙ ᣇ ࡔ ࡯ ࠻ ࡞ OQNOノᐲ ࠞ ࡦ ࠺ ࡜ Ფ ᐔ ᣇ ࡔ ࡯ ࠻ ࡞ EFOዮ᛬₸ 㧔ᢙ‫ߩޓ‬㧕‫ޓ‬㧝⴫㧡㧚5+ធ㗡⺆ਸ਼ᢙ‫ ޓ‬ធ㗡⺆⸥ภ ਸ਼ᢙ‫ ޓ‬ធ㗡⺆


࡛
࠲ 㨅
࠺
ࠪ




࠯
࠲
㨆
࠮ ࡦ ࠴




ࠛ ࠢ ࠨ 㧱
ࡒ
࡝




ࡍ
࠲ 㧼
ࡑࠗࠢࡠ




࠹
࡜ 㨀
࠽
ࡁ




ࠡ
ࠟ 㧳
ࡇ
ࠦ




ࡔ
ࠟ 㧹
ࡈࠚࡓ࠻




ࠠ
ࡠ 㨗
ࠕ
࠻




ࡋ ࠢ ࠻ 㨔
࠯ ࡊ ࠻




࠺
ࠞ FC
࡛ ࠢ ࠻



⸥ภ
F
E
O
v
P
R
H
C
\
[

⴫㧢㧚࿖㓙න૏♽ߣ૬↪ߐࠇࠆ߇࿖㓙න૏♽ߦዻߐߥ޿න૏

⴫㧟㧚࿕᦭ߩฬ⒓ߣߘߩ⁛⥄ߩ⸥ภߢ⴫ߐࠇࠆ5+⚵┙න૏
5+⚵┙න૏
⚵┙㊂
ઁߩ5+න૏ߦࠃࠆ 5+ၮᧄන૏ߦࠃࠆ
ฬ⒓
⸥ภ
⴫ߒᣇ
⴫ߒᣇ
C
TCF
ᐔ
㕙
ⷺ࡜ ࠫ ࠕ ࡦO࡮O D
C
E


┙
૕
ⷺ ࠬ࠹࡜ࠫࠕࡦUTO ࡮O  D
๟
ᵄ
ᢙࡋ
࡞
࠷
*\
U
ജ
࠾ ࡘ ࡯ ࠻ ࡦ
0O࡮MI࡮U࿶ ജ ᔕ ജࡄ ࠬ ࠞ ࡞
2C
0OO࡮MI࡮Uࠛࡀ࡞ࠡ࡯ ઀੐ ᾲ㊂ ࠫ ࡘ ࡯ ࡞,
0࡮O
O ࡮MI࡮UᎿ ₸  ᡼ ኿ ᧤ࡢ
࠶
࠻
9
,U
O࡮MI࡮U㔚 ⩄  㔚 ᳇ ㊂ࠢ ࡯ ࡠ ࡦ%U࡮#㔚૏Ꮕ㧔㔚࿶㧕 ⿠㔚ജ ࡏ
࡞
࠻
8
9#
O࡮MI࡮U࡮#㕒㔚ኈ㊂ࡈ ࠔ ࡜ ࠼(%8O࡮MI࡮U࡮#
㔚
᳇
ᛶ
᛫ࠝ
࡯
ࡓ
ǡ
8#
O࡮MI࡮U࡮#ࠦ ࡦ ࠳ ࠢ ࠲ ࡦ ࠬࠫ ࡯ ࡔ ࡦ ࠬ
5#8O࡮MI࡮U࡮#⏛᧤࠙ ࠛ ࡯ ࡃ
9D
8࡮UO࡮MI࡮U࡮#



⏛
᧤
ኒ
ᐲ࠹
ࠬ
࡜
6
9DO
MI࡮U ࡮#ࠗ ࡦ ࠳ ࠢ ࠲ ࡦ ࠬࡋ ࡦ ࡝ ࡯*
9D#O࡮MI࡮U࡮#࠮ ࡞ ࠪ ࠙ ࠬ ᷷ ᐲ ࠮࡞ࠪ࠙ࠬᐲ F͠శ᧤࡞ ࡯ ࡔ ࡦNO
EF࡮UT EO࡮O࡮EFEF
ᾖ
ᐲ࡞
ࠢ
ࠬ
NZ
NOO
O࡮O࡮EFO࡮EF㧔᡼኿ᕈᩭ⒳ߩ㧕᡼኿⢻ ࡌ ࠢ ࡟ ࡞$SUๆ ෼ ✢ ㊂  ⾰ ㊂ࠛ ࡀ ࡞
ࠣ
࡟
ࠗ
)[
,MI
O࡮Uࠡ ࡯ ಽ ਈ  ࠞ ࡯ ࡑ✢ ㊂ ᒰ ㊂  ๟ ㄝ✢ ㊂ ᒰ㊂ ᣇะᕈ✢㊂ᒰ㊂ ୘ ࠪ ࡯ ࡌ ࡞ ࠻
5X,MIO࡮U
ੱ✢㊂ᒰ㊂⚵❱✢㊂ᒰC ࡜ࠫࠕࡦ෸߮ࠬ࠹࡜ࠫࠕࡦߩ૶↪ߪ‫ޔ‬หߓᰴరߢ޽ߞߡ߽⇣ߥߞߚᕈ⾰ࠍ߽ߞߚ㊂ࠍ඙೎ߔࠆߣ߈ߩ⚵┙න૏ߩ⴫ߒᣇߣߒߡ೑ὐ߇޽ࠆ‫┙⚵ޕ‬න૏ࠍᒻ૞ࠆߣ߈ߩ޿ߊߟ߆ߩ
‫⴫ߪ଀↪ޓ‬㧠ߦ␜ߐࠇߡ޿ࠆ‫ޕ‬D ታ㓙ߦߪ‫⸥ߪߦᤨࠆߔ↪૶ޔ‬ภTCF෸߮UT߇↪޿ࠄࠇࠆ߇‫ޔ‬⠌ᘠߣߒߡ⚵┙න૏ߣߒߡ
ߩ⸥ภ̌㧝̍ߪ᣿␜ߐࠇߥ޿‫ޕ‬E ᷹శቇߢߪ‫ߩࡦࠕࠫ࡜࠹ࠬޔ‬ฬ⒓ߣ⸥ภUTࠍන૏ߩ⴫ߒᣇߩਛߦߘߩ߹߹⛽ᜬߒߡ޿ࠆ‫ޕ‬F ߎߩන૏ߪ‫ࠬ࠙ࠪ࡞࠮࡝ࡒߡߒߣ଀ޔ‬ᐲ㨙͠ߩࠃ߁ߦ5+ធ㗡⺆ࠍ઻ߞߡ↪޿ߡ߽⦟޿‫ޕ‬ᧄ࡟ࡐ࡯࠻ߪ⁛┙ⴕ᡽ᴺੱᣣᧄේሶജ⎇ⓥ㐿⊒ᯏ᭴߇ਇቯᦼߦ⊒ⴕߔࠆᚑᨐႎ๔ᦠߢߔ‫ޕ‬ᧄ࡟ࡐ࡯࠻ߩ౉ᚻਗ߮ߦ⪺૞ᮭ೑↪ߦ㑐ߔࠆ߅໧޿วࠊߖߪ‫ޔ‬ਅ⸥޽ߡߦ߅໧޿วࠊߖਅߐ޿‫ޕ‬ߥ߅‫ోߩ࠻࡯ࡐ࡟ᧄޔ‬ᢥߪᣣᧄේሶജ⎇ⓥ㐿⊒ᯏ᭴ࡎ࡯ࡓࡍ࡯ࠫ㧔http://www.jaea.go.jp㧕ࠃࠅ⊒ାߐࠇߡ޿߹ߔ‫ޕ‬⁛┙ⴕ᡽ᴺੱᣣᧄේሶജ⎇ⓥ㐿⊒ᯏ᭴ ⎇ⓥᛛⴚᖱႎㇱ ⎇ⓥᛛⴚᖱႎ⺖‫ޥ‬319-1195 ⨙ၔ⋵㇊⃗㇭᧲ᶏ᧛⊕ᣇ⊕ᩮ 2 ⇟࿾ 4
㔚⹤ 029-282-6387, Fax 029-282-5920, E-mail:ird-support@jaea.go.jp
This report is issued irregularly by Japan Atomic Energy Agency
Inquiries about availability and/or copyright of this report should be addressed toIntellectual Resources Section, Intellectual Resources Department,
Japan Atomic Energy Agency
2-4 Shirakata Shirane, Tokai-mura, Naka-gun, Ibaraki-ken 319-1195 Japan
Tel +81-29-282-6387, Fax +81-29-282-5920, E-mail:ird-support@jaea.go.jp© Japan Atomic Energy Agency, 2008⴫㧠㧚න૏ߩਛߦ࿕᦭ߩฬ⒓ߣߘߩ⁛⥄ߩ⸥ภࠍ฽߻5+⚵┙න૏ߩ଀
5+⚵┙න૏
ฬ⒓
⸥ภ
5+ၮᧄන૏ߦࠃࠆ⴫ߒᣇ☼ᐲࡄ ࠬ ࠞ ࡞ ⑽2C࡮UO ࡮MI࡮Uജ ߩ ࡕ ࡯ ࡔ ࡦ ࠻ ࠾ࡘ࡯࠻ࡦࡔ࡯࠻࡞0࡮OO ࡮MI࡮U⴫㕙ᒛജ ࠾ࡘ࡯࠻ࡦᲤࡔ࡯࠻࡞
0O
MI࡮U



ⷺ
ㅦ
ᐲ࡜ ࠫ ࠕ ࡦ Ფ ⑽TCFUO࡮O ࡮U Uⷺടㅦᐲ ࡜ ࠫ ࠕ ࡦ Ფ ᐔ ᣇ ⑽ TCFU O࡮O࡮UUᾲ ᵹ ኒ ᐲ  ᡼ ኿ ᾖ ᐲ ࡢ࠶࠻Ფᐔᣇࡔ࡯࠻࡞
9O
MI࡮U



ᾲ ኈ ㊂  ࠛ ࡦ ࠻ ࡠ ࡇ ࡯ ࠫࡘ࡯࡞Ფࠤ࡞ࡆࡦ
,O ࡮MI࡮U ࡮⾰ ㊂ ᾲ ኈ ㊂ 㧔 Ყ ᾲ ኈ ㊂ 㧕  ࠫࡘ࡯࡞Ფࠠࡠࠣ࡜ࡓ
, MI࡮- O࡮U࡮-
⾰ ㊂ ࠛ ࡦ ࠻ ࡠ ࡇ ࡯ Ფࠤ࡞ࡆࡦ
⾰ ㊂ ࠛ ࡀ ࡞ ࠡ ࡯ࠫࡘ࡯࡞Ფࠠࡠࠣ࡜ࡓ,MIO ࡮U ࡮㧔 Ყ ࠛ ࡀ ࡞ ࠡ ࡯ 㧕ࡢ࠶࠻Ფࡔ࡯࠻࡞Ფࠤ
ᾲવዉ₸
9 O࡮- O࡮MI࡮U࡮-࡞ࡆࡦࠫࡘ࡯࡞Ფ┙ᣇࡔ࡯࠻૕ Ⓧ ࠛ ࡀ ࡞ ࠡ ࡯
O ࡮MI࡮U
,O
࡞


㔚
⇇
ߩ
ᒝ
ߐࡏ࡞࠻Ფࡔ࡯࠻࡞
8OO࡮MI࡮U ࡮#ࠢ࡯ࡠࡦᲤ┙ᣇࡔ࡯࠻


૕
Ⓧ
㔚
⩄
O
࡮U࡮#
%O
࡞
ࠢ࡯ࡠࡦᲤᐔᣇࡔ࡯࠻


㔚
᳇
ᄌ
૏
O ࡮U࡮#
%O
࡞




⺃
㔚
₸ ࡈࠔ࡜࠼Ფࡔ࡯࠻࡞(OO ࡮MI ࡮U ࡮#


ㅘ
⏛
₸ ࡋࡦ࡝࡯Ფࡔ࡯࠻࡞*OO࡮MI࡮U ࡮#ࡕ ࡞ ࠛ ࡀ ࡞ ࠡ ࡯ࠫ ࡘ ࡯ ࡞ Ფ ࡕ ࡞,OQN
O ࡮MI࡮U ࡮OQN
ࡕ ࡞ ࠛ ࡦ ࠻ ࡠ ࡇ ࡯  ࠫࡘ࡯࡞Ფࡕ࡞Ფࠤ࡞, OQN࡮- O࡮MI࡮U࡮-࡮OQN
ࡕ
࡞
ᾲ
ኈ
㊂ ࡆࡦ

ᾖ ኿ ✢ ㊂ 㧔 㨄 ✢ ෸ ߮ ǫ ✢ 㧕 ࠢ࡯ࡠࡦᲤࠠࡠࠣ࡜ࡓ%MI
MI ࡮U࡮#


ๆ
෼
✢
㊂
₸ࠣ ࡟ ࠗ Ფ ⑽
)[U
O ࡮U





᡼
኿
ᒝ
ᐲ ࡢ࠶࠻Ფࠬ࠹࡜ࠫࠕࡦ9UTO ࡮O ࡮MI࡮U O ࡮MI࡮Uࡢ࠶࠻Ფᐔᣇࡔ࡯࠻࡞





᡼
኿
ノ
ᐲ
࡮UT
O
࡮O
࡮MI࡮U
MI࡮U
9
O
Ფࠬ࠹࡜ࠫࠕࡦ
⚵┙㊂

ฬ⒓
ಽ
ᤨ
ᣣ
ᐲ
ಽ
⑽
࡝࠶࠻࡞
࠻ࡦ
ࡀ࡯ࡄ
ࡌ࡞

⸥ภ
5+න૏ߦࠃࠆ୯
OKP OKPUJJOKPUF
FJUq q Ǹ TCF
̉ ̉  q Ǹ TCF
̍ ̍  ̉ Ǹ TCF
N‫ޔ‬. NFOOVV MI0R 0R$$  NP 0R

⴫㧣㧚࿖㓙න૏♽ߣ૬↪ߐࠇߎࠇߦዻߐߥ޿න૏ߢ
5+න૏ߢ⴫ߐࠇࠆᢙ୯߇ታ㛎⊛ߦᓧࠄࠇࠆ߽ߩ
ฬ⒓⸥ภ
5+න૏ߢ޽ࠄࠊߐࠇࠆᢙ୯㔚 ሶ ࡏ ࡞ ࠻ G8 G8  ,⛔৻ේሶ⾰㊂න૏ W W  MIᄤ ᢥ න ૏ WC WC  O⴫㧤㧚࿖㓙න૏♽ߦዻߐߥ޿߇࿖㓙න૏♽ߣ‫ޓޓޓ‬૬↪ߐࠇࠆߘߩઁߩන૏
ฬ⒓⸥ภ
5+න૏ߢ޽ࠄࠊߐࠇࠆᢙ୯
ᶏ㉿㧝ᶏ㉿Oࡁ࠶࠻㧝ࡁ࠶࠻㧝ᶏ㉿Ფᤨ  OUࠕ࡯࡞ C 㧝CFCOOࡋ ࠢ ࠲ ࡯ ࡞ JC 㧝JCJOOࡃ࡯࡞ DCT 㧝DCT/2CM2CJ2C2C
ࠝࡦࠣࠬ࠻ࡠ࡯ࡓ Έ 㧝ΈPOOࡃ࡯ࡦ D 㧝DHOO⴫㧥㧚࿕᦭ߩฬ⒓ࠍ฽߻%)5⚵┙න૏ฬ⒓⸥ภ 5+න૏ߢ޽ࠄࠊߐࠇࠆᢙ୯ࠛ࡞ࠣ GTI GTI ,࠳ࠗࡦ F[P F[P0

ࡐ
ࠕ
࠭
2
2F[P࡮UEO 2C࡮Uࠬ ࠻ ࡯ ࠢ ࠬ 5V 5VEOUOU

ࠟ
࠙
ࠬ
)
) 6ࠛ ࡞ ࠬ ࠹ ࠶ ࠼ 1G 1G Ǹ #Oࡑ ࠢ ࠬ ࠙ 㨴 ࡞ /Z /Z9Dࠬ࠴࡞ࡉ UD UDEFEOEFOࡎ࠻ RJ RJNZࠟ࡞ )CN )CNEOUOU

⴫㧚࿖㓙න૏ߦዻߐߥ޿ߘߩઁߩන૏ߩ଀
ฬ⒓⸥ภ
5+න૏ߢ޽ࠄࠊߐࠇࠆᢙ୯
ࠠࡘ࡝
࡯ %K %K$S࡟ ࡦ ࠻ ࠥ ࡦ 4 4%MI࡜࠼ TCF TCFE)[)[࡟
ࡓ TGO TGOE5X5X

:
✢
න
૏
:WPKV POࠟࡦࡑ ǫ ǫP66ࠫ ࡖ ࡦ ࠬ ࠠ ࡯ ,[ ,[9࡮Ox*\

ࡈ
ࠚ
࡞
ࡒ
HGTOKHO Oࡔ࡯࠻࡞♽ࠞ࡜࠶࠻
OGVTKEECTCVOI MI࠻࡞ 6QTT 6QTT  2C
ᮡ Ḱ ᄢ ᳇ ࿶ CVO CVO2C
ࠞ
ࡠ
࡝
࡯ ECN

ࡒ
ࠢ
ࡠ
ࡦ
v
vvO O㧔╙㧣㧘ᐕᡷ⸓㧕

JAEA-Data/Code 2008-014

Steam Explosion Simulation Code JASMINE v.3 User’s Guide
Kiyofumi MORIYAMA, Yu MARUYAMA and Hideo NAKAMURA
Reactor Safety Research Unit
Nuclear Safety Research Center
Japan Atomic Energy Agency
Tokai-mura, Naka-gun, Ibaraki-ken(Received April 28, 2008)

A steam explosion occurs when hot liquid contacts with cold volatile liquid. In this phenomenon,
ﬁne fragmentation of the hot liquid causes extremely rapid heat transfer from the hot liquid
to the cold volatile liquid, and explosive vaporization, bringing shock waves and destructive
forces. The steam explosion due to the contact of the molten core material and coolant water
during severe accidents of light water reactors has been regarded as a potential threat to theintegrity of the containment vessel. We developed a mechanistic steam explosion simulation
code, JASMINE, that is applicable to plant scale assessment of the steam explosion loads. This
document, as a manual for users of JASMINE code, describes the models, numerical solution
methods, and also some veriﬁcation and example calculations, as well as practical instructions
for input preparation and usage of the code.

Keywords: Severe Accident, Molten Core, Steam Explosion, Fuel-Coolant Interaction, Premixing, Propagation, Numerical Simulation, Safety Assessment, JASMINE

ii

JAEA-Data/Code 2008-014

水蒸気爆発解析コード JASMINE v.3 ユーザーズガイド
日本原子力研究開発機構 安全研究センター
原子炉施設安全評価研究ユニット
森山 清史・丸山 結・中村 秀夫(2008 年 4 月 28 日 受理)

水蒸気爆発は, 高温の液体が低温かつ揮発性の液体に接触するとき, 高温液体の細粒化により急激な伝熱と蒸気発生が生じ, 衝撃波や破壊力をもたらす現象である. 軽水炉シビアアクシデント時の炉心融体と冷却水の接触による水蒸気爆発は, 格納容器健全性への脅威となり得る現象として重要視されてきた. 著者らは水蒸気爆発の機構論的シミュレーションコード JASMINE を開発した.
JASMINE は実機規模の水蒸気爆発負荷の評価に適用可能なコードである. 本報告書は, JASMINE
コードのユーザーのためのマニュアルとして, モデル, 数値解法, コードの検証のための計算, その他の計算例, 及び, 実際的なコード使用に必要な入力データ作成と計算の実行方法等について解説したものである.

原子力科学研究所 (駐在)：〒 319-1195 茨城県那珂郡東海村白方白根 2-4iiii

JAEA-Data/Code 2008-014

Contents
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Background—Steam Explosion and Safety Assessment . . . . . . . . . . . . . . .
1.2 JASMINE Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Usage of This Manual . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
2
4

2 Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1 Modeling Framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Melt Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Melt Jet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 Melt Pool . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 Melt Particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.4 Extension of Melt Particle Model for Explosion Process . . . . . . . . . .
2.3 Two-Phase Flow Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Basic Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Constitutive Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5
5
6
6
8
11
18
21
21
22

3 Numerical Solution Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1 Coupling of Melt and Two-Phase Flow Models . . . . . . . . . . . . . . . . . . .
3.2 Numerical Solution Method for Melt Model . . . . . . . . . . . . . . . . . . . . .
3.2.1 Melt Jet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Melt Pool . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.3 Melt Particle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Numerical Solution Method for Two-Phase Flow . . . . . . . . . . . . . . . . . .
3.3.1 Overall Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.2 Finite Diﬀerence Form of Basic Equations . . . . . . . . . . . . . . . . . .
3.3.3 Boundary Condition Settings . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.4 Newtonian Iteration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.5 Organization of the Pressure Equations . . . . . . . . . . . . . . . . . . .
3.3.6 Elements of the Jacobian Matrix . . . . . . . . . . . . . . . . . . . . . . .

28
28
28
28
31
33
34
34
34
39
40
41
42

4 Veriﬁcation and Example Calculations . . . . . . . . . . . . . . . . . . . . . . . .
4.1 Check of the Numerical Behavior of the Models . . . . . . . . . . . . . . . . . . .
4.1.1 Melt Jet and Pool . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.2 Melt Particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.3 Shock Wave Propagation in Two-Phase Medium . . . . . . . . . . . . . .
4.2 Simulation of Premixing Experiments . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Simulation of Explosion Experiments . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Application for Plant Scale Simulation . . . . . . . . . . . . . . . . . . . . . . . .

47
47
47
51
55
60
64
75

5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

83

iii
iii

JAEA-Data/Code 2008-014

A Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

83

B Particle Heat Conduction Models . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.1 Surface Temperature Drop of Melt Particles . . . . . . . . . . . . . . . . . . . . .
B.2 Fast Transient Heat Release from Fine Fragments . . . . . . . . . . . . . . . . . .

86
86
87

C Brief Description of CIP Method . . . . . . . . . . . . . . . . . . . . . . . . . . .
C.1 Basic Concept of CIP Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C.2 Improvement for Steep Steps—CIP1 Scheme . . . . . . . . . . . . . . . . . . . . .
C.3 Non-Linear Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

90
90
91
91

D Steam Table Package WRSTEAMTAB . . . . . . . . . . . . . . . . . . . . . . .

93

E Input Preparation and Code Usage . . . . . . . . . . . . . . . . . . . . . . . . . . 97
E.1 Input File Format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
E.1.1 Preprocessing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
E.1.2 Input Variables and Their Ordering . . . . . . . . . . . . . . . . . . . . . 97
E.2 Running the Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
E.2.1 Normal Run . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
E.2.2 Output Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
E.2.3 Restart Run . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
E.2.4 Restart Calculation in “explosion mode” . . . . . . . . . . . . . . . . . . . 113
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

iv

iv

JAEA-Data/Code 2008-014

目次
1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 背景—水蒸気爆発と安全評価 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 JASMINE コード . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .2
1.3 本書の使用法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 モデル説明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1 モデルの枠組 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 メルトモデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2.1 メルトジェット . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2.2 メルトプール . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2.3 メルト粒子 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.4 爆発過程のためのメルト粒子モデル拡張 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3 二相流モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3.1 基礎式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3.2 構成式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3 数値解法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.1 メルトモデルと二相流モデルの連成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2 メルトモデルの数値解法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2.1 メルトジェット . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2.2 メルトプール . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2.3 メルト粒子 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.3 二相流モデルの数値解法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3.1 概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3.2 基礎式の差分形 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3.3 境界条件の設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.3.4 Newton 法の反復計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .40
3.3.5 圧力方程式系の構成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.3.6 Jacobi 行列の要素 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4 検証計算及びその他計算例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.1 各モデルの数値的挙動確認 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .47
4.1.1 メルトジェット及びプール . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.1.2 メルト粒子 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.1.3 二相媒体における衝撃波伝播 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.2 粗混合実験のシミュレーション . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.3 水蒸気爆発実験のシミュレーション . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.4 実機規模解析への適用例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5 まとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80v

v

JAEA-Data/Code 2008-014

謝辞 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
付録 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

A 記号表 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
B 粒子内の熱伝導モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
B.1 メルト粒子の表面温度降下 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
B.2 細粒からの急速放熱 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .87
C CIP 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
C.1 CIP 法の概念 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
C.2 急峻な変化がある場合の改良–CIP1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .91
C.3 非線型方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
D 蒸気表パッケージ WRSTEAMTAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
E 入力データ作成及びコード使用法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
E.1 入力データの形式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
E.1.1 プリプロセス . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
E.1.2 入力変数とその順序 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
E.2 コードの実行 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
E.2.1 通常の実行 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
E.2.2 出力ファイル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
E.2.3 リスタート計算の実行 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
E.2.4 “爆発モード” でのリスタート計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
参考文献 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

vi

vi

JAEA-Data/Code 2008-014

1 Introduction
1.1

Background—Steam Explosion and Safety Assessment

A steam explosion is a phenomenon in which high temperature liquid contacts with low temperature volatile liquid and transfer heat extremely quickly, causing rapid vaporization and shock
wave generation. In other words, a steam explosion is a thermodynamic process which converts
a part of the internal energy possessed by the high temperature liquid into mechanical energy,
which can produce mechanical damages to the surroundings.
In the ﬁeld of nuclear safety, the steam explosion due to the contact of the molten corematerial and coolant water has been recognized as a potential threat to the integrity of thecontainment vessel during severe accidents of light water reactors. The interaction between themolten fuel and coolant in nuclear power plants is also called fuel-coolant interaction (FCI) ina broad sense, which include both the explosive and mild interactions.
The steam explosion is phenomenologically understood as four stages[1, 2] as shown in
Fig. 1.1. The phenomena involved in each stage is described brieﬂy as follows.
• Premixing: coarse break-up and mixing of the high temperature liquid in the low temperature liquid (coolant) while thermally insulated by a vapor ﬁlm; time scale 0.1–several
seconds• Triggering: destabilization of the vapor ﬁlm and onset of rapid vaporization and fragmentation of the melt droplets• Propagation: ﬁne fragmentation of the melt droplets and rapid heat transfer associatedwith a propagating shock wave; time scale a few milliseconds• Expansion: vaporization of the coolant and expansion of the mixture potentially causing
damages on the surroundings
Vapor
film

Melt drop

1) Premixing

Pressure wave

2) Triggering

3) Propagation

4) Expansion

Figure 1.1: Phenomenological four stages of steam explosions.
Plant scale assessment on the steam explosion issue needs a certain analytical approachin order to extrapolate the knowledge obtained from experiments to the scale and conditions
relevant to power plant accidents.
Historically, a variety of models and methods were developed ranging from thermodynamic
models, which give the bounds of energy conversion eﬃciency, to multidimensional transient
thermohydraulic models, which mechanistically simulates the process.[2, 3]
According to the phenomenology of 4 stages described in the previous section, the assessmentof steam explosion energetics is usually done through the following two steps:−1−

1

JAEA-Data/Code 2008-014• Evaluation of the initial internal energy inventory in the high temperature liquid which is“premixed” with the coolant and ready to participate in the following explosive process,and• Evaluation of the mechanical energy output through the “propagation” and “expansion”(explosion) processes.
Note that the “triggering” is a tricky process from the view point of safety assessment.
It was observed, in some experimental conditions, that the triggering sometimes occurred but
sometimes not, depending on even a small uncontrollable change in the condition.[4] While,in the analytical assessment, assuming a triggering or not makes totally diﬀerent outputs, i.e.
explosion loads or nothing.
Considering that, we, at present, recommend to assume a strong enough triggering, andto do the explosion simulation to see the mechanical loads. Actually, it may happen that no
signiﬁcant explosive interaction occurs even if a strong triggering is assumed, depending on thepremixture condition, e.g. high void fraction or melt solidiﬁcation. In such a case, we may say,with some conﬁdence, that at most insigniﬁcant damage is predicted in the given condition.

1.2

JASMINE CodeWe developed a mechanistic FCI simulation code, JASMINE (JAEA Simulator for Multiphase
INteractions and Explosions).
The ﬁrst version of JASMINE developed in 1995 [5] was a three-ﬂuid Eulerian model withwater, vapor and melt ﬁelds aimed at the premixing simulation. This was then called JASMINEpre and used in the OECD International Standard Problem No.39 (ISP-39) exercise on a large
scale melt quenching experiment FARO-L14 [6].
The propagation model was separately developed as JASMINE-pro [7, 8] which was a ﬁveﬂuid model with additional two ﬁelds for the ﬁne fragments of melt generated in the propagation
stage and small amount of coolant which interacts with the ﬁne fragments.
A major modiﬁcation in JASMINE-pre was made in 1997. In the new version, JASMINE-prev.2, melt droplets generated by coarse break-up was modeled in Lagrangian framework. Finally,in 2003, the propagation related models in JASMINE-pro v.1.1 were merged into JASMINEpre v.2.2. This change produced JASMINE v.3, which can handle both the premixing andpropagation/expansion simulations. JASMINE-pre v.2 and JASMINE v.3 were used in theOECD’s international cooperative research program SERENA Phase-1[9, 10].1
The structure of the present code is illustrated in Fig. 1.2 in relation to the phenomena
expected in the fuel-coolant interaction in a nuclear power plant. JASMINE v.3 consists of
two separate parts for the modeling of the molten core behavior and the coolant multiphaseﬂow. The molten core model called MELT includes three sub-models for the melt jet, melt pooland melt particles. The multiphase ﬂow model, which handles the coolant thermohydraulics is a
modiﬁed version of ACE3D code developed at Japan Atomic Energy Research Institute (JAERI)
by Ohnuki et al. [11]. The melt jet and melt pool models are one-dimensional representations ofa molten core stream falling into a water pool and a continuous melt body agglomerated on thebottom, respectively. The melt particles generated by the melt jet break-up in the water pool
are modeled based on a Lagrangian grouped particle concept.
JASMINE simulates the whole process of the steam explosion in two steps. First, the premixing stage is simulated with the initial/boundary conditions for the ﬂow ﬁeld and melt inlet.
Then, the the explosion process—propagation and expansion—is simulated by feeding output

1

The latest version is 3.3a as of April, 2008. The basic structure and functions are the same for versions 3.x.

−2−

2

JAEA-Data/Code 2008-014

Ex-vessel FCI

Premixing mode
Melt jet breakup in
premixing process
Coolant: Two-phaseflow code ACE-3D

Melt jet
1D CIP

Melt droplets
flow and coolingGroup-particle model(Np particles in it)

Melt particles
Lagrangianrx

Debris bed
formation
Melt pool
1D FDMrz
Dp(xp, zp)(Vx,Vz)zx

Explosion mode
Particlegroupspremixing

Shock
front

Particlegroups"Fragmented"

Fragments
Particles

Figure 1.2: Concept and structure of JASMINE code.

3−−3

JAEA-Data/Code 2008-014

Table 1.1: Nature of physics in premixing and explosion processes and related terminology usedin this document.
Premixing
Explosion (propagation & expansion)
Physics
Dominating
physical gravity (buoyancy) driven con- shock wave propagation
process
vectionTime scale
0.1–10s
0.1–10ms
Melt size scale
1–10mm1–100µm
Terminology
Name of melt size (coarse) break-up of melt jet or (ﬁne) fragmentation of meltchange process
particles (droplets)
particles (droplets)
Name of melt particles
droplets (particles)(ﬁne) fragments
data from the premixing calculation at a selected time, as an initial condition, and by specifying
explosion model options in the input ﬁle.
Table 1.1 summarizes the individual features of the premixing and explosion processes andrelated terminology, that is used in this document.

1.3

Usage of This Manual

The present document describes the models, numerical solution methods, some veriﬁcation andexample calculations, as well as the practical usage of the code.
The instructions on the usage and input preparation, the minimum necessary information
to run the code quickly, are given in Appendix E.
It is recommended for the readers, at least, to brieﬂy go through Chapter 2 before runningthe code, so that they can grab the basic concept of models and deﬁnition and meaning of thetechnical terms.
Chapters 3 and 4 give an explanation of numerical solution methods and veriﬁcation/example
calculations, respectively. Notation of the variables is given in Appendix A. Appendix B–D
describes some modeling and numerical method details.
The summary, Chap. 5, includes the information on the limitation and problems in thepresent model.

4−−4

JAEA-Data/Code 2008-014

2 Model Description
2.1

Modeling Framework

JASMINE consists of a melt model and a two-phase ﬂow model as illustrated in Fig. 1.2. The
melt model consists of three components of sub-models: one-dimensional models for melt jet andpool, and Lagrangian grouped-particle model for melt particles. For the explosion calculation,an additional component is attached to the melt particle model for the ﬁne fragments generatedby shock wave–melt particle interaction. The coolant ﬂow involving water, vapor and noncondensible gases are handled by a modiﬁed version of two-phase ﬂow code ACE3D.
The calculation domain considered is a sector of a cylindrical volume with azimuthal width(angle) Θ. The volume is discretized x-z plain into a ﬁnite diﬀerence grid. The coolant ﬂow
model uses this grid. The melt jet is modeled as a vertical one-dimensional ﬂow accommodated in
the central column of the grid. The melt pool is modeled as a horizontal (radial) one-dimensional
ﬂow in the bottom layer of the grid. The melt particle groups, modeling the melt droplets, areproduced at the jet surface and migrate on the x-z grid as Lagrangian elements.
The basic equations and rules as well as constitutive models are described in the following
sections.
In evaluating the constitutive equations, it is often necessary to consider the multi-phaseﬂow regimes of the coolant ﬂow. We constructed the constitutive equations on a basic idea thatthe melt interacts with the liquid and gas components of the coolant according to the coolant
two-phase ﬂow regimes. The ﬂow regimes are deﬁned as follows.
• Bubbly (liquid-continuous): when the total void fraction α (vapor plus non-condensible
gases) is less than 0.3.
• Droplet (gas-continuous): when α is more than 0.75, and• Transition: when α is between 0.3 and 0.75.
Based on the above deﬁnition, normally the constitutive equations for melt-coolant interactions
are evaluated separately for the gas and liquid components, and averaged with the weighting
function:(0 ≤ α < 0.3) 0(α − 0.3)/0.45 (0.3 ≤ α ≤ 0.75) .
fα =
(2.1)

 1
(0.75 < α ≤ 1)

For example, variables Xliq and Xgas evaluated for the eﬀects of liquid and gas, respectively, areaveraged by
Xtwo-phase = (1 − fα )Xliq + fα Xgas .
(2.2)
Hereafter, this averaging scheme is applied if no speciﬁc description is given for the compositionof two-phase contributions.

−5−5

JAEA-Data/Code 2008-014

2.2

Melt Model

2.2.1

Melt Jet

Basic equations
One-dimensional conservation equations in z-direction for the melt jet are solved. The equations
for mass, internal energy and momentum are expressed as follows:�∂AJ ρJ∂vJ∂AJ ρJ+ vJ= − 2ΘAJ me − AJ ρJ
∂t
∂z
∂z
�

∂eJ
∂eJ
+ vJ
=−
∂t
∂z

2Θ q
AJ ρJ

for mass,

for energy, and

∂vJ
1 ∂pa
∂vJ
+ vJ
=−
+ Kf (va − vJ ) + g
∂t
∂z
ρJ ∂z

for momentum,

(2.3)
(2.4)
(2.5)

where AJ , ρJ , vJ and eJ denote the cross section, density, velocity and internal energy of themelt jet, respectively. Note that the area AJ is for the sector of angle Θ, i.e. AJ = ΘDJ2 /8
where DJ is the jet diameter. The pressure pa in the momentum equation is given by the outerﬂow ﬁeld. Some variables in the source terms: me , mass ﬂux due to droplet entrainment fromthe jet surface, q, heat ﬂux and Kf , friction factor, are to be given by constitutive equations.
Constitutive models
Jet break-up The constitutive equation for the mass entrainment from the jet, i.e. droplet
production from the jet surface, was deduced by relating the empirical correlation for jet breakup length and the mass balance of the jet. Figure 2.1 illustrates a situation that a melt jet iseroded in a water pool and reduces its diameter as it proceeds downward. The jet is totallybroken up at depth Lbrk , which is called the break-up length. Empirical correlations are available
for this break-up length. By considering the mass balance at a certain depth z, the change ofthe jet cross section dA, and the volume ﬂux for the droplet entrainment Ve [m/s] are related by
Ve = −VJ dA,
πDJ dz(2.6)

where A = πDJ2 /4. This equation connects the jet proﬁle, i.e. the cross section or the diameteras a function of the depth, and the entrainment volume ﬂux. The entrainment volume ﬂux isdetermined if a proﬁle is assumed. If it is assumed that the jet diameter decreases linearly withthe depth, i.e. dDJ /dz = const. = DJi /Lbrk , Eq. (2.6) givesVe =VJ DJi.
2 Lbrk(2.7)

If it is assumed that the jet cross section decreases linearly with the depth, i.e. dA/dz = const. =Ai /Lbrk , Eq. (2.6) gives
2
VJ DJi.
(2.8)
Ve =4DJ Lbrk
The variables DJi and Ai are the jet diameter and cross section at the water surface.
The present version of JASMINE uses Eq. (2.7) when the jet break-up length is shorter
than water pool depth Hpl . For�the case of a shallow water pool where Lbrk > Hpl (incomplete

break-up), a correction factor Hpl /Lbrk is attached to Eq. (2.7) based on a consideration on
the strength of the steam upward ﬂow surrounding the melt jet[12].
−6−6

JAEA-Data/Code 2008-014

AiMelt jet

DJi

Water pool

VJz

Steam flow

A
Ve

DJ

Lbrk

Dropletentrainment

Figure 2.1: Concept of the jet break-up model.
A set of empirical correlation for the jet break-up length was obtained by examination intoexperimental data by the authors and in the literature[13]. According to a jet diameter-based
Bond number,2
ρJ gDJi
,
(2.9)
BoJ =
σJ
the jet break-up length is given by a Taylor type correlation[14] when BoJ ≥ 50,
�

Lbrk
ρJ
= 10
DJi
ρl

�1/2

(2.10),or otherwise by Saito et al. correlation[15],
�

Lbrk
ρJ
= 2.1
DJi
ρl

�1/2 �

2
VJi
gDJi

�1/2

,

(2.11)

where ρJ and ρl are densities of the jet and water, VJi is the jet velocity at the water surface,and g is the gravitational acceleration.
The mass ﬂux of droplet entrainment me is given byme = Cent ρJ Ve ,(2.12)

where Cent is a constant. If the constant Cent is set to be unity, the jet breaks up at the length
determined by the empirical correlation. A larger value for Cent shortens the break-up length.
The diameter of the entrained droplet is, right now, left to be given as a constant value by
the user. Successive break-up of the droplets are considered in the particle model.
The x and z velocity components of the entrained droplet is given byvex =Cvx (2Ve ), andvez = Cvzwt vJ + (1 − Cvzwt )vc ,(2.13)(2.14)

where Cvx and Cvzwt are parameters which should have the values roughly ∼ 5 and ∼ 0.5,
respectively.
−7−

7

JAEA-Data/Code 2008-014

Friction The friction coeﬃcient is given by considering the relative ﬂow of the coolant in thesurrounding boundary layer,�ρcΘ
Kf = Cf J f
vr .
(2.15)
ρJ 2AJ
The constant Cf J is usually set to be unity. The friction factor f is given as follows:f= max(16/Red , 0.0791/Re0.25d ) ,(2.16)

Red = ρc dbnd vr /µc ,dbnd =(2.17)

0.37zbnd /Re0.2
z , and(2.18)

Rez = ρc zbnd vr /µc ,(2.19)

where Red and Rez are Reynolds numbers based on the boundary layer thickness dbnd and the
distance from the jet leading edge zbnd . Equation (2.18) is based on the development of turbulent
boundary layer on a ﬂat plate[16].
Heat transfer Heat transfer from the jet is currently not considered, because both the surface
area and the duration of travel in the coolant are much smaller than those of the melt particles
or the melt pool. Therefore, the heat ﬂux q is simply set 0.

2.2.2

Melt Pool

Basic equations
When the jet column or a particle group reaches the bottom, it may form a continuous melt body.
The melt pool model was introduced to contain such a mass. One-dimensional conservation
equations for the pool in the radial direction are solved. The conservation equations for themelt pool are given as follows:∂(whP ρP vP )∂whP ρP
+
= wms
∂t
∂x

for mass,

∂eP
∂eP
q
(ems − eP )max(ms , 0)
+ vP
=
+
∂t
∂x
h P ρP
hP ρP
∂vP
∂t

∂vP
1
+ vP
=−
∂x
ρP

�

(2.20)

for energy, and(2.21)

�∂pa g ∂(hP ρP )
+
∂x
2
∂x
(vms − vP ) max(ms , 0)
+ Kf a (va − vP ) − Kf w vP +h P ρPfor momentum,(2.22)

where hP , w(x)(≡ xΘ), ρP , eP and vP denote the pool height, the sector width of the analysis
domain at position x, density, internal energy and velocity, respectively. Variables ms , q, Kf aand Kf w denote the mass ﬂux of melt jet and particles falling on to the pool surface, the heat ﬂux
at the pool surface, the pool-ﬂuid friction factor and the pool-ﬂoor friction factor, respectively.
Constitutive models
Friction The friction is neglected (friction factors are set to be 0) at present. The friction
is not important relatively to the heat transfer and mass exchange for the primary purpose ofthis pool model, i.e. to accommodate the agglomerated mass of melt in a physically reasonable
manner.

−8−

8

JAEA-Data/Code 2008-014

Heat transfer The surface heat ﬂux is given by considering radiation, ﬁlm boiling, nucleateboiling and convection heat transfer, with a collection of available correlation equations describedin the following. Figure 2.2 shows the boiling curve considered in combining the eﬀects of various
modes of heat transfer. Boiling is considered when the temperature of the heat transfer surface
is higher than the saturation temperature of the coolant. Convection heat transfer to the coolant
liquid and/or gas is considered if it is suitable.
• Convection: natural convection over a horizontal ﬂat surface [17]
= N ucvf (λ/l)(Tw − T∞ ) ,qcvf(2.23)= 0.15Ra1/3 ,N ucvf(2.24)Ra = GrP r = {gβ(Tw − T∞ )l /ν }P r ,3(2.25)

2

where N ucvf denotes the average Nusselt number over the surface of characteristic scale l;Ra, Gr and P r denote Rayleigh, Grashof and Prandtl numbers. The melt pool diameteris taken for the length l. The correlation is valid roughly for N ucvf = 50 − 1000.
• Nucleate boiling: Kutateladze [18]7/6
qnb = 3.05 × 10−11 {λl (Tw − Tsat )}10/3 P rl�σlg(ρl − ρv )

�2/3 �

pρl
σl ∆hf g µl ρv

�7/3. (2.26)• Critical heat ﬂux (CHF): Zuber [19]qchf = 0.131∆hf g ρv.
2
{ρv /σl g(ρl − ρv )}1/4(2.27)

The corresponding temperature is found by substituting qchf for qnb in the nucleate boilingcorrelation, Eq. (2.26), and solving it for Tw .
• Film boiling: Berenson [20]= Nf bf (λvf /L)(Tw − Tl ) ,qf bf= 0.425

N uf bf�(2.28)
�1/4
2

gρvf (ρl − ρvf )L3 /µvf
cpv ∆Tsat /∆h�f g P rvf,(2.29)(2.30)

where the Laplace length L and modiﬁed latent heat ∆h�f g are deﬁned by
L =∆h�f g�σl /g(ρl − ρvf ) , and= ∆hf g + 0.5cpvf ∆Tsat .

(2.31)(2.32)

Vapor ﬁlm properties ρvf , µvf , cpvf and P rvf are evaluated at the ﬁlm temperature Tvf =0.5(Tw + Tv ).
• Minimum ﬁlm boiling (MFB) heat ﬂux: Berenson [20]qmf b = 0.09

∆hf g ρvf.
2
{ρvf /σl g(ρl − ρvf )}1/4 {(ρl − ρvf )/ρvf }1/2

(2.33)

The corresponding temperature is found by substituting qmf b for qf bf in the ﬁlm boilingcorrelation, Eq. (2.28), and solving it for Tw .
9−−9

JAEA-Data/Code 2008-014• Transition boiling: Linear interpolation is made between the CHF point (Tchf , qchf ) andthe MFB point (Tmf b , qmf b ) on the boiling curve plain. (Fig. 2.2)
• Radiation: Stefan-Boltzmann law
qrad = �m σSB (Tw4 − Tl4 ) .

(2.34)

The ﬁlm boiling and radiation models are combined by
qsum = qf bf + (7/8)qrad(2.35)

based on Bromley et al. [21] and Liu-Theofanous [22].
The void eﬀect on the boiling and radiation heat transfer is considered by attenuating theheat ﬂux when the ﬂow regime departs from liquid continuous bubbly ﬂow. Also, the boiling andradiation heat transfer is cut oﬀ at very high void regime α > 0.95 to avoid numerical problems.
This is made by multiplying the heat ﬂux by the functions deﬁned as follows.
• For nucleate boiling:• For ﬁlm boiling:



 1

(0 ≤ α < 0.3)
1 − (α − 0.3)/0.2 (0.3 ≤ α ≤ 0.5)
fα,nb =
 0(0.5 < α ≤ 1)

fα,f b = 1{(0.95 − α)/0.65}n 0(2.36).

(0 ≤ α < 0.3)
(0.3 ≤ α ≤ 0.95)(0.95 < α ≤ 1),

(2.37)(0 ≤ α < 0.3)
(0.3 ≤ α ≤ 0.95)(0.95 < α ≤ 1),

(2.38)

where n is a constant. We set n = 1 at present.
• For radiation:
fα,rd =


 1{(0.95 − α)/0.65}n 0where n is a constant. We set n = 1 at present.

Pool surface temperature The heat transfer at the pool surface depends on the pool surface
temperature. Oxide melts, a typical molten core material, has a relatively low thermal conductivity and the surface temperature can be signiﬁcantly lower than the average temperature. For
the melt pool, this aspect is considered by a simple model, based on an assumption of quadratic
temperature proﬁle across the melt pool depth. The surface temperature is expressed by

Tsf =

hP H
6λP
hP H
1+
6λP

Tav + Tc

(2.39)

where Tsf , Tav and Tc denote the pool surface, pool average and coolant temperatures, respectively; H denotes the heat transfer coeﬃcient on the pool surface. Since this model is quite crude,
the result is checked to ensure Tc ≤ Tsf ≤ Tav . If a strange value arises, Tsf = 0.5(Tav + Tc ) is
imposed.
− 10 −

10

JAEA-Data/Code 2008-014

qchf
Transition
BoilingHeat flux

Nucleate
Boiling

Film Boiling

Convection
q

Tsat

Tchf

mfb

Tmfb

Temperature of melt surface

Figure 2.2: Construction of heat transfer model with available heat transfer correlations invarious modes according to the boiling curve.
Mass addition When the jet leading edge reaches the bottom of the calculation domain, themass and energy of the jet is given to the central node of the pool. The momentum of the jet
is given to the pool by adding the stagnant impulse due to jet impingement on the central nodeof the pool. Modeling for the merge of settled particles will be described later in Section 2.2.3.
Freezing of the pool Freezing and melting of the pool is decided according to the averagetemperature of the whole body of the pool for simplicity. When the average temperature islower than the melting point, the pool is considered frozen (velocity is set 0 everywhere).

2.2.3

Melt Particles

Concept and basic equations
A concept of grouped particle was developed for two-dimensional cylindrical geometry. It enables
a practical use of the code for large systems such as plant scale simulation, in which the numberof the real debris particles could be too large to be simulated directly. The grouped particle
model has a scalability according to computational resources available, by choosing the size ofparticle groups adequately.
A “particle group” is a group of Np particles which are assumed to have the same properties
and occupy a ﬁnite space, 2rx by 2rz on the x-z plain. (see Fig. 1.2) This assumption of ﬁnite
size of groups is required to move the melt particle volume, force and heat source smoothly
among the two-phase ﬂow cells. If a group is assumed to be a point, the source terms for thetwo-phase ﬂow model would jump from a cell to another suddenly and might cause numericalproblems. Other attributes of the group, i.e. the realistic shape of the group, rotation, diﬀusive
feature etc., were stripped oﬀ to simplify the situation and to make this method feasible.
The motion and heat balance of a representative particle in a group are tracked by Lagrangian− 11 −

11

JAEA-Data/Code 2008-014

equations:∂v p∇pa F hy
=g−
+
∂t
ρp
mp

and

∂xp
= v p for motion, and
∂t
∂ep
q
=
for energy conservation,∂t
mp(2.40)(2.41)(2.42)

where v p , xp , mp and ep denote particle velocity (vpx , vpz ), position (xp , zp ), mass of a particle
and internal energy per unit mass, respectively. The hydrodynamic drag force F hy and heatexchange rate q are given by constitutive equations.
This method requires other rules to handle the generation and interferences of particle groups
as described in the followings.
Rules for generation of particle groups
When particles are generated on the jet surface, the mass, temperature, size and velocity ofparticles are given in each jet cell at each time step. However, releasing a particle group fromevery jet cell at every time step makes too many groups. Besides, the mass or number of particles
generated in one time step at one cell is usually very few. This is a matter of interfacing between
Eulerian and Lagrangian forms of discretization. To workaround this problem, we introduced
a concept of “pre-particle group” as a buﬀer. (Fig. 2.3) The particles just after the birth are
stored in pre-particle groups. They are released as real particle groups after they mature to haveenough mass or number of particles. Each pre-particle group is attached to each jet cell and has
the same vertical size as the jet cell. (2rz = ∆zjet ) The horizontal size of the pre-particle groupdevelops as the particles inside travels outward. Pre-particle groups do not exchange momentum
and heat with the coolant until they are released as real-particle groups. The criteria of release
are deﬁned as follows.
1. Total mass of pre-particle group is larger than 0, and
2. one of the conditions below is satisﬁed;• the number of particles Np > Npcr ,• x-direction size rx > rxcr (≡ ∆xmin /4),
• x outer edge position exceeds the ﬁrst cell boundary, or• the age exceeds a time limit (nhist time steps),
where Npcr and nhist are the criteria given by the user and ∆xmin is the minimum x-direction
grid width in the system.
When a pre-particle group is released as a real-particle group, the x-direction velocity thatwas given by Eq. (2.13) is multiplied by a stochastic function
fr = 0.02 + 0.98φr ,(2.43)

where φr is a random number distributed uniformly between 0 and 1. This stochastic function isa trial model to express stochastic nature of the initial velocity of melt droplets emitted throughturbulence surrounding the melt jet.
Rules for interference of particle groups with other particle groups or boundaries
Some rules are required to avoid the occurrence of physically unreasonable situation upon theinterference between particle groups or between a particle group and a boundary. The following
rules were deﬁned for this purpose.
− 12 −

12

JAEA-Data/Code 2008-014

Pre-particle-group
Jet

Particlegroup

zjet

2rz

2rx

Figure 2.3: Pre-particle group concept.

1. Bounce

1+

2>

pack

(=0.6)

2

1(damping factor 0.5)

2. Merge

if two groups have similar attributes(velocity, particle dia., temperature)

3. No interference

Figure 2.4: Rules for the interference between particle groups.

13 −
− 13

1. Complete
Bounce2. Damped
BounceSide wall

Central axis

JAEA-Data/Code 2008-014(damping factor 0.1)

3. Merge (melt)
or stack (solid)

BottomFigure 2.5: Rules for the interference between a particle group and boundary.
(a) Between particle groups (Fig. 2.4)
1. If the sum of the particle volume fractions of two groups in contact is more than the packing
limit αpack (= 0.6), the two groups bounce apart with a damping factor Cdamppar (= 0.5),
2. else if the two groups in contact have similar attributes: temperature, particle size andvelocity, matching within 20%, they merge into one group,
3. otherwise, the two groups have no interference, but just go through each other.
(b) With boundaries (Fig. 2.5)
1. If a group reaches the central boundary, it is bounced back inside.
2. If a group reaches the top or the side boundary, it is bounced back inside with a damping
factor Cdampwal (= 0.1).
3. If a group reaches the bottom, it merges into the pool if the pool or the particle groupis molten, otherwise the group piles up on the frozen pool making a debris bed. When agroup joins into the debris bed, the group is vertically collapsed so that its particle volume
fraction becomes the packing limit or the group height becomes the same as the particle
size (the height of a single particle layer).
Constitutive models
Hydrodynamic dragThe drag is evaluated by a correlation for a sphere subjected to ﬂow [23]:π1F hy = Cf p f Dp2 ρvr v r ,42��
f = max 24/Rep , 18.5/Re0.6
p , 0.44 ,Rep = ρDp vr /µ ,− 14 −

14(2.44)(2.45)(2.46)

JAEA-Data/Code 2008-014

where Cf p is a correction parameter to consider deformation and irregular shape of particles.
Variables Dp , v r , vr , ρ and µ are the melt particle diameter, relative velocity between the melt
particle and the coolant and its absolute value, density and viscosity of the coolant, respectively.
The drag force F hy is calculated for the liquid and gas components of the coolant separately,
and the eﬀects are combined with the method described by Eqs. (2.1) and (2.2).
Heat transfer The surface heat transfer coeﬃcients are evaluated by a set of correlations
for radiation, ﬁlm boiling, nucleate boiling and convection around a sphere as described in thefollowing. Those various modes of heat transfer are combined by considering the boiling regimes
as shown in Fig. 2.2. Boiling is considered when the temperature of the heat transfer surface ishigher than the saturation temperature of the coolant. Convection heat transfer to the coolant
liquid and/or gas is considered if it is suitable.
• Convection: forced convection around a sphere [24]qcvp = N ucvp (λ/Dp )(Tw − T∞ ) ,(2.47)

Rep = ρp Dp vr /µ ,(2.49)

N ucvp =1/3
2 + 0.6Re1/2,
p Pr(2.48)

where N ucvp denotes the average Nusselt number for a particle of diameter Dp exposed to
ﬂow with relative velocity vr ; Rep and P r denote Reynolds and Prandtl numbers.
• Nucleate boiling: Kutateladze [18], Eq.(2.26).
• Critical heat ﬂux (CHF): Zuber [19], Eq.(2.27).
• Film boiling: Liu-Theofanous [22]qf bp = N uf bp (λvf /Dp )(Tw − Tl ) ,N uf bp =Ff
N uf(2.51)

2vrl/gDp ,(2.53)= 1 − 0.2/{1 + (F rFr ==N us =(2.50){N u5p + (Ff N uf )5 }1/5 ,1/2− 1)} ,

(2.52)

1/2
N us + 0.072Re0.77
P rl (µl /µvf )(Sc� /Sp� ) ,l
1/2
0.5Rel (µl /µvf )(KR4 /Sp� )1/4 ,(2.54)(2.55)

Rel = ρl Dp vrl /µl ,(2.56)

� 1/4

N up /(1 + 2/N up ) = kc(Ar/Sp )

Mc

1/4

(2.57),

�� 2

M c = E /{(1 + E/Sp P rl )(RP rl Sp ) } ,3

E = (A + CB
A =

B =
C =

)

1/2 1/3

�3+ (A − CB
�

Sc

�

=

(2.59)

�2(2.60)

�2(2.61)(2.62)(2.63)

R = (µvf ρvf /µl ρl )D� =+ (1/3)Sc ,)(1/27)Sc + (1/3)R Sp P rl Sc + (1/4)R Sp P rl2 ,−(4/27)Sc�2 + (2/3)Sp� P rl Sc� − (32/27)R2 Sp� P rl+(1/4)Sp�2 P rl2 + (2/27)Sc�3 /R2 ,(1/2)R2 Sp� P rl ,2

K = ρl /ρvf ,Ar =(2.58)
�

1/2 1/3

1/2,

g(ρl − ρvf )Dp3 /(µ2vf /ρvf ) ,Dp /{σl /g(ρl − ρvf )}1/2 ,cpl ∆Tsub /∆h�f g P rl ,− 15 −

15(2.64)(2.65)(2.66)(2.67)

JAEA-Data/Code 2008-014

Sp� = cpvf ∆Tsup /∆h�f g P rvf ,∆h�f g(2.68)= ∆hf g + 0.5cpvf ∆Tsup ,
0.5/D�1/4


(2.69)(D� ≤ 0.14)

0.86/(1 + 0.28D� ) (0.14 < D� ≤ 1.25)
��� 2.4D /(1 + 3.0D ) (1.25 < D ≤ 6.6)


0.47D�1/4(6.6 < D� )

kc =

.

(2.70)The relative velocity vrl is the one between water and melt particles. The liquid temperature Tl is used in the superheat term, i.e. (Tw − Tl ), instead of the saturation temperature
because it is favorable for the numerical stability. The physical properties of vapor with
suﬃx vf is evaluated at the ﬁlm temperature Tvf = 0.5(Tw + Tv ). In the original correlation, the void eﬀect is considered in N us . However, we consider it by multiplying theheat ﬂux qf bp with a function as described later. Our method gives a result close to theoriginal expression.
• Minimum ﬁlm boiling (MFB) temperature: Kondo et al. [25]Tmf b = Tsat + C

�

�

27
Tcr − Tsat +
32

�

λl
λvf

�

N uc ∆Tsub,
Dp /δmin + N ur

(2.71)

where C is a constant (= 0.6); Tcr is the critical temperature of the coolant; δmin is
the minimum limit of the vapor ﬁlm thickness, practically 0.1mm. The Nusselt number
for convection around the vapor ﬁlm which envelopes the hot particle, N uc , is given by
Eq.(2.48). The Nusselt number for radiation heat transfer, N ur , is given based on the
Stefan-Boltzmann law Eq.(2.34) and by taking account of Tmf b � Tl ,3
N ur = �m σSB Tmfb (Dp /λvf ) .

(2.72)

The corresponding heat ﬂux qmf b is found by substituting Tmf b for Tw in the ﬁlm boilingcorrelation, Eq. (2.50).
• Transition boiling: linear interpolation is made between the CHF point (Tchf , qchf ) andthe MFB point (Tmf b , qmf b ) on the boiling curve plain. (Fig. 2.2)
• Radiation: Stefan-Boltzmann law, Eq.(2.34). It is combined with the ﬁlm boiling correlation by Eq. (2.35).
The void eﬀect on the boiling and radiation heat transfer is considered by attenuating theheat ﬂux when the ﬂow regime departs from liquid continuous bubbly ﬂow. Also, the boiling andradiation heat transfer is cut oﬀ at very high void regime α > 0.95 to avoid numerical problems.
The functions are the same as deﬁned for the melt pool model, Eq.(2.36)–(2.38), except thatthe power n is 0.3 for ﬁlm boiling to ﬁt the behavior of the original Liu-Theofanous correlation,and that n is left tunable through the user input for radiation because it may have a strong
inﬂuence and it is not easy to give a precise value with a simple model.
When the particles travel downward and ﬁnally settle on the ﬂoor making a debris bed, theparticle groups are vertically collapsed. In such a case, the heat transfer is degraded due to
ineﬀectiveness in the contact with the liquid coolant and in the emission of radiation. This heattransfer degradation is simulated by putting a reduction factor to the heat transfer coeﬃcient.
The reduction factor is given as an input parameter (typically ∼0.1). When the collapsed groupconsists of more than one particle layers, further reduction factor, Dp /(2rz ), is attached, withwhich groups containing more layers of particles have more ineﬀective heat transfer.

− 16 −

16

JAEA-Data/Code 2008-014

Secondary break-up of particles The secondary break-up of particles are modeled based
on total break-up time τbrk and maximum stable size of liquid droplets Dsmx evaluated by thecorrelations by Pilch & Erdman [26],Dsmx =
τbrk =
v∗ =

τ∗ =

Cd =

�

�

v ∗ −2
W ecr σd1−,v0
ρc v02D0
τ ∗ 0.5 ,v0 ���
0.75Cd τ ∗v0 �0.5,
1 + 0.75Cd �0.5 τ ∗
6(W e − 12)−0.25(12 ≤ W e < 18)




0.25

2.45(We−
12)(18 ≤ W e < 45)−0.25
14.1(W e − 12)(45 ≤ W e < 351)
, and


0.25 (351 ≤ W e < 2670)

0.766(W
e
−
12)



5.5
(2670 ≤ W e)4
 5.6 (Re < 4 × 10 )3.3 (4 × 104 ≤ Re < 105 ) , 1.4 (105 ≤ Re)(2.73)(2.74)(2.75)(2.76)(2.77)

Cd : friction factor,
D0 : original droplet diameter,
Re : Reynolds number (≡ ρc v0 D0 /µc ),v0 : original relative velocity between the droplet and coolant,v ∗ : relative velocity decrease during droplet break-up,
W e : Weber number (≡ ρc v02 D0 /σd ),W ecr : critical Weber number (≡ 12),
� : density ratio (≡ ρc /ρd ),µc : viscosity of the continuous ﬂuid (≡ 1/{α/µg + (1 − α)/µl }),ρc : density of the continuous ﬂuid (≡ αρg + (1 − α)ρl ),ρd : density of the droplet,σd : surface tension of droplet, andτ ∗ : non-dimensional break-up time.
The decrease rate of the droplet diameter is evaluated by

where

dDp
Dmed − Dp
=
,
dt
τbrk

(2.78)

where the mass median diameter of droplets after the break-up, Dmed , is simply considered to be
the half of the maximum stable diameter, i.e. Dmed = Dsmx /2. Note that the relative velocity
v0 is evaluated by using the two-phase mixture velocity for the coolant, vc = {αρg vg + (1 −α)ρl vl }/ρc .
Strictly according to the original meaning of Pilch-Erdman correlation, the relative velocity,v0 , should be the one before the break-up process begins. However, because of the diﬃcultyof ﬁnding the start and end of the break-up for each droplet (particle group) in the present
framework of the code, we use the relative velocity at the start of each time step to obtain the
diameter decrease rate due to break-up, dDp /dt, for that time step. This model is let in eﬀect
when the droplets are molten.
However, usage of this model is not recommended for practical purpose, i.e. safety evaluation,
because this model is strongly aﬀected by two-phase ﬂow velocity that sometimes shows unstablefeature in calculations. This model may cause unexpectedly small particle sizes, that signiﬁcantly
aﬀects the heat transfer, void generation, and all the phenomena calculated. So, we recommend
to turn oﬀ this model and do a parameter study on the droplet size given as an input constant.
− 17 −

17

JAEA-Data/Code 2008-014

Surface temperature of particles The surface temperature is required for the evaluation ofheat transfer and as a criterion of the droplet break-up. It is important especially for oxide melts
which have much smaller thermal conductivity than metals, and tend to have large diﬀerence
between the inside and surface temperatures.
It is not easy to fully model the temperature distribution inside the particle with consideration on the development of solid layer at the surface even if a temperature proﬁle is assumed.
Therefore, a simple approach was adopted, that is, a quadratic temperature proﬁle is assumed
inside a particle and the diﬀerence of surface and average temperatures were evaluated from thetime passed from the particle generation.
Based on the consideration detailed in Appendix B.1, the following formula is obtained.
Development of the thermal boundary layer thickness δ inside a particle surface is expressed by
�

3δ
1
δ 1−
+
4R 5

�

δ
R

�2 �

dδ = 3κdt(2.79)

where κ is the thermal diﬀusion coeﬃcient and R is the particle radius. At the ﬁrst time step (innumerical integration) after particle generation, δ is obtained by the ﬁrst order approximation√
δ = 6κ∆t.
(2.80)
After that, Eq. (2.79) is integrated numerically to get updated δ successively. The diﬀerence ofthe surface and average temperatures is obtained by
�

qδ
δ
1−
Tav − Tsf =
2λ
R

�

1δ
1
1−
+
2 R 10

�

δ
R

�2 ��

,

(2.81)

where q is the surface heat ﬂux evaluated by heat transfer models, λ is the thermal conductivityof melt.
Neglect of phase change (crust formation) does nothing until the surface temperature drops
to the melting point, though it will make separation from the reality after that. Therefore, thetiming of surface solidiﬁcation should be adequately predicted. Besides, the surface temperatureplateau near the phase change point is practically simulated through the average temperatureand the material property model that expresses the latent heat through an artiﬁcial speciﬁc heatpeak around the melting point.

2.2.4

Extension of Melt Particle Model for Explosion Process

Framework and basic equations
The explosion process is modeled based on the following concept on the primary physics included.
• When the melt particles are exposed to shock wave, i.e. an intense relative velocity of thesurrounding coolant, very ﬁne fragments are produced from the surface of the particles byhydrodynamic force.
• The ﬁne fragments have sizes of micrometer order and the heat release from them is
extremely quick, so that it vaporizes water quick enough to support high pressure topropel and escalate the shock wave. The heat release is dominated by the generation rateof ﬁne fragments.
An additional component of the melt model, “fragment group”, is attached to the particle
group model (Fig. 2.6) for implementation of the above described concept. A fragment groupis a swarm of fragments which is generated from a particle group. Fragment groups behave as
follows.
− 18 −

18

JAEA-Data/Code 2008-014

Particle group
Fragment
group
Figure 2.6: Particle group and fragment group.
• Generation: The fragment mass is separated from the “mother” particle group by hydrodynamic interaction. The mass conservation law is expressed bydmp= −ṁf , and(2.82)
dt
dMf= np ṁf(2.83)
dt
where mp is the mass of one particle in the particle group, Mf is the mass of fragmentsin the fragment group and np is the number of particles in the particle group. The massfragmentation rate for a particle, ṁf , is given by a constitutive model, Eq. (2.85).
• Heat transfer and energy conservation: The limiting process of the heat transfer is theheat conduction inside the fragment. The heat transfer coeﬃcient outside the fragment is
considered to be quite large and any steady state heat transfer correlation is not reliable inthis highly transient situation. With a given heat release rate per unit mass of fragmentsq̇f , the energy conservation of the fragments is expressed bydef
np ṁf
= ep
− q̇f
dt
Mf

(2.84)

where ef and ep are the speciﬁc internal energies of fragments and particles. The speciﬁc
internal energy of the mother particle is assumed constant during the explosion process.
• Size: Size of the fragment is given in the input as a constant, df . When the heat releasetime of the fragment is comparable to the numerical time step, the size does not aﬀect
much.
• Movement: No kinetic equation is solved for fragments. Instead, they are just assumed toaccompany the mother particle group. The primary function of the fragments is the quickheat release and the kinetic aspect is regarded not important.
Since the primary physical process relevant to the melt in the propagation phase is thegeneration of ﬁne fragments and quick heat release from them, there is no need of keeping themelt components of jet and pool. Therefore, the jet and pool components are re-casted into
particle groups with equivalent surface area and volume at the start of an explosion calculation.
Normally, the jet and pool have much less surface area to volume ratio than particles, andrelatively small contribution in the heat release in the explosion process.
Constitutive models for propagation simulation
Fine fragmentation model
is given byṁf

The mass fragmentation rate from one melt particle ṁf [kg/s]

1π= Cf rg ∗ Dp2 vr (ρc ρp )1/2 ,tb 6− 19 −

19(2.85)

JAEA-Data/Code 2008-014

t∗b

=

Bo =

�

1
: Caracharios et al. model[27], or
13.7/Bo1/4 : Yuen et al. model[28],3 CD ρc Dp vr2,16
σp(2.86)(2.87)

where t∗b is a non-dimensional time of break-up [26]; Dp is the particle diameter; vr is the relative
velocity between the particle and coolant; ρc and ρp denote the coolant and particle densities,
respectively; σp is the surface tension of the melt particle. For the drag coeﬃcient, CD , thevalue for Newtonian regime, 0.44, is used.[23] An empirical constant Cf rg is put in the equation
for tuning.
Equation (2.85) means a simple fact that the whole volume of a particle (π/6)Dp3 is strippedaway in a characteristic period of break-up
tb = t∗b�ρp D p.
ρc vr(2.88)

Caracharios et al.[27] model gives the non-dimensional time t∗b as a constant, and Yuen et al.[28]model includes the inﬂuence of Bond number Bo.
The coolant density, ρc , and the relative velocity, vr , can be evaluated in several ways. Thosein the Bond number, Bo, are evaluated simply by averaging for gas and liquid phases byρc = αρg + (1 − α)ρl ,αρg v g + (1 − α)ρl v l
vc =
and
ρc
vr = |v p − v c | .

(2.89)(2.90)(2.91)

Note that the Bond number is included with the power of 1/4 and the deﬁnition of thesevariables does not aﬀect the result strongly. Those included in the expression of ṁf , however,
have stronger inﬂuences and are given more carefully by the following three options. By default,
the liquid density and velocity is used,ρc = ρland vr = |v p − v l | ,(2.92)and the evaluated ṁf is attenuated in case of high void fractions with a factor deﬁned by
 1(0 ≤ α < 0.3)
(0.75 − α)/0.45 (0.3 ≤ α ≤ 0.75)
fα,f rg = 0(0.75 < α ≤ 1).

(2.93)The second option is to use liquid density and velocity as in Eq. (2.92) and put a factor (1 − α),and the last option is to use the average as in Eqs. (2.89)–(2.91).
There is another model parameter for the usage of this fragmentation model. It is theduration to keep the fragmentation model active after the passage of the shock front, τtriglif e .
So far, we tested the model with alumina steam explosion experiments KROTOS[29] and foundthat JASMINE reproduces the experimentally observed pressure pulse and kinetic energy fairly
well by Cf rg ∼ 0.35 with Caracharios et al.[27] model, with τtriglif e ∼ 1ms. (see Section 4.3)
Heat release modelby

Heat release rate from a unit mass of the fragment q̇f [J/kg] is given

q̇f

= Cqf

trel =

frel(Eini − Elow ) ,
trel

1 d2f
τf ,
4 κf
− 20 −

20(2.94)

JAEA-Data/Code 2008-014

where Eini and Elow denote the internal energy of fragments at the initial high and the ﬁnal lowtemperatures, κf is the thermal diﬀusion coeﬃcient of fragments. The two temperatures are
set to be the melt particle and the water temperatures, respectively. Equation (2.94) means thefraction frel of the heat possessed by the fragment is released in time trel . Based on an analysisof transient heat conduction inside a spherical fragment (see Appendix B.2 for the detail), theconstants are frel = 0.632 and τf = 0.045. An empirical constant Cqf is put for tuning.
An attenuation factor referring the void fractionf = (1 − α)0.2(2.95)

is put to Eq. (2.94) to avoid numerical diﬃculty which might be caused by depositing large
amount of heat to a cell which contains little water.
With typical oxide melt properties (κf ∼ 1 × 10−6 m2 /s) and the fragment diameter in theorder of 10 µm, the heat release time trel is in the order of micro seconds. This is practicallycomparable to the numerical time step for propagation calculations. Thus, conceptually instant release of heat from the fragments holds. Therefore, the fragment diameter df does not
signiﬁcantly aﬀect the result as long as a small enough value in the above mentioned range is
given.

2.3

Two-Phase Flow Model

2.3.1

Basic Equations

Extensions to the two-phase ﬂow model
The two-phase ﬂow code ACE3D developed by Ohnuki et al. [11], which solves a two-ﬂuid model
for steam-water or air-water system in three-dimensional coordinates, was extended to include
the following models required for the present purpose.
1. Convecting non-ﬂuid volume model to accommodate the melt model components among
the two-phase ﬂow cells
2. Momentum and heat exchange schemes between the two-phase ﬂow and the melt models
3. Non-condensible gas components to handle highly subcooled (non-condensible gas rich)conditions
Thus, the present coolant ﬂow model in JASMINE has steam, water and non-condensible
gas components. Mass and energy equations are solved individually for these components, while
momentum equations are solved for two-phases, gas mixture and water, based on an assumptionof mechanic equilibrium among gas components.
Deﬁnition of volume fractions
As usual in the conventional multi-ﬂuid modeling, we express the conservation equations withvolume fractions of each component. The relationship among the volume fractions: αs , steam;αa (a = a1 , a2 , . . . ang ), non-condensible gases; αl , water; α, total void fraction, is deﬁned by thefollowing equations.
α = αs +�

αl = 1 − ααa(2.96)(2.97)

This multi-ﬂuid modeling has a deviation from the reality that the gas components are actually
mixed together: each component occupies the total gas volume and has a partial pressure.
− 21 −

21

JAEA-Data/Code 2008-014

However, from the mass account point of view, the mass of each gas component are equivalent
between the model and the reality, as long as the ideal gas assumption is valid. This aspect is
discussed again later in Section 2.3.2
Note that the volume fractions of ﬂuids deﬁned above are based on the ﬂuid volume in cells
after subtracting the melt volume. This ﬂuid volume in a cell is denoted by V in the followings.
Conservation equations
Conservation equations for mass and energy of each component are formulated as follows:∂(V αk ρk ) + ∇ · (V αk ρk v k ) = Γk∂tfor mass, and�

(2.98)
�∂αk∂(V αk ρk ek ) + ∇ · (V αk ρk ek v k ) = −V p
+ ∇ · (αk v k )∂t∂t+ qik + qmk + qwk +�

qjk + Γk hk

for energy,(2.99)

j�=kwhere ρ, e, p and v denote the density, internal energy, pressure and velocity, respectively. The
subscript k speciﬁes a component s, l or a. The variable Γk denotes the mass generation rate ina cell. Source terms qik , qmk , qwk and qjk are the heat input from the gas-liquid interface (forwater and steam), from the melt, from the wall, and from other components, respectively.
The momentum equation is formulated by
Γ+
∂v k
1
Fmk
k
+ v k · ∇v k = − ∇p − [Kik (v k − v i )]i�=k −(v k − v i )i�=k ++ g . (2.100)∂tρk
V α k ρk
V α k ρk
The suﬃx k is either g (gas) or l (liquid water). The total void fraction α and averaged
gas properties are used for the phase g. The variables Kik and Fmk are the interface friction
coeﬃcient and the body force exerted by the melt, respectively. The phase generation term Γ+kequals Γk if Γk ≥ 0, or otherwise 0.

2.3.2

Constitutive Models

Phase equilibrium under the existence of non-condensible gases
As described above, non-condensible gases are handled in a multi-ﬂuid model, where each component is assumed to occupy a volume fraction αk , and to be compressed by the total pressure p.
While, in the physical reality, the gas components are mixed together extending to the total gas
volume fraction (void fraction) α, having lower density corresponding to the partial pressurespk (k = s, a1 , a2 , . . . ang ). This gap makes, however, no diﬀerence in the mass if the gases canbe treated as uniformly mixed ideal gases. Besides, the present method has an advantage ineasiness of coding because it is a natural extension of the two-ﬂuid model in the original ACE3D
code.
The saturation temperature of water is evaluated at the “pseudo partial pressure” of steam
deﬁned by
αs
ps = p .
(2.101)αAs a result, the saturation temperature becomes a function of αs and α as well as pressure.
Density and internal energy of vapor are also evaluated at the pseudo partial pressure of steam,and the density is converted to the value at the total pressure byρs (p) =pρs (ps ) .
ps− 22 −

22(2.102)

JAEA-Data/Code 2008-014

Evaporation and condensation
The amount of phase change in a cell is deﬁned by−(qis + qil ) + qint,
hs − hl= Ai his (Tsat − Ts ), andΓs = −Γl =(2.103)qis(2.104)

qil = Ai hil (Tsat − Tl ),(2.105)where qis and qil are heat exchange rates from the interface to steam and liquid water; qint
is the direct heat input to the interface from the melt; his and hil are interface heat transfer
coeﬃcients on the gas side and liquid side, respectively; Ai is the interface area in a cell. Notethat the saturation temperature Tsat is evaluated at the pseudo partial pressure of steam, ps .
Exchange among components
Dispersed ﬂow model The interface exchange terms are given based on the dispersed ﬂow
model in which either of gas or liquid is considered to be continuous and the other to be dispersed.
The regimes are selected according to the void fraction as described in Section 2.1, i.e. bubbly
for α < 0.3, droplet for α > 0.75 and transition in-between.
The averaging scheme for the transition zone, 0.3 ≥ α ≥ 0.75, is deﬁned by
Xtrans = (1 − fα )Xbubbly (α = 0.3) + fα Xdroplet (α = 0.75)(2.106)

with the interpolation function fα deﬁned by Eq.(2.1).
Following models were used for the exchange coeﬃcients.
• Friction coeﬃcient: dispersed ﬂow friction factor based on correlations used in TRAC-PF1
code[30]
3 αd Cd ρc |v k − v i |
,
4 dd
αk ρk

240
(Red ≤ 0.1031)


 24
0.687
(1 + 0.15Red ) (0.1031 ≤ Red < 989)
=

Re
d


0.44
(Red ≥ 989)
ρc d d Vr,=µc
W ecr σ=,ρc Vr2

Kik =

Cd

Red
dd

(2.107),(2.108)(2.109)(2.110)where suﬃxes c and d mean continuous and dispersed phases. The critical Weber number
to evaluate the diameter of the dispersed phase is set 4.0 for droplets and 7.5 for bubbles.
The relative velocity between the dispersed and continuous phases, Vr , that is speciﬁcallyused for evaluation of the friction factor, Cd , is evaluated either by directly applying therelative velocity |v k − v i |, or by force balance of buoyancy and drag on a bubble or droplet
in a steady state (at terminal velocity) given by
π 1
π 3
d (ρl − ρg )g = Cd d2d ρc Vr2 .
6 d
4 2(2.111)The latter is recommended to avoid troubles caused by numerical (not physical) disturbances in the velocities.
− 23 −

23

JAEA-Data/Code 2008-014• Heat transfer coeﬃcient: dispersed ﬂow heat transfer coeﬃcient based on correlations usedin TRAC-PF1 code[30]– Bubbly ﬂow, liquid side
hil = 0.02ρl cpl |v l − v g | ,– Bubbly ﬂow, gas side
his = 200

(2.112)λg,dd(2.113)– Droplet ﬂow, liquid side
hil = 0.02ρl cpl vcirc ,�

vcirc = min 0.5�µg
σ|v l − v g | , 1.4µg + µld d ρl�1/2 �

(2.114),(2.115)– Droplet ﬂow, gas side
1/2

hig = (2 + 0.74Red )Red =ρg dd |v l − v g |
.
µg

λg
,
dd

(2.116)(2.117)The size of the dispersed phase, dd , is determined by the same method as that used forthe friction factor, i.e. Eqs. (2.110) and (2.111). Once the size is determined, the interface
area is given by
6V αd
Ai =
.
(2.118)dd
Equilibration approach In some cases, heat transfer coeﬃcients to force the system to approach thermal equilibrium state in a certain period might be more eﬀective than the one
deduced from a geometrical model. This type of method is used for the heat exchange between
steam and non-condensible gases to keep them in close temperatures.
Heat exchange rate [W/K] is given byRsa =V (αρcp )m.
τrsa(2.119)The suﬃx m means one of two components having less internal energy deviation from theequilibrium state. With this formulation, excess heat in the component of less energy deviation
is given up to the other roughly in the order of time τrsa at temperature diﬀerence of 1 K. The
time constant is set τrsa ∼ 0.1s.
Diﬀusion of steam in gas phase Diﬀusion of steam in the gas phase is considered by an
auxiliary model. Because the basic conservation equations (Eq. (2.98)∼(2.100)) do not include
the diﬀusion term between gas components, it may happen that physically unrealistic gaps insteam concentration in the gas ﬁeld develops, which cause abrupt change of saturation temperature and aﬀect numerical stability. Thus, this diﬀusion model is implemented not to develop
such steam concentration gaps by enhancing intermixing of gas components. Diﬀusion of steam
against the mixture of non-condensible gases is expressed by∂V αs ρs+ ∇ · (V j s ) = 0∂t− 24 −

24(2.120)JAEA-Data/Code 2008-014

where j s denotes diﬀusional mass ﬂux of steam, and it is expressed by� 2�

cg
js = −
ρg�

αs
Ms Mnc D∇α�,(2.121)

where cg and ρg are total molar concentration and total mass concentration of gas mixture,
respectively; Ms and Mnc are molecular weight of steam and average molecular weight of noncondensible gases, respectively; D is the diﬀusion coeﬃcient. The concentrations and the averagemolecular weight are deﬁned by
cg =
ρg =
Mnc =

�

�

1 αs ρs � αa ρa
+
,
α Ms
Ma
�
1
(αs ρs +αa ρa ) , and
α
�
αa Ma
.
α − αs(2.122)(2.123)(2.124)The subscript a denotes non-condensible gas components, a = a1 , a2 , ..., ang . The diﬀusion
coeﬃcient is given with a large enough value, D = 0.01 − 0.1 [m2 /s], that ensures the eﬀect ofthe model.
This auxiliary model is solved outside of the solution of basic equations. Preceding thesolution of the basic equations, volume fraction of steam, αs , in each cell is modiﬁed by equations (2.120) and (2.121). Then, non-condensible gases are moved to compensate the volume(molarity) change of steam in each cell, to keep the consistency of total pressure and total gas
volume (molarity). The correction of non-condensible gas volume fractions are given by∂V αa ρa+ ∇ · (V j a ) = 0 , and∂t
ρa αaj a = −j s.
ρs α − αs(2.125)(2.126)

Exchange with melt
As described later in Section 3.1, melt calculation is done preceding the two-phase ﬂow calculation, and it gives the melt volume, heat exchanged and force exerted in each cell to the two-phaseﬂow model.
Volume The ﬂuid volume in each cell V are obtained as a result of melt calculation. The
mass and energy equations have dV /dt term derived from the transient term. The new time andthe old time values of V are passed to evaluate them.
Force The force exerted from the melt in a cell is distributed into steam, water and noncondensible gases so that all these components get a same acceleration. Though this treatment
is not precise in view of physics, it is a simple and numerically stable method that also gives aphysically reasonable result.
Heat The heat release from the melt in a cell is distributed into steam, water, steam-waterinterface and non-condensible gases. The heat deposited on the interface is used directly forevaporation. Rules for the partition of heat is as follows.
• Cells not including melt heat source (Fig. 2.7): Heat transfer between bulk water/gas isconsidered. The heat deposited/removed at the gas-liquid interface is used for evaporation/condensation. (Eq. (2.103))− 25 −

25

JAEA-Data/Code 2008-014

Water

Hil

Hig

Vapor
+
Gases

Interface(Evaporation or condensation)

Figure 2.7: Heat partition scheme for cells not including melt heat source.
• Cells including melt heat source (premixing) (Fig. 2.8):
1. Basic concept: A fraction Fhtint (∼ 0.02 is recommended at present) of the heat input
from the melt is deposited at the interface for evaporation. The rest of the heat is
distributed to water, steam and non-condensible gases so that every component hassame temperature rise. Normal gas-liquid interface heat transfer is suppressed withfactors fhigkill and fhilkill , typically 0.1 and 0.001, respectively, so that the eﬀect ofthe melt becomes dominant.
2. Water at Tsat : If water temperature reaches saturation temperature, the heat to bulk
water is assigned to the gas-liquid interface in stead.
3. No water: If there is no water (no latent heat capacity), the heat to the interface isassigned to the gas phase in stead.
• Cells including melt heat source (explosion) (Fig. 2.9): A fraction Kev (∼ 0.7 is recommended at present) of the heat from melt (ﬁne fragments) is deposited at the gas-liquidinterface for evaporation. The rest is assigned to bulk water. Normal heat transfer between the interface and bulk water/gas is suppressed, unless physically unstable conditionsof superheated water or supercooled steam, or too high gas temperature arises.
In the above scheme, the method distributing heat to water, steam and non-condensible
gases so that they have same temperature rise, is not a physically accurate one. However, we
adopted this method due to its simplicity and numerical stability, and experiences that we did
not have signiﬁcant problem with this.

− 26 −

26

JAEA-Data/Code 2008-014

1. Basic heat partition concept

Qmelt x (1-Fhtint)

Heat input
from melt

Qmelt

Water

Qmelt x Fhtint* distribute so that
all components havesame temperature rise

Hig
Vapor * normal interface
heat transfer is
x fhigkill +suppressed byGases
f**kill facters
Interface(Evaporation term)Hilx fhilkill

2. If water temperature reaches Tsat , the heat to "water" is assigned
to the "interface".
Qmelt x (1-Fhtint)

Qmelt

Water
=Tsat

Vapor
Hig
Hil
+
x
fhigkill
x fhilkill
Gases

Interface(Evaporation term)

Qmelt x Fhtint

3. If there is no water (no-evaporation heat capacity), the heat to "interface"is assigned to the "vapor+gases".
Qmelt x (1-Fhtint)

Vapor
+
Gases

Qmelt

Qmelt x Fhtint

Figure 2.8: Heat partition scheme for cells including melt heat source (premixing).

Heat input
from melt

Qmelt x (1-Kev)

Qmelt

Vapor
+
Gases

Water

Qmelt x Kev

Interface(Evaporation term)* normal interface heattransfer is allowed only
to avoid unstable conditions:- superheated water- supercooled steam
or too high gas temperature.

Figure 2.9: Heat partition scheme for cells including melt heat source (explosion).

27 −
− 27

JAEA-Data/Code 2008-014

3 Numerical Solution Method
3.1

Coupling of Melt and Two-Phase Flow Models

Coupling of the melt and two-phase models are made explicitly in time as illustrated in Fig. 3.1.
The melt model is called at the beginning of a new time step, referring the two-phase ﬂow
variables obtained in the previous step. After interpolating the two-phase ﬂow variables to themelt jet and melt pool grids, and to the location of each particle group, the three melt modelcomponents, i.e. jet, particles and pool, are calculated successively, with necessary informationexchanges among them: jet out-ﬂow to the pool, jet break-up to the particle, and particle
settlement to the pool.
The exchange terms from the melt models, i.e. melt volume, heat transfer and force, aresummed for two-phase ﬂow cells and given back to the two-phase ﬂow calculation for update to
the next time step.
Figure 3.2 shows the grid that is used for the numerical solution. Melt jet and pool grids
are accommodated in the two-phase grid on x-z plain.
In the following sections (Sections 3.2 and 3.3), solution methods for the sub-models aredescribed.

3.2

Numerical Solution Method for Melt Model

3.2.1

Melt Jet

The conservation equations, Eqs. (2.3)–(2.5), were solved by CIP (cubic-interpolated pseudoparticle) method [31], which is known as a high resolution scheme. Conceptually it is free fromnumerical diﬀusion that is caused by ﬁnite diﬀerence of the advection term. The detail of thismethod is described in Appendix C.
The grid for the melt jet is ﬁt in the z-direction grid of two-phase ﬂow model, with a mannerthat one two-phase ﬂow cell is subdivided into Nsub (∼ 5) jet cells to give the jet ﬁner resolution
than the coolant ﬂow ﬁeld. A staggered grid is adopted where scalar variables are deﬁned at thecell center and vector variables are deﬁned at the cell boundary. Figure 3.3 shows the geometryof the grid for the melt jet model.
In the CIP method, basic equations are solved simultaneously with their derivatives. A
general form of the conservation equations and their derivatives are expressed by
∂f
∂t
∂f �
∂t

∂f
= G , and
∂z
∂f �
∂vJ
+ vJ
= G� − f �
,
∂z
∂z

+ vJ

(3.1)
(3.2)

where f and G denote the advected variable and source term, respectively. For the mass, energy
and momentum equations, they are expressed as follows.
�

Mass : fi = [AJ ρJ ]i , Gi = −[ 2ΘAJ me ]i − [AJ ρJ ]i
Energy : fi = [eJ ]i , Gi = −��Momentum : fi+ 1 = [vJ ]i+ 1 ,2

2− 28 −

28

2Θ q
AJ ρJ

�

�

∂vJ
∂z

�

(3.3)

i

(3.4)
i

JAEA-Data/Code 2008-014

JAEA-Data/Code 2008-014

Time

Two-phase flow
ACE3D

Melt models
Particles

Jet

Pool

t Pressure,Volume fractions,Velocities,
Temperatures, Densities(Interpolate)(Interpolate)(Interpolate)(Calc.)
Particle gen.
(Calc.)

Flow into pool(Calc.)
Merge into pool

Melt Volume,
Momentum/Energy
Exchanges
Pressure, Volume fractions,(Calc.)

Jet Volume,
Particle position, Pool Volume,
Velocity, Temp. Velocity, Temp. Velocity, Temp.

t+ ∆ t Velocities, Temps., Densities(Cell average/sum)

Figure 3.1: Coupling algorithm of melt and two-phase ﬂow models.

Θ

Melt Jet
Grid

Two-phase
Flow Gridzx

Melt pool Grid

Figure 3.2: Grid for two-phase ﬂow, melt jet and melt pool.

– 29− –29 −

JAEA-Data/Code 2008-014

JAEA-Data/Code 2008-014

DJin, TJin
AJin, eJin, ρJin
njmax+1

njmax+1/2
njmax

vJin
njmax-1/2

i+1/2
i
i-1/2

z

njmin+1/2

njmin

njmin-1/2
njmin-1

njmin-3/2

DJout, TJout, vJout

Flow into melt pool

* In the code, cell boudnary index i-1/2, i+1/2, i+3/2 ...
are denoted by i-1, i+1, i+2 ....

Figure 3.3: Grid geometry and boundary conditions for melt jet model.

−–3030−–JAEA-Data/Code 2008-014

�

�

1 ∂pa
Gi+ 1 = −+ [Kf (va − vJ )]i+ 1 + gi+ 1
222
ρJ ∂z i+ 1

(3.5)

2

The suﬃx i or i + 12 denotes the location of the variable deﬁnition. The center of a cell i isrepresented by i, and its both ends are represented by i − 12 and i + 21 .
Because of this diﬀerence of the deﬁnition points between velocity and other variables weneed a scheme for averaging velocities. Interpolation of the velocity and its derivative is givenby the followings. It is third-order by taking advantage of having the values of derivatives.
[vJ ]i =
�

∂vJ
∂z

�

=

i

zi+ 1 − zi− 1
1
�
�
2
2
(vJ,i+ 1 + vJ,i− 1 ) −(vJ,i+) , and
1 − v
J,i− 12
2
2
2
2
8
3
1 �
�
(v
1 − v
J,i− 12 ) − (vJ,i+ 12 + vJ,i− 12 ).
2(zi+ 1 − zi− 1 ) J,i+ 2
4
2

(3.6)
(3.7)

2

In the momentum equation, the pressure gradient term was discretized as follows.
�

�

1 ∂pa
2
pa,i+1 − pa,i=ρJ ∂z i+ 1
ρi+1 − ρi zi+1 − zi

(3.8)

2

As described in Appendix C in detail, the conservation equations are solved in two steps,
i.e. advection and non-advection phases. Starting from the old time step value, f n , ﬁrst, theadvection phase is calculated with so called CIP1 scheme, and the intermediate value, f ∗ , isobtained. Then, the non-advection phase is solved by ordinary ﬁnite diﬀerence scheme and the
updated value, f n+1 , is obtained.
Boundary conditions are given as follows. (Fig. 3.3)
• At the inlet: velocity, diameter and temperature (internal energy) is given.
vJ,njmax+ 1 = vJin , DJ,njmax+1 = DJin , eJ,njmax+1 = eJin
2
Corresponding advected variables are calculated from them. Derivatives are set 0.
• At the bottom: Free ﬂow-out condition is given, i.e. vJ,njmin− 3 and fnjmin−1 (f = AJ ρJ or
2eJ ) are calculated from advection of linear interpolated proﬁles for [njmin − 32 , njmin − 21 ]and [njmin−1, njmin], respectively. Derivatives are also obtained from the linear proﬁles.

Evaluation of the ﬂow out mass needs high accuracy because a large part of the jet massis passed to the melt pool through it, and the accuracy is directly reﬂected on the total melt
mass conservation. So, the cell centered variables are 3rd order interpolated to the bottom end,
njmin − 12 , byfnjmin− 1

2

=

∆ =

1
∆ �
�
(fnjmin + fnjmin−1 ) + (fnjmin−1− fnjmin) , where
2
4
1
(z
1 − z
njmin− 12 ) .
2 njmin+ 2(3.9)(3.10)

Then, the variables passed to the pool, i.e. jet diameter DJout , velocity vJout and temperature
TJout , are evaluated by averaging the new and old time step values.

3.2.2

Melt Pool

The one-dimensional (radial direction) grid for the melt pool resides in a single layer of the twophase ﬂow grid on the bottom. The pool model uses the same x-direction grid with the two-phaseﬂow model. Figure 3.4 shows the grid geometry for the melt pool model. A conventional upwind scheme with a staggered mesh, i.e. velocities deﬁned at cell boundaries i − 21 and i + 12 , and− 31 −

31

JAEA-Data/Code 2008-014

npmin-1

npmin

x

npmax

i

i-1/2

npmin-1/2 npmin+1/2

i+1/2

npmax+1

npmax-1/2 npmax+1/2

Mirror boundary

Mirror boundary* In the code, cell boudnary index i-1/2, i+1/2, i+3/2 ...
are denoted by i-1, i+1, i+2 ....

Figure 3.4: Grid geometry and boundary conditions for melt pool model.
other variables at the cell center i, and the SIMPLE algorithm [32] were applied for numericalsolution.
The ﬁnite diﬀerence form of the conservation equations Eqs. (2.20)–(2.22) are expressed asfollows.
�1
n+1{max(vi+ 1 , 0)fin+1 + min(vi+ 1 , 0)fi+1}
2
2
xi+ 1 − xi− 1

fin+1 − fin
= −
∆t

2

2

2

2

for mass, where f = whP ρP ,n
en+1
P,i − eP,i�

= − max(vP,i , 0)

∆t

�

+

�

n+1
+ min(vi− 1 , 0)fin+1 } + [wms ]i{max(vi− 1 , 0)fi−1

−

n+1
en+1
P,i − eP,i−1+ min(vP,i , 0))xi − xi−1

n+1
en+1
P,i+1 − eP,ixi+1 − xi�(ems − eP )max(ms , 0)q+
h P ρPhP ρP
i(3.11)

�(3.12)

for energy, and
n+1
n
vP,i+
1 − v
P,i+ 12∆t2

= − max(vP,i+ 1 , 0)
−

+

2

1

ρP,i+ 1

�

2

�

n+1
n+1
vP,i+
1 − v
P,i− 1
2

2

xi+ 1 − xi− 12+ min(vP,i+ 1 , 0)2pa,i+1 − pa,i g [hP ρP ]i+1 − [hP ρP ]i+
xi+1 − xi2xi+1 − xi

Kf a (va − vP ) − Kf w vP +

n+1
n+1
vP,i+
3 − v
P,i+ 1

2

�

2
2

�

2

xi+ 3 − xi+ 1(vms − vP ) max(ms , 0)
hP ρP
i+ 1

2




(3.13)

2for momentum.
The mass and momentum equations are solved by the SIMPLE method, ﬁrst, to obtain the
updated velocity and mass ﬁelds. Then, the energy equation is solved.
At the center of the melt pool where the melt jet drops into, the pressure by jet impingement,
pJin =

1 ρJ vJ2 DJ2
2 4x21(3.14)

is added to the pool surface pressure pa , where x1 is the size of the central cell for the melt pool.

− 32 −

32

JAEA-Data/Code 2008-014

When the melt jet comes into the center of the melt pool, it gives sudden increase of thepool height and velocity, and easily causes numerical instability. To suppress this problem, an
artiﬁcial viscosity termµa ∂ 2 vP(3.15)
ρP ∂x2is added to the right hand side of the momentum equation. The artiﬁcial viscosity by Fletcher[33]��� ∂vP ��µa(b0 ∆x)2 ��
=
∂x �

ρP
0

∂vP
vP < 0, or∂xotherwiseif(3.16)

is used with the constant b0 ∼ 10. The discretized form of the artiﬁcial viscosity term, when it
is not zero, is given by�µa ∂ 2 vPρP ∂x2

�

i+ 21

�
�
� ∂vP �2�
�= b0 (xi+1 − xi ) �
∂x �

 n+1
n+1
 vP,i+ 3 − vP,i+ 1
2

DW 

2

xi+ 3 − xi+ 1
2

2

−



n+1
n+1
vP,i+

1 − v
P,i− 1
2

2

xi+ 1 − xi− 1 
2,(3.17)

2

where the velocity gradient with suﬃx DW means down-wind ﬁnite diﬀerence, given by� ��
vP,i+ 3 − vP,i+ 1 ��

�


2
2 �
�


� x 3 −x 1 �
�
�


�
� ∂vP �i+ 2
i+ 2 �����=� ∂x �
�v
�

1 �
1 − v

DW
�
P,i−
P,i+


2
2 �
�


 �� x 1 − x 1 ��
i+ 2
i− 2

if vP,i+ 1 > 0 , or2

(3.18)
otherwise.

Boundary conditions at both ends of the grid are given as mirror conditions. (Fig. 3.4) The
velocity is set 0 at npmin − 12 and npmax + 21 . Cell centered variables, i.e. pool hight andinternal energy, in the boundary cells are given by
wnpmin−1
fnpmin−1 = fnpmin,(3.19)
wnpminwnpmax+1
fnpmax+1 = fnpmax.
(3.20)
wnpmax3.2.3

Melt Particle

The kinetic and energy conservation equations for a particle in a particle group, Eq. (2.40)–(2.42), are solved by a simple explicit scheme for time advancement,− v np
v n+1
p
∆t
n+1
xp − xnp
∆t
en+1
− enp
p
∆t

= g−

∇pa F hy+,ρp
mp(3.21)= v n+1, andp(3.22)

q.
mp(3.23)=

Additional equations for the explosion process, Eqs. (2.82)–(2.84), are also solved in the same
way.
mn+1
− mnp
p
∆t
Mfn+1 − Mfn
∆t

en+1
− enf
f
∆t

= −ṁf(3.24)= np ṁf , and(3.25)= ep− 33 −

33

np ṁf− q̇f .
Mf

(3.26)

JAEA-Data/Code 2008-014

Evaluation of the variables ∇pa , F hy , q, ṁf and q̇f needs interpolation of two-phase ﬂow variables to the location of particle groups. For this purpose, the indexes, i and k, of two-phase ﬂow
cells in which particle groups exist are searched by bi-section algorithm, and the interpolationof two-phase ﬂow variables is done at the start of the calculation of the melt particle model.
During the calculation, the melt volume, force and heat exchanged with the two-phase ﬂow
are summed and stored for every two-phase ﬂow cell and passed to the two-phase ﬂow model.

3.3

Numerical Solution Method for Two-Phase Flow

3.3.1

Overall Scope

Figure 3.5 shows the geometry of the grid for two-phase ﬂow model. The ACE3D code isoriginally a 3D two-phase ﬂow code, capable of solving a 3D basic equation set on Cartesian or
cylindrical coordinate. JASMINE uses the code, however, with limitation of the geometry to
2D cylindrical coordinate. In the followings, equations are expressed in 3D Cartesian system forbrevity.
Numerical solutions are obtained by a conventional ﬁnite diﬀerence method, i.e. semi-implicit
scheme for time advancement, staggered grid and upwind scheme for spatial discretization.
Practically, the following 10 + 2ng discretized equations are solved.
• Mass: steam, steam+water, non-condensible gases (2 + ng )• Energy: steam, steam+water, non-condensible gases (2 + ng )• Momentum: gas (steam and non-condensible gases), water (2 phases × 3 directions = 6)
Independent variables solved for are the total pressure, 2 + ng temperatures, 1 + ng volume
fractions (water volume fraction is known if others are obtained), and 3 velocity components for2 phases—10 + 2ng variables in total.
The mass and energy equations of the water-steam mixture are included so that at least oneof the two equations is kept signiﬁcant when either of water or steam is vanished. Also, a variable
conversion, T → αT , is done in constructing the innermost system equation set which is solved
by matrix inversion and Newtonian iteration, not to have singularity for vanished components.
The overall solution method is brieﬂy as follows. The velocities at the new time step canbe expressed by pressures in the adjacent cells by transforming momentum equations. By using
these expressions, velocities are eliminated from mass and energy equations. The mass andenergy equations are, then, linearized for the variance of independent variables, i.e. pressure,volume fractions of steam and non-condensible gases, and αT of all the components. Thus, asystem of linear equations for Newtonian iteration is obtained. A set of pressure equation, that
only includes pressure at adjacent cells, is picked from the linear system. In solving the linear
system, pressure correction is obtained ﬁrst by the pressure equation, then correction of othervariables are obtained by substituting the pressure correction into the rest of the equations in
the linear system. This process is iterated until residual becomes small enough, then all thevariables except velocities at the new time step are obtained. Finally, velocities are calculated
from pressure distribution by the momentum equations.

3.3.2

Finite Diﬀerence Form of Basic Equations

In the following expressions, suﬃx i, j and k denote x(r), y(θ) and z direction indexes of thecoordinates, respectively. When the coordinate index is omitted, that means the location where
the equation is deﬁned. (e.g. pi,j,k → p, ui+ 1 ,j,k → u) Suﬃx m denotes a component either s2(steam), l (water) or a(= a1 , a2 , ..., ang ) (non-condensible gas) for the mass and energy equations.
− 34 −

34

JAEA-Data/Code 2008-014

C.L.
nz+1
nz+1/2nz
nz-1/2

nz-1
.
.
.

k+1/2

k
k-1/2

.
.
.

3

z

2
3/2

1
1/2

k=0

i=0

1/2

1

2

......

3/2

i-1/2

i

i+1/2

nx+1
...... nx-1 nx
nx-1/2 nx+1/2

x* Only x-z plain is shown. JASMINE uses only j=1 layer iny(theta)-direction of 3D field supported by ACE3D.
* In the code, cell boudnary index i-1/2, i+1/2, i+3/2 ...
are denoted by i-1, i+1, i+2 ....

Figure 3.5: Grid geometry for two-phase ﬂow model.

35 −
− 35

JAEA-Data/Code 2008-014

In the momentum equation, m denotes either of g (gas mixture) or l (water). Superscripts n
and n + 1 denote the old and new time step values. An expression v �X� means that the variable
X is taken at the upwind side of the velocity deﬁnition point for v. For example, if ui+ 1 > 0,2
u �X� = ui+ 1 Xi .
2

Mass at (i, j, k):dV(αm ρm )ndt− �αm ρm �ni− 1 ,j,k un+1
�αm ρm �ni+ 1 ,j,k un+1m,i+ 1 ,j,km,i− 1 ,j,k

Mm = V {(αm ρm )n+1 − (αm ρm )n } + ∆t
+∆tV

2

2

2

2

+∆tV

2

n+1
n+1
�αm ρm �ni,j+ 1 ,k vm,i,j+− �αm ρm �ni,j− 1 ,k vm,i,j−
1
1
,k
,k
2

2

2

2

yj+ 1 − yj− 1
2

+∆tV

2

xi+ 1 − xi− 1

2

n
n+1
n+1
�αm ρm �ni,j,k+ 1 wm,i,j,k+
1 − �αm ρm �i,j,k− 1 w
m,i,j,k− 1
2

2

2

2

zk+ 1 − zk− 12−∆tΓn+1=0m

2(3.27)

Energy at (i, j, k):dV(αm ρm em )ndt
�αm ρm em �ni+ 1 ,j,k un+1
− �αm ρm em �ni− 1 ,j,k un+1m,i+ 1 ,j,km,i− 1 ,j,k

Em = V {(αm ρm em )n+1 − (αm ρm em )n } + ∆t
+∆tV

2

2

2

2

+∆tV

2

n+1
n+1
�αm ρm em �ni,j+ 1 ,k vm,i,j+− �αm ρm em �ni,j− 1 ,k vm,i,j−
1
1
,k
,k
2

2

2

2

+∆tV p 2

2

2

zk+ 1 − zk− 12

n+1 − αn
αm
m
+
∆t

2

�αm �ni+ 1 ,j,k un+1
− �αm �ni− 1 ,j,k un+1m,i+ 1 ,j,km,i− 1 ,j,k
2

2

2

2

xi+ 1 − xi− 1
2

2

n+1
n+1
�αm �ni,j+ 1 ,k vm,i,j+− �αm �ni,j− 1 ,k vm,i,j−
1
1
,k
,k
2
2
2

2

yj+ 1 − yj− 1
2

+

2

n
n+1
n+1
�αm ρm em �ni,j,k+ 1 wm,i,j,k+
1 − �αm ρm em �i,j,k− 1 w
m,i,j,k− 1



+

2

yj+ 1 − yj− 1
2

+∆tV

2

xi+ 1 − xi− 1

2


n
n+1
n+1
�αm �ni,j,k+ 1 wm,i,j,k+
1 − �αm �i,j,k− 1 w
m,i,j,k− 1
2
2
2

2

zk+ 1 − zk− 1
2




n+1
n
n
−∆t qim
+ qm
+ qw
+



2

�

n+1
n+1
qjm
+ Γn+1m hm

j�=mMomentum, x-direction at (i + 21 , j, k):nun+1
m − um
∆t

+
+

unm

�

∂um
∂x

�

i+ 12 ,j,k

n
+ vm

�

�

∂um
∂y

�





i+ 12 ,j,k

�



=0

n
+ wm

(3.28)

�

∂um
∂z

�

i+ 21 ,j,kn+1
Γ+
1 pn+1i+1 − pim
n+1
+ Kim +(un+1)i�=mm − uiρm xi+1 − xiV αm ρm− 36 −

36

JAEA-Data/Code 2008-014−

Fmk− gx = 0V α m ρm

(3.29)The momentum equations constitute simultaneous equations of ug and ul . By solving them
analytically, we obtain expressions of ug and ul at the new time step in terms of pressure at new
time step,
un+1
g,i+ 1n= Cgx,i+ 1 + Dgx,i+ 1 (pn+1i+1 − pi ) ,(3.30)

un+1
l,i+ 1n= Clx,i+ 1 + Dlx,i+ 1 (pn+1i+1 − pi ) ,(3.31)

2

2

2

2

2

2

where Cgx , Dgx , Clx , Dlx are constants (including only old time step variables). In the sameway, expression of every component of gas and liquid velocities in terms of pressure is obtained.
n+1
vg,j+
1

n
= Cgy,j+ 1 + Dgy,j+ 1 (pn+1
j+1 − pj ) ,(3.32)

n+1
vl,j+
1

n
= Cly,j+ 1 + Dly,j+ 1 (pn+1
j+1 − pj ) ,(3.33)

n+1
wg,k+
1

n
= Cgz,k+ 1 + Dgz,k+ 1 (pn+1
k+1 − pk ) ,(3.34)

n+1
wl,k+
1

n
= Clz,k+ 1 + Dlz,k+ 1 (pn+1
k+1 − pk ) .

(3.35)

2

2

2

2

2

2

2

2

2

2

2

2

These equations are substituted into mass and energy equations to eliminate velocities at the
new time step.
After eliminating velocities, Eqs. (3.27) and (3.28) are expressed in the following form,
n+1
n+1
Mm (pn+1 , pn+1
adj , αm , Tm )= V (αm ρm )n+1 + V

�
adj

n+1
� m,adj (pn+1m) − ∆tΓn+1m
adj − pdV�m = 0∆t − V M
dtn+1
n+1
Em (pn+1 , pn+1
adj , αm , Tm )+(αm ρm )n= V (αm ρm em )n+1 + V

�
adj

(1)

(3.36)(2)

n+1(e�m,adj + e�m,adj pn+1 )(pn+1)
adj − pn+1
� m − ∆t(q n+1 + Γn+1 hn+1 ) − ∆t+V pn+1 αm+ V pn+1 Dmmim�m + (αm ρm em )n dV ∆t = 0 .
−V E
dt

�

n+1
qjm

j�=m

(3.37)The suﬃx adj means adjacent cells; m denotes either of s (steam), l (water) or ai (i-th element
�m , e�(n) , D
� m and E
�m are constants including
� m,adj , M
of non-condensible gases). Symbols mm,adj
only old time values, and are derived from the advection terms. The summation in the second
term of the right hand side of mass equation, Eq. (3.36), is expanded as
�
adj

n+1
� m,adj (pn+1m) =adj − p

+
+
+

∆t
n+1
n+1
)
D
1 �αm ρm �
i+ 21 (pi+1 − p
xi+ 1 − xi− 1 mx,i+ 2
2

2

∆t
n+1
Dmx,i− 1 �αm ρm �i− 1 (pn+1)
i−1 − p2
2
xi+ 1 − xi− 1
2

2

2

2

2

2

∆t
n+1
n+1
D
)
1 �αm ρm �
j+ 21 (pj+1 − p
yj+ 1 − yj− 1 my,j+ 2
∆t
n+1
n+1
)
D
1 �αm ρm �
j− 21 (pj−1 − p
yj+ 1 − yj− 1 my,j− 2− 37 −

37

JAEA-Data/Code 2008-014

+

∆t
n+1
n+1
D
)
1 �αm ρm �
k+ 21 (pk+1 − p
zk+ 1 − zk− 1 mz,k+ 2
2

+

2

∆t
n+1
n+1
D
).
1 �αm ρm �
k− 21 (pk−1 − p
zk+ 1 − zk− 1 mz,k− 22(3.38)

2

�m is expressed by
The constant M

�m = αm ρmM
�αm ρm �i+ 1 Cmx,i+ 1 − �αm ρm �i− 1 Cmx,i− 1
2
2
2
2
+ ∆t
xi+ 1 − xi− 12(3.39)2+ ∆t

�αm ρm �j+ 1 Cmy,j+ 1 − �αm ρm �j− 1 Cmy,j− 1+ ∆t

�αm ρm �k+ 1 Cmz,k+ 1 − �αm ρm �k− 1 Cmz,k− 1

2

2

2

2

yj+ 1 − yj− 1
2

2

2

2

2

2

zk+ 1 − zk− 1
2.

(3.40)

2

The summation in the second term of the right hand side of energy equation, Eq. (3.37), is
expanded as
�
adj

(1)

(2)

n+1
(e�m,adj + e�m,adj pn+1 )(pn+1)=
adj − p

∆t
n+1
n+1
D
)(pn+1)1 (�αm ρm em �
i+1 − pi+ 21 + �αm �i+ 21 p
xi+ 1 − xi− 1 mx,i+ 2
2

2

2

2

2

2

2

2

2

2

2

2

∆t
n+1
n+1
D
)(pn+1)1 (�αm ρm em �i−1 − pi− 21 + �αm �i− 21 p
xi+ 1 − xi− 1 mx,i− 2

∆t
n+1
n+1
D
)(pn+1)1 (�αm ρm em �
j+1 − pj+ 21 + �αm �j+ 21 p
yj+ 1 − yj− 1 my,j+ 2
∆t
n+1
n+1
D
)(pn+1)1 (�αm ρm em �
j−1 − pj− 21 + �αm �j− 21 p
yj+ 1 − yj− 1 my,j− 2

∆t
n+1
n+1
D
)(pn+1)1 (�αm ρm em �
k+ 21 + �αm �k+ 21 p
k+1 − pzk+ 1 − zk− 1 mz,k+ 2

∆t
n+1
n+1
D
)(pn+1).
1 (�αm ρm em �
k− 21 + �αm �k− 21 p
k−1 − pzk+ 1 − zk− 1 mz,k− 2(3.41)

� m and E
�m are expressed by
The constant D� m = −αm
D

+

+
+

∆t
(�αm �i+ 1 Cmx,i+ 1 − �αm �i− 1 Cmx,i− 1 )22
22
xi+ 1 − xi− 1
2

2

2

2

2

2

∆t
(�αm �j+ 1 Cmy,j+ 1 − �αm �j− 1 Cmy,j− 1 )2
2
22
yj+ 1 − yj− 1∆t(�αm �k+ 1 Cmz,k+ 1 − �αm �k− 1 Cmz,k− 1 ) , and22
2
2
zk+ 1 − zk− 1

�m = αm ρm em + ∆tEqmelt,mV− 38 −

38(3.42)

JAEA-Data/Code 2008-014

−
−
−

∆t
(�αm ρm em �i+ 1 Cmx,i+ 1 − �αm ρm em �i− 1 Cmx,i− 1 )22
22
xi+ 1 − xi− 1
2

2

2

2

∆t
(�αm ρm em �j+ 1 Cmy,j+ 1 − �αm ρm em �j− 1 Cmy,j− 1 )2
2
22
yj+ 1 − yj− 1∆t(�αm ρm em �k+ 1 Cmz,k+ 1 − �αm ρm em �k− 1 Cmz,k− 1 ) .
2
2
2
2
zk+ 1 − zk− 1
2(3.43)

2

Now, the mass and energy equations at a cell (i, j, k) include pressure, volume fractions andtemperatures at the local cell and pressure at adjacent cells (i ± 1, j ± 1 and k ± 1) as unknown
variables. In other words, the equations for a cell is connected with those of adjacent cells only
through pressure, and the dependence is linear. The dependence on the variables in the local
cell is, however, non-linear, through the phase change term and physical properties.
So, the mass and energy equations are solved by Newtonian iteration as described later in
Section 3.3.4.

3.3.3

Boundary Condition Settings

Boundary condition for velocities
Boundary conditions for velocities are categorized and given by the followings. Expressions are
given at the cell (i, 1, 1) as an example (See Fig. 3.5).
Normal velocity The normal velocity is deﬁned on the wall surface. Therefore, the conditionis given by setting it adequately through Eqs. (3.30)∼(3.35).
• Fill: given velocity at the wall surface; wi,1, 1 = win , that is satisﬁed by setting Cmz, 1 = win
22
and Dmz, 1 = 0 (m = g or l)2• Break: given pressure at the boundary cell (pi,1,0 ) (outside of the system), solution of the
momentum equation to obtain wi,1, 1
2

Tangential velocity Since the tangential velocity is not deﬁned just on the wall surface, thevelocity in the boundary cell ui+ 1 ,1,0 (outside of the system) is given in terms of that inside of2
the system ui+ 1 ,1,1 so that the condition on the surface is satisﬁed.
2• Slip: zero shear force at wall surface, i.e. derivative of the tangential velocity in the normal
direction (∂u/∂z) = 0;
ui+ 1 ,1,0 = ui+ 1 ,1,1
2

2

• Non-slip: zero tangential velocity at wall surface,∆z1 ui+ 1 ,1,0 + ∆z0 ui+ 1 ,1,1
2
2
ui+ 1 ,1, 1 �= 0 (∆z0 = z 1 − z− 1 , ∆z1 = z 3 − z 1 );
2
2
2
2
2
2
∆z0 + ∆z1
∆z0
therefore, ui+ 1 ,1,0 = −u 1
2∆z1 i+ 2 ,1,1
Velocity boundary conditions on other wall surfaces are given in the same manner.
Boundary condition for cell-centered variables
Cell-centered variables, pressure, volume fractions and temperatures, should be given at theboundary with “ﬁll” condition, where inlet velocity is given. Values of those variables are set atthe boundary cell so that they are used for the ﬂuid ﬂowing into the system. As an example,
when k = 1/2 boundary is a “ﬁll” boundary, pi,1,0 , αm,i,1,0 and Tm,i,1,0 should be given (suﬃxm denotes a component).
− 39 −

39

JAEA-Data/Code 2008-014

3.3.4

Newtonian Iteration

General method
The solution of a set of non-linear equations
Fk (x1 , x2 , ..., xi , ..., xn ) = 0(3.44)

where k = 1, ..., n is obtained iteratively as follows. Assume that the unknown variables xi has
initial values that give non-zero values of Fk (residual errors). The desired correction of Fk isδFk = −Fk .

(3.45)On the other hand, the variation of Fk is approximated in terms of variation of xi by Taylor
expansion
∂Fk
∂Fk
∂Fk
δx1 +
δx2 + · · · +δxn .
(3.46)
δFk �
∂x1
∂x2
∂xn
Equations (3.45) and (3.46) gives a linear equation set to give correction values of xi ,








∂F1 /∂x1
 ∂F2 /∂x1
J =
..

.


(3.47)

Fnδxn

where J denote Jacobian matrix,



δx1
F1
 .. 
 .. 
J  .  = − .  ,∂F1 /∂x2∂F2 /∂x2

· · · ∂F1 /∂xn
· · · ∂F2 /∂xn 

 .
..

.


(3.48)∂Fn /∂x1 ∂Fn /∂x2 · · · ∂Fn /∂xn

Obtaining the correction values δxi by solving Eq. (3.47) and doing correction byxnext= xcurrent
+ δxi
i
i

(3.49)

until the residual errors Fk become close enough to zero makes conversion of the solution(x1 , · · · , xn ).
Application to the present case
The above method is applied for the mass and energy equations. However, the pressure ﬁeld
is treated diﬀerently as follows because it includes connection to adjacent cells. In the present
case, the equation for correction, Eq. (3.47), for a cell (i, j, k) is given by


δp
 δα
s


 δβs

 δβl


J  δαa1
 δβ

a1

..

.




Msl
 M

s




 Es



 Esl





 = −  M a1
 E



a1


..


.





 δαang 

δβang







∂Msl /∂Padj
 ∂M /∂P

s
adj




 ∂Es /∂Padj





 �  ∂Esl /∂Padj
 ∂Ma /∂Padj

1

−


 adj  ∂Ea1 /∂Padj


..


.











 δPadj ,







 ∂Mang /∂Padj  Mang Eang

∂Eang /∂Padj− 40 −

40(3.50)

JAEA-Data/Code 2008-014

where suﬃx sl means sum of the steam and water equations, i.e. Msl = Ms + Ml , Esl = Es + El ;βm ≡ αm Tm is used as an independent variable instead of Tm to avoid trouble in getting Tm
for missing components (αm = 0) . The Jacobian matrix J is, for the present case, a matrix ofpartial derivatives of the variables in the ﬁrst term in the right hand side by p, αs , βs , βl , αa1 ,
βa1 , . . . , αang and βang .
Note that the pressures in adjacent cells are kept out of the Jacobian matrix and they aresolved separately together with all the pressures in the system. By matrix inversion, Eq. (3.50)
becomes∂Msl /∂padj
Msl
δp ∂M /∂p M  δα 
s
s 
s 
adj 









 ∂Es /∂padj  Es 
 δβs  ∂Esl /∂padj  Esl  δβl 






�


 δαa −1
−1  Ma1  − J
1  = −J(3.51)
 ∂Ma1 /∂padj  δpadj .
 ∂E /∂p


 E

 δβ


a1
a1 
a1 
adj 
adj 





..
.. 
..






.
. 
.







 δαang This equation can be rewritten by
δp
δαs
δβs
δβl
..
.

∂Eang /∂padj

Eang

δβang




 ∂Mang /∂padj  Mang 







f1
f2
f3
f4
..
.



g1
g2
g3
g4
..
.
























 �








 = −− δpadj .









 adj 













 δαang  f4+2ng −1 
 g4+2ng −1 

δβang

f4+2ng(3.52)

g4+2ng

The ﬁrst line of this equation

δp = f1 −�

g1 δpadj(3.53)

adj

includes only pressures of the cell and adjacent cells as unknown variables, and can be solved asa system wide pressure equation set. After obtaining the correction values of all the pressuresin the system by Eq. (3.53), correction values for other variables are obtained cell by cell, byapplying the rest of the lines in Eq. (3.52),


δαs
δβs
δβl
..
.





f2
f3
f4
..
.

g2
g3
g4
..
.
















 �









 δpadj .
− = −








 adj 









 g4+2ng −1 
 f4+2ng −1  δαang 

δβang

3.3.5

(3.54)

g4+2ngf4+2ngOrganization of the Pressure Equations

The pressure equation, Eq. (3.53) include the pressures in cell (i, j, k) and six adjacent cells. To
solve this system wide set of equations, it is convenient to have a one dimensional index (onedimensional array), and it is deﬁned bym(i, j, k) = i + nx (j − 1) + nx ny (k − 1) ,− 41 −

41(3.55)

JAEA-Data/Code 2008-014

where nx and ny denote the number of cells in x(r) and y(θ) directions. The matrix for thecoeﬃcients of the system wide pressure equation set becomes a multi-band coarse matrix ofnx ny nz dimensions which has non-zero elements in the diagonal (i, j, k), both sides of the diagonal (i ± 1, j, k), and in the distances nx (i, j ± 1, k) and nx ny (i, j, k ± 1) from the diagonal, i.e.
seven non-zero bands. The mapping scheme by Eq. (3.55) is advantageous when the number ofcells in z direction is larger than that in x or y directions.
The pressure equation is solved by ILUBCG (Incomplete LU decomposition BiConjugate
Gradient) method[34].

3.3.6

Elements of the Jacobian Matrix

Full description of the mass and energy equations
A full description of the mass and energy equations that is solved, and derivatives for theJacobian matrix are given as follows. Note that the derivatives are given here in terms of theoriginal independent variables p, αm and Tm . A variable transformation to avoid numericalproblem with temperatures of missing components is described later.
Mass (steam):�Ms = V (αs ρs )n+1 + V

adj

n+1
� s,adj (pn+1) − ∆tΓsm
adj − pdV�s+(αs ρs )n∆t − V M
dt

�
∂Γs
∂ρs
− ∆t
−V
ms,adj
∂p
∂p
adj

∂Ms
∂p

= V αs

∂Ms
∂αs
∂Ms
∂Ts
∂Ms
∂Tl
∂Ms
∂αa
∂Ms
∂Ta

= V ρs − ∆t

∂Γs
∂αs

(3.56)(3.57)(3.58)

∂Γs
∂ρs
− ∆t
∂Ts
∂Ts
∂Γs
= −∆t
∂Tl
∂Γs
= −∆t
∂αa

(3.59)= 0(3.62)= V αs(3.60)(3.61)

Mass (steam+water):
Msl = V��αs ρs + 1 − αs −+(αs ρs + αl ρl )n
�

�
a

�

αa ρl

+V

�

= V

∂Msl
∂αs
∂Msl
∂Ts

= V (ρs − ρl )αs ρs + 1 − αs −

�

�

αa ρl

a

∂ρs
∂Ts

�
adj

dV
�s + M�l )∆t − V (M
dt

∂Msl
∂p

= V αs�n+1

n+1
� s,adj + m
� l,adj )(pn+1(m)
adj − p(3.63)

�

−V

�
adj

� s,adj + m
� l,adj )(m(3.64)(3.65)(3.66)− 42 −

42

JAEA-Data/Code 2008-014

∂Msl
∂Tl
∂Msl
∂αa
∂Msl
∂Ta

= V αl∂ρl∂Tl(3.67)= −V ρl(3.68)= 0(3.69)

Mass (non-condensible gases):
Ma = V (αa ρa )n+1 + V

�
adj

n+1
� a,adj (pn+1m)
adj − pdV�a+(αa ρa )n∆t − V M
dt(3.70)∂Ma∂p= V αa(3.71)

∂Ma
∂αs
∂Ma
∂Ts
∂Ma
∂Tl
∂Ma
∂αa
∂Ma
∂Ta

= 0

�
∂ρa
� a,adj
m
−V
∂p
adj

(3.72)= 0(3.73)= 0(3.74)= V ρa(3.75)= V αa∂ρa∂Ta(3.76)

Energy (steam):
Es = V (αs ρs es )n+1 + V

�
adj

(1)

(2)

n+1
(e�s,adj + e�s,adj pn+1 )(pn+1)
adj − p� s − ∆t(q n+1 + Γn+1 hn+1 ) − ∆t+V pn+1 αsn+1 + V pn+1 Dssis�s + (αs ρs es )n dV ∆t
−V E
dt

∂Es
∂p

∂Es
∂αs
∂Es
∂Ts
∂Es
∂Tl
∂Es
∂αa
∂Es
∂Ta

= V αs

=
=
=
=

�

∂ρs
∂es
es + ρs
∂p
∂p

�

�

−V

�
adj

(1)

�

n+1
qsa

a

(2)

{e�s,adj + e�s,adj (2p − padj )}�� s − ∆t ∂qis + ∂Γs hs + Γs ∂hs+V αs + V D
∂p
∂p
∂p
�
�
∂qis ∂Γs
V ρs es + V p − ∆t
+
hs
∂αs
∂αs
�
�
�
�
�
∂ρs
∂es
∂qis ∂Γs∂hs
V αses + ρs− ∆t+
h s + Γs
− ∆t
Rsa
∂Ts
∂Ts
∂Ts
∂Ts
∂Ts
a
∂Γs
−∆t
hs
∂T
� l
�
∂qis
∂Γs
−∆t
+
hs
∂αa ∂αa= ∆tRsa(3.77)(3.78)(3.79)(3.80)(3.81)(3.82)(3.83)− 43 −

43

JAEA-Data/Code 2008-014

Energy (steam+water):
Esl = V��αs ρs es + 1 − αs −

�

�

�n+1

αa ρl el
a
� (1)
(1)
(2)
(2)
n+1
+V
(e�s,adj + e�l,adj + e�s,adj pn+1 + e�l,adj pn+1 )(pn+1)
adj − p
adj
�
�n+1
�
�
�
�
n+1
n+1 �
n+1
�
+V p
1−
αa
+Vp
(Ds + Dl ) − ∆tqsa + qint
a
a
�
�
�
�n
�

�s + E�l ) + αs ρs es + 1 − αs −
−V (E

∂Esl
∂p

= V

�

αs

−V

�

�
adj

∂ρs
∂es
es + ρs
∂p
∂p
(1)

(1)

�

αa ρl el

a

+ αl
(2)

�

dV
∆t
dt

∂ρl
∂el
el + ρl
∂p
∂p
(2)

��

{e�s,adj + e�l,adj + (e�s,adj + e�l,adj )(2p − padj )}�s + D� l)+V (αs + αl ) + V (D∂Esl
= V (ρs es − ρl el )
∂αs
�
�
�
∂ρs
∂es
∂Esl
− ∆t
Rsa
= V αses + ρs
∂Ts
∂Ts
∂Ts
a
�
�
∂Esl
∂ρl
∂el
= V αlel + ρl
∂Tl
∂Tl
∂Tl
∂Esl
= −V ρl el − V p∂αa∂Esl= ∆tRsa∂Ta
Energy (non-condensible gases):
Ea = V (αa ρa ea )n+1 + V

�
adj

(1)

∂Ea
∂αs
∂Ea
∂Ts
∂Ea
∂Tl
∂Ea
∂αa
∂Ea
∂Ta

�a + (αa ρa ea )n dV ∆t−V Edt= V αa�∂ρa∂eaea + ρa∂p∂p(3.85)(3.86)(3.87)(3.88)(3.89)(3.90)(2)

n+1(e�a,adj + e�a,adj pn+1 )(pn+1)
adj − p

� a + ∆t
+V pn+1 αan+1 + V pn+1 D∂Ea∂p(3.84)

�−V�a+V αa + V D

�
adj

�

n+1
qsa

a

(3.91)(1)(2){e�a,adj + e�a,adj (2p − padj )}(3.92)= 0(3.93)= ∆tRsa(3.94)= 0(3.95)= V ρa ea − V p(3.96)= V αa�∂ρa∂eaea + ρa∂Ta∂Ta�− ∆tRsa− 44 −

44(3.97)

JAEA-Data/Code 2008-014

Derivatives of the phase change and heat transfer terms
The phase change (evaporation, condensation) term given by Eq. (2.103) is dependent on steam
and water temperatures and also on the pressure and volume fractions through the saturation
temperature at steam partial pressure. Derivatives for those variables necessary for the Jacobianmatrix are given as follows.
∂Γs
∂p
∂Γs
∂αs
∂Γs
∂Ts
∂Γs
∂Tl
∂Γs
∂αa
∂qis
∂p
∂qis
∂αs
∂qis
∂Ts
∂qis
∂αa

=
=
=
=
=
=
=

his + hil ∂Tsat ∂ps
Γs
−Ai
− Ai
hs − hl ∂ps ∂phs − hl
his + hil ∂Tsat ∂ps−Ai
hs − hl ∂ps ∂αshisΓs ∂hsAi− Ai
hs − hl
hs − hl ∂TshilΓs ∂hlAi+ Ai
hs − hl
hs − hl ∂Tl
his + hil ∂Tsat ∂ps−Ai
hs − hl ∂ps ∂αa∂Tsat ∂ps
Ai his∂ps ∂p∂Tsat ∂ps
Ai his∂ps ∂αs�∂hs ∂hl−∂p∂p�= −Ai his
= Ai his(3.98)(3.99)(3.100)(3.101)(3.102)(3.103)(3.104)(3.105)∂Tsat ∂ps∂ps ∂αa(3.106)The derivatives of pseudo partial pressure of steam (Eq. (2.101)) are given as follows.
∂ps
∂p
∂ps
∂αs
∂ps
∂αa

αs
α
α − αs= pα2αs= −p 2α=(3.107)(3.108)(3.109)Variable transformation
To avoid the singularity for Tm of missing components, a variable transformation T → αT isdone before constructing the Jacobian matrix. The transformation is expressed byf (p, αs , Ts , αa , Ta (a = a1 , ..., ang )) → f (P, As , βs , Aa , βa (a = a1 , ..., ang ))(3.110)where the new set of independent variables are deﬁned as follows.
= p(3.111)

As = αs(3.112)Pβs = αs Ts
βl = αl Tl = (1 − αs −Aa = αa ,βa = αa Ta�(3.113)αai )Tl(3.114)(3.115)(3.116)− 45 −

45

JAEA-Data/Code 2008-014

The relation between partial derivatives are given as follows.
∂f
∂P
∂f
∂As
∂f
∂βs
∂f
∂βl
∂f
∂Aa
∂f
∂βa

=
=
=
=
=
=

∂f
∂p
∂f
∂f Ts
∂f Tl
−
+
∂αs ∂Ts αs ∂Tl αl∂f 1∂Ts αs∂f 1∂Tl αl
∂f
∂f Tl
∂f Ta
+
−
∂αa ∂Tl αl ∂Ta αa∂f 1∂Ta αa− 46 −

46(3.117)(3.118)(3.119)(3.120)(3.121)(3.122)

JAEA-Data/Code 2008-014

4 Veriﬁcation and Example Calculations
4.1

Check of the Numerical Behavior of the Models

4.1.1

Melt Jet and Pool

PurposeFree fall and spreading of a column of liquid was simulated for conﬁrmation of physically proper
behavior of the melt jet and pool models and their mass conservation accuracy.
Calculation conditions
Figure 4.1 shows the geometric condition. A liquid column, diameter 90mm and initial velocity
0.5m/s at inlet, ﬂows into a closed chamber of radius 0.2m and height 1m ﬁlled with gas. The
liquid column is subjected to the gravitational acceleration. The chamber is discretized into
0.1m cells in vertical direction, and each vertical cell is sub-divided into ﬁner cells for the melt
jet. The base case condition and some variations are deﬁned as follows.
• Base case: jet cell size 2cm, time step 1ms• “dz” case: jet cell size 5mm (×0.25)• “dt” case: time step 0.1ms (×0.1)
Also, the result is compared with theoretical solution (“Theoretical” in ﬁgures) and the up-windﬁnite diﬀerence solution (“FDM up-wind” in ﬁgures).
Calculation results
The leading edge proﬁle and trailing edge proﬁle are shown in Figs. 4.2 and 4.3, respectively.
The base case solution shows a reasonable agreement with the theoretical solution, and thatis much better than the up-wind ﬁnite diﬀerence because of the usage of CIP method. The
calculation with a ﬁner mesh showed a result closer to the theoretical solution. The case with asmaller time step showed a result almost identical to the base case.
Figure 4.4 shows the comparison of the results with the free fall curve. Admitting that theplotted leading edge position by calculation is aﬀected by the criterion with which the existenceof the jet is detected in cells, the leading edge position agrees with the progress by free fall.
Figure 4.5 shows the evolution of jet and pool proﬁles in the base case. The mass conservation
accuracy is checked during the mass transfer from the jet model to the pool model. Figure 4.6
shows the model has a good mass conservation. The total liquid mass conservation error throughthe simulation was < 0.2%.

− 47 −

47

JAEA-Data/Code 2008-014

C.L.

Liquid column
inlet D0, V0

1.0m

Closed space
filled with gas

0.2m

Figure 4.1: Melt jet and pool test calculation geometry.

1

time=0.3 s

Height (m)

0.8

Base
dz x0.25
dt x0.1
FDM upwind
Theoretical

0.6

0.4

0.2

0

0

0.02

0.06
0.04
Radius (m)Figure 4.2: The jet leading edge proﬁle at 0.3s.

48 −
− 48

JAEA-Data/Code 2008-014

1

time=1.0 s

Height (m)

0.8

Base
dz x0.25
dt x0.1
FDM upwind
Theoretical

0.6

0.4

0.2

0

0

0.02

0.04

0.06

Radius (m)Figure 4.3: The jet trailing edge proﬁle at 1.0s.

1Basedz x0.25
dt x0.1
Free fall

Height (m)

0.8

0.6

0.4

0.2

0

0

0.1

0.2

0.3

0.4

0.5

Time (s)Figure 4.4: Comparison of the jet leading edge progress with the free fall curve.

49 −
− 49

JAEA-Data/Code 2008-0141

t=0.1 st=0.2 st=0.3 st=0.4 s

0.1

0.1

0.1

t=0.5 sHeight (m)

0.8
0.6
0.4
0.2
0

0

0.1

0.2 0

0.2 0

0.2 0

0.2 0

0.1

0.2

1

t=1.0 st=0.8 st=0.6 st=1.4 st=1.2 sHeight (m)

0.8

Jet

0.6
0.4

Pool

0.2
0

0.1

0

Radius (m)

0.2 0

0.1

Radius (m)

0.2 0

0.1

Radius (m)

0.2 0

0.1

0.2 0

Radius (m)

0.1

0.2

Radius (m)Figure 4.5: Overall proﬁles of the melt jet/pool evolution.

20

18

16

Mass (kg)

14

12

10

8

6

4
Jet mass
Pool mass
Total mass

2

0

0

0.5

1

1.5

2

Time (s)Figure 4.6: Plot of the jet, pool and total melt mass for the check of mass conservation.

50 −
− 50

JAEA-Data/Code 2008-014

4.1.2

Melt Particles

PurposeFree fall of the melt particles and their merging into the melt pool, or stacking on the ﬂoor was
simulated for conﬁrmation of physically proper behavior of the melt particles and pool models
and their mass conservation accuracy.
Calculation conditions
Figure 4.7 shows the geometric condition. Melt particle groups are initially distributed in theupper area of the calculation ﬁeld that is a cylinder of radius 0.2m and height 1m ﬁlled withsteam, with constant radial and vertical intervals. The diameter of the individual particle in
the groups is 1mm and the number of the particles in each group is given so that their particle
volume fraction becomes 0.2. The initial vertical velocity of the particle groups is set 0 andthey are subjected to the free fall. The base case condition and some variations are deﬁned asfollows.
• Base case: particles are molten, initial x-direction velocity 0, time step 1ms• “dt” case: time step 0.1ms (×0.1)• “Slid” case: particles are solid• “Solid-rand” case: particles are solid and have initial x-direction velocity randomly givenin the range −0.2 < vx < 0.2m/s
Calculation results
The evolution of the particle group distribution is shown in Figs 4.8–4.10, by the positions of
the centers of particle groups. In the base case, the particle groups fall on to the ﬂoor and makea melt pool. In “Solid” and “Solid-rand” cases, the particle groups stack on the ﬂoor and makea debris bed. In the “Solid-rand” case, the randomly given horizontal initial velocities make adiﬀerence from the “Solid” case.
Figure 4.11 shows a comparison between the leading edge progress obtained in the calculationand the free fall curve. They agreed well.
Figure 4.12 shows that a good mass conservation accuracy is kept during the mass transferfrom the particle group model to the pool model. The total melt mass conservation error throughthe simulation was < 0.01%.

− 51 −

51

JAEA-Data/Code 2008-014

C.L.

1.0m

ParticlegroupsClosed space
filled with gas

0.2m

Figure 4.7: Melt particles and pool calculation geometry.

1

t=0.0 st=0.2 st=0.6 st=0.4 sHeight (m)

0.8
0.6
0.4

Particle groups

0.2

Pool
00

0.1
0.2 0
Radius (m)

0.1
0.2 0
Radius (m)

0.1
0.2 0
Radius (m)0.1
0.2Radius (m)Figure 4.8: Evolution of particle group distribution: Base case.

52 −
− 52

JAEA-Data/Code 2008-014

1

t=0.0 st=0.2 st=0.6 st=0.4 sHeight (m)

0.8
0.6

Particle groups
stacked on the
floor

0.4
0.2
00

0.1
0.2 0
Radius (m)

0.1
0.2 0
Radius (m)

0.1
0.2 0
Radius (m)0.1
0.2Radius (m)Figure 4.9: Evolution of particle group distribution: “Solid” case.

1

t=0.0 st=0.2 st=0.4 st=0.6 sHeight (m)

0.8
0.6

Particle groups
stacked on the
floor

0.4
0.2
00

0.1
0.2 0
Radius (m)

0.1
0.2 0
Radius (m)

0.1
0.2 0
Radius (m)0.1
0.2Radius (m)Figure 4.10: Evolution of particle group distribution: “Solid-rand” case.

53 −
− 53

JAEA-Data/Code 2008-014

0.8

Base
dt x0.1
Solid
Theoretical

0.7

0.6

Height (m)

0.5

0.4

0.3

0.2

0.1

0

0

0.1

0.2

0.3

0.4

0.5

Time (s)Figure 4.11: Comparison of the particle leading edge progress with the free fall curve.

90

80

70

Mass (kg)

60

50

40

30

20
Droplet mass
Pool mass
Total mass

10

0

0

0.1

0.2

0.3

0.4

0.5

0.6

Time (s)Figure 4.12: Plot of the droplet, pool and total melt mass for the check of mass conservation.

54 −
− 54

JAEA-Data/Code 2008-014

4.1.3

Shock Wave Propagation in Two-Phase Medium

PurposeBecause the shock wave propagation in a two-phase medium is a primary physical process to
be captured by the two-phase ﬂow model in a steam explosion simulation, the ability of thetwo-phase ﬂow model in this aspect was checked by simple simulations.
First, dependence of the sonic velocity on the void fraction was checked by a nearly onedimensional geometry. Then, spherical wave propagation was tested in a two-dimensional (cylindrical) geometry.
Two-phase sonic velocity
Calculation conditions Figure 4.13 shows the geometric condition. A cylindrical chamberof radius 0.2m and height 2.5m, ﬁlled with water up to the level 2.05m was considered. A
pressure source was initially placed in the center bottom cell, where void fraction 0.5 and pressure
0.12MPa is given. The initial system pressure was set 0.1MPa.
The sonic velocity was evaluated by tracing the propagation of the pressure front, andcompared with theoretical values evaluated for homogeneous bubbly medium given byctphm (α) =

�

�

ρl
ρg
α + α(1 − α) + (1 − α)2 + α(1 − α)
ρg
ρl
2

��

cg
cl

�2 �−1/2.

(4.1)

The void fraction in the water column was changed as a parameter.
Calculation results Figures 4.14–4.16 show the pressure histories at height 0.475m, 0.975m
and 1.95m with void fractions 10−5 , 10−3 and 0.1 as a parameter. The time of the arrival of theshock front was read at the half value of the ﬁrst peak. It was observed that the shock front
propagation delayed with increase of the void fraction, and that the shock wave was rapidly
damped when the void fraction was 0.1.
The sonic velocity was calculated from the shock arrival time, and plotted with the theoreticalvalue (Eq. (4.1)) in Fig. 4.17. The sonic velocity evaluated by the calculation results agreed wellwith the theoretical values.
Spherical shock wave propagation
Calculation conditions Figure 4.18 shows the geometric condition. A cylindrical chamberof radius 1m and height 2m, ﬁlled with water up to the level 1.78m was considered. A pressuresource was initially placed in the center bottom cell, where void fraction 0.5 and pressure 1MPa
is given. The initial system pressure was set 0.1MPa. The void fraction in the water pool was
set 10−4 .
Calculation results Figure 4.19 shows the pressure proﬁles obtained by the calculation. The
plots show that the shock front propagates spherically from the pressure source, and reﬂects
at the side wall and the ﬂoor. Thus, it is demonstrated that the wave dynamics is reasonably
captured by the present two-phase ﬂow model.

− 55 −

55

JAEA-Data/Code 2008-014

C.L.
2.5

Air
2

Water

1.5

Cells wherepressure riseis detected
1

0.5

0

Pressuresource(Air 0.12MPa
Void frac. 0.5)
0

0.2

Figure 4.13: Geometry for the shock wave propagation test calculation.

56 −
− 56

JAEA-Data/Code 2008-014

0.135

0.13

=1e-5
=1e-3
=0.1

At z=0.475m

Pressure (MPa)

0.125

0.12

0.115

0.11

0.105

0.1

0.095
1e-04

0.001

0.01
Time (s)Figure 4.14: Pressure histories at z = 0.475m with various void fractions.

0.13

0.125

=1e-5
=1e-3
=0.1

At z=0.975m

Pressure (MPa)

0.12

0.115

0.11

0.105

0.1

0.095
1e-04

0.001

0.01
Time (s)Figure 4.15: Pressure histories at z = 0.975m with various void fractions.

57 −
− 57

JAEA-Data/Code 2008-014

0.12

=1e-5
=1e-3
=0.1

At z=1.95m

Pressure (MPa)

0.115

0.11

0.105

0.1

0.095
1e-04

0.001

0.01
Time (s)Figure 4.16: Pressure histories at z = 1.95m with various void fractions.

Sonic velocity ratio (=chom/cg)

10

TheoryCalc (z=0.475m)
Calc (z=0.975m)
Calc (z=1.95m)

1

0.1

0.01

1e-05

1e-04

0.001

0.01

0.1

1

Void fraction

Figure 4.17: Dependence of two-phase sonic velocity on void fraction—comparison of homogeneous two-phase medium theory and JASMINE calculation.

58 −
− 58

JAEA-Data/Code 2008-014

C.L.

Void
frac.
1

0.8

2m

0.6

0.4

Pressuresource1MPa Air
Void frac.=0.5

0.2

0

z

1m

r

Figure 4.18: Geometry for the two-dimensional shock wave propagation test calculation.

1.2e+06

t= 0.0000

t= 0.0005

t= 0.0010

t= 0.0020

t= 0.0025

1e+06
800000600000400000r

2000000z
1.2e+06

t= 0.0015

1e+06

Pressure
scale (Pa)
1.2e+06
1e+06

8000008000006000006000004000004000002000002000000

0Figure 4.19: Pressure proﬁles obtained by the calculation, showing spherical propagation of theshock front and its interference with the reﬂected wave from the side wall.

59 −
− 59

JAEA-Data/Code 2008-014

Table 4.1: Conditions and results of FARO melt jet quenching experiments, and analyticalconditions.
L14
Melt
Material

L28

L31

Corium (Tsol /Tliq = 2830/2850(K),80wt%UO2 -20wt%ZrO2 )
125
175
92
3073
3053
2990
92
44
48

Released mass (kg)Temperature (K)Jet diameter at inlet (mm)Water
Temperature (Subcool) (K)537 (1)424 (1) 291 (104)Pool depth (m)
2.05
1.44
1.45
Water vessel diameter (m)
0.71
0.71
0.71
Cover gas
Atmosphere
Steam
Steam
Argon
Initial pressure (MPa)
5.0
0.51
0.22
1.26
3.53
3.49
Cover gas volume(m3 )Results
Pressure rise at plateau (MPa)
2.5
1.2
0.04
Level swell (m)
∼1.4
∼0.6
∼0.3
Agglomerated mass (kg)
20 (16%) 77 (44%)
0
5.0
3.0
3.3
Debris mean dia.a) (mm)Analytical
Grid (Water pool zone)
6×37(6×20)
8×25(6×14)
Central column radius (m)
0.15
0.10
0.10
Time step (ms)∼0.5
Jet initial velocity(m/s)
3.0
∼3.0
∼2.7
Jet ﬂow-in duration(s)
0.80
5.2
2.5
Melt particle dia.(const.)(mm)
5.0
3.0
3.0
Taylor type, Cent = 1
Jet break-up length modelb)
Surface temperature drop model
Use
Npcr c)1000
nhist c)1000
Fhtint d)
0.02a) mass median diameter b) see 2.2.1 c) see 2.2.3 d) see 2.3.2

4.2

Simulation of Premixing Experiments

PurposeFor the veriﬁcation of the premixing-related model functions and the parameter tuning, we
performed simulation of FARO experiments that were performed at JRC Ispra of EU[35, 36].
In those experiments, ∼100kg of molten corium (mixture of UO2 and ZrO2 ) was dropped into awater pool. Data for the pressurization and water level swelling by steam generation during themelt jet break-up and quenching, particle size distribution of the quenched debris are available.
We referred three experiments, L14, L28 and L31 with diﬀerent subcools, pressures and jet sizes,
and compared the calculation and experimental results.
Calculation conditions
The experimental conditions of the three experiments, and analytical conditions for the simulations are summarized in Table 4.1. Figure 4.20 shows the grids to simulate the FARO
experimental geometry.

60 −
− 60

JAEA-Data/Code 2008-014

C.L.

5

C.L.
3

4
z (m)

Melt jet
inlet

2.5

Melt jet
inlet

Water
vessel

2

3
z (m)

1.5

2
1

1
0.5

0

0

0 0.35
r (m)

0

0.3
r (m)

0.6

Figure 4.20: Grid geometry to simulate the experimental system of FARO-L14 (left), L28 andL31 (right).

61 −
− 61

JAEA-Data/Code 2008-014

L14 exp.
L14 cal.
L28 exp.
L28 cal.
L31 exp.
L31 cal.

Pressure (Pa)

107

10

6

0

2

6

4

8

10

Time (s)Figure 4.21: Comparison of calculation and experimental results: Pressure history.
Calculation results
Calculation results are compared with experimental data in Figs. 4.21–4.23. The pressure histories (Fig. 4.21) and the level swellings (Fig. 4.22) were roughly in reasonable agreement. In
detail, we can see a trend of underestimating the pressure rise while overestimating the levelswelling, especially for low pressure conditions (L28 and L31). This trend depicts the present
limitation in the premixing modeling.
We consider that steam bubbles generated in the premixing zone escapes faster in reality
than in the simulation especially in low pressure conditions. It seems that we need considerable
work in two-phase ﬂow constitutive models if we want to improve this situation.
Figure 4.23 shows the comparison on the mass fraction of agglomerated melt, that was foundas a continuous lump at the bottom. The lump can be made either by direct arrival of the jetcolumn without fragmentation during the fall or by re-agglomeration of once-fragmented melt
droplets. The simulation results did not agree with experimental results.
These results show that ability of the present model in simulating the complicated premixing
phenomena and long term cooling behavior is not in a satisfactory level. Nevertheless, we canuse the present model with a certain level of conﬁdence for the purpose of steam explosion
simulations, in which mainly the “premixed” melt mass in the initial phase of premixing isneeded. (see the next section.)

62 −
− 62

JAEA-Data/Code 2008-014

3.5

L14 exp.
L14 cal.
L28 exp.
L28 cal.
L31 exp.
L31 cal.

Cross sectional average void=0.8 for cal.

Water level (m)

3

2.5

2

1.5

1

0

2

4

6

8

10

Time (s)Figure 4.22: Comparison of calculation and experimental results: Water level swelling.

Agglomerated mass fraction (cake or pool)

0.8

L14 exp.
L14 cal.
L28 exp.
L28 cal.
L31 exp.
L31 cal.

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0

0

2

4

6

8

10

Time (s)Figure 4.23: Comparison of calculation and experimental results: Agglomerated melt massfraction.

63 −
− 63

JAEA-Data/Code 2008-014

4.3

Simulation of Explosion Experiments

PurposeFor the veriﬁcation of the explosion-related model functions and the parameter tuning, we
performed simulation of KROTOS experiments and a FARO experiment that were performed
at JRC Ispra of EU[37, 38]. In KROTOS experiments, several kilograms of molten alumina or
corium (mixture of UO2 and ZrO2 ) was dropped into a column of water, and a triggering for asteam explosion was given by a pressure pulse generated by releasing compressed argon gas. We
referred three cases, 44, 42 and 37 from the KROTOS series experiments. We referred anotherexperiment, FARO-L33, a steam explosion experiment in a larger scale, in which about 100kg
of corium melt was dropped into a water pool, and a steam explosion was triggered by pressure
pulses from a detonator. For those experiments, data for transient pressure histories at several
locations showing shock wave propagation and escalation behavior, and particle size distributionof the ﬁne fragments produced by the steam explosion are available.
Calculation conditions
In the calculation, we simulated the experimental condition as summarized in the Table 4.2.
Figure 4.24 shows the grids to simulate the KROTOS and FARO experimental systems.
The calculation of the steam explosion needs two steps, premixing and explosion calculations.
In the premixing calculation, the model was tuned so that experimentally observed premixing
behavior was reasonably simulated. Namely, the jet break-up length correlation and its tuning
factor were selected so that jet break-up length observed in a similar condition but without
explosion was reproduced in the simulation. The melt particle diameter also was given in asimilar way. For FARO-L33 experiment, the Taylor type jet break-up length correlation (seeSection 2.2.1) should be valid, and it was applied. While, in KROTOS experiments, it is not
the case due to a transient nature of the melt delivery.[37]After conﬁrming that reasonable premixing calculation results, i.e. agreement with experimentally observed overall void fraction at the time of triggering, were obtained, we used thoseresults as the initial conditions for the explosion calculations.
In the explosion calculations, on the other hand, we set a consistent set of model parameters
for all the cases independent of the melt material and other conditions, i.e. Cf rg = 0.35,
df rg = 50(µm), ttriglif e = 1 (ms) and Kev = 0.7. It was demonstrated that as long as thepremixing condition was given reasonably, the explosion results were simulated reasonably withthe same set of the explosion model parameters. This set of explosion model parameters was
obtained by tuning them to have agreement with KROTOS alumina experiments on the pressure
pulses and the mass fractions of ﬁne debris produced by the steam explosion.
A rationale of this tuning method is as follows. Figure 4.25 shows the diﬀerence between thealumina and corium cases of KROTOS experiments on the total, molten (not frozen), molten
and premixed (in the zone where void fraction is below 0.75) masses of corium obtained by thepremixing calculations. Alumina melt is kept molten through out the period of the premixing
stage and does not cause too much steam generation in the KROTOS condition. Thus, an ideal
premixing condition, i.e. most of the melt is molten and well-mixed with water, is realized.
On the other hand, with corium melt, considerable amount of the melt is already frozen or
enveloped in a highly voided zone at the time of triggering, and, only a small part of the melt
can actually participate in the explosion process. This diﬀerence is caused by material physicalproperties and the initial temperature. Alumina has lower melting and initial temperature,
lower density and a larger latent heat than corium. Then, it makes larger droplets and takes
much longer to freeze than corium. Also, heat transfer to water is much slower than corium
due to the lower temperature as well as smaller surface area per volume. This explains the well− 64 −64

JAEA-Data/Code 2008-014

C.L.
3.0

C.L.

2.5

1.5
Jet inlet

Jet
Inlet

Water vessel(inner wall)

Initialwater level
K5

1.0

2.0

Initial water
level

z (m)

1.5

Water vessel(inner wall)

z (m)

K4
K3

0.5

K2

Triggersource

K1
K0

0.0

1.0

P-measurement
positions in theexperimentP-measurement
positions in theexperiment
0.5
Cells where
the pressure is Trigger
comparedsourceCells wherepressure iscompared

0.0

0.1
r (m)

0.355
r (m)

Figure 4.24: Grid geometry to simulate steam explosion experiments KROTOS series (left) and
FARO-L33 (right).
known experimental observation that corium melt usually does not cause as strong explosions
as alumina melt does. [39]Considering the above, the alumina experiments are more suitable as references for thetuning of fragmentation model parameters due to the simple and ideal premixing condition.
Application of the same parameter setting for corium cases should reasonably estimate weaker
explosions if the diﬀerence in the premixing condition, i.e. more fraction of frozen melt andhigher void fraction, is adequately simulated.
External triggering was simulated by placing a certain amount of high pressure gas in thecenter bottom cell. For the KROTOS experiments, the gas chamber used in the experiment
was simulated by its speciﬁcation (15cm3 , 14MPa, pV ∼ 220kJ). For FARO-L33, an equivalent
volume and pressure of the gas (28.6cm3 , 35MPa) was deduced from the detonator energy ∼1kJ.
Figure 4.26 shows the calculation result for propagation of the pressure pulse by the triggering
device in pure water in the FARO-L33 geometry. It was conﬁrmed that the calculated trigger
pulse is in a reasonable agreement with experimental data[40].
Calculation results
Figures 4.27 and 4.28 show the comparison of calculated and experimentally observed pressure
histories at various vertical locations for KROTOS alumina experiments. With the tuned explosion model parameters, the calculation results reasonably agreed with the experimental data.
The inﬂuence of the water subcool was not signiﬁcant. Figure 4.29 shows the calculation results for KROTOS-37, the corium experiment. It was reported that no propagating energetic− 65 −

65

JAEA-Data/Code 2008-014

Table 4.2: Conditions and results of KROTOS-44, 42, 37 and FARO-L33 experiments, andanalytical conditions (“K#” denote run numbers of KROTOS series).
K44
K42
K37
L33
Material
Aluminaa)
Coriumb)
Released mass (kg)
1.50
1.54
3.22
100
Temperature(K)
2673
2465
3018
3070
Jet release dia. (mm)3048Free fall height(m)
0.44
0.77
Water
Temperature(Subcool) (K)363(10)293(80)296(77)294(124)
Depth (m)
1.105
1.62
Vessel diameter (m)
0.20
0.71
Cover gas Atmosphere
Steam
Argon
Pressure (MPa)
0.10
0.41
Cover gas volume(m3 )
0.290
3.496
Premixing Level swell (cm)
12
3
30
9
results
Jet break-up length (m)
∼0.3
∼0.8
∼1.1
Explosion Pressure peaks (MPa)
∼50
∼50
NAc)
∼6
results
Half height width (ms)
∼1.5
∼1
NAc)
∼3
Debris <0.106mm (%)
47
31
1.4
∼8d)
Analytical Time step (µs)∼2conditions Jet initial velocity(m/s)
3.0∼2.9
Jet ﬂow-in duration(s)
0.270
0.255
0.191
2.6
Melt Particle dia.(mm)102
3
Jet break-up length modele)
S: Cent = 2.5
S: Cent = 1.7 T: Cent = 1
Surface temperatureUsedrop model
Npcr f )100
1000
nhist f )500
1000
External trigger
14.8MPa-15cm3
35MPa-29cm3
Trigger time (s)
0.9
0.5
1.12
g)
Cf rg
0.35
df rg (µm) g)
50
ttriglif e (ms) g)
1.0
ptrig (MPa) g)
0.2
0.5
Fragmentation criterion
Tav ≥ Tmp
Kev h)
0.7a) Tmelt = 2300(K) b) UO2 (80)-ZrO2 (20) (wt%), Tsol /Tliq = 2830/2850(K)c) Reportedly “no propagating energetic explosions”[37].
d) Extrapolated from the size distribution reported in [41]. The ratio to themelt mass in the system at the triggering time, ∼40kg, is ∼20%.
e) S: Saito et al. correlation, T: Taylor type correlation (see 2.2.1)f ) see 2.2.3 g) see 2.2.4 h) see 2.3.2
Melt

66 −
− 66

JAEA-Data/Code 2008-014

K42 Molten, <0.75

K42 Total mass
K42 Molten

Mass (kg)

1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

Triggering

1

0.8

0.6

0.4

0.2

0

Time (s)K37 Molten, <0.75

K37 Total mass
K37 Molten

3.5

3Mass (kg)2.5

Triggering

2

1.5

1

0.5

0

1

0.8

0.6

0.4

0.2

0

Time (s)Figure 4.25: Diﬀerence of the premixing melt mass due to freezing and void in KROTOS alumina(K42, above) and corium (K37,below) experiments.

8

490mm
715mm
940mm
1165mm
1390mm7

Pressure (MPa)

6
5
4
3
2
1
0

0

0.5

1

1.5
Time (ms)

2

2.5

3

Figure 4.26: Simulation of FARO-L33 external trigger (propagation of the pressure pulse in purewater).

67 −− 67

JAEA-Data/Code 2008-014

explosion was observed in the experiment.[37] In the calculation, the triggering pressure pulse
escalated only near the bottom, K1–K2, and damped out later.
Figure 4.30 shows the comparison of calculated and experimentally observed pressure histories at various vertical locations for FARO-L33. The scale of the pressure pulse is reasonably
simulated by the calculation.
Figures 4.31 and 4.32 show the evolution of the total ﬂuid kinetic energy in the calculationin comparison with experimentally estimated kinetic energies. A good agreement was obtainedfor KROTOS alumina cases and FARO-L33. The calculation for KROTOS-37 showed very low
kinetic energy, in agreement with the fact of no energetic explosion in the experiment.
Generally, very ﬁne debris of the size smaller than 0.1 mm is produced by steam explosions.
In the calculation, the ﬁne debris produced by the explosion is modeled by ﬁne fragments of
uniform size, 50µm. Figures 4.33 and 4.34 show comparison of the mass fraction of ﬁne debris inexperiments and in the calculations. Though rigorous comparison is not allowed due to diﬀerencein the classiﬁcation of the debris size, the calculated debris mass fractions are comparable tothose observed in the experiments.

− 68 −

68

JAEA-Data/Code 2008-014

Pressure (MPa)80

Pressure (MPa)
Pressure (MPa)

K4

EXP
CAL

K3

EXP
CAL

K2

EXP
CAL

K1

EXP
CAL

40
20

60
40
20
0
80

60
40
20
0
80

Pressure (MPa)

EXP
CAL

60

0
80

60
40
20
0
80

Pressure (MPa)

K5

60
40
20
0

0

1

2

3

4
Time (ms)

5

6

7

8

Figure 4.27: Comparison of calculation and experimental results: KROTOS 44 (alumina, saturation temperature) pressure histories at diﬀerent vertical locations (K1–K5 refers locations ofmeasurement, see Fig. 4.24).

69 −
− 69

JAEA-Data/Code 2008-014

Pressure (MPa)80

Pressure (MPa)
Pressure (MPa)

K4

EXP
CAL

K3

EXP
CAL

K2

EXP
CAL

K1

EXP
CAL

40
20

60
40
20
0
80

60
40
20
0
80

Pressure (MPa)

EXP
CAL

60

0
80

60
40
20
0
80

Pressure (MPa)

K5

60
40
20
0

0

1

2

3

4
Time (ms)

5

6

7

8

Figure 4.28: Comparison of calculation and experimental results: KROTOS 42 (alumina, subcooled condition) pressure histories at diﬀerent vertical locations (K1–K5 refers locations ofmeasurement, see Fig. 4.24).

70 −
− 70

JAEA-Data/Code 2008-014

30

CAL K1CAL K2CAL K3
CAL K4
CAL K525

Pressure (MPa)

20

15

10

5

0
0

1

2

3

4

5

6

7

8

Time (ms)Figure 4.29: Calculated pressure histories for KROTOS 37 (corium, subcooled condition) at
various vertical locations.

71 −
− 71

Pressure (MPa)JAEA-Data/Code 2008-014

10z=1615mm

EXP
CAL

z=1390mm

EXP
CAL

z=1165mm

EXP
CAL

z= 940mm

EXP
CAL

z= 715mm

EXP
CAL

z= 490mm

EXP
CAL

8
6
4
2

Pressure (MPa)

0

10
8
6
4
2

Pressure (MPa)

0

10
8
6
4
2

Pressure (MPa)

0

10
8
6
4
2

Pressure (MPa)

0

10
8
6
4
2

Pressure (MPa)

0

10
8
6
4
2
0

0

2

4

6

8

12
10
Time (ms)

14

16

18

20

Figure 4.30: Comparison of calculation and experimental results: FARO-L33 (corium, subcooled,large scale) pressure histories at diﬀerent vertical locations.
72 −
− 72

JAEA-Data/Code 2008-014

0.2

KROTOS-44
KROTOS-44 exp.
KROTOS-42
KROTOS-42 exp.
KROTOS-37

Fluid kinetic energy (MJ)

0.15

0.1

0.05

0

0

1

2

3

4
Time (ms)

5

6

7

8

Figure 4.31: Comparison of calculation and experimental results: KROTOS total ﬂuid kinetic
energy (Experimental data is available for the ﬁnal value; for KROTOS 37, only “no explosion”
was reported).
0.3

FARO-L33
FARO-L33 exp.

Fluid kinetic energy (MJ)

0.25

0.2

0.15

0.1

0.05

0

0

5

10

15

20

Time (ms)Figure 4.32: Comparison of calculation and experimental results: FARO-L33 total ﬂuid kinetic
energy (Experimental data is available for the ﬁnal value).

73 −
− 73

JAEA-Data/Code 2008-014

0.6

KROTOS-44
KROTOS-44 exp.<0.106mm
KROTOS-42
KROTOS-42 exp.<0.106mm

Fragmented mass fraction

0.5

0.4

0.3

0.2

0.1

0

0

1

2

3

4
Time (ms)

5

6

7

8

Figure 4.33: Comparison of calculation and experimental results: KROTOS fragmented mass
(Experimental data is available for the ﬁnal value).

0.2

FARO-L33
FARO-L33 exp.<0.25mm

Fragmented mass fraction

0.15

0.1

0.05

0

0

5

10

15

20

Time (ms)Figure 4.34: Comparison of calculation and experimental results: FARO-L33 fragmented mass
(Experimental data is available for the ﬁnal value).

74 −
− 74

JAEA-Data/Code 2008-014

4.4

Application for Plant Scale Simulation

PurposeAs an example for the application of JASMINE in plant scale, a simulation of PWR ex-vessel
steam explosion is demonstrated.
Calculation conditions
Figure 4.35 shows the geometric conditions and the grid for the present analysis assuming atypical PWR cavity. The calculation condition is as follows.
• Water pool size: 4.1m deep, 2.75m radius• Pressure vessel lower head: 2.2m radius, bottom position z=5.0m (melt free fall height
0.9m)• Pressure/water temperature: 0.2MPa, 342K (subcool 50K)

• Boundary condition: ﬂow-out boundary at the top (p =0.2MPa, constant), Slip walls atbottom and side
• Melt material: UO2 (80wt%)-ZrO2 (20wt%) corium, Tsol /Tliq = 2830/2850K, Tini = 2950K(superheat 100K)• Melt jet diameter: 0.5m

• Melt release velocity: ∼5.5m/s• Triggering pressure source: bottom center, p =10MPa, pV =1MJ

• Triggering time: at the ﬁrst peak of premixing mass, 0.7s

• Premixing model parameters: jet break-up length by Taylor type correlation with Cent = 1,
droplet diameter 5mm• Explosion model parameters: Cf rg = 0.35, df rg = 50(µm), ttriglif e = 1 (ms) and Kev = 0.7(as tuned in FARO/KROTOS simulation in the previous section)
Premixing calculation result and triggering time
The triggering time was chosen according to the premixing behavior. Figure 4.36 shows voidand melt distribution proﬁles at diﬀerent times in the premixing calculation. An oscillatory
nature is observed, i.e. penetration and break-up of the melt jet, large vapor pocket generation
causing worse contact of the water and melt, escape of the vapor and re-contact of the waterwith the melt.
The melt mass that eﬀectively participates in the explosion process is the mass that is kept
molten and resides where void fraction is not too high. So, the mass of the melt jet and particles(melt pool is excluded) whose temperature is below the melting point, and that exist in cells
where void fraction is less than 0.75 is deﬁned as “premixed mass”.
Figure 4.37 shows the history of the total mass, the molten mass and the premixed mass. The
molten mass of the jet and particles are kept at almost constant level, about 5300kg, due to thebalance of supply, freezing and merge to the pool. The premixed mass is oscillating reﬂecting theabove mentioned oscillation of the void generation and escape. From our calculation experiences,the most energetic explosion is obtained by applying a triggering at the time when the premixed
mass takes the ﬁrst peak, about 0.7s from Fig. 4.37.
− 75 −

75

JAEA-Data/Code 2008-014

C.L.

RPV
lower head

4.1m

5m

Triggersource2.75m
Figure 4.35: Grid geometry for the simulation of PWR ex-vessel steam explosion.

76 −
− 76

JAEA-Data/Code 2008-014

t= 0.200st= 0.500s

t= 1.000s

t= 2.000s

Void
scale
1.0
Melt
particles

Melt jet

0.8

0.6

0.4
Frozenparticles

0.2

Melt
pool

0.0

Figure 4.36: Proﬁles of the melt components (jet, pool and particles) and the void fraction in
the premixing stage of PWR ex-vessel steam explosion simulation.

Melt total

MoltenMolten, <0.75

12000

10000

Mass (kg)

8000

6000

4000

2000

0

0

0.5

1

1.5

2

2.5

3

Time (s)Figure 4.37: Calculated melt mass (total, molten and molten-premixed) in the premixing of
PWR ex-vessel steam explosion simulation.

77 −
− 77

JAEA-Data/Code 2008-014

1.8e+08

1.6e+08t= 0.0008s

t= 0.0000s

t= 0.0016s

1.4e+08

1.2e+08

1e+08

8e+07

6e+07

4e+07

r

2e+07

0

z
1.8e+08

1.6e+08t= 0.0032s

Pressure
scale (Pa)t= 0.0048s

t= 0.0064s

1.4e+08

1.2e+08

1e+08

8e+07

6e+07

1.8e+08
1.6e+08
1.4e+08
1.2e+08
1e+08
8e+07
6e+07
4e+07
2e+07
0

4e+07

2e+07

0

Figure 4.38: Proﬁles of the pressure wave propagation in the explosion stage of PWR ex-vessel
steam explosion simulation.
Explosion calculation result
Figure 4.38 shows the pressure wave propagation obtained by the explosion calculation. Figure 4.39 shows the evolution of the total ﬂuid kinetic energy. In this case, about 230MJ of the
kinetic energy is reached at about 13ms after the triggering.
As output information of the explosion calculation, the impulse on wall structures can bealso obtained by integrating the pressure in the cells adjacent to the wall. The impulse would
be more useful for the evaluation of the integrity or damage of the walls than the kinetic energyof the ﬂuid.[42, 43]− 78 −

78

JAEA-Data/Code 2008-014

250

Kinetic energy (MJ)

200

150

100

50

0

0

0.01

0.02

0.03

0.04

0.05

Time (s)Figure 4.39: Calculated total kinetic energy of ﬂuids in the explosion stage of PWR ex-vessel
steam explosion simulation.

79 −
− 79

JAEA-Data/Code 2008-014

5 Summary
An FCI simulation code JASMINE was developed. The concept, models, numerical solution
methods, veriﬁcation and example calculations were described in this report. Also, some model
details and practical information are collected in appendixes.
Some papers have already been published on application of JASMINE code for simulationof experiments[39], risk assessment for LWRs[42], and on a strategy for the application of FCI
simulation codes for plant risk assessment[43]. Those papers also include helpful information
for the usage of JASMINE code. Some of the model parameter setting have been modiﬁed fromprevious works. The setting described in this report is the latest one compatible with JASMINE
version 3.3a.
Present limitations of the models and calculation methods, and known practical problems
are summarized below.
Limitations, cautions and known problems:• JASMINE code is, at present, applicable only for the assessment of steam explosion loads.
The ability to simulate debris bed formation process and to assess long term coolability ofthe debris bed or melt pool is still poor.
• The present modeling and veriﬁcation concept relies on the idea of limitation of the steam
explosion energetics due to void and freezing. So, the status and improvement of modelingtechniques and experimental data on those aspects should be watched carefully.
• Void fraction issue– The current ability of JASMINE code in evaluating void distribution in the premixing
stage for long term is not good. We tried tuning related parameters conservatively,
not to overestimate overall void fraction. To get more reliable simulation of void
distribution in wide range of situations, it seems to need a fundamental improvementof 3D multiphase mixing modeling.
– Inﬂuence of the void in the explosion process is considered with an intuitive attenuation factor, i.e. 1 → 0 for α = 0.3 → 0.75. However, such an eﬀect of void on the ﬁne
fragmentation of the melt has not been experimentally conﬁrmed at a fundamental
level. Relevant experimental work is desirable for providing a reliable basis.
• Limitation by the numerical framework– As long as we use the conventional multiphase ﬂow simulation method with spatiallyaveraged volume fractions in relatively large cells, the void fraction in cells depends onthe cell size. Water in small cells are more quickly heated-up and produce void thanthat in large cells due to the diﬀerence of heat capacity. It aﬀects the attenuation offragmentation by void.
– At the prototypic range of melt temperature, ∼3000K, radiation dominates the heatrelease and freezing of the melt. The radiation heat transfer may extend beyond the
local cell when void fraction in the cell is large.1 The present modeling of radiation
1
Thermal radiation is absorbed by water mostly in several millimeters of depth. However, steam is transparent
to the radiation. Thus, the radiation from melt droplets in a cell ﬁlled with steam should travel to adjacent cells
where is water.

− 80 −

80

JAEA-Data/Code 2008-014

is limited within a single cell. Better radiation model considering heat deposition tosurrounding cells would improve the results.
– Geometrical arrangement: The melt jet is only available in the central column of theanalytical space. The one-dimensional modeling of the melt jet dose not allow strong
deformation of the jet leading edge (too much expansion in lateral direction causes a
numerical problem).
• Factors not modeled and that is relevant to the steam explosion process– Three dimensional eﬀects: an explosion initiated not at the center of the pool (3Dpressure wave dynamics)– Chemical eﬀects: hydrogen generation and chemical reaction heat due to oxidationof metallic components in the melt[35, 44]– Leading edge eﬀects in jet break-up• A care should be taken about extrapolation in scale, i.e. an application for plant scale(∼10 tons) means an extrapolation by two orders of magnitude in the melt mass from theveriﬁcation basis (∼100kg, in the FARO experiment). There can be unknown scale factor
that causes additional uncertainty in such extrapolation.
• Practical limitations– Inﬂexible physical property packages: The water and melt physical property packages
are not ﬂexible. The steam table package is totally hard coded; thus, usage of otherkind of coolant ﬂuid needs considerable code modiﬁcation. Several kinds of melt
property packages are already made available, as separate source ﬁles. Change ofthe melt property packages needs re-linking of the object ﬁle of the melt property
package.
– Redundancy in the output data: The binary output ﬁles, i.e. plot ﬁle for two-phaseﬂow and melt dump ﬁle for the melt, include most of useful information. Other
text format ﬁles are redundant though they may be convenient for a quick check of
calculation results.
– Restart function is only available when the melt model is used.
– Coding style is old. (Fortran 77)− 81 −

81

JAEA-Data/Code 2008-014

Acknowledgments
The authors acknowledge valuable suggestions and discussions of the late Mr. N. Yamano, supports by Dr. A. Ohnuki and Dr. H. Akimoto who developed and provided two-phase ﬂow code
ACE3D, an excellent work in modifying ACE-3D code for JASMINE by Mr. H. Kamo, andvaluable discussions on FCI modeling and its applications with OECD/SERENA Phase-1 members. Part of the FARO/KROTOS experimental data was provided by European Commission
for SERENA Phase-1 participants.

82− 82 −

JAEA-Data/Code 2008-014

Appendix A Notation
Notation of major variables in this report is summarized here.
Alphabets
A
Bo
cp
C
CdampparCdampwal
Cent
Cf rg
Cvx
Cvzwt

: cross section, surface or interface area
: Bond number: speciﬁc heat: constant, tuning parameter: damping factor for particle-particle collision: damping factor for particle-wall collision: correction factor on the mass ﬂux of droplet entrainment from the melt jet: tuning parameter for the ﬁne fragmentation model: empirical constant on the entrainment x-velocity: empirical constant for the weight in averaging melt jet and ambient ﬂuid velocities to give the entrainment z-velocity
D, d: diameterD: diﬀusion coeﬃcientdbnd: boundary layer thickness around the melt jetDe: diameter of melt droplets entrained from the melt jet surface
df: ﬁne fragment diameter
Dmed: mass median diameter of melt droplets
Dsmx: maximum stable diameter of melt droplets
Dp: melt particle (droplet) diametere: internal energyf: friction factor, function
F: force, function
Fhtint , fhigkill , fhilkill : parameters for control of coolant evaporation in premixing stage
F hy: hydrodynamic drag force on a particle
fα: function of the void fractionFr: Froude numberg, g: gravitational acceleration
Gr: Grashof numberh: heat transfer coeﬃcient, heat transfer rate, pool hightH: height of the calculation domain, heat transfer coeﬃcientHpl: water pool depthj: diﬀusional mass ﬂux
K: friction factor
Kev: parameter for control of coolant evaporation in explosion stagel: length
Lbrk: jet breakup lengthM: molecular weightme: mass ﬂux of the entrainment from jet surface
mp: mass of a particle in a particle groupms: mass ﬂux falling on the pool surfacen, n + 1: (superscript) old or new time stepng
: number of non-condensible gas components
nhist: step number criterion for release of real particle groups− 83 −

83

JAEA-Data/Code 2008-014

Np
Npcr
Nu
p
ps
Pr
q
Q
q̇
r x , rz
rxcr
R
Ra
Re
t
T
u, v, w
v, v
V
w(x)
We
xp
x, z
zbnd

: number of particles in a particle group: criterion in number of particles for real particle-group generation: Nusselt number: pressure: pseudo partial pressure of steam: Prandtl number: surface heat ﬂux, heat input in a cell: heat: heat release rate
: half width of a particle group in the x and z direction: criterion in x-direction size for real particle-group generation: radius of the calculation domain, particle radius: Rayleigh number: Reynolds number: time: temperature: velocity components in x (r), y (θ) and z direction: velocity: ﬂuid volume in a cell, velocity, volume ﬂux: width of the calculation domain at position x (≡ xΘ): Weber number: position of the center of particle group (≡ (xp , zp )): horizontal and vertical coordinates: distance from the leading edge of the melt jet

Greek symbols
α
αpack
β
Γ
δ
∆hf g
∆t
∆x, ∆z
∆xmin
�
Θ
κ
λ
µ
ν
π
ρ
σ
σSB
τ

: void fraction, volume fraction: maximum packing ratio for particle (∼ 0.6)
: volume expansion coeﬃcient, an alternative independent variable (≡ αT ): mass generation rate
: boundary layer thickness, variance: latent heat of evaporation: time step: grid size for x, z coordinate: minimum x-direction grid size in the system: density ratio, emissivity: Azimuthal width (angle) of the calculation domain: thermal diﬀusion coeﬃcient: thermal conductivity: viscosity: kinetic viscosity: Ludolphian number: density: surface tension: Stefan-Boltzmann constant: time, time constantSuﬃxesaav: non-condensible gas, ambient ﬂuid
: average− 84 −

84

JAEA-Data/Code 2008-014

bnd
c
cr
d
e
f
i
i, j, k
J
l
liq
m
mp
p
P
r
rad
s
sat
sf
sol
sub
sup
v
w
∞

: boundary layer: continuous phase, core, center: critical, critical point: dispersed phase: entrainment: ﬁne fragment: initial value: dummy index for coordinates or components: melt jet: liquid water: liquidus point: melt, dummy index for components: melting point: melt particle: melt pool
: relative, e.g. between melt and coolant: radiation: steam: saturation: surface: solidus point: subcool: superheat: vapor: wall: inﬁnity, far enough distance− 85 −

85

JAEA-Data/Code 2008-014

Appendix B Particle Heat Conduction
Models
B.1

Surface Temperature Drop of Melt Particles

The concept of the model is illustrated in Fig. B.1. A particle with uniform initial temperature
Tc cooled at the surface is considered. A temperature proﬁle of quadratic function is assumed
in the thermal boundary layer,
Tc − Tsf{r − (R − δ)}2 ,δ2T (r) = Tc −(B.1)

where Tsf denotes the surface temperature; δ denotes the boundary layer thickness.
The heat loss corresponding to the temperature drop in the boundary layer, Q, is expressed
by
�
R

Q=

R−δ

4πr2 cρ(Tc − T )dr ,(B.2)

where c and ρ denote the speciﬁc heat and density, respectively. Substituting T with the assumed
temperature proﬁle, Eq. (B.1), gives4δQ = πR3 cρ(Tc − Tsf )
3
R

�

1
1−
2

�

δ
R

�

1
+
10

�

δ
R

�2 �

.

(B.3)

The surface heat ﬂux, q, should satisfy the continuity,q = −λ

�

dT
dr

�

= 2λ
r=R

Tc − Tsf
,
δ

(B.4)

and also the heat balance

dQ
.
(B.5)
dt
Assume that the surface heat ﬂux q is kept constant during the period the surface temperature drops, i.e. both δ and Tsf in Eq. (B.3) change. By substituting
4πR2 q =Tc − Tsf =

δq
2λ

(B.6)

obtained from Eq. (B.4) for Tc − Tsf in Eq. (B.3), and by taking derivative of Eq. (B.3) in terms
of t, we have
�
� �
� � �
dQ
4 2q
3 δ
1 δ 2 dδ= πR δ 1 −
+
,
(B.7)
dt
3
κ
4 R
5 R
dt
where κ = λ/(cρ) is the thermal diﬀusion coeﬃcient. By connecting Eq. (B.7) and Eq. (B.5),
we have a diﬀerential equation for the evolution of thermal boundary layer,
�

3
δ 1−
4

�

δ
R

�

1
+
5

�

δ
R

�2 �

dδ = 3κdt .

(B.8)

We can get Tc − Tsf by Eq. (B.6) with δ obtained by solving Eq. (B.8).
An expression for Tav − Tsf is obtained by the deﬁnition of average temperature,4
4 3
πR cρTav = πR3 cρTc − Q .
33− 86 −86(B.9)

JAEA-Data/Code 2008-014

Equations (B.3) and (B.9) gives
Tc =Tav − Fδ Tsf,
1 − Fδ

��where Fδ is a function of δ deﬁned by
δ
Fδ =
R

1
1−
2

δ
R

�

1
+
10

�

(B.10)

δR

�2 �(B.11).

Equations (B.6) and (B.10) lead an expression for the diﬀerence of the average and surface
temperatures,δq
Tav − Tsf =(B.12)(1 − Fδ ) .
2λParticle

Thermalboundary layer

Tc
Tav
Tsf

R

Figure B.1: The surface temperature drop model.

B.2

Fast Transient Heat Release from Fine Fragments

The concept of the model illustrated in Fig. B.2 is similar to the one for the surface temperature
consideration of the melt particles (Appendix B.1). However, in the present model, the global
cooling phase which comes after full development of the thermal boundary layer is also in thescope. Therefore, the present model consists of the following two phases.
• Phase I: Boundary layer development phase.
• Phase II: Global cooling phase.

In phase I, the model is basically the same as the one for the surface temperature model for
melt particles. The heat loss is expressed by the boundary layer thickness, δ, byδ
4
Q = πR3 cρ(Tc0 − Tsf )
3
R

�

1
1−
2

�

δ
R

�

1
+
10

�

δ
R

�2 �

.

(B.13)

The assumed temperature proﬁle in the boundary layer,
T (r) = Tc0 −Tc0 − Tsf{r − (R − δ)}2 ,δ2− 87 −

87(B.14)

JAEA-Data/Code 2008-014

gives the heat ﬂux at the surface by
q = 2λTc0 − Tsf.
δ

(B.15)In the present model, we assume constant surface temperature, Tsf = const., instead of constant surface heat ﬂux as we did in the previous section B.1. The constant surface temperature
is valid typically for the interfacial heat conduction of two objects in sudden contact.
From the heat balance at the surface,
4πR2 q =dQ,
dt(B.16)and Eqs. (B.13) and (B.15), we obtain the diﬀerential equation for the boundary layer development,
�
� � �
δ
3 δ 2
δ 1− +dδ = 6κdt .
(B.17)
R 10 R
Non-dimensional form of the equations are obtained by normalizing δ, t and Q by
x=

κt
δ
Q
,
, τ = 2 , and β =RRQ0(B.18)

where Q0 is the total heat inventory4
Q0 = πR3 cρ(Tc0 − Tsf ) .
3(B.19)The non-dimensional forms of Eqs. (B.13) and (B.17) are�

11
β = x 1 − x + x2
2
10

�

, and

(B.20)

3x 1 − x + x2 dx = 6dτ .
10(B.21)

��

√
We can solve x by numerical integration starting with initial value x0 = 12τ0 based on the ﬁrst
order analytical solution. During phase I process, the change of the non-dimensional parameters
is x = 0 ∼ 1, β = 0 ∼ β1 (≡ 3/5) and τ = 0 ∼ τ1 . The non-dimensional time of the end of phaseI, τ1 , is 0.040 according to the numerical solution.
During phase II, the center temperature, Tc decreases from Tc0 , while the boundary layer
thickness is ﬁxed at R. The additional heat loss during phase II is∆Q = Q − Q1 =8
πR3 cρ(Tc0 − Tc ) ,15(B.22)

where Q1 is the heat loss at the end of phase I, which corresponds to β1 . In this case, the surface
heat ﬂux balance, Eq. (B.16), leads the time evolution equation,dTc = −15(Tc − Tsf )

κ
dt .
R2

(B.23)

By normalizing the center temperature byθ=Tc − Tsf,
Tc0 − Tsf

(B.24)

we have normalized form of Eqs. (B.22) and (B.23),2β = β1 + (1 − θ) , and5dθ = −15θdτ .
− 88 −

88(B.25)(B.26)

JAEA-Data/Code 2008-014

Equation (B.26) is solved analytically, with the initial condition τ = τ1 and θ = 1.
θ = e−15(τ −τ1 ) .

(B.27)The non-dimensional parameters change by θ = 1 → 0, β = β1 → 1 and τ = τ1 → ∞ duringphase II.
Let us deﬁne the non-dimensional time constant for the heat release from a fragment, τf ,as non-dimensional time, τ , at which the heat loss fraction ,β, equals 1 − e−1 � 0.632. The
non-dimensional time constant, τf , is in phase II, and the value is about 0.046.
Fine fragment

Thermalboundary layer

Tc0

Tc0
Tc

Tsf

Tsf

R
R
Phase I: Boundary layer development Phase II: Global cooling

Figure B.2: The heat release model for ﬁne fragments.

− 89 −

89

JAEA-Data/Code 2008-014

Appendix C Brief Description of CIP
MethodC.1

Basic Concept of CIP Method

CIP—cubic interpolated pseudo particle—method was developed by Yabe and Aoki[31]. The
basic concept of this scheme is cubic interpolation of discrete data and handling the advection
term by moving the proﬁle like in Lagrange method. Consider a one-dimensional hyperbolic
equation
∂f
∂f
+v
(C.1)
= 0,
∂t
∂z
where advection velocity v is constant. Given a spatial proﬁle of the function f and the velocity
v at a certain time t, the value of f at a certain grid point zi after a short moment ∆t can beapproximated as a result of advection.
f (zi , t + ∆t) � f (zi − v∆t, t)(C.2)

In the CIP method, the values of f between grid points are interpolated with a cubic polynomial
Fi (z) = {(ai X + bi )X + fi� }X + fi ,(C.3)

where X = z − zi , and fi� is the spatial derivative of f at the point zi . The constants ai and bi
are derived by considering the continuity of f and f � at grid points i and i + 1 or i − 1 dependingon the direction of the velocity,
ai =
bi =
φ =

−2φ + fi�2 + fi�
,∆z 2
3φ − fi�2 − 2fi�,∆zfi2 − fi
,
∆z

(C.4)
(C.5)
(C.6)

where i2 = i + 1 when v ≤ 0 or i2 = i − 1 otherwise, and ∆z = zi2 − zi . Note that the grid is
deﬁned such as . . . zi−1 < zi < zi+1 . . ..
The controlling equation for the spatial derivative f � can be deduced from the original equation,
∂f �
∂v
∂f �
+v
= −f � .
(C.7)
∂t
∂z
∂z
In linear case where v is constant, the right hand side vanishes and it becomes identical withthe original equation.
Therefore, once the proﬁle of the function and its derivative is given at t as f and f � , theproﬁle at the next time step t + ∆t can be calculated using the constants ai and bi byfinew = Fi (zi − v∆t) = {(ai ξ + bi )ξ + fi� }ξ + fi ,fi�new=Fi� (zi − v∆t) = (3ai ξ + 2bi )ξ + fi� ,where ξ = −v∆t.

− 90 −

90

(C.8)
(C.9)

JAEA-Data/Code 2008-014

C.2

Improvement for Steep Steps—CIP1 Scheme

A steep edge of function f can be encountered at places such as leading or trailing edge of amelt jet, where the material suddenly begins or quits to exist. The basic scheme brings overor under-shoot at such places, due to the feature of cubic polynomial interpolation. A simple
counter measure for it is given in Ref.[31].
Assume that a certain grid point is located at the edge of the proﬁle, i.e. one side is almost
ﬂat and the other side has a steep change. The interpolated function proﬁle at the ﬂat side issigniﬁcantly diﬀerent from the reality and causes unphysical over- or under-shoot. To prevent
this feature, two diﬀerent derivatives are considered on the front and the back side of such an“edge point”, and the use of derivative based on the interpolated proﬁle is given up on the ﬂatside, but give the derivative by ﬁnite diﬀerence between two points instead. This scheme iscalled CIP1[31]. An index np is introduced to put a ﬂag on the point at the edge.
|fi − fi−1 ||fi−1 − fi−2 |< � and< � (gap on the front side) or|fi+1 − fi ||fi+1 − fi ||fi+1 − fi ||fi+2 − fi+1 |if< � and< � (gap on the back side) or|fi − fi−1 ||fi − fi−1 |otherwise.

np = −1 if(C.10)

np = 1(C.11)

np = 0(C.12)

The constant � should be a small value like 0.05 . The derivatives on both sides of the point i
� (front) and f � (back). They are initially set equal to f � and changed
is taken diﬀerently as fi,1ii,0in case an edge is detected as follows.
fi − fi−1
zi − zi−1f
i+1 − fi�
fi,1=zi+1 − zi
�
fi,0
=

if

np = −1,(C.13)

if

np = 1.

(C.14)(C.15)The improved scheme becomes as follows.
if ξ ≥ 0, i2 = i + 1

and i3 = 1,

(C.16)and i3 = 0.

(C.17)∆z = zi2 − zi ,fi2 − fi,φ =∆z
�−2φ + fi�2 ,1−i3 + fi,i
3
ai =
,
∆z 2
�
3φ − fi�2 ,1−i3 − 2fi,i
3
,
bi =
∆z
�
finew = {(ai ξ + bi )ξ + fi,i}ξ + fi ,3(C.18)

otherwisefi�new=

i2 = i − 1�(3ai ξ + 2bi )ξ + fi,i.
3

(C.19)(C.20)(C.21)(C.22)(C.23)

This CIP1 scheme is used in the melt jet model of JASMINE.

C.3

Non-Linear Equation

A non-linear hyperbolic equation with a source term g
∂f v
∂f
+
=g
∂t
∂z
− 91 −

91(C.24)

JAEA-Data/Code 2008-014

is transformed into a non-conservation form
∂f
∂v
∂f
+v
=g−f
≡ G.
∂t
∂z
∂z

(C.25)

Corresponding equation for the derivative is
∂f �
∂f �
∂v
+v
= G� − f � .
∂t
∂z
∂z

(C.26)

These equations are solved in two steps, so called advection and non-advection phases.
• Advection phase:

∂f
∂f
+v
=0
∂t
∂z
∂f �
∂f �
+v
=0
∂t
∂z

(C.27)∂f=G∂t(C.29)∂v∂f �= G� − f �∂t∂z(C.30)• Non-advection phase:(C.28)Starting from values at the old time step, ﬁrst, the advection phase is calculated by the method
described in the previous section, and an intermediate result is obtained. Then, the nonadvection phase is solved in any method such as ﬁnite-diﬀerence, and the ﬁnal result for the newtime step is obtained. The following equations are ﬁnite diﬀerence form of the non-advection
phase, with intermediate values denoted by asterisk * and new time step values denoted bysuperscript new.
finew − fi∗= Gi(C.31)∆tfi�new − fi�∗
∆t

=

�

�

∂v
∂z i
� �
new − f ∗ )∗new(fi+1 − fi+1 ) − (fi−1
i−1
�∗ ∂v
− fi
∆t(zi+1 − zi−1 )∂z i= G�i − fi�∗− 92 −

92(C.32)

JAEA-Data/Code 2008-014

Appendix D Steam Table Package
WRSTEAMTABWRSTEAMTAB is a rapid running steam table package that has table data generated by theformulas published in JSME (The Japan Society of Mechanical Engineers) steam table (1980
edition)[45, 46] and calculates physical properties at a given set of pressure and temperature
by interpolation of the table. The calculation for the interpolation is much faster than direct
calculation of the original formulas.
The steam table data consist of the values of density and enthalpy, and their derivatives
in terms of temperature and pressure. The same cubic polynomial interpolation scheme as in
the CIP method (Appendix. C) is used on the pressure-temperature plane. The table points
are shown in Fig. D.1. Because of the utilization of the derivatives at each point, even a
relatively coarse grid gives an adequate level of accuracy, Also, adoption of this interpolation
method automatically brings the consistency of the derivatives that is necessary for the numericalsolution of the two-phase ﬂow model.
The table data covers the temperature range 274 < Tl < 787(K) for water and 281 < Tg <999(K) for steam, and pressure range 103 < P < 108 (Pa). Extrapolation is provided at theedges of the table range either by linear extrapolation or ideal gas approximation to cover a
wider range, P ∼ 2 × 108 (Pa), Tl ∼ 800(K) and Tg ∼ 3000(K). Especially the high temperature
range is needed for steam explosion modeling, which involves high temperature melt.
The two-phase ﬂow model in JASMINE considers a two-phase medium even in the supercritical conditions. In other words, the conventional two-phase ﬂow model is extended into thesupercritical regime to cover the situation that may arise in the steam explosion simulation.
To support this treatment, the saturation line is extended into the supercritical regime, andliquid/steam properties are given keeping a consistency, i.e. slightly smaller density and slightly
higher enthalpy for steam than water. Both water and steam properties are made available
near the saturation line so that superheated water (superheat ∼ 10K) and supercooled steam(supercool ∼10K) is managed numerically.
The accuracy of the saturation temperature, liquid density, liquid enthalpy, steam density
and steam enthalpy are shown in Figs D.2–D.6. The very steep variation near the critical point
is intentionally smoothed because it can be harmful for stable numerical solution. Then, theaccuracy around that area is not good. However, an adequate accuracy is obtained in otherareas.

− 93 −

93

JAEA-Data/Code 2008-014

1100

Steam
Water
Sat.
Crit.

1000

900
Temperature (K)

800

700

600

500

400

300

200
1e+02 1e+03 1e+04 1e+05 1e+06 1e+07 1e+08
Pressure (Pa)Deviation in Tsat (K)Figure D.1: The table data points for WRSTEAMTAB.

4

2

0

-2

-4
10

3

10

4

5

10
Pressure (Pa)

10

6

10

7

Figure D.2: Deviation of the calculated saturation temperature by WRSTEAMTAB from the
JSME formula, Tsat (WRSTEAMTAB) − Tsat (JSME).

94 −
− 94

JAEA-Data/Code 2008-014

0.03

0.02

0.01

0.4

0

0.2

-0.01

0

-0.02

-0.2

-0.03

-0.4

700
3

10

600

104

5

10
6
Pressure (Pa) 10

400

10

7

10

8

500Temperature (K)300

Figure D.3: Relative deviation of the calculated water density by WRSTEAMTAB from the
JSME formula, (ρl (WRSTEAMTAB) − ρl (JSME))/ρl (JSME).

0.03

0.02

0.01

0.4

0

0.2

-0.01

0

-0.02

-0.2

-0.03

-0.4

700
3

10

600

104

5

10
6
Pressure (Pa) 10

400

10

7

10

8

500Temperature (K)300

Figure D.4: Relative deviation of the calculated water enthalpy by WRSTEAMTAB from the
JSME formula, (hl (WRSTEAMTAB) − hl (JSME))/hl (JSME).

95 −
− 95

JAEA-Data/Code 2008-014

0.03

0.02

0.01

0.4

0

0.2

-0.01

0

-0.02

-0.2

-0.4
3

10

104

5

10
6
Pressure (Pa) 10

10

7

10

8

-0.03
1000
900
800
700
600
500 Temperature (K)400300

Figure D.5: Relative deviation of the calculated steam density by WRSTEAMTAB from the
JSME formula, (ρg (WRSTEAMTAB) − ρg (JSME))/ρg (JSME).

0.03

0.02

0.01

0.4

0

0.2

-0.01

0

-0.02

-0.2

-0.4
3

10

104

5

10
6
Pressure (Pa) 10

10

7

10

8

-0.03
1000
900
800
700
600
500 Temperature (K)400300

Figure D.6: Relative deviation of the calculated steam enthalpy by WRSTEAMTAB from the
JSME formula, (hg (WRSTEAMTAB) − hg (JSME))/hg (JSME).

96 −
− 96

JAEA-Data/Code 2008-014

Appendix E Input Preparation and Code
Usage
E.1

Input File Format

The input ﬁle is prepared as a plain text ﬁle and read by the code line by line, then interpreted
according to the order of the data. The interpreter has preprocessing functions to do withcomments, inclusion of other ﬁles and echoing. In the following sections, the preprocessing
functions and order of the data in the input ﬁle are described.

E.1.1

Preprocessing

Line length Maximum length of lines in the input ﬁle is 140 characters. The ﬁrst 132 characters are interpreted as data and rest of the line is neglected as an ID ﬁeld.
Comments Comment lines can be included in the input ﬁle. Any characters following double
slashes “//” are ignored as comments.
//

this is a comment line.
35
12.6// anything after "//" is ignored.

Also, C-style commenting is available. That is, characters embraced by a pair of “/*” and “*/”is interpreted as a comment.
/* node */ 15 /* element */ 6
23.8 /*
this is a comment.
*/
56.7
Nesting of the “/* */” pair, as in the next example, is not allowed. The job is terminated if it
is detected.
/*/* this is an example of nesting. */ */Directives There are directives that give instructions to the code about the processing of the
input ﬁle. A directive is recognized by one of the following keywords at the head of a line.
• #include <file name> : To include a ﬁle at the position this directive appears. Nestingof including is not allowed.
• #listoff : To stop echoing the input data. (The input data is echoed into the list output
unit.)• #liston : To resume echoing the input data.

E.1.2

Input Variables and Their Ordering

An input ﬁle consists of following sections. The data are recognized according to the position(order) that they appears. In the followings, variable names, types, meanings, values to be given
etc. are listed. Values given in ( ) are default or recommended values.
− 97 −

97

JAEA-Data/Code 2008-014

Flag for the namelist options
The ﬁrst section is a ﬂag to instruct whether a namelist options (see below) are given or not.
Variable
iopt

Type
Int

Description
Flag on the usage of namelist options
0 = No use of namelist options
1 = Use of namelist options

Title
Title of the case is read as a 132-character string.
Variable
title

Type
Char

Description
Title of the case in a132 format

Namelist options
If iopt=1 is speciﬁed in the ﬁrst section, a namelist “option” is read. The input format inthis section conforms to the Fortran namelist input syntax. Variables deﬁned in the namelist“option” is as follows.
Variable
cvm
sdtim
sdstp
editop

Type
Real
Real
Int
Int

idebugInt

idbeosInt

istdiﬀInt

ddst
fhtintRealReal

tlsupqcut

RealfhigkillRealfhilkill

Real

Description
Coeﬃcient of the virtual mass term for two-phase ﬂow (0.0)
Time step for message (stdout) output (0.1[s])Step interval for message output (1000)
List output selection (0)0 = Primary variables only
1 = More variables
2 = All the available variables
Flag to control debug information output (0)0 = No output
1 = Output
Flag for EOS data output at bad convergence (0)0 = No output
1 = Output (A ﬁle “out.db” is written.)
Flag to control the steam pseudo-diﬀusion model for smoothing ofsteam partial pressure ﬁeld (1)0 = No use
1 = UseCoeﬃcient of steam pseudo-diﬀusion (0.1)
Fraction of the heat input from the melt that is directly deposited
to gas-liquid interface for evaporation (0.02)
Superheat of water at which deposition of heat from melt to bulk
water is cut-oﬀ (0.5[K])Attenuation factor for the heat transfer coeﬃcient between superheated steam and gas-liquid interface for cells having heat input
from melt (0.1)
Attenuation factor for the heat transfer coeﬃcient between subcooled water and gas-liquid interface for cells having heat input
from melt (0.001)− 98 −

98

JAEA-Data/Code 2008-014

fcondresReal

fevapresRealdtminqvbd

Realihighil0Int

taueqxplg

Realirsarla0Int

ici0
fcfmistfcfbubb
ieqdia

Int
Real
Real
Int

nftran

Real

icimj

Int

Attenuation factor for the gas-liquid interface heat transfer
when condensation takes place and the gas phase involves noncondensible gas, i.e. consideration on non-condensible gas accumulation near the interface (0.02)
Attenuation factor for the gas-liquid interface heat transfer whenevaporation takes place and the gas phase involves non-condensible
gas, i.e. consideration on non-condensible gas dilution near theinterface (0.5)
Lower limit for the denominator of qvbd/delt in the evaporation
term, where qvbd[J] is the heat input from melt in a cell for one
time step (1e-6[s]) (When time step is decreased in case of failure inconversion of Newton iteration, qvbd/delt increases and may makethe situation worse. This is a measure to limit the increase of thesource term and continue the calculation. Energy conservation is
maintained.)
Flag to control the heat transfer coeﬃcient between water and gasliquid interface (0)0: Default; Recommended mode for premixing calculation
1: hig=hil=0 (no heat transfer)
2: Recommended mode for explosion calculation; Limited heattransfer, only to relax supercooled steam or superheated water (unstable state) or too high steam superheat (Tg > (Tm + Tsat )/2)
Relaxation time of highly super heated steam when ihighil0=2 (1e3[s])1: Set rsa=rla=0 (no heat transfer between non-condensible gas
and water or steam) (0)1: Set cix=ciy=ciz=0 (no interface friction)(0)Attenuation factor of two-phase friction for droplet ﬂow (1.0)
Attenuation factor of two-phase friction for bubbly ﬂow (1.0)
Control of bubble/droplet size model in two-phase ﬂow (1)0: Evaluation by local relative velocity and W e (sensitive to velocity variation)
1: Use equilibrium size by balance of drag and gravity/buoyancy
Power of interpolation function for two-phase friction in the transition regime, α = 0.3 − 0.75; a larger value causes faster escape of
bubbles from the pool (1.0–4.0)(1.0)
n(> 0): Use separated (annular) ﬂow model for columns 1–n fromthe center (0)Time step section
Time step control variables are given in this section.
Variable
sttimeendtimdeltmxdeltmn

Type
Real
Real
Real
Real

Description
Start time [s]End time [s]Maximum time step [s]Minimum time step [s]− 99 −

99

JAEA-Data/Code 2008-014

Variable
edtim
pltim
httim

Type
Real
Real
Real

Description
List output time step [s]Plot (binary) output time step [s]History output time step [s]Next, parameters for the time step control is given. The time step is increased or decreased
according to the number of iterations necessary for conversion of the solution. The Courant
limit is also checked and the time step is kept smaller than it.
VariableitrincTypeInt

itrredIntfdtincfdtredReal
Real

Description
Number of iterations; time step is increased when conversion isreached earlier than this number (5)Number of iterations; time step is decreased when conversion is notreached at this number (10)Increase factor of time step (1.20)
Decrease factor of time step (0.95)

Conversion check condition
The maximum iteration number and convergence criterion are given in this section. If the
maximum iteration number for the pressure convergence is exceeded, the iteration is quited, thetime step is reduced to 80% and the pressure convergence is re-tried. When the relative error of
the pressure equation set becomes smaller than given criterion, it is considered that convergenceis achieved.
Variable
maxitrepscnv

Type
Int
Real

Description
Maximum limit iteration number (20)Conversion criterion by the relative residual (1e-4)

History output section
Number of cells for which history output is made, nhstcell, is speciﬁed ﬁrst. Then, the x(r)and z direction indexes for those cells, (ihst(i),khst(i)), are given. The indexes are read by thefollowing syntax.
read(iwkunt,*) (ihst(i),khst(i), i=1,nhstcell)
Variable
nhstcell

TypeInt

Description
Number of cells for which history output is made (max. 50)

Variableihst(i)khst(i)

Type
Int
Int

Descriptionx(r) direction index of history output cells
z direction index of history output cells

Water level output section
Number of x(r) positions where water level is detected and written in the history output, nwlevel,is speciﬁed ﬁrst. Then, x(r) indexes for those positions are given. The indexes, iwlev(i), areread by the following syntax.
read(iwkunt,*) (iwlev(i), i=1,nwlevel)
VariablenwlevelType
Int

Description
Number of x(r) positions where water level is detected and written
in the history output (max. 20)− 100 −

100

JAEA-Data/Code 2008-014

Variableiwlev(i)TypeInt

Descriptionx(r) indexes where water level output is madeThe water level detection is made by the following two methods, and written as leva and levb,
respectively.
(a) Seek cells from the top at speciﬁed x(r) position, and ﬁnd a cell where α < αlev is ﬁrst
attained. Assume the z index of this cell is k1 . Find the z position where α = αlev , bylinear interpolation between k1 and k1 + 1(b) Seek cells from the top at speciﬁed x(r) position, and ﬁnd a cell where α < αlev is ﬁrst
attained. Assume that there is the water surface in the cell, and ﬁnd the z position of thewater surface from α, the top and bottom positions of the cell.
The variable αlev (≡ 0.8) is the criterion to ﬁnd the water surface.
System deﬁnition section
Deﬁnition of the analytical system is made in this section. The mesh number in every direction,
number of non-condensible gas components and their indexes, usage of the melt model, coordinate system, boundary conditions and gravity control are given. For the boundary condition,wall, inlet with given velocity, or ﬂow out with given pressure can be speciﬁed. Or, all theboundaries without speciﬁc deﬁnition are regarded as slip walls.
Variable
nx
ny
nz

Type
Int
Int
Int

Description
x (r) direction mesh numbery (θ) direction mesh number (1)z direction mesh number

Variable
ngas

Type
Int

Description
Number of the non-condensible gas components(max. 4)

if ngas ≥ 1, the kinds of the gas components are given.
Variable Typekindgs(1) Int
kindgs(2) Int···

Description
The 1st gas component
The 2nd gas component(Repeat this ngas times.)
Indexes for gas kinds:
1 = air
2 = hydrogen3 = helium
4 = nitrogen
5 = argon

Variable
igeom

Type
Int

cyclicInt

nvgrav

Int

Description
Selection of the coordinate (0)0 = Cylindrical coordinate
1 = Cartesian coordinatey (θ) boundary condition control (0)0 = wall
1 = cyclic condition
Gravity control
0 = z direction gz = −9.807(m/s2 )1 = gx , gy and gz are given in gravity components section later− 101 −

101

JAEA-Data/Code 2008-014

Variable
igbcz(1)TypeInt

igbcz(2)Int

Description
z direction boundary condition at the bottom (k=0)
0 = wall1 = inlet condition
2 = ﬂow out conditionz direction boundary condition at the top (k=nz+1)
0 = wall1 = inlet condition
2 = ﬂow out condition

For the inlet boundary, constant velocity is given; for the ﬂow out boundary, constant pressure
is given.
The next section speciﬁes walls deﬁned inside the calculation domain. It is used to stop
coolant ﬂow in horizontal or vertical direction, useful for simulating vessel walls or structural
walls. The number of the walls are ﬁrst given and deﬁnition of the positions are to follow.
Note that the melt model does not know the walls deﬁned here. Conﬁnement of the zone
where the melt can move is set later in the melt section.
Variable
nxwallVariable
ixwall(i)kxwall(1,i)
kxwall(2,i)
Variable
nzwall

Type
Int
Type
Int
Int
Int
Type
Int

Variablekzwall(i)izwall(1,i)

Type
Int
Int

izwall(2,i)

Int

Description
Number of vertical internal walls (x(r) direction stops) (max. 20)
Descriptionx(r) index of internal wall i (cell boundary index)
z index of the bottom of internal wall i (cell center index)
z index of the top of internal wall i (cell center index)
Description
Number of horizontal internal walls (z direction stops)(max. 20)
Descriptionz index of internal wall i (cell boundary index)
x(r) index of the inner (left) end of internal wall i (cell center
index)
x(r) index of the outer (right) end of internal wall i (cell center
index)

Mesh data section
Cell boundary positions are given for z, x (r) and y (θ) coordinates. Though JASMINE onlyuses 2D cylindrical coordinate (r-z), θ direction cell size (of the only one grid) is needed. It
deﬁnes the azimuthal slice size of the system.
Variable
z

Type
Real

Description
z direction cell boundary positions (nz+1 data)

Variable
rad

Type
Real

Description
x (r) direction cell boundary positions (nx+1 data)

Variable
th

Type
Real

Descriptiony (θ) direction cell boundary positions (ny+1 data)− 102 −

102

JAEA-Data/Code 2008-014

Gravity components section
This section is needed when nvgrav = 1 is speciﬁed.
Variable
gc
gravx
gravy
gravz

Type
Real
Real
Real
Real

Description
Absolute value of the gravitational accelerationx (r) directional cosiney (θ) directional cosine
z directional cosine

Cell variable data section
Fluid volume in cells and area of cell surfaces are initially set equal to the geometrical volume
and area. (No obstacle objects) Volume fractions, pressures, temperatures and velocities in thecells are initially set uniform for the entire system, and modiﬁcations are made for speciﬁed
rectangular zones later.
Variable
alpn0
pn0
tvn0
tln0

Type
Real
Real
Real
Real

Description
Initial value of steam volume fraction
Initial value of pressure
Initial value of steam temperature
Initial value of water temperature

Variable
vvxr0
vvyt0
vvz0

Type
Real
Real
Real

Description
Initial value of x (r) direction gas velocity
Initial value of y (θ) direction gas velocity
Initial value of z direction gas velocity

Variable
vlxr0
vlyt0
vlz0

Type
Real
Real
Real

Description
Initial value of x (r) direction water velocity
Initial value of y (θ) direction water velocity
Initial value of z direction water velocity

Volume fractions and temperatures of non-condensible gases are given if ngas ≥ 1. The
following data is repeated for ngas times.
Variable
alp0a
ta0

Type
Real
Real

Description
Volume fraction of a non-condensible gas component
Temperature of a non-condensible gas component

The partial pressure of a non-condensible gas component, Pa , is expressed by
Pa =

αa
�
P
(αs + αa )(E.1)

where αs is the steam volume fraction, αa is the volume fraction of non-condensible gas component a, P is the total pressure.
Modiﬁcation of the initial values for material distribution in cells can be made as follows,
for arbitrarily deﬁned rectangular zones.
Variable Typenmatarea Int

Description
Number of rectangular zones for which material distribution is
modiﬁedIf nmatarea≥1, the following speciﬁcation is given for every rectangular zone.
− 103 −

103

JAEA-Data/Code 2008-014

Variable
is

Type
Int

Descriptionx(r) index for the cell in one end of the diagonal of the rectangularzone (start point)z index for the cell in one end of the diagonal of the rectangularzone (start point)x(r) index for the cell in another end of the diagonal of the rectangular zone (end point)z index for the cell in another end of the diagonal of the rectangularzone (end point)

ks

Int

ie

Int

ke

Int

Variable
alpst
tv
tl

Type
Real
Real
Real

Description
Steam volume fraction in the zone
Steam temperature in the zone
Water temperature in the zone

Variable
alpgs(i)tvgs(i)

Type
Real
Real

Description
Volume fraction of non-condensible gas component i in the zone
Temperature of non-condensible gas component i in the zone

These data for non-condensible gas components are given ngas sets in the order of the gas index.
Boundary condition section
If igbcz is not 0 in the system section, boundary conditions—inlet or ﬂow out—are speciﬁed at
the bottom and top of the analytical domain. The bottom boundary is speciﬁed ﬁrst and the
top boundary follows.
The method of boundary condition speciﬁcation is as follows, and common for both thebottom and top.
The next variable is needed if igbcz= 2 (ﬂow out condition). The size and volume of theboundary cell is calculated by this number. If ibcgz= 1 (inlet condition), this data is not read,and the boundary cell size is set as same as the adjacent cell inside of the system.
Variable
dz

Type
Real

Description
z direction size of the boundary cell

The following volume fractions, temperatures and velocities are necessary regardless of thevalue of igbcz.
If positive values are given for the volume fraction, pressure and temperatures, those values
are set for all the cells along the boundary. If negative value(s) are given for one or morevariable(s), cell dependent values are read from the next section for those variables. Velocities
are always regarded cell dependent.
Variable
alpha
press
tempv
templ

Type
Real
Real
Real
Real

Description
Steam volume fraction
Pressure
Steam temperature
Water temperature

If a negative value is set for alpha, steam volume fractions for every cell along the boundary
is read by the following. The reading format is the same as that for velocities.
Variable
alpha

Type
Real

Description
Steam volume fraction− 104 −

104

JAEA-Data/Code 2008-014

It follows in the same way for pressure and temperatures.
Variable
press

Type
Real

Description
Pressure

Variable
tempv

Type
Real

Description
Steam temperature

Variable
templ

Type
Real

Description
Water temperature

Velocities are read as 2D data by the following format.

100

do 100 j = 1 , nyread(iwkunt,*) ( vel(i,j,k) , i = 1 , nx )continue

If the boundary is the inlet (bottom), k = 0, or if it is the outlet (top), k = nz. The velocities
read here override those read in the cell data section. In case of inlet boundary, the velocities
read here become the inlet velocities.
Variable
vvz

Type
Real

Description
Gas velocity

Variable
vlz

Type
Real

Description
Water velocity

If non-condensible gas exists, i.e. ngas≥ 1, their volume fractions and temperatures are alsoneeded for the boundary condition. If positive values are given for them, they are regarded
uniform along the boundary. If negative values are given, it is regarded cell dependent and readfrom the following sections.
Variable
alp0a
ta0

Type
Real
Real

Description
Gas volume fraction
Gas temperatureThis section is repeated for ngas times.
Melt input section
The next part is the input data for the melt model. The ﬁrst is the switching ﬂag of the usageof the melt model.
Variable
imelt

Type
Int

Description
Flag for the usage of the melt model (1)0 = Not use
1= Use

The following melt model related input data are necessary to be ﬁlled even though imelt=0
is speciﬁed. The ﬁrst section is about the mesh and boundary.

− 105 −

105

JAEA-Data/Code 2008-014

Variableprmpp iedg

Type
Int

prmj ktop

Intprmj nsubInt

Descriptionx(r) index (coolant ﬂow cell boundary index) of the right(outer) end of the zone in which the melt components(pool and particles) movez index (coolant ﬂow cell boundary index) of the melt jet
inlet position, that is also the top end of the zone of melt
components movement
Number of sub-division of the z direction coolant ﬂow cells
for the ﬁner grid for melt jet (4)The next section speciﬁes the melt jet inlet condition.
Variableprmj tempini
prmj nini

Type
Real
Int

Description
Initial melt jet temperature (const.)
Number of the time dependent data for the velocity anddiameter of the melt jet (max. mjetini=500)

Time dependent melt jet velocity and diameter data follows. They are read by the following
format.
read(iwkunt,*) (prmj_ini(i), i=1,prmj_nini)So, prmj nini pieces of values are necessary for each of the time, velocity and diameter.
Variableprmj tini(i)

Type
Real

Description
Time (s) at which the melt jet inlet velocity and diameteris given

Variableprmj vini(i)

Type
Real

Description
Melt jet inlet velocity at times speciﬁed in prmj tini(i)Variableprmj dini(i)

Type
Real

Description
Melt jet inlet diameter at times speciﬁed in prmj tini(i)The melt jet inlet velocity and diameter are interpolated or extrapolated for time by thefollowing scheme.
• For time between prmj tini(1) and prmj tini(prmj nini) : linear interpolation between thepair of data just before and after the time wanted• For time before prmj tini(1): the value at prmj tini(1)• For time after prmj tini(prmj nini): the value at prmj tini(prmj nini)
A function to put melt particle groups at arbitrary positions at arbitrary times is equipped.
One data set for this purpose includes position, attributes, the time to start throwing-in, thefrequency and interval of throwing-in. The number of the data sets is to be given ﬁrst.
Variable
nparinType
Int

Description
Number of the data sets for arbitrary particle group
throwing-in

Then, nparin sets of data are given. One data set includes data for the following variables. Careshould be taken not to give physically inconsistent particle group size, particle diameter andnumber of particles in the group, e.g. too densely packed group.
− 106 −

106

JAEA-Data/Code 2008-014

Variable
x
z
vx
vz
rgx
rgz
dp
tav

Type
Real
Real
Real
Real
Real
Real
Real
Real

np
tstartfrequency
tint

Real
Real
Int
Real

Description
x (r) position (coordinate) of the center of the groupz position (coordinate) of the center of the groupx (r) direction velocity
z direction velocity
1/2 size of the group in x (r) direction
1/2 size of the group in z direction
Particle diameter
Average temperature of the particle (the same temperature is assumed for the particle surface)
Number of particles in the group
Time to start throwing-in the particle group
Total frequency of throwing-in the particle group
Interval of throwing-in

Next, parameters for the constitutive equations in the melt model are given.
Variableprmj brkmdl

Type
Int

prmj sacvx
prmj sacvzwt
prmj sacentprmj saedia

Real
Real
Real
Real

prmj sacfrc

Real

Variable
prmpa mergecrt

Type
Real

prmpa cdamppar
prmpa cdampwal

Real
Real

prmpa crnump

Real

Description
Selection of the jet break-up length model (0=Saito etal. (local water density, in eﬀect under the water level),
1=Saito-type2 (cross section averaged two-phase density
for ambient ﬂuid density, cut-oﬀ if it is less than 1/1000
of melt density), 2=Saito-type3 (ambient ﬂuid density deﬁned according to cross section averaged void fraction αav :max. two-phase density in the radial proﬁle (∼water density) if αav < 0.5; cross section averaged two-phase density
if 0.75 < αav with cut-oﬀ if it is less than 1/1000 of melt
density; linear interpolation between them), 3=JAERI
model (the same ambient ﬂuid density with Saito-type3,
use Taylor type or Saito type according to BoJ , with limitation in the valid range for Saito type correlation) (3)Cvx (5.0)
Cvzwt (0.5)
Cent (1.0)
De 2–6mm for corium, ∼10mm for alumina; should begiven according to experimental results, or a certain range
is surveyed parametrically (involving values causing maximum explosion load)
Cf rc (1.0)

Description
Criterion for matching (relative diﬀerence) for merge ofparticle groups (0.2)
Damping factor for collision of two particle groups (0.5)
Damping factor for collision of a particle group and a wall
(0.1)
Criterion for the number of particles for production of aparticle group (10–5000)− 107 −

107

JAEA-Data/Code 2008-014

prmpa crhistIntprmpa ccfrc

Realprmpa chtc

Realprmpa nattrad

Realprmpa coxd

Realprmpa mh2max

Realprmpa fhtstl

Realprmpa dminparRealprmpa ivxranIntprmpa ibrkcriIntprmpa inotdrp

Intprmpa imrgtsf

Intprmpa ihtpackIntprmpa inomerge

Int

Variable
prmp inoht

Type
Int

prmp inotdrp

Intprmp inopamb

Int

Criterion for the number of steps for production of a particle group (1000)
Correction factor for particle-ﬂuid friction factor based on
the solid sphere correlation (1.0)
Correction factor for particle-ﬂuid heat transfer coeﬃcient(also used for the melt pool heat transfer) (1.0)
Exponent of the void attenuation factor for radiation heattransfer from melt particles; smaller values (< 1) let theradiation eﬀect persistent till higher void fractions (1.0)
Correction factor for particle oxidation model (1.0,
presently not implemented)
Hydrogen production for unit mass of melt [kg/kg] (2.4e-3,
presently not implemented)
Heat transfer degradation factor for settled particle groups(∼0.1)The min. limit particle size in the secondary break-up ofparticles in the premixing stage (1e-3[m])Flag for the random factor for x(r) velocity of particle
groups at release (0=no random factor, 1=enable random
factor) (1)Selection of temperature criterion for secondary break-up
of particles (−1)
0: Tav > Tm.p. (average temp. higher than melting point)
1: Tsf > Tsol. (surface temp. higher than solidus point)
2: Tsf > Tsol. or Tav > Tliq.
−1: Suppress secondary break-up
Flag for suppression of the particle surface temperaturedrop model (1=suppress) (0)Flag for usage of surface temperature criterion for particlepool merging; by enabling this, merge becomes less frequent and the surface area persists; if disabled, particle
average temperature criterion is used (1=use) (0)Flag for the collapse and heat transfer degradation forsettled particle groups (see prmpa fhtstl). (1=use) (1)
Flag for disabling merge between particle groups (e.g.
when tracking of a group is wanted) (1=disable) (0)Description
Flag for disabling melt pool heat transfer model (1=disable) (0)Flag for disabling melt pool surface temperature drop
model (1=disable) (0)Flag for disabling the ambient (ﬂuid) pressure term in themelt pool momentum equation (for stability) (1=disable)(1)

The followings are input data required for “explosion mode” calculations. When “explosionmode” is speciﬁed in a restart run, melt jet and melt pool components are converted into particle
groups, and explosion related constitutive models are enabled.
− 108 −

108

JAEA-Data/Code 2008-014

Variableprmf ixpl

Type
Int

Description
Flag for “explosion mode”
0 : OFF
1 : ON
2 : ON (conversion of jet/pool into particle groups is suppressed) (restart after an explosion mode run does nothave jet/pool components and this does not have eﬀect)

The following data are read only when prmf ixpl> 0 is speciﬁed.
Variableprmf ifrgmdl

Type
Int

prmf cfrg
prmf ifrgcri

Real
Int

pmrf ifrgattIntprmf ikevmdl

Intprmf ckevfrg

Realprmf cvicfrg

Realprmf cqfrg
prmf ttriglife

Real
Real

prmf ptrig

Realprmf nmtrig

Intprmf dfrg
prmf dmminfrg

Real
Real

Description
Selection of fragmentation model (0)0 : Caracharios(1983) model
1 : Yuen(1994) model
Tuning factor for fragmentation model (0.35)
Flag for temperature criterion for fragmentation (0)0: Tav > Tm.p. (average temp. higher than melting point)
1: Tsf > Tsol. (surface temp. higher than solidus point)
2: Tsf > Tsol. or Tav > Tliq.
Flag for the void attenuation factor for fragmentation (2)0: Two-phase ﬂuid-particle interaction scheme
1: Interaction with water, attenuation by 1 − α
2: Interaction with water, cut-oﬀ in the range α = 0.3 −0.75
Flag for the model for partition of heat released from ﬁne
fragments to evaporation/bulk heating (0)0 : Given by constant (=prmf kevfrg)
1 : Original model (not tested)
Fraction of heat from ﬁne fragments assigned for evaporation (0.7) (used if prmf ikevmdl=0)
Constant for the original model (not tested) (used if
prmf ikevmdl=1)
Factor on the heat release rate of ﬁne fragments (1.0)
Period during which particle fragmentation is enabled after the local triggering (passage of pressure front) (1e-3[s])Threshold for local triggering (cell pressure) (5e5–10e6[Pa], should be set higher than initial pressure andlower enough than pressure pulses by explosion)
Number of possible triggering for a particle group (fragmentation is enabled for a given period after triggering,and this sequence is repeated for the frequency given here)(< 5) (1)Size of ﬁne fragments (50e-6[m])The min. size of particles that can undergo ﬁne fragmentation (too small particles are not fragmented further) (if
0 is speciﬁed, the same as fragment size is set) (0)− 109 −

109

JAEA-Data/Code 2008-014

The recommended values for explosion model parameters given above are found as a resultof test calculations. For users who try this tuning by themselves, some generic features of theparameters summarized below would give a perspective.
• prmf cfrg : The peak of the pressure pulse changes.
• prmf ttriglife : The width of the pressure pulse changes.
• prmf ckevfrg : Smaller values for this parameter make less part of heat released from ﬁne
fragments goes to the evaporation, and pressurization. Thus, smaller pressure pulses andmechanical energy outputs arise for production of unit mass of ﬁne fragments.
To ﬁt the pressure pulse intensity and width (or mechanical energy output), tune prmf cfrg andprmf ttriglife to get better ﬁt. When good agreement is obtained for the pressure pulse while
the comparison of the ﬁne debris mass (e.g. mass of debris ﬁner than ∼0.1mm in experiments)
and the calculated fragment mass is not good, tuning of prmf ckevfrg should be tried. To keep
the same size of the pressure pulse and to get more fragment mass, a larger value for prmf cfrg
and a smaller value for prmf ckevfrg should be tried.

E.2

Running the Code

E.2.1

Normal Run

Instructions for arrangement of the source tree and compilation are given in electronic texts
included in the distributed source package. Execution of JASMINE on a Unix-like enviroment
is as follows.
jasmine.{PROP} -i input_file -o output_file -p plot_file -h hist_file \[-r restart_info_file]
We recommend to name the load module (executable ﬁle) “jasmine.{PROP}” where {PROP} isthe name of the melt physical property package. This is a counter measure of the inﬂexible
handling of the physical property packages, i.e. we need to make separate load modules fordiﬀerent melt physical property packages. The ﬁles speciﬁed in the command line are as follows.
input ﬁle Input data ﬁleoutput ﬁle List output ﬁle for coolant two-phase ﬂow
plot ﬁle Plot (binary) output ﬁle for coolant two-phase ﬂow
hist ﬁle History output ﬁle for coolant two-phase ﬂow
The input ﬁle can not be omitted. If the list output ﬁle is omitted (including -o), ./work/outlist
is used as default. If the plot output ﬁle is omitted (including -p), ./work/plotfile is usedas default. If the history output ﬁle is omitted (including -h), ./work/histfile is used as
default. To use these default output ﬁle names, it is necessary to have a work directory ./workin advance.
The option -r speciﬁes a restart calculation, as explained later.

− 110 −

110

JAEA-Data/Code 2008-014

E.2.2

Output FilesThe output ﬁles include both the coolant two-phase ﬂow data (output from ACE3D) and melt
data (output from the melt package). The ACD3D output ﬁle names can be given by thecommand line options as described above. Their contents are as follows.
List output (output ﬁle) Cell variables, error messages and so on, are written in plain text.
The format of the list is an x(r) × z table for every variable.
Plot output (plot ﬁle) Major variables for the two-phase ﬂow are written in a binary (unformatted) ﬁle. Post process programs included in the JASMINE package, i.e. readplot
and mkrsdat, can read out the data from the plot ﬁle for speciﬁed time or step. The datain the plot ﬁle can be cut out for visualization by the post process tools or for restartcalculations.
History output (hist ﬁle) History of the total mass and energy of every component, waterlevel, etc. and major cell variables for cells speciﬁed in the input ﬁle are written in acolumn formatted text ﬁle., i.e. one line for one output step. Lines other than data areheaded by #. The data in this ﬁle can be read and plotted by any plotting tool like
gnuplot[47]. Note that this ﬁle has large number of columns, i.e. ∼ 22 + 2Nwatlev + 12Nhist(Nwatlev and Nhist are the number of columns for water level output and the number ofcells for history output, respectively); thus, some software may fail in reading and plotting
the data.
Message output Short messages to show progress of the calculation or warnings are printed
on the screen if not redirected into a ﬁle.
The outputs from the melt package are written in ﬁles whose names are hard corded, as
follows.
Jet list output (out.j) List of melt jet variables in a text ﬁle, written at the same timing withthe two-phase ﬂow list output.
Pool list output (out.p) List of melt pool variables in a text ﬁle, written at the same timingwith the two-phase ﬂow list output.
Particles list output (out.par) List of particle group variables in a text ﬁle, written at thesame timing with the two-phase ﬂow list output.
Fragments list output (out.frg) List of fragment group variables in a text ﬁle, written at
the same timing with the two-phase ﬂow list output.
Melt history output (out.mlt) History of general information for the melt model, i.e. total
melt mass, internal energy, mass fractions of jet, pool and particles, etc. in a column
formatted text ﬁle. The output timing is the same with two-phase ﬂow history output.
Melt-two-phase ﬂow exchange output (out.m2f ) List of data for exchanges between melt
and two-phase ﬂow models (volume, force and heat) in a text ﬁle, written at the sametiming with the two-phase ﬂow list output.
Melt dump (out.mdp) Major variables for the melt model are written in a binary (unformatted) ﬁle, written at the same timing with the two-phase ﬂow plot output. This ﬁle isprocessed by the post process program, mkrsdat, for the preparation of a restart data forthe melt model.

− 111 −

111

JAEA-Data/Code 2008-014

E.2.3

Restart Run

When a calculation is failed (mainly by getting into a too small time step due to bad convergence), the calculation can be restarted at some steps before the failure, possibly with modiﬁed
input parameters or ﬁeld variables.
For a restart run, a restart data ﬁle should be made from the two-phase ﬂow plot ﬁle andthe melt dump ﬁle generated by the previous calculation. The restart data ﬁle is read as an
initial condition for the restart run. Also, arbitrary change of cell variables can be made at thestart of the restart run. The procedure is described below.
1. Make a new directory for the restart calculation data.
It is because a restart run writes a new set of output ﬁles. If a restart run is executed in
the same directory with the original calculation, some ﬁles (melt outputs) are overwritten
and lost.
2. Copy the original input ﬁle into the new directory.
This becomes the input ﬁle for the restart run, and it may be modiﬁed as long as theconsistency of the physical conditions are kept. Assume the ﬁle name is in.
3. Make a restart data ﬁle.
A post process program, mkrsdat, cuts data out of the two-phase ﬂow plot and the melt
dump ﬁles, and produces a “restart data ﬁle”. Assume that we execute the restart run in
the current directory, the previous calculation data are in another directory named ../a1,
and that the names of the two-phase ﬂow plot ﬁle and melt dump ﬁle are pl and out.mdp,
respectively. First, do the following to see the steps of the data included in the plot andmelt dump ﬁles.
mkrsdat ../a1/pl ../a1/out.mdp
Assume that we want the data at the step 1234. Then, -s option speciﬁes that, as follows.
mkrsdat ../a1/pl ../a1/out.mdp -s 1234
The restart data ﬁle is a binary ﬁle and the name is set rs.dat by default. The detailed
usage of mkrsdat command is printed by calling it without options or arguments.
4. Make a restart instruction ﬁle.
A restart instruction ﬁle is a text ﬁle where instructions needed in performing a restartrun is written. The format is as follows.
//// restart instruction file: example////
Description//--------------------------------//// makes comments that is ignored// restart data filers.datrestart data file name// modifications// ndatmod
8
number of data modification items// namei k value
instructions for data modifications− 112 −

112

JAEA-Data/Code 2008-014

pn
alpstnalpgsn1//
tvn
tln
tvgsn1
//
//
time
step

3
3
3

3
3
3

2.4e5
0.100
0.100

3
3
3

3
3
3

410.
s
410.

0.
1

variable name, cell index(i,k), value
volume fraction of a non-condensible gas
alpgsn#, # is one digit int for gas no.
"s" means the saturation temperature
non-condensible gas 1 temperature
the following data do not have cell no.
time given at the start of restart run
step no. at the start of restart run

In case that we need to set the temperature in a cell at the saturation point, and, at thesame time, we need to modify the pressure, volume fraction for steam or non-condensible
gas in the cell, the modiﬁcation of pressure or volume fraction should be made before thespeciﬁcation of saturation temperature. It is because the calculation of the saturation
temperature needs those information. Assume the name of the ﬁle is rsi.
5. Execute the restart run.
Run JASMINE code by specifying the restart instruction ﬁle with -r option.
jasmine.{PROP} -i in -o output_file -p plot_file -h hist_file -r rsi
This restart function is helpful for eﬀective debugging. When a calculation is failed, we cantry a restart calculation from a few steps before the failure and with a smaller data output time
step to ﬁnd out the reason of the failure.

E.2.4

Restart Calculation in “explosion mode”

The basic procedure is the same as described in the previous section, E.2.3. Specifying prmf ixpl=1
in the input ﬁle makes the restart run an explosion mode calculation, that diﬀers from the normal
one in the following points.
• The melt jet and pool data are converted into particle groups. If we want another explosion
mode calculation after an explosion model calculation, this conversion is not in eﬀect
because there is no melt jet or pool data. (prmf ixpl=2 explicitly speciﬁes suppression ofthis conversion)• Constitutive models for the explosion mode are used in the melt model calculation. Those
are as follows.
– No heat transfer, break-up or merge for the particle groups are considered.
– Fine fragments are generated by hydrodynamic interaction of the melt particles andcoolant, and the quick heat release from the ﬁne fragments dominates the evaporationand heat-up of the surrounding coolant.
Additionally, the following modiﬁcation in the input ﬁle and the cell data is necessary.
• Modify the followings in the input ﬁle: specify ihighil0=2 (suppress evaporation or condensation due to the thermal equilibration between the gas and water); give a small enoughvalue for dtminqvbd, e.q. 1e−8s � normal dt.

− 113 −

113

JAEA-Data/Code 2008-014• Set the external triggering in the restart instruction ﬁle: set pressure and non-condensible
gas volume fraction for a certain cell where a high pressure bubble is assumed, with pVvalue matched with an experimental condition or in an adequate range. Note that weshould also assume a high temperature for the gas not to cause troubles due to temperaturedrop of the expanding gas. If the triggering by pressurized non-condensible gas is planned,
the non-condensible gas should be included from the premixing calculation.
• Set time steps adequate for the explosion mode: ∼ 2µs for the maximum time step for the
calculation, ∼ 0.1ms for the output time step.

− 114 −

114

JAEA-Data/Code 2008-014

References[1] A.W. Cronenberg and R. Benz. Vapor explosion phenomena with respect to nuclear reactor
safety assessment. NUREG/CR-0245, TREE-1242, Idaho National Engineering Laboratories, 1978.
[2] M.L. Corradini, B.J. Kim, and M.D. Oh. Vapor explosion in light water reactors: a reviewof theory and modeling. Progress in Nuclear Energy, 22(1), pp.1–117, 1988.
[3] R. Meignen, D. Magallon, K-H. Bang, G. Berthoud, S. Basu, M. Bürger, M. Buck, M.L.
Corradini, H. Jacobs, O. Melikhov, M. Naitoh, K. Moriyama, R. Sairanen, J-H. Song,N. Suh, and T.G. Thefanous. Comparative review of FCI computer models used in the
OECD-SEREBA program. Proceedings of International Congress on Advances in Nuclear
Power Plants (ICAPP’05), Seoul, Korea, 2005. (CD-ROM:Paper 5087).
[4] N. Yamano, Y. Maruyama, T. Kudo, A. Hidaka, and J. Sugimoto. Phenomenological studies
on melt-coolant interactions in the ALPHA program. Nuclear Engineering and Design, 155,
pp.369–389, 1995.
[5] K. Moriyama, N. Yamano, Y. Maruyama, T. Kudo, and J. Sugimoto. Study of premixing phase of steam explosion with jasmine code in ALPHA program. Proceedings of 4thInternational Conference on Nuclear Engineering, New Orleans, volume 1B, pp.903–915,
1996.
[6] A. Annunziato, C. Addabo, A. Yerkess, R. Silverii, W. Brewka, and G. Leva. FARO test
L-14 on fuel-coolant interaction and quenching– comparison report, volume I: Analysis of
the results – OECD/CSNI international standerd problem 39. NEA/CSNI/R(97)31/PartI, OECD/NEA/CSNI, 1998.
[7] Y. Yang, K. Moriyama, Y. Maruyama, H-S. Park, and J. Sugimoto. Propagation calculation for steam explosion with JASMINE-pro code. Proceedings of the 7th International
Conference on Nuclear Engineering, Tokyo (ICONE7), 1999. ICONE-7235.
[8] Y. Yang, S. Nilsuwankosit, K. Moriyama, Y. Maruyama, N. Nakamura, and K. Hashimoto.
JASMINE-pro: A computer code for the analysis of the propagation process in a steam
explosion–user’s manual. JAERI-Data/Code 2000-035, Japan Atomic Energy ResearchInstitute, 2000.
[9] D. Magallon, K-H. Bang, S. Basu, G. Berthoud, M. Bürger, M.L. Corradini, H. Jacobs,
R. Meignen, O. Melikhov, K. Moriyama, M. Naitoh, J-H. Song, N. Suh, and T.G. Thefanous.
Status of international programme SERENA on fuel-coolant interaction. Proceedings ofInternational Congress on Advances in Nuclear Power Plants (ICAPP’05), Seoul, Korea,
2005. (CD-ROM:Paper 5382).
[10] D. Magallon, K-H. Bang, S. Basu, G. Berthoud, M. Bürger, M.L. Corradini, H. Jacobs,
R. Meignen, O. Melikhov, K. Moriyama, M. Naitoh, J-H. Song, N. Suh, and T.G. Thefanous.
115− 115 −

JAEA-Data/Code 2008-014

Results pf phase-1 of OECD programme SERENA on fuel-coolant interaction. The 1st
European Review Meeting on Severe Accident Research (ERMSAR-2005), Aix-en-Provence,France, 2005.
[11] A. Ohnuki, H. Kamo, and H. Akimoto. Development of multidimensional two-ﬂuid model
code ACE-3D for evaluation of constitutive equations. JAERI-Data/Code 96-033, Japan
Atomic Energy Research Institute, 1996. [in Japanese].
[12] K. Moriyama, H. Nakamura, and Y. Maruyama. Analytical tool development for coarse
break-up of a molten jet in a deep water pool. Nuclear Engineering and Design, 236,
pp.2010–2025, 2006.
[13] K. Moriyama, Y. Maruyama, T. Usami, and H. Nakamura. Coarse break-up of a streamof oxide and steel melt in a water pool. JAERI-Research 2005-017, Japan Atomic Energy
Research Institute, 2005.
[14] G.I. Taylor. The dispersion of jets of metals of low melting point in water. G.K. Batchelor,
editor, The scientiﬁc papers of Sir Geoﬀrey Ingram Taylor, vol.3 Aerodynamics and the mechanics of projectiles and explosions, pp.304–305. Cambridge University Press, Cambridge,
1963.
[15] M. Saito, K. Sato, and S. Imahori. Experimental study on penetration behaviors of water jet
into freon-11 and liquid nitrogen. ANS Proc. National Heat Transfer Conference, Houston,
U.S., volume 3, pp.173–183, 1988.
[16] H. Schlichting. Boundary layer theory, chapter XXI. McGraw Hill, 7 edition, 1979.
[17] JSME, editor. JSME data book: Heat Transfer (4th edition), pp.68–70. JSME, 1994. [inJapanese].
[18] S.S. Kutateladze. Heat transfer in condensation and boiling. ACE-tr-3770, US AEC, 1952.
[19] N. Zuber et al. Int. Development in Heat Transfer, Part II, p.230. ASME, 1961. paper 27.
[20] P.J. Berenson. Film boiling heat transfer from a horizontal surface. ASME Journal of Heat
Transfer, 83, pp.351–358, 1961.
[21] L.A. Bromley, N.R. Leroy, and J.A. Robbers. Heat transfer in forced convection ﬁlm boiling.
Industrial and Engineering Chemistry, 45, pp.2639–2646, 1953.
[22] C. Liu and T.G. Theofanous. Film boiling on spheres in single- and two-phase ﬂows. part
I: Experimental studies. ANS Proceedings of National Heat Transfer Conference, Portland,
U.S., pp.34–47, 1995.
[23] R.B. Bird, W.E. Stewart, and E.N. Lightfoot. Transport Phenomena, chapter 6.3. John
Wiley & Sons, 1960.
[24] W.E. Ranz and W.R. Marshall. Evaporation from drops: part I. Chemical Engineering
Progress, 48, pp.141–146, 1952.
[25] Sa. Kondo, K. Konishi, K. Morita, N. Shirakawa, A. Furutani, and D. J. Brear. Experimental study on simulated molten jet-coolant interactions. Nuclear Engineering and Design,
155, pp.73–84, 1995.
[26] M. Pilch and C. A. Erdman. Use of breakup time data and velocity history data to predictthe maximum size of stable fragments for acceleration-induced breakup of a liquid drop.
International Journal of Multiphase Flow, 13(6), pp.741–757, 1987.
− 116 −

116

JAEA-Data/Code 2008-014[27] C. Caracharios, M. Burger, and H. Unger. A transient two-phase model to describe thermal
detonations based on hydrodynamic fragmentation. Proc. Int. Meeting on LWR Severe
Accident Evaluation, Cambridge, 1983.
[28] W.W. Yuen, X. Chen, and T.G. Theofanous. On the fundamental microinteractions that
support the propagation of steam explosions. Nuclear Engineering and Design, 146, pp.133–
146, 1994.
[29] I. Huhtiniemi, H. Hohmann, R. Faraoni, M. Field, R. Gambaretti, and K. Klein. KROTOS
38 to 44: data report. TN I.96.37, European Commission Joint Research Centre, Italy,
1996.
[30] D.R. Liles and J.H. Mahaﬀy. TRAC-PF1/MOD1: an advanced best-estimate computer
program for pressurized water reactor thermal-hydraulic analysis. NUREG/CR-3858, LA10157-MS, U. S. Nuclear Regulatory Commission, 1994.
[31] T. Yabe and T. Aoki. A universal solver for hyperbolic equations by cubic-polynomialinterpolation, I. one-dimensional solver. Computer Physics Communications, 66, pp.219–
232, 1991.
[32] S.V. Patankar. Numerical heat transfer and ﬂuid ﬂow, chapter 6. Hemisphere publishingcorporation, 1980.
[33] D.F. Fletcher. An improved mathematical model of melt/water detonations – i. model
formulation and example results. International Journal of Heat and Mass Transfer, 34(10),pp.2435–2448, 1991.
[34] R. Barrett, M. Berry, T.F. Chan, et al. Templates for the solution of linear systems:
Building blocks for iterative methods. Society for Industrial and Applied Mechanics (SIAM),
1993.
[35] D. Magallon, I. Huhtiniemi, and H. Hohmann. Lessons learnt from FARO/TERMOS corium
melt quenching experiments. Nuclear Engineering and Design, 189, pp.223–238, 1999.
[36] D. Magallon and I. Huhtiniemi. Corium melt quenching tests at low pressure and subcooled
water in FARO. Nuclear Engineering and Design, 204, pp.369–376, 2001.
[37] I. Huhtiniemi, D. Magallon, and H. Hohmann. Results of recent KROTOS FCI tests:
alumina versus corium melts. Nuclear Engineering and Design, 189, pp.379–389, 1999.
[38] D. Magallon and I. Huhtiniemi. Energetic event in fuel-coolant interaction test FARO L33. Proc. 9th International Conference on Nuclear Engineering (ICONE-9), Nice Acropolis,
France, 2001. (Paper no.285).
[39] K. Moriyama, H. Nakamura, and Y. Maruyama. Simulation of alumina and corium steam
explosion experiments with JASMINE v.3. Proceedings of the 6th International Conferenceon Nuclear Reactor Thermal Hydraulics, Operations and Safety (NUTHOS-6), Nara, Japan,
2004. (CD-ROM: Paper ID.264).
[40] R. Silverii and D. Magallon. FARO LWR programme: test L-33 data report. Technical
Note No.I.00.124, EU JRC Ispra, 2000.
[41] A. Annunziato, C. Addabbo, and D. Magallon. FARO test L-33 quick look report. TN
I.00.111, European Commission Joint Research Centre, Italy, 2000.

− 117 −

117

JAEA-Data/Code 2008-014[42] K. Moriyama, S. Takagi, K. Muramatsu, H. Nakamura, and Y. Maruyama. Evaluationof containment failure probability by ex-vessel steam explosion in Japanese LWR plants.
Journal of Nuclear Science and Technology, 43(7), pp.774–784, 2006.
[43] K. Moriyama and H. Nakamura. A strategy for the application of steam explosion codes toreactor analysis. Proc. Technical Meeting on Severe Accident and Accident Management,
Tokyo, Japan, 2006. (CD-ROM).
[44] D.H. Cho, D.R. Armstrong, and W.H. Gunther. Experiments on interactions between
zirconium-containing melt and water (ZREX): hydrogen generation and chemical augmentation of energetics. Proc. OECD/CSNI specialists Meeting on Fuel-Coolant Interactions, Tokai-mura, Japan (JAERI-Conf 97-011, NEA/CSNI/R(97)26)(Part II), pp.595–
607, 1997.
[45] The Japan Society of Mechanical Engineers, editor. JSME Steam Tables. The Japan Societyof Mechanical Engineers, 1981.
[46] K. Kobayashi, K. Namatame, M. Akimoto, M. Okazaki, S. Sasaki, K. Soda, and
K. Sanokawa. Digital computer subroutine STEAM for JSME steam table. JAERI-M
6967, Japan Atomic Energy Research Institute, 1977.
[47] http://www.gnuplot.info. Gnuplot is a portable open source plotting software.

− 118 −

118

࿖㓙න૏♽㧔5+㧕

⴫㧝㧚5+ၮᧄන૏
5+ၮᧄන૏
ၮᧄ㊂ฬ⒓⸥ภ㐳ߐࡔ ࡯ ࠻ ࡞ O
⾰㊂ ࠠࡠࠣ࡜ࡓ MI
ᤨ
㑆
⑽
U
㔚
ᵹࠕ ࡦ ࡍ ࠕ #ᾲജቇ᷷ᐲ ࠤ ࡞ ࡆ ࡦ ‛ ⾰ ㊂ࡕ࡞ OQNశᐲ ࠞ ࡦ ࠺ ࡜ EF⴫㧞㧚ၮᧄන૏ࠍ↪޿ߡ⴫ߐࠇࠆ5+⚵┙න૏ߩ଀
5+ၮᧄන૏
⚵┙㊂ฬ⒓⸥ภ㕙Ⓧᐔ ᣇ ࡔ ࡯ ࠻ ࡞O૕Ⓧ┙ ᴺ ࡔ ࡯ ࠻ ࡞Oㅦ ߐ 㧘 ㅦ ᐲࡔ ࡯ ࠻ ࡞ Ფ ⑽OUടㅦᐲࡔ ࡯ ࠻ ࡞ Ფ ⑽ Ფ ⑽
OU
ᵄ
ᢙᲤ
ࡔ
࡯
࠻
࡞
O
ኒ ᐲ 㧔 ⾰ ㊂ ኒ ᐲ 㧕 ࠠࡠࠣ࡜ࡓᲤ┙ᴺࡔ࡯࠻࡞ MIO
⾰㊂૕Ⓧ㧔Ყ૕Ⓧ㧕 ┙ᴺࡔ࡯࠻࡞Ფࠠࡠࠣ࡜ࡓ OMI㔚ᵹኒᐲ ࠕࡦࡍࠕᲤᐔᣇࡔ࡯࠻࡞
#O⏛ ⇇ ߩ ᒝ ߐࠕ ࡦ ࡍ ࠕ Ფ ࡔ ࡯ ࠻ ࡞#O㧔 ‛ ⾰ ㊂ ߩ 㧕 Ớ ᐲ ࡕ ࡞ Ფ ┙ ᣇ ࡔ ࡯ ࠻ ࡞ OQNOノᐲ ࠞ ࡦ ࠺ ࡜ Ფ ᐔ ᣇ ࡔ ࡯ ࠻ ࡞ EFOዮ᛬₸ 㧔ᢙ‫ߩޓ‬㧕‫ޓ‬㧝⴫㧡㧚5+ធ㗡⺆ਸ਼ᢙ‫ ޓ‬ធ㗡⺆⸥ภ ਸ਼ᢙ‫ ޓ‬ធ㗡⺆


࡛
࠲ 㨅
࠺
ࠪ




࠯
࠲
㨆
࠮ ࡦ ࠴




ࠛ ࠢ ࠨ 㧱
ࡒ
࡝




ࡍ
࠲ 㧼
ࡑࠗࠢࡠ




࠹
࡜ 㨀
࠽
ࡁ




ࠡ
ࠟ 㧳
ࡇ
ࠦ




ࡔ
ࠟ 㧹
ࡈࠚࡓ࠻




ࠠ
ࡠ 㨗
ࠕ
࠻




ࡋ ࠢ ࠻ 㨔
࠯ ࡊ ࠻




࠺
ࠞ FC
࡛ ࠢ ࠻



⸥ภ
F
E
O
v
P
R
H
C
\
[

⴫㧢㧚࿖㓙න૏♽ߣ૬↪ߐࠇࠆ߇࿖㓙න૏♽ߦዻߐߥ޿න૏

⴫㧟㧚࿕᦭ߩฬ⒓ߣߘߩ⁛⥄ߩ⸥ภߢ⴫ߐࠇࠆ5+⚵┙න૏
5+⚵┙න૏
⚵┙㊂
ઁߩ5+න૏ߦࠃࠆ 5+ၮᧄන૏ߦࠃࠆ
ฬ⒓
⸥ภ
⴫ߒᣇ
⴫ߒᣇ
C
TCF
ᐔ
㕙
ⷺ࡜ ࠫ ࠕ ࡦO࡮O D
C
E


┙
૕
ⷺ ࠬ࠹࡜ࠫࠕࡦUTO ࡮O  D
๟
ᵄ
ᢙࡋ
࡞
࠷
*\
U
ജ
࠾ ࡘ ࡯ ࠻ ࡦ
0O࡮MI࡮U࿶ ജ ᔕ ജࡄ ࠬ ࠞ ࡞
2C
0OO࡮MI࡮Uࠛࡀ࡞ࠡ࡯ ઀੐ ᾲ㊂ ࠫ ࡘ ࡯ ࡞,
0࡮O
O ࡮MI࡮UᎿ ₸  ᡼ ኿ ᧤ࡢ
࠶
࠻
9
,U
O࡮MI࡮U㔚 ⩄  㔚 ᳇ ㊂ࠢ ࡯ ࡠ ࡦ%U࡮#㔚૏Ꮕ㧔㔚࿶㧕 ⿠㔚ജ ࡏ
࡞
࠻
8
9#
O࡮MI࡮U࡮#㕒㔚ኈ㊂ࡈ ࠔ ࡜ ࠼(%8O࡮MI࡮U࡮#
㔚
᳇
ᛶ
᛫ࠝ
࡯
ࡓ
ǡ
8#
O࡮MI࡮U࡮#ࠦ ࡦ ࠳ ࠢ ࠲ ࡦ ࠬࠫ ࡯ ࡔ ࡦ ࠬ
5#8O࡮MI࡮U࡮#⏛᧤࠙ ࠛ ࡯ ࡃ
9D
8࡮UO࡮MI࡮U࡮#



⏛
᧤
ኒ
ᐲ࠹
ࠬ
࡜
6
9DO
MI࡮U ࡮#ࠗ ࡦ ࠳ ࠢ ࠲ ࡦ ࠬࡋ ࡦ ࡝ ࡯*
9D#O࡮MI࡮U࡮#࠮ ࡞ ࠪ ࠙ ࠬ ᷷ ᐲ ࠮࡞ࠪ࠙ࠬᐲ F͠శ᧤࡞ ࡯ ࡔ ࡦNO
EF࡮UT EO࡮O࡮EFEF
ᾖ
ᐲ࡞
ࠢ
ࠬ
NZ
NOO
O࡮O࡮EFO࡮EF㧔᡼኿ᕈᩭ⒳ߩ㧕᡼኿⢻ ࡌ ࠢ ࡟ ࡞$SUๆ ෼ ✢ ㊂  ⾰ ㊂ࠛ ࡀ ࡞
ࠣ
࡟
ࠗ
)[
,MI
O࡮Uࠡ ࡯ ಽ ਈ  ࠞ ࡯ ࡑ✢ ㊂ ᒰ ㊂  ๟ ㄝ✢ ㊂ ᒰ㊂ ᣇะᕈ✢㊂ᒰ㊂ ୘ ࠪ ࡯ ࡌ ࡞ ࠻
5X,MIO࡮U
ੱ✢㊂ᒰ㊂⚵❱✢㊂ᒰC ࡜ࠫࠕࡦ෸߮ࠬ࠹࡜ࠫࠕࡦߩ૶↪ߪ‫ޔ‬หߓᰴరߢ޽ߞߡ߽⇣ߥߞߚᕈ⾰ࠍ߽ߞߚ㊂ࠍ඙೎ߔࠆߣ߈ߩ⚵┙න૏ߩ⴫ߒᣇߣߒߡ೑ὐ߇޽ࠆ‫┙⚵ޕ‬න૏ࠍᒻ૞ࠆߣ߈ߩ޿ߊߟ߆ߩ
‫⴫ߪ଀↪ޓ‬㧠ߦ␜ߐࠇߡ޿ࠆ‫ޕ‬D ታ㓙ߦߪ‫⸥ߪߦᤨࠆߔ↪૶ޔ‬ภTCF෸߮UT߇↪޿ࠄࠇࠆ߇‫ޔ‬⠌ᘠߣߒߡ⚵┙න૏ߣߒߡ
ߩ⸥ภ̌㧝̍ߪ᣿␜ߐࠇߥ޿‫ޕ‬E ᷹శቇߢߪ‫ߩࡦࠕࠫ࡜࠹ࠬޔ‬ฬ⒓ߣ⸥ภUTࠍන૏ߩ⴫ߒᣇߩਛߦߘߩ߹߹⛽ᜬߒߡ޿ࠆ‫ޕ‬F ߎߩන૏ߪ‫ࠬ࠙ࠪ࡞࠮࡝ࡒߡߒߣ଀ޔ‬ᐲ㨙͠ߩࠃ߁ߦ5+ធ㗡⺆ࠍ઻ߞߡ↪޿ߡ߽⦟޿‫ޕ‬ᧄ࡟ࡐ࡯࠻ߪ⁛┙ⴕ᡽ᴺੱᣣᧄේሶജ⎇ⓥ㐿⊒ᯏ᭴߇ਇቯᦼߦ⊒ⴕߔࠆᚑᨐႎ๔ᦠߢߔ‫ޕ‬ᧄ࡟ࡐ࡯࠻ߩ౉ᚻਗ߮ߦ⪺૞ᮭ೑↪ߦ㑐ߔࠆ߅໧޿วࠊߖߪ‫ޔ‬ਅ⸥޽ߡߦ߅໧޿วࠊߖਅߐ޿‫ޕ‬ߥ߅‫ోߩ࠻࡯ࡐ࡟ᧄޔ‬ᢥߪᣣᧄේሶജ⎇ⓥ㐿⊒ᯏ᭴ࡎ࡯ࡓࡍ࡯ࠫ㧔http://www.jaea.go.jp㧕ࠃࠅ⊒ାߐࠇߡ޿߹ߔ‫ޕ‬⁛┙ⴕ᡽ᴺੱᣣᧄේሶജ⎇ⓥ㐿⊒ᯏ᭴ ⎇ⓥᛛⴚᖱႎㇱ ⎇ⓥᛛⴚᖱႎ⺖‫ޥ‬319-1195 ⨙ၔ⋵㇊⃗㇭᧲ᶏ᧛⊕ᣇ⊕ᩮ 2 ⇟࿾ 4
㔚⹤ 029-282-6387, Fax 029-282-5920, E-mail:ird-support@jaea.go.jp
This report is issued irregularly by Japan Atomic Energy Agency
Inquiries about availability and/or copyright of this report should be addressed toIntellectual Resources Section, Intellectual Resources Department,
Japan Atomic Energy Agency
2-4 Shirakata Shirane, Tokai-mura, Naka-gun, Ibaraki-ken 319-1195 Japan
Tel +81-29-282-6387, Fax +81-29-282-5920, E-mail:ird-support@jaea.go.jp© Japan Atomic Energy Agency, 2008⴫㧠㧚න૏ߩਛߦ࿕᦭ߩฬ⒓ߣߘߩ⁛⥄ߩ⸥ภࠍ฽߻5+⚵┙න૏ߩ଀
5+⚵┙න૏
ฬ⒓
⸥ภ
5+ၮᧄන૏ߦࠃࠆ⴫ߒᣇ☼ᐲࡄ ࠬ ࠞ ࡞ ⑽2C࡮UO ࡮MI࡮Uജ ߩ ࡕ ࡯ ࡔ ࡦ ࠻ ࠾ࡘ࡯࠻ࡦࡔ࡯࠻࡞0࡮OO ࡮MI࡮U⴫㕙ᒛജ ࠾ࡘ࡯࠻ࡦᲤࡔ࡯࠻࡞
0O
MI࡮U



ⷺ
ㅦ
ᐲ࡜ ࠫ ࠕ ࡦ Ფ ⑽TCFUO࡮O ࡮U Uⷺടㅦᐲ ࡜ ࠫ ࠕ ࡦ Ფ ᐔ ᣇ ⑽ TCFU O࡮O࡮UUᾲ ᵹ ኒ ᐲ  ᡼ ኿ ᾖ ᐲ ࡢ࠶࠻Ფᐔᣇࡔ࡯࠻࡞
9O
MI࡮U



ᾲ ኈ ㊂  ࠛ ࡦ ࠻ ࡠ ࡇ ࡯ ࠫࡘ࡯࡞Ფࠤ࡞ࡆࡦ
,O ࡮MI࡮U ࡮⾰ ㊂ ᾲ ኈ ㊂ 㧔 Ყ ᾲ ኈ ㊂ 㧕  ࠫࡘ࡯࡞Ფࠠࡠࠣ࡜ࡓ
, MI࡮- O࡮U࡮-
⾰ ㊂ ࠛ ࡦ ࠻ ࡠ ࡇ ࡯ Ფࠤ࡞ࡆࡦ
⾰ ㊂ ࠛ ࡀ ࡞ ࠡ ࡯ࠫࡘ࡯࡞Ფࠠࡠࠣ࡜ࡓ,MIO ࡮U ࡮㧔 Ყ ࠛ ࡀ ࡞ ࠡ ࡯ 㧕ࡢ࠶࠻Ფࡔ࡯࠻࡞Ფࠤ
ᾲવዉ₸
9 O࡮- O࡮MI࡮U࡮-࡞ࡆࡦࠫࡘ࡯࡞Ფ┙ᣇࡔ࡯࠻૕ Ⓧ ࠛ ࡀ ࡞ ࠡ ࡯
O ࡮MI࡮U
,O
࡞


㔚
⇇
ߩ
ᒝ
ߐࡏ࡞࠻Ფࡔ࡯࠻࡞
8OO࡮MI࡮U ࡮#ࠢ࡯ࡠࡦᲤ┙ᣇࡔ࡯࠻


૕
Ⓧ
㔚
⩄
O
࡮U࡮#
%O
࡞
ࠢ࡯ࡠࡦᲤᐔᣇࡔ࡯࠻


㔚
᳇
ᄌ
૏
O ࡮U࡮#
%O
࡞




⺃
㔚
₸ ࡈࠔ࡜࠼Ფࡔ࡯࠻࡞(OO ࡮MI ࡮U ࡮#


ㅘ
⏛
₸ ࡋࡦ࡝࡯Ფࡔ࡯࠻࡞*OO࡮MI࡮U ࡮#ࡕ ࡞ ࠛ ࡀ ࡞ ࠡ ࡯ࠫ ࡘ ࡯ ࡞ Ფ ࡕ ࡞,OQN
O ࡮MI࡮U ࡮OQN
ࡕ ࡞ ࠛ ࡦ ࠻ ࡠ ࡇ ࡯  ࠫࡘ࡯࡞Ფࡕ࡞Ფࠤ࡞, OQN࡮- O࡮MI࡮U࡮-࡮OQN
ࡕ
࡞
ᾲ
ኈ
㊂ ࡆࡦ

ᾖ ኿ ✢ ㊂ 㧔 㨄 ✢ ෸ ߮ ǫ ✢ 㧕 ࠢ࡯ࡠࡦᲤࠠࡠࠣ࡜ࡓ%MI
MI ࡮U࡮#


ๆ
෼
✢
㊂
₸ࠣ ࡟ ࠗ Ფ ⑽
)[U
O ࡮U





᡼
኿
ᒝ
ᐲ ࡢ࠶࠻Ფࠬ࠹࡜ࠫࠕࡦ9UTO ࡮O ࡮MI࡮U O ࡮MI࡮Uࡢ࠶࠻Ფᐔᣇࡔ࡯࠻࡞





᡼
኿
ノ
ᐲ
࡮UT
O
࡮O
࡮MI࡮U
MI࡮U
9
O
Ფࠬ࠹࡜ࠫࠕࡦ
⚵┙㊂

ฬ⒓
ಽ
ᤨ
ᣣ
ᐲ
ಽ
⑽
࡝࠶࠻࡞
࠻ࡦ
ࡀ࡯ࡄ
ࡌ࡞

⸥ภ
5+න૏ߦࠃࠆ୯
OKP OKPUJJOKPUF
FJUq q Ǹ TCF
̉ ̉  q Ǹ TCF
̍ ̍  ̉ Ǹ TCF
N‫ޔ‬. NFOOVV MI0R 0R$$  NP 0R

⴫㧣㧚࿖㓙න૏♽ߣ૬↪ߐࠇߎࠇߦዻߐߥ޿න૏ߢ
5+න૏ߢ⴫ߐࠇࠆᢙ୯߇ታ㛎⊛ߦᓧࠄࠇࠆ߽ߩ
ฬ⒓⸥ภ
5+න૏ߢ޽ࠄࠊߐࠇࠆᢙ୯㔚 ሶ ࡏ ࡞ ࠻ G8 G8  ,⛔৻ේሶ⾰㊂න૏ W W  MIᄤ ᢥ න ૏ WC WC  O⴫㧤㧚࿖㓙න૏♽ߦዻߐߥ޿߇࿖㓙න૏♽ߣ‫ޓޓޓ‬૬↪ߐࠇࠆߘߩઁߩන૏
ฬ⒓⸥ภ
5+න૏ߢ޽ࠄࠊߐࠇࠆᢙ୯
ᶏ㉿㧝ᶏ㉿Oࡁ࠶࠻㧝ࡁ࠶࠻㧝ᶏ㉿Ფᤨ  OUࠕ࡯࡞ C 㧝CFCOOࡋ ࠢ ࠲ ࡯ ࡞ JC 㧝JCJOOࡃ࡯࡞ DCT 㧝DCT/2CM2CJ2C2C
ࠝࡦࠣࠬ࠻ࡠ࡯ࡓ Έ 㧝ΈPOOࡃ࡯ࡦ D 㧝DHOO⴫㧥㧚࿕᦭ߩฬ⒓ࠍ฽߻%)5⚵┙න૏ฬ⒓⸥ภ 5+න૏ߢ޽ࠄࠊߐࠇࠆᢙ୯ࠛ࡞ࠣ GTI GTI ,࠳ࠗࡦ F[P F[P0

ࡐ
ࠕ
࠭
2
2F[P࡮UEO 2C࡮Uࠬ ࠻ ࡯ ࠢ ࠬ 5V 5VEOUOU

ࠟ
࠙
ࠬ
)
) 6ࠛ ࡞ ࠬ ࠹ ࠶ ࠼ 1G 1G Ǹ #Oࡑ ࠢ ࠬ ࠙ 㨴 ࡞ /Z /Z9Dࠬ࠴࡞ࡉ UD UDEFEOEFOࡎ࠻ RJ RJNZࠟ࡞ )CN )CNEOUOU

⴫㧚࿖㓙න૏ߦዻߐߥ޿ߘߩઁߩන૏ߩ଀
ฬ⒓⸥ภ
5+න૏ߢ޽ࠄࠊߐࠇࠆᢙ୯
ࠠࡘ࡝
࡯ %K %K$S࡟ ࡦ ࠻ ࠥ ࡦ 4 4%MI࡜࠼ TCF TCFE)[)[࡟
ࡓ TGO TGOE5X5X

:
✢
න
૏
:WPKV POࠟࡦࡑ ǫ ǫP66ࠫ ࡖ ࡦ ࠬ ࠠ ࡯ ,[ ,[9࡮Ox*\

ࡈ
ࠚ
࡞
ࡒ
HGTOKHO Oࡔ࡯࠻࡞♽ࠞ࡜࠶࠻
OGVTKEECTCVOI MI࠻࡞ 6QTT 6QTT  2C
ᮡ Ḱ ᄢ ᳇ ࿶ CVO CVO2C
ࠞ
ࡠ
࡝
࡯ ECN

ࡒ
ࠢ
ࡠ
ࡦ
v
vvO O㧔╙㧣㧘ᐕᡷ⸓㧕









 


































































































































 
















 





    




  

















  















































