CN114169203A - Two-phase flow total-hidden numerical method for transient safety analysis of nuclear power plant - Google Patents

Two-phase flow total-hidden numerical method for transient safety analysis of nuclear power plant Download PDF

Info

Publication number
CN114169203A
CN114169203A CN202111501842.1A CN202111501842A CN114169203A CN 114169203 A CN114169203 A CN 114169203A CN 202111501842 A CN202111501842 A CN 202111501842A CN 114169203 A CN114169203 A CN 114169203A
Authority
CN
China
Prior art keywords
equation
phase
implicit
kth
iteration
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.)
Granted
Application number
CN202111501842.1A
Other languages
Chinese (zh)
Other versions
CN114169203B (en
Inventor
苟军利
樊杰
单建强
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN202111501842.1A priority Critical patent/CN114169203B/en
Publication of CN114169203A publication Critical patent/CN114169203A/en
Application granted granted Critical
Publication of CN114169203B publication Critical patent/CN114169203B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR 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; CALCULATING OR 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; CALCULATING OR 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; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR 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; CALCULATING OR 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

Abstract

The invention discloses a two-phase flow total-hidden numerical method for transient safety analysis of a nuclear power plant. Compared with the prior art, the method solves the two fluid models by using the full-implicit numerical algorithm based on JFNK, reduces the difficulty of full-implicit numerical solution of the traditional two fluid models, improves the stability and the calculation efficiency of the numerical algorithm, improves the calculation efficiency by using the preprocessing technology of the Jacobian matrix based on semi-implicit difference, and can more accurately and efficiently simulate and analyze the normal operation and accident transient phenomena of the nuclear reactor.

Description

Two-phase flow total-hidden numerical method for transient safety analysis of nuclear power plant
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 of self-evident importance as a primary means and tool for understanding and studying nuclear reactor normal operation and accident transients, and system design verification. From the last 60 s until now, there have been many well-known reactor transient safety analysis programs developed, such as: RELAP5, TRAC, ATHLET, CATARE, 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.
The current widely used reactor transient safety analysis procedures are as follows: the RELAP5 or TRACE program uses two fluid models to simulate two-phase flow, the numerical solution adopts a semi-implicit solution, the precision is relatively low, and the numerical stability is influenced by the sound velocity coulomb limit value of the material, so the time step length cannot be very large, and the efficiency is low when the program simulates some problems of long transient state;
for some advanced programs that use implicit numerical algorithms, such as: the CARHARE2 program, the numerical solution of which adopts Newton method, because two fluid models contain a large number of constitutive models, the difficulty is very high when constructing 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 Jacobian matrix of the traditional full-implicit solution method, researchers in recent years propose a Newton-Krylov subspace iteration method (JFNK method) without the Jacobian matrix. At present, in the field of nuclear reactor transient safety analysis programs, the JFNK algorithm is mostly applied to the solution of a drift flow model, when the JFNK algorithm is applied to the solution of two fluid models, the preprocessing technology adopts some pure mathematical preprocessing technologies, a Jacobian matrix still needs to be formed, and the difficulty of constructing the Jacobian matrix still exists.
In summary, a semi-implicit solution adopted in a currently common nuclear reactor transient safety analysis program has certain disadvantages, a traditional fully-implicit algorithm has certain difficulties, and although some researches apply the JFNK algorithm to the solution of a drift flow model and a two-fluid model, the preprocessing technology is not concise, so that the model and the algorithm need to be improved to make the nuclear reactor transient safety analysis program 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, a full-implicit solution of two fluid models based on a JFNK algorithm is realized, so that a numerical algorithm is not limited by a material sound velocity coulomb value; secondly, a semi-implicit difference-based Jacobian matrix preprocessing technology is established, and the calculation efficiency is improved; finally, a two-phase flow full-hidden numerical method for the transient safety analysis of the nuclear power plant is obtained.
The invention is realized by adopting the following technical scheme:
a two-phase flow full-hidden numerical method for transient safety analysis of a nuclear power plant comprises the following steps:
the first step is as follows: dividing control bodies according to the operation parameters and the pipeline structure parameters of the nuclear power plant system, and establishing a nonlinear equation set of two-fluid two-phase flow models of all the control bodies:
Figure BDA0003401989850000021
wherein the content of the first and second substances,
Figure BDA0003401989850000022
is the residual form of the vapor phase momentum equation at the inlet boundary i-1/2 of control volume i;
Figure BDA0003401989850000023
is the residual form of the liquid phase momentum equation at the inlet boundary i-1/2 of the control volume i;
Figure BDA0003401989850000024
is the residual error form of the vapor phase energy equation of the control body i;
Figure BDA0003401989850000025
controlling a residual form of a liquid phase energy equation of the body i;
Figure BDA0003401989850000026
controlling a residual error form of a vapor phase mass conservation equation of the body i; fP,iControlling a residual form of a mixed phase mass conservation equation of the body i;
Figure BDA0003401989850000027
is the residual form of the vapor phase momentum equation at the exit boundary i +1/2 of the control volume i;
Figure BDA0003401989850000031
is the residual form of the liquid phase momentum equation at the exit boundary i +1/2 of the control volume i;
the second step is that: solving a nonlinear equation set established in the kth Newton iteration step at the (n + 1) th time step by adopting a JFNK method to obtain the gas-phase internal energy, the liquid-phase internal energy, the void fraction, the pressure, the gas-phase speed and the liquid-phase speed of all control bodies at the (n + 1) th time step, wherein the JFNK method comprises the following steps: a preprocessing matrix M, Krylov subspace iteration and a newton iteration are constructed.
The further improvement of the invention is that the conservation equation of the two-fluid two-phase flow model in the first step is dispersed based on the staggered grid, the time backward difference and the space first-order windward difference, and the obtained full-implicit discrete equation of the control body i and the outlet boundary i +1/2 is as follows:
the k-phase energy equation of the control body i, wherein k ═ g represents the vapor phase and k ═ l represents the liquid phase:
Figure BDA0003401989850000032
vapor phase mass equation of control volume i
Figure BDA0003401989850000033
Mixed phase mass equation of control volume i
Figure BDA0003401989850000034
Equation of k-phase momentum for the control volume i exit boundary i + 1/2:
Figure BDA0003401989850000035
in the formula: Δ x is the control body length; Δ t is the time step of transient analysis; a is the sectional area of the control body; p is the control body pressure; alpha is alphak、ρk、UkRespectively representing the volume fraction, density and internal energy of the k phase; gamma-shapedi,k、Qw,k、Qi,k、Ew,k、Ei,kRespectively representing the heat transfer mass of the k-phase unit volume, the heat exchange quantity of the control body and the wall surface, the gas-liquid interphase heat exchange quantity of the control body, the heat transfer quantity caused by wall surface mass transfer and the heat transfer quantity caused by interphase mass transfer; vk、Fg,k、Fw,k、Fc,k、Fi,k、MkRespectively representing the speed of a k phase, wall gravity, wall friction, local resistance, interphase friction and momentum transfer caused by mass transfer;
Figure BDA0003401989850000041
the amount of the boundary of the equal control body is assigned through a first-order windward difference; k ═ g denotes the vapor phase and k ═ l denotes the liquid phase.
In a further development of the invention, in a second step of the method, the system of nonlinear equations established in the first step is first linearized:
J(xk)·δxk＝-F(xk)
wherein, J (x)k) Jacobian matrix for the kth Newton iteration, deltaxkFor the kth Newton iterationSolving for delta, F (x)k) Set of non-linear equations for the kth Newton iteration, xkIndependent variable for k-th Newton iteration
Figure BDA0003401989850000042
Wherein the content of the first and second substances,
Figure BDA0003401989850000043
controlling the vapor velocity at the inlet boundary i-1/2 of volume i for the kth Newton iteration;
Figure BDA0003401989850000044
controlling the liquid phase velocity at the inlet boundary i-1/2 of volume i for the kth Newton iteration;
Figure BDA0003401989850000045
controlling the vapor phase internal energy of the body i for the kth Newton iteration;
Figure BDA0003401989850000046
controlling the liquid phase internal energy of the body i for the kth Newton iteration;
Figure BDA0003401989850000047
controlling the vapor phase void fraction of the body i for the kth Newton iteration; pi kControlling the pressure of the volume i for the kth Newton iteration;
Figure BDA0003401989850000048
controlling the vapor phase velocity at the exit boundary i +1/2 of volume i for the kth Newton iteration;
Figure BDA0003401989850000049
the liquid phase velocity at the exit boundary i +1/2 of volume i is controlled for the kth newton iteration.
The further improvement of the invention is that the linear equation set is preprocessed, and the linear equation set becomes:
(J(xk)M-1)·(Mδxk)＝-F(xk)
the construction method of the preprocessing matrix M comprises the following steps:
based on staggered grids and finite volume difference, semi-implicit discrete is carried out on two-phase flow equations of two fluids to obtain a nonlinear equation set
Figure BDA0003401989850000051
Wherein the content of the first and second substances,
Figure BDA0003401989850000052
the residual error form of the semi-implicit discrete vapor phase momentum equation at the inlet boundary i-1/2 of the control body i;
Figure BDA0003401989850000053
the residual error form of the liquid phase momentum equation at the inlet boundary i-1/2 of the control body i after semi-implicit dispersion;
Figure BDA0003401989850000054
the method is a residual error form after semi-implicit dispersion of a vapor phase energy equation of a control body i;
Figure BDA0003401989850000055
the method is a residual error form after semi-implicit dispersion of a liquid phase energy equation of a control body i;
Figure BDA0003401989850000056
the residual error form is a residual error form after semi-implicit dispersion of a vapor phase mass conservation equation of a control body i; fP′,iThe residual error is a residual error form after semi-implicit dispersion of a mixed phase mass conservation equation of a control body i;
Figure BDA0003401989850000057
is a residual error form after semi-implicit dispersion of a vapor phase momentum equation at an outlet boundary i +1/2 of the control body i;
Figure BDA0003401989850000058
is the exit boundary i + of control body i1/2, a residual form of the liquid phase momentum equation after semi-implicit discretization;
the elements of the d-th row and the e-th column of the preprocessing matrix M are as follows:
Figure BDA0003401989850000059
wherein F' (x)k)dFor the d-th equation in the system of linear equations F' (x),
Figure BDA00034019898500000510
to solve for variable xkThe e-th solution variable in (1).
The further improvement of the invention is that the linear equation system after pretreatment is solved by using a Krylov subspace iteration method, which specifically comprises the following steps:
firstly, establishing a standard orthogonal basis vector Vm and a Shanghaineberg matrix of a Krylov subspace Km with the dimension of m through an Arnold process
Figure BDA00034019898500000511
Namely:
Vm＝(v1,v2,...,vm)
Figure BDA0003401989850000061
Figure BDA0003401989850000062
Figure BDA0003401989850000063
second, the preprocessed set of equations (J (x))k)M-1)·(Mδxk)＝-F(xk) Conversion into a new system of equations
Figure BDA0003401989850000064
Wherein r is(0)＝-F(xk)-A·xk，A＝J(xk)M-1M is the dimension of the Krylov subspace, zm∈RnTo solve for the variables, z(m)＝Vm·Y＝(v1,v2,...,vm)·(y1,y2,...,ym)T，yi(i is 1, m) is a base vector viCorresponding base, β ═ r(0)||，
Figure BDA0003401989850000065
New system of equations
Figure BDA0003401989850000066
Is an overdetermined equation set with m unknowns and m +1 equations, Y can be obtained by solving by using QR decomposition based on HausHold transform, and z is further calculatedm。
The further improvement of the invention is that the calculation of the dot product of the Jacobian matrix and the vector involved in the Krylov subspace iteration method is totally calculated as follows through difference approximation:
Figure BDA0003401989850000067
wherein epsilon is a difference step length or a disturbance parameter, and x is solved according to the kth stepkVector vkAnd machine accuracy εmachDetermining
Figure BDA0003401989850000068
The invention is further improved in that the convergence conditions of the Krylov subspace iteration are two:
1)
Figure BDA0003401989850000069
2)
Figure BDA0003401989850000071
if either convergence condition is satisfied, the Krylov subspace iteration converges, the system of equations (J (x)k)M-1)·(Mδxk)＝-F(xk) Is solved as δ xk＝M-1(x(k)+zm) (ii) a Otherwise, the dimension m of the Krylov subspace becomes m +1, and the equation set (J (x) is solved againk)M-1)·(Mδxk)＝-F(xk)。
A further improvement of the invention is that the solution δ x of the system of linear equations resulting from the iterative convergence of the Krylov subspacekUpdate the solution x for the kth Newton iterationk
xk+1＝xk+δxk
Judging whether Newton iteration converges, wherein the Newton iteration convergence conditions have three conditions:
1)||F(xk+1)||≤1.0-3；
2)
Figure BDA0003401989850000072
3)
Figure BDA0003401989850000073
wherein, nvol is the total number of control bodies in the nuclear reactor system, and njun is the total number of takeover tubes in the nuclear reactor system;
if the k-th Newton iteration does not satisfy all the above conditions, the Newton iteration is considered not converged, k is k +1, and x is used ask+1As xkPerforming Krylov subspace iteration again; otherwise Newton's iteration converges, xk+1For solutions at the n +1 th time step, according to xk +1The vapor phase internal energy, liquid phase internal energy, void fraction, pressure, vapor phase velocity and liquid phase velocity of all the control bodies in the (n + 1) th time step can be obtained.
Compared with the prior art, the invention has at least the following beneficial technical effects:
1: because a semi-implicit solution method is adopted in a common nuclear reactor transient safety analysis program, the stability of a numerical method is limited by the sound velocity coulomb value of a material, the time step cannot be very large, and the efficiency is low when a slow transient problem is calculated. The invention uses a totally implicit numerical method of two fluid models based on JFNK algorithm, the method can accurately and efficiently simulate the problem of two-phase flow in a nuclear reactor system, and compared with the traditional Newton method, the method can save the difficulty of writing Jacobian matrix.
2: because the JFNK algorithm is applied to the research of two fluid models at present, the preprocessing technology mostly adopts the preprocessing technology based on the pure mathematics built in the PETSc program package, and a Jacobian matrix based on a total-hidden numerical method needs to be written out. The invention uses the semi-implicit difference-based Jacobian matrix preprocessing technology, only needs to write out the Jacobian matrix based on the semi-implicit numerical method, greatly reduces the difficulty of the numerical method, and the preprocessing matrix has physical significance.
Drawings
Fig. 1 is a schematic flow chart of a two-phase flow implicit numerical method for transient safety analysis of a nuclear power plant according to the present invention.
Fig. 2 is a comparison graph of the void fraction distribution calculated by the total hidden numerical method and the results of the experiment value and the semi-hidden numerical method, where 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 the embodiments and features of the embodiments may be combined with each other without conflict.
The following takes a one-dimensional vertical pipeline flow heat transfer problem similar to the flow heat transfer of a nuclear power plant reactor core as an example, and the implementation steps of the two-phase flow total-hidden numerical method for the transient safety analysis of the nuclear power plant provided by the invention are described in detail with reference to fig. 1.
The first step is as follows: 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
Figure BDA0003401989850000091
Wherein the content of the first and second substances,
Figure BDA0003401989850000092
is the residual form of the vapor phase momentum equation at the inlet boundary i-1/2 of control volume i;
Figure BDA0003401989850000093
is the residual form of the liquid phase momentum equation at the inlet boundary i-1/2 of the control volume i;
Figure BDA0003401989850000094
is the residual error form of the vapor phase energy equation of the control body i;
Figure BDA0003401989850000095
controlling a residual form of a liquid phase energy equation of the body i;
Figure BDA0003401989850000096
controlling a residual error form of a vapor phase mass conservation equation of the body i; fP,iThe residual form of the mixed (vapor + liquid) mass conservation equation for control volume i;
Figure BDA0003401989850000097
is the residual form of the vapor phase momentum equation at the exit boundary i +1/2 of the control volume i;
Figure BDA0003401989850000098
to control the liquid phase momentum at the outlet boundary i +1/2 of volume iResidual form of the program;
the second step is that: discretizing the conservation equation of the two-fluid two-phase flow model in the first step based on the staggered grids, the time backward difference and the 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:
the energy equation of the k (k ═ g represents the vapor phase and k ═ l represents the liquid phase) phase of the control body i:
Figure BDA0003401989850000099
vapor phase mass equation of control volume i
Figure BDA00034019898500000910
Mass equation of mixed phase (vapor phase + liquid phase) of control body i
Figure BDA00034019898500000911
The equation for the phase momentum of k (k ═ g for the vapor phase and k ═ l for the liquid phase) at the outlet boundary i +1/2 of control body i:
Figure BDA0003401989850000101
in the formula: Δ x is the control body length; Δ t is the time step of transient analysis; a is the sectional area of the control body; p is the control body pressure; alpha is alphak、ρk、UkRespectively representing the volume fraction, density and internal energy of the k phase; gamma-shapedi,k、Qw,k、Qi,k、Ew,k、Ei,kRespectively representing the heat transfer mass of the k-phase unit volume, the heat exchange quantity of the control body and the wall surface, the gas-liquid interphase heat exchange quantity of the control body, the heat transfer quantity caused by wall surface mass transfer and the heat transfer quantity caused by interphase mass transfer; vk、Fg,k、Fw,k、Fc,k、Fi,k、MkRespectively representing the speed of a k phase, wall gravity, wall friction, local resistance, interphase friction and momentum transfer caused by mass transfer;
Figure BDA0003401989850000102
the amount of the boundary of the equal control body is assigned through a first-order windward difference; k ═ g denotes the vapor phase and k ═ l denotes the liquid phase.
And forming a nonlinear equation system by using the 20 control bodies and 120 equations of the inlet-outlet boundary of each control body:
Figure BDA0003401989850000103
the solution variables are:
x＝(…,Vg,i-1/2,Vl,i-1/2,Ug,i,Ul,i,αg,i,Pi,Vg,i+1/2,Vl,i+1/2,…)T
wherein, Vg,i-1/2To control the vapor velocity at the inlet boundary i-1/2 of volume i; vl,i-1/2To control the liquid phase velocity at the inlet boundary i-1/2 of volume i; u shapeg,iTo control the gas phase internal energy of the body i; u shapel,iTo control the liquid phase internal energy of the body i; alpha is alphag,iControlling the gas phase void fraction of the body i; piTo control the pressure of the body i; vg,i+1/2To control the vapor velocity at the exit boundary i +1/2 of volume i; vl,i+1/2To control the liquid phase velocity at the outlet boundary i +1/2 of volume i;
the third step: linearizing the set of nonlinear equations established in the second step:
J(xk)·δxk＝-F(xk)
wherein, J (x)k) Jacobian matrix for the kth Newton iteration, deltaxkIncremental solution variables, F (x), for the kth Newton iterationk) Set of non-linear equations for the kth Newton iteration, xkIndependent variable for k-th Newton iteration
The fourth step: preprocessing the linear equation set formed in the third step, wherein the construction method of the preprocessing matrix M comprises the following steps:
based on staggered grids and finite volume difference, semi-implicit discrete is carried out on two-phase flow equations of two fluids to obtain a nonlinear equation set
Figure BDA0003401989850000111
Wherein the content of the first and second substances,
Figure BDA0003401989850000112
the residual error form of the semi-implicit discrete vapor phase momentum equation at the inlet boundary i-1/2 of the control body i;
Figure BDA0003401989850000113
the residual error form of the liquid phase momentum equation at the inlet boundary i-1/2 of the control body i after semi-implicit dispersion;
Figure BDA0003401989850000114
the method is a residual error form after semi-implicit dispersion of a vapor phase energy equation of a control body i;
Figure BDA0003401989850000115
the method is a residual error form after semi-implicit dispersion of a liquid phase energy equation of a control body i;
Figure BDA0003401989850000116
the residual error form is a residual error form after semi-implicit dispersion of a vapor phase mass conservation equation of a control body i; fP′,iThe residual error form is a residual error form after semi-implicit dispersion of a mixed (vapor phase + liquid phase) mass conservation equation of a control body i;
Figure BDA0003401989850000117
is a residual error form after semi-implicit dispersion of a vapor phase momentum equation at an outlet boundary i +1/2 of the control body i;
Figure BDA0003401989850000118
residual after semi-implicit discretization of liquid phase momentum equation at exit boundary i +1/2 for control volume iA difference form;
the elements of the d-th row and the e-th column of the preprocessing matrix M are as follows:
Figure BDA0003401989850000119
wherein F' (x)k)dFor the d-th equation in the system of linear equations F' (x),
Figure BDA00034019898500001110
to solve for variable xkThe e-th solution variable in (1).
After preprocessing, the linear equation system becomes:
(J(xk)M-1)·(Mδxk)＝-F(xk)
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:
firstly, establishing a standard orthogonal basis vector Vm and a Shanghaineberg matrix of a Krylov subspace Km with the dimension of m (m is 1 in the first iteration) through an Arnold process
Figure BDA0003401989850000121
Namely:
Vm＝(v1,v2,...,vm)
Figure BDA0003401989850000122
Figure BDA0003401989850000123
Figure BDA0003401989850000124
second, the preprocessed set of equations (J (x))k)M-1)·(Mδxk)＝-F(xk) Conversion into a new system of equations
Figure BDA0003401989850000125
Wherein r is(0)＝-F(xk)-A·xk，A＝J(xk)M-1M is the dimension of the Krylov subspace, zm∈RnTo solve for the variables, z(m)＝Vm·Y＝(v1,v2,...,vm)·(y1,y2,...,ym)T，yi(i is 1, m) is a base vector viCorresponding base, β ═ r(0)||，
Figure BDA0003401989850000126
New system of equations
Figure BDA0003401989850000127
Is an overdetermined equation set with m unknowns and m +1 equations, Y can be obtained by solving by using QR decomposition based on HausHold transform, and z is further calculatedm。
The dot product calculation of the Jacobian matrix and the vector involved in the fifth step of the Krylov subspace iteration method is all calculated by difference approximation as follows:
Figure BDA0003401989850000128
wherein epsilon is a difference step length or a disturbance parameter, and x is solved according to the kth stepkVector vkAnd machine accuracy εmachDetermining
Figure BDA0003401989850000131
For a typical 64-bit computer,. epsilonmach＝10-8。
And a sixth step: for the fifth oneZ of the band of step solutionmJudging whether the Krylov subspace iteration converges, wherein two convergence conditions are as follows:
1)
Figure BDA0003401989850000132
2)
Figure BDA0003401989850000133
if either convergence condition is satisfied, the Krylov subspace iteration converges, the system of equations (J (x)k)M-1)·(Mδxk)＝-F(xk) Is solved as δ xk＝M-1(x(k)+zm) (ii) a Otherwise, changing the dimension m of the Krylov subspace to m +1, returning to the fifth step again, and solving the equation set (J (x)k)M-1)·(Mδxk)＝-F(xk) Solution z ofm+1；
The seventh step, the solution delta x of the linear equation set obtained by the iterative convergence of the Krylov subspace in the sixth stepkUpdate the solution x for the kth Newton iterationk+1＝xk+δxkAnd judging whether Newton iteration converges. The newton iterative convergence condition has three:
1)||F(xk+1)||≤1.0-3；
2)
Figure BDA0003401989850000134
3)
Figure BDA0003401989850000135
if the k-th Newton iteration does not satisfy all the above conditions, the Newton iteration is considered not converged, k is k +1, and x is used ask+1As xkRepeating the fifth, sixth and seventh steps; otherwise Newton's iteration converges, xk+1For solutions at the n +1 th time step, according to xk+1The vapor phase internal energy, liquid phase internal energy, void fraction, pressure, vapor phase velocity and liquid phase velocity of all the control bodies in the (n + 1) th time step can be obtained.
Through the above process, the numerical calculation result is shown in fig. 2. The comparison of the calculation results of a plurality of groups of calculation conditions shows that the calculation result obtained by the total-hidden numerical algorithm in the invention is better in accordance with the experimental result and is close to the calculation result of the semi-hidden numerical algorithm in precision. In the numerical simulation of the experimental conditions, limited by the Counrant value, in order to obtain a stable time value calculation result, the time step of the explicit or semi-implicit numerical algorithm of the experimental conditions does not exceed 0.008s, generally goes 0.005s or less, and the time step Δ t in the fully-implicit numerical method is 0.01 s. Therefore, the time step of the fully-implicit numerical method is not limited by the CoUnrant value, and the numerical calculation efficiency of the transient safety analysis of the nuclear power plant can be improved.
Although the invention has been described in detail hereinabove with respect to a general description and specific embodiments thereof, it will be apparent to those skilled in the art that modifications or improvements may be made thereto based on the invention. Accordingly, such modifications and improvements are intended to be within the scope of the invention as claimed.

