1. Introduction

Large break loss-of-coolant accident (LBLOCA) is a type of pipeline rupture accident that belongs to the design-basis accident category in light water reactors (LWRs). In the early phase of a LOCA blowdown, a rapid depressurization induces pressure waves that propagate and reverberate within the reactor loop and core for a short interval. A LOCA induced propagation of depressurization waves exerts hydrodynamic loads on internal structures, causing loosening of bolts, contact-induced wear, loss of coolable geometry, and excessive localized stress. LOCA analysis is thus vital for reactor safety analysis, and it is of utmost importance to develop an efficient and reliable method to predict the LOCA dynamic response.

Dynamic analysis of LOCA induced structural response was initiated during 1960-70. Krieg et al. [1] summarized the state-of-the-art approaches and developed codes for analyzing the dynamic response of vessel internals. Griggs et al. [2] used a 1D analysis code RELAP presented a solution for forces and deformations of core structures during LOCA. A 1200 MW pressurized water reactor (PWR) and the German Heissdampf Reactor (HDR) were chosen as examples to validate their analysis methods. Belytschko and Schumann [3] summarized the LOCA dynamic analysis reports focusing on the importance of the LOCA phenomenon, appropriate computational methods, and experiments. Wolf [4,5] presented the first systematic experiments and test series from V31.2 to V34 on the HDR facility. Wolf’s experimental results have triggered a bundle of efforts to model FSI in LOCA, and have served as a benchmark for validation of numerical simulations. Andersson et al. [6] presented a numerical simulation of LOCA response by using the ADINA FSI capability. HDR V31.1, boiling water reactor (BWR) streamline break and BWR waterline break were successfully simulated and validated. The effects of FSI on LOCA response should be accounted for. Robbe et al. [7] used the EUROPLEXUS code to model the depressurization wave occurring at the beginning of a LOCA in a 4-Loop PWR. The primary circuit was simplified and represented with a pipe model, enabling the simulation of the entire circuit. Casadei and Potapov [8] described the development of non-conforming interface algorithms for a family of FSI problems. The algorithms circumvent the difficulties of generating a conforming fluid-structure mesh along the interface, a prerequisite for previously developed FSI algorithms, and can substantially increase the flexibility and overall efficiency of the numerical simulation. Brandt et al. [9] developed a fully coupled simulation of the HDR experiments by combining the computational fluid dynamic (CFD) software Fluent and Star-CD with a finite element software ABAQUS. MpCCI was used for two-way coupling and data transfer between the CFD solver and the structural solver. Comparisons between numerical simulations and experiments were made for pressure, displacement, and strain, and good agreements were obtained. Hermansky et al. [10] investigated the dynamic response of the WWER440/V213 reactor vessel internals to large-break LOCA. They used the MSC.Dytran code and the arbitrary Lagrangian Eulerian (ALE) coupling for FSI, with Lagrangian solid element for the internals and Eulerian elements for the water coolant. Faucher et al. [11] described a follow-up work by using the EUROPLEXUS software for a full-scale simulation on LOCA. The complete primary circuit was modeled using a coupled 1D/3D strategy, and the strategy was demonstrated for a French 900 MW PWR. Sheykhi et al. [12-13] compiled a calculation program and combined a finite element software ABAQUS to predict the pressure extremes and stress distributions of different reactor pressure vessels under the LOCA accident.

In a parallel study with the aforementioned efforts where coolants were modeled as viscous or nonviscous, compressible or incompressible fluids, and treated mainly in a CFD-based approach, other alternative methods such as the potential-flow approximation and the acoustic finite element method were also introduced to model the induced pressure wave during LOCA. The basic assumption in this approach lies in the fact that the induced pressure is small, and the fluids can be simplified as inviscid flows to some degree. Ludwig and Schumann [14] adopted this potential flow assumption and modeled the LOCA dynamic response by using the FLUX code together with some in-house (finite element method) FEM codes. Au-Yang et al. [15,16] developed a structural priority approach for weakly coupled simulation of LOCA, where the induced pressure was governed by the acoustic equation, and the FSI could be expressed as an added time-invariant hydrodynamic mass term. In Au-Yang’s model, a PWR was simplified as a 1D lumped mass beam model. However, the structural priority approach is restricted to small structural displacements with cylindrically symmetric fluid-structure boundaries. A subsonic potential-based fluid formulation was implemented in ADINA, and Sussman and Sundqvist [17] used this program to perform FSI analysis of the HDR blowdown experiment that demonstrated the efficiency and accuracy of the formulation. Sommerville [18] and Luke [19] demonstrated the use of ANSYS acoustic elements and benchmark of the application for LOCA analysis. Compared with accurate CFD calculation, the potential-based or acoustic method only needs to solve the acoustic wave equation without solving the complicated Navier-Stokes (N-S) equation. This simplification can dramatically reduce the computation cost. In addition, the calculation efficiency is slightly higher than the CFD fluid-structure coupling method. The use of commercially available acoustic finite element software enables a more cost-effective and transparent analysis of LOCA responses. Even after being established as a simple, efficient, and rational method with enough engineering accuracy, the potential-based or acoustic finite element method is criticized for excluding fluid viscosity and for having no convincing explanation of the damping mechanism. Using this method will over-predict the motion of the core barrel in the late phase of LOCA. The potential-based or acoustic method maybe suitable for early-stage short duration LOCA analysis, where the time elapsed is too short for the effect of damping and response attenuation to become predominant.

