平成 27 年度 修士論文

原子炉反応度事故解析の高度化に関する研究

原子炉工学研究室
松浦健太

北海道大学大学院 工学院
エネルギー環境システム専攻

目次
第１章

序論 ........................................................ 2

1.１

緒言 .............................................................. 2

1.2

国内における反応度投入事象 ........................................ 41.3

原子炉における各種フィードバック .................................. 7

1.4

研究目的 .......................................................... 8

第2章

解析手法 .................................................... 92.1

CBZ 概要 ........................................................... 9

2.2

EFP モデル ........................................................ 11

2.3

TRAC 概要 ......................................................... 13

2.4

結合計算概要 ..................................................... 16

第3章

燃料ピンセル過渡応答解析 ................................... 17

3.1

解析条件 ......................................................... 17

3.2

EFP モデルにおける解析 ............................................ 21

3.2.１

即発臨界未満の正の反応度を印加 ............................... 223.2.２

負の反応度を印加............................................. 25

3.2.３

未臨界体系への中性子パルス照射 ............................. 27

3.3

熱水力結合解析における解析 ....................................... 29

3.３.１

第4章

結論 ....................................................... 34

参考文献
付録

即発臨界以上の正の反応度印加 ................................. 30

37
39

1

第1章

序論

1.1 緒言世界では、石油枯渇への懸念や地球温暖化防止の観点から低炭素社会を目指した取り組
みが数多く行われている。その一つの選択肢として、二酸化炭素の排出が少ない原子力発電が注目され、エネルギー需要の著しいアジア諸国を中心に各国で新規建設・計画が行われている。2014 年 8 月時点では、全世界 31 ヶ国で合計 435 基の原子力発電所が稼働中であり、13 ヶ国で 72 基が建設中となっている[1]。現在日本では 48 基の原子力発電所があるものの、東日本大震災以降、多くが定期検査により停止中となっている。日本はエネルギー
資源に乏しく輸入に頼らざるを得ない。近年の中東情勢ならびに米国におけるシェールガ
ス革命、原油輸出解禁などにより火力発電燃料価格の不安定かつ予測困難な変動やエネル
ギーセキュリティーの問題が生じている。また、近年では再生可能エネルギーを用いた発電が実用化に向けて着々と計画されているものの、出力変動や発電特性から原子力に変わる代替エネルギー源とするには困難を伴う事が予測される。その為、2014 年 4 月に閣議決定されたエネルギー基本計画で述べられているように、原子力発電は「重要なベースロー
ド電源」として必要であると言える。
2011 年 3 月 11 日、東日本大震災による地震と津波により東京電力株式会社福島第一原子力発電所事故が発生した。この事故以降、事故以前にはシビアアクシデント対策が規制の対象とされず十分な備えがなかった事に加え、バックフィットを適用する法的仕組みがな
く、海外の知見を導入し不確実なリスクに対応して安全の向上を図る事はあくまで事業者の自主性に任されていたといった指摘を受け[2]、新たに新規制基準を策定する要求が高まってきた。その結果、福島第一原子力発電所の事故における教訓や、これまでは要求されてこなかった自然災害への対策や火災・内部溢水対策、シビアアクシデント対策やテロへの備えなど、様々な事項が追加された新規制基準が策定され、2013 年 7 月 8 日より施行された。これに伴い、各原子力発電所ではさらなる安全性確保の為に様々な過酷事故対策が行われ、事故時における安全性・信頼性は向上することとなった。
新規制基準で主に見直された事柄としてはシビアアクシデントに至る地震・津波の厳格な評価やそれらへの対策、火災対策や電源の多重化をはじめ、シビアアクシデントに対処する為の機能が挙げられる(図 1.1 参照)。しかしながら、新規制基準で見直しが行われた箇所以外にも事故リスクは存在するため、種々の事故リスクの適切な評価による安全性・信頼性の向上や、保全作業上の安全対策等の弛まぬ改善は必要であり、今後も多くの検討が必要であると言える。

2

図 1.1 従来の規制基準と新規制基準の比較[3]31.2 国内における反応度投入事象
原子炉施設の安全設計の妥当性を確認するため、「発電用軽水型原子炉施設の安全評価に関する審査指針」では、異常状態すなわち「運転時の異常な過渡変化」および「事故」について解析し、評価を行うことを求めている[4]。運転時の異常な過渡変化では「原子炉起動時における制御棒の異常な引き抜き」や「出力運転中の制御棒の異常な引き抜き」、事故においては「制御棒飛び出し」(PWR)、「制御棒落下」(BWR)といった制御棒引き抜きに関する項目について触れられており、原子炉に反応度が投入される事象である反応度事故（RIA :Reactivity Initiated Accident）についての評価が行われている。
一方で、1999 年における志賀原子力発電所 1 号機における臨界事故を始め国内原子力発
電所では 10 件の制御棒引き抜け事象が発生している(図 1.2 参照)。そのうち 8 件が上記審査指針で想定されている起動・運転状態ではなく定期検査中に発生している。また、8 件の定期検査中における制御棒引き抜け事象のうち福島第一原子力発電所 3 号機ならびに志賀原子力発電所 1 号機の 2 件では臨界状態となったという報告がなされている[5][6][7]。これらは制御棒駆動機構（以下 CRD）水圧系の水圧制御ユニット（以下 HCU）隔離作業中に発生しており、その発生メカニズムを、志賀原子力発電所 1 号機の例を元に説明する。
図 1.3、1.4 に HCU ならびに CRD の概略図を示す。CRD は 101 弁ならびに 102 弁により HCUと接続されている。101 弁側より水圧をかける事により制御棒は炉内へと挿入され、102 弁側より水圧をかける事により制御棒を固定しているコレットフェインガが解放され、制御棒の引き抜きが生じる。定期検査中の HCU 隔離作業においては、CRD の弁を順次閉止していくことにより CRD を HCU より隔離していく。しかし、CRD の隔離数が増える事により CRD ポ
ンプの流量が減るため、Q-H 曲線に従い HCU の圧力は上昇していく。その結果、HCU に接続されている CRD では 101 弁ならびに 102 弁に水圧が掛けられた状態となる。この時、101 弁側への加圧により制御棒が過剰に挿入され、制御棒を固定しているコレットフィンガが解放される。そして、102 弁側への加圧によりコレットフィンガが開放状態で固定される。志
賀における事故では、弁閉止手順を間違えた事により挿入側圧力が消失したため制御棒の引き抜け事象が発生した。
志賀の事例では炉内中性子高の警報によりスクラム信号が発信されたため引き抜けは停止したものの、緊急挿入用制御ユニットアキュムレータの充填圧力が無かったためスクラ
ムは行われず、手動にて弁を開放するまで臨界は継続した。この時、中性子束レベルの記録より 15 分間臨界が継続した事が分かっている。
このように、RIA は潜在的な危険として確かに存在しており、そのリスク評価ならびに安全対策は必要であると言える。

4

図 1.2 国内における制御棒引き抜け事象