Claims (8)

1. A two-phase flow full-hidden numerical method for transient safety analysis of a nuclear power plant is characterized by comprising the following steps:
the first step is as follows: dividing control bodies according to the operation parameters and the pipeline structure parameters of the nuclear power plant system, and establishing a nonlinear equation set of two-fluid two-phase flow models of all the control bodies:
Figure FDA0003401989840000011
wherein the content of the first and second substances,
Figure FDA0003401989840000012
is the residual form of the vapor phase momentum equation at the inlet boundary i-1/2 of control volume i;
Figure FDA0003401989840000013
is the residual form of the liquid phase momentum equation at the inlet boundary i-1/2 of the control volume i;
Figure FDA0003401989840000014
is the residual error form of the vapor phase energy equation of the control body i;
Figure FDA0003401989840000015
controlling a residual form of a liquid phase energy equation of the body i;
Figure FDA0003401989840000016
controlling a residual error form of a vapor phase mass conservation equation of the body i; fP,iControlling a residual form of a mixed phase mass conservation equation of the body i;
Figure FDA0003401989840000017
is the residual form of the vapor phase momentum equation at the exit boundary i +1/2 of the control volume i;
Figure FDA0003401989840000018
is the residual form of the liquid phase momentum equation at the exit boundary i +1/2 of the control volume i;
the second step is that: solving a nonlinear equation set established in the kth Newton iteration step at the (n + 1) th time step by adopting a JFNK method to obtain the gas-phase internal energy, the liquid-phase internal energy, the void fraction, the pressure, the gas-phase speed and the liquid-phase speed of all control bodies at the (n + 1) th time step, wherein the JFNK method comprises the following steps: a preprocessing matrix M, Krylov subspace iteration and a newton iteration are constructed.
2. The two-phase flow total-implicit numerical method for the transient safety analysis of the nuclear power plant according to claim 1, wherein the conservation equation of the two-phase flow model in the first step is discretized based on the staggered grid, the time backward difference and the first-order spatial windward difference, and the total-implicit discrete equation of the control body i and the outlet boundary i +1/2 is obtained as follows:
the k-phase energy equation of the control body i, wherein k ═ g represents the vapor phase and k ═ l represents the liquid phase:
Figure FDA0003401989840000021
vapor phase mass equation of control volume i
Figure FDA0003401989840000022
Mixed phase mass equation of control volume i
Figure FDA0003401989840000023
Equation of k-phase momentum for the control volume i exit boundary i + 1/2:
Figure FDA0003401989840000024
in the formula: Δ x is the control body length; Δ t is the time step of transient analysis; a is the sectional area of the control body; p is the control body pressure; alpha is alphak、ρk、UkRespectively representing the volume fraction, density and internal energy of the k phase; gamma-shapedi,k、Qw,k、Qi,k、Ew,k、Ei,kRespectively representing the heat transfer mass of the k-phase unit volume, the heat exchange quantity of the control body and the wall surface, the gas-liquid interphase heat exchange quantity of the control body, the heat transfer quantity caused by wall surface mass transfer and the heat transfer quantity caused by interphase mass transfer; vk、Fg,k、Fw,k、Fc,k、Fi,k、MkRespectively representing the speed of a k phase, wall gravity, wall friction, local resistance, interphase friction and momentum transfer caused by mass transfer;
Figure FDA0003401989840000025
the amount of the boundary of the equal control body is assigned through a first-order windward difference; k ═g represents a vapor phase and k ═ l represents a liquid phase.
3. The two-phase flow total-hidden numerical method for the transient safety analysis of the nuclear power plant as claimed in claim 1, wherein the second step of the method is to linearize the nonlinear equation set established in the first step:
J(xk)·δxk＝-F(xk)
wherein, J (x)k) Jacobian matrix for the kth Newton iteration, deltaxkIncremental solution variables, F (x), for the kth Newton iterationk) Set of non-linear equations for the kth Newton iteration, xkIndependent variable for k-th Newton iteration
Figure FDA0003401989840000031
Wherein the content of the first and second substances,
Figure FDA0003401989840000032
controlling the vapor velocity at the inlet boundary i-1/2 of volume i for the kth Newton iteration;
Figure FDA0003401989840000033
controlling the liquid phase velocity at the inlet boundary i-1/2 of volume i for the kth Newton iteration;
Figure FDA0003401989840000034
controlling the vapor phase internal energy of the body i for the kth Newton iteration;
Figure FDA0003401989840000035
controlling the liquid phase internal energy of the body i for the kth Newton iteration;
Figure FDA0003401989840000036
controlling the vapor phase void fraction of the body i for the kth Newton iteration; pi kControlling the pressure of the volume i for the kth Newton iteration;
Figure FDA0003401989840000037
controlling the vapor phase velocity at the exit boundary i +1/2 of volume i for the kth Newton iteration;
Figure FDA0003401989840000038
the liquid phase velocity at the exit boundary i +1/2 of volume i is controlled for the kth newton iteration.
4. The two-phase flow total-hidden numerical method for the transient safety analysis of the nuclear power plant as claimed in claim 3, wherein the linear equation set is preprocessed, and the linear equation set becomes:
(J(xk)M-1)·(Mδxk)＝-F(xk)
the construction method of the preprocessing matrix M comprises the following steps:
based on staggered grids and finite volume difference, semi-implicit discrete is carried out on two-phase flow equations of two fluids to obtain a nonlinear equation set
Figure FDA0003401989840000039
Wherein the content of the first and second substances,
Figure FDA00034019898400000310
the residual error form of the semi-implicit discrete vapor phase momentum equation at the inlet boundary i-1/2 of the control body i;
Figure FDA00034019898400000311
the residual error form of the liquid phase momentum equation at the inlet boundary i-1/2 of the control body i after semi-implicit dispersion;
Figure FDA00034019898400000312
the method is a residual error form after semi-implicit dispersion of a vapor phase energy equation of a control body i;
Figure FDA00034019898400000313
the method is a residual error form after semi-implicit dispersion of a liquid phase energy equation of a control body i;
Figure FDA00034019898400000314
the residual error form is a residual error form after semi-implicit dispersion of a vapor phase mass conservation equation of a control body i; f'P,iThe residual error is a residual error form after semi-implicit dispersion of a mixed phase mass conservation equation of a control body i;
Figure FDA0003401989840000041
is a residual error form after semi-implicit dispersion of a vapor phase momentum equation at an outlet boundary i +1/2 of the control body i;
Figure FDA0003401989840000042
the residual error form of the liquid phase momentum equation at the outlet boundary i +1/2 of the control body i after semi-implicit dispersion;
the elements of the d-th row and the e-th column of the preprocessing matrix M are as follows:
Figure FDA0003401989840000043
wherein F' (x)k)dFor the d-th equation in the system of linear equations F' (x),
Figure FDA0003401989840000044
to solve for variable xkThe e-th solution variable in (1).
5. The two-phase flow total-hidden numerical method for the transient safety analysis of the nuclear power plant as claimed in claim 4, wherein the preprocessed linear equation set is solved by using a Krylov subspace iteration method, which specifically comprises:
firstly, establishing a standard orthogonal basis vector Vm and a Shanghaineberg matrix of a Krylov subspace Km with the dimension of m through an Arnold process
Figure FDA0003401989840000045
Namely:
Vm＝(v1,v2,...,vm)
Figure FDA0003401989840000046
Figure FDA0003401989840000047
Figure FDA0003401989840000048
second, the preprocessed set of equations (J (x))k)M-1)·(Mδxk)＝-F(xk) Conversion into a new system of equations
Figure FDA0003401989840000049
Wherein r is(0)＝-F(xk)-A·xk，A＝J(xk)M-1M is the dimension of the Krylov subspace, zm∈RnTo solve for the variables, z(m)＝Vm·Y＝(v1,v2,...,vm)·(y1,y2,...,ym)T，yi(i is 1, m) is a base vector viCorresponding substrate, β ═ r: (a), (b), and (c)0)||，
Figure FDA0003401989840000051
New system of equations
Figure FDA0003401989840000052
Is an overdetermined equation with m unknowns and m +1 equationsSet of Y obtained by solving using HausHold transform based QR decomposition, and then calculating zm。
6. The two-phase flow total-hidden numerical method for the transient safety analysis of the nuclear power plant according to claim 5, wherein the calculation of the dot product of the Jacobian matrix and the vector involved in the Krylov subspace iteration method is all calculated by difference approximation as:
Figure FDA0003401989840000053
wherein epsilon is a difference step length or a disturbance parameter, and x is solved according to the kth stepkVector vkAnd machine accuracy εmachDetermining
Figure FDA0003401989840000054
7. The two-phase flow total-hidden numerical method for the transient safety analysis of the nuclear power plant according to claim 5, wherein the convergence condition of the Krylov subspace iteration is two:
1)
Figure FDA0003401989840000055
2)
Figure FDA0003401989840000056
if either convergence condition is satisfied, the Krylov subspace iteration converges, the system of equations (J (x)k)M-1)·(Mδxk)＝-F(xk) Is solved as δ xk＝M-1(x(k)+zm) (ii) a Otherwise, the dimension m of the Krylov subspace becomes m +1, and the equation set (J (x) is solved againk)M-1)·(Mδxk)＝-F(xk)。
8. The two-phase flow total-hidden numerical method for the transient safety analysis of the nuclear power plant as claimed in claim 7, wherein the solution δ x of the linear equation system obtained by the iterative convergence of the Krylov subspacekUpdate the solution x for the kth Newton iterationk
xk+1＝xk+δxk
Judging whether Newton iteration converges, wherein the Newton iteration convergence conditions have three conditions:
1)||F(xk+1)||≤1.0-3；
2)
Figure FDA0003401989840000061
3)
Figure FDA0003401989840000062
wherein, nvol is the total number of control bodies in the nuclear reactor system, and njun is the total number of takeover tubes in the nuclear reactor system;
if the k-th Newton iteration does not satisfy all the above conditions, the Newton iteration is considered not converged, k is k +1, and x is used ask+1As xkPerforming Krylov subspace iteration again; otherwise Newton's iteration converges, xk+1For solutions at the n +1 th time step, according to xk+1The vapor phase internal energy, liquid phase internal energy, void fraction, pressure, vapor phase velocity and liquid phase velocity of all the control bodies in the (n + 1) th time step can be obtained.
CN202111501842.1A 2021-12-09 2021-12-09 Two-phase flow total hidden numerical method for transient safety analysis of nuclear power plant Active CN114169203B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111501842.1A CN114169203B (en) 2021-12-09 2021-12-09 Two-phase flow total hidden numerical method for transient safety analysis of nuclear power plant

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111501842.1A CN114169203B (en) 2021-12-09 2021-12-09 Two-phase flow total hidden numerical method for transient safety analysis of nuclear power plant

