DEVELOPMENT OF AN IMPLICIT NUMERICAL METHOD FOR TWO-PHASE FLOW ANALYSIS OF PRESSURIZED WATER REACTORS

P. Jobczyk; H.Y. Yoon

doi:10.6112/kscfe.2025.30.3.102

Preview

Original Article

Journal of Computational Fluids Engineering. 30 September 2025. 102-122
https://doi.org/10.6112/kscfe.2025.30.3.102

DEVELOPMENT OF AN IMPLICIT NUMERICAL METHOD FOR TWO-PHASE FLOW ANALYSIS OF PRESSURIZED WATER REACTORS

가압경수로 2상유동 해석을 위한 내재적 수치해법 개발

P. Jobczyk¹

H.Y. Yoon¹^†

Jobczyk Patryk¹

윤 한영¹^†

¹Department of Nuclear Power Plant Engineering, KEPCO International Nuclear Graduate School

¹한국전력국제원자력대학원대학교 원자력산업학과

^{*Corresponding Author}

License (open-access, http://creativecommons.org/licenses/by-nc/4.0):

This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0) which permits unrestricted noncommercial use, distribution, and reproduction in anymedium, provided the original work is properly cited.

ABSTRACT

This study presents the development of an Implicit numerical method implemented for two-phase flows in Pressurized Water Reactors, to improve stability and decrease simulation time. It focuses on solving governing equations for the two-fluid model by treating the phase and spatial terms implicitly using implementation in CUPID. The verification of this method is performed on a blowdown and refill scenario, to verify the ability of numerical scheme to simulate rapid pressure decrease and phase change phenomena. The validation is performed against FLECHT-SEASET reflood experiment. showing accurate results and reduced calculation time. The implicit method offers precise solution with efficient simulation with high CFL number and has potential for further use and development in nuclear thermal hydraulic field.

Keywords

Two-phase flow

Two-fluid model

Nuclear thermal hydraulics

Implicit numerical method

Computational fluid dynamics

CFD

키워드

2상유동

2유체모델

원자로열수력

내재적 수치해법

전산유체역학

MAIN

Nomenclature
1. 개 요
2. 수치해석 방법
2.1 지배 방정식
2.2 반내재적(semi-implicit) 수치해법
2.3 내재적(implicit) 수치해법
3. 냉각재상실사고(LOCA) 개념 문제 확인(Verification) 해석
4. 재관수(Reflood) 실험에 대한 검증(Validation)
5. 결 론
Appendix

Nomenclature

$A$ 압력방정식 주대각 요소 계수

$B$ 압력방정식 원천항

$C$ 상수 계수

$d \vec{r}$ 거리 벡터

$e$ 내부에너지

$\vec{F}$ 운동량 벡터

$G$ 압력방정식 비대각 요소 중 기체 계수

$H$ 열전달계수

$h$ 엔탈피

$k$ 열전도도

$L$ 압력방정식 비대각 요소 중 액체 계수

$m$ 운동량방정식 행렬 계수

$\vec{n}$ 수직 벡터

$o$ 에너지방정식 행렬 계수

$P$ 압력

$\vec{q}$ 열속 벡터

