Parallel computing method suitable for two-fluid six-equation system

Abstract

The invention belongs to the technical field of fluid calculation, and particularly relates to a parallel calculation method suitable for a two-fluid six-equation system. Dynamically adjusting the time step length based on the Languang number check to ensure the stability; parallel computing generalized source items and storing the generalized source items to a node data container; assembling a mass and energy conservation equation coefficient matrix; partitioning the sparse matrix and inversing, and establishing a partitioning scalar and vector matrix calculation formula; combining the constraint equation, constructing a pressure matrix equation, and numbering loops and nodes in a layered manner by using breadth-first traversal; solving a pressure matrix in parallel; calculating a mass error through the mixed density deviation and verifying the mass error; and carrying out loop iteration on the steps until the set simulation time is reached. According to the method, the solving efficiency of the complex system is improved through the block matrix, parallel computing and error control. The method has the beneficial effects that a two-fluid six-equation parallel computing process framework is established, the computing process is complete and effective, and the solving efficiency is greatly improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of fluid calculation, and particularly relates to a parallel calculation method applicable to a two-fluid six-equation system. Background Art

[0002] The two-fluid six-equation system is a typical non-linear equation system with characteristics such as time transient, spatial coupling, and the combined action of multiple physical factors. In the process of solving, it is necessary to consider the advancement based on the simulation step size in the time dimension, the transport and conservation of mass flow, momentum flow, and energy flow in the spatial dimension, and the influence of physical process phenomena on the equation coefficients. Therefore, the solution calculation process of the two-fluid six-equation system is complex and time-consuming. Summary of the Invention

[0003] The purpose of the present invention is to provide a parallel calculation method applicable to a two-fluid six-equation system, which can realize block parallel calculation of matrix equations, parallel calculation of source term functions, automatic scaling of time steps, and efficient management of the solution memory, accelerate the solution efficiency of the two-fluid six-equation system, and provide technical support for the rapid design verification of thermal-hydraulic systems and the real-time simulation of digital twin models.

[0004] The technical solution of the present invention is as follows: A parallel calculation method applicable to a two-fluid six-equation system includes the following steps:

[0005] Step 1: Conduct Courant number check. Based on the state variable values of each control volume and connecting pipe of the thermal-hydraulic system at the current moment, calculate the Courant limit of the current two-fluid six-equation system. If the current time step is greater than the Courant limit, adjust the current time step to the Courant limit; if the current time step is less than or equal to the Courant limit, do not adjust the current time step.

[0006] Step 2: Based on the current time step and the current state variable values of the control volume and connecting pipe, conduct parallel calculation of the generalized source term function of the basic node unit to obtain the complete data information required for assembling the coefficient matrix of the two-fluid six-equation, and store it in the corresponding data container of the basic node unit.

[0007] Step 3: Read the data information stored in the data container of the basic node unit, and assemble the coefficient matrix of the mass and energy conservation equations and the coefficient matrix of the momentum conservation equation according to the topological structure of the thermal-hydraulic system.

[0008] Step 4: Conduct a block operation on the coefficient matrix of the mass and energy conservation equations and the coefficient matrix of the momentum conservation equation of the two-fluid six-equation system, convert the coefficient matrix with the characteristics of a large sparse matrix into a combination of several small matrices, and further conduct an inversion operation on the coefficient matrix to establish a scalar array block matrix calculation formula for the mass and energy conservation equations based on the control volume and a vector array block matrix calculation formula for the momentum conservation equation based on the connecting pipe.

[0009] Step 5: Combine the constraint equations of the control body pressure and the pipe vapor-liquid flow rate in the scalar array block matrix calculation formula and the constraint equations of the pipe vapor-liquid flow rate and the control body pressure in the vector array block matrix calculation formula to form a pressure matrix equation with only the control body pressure as the solution variable. Furthermore, in the process of assembling the pressure matrix equation, a breadth-first traversal algorithm is used to perform two-level numbering on the thermal hydraulic circuits and basic node units, and the pressure matrix equations of each independent circuit are established according to the circuit number;

[0010] Step 6: Based on the direct solution method, the pressure matrix equations of each independent loop are solved and calculated in parallel to obtain the pressure of each control body of the thermal hydraulic system, and the pressure is stored and updated in the data container of each control body;

[0011] Step 7: Call the vector array block matrix calculation formula of the momentum conservation equation based on the nozzle, take the control body pressure as input, calculate the steam and liquid flow rate of each nozzle of the thermal hydraulic system, and store and update it in the data container of each nozzle;

[0012] Step 8: Call the block matrix calculation formula based on the mass and energy conservation equations of the control body, take the steam-liquid flow rate of each pipe of the thermal-hydraulic system as input, calculate the scalars such as the cavitation fraction, steam phase internal energy, liquid phase internal energy, and non-condensable gas fraction of each control body of the thermal-hydraulic system, and store and update them in the data container of each control body;

[0013] Step 9: Calculate the mass error of the solution of the two-fluid six-equation system, where the mass error is the deviation between the vapor-liquid mixture density obtained by the first-order Taylor expansion and the vapor-liquid mixture density calculated by the physical property state equation, and determine the numerical value of the current mass error and the set mass error upper limit and mass error lower limit;

[0014] Step 10: The simulation time is advanced along the timestamp with a time step length. Steps 1 to 9 are repeated until the current simulation time is greater than or equal to the set simulation time.

[0015] The Courant number check described in step 1 includes the following:

[0016] The Courant number Co is the ratio of the flow distance of the fluid microelement to the length of the local grid unit under the simulation time step, which is defined as:

[0017]

[0018] Where: V is the flow velocity, Δt is the time step, and Δx is the grid unit length;

[0019] When the Courant number Co is less than 1, in one time step, the fluid element cannot directly "cross" the grid cell. When the Courant number Co is greater than 1, in one time step, the fluid element can directly "cross" the grid cell. To ensure the stability and reliability of fluid simulation calculations, the Courant number Co of all grid cells should be less than 1. To prevent the Courant number from being greater than 1 during the fluid simulation process, before each time advancement, based on the condition that the Courant number is not greater than 1 and combined with the flow velocity V at the previous time step n calculate the Courant time of each basic node unit. The Courant time is the time step corresponding to when the Courant number of the basic node unit is equal to 1. When the pre-used time step is greater than the minimum value of the Courant times of each basic node unit in the thermal-hydraulic system, that is, the Courant limit, then change the time step to the Courant limit. The calculation method of the Courant limit is:

[0020]

[0021] where M is the total number of non-boundary nozzles, Δt upc，m and Δt downc，m are the upstream and downstream upwind Courant times of the m-th nozzle, V K，m is the volume of the upstream control volume of the m-th nozzle, A m is the flow area of the m-th nozzle, v g，m and v f，m are the vapor-phase and liquid-phase flow velocities of the m-th nozzle, Δx K，m and Δx L，m are the lengths of the upstream and downstream control volumes of the m-th nozzle; v f，L，m and v g，L，m are the liquid-phase and vapor-phase flow velocities at the center of the downstream control volume of the m-th nozzle.

[0022] The generalized source term function in the parallel calculation of the generalized source term function described in step 2 refers to the source terms in physics such as wall heat transfer calculation, interphase heat transfer calculation, interphase friction calculation, wall friction calculation, and local resistance calculation, the solution variable coefficients and constant terms of the mass, energy, and momentum conservation discrete equations, and the calculation preparation parts required for calculating these source terms, solution variable coefficients, constant terms, etc., including property calculation, nozzle variable calculation based on control volume variables, and control volume variable calculation based on nozzle variables.

[0023] The process of parallel calculation of the generalized source term function includes:

[0024] Step 21: Define the data container corresponding to the basic node unit to store the data information required for the assembly and solution of the two flow velocity and six equation coefficient matrices. The data information is stored and read and written using a certain data structure;

[0025] Step 22: Define the data reading Get(id, VarType) function and data writing Set(id, VarType, inputVal) function for the data container of the basic node unit. The Get function and Set function read and write the data information at the corresponding position of the data container data structure according to the basic node unit number id and variable name VarType;

[0026] Step 23: Use the Get(id, VarType) function and Set(id, VarType, inputVal) function to convert each source term function into a source term function model with the basic node unit number id as the only input variable. The source term function model reads the data stored in the data container of the basic node unit numbered id as the input variable of the source term function through the Get function, calculates the function output result using the mathematical principle formula of the source term function, and then uses the Set function to write the function output result into the memory space corresponding to the source term function output variable in the data container of the basic node unit numbered id;

[0027] Step 24: Define a source term function enumeration table. The source term function enumeration table encompasses all the generalized source term function types for the two-fluid six-equation system calculation. The position serial number in the enumeration table represents the execution order of the source term functions, expressing the mathematical logic of the two-fluid six-equation system source term function calculation. To reduce the system traversal times and the communication times between sub-threads and the main thread during the parallel calculation of the source term function model, introduce source term function major categories. The source term function major category refers to a set of source term functions that are adjacent in execution order in the source term function enumeration table, belong to the same type of basic node unit, and are spatially independent of each other;

[0028] Step 25: According to the classification principle and execution order of the source term function major categories, traverse the basic node units of the thermal-hydraulic system under the same source term function major category, and carry out the source term function model calculation of each basic node unit using the multi-basic node unit parallel calculation method.

[0029] The said Step 25 includes:

[0030] Step 251: Based on the storage and reading / writing of the basic node unit data information in the data container

[0031] The data container refers to the memory space allocated according to a certain data structure. Use a variable enumeration table to define the storage data structure of the parameters, variables, and source term functions of the basic nodes. The variable enumeration table is an array in which all the input and output variables and parameters of all the generalized source term functions of the basic node unit are stored according to a certain data structure. The variable enumeration table contains all the input and output data required for the calculation of the generalized source term functions of the basic node unit;

[0032] After defining the variable enumeration table, define the read function Get(id, VarType) and the write function Set(id, VarType, inputVal). The former reads the variable VarType stored in the data container of the basic node unit with the number id, and the latter writes the externally obtained parameter inputVal to the memory location corresponding to the variable VarType in the data container of the basic node unit with the number id. The variable VarType must be defined in the variable enumeration table. Through the Get(id, VarType) function and the write Set(id, VarType, inputVal) function, direct operation on the memory space corresponding to the parameter variable is achieved;

[0033] Step 252: Definition of the source term function model with the basic node unit number id as the only input variable

[0034] According to the mathematical principle formula of the source term function, using the Get(id, VarType) function and the Set(id, VarType, inputVal) function, convert each source term function into a source term function model with the basic node unit number id as the only input variable;

[0035] Step 253: Multi-threaded parallel calculation based on the major categories of source term functions

[0036] Define the source term function enumeration table. The source term function enumeration table is an array that stores all the generalized source term function types of the two-fluid six-equation system according to a certain data structure. The position serial number in the enumeration table represents the execution order of the source term function. The execution order of the source term function is defined in the source term function enumeration table, expressing the mathematical logic of the source term function calculation of the two-fluid six-equation system. Under the same source term function, traverse the basic node units of the thermal-hydraulic system, and use the multi-basic node unit parallel calculation method to carry out the calculation of the source term function model of each basic node unit.

[0037] The said step 252 includes:

[0038] Step 2521: According to the source term function name, define the input variables, output variables and process variables of the source term function;

[0039] Step 2522: Use Get(id, VarType) to read the data required for the input variables from the data container of the basic node unit with the number id;

