Evaluation of different Krylov subspace methods for simulation of the water faucet problem

NUCLEAR ENERGY SCIENCE AND ENGINEERING

Evaluation of different Krylov subspace methods for simulation of the water faucet problem

Hong-Yang Wei ，

Kevin Briggs，

Victor Quintanilla，

Yi-Tung Chen

Nuclear Science and Techniques

Vol.32, No.5

Article number 44

Published in print 01 May 2021

Available online 07 May 2021

DOI：10.1007/s41365-021-00888-w

43600

In this study, a one-dimensional two-phase flow four-equation model was developed to simulate the water faucet problem. The performance of six different Krylov subspace methods, namely the generalized minimal residual (GMRES), transpose-free quasi-minimal residual, quasi-minimal residual, conjugate gradient squared, biconjugate gradient stabilized, and biconjugate gradient, were evaluated with and without the application of an incomplete LU (ILU) factorization preconditioner for solving the water faucet problem. The simulation results indicate that using the ILU preconditioner with the Krylov subspace methods produces better convergence performance than that without the ILU preconditioner. Only the GMRES demonstrated an acceptable convergence performance under the Krylov subspace methods without the preconditioner. The velocity and pressure distribution in the water faucet problem could be determined using the Krylov subspace methods with an ILU preconditioner, while GMRES could determine it without the need for a preconditioner. However, there are significant advantages of using an ILU preconditioner with the GMRES in terms of efficiency. The different Krylov subspace methods showed similar performance in terms of computational efficiency under the application of the ILU preconditioner.

Water faucet problemKrylov subspace methodsILU preconditioner

1 Introduction

Two-phase flow behavior is one of the most important transport phenomena in industry and is found in several systems, including nuclear power plants. Under normal operating conditions, two-phase flow occurs in the primary loop of the boiling water reactor and the secondary loop of the pressurized water reactor. The evaluation of two-phase flow behavior is related to the design of nuclear reactors in both steady and transient states, safety evaluation, and analysis of events such as loss-of-coolant accidents (LOCA) and reflooding [1, 2, 3]. Two-phase flow models are often used to evaluate complex thermal–hydraulic phenomena in facilities, such as the steam–gas pressurizers in nuclear power plants [4]. One-dimensional two-phase flow behavior simulation is widely employed owing to its efficiency. Numerical model system codes such as RELAP5, RELAP7, and CATHAR are widely used for evaluating the safety of nuclear reactors, and the numerical solvers that are used to calculate the behavior of one-dimensional two-phase flow play an important role in terms of code performance. The Gauss elimination linear equation solver is used in TRACE, while CATHAR utilizes the Newton-Raphson method. Recently modified or developed system codes also use high-efficiency matrix solvers or numerical methods. For example, border-profile lower upper (BPLU), which is now used with RELAP5, has a significantly better performance than previous numerical methods [1, 5]. The nonlinear numerical solver JFNK is used for simulations that involve RELAP7 [1, 6]. Hajizadeh et al. also developed a new two-phase flow simulation model in which the SIMPLE algorithm was applied to the drift flux model [7].

The Krylov subspace methods are widely used in numerical simulations including the simulation of two-phase flow. The Krylov subspace methods are considerably more efficient than traditional numerical solvers such as Gaussian elimination or row reduction. Dawson et al. conducted a numerical investigation of three-dimensional two-phase flow behavior [8], in which the finite difference method was used for space discretization, and the backward Euler method was applied for time discretization. The generalized minimal residual (GMRES) method was used, and the efficiency of the Krylov subspace method was evaluated. Nordsveen et al. developed a two-phase model for the investigation of gas-liquid flow behavior using a semi-implicit time-discretization method [9]. Staggered and finite difference methods were used for space discretization. Both the GMRES and the biconjugate gradient stabilized (BICGSTAB) Krylov subspace methods were used in the evaluation and comparison with the Gaussian band solver. The BILU3D preconditioned BICGSTAB method was used by Downar et al. for the numerical solution of multi-dimensional two-fluid equations in the EPRI code VIPRE-02 to investigate computational efficiency [10]. The implemented numerical method was used to simulate a PWR main steam line break problem. The results of simulation with BICGSTAB indicate that the average number of inner iterations required per outer iteration is approximately 50 % less than that via application of the ADI method. To increase the computational efficiency of the two-phase flow simulation in porous media, Bergamaschi et al. evaluated the effects of using a preconditioner with the Krylov subspace method. A finite element method was used for the discretization of the two-phase flow models [11]. Both 2D and 3D numerical tests were conducted to evaluate the computational performance. Based on the simulation results, the number of linear iterations required was reduced using the proposed acceleration in terms of convergence. Mohitpour et al. evaluated different Krylov subspace methods (GMRES, flexible GMRES (FGMRES), direct quasi-GMRES (DQGMRES), conjugate gradient on the normal equations, biconjugate gradient (BCG), and transpose-free quasi-minimal residual (TFQMR)) for the simulation of two-phase flow in pressurized water reactor (PWR) fuel bundles [12]. The two-phase flow behavior was simulated using a three-dimensional two-fluid one-pressure model. The finite difference method was used for space discretization with semi-implicit time discretization. No significant deviation in the simulation of two-phase flow behavior was predicted by these Krylov subspace methods. In addition, GMRES, FGMRES, and DQGMRES exhibited better performance in terms of efficiency and stability than the other three Krylov subspace methods. The Krylov subspace method has also been used in other aspects of nuclear engineering simulation. For example, Liu et al. used the equivalent low-order angular flux nonlinear finite difference method with right preconditioned GMRES in a method of characteristics transport calculation [13].

