SPH-based Multiphysics Simulation Method for Thermo-Fluid-Structure Interaction in Severe Accidents = 중대사고 조건의 열·유체·구조 상호작용 해석을 위한 SPH 기반 다물리 시뮬레이션 기법 개발

국문 초록 (Abstract)

원자력발전소의 중대사고란 노심 손상이 발생하고, 그 결과로 방사성 물질이 격납계통 밖으로 방출될 가능성이 있는 사고 시나리오를 의미한다. 이러한 사고는 다중 안전계통에 의해 발생 확률이 매우 낮게 유지되지만, 일단 발생할 경우 원전의 안전성은 물론 사회적 신뢰에 심각한 영향을 미친다. 실제로 스리마일섬(1979), 체르노빌(1986), 후쿠시마(2011) 사고는 중대사고가 원자력 산업 전체에 미치는 기술적·사회적 파급효과를 여실히 보여준 사례이다.

따라서 중대사고를 정량적으로 예측하고 정확하게 해석하는 것은 원전의 안전성 향상을 위한 핵심 요소이다. 이는 사고 시 핵연료와 구조물의 열적·기계적 거동을 파악하고, 효과적인 완화 전략을 수립하며, 격납계통의 구조 건전성을 확보하고, 방사성 물질의 환경 유출을 최소화하는 데 필수적인 기반을 제공한다.

중대사고 해석은 고온·고점성 유체, 자유표면 유동, 고체-액체 상변화, 탄소성·크리프·파손 거동, 열·유체·구조 상호작용 등 다양한 물리현상이 복합적으로 작용하는 매우 복잡한 문제이다. 이러한 다물리적 현상들은 상호 밀접히 연결되어 있으므로, 신뢰성 있는 안전성 평가를 위해서는 하나의 통합된 해석 틀 내에서 열·유체·구조 상호작용 현상를 정밀하게 예측할 필요가 있다. 이러한 다물리 해석을 위한 수치기법은 정확성, 안정성, 계산 효율성, 구현 용이성을 갖추어야 한다. 입자기반 무격자 해석기법인 완화입자동역학(SPH)는 자유표면 유동, 상변화, 다물체 충돌, 구조물 손상 및 파괴 해석에 유리한 장점을 제공한다. 이에 본 연구에서는 열·유체·구조 상호작용 해석 SPH 기반의 통합 모델을 개발하여 중대사고 해석에 요구되는 조건을 충족하도록 하였다.

구조해석을 위해 수치안전성이 뛰어난 전라그랑주 SPH(TLSPH)에서 사용할 수 있는 초탄성(Hyperelastic)·탄소성(Elastoplasticity) 모델을 개발하였다. 해당 구성 모델은 탄성영역 내 에너지 소산 문제가 없으므로 탄성 거동에서 탄소성 거동까지 연속적으로 해석 가능하여 범용성을 갖추고 있다. 로그 변형률을 도입함으로써 탄성예측·회귀착점(Elastic Predictor·Return Mapping) 알고리즘 이 단순화되며, 이는 대변형 조건에서도 소변형 이론에서의 알고리즘과 유사한 형태로 구현된다. SPH의 공점(Collocation) 성질로 인해 발생하는 고유의 모래시계 불안정성(Hourglass Instability)을 해상도에 적응하는 비등방성 커널과 JST 기반 안정화 기법을 병용하여 효과적으로 완화하였다. 제안된 구조 모델은 외팔보의 탄소성 굽힘, 투사체 충돌에 의한 탄소성 굽힘, 금속 원기둥의 네킹, 금속 빌렛의 압축열간단조 등 다양한 검증 사례를 통해 기존 실험 및 수치 해석 결과와의 일치성을 확인하였다.