Publications (2)

Publication Number Publication Date
CN114169203A true CN114169203A (en) 2022-03-11
CN114169203B CN114169203B (en) 2024-01-16

Family

ID=80485086

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111501842.1A Active CN114169203B (en) 2021-12-09 2021-12-09 Two-phase flow total hidden numerical method for transient safety analysis of nuclear power plant

Country Status (1)

Country Link
CN (1) CN114169203B (en)

Citations (4)

* Cited by examiner, † Cited by third party
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
US20200194140A1 (en) * 2018-12-18 2020-06-18 Deep Isolation, Inc. Radioactive waste repository systems and methods
CN112906271A (en) * 2021-02-22 2021-06-04 中国核动力研究设计院 Reactor transient physical thermal full-coupling fine numerical simulation method and system

Patent Citations (4)

* Cited by examiner, † Cited by third party
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
US20200194140A1 (en) * 2018-12-18 2020-06-18 Deep Isolation, Inc. Radioactive waste repository systems and methods
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)

* Cited by examiner, † Cited by third party
Title
孙涛;徐钊;: "核电厂非能动安全壳热量导出系统瞬态模拟程序开发研究", 核动力工程, no. 03 *

Also Published As

Publication number Publication date
CN114169203B (en) 2024-01-16

Similar Documents