The water faucet problem was first proposed by Ranson for the investigation of two-phase flow behavior [14] and has been widely used in evaluating the performance of two-phase flow numerical simulations. However, the analytical solution has only considered the void fraction and liquid phase velocity, and this has yielded limitations in terms of validating the numerical simulation of the gas flow behavior and pressure distribution that occur in water faucet two-phase flow. Zou et al. obtained an analytical solution for both liquid and gas using the water faucet problem [15]. The liquid and gas pressures were assumed to differ, and the analysis led to a new solution for the gas velocity and pressure distribution, which could be beneficial for the further verification of two-phase numerical simulations that are based on the water faucet problem.

Nourgaliev et al. used a characteristics-based method for a two-fluid model numerical solution [16]. The primary step in this method was the separation of the conservative and non-conservative parts in the governing equations. The non-conservative part was implemented into the source term, and the conservative part was treated using the characteristics-based method. The water faucet problem was used as a benchmark. Various aspects, such as the time discretization, treatment of conservative terms, treatment of non-conservative terms, and convergence of the numerical solutions, were then evaluated. Gallouet et al. conducted a numerical investigation of two-phase flow behavior based on a two-fluid two-pressure model. The convective, source, and diffusive terms were considered in the governing equations [17]. The finite volume method was used as the numerical solution method, and one-dimensional conditions were considered. Both the Rusanov scheme and VFRoe-ncv scheme, which is an approximate Godunov scheme, were used. Numerical tests, such as the shock tube test, water faucet problem test, and sedimentation test, were then conducted, including the VFRoe-ncv scheme test, which was conducted using the water faucet problem to obtain results including the air volume fraction and water velocity. Morin et al. developed a Roe scheme for the solution of a six-equation two-fluid model [18], in which the water faucet problem was used for the numerical experiments and validation. Fullmer et al. conducted a numerical investigation of one-dimensional two-phase flow by considering a higher-order finite difference method [19]. Owing to the diffusion problem that results from the first-order method, a higher-order method was used and evaluated. The water faucet problem was again used for benchmark purposes. Delchini et al. used RELAP-7 to simulate the water faucet problem [20]. The simulation of two-phase flow and heat transfer with RELAP-7 was based on the two-pressure seven-equation model. RELAP-7 applies the all-Mach flow entropy viscosity method (EVM) for numerical stabilization. The water faucet problem was used as the benchmark to investigate both the two-phase model and the numerical stabilization ability of RELAP-7. Chen et al. conducted a numerical investigation of the two-phase flow using a one-step coupled solution method [21]. The water faucet problem and subcooled boiling problem were used to validate the developed code, and the effect of different nodes on the gas void fraction was investigated with the water faucet problem as a benchmark. In addition, a comparison of the numerical and theoretical results in terms of the void fraction, liquid velocity, and gas velocity was conducted, and the numerical results were consistent with the theoretical results.

Owing to the complexity of two-phase flow behavior, further evaluation of the performance of the Krylov subspace method in two-phase flow simulation is still necessary. Therefore, it is important to evaluate the different Krylov subspace methods that are used to simulate two-phase flow behavior. The water faucet problem is selected as the benchmark for these simulations because it is widely used as a benchmark for the numerical simulation of two-phase flow. The results of this study can be used as a reference for the numerical simulation of one-dimensional two-phase flow in terms of the selection of numerical solvers.

2 Simulation methods

2.1 Governing equations

The governing equations that are used to simulate water faucet behavior are the simplifications of RELAP5’s mass and momentum conservation equations [22]. The terms describing phase change, interfacial heat transfer, and interfacial friction are ignored in tests associated with the water faucet problem.

The gas mass conservation equation is expressed as:

\frac{\partial}{\partial t} (α_{g} ρ_{g}) + \frac{\partial}{\partial x} (α_{g} ρ_{g} u_{g}) = 0

(1)

The liquid mass conservation equation is expressed as:

\frac{\partial}{\partial t} (α_{f} ρ_{f}) + \frac{\partial}{\partial x} (α_{f} ρ_{f} u_{f}) = 0.

(2)

The gas momentum conservation equation is expressed as:

α_{g} ρ_{g} \frac{\partial u_{g}}{\partial t} + \frac{1}{2} α_{g} ρ_{g} \frac{\partial u_{g}^{2}}{\partial x} = - α_{g} \frac{\partial P}{\partial x} + α_{g} ρ_{g} g_{x} .

(3)

The liquid momentum conservation equation is expressed as:

α_{f} ρ_{f} \frac{\partial u_{f}}{\partial t} + \frac{1}{2} α_{f} ρ_{f} \frac{\partial u_{f}^{2}}{\partial x} = - α_{f} \frac{\partial P}{\partial x} + α_{f} ρ_{f} g_{x} .

(4)

2.2 Discretization methods

A fully implicit method was used for time discretization, and the first-order space discretization method was coupled with the staggered mesh method for space discretization. Based on the analytical solution obtained by Zou et al. [15], the direction of the gas velocity will be opposite to the direction of the liquid velocity. Therefore, the developed code should have the capacity to capture the opposite flow directions of liquid and gas.

The gas mass conservation equation can be discretized as follows:

(1) If $u_{g, j + \frac{1}{2}} > 0$ ,

\frac{α_{g, j} ρ_{g, j} - α_{g, j}^{n} ρ_{g, j}^{n}}{Δ t} + \frac{α_{g, j} ρ_{g, j} u_{g, j + \frac{1}{2}} - α_{g, j - 1} ρ_{g, j - 1} u_{g, j - \frac{1}{2}}}{Δ x} = 0.

(5)

(2) If $u_{g, j + \frac{1}{2}} > 0$ ,

\frac{α_{g, j} ρ_{g, j} - α_{g, j}^{n} ρ_{g, j}^{n}}{Δ t} + \frac{α_{g, j + 1} ρ_{g, j + 1} u_{g, j + \frac{1}{2}} - α_{g, j} ρ_{g, j} u_{g, j - \frac{1}{2}}}{Δ x} = 0.

(6)

The liquid mass conservation equation can be discretized as follows:

\frac{α_{f, j} ρ_{f, j} - α_{f, j}^{n} ρ_{f, j}^{n}}{Δ t} + \frac{α_{f, j} ρ_{f, j} u_{f, j + \frac{1}{2}} - α_{f, j - 1} ρ_{f, j - 1} u_{f, j - \frac{1}{2}}}{Δ x} = 0.

(7)

The gas momentum conservation equation can be discretized as follows:

(1) If $u_{g, j + \frac{1}{2}} > 0$ ,

\begin{matrix} \frac{\bar{{(α_{g} ρ_{g})}_{J + \frac{1}{2}}} u_{g, j + \frac{1}{2}} - \bar{{(α_{g} ρ_{g})}_{J + \frac{1}{2}}} {(u_{g, j + \frac{1}{2}})}^{n}}{Δ t} + \frac{1}{2} {(α_{g} ρ_{g})}_{j} \frac{{(\bar{u_{g, J + 1}})}^{2} - {(\bar{u_{g, J}})}^{2}}{Δ x} \\ = - \bar{{(α_{g})}_{J + \frac{1}{2}}} \frac{P_{j + 1} - P_{j}}{Δ x} + \bar{{(α_{g} ρ_{g})}_{J + \frac{1}{2}}} g_{x} . \end{matrix}

(8)

As an example $\bar{{(α_{g} ρ_{g})}_{J + \frac{1}{2}}}$ , can be expressed as $\frac{\bar{{(α_{g} ρ_{g})}_{J}} + \bar{{(α_{g} ρ_{g})}_{J + 1}}}{2}$ .

(2) If $u_{g, j + \frac{1}{2}} \leq 0$ ,

\begin{matrix} \frac{\bar{{(α_{g} ρ_{g})}_{J + \frac{1}{2}}} u_{g, j + \frac{1}{2}} - \bar{{(α_{g} ρ_{g})}_{J + \frac{1}{2}}} {(u_{g, j + \frac{1}{2}})}^{n}}{Δ t} + \frac{1}{2} {(α_{g} ρ_{g})}_{j + 1} \frac{{(\bar{u_{g, J + 1}})}^{2} - {(\bar{u_{g, J}})}^{2}}{Δ x} \\ = - \bar{{(α_{g})}_{J + \frac{1}{2}}} \frac{P_{j + 1} - P_{j}}{Δ x} + \bar{{(α_{g} ρ_{g})}_{J + \frac{1}{2}}} g_{x} . \end{matrix}

(9)

The liquid momentum conservation equation can be discretized as follows:

\begin{matrix} \frac{\bar{{(α_{f} ρ_{f})}_{J + \frac{1}{2}}} u_{f, j + \frac{1}{2}} - \bar{{(α_{f} ρ_{f})}_{J + \frac{1}{2}}} {(u_{f, j + \frac{1}{2}})}^{n}}{Δ t} + \frac{1}{2} {(α_{f} ρ_{f})}_{j} \frac{{(\bar{u_{f, J + 1}})}^{2} - {(\bar{u_{f, J}})}^{2}}{Δ x} \\ = - \bar{{(α_{f})}_{J + \frac{1}{2}}} \frac{P_{j + 1} - P_{j}}{Δ x} + \bar{{(α_{f} ρ_{f})}_{J + \frac{1}{2}}} g_{x} . \end{matrix}

(10)

2.3 Krylov subspace method [23, 24]

The Petrov-Galerkin condition can be applied to obtain the solution for the linear system $A x = b$ via the following expression:

b - A x_{m} ⊥ L_{m},

(11)