준취성 손상에 대해서는 등방성 손상모델을 적용하였다. 손상변수는 입자 크기, 파괴 에너지, 재료 특성 등을 반영하는 응집력 기반 모델로 공식화되어, 입자 해상도 및 분포에 대한 민감도를 줄이고 수치적 수렴성을 향상시킨다. 손상변수는 단일 스칼라 변수로 정의되며, 명시적인 균열 추적이나 자유도 보강법(Enrichment) 없이 구현 가능하여 계산 알고리즘을 간소화한다. 입자 간 평균 손상도가 임계값을 초과하면 입자간 상호작용을 비활성화함으로써 취성 파손를 모사한다. 해당 손상 및 파손 모델은 모드 I 및 혼합모드 3점 굽힘 실험, Kalthoff–Winkler 실험의 전단 하중, 압축 인장하중조건의 동적균열, 동적 균열 분기 등 다양한 검증 사례를 통해 실험 및 기존 수치 결과와의 정합성을 보였다.

유체 해석기법으로는 점성 항에 대해 암시적 접근법을 사용하는 명시적 비압축성 SPH(EISPH)을 사용하며, 이를 구조 해석기법과 연계하였다. 일반적인 비압축성 유동 해석에서는 압력 푸아송 방정식(PPE)을 전역적으로 풀어 압력을 계산하지만, EISPH는 전역 연립방정식의 해를 요구하지 않고 국소적인 압력 계산을 통해 계산 비용을 절감하면서도 높은 수치 정확도를 유지할 수 있다.

유체·구조 상호작용을 해석하기 위해, 비슬립 조건과 힘의 평형 조건을 만족하는 인터페이스 모델을 구현하였다. 계산 효율을 높이기 위해, 유체 영역과 구조 영역을 서로 다른 시간 간격으로 적분하는 순차적 시차(Sequential Staggered) 시간 적분 기법을 도입하였다. 이 방식에서는 유체 영역을 고유 시간 간격으로 먼저 계산하고, 구조 영역은 더 짧은 시간 간격으로 반복 적분한 뒤, 양 영역을 동일한 물리 시간에서 동기화함으로써 계산의 안정성과 효율성을 동시에 확보할 수 있다.

유체·구조 상호작용 모델은 정수압 하중에 의한 플레이트 변형, 충돌에 의한 쐐기 변형, 슬로싱에 의한 보의 진동 등 다양한 검증 사례에 적용되었다. 해석 결과는 이론 해석, 실험 결과 및 기존 수치 해석 결과와 우수한 정합성을 보여, 사용된 모델이 선행 연구를 정확하게 재현하고 있음을 확인하였다.

격자형(Cartesian-type) 입자 배치는 실제 형상과 불일치하는 인공적인 노치(Notch)를 형성할 수 있으며, 이로 인해 비물리적인 응력 집중과 과도한 영구 변형 및 손상이 발생할 수 있다. 따라서 복잡하거나 곡면을 포함하는 형상을 정밀하게 반영하는 입자 초기분포가 필수적이다. 본 연구에서는 CAD 기반 형상을 바탕으로 입자 분포를 생성하는 입자 생성기를 개발하였다. CAD 파일에서 추출한 형상 정보를 바탕으로, 표면 부근에서는 정확한 형상을 따르도록 입자를 재배치하는 표면 인자 입자 이동 기법(Surface-aware PST)을 적용하고, 내부에서는 균일한 분포를 유지한다. 결과적으로, 이러한 고품질 입자 분포는 응력장 및 손상장의 비정상적 발달을 완화시키고, 전체 해석의 신뢰도를 높이는 데 기여한다.

개발된 SPH-FSI 통합 해석기법은 실제 원전 규모의 중대사고 시나리오에 적용되어 그 활용 가능성을 입증하였다. 노심 용융물의 RPV 내 유지(In-vessel Retention) 단계에서는 산화물 및 금속 용융물의 자연대류에 의한 대류 열전달이 압력용기(RPV) 하부 헤드의 상부에 집중되며, 이로 인해 해당 영역은 열침식에 취약해지고 전반적인 용기 파손으로 이어질 수 있다. 본 연구에서는 RPV 하부 헤드의 단면 파단(Guillotine-type Break)을 가정하여 RPV와 용융물이 동시에 외부 코어캐처에 충돌하는 시나리오를 설정하였다. 이 사건의 발생 가능성은 극히 낮지만, 코어캐처의 구조적 건전성과 기능적 성능에 가장 큰 위험도를 가지므로 안전성 평가가 필요하다.