図 1.3 HCU 概略図5図 1.3 CRD 構造概略図[8]61.3 原子炉における各種フィードバック
原子炉動特性解析について述べる前に原子炉における各種フィードバックについて述べる。通常、原子炉の出力は制御棒の挿入や再循環ラインの流量により制御されるが、出力上昇が速く能動的な制御が追従できないような状態が生じても、適切な設計を行うことで核分裂連鎖反応を自然に制御することが可能となる。これを原子炉の自己制御性と呼ぶ(図1.2 参照)。自己制御性の要因としては、複数の負のフィードバックが存在しており、ドップラー効果や密度効果（ボイド効果）がそれにあたる。
原子炉の燃料温度が上昇した時、燃料の原子は熱運動が激しくなり、ある特定のエネル
ギーを持つ中性子を対象に考えた場合相対速度のばらつきが大きくなる。その時、鋭いピークを持つ共鳴吸収断面積を考えると、その共鳴吸収断面積は実効的に潰れたような形となり広がりを持つようになる。この時、共鳴吸収断面積の幅が広がる事により共鳴の自己遮蔽効果が低下し、中性子の吸収量が増加する。軽水炉の場合、燃料には濃縮度が 2~4%の低濃縮ウランが用いられ、多量の 238U が含まれている。238U は多数の鋭い共鳴ピークを持つ
為、-燃料温度の増加により反応度が低下する。
また、減速材温度が上昇する時、減速材は膨張し密度が低下する。その結果として減速材による中性子の減速能力が低下するため反応度が低下する。BWR の場合、減速材温度の上昇により生成されるボイド量に変化が生じ、減速材の密度変化が発生する。

図 1.4 原子炉の自己制御性[9]71.4 研究目的
現在、原子力発電所では福島第一原子力発電所の事故を経て更なる安全性向上の為に各種対策が行われている。しかしながら、現在の新規制基準では、重大事故対策や耐震・津波対策、そして自然現象・火災に対する考慮など、原子力発電所に対する外部からの影響を遮断もしくは緩和する項目が中心であるため、反応度事故などの他の事故のリスクも適切に評価していく必要がある。本研究の目的は、更なる安全性の向上の為に RIA を対象とした解析の高度化を行うことである。
RIA 解析の高度化にあたり、炉物理モデルの詳細化ならびに熱水力解析との統合による核熱結合計算を行うアプローチを本研究では採用した。炉物理計算においては、中性子・核反応による物理現象の忠実なモデル化を行うことにより、RIA 時における熱の発生とそれに伴う反応度フィードバックをより正確に考慮する事が可能となる。そのため、核分裂生成物による遅発中性子生成のメカニズムを従来のモデルより厳密に取り扱う解析モデルの開発を行った。また、熱水力解析では熱水力解析コード TRAC を用いる事により二流体モデルを適用し、汎用炉物理コードと結合する事によって RIA 時の高精度な過渡応答解析を行い、原子炉過渡応答解析の高度化を図った。

8

第2章

解析手法

2.1 CBZ 概要
本研究では、汎用炉物理解析コードシステム CBZ を用いて核計算を行った。CBZ は核分裂性物質を含む系における核分裂連鎖反応や原子炉燃料の燃焼（崩壊をはじめ、中性子核反応による原子核の転換等）、放射線の遮蔽といった、原子炉内とその周囲における中性子・ガンマ線の輸送、それらの原子核・原子との相互作用に関わる物理現象を数値的に模擬するためのコードシステムである。
CBZ は、日本原子力研究開発機構で開発されたコードシステム CBG をベースとしており、2012 年より本研究室で開発が行われている。コンピュータ言語 C++で記述されており、炉物理の計算に用いるための情報（炉心構造や燃料集合体を表現するための幾何形状や媒質、計算手法や収束条件等）は全て C++の「クラス」として定義されている。また、中性子やガンマ線の挙動を記述する輸送方程式、拡散方程式を解くソルバーが複数実装されており、それらソルバーはそれぞれ関連する複数のクラスで構成されている。そして、ソルバー間のデータ（中性子反応断面積データ等）のやり取りはクラスのインスタンスを通して行われるため、複数のソルバーを組み合わせた多種多様な炉物理計算が容易に実現可能であることが特徴である。
CBZ は決定論的手法に基づいたコードシステムであり、エネルギーは離散化して取り扱っている。そのため、共鳴自己遮蔽効果を考慮して実効群断面積を計算する必要がある。その際に用いる多群ライブラリとして、CBZLIB が実装されている。CBZLIB は Bondarenko 型の炉定数であり、無限希釈断面積、非等方散乱の高次のルジャンドル係数を含む散乱マト
リクスに加えて、共鳴自己遮蔽因子を有している。ライブラリは NJOY-99 コード等を用いて評価済み核データファイルを処理することにより作成されており、任意のエネルギー群構造を使用することが可能である。CBZ を用いた一般的な解析では、軽水炉を対象とした場
合には、SRAC コードシステムと同様の群構造であるエネルギー107 群のライブラリを使用する。
動特性計算は、以前開発した衝突確率法による動特性解析手法を、新たに CBZ の機能として追加して用いており[10]、以下の式(2.1)のように定式化されている。なお、式(2.1)
ではエネルギー群の記載は省いている。

9

1

←
∈

Σtj:領域 j における全反応断面積[1/cm]
Σsj:領域 j における散乱断面積[1/cm]
Vj:領域 j における体積[cm3]φj:領域 j における中性子束[1/cm2/s]
si:領域 i における中性子生成量[1/s/cm3]

v:中性子速度[cm/s]
P(j←i):領域 i で発生した中性子が領域 j で初めて反応する確率[-]

10

2.1

2.2 EFP モデル
原子炉の動特性解析を行う上で、遅発中性子は重要な役割を持つ。遅発中性子は、主として核分裂生成物（以下 FP）の崩壊により生成される中性子であり、遅発中性子を生成す
る FP を遅発中性子先行核と呼ぶ。遅発中性子先行核としては、87Br, 88Br, 137I, 138I, 139I など多数の核種が存在し、異なる半減期で崩壊し遅発中性子を生成する。その為、異なる半減期により遅発中性子が生成される状態を再現するために、通常原子炉動特性解析では遅発中性子先行核を半減期ごとに 6 個のグループに分類し取り扱う（表 2.1 に一例を示す）。
核分裂により生成される全中性子に対して、遅発中性子の占める割合は 1%未満(表 2.1 に示す例では 0.65%)と非常に少ない。しかし、原子炉の過渡応答を緩やかにする為、原子炉の制御という面では極めて重要であると言える。

表 2.1 235U 熱中性子核分裂による遅発中性子データ(JENDL-4.0 より)
グループ

半減期[s]

崩壊定数λi[s-1]先行核生成割合 ai

1

55.72

0.01240.033072

22.72

0.03050.21901

3

6.22

0.111

0.19594

4

2.30

0.301

0.39496

5

0.610

1.14

0.11504

6

0.230

3.01

0.04199