The development of computer codes and of commercial software in particular has invoked interest in LOCA and other similar FSI issues. Here, we present a coupling analysis strategy by integrating a CFD solver Fluent, and a structural solver ANSYS, both of which are from the same vendor ANSYS. The motivation in doing so is because Fluent can solve continuity equations, momentum equations, and turbulence equations in a separate solver and ANSYS Mechanical solves deformation equations in a segregated solver. The system coupling under the Workbench framework couples the ANSYS Mechanical and Fluent solvers by running a co-simulation to exchange data and achieve one-way or two-way FSI in a straightforward and seamless way. However, selecting these two ideal candidates poses a difficulty in performing transient FSI analysis of prestressed structures that are not supported by the ANSYS Workbench. The initial prestressed state is essential for a pressurized PWR, therefore, for the LOCA analysis of PWR, the prestressed state must be tackled properly. Here, we report a numerical coupling strategy consisting of an initialization step to attain the prestressed state, and a subsequent restart process to model transient FSI problems. It is implemented via the system coupling in the ANSYS Workbench, and is ideally suitable for pressurized PWRs. We discuss the implementations and the FSI coupling scheme in detail in this work. Rigorous comparison with reported HDR experiments is made for the pressures, displacements, and strains for a number of different positions that validate the correctness and accuracy of the proposed approach. We finally give some concluding remarks and future outlook on this topic.

2. Model parameters and boundary conditions

The model shown in Fig. 1 is a half-model with symmetric boundary conditions. Figure 1 (a) is the CAD drawing of the HDR reactor. A core barrel with a lumped mass ring at the bottom represents the simplified internal structures. The bottom center of the vessel was chosen as the origin of the coordinate system. The blowdown nozzle is located at Z=8.85 m with the length of break nozzle 1.369 m and break diameter 0.2 m. The break nozzle is the only outlet and all other outlets are closed in accordance with the boundary conditions adopted in the model experiment. The total liquid volume of the system is approximately 33.36 kiloliters. Figure 1 (b) shows the fluid (left) and the structural (right) meshes used for simulation. The fluid mesh consists of approximately 290,000 hexahedral grids with a maximum mesh size of 150 mm, and the solid mesh has about 60,000 8-node hexahedral elements with a maximum mesh size of 160 mm. The bottom of the vessel was fixed.

Fig. 1.Full-screen View

(Color online) CAD drawing of the (a) HDR reactor, and (b) fluid (left) and structural (right) meshes

Table 1 summarizes the test conditions of the HDR blowdown experiment V32 chosen for this simulation. In the test, the downcomer temperature was 240 °C, and the core temperature was varied axially from 283 °C to 308 °C. An average temperature of 274 °C was selected for the simulation. The dynamic viscosity of the fluid is approximately μ= 0.0001 Pa∙s in the concerned temperature range.

Figure 2 is a plot of the pressure history of the break nozzle up to 10 ms. After 10 ms, the pressure at the break location attains a constant value. For LOCA, which typically has a very short duration in the order of milliseconds, the break pressure condition is crucial for dynamic response analysis. The pressure drop at the break can be monitored in the experiment either by pressure sensors, or estimated by some system codes such as RELAP or APROS [9]. We are concerned with the numerical strategy for LOCA dynamic analysis, and thus we adopt a simplified pressure condition (the solid red line) at the end of the break nozzle as a pressure boundary condition for the Fluent simulation. We thus obtain the results presented in Figs. 4-7. Also shown in Fig. 2 are plots of the experimental measurements, and we used these as the pressure boundary condition to calculate and obtain the results of Fig. 8. Figure 2 indicates that our simplified pressure curve captures the overall trend of the pressure drop during LOCA.