RPV의 충돌각, RPV 초기 온도, RPV 전단 파단 위치 등을 변수로 민감도 분석을 수행하였다. 주로 용융물의 재배치 거동과 RPV 및 코어캐처의 구조 응답을 정량화하는 데 초점을 두었다. RPV와 용융물의 총 낙하질량이 매우 크므로 충돌영역 근처에서 희생 콘크리트의 국소적인 관통이 발생했지만, 하부의 탄소강층은 구조 건전성을 유지하였다. 이는 고온 상태의 RPV에 비해 콘크리트의 강도가 낮고, 상온의 탄소강이 매우 높은 강도를 갖기 때문이다. 탄소강의 건전성이 유지된다는 점은 하부에 설치된 피동 노외 노심용융물냉각계통(PECS)의 기능이 보존될 가능성이 높음을 시사한다. RPV의 온도가 높을수록 RPV 내에서의 소성 에너지 소산이 증가하며, 그 결과 용융물의 분산이 줄어들었다. 반면, RPV의 충돌각이 클수록 용융물의 분산이 증가하는 경향을 보였다. 이러한 결과는 코어캐처 설계 여유를 설정하고, 규제기관의 안전성 평가 기준 수립에 중요한 정보를 제공할 것으로 기대된다.

이상의 결과를 종합하면, 본 연구에서 제안된 SPH-FSI 통합 해석기법은 원자력 분야 중대사고 해석을 위한 실용적이고 고신뢰도의 수치해석 도구로서의 타당성을 확보하였다.

다국어 초록 (Multilingual Abstract)

Severe accidents in nuclear power plants refer to scenarios in which core damage occurs, potentially leading to the release of radioactive materials. Although the probability of such events is extremely low, their consequences can be catastrophic, with significant impacts on plant safety and public trust. Historical events such as the accidents at Three Mile Island (1979), Chernobyl (1986), and Fukushima (2011) have clearly demonstrated the far-reaching implications of severe accidents on the nuclear industry.

Accurate prediction and analysis of severe accidents are therefore essential for improving nuclear safety. They enable the development of effective mitigation strategies, ensure the structural integrity of containment systems, and help minimize the release of radioactive materials to the environment.

Severe accident analysis involves the complex interplay of various physical phenomena: high-temperature, highly viscous fluids, free-surface flows, solid-liquid phase change, elastoplasticity-creep-fracture behaviors, and thermo-fluid-structure interaction (thermo-FSI). Because each mechanism feeds back into the others, a credible safety assessment should resolve thermo-FSI within one tightly coupled framework. Key features for numerical framework considering the above-mentioned multiphysics are accuracy, stability, efficiency, and straightforward implementation. As a mesh-free solver, smoothed particle hydrodynamics (SPH) offers great advantages to treat free-surface flows, phase change, collision of multiple complex bodies, and structural damage and failure. In this respect, this study presents a thermos-FSI model based exclusively on SPH to meet these requirements.

A finite-strain hyperelastic-plastic model is formulated in total Lagrangian SPH (TLSPH)-compatible manner, which enables purely elastic to elastoplastic analysis without elastic energy dissipation issue. By using the logarithmic strain, the elastic predictor/return mapping algorithm is significantly simplified, which closely mirrors the formulation used in infinitesimal-strain plasticity. Spurious hourglass instability, induced by the collocated nature of SPH, is effectively mitigated by a resolution-adaptive anisotropic kernel combined with a JST-based stabilization scheme. The hyperelastic-plastic model is validated through benchmarks including elastoplastic bending of a cantilever beam, elastoplastic bending by projectile impact, necking of a cylindrical bar, and thermo-upsetting of a metallic billet, which shows good consistency with experimental and numerical data reported previously.