Where $x_{m}$ is in the subspace $x_{0} + K_{m}$ with m dimensions, and $x_{0}$ is the initial guess. The method is considered to be a Krylov subspace method when subspace $K_{m}$ is the Krylov subspace.

Generally, there are two types of Krylov subspace methods: one in which $L_{m} = K_{m}$ , and the other where $L_{m} = K_{m} (A^{T}, r_{0})$ . Several different Krylov subspace methods are used, of which GMRES, BICGSTAB, quasi-minimal residual (QMR), TFQMR, conjugate gradient squared (CGS), and biconjugate gradient (BICG) were used in the simulation. However, when compared to direct numerical solvers, the main disadvantage of iterative solvers such as the Krylov subspace methods is their robustness. The use of a preconditioner renders the Krylov subspace method advantageous in terms of both efficiency and robustness. Therefore, the incomplete LU factorization (ILU) preconditioner is also investigated in this work. For sparse matrix A, the ILU factorization process calculates a sparse upper triangular matrix L and a sparse lower triangular matrix U with certain constraints for application in the residual matrix R=LU-A. The numerical programming code was developed using MATLAB 2019, and the numerical solvers used in this code were solvers that are implemented in MATLAB 2019.

3 Simulation results analysis

The residual analysis, the simulation results of the liquid and gas flow behavior in the water faucet, and the solution time were obtained to analyze the numerical performance of the different Krylov subspace methods with and without an ILU preconditioner. The effects of the timestep and node number were also noted during evaluation.

3.1 Residual analysis

Figure 1 shows the relative residuals for different Krylov subspace methods using different node numbers with and without the application of the ILU preconditioner. The timestep was set to 1.0×10^-3 s. The results indicate that the convergence performance of the Krylov subspace methods is much better with an ILU preconditioner than without it. The various Krylov subspace methods demonstrate differing performances when run without the ILU preconditioner. The relative residual is much higher when using CGS and BICG than it is with the other Krylov subspace methods, which indicates that CGS and BICG have a serious robustness issue when simulating the water faucet problem. The relative residual was observed to reach a level of approximately 1.0×10^-10 under the GMRES, which means that GMRES may perform well even without the application of the ILU preconditioner. For cases in which the ILU preconditioner was not used, the calculation performances of TFQMR, QMR, and BICGSTAB were similar, with better results than those obtained with CGS and BICG, and worse than the calculation performance of GMRES. The performance of the GMRES method was quite different from that of the other methods when an ILU preconditioner was used. This is because the other Krylov subspace methods reached the required relative residual after only a few iterations; however, GMRES could not reach the required relative residual when the ILU preconditioner was utilized. This is mainly because the results were the same for two consecutive iterations, after which GMRES stagnated. GMRES therefore appears to have irregular or divergent performance under certain conditions, and similar behavior has also been observed in other studies [12]. CGS exhibited a better calculation performance than the other Krylov subspace methods, regardless of the node number. TFQMR and BICGSTAB demonstrated similar calculation performances, as did QMR and BICG; however, no significant difference was observed among these methods. No significant effect on the calculation performance was associated with the node number used with the different Krylov subspace methods.

Fig. 1

(Color online) Relative residual for different Krylov subspace methods with and without the application of an ILU preconditioner with different node numbers

Figure 2 shows the relative residual for different Krylov subspace methods with and without the ILU preconditioner at different time steps. The convergence performance of the Krylov subspace methods with an ILU preconditioner is much better than that without an ILU preconditioner, regardless of the time step investigated. The GMRES generally has a better performance than the other Krylov subspace methods without the application of the ILU preconditioner in this study. The relative residual was observed to reach 1.0×10^-10 when the GMRES method was used without the ILU preconditioner. The calculation performances of CGS and BICG are much worse than those of the other Krylov subspace methods, while TFQMR, QMR, and BICGSTAB have similar calculation performances. However, GMRES could not reach the required relative residual with the ILU preconditioner, and the CGS method showed better calculation performance than the other Krylov subspace methods when ILU was used. TFQMR, QMR, BICGSTAB, and BICG have similar calculation performances. The timestep does not have a significant effect on the calculation performance in terms of the different types of Krylov subspace methods.

Fig. 2

(Color online) Relative residuals for the different Krylov subspace methods with and without the application of the ILU preconditioner at different timesteps

Figure 3 shows the error tolerance (Newton-Raphson tolerance) for the first timestep using different Krylov subspace methods with the ILU preconditioner and different node numbers. The time step was set at 1.0×10^-3 s; the Newton-Raphson error tolerance was set at 1.0×10^-12; the maximum iterations per timestep was set to 1000, and the Krylov subspace method tolerance was set at 1.0×10^-12. The selected Newton-Raphson error tolerance setting was used to investigate the smallest error tolerance that the developed code could reach under these conditions. The error tolerance varies with the Newton-Raphson iteration. As shown in the figure, the error tolerance in the Newton-Raphson iteration is mainly between 1.0×10^-8 and 1.0×10^-9 when the node number is 10, mainly between 1.0×10^-9 and 1.0×10^-10 when the node number is 100, and mainly between 1.0×10^-9 and 1.0×10^-10 when the node number is 1000. The error tolerance could reach a low value in a few steps, regardless of the Krylov subspace method used with the ILU preconditioner. In addition, there is no significant difference in the results obtained with the different Krylov subspace methods in terms of the Newton-Raphson tolerance. This is mainly because the low relative residual could be reached using all the different Krylov subspace methods when the ILU preconditioner was applied.