上で述べた 6 群モデルは極めて実用的であり、原子力発電所の設計・解析に耐え得る方法として広く認知されている。しかし、この 6 群モデルは極めて有用であるものの、多くの物理現象を簡略化してモデル化されている。一方で、FP には 1000 種類以上の核種が存在し、遅発中性子先行核の種類も 100 種以上が存在している。また、FP の崩壊により他の FPが生成される事もある。これらのメカニズムを厳密に取り扱うために、全ての遅発中性子
先行核を含む、全ての核分裂生成物を陽的に取り扱う原子炉過渡解析モデルとして EFP モ
デル（Explicit representation of Fission Products model）の開発を行った。EFP モデ
ルを用いる事で、6 群モデルでは簡略化のために無視せざるを得なかった種々の物理現象を再現することが可能となる。その一例として、FP の崩壊による他の FP の生成が挙げられる。
この効果を、これ以降チェーン効果と呼称する。また、FP の崩壊を厳密に取り扱うことにより、核分裂により放出される即発γ線に加え、FP の崩壊により放出される遅発γ線をより厳密に考慮することが可能となる。それにより、ゆくゆくはγ線・中性子線により減速

11

材が直接加熱される効果である PMH 効果(Prompt Moderator Heating Effect)の高精度な評価を行う事が出来る。
図 2.2 に従来の 6 群モデルと EFP モデルの概略図を示す。従来の 6 群モデルでは先行核生成割合に従い遅発中性子先行核グループの生成を行っている。その後、各グループの半減期に従い崩壊し遅発中性子を生成する。一方 EFP モデルでは、核分裂収率を用いる事でFP の生成量を決定している。ここで、チェーン効果を考慮する場合は独立核分裂収率を用いる事で 1 回の核分裂で生成される FP の生成量を計算し、チェーン効果を考慮しない場合
は累積核分裂収率を用いる事により 1 回の核分裂により生成される最終的な FP の生成量を計算している。各 FP は収率に従い生成された後に、それぞれの半減期に従い崩壊する。こ
こで、チェーン効果を考慮する場合は FP の崩壊時の分岐比に従い他の核種の生成を行う事
で崩壊による FP の生成を行っている。

図 2.2 EFP モデル概略図

12

2.3 TRAC 概要
本研究では解析コード TRAC を用いる事により、熱伝導・熱伝達計算並びに減速材領域におけるボイド率の計算をおこなっている。TRAC は二流体モデルを用いた二相流解析コード
であり、その特徴を以下に示す。
・圧力、温度、ボイド率をセル中心、流速をセル境界で評価する。
・流動解析の基礎方程式ならびに構成方程式は差分法で解く。
・解析手法は J.Mahaffy により提唱された SETS 法(stability-enhancing two-step method)を用いて解いている。SETS 法では半陰解法の前後に安定化ステップを設けており、質
量・エネルギー・運動量の伝播を陰的に取り扱っている。それにより、CFL 条件（Courant-Friedrichs-Lewy Condition）の制約を緩め、かつ方程式における各項の取り扱いを最小化している。
・各構成方程式はモジュール化されており、プログラムの変更・保守管理が容易に行えるようになっている。
・水・水蒸気系あるいは水・空気系に対する物性値計算が標準で組み込まれているが、流体の物性値計算箇所もモジュール化されている為物性値を変更することにより他の流体の解析を行うことも出来る。
ここで、TRAC コードにおける管内二相流解析モジュールである PIPE コンポーネントについて触れる。PIPE コンポーネントは一次元モデルにより解析を行っており、このコンポー
ネントでは気相・液相に対して連続の式、運動方程式、エネルギー式を用いる事により解析を行っている。それぞれの式を、次の微分方程式(2.2)~(2.7)に示す。また、この時の添字 g,l はそれぞれ気相と液相である事を示す。

気相の連続の式
2.2
液相連続の式

11

2.3

気相の運動方程式

cos

13

2.4

液相の運動方程式
1

1

| |

1
cos

2.5

気相のエネルギー式
2.6
二相のエネルギー式
1

1

1

2.7

α:ボイド率[-]ρg(l):気相（液相）密度[kg/m3]
ug(l):気相（液相）流速[m/s]Γ:相変化に伴う蒸気生成速度[kg/m3/s]
P:圧力[Pa]g:重力加速度[m/s2]
Ci:気液界面におけるせん断応力係数[kg/m4]
Cwg(wl):気相(液相)の壁面せん断応力係数[kg/m4]
qig(il):気相から液相へ（液相から気相へ）流入する熱量[kg/m/s3]
qwg(wl):壁面から気相(液相)へ流入する熱量[kg/m/s3]
hlg:蒸発潜熱[m2/s2]
θ:流路の水平面に対する角度[°]以上 6 本の基礎方程式に加え、次の(2.8)~(2.11)に示す気液密度と内部エネルギーについての状態方程式が用いられている。
,

2.8

,

2.9

,

2.10

,

2.11

14

また、TRAC コードでは、熱伝導体における熱伝導、熱伝導体と流体間における熱伝達計
算を HEAT CONDUCTOR コンポーネントで行っている。HEAT CONDUCTOR コンポーネントでは二次元熱伝導モデルにより解析を行っており以下に特徴を示す。
・軸方向、径方向の二次元熱伝導を扱い、オプションにより軸方向熱伝達がない場合の計算も行う事が可能である。
・円筒形・矩形形状の熱伝導体の計算を行うことが出来る。
・ファインメッシュ計算は行わない。
・ギャップ計算は省略する。
・ジルカロイ-水反応による効果は省略する。
・熱伝導体の上下端においては断熱であると仮定する。
円筒体系において HEAT CONDUCTOR コンポーネントで用いられている 2 次元非定常熱伝導
方程式の基礎式を(2.12)に示す。
1′′

Cp:定圧比熱[m2/s2/K]
r:半径[m]ρ:密度[kg/m3]q’’:体積発熱[kg/s2/m]
k(K):径方向（軸方向）に対する熱伝達率[kg/s3/K]

15

2.12

2.4 結合計算概要
本研究では CBZ コードと TRAC コードを統合することにより、核熱結合解析を行った。
図 2.3 に、解析フローを示すとともに、以下に統合時に行った変更を記載する。
・TRAC コードでは二相流の状態によって解析の時間幅を調整することにより、解析の高速化が図られている。しかし、結合計算を行う上では CBZ 側の時間幅と一致させる必要があるため当該機能を停止した。
・TRAC コードに新規制御モジュールを作成し、CBZ コード側からの発熱量受け渡し並びに TRAC コードからの減速材・燃料ピン温度受け渡し機能を実装した。
・CBZ コードに、温度に対応して減速材密度ならびに微視的核反応断面積を調整する機能の追加をした。
・入力データは CBZ コードならびに TRAC コードで個別に用意し、熱水力解析や核計算に用いるパラメータはそれぞれのモジュールにより取り込みを行った。

16

図 2.3 結合計算時解析フロー

17

第3章

燃料ピンセル過渡応答解析

3.1 解析条件
解析体系は PWR の燃料ピンを想定したピンセル体系として、解析体系は燃料と減速材からなる簡易的な体系とし、境界条件を等方反射条件とした。概略図と諸元を図 3.1 ならびに表 3.1 に示す。フィードバックを考慮しない解析において、断面積ライブラリ作成に用いた温度も併せて示す。燃料組成は濃縮度 3.4%の低濃縮ウランを想定した。核反応断面積の計算に用いた原子の数密度を表 3.2,3.3 に示す。ここで、表 3.2,3.3 で示す数密度は燃
料温度 800K、減速材温度 600K のものを示している事に注意が必要である。