Fig. 2.Full-screen View

LOCA induced pressure drop condition at the break nozzle obtained from experiment with linear fitting of the experiment [4].

Fig. 4.Full-screen View

(Color online) Plots of pressure for (a) von Mises stress (b) contours at selected instants of time.

Fig. 8.Full-screen View

Effect of pressure boundary conditions. Pressure on the downcomer wall at Z =8.85 m, φ=90° (left) and Relative radial displacement between the core barrel and RPV at Z =7.15 m, φ=90° (right).

The initial pressure of the coolant of 11 MPa renders the vessel in a pressurized state with initial deformation and prestress. For a transient FSI analysis, this initial displacement and the prestress must be tackled properly. With regard to the numerical analysis of prestressed structures, Chen et al. [20] have developed a simplified technique for modal analysis of symmetrical prestressed structures using group theory, which has been verified in 2D and 3D structures. In Timperi's previous work [21] within the framework of ABAQUS and Fluent with data transfer managed by MpCCI, a separate structural static load analysis step was needed to calculate the equilibrium prestressed state of the vessel, after the static load step. The applied pressure was provided by CFD code, and then a transient FSI coupling process followed. The introduction of a data transfer code makes the implementation and debugging of FSI more complicated. It is desirable to utilize an integration of Fluent and ANSYS for coupled FSI analysis for a variety of engineering problems, since both codes belong to one single company ANSYS. Unfortunately, System Coupling within the ANSYS Workbench framework cannot specify multiple load steps in ANSYS Mechanical, which means that System Coupling cannot provide a pre-stressed structural model before FSI coupling as in the above ABAQUS-MpCCI-Fluent framework. A direct two-way FSI analysis with an initial pressure of 11 MPa results in an unphysical pressure drop when a compressible fluid model is used. The unphysical pressure drop in the earlier steps right after the initiation of computation is due to the initial deformation of the vessel wall, and the resultant gap between the fluid and solid interface. The unphysical pressure drop due to this numerical artifact would ruin the following FSI analysis, and this artifact must be ruled out.

We propose a numerical strategy shown schematically in Fig. 3 to solve this issue within the framework of ANSYS Workbench. The numerical procedure consists of two transient analysis processes. The first process, i.e., an initialization process, initializes a Fluent pressure P = 11 MPa, runs transient analysis with only force data transfer from fluid to solid, whereas the displacement data transfer from solid to fluid is suppressed. The first transient analysis terminates after enough numbers of iterations, after an equilibrium state is achieved. This one-way transient analysis attains an equilibrium state with unchanged fluid pressure and a prestressed structure. Then, the second transient analysis is started with the prestressed state as the initial condition, and loaded pressure boundary condition. This runs with two-way data transfer in a conventional manner. Thus, the structure deforms because of the propagation of pressure waves, and the fluid pressure drops further because of the deformation of the interface. Using this initialization and restart scheme, the prestressed state of the vessel is properly handled, and the FSI analysis can be performed seamlessly by virtue of the ANSYS Workbench.

Fig. 3.Full-screen View

(Color online) An initialization and restart strategy for transient FSI analysis of a prestressed structure.

The simulations were run on a PC with Intel(R) Xeon(R) Silver 4114 CPU 2.20 GHz processor in the Windows environment. A time step $�$t = 1×10^-3 s was used, and the first transient analysis ran 4,000 steps to guarantee equilibrium. The second transient analysis was run with time step $�$t = 5×10^-5 s for 0.2 s. The whole simulation took approximately 123 hours to complete.

The simulation parameters used are summarized as follows. The coolant was treated as a compressible fluid, and the following pressure-density relation was adopted:

where ρ₀ = 758.5 kg/m³ and p₀ =11 MPa are the reference density and pressure, respectively. The bulk modulus was set as κ = 866.6 MPa. The fluid was assumed to be isothermal, and 274 °C was chosen as the average temperature for this simulation. The dynamic viscosity was approximately μ = 0.0001 Pa·s for this average temperature. The commercial CFD software Fluent was used for the simulation with the standard large Reynolds number k-$\epsilon$ model. The PISO method was used for pressure correction. The convection term in the momentum and turbulent flow equations was discretized using Fluent’s third-order MUSCL scheme, and a diffusion term by a second-order central difference scheme. ANSYS Mechanical was used for finite element simulation. Table 2 lists the material parameters used in the simulation. For the core barrel, isotropic linear elastic properties with Young’s modulus E = 175 GPa, Poisson’s ratio ν = 0.3, and density ρ = 7,900 kg/m³ were adopted. For the RPV, the corresponding material parameters were E = 190 GPa, ν = 0.3, and ρ = 7,850 kg/m³. A stiffness propor tional Rayleigh damping was used for structural dynamics modeling with a damping coefficient β = 6.4×10^-5, which yields nearly 2 % model damping at a frequency of 100 Hz.

3. Results and discussions

Figure 4 shows the contours of pressure at (a) von-Mises stress, and (b) at selected instants of time. From Fig. 4 (a), we see that a depressurization wave initiates at the break nozzle and propagates along the break pipe and the downcomer region, and eventually spreads out over the whole domain of the coolant. Figure 4 (b) is a plot of the stress contour of the internal structures at various instants of time, from t = 0 ms to 115 ms. It is clear that the vessel undergoes an initial deformation and attains a prestress of about 160 MPa, while the core barrel is immersed in the coolant and is stress-free. The pressure of the fluid enclosed by the vessel decreases with the development and propagation of the depressurization wave, and thus, the stress level of the vessel decreases. Meanwhile, the core barrel, which is suspended inside the pressure vessel, pendulates because of the propagation of the pressure wave. This pendulation induces stress in the upper portion of the core barrel, in particular at the clamped end.

To further validate our numerical strategy and the two-way FSI results, we made a systematic comparison between the two-way FSI results with the one-way FSI results, and the experimental measurement reported by Wolf et al. [4,5]. For the one-way FSI simulation, a transient CFD simulation was carried out with the solid structure being treated as a rigid wall, and then the time-dependent fluid pressure was exerted to the elastic structure to calculate the deformation and the stress distribution.

Figure 5 (a) shows a plot of the pressure on the downcomer wall at Z=8.85 m, φ=90° (left), and at Z=2.3 m, φ=90° (right). Figure 5 (b) shows a plot of the pressure on the downcomer wall at Z = 7.78 m, φ = 270° (left), and on the core axis Z = 5.05 m (right). The inset picture shows the schematic of the reactor and the location where the pressure is compared is marked by a red circle. Good agreement is achieved as compared to measurements, especially in the two-way FSI up to 0.1 s, when the pressure drops from the initial 11 MPa to about 8 MPa for the position of interest. Afterwards, there is a noticeable discrepancy between the simulation and measurements. This is due to the limitation of the single-phase model used in the simulation. After t = 0.1 s, flash evaporation occurs. This is a physical phenomenon that when the ambient pressure suddenly drops below the initial saturation pressure of the liquid, the liquid changes from the initial equilibrium state to a superheated state. Due to the rapid drop in pressure, the whole energy cannot be contained in the liquid as sensible heat, and the heat surplus is converted into latent heat of vaporization. During this process, violent phase changes can be observed. In this state, the fluid temperature is higher than the boiling point at the particular pressure, and it rapidly boils and vaporizes to form a two-phase flow. Note that the critical phase-transition pressure of water at 274 °C is exactly about 8 MPa, implying that part of the coolant transfers into the two-phase system, and the disagreement after 0.1 s between model prediction and experiment is a rational result as expected.

Fig. 5.Full-screen View

(Color online) Validation of the two-way FSI analysis by comparing the experimental results with the one-way FSI results. (a) Pressure on the downcomer wall at Z =8.85 m, φ=90° (left) and at Z =2.3 m, φ=90° (right). (b) Pressure on the inner downcomer wall at Z = 7.78 m, φ = 270° (left) and on the core axis Z = 5.05 m, φ = 270° (right). (c) Relative radial displacement between the core barrel and RPV at Z =7.15 m, φ=90° (left) and at Z =2.3 m, φ=90° (right). (d) Pressure difference across the core barrel wall at Z =8.85 m, φ=90° (left) and at Z =5.55 m, φ=90° (right).