$q^{'''}$ 단위체적당 열생성율

$r$ 에너지방정식 원천항

$S$ 계산 셀 면적

$\vec{s}$ 운동량방정식 원천항 벡터

$T$ 온도

$t$ 시간

$\vec{u}$ 속도 벡터

$V$ 계산 셀 체적

Greek Letter

𝛼 체적율

𝛽 압력구배항 계수

𝛿 변화량

𝜙 종속 변수

𝜌 밀도

𝜇 점성

𝜓 체적 유량

𝛤 상변화율

Superscripts

$(n)$ n 번째 반복 계산

$n$ 이전 차분 시간

$n$ +1 새로운 차분 시간

$s$ 포화 상태

Subscripts

$d$ 계면 마찰

$f$ 계산 셀 면

$g$ 기체

$g l$ 액체외 비응축성기체 사이의 계면

$i$ i 번째 계산 셀

$i j$ i 번째 계산 셀과 j 번째 계산 셀 사이

$i g$ 기체 계면

$i l$ 액체 계면

$i n t$ 기체와 액체 사이의 계면

$j$ j 번째 계산 셀

$k$ 액체 또는 기체

$v$ 증기

$v m$ 가상 질량

1. 개 요

가압경수로(PWR)의 열수력 해석은 복잡한 기체-액체 계면 구조와 대규모 상변화를 포함한 2상(two-phase) 유동 현상을 다루기 때문에 많은 기술적 난제들을 동반한다. 예를 들어, 원자로 안전성 평가를 위한 대표적 가상사고 중 하나인 냉각재 상실 사고(LOCA) 조건에서는 원자로 냉각재 시스템 압력이 150기압에서 대기압으로 몇 분 이내로 급강하하며, 이로 인해 급격한 상변화(phase change)가 발생한다. 또한, LOCA 후반부에서는 핵연료봉을 냉각하기 위해 안전 주입수(safety injection water)가 원자로 압력 용기에 도달하고, 이때 과열된 핵연료봉 표면에서 대량의 벽면 비등과 증기 응축이 발생한다. 이와 같은 원자로 2상 유동 현상에 대한 정확한 수치해석은 매우 어려운 분야이다.

신뢰할 수 있는 원자로 열수력 해석은 APR1400[1]과 같은 PWR의 성능과 안전을 확보하기 위하여 필수적 요소이다. PWR의 2상 유동 해석에는 일반적으로 증기와 액체로 구성된 2유체(two-fluid) 모델[2]을 적용한다. 현재 상업용 PWR의 안전과 성능 분석을 위해 사용되는 해석 코드로는 RELAP5[3], CATHARE[4], MARS[5] 등이 활용된다. 이들 2상 유동 해석 코드들은 원자로 냉각재 계통의 복잡성과 계산 시간을 고려하여 단순화된 1차원 모델을 기반으로 하고 있는데, 모델의 예측 불확실도를 고려하기 위하여 보수적 해석 가정 및 입력 모델이 적용된다. 1차원 해석 코드의 경우 대류항 및 확산항 등 공간 미분항의 역할이 매우 제한적으로 , 상변화에 대하여 내재적 해석 방법을 적용하고 공간 미분에 대해서는 양해법을 적용하고 있다. 이에 따라서, 수치해석 방법으로는 반내재적(semi-implicit) 방법이 적용되며, 차분 시 구간의 크기는 CFL(Courant Friedrichs Levy) 수에 의해 제한된다. 1차원 안전 및 성능 해석의 경우 제어 체적(계산 격자) 수가 수백 개 수준으로 적기 때문에 차분 시 구간의 크기 제한은 큰 문제가 되지 않는다. 그러나, 최근 원자로의 3차원 해석 열수력 해석이 확대됨에 따라서 반내재적 방법 적용으로 인한 차분 시 구간 크기 제한에 대한 개선이 필요한 실정이다.

원자로 3차원 2상 유동 해석을 위하여 기존 1차원 해석 코드를 확장하여, RELAP5-3D[6], CATHARE-3D[7], MARS- 3D[8] 등이 개발되었으나, 이들 해석 코드에서는 계산 격자수가 약 10,000개 수준으로 해석 해상도(resolution)를 크게 높이기에는 한계가 있다. 더욱 정교한 3차원 2상 유동 해석을 위하여 FLUENT[9], CFX[10], STAR-CCM[11] 등의 상용 CFD 코드의 적용이 가능하다. 다만, 상용 CFD 코드들은 대류항 및 확산항 등에 대해서는 내재적인 해석 방법을 적용하지만 상변화는 양해법을 적용하고 있어서, PWR LOCA에서와 같이 대규모 상변화 현상을 동반한 2상 유동 해석에 적용하기에는 한계가 있다.

비교적 최근에, NEPTUNE-CFD[12], CUPID[13,14,15]와 같은 다상(multi-phase) 유동 CFD 해석 코드들이 PWR의 3차원 2상 유동 해석을 위해 개발되었다. CUPID는 3차원 2유체(two-fluid) 모델을 채택하였고, 지배 방정식은 비정형 격자(unstructured mesh)를 사용하는 유한체적법(Finite Volume Method, FVM)으로 이산화되었다. CUPID에는 기존 원자로 안전 해석 코드에서 사용된 반내재적 방법이 적용되었으며, 이로 인하여 미세한 격자를 사용할 경우 CFL 수 제한으로 인하여 차분 시 구간 크기가 매우 작게 설정되는 결과를 초래한다.

내재적 수치해석 방법(implicit numerical solution scheme)은 CFL 조건을 초과하는 큰 시 구간 크기를 허용함으로써, PWR의 3차원 2상 유동 해석을 실용적인 계산 시간 내에 수행하는 것이 가능하다. 또한, 1차원 안전 해석의 경우, 내재적 수치해석 방법은 많은 입력 조건에 대한 대량의 해석 경우들이 요구되는 불확실성 분석(uncertainty analysis)에 소요 되는 시간을 크게 단축할 수 있다.

본 연구에서는 3차원 2유체 모델에 적용이 가능한 내재적 수치해석 방법을 제시하고 이를 CUPID 코드에 구현하였다. 상변화 모델은 운동량 방정식으로부터 유도된 압력 방정식과 결합하여 내재적으로 계산되며, 동시에 지배 방정식의 대류항(convection) 및 확산항(diffusion) 또한 내재적으로 처리된다. 시스템 행렬의 크기를 줄이기 위하여 계산 절차를 두 단계로 분리하였는데, 이는 각각 상-연계(phase-link)와 공간-연계(space-link) 단계로 구성된다. 계산 단계 분리에서 발생할 수 있는 수치 해의 차이를 최소화하기 위하여 반복 계산 기법을 적용한다.

제안된 내재적 수치해법의 확인(verification)을 위하여 PWR의 LOCA에서 발생하는 대규모 상변화를 모사하는 개념적 2상 유동 문제를 설정하였다. 개념 문제의 과도현상은 고압(150기압)의 단상 액체가 급격히 감압 되어 플래싱(flashing) 현상에 의한 대규모 상변화가 발생하는 것으로 시작된다. 압력이 특정 설정값 이하로 감소하면, 냉각수 주입(cold water injection)이 시작되어 수위(water level)가 회복된다. 본 논문에서는, 위의 시나리오를 바탕으로 반내재적 방법과 내재적 방법을 각각 적용한 계산을 수행하고 그 계산 결과를 비교 분석하였다.

다음으로, 내재적 수치해법의 검증(validation)을 위하여 PWR LOCA의 재관수(reflood) 현상을 규명하기 위한 실험 FLECHT- SEASET[16]에 대한 비교 해석을 수행하였다. 수치 해의 정확성 및 해석 시간 비교를 위해 1차원 및 3차원 격자를 모두 사용하였다.

제2장에서는 지배 방정식, 반내재적 수치해법, 내재적 수치해법을 포함한 수치 모델을 설명한다. 제3장에서는 급격한 감압 현상을 동반한 개념 문제에 관한 확인 해석 결과를 제시하고, 제4장에서는 LOCA 재관수 실험에 대한 검증 해석 결과를 설명한다. 제5장에서는 본 연구의 요약과 결론을 제시하였다.

2. 수치해석 방법

2.1 지배 방정식

원자로 내 2상 유동 해석에서는 일반적으로 2유체 모델이 적용되는데, 이 모델에서는 기상과 액상에 대하여 질량, 운동량, 에너지 보존 방정식이 각각 정의된다. 기상 및 액상에 대한 질량 보존방정식은 아래와 같이 표현된다.

(1)

\frac{\partial}{\partial t} (α_{k} ρ_{k}) + \nabla ∙ (α_{k} ρ_{k} {\vec{u}}_{k}) = C_{k} ∙ Γ_{v}

위의 방정식에서, $k = g$ 는 기체 상태이고 $k = l$ 은 액체 상태를 나타낸다. 아래 첨자 $v$ 는 증기 상태를 나타낸다. 기상의 경우 $C_{g}$ =1 이며, 액상의 경우 $C_{l}$ -1 이다. 기상은 증기와 공기, 수소와 같은 비응축성(non-condensable) 기체를 포함한다. $α_{k}$ , $ρ_{k}$ , $u_{k}$ 는 $k$ -상의 체적 분율, 밀도, 속도를 나타낸다. $Γ_{v}$ 는 단위 부피당 상 변화율로 정의하고 아래와 같은 식으로 표현된다.

(2)

Γ_{v} = - \frac{H_{i g} [T^{s} (P_{v}) - T_{g}] + H_{i l} [T^{s} (P_{v}) - T_{l}]}{h_{g}^{*} - h_{l}^{*}},

여기서 $P_{v}$ 및 $T^{s} (P_{v})$ 는 각각 증기 압력 및 $P_{v}$ 에서의 포화 온도이고, $H_{i k}$ 및 $T_{k}$ 는 $k$ -상의 계면 열전달 계수 및 온도이다. $h_{k}^{*}$ 는 $k$ -상의 엔탈피이며, 기화의 경우 ( $Γ_{v}$ ≥0), $(h_{g}^{*}, h_{l}^{*}) = (h_{g}^{s}, h_{l})$ 로 정의되고 응축의 경우( $Γ_{v}$ <0), $(h_{g}^{*}, h_{l}^{*}) = (h_{g}, h_{l}^{s})$ 로 정의된다. 위 첨자 “ $s$ ” 는 포화 조건에서의 값을 나타낸다. 상 변화율( $Γ_{v}$ )은 압력 및 온도의 함수이므로 내재적으로 계산되어야 한다.

기체 상태는 증기와 비응축성 기체의 혼합물이고, 비응축성 기체 거동 예측을 위하여 추가적인 질량 보존 방정식이 다음과 같이 정의된다. 여기서 $X_{n}$ 은 비응축성 기체의 질량 분율이다.

(3)

\frac{\partial}{\partial t} (α_{g} ρ_{g} X_{n}) + \nabla ∙ (α_{g} ρ_{g} X_{n} {\vec{u}}_{g}) = 0

수치해석의 편의를 위하여 식 (3)의 좌변은 기체 상태 질량 보존 방정식(1)을 이용하여 비보존(non-conservaative) 형태로 전개된다.

(4)

α_{g} ρ_{g} \frac{\partial X_{n}}{\partial t} + \nabla ∙ (α_{g} ρ_{g} X_{n} \vec{u_{g}}) - X_{n} \nabla ∙ (α_{g} ρ_{g} \vec{u_{g}}) = - X_{n} Γ_{v}

운동량 보존 방정식은 기체 혼합물과 액체 상태에 대하여 각각 독립적으로 정의된다. 비응축성 기체의 운동량 방정식은 기체 속도와 같다고 가정하여 별도로 고려하지 않는다.

(5)

\frac{\partial}{\partial t} (α_{k} ρ_{k} {\vec{u}}_{k}) + \nabla ∙ (α_{k} ρ_{k} {\vec{u}}_{k} {\vec{u}}_{k}) = - α_{k} \nabla P + \nabla ∙ (α_{k} μ_{k} \nabla \vec{u_{k}}) + α_{k} ρ_{k} \vec{g} + C_{k} (\vec{F^{i n t}} + \vec{F^{i m}})

위의 식에서, $\vec{F^{i n t}}$ 와 ${\vec{F}}^{v m}$ 은 각각 계면 운동력(interfacial force)와 가상 질량력(virtual mass force)을 의미한다. $\vec{F^{i n t}}$ 는 계면 항력, 난류 확산력, 벽면 윤활력 등을 포함하는데, 이 중에서 계면 항력은 수치 안정성 확보를 위하여 반드시 내재적으로 계산되어야 한다. 계면 항력은 일반적으로 다음과 같은 식으로 정의할 수 있는데 여기서 $C_{d}$ 는 항력 계수이다.

(6)

{\vec{F}}^{i n t} = C_{d} ({\vec{u}}_{l} - {\vec{u}}_{g})

가상 질량력 또한 내재적으로 계산되어야 하며, 다음과 같이 표현된다. 여기서 $C_{v m}$ 은 가상 질량 계수(virtual mass coefficient)이다. 가상질량력은 기상-액상 혼합 유체가 가속 또는 감속되는 상황에서 수치 안정성을 확보하는 역할을 하므로, 시간 미분항이 주요한 역할을 하고 있으며 공간 미분을 고려한 경우와의 차이는 매우 작다. 따라서, 대부분의 2유체 모델 기반 코드에서는 시간 미분항만을 고려하고 있다.

(7)

{\vec{F}}^{v m} = C_{v m} \frac{\partial ({\vec{u}}_{l} - {\vec{u}}_{g})}{\partial t}

식 (4)와 같은 방식으로, 식 (5)의 좌변도 식 (1)의 질량 보존 방정식을 이용하여 아래와 같이 비보존 형태로 전개할 수 있다.

(8)

\frac{\partial}{\partial t} (α_{k} ρ_{k} {\vec{u}}_{k}) + \nabla ∙ (α_{k} ρ_{k} {\vec{u}}_{k} {\vec{u}}_{k}) = α_{k} ρ_{k} \frac{\partial \vec{u_{k}}}{\partial t} + \nabla ∙ (α_{k} ρ_{k} {\vec{u}}_{k} {\vec{u}}_{k}) + \vec{u_{k}} \frac{\partial (α_{k} ρ_{k})}{\partial t} = α_{k} ρ_{k} \frac{\partial {\vec{u}}_{k}}{\partial t} + \nabla ∙ (α_{k} ρ_{k} {\vec{u}}_{k} {\vec{u}}_{k}) - \vec{u_{k}} \nabla ∙ (α_{k} ρ_{k} {\vec{u}}_{k}) + \vec{u_{k}} Γ_{k}

식 (6), (7), (8)을 식 (5)에 대입하면, 운동량 방정식은 다음과 같이 표현된다.

(9)

α_{k} ρ_{k} \frac{\partial \vec{u_{k}}}{\partial t} + \nabla ∙ (α_{k} ρ_{k} {\vec{u}}_{k} {\vec{u}}_{k}) - {\vec{u}}_{k} \nabla ∙ (α_{k} ρ_{k} {\vec{u}}_{k}) = - α_{k} \nabla P + \nabla ∙ (α_{k} μ_{k} \nabla {\vec{u}}_{k}) + α_{k} ρ_{k} \vec{g} + C_{k} ∙ [C_{d} (\vec{u_{l}} - \vec{u_{g}}) + C_{v m} \frac{\partial (\vec{u_{l}} - \vec{u_{g}})}{\partial t} - \vec{u_{k}} Γ_{v}]

에너지 보존 방정식은 증기 및 비응축성 기체 사이의 열적 평형(thermal equilibrium)을 가정하여 기체 혼합물과 액체에 대하여 두 개의 방정식을 적용한다.

(10)

\frac{\partial}{\partial t} (α_{k} ρ_{k} e_{k}) + \nabla ∙ (α_{k} ρ_{k} e_{k} \vec{u_{k}}) = - P \frac{\partial α_{k}}{\partial t} - P \nabla ∙ (α_{k} \vec{u_{k}}) + \nabla ∙ (α_{k} q_{k}) + q_{k}^{'''} + H_{i k} (T^{s} - T_{k}) + C_{k} ∙ [Γ_{v} h_{k}^{*} - H_{g l} (T_{g} - T_{l})]

