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

Evaluation of different Krylov subspace methods for simulation of the water faucet problem 1 corresp first-author Corresponding author, hongyang.wei@unlv.edu Corresponding author, hongyang.wei@unlv.edu hongyang.wei@unlv.edu https://orcid.org/0000-0002-5359-2121 1 1 1 Department of Mechanical Engineering, University of Nevada, Las Vegas, NV89154, U.S.A. Department of Mechanical Engineering, University of Nevada, Las Vegas, NV89154, U.S.A. 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 problem Krylov subspace methods ILU preconditioner 1 Introduction 1 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 1 , 2 2 , 3 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 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 1 , 5 5 ]. The nonlinear numerical solver JFNK is used for simulations that involve RELAP7 [ 1 1 , 6 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 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 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 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 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 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 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 13 ]. The water faucet problem was first proposed by Ranson for the investigation of two-phase flow behavior [ 14 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 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 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 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 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 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 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 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 1 2.1 Governing equations 2 The governing equations that are used to simulate water faucet behavior are the simplifications of RELAP5’s mass and momentum conservation equations [ 22 22 ]. The terms describing phase change, interfacial heat transfer, a nd interfacial friction are ignored in tests associated with the water faucet problem. The gas mass conservation equation is expressed as: <math display="block" id="M1"><mrow><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow><mo>+</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub><msub><mi>u</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow><mo>=</mo><mn>0</mn></mrow></math> The liquid mass conservation equation is expressed as: <math display="block" id="M2"><mrow><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>f</mtext></msub><msub><mi>ρ</mi><mtext>f</mtext></msub></mrow><mo>)</mo></mrow><mo>+</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>f</mtext></msub><msub><mi>ρ</mi><mtext>f</mtext></msub><msub><mi>u</mi><mtext>f</mtext></msub></mrow><mo>)</mo></mrow><mo>=</mo><mn>0.</mn></mrow></math> The gas momentum conservation equation is expressed as: <math display="block" id="M3"><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub><mfrac><mrow><mo>∂</mo><msub><mi>u</mi><mtext>g</mtext></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub><mfrac><mrow><mo>∂</mo><msub><mi>u</mi><mtext>g</mtext></msub><msup><mrow/><mn>2</mn></msup></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>=</mo><mo>−</mo><msub><mi>α</mi><mtext>g</mtext></msub><mfrac><mrow><mo>∂</mo><mi>P</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>+</mo><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub><msub><mi>g</mi><mi>x</mi></msub><mo>.</mo></mrow></math> The liquid momentum conservation equation is expressed as: <math display="block" id="M4"><mrow><msub><mi>α</mi><mtext>f</mtext></msub><msub><mi>ρ</mi><mtext>f</mtext></msub><mfrac><mrow><mo>∂</mo><msub><mi>u</mi><mtext>f</mtext></msub></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><msub><mi>α</mi><mtext>f</mtext></msub><msub><mi>ρ</mi><mtext>f</mtext></msub><mfrac><mrow><mo>∂</mo><msub><mi>u</mi><mtext>f</mtext></msub><msup><mrow/><mn>2</mn></msup></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>=</mo><mo>−</mo><msub><mi>α</mi><mtext>f</mtext></msub><mfrac><mrow><mo>∂</mo><mi>P</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>+</mo><msub><mi>α</mi><mtext>f</mtext></msub><msub><mi>ρ</mi><mtext>f</mtext></msub><msub><mi>g</mi><mi>x</mi></msub><mo>.</mo></mrow></math> 2.2 Discretization methods 2 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 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 <math id="M5"><mrow><msub><mi>u</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>></mo><mn>0</mn></mrow></math> , <math display="block" id="M6"><mrow><mfrac><mrow><msub><mi>α</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi></mrow></msub><msub><mi>ρ</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi></mrow></msub><mo>−</mo><msubsup><mi>α</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi></mrow><mi>n</mi></msubsup><msubsup><mi>ρ</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi></mrow><mi>n</mi></msubsup></mrow><mrow><mi>Δ</mi><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><msub><mi>α</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi></mrow></msub><msub><mi>ρ</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi></mrow></msub><msub><mi>u</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>−</mo><msub><mi>α</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi><mo>−</mo><mn>1</mn></mrow></msub><msub><mi>ρ</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi><mo>−</mo><mn>1</mn></mrow></msub><msub><mi>u</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mrow><mi>Δ</mi><mi>x</mi></mrow></mfrac><mo>=</mo><mn>0.</mn></mrow></math> (2) If <math id="M7"><mrow><msub><mi>u</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>></mo><mn>0</mn></mrow></math> , <math display="block" id="M8"><mrow><mfrac><mrow><msub><mi>α</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi></mrow></msub><msub><mi>ρ</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi></mrow></msub><mo>−</mo><msubsup><mi>α</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi></mrow><mi>n</mi></msubsup><msubsup><mi>ρ</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi></mrow><mi>n</mi></msubsup></mrow><mrow><mi>Δ</mi><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><msub><mi>α</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi><mo>+</mo><mn>1</mn></mrow></msub><msub><mi>ρ</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi><mo>+</mo><mn>1</mn></mrow></msub><msub><mi>u</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>−</mo><msub><mi>α</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi></mrow></msub><msub><mi>ρ</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi></mrow></msub><msub><mi>u</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mrow><mi>Δ</mi><mi>x</mi></mrow></mfrac><mo>=</mo><mn>0.</mn></mrow></math> The liquid mass conservation equation can be discretized as follows: <math display="block" id="M9"><mrow><mfrac><mrow><msub><mi>α</mi><mrow><mtext>f</mtext><mo>,</mo><mi>j</mi></mrow></msub><msub><mi>ρ</mi><mrow><mtext>f</mtext><mo>,</mo><mi>j</mi></mrow></msub><mo>−</mo><msubsup><mi>α</mi><mrow><mtext>f</mtext><mo>,</mo><mi>j</mi></mrow><mi>n</mi></msubsup><msubsup><mi>ρ</mi><mrow><mtext>f</mtext><mo>,</mo><mi>j</mi></mrow><mi>n</mi></msubsup></mrow><mrow><mi>Δ</mi><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><msub><mi>α</mi><mrow><mtext>f</mtext><mo>,</mo><mi>j</mi></mrow></msub><msub><mi>ρ</mi><mrow><mtext>f</mtext><mo>,</mo><mi>j</mi></mrow></msub><msub><mi>u</mi><mrow><mtext>f</mtext><mo>,</mo><mi>j</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>−</mo><msub><mi>α</mi><mrow><mtext>f</mtext><mo>,</mo><mi>j</mi><mo>−</mo><mn>1</mn></mrow></msub><msub><mi>ρ</mi><mrow><mtext>f</mtext><mo>,</mo><mi>j</mi><mo>−</mo><mn>1</mn></mrow></msub><msub><mi>u</mi><mrow><mtext>f</mtext><mo>,</mo><mi>j</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mrow><mi>Δ</mi><mi>x</mi></mrow></mfrac><mo>=</mo><mn>0.</mn></mrow></math> The gas momentum conservation equation can be discretized as follows: (1) If <math id="M10"><mrow><msub><mi>u</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>></mo><mn>0</mn></mrow></math> , <math display="block" id="M11"><mtable><mtr><mtd><mfrac><mrow><mover accent="true"><mrow><msub><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>J</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo stretchy="true">¯</mo></mover><msub><mi>u</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>−</mo><mover accent="true"><mrow><msub><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>J</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo stretchy="true">¯</mo></mover><msup><mrow><mrow><mo>(</mo><mrow><msub><mi>u</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo>)</mo></mrow></mrow><mi>n</mi></msup></mrow><mrow><mi>Δ</mi><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><msub><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow><mi>j</mi></msub><mfrac><mrow><msup><mrow><mrow><mo>(</mo><mrow><mover accent="true"><mrow><msub><mi>u</mi><mrow><mtext>g</mtext><mo>,</mo><mi>J</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow><mo stretchy="true">¯</mo></mover></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mover accent="true"><mrow><msub><mi>u</mi><mrow><mtext>g</mtext><mo>,</mo><mi>J</mi></mrow></msub></mrow><mo stretchy="true">¯</mo></mover></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mrow><mi>Δ</mi><mi>x</mi></mrow></mfrac></mtd></mtr><mtr><mtd><mo>=</mo><mo>−</mo><mover accent="true"><mrow><msub><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>J</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo stretchy="true">¯</mo></mover><mfrac><mrow><msub><mi>P</mi><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mi>P</mi><mi>j</mi></msub></mrow><mrow><mi>Δ</mi><mi>x</mi></mrow></mfrac><mo>+</mo><mover accent="true"><mrow><msub><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>J</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo stretchy="true">¯</mo></mover><msub><mi>g</mi><mi>x</mi></msub><mo>.</mo></mtd></mtr></mtable></math> As an example <math id="M12"><mrow><mover accent="true"><mrow><msub><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>J</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo stretchy="true">¯</mo></mover></mrow></math> , can be expressed as <math id="M13"><mrow><mfrac><mrow><mover accent="true"><mrow><msub><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow><mi>J</mi></msub></mrow><mo stretchy="true">¯</mo></mover><mo>+</mo><mover accent="true"><mrow><msub><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>J</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow><mo stretchy="true">¯</mo></mover></mrow><mn>2</mn></mfrac></mrow></math> . (2) If <math id="M14"><mrow><msub><mi>u</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>≤</mo><mn>0</mn></mrow></math> , <math display="block" id="M15"><mtable><mtr><mtd><mfrac><mrow><mover accent="true"><mrow><msub><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>J</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo stretchy="true">¯</mo></mover><msub><mi>u</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>−</mo><mover accent="true"><mrow><msub><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>J</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo stretchy="true">¯</mo></mover><msup><mrow><mrow><mo>(</mo><mrow><msub><mi>u</mi><mrow><mtext>g</mtext><mo>,</mo><mi>j</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo>)</mo></mrow></mrow><mi>n</mi></msup></mrow><mrow><mi>Δ</mi><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><msub><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow></msub><mfrac><mrow><msup><mrow><mrow><mo>(</mo><mrow><mover accent="true"><mrow><msub><mi>u</mi><mrow><mtext>g</mtext><mo>,</mo><mi>J</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow><mo stretchy="true">¯</mo></mover></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mover accent="true"><mrow><msub><mi>u</mi><mrow><mtext>g</mtext><mo>,</mo><mi>J</mi></mrow></msub></mrow><mo stretchy="true">¯</mo></mover></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mrow><mi>Δ</mi><mi>x</mi></mrow></mfrac></mtd></mtr><mtr><mtd><mo>=</mo><mo>−</mo><mover accent="true"><mrow><msub><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>J</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo stretchy="true">¯</mo></mover><mfrac><mrow><msub><mi>P</mi><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mi>P</mi><mi>j</mi></msub></mrow><mrow><mi>Δ</mi><mi>x</mi></mrow></mfrac><mo>+</mo><mover accent="true"><mrow><msub><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>g</mtext></msub><msub><mi>ρ</mi><mtext>g</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>J</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo stretchy="true">¯</mo></mover><msub><mi>g</mi><mi>x</mi></msub><mo>.</mo></mtd></mtr></mtable></math> The liquid momentum conservation equation can be discretized as follows: <math display="block" id="M16"><mtable><mtr><mtd><mfrac><mrow><mover accent="true"><mrow><msub><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>f</mtext></msub><msub><mi>ρ</mi><mtext>f</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>J</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo stretchy="true">¯</mo></mover><msub><mi>u</mi><mrow><mtext>f</mtext><mo>,</mo><mi>j</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>−</mo><mover accent="true"><mrow><msub><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>f</mtext></msub><msub><mi>ρ</mi><mtext>f</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>J</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo stretchy="true">¯</mo></mover><msup><mrow><mrow><mo>(</mo><mrow><msub><mi>u</mi><mrow><mtext>f</mtext><mo>,</mo><mi>j</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo>)</mo></mrow></mrow><mi>n</mi></msup></mrow><mrow><mi>Δ</mi><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><msub><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>f</mtext></msub><msub><mi>ρ</mi><mtext>f</mtext></msub></mrow><mo>)</mo></mrow><mi>j</mi></msub><mfrac><mrow><msup><mrow><mrow><mo>(</mo><mrow><mover accent="true"><mrow><msub><mi>u</mi><mrow><mtext>f</mtext><mo>,</mo><mi>J</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow><mo stretchy="true">¯</mo></mover></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mover accent="true"><mrow><msub><mi>u</mi><mrow><mtext>f</mtext><mo>,</mo><mi>J</mi></mrow></msub></mrow><mo stretchy="true">¯</mo></mover></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mrow><mi>Δ</mi><mi>x</mi></mrow></mfrac></mtd></mtr><mtr><mtd><mo>=</mo><mo>−</mo><mover accent="true"><mrow><msub><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>f</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>J</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo stretchy="true">¯</mo></mover><mfrac><mrow><msub><mi>P</mi><mrow><mi>j</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mi>P</mi><mi>j</mi></msub></mrow><mrow><mi>Δ</mi><mi>x</mi></mrow></mfrac><mo>+</mo><mover accent="true"><mrow><msub><mrow><mrow><mo>(</mo><mrow><msub><mi>α</mi><mtext>f</mtext></msub><msub><mi>ρ</mi><mtext>f</mtext></msub></mrow><mo>)</mo></mrow></mrow><mrow><mi>J</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo stretchy="true">¯</mo></mover><msub><mi>g</mi><mi>x</mi></msub><mo>.</mo></mtd></mtr></mtable></math> 2.3 Krylov subspace method [ 23 23 , 24 24 ] 2 The Petrov-Galerkin condition can be applied to obtain the solution for the linear system <math id="M17"><mrow><mi>A</mi><mi>x</mi><mo>=</mo><mi>b</mi></mrow></math> via the following expression: <math display="block" id="M18"><mrow><mi>b</mi><mo>−</mo><mi>A</mi><msub><mi>x</mi><mi>m</mi></msub><mo>⊥</mo><msub><mi mathvariant="script">L</mi><mi>m</mi></msub><mo>,</mo></mrow></math> Where <math id="M19"><mrow><msub><mi>x</mi><mi>m</mi></msub></mrow></math> is in the subspace <math id="M20"><mrow><msub><mi>x</mi><mn>0</mn></msub><mo>+</mo><msub><mi mathvariant="script">K</mi><mi>m</mi></msub></mrow></math> with m dimensions, and <math id="M21"><mrow><msub><mi>x</mi><mn>0</mn></msub></mrow></math> is the initial guess. The method is considered to be a Krylov subspace method when subspace <math id="M22"><mrow><msub><mi mathvariant="script">K</mi><mi>m</mi></msub></mrow></math> is the Krylov subspace. Generally, there are two types of Krylov subspace methods: one in which <math id="M23"><mrow><msub><mi mathvariant="script">L</mi><mi>m</mi></msub><mo>=</mo><msub><mi mathvariant="script">K</mi><mi>m</mi></msub></mrow></math> , and the other where <math id="M24"><mrow><msub><mi mathvariant="script">L</mi><mi>m</mi></msub><mo>=</mo><msub><mi mathvariant="script">K</mi><mi>m</mi></msub><mo stretchy="false">(</mo><msup><mi>A</mi><mi>T</mi></msup><mo>,</mo><msub><mi>r</mi><mn>0</mn></msub><mo stretchy="false">)</mo></mrow></math> . 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 1 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 2 Figure 1 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 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. Figure 2 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. Figure 3 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. Figure 4 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 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. Figure 5 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) 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. 3.2 Simulation results 2 Figure 6 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 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 14 ]. Zou et al. extended Ranson’s work and provided theoretical solutions for the gas velocity and pressure distribution [ 15 15 ]. The discontinuity location can be expressed as <math display="block" id="M25"><mrow><msub><mi>x</mi><mi>d</mi></msub><mo>=</mo><msqrt><mrow><msubsup><mi>u</mi><mrow><mtext>f,int</mtext></mrow><mn>2</mn></msubsup><mo>+</mo><mn>2</mn><mi>g</mi><mi>x</mi></mrow></msqrt><mo>.</mo></mrow></math> The liquid velocity can be expressed as <math display="block" id="M26"><mrow><msub><mi>u</mi><mtext>f</mtext></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><msqrt><mrow><msubsup><mi>u</mi><mrow><mtext>f,int</mtext></mrow><mn>2</mn></msubsup><mo>+</mo><mn>2</mn><mi>g</mi><mi>x</mi></mrow></msqrt></mrow></mtd><mtd columnalign="left"><mrow><mtext>if</mtext><mi>x</mi><mo>≤</mo><msub><mi>x</mi><mtext>d</mtext></msub></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>u</mi><mrow><mtext>f,int</mtext></mrow></msub><mo>+</mo><mi>g</mi><mi>t</mi></mrow></mtd><mtd columnalign="left"><mrow><mtext>otherwise</mtext></mrow></mtd></mtr></mtable></mrow></mrow><mo>.</mo></mrow></math> The gas velocity can be expressed as <math display="block" id="M27"><mrow><msub><mi>u</mi><mtext>g</mtext></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mn>0</mn></mtd><mtd columnalign="left"><mrow><mtext>if</mtext><mi>x</mi><mo>≤</mo><msub><mi>x</mi><mtext>d</mtext></msub></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><mo>−</mo><mfrac><mrow><mn>1</mn><mo>−</mo><msub><mi>α</mi><mrow><mtext>g,int</mtext></mrow></msub></mrow><mrow><msub><mi>α</mi><mrow><mtext>g,int</mtext></mrow></msub></mrow></mfrac></mrow></mtd><mtd columnalign="left"><mrow><mtext>otherwise</mtext></mrow></mtd></mtr></mtable></mrow></mrow><mo>.</mo></mrow></math> The gas void fraction can be expressed as <math display="block" id="M28"><mrow><msub><mi>α</mi><mtext>g</mtext></msub><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mn>1</mn><mo>−</mo><mfrac><mrow><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><msub><mi>α</mi><mrow><mtext>g,int</mtext></mrow></msub><mo stretchy="false">)</mo><msub><mi>u</mi><mrow><mtext>l,int</mtext></mrow></msub></mrow><mrow><msqrt><mrow><msubsup><mi>u</mi><mrow><mtext>f,int</mtext></mrow><mn>2</mn></msubsup><mo>+</mo><mn>2</mn><mi>g</mi><mi>x</mi></mrow></msqrt></mrow></mfrac></mrow></mtd><mtd columnalign="left"><mrow><mtext>if</mtext><mi>x</mi><mo>≤</mo><msub><mi>x</mi><mtext>d</mtext></msub></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>α</mi><mrow><mtext>g,int</mtext></mrow></msub></mrow></mtd><mtd columnalign="left"><mrow><mtext>otherwise</mtext></mrow></mtd></mtr></mtable></mrow></mrow><mo>.</mo></mrow></math> The pressure distribution can be expressed as <math display="block" id="M29"><mrow><mi>p</mi><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mtable columnalign="left"><mtr><mtd><msub><mi>p</mi><mrow><mtext>outlet</mtext></mrow></msub><mo>−</mo><mfrac><mrow><msub><mi>ρ</mi><mtext>g</mtext></msub><mi>g</mi><mo stretchy="false">(</mo><mi>L</mi><mo>−</mo><mi>x</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mi>a</mi><mrow><mtext>g,int</mtext></mrow></msub></mrow></mfrac><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><msub><mi>ρ</mi><mtext>g</mtext></msub><msup><mrow><mo>(</mo><mrow><mfrac><mrow><mn>1</mn><mo>−</mo><msub><mi>α</mi><mrow><mtext>g,int</mtext></mrow></msub></mrow><mrow><msub><mi>α</mi><mrow><mtext>g,int</mtext></mrow></msub></mrow></mfrac><mi>g</mi><mi>t</mi></mrow><mo>)</mo></mrow><mn>2</mn></msup><mo>+</mo></mtd></mtr><mtr><mtd><msub><mi>ρ</mi><mtext>g</mtext></msub><mfrac><mrow><mn>1</mn><mo>−</mo><msub><mi>α</mi><mrow><mtext>g,int</mtext></mrow></msub></mrow><mrow><msub><mi>α</mi><mrow><mtext>g,int</mtext></mrow></msub></mrow></mfrac><mi>g</mi><mi>t</mi><mo stretchy="false">(</mo><msub><mi>u</mi><mrow><mtext>l,int</mtext></mrow></msub><mo>+</mo><mi>g</mi><mi>t</mi><mo stretchy="false">)</mo><mo>−</mo><msub><mi>ρ</mi><mtext>g</mtext></msub><mi>g</mi><mo stretchy="false">(</mo><msub><mi>x</mi><mtext>d</mtext></msub><mo>−</mo><mi>x</mi><mo stretchy="false">)</mo></mtd></mtr></mtable></mtd><mtd><mrow><mtext>if</mtext><mi>x</mi><mo>≤</mo><msub><mi>x</mi><mtext>d</mtext></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>p</mi><mrow><mtext>outlet</mtext></mrow></msub><mo>−</mo><mfrac><mrow><msub><mi>ρ</mi><mtext>g</mtext></msub><mi>g</mi><mo stretchy="false">(</mo><mi>L</mi><mo>−</mo><mi>x</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mi>a</mi><mrow><mtext>g,int</mtext></mrow></msub></mrow></mfrac></mrow></mtd><mtd><mrow><mtext>otherwise</mtext></mrow></mtd></mtr></mtable></mrow></mrow><mo>.</mo></mrow></math> Figure 7 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 Figure 5 . Although the relative residual of the GMRES with ILU could not reach 1.0×10 -12 as shown in Figure 1 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. 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 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 Figs. 1 and 2 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 Figs. 1 and 2 2 . Figure 9 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. Figure 10 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. 3.3 Simulation time 2 Figure 11 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 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. Figure 12 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 Figs. 1 and 2 2 . CGS and BICGSTAB are slightly better than other Krylov subspace methods in terms of efficiency. Table 1 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. 1 1   1 1   1 1 Relative residual 4 1 Water faucet problem simulation results 1 1 Computational efficiency  1 1   1 1   1 1   1 1 Liquid velocity 1 1 Gas velocity 1 1 Gas void fraction 1 1 Pressure 1 1   1 2 GMRES 1 1 With ILU 1 1 + 1 1 + 1 1 + 1 1 + 1 1 + 1 1 + 1 1 Without ILU 1 1 + 1 1 + 1 1 + 1 1 + 1 1 + 1 1 - 1 2 TFQMR 1 1 With ILU 1 1 + 1 1 + 1 1 + 1 1 + 1 1 + 1 1 + 1 1 Without ILU 1 1 - 1 1 + 1 1 - 1 1 + 1 1 - 1 1 N/A 1 2 QMR 1 1 With ILU 1 1 + 1 1 + 1 1 + 1 1 + 1 1 + 1 1 + 1 1 Without ILU 1 1 - 1 1 + 1 1 - 1 1 + 1 1 - 1 1 N/A 1 2 CGS 1 1 With ILU 1 1 + 1 1 + 1 1 + 1 1 + 1 1 + 1 1 + 1 1 Without ILU 1 1 - 1 1 + 1 1 - 1 1 + 1 1 - 1 1 N/A 1 2 BICGSTAB 1 1 With ILU 1 1 + 1 1 + 1 1 + 1 1 + 1 1 + 1 1 + 1 1 Without ILU 1 1 - 1 1 + 1 1 - 1 1 + 1 1 - 1 1 N/A 1 2 BICG 1 1 With ILU 1 1 + 1 1 + 1 1 + 1 1 + 1 1 + 1 1 + 1 1 Without ILU 1 1 - 1 1 + 1 1 - 1 1 + 1 1 - 1 1 N/A 4 Summary and future work 1 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 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 http://doi.org/10.1016/j.pnucene.2014.05.011 Theory and implementation of nuclear safety system codes-Part I: Conservation equations, flow regimes, numerics and significant assumptions G. A. Roth , F. Aydogan , Theory and implementation of nuclear safety system codes-Part II: System code closure relations, validation, and limitations . Prog. Nucl. Energy 76 , 55 - 72 ( 2014 ). doi: 10.1016/j.pnucene.2014.05.003 http://doi.org/10.1016/j.pnucene.2014.05.003 Theory and implementation of nuclear safety system codes-Part II: System code closure relations, validation, and limitations Y. P. Zhang , S. P. Niu , L. T. Zhang , et al ., A review on analysis of LWR severe accident . J. Nucl. Eng. Radiat. Sci . 1 ( 4 ), 041018 ( 2015 ). doi: 10.1115/1.4030364 http://doi.org/10.1115/1.4030364 A review on analysis of LWR severe accident X.Z. Ma , H.J. Jia , Y. Liu , Research on the steam-gas pressurizer model with Relap5 code . Nucl. Sci. Tech . 28 , 58 ( 2017 ). doi: 10.1007/s41365-017-0223-x http://doi.org/10.1007/s41365-017-0223-x Research on the steam-gas pressurizer model with Relap5 code F. A. Braz Filho , G. Sabundjian , G. B. Ribeiro , et al ., Assessment of RELAP5 matrix solvers for a two-phase natural circulation loop . Ann. Nucl. Energy 105 , 249 - 258 ( 2017 ). doi: 10.1016/j.anucene.2017.03.023 http://doi.org/10.1016/j.anucene.2017.03.023 Assessment of RELAP5 matrix solvers for a two-phase natural circulation loop R. A. Berry , J. W. Peterson , H. Zhang , et al ., Relap-7 theory manual (No. INL/EXT-14-31366-Rev003) . Idaho National Lab. (INL) , Idaho Falls, ID (United States) , 2018 . Relap-7 theory manual (No. INL/EXT-14-31366-Rev003) Hajizadeh, A. , Kazeminejad, H. & Talebi, S . High-order fully implicit SIMPLE-based model for fully implicit simulation of upward two-phase flow . Nucl. Sci. Tech . 29 , 114 ( 2018 ). doi: 10.1007/s41365-018-0456-3 http://doi.org/10.1007/s41365-018-0456-3 High-order fully implicit SIMPLE-based model for fully implicit simulation of upward two-phase flow C. N. Dawson , H. Klíe , M. F. Wheeler , et al ., A parallel, implicit, cell‐centered method for two‐phase flow with a preconditioned Newton-Krylov solver . Comput. Geosci . 1 , 215 - 249 ( 1997 ). doi: 10.1023/A:1011521413158 http://doi.org/10.1023/A:1011521413158 A parallel, implicit, cell‐centered method for two‐phase flow with a preconditioned Newton-Krylov solver M. Nordsveen , R. Moe , Preconditioned Krylov subspace methods used in solving two‐dimensional transient two‐phase flows . Int. J. Numer. Methods Fluids 31 ( 7 ), 1141 - 1156 ( 1999 ). doi: https://doi.org/10.1002/(SICI)1097-0363(19991215)31:7%3c1141::AIDFLD916%3e3.0.CO;2-G https://doi.org/10.1002/(SICI)1097-0363(19991215)31:7%3c1141::AIDFLD916%3e3.0.CO;2-G Preconditioned Krylov subspace methods used in solving two‐dimensional transient two‐phase flows T. J. Downar , H. G. Joo , A preconditioned Krylov method for solution of the multi-dimensional, two fluid hydrodynamics equations . Ann. Nucl. Energy 28 ( 12 ), 1251 - 1267 ( 2001 ). doi: 10.1016/S0306-4549(00)00124-9 http://doi.org/10.1016/S0306-4549(00)00124-9 A preconditioned Krylov method for solution of the multi-dimensional, two fluid hydrodynamics equations L. Bergamaschi , R. Bru , A. Martı , et al ., Quasi-Newton acceleration of ILU preconditioners for nonlinear two-phase flow equations in porous media . Adv. Eng. Softw . 46 ( 1 ), 63 - 68 ( 2012 ). doi: 10.1016/j.advengsoft.2010.10.011 http://doi.org/10.1016/j.advengsoft.2010.10.011 Quasi-Newton acceleration of ILU preconditioners for nonlinear two-phase flow equations in porous media M. Mohitpour , G. Jahanfarnia , M. Shams , An advancement in iterative solution schemes for three-dimensional, two-fluid modeling of two-phase flow in PWR fuel bundles . Ann. Nucl. Energy 63 , 83 - 99 ( 2014 ). doi: 10.1016/j.anucene.2013.07.007 http://doi.org/10.1016/j.anucene.2013.07.007 An advancement in iterative solution schemes for three-dimensional, two-fluid modeling of two-phase flow in PWR fuel bundles Liu, LX. , Hao, C. & Xu, YL . Equivalent low-order angular flux nonlinear finite difference equation of MOC transport calculation . NUCL SCI TECH 31 , 125 ( 2020 ). doi: 10.1007/s41365-020-00834-2 http://doi.org/10.1007/s41365-020-00834-2 Equivalent low-order angular flux nonlinear finite difference equation of MOC transport calculation V. H. Ransom , Numerical benchmark test No. 2.1: Faucet flow . Multiph. Sci. Technol . 3 ( 1-4 ), 465 - 467 ( 1987 ). DOI: 10.1615/MultScienTechn.v3.i1-4.280 http://doi.org/10.1615/MultScienTechn.v3.i1-4.280 Numerical benchmark test No. 2.1: Faucet flow L. Zou , H. Zhao , H. Zhang , New analytical solutions to the two-phase water faucet problem . Prog. Nucl. Energy 91 , 389 - 398 ( 2016 ). doi: 10.1016/j.pnucene.2016.05.013 http://doi.org/10.1016/j.pnucene.2016.05.013 New analytical solutions to the two-phase water faucet problem R. Nourgaliev , N. Dinh , T. Theofanous , A characteristics-based approach to the numerical solution of the two-fluid model . In Fluids Engineering Division Summer Meeting , Vol. 36967 , 1729 - 1748 ( 2003 ). https://doi.org/10.1115/FEDSM2003-45551 https://doi.org/10.1115/FEDSM2003-45551 A characteristics-based approach to the numerical solution of the two-fluid model T. Gallouët , J. M. Hérard , N. Seguin , Numerical modeling of two-phase flows using the two-fluid two-pressure approach . Math. Models Methods Appl. Sci . 14 (0 5 ), 663 - 700 ( 2004 ). doi: 10.1142/S0218202504003404 http://doi.org/10.1142/S0218202504003404 Numerical modeling of two-phase flows using the two-fluid two-pressure approach A. Morin , T. Flåtten , S. T. Munkejord , A Roe scheme for a compressible six‐equation two‐fluid model . Int. J. Numer. Methods Fluids 72 ( 4 ), 478 - 504 ( 2013 ). doi: 10.1002/fld.3752 http://doi.org/10.1002/fld.3752 A Roe scheme for a compressible six‐equation two‐fluid model W. D. Fullmer , M. A. L. de Bertodano , X. Zhang , Verification of a higher-order finite difference scheme for the one-dimensional two-fluid model . J. Comput. Multiph. Flows 5 ( 2 ), 139 - 155 ( 2013 ). doi: 10.1260/1757-482X.5.2.139 http://doi.org/10.1260/1757-482X.5.2.139 Verification of a higher-order finite difference scheme for the one-dimensional two-fluid model M. O. Delchini , J. C. Ragusa , R. A. Berry , Simulations of single-and two-phase shock tubes and gravity-driven wave problems with the RELAP-7 nuclear reactor system analysis code . Nucl. Eng. Des . 319 , 106 - 116 ( 2017 ). doi: 10.1016/j.nucengdes.2017.04.034 http://doi.org/10.1016/j.nucengdes.2017.04.034 Simulations of single-and two-phase shock tubes and gravity-driven wave problems with the RELAP-7 nuclear reactor system analysis code J. Chen , H. Chen , X. Zhang , Implementation and validation of a one-step coupled solution method for the two-fluid model . Nucl. Eng. Des . 348 , 56 - 64 ( 2019 ). doi: 10.1016/j.nucengdes.2019.04.008 http://doi.org/10.1016/j.nucengdes.2019.04.008 Implementation and validation of a one-step coupled solution method for the two-fluid model Information Systems Laboratories, Inc . RELAP5/MOD3.3 code manual volume i: code structure, system models, and solution methods , ( Idaho Falls, Idaho , 2001 ) RELAP5/MOD3.3 code manual volume i: code structure, system models, and solution methods R. Barrett , M. Berry , T. F. Chan , et al ., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods . ( Society for Industrial and Applied Mathematics , 1994 ) Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods Y. Saad , Iterative methods for sparse linear systems . ( SIAM , 2003 ) Iterative methods for sparse linear systems All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Hong-Yang Wei, Kevin Briggs, Victor Quintanilla and Yi-Tung Chen. The first draft of the manuscript was written by Hong-Yang Wei and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. 10.1007/s41365-021-00888-w 2021-02-15 2021-03-21 2021-03-29 2021-05-07 2021-05-01 2021-12-20 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2021 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-021-00888-w http://dx.doi.org/10.1007/s41365-021-00888-w . 2021 Nuclear Science and Techniques 05 32