[0040] Step 2523: According to the mathematical principle formula of the source term function, obtain the process variable data and output variable data of the function calculation through the input variable data;

[0041] Step 2524: Write the obtained output variable data into the memory location corresponding to the output variable of the source term function in the basic node unit data container numbered id using Set(id, VarType, inputVal).

[0042] The assembly process of the coefficient matrices of the mass and energy conservation equations and the momentum conservation equation in Step 3 includes that in the thermal-hydraulic system model, assuming there are N non-boundary control volumes and M non-boundary nozzles, then there are 5N control volume variables, including pressure, void fraction, specific internal energy of the vapor phase, specific internal energy of the liquid phase, non-condensable gas fraction, and 2M nozzle variables, including vapor phase flow rate and liquid phase flow rate. The discrete forms of the scalar equations such as the mass conservation equation, mass conservation difference equation, vapor phase energy conservation equation, liquid phase energy conservation equation, and non-condensable gas transport equation for N control volumes are sorted out in the following form:

[0043]

[0044] where, a 1 to a 5 are respectively the variable coefficients in the scalar equation, s 1 to s 4 are respectively the variable coefficients in the scalar equation, b is the constant term of the equation. The subscript g represents the vapor phase, the subscript f represents the liquid phase, the subscript L represents the current control volume, the subscript j represents the nozzle number, j + 1 represents the downstream nozzle, the superscript n represents the nth time moment, the superscript n + 1 is the (n + 1)th moment, a 1 to a 5 , s 1 to s 4 , and b are all obtained by calculating the source term function and stored in the data containers of the corresponding control volumes; the scalar equations for N control volumes are combined into matrix form:

[0045] AX = SV + B (6)

[0046] where,

[0047]

[0048] A is a 5N×5N coefficient matrix for solving scalars such as control volume pressure, void fraction, specific internal energy of the vapor phase, specific internal energy of the liquid phase, and non-condensable gas fraction, X is an array of 5N control volume variables, S is a 5N×2M coefficient matrix for solving vectors such as vapor phase flow rate and liquid phase flow rate of the upstream and downstream nozzles of the control volume, V is an array of 2M nozzle variables, and B is an array of 5N constant terms;

[0049] According to the topological structure of the thermal-hydraulic system, the momentum conservation equations and momentum conservation difference equations of the M nozzles are sorted in the following form:

[0050]

[0051] where g 1 、g 2 、w 1 、w 2 are the variable coefficients in the equations respectively, f is the constant term of the equations, and g 、g 1 、g 2 、w 1 、w 2 、f are all obtained by calculating the source term function and stored in the data containers of the corresponding nozzles;

[0052] The momentum conservation equations of the M nozzles are combined into a matrix form:

[0053] GV = WP + F (8)

[0054] where,

[0055] G is a 2M×2M matrix, V is an array of vapor-phase and liquid-phase flow velocity variables of 2M nozzles, W is a 2M×N matrix, P is an array of pressure variables of N control volumes. The coefficient matrices A, S, G, and W are sparse matrices. The row elements of the matrix correspond to the coefficients of all solution variables in a certain conservation equation. When a certain variable is not in this equation, its coefficient is 0. According to the discrete format of the two-fluid six equations, the scalar conservation discrete equations such as the mass, energy, and non-condensable gas of the control volume only contain the variables and parameters of this control volume and the upstream and downstream nozzles, while the momentum conservation discrete equation of the nozzle only contains the variables and parameters of this nozzle and the upstream and downstream control volumes.

[0056] The two-fluid six-equation coefficient matrix blockification and inversion operation in step 4 includes adjusting the variable array to make the variables of the same basic node unit adjacent in the variable array, and performing non-zero element concentration processing on the coefficient matrix of the scalar equation, that is:

[0057]

[0058] (11) can be simplified to:

[0059]

[0060] where: A 1 and A 2 are 5×5 matrices, X 1 and X 2 are basic node 5×1 scalar arrays, S 1and S 2 is a 5×4 matrix, V 1 and V 2 is an array of 2×1 vectors of basic nodes, B 1 and B 2 is an array of 5×1 constants.

[0061] Performing an inverse operation on equation (12), then:

[0062]

[0063] Without loss of generality, under N non-boundary control volumes and M non-boundary takeovers, equation (13) can be written as:

[0064]

[0065] where, X 1 to X N is an array of 5×1 scalars of basic nodes, A 1 to A N is a 5×5 matrix, S 1 to S N is a 5×2M matrix, V 1 to V M is an array of 2×1 vectors of basic nodes, B 1 to B N is an array of 5×1 constants;

[0066] It can be seen from equation (14) that after block processing the scalar equation AX = SV + B, its time complexity changes from O(125N 3 ) to O(125N), and the theoretical calculation cost increases linearly with the number of control volumes;

[0067] Since the element values of the A i (i = 1,..., N) matrices and B i (i = 1,..., N) can all be calculated from the parameter variables stored in the data containers corresponding to the current control volume, while only the columns corresponding to the takeovers connected to the current control volume in S i (i = 1,..., N) have non-zero elements, so and 's calculation is only related to the current control volume and its connected takeovers, meeting the conditions for parallel calculation; at the same time, since A i is a 5×5 matrix, its inverse matrix is easily obtained, and the positions of the non-zero element columns can be quickly confirmed through the topological structure information of the current control volume and the takeovers, further improving the calculation efficiency;

[0068] ​Similarly, by adjusting the variable array to make the variables of the same basic node unit adjacent in the variable array, the non-zero elements of the coefficient matrix of the scalar equation are centralized:

[0069]

[0070] Equation (15) can be simplified to:

[0071]

[0072] Where: G 1 and G 2 are 2×2 matrices, X 1 and X 2 are 2×1 vector arrays of basic node units, W 1 and W 2 are 2×2 matrices, P 1 and P 2 are the pressures of basic node units, F 1 and F 2 are 2×1 constant arrays;

[0073] Performing an inversion operation on Equation (16), then:

[0074]

[0075] Without loss of generality, under N non-boundary control volumes and M non-boundary nozzles, Equation (17) can be written as:

[0076]

[0077] Where, V 1 to V M are 2×1 vector arrays of basic nodes, G 1 to G M are 2×2 matrices, W 1 to W M are 2×N matrices, F 1 to F N are 2×1 constant arrays;

[0078] Similarly, by block-processing the vector matrix GV = WP + F, its time complexity changes from O(4M 2 ) to O(4M), and the computational cost increases linearly with the number of nozzles. and The calculations of are only related to the current nozzle and its upstream and downstream control volumes, meeting the conditions for parallel computing.

[0079] The generation of the pressure matrix equation in step 5 includes, after the two-fluid six-equation coefficient matrix is block-partitioned and inverted, extracting all the matrix element expressions regarding the control volume pressure in equation (14), i.e., X i , and 's first row, and determining the non-zero element distribution of the first row of through the connection topology information between the control volume and the nozzles. The matrix element expression of equation (14) regarding pressure is in the following form:

[0080]

[0081] where j is the inlet nozzle number, m is the total number of inlet nozzles, jp is the outlet nozzle number, n is the total number of outlet nozzles, b i is 's constant term regarding pressure P i , C P，g，j , C P，f，j , C P，f，jp , C P，f，jp are the non-zero elements of the first row of ;

[0082] Furthermore, extract the matrix element expressions regarding the flow velocities of the inlet and outlet nozzles of the i-th control volume in equation (18). Taking the j-th inlet and jp-th outlet nozzles as an example, the following form is obtained:

[0083] V g，j = f V，g，j + C V，g，K，j P K，j + C V，g，i，j P i (20)

[0084] V f，j = f V，f，j + C V，f，K，j P K，j + C V，f，i，j P i (21)

[0085] V g，jp = f V，g，jp + C V，g，i，jp P i + C V，g，L，jp P L，jp (22)

[0086] V f，jp = f V，f，jp + C V，f，i，jp P i + C V，f，L，jp P L，jp (23)

[0087] Among them, f V，g，j and f V，f，j are respectively the first - row and second - row constant terms of V，g，K，j and C V，f，K，j are respectively the first - row and second - row matrix elements of V，g，i，j with respect to the pressure of the upstream control volume of the j - th take - over, C V,f,i,j and C are respectively V,g,jp the first - row and second - row matrix elements of V,f,jp with respect to the pressure of the i - th control volume, f and f V,g,i,jp are respectively V，f，i，jp the first - row and second - row constant terms of and C V，g，L，jp and C V,f,L,jp are respectively the first - row and second - row matrix elements of

[0088] Substituting equations (20) to (23) into equation (19), we can obtain:

[0089]

[0090] The coefficients and constants of equation (24) can all be obtained through the matrix of . Further, traversing by the control - volume number, equation (24) is converted into a matrix equation

[0091] A p P = b p (25)

[0092] A p is a large sparse matrix. The non - zero elements of the i - th (i = 1, …, N) row of the A p matrix only include the i - th column and the column numbers corresponding to the pressures of the non - boundary control volumes that are connected by take - overs to the i - th control volume. The number of its non - zero elements is m + n + 1, where m is the total number of inlet take - overs and n is the total number of outlet take - overs, which is much smaller than the number of control volumes N. The pressures P n+1 of each control volume in the thermal - hydraulic system can be solved through equation (25), and then V n+1 is calculated through equation (18), and finally other scalars except pressure are calculated through equation (14)

[0093] In step 5 described above, it also includes establishing a pressure equation for the thermal - hydraulic system, and its structure is as shown in (26)

[0094]

[0095] Perform a linear transformation on equation (26), that is, adjust the order of the variables in the variable array, and group the control volume variables of the same loop together, as shown in equation (27).

[0096]

[0097] Equation (27) is simplified to:

[0098]

[0099] Where: and are the block matrices of the pressure equations corresponding to Loop 1 and Loop 2 respectively; and are the arrays of pressure variables for Loop 1 and Loop 2; and are the arrays of constant terms for Loop 1 and Loop 2;

[0100] Considering k independent loops, we get:

[0101]

[0102] That is:

[0103] The step 5 includes the following steps:

[0104] Step 51: Use the breadth-first traversal algorithm to number the basic node units of the thermal-hydraulic system, generating a two-layer numbering including loop numbers and node numbers;

[0105] Step 52: Without distinguishing the loop numbers of the basic node units, use the parallel computing method to traverse the node numbers of the control volume and perform the coefficient calculations required for assembling the coefficient matrix of the pressure equation;

[0106] Step 53: Sequentially according to the loop numbers of the basic node units, traverse the node numbers of the basic node units under the loop number, obtain the pressure matrix coefficients corresponding to the basic node units of the control volume type, and assemble the block coefficient matrix of the pressure equation for the loop number.

[0107] The step 51 includes the following steps:

[0108] Step 511: Change the directed edges of the flattened topological structure diagram of the thermal-hydraulic system to undirected edges. The flattened topological structure diagram is a system model that restores the thermal-hydraulic system model constructed by the user to a system model completely composed of basic node unit models and connected in an orderly manner of control volume - takeover - control volume;

[0109] Step 512: Establish a data container A for the graph vertices that includes the undirected graph of the flattened topology of the thermohydraulic system. The data container A stores the basic node unit names of all graph vertices and the names of the model components where the basic node units are located in a queue manner;

[0110] Step 513: Obtain the first graph vertex stored in the data container. The loop number of this graph vertex is marked as 1 and stored in a temporary data container B. The temporary data container B expresses the order of the graph vertex elements stored in a queue form;