위 식에서, $e_{k}$ , $h_{k}$ 와 $q_{k}^{'''}$ 는 $k$ -상의 내부에너지, 엔탈피, 체적 열원(volumetric heat source)을 의미한다. ${\vec{q}}_{k}$ 는 $k$ -상의 열속으로, ${\vec{q}}_{k} = k_{k} \nabla T_{k}$ 와 같이 정의된다. $H_{g l}$ 은 비응축성 기체와 액체 사이의 열전달 계수를 의미한다. 식 (10)의 좌변에 식 (1)을 적용하면, 다음과 같이 비보존 형태로 전개된다.

(11)

\frac{\partial}{\partial t} (α_{k} ρ_{k} e_{k}) + \nabla ∙ (α_{k} ρ_{k} e_{k} {\vec{u}}_{k}) = α_{k} ρ_{k} \frac{\partial e_{k}}{\partial t} + \nabla ∙ (α_{k} ρ_{k} e_{k} {\vec{u}}_{k}) - e_{k} \nabla ∙ (α_{k} ρ_{k} {\vec{u}}_{k}) + e_{k} Γ_{k}

마찬가지로, 식 (11)의 우변에 있는 첫 번째 항과 두 번째 항도 식 (1)의 질량 보존 방정식을 이용하면, 다음과 같이 표현 할 수 있다.

(12)

- P \frac{\partial α_{k}}{\partial t} - P \nabla ∙ (α_{k} \vec{u_{k}}) = - \frac{P}{ρ_{k}} [Γ_{k} - α_{k} \frac{\partial ρ_{k}}{\partial t} - \nabla ∙ (α_{k} ρ_{k} \vec{u_{k}}) + ρ_{k} \nabla ∙ (α_{k} \vec{u_{k}})]

식 (11)과 (12)를 식 (10)에 대입하면, 최종 에너지 보존 방정식은 아래와 같이 표현된다.

(13)

α_{k} ρ_{k} \frac{\partial e_{k}}{\partial t} + \nabla ∙ (α_{k} ρ_{k} e_{k} {\vec{u}}_{k}) - e_{k} \nabla ∙ (α_{k} ρ_{k} {\vec{u}}_{k}) = - \frac{P α_{k}}{ρ_{k}} \frac{\partial ρ_{k}}{\partial t} + \frac{P}{ρ_{k}} [\nabla ∙ (α_{k} ρ_{k} {\vec{u}}_{k}) - ρ_{k} \nabla ∙ (α_{k} {\vec{u}}_{k})] + H_{i k} (T^{s} - T_{k}) + C_{k} ∙ [Γ_{v} (h_{k}^{*} - e_{k} - \frac{P}{ρ_{k}}) - H_{g l} (T_{g} - T_{l})] + \nabla ∙ (α_{k} q_{k}) + q_{k}^{'''}

위에서 기술한 지배 방정식들은 Fig. 1에 나타낸 비정형 격자를 이용하여 유한체적법(FVM)으로 이산화된다. 모든 주요 변수들은 ( $α_{g}, α_{l}, \vec{u_{g}}, \vec{u_{l}}, e_{g}, e_{l}, X_{n}$ ) 셀 중심에서 정의된다.

https://cdn.apub.kr/journalsite/sites/kscfe/2025-030-03/N0500300307/images/jkscfe_2025_303_102_F1.jpg

Fig. 1.

Unstructured mesh and geometric parameters used for discretization

지배 방정식에 나타나는 $ρ_{k}$ , $T_{k}$ , $T^{s}$ 와 같은 종속 변수들은 $P$ , $e_{k}$ , $X_{n}$ 와 같은 독립 변수들의 변화량에 대해 선형화하여 표현된다.

(14)

ϕ_{g}^{n + 1} = ϕ_{g}^{n} + {(\frac{\partial ϕ_{g}}{\partial P})}^{n} (P^{n + 1} - P^{n}) + {(\frac{\partial ϕ_{g}}{\partial e_{g}})}^{n} (e_{g}^{n + 1} - e_{g}^{n}) + {(\frac{\partial ϕ_{g}}{\partial X_{n}})}^{n} (X_{n}^{n + 1} - X_{n}^{n})

(15)

ϕ_{l}^{n + 1} = ϕ_{l}^{n} + {(\frac{\partial ϕ_{l}}{\partial P})}^{n} (P^{n + 1} - P^{n}) + {(\frac{\partial ϕ_{l}}{\partial e_{l}})}^{n} (e_{l}^{n + 1} - e_{l}^{n})

$ϕ = ρ, T, T^{s}$ 는 밀도, 온도 및 포화 온도를 나타내며, 위 첨자 $n, n + 1$ 은 차분 시간 단계( $t^{n + 1} = t^{n} + δ t$ )를 의미한다. 시간 이산화에는 1차 양해법(explicit Euler method)이 적용된다.

2.2 반내재적(semi-implicit) 수치해법

반내재적 수치해법의 첫 번째 계산 단계에서는 운동량 방정식을 양해법으로 계산한다. 식 (8)의 이산화된 형태는 $i - t h$ 번째 셀에 대해 다음과 같이 표현될 수 있다.

(16)

{[\begin{matrix} m_{g g} m_{g l} \\ m_{l g} m_{l l} \end{matrix}]}_{i} {[\begin{matrix} {\vec{u}}_{g}^{*} \\ {\vec{u}}_{l}^{*} \end{matrix}]}_{i} = - {[\begin{matrix} V α_{g} \\ V α_{l} \end{matrix}]}_{i} \nabla P_{i} + {[\begin{matrix} {\vec{s}}_{g} \\ {\vec{s}}_{l} \end{matrix}]}_{i}

위 식에서 위 첨자 “*” 는 양해법에 따른 운동량 계산 후 업데이트된 값을 나타낸다. 위 첨자가 표시되지 않은 경우의 해당 값은 $t^{n}$ 에서의 값을 나타낸다. $V_{i}$ 는 $i - t h$ 번째 셀의 체적을 나타내며, 행렬 요소 $m_{11}$ , $m_{12}$ , $m_{21}$ , $m_{22}$ , ${\vec{s}}_{g}$ , ${\vec{s}}_{l}$ 은 다음과 같이 정의된다.

(17)

m_{g g, i} = {(α_{g} ρ_{g} + C_{v m})}_{i} \frac{V_{i}}{Δ t} + V_{i} C_{d}, m_{g l, i} = - C_{v m, i} \frac{V_{i}}{Δ t} - V_{i} C_{d} m_{l g, i} = - C_{v m, i} \frac{V_{i}}{Δ t} - V_{i} C_{d}, m_{l l, i} = {(α_{l} ρ_{l} + C_{v m})}_{i} \frac{V_{i}}{Δ t} + V_{i} C_{d}

(18)

{\vec{s}}_{k, i} = {(α_{k} ρ_{k} {\vec{u}}_{k})}_{i} \frac{V_{i}}{Δ t} + \sum_{f = 1}^{N_{neigh}} {(α_{k} ρ_{k})}_{f} \frac{{\vec{u}}_{k, j} - {\vec{u}}_{k, i}}{\vec{d r_{i j}} ∙ {\vec{n}}_{f}} S_{f} - \sum_{f = 1}^{N_{neigh}} {(α_{k} ρ_{k})}_{f} ({\vec{u}}_{k, f} - {\vec{u}}_{k, i}) ψ_{k, f} - C_{k} ∙ V_{i} {[C_{v m} (\vec{u_{l}} - \vec{u_{g}}) \frac{1}{Δ t} + \vec{u_{k}} Γ_{v}]}_{i} + V_{i} α_{k} ρ_{k} \vec{g}

식 (18) 우변의 두 번째 및 세 번째 항은 셀 면들에 대한 합으로, $f$ =1 부터 $f = N_{n e i g h}$ 까지의 합을 나타낸다. 여기서, $N_{neigh}$ 는 $i - t h$ 번째 셀이 소유한 면의 개수를 의미한다. 예를 들면, Fig. 1에서 $N_{neigh}$ 은 3이다. $ψ_{k, f}$ 는 표면 벡터 ${\vec{S}}_{f}$ 와 표면 속도 ${\vec{u}}_{k, f}$ 의 내적으로 정의된 체적 유량( $ψ_{k, f} = {\vec{u}}_{k, f} ∙ {\vec{S}}_{f}$ ) 이다. ${\vec{u}}_{k, f}$ 는 셀 중심 값을 이용하여 내삽된 값으로 ${\vec{u}}_{k, f} = w ∙ {\vec{u}}_{k, i} + (1 - w) ∙ {\vec{u}}_{k, j}$ 로 정의된다. 여기서 $w$ 는 인접한 두 셀 사이의 거리에 대한 가중치이다.

식 (16)을 이용하여 양해법에 따른 운동량 계산 후, 새로운 시간( $t^{n + 1} = t^{n} + δ t$ )에서의 속도 벡터는 아래와 같이 결정된다.

(19)

{[\begin{array}{l} m_{g g} m_{g l} \\ m_{l g} m_{l l} \end{array}]}_{i} {[\begin{array}{l} {\vec{u}}_{g}^{n + 1} \\ {\vec{u}}_{l}^{n + 1} \end{array}]}_{i} = {[\begin{array}{l} m_{g g} m_{g l} \\ m_{l g} m_{l l} \end{array}]}_{i} {[\begin{array}{l} {\vec{u}}_{g}^{*} \\ {\vec{u}}_{l}^{*} \end{array}]}_{i} - {[\begin{array}{l} V α_{g, i} \\ V α_{l, i} \end{array}]}_{i} \nabla δ P_{i}

여기서 $δ P_{i}$ 는 $δ P_{i} = P_{i}^{n + 1} - P_{i}^{n}$ 로 정의된다. 식 (19)는 아래와 같이 간단히 표현될 수 있다.

(20)

{\vec{u}}_{k, i}^{n + 1} = {\vec{u}}_{k, i}^{*} - β_{k, i} \nabla δ P_{i}

여기서, $β_{k, i}$ .는 $β_{k, i} = V_{i} \sum_{m = g, l} m_{m k}^{- 1} α_{m, i}$ 로 정의된다. 식 (20)의 양변에 발산( $\nabla ∙$ )연산자를 적용하고 제어 체적에 대해 적분한 후, 결과를 이산화하면 다음과 같은 방정식이 도출된다.

(21)

\sum_{f = 1}^{N_{reigh}} ψ_{k, f}^{n + 1} = \sum_{f = 1}^{N_{neigh}} ψ_{k, f}^{*} - \sum_{f = 1}^{N_{reigh}} β_{k, f} \frac{δ P_{j} - δ P_{i}}{\vec{d r_{i j}} ∙ \vec{n_{f}}} S_{f}

식 (21)에는 2개의 미지수( $ψ_{k, f}^{n + 1}$ , $δ P_{i}$ )가 포함되어 있으므로, 질량 보존 방정식을 활용하여 $ψ_{k, f}^{n + 1}$ 소거한다. 이를 위하여, 기체와 액체의 질량 방정식을 합산하면 다음과 같은 방정식을 얻는다.