For quasi-brittle damage, an isotropic damage model is employed. The damage evolution is governed by a cohesive force-based law, which incorporates the effects of particle size, fracture energy, and material properties. This formulation reduces sensitivity to particle resolution and distribution, thereby enhancing numerical convergence. The model requires only a single damage-related scalar parameter and avoids the need for explicit crack tracking or kinematic enrichment, significantly simplifying the computational algorithm. Quasi-brittle fracture is represented by deactivating particle interactions when the average damage between interacting particles exceeds a threshold (~1.0). The damage and fracture models are validated against mode-I and mixed-mode three-point bending tests, dynamic shear loading, dynamic fracture of concrete under compact tension, dynamic crack branching scenarios, showing strong agreement with experimental and numerical results reported in the literature.

For a fluid solver, an explicit incompressible SPH (EISPH) with an implicit viscosity solver is coupled with the structural solver. In incompressible-flow approaches, the flow pressure is generally computed by solving the pressure Poisson equation (PPE). EISPH eliminates the need for solving a global system of the PPE, allowing for a fully explicit time integration while maintaining numerical accuracy. For FSI, an interface model, satisfying no-slip condition and force equilibrium, is implemented. To efficiently reduce the computational cost of the fluid domain, a sequential staggered time integration scheme is adopted. The core idea is to advance the fluid domain using its own time step, while sub-iterating the structural domain with a smaller time step until both domains are synchronized at the same physical time. To verify the implementation, the FSI model is applied to benchmark problems including hydrostatic water loading, wedge slamming, and sloshing-induced beam oscillations. The model results show good consistency with theoretical solutions, experimental measurements, and numerical simulations reported in the literature.

Regular Cartesian-type particle distributions may introduce artificial surface notches that are inconsistent with the actual geometry. These geometric inaccuracies can lead to unphysical stress concentrations, which in turn cause excessive permanent deformation and damage. This limitation underscores the importance of generating particle distributions that more accurately capture curved or complex boundaries in meshless methods. To address this issue, a versatile particle generator is developed to produce body-fitted, notch-free initial particle distributions for complex geometries. Based on surface data extracted from computer-aided design (CAD) files, a surface-aware particle shifting technique (PST) is employed to redistribute particles near the surface, ensuring closer conformity to the exact geometry while preserving uniform spacing within the interior. The resulting high-quality particle distributions lead to smoother stress and damage fields.

The integrated SPH-FSI framework is applied to a plant-scale severe accident scenario to demonstrate its applicability. During the in-vessel retention phase, convective heat transfer, driven by natural convection of oxidic and metallic coriums, may concentrate in the upper part of the reactor pressure vessel (RPV) lower head. The so-called heat focusing makes the upper part vulnerable to thermal ablation, potentially leading to global fracture. Among various global fracture modes, the selected accident scenario assumes that the RPV and entrained corium simultaneously impact an ex-vessel core catcher due to a guillotine-type break of the RPV lower head. Although the probability of this event is extremely low, it poses significant risk to the structural integrity and performance of the core catcher. Therefore, a performance and safety assessment of this scenario is considered worthwhile.

Parametric studies are carried out to evaluate the effects of impact angle (0° to 60°), initial RPV temperature (900 to 1500 °C), and the location of the guillotine break. These studies focus on quantifying corium relocation behavior and the structural response of the RPV and core catcher. While local penetration of the sacrificial concrete is observed due to the large mass of the RPV and corium, the overall structural integrity is preserved by the bottom carbon steel layer. The preserved integrity of the bottom carbon steel layer suggests that the underlying PECS is likely to remain functional. Higher RPV temperatures lead to more plastic dissipation within the RPV, reducing corium dispersion. Larger impact angles, on the other hand, intensify the corium dispersion. These findings contribute to defining design margins and support regulatory evaluation of ex-vessel core catcher strategies.

Collectively, the advances presented herein establish the integrated SPH-FSI framework as a practical, high-fidelity tool for severe-accident analysis in nuclear engineering.

목차 (Table of Contents)