図 3.1 解析体系概略図

表 3.1 解析体系諸元

燃料領域

減速材

メッシュ数

2

5

半径(cm)

0.412

0.7137

温度(K)

800

600

表 3.2 燃料領域原子数密度(1/barn/cm)
235

U

7.753E-4

238

U

2.175E-2

O

4.505E-2

1617

表 3.3 減速材領域原子数密度(1/barn/cm)
1

H

16

5.572E-2O

10

2.786E-2

B

4.592E-6

Ni

3.688E-4

Cr

1.609E-4

Fe

1.306E-4

中性子のエネルギーは核計算コード SRAC と同一のエネルギー構造である 107 群構造を使用した。遅発中性子先行核 6 群モデルにおける各種定数や分岐比ならびに遅発中性子ス
ペクトル、核反応断面積の計算には JENDL-4.0 の核データを使用した。EFP モデルで必要とされる核分裂収率データには JENDL FP Fission Yields Data File-2011 を使用し、FP の分
岐比・半減期といった各種崩壊データは JENDL FP Decay Data File-2011 のものを使用した。そして、遅発中性子先行核より放出される中性子のエネルギースペクトルに関しては
ENDF/B- .1 のデータを使用した。また、断面積データの 107 群への縮約には CBZ コードの自己遮蔽計算クラスを利用した。実効増倍率の調整や反応度の印加は体系のバックリン
グを調整することにより行っている。遅発中性子先行核 6 群モデルと EFP モデル間における遅発中性子割合βに差が生じないよう、EFP モデルにおける 1 回の核分裂により生成される遅発中性子量を計算し、EFP モデルにおけるβと 6 群モデルにおけるβが一致するよ
うに EFP モデルの収率に係数をかける事で、FP 全体の生成量を調整した。また、今回の体系では、燃料領域には 235U ならびに 238U が存在しており、238U の高速核分裂による FPの生成を考慮する必要もある。その為、燃料領域では巨視的反応断面積により核分裂反応を計算せず、
微視的反応断面積と数密度により 235U と 238U の反応率を計算している。また、今回の解析においては FP の中性子による核反応は考慮していない。
次に、熱水力解析の諸条件について述べる。解析体系概略図並びに PIPE コンポーネントの構成を図 3.2,3.3 に、諸元を表 3.4 に示す。境界条件は圧力で規定し、低圧側・高圧側を共に 15.5MPa とした。燃料ペレットの物性値は、今回は簡便のため温度依存性を考慮せず
表 3.5 の値を用いて計算を行った。

18

表 3.4 熱水力解析体系諸元
流路面積(m2)

1.47E-4

等価水力直径(m)1.47E-4

水温(K)

574.2

圧力(Pa)

15.5E+6

表 3.5 燃料ペレット物性値
密度(kg/m3)

10.42

比熱(J/K/kg)

300

熱伝導率(W/m/K)

4.6

図 3.2 熱水力解析体系概略図19

図 3.3 PIPE コンポーネントにおける解析体系203.2 EFP モデルにおける解析
この節では、従来の 6 群モデルと EFP モデル過渡応答解析を行うことにより、EFP モデルによる影響の確認を行った。解析体系は前節で述べたものを使用する。また、この節の解析では、従来の遅発中性子先行核 6 群モデルでは 1 ケース、EFP モデルではチェーン効果の有無により 2 ケース、合計 3 ケースの解析を行った。初期増倍率ならびに投入反応度は各解析で述べる。また、この節の解析では EFP モデルによる効果に着目するために熱的・核的フィードバックは考慮していない。

21

3.2.1 即発臨界未満の正の反応度を印加
本解析では、臨界状態の体系に即発臨界未満の正の反応度を印加した時の過渡応答解析を行った。初期状態を臨界として計算した後に、0.1 秒において 60¢の反応度を印加した。
出力変動を図 3.4 に示す。

Output(arb. unit)

6

6group
EFPmodelEFPmodel(chain effect)

5
4
3
2
1

0.0

0.4

0.8

1.2

1.6

2.0

Time(sec)
図 3.4 臨界体系に 60¢の反応度を印加した場合の出力推移

この解析結果では、EFP モデルと既存の 6 群モデルとでは出力推移に差が生じている事が分かる。この理由としては、EFP モデルと 6 群モデル間における遅発中性子の時間依存の放出率が関係していると考えられる。図 3.5 に熱中性子照射時における U-235 の遅発中性子放出率のグラフを示す。このグラフは、一回の核分裂で生成される遅発中性子の放出率の時間推移を示しており、JENDL-4.0 における 6 群モデルにおける放出率ならびに全てのFP を取り扱った総和計算による放出率を図示している。なお、このグラフでは全遅発中性子の生成量の規格化がなされていない為、6 群モデルと比較し総和計算では遅発中性子の放出量が 1 割程度多い事に留意が必要である。それを考慮しても、このグラフから 6 群モデルと比較して総和計算ではより早く放出率のピークが生じている事が分かる。これは 6 群モデルと比較し総和計算ではより早いタイミングで遅発中性子が生成される事を表しており、即ち 6 群モデルと比較して EFP モデルでは中性子の増倍が早くなるということを意味
している。

22

図 3.5 235U バーストフィッション時の遅発中性子放出率[11]この原因として考えられるのは、評価済み核データが挙げられる。EFP モデルでは全てのFP を取り扱っており、取り扱う核データとしては核分裂収率や、どのように崩壊するかの
崩壊様式、分岐比、崩壊定数、崩壊時における中性子スペクトルといった崩壊データが挙げられる。しかし、取り扱っている FP の中には物理定数の測定が困難な短半減期核種も存在しており、その多くはシステマティックスあるいは核模型理論により推定されている。
そのため、真値と比較して大きく違う事も考えられる。このような影響により、6 群モデルと EFP モデル間に差が生じたと考えられる。6 群モデルは直接測定値を元に作成されており、精度も良い。その為、EFP モデルで 6 群モデルの結果を再現する事が今後の課題としては挙げられる。
また、EFP モデル間でもチェーン効果の有無により増倍に差が生じている事が図 3.4 から分かる。この原因としては、チェーン効果を考慮することにより遅発中性子先行核の実効的な半減期が長くなっているという事が考えられる。この効果の概略図を図 3.6 に示す。
EFP モデルでは、チェーン効果を考慮する場合には独立核分裂収率を、考慮しない場合には累積核分裂収率を使用する事を 2.2 節で述べた。独立核分裂収率は、1 回の核分裂で生成される核種の生成率を表している。その為、独立核分裂の総和は 2 となる。一方、累積核

23図 3.6 チェーン効果による実効的な半減期への寄与