Figure 5 (c) shows a plot of the relative radial displacement between the core barrel and the vessel at Z=7.15 m, φ=90° (left), and at Z=2.3 m, φ=90° (right). Regardless of the limitation of the single-phase model for pressure prediction after phase-transformation, the proposed two-way FSI simulation gives a surprisingly good correlation with measurements, while the one-way FSI simulation exhibits very wavy non-smooth results. The smooth two-way FSI results and the non-smooth one-way FSI results obtained in this study are similar to those obtained by Casadei [8] and Brandt [9] using different codes and coupling strategies. In a fully coupled simulation, the wall movement damps the pressure oscillation and gives smooth results. However, in the one-way FSI simulation, the fluid pressure does not react to wall movement and gives unrealistic modulation because of the sensitivity of the fluid pressure to fluid volume change. These conclusions are also reflected in the pressure difference results because an important indicator for measuring fluid-solid interaction in the HDR blowdown test is the pressure difference across the core barrel, as it creates a hydrodynamic load in the structure.

Figure 5 (d) is a plot of the pressure difference across the core barrel at Z=8.85 m, φ=90° (left), and at Z=5.55 m, φ=90° (right). It should be noted that the one-way FSI simulation not only overestimates the amplitude of the hydrodynamic load, but also predicts a completely wrong time history. The calculation results of the two-way FSI simulation agree well with the experiment during the whole period, that is, the pressure difference is basically unaffected by the phase change.

Figure 6 compares the relative radial displacement between the core barrel and the RPV at various locations. Positive relative displacement represents a reduction in distance between the walls and vice versa.

Fig. 6.Full-screen View

Comparison between predicted displacements and experiments for various positions. (a) Relative radial displacement between the core barrel and RPV at Z =7.15 m, φ=90° (left) and φ=270° (right). (b) Relative radial displacement between the core barrel and RPV at Z =5.55 m, φ=90° (left) and φ=270° (right). (c) Relative radial displacement between the core barrel and RPV at Z =2.3 m, φ=90° (left) and φ=270° (right).

Figure 7 shows a comparison of the strains between the numerical simulations and the experimental measurements at different locations. The strain oscillates at a frequency of 200 Hz, which is consistent with the frequency of the pressure wave propagating back and forth in the nozzle. This oscillation is not shown in the experiment because of the boiling inside and outside of the nozzle, which makes the water "soft" in this region.

Fig. 7.Full-screen View

Comparison between predicted strains and experiments for different positions. (a) Hoop strain on the inner surface of the core barrel wall Z =5.45 m, φ=135° (left) and at Z =5.45 m, φ=225° (right). (b) Hoop (left) and axial (right) strain on the outer surface of the core barrel at Z = 8.85 m, φ=90°.

In Figs. 6 and 7, only the two-way FSI results are presented and compared with experiments. The overall behavior of the calculated displacements and strains correspond quite well with the experiment.

Figure 8 shows plots of the pressure field and displacement of a two-way FSI simulation under linear fitting and experimentally measured pressure boundary conditions. Changing the input conditions has little effect on the resulting pressure field and displacement.

4. Concluding Remarks

LOCA analysis tools and codes have evolved from in-house developed codes to integration of commercial software. One particular interest is to use an integrated FSI environment such as workbench for transient analysis of pressurized structures. We describe here a numerical strategy based solely on the system coupling in the ANSYS workbench toward this goal. The strategy consists of two transient FSI analyses, one with partial data transfer being suppressed to achieve the equilibrium prestressed state, and the other with a conventional two-way data transfer between the solid and the fluid, and restarted after the initialization process. We explained in detail the implementation of the numerical strategy and the data transfer during the two transient processes.

To validate the correctness and feasibility of the proposed strategy, we analyzed the German HDR test V32 since we could access experimental data reported in literature. We compared the pressure history, relative displacement between the core barrel and the vessel, and strains at different locations taken from experiments and given by the two-way FSI simulation. The results of the one-way FSI analysis are also given for comparison. Excellent agreement is achieved between simulation and experiment. The obtained results were not sensitive to pressure boundary conditions. Because of the good popularity and acknowledgment of ANSYS and Fluent in both academia and industry, our strategy is an appealing candidate for efficient and reliable LOCA analysis. The strategy can also be used straightforwardly in the dynamic analysis of any prestressed structures, including modal analysis, transient analysis, and other flow-induced vibration problems. Future work may include extension of the present single-phase model to a two-phase model, and applying the method to analyze dynamics of internal structures that are not fixed rigidly, but are joined by bolts with contact and friction.