CN114169203B - A fully implicit numerical method for two-phase flow for transient safety analysis of nuclear power plants - Google Patents
A fully implicit numerical method for two-phase flow for transient safety analysis of nuclear power plants Download PDFInfo
- Publication number
- CN114169203B CN114169203B CN202111501842.1A CN202111501842A CN114169203B CN 114169203 B CN114169203 B CN 114169203B CN 202111501842 A CN202111501842 A CN 202111501842A CN 114169203 B CN114169203 B CN 114169203B
- Authority
- CN
- China
- Prior art keywords
- equation
- phase
- control
- kth
- vapor phase
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/08—Thermal analysis or thermal optimisation
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Data Mining & Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Operations Research (AREA)
- Computer Hardware Design (AREA)
- Computing Systems (AREA)
- Monitoring And Testing Of Nuclear Reactors (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种用于核电厂瞬态安全分析的两相流全隐数值方法，该方法首先将两流体的两相流模型进行全隐式离散形成方程组，然后基于半隐离散的雅克比矩阵构造预处理矩阵，使用JFNK方法迭代求解全隐式离散的线性方程组直至收敛，输出用于核电厂瞬态安全分析的参数。与现有技术相比，本发明方法使用基于JFNK的全隐式数值算法求解两流体模型，降低了传统的两流体模型全隐式数值求解的难度，提高了数值算法的稳定性和计算效率，并使用了基于半隐式差分的雅克比矩阵的预处理技术，提高了计算效率，能够更加精确且高效地模拟和分析核反应堆正常运行和事故瞬变现象。
The invention discloses a fully implicit numerical method for two-phase flow used for transient safety analysis of nuclear power plants. This method first performs fully implicit discretization of the two-phase flow model of two fluids to form a system of equations, and then based on the semi-implicit discretization of Jacques Construct a preprocessing matrix using the ratio matrix, use the JFNK method to iteratively solve the fully implicit discrete linear equations until convergence, and output parameters for transient safety analysis of nuclear power plants. Compared with the existing technology, the method of the present invention uses a fully implicit numerical algorithm based on JFNK to solve the two-fluid model, which reduces the difficulty of the traditional fully implicit numerical solution of the two-fluid model and improves the stability and calculation efficiency of the numerical algorithm. It also uses preprocessing technology based on semi-implicit differential Jacobian matrix to improve calculation efficiency and can more accurately and efficiently simulate and analyze normal operation and accident transient phenomena of nuclear reactors.
Description
Technical Field
The invention belongs to the technical field of transient safety analysis of nuclear power plants, and particularly relates to a two-phase flow total-hidden numerical method for transient safety analysis of a nuclear power plant.
Background
The nuclear reactor transient safety analysis program is used as a main means and tool for recognizing and researching the normal operation and accident transient phenomenon of the nuclear reactor and verifying the system design, and the importance of the nuclear reactor transient safety analysis program is self-evident. From the last 60 th century until now, a number of well-known reactor transient safety analysis programs have been developed, such as: RELAP5, TRAC, ATHLET, CATHARE, etc. have been widely used in system design and accident safety analysis of various nuclear power plants or nuclear reactors. There are still some problems to be solved and optimized.
Currently widely used procedures for analyzing reactor transient safety include: RELAP5 or TRACE program, it uses two fluid models to simulate two-phase flow, the numerical solution adopts semi-implicit solution, the precision is relatively low, and the numerical stability is affected by the coulomb limit of sound velocity of the material, so the time step cannot be taken very much, the efficiency is low when the program simulates some long transient problems;
for some advanced programs that employ the fully hidden numerical algorithm, such as: the CARHARE2 program adopts Newton method, and because the two fluid models contain a large number of constitutive models, the difficulty is great in constructing the Jacobian matrix of the fully implicit discrete equation;
in order to solve the problems of low semi-implicit numerical efficiency and difficult construction of the traditional fully implicit solution Jacobian matrix, researchers in recent years propose a Newton-Krylov subspace iteration method (JFNK method) without Jacobian matrix. At present, in the field of nuclear reactor transient safety analysis programs, the JFNK algorithm is mostly applied to solving of drift flow models, and when the JFNK algorithm is applied to solving of two-fluid models, the preprocessing technology adopts some pure mathematical preprocessing technologies, and the Jacobian matrix still needs to be formed, so that the difficulty of constructing the Jacobian matrix still exists.
In summary, the semi-implicit solution adopted in the currently-used nuclear reactor transient safety analysis program has certain disadvantages, the traditional total hidden algorithm has certain difficulties, and although some researches apply the JFNK algorithm to the solution of the drift flow model and the two-fluid model, the preprocessing technology is not concise enough, so that the model and the algorithm need to be improved to enable the nuclear reactor transient safety analysis program to be more accurate and efficient.
Disclosure of Invention
The invention provides a two-phase flow total hidden numerical method for transient safety analysis of a nuclear power plant, which improves the instability of a numerical algorithm of a transient safety analysis program of an existing nuclear reactor. Firstly, realizing a fully implicit solution of two fluid models based on an JFNK algorithm, so that a numerical algorithm is not limited by the coulomb value of the sound velocity of a material; secondly, a pre-processing technology of the Jacobian matrix based on the semi-hidden difference is established, and the calculation efficiency is improved; finally, a two-phase flow total hidden numerical method for transient safety analysis of the nuclear power plant is obtained.
The invention is realized by adopting the following technical scheme:
a two-phase flow total hidden numerical method for transient safety analysis of a nuclear power plant, comprising:
the first step: dividing control bodies according to the operation parameters and pipeline structure parameters of the nuclear power plant system, and establishing a nonlinear equation set of a two-fluid two-phase flow model of all the control bodies:
wherein,residual form of the vapor phase momentum equation at the inlet boundary i-1/2 for control volume i;Residual form of the liquid phase momentum equation at the inlet boundary i-1/2 of the control volume i;Residual form of the vapor phase energy equation for control volume i;Controlling the residual form of the liquid phase energy equation of the body i;Controlling the residual form of the vapor phase mass conservation equation of the body i; f (F) P,i Controlling the residual form of the mixed phase mass conservation equation of the body i;Residual form of the vapor phase momentum equation at the outlet boundary i+1/2 of the control volume i;Residual form of the liquid phase momentum equation at the outlet boundary i+1/2 of the control volume i;
and a second step of: and solving a nonlinear equation set established in the (n+1) th time step and the (k) th Newton iteration step by adopting a JFNK method to obtain the vapor phase internal energy, the liquid phase internal energy, the cavitation share, the pressure, the vapor phase speed and the liquid phase speed of all control bodies in the (n+1) th time step, wherein the JFNK method comprises the following steps: the preprocessing matrix M, krylov subspace iteration and newton iteration are constructed.
The invention is further improved in that the conservation equation of the two-fluid two-phase flow model in the first step is discretized based on staggered grids and time backward difference and space first-order windward difference, and the obtained fully implicit discrete equation of the control body i and the outlet boundary i+1/2 is as follows:
k-phase energy equation for control i, k=g for vapor phase and k=l for liquid phase:
vapor phase mass equation for control volume i
Mixed phase mass equation for control volume i
K-phase momentum equation controlling the bulk i exit boundary i+1/2:
wherein: Δx is the control body length; Δt is the time step of transient analysis; a is the cross section of the control body; p is the control body pressure; alpha k 、ρ k 、U k Respectively representing the volume fraction, density and internal energy of k phase; Γ -shaped structure i,k 、Q w,k 、Q i,k 、E w,k 、E i,k Respectively representing the transmission quantity of k-phase unit volume, the heat exchange quantity of a control body and a wall surface, the heat exchange quantity between a control body gas phase and a liquid phase, the heat transfer quantity caused by wall surface mass transfer and the heat transfer quantity caused by interphase mass transfer; v (V) k 、F g,k 、F w,k 、F c,k 、F i,k 、M k Respectively representing the velocity of k phase, wall gravity, wall friction, local resistance, interphase friction and momentum transfer caused by mass transfer;the quantity of the boundary of the control body is assigned through first-order windward differential assignment; k=g represents a vapor phase, and k=l represents a liquid phase.
The invention is further improved in that in the second step of the method, firstly, the nonlinear equation set established in the first step is linearized:
J(x k )·δx k ＝-F(x k )
wherein J (x) k ) Jacobian matrix, δx, for the kth Newton iteration k Solution variable increment for the kth Newton iteration, F (x k ) Nonlinear equation set, x, for the kth Newton iteration k Independent variable for the kth Newton iteration
Wherein,controlling the vapor phase velocity at the inlet boundary i-1/2 of volume i for the kth Newton iteration;Controlling the liquid phase velocity at the inlet boundary i-1/2 of the body i for the kth Newton iteration;Controlling the vapor phase internal energy of the body i for the kth Newton iteration;Controlling the liquid phase internal energy of the body i for the kth Newton iteration;Controlling the fraction of vapor phase cavitation of volume i for the kth newton iteration; p (P) i k Controlling the pressure of the body i for the kth newton iteration;Controlling the vapor phase velocity at the exit boundary i+1/2 of the volume i for the kth Newton iteration;The liquid phase velocity at the exit boundary i+1/2 of volume i is controlled for the kth Newton iteration.
The invention is further improved in that the linear equation set is preprocessed, and becomes:
(J(x k )M -1 )·(Mδx k )＝-F(x k )
the construction method of the pretreatment matrix M comprises the following steps:
based on the staggered grid and the finite volume difference, semi-implicit discrete is carried out on the two-fluid two-phase flow equation to obtain a nonlinear equation set
Wherein,a residual form after semi-implicit discrete of a vapor phase momentum equation at an inlet boundary i-1/2 of a control body i;The residual form after the semi-implicit discrete of the liquid phase momentum equation at the inlet boundary i-1/2 of the control body i;a residual form after semi-implicit discrete of a vapor phase energy equation for the control volume i;The residual form after the semi-implicit discrete of the liquid phase energy equation of the control body i;A residual form after semi-implicit discrete of a vapor phase mass conservation equation of the control body i; f (F) P ′ ,i A residual form after semi-implicit discrete of a mixed phase mass conservation equation of the control body i;A residual form after semi-implicit discrete of a vapor phase momentum equation at an outlet boundary i+1/2 of the control body i;The residual form after the semi-implicit discrete of the liquid phase momentum equation at the outlet boundary i+1/2 of the control body i;
the elements of the d-th row and e-th column of the pretreatment matrix M are:
wherein F' (x) k ) d Is the d-th equation in the linear equation set F' (x),to solve for the variable x k The e-th solution variable in (a).
The invention further improves that the linear equation system after pretreatment is solved by using a Krylov subspace iteration method, which concretely comprises the following steps:
first, a standard orthogonal basis vector Vm and an epinastine matrix of a Krylov subspace Km with dimension m are established through an Arnoldi processNamely:
Vm＝(v 1 ,v 2 ,...,v m )
next, the preprocessed equation set (J (x k )M -1 )·(Mδx k )＝-F(x k ) Conversion to a new system of equations
Wherein r is (0) ＝-F(x k )-A·x k ，A＝J(x k )M -1 M is the dimension of the Krylov subspace, z m ∈R n To solve for the variables, z (m) ＝V m ·Y＝(v 1 ,v 2 ,...,v m )·(y 1 ,y 2 ,...,y m ) T ，y i (i=1, m) is the basis vector v i Corresponding substrate, β= ||r (0) ||，
New system of equationsIs an overdetermined equation set with m unknowns and m+1 equations, and can be solved by using QR decomposition based on Hastelloy transformation to obtain Y, so as to calculate z m 。
The invention is further improved in that the dot product calculation of the Jacobian matrix and the vector involved in the Krylov subspace iteration method is all calculated by differential approximation as follows:
wherein epsilon is a differential step length or a disturbance parameter, and x is solved according to the kth step k Vector v k Machine accuracy epsilon mach Determination of
The invention is further improved in that the convergence condition of the Krylov subspace iteration is two:
1)
2)
if either convergence condition is satisfied, the Krylov subspace iteratively converges, a system of equations (J (x k )M -1 )·(Mδx k )＝-F(x k ) Solution of δx k ＝M -1 (x (k) +z m ) The method comprises the steps of carrying out a first treatment on the surface of the Otherwise the dimension m of the Krylov subspace becomes m+1, and the system of equations is solved again (J (x k )M -1 )·(Mδx k )＝-F(x k )。
The invention is further improved in that the solution delta x of the linear equation system obtained by iterative convergence of the Krylov subspace k Updating solution x for the kth Newton iteration k
x k+1 ＝x k +δx k
Judging whether Newton iteration converges or not, wherein the Newton iteration converges conditions are three:
1)||F(x k+1 )||≤1.0 -3 ；
2)
3)
wherein nvol is the total number of control bodies in the nuclear reactor system, njun is the total number of takeover in the nuclear reactor system;
if the kth newton iteration does not meet all of the above conditions, then the newton iteration is considered non-converging, k=k+1, with x k+1 As x k Re-performing Krylov subspace iteration; otherwise Newton iteration converges, x k+1 For the solution of the n+1th time step, according to x k +1 The vapor phase internal energy, liquid phase internal energy, void fraction, pressure, vapor phase velocity, and liquid phase velocity of all control bodies can be obtained for the n+1th time step.
Compared with the prior art, the invention has at least the following beneficial technical effects:
1: because the conventional nuclear reactor transient safety analysis program adopts a semi-implicit solution, the stability of the numerical method is limited by the coulomb value of the sound velocity of the material, and the time step cannot be taken very much, so that the efficiency is lower when the problem of slow transient is calculated. The method uses the fully hidden numerical method of the two-fluid model based on the JFNK algorithm, can accurately and efficiently simulate the two-phase flow problem in a nuclear reactor system, and can save the difficulty of writing the Jacobian matrix compared with the traditional Newton method.
2: because the JFNK algorithm is currently applied to the research of two-fluid models, the preprocessing technology mostly adopts the preprocessing technology based on the pure numerology built in the PETSc program package, and the Jacobian matrix based on the fully hidden numerical method needs to be written. According to the invention, the pretreatment technology of the Jacobian matrix based on the semi-hidden difference is used, only the Jacobian matrix based on the semi-hidden numerical method is needed to be written, the difficulty of the numerical method is greatly reduced, and the pretreatment matrix has physical significance.
Drawings
Fig. 1 is a flow chart of the fully hidden numerical method for the two-phase flow for transient safety analysis of the nuclear power plant.
Fig. 2 is a graph comparing cavitation share distribution calculated by the total hidden numerical method with experimental values and calculation results of the semi-hidden numerical method, wherein fig. 2 (a) is a working condition 1, fig. 2 (b) is a working condition 2, fig. 2 (c) is a working condition 3, and fig. 2 (d) is a working condition 4.
Detailed Description
Exemplary embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be embodied in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art. It should be noted that, without conflict, the embodiments of the present invention and features of the embodiments may be combined with each other.
The implementation steps of the two-phase flow total hidden numerical method for transient safety analysis of a nuclear power plant proposed in the present invention will be described in detail with reference to fig. 1, taking a one-dimensional vertical pipe flow heat transfer problem similar to the flow heat transfer of the reactor core of a nuclear power plant as an example.
The first step: dividing a one-dimensional reactor core into 20 control bodies from bottom to top, and establishing a two-fluid two-phase flow equation of each control body
Wherein,to control the vapor phase at the inlet boundary i-1/2 of the body iResidual form of momentum equation;Residual form of the liquid phase momentum equation at the inlet boundary i-1/2 of the control volume i;Residual form of the vapor phase energy equation for control volume i;Controlling the residual form of the liquid phase energy equation of the body i;Controlling the residual form of the vapor phase mass conservation equation of the body i; f (F) P,i Controlling the residual form of the mixed (vapor phase + liquid phase) mass conservation equation of the body i;Residual form of the vapor phase momentum equation at the outlet boundary i+1/2 of the control volume i;Residual form of the liquid phase momentum equation at the outlet boundary i+1/2 of the control volume i;
and a second step of: dispersing the conservation equation of the two-fluid two-phase flow model in the first step based on staggered grids and time backward difference and space first-order windward difference to obtain a fully implicit discrete equation of the control body i and the outlet boundary i+1/2 as follows:
k (k=g represents vapor phase, k=l represents liquid phase) phase energy equation of control body i:
vapor phase mass equation for control volume i
Mixed phase (vapor phase + liquid phase) mass equation for control volume i
K (k=g represents vapor phase, k=l represents liquid phase) phase momentum equation for control volume i exit boundary i+1/2:
wherein: Δx is the control body length; Δt is the time step of transient analysis; a is the cross section of the control body; p is the control body pressure; alpha k 、ρ k 、U k Respectively representing the volume fraction, density and internal energy of k phase; Γ -shaped structure i,k 、Q w,k 、Q i,k 、E w,k 、E i,k Respectively representing the transmission quantity of k-phase unit volume, the heat exchange quantity of a control body and a wall surface, the heat exchange quantity between a control body gas phase and a liquid phase, the heat transfer quantity caused by wall surface mass transfer and the heat transfer quantity caused by interphase mass transfer; v (V) k 、F g,k 、F w,k 、F c,k 、F i,k 、M k Respectively representing the velocity of k phase, wall gravity, wall friction, local resistance, interphase friction and momentum transfer caused by mass transfer;the quantity of the boundary of the control body is assigned through first-order windward differential assignment; k=g represents a vapor phase, and k=l represents a liquid phase.
Forming a nonlinear equation set by the total 120 equations of the 20 control bodies and the inlet-outlet boundaries of each control body:
the solution variables are:
x＝(…,V g,i-1/2 ,V l,i-1/2 ,U g,i ,U l,i ,α g,i ,P i ,V g,i+1/2 ,V l,i+1/2 ,…) T
wherein V is g,i-1/2 To control the vapor phase velocity at the inlet boundary i-1/2 of the body i; v (V) l,i-1/2 To control the liquid phase velocity at the inlet boundary i-1/2 of the body i; u (U) g,i To control the vapor phase internal energy of the body i; u (U) l,i To control the liquid phase internal energy of the body i; alpha g,i To control the fraction of vapor phase cavitation of the body i; p (P) i For controlling the pressure of the body i; v (V) g,i+1/2 To control the vapor phase velocity at the outlet boundary i+1/2 of the body i; v (V) l,i+1/2 To control the liquid phase velocity at the outlet boundary i+1/2 of the body i;
and a third step of: linearizing the nonlinear equation set established in the second step:
J(x k )·δx k ＝-F(x k )
wherein J (x) k ) Jacobian matrix, δx, for the kth Newton iteration k Solution variable increment for the kth Newton iteration, F (x k ) Nonlinear equation set, x, for the kth Newton iteration k Independent variable for the kth Newton iteration
Fourth step: the linear equation set formed in the third step is preprocessed, and the construction method of the preprocessing matrix M comprises the following steps:
based on the staggered grid and the finite volume difference, semi-implicit discrete is carried out on the two-fluid two-phase flow equation to obtain a nonlinear equation set
Wherein,a residual form after semi-implicit discrete of a vapor phase momentum equation at an inlet boundary i-1/2 of a control body i;The residual form after the semi-implicit discrete of the liquid phase momentum equation at the inlet boundary i-1/2 of the control body i;a residual form after semi-implicit discrete of a vapor phase energy equation for the control volume i;The residual form after the semi-implicit discrete of the liquid phase energy equation of the control body i;A residual form after semi-implicit discrete of a vapor phase mass conservation equation of the control body i; f (F) P ′ ,i A residual form after semi-implicit discretization of a mixed (vapor phase+liquid phase) mass conservation equation for the control body i;A residual form after semi-implicit discrete of a vapor phase momentum equation at an outlet boundary i+1/2 of the control body i;The residual form after the semi-implicit discrete of the liquid phase momentum equation at the outlet boundary i+1/2 of the control body i;
the elements of the d-th row and e-th column of the pretreatment matrix M are:
wherein F' (x) k ) d Is the d-th equation in the linear equation set F' (x),to solve for the variable x k The e-th solution variable in (a).
The linear equation set becomes after pretreatment:
(J(x k )M -1 )·(Mδx k )＝-F(x k )
and fifthly, solving the preprocessed linear equation set formed in the fourth step by using a Krylov subspace iteration method, wherein the method specifically comprises the following steps:
first, a standard orthogonal basis vector Vm and a Shandong Berger matrix of a Krylov subspace Km of dimension m (first iteration m is 1) are established by an Arnoldi processNamely:
Vm＝(v 1 ,v 2 ,...,v m )
next, the preprocessed equation set (J (x k )M -1 )·(Mδx k )＝-F(x k ) Conversion to a new system of equations
Wherein r is (0) ＝-F(x k )-A·x k ，A＝J(x k )M -1 M is the dimension of the Krylov subspace, z m ∈R n To solve for the variables, z (m) ＝V m ·Y＝(v 1 ,v 2 ,...,v m )·(y 1 ,y 2 ,...,y m ) T ，y i (i=1, m) is the basis vector v i Corresponding substrate, β= ||r (0) ||，
New system of equationsIs an overdetermined equation set with m unknowns and m+1 equations, and can be solved by using QR decomposition based on Hastelloy transformation to obtain Y, so as to calculate z m 。
The dot product calculation of the jacobian matrix and the vector involved in the fifth step Krylov subspace iteration method is all calculated by differential approximation as follows:
wherein epsilon is a differential step length or a disturbance parameter, and x is solved according to the kth step k Vector v k Machine accuracy epsilon mach Determination of
For a commonly used 64-bit computer ε mach ＝10 -8 。
Sixth step: z for the band solved in the fifth step m Judging whether the Krylov subspace iteration converges or not, wherein two convergence conditions are provided:
1)
2)
if either convergence condition is satisfied, the Krylov subspace iteratively converges, a system of equations (J (x k )M -1 )·(Mδx k )＝-F(x k ) Solution of δx k ＝M -1 (x (k) +z m ) The method comprises the steps of carrying out a first treatment on the surface of the Otherwise Krylov subspaceThe dimension m of (2) becomes m+1, and returns to the fifth step again, and the solution equation set (J (x k )M -1 )·(Mδx k )＝-F(x k ) Solution z of (2) m+1 ；
Seventh, solution delta x of linear equation set obtained by utilizing Krylov subspace iteration convergence in sixth step k Updating solution x for the kth Newton iteration k+1 ＝x k +δx k And judging whether Newton iteration converges or not. The newton iteration convergence condition has three:
1)||F(x k+1 )||≤1.0 -3 ；
2)
3)
if the kth newton iteration does not meet all of the above conditions, then the newton iteration is considered non-converging, k=k+1, with x k+1 As x k The fifth, sixth and seventh steps are carried out again; otherwise Newton iteration converges, x k+1 For the solution of the n+1th time step, according to x k+1 The vapor phase internal energy, liquid phase internal energy, void fraction, pressure, vapor phase velocity, and liquid phase velocity of all control bodies can be obtained for the n+1th time step.
Through the above-described procedure, the numerical calculation results are shown in fig. 2. The calculation results obtained by the fully hidden numerical algorithm in the invention are better in accordance with the experimental results and are similar to the calculation results of the semi-hidden numerical algorithm in precision. In numerical simulation of the several groups of experimental conditions, limited by the count value, in order to obtain stable time numerical calculation results, the time step of the explicit or semi-implicit numerical algorithm of the several groups of conditions is not more than 0.008s, typically 0.005s or less, while in the fully implicit numerical method, the time step Δt=0.01 s. Therefore, the time step of the fully implicit numerical method is not influenced by the count value, and the numerical calculation efficiency of the transient safety analysis of the nuclear power plant can be improved.
While the invention has been described in detail in the foregoing general description and with reference to specific embodiments thereof, it will be apparent to one skilled in the art that modifications and improvements can be made thereto. Accordingly, such modifications or improvements may be made without departing from the spirit of the invention and are intended to be within the scope of the invention as claimed.
Claims (7)
1. A two-phase flow total hidden numerical method for transient safety analysis of a nuclear power plant, comprising:
the first step: dividing control bodies according to the operation parameters and pipeline structure parameters of the nuclear power plant system, and establishing a nonlinear equation set of a two-fluid two-phase flow model of all the control bodies:
wherein,residual form of the vapor phase momentum equation at the inlet boundary i-1/2 for control volume i;Residual form of the liquid phase momentum equation at the inlet boundary i-1/2 of the control volume i;Residual form of the vapor phase energy equation for control volume i;Controlling the residual form of the liquid phase energy equation of the body i;Controlling the residual form of the vapor phase mass conservation equation of the body i; f (F) P,i Control bodyi residual form of mixed phase conservation of mass equation;Residual form of the vapor phase momentum equation at the outlet boundary i+1/2 of the control volume i;Residual form of the liquid phase momentum equation at the outlet boundary i+1/2 of the control volume i;
and a second step of: and solving a nonlinear equation set established in the (n+1) th time step and the (k) th Newton iteration step by adopting a JFNK method to obtain the vapor phase internal energy, the liquid phase internal energy, the cavitation share, the pressure, the vapor phase speed and the liquid phase speed of all control bodies in the (n+1) th time step, wherein the JFNK method comprises the following steps: constructing a preprocessing matrix M, krylov subspace iteration and a Newton iteration; firstly, linearizing the nonlinear equation set established in the first step:
J(x k )·δx k ＝-F(x k )
wherein J (x) k ) Jacobian matrix, δx, for the kth Newton iteration k Solution variable increment for the kth Newton iteration, F (x k ) Nonlinear equation set, x, for the kth Newton iteration k Independent variable for the kth Newton iteration
Wherein,controlling the vapor phase velocity at the inlet boundary i-1/2 of volume i for the kth Newton iteration;Controlling the liquid phase velocity at the inlet boundary i-1/2 of the body i for the kth Newton iteration;Controlling the vapor phase internal energy of the body i for the kth Newton iteration;Controlling the liquid phase internal energy of the body i for the kth Newton iteration;Controlling the fraction of vapor phase cavitation of volume i for the kth newton iteration; p (P) i k Controlling the pressure of the body i for the kth newton iteration;Controlling the vapor phase velocity at the exit boundary i+1/2 of the volume i for the kth Newton iteration;The liquid phase velocity at the exit boundary i+1/2 of volume i is controlled for the kth Newton iteration.
2. The two-phase flow total hidden numerical method for transient safety analysis of a nuclear power plant according to claim 1, wherein the conservation equation of the two-fluid two-phase flow model in the first step is discretized based on a staggered grid and a time backward differential, space first-order windward differential, and the obtained total hidden discrete equation of the control body i and the outlet boundary i+1/2 is as follows:
k-phase energy equation for control i, k=g for vapor phase and k=l for liquid phase:
vapor phase mass equation for control volume i
Mixed phase mass equation for control volume i
K-phase momentum equation controlling the bulk i exit boundary i+1/2:
wherein: Δx is the control body length; Δt is the time step of transient analysis; a is the cross section of the control body; p is the control body pressure; alpha k 、ρ k 、U k Respectively representing the volume fraction, density and internal energy of k phase; Γ -shaped structure i,k 、Q w,k 、Q i,k 、E w,k 、E i,k Respectively representing the transmission quantity of k-phase unit volume, the heat exchange quantity of a control body and a wall surface, the heat exchange quantity between a control body gas phase and a liquid phase, the heat transfer quantity caused by wall surface mass transfer and the heat transfer quantity caused by interphase mass transfer; v (V) k 、F g,k 、F w,k 、F c,k 、F i,k 、M k Respectively representing the velocity of k phase, wall gravity, wall friction, local resistance, interphase friction and momentum transfer caused by mass transfer;the quantity of the boundary of the control body is assigned through first-order windward differential assignment; k=g represents a vapor phase, and k=l represents a liquid phase.
3. The two-phase flow total hidden numerical method for transient safety analysis of a nuclear power plant according to claim 1, wherein the linear equation set is preprocessed, and becomes:
(J(x k )M -1 )·(Mδx k )＝-F(x k )
the construction method of the pretreatment matrix M comprises the following steps:
based on the staggered grid and the finite volume difference, semi-implicit discrete is carried out on the two-fluid two-phase flow equation to obtain a nonlinear equation set
Wherein,a residual form after semi-implicit discrete of a vapor phase momentum equation at an inlet boundary i-1/2 of a control body i;the residual form after the semi-implicit discrete of the liquid phase momentum equation at the inlet boundary i-1/2 of the control body i;A residual form after semi-implicit discrete of a vapor phase energy equation for the control volume i;The residual form after the semi-implicit discrete of the liquid phase energy equation of the control body i;A residual form after semi-implicit discrete of a vapor phase mass conservation equation of the control body i; f'. P,i A residual form after semi-implicit discrete of a mixed phase mass conservation equation of the control body i;A residual form after semi-implicit discrete of a vapor phase momentum equation at an outlet boundary i+1/2 of the control body i;To controlA residual form after semi-implicit discretization of a liquid phase momentum equation at an outlet boundary i+1/2 of the body i;
the elements of the d-th row and e-th column of the pretreatment matrix M are:
wherein F' (x) k ) d Is the d-th equation in the linear equation set F' (x),to solve for the variable x k The e-th solution variable in (a).
4. The two-phase flow total hidden numerical method for transient safety analysis of a nuclear power plant according to claim 1, wherein the preprocessed linear equation set is solved by using a Krylov subspace iteration method, and the method specifically comprises the following steps:
first, a standard orthogonal basis vector Vm and an epinastine matrix of a Krylov subspace Km with dimension m are established through an Arnoldi processNamely:
Vm=(v 1 ,v 2 ,...,v m )
next, willThe preprocessed system of equations (J (x k )M -1 )·(Mδx k )＝-F(x k ) Conversion to a new system of equations
Wherein r is (0) ＝-F(x k )-A·x k ，A＝J(x k )M -1 M is the dimension of the Krylov subspace, z m ∈R n To solve for the variables, z (m) ＝V m ·Y＝(v 1 ,v 2 ,...,v m )·(y 1 ,y 2 ,...,y m ) T ，y i (i=1, m) is the basis vector v i Corresponding substrate, β= ||r (0 )||，
New system of equationsIs an overdetermined equation set with m unknowns and m+1 equations, and can be solved by using QR decomposition based on Hastelloy transformation to obtain Y, so as to calculate z m 。
5. The two-phase flow total hidden numerical method for transient safety analysis of a nuclear power plant according to claim 4, wherein the dot product calculation of the jacobian matrix and the vector involved in the Krylov subspace iteration method is all calculated by differential approximation:
wherein epsilon is a differential step length or a disturbance parameter, and x is solved according to the kth step k Vector v k Machine accuracy epsilon mach Determination of
6. The method for two-phase flow total hidden numerical value for transient safety analysis of nuclear power plant according to claim 4, wherein the convergence condition of Krylov subspace iteration has two:
1)
2)
if either convergence condition is satisfied, the Krylov subspace iteratively converges, a system of equations (J (x k )M -1 )·(Mδx k )＝-F(x k ) Solution of δx k ＝M -1 (x (k) +z m ) The method comprises the steps of carrying out a first treatment on the surface of the Otherwise the dimension m of the Krylov subspace becomes m+1, and the system of equations is solved again (J (x k )M -1 )·(Mδx k )＝-F(x k )。
7. The method for two-phase flow total hidden numerical value for transient safety analysis of nuclear power plant according to claim 6, wherein the solution δx of a linear equation set obtained by iterative convergence of Krylov subspace k Updating solution x for the kth Newton iteration k
x k+1 ＝x k +δx k
Judging whether Newton iteration converges or not, wherein the Newton iteration converges conditions are three:
1)||F(x k+1 )||≤1.0 -3 ；
2)
3)
wherein nvol is the total number of control bodies in the nuclear reactor system, njun is the total number of takeover in the nuclear reactor system;
if the kth newton iteration does not meet all of the above conditions, then the newton iteration is considered non-converging, k=k+1, with x k+1 As x k Re-performing Krylov subspace iteration; otherwise Newton iteration converges, x k+1 For the solution of the n+1th time step, according to x k+1 The vapor phase internal energy, liquid phase internal energy, void fraction, pressure, vapor phase velocity, and liquid phase velocity of all control bodies can be obtained for the n+1th time step.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111501842.1A CN114169203B (en) | 2021-12-09 | 2021-12-09 | A fully implicit numerical method for two-phase flow for transient safety analysis of nuclear power plants |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111501842.1A CN114169203B (en) | 2021-12-09 | 2021-12-09 | A fully implicit numerical method for two-phase flow for transient safety analysis of nuclear power plants |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114169203A CN114169203A (en) | 2022-03-11 |
| CN114169203B true CN114169203B (en) | 2024-01-16 |
Family
ID=80485086
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111501842.1A Expired - Fee Related CN114169203B (en) | 2021-12-09 | 2021-12-09 | A fully implicit numerical method for two-phase flow for transient safety analysis of nuclear power plants |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114169203B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117290642B (en) * | 2022-10-28 | 2024-09-17 | 国家电投集团科学技术研究院有限公司 | Coupling method, device and equipment of thermodynamic and hydraulic model based on Newton-Raphson solver |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104091085A (en) * | 2014-07-18 | 2014-10-08 | 安徽工业大学 | Cavitation noise feature estimation method based on propeller wake flow pressure fluctuation computing |
| CN111125972A (en) * | 2019-12-26 | 2020-05-08 | 西安交通大学 | Hydraulic load analysis method for water loss accident of break of nuclear power plant |
| CN112906271A (en) * | 2021-02-22 | 2021-06-04 | 中国核动力研究设计院 | Reactor transient physical thermal full-coupling fine numerical simulation method and system |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| TW202036599A (en) * | 2018-12-18 | 2020-10-01 | 美商深絕公司 | Radioactive waste repository systems and methods |
-
2021
- 2021-12-09 CN CN202111501842.1A patent/CN114169203B/en not_active Expired - Fee Related
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104091085A (en) * | 2014-07-18 | 2014-10-08 | 安徽工业大学 | Cavitation noise feature estimation method based on propeller wake flow pressure fluctuation computing |
| CN111125972A (en) * | 2019-12-26 | 2020-05-08 | 西安交通大学 | Hydraulic load analysis method for water loss accident of break of nuclear power plant |
| CN112906271A (en) * | 2021-02-22 | 2021-06-04 | 中国核动力研究设计院 | Reactor transient physical thermal full-coupling fine numerical simulation method and system |
Non-Patent Citations (1)
| Title |
|---|
| 核电厂非能动安全壳热量导出系统瞬态模拟程序开发研究;孙涛;徐钊;;核动力工程(03);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114169203A (en) | 2022-03-11 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107066745B (en) | A method for obtaining the three-dimensional neutron flux density distribution of fast neutron reactor core transient process | |
| Owens et al. | Computational fluid dynamics simulation of structured packing | |
| CN114444413B (en) | A sub-channel level three-dimensional thermal hydraulic analysis method for plate fuel reactor core | |
| Ajmani et al. | Preconditioned conjugate gradient methods for the Navier-Stokes equations | |
| CN116956770B (en) | A multi-physics field coupling method for heat pipe reactor core | |
| CN107122331A (en) | A kind of coupling of multiple physics method in presurized water reactor transient state calculating | |
| US20260018310A1 (en) | Three-dimensional thermal-hydraulic analysis method and system for reactor core | |
| CN114169203B (en) | A fully implicit numerical method for two-phase flow for transient safety analysis of nuclear power plants | |
| CN117672430B (en) | Second moment calculation method for turbulent heat flux of liquid metal | |
| Ishii et al. | Interfacial area transport equation and implementation into two-fluid model | |
| CN110543705B (en) | Boiling simulation solving acceleration method in typical channel of nuclear reactor | |
| CN117709212A (en) | Transient safety numerical analysis method and system for gas cooled reactor system | |
| Liu et al. | Numerical simulation of shock wave problems with the two-phase two-fluid model | |
| CN118520638A (en) | Fully hidden transient analysis method and system for thermal hydraulic system of gas cooled reactor | |
| Tentner et al. | Modeling of two-phase flow in a BWR fuel assembly with a highly-scalable CFD code | |
| CN111680405A (en) | A Calculation Method of Hydraulic Characteristics of Natural Circulation Capability | |
| Do et al. | Highly scalable meshless multigrid solver for 3D thermal-hydraulic analysis of nuclear reactors | |
| CN113268828B (en) | Optimal Design Method for Asymmetric Leading Edge Structure of Turbine Blade Based on Uniform Design | |
| Sun et al. | Development and validation of a turbulence model based on drift velocity | |
| Sung et al. | Efficient aerodynamic design method using a tightly coupled algorithm | |
| Guo et al. | A parallel strategy applied to the simplified thermal-hydraulic part of system program | |
| Choi et al. | Development of multi-scale thermal hydraulic analysis code, CUPID/SPACE based on implicit code coupling scheme | |
| CN119720671B (en) | A method for analyzing nuclear reactor core channels based on a two-fluid model | |
| Douglas et al. | Multigrid solution of flame sheet problems on serial and parallel computers | |
| CN119207634A (en) | FVM-HBM Coupling Algorithm for Solving Flow and Mass Transfer Problems |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20240116 |