分裂収率では 1 回の核分裂で最終的に生成される核種の生成率を表している。
このように、累積核分裂収率は崩壊プロセスを経て最終的に生成される FP の量を表しているため、FP の崩壊により遅発中性子が生成される事を考慮しない分、チェーン効果を考慮した場合と比較して遅発中性子先行核の生成が早くなると考えられる（図 3.6 における核種 B が相当する）。その為、遅発中性子先行核から放出される遅発中性子に着目した時、チェーン効果を考慮した場合は実効的な半減期が延び遅発中性子の放出が遅くなると言える。
この効果により、チェーン効果の有無により出力推移に差が生じたと考えられる。

24

3.2.2 負の反応度を印加
本解析では、臨界状態の体系に負の反応度を印加した時の過渡応答解析を行った。初期状態を臨界として計算した後に、0.1 秒において-25¢の反応度を印加した。出力変動を図 3.7に示す。
今回は、印加反応度が小さい事による影響か、チェーン効果の有無による差は非常に小さなものとなった。3.2.1 節では、6 群モデルと比較して EFP モデルでは遅発中性子の生成が早い為、出力の増倍が早くなる事を述べた。今回の計算結果では 6 群モデルと比較して、EFP モデルにおける出力減衰の方が早い結果となった。この要因としては、図 3.5 で示した遅発中性子生成率が挙げられる。EFP モデルでは、遅発中性子生成率がより早くピークに到達する。これは、遅発中性子先行核の崩壊が 6 群モデルと比較して EFP モデルでは早い事を意味している。その為、未臨界状態の体系では遅発中性子先行核の減衰が早くなり 6群モデルと比較してより早く減衰する事が分かる。
チェーン効果の有無による効果も同様であり、実効的な半減期が短いチェーン効果を考慮しないモデルではチェーン効果を考慮した場合と比較して遅発中性子先行核の減衰が早い為、より早く減衰していると考えられる。

Output(arb. unit)

1.0

6group
EFPmodelEFPmodel(chain effect)

0.8
0.6
0.4
0.2

0.0

0.4

0.8

1.2

1.6

Time(sec)
図 3.7 臨界体系に-25¢の反応度を印加した場合の出力推移

25

2.0

また、図 3.8 に本解析における 0.1 秒以降の出力推移を対数グラフで表したグラフを示す。
当初、本解析では臨界体系に負の反応度を印加することにより、パルス中性子照射時の挙動を模擬することを目的としていた。そのため、反応度印加時の即発跳躍では単一の減衰定数による指数関数的な減衰が生じると考えていた。しかし、結果は図 3.8 に示すようになり、対数グラフ上で確認すると直線的ではなく連続的に変化する曲線状の挙動となった。
すなわち、減衰定数は連続的に変化しており、当初の予想とは異なる結果となった。
この原因としては、遅発中性子先行核による影響が考えられる。通常、パルス中性子照
射体系は未臨界状態であり、遅発中性子先行核を始めとした FP は存在していない、あるいは非常に少ない状態である。その為、パルス中性子照射後の核分裂により FP が生成され遅発中性子が発生するまで、体系は即発中性子のみで増倍が行われる事となり特定の減衰定数を持つこととなる。一方、臨界体系に対して負の反応度を印加した場合、体系には既に十分な量の遅発中性子先行核が存在している。その為、反応度を印加した直後の即発跳躍には即発中性子による増倍に加え、既に存在している遅発中性子先行核からの遅発中性子による影響も出ている事が考えられる。その結果として、図 3.8 に示す出力推移となった事が推察できる。

1

Output(arb. unit)

6group
EFPmodelEFPmodel(chain effect)

0.1

0.101

0.102

0.103

Time(sec)
図 3.8 反応度印加に伴う即発跳躍時の出力推移

26

0.104

3.2.3 未臨界体系への中性子パルス照射
本解析では、パルス中性子法による実験を模擬した解析として、未臨界状態の体系に負の反応度を印加した時の過渡応答解析を行った。初期状態を実効増倍率 0.9900 の未臨界体系として 0.1 秒において燃料領域中心に中性子パルスを照射した。今回は、パルス照射により生成された FP による出力変動に着目するため、0.1 秒以降の出力変動に着目した。この時の出力変動のグラフを図 3.9 に、6 群モデルと EFP モデル間の差異に着目するため詳細部に着目したグラフを図 3.10 に示す。パルス中性子法では、中性子束が指数関数的に減衰することが広く知られており、その減衰定数が重要な役割を持つ。その為、本解析における出力変動は対数グラフで示す事に留意が必要である。
図 3.9,3.10 では出力推移の傾きより、6 群モデルが 0.107 秒段階で減衰定数に変化が生じている事と比較して EFP モデルではより早い 0.105 秒段階で減衰定数が変わっていることが確認できる。これは、EFP モデルでは遅発中性子先行核を陽的に取り扱っているため超短半減期核種も考慮している影響が考えられる。また、6 群モデルと比較して EFP モデルの最終的な出力レベルが高くなっている要因としては、3.2.1 でも触れているように遅発中性子放出率による影響が考えられる。EFP モデルでは遅発中性子に放出が早い為、6 群モデルと比較して早い段階で遅発中性子による影響が生じ始める。その結果として、今回の出力推移に繋がっている事と考えられる。しかし一方で、理論上では最終的な出力レベルは遅発中性子割合βならびに印加した反応度に依存する。そのため、遅発中性子割合を同一となるように調整している今回の解析でこのような差異が生じた原因は不明な点も多いと言えるためより多くの検討が必要であると考えられる。

27

Output(arb. unit)

10-2

6group
EFPmodelEFPmodel(chain effect)

10-3
10-4
10

-5

10

-6

10

-7

0.100

0.104

0.108

0.112

0.116

0.120

Time(sec)
図 3.9 未臨界体系におけるパルス中性子照射時の出力推移

Output(arb. unit)

10-4

10

-5

10

-6

0.104

6group
EFPmodelEFPmodel(chain effect)

0.106

0.108

0.110

Time(sec)
図 3.10 未臨界体系におけるパルス中性子照射時の出力推移（詳細）283.3 熱水力結合解析における解析
この節では、TRAC コードと CBZ コードの結合解析を行う事により、核的フィードバッ
クと熱水力計算による減速材温度フィードバックによる効果を確認した。解析体系は 3.1
で述べたものを使用する。また、燃料棒の出力は一般的な PWR 燃料棒の線出力密度を参考に 20kW/m として解析を行った。そして、今回の解析では時間幅を 1.0E-04 秒として計算を行った。本解析では、TRAC コードによる熱伝導解析モジュールにより燃料温度を計算
し、CBZ コードにおいて燃料領域の核反応断面積データを変更する事により核的フィード
バックを考慮した。なお、燃料領域の核反応断面積は燃料温度に体積により重み付けをした平均温度を用いて決定した。また、減速材温度フィードバックは TRAC コードにおける熱伝達計算の結果により減速材領域における核反応断面積ならびに減速材密度を変更する
ことにより考慮した。ここで、減速材密度は基準状態からの比をとっている。この解析では、初期条件として燃料ペレット内温度を 800K、減速材温度を 575K として断面積データを設定しており、減速材密度変化を考慮する上での基準状態も 575K とした。解析は、核的フィードバックを考慮したケースならびに核的フィードバックに加え、減速材温度フィー
ドバックを考慮したケースの合計 2 ケースの解析を行った。この節の解析では熱的・核的フィードバックによる効果に着目するために EFP モデルは使用せず 6 群モデルでの解析を
行った。