[0111] Step 514: Obtain the first graph vertex element in the temporary data container B and traverse all the connected graph vertex elements. It is divided into the following situations:

[0112] Situation 1: The connected graph vertex elements have been numbered, continue traversing until the traversal ends;

[0113] Situation 2: The connected graph vertex elements have not been numbered, store them in the temporary data container B in a queue and continue traversing until the traversal ends.

[0114] Step 515: After the traversal of the connected graph vertex elements ends, delete the first graph vertex element from the temporary data container B and store it in the result data container C;

[0115] Step 516: Record the type of this graph vertex element and the sequence number of the graph vertex element stored in the result data container C under different types. The types of graph vertex elements include control volumes and nozzles. The sequence number of the control volume type graph vertex element stored in the result data container C is used as the node number of the control volume, and the sequence number of the nozzle type graph vertex element stored in the result data container C is used as the node number of the nozzle;

[0116] Step 517: Judge whether this graph vertex element has been marked with a loop number, and it is divided into the following situations:

[0117] Situation 1: This graph vertex element has not been marked with a loop number. The loop number of this graph vertex is marked as the loop number of the graph vertex with the previous sequence number in the result data container C, and then continue with the next step operation;

[0118] Situation 2: This graph vertex element has been marked with a loop number, continue with the next step operation;

[0119] Step 518: Judge whether the temporary data container B is empty, and it is divided into the following situations:

[0120] Situation 1: The temporary data container B is not empty, jump to Step 514;

[0121] Situation 2: The temporary data container B is empty, continue with the next step operation;

[0122] Step 519: Traverse the graph vertex elements of data container A, and determine whether there are graph vertex elements in data container A that have not been numbered by nodes, which are divided into the following two cases:

[0123] Case 1: There are graph vertex elements that have not been numbered by nodes. Arbitrarily select a graph vertex element that has not been numbered by nodes and store it in the temporary data container B. The loop number of this graph vertex element is marked as the loop number of the graph vertex element finally stored in the result data container C + 1, and jump to step 514;

[0124] Case 2: There are no graph vertex elements that have not been numbered by nodes, then the basic node unit numbering ends.

[0125] The specific calculation process in step 6 is as follows:

[0126] Suppose the independent loop has N control volumes. The non-zero elements in the i-th (i = 1,..., N) row of the pressure matrix only include the i-th column and the column numbers corresponding to the pressures of non-boundary control volumes that are connected by pipelines to the i-th control volume. The direct solution method based on the sparse LU decomposition method is used as the solution algorithm for the pressure matrix equation. The sparse LU decomposition method aims to obtain the elementary transformation matrix Q, the lower triangular matrix L, and the upper triangular matrix U, such that:

[0127] QA = LU (31)

[0128] Among them, in the LU decomposition process, the elementary transformation matrix Q extracts the element with the largest absolute value in each column through elementary row transformation as the pivot element for elimination. The pressure equation block matrix (30) is converted to:

[0129]

[0130] is the elementary transformation matrix, is to reorder the elements of the array and then, according to the structural characteristics of the lower triangular matrix and the upper triangular matrix , obtain the independent loop pressure distribution

[0131] The calculation of the vapor-liquid flow rate in the pipeline in step 7 is as follows: Call the formula (18) for block matrix calculation of the vector array of the momentum conservation equation based on the pipeline, and use the control volume pressure P n+1 as the input to calculate the flow rate of each pipeline in the thermal-hydraulic system Store and update it in the data container of each pipeline;

[0132]

[0133] Wherein: the superscript n+1 represents the calculated value at the current moment, n represents the calculated value at the previous moment, and M is the number of non-boundary takeovers of the thermal-hydraulic system.

[0134] In step 8 described above, the void fraction, vapor-liquid specific internal energy, and non-condensable gas fraction of the control volume are calculated as follows: Call the scalar array block matrix calculation formula (14) of the mass and energy conservation equations based on the control volume, and use the vapor-liquid flow velocity V of the thermal-hydraulic system takeovers n+1 as the input to calculate the void fraction of the control volume of the thermal-hydraulic system vapor specific internal energy liquid specific internal energy non-condensable gas fraction Store the updated data in the data containers of each control volume;

[0135]

[0136] In step 9 described above, the mass error of the two-fluid six-equation system solution is the deviation between the vapor-liquid mixture density obtained through first-order Taylor expansion and the vapor-liquid mixture density calculated through the physical property state equation. Its expression is:

[0137] ε m =max(ε mlocal , ε mtotal ) (35)

[0138]

[0139] ρ i =α g，i ρ g，i +(1-α g，i )ρ f，i (39)

[0140]

[0141] Wherein, and are the vapor density and liquid density obtained through first-order Taylor expansion for the i-th control volume; ρ g，i and ρ f，i are the vapor density and liquid density obtained through the physical property state equation or interpolation table for the i-th control volume; k is a factor greater than 1; N is the number of non-boundary control volumes;

[0142] Set the mass error upper limit MassError_upBond and MassError_lowBond, and judge the numerical size of the current mass error with the mass error upper limit and mass error lower limit, which is divided into the following situations:

[0143] Case 1: If the current mass error is greater than the set upper limit of the mass error, reduce the time step, roll back the simulation time to the previous moment, and the state variable values stored in each control volume and the takeover data container are rolled back to the values at the previous moment, that is, rolled back to the values before being updated, and re - perform the calculation and solution of the two - fluid six - equation system;

[0144] Case 2: If the current mass error is less than the set lower limit of the mass error, increase the time step and continue the time advancement;

[0145] Case 3: If the current mass error is between the set upper and lower limits of the mass error, do not adjust the time step and continue the time advancement.

[0146] The beneficial effects of the present invention are as follows: (1) The present invention solves the problem of calculating and solving the two - fluid six - equation based on parallel multi - threads, establishes a parallel calculation process framework for the two - fluid six - equation, supports the parallel calculation of multi - basic node units for the generalized source term function of the two - fluid six - equation, the parallel calculation of multi - basic node units for the block - inverse of the coefficient matrix of the two - fluid six - equation, the parallel calculation of multi - loop for the pressure matrix equation, the parallel calculation of multi - basic node units for the takeover flow velocity based on the control - volume pressure, and the parallel calculation of multi - basic node units for solving scalars such as the void fraction of the control volume based on the takeover flow velocity. The calculation process is complete and effective, and the solving efficiency is greatly improved. The following table shows the comparison of the parallel calculation efficiency of the two - fluid six - equation after adopting the present invention:

[0147] Table 1 Comparison of Parallel Calculation Efficiency of Typical Thermal - Hydraulic System Models

[0148]

[0149] (2) The present invention establishes a framework for storing, reading, writing, and multi - node parallel calculation of data information of basic node units in the thermal - hydraulic system based on data containers. There is a one - to - one correspondence between the data container and the basic node unit. It can not only store the data information of the basic node units required for solving the two - velocity six - equation coefficient system, realize the fast retrieval and positioning of parameter variables based on the basic node unit number; but also create data reading and writing functions based on the basic node unit number to directly operate on the memory space corresponding to the parameter variables, avoiding a large number of intermediate variables and value - passing functions; and can also use the data container as a calculation carrier to realize the simultaneous calculation of the generalized source term function, matrix coefficient calculation, or parameter variable calculation for multiple basic node units. During the multi - node parallel calculation process, the data storage and reading and writing related to each basic node unit are only executed in its corresponding data container, avoiding the mutual interference of parameter variables of each basic node unit during the multi - node parallel calculation.

[0150] (3) The present invention provides a two-layer numbering method for the loop numbers and node numbers of the basic node units in the thermal-hydraulic system, which not only supports assembling the pressure matrix equations of each independent loop according to the loop numbers, that is, traversing the pressure matrix coefficients corresponding to the basic node units of the control volume type under the same loop number to assemble the pressure matrix equation of the independent loop corresponding to the loop number, and realizing the parallel solution of the global pressure equation by block according to the loop; but also supports carrying out the calculation of the pressure matrix coefficients of each control volume by multiple basic node units according to the node numbers, realizing the flattening process of the assembly calculation of the coefficient matrices of the pressure equations of different independent loops, ensuring the uniform distribution of the load tasks of each thread in the parallel calculation, and improving the efficiency of the parallel calculation. Description of the Drawings

[0151] Figure 1 It is the classification principle and execution order of the large category of the source term function for the two-fluid six-equation system;

[0152] Figure 2 It is an example of the thermal-hydraulic system model;

[0153] Figure 3 It is an example of the multi-thread parallel calculation process based on the large category of the source term function;

[0154] Figure 4 It is the topological structure of the thermal-hydraulic system;

[0155] Figure 5 It is an example of the topological structure of the thermal-hydraulic system. Detailed Embodiment

[0156] The present invention will be further described in detail below with reference to the drawings and specific embodiments.

[0157] A parallel calculation method applicable to the two-fluid six-equation system includes the following steps:

[0158] Step 1: Conduct a Courant number check. Based on the state variable values of each control volume and pipe connection in the thermal-hydraulic system at the current moment, calculate the Courant limit of the current two-fluid six-equation system. If the current time step is greater than the Courant limit, adjust the current time step to the Courant limit; if the current time step is less than or equal to the Courant limit, do not adjust the current time step.

[0159] The Courant number check includes the following:

[0160] The Courant number Co is a ratio that measures the flow distance of the fluid microelement to the local grid cell length under the simulation time step, and its definition is:

[0161]

[0162] Among them: V is the flow velocity, Δt is the time step, and Δx is the grid cell length.

[0163] When the Courant number Co is less than 1, within one time step, a fluid element cannot directly "cross" a grid cell, that is, the grid cell "captures" the fluid element within the time step, and the structural parameters, state parameters, source term functions, etc. of the grid cell can act on the fluid element. When the Courant number Co is greater than 1, within one time step, the fluid element can directly "cross" the grid cell, that is, the grid cell cannot "capture" the fluid element within the time step, and various parameters and source term functions of the grid cell cannot act on the fluid element, resulting in the lack of flow condition information of the fluid element during the simulation calculation process. Therefore, to ensure the stability and reliability of fluid simulation calculations, the Courant number Co of all grid cells should be less than 1.

[0164] To prevent the Courant number from being greater than 1 during the fluid simulation process, before each time advancement, based on the limitation that the Courant number is not greater than 1 and combined with the flow velocity V of the previous time step n the Courant time of each basic node unit is calculated. The Courant time is the time step corresponding to when the Courant number of the basic node unit is equal to 1. When the pre - used time step is greater than the minimum Courant time of each basic node unit in the thermal - hydraulic system, that is, the Courant limit, then the time step is changed to the Courant limit. The calculation method of the Courant limit is:

[0165]

[0166] where M is the total number of non - boundary nozzles, Δt upc，m and Δt downc，m are the upstream and downstream upwind Courant times of the m - th nozzle, V K，m is the volume of the upstream control volume of the m - th nozzle, A m is the flow area of the m - th nozzle, v g，m and v f，m are the vapor - phase and liquid - phase flow velocities of the m - th nozzle, Δx K，m and Δx L，m are the lengths of the upstream and downstream control volumes of the m - th nozzle; v f，L，m and v g，L，m are the liquid - phase and vapor - phase flow velocities at the center of the downstream control volume of the m - th nozzle.