Publication Publication Date Title
CN107066745B (en) Method for obtaining three-dimensional neutron flux density distribution in fast neutron reactor core transient process
CN107871046B (en) Simulation method for spray blending in low-temperature propellant storage tank
Zhang et al. A class of DG/FV hybrid schemes for conservation law IV: 2D viscous flows and implicit algorithm for steady cases
CN114065662B (en) Method suitable for rapidly predicting airfoil flow field with variable grid topology
Salko Jr et al. CTF 4.0 Theory Manual
Rumpfkeil et al. Efficient hessian calculations using automatic differentiation and the adjoint method with applications
CN114707394A (en) Reactor transient neutron flux simulation method considering deformation effect under fixed grid
Chaudri et al. Development and validation of a fast sub-channel code for LWR multi-physics analyses
CN114169203B (en) Two-phase flow total hidden numerical method for transient safety analysis of nuclear power plant
CN110705184B (en) Virtual volume momentum source method for fine numerical solution of reactor core
Yoon et al. A multiscale and multiphysics PWR safety analysis at a subchannel scale
Jeong et al. A semi-implicit numerical scheme for a transient two-fluid three-field model on an unstructured grid
Celik et al. Numerical simulation of circulation in gas-liquid column reactors: isothermal, bubbly, laminar flow
Gallorini et al. An adjoint‐based solver with adaptive mesh refinement for efficient design of coupled thermal‐fluid systems
Tentner et al. Modeling of two-phase flow in a BWR fuel assembly with a highly-scalable CFD code
Cai et al. Derivative-free level-set-based multi-objective topology optimization of flow channel designs using lattice Boltzmann method
Perez-Castro et al. PhotoBioLib: A Modelica library for modeling and simulation of large-scale photobioreactors
Ge et al. Development and application of a new high-efficiency sparse linear system solver in the thermal-hydraulic system analysis code
CN116579205A (en) Calculation method for pressurized water reactor nuclear thermal coupling
Bullerwell A Multiphysics Coupling Scheme Between OpenMC and OpenFOAM for High-Fidelity, Flexible Multi-Physics Analysis of Nuclear Reactors
Douglas et al. Multigrid solution of flame sheet problems on serial and parallel computers
Lee et al. Multi-physics approach for nuclear reactor analysis using thermal-hydraulics and neutron kinetics coupling methodology
Rumpfkeil et al. Efficient hessian calculations using automatic differentiation and the adjoint method
Cadinu et al. Relating system-to-CFD coupled code analyses to theoretical framework of a multiscale method
Amoignon Adjoint-based aerodynamic shape optimization

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