29

3.3.1 即発臨界以上の正の反応度印加
本解析では、臨界状態の体系に即発臨界以上の正の反応度を印加した時の過渡応答解析を行った。初期状態を臨界として計算した後、同一のバックリング調整を行うことにより 0秒段階で約 2＄の反応度を印加した。出力変動を図 3.11 に示し、-0.02 秒から 0.1 秒までの出力変動を図 3.12 に示す。また、核反応断面積の決定に用いた燃料平均温度を図 3.13 に示しその詳細を図 3.14 に、減速材温度を図 3.15 に示し、減速材のボイド率推移を図 3.16 に示す。
図 3.11 より、減速材フィードバックによる効果が反応度投入以降の出力変動に大きな影響をもたらしている事が確認できる。また、図 3.12 に示しているように、ドップラーフィ
ードバックにより極めて早く出力が抑制されている事が分かる。また、図 3.12 や図 3.13、
図 3.15 より燃料温度の減少に伴いドップラーフィードバックが減少するも、減速材温度の上昇による減速材フィードバックにより十分に補償されている事が分かる。
また、図 3.15 では減速材フィードバックを考慮しない解析において 0.8 秒近辺より減速材温度がほぼ一定に到達している。この理由としては、減速材が飽和温度に到達した事が挙げられる。15MPa における飽和温度は 615K であり、図 3.16 に示すように 0.8 秒近辺よりボイド率が上昇している事からも減速材が飽和温度に到達した事が伺える。
この解析では、図 3.12 に示されているように減速材温度フィードバック考慮の有無により最大出力に違いが生じている。この原因としては解析条件による減速材温度の固定が挙げられる。減速材温度フィードバックを考慮しない解析では、減速材領域における断面積
データの更新ならびに密度変化の考慮時に用いる減速材温度を 575K と固定して解析を行っている。一方、考慮を行う計算では熱水力解析により算出されている温度を使用し解析を行っている。その為、2 つの体系は異なる増倍率を持つこととなる。結果として、今回の解析のように同一の外乱を与えた計算では印加される反応度に差異が生じ、最大出力の差が生じたものと思われる。

30

4

without coolant feedback
with coolant feedback

Output(arb. unit)

3.5
3
2.5
2
1.5
1
0.5

0

0.2

0.4

0.6

0.8

Time(sec)
図 3.11 臨界体系に 2.5＄の反応度印加時の過渡応答

4

Output(arb. unit)

3.5

without coolant feedback
with coolant feedback

3
2.5
2
1.5
1
0.5
-0.02

0

0.02

0.04

0.06

0.08

Time(sec)
図 3.12 臨界体系に 2.5＄の反応度印加時の過渡応答（詳細）

31

0.1

1600

without coolant feedback
with coolant feedback

Temperature(K)

1500
1400
1300
1200
1100
1000
900
800
700

0

0.2

0.4

0.6

0.8

Time(sec)
図 3.13 燃料温度時間推移

1500

Temperature(K)

1400
1300
1200
1100
1000
900
without coolant feedback
with coolant feedback

800
700

0

0.02

0.04

0.06

Time(sec)
図 3.14 燃料温度時間推移（詳細）

32

0.08

Temperature(K)

625
620
615
610
605
600
595
590
585
580
575
570

without coolant feedback
with coolant feedback
0

0.2

0.4

0.6

0.8

Time(sec)
図 3.15 減速材温度推移

Vapor fraction

0.4

without coolant feedback
with coolant feedback

0.3
0.2
0.1
0

0

0.2

0.4
Time(sec)

図 3.16 ボイド率推移

33

0.6

0.8

第4章

結論

本研究では遅発中性子先行核を陽的に取り扱うモデルの開発を行い、核分裂生成物による遅発中性子生成のメカニズムをより厳密に取り扱う解析モデルの開発を行った。そして、
衝突確率法による動特性解析手法を用いる事により、臨界体系への反応度印加、負の反応度を印加、未臨界体系へのパルス中性子照射の解析を行った。その結果、EFP モデルでは遅発中性子放出率が既存の 6 群モデルと比較し早くピークを迎える事により、出力の増倍が促進される方向に誤差が生じる事を確認した。また、核分裂生成物の崩壊による他核種の生成をモデル化することにより、チェーン効果を考慮することで遅発中性子の放出に時間遅れが生じるという知見を得た。
そして、汎用炉物理解析コード CBZ と二相流解析コード TRAC との結合計算を行うことにより核的フィードバックならびに減速材温度フィードバックを考慮した解析を行い、RIA 時
における過渡応答解析を核的ならびに熱水力的側面より検討する事が可能となった。
今後の課題としては、EFP モデルの改良ならびに調整が挙げられる。EFP モデルを用いる
事により、FP の崩壊を厳密に取り扱う事が可能となった。その為、現在の EFP モデルを拡張し、中性子輸送計算に加えγ線輸送計算を実装することにより、減速材直接加熱（PMH）
効果を厳密に取り扱うことが可能となる。そして、γ線や中性子による被覆管・減速材の加熱効果を体積発熱として熱水力解析に使用することで、より物理現象に即した解析を行うことが出来るようになると考えられる。
また、FP を厳密に取り扱う事が可能となったため FP の中性子捕獲反応といった FP による核反応を考慮する事も課題として挙げられる。通常、遅発中性子先行核の核反応断面積
は崩壊定数と比較して非常に小さい。その為、通常の解析においては無視出来るとかんがえられる。しかし、RIA といった中性子束レベルが桁オーダーで上昇するばあいにはその効果を確認することが予測される。また、加速器駆動未臨界炉(ADS:accelerator-driven
subcritical reactor)といった FP の燃焼を取り扱う系における過渡応答解析にも効果を発揮するものと考えられる。
他にも、現在の EFP モデルは従来の 6 群モデルと比較して計算結果に大きな差が生じている事が分かる。これは、遅発中性子先行核の多くは短半減期核種であり、収率や各構造が実験的にはあまり知られていない FP 核種である事が挙げられる。また、遅発中性子先行
核として寄与する核種数が比較的少ない為、それぞれの核種のデータの精度が遅発中性子の総和計算に影響を及ぼしているとも考えられる。その為、6 群モデルの定数、あるいはEFP モデルに用いられている炉定数を調整することによりモデル間の差異を少なくし、モデ
ル間のメカニズムによる挙動の違いを見る事も課題として挙げられる。

34

一方、CBZ コードの過渡応答解析ならびに TRAC コードの熱水力解析自体は検証がなされているものの、結合解析としての検証は未だ不十分である。その為、実験と解析を行うことにより検証を行う必要があると考えられる。そして、今回の熱水力解析に用いたモデル
は簡略化の為に物性値の温度変化を省いた点や、セルの分割が不十分、そして発熱分布を一定としているなど改善の余地が多数挙げられる。今後は、解析体系の詳細化を行うこと
により、より良い解析結果が出せると考えられる。

35