[0167] Step 2: Based on the current time step and the current values of the control volume and takeover state variables, perform parallel calculations of the generalized source term functions for the basic node elements to obtain the complete data information required for assembling the two-fluid six-equation coefficient matrix, and store it in the corresponding data container of the basic node elements. The generalized source term functions refer to the physical source terms such as wall heat transfer calculation, interphase heat transfer calculation, interphase friction calculation, wall friction calculation, local resistance calculation, etc., the solution variable coefficients and constant terms of the mass, energy, and momentum conservation discrete equations, and the calculation preparation parts required for calculating these source terms, solution variable coefficients, constant terms, etc., including physical property calculation, takeover variable calculation based on control volume variables, control volume variable calculation based on takeover variables, etc.

[0168] The generalized source term functions in the parallel calculation of the generalized source term functions refer to the physical source terms such as wall heat transfer calculation, interphase heat transfer calculation, interphase friction calculation, wall friction calculation, local resistance calculation, etc., the solution variable coefficients and constant terms of the mass, energy, and momentum conservation discrete equations, and the calculation preparation parts required for calculating these source terms, solution variable coefficients, constant terms, etc., including physical property calculation, takeover variable calculation based on control volume variables, control volume variable calculation based on takeover variables, etc.

[0169] The process of parallel calculation of the generalized source term functions includes:

[0170] Step 21: Define the data container corresponding to the basic node element, and store the data information required for assembling and solving the two-fluid velocity six-equation coefficient matrix. The data information is stored and read using a certain data structure.

[0171] Step 22: Define the data reading Get(id, VarType) function and data writing Set(id, VarType, inputVal) function for the data container of the basic node element. The Get function and Set function read and write the data information at the corresponding position of the data container data structure according to the basic node element number id and the variable name VarType.

[0172] Step 23: Use the Get(id, VarType) function and Set(id, VarType, inputVal) function to convert each source term function into a source term function model with the basic node element number id as the only input variable. The source term function model reads the data stored in the data container of the basic node element numbered id as the input variable of the source term function, calculates the function output result using the mathematical principle formula of the source term function, and then uses the Set function to write the function output result into the memory space corresponding to the source term function output variable in the data container of the basic node element numbered id.

[0173] Step 24: Define the source term function enumeration table. The source term function enumeration table encompasses all generalized source term function types for the two-fluid six-equation system calculation. The position serial number in the enumeration table represents the execution order of the source term functions, expressing the mathematical logic of the source term function calculation in the two-fluid six-equation system. Further, to reduce the system traversal times and the communication times between the sub-threads and the main thread during the parallel calculation of the source term function model, a major category of source term functions is introduced. The major category of source term functions refers to a set of source term functions that are adjacent in execution order in the source term function enumeration table, belong to the same type of basic node unit, and are spatially independent of each other. For example, physical property source term functions such as density, viscosity coefficient, thermal conductivity, and specific heat at constant pressure belong to the same major category.

[0174] Step 25: According to the classification principle and execution order of the major category of source term functions, traverse the basic node units of the thermohydraulic system under the same major category of source term functions, and carry out the source term function model calculation of each basic node unit by means of parallel calculation of multiple basic node units.

[0175] Step 251: Storage and reading / writing of basic node unit data information based on the data container

[0176] The data container refers to the memory space allocated according to a certain data structure. The storage data structure of the parameters, variables, and source term functions of the basic node is defined using a variable enumeration table. The variable enumeration table refers to an array in which all input and output variables and parameters of the generalized source term functions of the basic node unit are stored according to a certain data structure. The variable enumeration table contains all input and output data required for the calculation of the generalized source term functions of the basic node unit.

[0177] After defining the variable enumeration table, define the reading function Get(id, VarType) and the writing function Set(id, VarType, inputVal). The former reads the variable VarType stored in the data container of the basic node unit numbered id, and the latter writes the externally obtained parameter inputVal to the memory location corresponding to the variable VarType in the data container of the basic node unit numbered id. The variable VarType must be defined in the variable enumeration table. Through the Get(id, VarType) function and the writing Set(id, VarType, inputVal) function, direct operations on the memory space corresponding to the parameter variables are realized, avoiding a large number of equivalent or value-passing functions.

[0178] Step 252: Definition of the source term function model with the basic node unit number id as the only input variable

[0179] According to the mathematical principle formula of the source term function, the Get(id, VarType) function and the Set(id, VarType, inputVal) function are used to convert each source term function into a source term function model with the basic node unit number id as the only input variable. The specific method is as follows:

[0180] Step 2521: Define the input variables, output variables, and process variables of the source term function according to the source term function name.

[0181] Step 2522: Use Get(id, VarType) to read the data required for the input variables from the data container of the basic node unit with number id.

[0182] Step 2523: According to the mathematical principle formula of the source term function, obtain the process variable data and output variable data of the function calculation through the input variable data.

[0183] Step 2524: Write the obtained output variable data into the memory location corresponding to the output variable of the source term function in the data container of the basic node unit with number id using Set(id, VarType, inputVal).

[0184] Step 253: Multithreaded parallel calculation based on the large category of source term functions

[0185] Define a source term function enumeration table, which is an array that stores all the generalized source term function types of the two-fluid six-equation system according to a certain data structure. The position sequence number in the enumeration table represents the execution order of the source term function.

[0186] The execution order of the source term functions is defined in the source term function enumeration table, which expresses the mathematical logic of the source term function calculation in the two-fluid six-equation system. Under the same source term function, the basic node units of the thermohydraulic system can be traversed, and the calculation of the source term function model for each basic node unit can be carried out in a parallel calculation mode with multiple basic node units. However, the calculation of each source term function model requires traversing the basic node units of the thermohydraulic system. To reduce the number of traversals of the basic node units of the thermohydraulic system, a parallel calculation method based on the major categories of source term functions is adopted for the traversal calculation of the basic node units of the thermohydraulic system, that is, under the same major category of source term functions, the number of traversals of the basic node units of the thermohydraulic system is 1, rather than 1 traversal for each individual source term function. The major categories of source term functions are those source term functions that are adjacent in execution order in the source term function enumeration table, belong to the same type of basic node units, and are not dependent on each other spatially. For example, physical property source term functions such as density, viscosity coefficient, thermal conductivity, and specific heat at constant pressure belong to the same major category. At the same time, adopting the parallel calculation method based on the major categories of source term functions avoids the situation where the sub-thread must return to the main thread to receive new source term function calculation instructions after completing the calculation of a single source term function during the parallel calculation process, effectively reducing the number of data communications between the sub-thread and the main thread during the parallel calculation process.

[0187] The classification principle and execution order of the major categories of source term functions in the two-fluid six-equation system are as Figure 1 shown, and successively include the major category of source term functions for control volume fluid property calculation, the major category of source term functions for header average quantity calculation, the major category of source term functions for control volume average quantity calculation, the major category of source term functions for wall heat transfer, the major category of source term functions for two-phase flow pattern judgment, the major category of source term functions for inter-phase heat transfer, the major category of source term functions for inter-phase and wall friction, the major category of source term functions for HLOSS item calculation, the major category of source term functions for dissipation item calculation, the major category of source term functions for display velocity calculation, the major category of source term functions for critical flow, and the major category of source term functions for CCFL.

[0188] Taking Figure 2 the thermohydraulic system model as an example, the thermohydraulic system model includes 5 non-boundary control volumes and 6 non-boundary headers. The multi-thread parallel calculation process based on the major categories of source term functions is as Figure 3 shown.

[0189] Step 3: Read the data information stored in the basic node unit data container, and assemble the coefficient matrices of the mass and energy conservation equations and the coefficient matrix of the momentum conservation equation according to the thermohydraulic system topology.

[0190] The assembly process of the coefficient matrices of the mass and energy conservation equations and the coefficient matrix of the momentum conservation equation includes that in the thermal-hydraulic system model, assuming there are N non-boundary control volumes and M non-boundary nozzles, then there are 5N control volume variables (pressure, void fraction, specific internal energy of the vapor phase, specific internal energy of the liquid phase, non-condensable gas fraction) and 2M nozzle variables (vapor phase flow rate, liquid phase flow rate). The discrete forms of the scalar equations such as the mass conservation and equation, mass conservation difference equation, vapor phase energy conservation equation, liquid phase energy conservation equation, non-condensable gas transport equation, etc. for the N control volumes are sorted out in the following form,

[0191]

[0192] where P is the fluid pressure, a is the void fraction, U is the specific internal energy, X is the non-condensable gas fraction, v is the flow rate, a 1 to a 5 are respectively the variable coefficients in the scalar equation, s 1 to s 4 are respectively the variable coefficients in the scalar equation, b is the constant term of the equation. The subscript g represents the vapor phase, the subscript f represents the liquid phase, the subscript L represents the current control volume, the subscript j represents the nozzle number, j + 1 represents the downstream nozzle, the superscript n represents the nth time moment, and the superscript n + 1 represents the (n + 1)th moment. a 1 to a 5 , s 1 to s 4 , and b are all obtained by calculating the source term function and stored in the data containers of the corresponding control volumes. The scalar equations for the N control volumes are combined into matrix form:

[0193] AX = SV + B (6)

[0194] where,

[0195]

[0196] A is the 5N×5N coefficient matrix for solving the scalars such as the control volume pressure, void fraction, specific internal energy of the vapor phase, specific internal energy of the liquid phase, non-condensable gas fraction, etc., X is the 5N control volume variable array, S is the 5N×2M coefficient matrix for solving the vectors such as the vapor phase flow rate and liquid phase flow rate of the control volume upstream and downstream nozzles, V is the 2M nozzle variable array, and B is the 5N constant term array.

[0197] According to the topological structure of the thermal-hydraulic system, the momentum conservation and equation and momentum conservation difference equation for the M nozzles are sorted out in the following form,

[0198]

[0199] where, g1 , g 2 , w 1 , w 2 are respectively P L n+1 , P K n+1 the variable coefficients in the equation, and f is the constant term of the equation. g 1 , g 2 , w 1 , w 2 , f are all obtained by calculating the source term function and stored in the corresponding data containers of the takeover.

[0200] The momentum conservation equations for M takeovers are combined into matrix form:

[0201] GV = WP + F (8)

[0202] where

[0203] G is a 2M×2M matrix, V is an array of vapor and liquid flow velocity variables for 2M takeovers, W is a 2M×N matrix, P is an array of pressure variables for N control volumes, and F is an array of 2M constant terms. The coefficient matrices A, S, G, and W are sparse matrices. The row elements of the matrix correspond to the coefficients of all solution variables in a certain conservation equation. When a certain variable is not in this equation, its coefficient is 0. According to the discrete format of the two-fluid six-equation, the scalar conservation discrete equations of mass, energy, non-condensable gas, etc. of the control volume only contain the variables and parameters of this control volume and the upstream and downstream takeovers, while the momentum conservation discrete equation of the takeover only contains the variables and parameters of this takeover and the upstream and downstream control volumes. The coefficients of the variables of the nodes not adjacent to this control volume or takeover in the matrix row elements must be 0, so there will be a large number of 0 elements in the matrix.

[0204] For example, for the following thermal-hydraulic system model, which includes two non-boundary control volumes V1 and V2 and 2 non-boundary takeovers. According to the topology of the thermal-hydraulic system, the coefficient matrix of the mass, energy, and non-condensable gas conservation equations AX = SV + B is:

[0205]

[0206] The coefficient matrix of the momentum conservation equation GV = WP + F is:

[0207]