Fig. 3

(Color online) The error tolerance (Newton-Raphson tolerance) for the first-timestep for different Krylov subspace methods with consideration of ILU preconditioner for different node numbers

Figure 4 shows the iteration (Newton-Raphson iteration) number for each time step during the simulation using different Krylov subspace methods and the ILU preconditioner. The time step was set at 1.0×10^-3 s, the simulation time was set at 0.1 s, the Newton-Raphson error tolerance was set to 1.0×10^-6, the maximum iterations per timestep was set to 100, and the Krylov subspace method tolerance was set to 1.0×10^-6. The Newton-Raphson iteration number is small in the different Krylov subspace methods when the ILU preconditioner is applied. Based on the simulation results, few Newton-Raphson iterations are required to reach the setting tolerance. In addition, the node number does not have a significant effect on the calculation performance of the different Krylov subspace methods examined. However, the number of Newton-Raphson iterations reached 100 when the node number is 100 under the BICG and BICGSTAB methods, indicating that these methods have issues in terms of convergence under some conditions. BICG has an irregular convergence performance [20]. The restarted GMRES method is used by BICGSTAB for the improvement of BICG, but BICGSTAB may have a stagnation issue. This may be the main reason for the convergence problems observed when using BICG and BICGSTAB at 100 nodes.

Fig. 4

(Color online) Iteration (Newton-Raphson iteration) number at each timestep during the simulation using different Krylov subspace methods with the ILU preconditioner. (a) Node number: 10; (b) Node number: 100.

Figure 5 shows the error tolerance (Newton-Raphson iteration) at each time step during the simulation process using different Krylov subspace methods with the ILU preconditioner. The time step was set at 1.0×10^-3 s; the simulation time was set at 0.1 s; the Newton-Raphson error tolerance was set to 1.0×10^-6; the maximum iterations per timestep was set at 100, and the Krylov subspace method tolerance was set to 1.0×10^-6. When the node number is 10, all the Krylov subspace methods could reach the set Newton-Raphson error tolerance requirement with the ILU preconditioner. However, when the node number is 100, the BICG and BICGSTAB methods had an issue with respect to the Newton-Raphson error tolerance, which reaches approximately 5 around 80 timestep when using the two methods. As shown in Figure 4 (b), BICG and BICGSTAB have convergence problems when the node number is 100, which is probably a result of the convergence issues that are inherent in BICG and BICGSTAB.

Fig. 5

(Color online) Error tolerance (Newton-Raphson iteration) at each timestep during the simulation process using different Krylov subspace methods with the ILU preconditioner. (a) Node number: 10; (b) Node number: 100.

3.2 Simulation results

Figure 6 shows a schematic of the flow behavior in the water faucet problem. A vertical pipe is initially filled with both liquid and air. The inner part of the uniform column part is filled with liquid at an initial velocity of 10.0 m/s. Air at an initial velocity of 0.0 m/s surrounds the liquid. The initial void fraction of the liquid is 0.8. At the top of the pipe, the liquid inlet velocity is 10.0 m/s, and the air inlet velocity is 0.0 m/s. The bottom of the pipe is left open at a constant pressure. The liquid rapidly starts to fall under gravity. The acceleration in the flow that results from gravity leads to discontinuity at a particular location in the falling process, as can be observed in Fig. 6, after which the system eventually reaches a steady state. The theoretical solutions for the liquid velocity and gas void fraction were first proposed by Ranson [14]. Zou et al. extended Ranson’s work and provided theoretical solutions for the gas velocity and pressure distribution [15].

Fig. 6

Schematic of the flow behavior in the water faucet problem [15]

The discontinuity location can be expressed as

x_{d} = \sqrt{u_{f,int}^{2} + 2 g x} .

(12)

The liquid velocity can be expressed as