(22)

\frac{1}{ρ_{g}} \nabla ∙ (α_{g} ρ_{g} {\vec{u}}_{g}) + \frac{1}{ρ_{l}} \nabla ∙ (α_{l} ρ_{l} {\vec{u}}_{l}) = Γ_{v} (\frac{1}{ρ_{g}} - \frac{1}{ρ_{l}}) - \frac{α_{g}}{ρ_{g}} \frac{\partial ρ_{g}}{\partial t} - \frac{α_{l}}{ρ_{l}} \frac{\partial ρ_{l}}{\partial t}

위의 방정식을 적분하고 이산화하면 다음과 같은 방정식이 된다.

(23)

\frac{1}{V_{i}} \sum_{k = g, l} [\frac{1}{ρ_{k}} \sum_{f = 1}^{N_{neigh}} {(α_{k} ρ_{k})}_{f} ψ_{k, f}^{n + 1}] = Γ_{v, i}^{n + 1} {(\frac{1}{ρ_{g}} - \frac{1}{ρ_{l}})}_{i} - {(\frac{α_{g}}{ρ_{g}} \frac{δ ρ_{g}}{δ t} + \frac{α_{l}}{ρ_{l}} \frac{δ ρ_{l}}{δ t})}_{i}

위 식에서, $Γ_{v}^{n + 1}$ 와 $δ ρ_{k}$ 는 $δ P$ 에 대하여 다음과 같이 선형화한다.

(24)

Γ_{v}^{n + 1} = - \frac{H_{i g} [T^{s, n + 1} - T_{g}^{n + 1}] + H_{i l} [T^{s, n + 1} - T_{l}^{n + 1}]}{h_{g}^{*} - h_{l}^{*}} = - \frac{H_{i g} [T^{s} - T_{g}] + H_{i l} [T^{s} - T_{l}]}{h_{g}^{*} - h_{l}^{*}} - \frac{1}{h_{g}^{*} - h_{l}^{*}} [H_{i g} (\frac{\partial T^{s}}{\partial P} - \frac{\partial T_{g}}{\partial P}) + H_{i l} (\frac{\partial T^{s}}{\partial P} - \frac{\partial T_{l}}{\partial P})] δ P

(25)

δ ρ_{k} = (\frac{\partial ρ_{k}}{\partial P}) δ P

식 (21), (24), (25)를 식 (23)에 대입한 후 정리하면 다음과 같은 방정식을 얻는다.

(26)

A δ P_{i} = \sum_{f = 1}^{N_{neigh}} G [ψ_{g, f}^{*} - β_{g, f} \frac{δ P_{j} - δ P_{i}}{\vec{d r_{i j}} ∙ \vec{n_{f}}} S_{f}] + \sum_{f = 1}^{N_{neigh}} L [ψ_{l, f}^{*} - β_{l, f} \frac{δ P_{j} - δ P_{i}}{\vec{d r_{i j}} ∙ \vec{n_{f}}} S_{f}] + B

위 식에서, $A$ , $B$ , $G$ , $L$ 계수들은 다음과 같이 정의된다.

(27)

A = - \frac{V_{i}}{h_{g}^{*} - h_{l}^{*}} (\frac{1}{ρ_{g}} - \frac{1}{ρ_{l}}) [H_{i g} (\frac{\partial T^{s}}{\partial P} - \frac{\partial T_{g}}{\partial P}) + H_{i l} (\frac{\partial T^{s}}{\partial P} - \frac{\partial T_{l}}{\partial P})] - \frac{V_{i} α_{g}}{ρ_{g} δ t} (\frac{\partial ρ_{g}}{\partial P}) - \frac{V_{i} α_{l}}{ρ_{l} δ t} (\frac{\partial ρ_{l}}{\partial P}) B = \frac{V_{i}}{h_{g}^{*} - h_{l}^{*}} (\frac{1}{ρ_{g}} - \frac{1}{ρ_{l}}) [H_{i g} (T^{s} - T_{g}) + H_{i l} (T^{s} - T_{l})] G = {(α_{g} ρ_{g})}_{f}, L = {(α_{l} ρ_{l})}_{f}

식 (26)은 아래와 같이 좀 더 간소화된 형태로 표현될 수 있다.

(28)

[1 + \sum_{f = 1}^{N_{neigh}} C C] δ P_{i} - \sum_{f = 1}^{N_{neigh}} C C δ P_{j} = B B

위에서, $B B = A^{- 1} B + A^{- 1} \sum_{f = 1}^{N_{neigh}} (G ψ_{g, f}^{*} + L ψ_{l, f}^{*})$ , $C C = - A^{- 1} (G β_{g, f} + L β_{l, f}) \frac{S_{f}}{\vec{d r_{i j}} ∙ \vec{n_{f}}}$ 이다. 식 (28)은 선형화된 압력 방정식을 나타내는데, BiCGSTAB[18] 반복 솔버(solver)를 적용하여 계산한다.

식 (28)에 의한 압력 계산 후, 속도 벡터 ( ${\vec{u}}_{k}^{n + 1}$ )와 체적 유량 ( $ψ_{k}^{n + 1}$ )은 식 (20)과 (21)을 사용하여 계산된다. 다음 단계로, 기포율( $α_{g}^{n + 1}$ ), 에너지( $e_{k}^{n + 1}$ ), 비응축성 기체 질량 분율( $X_{n}^{n + 1}$ )을 결정하기 위해 질량 보존 방정식과 에너지 보존 방정식을 활용한다. 먼저, 기체와 액체의 질량 보존 방정식이 활용하여 기체 체적 분율( $α_{g}^{n + 1}$ )을 계산한다.

(29)

\frac{\partial (α_{g} ρ_{g})}{\partial t} + \frac{\partial (α_{l} ρ_{l})}{\partial t} + \nabla ∙ (α_{g} ρ_{g} {\vec{u}}_{g}) + \nabla ∙ (α_{l} {\vec{ρ}}_{l} {\vec{u}}_{l}) = 0

위의 방정식을 이산화하면, 다음과 같다.

(30)

V_{i} α_{g, i}^{n + 1} (ρ_{g}^{n + 1} - ρ_{g, i}^{n + 1}) δ t^{- 1} = V_{i} [{(α_{g} ρ_{g})}_{i}^{n} - {(α_{l} ρ_{l})}_{i}^{n} - ρ_{g, i}^{n + 1}] δ t^{- 1} - \sum_{f = 1}^{N_{neigh}} {(α_{g} ρ_{g})}_{f}^{n} ψ_{g, f}^{n + 1} - \sum_{f = 1}^{N_{neigh}} {(α_{l} ρ_{l})}_{f}^{n} ψ_{l, f}^{n + 1}

식 (30)을 계산하여 기체 체적 분율 ( $α_{g}^{n + 1}$ )을 계산하면, 액체 체적 분율은 $α_{l}^{n + 1} = 1 - α_{g}^{n + 1}$ 에 의해 얻어진다.

다음으로, 비응축성 기체의 질량 보존 방정식과 에너지 보존 방정식(식 (4)및 (13))을 이산화하면 3개의 미지수 ( $δ e_{g}, δ e_{l}, δ X_{n}$ )에 대해 다음과 같이 정리할 수 있다.

(31)

\bar{O} \bar{X} = \bar{R}

위 방정식에 사용된 행렬은 다음과 같으며, 계수들은 부록에 정리되어 있다.

(32)

\bar{O} = [\begin{array}{l} o_{11} & o_{12} & o_{13} \\ o_{21} & o_{22} & o_{23} \\ o_{31} & o_{32} & o_{33} \end{array}], \bar{X} = [\begin{matrix} δ e_{g} \\ δ e_{l} \\ δ X_{n} \end{matrix}], \bar{R} = [\begin{array}{l} r_{1} \\ r_{2} \\ r_{3} \end{array}]

식 (32)를 계산하면, 새로운 시간( $t^{n + 1} = t^{n} + δ t$ )에서의 에너지 및 비응축성 기체 질량 분율은 아래와 같이 계산된다.

(33)

e_{g}^{n + 1} = e_{g}^{n} + δ e_{g}, e_{l}^{n + 1} = e_{l}^{n} + δ e_{l}, X_{n}^{n + 1} = X_{n}^{n + 1} + δ X_{n}

마지막 단계로, 온도와 밀도 및 기타 물성 등의 종속 변수들은 상태 방정식을 이용하여 계산한다.

2.3 내재적(implicit) 수치해법

CFL 수에 의하여 차분 시 구간 크기가 제한되는 반내재적 방법의 계산 시간을 단축하기 위하여 내재적 수치해석 방법이 개발되었다. 2유체 모델은 본질적으로 기체와 액체의 “상-연계(phase link)”에 대해 내재적 계산이 필수적이다. 예를 들어, 기체와 액체의 속도는 계면 운동량 전달을 통해 서로 영향을 미치므로 동시에 내재적으로 계산되어야 한다. 이러한 상-연결의 내재적 계산은 2.2절에서 설명한 반내재적 계산에서도 필요하다. 본 연구에서는 상-연계의 내재적 계산뿐만 아니라, 지배 방정식의 대류 및 확산항도 내재적으로 계산한다. 이는 모든 격자 지점의 값을 동시에 계산한다는 점에서 “공간-연계(space link)” 계산이라 할 수 있다. 식 (16)에 제시된 운동량 방정식의 대류항 및 확산항에 내재적 이산화를 적용하면, 아래와 같은 방정식을 얻을 수 있다.

(34)

[\begin{matrix} m_{g g} m_{g l} \\ m_{l g} m_{l l} \end{matrix}] [\begin{array}{l} {\vec{u}}_{g}^{n + 1} - {\vec{u}}_{g}^{n} \\ {\vec{u}}_{l}^{n + 1} - {\vec{u}}_{l}^{n} \end{array}] = - [\begin{matrix} V α_{g} \\ V α_{l} \end{matrix}] \nabla P_{i}^{n + 1} + [\begin{array}{l} - conv ({\vec{u}}_{g}^{n + 1}) + diff ({\vec{u}}_{g}^{n + 1}) \\ - conv ({\vec{u}}_{l}^{n + 1}) + diff ({\vec{u}}_{l}^{n + 1}) \end{array}] + [\begin{array}{l} \vec{s i_{g}} \\ \vec{s i_{l}} \end{array}]