• Chapter 1. Introduction 15
• 1.1. Background and motivation 15
• 1.2. Literature review on an SPH-based FSI model 16
• 1.2.1. Large-strain thermo-elastoplasticity model 16
• 1.2.2. Quasi-brittle damage and fracture model 17
• 1.2.3. Thermo-fluid-structure interaction model 19
• 1.3. Purpose and objectives 19
• 1.4. Research scope 20
• Chapter 2. Fluid model 23
• 2.1. Smoothed Particle Hydrodynamics 23
• 2.2. Fluid solver 24
• 2.2.1. Governing equations 24
• 2.2.2. Weakly-compressible SPH 25
• 2.2.3. Explicit incompressible SPH 29
• 2.3. Implicit viscosity solver 30
• Chapter 3. Large-strain elastoplasticity model 38
• 3.1. Governing equations 38
• 3.2. Strain measure 39
• 3.3. Stress update algorithm for elastoplasticity 42
• 3.3.1. Hypoelastic-based approach 42
• 3.3.2. Hyperelastic-based approach 43
• 3.4. SPH approximation 45
• 3.4.1. Anisotropic kernel function 45
• 3.4.2. Spatial discretization 47
• 3.4.3. Conservation of energy and linear/angular momenta 49
• 3.5. Additional numerical schemes 52
• 3.5.1. Hourglass control scheme 52
• 3.5.2. Continuum collision model 53
• 3.5.3. Time integration scheme 55
• 3.6. Verification and validation of the elastoplasticity model 56
• 3.6.1. Cantilever beam under bending moments at free end 56
• 3.6.2. Elastoplastic bending of a cantilever beam 58
• 3.6.3. Elastoplastic bending of a metal beam by projectile impact 62
• 3.6.4. Necking of a cylindrical metal bar 63
• 3.6.5. Thermo-upsetting of a metallic billet 65
• Chapter 4. Quasi-brittle damage and fracture models 91
• 4.1. Continuum damage mechanics 91
• 4.2. Total Lagrangian SPH with damage mechanics 95
• 4.3. Improved non-local formulation 97
• 4.4. Verification and validation 100
• 4.4.1. Mode-I fracture by three-point bending tests 100
• 4.4.2. Mixed mode fracture by three-point bending tests 102
• 4.4.3. Dynamic shear loading of Kalthoff-Winkler experiment 103
• 4.4.4. Dynamic fracture of concrete subjected to compact tension 105
• 4.4.5. Dynamic crack branching 106
• Chapter 5. Thermo-fluid-structure interaction model 123
• 5.1. Boundary model at fluid-structure interface 123
• 5.1.1. Kinematic boundary model 123
• 5.1.2. Dynamic boundary model 124
• 5.1.3. Thermal boundary model 126
• 5.2. Time integration scheme 127
• 5.3. Verification and validation 128
• Chapter 6. CAD-compatible particle generator 140
• 6.1. Working principle of particle generator 140
• 6.2. Application to complex 3-dimensional geometries 142
• Chapter 7. Core catcher integrity & corium relocation under guillotine-type RPV failure 146
• 7.1. Introduction 146
• 7.1.1. Literature review of an ex-vessel core catcher 146
• 7.1.2. Targeted severe accident scenario 148
• 7.2. Simulation setup 149
• 7.2.1. Constitutive model of molten coriums 149
• 7.2.2. Constitutive model of structures 150
• 7.2.3. Simulation configuration and boundary conditions 154
• 7.2.4. Simulation parameters 155
• 7.3. Simulation results and discussion 156
• 7.3.1. Particle resolution study 156
• 7.3.2. Effect of initial temperature of the RPV 157
• 7.3.3. Effect of impact angle of the RPV 159
• 7.3.4. Effect of global failure location of the RPV 159
• 7.3.5. Effect of ICI nozzles in the lower head of the RPV 159
• Chapter 8. Conclusion 181
• 8.1. Summary 181
• 8.2. Implications in severe accident analysis 182
• 8.3. Recommendations of future study 184
• References 186
• 국문 초록 197