[0208] Step 4: Perform a block operation on the coefficient matrices of the mass and energy conservation equations and the momentum conservation equation of the two-fluid six-equation system, convert the coefficient matrix with the characteristics of a large sparse matrix into a combination of several small matrices, and further perform an inversion operation on the coefficient matrix. Establish a block matrix calculation formula for the scalar array of the mass and energy conservation equations based on the control volume and a block matrix calculation formula for the vector array of the momentum conservation equation based on the takeover, significantly reducing the time complexity of the coefficient matrix.

[0209] Among them, the block operation and inversion operation of the two-fluid six-equation coefficient matrix

[0210] By adjusting the variable array in equation (9), make the variables of the same basic node unit adjacent in the variable array, and perform a non-zero element concentration process on the coefficient matrix of the scalar equation, that is:

[0211]

[0212] Equation (11) can be simplified to:

[0213]

[0214] Where: A 1 and A 2 are 5×5 matrices, X 1 and X 2 are 5×1 scalar arrays of the basic node, S 1 and S 2 are 5×4 matrices, V 1 and V 2 are 2×1 vector arrays of the basic node, B 1 and B 2 are 5×1 constant arrays.

[0215] Perform an inversion operation on equation (12), then:

[0216]

[0217] Without loss of generality, under N non-boundary control volumes and M non-boundary takeovers, equation (13) can be written as:

[0218]

[0219] Where, X 1 to X N are 5×1 scalar arrays of the basic node, A 1 to A N are 5×5 matrices, S 1 to S N are 5×2M matrices, V 1 to V Mis a 2×1 vector array of basic nodes, B 1 to B N is a 5×1 constant array.

[0220] It can be seen from equation (14) that after the scalar equation AX = SV + B is block-processed, its time complexity changes from O(125N 3 ) to O(125N), and the theoretical calculation cost increases linearly with the number of control volumes.

[0221] Since the elements of the A i (i = 1, …, N) matrix and B i (i = 1, …, N) can all be calculated from the parameter variables stored in the data container corresponding to the current control volume, while only the columns corresponding to the pipes connected to the current control volume in S i (i = 1, …, N) have non-zero elements. Therefore the calculation of only relates to the current control volume and its connected pipes and meets the conditions for parallel calculation. At the same time, since A i is a 5×5 matrix, its inverse matrix is easily obtained, and the positions of the non-zero element columns can be quickly confirmed through the topological structure information of the current control volume and the pipes, further improving the calculation efficiency.

[0222] Similarly, by adjusting the variable array in equation (10) to make the variables of the same basic node unit adjacent in the variable array, the non-zero elements of the coefficient matrix of the scalar equation are centralized:

[0223]

[0224] Equation (15) can be simplified to:

[0225]

[0226] where: C 1 and C 2 are 2×2 matrices, X 1 and X 2 are 2×1 vector arrays of basic node units, W 1 and W 2 are 2×2 matrices, P 1 and P 2 are the pressures of the basic node units, and F 1 and F 2 are 2×1 constant arrays.

[0227] Performing an inverse operation on equation (16), then:

[0228] ​

[0229] Without loss of generality, under N non-boundary control volumes and M non-boundary takeovers, Equation (17) can be written as:

[0230]

[0231] where V 1 to V M is a 2×1 vector array of basic nodes, G 1 to G M is a 2×2 matrix, W 1 to W M is a 2×N matrix, and F 1 to F N is a 2×1 constant array.

[0232] Similarly, by block processing the vector matrix GV = WP + F, its time complexity changes from O(4M 2 ) to O(4M), and the computational cost increases linearly with the number of takeovers. and are calculated only related to the current takeover and its upstream and downstream control volumes, meeting the conditions for parallel computing. Since G j is a 2×2 matrix, its inverse matrix is easily obtained, and the positions of the non-zero element columns of can be quickly confirmed through the topological structure information of the current takeover and control volume, further improving the computational efficiency.

[0233] Step 5: Combine the constraint equations of the control volume pressure and the vapor-liquid flow rate of the takeover in the scalar array block matrix calculation formula and the constraint equations of the vapor-liquid flow rate of the takeover and the control volume pressure in the vector array block matrix calculation formula to form a pressure matrix equation with only the control volume pressure as the solution variable. Further, during the assembly process of the pressure matrix equation, the thermal-hydraulic loop and the basic node unit are numbered in two layers using the breadth-first traversal algorithm, and the pressure matrix equations of each independent loop are established according to the loop numbers.

[0234] Among them, the generation of the pressure matrix equation includes, after block processing and inversion of the two-fluid six-equation coefficient matrix, extracting all the matrix element expressions of the control volume pressure in Equation (14), that is, the first rows of X i , and , and determining the distribution of the non-zero elements in the first row of through the connection topological information between the control volume and the takeover. Then, the matrix element expression of Equation (14) regarding the pressure can be written in the following form:

[0235]

[0236] where j is the inlet nozzle number, m is the total number of inlet nozzles, jp is the outlet nozzle number, n is the total number of outlet nozzles, and b i is the constant term with respect to pressure P i , C P,g,j , C P,f,j , C P,g,jp , C P,f,jp is the non-zero elements of the first row.

[0237] Furthermore, extract the matrix element expressions for the flow velocities of the inlet and outlet nozzles of the i-th control volume in equation (18). Taking the j-th inlet and jp-th outlet nozzles as an example, the following form can be obtained:

[0238] V g,j = f V,g,j + C V,g,K,j P K,j + C V,g,i,j P i (20)

[0239] V f，j = f V，f，j + C V，f，K，j P K，j + C V，f，i，j P i (21)

[0240] V g，jp = f V，g，jp + C V，g，i，jp P i + C V，g，L，jp P L，jp (22)

[0241] V f，jp = f V,f,jp + C V，f，jp P i + C V，f，L，jp P L，jp (23)

[0242] where f V，g，j and f V,f，j are respectively the constant terms of the first row and the second row, C V，g，K，j and C V，f，K，j are respectively the matrix elements of the first row and the second row with respect to the pressure of the upstream control volume of the j-th nozzle, C V，g，i，j and C V，f，i，j are respectively the matrix elements of the first row and the second row with respect to the pressure of the i-th control volume, f V，g，jp and f V,f,jp are respectively The constant terms in the first and second rows of, C V，g，i，jp and C V，f，i，jp are respectively The matrix elements in the first and second rows of with respect to the pressure of the i-th control volume, C V，g，L，jp and C V，f，L，jp are respectively The matrix elements in the first and second rows of with respect to the pressure of the control volume downstream of the jp-th pipe connection.

[0243] Substituting equations (20) to (23) into equation (19), we can obtain:

[0244]

[0245] The coefficients and constants of each term in equation (24) can all be obtained through and other matrices. Further, by traversing according to the control volume numbers, equation (24) is converted into a matrix equation,

[0246] A p P = b p (25)

[0247] A p is a large sparse matrix. The non-zero elements in the i-th (i = 1,..., N) row of the A p matrix only include the i-th column and the column numbers corresponding to the pressures of the non-boundary control volumes that are connected by pipe connections to the i-th control volume. The number of its non-zero elements is m + n + 1 (m is the total number of inlet pipe connections, n is the total number of outlet pipe connections), which is much smaller than the number of control volumes N. The pressures P of each control volume in the thermal-hydraulic system can be solved through equation (25) n+1 , and then V can be calculated through equation (18) n+1 , and finally other scalars except pressure can be calculated through equation (14)

[0248] In the process design of an actual thermal-hydraulic system, there may be several loops that are not fluidly connected to each other. These loops cannot affect each other's pressure fields through fluid flow. It is manifested that the non-zero elements in the row corresponding to the control volume of a certain loop in the pressure equation are only related to that loop, and the coefficients of other loop control volumes in that row are 0. Then, through a certain linear transformation, the pressure equation can be converted into several block matrices, and each block matrix corresponds to each independent loop.

[0249] Taking Figure 5 the topological structure of the thermal-hydraulic system as an example, this thermal-hydraulic system contains 2 loops that are not fluidly connected to each other, and heat interaction is carried out between the loops through heat components.

[0250] A pressure equation is established for this thermal-hydraulic system, and its structure is as shown in (26).

[0251]

[0252] Perform a linear transformation on equation (26), that is, adjust the order of the variables in the variable array, and group the control volume variables of the same loop together, as shown in equation (27).

[0253]

[0254] Equation (27) can be simplified to:

[0255]

[0256] Where: and are the block matrices of the pressure equations corresponding to Loop 1 and Loop 2 respectively; and are the arrays of pressure variables for Loop 1 and Loop 2; and are the arrays of constant terms for Loop 1 and Loop 2.

[0257] Without loss of generality, considering k independent loops, we get:

[0258]

[0259] That is:

[0260]

[0261] This shows that each independent loop has an independent pressure equation and can be solved separately, meeting the conditions for parallel solution. To quickly assemble the coefficient matrix of the pressure equation in equation (30), eliminate the linear transformation from equation (27) to equation (28), and avoid the problem of uneven parallel computing resource allocation due to different numbers of control volumes in each independent loop during the block assembly of the coefficient matrix of the pressure equation based on loops, the following strategy is adopted:

[0262] Step 51: Use the breadth-first search algorithm to number the basic node units of the thermal-hydraulic system, generating a two-layer numbering including loop numbers and node numbers, including the following sub-steps:

[0263] Step 511: Change the directed edges of the flattened topological structure diagram of the thermal-hydraulic system to undirected edges. The flattened topological structure diagram is a system model that restores the thermal-hydraulic system model constructed by the user to a system model completely composed of basic node unit models and connected in an orderly manner according to control volume - takeover - control volume (V-J-V).

[0264] Step 512: Establish a data container A for the graph vertices that contains the undirected graph of the flattened topology of the thermal-hydraulic system. The data container A stores the basic node unit names of all graph vertices and the names of the model components where the basic node units are located in a queue manner.

[0265] Step 513: Obtain the first graph vertex stored in the data container. The loop number of this graph vertex is marked as 1 and is stored in a temporary data container B. The temporary data container B expresses the order of the graph vertex elements stored in a queue form.

[0266] Step 514: Obtain the first graph vertex element in the temporary data container B and traverse all the connected graph vertex elements of this graph vertex element, which is divided into the following situations:

[0267] Situation 1: The connected graph vertex element has been numbered, continue traversing until the traversal ends;

[0268] Situation 2: The connected graph vertex element has not been numbered, store it in the temporary data container B in a queue, and continue traversing until the traversal ends.

[0269] Step 515: After the traversal of the connected graph vertex elements ends, delete the first graph vertex element from the temporary data container B and store it in the result data container C;

[0270] Step 516: Record the type of this graph vertex element and the sequence number of the graph vertex element stored in the result data container C under different types. The types of the graph vertex elements include control volumes and nozzles. Take the sequence number of the control volume type graph vertex element stored in the result data container C as the node number of the control volume, and take the sequence number of the nozzle type graph vertex element stored in the result data container C as the node number of the nozzle.

[0271] Step 517: Judge whether this graph vertex element has been marked with a loop number, which is divided into the following situations:

[0272] Situation 1: This graph vertex element has not been marked with a loop number. The loop number of this graph vertex is marked as the loop number of the graph vertex with the previous sequence number in the result data container C, and then continue with the next sub-step operation;

[0273] Situation 2: This graph vertex element has been marked with a loop number, continue with the next sub-step operation.

[0274] Step 518: Judge whether the temporary data container B is empty, which is divided into the following situations:

[0275] Situation 1: The temporary data container B is not empty, jump to Step 514;

[0276] Situation 2: The temporary data container B is empty, continue with the next sub-step operation.

[0277] Step 519: traverse the graph vertex elements of data container A to determine whether there are graph vertex elements in data container A that are not node-numbered, which is divided into the following two cases:

[0278] Case 1: There are graph vertex elements that are not node-numbered. Randomly select a graph vertex element that is not node-numbered and store it in the temporary data container B. The loop number of the graph vertex element is marked as the loop number of the graph vertex element that is currently stored in the result data container C last + 1, and jump to step 514;

[0279] Case 2: If there are no graph vertex elements that are not node-numbered, the basic node element numbering ends.

[0280] Step 52: Without distinguishing the loop numbers of the basic node units, the node numbers of the control body are traversed in parallel computing to calculate the coefficients required for assembling the pressure equation coefficient matrix, as shown in formula (100). Since the coefficient calculation required for assembling the pressure equation coefficient matrix is ​​based on the node number of the control body as the basis for parallel multi-threaded task allocation instead of the loop number, the flattening of the pressure equation coefficient matrix assembly calculation of different independent loops is achieved, ensuring the uniform distribution of the load and tasks of each thread of parallel computing, and improving the efficiency of parallel computing.

[0281] Step 53: According to the loop number of the basic node unit, traverse the node number of the basic node unit under the loop number in turn, obtain the pressure matrix coefficient corresponding to the basic node unit of the control body class, and assemble the block coefficient matrix of the pressure equation of the loop number, see formula (30), to create conditions for parallel calculation of multi-loop pressure equations.

[0282] Step 6: Based on the direct solution method, the pressure matrix equations of each independent loop are solved and calculated in parallel to obtain the pressure of each control body of the thermal hydraulic system, and the pressure is stored and updated in the data container of each control body, as follows:

[0283] Assuming that the independent loop has N control bodies, the non-zero elements in the i-th (i=1, ..., N) row of the pressure matrix only include the i-th column and the number of columns corresponding to the pressure of non-boundary control bodies connected to the i-th control body by a connecting pipe. The number of non-zero elements is much smaller than the number of control bodies N, so the pressure matrix is ​​a large sparse matrix.

[0284] Since the number of connected control bodies of each control body in the thermal-hydraulic system is different, and the values ​​of the non-zero elements of the pressure equation coefficient matrix vary greatly during the transient operation, the diagonal dominance, positive definiteness and symmetry of the matrix cannot be guaranteed, and the influence of the matrix iteration residual convergence condition on the calculation accuracy cannot be effectively evaluated. Therefore, this patent of the present invention adopts a direct solution method based on the sparse LU decomposition method as the solution algorithm for the pressure matrix equation.

[0285] The sparse LU decomposition method aims to obtain an elementary change matrix Q, a lower triangular matrix L (all the main diagonal elements are 1) and an upper triangular matrix U, such that:

[0286] QA=LU (31)

[0287] Among them, the elementary change matrix Q extracts the absolute maximum element of each column as the elimination principal element through basic row transformation during the LU decomposition process, reducing the condition number of the matrix to ensure the numerical stability of the calculation. Then, the block matrix (30) of the pressure equation is converted to:

[0288]

[0289] because is an elementary transformation matrix, then It is an array The elements are reordered and then the lower triangular matrix and the upper triangular matrix Structural features to obtain independent circuit pressure distribution

[0290] Step 7: Call the vector array block matrix calculation formula of the momentum conservation equation based on the pipe, take the control body pressure as input, calculate the vapor-liquid flow rate of each pipe of the thermal hydraulic system, and store and update it in the data container of each pipe. The vapor-liquid flow rate of the pipe is calculated as follows:

[0291] Call the momentum conservation equation vector array block matrix calculation formula (18) based on the takeover to control the body pressure P n+1 As input, calculate the flow rate of each pipe in the thermal hydraulic system The updates are stored in the data containers of each takeover.

[0292]

[0293] Where: the superscript n+1 represents the calculated value at the current moment, n represents the calculated value at the previous moment, and M is the number of non-boundary pipes in the thermal-hydraulic system.

[0294] From formula (33), it can be seen that the flow rate of each pipe is Only for takeover W j and F j It is related to the current flow rate of other pipes, so (33) can be calculated in parallel using multiple threads. j The non-zero elements of are only related to the upstream and downstream control bodies of the jth pipe, so the topological structure information of the thermal hydraulic system can be used to quickly locate W j The non-zero elements of provide efficiency in matrix calculations.

[0295] Step 8: Call the scalar array block matrix calculation formula of the mass and energy conservation equations based on the control volume, take the vapor-liquid flow velocities of each pipe in the thermal-hydraulic system as inputs, calculate the scalar quantities such as the void fraction, vapor specific internal energy, liquid specific internal energy, and non-condensable gas fraction of each control volume in the thermal-hydraulic system, and store and update them in the data containers of each control volume.

[0296] The calculation of the void fraction, vapor-liquid specific internal energy, and non-condensable gas fraction of the control volume is as follows:

[0297] Call the scalar array block matrix calculation formula (Equation (14)) of the mass and energy conservation equations based on the control volume, take the vapor-liquid flow velocity V of the pipes in the thermal-hydraulic system as the input, and calculate the void fraction of the control volume in the thermal-hydraulic system. n+1 For input, calculate the void fraction of the control volume in the thermal-hydraulic system Vapor specific internal energy Liquid specific internal energy Non-condensable gas fraction And other scalar quantities, and store and update them in the data containers of each control volume.

[0298]

[0299] Among them, since the pressure P of each control volume n+1 has been obtained, the first row calculation of the scalar array block matrix for each control volume is not performed in Equation (34), and only the second row and the fifth row are calculated. It can be seen from Equation (34) that the of each control volume only depends on the S i and B i related to it. Therefore, Equation (34) can adopt a multi-threaded parallel calculation method. At the same time, since the non-zero elements of S i only relate to the upstream and downstream pipes of the i-th control volume, the non-zero elements of S i can be quickly located through the topological structure information of the thermal-hydraulic system, improving the efficiency of matrix calculation.

[0300] Step 9: Calculate the mass error of the solution of the two-fluid six-equation system. The mass error is the deviation between the vapor-liquid mixture density obtained through the first-order Taylor expansion and the vapor-liquid mixture density calculated through the physical property state equation, and judge the numerical magnitudes of the current mass error and the set mass error upper limit and mass error lower limit. It is divided into the following situations:

[0301] Situation 1: If the current mass error is greater than the set mass error upper limit, the time step becomes 1 / 2 of the original time step, the simulation time rolls back to the previous moment, and the values of the state variables stored in the data containers of each control volume and pipe roll back to the values at the previous moment, that is, roll back to the values before being updated, and jump to Step 2.

[0302] Case 2: If the current mass error is less than the set lower limit of the mass error, change the time step to 1.1 times the original time step and continue the time advancement.

[0303] Case 3: If the current mass error is between the set upper and lower limits of the mass error, do not adjust the time step and continue the time advancement.

[0304] The mass error in the solution of the two-fluid six-equation system is the deviation between the vapor-liquid mixture density obtained through the first-order Taylor expansion and the vapor-liquid mixture density calculated through the physical property state equation, and its expression is:

[0305] ε m = max(ε mlocal , ε mtotal ) (35)

[0306]

[0307] ρ i = α g，i ρ g，i + (1 - α g，i )ρ f，i (39)

[0308]

[0309] Where and are the vapor phase density and liquid phase density obtained through the first-order Taylor expansion for the i-th control volume; ρ g，i and ρ f，i are the vapor phase density and liquid phase density obtained through the physical property state equation or interpolation table for the i-th control volume; k is a factor greater than 1 (recommended to be 1.41); N is the number of non-boundary control volumes.

[0310] Set the upper limit of the mass error MassError_upBond and the lower limit of the mass error MassError_lowBond, and judge the numerical magnitudes of the current mass error, the upper limit of the mass error, and the lower limit of the mass error. It is divided into the following cases:

[0311] Case 1: If the current mass error is greater than the set upper limit of the mass error, reduce the time step, roll back the simulation time to the previous moment, and the values of the state variables stored in each control volume and the takeover data container are rolled back to the values at the previous moment, that is, rolled back to the values before being updated, and re-perform the calculation and solution of the two-fluid six-equation system.

[0312] Case 2: If the current mass error is less than the set lower limit of the mass error, increase the time step and continue the time advancement.

[0313] Case 3: If the current mass error is between the set upper and lower limits of the mass error, the time step is not adjusted and the time advancement continues.

[0314] Step 10: The simulation time advances along the time stamp with a step size of one time step. Repeat steps 1 to 9 until the current simulation time is greater than or equal to the set simulation time.

Claims

1. A parallel computing method for a two-fluid six-equation system, characterized in that: The steps include: Step 1: Carry out Courant number check. Based on the state variable values ​​of each control body and nozzle of the thermal hydraulic system at the current moment, calculate the Courant limit of the current two-fluid six-equation system. If the current time step is greater than the Courant limit, adjust the current time step to the Courant limit. If the current time step is less than or equal to the Courant limit, do not adjust the current time step. Step 2: Based on the current time step and the current control volume and takeover state variable values, carry out parallel calculation of the generalized source term function of the basic node unit, obtain the complete data information required for assembling the coefficient matrix of the two-fluid six equations, and store it in the data container corresponding to the basic node unit; Step 3: Read the data information stored in the basic node unit data container, and assemble the coefficient matrix of the mass and energy conservation equations and the coefficient matrix of the momentum conservation equation according to the topological structure of the thermal hydraulic system; Step 4: Perform block operations on the coefficient matrices of the mass and energy conservation equations and the momentum conservation equation of the two-fluid six-equation system, convert the coefficient matrices with large sparse matrix characteristics into a combination of several small matrices, and further perform inversion operations on the coefficient matrices to establish the scalar array block matrix calculation formula for the mass and energy conservation equations based on the control body and the vector array block matrix calculation formula for the momentum conservation equation based on the takeover; Step 5: Combine the constraint equations of the control body pressure and the pipe vapor-liquid flow rate in the scalar array block matrix calculation formula and the constraint equations of the pipe vapor-liquid flow rate and the control body pressure in the vector array block matrix calculation formula to form a pressure matrix equation with only the control body pressure as the solution variable. Furthermore, in the process of assembling the pressure matrix equation, a breadth-first traversal algorithm is used to perform two-level numbering on the thermal hydraulic circuits and basic node units, and the pressure matrix equations of each independent circuit are established according to the circuit number; Step 6: Based on the direct solution method, the pressure matrix equations of each independent loop are solved and calculated in parallel to obtain the pressure of each control body of the thermal hydraulic system, and the pressure is stored and updated in the data container of each control body; Step 7: Call the vector array block matrix calculation formula of the momentum conservation equation based on the nozzle, take the control body pressure as input, calculate the steam and liquid flow rate of each nozzle of the thermal hydraulic system, and store and update it in the data container of each nozzle; Step 8: Call the block matrix calculation formula based on the mass and energy conservation equations of the control body, take the steam-liquid flow rate of each pipe of the thermal-hydraulic system as input, calculate the scalars such as the cavitation fraction, steam phase internal energy, liquid phase internal energy, and non-condensable gas fraction of each control body of the thermal-hydraulic system, and store and update them in the data container of each control body; Step 9: Calculate the mass error of the solution of the two-fluid six-equation system, where the mass error is the deviation between the vapor-liquid mixture density obtained by the first-order Taylor expansion and the vapor-liquid mixture density calculated by the physical property state equation, and determine the numerical value of the current mass error and the set mass error upper limit and mass error lower limit; Step 10: The simulation time is advanced along the timestamp with a time step length. Steps 1 to 9 are repeated until the current simulation time is greater than or equal to the set simulation time.