위 식에서, $conv [{\vec{u}}_{k}^{n + 1}] = \sum_{f = 1}^{n_{neigh}} {(α_{k} ρ_{k})}_{f} ({\vec{u}}_{k, f}^{n + 1} - {\vec{u}}_{k, i}^{n + 1}) ψ_{k, f}^{n + 1}$ , $diff [{\vec{u}}_{k}^{n + 1}] = \sum_{f = 1}^{n_{neigh}} {(α_{k} ρ_{k})}_{f} \frac{{\vec{u}}_{k, j}^{n + 1} - {\vec{u}}_{k, i}^{n + 1}}{\vec{d r_{i j} ∙} {\vec{n}}_{f}} S_{f}$ 로 정의된다. 식 (34)에서, 행렬 요소 $(m_{g g}, m_{g l}, m_{l g}, m_{l l})$ 은 식 (17)에서 정의된 것과 같으며, 소스 벡터인 $({\vec{s i}}_{g}, {\vec{s i}}_{l})$ 은 식 (18)에서 대류항과 확산항을 제외한 값으로 얻어진다. 만약 상-연계와 공간-연계를 동시에 계산하게 된다면, 식 (34)의 계산 비용은 매우 커지게 되므로, 계산을 두 단계로 분리하여 수행한다.

첫 번째 단계에서는 상-연계에 대한 내재적 계산을 수행한다. 식 (34)의 운동량 방정식에서 상-연계에 해당하는 항은 계면 항력과 가상 질량력항 이며, 이들에 대한 내재적 표현은 아래와 같다.

(35)

[\begin{matrix} m_{g g} m_{g l} \\ m_{l g} m_{l l} \end{matrix}] [\begin{array}{lr} {\vec{u}}_{g}^{(1)} - {\vec{u}}_{g}^{n} \\ {\vec{u}}_{l}^{(1)} - {\vec{u}}_{l}^{n} \end{array}] = - [\begin{matrix} V α_{g} \\ V α_{l} \end{matrix}] \nabla P_{i}^{n} + [\begin{array}{l} {\vec{s i}}_{g} \\ {\vec{s i}}_{l} \end{array}]

첫 번째 단계 계산 후, $k$ -상 속도는 ${\vec{u}}_{k}^{(1)}$ 로 저장된다.

두 번째 단계에서는 대류항과 확산항을 내재적으로 계산하며, 업데이트된 속도 벡터는 ${\vec{u}}_{k}^{(2)}$ 로 저장된다.

(36)

m_{g g} [{\vec{u}}_{g}^{(2)} - {\vec{u}}_{g}^{(1)}] = - conv ({\vec{u}}_{g}^{(2)}) + diff ({\vec{u}}_{g}^{(2)}) m_{l l} [{\vec{u}}_{l}^{(2)} - {\vec{u}}_{l}^{(1)}] = - conv ({\vec{u}}_{l}^{(2)}) + diff ({\vec{u}}_{l}^{(2)})

식 (36)의 계산이 완료된 후, 식 (28)과 같은 방식으로 압력 방정식을 계산하면 새로운 시간의 속도 벡터는 다음과 같이 표현될 수 있다.

(37)

[\begin{matrix} m_{g g} m_{g l} \\ m_{l g} m_{l l} \end{matrix}] [\begin{matrix} {\vec{u}}_{g}^{n + 1} - {\vec{u}}_{g}^{(2)} \\ {\vec{u}}_{l}^{n + 1} - {\vec{u}}_{l}^{(2)} \end{matrix}] = - [\begin{matrix} V α_{g} \\ V α_{l} \end{matrix}] \nabla δ P_{i}^{n}

식 (35), (36), (37)을 합하면, 두 계산 단계를 통한 내재적 운동량 방정식은 다음과 같은 식으로 도출된다.

(38)

[\begin{matrix} m_{g g} m_{g l} \\ m_{l g} m_{l l} \end{matrix}] [\begin{matrix} {\vec{u}}_{g}^{n + 1} - {\vec{u}}_{g}^{n} \\ {\vec{u}}_{l}^{n + 1} - {\vec{u}}_{l}^{n} \end{matrix}] = - [\begin{matrix} V α_{g} \\ V α_{l} \end{matrix}] \nabla P_{i}^{n + 1} + [\begin{array}{l} - conv ({\vec{u}}_{g}^{(2)}) + diff ({\vec{u}}_{g}^{(2)}) \\ - conv ({\vec{u}}_{l}^{(2)}) + diff ({\vec{u}}_{l}^{(2)}) \end{array}] + [\begin{array}{l} {\vec{s i}}_{g} \\ {\vec{s i}}_{l} \end{array}]

식 (38), (34)를 비교하면, 계산 단계를 분리함으로써 생기는 오차는 $|{\vec{u}}_{k}^{n + 1} - {\vec{u}}_{k}^{(2)}|$ 에 비례한다. 이 오차는 상대적으로 작고, 일반적으로 약 2회의 반복 계산(iteration)을 통해 수렴된 수치 해를 얻을 수 있다.

운동량 방정식에서와 같은 방법으로, 에너지 및 비응축성 기체 방정식(식 (13), (4))에 내재적 이산화를 적용하면 다음과 같은 방정식을 얻을 수 있다.

(39)

[\begin{array}{l} o_{11} & o_{12} & o_{13} \\ o_{21} & o_{22} & o_{23} \\ o_{31} & o_{32} & o_{33} \end{array}] [\begin{matrix} e_{g}^{n + 1} - e_{g}^{n} \\ e_{l}^{n + 1} - e_{l}^{n} \\ X_{n}^{n + 1} - X_{n}^{n} \end{matrix}] = [\begin{matrix} - conv [e_{g}^{n + 1}] + diff [e_{g}^{n + 1}] \\ - conv [e_{l}^{n + 1}] + diff [e_{l}^{n + 1}] \\ - conv [X_{n}^{n + 1}] + conv [X_{n}^{n + 1}] \end{matrix}] + [\begin{array}{l} r i_{1} \\ r i_{2} \\ r i_{3} \end{array}]

여기서, $conv [e_{k}^{n + 1}] = \frac{δ t}{V_{i}} \sum_{f = 1}^{n_{neigh}} {(α_{k} ρ_{k})}_{f} (e_{k, f}^{n + 1} - e_{k, i}^{n + 1}) ψ_{k, f}^{n + 1}$ , $conv [X_{n}^{n + 1}] = \frac{δ t}{V_{i}} \sum_{f = 1}^{n_{neigh}} {(α_{g} ρ_{g})}_{f} (X_{n, f}^{n + 1} - X_{n, i}^{n + 1}) ψ_{g, f}^{n + 1}$ , $diff [e_{k, f}^{n + 1}] = \frac{δ t}{V_{i}} \sum_{f = 1}^{n_{neigh}} {(α_{k} k_{k})}_{f} \{e_{k, f}^{n} + {(\frac{\partial T_{k}}{\partial e_{k}})}_{j} δ e_{k, j} - e_{k, i}^{n} - {(\frac{\partial T_{k}}{\partial e_{k}})}_{i} δ e_{k, i}\} \frac{S_{f}}{\vec{d r_{i j}} ∙ {\vec{n}}_{f}}$ 로 정의 된다. 식 (39)에서 행렬 요소 $(o_{11}, . . ., o_{33})$ 는 식 (34)에서 정의된 것과 같으며, 소스 벡터 $(r i_{1}, r i_{2}, r i_{3})$ 는 식 (32)의 $(r_{1}, r_{2}, r_{3})$ 에서 대류항과 확산항을 제외한 값으로 얻어진다. 운동량 방정식에서의 계산과 같은 방식으로, 식 (39)의 계산도 공간-연계와 상-연계 단계로 분리하여 수행한다. 첫 번째 단계는 내재적 공간-연계 계산에 해당한다.

(40)

o_{11} [e_{g}^{(1)} - e_{g}^{n}] = - conv [e_{g}^{(1)}] + diff [e_{g}^{(1)}] o_{22} [e_{l}^{(1)} - e_{l}^{n}] = - conv [e_{l}^{(1)}] + diff [e_{g}^{(1)}] o_{33} [X_{n}^{(1)} - X_{n}^{n}] = - conv [X_{n}^{(1)}] + diff [X_{n}^{(1)}]

여기서 $e_{k}^{(1)}$ 는 첫 번째 단계 이후의 $k$ -상 내부 내부에너지를 나타낸다. 식 (40)을 계산한 후, 다음 단계로 내재적 상-연계 계산을 수행하는데, 이때 아래와 같이 상-계면에서의 내재적 열전달 계산을 고려한다.

(41)

[\begin{array}{l} o_{11} & o_{12} & o_{13} \\ o_{21} & o_{22} & o_{23} \\ o_{31} & o_{32} & o_{33} \end{array}] [\begin{matrix} e_{g}^{(2)} - e_{g}^{n} \\ e_{l}^{(2)} - e_{l}^{n} \\ X_{n}^{(2)} - X_{n}^{n} \end{matrix}] = [\begin{matrix} o_{11} (e_{g}^{(1)} - e_{g}^{n}) + r i_{1} \\ o_{22} (e_{l}^{(1)} - e_{l}^{n}) + r i_{2} \\ o_{33} (X_{n}^{(1)} - X_{n}^{n}) + r i_{3} \end{matrix}]

식 (40)와 (41)을 합하면, 두 계산 단계의 내재적 에너지 계산은 다음과 같다.

(42)

[\begin{array}{l} o_{11} & o_{12} & o_{13} \\ o_{21} & o_{22} & o_{23} \\ o_{31} & o_{32} & o_{33} \end{array}] [\begin{matrix} e_{g}^{(2)} - e_{g}^{n} \\ e_{l}^{(2)} - e_{l}^{n} \\ X_{n}^{(2)} - X_{n}^{n} \end{matrix}] = [\begin{matrix} - conv [e_{g}^{(1)}] + diff [e_{g}^{(1)}] \\ - conv [e_{l}^{(1)}] + diff [e_{l}^{(1)}] \\ - conv [X_{n}^{(1)}] \end{matrix}] + [\begin{array}{l} r i_{1} \\ r i_{2} \\ r i_{3} \end{array}]