謝辞
本研究を進めるにあたり、指導教官である奈良林直教授ならびに、千葉豪准教授に深く感謝いたします。特に、CBZ の開発者である千葉准教授には数多くのご指導を頂き心より感謝いたします。また、研究に有意義な指摘をしていただいた同期の梶原氏にも感謝を、
そして今後のご活躍をお祈りさせて頂きます。また、大学での生活を支えてくださった親に最大の感謝を。

36

参考文献[1]World Nuclear Association , Information Library , Nuclear share figures
http://www.world-nuclear.org/info/Facts-and-Figures/Nuclear-generation-by-country/（2016/1/18 アクセス）[2]国会事故調査委員 報告書
http://warp.da.ndl.go.jp/info:ndljp/pid/3856371/naiic.go.jp/（2016/1/18 アクセス）[3]原子力規制委員会 実用発電用原子炉及び核燃料施設等に係る新規制基準について
http://www.nsr.go.jp/data/000070101.pdf（2016/1/18 アクセス）[4]文部科学省 告示・通達

発電用軽水型原子炉施設の安全評価に関する審査指針についてhttp://www.mext.go.jp/b_menu/hakusho/nc/t19900830003/t19900830003.html（2016/1/18 アクセス）[5] 原子力保安院,(2007),北陸電力株式会社志賀原子力発電所１号機
における平成 11 年の臨界事故及びその他の原子炉停止中の想定外の制御棒の引き抜け事象に関する調査報告書（概要）[6] 原子力保安院,(2007),北陸電力株式会社志賀原子力発電所１号機
における平成 11 年の臨界事故及びその他の原子炉停止中の想定外の制御棒の引き抜け事象に関する調査報告書（概要）[7] 日本原子力学会,(2007),制御棒引き抜け事象調査委員会 (3)臨界事象 WG 報告[8] 失敗知識データベース・志賀原発、臨界事故発覚
http://www.sozogaku.com/fkd/cf/CZ0200701.html（2016/1/30 アクセス）[9] 電気事業連合会,(2015),原子力・エネルギー図面集 2015[10]松浦健太,(2013), 衝突確率法による燃料ピン内温度分布を考慮した反応度フィード
バック動特性解析[11] Go Chiba, Tadashi Narabayashi, (2015) ,Uncertainty quantification of total
delayed neutron yields and time-dependent delayed neutron emission rates in frameof summation calculations Ann. Nucl. Energy 85, pp846-855[12] Y.Kawamoto, G.Chiba, M.Tsuji, T.Narabayashi, (2015), Numerical solution of
matrix exponential in burn-up equation using mini-max polynomial approximation ,Ann. Nucl. Energy , 80, pp219-224

37

[13]日本機会学会,(1999),蒸気表[14] 山根 祐一, 酒井 幹夫, 阿部 仁, 山本 俊弘, 奥野 浩,三好 慶典, (2006), MOX 粉体系 の 過 渡 臨 界 事 象 評 価 の た め の 物 性 値 デ ー タ , JAEA-Data/Code 2006-021; October
2006,75p.
[15]小林 啓祐,(1996), 原子炉物理 初版, コロナ社
[16] 岡嶋 成晃, 久語 輝彦, 森 貴正, (2012), 原子力教科書 原子炉物理学第 1 版, オー
ム社[17] 岡 芳明, 鈴木 勝男, (2008), 原子力教科書 原子炉動特性とプラント制御 第 1 版,
オーム社
[18]原子力安全研究協会実務テキスト編集委員会, (2008), 軽水炉発電のあらまし 第 3 版,
原子力安全研究協会[19]日本原子力学会炉物理部会, (2014), 第 46 回炉物理夏期セミナーテキスト 炉物理と原子力安全-未臨界と動特性の基礎と応用-, 日本原子力学会

38

付録
TRAC-CBZ 結合クラスソースコード-CBZ 側（使用言語：C++）
trac_c.h#ifndef TRAC
#define TRAC

#include <stdio.h>#include <iostream>#include <fstream>#include <vector>#include <cstdio>#include <cmath>#include <stdlib.h>#include <math.h>using namespace std;
extern "C" {extern void mini_();}
extern "C" {extern void cbz_();}
extern "C" {extern void cbza_();}
extern "C" {extern void cinput_();}
extern "C" {extern void cbztrans_();}
extern "C" {extern void cbzheatinput_( double *, double *, int *);}
extern "C" {extern void cbztemp_( double *, int *, int *,double *);};
using namespace std;
typedef double real;
typedef vector< vector<real> > dim2d;
typedef vector< vector< vector<real> > > dim3d;
typedef vector< real > dim1d;
typedef vector< vector<int> > dim2di;typedef vector< vector< vector<int> > > dim3di;typedef vector<int> dim1di;39