2. A parallel computing method for a two-fluid six-equation system as claimed in claim 1, characterized in that: The Courant number check described in step 1 includes the following: The Courant number Co is the ratio of the flow distance of the fluid microelement to the length of the local grid unit under the simulation time step, which is defined as: Where: V is the flow velocity, Δt is the time step, and Δx is the grid unit length; When the Courant number Co is less than 1, the fluid element cannot directly "cross" the grid unit in one time step. When the Courant number Co is greater than 1, the fluid element can directly "cross" the grid unit in one time step. To ensure the stability and reliability of fluid simulation calculations, the Courant number Co of all grid units should be less than 1. To prevent the Courant number from being greater than 1 during fluid simulation, before each time advancement, the Courant number is not greater than 1 and combined with the flow velocity V of the previous time step. n Calculate the Courant time of each basic node unit. The Courant time is the time step corresponding to the Courant number of the basic node unit equal to 1. When the pre-used time step is greater than the minimum value of the Courant time of each basic node unit of the thermal hydraulic system, that is, the Courant limit, then the time step is changed to the Courant limit. The calculation method of the Courant limit is: Where M is the total number of non-boundary takeovers, Δt upc，m and Δt downc，m V is the upstream upwind Courant time and downstream upwind Courant time of the mth connection, K，m is the volume of the upstream control volume of the mth takeover, A m is the flow channel area of ​​the mth connection, v g，m and v f，m is the vapor and liquid flow rate of the mth connection, Δx K，m and Δx L，m is the length of the upstream and downstream control bodies of the mth takeover; v f，L，m and v g，L，m are the liquid and vapor flow rates at the center of the downstream control body of the mth pipe.

3. A parallel computing method applicable to a two-fluid six-equation system as claimed in claim 1, characterized in that: The generalized source term function in the parallel calculation of the generalized source term function described in step 2 refers to physical source terms such as wall heat transfer calculation, phase heat transfer calculation, phase friction calculation, wall friction calculation, local resistance calculation, etc., and the variable coefficients and constant terms of the discrete equations of conservation of mass, energy, and momentum, as well as the calculation preparation parts required for calculating these source terms, solving variable coefficients, constant terms, etc., including physical property calculation, takeover variable calculation based on control body variables, and control body variable calculation based on takeover variables.

4. A parallel computing method applicable to a two-fluid six-equation system as claimed in claim 3, characterized in that: The parallel calculation process of generalized source term function includes: Step 21: define the data container corresponding to the basic node unit, store the data information required for assembling and solving the coefficient matrix of the two flow velocity six equations, and use a certain data structure to store and read and write the data information; Step 22: define a data reading Get(id, VarType) function and a data writing Set(id, VarType, inputVal) function for a basic node unit data container, wherein the Get function and the Set function read and write data information at corresponding positions of the data container data structure according to the basic node unit number id and the variable name VarType; Step 23: Use Get(id, VarType) function and Set(id, VarType, inputVal) function to convert each source item function into a source item function model with the basic node unit number id as the only input variable, the source item function model reads the data stored in the basic node unit data container numbered id as the input variable of the source item function through the Get function, calculates the function output result using the mathematical principle formula of the source item function, and then uses the Set function to write the function output result into the memory space corresponding to the source item function output variable in the basic node unit data container numbered id; Step 24: define a source term function enumeration table, which includes all generalized source term function types for calculating the two-fluid six-equation system. The position number in the enumeration table represents the execution order of the source term function, and expresses the mathematical logic of the source term function calculation of the two-fluid six-equation system. In order to reduce the number of system traversals and the number of communications between the sub-thread and the main thread in the parallel calculation process of the source term function model, a source term function category is introduced. The source term function category refers to a set of source term functions that are adjacent in execution order in the source term function enumeration table and belong to the same type of basic node units and have no spatial dependence on each other; Step 25: According to the classification principle and execution order of the source term function category, under the same source term function category, traverse the basic node units of the thermal-hydraulic system, and use multi-basic node unit parallel calculation method to carry out the source term function model calculation of each basic node unit.

5. A parallel computing method applicable to a two-fluid six-equation system as claimed in claim 4, characterized in that: The step 25 comprises: Step 251: Data container-based basic node unit data information storage and reading and writing The data container refers to a memory space allocated according to a certain data structure, and a variable enumeration table is used to define the storage data structure of the parameters, variables, and source term functions of the basic node. The variable enumeration table refers to an array in which the input and output variables and parameters of all generalized source term functions of the basic node unit are stored according to a certain data structure. The variable enumeration table contains all the input and output data required for the calculation of the generalized source term function of the basic node unit; After defining the variable enumeration table, define the Get(id, VarType) function and the Set(id, VarType, inputVal) function. The former reads the variable VarType stored in the basic node unit data container numbered id, and the latter writes the externally obtained parameter inputVal to the memory location corresponding to the basic node unit data container variable VarType numbered id. The variable VarType must be defined in the variable enumeration table. Through the Get(id, VarType) function and the Set(id, VarType, inputVal) function, direct operation of the memory space corresponding to the parameter variable is realized; Step 252: Define the source term function model with the basic node unit number id as the only input variable According to the mathematical principle formula of the source term function, the Get(id, VarType) function and the Set(id, VarType, inputVal) function are used to convert each source term function into a source term function model with the basic node unit number id as the only input variable; Step 253: Multi-threaded parallel computing based on the source function category A source term function enumeration table is defined, wherein the source term function enumeration table is an array that stores all generalized source term function types of a two-fluid six-equation system according to a certain data structure. The position numbers in the enumeration table represent the execution order of the source term functions. The execution order of the source term functions is defined in the source term function enumeration table, and the mathematical logic of the source term function calculation of the two-fluid six-equation system is expressed. Under the same source term function, the basic node units of the thermal-hydraulic system are traversed, and the source term function model calculation of each basic node unit is carried out using a multi-basic node unit parallel calculation method.

6. A parallel computing method for a two-fluid six-equation system as claimed in claim 5, characterized in that: The step 252 includes: Step 2521: According to the source function name, define the input variables, output variables and process variables of the source function; Step 2522: Use Get(id, VarType) to read the data required for the input variable from the basic node unit data container numbered id; Step 2523: According to the mathematical principle formula of the source term function, the process variable data and output variable data of the function calculation are obtained through the input variable data; Step 2524: Use Set(id, VarType, inputVal) to write the obtained output variable data into the memory location corresponding to the output variable of the source item function in the basic node unit data container numbered id.

7. A parallel computing method applicable to a two-fluid six-equation system as claimed in claim 1, characterized in that: The process of assembling the coefficient matrix of the mass and energy conservation equations and the coefficient matrix of the momentum conservation equation in step 3 includes: in the thermal hydraulic system model, assuming that there are N non-boundary control bodies and M non-boundary nozzles, then there are 5N control body variables, including pressure, cavitation fraction, vapor phase internal energy, liquid phase internal energy, non-condensable gas fraction and 2M nozzle variables, including vapor phase flow rate and liquid phase flow rate, and the discrete formats of scalar equations such as mass conservation and equation, mass conservation difference equation, vapor phase energy conservation equation, liquid phase energy conservation equation, non-condensable gas transport equation of the N control bodies are respectively organized in the following form: Among them, a1 to a5 are The variable coefficients in the scalar equation, s1 to s4, are The variable coefficient in the scalar equation, b is the constant term of the equation. Subscript g represents the vapor phase, subscript f represents the liquid phase, subscript L represents the current control body, subscript j represents the pipe number, j+1 represents the downstream pipe, superscript n represents the nth time moment, superscript n+1 represents the n+1th time moment, a1 to a5, s1 to s4, and b are all calculated by the source term function and stored in the data container of the corresponding control body; the scalar equations of N control bodies are combined into a matrix form: AX=SV+B (6) in, A is a 5N×5N coefficient matrix for solving scalars such as control volume pressure, cavitation fraction, vapor phase internal energy, liquid phase internal energy, and non-condensable gas fraction. X is an array of 5N control volume variables. S is a 5N×2M coefficient matrix for solving vectors such as vapor and liquid flow rates of upstream and downstream nozzles of the control volume. V is an array of 2M nozzle variables. B is an array of 5N constant terms. According to the topological structure of the thermal hydraulic system, the momentum conservation sum equation and momentum conservation difference equation of M nozzles are arranged in the following form: Among them, g1, g2, w1, and w2 are The variable coefficients in the equation, f is the constant term of the equation, g1, g2, w1, w1, f are all calculated by the source term function and stored in the data container of the corresponding receiver; The momentum conservation equations for the M nozzles are combined into a matrix form: GV=WP+F (8) in, G is a 2M×2M matrix, V is an array of vapor and liquid flow rate variables of 2M nozzles, W is a 2M×N matrix, P is an array of pressure variables of N control bodies, and the coefficient matrices A, S, G and W are sparse matrices. The row elements of the matrices correspond to the coefficients of all solved variables in a certain conservation equation. When a variable is not in this equation, its coefficient is 0. According to the discrete format of the two-fluid six equations, the scalar conservation discrete equations of the control body's mass, energy, non-condensable gas, etc. only contain the variables and parameters of the control body and the upstream and downstream nozzles, while the momentum conservation discrete equation of the nozzle only contains the variables and parameters of the nozzle and the upstream and downstream control bodies.

8. A parallel computing method applicable to a two-fluid six-equation system as claimed in claim 1, characterized in that: The block division and inversion operation of the coefficient matrix of the two-fluid six equations in step 4 includes adjusting the variable array so that the variables of the same basic node unit are adjacent in the variable array, and centralizing the non-zero elements of the coefficient matrix of the scalar equation, that is: Formula (11) can be simplified as: Where: A1 and A2 are 5×5 matrices, X1 and X2 are basic node 5×1 scalar arrays, S1 and S2 are 5×4 matrices, V1 and V2 are basic node 2×1 vector arrays, and B1 and B2 are 5×1 constant arrays. By inverting equation (12), we can obtain: Without loss of generality, under the condition of N non-boundary control volumes and M non-boundary takeovers, equation (13) can be written as: Among them, X1 to X N A 5×1 scalar array for the base node, A1 to A N is a 5×5 matrix, S1 to S N is a 5×2M matrix, V1 to V M A 2×1 vector array of basic nodes, B1 to B N is a 5×1 constant array; From formula (14), we can see that after the scalar equation AX = SV + B is divided into blocks, its time complexity is reduced from O(125N 3 ) becomes O(125N), and the theoretical computational cost increases linearly with the number of control bodies; Because A i (i=1,…,N) matrix element values ​​and B i (i=1, ..., N) can be calculated by the parameter variables stored in the data container corresponding to the current control body, and S i In (i=1, ..., N), only the columns corresponding to the pipes connected to the current control body have non-zero elements, so and The calculation of is only related to the current control body and its connected takeover, which meets the conditions for parallel calculation. i is a 5×5 matrix, and its inverse matrix It is easy to obtain and can be quickly confirmed through the topological structure information of the current control body and the takeover The position of the non-zero element column further improves the computational efficiency; Similarly, by adjusting the variable array, the variables of the same basic node unit are placed adjacent to each other in the variable array, and the non-zero elements of the coefficient matrix of the scalar equation are concentrated: Formula (15) can be simplified as: Among them: G1 and G2 are 2×2 matrices, X1 and X2 are 2×1 vector arrays of basic node elements, W1 and W2 are 2×2 matrices, P1 and P2 are basic node element pressures, and F1 and F2 are 2×1 constant arrays; Perform the inverse operation on equation (16), then: Without loss of generality, under the condition of N non-boundary control volumes and M non-boundary takeovers, equation (17) can be written as: Among them, V1 to V M A 2×1 vector array of basic nodes, G1 to G M is a 2×2 matrix, W1 to W M is a 2×N matrix, F1 to F N is a 2×1 constant array; Similarly, by dividing the vector matrix GV=WP+F into blocks, the time complexity is reduced from O(4M 2 ) becomes O(4M), and the computational cost increases linearly with the number of takeovers. and The calculation is only related to the current takeover and its upstream and downstream control bodies, and has the conditions for parallel calculation.