내재적 운동방정식 계산에서와 유사하게, 계산 단계를 분리함으로써 생기는 오차는 $|e_{k}^{(2)} - e_{k}^{(1)}|$ 와 $|X_{n}^{(2)} - X_{n}^{(1)}|$ 에 비례하며, 약 2회의 반복 계산을 통해 최소화할 수 있다.

3. 냉각재상실사고(LOCA) 개념 문제 확인(Verification) 해석

제안된 내재적 방법의 수치 안정성 및 정확성을 확인하기 위한 해석을 수행하였다. 본 확인 해석의 주요 관점은 PWR의 가상사고 조건에서 발생하는 급격한 상변화를 동반하는 문제에 대해, 제안된 내재적 수치해석 방법의 수치적 안정성 및 정확성을 평가하는 데 있다. 이 해석 문제는 Fig. 1에 나타난 바와 같이 총 1,250개(25×50)의 계산 셀이 사용되는 2차원 개념적 문제이며, 초기 조건으로 왼쪽 하부에 20 MW의 균일한 열생성원과 오른쪽 상부에 같은 크기의 열 흡수원(sink)을 설정하여 단상 액체(single phase liquid)의 자연 순환 유동을 형성한다. 단상 액체인 물은 150기압에서 300°C의 과냉각(subcooled) 상태이다. 이후, 300초 시점에 Fig. 2와 같이 왼쪽 상부 경계에서 액체 블로우다운(blowdown) 유동 경계 조건이 적용되고, 이로 인해 급격한 압력 강하에 따른 플래싱(flashing)현상이 발생하여 대량의 증기가 생성된다. 플래싱으로 인해 수위는 계속해서 감소하다가, 시스템 압력이 10 MPa 이하로 떨어지면 오른쪽 상단에서 100°C의 냉각수가 주입되기 시작한다. 이는 PWR의 안전 주입(safey injection) 유동을 모사한 것으로, 이 단계에서는 과냉각수 주입에 따른 증기 응축(condensation)을 동반한 수위 증가가 이루어지며, 약 2,000초 부근에서 초기 상태의 수위를 회복한다.

https://cdn.apub.kr/journalsite/sites/kscfe/2025-030-03/N0500300307/images/jkscfe_2025_303_102_F2.jpg

Fig. 2.

Schematic of pressurized water blowdown and refill process

이 개념적 문제를 대상으로 CUPID의 내재적 수치해석 방법을 적용한 해석을 수행하였으며, 그 결과를 반내재적 수치해석 방법의 계산 결과와 비교하였다. Fig. 3은 내재적 방법을 적용하여 CFL 수를 0.1에서 5로 증가시키며 계산한 수위 변화 양상을 반내재적 방법을 적용한 계산 결과와 비교한 것이다. 초기 300초까지 물의 수위가 50 m³로 유지되다가, 이후 왼쪽 상부에서 블로우다운(blowdown)이 시작되면 압력 강하로 인한 플래싱으로 급격한 기화가 발생하고, 약 846초에 수위가 최솟값인 36 m³에 도달하게 되며, 이후 오른쪽 상부에서 냉각수 주입이 시작되면 수위는 점차 회복된다. 내재적 방법을 적용한 계산 결과는 모든 CFL 수에 대하여 반내재적 방법을 적용한 계산 결과와 매우 유사한 경향을 보인다. CFL은 $u_{\max} \frac{Δ t}{Δ x}$ 에 의하여 정의되는데, $u_{\max}$ 는 전 유동 영역에서의 최대 속도를 적용하였다. CFL=5을 사용한 내재적 방법에서 최소 수위가 다소 낮게 나타났으나 그 차이는 크지 않다. 계산 시간 측면에서는, Table 1에 나타난 바와 같이 내재적 방법의 경우(CFL=5) 67초이며, 반내재적 방법(CFL=0.1)은 약 1,035초가 소요되었다. 결과적으로, 내재적 방법을 적용하였을 때 계산 정확도를 유지하면서도 약 15배의 계산 속도 향상이 달성되었음을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/kscfe/2025-030-03/N0500300307/images/jkscfe_2025_303_102_F3.jpg

Fig. 3.

Water volume transient during blowdown and refill

Table 1.

Comparison of calculation times for the conceptual problem

Numerical scheme	CFL number	Calculation time(s)
Explicit	0.1	1035
Implicit	0.1	1806
Implicit	1.0	209
Implicit	2.0	115
Implicit	3.0	84
Implicit	5.0	67

4. 재관수(Reflood) 실험에 대한 검증(Validation)

내재적 수치해석 방법의 실제 현상 예측에 대한 검증(validation)을 위하여 PWR의 주요 고려 사건 중 하나인 LOCA 발생 시, 핵연료봉 거동과 관련된 중요한 안전성 평가 실험인 재관수(reflood) 실험에 대한 해석을 수행하였다. LOCA 발생 초기 냉각재 손실로 인해 핵연료봉이 과열되며 이후 노심에 냉각수가 재주입되면서 과열된 연료봉이 다시 냉각되는데, 재관수 실험은 냉각수가 재주입 되는 상황에서 가열된 핵연료봉의 냉각 과정을 모의한다.

핵연료봉 실제 길이를 대상으로 수행된 재관수 실험(FLECHT-SEASET)에 대하여 내재적 및 반내재적 방법을 적용한 CUPID 코드 해석을 수행하였다. FLECHT-SEASET 실험은 Westinghouse 17 × 17 연료봉 집합체를 모사하며, 가열 길이는 기존 PWR과 동일한 3.66 m이고, 수로 직경은 0.194 m이다. Fig. 4는 수로의 수평 단면 및 수직 형상을 나타낸다. 총 161개의 가열된 봉(heated rod)과 16개의 가열되지 않은 봉(unheated rod)이 좌측 그림과 같이 배열되어 있다. Fig. 4의 우측에 표시된 화살표는 벽면 온도 측정 지점의 높이를 나타낸다. 검증 대상 실험 조건으로는 31203, 31302, 31701 세 가지가 선택되었고, 이들에 대한 실험 조건은 Table 2에 요약되어 있다. Fig. 5에서는 가열봉에 적용되는 코사인(cosine) 형태 축 방향 출력(power) 분포를 나타내었다. 1차원 및 3차원 해석을 위한 두 종류의 격자가 생성되었으며, 이는 Fig. 6에 나타나 있다. 격자수는 1차원 해석에서 20개, 3차원 해석에서 800개로 구성되어 있다. 해당 격자는 가열봉, 가열되지 않은 봉, 및 내부 구조를 고려한 다공성 매질을 적용하였다. 총 161개의 가열봉의 열전도 계산을 위하여 반경 방향으로 7개의 노드가 있는 2차원 열 구조 해석 모델을 적용하였다. 가열봉은 반경 위치에 따라 보론 나이트라이드(Boron Nitride), 인코넬(Inconel), 모넬 K-600(Monel K-600)의 세 가지 재료로 구성되어 있다.

https://cdn.apub.kr/journalsite/sites/kscfe/2025-030-03/N0500300307/images/jkscfe_2025_303_102_F4.jpg

Fig. 4.

Test channel of the FLECHT-SEASET experiment

Table 2.

Test conditions of the FLECHT-SEASET experiment

Test number	Flooding rate (m/s)	Pressure (MPa)	Initial rod temperature(°C)	Coolant temperature(°C)	Peak power (kW/m)
31203	0.0384	0.28	872	52	2.3
31701	0.1550	0.28	872	53	2.3
31302	0.0765	0.28	869	52	2.3

https://cdn.apub.kr/journalsite/sites/kscfe/2025-030-03/N0500300307/images/jkscfe_2025_303_102_F5.jpg

Fig. 5.

Axial power profile of the FLECHT-SEASET experiment

https://cdn.apub.kr/journalsite/sites/kscfe/2025-030-03/N0500300307/images/jkscfe_2025_303_102_F6.jpg

Fig. 6.

Mesh configuration of 1D and 3D calculations

CUPID를 이용한 가열봉의 벽면 온도 해석 결과는 Fig. 4에서 제시한 바와 같이 하단(61 mm), 중간(1,820 mm), 상단(3,050 mm)의 세 위치에서 실험 결과와 비교하였다. 모든 계산은, 1차원 및 3차원 계산 결과를 동일 높이에서 비교하였다. 3차원 계산 데이터는 1차원 결과와 비교하기 위해 각 높이에서 수로 단면적으로 평균하였다. 내재적 방법을 적용한 계산 시간 단축 효과를 평가하기 위하여, CFL 수를 증가시키면서 계산이 수행하여 그 결과를 비교하였다. 특히, 계산 시간이 훨씬 긴 3차원 계산에는 내재적 방법의 효과를 평가하기 위하여 CFL 수를 1.0 보다 큰 값으로 설정하였다. 모든 계산은 MPI 라이브러리를[17] 사용한 병렬 계산이 적용되었으며, 4개의 프로세서 적용을 위하여 계산 격자에 대하여 4개의 계산 도메인을 설정하였다.

과도현상 초기에는 유동 채널이 증기로 가득 차 있어, 가열봉 내부의 열원으로 인하여 벽면 온도가 상승한다. 연료봉 피복관의 최대 온도(PCT, Peak Cladding Temperature)는 안전 한계 온도(1,477 K)를 초과 시 연료봉 용융을 초래할 수 있어, 안전 평가를 위한 주요 지표로 간주 된다. PCT는 축 방향 발열량이 최대인 중간 높이(1,820 mm)에서 나타났다. 이후 채널 하부로 냉각수가 주입되면서 벽면 온도가 감소하기 시작한다. 수면이 연료봉 가열 표면을 덮는 시점부터 벽면 온도는 약 400 K 수준으로 급격히 하강한다. 급랭 전선(quench front)이라 불리는 수위는 시간이 지남에 따라서 상단으로 이동하는데, 특정 높이에서 수면이 연료봉 가열 표면을 덮는 시점을 급랭 시간(quenching time)이라 정의한다.