class trac_c{
private:public:void run_trac();void run_trac_main();void setup_for_dyn();void run_trac_onestep();
void put_output(double *output,real *out_dis, int *pincells);
void get_temp(double *temper,int *pincells,int *cell,double *alpha);//temper};#endiftrac_c.cxx#include <cstdlib>#include "trac_c.h"using namespace std;
void trac_c::run_trac(){mini_();}
void trac_c::run_trac_main(){ cbz_();}
void trac_c::setup_for_dyn(){cbza_();}
void trac_c::run_trac_onestep(){cbztrans_();}
void trac_c::put_output(double *output,double *out_dis, int *pincells){
cbzheatinput_( output, out_dis, pincells);}
void trac_c::get_temp(double *temper,int *pincells,int *cell,double *alpha){
cbztemp_( temper,pincells,cell,alpha);}

40

TRAC-CBZ 結合クラスソースコード-TRAC 側（使用言語：Fortran）
cbz.forsubroutine cbz
c this is made to use miniTRAC in CBZimplicit real*8 ( a-h , o-z )call mini
end

mini.forsubroutine minic subroutine for using TRAC by subroutine .
implicit real*8 ( a-h , o-z )include 'setprm.inc'include 'bdarry.inc'include 'brktbl.inc'include 'contrl.inc'include 'ctitle.inc'include 'cconv.inc'include 'cvoid.inc'include 'filtbl.inc'include 'fltab.inc'include 'iounit.inc'include 's1dptb.inc'include 's1dptr.inc'include 'tsatcn.inc'common / cversn /character*8 sttime , entime
call time ( sttime )varerm = 0.0d0
varer = 0.0d0
dtlu= 0.0d0

dtll= 0.0d0dtvu= 0.0d041

dtvl= 0.0d0darv= 0.0d0darl= 0.0d0tlmx= 0.0d0tvmx= 0.0d0pmx= 0.0d0dau= 0.0d0dal= 0.0d0oau= 0.0d0oal= 0.0d0

delt= 1.0d-7

odelt = deltcall input
call init
call transcall time ( entime )write(iout,1000) sttime , entime
write(*,1000) sttime , entime
stop
1000 format(1h0,'**********************'/*

1h ,'* problem statistics *'/*1h ,'**********************'/*

1h0,10x,'start time : ', a8 /*1h ,10x,' end time : ', a8 )

end

42

cinput.forsubroutine cinput
implicit real*8 ( a-h,o-z )c input routine for coupling analysis
include 'setprm.inc'include 'iounit.inc'include 'ctitle.inc'include 'contrl.inc'include 'tsatcn.inc'include 'bdarry.inc'include 'fltabs.inc'include 'heatpt.inc'common / cversn /character*60 filnam(3)write(*,'(a)') ' enter file names..'
write(*,'(a)') ' input file name= ./'filnam(1) = "input.in"
write(*,'(a)') ' output file name = cal.out'
filnam(2) = "cal.out"
write(*,'(a)') ' graphic file name = cal.grp'
filnam(3) = "cal.grp"
open ( unit=irdcrd , file = filnam(1) )open ( unit=iout , file = filnam(2) )open ( unit=igrf , file = filnam(3) ,*access='sequential' , form='formatted')

open ( unit=in , status='scratch' ,*access='sequential' , form='formatted')

open ( unit=iwkfl , status='scratch' ,*access='sequential' , form='formatted')call rdcard ( irdcrd , in , iout , 'minitrac' )rewind ( in )read(in,'(a72)') title
write(*,'(a)') 'end title'43

read( in , * ) timet
if ( timet .lt. 0.0d0 ) timet = 0.0d0
write(*,'(a)') 'end timet'
read( in , * ) ieos , ipak
if ( ieos .ne. 1 ) ieos = 0
if ( ipak .ne. 0 ) ipak = 1
write(*,'(a)') 'end ieos,ipak'
read ( in , * ) epso , oitmax
write(*,'(a)') 'end epso oitmax'
read ( in , * ) nfb
if ( nfb .lt. 0 ) nfb = 0
write(*,'(a)') 'end nfb'
read ( in , * ) ihtcnd
if ( ihtcnd.lt.0 .and. ihtcnd.gt.2 ) ihtcnd = 0
write(*,'(a)') 'end ihtcnd'
write( iout , 5100 )write( iout , 6100 ) timet , ieos , ipak , epso , oitmax , nfb ,*ihtcndcall seteos
call rpipe
read( in , * ) islb , isrb
write( iout , 5200 )write( iout , 6200 ) islb , isrb
if ( islb .eq. 0 .or. islb .eq. 2 ) go to 10
call error ( 1 , '*input* left boundary mismatch' , 30 )10 continue
if ( isrb .eq. 0 .or. isrb .eq. 2 ) go to 20
call error ( 1 , '*input* right boundary mismatch' , 31 )20 continue

if ( islb .eq. 0 ) thencall rfill(1)else if ( islb .eq. 2 ) then44

call rbreak(1)end ifif ( isrb .eq. 0 ) thencall rfill(2)else if ( isrb .eq. 2 ) thencall rbreak(2)end ifif( nfb .gt. 0 ) call rfb
ihtprt = 0
if( ihtcnd .gt. 0 ) call rcore
return5100 format(1h0,'********************'/*

1h ,'* main data card *'/*1h ,'********************')

5200 format(1h0,'********************'/*

1h ,'* boundary cards *'/*1h ,'********************')

6100 format(1h0,'start time (timet)= ',1pd15.3/*

1h ,'gas option (ieos) (0 = steam/ 1 = air) =',i15/*

1h ,'water packer (ipak) (0 = off/ 1 = on) =',i15/*

1h ,'convergence criterion (epso)*

1h ,'maximum number of iteration (oitmax)*

1h ,'number of side fill and break (nfb)=',i15/*

1h ,'type of heat conductor (ihtcnd)=',i15/

*

1h ,'

')

= ',1pd15.3/

0 = no/ 1 = cyl/ 2 = slab

6200 format(1h0,'left boundary condition (islb)*

1h ,'right boundary condition (isrb)

*

1h ,'

0 = fill/ 2 = break= ',i5/= ',i5/
')

end

45

=',i15/

cbza.forsubroutine cbza
c initialize routine for coupling analysis
implicit real*8 ( a-h , o-z )include 'setprm.inc'include 'bdarry.inc'include 'brktbl.inc'include 'contrl.inc'include 'ctitle.inc'include 'cconv.inc'include 'cvoid.inc'include 'filtbl.inc'include 'fltab.inc'include 'iounit.inc'include 's1dptb.inc'include 's1dptr.inc'include 'tsatcn.inc'common / cversn /character*8 sttime , entime
varerm = 0.0d0
varer = 0.0d0
dtlu= 0.0d0

dtll= 0.0d0dtvu= 0.0d0dtvl= 0.0d0darv= 0.0d0darl= 0.0d0tlmx= 0.0d0tvmx= 0.0d0pmx= 0.0d0dau= 0.0d0dal= 0.0d0oau= 0.0d046

oal= 0.0d0

delt= 1.0d-7

odelt = deltcall cinput
call init
nstep = 0iofail = 0
nmfail = 0
oitno = 0
iendp = 0

end

47

cbztrans.forsubroutine cbztransimplicit real*8(a-h,o-z)c------------------------------------------------------------------c

subroutine cbz transcdo overall calculation for one step.

c------------------------------------------------------------------include 'contrl.inc'include 'cconv.inc'include 'iounit.inc'
print *, "in the trans"
10 continue
odelt = delt
nmfail = 0
20 continue
call timstp( iendp, iofail, nmfail )dto = odelt
im100 = oitno
im100x = im100
if ( nstep .eq. 0 ) im100x = 0
oitno = 0iofail = 0
if ( nstep .eq. 0 ) thencall edit
call graf
end ifcall prep
ipakon = 030 continue
call hout( iofail , nmfail )if ( oitno .ne. 0 ) go to 31
iofail = 0
go to 3048

31 continueif( iofail .ne. 0 ) go to 20
it= oitnooitno = -1
lbkpst = .false.
call post
if ( .not.lbkpst ) go to 32
iofail = 1
nmfail = nmfail + 1
oitno = -100
write(iout,2001) jbkpst , nstep
call post
go to 2032 continue
oitno = itnstep = nstep + 1
timet = timet + deltcall pstepq
call edit
return2001 format(1h0,'***** back up forced in post at cell ',i5,*' time step',i5,'.')

1000 format(15h end of problem )end49

cbzheatinput.forsubroutine cbzheatinput( pow, dist, pinreg )c subroutine to input the data to TRAC
implicit real*8 ( a-h , o-z )include 'setprm.inc'include 'iounit.inc'include 'fltab.inc'include 'heatpt.inc'integer pinreg
double precision pow
double precision dist( pinreg )
integer i
rpwtb(1,1) = 0.0
rpwtb(2,1) = pow
rpwtb(1,2) = 99
rpwtb(2,2) = pow
end50cbztemp.for
subroutine cbztemp( temp,pinreg,cells,alpha )c subroutine to putout the data to CBZ
implicit real*8 ( a-h , o-z )include 'setprm.inc'include 'iounit.inc'include 'fltab.inc'include 'heatpt.inc'include 's1dptr.inc'integer pinreg
integer ftemp
integer cells
double precision alpher
double precision temp( pinreg+1 )integer i
ftemp = cells - 1
do i = 1, pinregctemp(i) = rftn( i + ftemp*pinreg )temp(i) = rftn( i )end dotemp(pinreg+1) = tln(cells)
alpha = alpn(cells)
end51