u_{f} (x, t) = {\begin{array}{l} \sqrt{u_{f,int}^{2} + 2 g x} & if x \leq x_{d} \\ u_{f,int} + g t & otherwise \end{array} .

(13)

The gas velocity can be expressed as

u_{g} (x, t) = {\begin{array}{l} 0 & if x \leq x_{d} \\ - \frac{1 - α_{g,int}}{α_{g,int}} & otherwise \end{array} .

(14)

The gas void fraction can be expressed as

α_{g} (x, t) = {\begin{array}{l} 1 - \frac{(1 - α_{g,int}) u_{l,int}}{\sqrt{u_{f,int}^{2} + 2 g x}} & if x \leq x_{d} \\ α_{g,int} & otherwise \end{array} .

(15)

The pressure distribution can be expressed as

p (x, t) = {\begin{matrix} \begin{array}{l} p_{outlet} - \frac{ρ_{g} g (L - x)}{a_{g,int}} + \frac{1}{2} ρ_{g} {(\frac{1 - α_{g,int}}{α_{g,int}} g t)}^{2} + \\ ρ_{g} \frac{1 - α_{g,int}}{α_{g,int}} g t (u_{l,int} + g t) - ρ_{g} g (x_{d} - x) \end{array} & if x \leq x_{d} \\ p_{outlet} - \frac{ρ_{g} g (L - x)}{a_{g,int}} & otherwise \end{matrix} .

(16)

Figure 7 shows the flow behavior at 0.1 s for the different Krylov subspace methods with the ILU preconditioner. The node number was set at 100; the time step was set at 1.0×10^-3 s; the simulation time was set to 0.1 s; the Newton-Raphson error tolerance was set to 1.0×10^-6, the maximum iterations per timestep was set at 100, and the Krylov subspace method tolerance was set to 1.0×10^-6. Based on the simulation results, regardless of the Krylov subspace method, the simulation results are generally consistent with the analytical solution. The discontinuous behavior was reasonably captured by the simulation code. This is mainly because the application of the ILU preconditioner means that the error tolerance (Newton-Raphson iteration) can reach a low value in most cases, as shown in Figure 5. Although the relative residual of the GMRES with ILU could not reach 1.0×10^-12 as shown in Figure 1, it could reach approximately 1.0×10^-10, which is less than the current tolerance setting used with the Krylov subspace method, and it can therefore be considered able to capture the water faucet behavior based on the simulation results.

Fig. 7

(Color online) Flow behavior at 0.1 s for the different Krylov subspace methods with the ILU preconditioner. (a) Liquid velocity; (b) Gas velocity; (c) Gas void fraction; (d) Pressure.

For the liquid flow behavior, the inlet liquid velocity was fixed at 10.0 m/s in the simulation settings. In the analytical solution, the liquid velocity increases from 10.0 m/s to approximately 11.0 m/s as the liquid flows to a position 1.0 m away from the inlet. The velocity of the liquid thereafter remained the same. The liquid velocity distribution obtained by the numerical simulation was consistent with the analytical solution. The liquid void fraction in the simulation domain decreases over time before steady state is achieved, while the gas needs to flow into the simulation domain to fill the space. Because the gas velocity at the inlet position is set to 0.0 m/s, gas flows into the simulation domain from the outlet position. In the analytical solution, the gas velocity remained at 0.0 m/s from the inlet position to a distance of approximately 1.0 m and reached approximately -4.0 m/s at distances greater than approximately 1.0 m from the inlet. There is an obvious discontinuity in the gas velocity distribution at approximately 1.0 m. The gas velocity obtained by the simulation changes smoothly from 0.0 m/s to approximately -4.0 m/s at around 1.0 m from the inlet. These simulation results are mainly due to numerical diffusion. In the analytical solution, the gas void fraction increases from 0.2 at the inlet to approximately 0.27 at a position approximately 1.0 m from the inlet. The gas void fraction is 0.2, and this remains the same at distances greater than 1.0 m. The maximum gas void fraction is approximately 0.27, which occurs at the discontinuity. In the numerical simulation, the gas void fraction is also 0.2 at distances greater than 1.0 m. The maximum gas void fraction obtained by the numerical simulation is between 0.24 and 0.25, which is smaller than that found via analysis. There is also an obvious discontinuity in the pressure distribution, which lies approximately 1.0 m from the inlet. The pressure distribution along the distance obtained by numerical simulation matches the analytical solution reasonably well.

Figure 8 shows the flow behavior at 0.1 s for different Krylov subspace methods without the ILU preconditioner. The node number was set to 100; the timestep was set at 1.0×10^-3 s, and the simulation time was set to 0.1 s. The Newton-Raphson error tolerance was set to 1.0×10^-6; the maximum iterations per timestep was set at 100, and the Krylov subspace method tolerance was set to 1.0×10^-6. Based on the simulation results, the developed code was found to capture the liquid velocity and liquid void fraction regardless of the Krylov subspace method used. However, except for the GMRES method, all the Krylov subspace methods have issues in capturing the gas velocity and pressure distribution. This is mainly because the liquid in the system is mainly driven by gravity, and the pressure difference does not have a significant effect on the flow behavior of the liquid. However, the gas velocity is based on the conditions under which the liquid flows and the pressure distribution. If the code cannot capture the pressure distribution properly, then it also has problems in terms of capturing the gas velocity. The GMRES could capture the pressure distribution and gas velocity accurately; this is mainly because it could also obtain a small relative residual even without the ILU preconditioner, as shown in Figs. 1 and 2. However, the other Krylov subspace methods used in this study could not reach a small relative residual without the ILU preconditioner, as shown in Figs. 1 and 2.

Fig. 8

(Color online) Flow behavior at 0.1 s for different Krylov subspace methods without the ILU preconditioner. (a) Liquid velocity; (b) Gas velocity; (c) Gas void fraction; (d) Pressure.

Figure 9 shows the flow behavior at 0.5 s for the different Krylov subspace methods with the ILU preconditioner. The node number was set at 100; the timestep was set at 1.0×10^-3 s; the Newton-Raphson error tolerance was set to 1.0×10^-6; the maximum iterations per timestep was set to 100, and the Krylov subspace method tolerance was set at 1.0×10^-6. Based on the simulation results, the velocity flow behavior of both liquid and gas can be captured, regardless of the Krylov subspace method used. Compared to the flow behavior at 0.2 s, the discontinuous point propagates to approximately 6.0 m at 0.5 s. In the analytical solution, the liquid velocity is approximately 15 m/s following a discontinuity after increasing from 10.0 m/s to approximately 15.0 m/s. The liquid velocity obtained using the numerical solution matched that obtained using the analytical solution. In the analytical solution, the gas velocity is approximately -20 m/s after the discontinuity. The gas velocity distribution obtained by numerical simulation therefore matches that found via the analytical solution well, except around the point of discontinuity. The diffusion behavior caused by the numerical method around the discontinuity can also be observed in the gas velocity distribution. The analytical solution suggests that the maximum gas void fraction occurs at the point of discontinuity, with a value of approximately 0.45. The maximum gas void fraction in the numerical simulation is at around 5.0 m, which exceeds the discontinuous point predicted by the analytical solution. The maximum gas void fraction obtained by numerical simulation is less than that predicted by the analytical solution, and this is mainly caused by the numerical diffusion in the simulation; however, errors may have occurred in the pressure distribution captured before the discontinuous point. The main reason for this is that the single-pressure model is used in the simulation code. Different pressure conditions were also considered for the liquid and gas in the analytical solution.

Fig. 9

(Color online) Flow behavior at 0.5 s for different Krylov subspace methods with ILU preconditioner. (a) Liquid velocity; (b) Gas velocity; (c) Gas void fraction; (d) Pressure.

Figure 10 shows the flow behavior at 1.2 s for the different Krylov subspace methods with the ILU preconditioner. The node number was set at 100, and the timestep was set at 1.0×10^-3 s. The Newton-Raphson error tolerance was set to 1.0×10^-6; the maximum iterations per timestep was set at 100, and the Krylov subspace method tolerance was set to 1.0×10^-6. Regardless of the Krylov subspace method used, both the liquid and the gas velocity flow behavior could be captured. The analytical solution indicates that the liquid velocity increases from 10.0 m/s to approximately 18.5 m/s along the distance travelled. There is no discontinuity in the liquid velocity along the distance in a steady state. The liquid velocity obtained by the numerical simulation matches well with the analytical solution, which gives a gas velocity of 0.0 m/s along the distance. This is mainly because the gas void fraction does not change over time, and gas injection was not necessary. A gas velocity of approximately 0.0 m/s was also obtained along the distance by numerical simulation. The maximum absolute gas velocity is approximately 4.0×10^-3 m/s, which is very close to that obtained via the analytical solution. The gas void fraction increases along the distance from 0.2 to more than 0.55 in the analytical solution. The gas void fraction obtained by the numerical simulation matches the analytical solution well. The pressure increases linearly along the distance in the analytical solution, with a similar result obtained via numerical simulation. However, there are differences between the results of the analytical solution and the numerical simulation. The main reason for this is that the two-pressure model was applied in the analytical solution, whereas the single-pressure model was applied in the numerical simulation.

Fig. 10

(Color online) Flow behavior at 1.2 s using different Krylov subspace methods with the ILU preconditioner. (a) Liquid velocity; (b) Gas velocity; (c) Gas void fraction; (d) Pressure.

3.3 Simulation time

Figure 11 shows the computational time taken by GMRES with and without application of the ILU preconditioner. The timestep was set at 1.0×10^-3 s. The total simulation time was 0.1 s. The Newton-Raphson error tolerance was set to 1.0×10^-6; the maximum iterations per timestep was set at 100, and the Krylov subspace method tolerance was set to 1.0×10^-6. As shown in Fig. 8, only the GMRES without the ILU preconditioner could sufficiently capture the pressure distribution and gas velocity. Therefore, the GMRES method was selected for the comparison of computational time with and without the application of the ILU preconditioner. Based on the simulation results, there is a significant difference in the computational time taken for cases with and without the application of the ILU preconditioner, especially when the node number is large. When the node number is less than 50, there is no significant difference between the two conditions. However, increasing the number of nodes led to an increase in the computational time. This is mainly because without the ILU preconditioner, the GMRES solver requires a large number of iterations to reach or become close to the set tolerance, while only few steps are required to reach or become close to the setting tolerance when the ILU preconditioner is used with the GMRES solver. In addition, the computational time taken without the ILU preconditioner has a nearly cubic relationship with the node number. The computational time for the GMRES with the ILU preconditioner is therefore much smaller than that taken when the ILU preconditioner was not applied, for which a large number of nodes is required.

Fig. 11

(Color online) Computational time for GMRES with and without application of the ILU preconditioner

Figure 12 shows the computational time taken for different Krylov subspace methods with the ILU preconditioner. The time step was set at 1.0×10^-3 s, and the total simulation time was set at 0.1 s. The Newton-Raphson error tolerance was set as 1.0×10^-6; the maximum iterations per timestep was set to 100, and the Krylov subspace method tolerance was set to 1.0×10^-6. The node numbers were set at 25, 50, 100, 200, and 300 for the different tests. The CPU calculation time taken for cases with 25 nodes was similar to that with 50 nodes. The CPU calculation time increased almost linearly when the number of nodes increased from 50 to 200. Based on the simulation time, there is no significant difference in the computational time taken by the different Krylov subspace methods when the ILU preconditioner was used. This is mainly because the Krylov subspace methods did not require many iterations to reach or approach the required tolerance when the ILU preconditioner was used, as shown in Figs. 1 and 2. CGS and BICGSTAB are slightly better than other Krylov subspace methods in terms of efficiency.

Fig. 12

(Color online) Computational time taken by the different Krylov subspace methods under application of ILU preconditioner

Table 1 summarizes the simulated performances of the six Krylov subspace methods in terms of relative residuals, simulation results, and computational efficiency. The symbol ‘+’ in the table indicates that the simulated performance is acceptable for the relative residual and the results of the water faucet problem simulation. The symbol ‘+’ in the computational efficiency column indicates the relative efficiency performance. For the relative residual, all six Krylov subspace methods with the ILU preconditioner showed a relatively good performance, whereas only GMRES could provide a relatively good simulation performance without the ILU preconditioner. Although the relative residual obtained with GMRES could not reach the level of 1.0×10^-12, it could reach an acceptable level for the water faucet problem simulation. Regarding the results of the water faucet problem simulation, all six Krylov subspace methods were found able to capture the liquid velocity and gas void fraction with an ILU preconditioner, whereas only GMRES could capture the gas velocity and pressure without an ILU preconditioner. Regarding the computational efficiency, the time cost of using GMRES without an ILU preconditioner was much greater than that using GMRES with an ILU preconditioner. All six Krylov subspace methods demonstrated a similar performance with an ILU preconditioner in terms of computational efficiency.

Simulated performances of different Krylov subspace methods

		Relative residual	Water faucet problem simulation results				Computational efficiency
			Liquid velocity	Gas velocity	Gas void fraction	Pressure
GMRES	With ILU	+	+	+	+	+	+
GMRES	Without ILU	+	+	+	+	+	-
TFQMR	With ILU	+	+	+	+	+	+
TFQMR	Without ILU	-	+	-	+	-	N/A
QMR	With ILU	+	+	+	+	+	+
QMR	Without ILU	-	+	-	+	-	N/A
CGS	With ILU	+	+	+	+	+	+
CGS	Without ILU	-	+	-	+	-	N/A
BICGSTAB	With ILU	+	+	+	+	+	+
BICGSTAB	Without ILU	-	+	-	+	-	N/A
BICG	With ILU	+	+	+	+	+	+
BICG	Without ILU	-	+	-	+	-	N/A

4 Summary and future work

In this study, a one-dimensional two-phase flow four-equation model was developed to simulate the water faucet problem. An evaluation of different Krylov subspace methods with and without the application of an ILU preconditioner was conducted, and the simulation results were obtained and analyzed.

(1) The convergence performance of Krylov subspace methods with an ILU preconditioner is much better than that of the Krylov subspace methods without a preconditioner. The analysis of the relative residuals for different Krylov subspace methods was performed with and without the application of the ILU preconditioner; with different node numbers; and the relative residual for different Krylov subspace methods with and without consideration of the ILU preconditioner at different time steps. Based on the analysis results, GMRES was found to have better calculation performance than the other Krylov subspace methods when the ILU preconditioner was not applied. The CGS method had a better calculation performance than the other Krylov subspace methods when the ILU preconditioner was used, but there was no significant difference among the Krylov subspace methods.

(2) Based on the simulation results and comparison with the analytical solution, the Krylov subspace methods used in this study could capture the water faucet behavior in terms of liquid velocity, gas velocity, gas void fraction, and pressure distribution when the ILU preconditioner was applied. Differences between the analytical solution and numerical simulation were observed in terms of predicating the pressure distribution. This is mainly because the single-pressure model was used in the numerical simulation. All the Krylov subspace methods except GMRES could capture the liquid velocity and gas void fraction without the ILU preconditioner; however, they could not capture the behavior in terms of gas velocity and pressure distribution. GMRES could capture both the gas velocity and the pressure distribution without the need for an ILU preconditioner.

(3) The computational time taken by GMRES with the ILU preconditioner was much less than that taken when the ILU preconditioner was not applied. There was no significant difference in the computational time taken by the different Krylov subspace methods with the ILU preconditioner. The CPU calculation times taken by CGS and BICGSTAB were slightly shorter than the CPU calculation time taken under the other Krylov subspace methods.

In the future, evaluating the Krylov subspace methods with the two-pressure model and the energy conservation equations in the simulation of two-phase flow and heat transfer behavior would be beneficial to better understand their application in terms of performance.

References:

[1]

G. A. Roth, F. Aydogan,

Theory and implementation of nuclear safety system codes-Part I: Conservation equations, flow regimes, numerics and significant assumptions

. Prog. Nucl. Energy 76, 160-182 (2014). doi: 10.1016/j.pnucene.2014.05.011