Figs. 7, 8, 9에는 실험 조건 31203, 31302, 31701에 대한 1차원 계산 결과를 각각 비교하여 나타내었는데, 왼쪽 그래프는 CFL 0.1에서 내재적 방법과 반내재적 방법에 의한 계산 결과를 실험 데이터와 비교한 것이다. 내재적 계산 방법의 계산 결과는 반내재적 계산 방법의 계산 결과와 거의 일치하였고, 실험에서 측정된 피복관 최대 온도(PCT) 또한 잘 예측되었다. 다만, 상단 위치에서의 급랭 시간은 계산 결과가 실험값에 비하여 다소 빠르게 예측되었다. 우측 그래프는 내재적 방법을 적용하였을 때, CFL 수를 0.1에서 1.0까지 증가시키면서 계산한 결과를 비교한 것이다. 내재적 계산에 사용된 CFL 수를 반내재적 계산에서 사용된 값(0.1)의 약 5~10배까지 증가시킨 경우의 계산 결과가 거의 일치하는 것으로 나타났다. 특히, 해석의 주요 관심 인자인 PCT 예측은 모두 일치하는 것으로 나타났다.

https://cdn.apub.kr/journalsite/sites/kscfe/2025-030-03/N0500300307/images/jkscfe_2025_303_102_F7.jpg

Fig. 7.

Comparison of 1D calculation results for case 31203

https://cdn.apub.kr/journalsite/sites/kscfe/2025-030-03/N0500300307/images/jkscfe_2025_303_102_F8.jpg

Fig. 8.

Comparison of 1D calculation results for case 31302

https://cdn.apub.kr/journalsite/sites/kscfe/2025-030-03/N0500300307/images/jkscfe_2025_303_102_F9.jpg

Fig. 9.

Comparison of 1D calculation results for case 31701

Figs. 10, 11, 12에는 실험 조건 31203, 31302, 31701에 대한 3차원 계산 결과를 비교하여 나타내었다. 3차원 계산은 1차원 계산보다 소요 시간이 훨씬 오래 걸리기 때문에 더 큰 CFL 수가 사용되었다. 해당 실험은 1차원 현상에 중점을 두고 설계되었기 때문에, 3차원 계산 결과도 1차원 결과와 유사하게 나타났다. 좌측 그래프에서는 CFL 수가 약 0.6일 때 내재적 방법 및 반내재적 방법의 계산 결과가 거의 일치하였으며, 실험에서 측정된 PCT도 정확히 예측되었다. 우측 그래프에서는 다양한 CFL 수(허용 가능한 최댓값은 6.0~7.2)를 적용한 계산 결과를 비교하였다. 1차원 계산과 마찬가지로, PCT 예측에서의 차이는 매우 작아 무시할 수 있는 수준이었다.

https://cdn.apub.kr/journalsite/sites/kscfe/2025-030-03/N0500300307/images/jkscfe_2025_303_102_F10.jpg

Fig. 10.

Comparison of 3D calculation results for case 31203

https://cdn.apub.kr/journalsite/sites/kscfe/2025-030-03/N0500300307/images/jkscfe_2025_303_102_F11.jpg

Fig. 11.

Comparison of 3D calculation results for case 31302

https://cdn.apub.kr/journalsite/sites/kscfe/2025-030-03/N0500300307/images/jkscfe_2025_303_102_F12.jpg

Fig. 12.

Comparison of 3D calculation results for case 31701

실험 조건 31203, 31302, 31701에 대한 내재적 방법과 반내재적 방법의 계산 시간을 Tables 3, 4, 5에 정리하여 비교하였다. 동일한 CFL 수를 사용할 경우, 내재적 방법의 계산 시간은 반내재적 방법의 계산 시간보다 더 오래 소요된다. 이는 내재적 방법을 적용할 때 대류 및 확산 항의 내재적 계산을 위하여 추가적인 계산이 필요하기 때문이다. 내재적 방법을 적용한 경우, 계산 시간은 CFL 수의 역수에 비례하여 감소하였는데, 계산 시간은 반내재적 방법을 적용한 계산 대비 약 6.0배에서 8.1배 단축되었다.

Table 3.

Comparison of calculation times for case 31203

Numerical scheme	1D		3D
Numerical scheme	CFL number	Calculation time(s)	CFL number	Calculation time(s)
Explicit	0.1	63	0.8	2166
Implicit	0.1	102	0.8	2316
Implicit	0.5	27	3.7	595
Implicit	0.8	13	7.2	342

Table 4.

Comparison of calculation times for case 31302

Numerical scheme	1D		3D
Numerical scheme	CFL number	Calculation time(s)	CFL number	Calculation time(s)
Explicit	0.1	125	0.6	1505
Implicit	0.1	202	0.6	1715
Implicit	0.4	47	3.2	448
Implicit	1.0	26	6.1	252

Table 5.

Comparison of calculation times for case 31701

Numerical scheme	1D		3D
Numerical scheme	CFL number	Calculation time(s)	CFL number	Calculation time(s)
Explicit	0.1	32	0.6	877
Implicit	0.1	58	0.6	910
Implicit	0.5	13	3.0	209
Implicit	0.9	10	6.0	132

5. 결 론

본 연구에서는 대규모 상변화가 수반되는 PWR의 2상 유동 해석을 위한 내재적 수치해석 방법을 개발하였다. 상변화 모델은 운동량 방정식에서 유도된 압력 방정식과 연계하여 내재적으로 계산되었다. 동시에, 지배 방정식의 대류항과 확산항도 내재적으로 계산되었다. 시스템 행렬의 크기를 줄이기 위해 계산 절차는 “상-연계(phase-link)”와 “공간-연계(space-link)”의 두 단계로 분리되었다.

개발된 내재적 수치해석 방법은 CUPID 코드에서 구현되었으며, LOCA 중 발생하는 2상 유동 블로우다운(blowdown)과 재급수(refill)현상을 모사한 개념적 문제를 통해 확인 해석을 수행하였다. 블로우다운과 재급수로 인해 대규모 상변화가 발생하는 급격한 과도 현상은 내재적 방법을 적용하여 성공적으로 계산되었으며, 계산 결과는 반내재적 방법과 동일하였다. 내재적 방법을 적용한 경우, 큰 CFL 수를 적용함으로써 계산 결과는 거의 동일한 수준을 유지하면서 약 15배의 계산 시간 단축 효과를 얻을 수 있었다.

또한 개발된 내재적 수치해석 방법의 검증을 위해, LOCA의 재급수(refill) 현상을 모사한 FLECHT-SEASET 실험에 대하여 내재적 방법을 적용하여 해석을 수행하였다. 내재적 방법과 반내재적 방법의 계산 결과는 거의 일치하였으며, 실험에서 측정된 연료봉 피복관 최대 온도(PCT)도 잘 예측되었다. 내재적 방법에 의한 계산의 경우, 3차원 계산에서 최대 7.2까지의 큰 CFL 수가 허용되었으며, 그 결과 계산 시간은 반내재적 방법의 계산 대비 약 6.0배에서 8.1배 단축되었다.

제안된 내재적 수치해석 방법의 정확성은 일련의 테스트 문제에서 반내재적 수치해석 방법과 동일한 수치 해를 도출함으로써 입증되었다. 내재적 방법의 경우 큰 CFL 수를 적용한 계산이 안정적으로 계산되었고, 작은 CFL 수를 사용할 때의 결과와의 차이가 미미하여 수치해석의 강건성과 효율성이 확인되었다. 따라서 본 연구에서 개발된 내재적 수치해석 방법은 원자로 해석을 위한 2유체 모델의 신속한 계산에 매우 유용하게 활용될 수 있다.

Acknowledgements

본 연구는 한국전력국제원자력대학원대학교의 연구비 지원(2025~2026)에 의해 수행되었다.

Appendix

The coefficients in Eq. (43) are as follows.

(43)