9. A parallel computing method applicable to a two-fluid six-equation system as claimed in claim 1, characterized in that: The pressure matrix equation generation in step 5 includes extracting all matrix element expressions related to the control volume pressure in formula (14) after dividing and inverting the coefficient matrix of the six equations of the two fluids, that is, X i , and The first line of the control body and the connection topology information of the takeover are used to determine The first row of non-zero elements is distributed, and the matrix elements of (14) with respect to pressure are expressed as follows: Where j is the inlet pipe number, m is the total number of inlet pipes, jp is the outlet pipe number, n is the total number of outlet pipes, and b i for About Pressure P i The constant term, C P，g，j , C P，f，j , C P，g，jp , C P,f,jp for Non-zero elements in the first row; Furthermore, the matrix element expression of the flow rate of the inlet and outlet pipes of the i-th control body in formula (18) is extracted, and the following form is obtained by taking the inlet and outlet pipes of the j-th and jp-th as examples: V g，j =f V，g，j +C V，g，K，j P K，j +C V，g，i，j P i (20) V f，j =f V，f，j +C V,f，K,j P K，j +C V，g，i，j P i (21) V g，jp =f V，g，jp +C V，g，i，jp P i +C V，g，L，jp P L，jp (22) V f，jp =f V，f，jp +C V，f，i，jp P i +C V，f，L，jp P L，jp (23) Among them, f V，g，j and f V，f，j They are The first and second rows of constant terms, C V，g，K，j and C V，f，K，j They are The first and second rows of the matrix elements of the pressure of the upstream control volume of the jth nozzle, C V，g，i，j and C V，f，i，j They are The first and second rows of the matrix elements of the pressure of the i-th control volume, f V,g,jp and C V，f，jp They are The first and second rows of constant terms, C V，g，i，jp and C V，f，i，jp They are The first and second rows of the matrix elements of the pressure of the i-th control volume, C V，g，L，jp and C V，f，L，jp They are The first and second rows of the matrix elements are about the pressure of the downstream control body of the jpth takeover; Substituting equations (20) to (23) into equation (19), we can obtain: All coefficients and constants in equation (24) can be obtained by The matrix is ​​obtained. Further, by traversing according to the control body number, equation (24) is converted into a matrix equation. A p P=b p (25) A p is a large sparse matrix, A p The non-zero elements of the i-th (i=1, ..., N) row of the matrix only include the i-th column and the number of columns corresponding to the pressure of the non-boundary control body connected to the i-th control body by the pipe. The number of non-zero elements is m+n+1, where m is the total number of inlet pipes and n is the total number of outlet pipes, which is much smaller than the number of control bodies N. The pressure P of each control body in the thermal hydraulic system can be obtained by solving formula (25): n+1 , and then calculate V by formula (18) n+1 Finally, other scalars except pressure are calculated by equation (14) 10. A parallel computing method applicable to a two-fluid six-equation system as claimed in claim 9, characterized in that: The step 5 also includes establishing a pressure equation for the thermal hydraulic system, the structure of which is shown in (26), Perform a linear transformation on equation (26), that is, adjust the order of the variables in the variable array and bring the control volume variables of the same loop together, as shown in equation (27). Formula (27) is simplified to: in: and They are the block matrices of the pressure equation corresponding to loop 1 and loop 2 respectively; and The pressure variable array for loop 1 and loop 2; and is the constant term array of loop 1 and loop 2; Considering k independent loops, we get: Right now:

11. A parallel computing method applicable to a two-fluid six-equation system as claimed in claim 10, characterized in that: The step 5 includes the following steps: Step 51: Use a breadth-first traversal algorithm to number the basic node units of the thermal hydraulic system, and generate a two-layer number including a loop number and a node number; Step 52: Without distinguishing the loop numbers of the basic node units, a parallel calculation method is adopted to traverse the node numbers of the control body and calculate the coefficients required for assembling the pressure equation coefficient matrix; Step 53: According to the loop number of the basic node unit, traverse the node number of the basic node unit under the loop number in turn, obtain the pressure matrix coefficient corresponding to the basic node unit of the control body class, and assemble the block coefficient matrix of the pressure equation of the loop number.

12. A parallel computing method applicable to a two-fluid six-equation system as claimed in claim 11, characterized in that: The step 51 includes the following steps: Step 511: changing the directed edges of the flattened topological structure graph of the thermal hydraulic system into undirected edges, wherein the flattened topological structure graph is a system model that restores the thermal hydraulic system model constructed by the user into a system model that is completely composed of basic node unit models and is connected in an orderly manner of control body-takeover-control body; Step 512: establishing a data container A containing graph vertices of an undirected graph of a flattened topological structure of a thermal hydraulic system, wherein the data container A stores the names of basic node units of all graph vertices and the names of the model components to which the basic node units belong in a queue manner; Step 513: Get the first graph vertex stored in the data container, the loop number of the graph vertex is marked as 1, and store it in a temporary data container B. The temporary data container B expresses the order in which the graph vertex elements are stored in the form of a queue; Step 514: Get the first graph vertex element of the temporary data container B, and traverse all graph vertex elements connected to the graph vertex element, which are divided into the following cases: Case 1: The connected graph vertex elements have been numbered, and the traversal continues until the traversal ends; Case 2: The connected graph vertex elements are not numbered and are queued and stored in a temporary data container B. The traversal continues until the traversal is completed. Step 515: After the traversal of the connected graph vertex elements is completed, the first graph vertex element is deleted from the temporary data container B and stored in the result data container C; Step 516: Record the type of the graph vertex element and the sequence number of the graph vertex elements of different types stored in the result data container C, wherein the graph vertex element types include control bodies and takeovers, and use the sequence number of the control body class graph vertex element stored in the result data container C as the node number of the control body, and use the sequence number of the takeover class graph vertex element stored in the result data container C as the node number of the takeover; Step 517: Determine whether the graph vertex element has been marked with a loop number, which can be divided into the following cases: Case 1: The graph vertex element is not marked with a loop number, and the loop number of the graph vertex is marked with the loop number of the graph vertex with the previous sequence number of the result data container C, and then continue to the next step; Case 2: The vertex elements of the graph have been marked with loop numbers, and the next step is continued; Step 518: Determine whether the temporary data container B is empty, which can be divided into the following cases: Case 1: Temporary data container B is not empty, jump to step 514; Case 2: Temporary data container B is empty, and the next step is continued; Step 519: traverse the graph vertex elements of data container A to determine whether there are graph vertex elements in data container A that are not node-numbered, which is divided into the following two cases: Case 1: There are graph vertex elements that are not node-numbered. Randomly select a graph vertex element that is not node-numbered and store it in the temporary data container B. The loop number of the graph vertex element is marked as the loop number of the graph vertex element that is currently stored in the result data container C last + 1, and jump to step 514; Case 2: If there are no graph vertex elements that are not node-numbered, the basic node element numbering ends.

13. A parallel computing method applicable to a two-fluid six-equation system as claimed in claim 1, characterized in that: The specific calculation process in step 6 is as follows: Assuming that the independent loop has N control bodies, the non-zero elements of the i-th (i=1, ..., N) row of the pressure matrix only include the i-th column and the number of columns corresponding to the pressure of the non-boundary control body connected to the i-th control body by a pipe, and a direct solution method based on the sparse LU decomposition method is used as the solution algorithm of the pressure matrix equation. The sparse LU decomposition method aims to obtain the elementary change matrix Q, the lower triangular matrix L and the upper triangular matrix U, so that: QA=LU (31) Among them, the elementary change matrix Q extracts the element with the maximum absolute value of each column as the elimination principal element through basic row transformation during the LU decomposition process, and the pressure equation block matrix (30) is converted into: is the elementary transformation matrix, It is an array The elements are reordered and then the lower triangular matrix and the upper triangular matrix Structural features to obtain independent circuit pressure distribution 14. A parallel computing method for a two-fluid six-equation system as claimed in claim 1, characterized in that: The vapor-liquid flow rate of the pipe in step 7 is calculated as follows: The momentum conservation equation vector array block matrix calculation formula (18) based on the pipe is called to control the body pressure P n+1 As input, calculate the flow rate of each pipe in the thermal hydraulic system Store updates to each taken-over data container; Where: the superscript n+1 represents the calculated value at the current moment, n represents the calculated value at the previous moment, and M is the number of non-boundary pipes in the thermal-hydraulic system.

15. A parallel computing method applicable to a two-fluid six-equation system as claimed in claim 1, characterized in that: In step 8, the cavitation fraction, vapor-liquid ratio internal energy and non-condensable gas fraction of the control volume are calculated as follows: The scalar array block matrix calculation formula (14) based on the mass and energy conservation equation of the control volume is called, and the vapor-liquid flow rate V of the thermal hydraulic system is taken over. n+1 As input, calculate the cavitation fraction of the control volume of the thermal hydraulic system Steam internal energy Liquid phase internal energy Non-condensable gas fraction Store updates into the data container of each control body; 16. A parallel computing method applicable to a two-fluid six-equation system as claimed in claim 1, characterized in that: The mass error of the solution of the two-fluid six-equation system in step 9 is the deviation between the vapor-liquid mixture density obtained by the first-order Taylor expansion and the vapor-liquid mixture density calculated by the physical property state equation, and its expression is: e m =max(e mlocal ,he mtotal ) (35) r i =a g，i r g，i +(1-a g，i )r f，i (39) in, and is the vapor phase density and liquid phase density obtained by the first-order Taylor expansion of the ith control volume; ρ g，i and ρ f，i is the vapor phase density and liquid phase density obtained by the physical property state equation or interpolation table of the ith control volume; k is a factor greater than 1; N is the number of non-boundary control volumes; Set the upper limit of mass error MassError_upBond and MassError_lowBond to judge the value of the current mass error and the upper limit of mass error and the lower limit of mass error, which are divided into the following situations: Case 1: If the current mass error is greater than the set upper limit of the mass error, the time step is reduced, the simulation time is rolled back to the previous moment, and the state variable values ​​stored in each control body and the takeover data container are rolled back to the values ​​of the previous moment, that is, rolled back to the values ​​before the update, and the calculation and solution of the two-fluid six-equation system are restarted; Case 2: If the current mass error is less than the set mass error lower limit, increase the time step and continue time advancement; Case 3: If the current mass error is between the set upper and lower limits of the mass error, the time step is not adjusted and time advancement continues.

Images