o_{11} = α_{g}^{n} ρ_{g}^{n} - \frac{α_{g}^{n} P^{n}}{ρ_{g}^{n}} \frac{\partial ρ_{g}}{\partial e_{g}} - δ t H_{i g} (1 - \frac{Δ h_{g}^{*}}{h_{g l}^{*}}) (\frac{\partial T^{s}}{\partial e_{g}} - \frac{\partial T_{g}}{\partial e_{g}}) + δ t H_{g l} \frac{\partial T_{g}}{\partial e_{g}} + δ t \frac{Δ h_{g}^{*}}{h_{g l}^{*}} H_{i l} \frac{\partial T^{s}}{\partial e_{g}} o_{12} = - δ t H_{g l} \frac{\partial T_{l}}{\partial e_{l}} - δ t \frac{Δ h_{g}^{*}}{h_{g l}^{*}} H_{i l} \frac{\partial T_{l}}{\partial e_{l}} o_{13} = - \frac{α_{g}^{n} P^{n}}{ρ_{g}^{n}} \frac{\partial ρ_{g}}{\partial X_{n}} - δ t H_{i g} (1 - \frac{Δ h_{g}^{*}}{h_{g l}^{*}}) (\frac{\partial T^{s}}{\partial X_{n}} - \frac{\partial T_{g}}{\partial X_{n}}) + δ t H_{g l} \frac{\partial T_{g}}{\partial X_{n}} + δ t \frac{Δ h_{g}^{*}}{h_{g l}^{*}} H_{i l} \frac{\partial T^{s}}{\partial X_{n}} r_{1} = - \frac{δ t}{V_{i}} \sum_{f = 1}^{N_{n e i g h}} \{[{(α_{g} ρ_{g} e_{g})}_{f}^{n} - e_{g, i}^{n} {(α_{g} ρ_{g})}_{f}^{n}] ψ_{g, f}^{n + 1} - [{(α_{g} ρ_{g})}_{f}^{n} - ρ_{g, i}^{n} α_{g, f}^{n}] ψ_{g, f}^{n + 1} - {(α_{g} {\vec{q}}_{g})}_{f}^{n} ∙ {\vec{S}}_{f}\} + δ t H_{i g} (1 - \frac{Δ h_{g}^{*}}{h_{g l}^{*}}) (T^{s, n + 1} - T_{g}^{n + 1}) - δ t H_{g l} (T_{g}^{n + 1} - T_{l}^{n + 1}) - δ t \frac{Δ h_{g}^{*}}{h_{g l}^{*}} H_{i l} (T^{s, n + 1} - T_{l}^{n + 1}) o_{21} = - δ t H_{g l} \frac{\partial T_{g}}{\partial e_{g}} - δ t H_{i l} (1 + \frac{Δ h_{l}^{*}}{h_{g l}^{*}}) \frac{\partial T^{s}}{\partial e_{g}} - δ t \frac{Δ h_{l}^{*}}{h_{g l}^{*}} H_{i g} (\frac{\partial T^{s}}{\partial e_{g}} - \frac{\partial T_{g}}{\partial e_{g}}) o_{22} = α_{l}^{n} ρ_{l}^{n} - \frac{α_{l}^{n} P^{n}}{ρ_{l}^{n}} \frac{\partial ρ_{l}}{\partial e_{l}} + δ t H_{i l} \frac{\partial T_{l}}{\partial e_{l}} + δ t H_{g l} \frac{\partial T_{l}}{\partial e_{l}} o_{23} = - δ t H_{i l} (1 + \frac{Δ h_{l}^{*}}{h_{g l}^{*}}) \frac{\partial T^{s}}{\partial X_{n}} - δ t H_{g l} \frac{\partial T_{g}}{\partial X_{n}} - δ t \frac{Δ h_{l}^{*}}{h_{g l}^{*}} H_{i g} (\frac{\partial T^{s}}{\partial X_{n}} - \frac{\partial T_{g}}{\partial X_{n}}) r_{2} = - \frac{δ t}{V_{i}} \sum_{f = 1}^{N_{n e i g h}} \{[{(α_{l} ρ_{l} e_{l})}_{f}^{n} - e_{l, i}^{n} {(α_{l} ρ_{l})}_{f}^{n}] ψ_{l, f}^{n + 1} - [{(α_{l} ρ_{l})}_{f}^{n} - ρ_{l, i}^{n} α_{l, f}^{n}] ψ_{l, f}^{n + 1} - {(α_{l} {\vec{q}}_{l})}_{f}^{n} ∙ {\vec{S}}_{f}\} + δ t H_{i l} (1 + \frac{Δ h_{l}^{*}}{h_{g l}^{*}}) (T^{s, n + 1} - T_{l}^{n + 1}) + δ t H_{g l} (T_{g}^{n + 1} - T_{l}^{n + 1}) + δ t \frac{Δ h_{l}^{*}}{h_{g l}^{*}} H_{i g} (T^{s, n + 1} - T_{g}^{n + 1}) o_{31} = - \frac{δ t X_{n}}{h_{g l}^{*}} [H_{i g} (\frac{\partial T^{s}}{\partial e_{g}} - \frac{\partial T_{g}}{\partial e_{g}}) + H_{i l} \frac{\partial T^{s}}{\partial e_{g}}] o_{32} = \frac{δ t X_{n}}{h_{g l}^{*}} H_{i l} \frac{\partial T_{l}}{\partial e_{l}} o_{33} = α_{g}^{n} ρ_{g}^{n} - \frac{δ t X_{n}}{h_{g l}^{*}} [H_{i g} (\frac{\partial T^{s}}{\partial X_{n}} - \frac{\partial T_{g}}{\partial X_{n}}) + H_{i l} \frac{\partial T^{s}}{\partial X_{n}}] r_{3} = - \frac{δ t}{V_{i}} \sum_{f = 1}^{N_{neigh}} [{(α_{g} ρ_{g} X_{n})}_{f}^{n} - X_{n, i}^{n} {(α_{g} ρ_{g})}_{f}^{n}] ψ_{g, f}^{n + 1} + \frac{δ t X_{n}}{h_{g l}^{*}} [H_{i g} (T^{s, n + 1} - T_{g}^{n + 1}) + H_{i l} (T^{s, n + 1} - T_{l}^{n + 1})] > where Δ h_{g}^{*} = h_{g}^{*} - e_{g} - \frac{P}{ρ_{g}}, Δ h_{l}^{*} = h_{l}^{*} - e_{l} - \frac{P}{ρ_{l}}, and h_{g l}^{*} = h_{g}^{*} - h_{l}^{*} .

References

2018, United States Nuclear Regulatory Commission, “APR1400 Final Safety Evaluation,” Available online: https://www.nrc.gov/reactors/new-reactors/large-lwr/design-cert/apr1400/ser-final.html.

1982, Ishii, M. and Mishima, K., “Two-fluid model and hydrodynamic constitutive relations,” Nucl. Eng. Des., Vol.82, pp.107-126.

10.1016/0029-5493(84)90207-3

2001, United States Nuclear Regulatory Commission, “RELAP5/MOD3.3 Code Manuals,” Idaho National Laboratory.

2003, Robert, M., Farvacque, M., Parent, M. and Faydide, B., “CATHARE 2 V2.5: A Fully Validated CATHARE Version for Various Applications,” Proc. 10th Int. Topl. Mtg. Nuclear Reactor Thermal Hydraulics (NURETH 10), Seoul, Korea, Oct. 5-9.

1999, Jeong, J.J., Ha, K.S., Chung, B.D. and Lee, W.J., “Development of a multi-dimensional thermal-hydraulic system code, MARS 1.3.1,” Ann. Nucl. Energy, Vol.26, No.18, pp.1611-1642.

10.1016/S0306-4549(99)00039-0

2005, The RELAP5-3D Code Development Team, “RELAP5-3D Code Manula Volume I: Code Structure, System Models and Solution Method,” Idaho National Laboratory.

2011, Emonot, P., Souyri, A., Gandrille, J.L. and Barré, F., “CATHARE-3: A new system code for thermal-hydraulics in the context of the NEPTUNE project,” Nucl. Eng. Des., Vol.241, No.11, pp.4476-4481.

10.1016/j.nucengdes.2011.04.049

1999, Jeong, J.-J., Ha, K.S., Chung, B.D. and Lee, W.J., “Development of a multi-dimensional thermal-hydraulic system code, MARS 1.3.1,” Ann. Nucl. Energy, Vol.26, No.18, pp.1611-1642.

10.1016/S0306-4549(99)00039-0

2025, ANSYS, Inc, “Ansys CFX,” Available online: https://www.ansys.com/products/fluids/ansys-cfx.

2025, ANSYS, Inc, “Ansys Fluent,” Available online: https://www.ansys.com/products/fluids/ansys-fluent.

2025, SIEMENS Digital Industries Software, Inc, “Simcenter STAR-CCM+ CFD software,” Available online: https://plm.sw.siemens.com/en-US/simcenter/fluids-thermal-simulation/star-ccm/.

2007, Guelfi, A., Bestion, D., Boucker, M., Boudier, P., Fillion, P., Grandotto, M., Hérard, J.-M., Hervieu, E. and Péturaud, P., “NEPTUNE: a new software platform for advanced nuclear thermal hydraulics,” Nucl. Sci. Eng., Vol.156, No.3, pp.281-324.

10.13182/NSE05-98

2010, Jeong, J.J., Yoon, H.Y., Park, I.K. and Cho, H.K., “The CUPID code development and assessment strategy,” Nucl. Eng. Technol., Vol.42, No.6, pp.636-655.

10.5516/NET.2010.42.6.636

2014, Yoon, H.Y., Lee, J.R., Kim, H., Park, I.K., Song, C.H., Cho, H.K. and Jeong, J.J., “Recent improvements in the CUPID code for a multi-dimensional two-phase flow analysis of nuclear reactor components,” Nucl. Eng. Technol, Vol.46, No.5, p.655.

10.5516/NET.02.2014.023

2020, Yoon, H.Y., Park, I.K., Lee, J.R., Lee, S.J., Cho, Y.J., Do, S.J., Cho, H.K. and Jeong, J.J., “A multiscale and multiphysics PWR safety analysis at a subchannel scale,” Nucl. Sci. Eng., Vol.194, pp.633-649.

10.1080/00295639.2020.1727698

1984, Hochreiter, L.E., Rupprecht, S.D., Dederer, J.T., Wong, S., Rosal, E.R., Sinwell, B.R., Kovdalski, F.J. and Quaglia, R., “PWR FLECHT SEASET Systems and reflux condensation data evaluation and analysis report,” NRC/EPRI/Westinghouse Report No. 14, NUREG/CR-3654.

2012, Message Passing Interface Forum, “MPI: a message passing interface standard,” High-Performance Computing Center, University of Stuttgart.

1992, Henk A. and van der Vorst, “Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems,” SIAM J. Sci. Comput., Vol.13, No.2, pp.631-644.

10.1137/0913035

Journal of Computational Fluids Engineering ISSN:1598-6071(Print) 3022-6252(Online) 한국전산유체공학회지

Preview

DEVELOPMENT OF AN IMPLICIT NUMERICAL METHOD FOR TWO-PHASE FLOW ANALYSIS OF PRESSURIZED WATER REACTORS

ABSTRACT

MAIN

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

Fig. 1.

Unstructured mesh and geometric parameters used for discretization

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

(34)

(35)

(36)

(37)

(38)

(39)

(40)

(41)

(42)

Fig. 2.

Schematic of pressurized water blowdown and refill process

Fig. 3.

Water volume transient during blowdown and refill

Table 1.

Comparison of calculation times for the conceptual problem

Fig. 4.

Test channel of the FLECHT-SEASET experiment

Table 2.

Test conditions of the FLECHT-SEASET experiment

Fig. 5.

Axial power profile of the FLECHT-SEASET experiment

Fig. 6.

Mesh configuration of 1D and 3D calculations

Fig. 7.

Comparison of 1D calculation results for case 31203

Fig. 8.

Comparison of 1D calculation results for case 31302

Fig. 9.

Comparison of 1D calculation results for case 31701

Fig. 10.

Comparison of 3D calculation results for case 31203

Fig. 11.

Comparison of 3D calculation results for case 31302

Fig. 12.

Comparison of 3D calculation results for case 31701

Table 3.

Comparison of calculation times for case 31203

Table 4.

Comparison of calculation times for case 31302

Table 5.