Two-fluid six-equation coefficient matrix assembly and storage method suitable for Modelica

Abstract

The invention belongs to the technical field of thermal hydraulic power, and particularly relates to a two-fluid six-equation coefficient matrix assembly and storage method suitable for Modelica. Comprising the following steps: converting a hierarchical thermal hydraulic system model into a planarization topological connection diagram; traversing a topological structure, and uniformly numbering control bodies and takeover nodes; calculating parameters required by an assembly coefficient matrix based on the node data; a sparse coefficient matrix of the mass / energy conservation equation is assembled according to equation and variable sorting, and is stored in a row compression format; and aiming at the takeover nodes, assembling a sparse coefficient matrix of the momentum conservation equation, and storing the sparse coefficient matrix in a row compression format. Through topology planarization, unified numbering and block sparse storage, efficient assembly and optimal management of complex system equations are realized. The method has the beneficial effects that the problem of automatic generation of the planarization topological connection general drawing of the thermal hydraulic system model under the unified modeling framework based on Modelica is solved, and complete input is provided for two-fluid six-equation coefficient matrix assembly.

Description

Technical Field

[0001] The present invention belongs to the technical field of thermal hydraulics, and in particular relates to a coefficient matrix assembly and storage method for two-fluid six-equation equations applicable to Modelica. Background Art

[0002] In order to overcome the limitations of native Modelica in expressing field equations represented by two-fluid six-equation system, the FMU model framework is introduced and the syntax and semantics of Modelica are expanded. In addition, a thermal-hydraulic hierarchical modeling method based on the basic node unit model is established within the FMU model framework. This not only solves the problems of cross-merger of algebraic loops and complicated coupling logic in the unified compilation and analysis process of the two-fluid six-equation system and the differential-algebraic equation system, but also solves the problem that the modeling method based on native Modelica is difficult to achieve reconfiguration and reuse of thermal-hydraulic models, and that a large number of variable transfer equations must be defined in the fluid interface of the thermal-hydraulic model to realize the flow constraint relationship between adjacent nodes. Summary of the invention

[0003] The purpose of the present invention is to provide a coefficient matrix assembly and storage method for two-fluid six-equation suitable for Modelica, which can solve the following problems: first, how to automatically generate a flattened topological connection diagram composed of basic node unit models according to the connection relationship of thermal-hydraulic component models, and automatically number each node in the flattened topological connection diagram, and record the structural parameters, initial parameters and the type and name of the source term closed function of each node according to the node number, so as to provide complete input for the coefficient matrix assembly of the two-fluid six-equation; second, how to automatically assemble the coefficient matrix of the two-fluid six-equation according to the flattened topological structure composed of the basic node unit models, and automatically establish the association relationship between the coefficients of each unit of the matrix and the structural parameters, characteristic parameters, initial parameters, state parameters, generalized source term functions, conservation equation variable coefficients and other data information stored in each basic node data container, so as to achieve the completeness, closure and real-time performance of the coefficient matrix and ensure the efficient storage of the coefficient matrix.

[0004] The technical solution of the present invention is as follows: a method for assembling and storing coefficient matrices of two-fluid six-equation equations applicable to Modelica, comprising the following steps:

[0005] Step 1: Automatically generate the topology of the thermal-hydraulic system, and restore the thermal-hydraulic system model with hierarchical encapsulation characteristics to a flattened topological connection diagram consisting entirely of orderly connected basic node unit models;

[0006] Step 2: Traverse the basic node units of the flattened topological structure connection diagram corresponding to the thermal-hydraulic system model, and uniformly number the control body type basic node units and the pipe type basic node units according to the direction, branches, and loops of the orderly connection of the topological structure;

[0007] Step 3: Based on the flattened topology and basic node unit numbers corresponding to the thermal-hydraulic system model, the data information of the control body basic node units and the nozzle basic node units are processed and calculated to obtain the complete data information required for assembling the coefficient matrix of the two-fluid six equations;

[0008] Step 4: Based on the basic node units of the control body, according to the topological structure and the set equation and variable sorting, read the data information of the basic node units, and assemble the coefficient matrix of the discrete equations of mass and energy conservation. The coefficient matrix consists of the coefficient matrix of scalar solutions such as the pressure of the control body, the cavity fraction, the vapor phase internal energy, and the liquid phase internal energy; the coefficient matrix of vector solutions such as the vapor phase flow rate and liquid phase flow rate of the upstream and downstream pipes of the control body; and the constant term. The scalar coefficient matrix and the vector coefficient matrix are sparse matrices and are stored in a row-compressed format.

[0009] Step 5: Based on the basic node units of the takeover type, according to the topological structure and the set equation sorting and variable sorting, read the data information of the basic node units, and assemble the coefficient matrix of the momentum conservation discrete equation. The coefficient matrix is composed of the coefficient matrix of the solution vectors such as the vapor phase flow rate and liquid phase flow rate of the takeover, the coefficient matrix of the upstream and downstream control body pressures of the takeover, and the constant term. The solution vector coefficient matrix and pressure coefficient matrix are sparse matrices and are stored in row compression format.

[0010] The method for automatically generating the topological structure of the thermal hydraulic system model in step 1 comprises the following steps:

[0011] Step 11: Identify the Modelica-based thermal-hydraulic model through model annotation keywords. The model annotation keywords are divided into basic node unit model identification keywords and thermal-hydraulic component model identification keywords. The model annotation keywords belong to the expanded Modelia syntax semantics, that is, by introducing the FMU model framework to support thermal-hydraulic hierarchical modeling based on basic node units; the basic node unit model includes the control body model v, the takeover model j, the time-dependent control body model tdv, and the time-dependent takeover model tdj. The thermal-hydraulic component model refers to the thermal-hydraulic model of pipes, containers and other internally defined discrete nodes, as well as the higher-level thermal-hydraulic model formed by connecting and encapsulating the thermal-hydraulic component models;

[0012] Step 12: Thermal-hydraulic component model flattening. The thermal-hydraulic component model is converted into a component model composed of control volumes, nozzles, flow interfaces, and directed connection edges according to the discrete node number and hierarchical encapsulation information.

[0013] For the pipe model, according to the discrete node number N defined by the pipe model, the pipe model is converted into a component model with N control volumes, N - 1 nozzles, 4N - 2 basic node unit flow interfaces, 2 component model flow interfaces, and 2N directed connection edges connected by connectors.

[0014] For the pipe model and the thermal-hydraulic model after pipe model encapsulation, the model is decomposed into pipe models along the hierarchical encapsulation information, and then the pipe is decomposed into a model composed of control volumes, nozzles, component interfaces, and directed connection edges.

[0015] Step 13: According to the model interface connection relationship in the thermal-hydraulic system modeling process, the compiler connects the basic node unit model and the flattened thermal-hydraulic component model to form a general model connection diagram. The vertices of the diagram include the basic node unit model and component interfaces. Among them, the component interfaces are the inlet flow interface and outlet flow interface of the thermal-hydraulic component model. The vertices of the diagram are stored using the data container NumapNode, and each vertex corresponds to a data container. The data container NumapNode corresponding to the basic node unit model stores the data information required for the solution of the two-fluid six-equation, including the model name, structural parameters, characteristic parameters, initial parameters, state parameters, upstream and downstream node information, and generalized source terms where the vertex is located. The data container NumapNode corresponding to the component interface only stores the thermal-hydraulic component model information, including the component name and interface type, and does not store the data information required for the solution of the two-fluid six-equation.

[0016] Step 14: Traverse the vertices of the general model connection diagram to obtain the type and information of the vertex NumapNode, as well as the number and information of the upstream and downstream vertices. According to the connection rules of the thermal-fluid model interface, when the vertex is a component interface, check the number of upstream and downstream vertices to confirm whether both are 1. If not, generate an error message and jump out of the topological structure generation logic. If so, delete the vertex and the directed connection edges connecting it to the upstream and downstream vertices, and then generate a directed connection edge from the upstream vertex K_Node to the downstream vertex L_Node.

[0017] Step 15: Traverse the directed connection edges of the flattened topology graph. If the two graph vertices NumapNode types connected by the directed connection edge are both takeover graph vertices, including takeover j or time-related takeover tdj, then output an error message; traverse the takeover graph vertices of the flattened topology graph. If the number of directed connection edges at the import or export of the takeover graph vertex is less than 0.00, then output an error message; if the two graph vertices K_Node and L_Node types connected by the directed connection edge are both control body v or time-related control body tdv, then delete the directed connection between the graph vertices K_Node and L_Node. Connecting edges: add a graph vertex of type pipe j between two graph vertices, and then establish graph vertices K_Node and L_Node points and directed connection edges between the newly added graph vertices based on the upstream and downstream relationships. The flow area and hydraulic equivalent diameter of the added graph vertex of type pipe j are equal to the flow area and hydraulic equivalent diameter of the graph vertex with the smallest flow area among graph vertices K_Node and L_Node. The length of pipe j is equal to half of the sum of the lengths of graph vertices K_Node and L_Node. By traversing the directed connection edges of the flattened topology graph, a flattened topology graph in the form of "VJV" is formed.

[0018] The specific steps of numbering the basic node units in step 2 are as follows:

[0019] Step 21: Traverse the graph vertices of the flattened topology connection graph, obtain the loop import graph vertices and their number, and select a loop import graph vertex as the first numbering starting point; if there is no loop import graph vertex, randomly select a control body class graph vertex as the first numbering starting point, and the first numbering starting point is recorded as V1. The loop import graph vertex refers to the control body class graph vertex with only outgoing directed connection edges. The numbering starting point refers to the initial starting point when the numbering is advanced along the topology connection graph, that is, there is no upstream graph vertex in the numbering advancement route;

[0020] Step 22: Using the first numbered starting point as the path starting point, determine the first path for numbered advancement. The first junction class graph vertex of the path is denoted as J1. If there is a branch downstream of the first numbered starting point, select a subpath of the numbered starting point as the first path for numbered advancement. If there is no branch downstream of the first numbered starting point, select the path where the outlet junction of the numbered starting point is located as the first path for numbered advancement. The path refers to a single line without branches in the topological structure connection diagram that uses a VJV directed connection, expressing the flow route of the fluid from upstream to downstream in a single pipeline in the thermal-hydraulic system.

[0021] Step 23: Advance along the VJV directed connection of the path, sequentially number the control body class graph vertices and the takeover class graph vertices respectively until the end of the path, and return to the starting point of the path. The end of the path is the control body class graph vertex, including branch points, numbered graph vertices and outlet graph vertices. The branch point refers to the control body class graph vertex with two or more outlet connection takeovers, and the outlet graph vertex is the control body class graph vertex with only inlet directed connection edges. The above logic of numbering and advancing along the path is defined as the path advancement numbering logic S1, which is applicable to the numbering advancement of the first path and all subsequent paths;

[0022] Step 24: After returning to the path starting point, determine whether the path starting point is a branch point. If not, return directly to the numbering starting point. At this time, the numbering starting point and the path starting point are the same graph vertex. If so, it is necessary to perform path advancement numbering along other subpaths of the path starting point. The execution method adopts the path advancement numbering logic S1 until all subpaths of the path starting point have performed path advancement numbering. Then return to the first branch point upstream of the path starting point, use this branch point as the path starting point, perform path advancement numbering on all subpaths of this branch point, and then return to the branch point further upstream to perform the above path advancement numbering operation until there are no upstream branch points and return to the numbering starting point. The above numbering advancement logic that incorporates the influence of branch points is defined as branch coverage logic S2. Logic S2 includes path advancement numbering logic S1 and is applicable to the numbering advancement of vertices in the topological structure graph with the numbering starting point as the upstream initial point.

[0023] Step 25: After returning to the numbering starting point, select other loop import graph vertices as the next numbering starting point to continue numbering the graph vertices. The numbering advancement logic adopts branch covering logic S2 until all loop import graph vertices are traversed as numbering starting points. If all loop import graph vertices have been selected as numbering starting points, or the flattened topology connection graph does not contain loop import graph vertices that have not been selected as numbering starting points, then the graph vertices of the flattened topology connection graph are traversed to determine whether all graph vertices have been numbered. If so, the numbering ends; if not, it means that there are multiple loops, then the branch covering logic S2 is executed, and the graph vertices of the flattened topology connection graph are traversed again and the branch covering logic S2 is executed until all graph vertices have been numbered. The above-mentioned numbering advancement logic that incorporates the influence of loop import graph vertices and multiple loops is defined as import graph vertex and loop covering logic S3. Logic S2 includes path advancement numbering logic S1 and branch covering logic S2, which is applicable to the numbering advancement of the entire flattened topology structure connection graph.

[0024] In step 3, the assembly of the two-fluid six-equation coefficient components is carried out based on the discrete forms of the mass conservation sum equation, the mass conservation difference equation, the vapor-phase energy conservation equation, the liquid-phase energy conservation equation, the momentum conservation sum equation, and the momentum conservation difference equation.

[0025] The discrete forms of the equations in step 3 are as follows:

[0026] (1) Mass conservation sum equation

[0027]

[0028] (2) Mass conservation difference equation

[0029]

[0030] (3) Vapor-phase energy conservation equation

[0031]

[0032] (4) Liquid-phase energy conservation equation

[0033]

[0034] (5) Momentum conservation sum equation

[0035]

[0036] (6) Momentum conservation difference equation

[0037]

[0038] indicates that the variable is a temporary result, indicates the use of the upwind scheme, that is:

[0039] where V is the control volume, A is the cross-sectional area of the control volume, P is the pressure, T is the temperature, U is the specific internal energy, H is the heat transfer coefficient, ρ is the density, α is the void fraction, v is the flow velocity, ∈ indicates whether it is wall boiling or wall condensation, Δt is the time step, Δx is the control volume length, the subscript g represents the vapor phase, the subscript f represents the liquid phase, the subscript L represents the downstream control volume, the subscript K represents the upstream control volume, the subscript j represents the nozzle number, j - 1 represents the upstream nozzle, 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;

[0040] After semi-implicit discretization, the unknowns in equations (1) to (4) include those at the current (n + 1)th moment Among them, the first five variables do not belong to the solution variables of the two-fluid six equations and need to be expanded using Taylor series, i.e.:

[0041]

[0042]

[0043] Substitute equations (7) to (11) into equations (1) to (4), and through variable coefficient arrangement and like-term combination, equations (1) to (4) are uniformly expressed in the following form:

[0044]

[0045] Among them, a1 to a4 are respectively the variable coefficients in the equation, s1 to s4 are respectively the variable coefficients in the equation, and b is the constant term of the equation. Among them, s1 to s4 are the convective transport coefficients of the mass or energy of the pipe connection, which are related to the control volume number, pipe connection number, and equation number, and can be written as: s n，j，q，k , where n is the control volume number, j is the pipe connection number, q is the equation number, and k is 1 for the vapor phase g or 2 for the liquid phase f;

[0046] In the momentum conservation equation (5) and the momentum conservation difference equation (6), the unknown variables only include those at the current n + 1 moment Through variable coefficient arrangement and like-term combination, equations (5) and (6) are uniformly expressed in the following form,

[0047]

[0048] Among them, g1, g2, w1, and w2 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; To provide the complete data information required for the assembly of the two-fluid six equation coefficient matrix, combining equations (12) and (13), the content stored in the data container corresponding to the basic node element includes the model name where the basic node element is located, the type of the basic node element, the number of the basic node element, the upstream and downstream node information, the structural parameters, the characteristic parameters, the initial parameters, and the generalized source term. The generalized source term refers to the physical source terms 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 physical property calculation, pipe connection variable calculation based on control volume variables, and control volume variable calculation based on pipe connection variables.

[0049] In step 3, according to the flattened topological structure corresponding to the thermal-hydraulic system model and the basic node unit numbers, calculate and update in real time the data information of the basic node units of the control volume type and the basic node units of the pipe connection type, and obtain the complete data information required for assembling the two-fluid six-equation coefficient matrix. The specific calculation steps are as follows:

[0050] Step 31: Obtain the upstream and downstream pipe connection numbers of the basic node units of the control volume type according to the overall diagram of the flattened topological structure, so as to establish the constraint relationship between the solution variables of pressure, void fraction, vapor specific internal energy, and liquid specific internal energy based on the control volume number and the solution variables of vapor velocity and liquid velocity based on the pipe connection number in the mass sum-difference conservation discrete equation and the energy conservation discrete equation;

[0051] Step 32: Obtain the upstream and downstream control volume numbers of the basic node units of the pipe connection type according to the overall diagram of the flattened topological structure, so as to establish the constraint relationship between the solution variables such as vapor velocity and liquid velocity based on the pipe connection number and the pressure variable based on the control volume number in the momentum sum-difference conservation discrete equation;

[0052] Step 33: The vapor velocity and liquid velocity of the control volume required by the momentum conservation equation are obtained by averaging the vapor-liquid velocities of the upstream and downstream pipe connections and stored in the corresponding data container of the basic node unit of the control volume type;

[0053] Step 34: The physical property parameters such as density, viscosity coefficient, thermal conductivity, and specific heat at constant pressure of the pipe connection, as well as the scalar quantities such as void fraction, vapor specific internal energy, and vapor specific internal energy, are all the physical property parameters and scalar quantities of the incoming flow control volume and are stored in the corresponding data container of the basic node unit of the pipe connection type and directly used in the calculation of the source term closure relation of the momentum conservation discrete equation;

[0054] Step 35: For the basic node units of the control volume type, calculate the generalized source terms of the mass and energy conservation discrete equations, and store the calculation results in the data container of the control volume. The calculation of the generalized source terms includes the calculation of physical properties, interphase mass transfer, interphase heat transfer, wall heat transfer, energy dissipation and other source terms. And according to the upstream and downstream pipe connection numbers of the control volume, calculate the coefficients and constant terms of each solution variable of the mass and energy conservation discrete equations sorted by formula (12), including a1 to a4, s1 to s4, b, which can be directly used when assembling the coefficient matrix. Among them, s1 to s4 need to correspond to the upstream and downstream pipe connection numbers of the control volume, the equation number, and the type of vapor-liquid velocity;

[0055] Step 36: For the basic node units of the takeover type, calculate the generalized source term of the discrete equation of momentum conservation, and store the calculation results in the data container of the control volume. The calculation of the generalized source term includes the calculation of wall friction, interfacial friction, local resistance, artificial viscosity terms, etc., and the treatment of special pumps, critical flows, etc. According to the control volume numbers upstream and downstream of the takeover, calculate the coefficient and constant term of each solution variable of the momentum conservation equation and the difference of the momentum conservation equation sorted by Equation (13), including g1, g2, w1, w2, f, which can be directly used when assembling the coefficient matrix;

[0056] Step 37: There may be multiple upstream and downstream takeovers for one control volume, but each takeover has only one upstream control volume and one downstream control volume. When dealing with the multi-pass situation, it is expressed by Equation (13). When there is one more upstream takeover for the control volume, two terms are added to the mass conservation equation of the control volume and Sorted according to Equation (12), that is:

[0057]

[0058] It can be seen from Equation (14) that the left scalar array, the variable coefficient array of the scalar, and the constant term remain unchanged, while the right vector array increases and The corresponding variable coefficient array of the right vector increases and

[0059] The coefficient matrix of the two-fluid six equations in Step 4 is obtained by sorting according to Equations (12) and (13). The j-th element a ij in the i-th row of the matrix corresponds to the coefficient of the j-th variable in the i-th equation. The j-th variable refers to the variable ranked in the j-th row of the variable array. The row elements of the coefficient matrix correspond to the coefficients of all solution variables. When a certain variable is not in this equation, its coefficient is 0. The large sparse matrix is stored in row-compressed format. In the thermal-hydraulic system model, assume that there are N non-boundary control volumes and M non-boundary takeovers. There are 4N control volume variables, including pressure, void fraction, specific internal energy of the vapor phase, specific internal energy of the liquid phase, and 2M takeover variables, including vapor phase velocity and liquid phase velocity. Sort the mass conservation equation, mass conservation difference equation, vapor phase energy conservation equation, and liquid phase energy conservation equation of N control volumes according to Equation (12), and merge them into matrix form:

[0060] AX = SV + B (20)

[0061] Among them,

[0062]

[0063] A is a 4N×4N coefficient matrix for solving scalar quantities such as control volume pressure, void fraction, vapor specific internal energy, and liquid specific internal energy. X is an array of 4N control volume variables. S is a 4N×2M coefficient matrix for solving vectors such as vapor velocity and liquid velocity at the upstream and downstream nozzles of the control volume. V is an array of 2M nozzle variables. B is an array of 4N constant terms.

[0064] The steps for assembling the coefficient matrix of the mass and energy conservation discrete equations in step 4 are as follows:

[0065] Step 41: Obtain the number of non-boundary control volumes N, the number of non-boundary nozzles M, the number of flow boundary nozzles M_TDJ, and the number of multi-pass extra nozzles M_ExtraConnectJun according to the topological structure information.

[0066] Step 42: Calculate the number of non-zero elements of matrix A as 16N, and the number of non-zero elements of matrix S as 12N + 6.

[0067] Step 43: Apply for the memory of the row-compressed format arrays of matrix A and matrix S according to the number of non-zero elements of matrix A and matrix S, including: A_data, A_colindex, A_rowindex and S_data, S_colindex, S_rowindex, and apply for the memory of the constant term value B.

[0068] Step 44: Set the mass and energy conservation equations and variable sorting as:

[0069]

[0070] Step 45: In the range of non-boundary control volumes, sequentially select control volumes from i = 1, where i is a natural number from 1 to N, obtain the data container of the i-th control volume, and read the structural parameters, characteristic parameters, initial parameters, state parameters, physical property parameters, generalized source term calculation results stored in the data container, and the coefficients and constants of each solution variable of the mass and energy conservation discrete equations sorted according to equation (12), including variable coefficients a1 to a4, s1 to s4, b. The variable coefficients and constants are different in different conservation discrete equations, and s1 to s4 need to correspond to the control volume number, nozzle number, and equation number. s1 to s4 can all be rewritten as: s n，j，q，k , where n is the control volume number, j is the nozzle number, q is the equation number, k is 1 for vapor g or 2 for liquid f. If there is a multi-pass, the number of variable coefficients of the solution vector will increase. Assuming an increase of m upstream or downstream nozzles, then the variable coefficients of the solution vector of each type of equation will increase from s1 to s4 to s1 to s 4+2m ;

[0071] Step 46: Calculate the number of rows and columns of the variable coefficients of the solution scalar in matrix A according to the set equations and variable sorting.

[0072] Step 47: Store the variable coefficients in the row-compressed format array of the A matrix according to the number of rows and columns of the A matrix for the variable coefficients of the solved scalar.

[0073] Step 48: Obtain the upstream and downstream connection information of the i-th control volume according to the topological structure, clarify the numbers and data information of each upstream and downstream connection. Assume the number of the connection is j, and confirm whether the connection numbered j is a flow boundary connection, which is divided into the following situations:

[0074] Situation 1: If the connection numbered j is a flow boundary connection, then read the physical property parameters, structural parameters, gas-liquid flow velocities and other data information of the data container of the connection numbered j, as well as the variable coefficients of the gas-liquid flow velocities of the connection numbered j in the mass conservation sum equation, mass conservation difference equation, gas-phase energy conservation equation and liquid-phase mass conservation equation in the data container of the i-th control volume, calculate the expression term of the gas-liquid flow velocity of the connection numbered j in the equation. This expression term is equal to the gas-liquid flow velocity multiplied by its corresponding variable coefficient in the equation, and add it to the constant term of the corresponding equation of the i-th control volume, update the constant term of the data container, and jump to Step 410;

[0075] Situation 2: If the connection numbered j is a non-boundary connection, then calculate the number of rows and columns of the S matrix for the variable coefficients of the solved vector according to the equation sorting of formula (21) and the variable sorting of formula (22). Taking the mass conservation difference equation as an example, the number of rows of the variable coefficient corresponding to the gas-phase flow velocity of the connection numbered j in the mass conservation difference equation of the i-th control volume in the S matrix is 4i - 2, and the number of columns is j - m_tdj, where m_tdj is the number of flow boundary connections among the connections numbered from 1 to j, and its value is the variable coefficient s of the gas-phase flow velocity of the connection numbered j in the mass conservation difference equation stored in the data container of the i-th control volume n，j，2，1 where: n is the number of the i-th control volume, and j is the number of the connection.

[0076] Step 49: Store the variable coefficients in the row-compressed format array of the S matrix according to the number of rows and columns of the S matrix for the variable coefficients of the solved vector.

[0077] Step 410: Store the constant terms b of the mass conservation sum equation, mass conservation difference equation, gas-phase energy conservation equation and liquid-phase mass conservation equation in the data container of the i-th control volume into the corresponding positions of the constant term array B according to the mass and energy conservation equation sorting.

[0078] Step 411: Repeat Steps 45 to 410, and successively complete the storage of the variable coefficients and constant terms for the control volume numbers i from 1 to N, then the A matrix, S matrix and constant term array B can be assembled and completed.

[0079] The assembly of the coefficient matrix of the momentum conservation discrete equation in step 5 is as follows:

[0080] In the thermal-hydraulic system model, assuming there are N non-boundary control volumes and M non-boundary nozzles, then there are 4N control volume variables, including pressure, void fraction, specific internal energy of the vapor phase, specific internal energy of the liquid phase, and 2M nozzle variables, including vapor velocity and liquid velocity. The momentum conservation equations and momentum conservation difference equations of the M nozzles are sorted according to equation (13) and combined into matrix form:

[0081] GV = WP + F (42)

[0082] Where,

[0083] G is a 2M×2M matrix, V is an array of 2M vapor-phase and liquid-phase velocity variables of the nozzles, W is a 2M×N matrix, and P is an array of pressure variables of the N control volumes.

[0084] The assembly of the coefficient matrix of the momentum conservation discrete equation in step 5 is as follows:

[0085] Step 51: Obtain the number N of non-boundary control volumes, the number M of non-boundary nozzles, and the number N_P of pressure boundary TDVs according to the topological structure information;

[0086] Step 52: Calculate the number of non-zero elements 4M of the G matrix and the number of non-zero elements 4M - 2N_P of the W matrix;

[0087] Step 53: Apply for the memory of the row-compressed format arrays of the G matrix and the W matrix according to the number of non-zero elements of the G matrix and the W matrix, including: G_data, G_colindex, G_rowindex, and W_data, W_colindex, W_rowindex, and apply for the memory of the constant term value F;

[0088] Step 54: Set the momentum conservation equation and variable sorting as:

[0089]

[0090] Step 55: Among the number M of non-boundary nozzles, sequentially select the nozzles with j = 1, where j is a natural number from 1 to M according to the nozzle number, obtain the data container of the j-th nozzle, and read the structural parameters, characteristic parameters, initial parameters, state parameters, physical property parameters, generalized source term calculation results stored in the data container, and the coefficients and constant terms of each solution variable of the momentum conservation equation sorted according to equation (13), including g1, g2, w1,;

[0091] Step 56: Calculate and solve for the number of rows and columns of the variable coefficients of the vector in the G matrix according to the momentum conservation equation sorting of Equation (43) and the variable sorting of Equation (14).

[0092] Step 57: Store the variable coefficients in the row-compressed format array of the G matrix according to the number of rows and columns of the variable coefficients of the solution vector in the G matrix.

[0093] Step 58: According to the topological structure, obtain the upstream control volume and downstream control volume of the j-th nozzle, clarify the numbers and data information of the upstream control volume and downstream control volume. Assume the number of the upstream control volume or downstream control volume is n, and confirm whether the control volume numbered n is a pressure boundary control volume. It is divided into the following situations:

[0094] Situation 1: Assume that the control volume numbered n is a pressure boundary control volume. Then read the physical property parameters, structural parameters, pressure and other data information of the data container of the control volume numbered n, as well as the variable coefficients of the pressure of the control volume numbered n in the momentum conservation sum equation and momentum conservation difference equation in the data container of the j-th nozzle. Calculate the expression term for the pressure of the control volume numbered n. This expression term is equal to the control volume pressure multiplied by its corresponding variable coefficient in the equation, and add it to the constant term of the corresponding equation of the j-th nozzle. Update the constant term of the data container, and jump to Step 510.

[0095] Situation 2: Assume that the control volume numbered n is not a pressure boundary control volume. According to the momentum conservation equation sorting of Equation (43) and the variable sorting of Equation (44), calculate the number of rows and columns of the variable coefficients of the control volume pressure in the W matrix.

[0096] Step 59: Store the variable coefficients in the row-compressed format array of the W matrix according to the number of rows and columns of the variable coefficients of the solution vector in the W matrix.

[0097] Step 510: According to the momentum conservation equation sorting of Equation (43), store the constant terms of the momentum conservation sum equation and momentum conservation difference equation of the data container of the j-th nozzle into the corresponding positions of the constant term array F.

[0098] Step 511: Repeat Steps 55 to 510, and sequentially complete the storage of the variable coefficients and constant terms for the control volume numbers j from 1 to M, and assemble the G matrix, W matrix, and constant term array F.

[0099] The beneficial effects of the present invention are: (1) the present invention solves the problem of automatically generating a flattened topological connection diagram of a thermal-hydraulic system model under a unified modeling framework based on Modelica, that is, automatically generating a flattened topological connection diagram consisting of basic node unit models according to the component connection relationship of the thermal-hydraulic system model, and automatically numbering each node in the flattened topological connection diagram, and recording the type, structural parameters, characteristic parameters, initial parameters, generalized source term function, conservation equation variable coefficients, etc. of each basic node unit according to the node number, providing complete input for the coefficient matrix assembly of the two-fluid six equations.

[0100] (2) The method for automatically generating a topological connection diagram of the present invention can adapt to the generation of topological structures of thermal-hydraulic system models with multiple loops, multiple thermal-hydraulic boundaries, and multiple branch structures. The basic node unit numbering of the present invention is unique and fully covered, and the numbering advancement path is well consistent with the fluid flow direction. It can quickly retrieve and traverse the numbering of the basic node units according to the flow direction, thereby improving the efficiency of vertex traversal of the flattened topological connection diagram of the thermal-hydraulic system.

[0101] (3) The present invention solves the problem of automatically assembling the coefficient matrix of the two-fluid six-equation based on the flattened topological connection diagram of the thermal-hydraulic system. According to the mathematical form of the discrete format of the two-fluid six-equation, the present invention automatically establishes the correlation relationship between the coefficients of each unit of the matrix and the structural parameters, initial parameters, state parameters, generalized source term functions, conservation equation variable coefficients and other data information stored in each basic node data container, thereby realizing the closedness, completeness and real-time performance of the coefficient matrix, supporting parallel calculation of the coefficient matrix assembly, and supporting efficient storage of the coefficient matrix. BRIEF DESCRIPTION OF THE DRAWINGS

[0102] Figure 1 This is an example of flattening the pipeline model;

[0103] Figure 2 This is an example of flattening the pipeline model and the thermal hydraulic model after encapsulation of the pipeline model;

[0104] Figure 3 Generate logic for flattening the topology of thermal hydraulic models;

[0105] Figure 4 Generate logic for flattening the topology of thermal hydraulic models;

[0106] Figure 5 An example of the effect of generating a flattened topology diagram for a thermal hydraulic model;

[0107] Figure 6 Select the logic diagram for the numbering starting point;

[0108] Figure 7 Select the logic diagram for the first path of number advancement;

[0109] Figure 8 is the path propulsion logic diagram S1;

[0110] Figure 9 is the branch coverage logic diagram S2;

[0111] Figure 10 is the inlet diagram vertex and loop coverage logic diagram S3;

[0112] Figure 11 is an example of the basic node unit numbering of the flattened topological connection general diagram under multiple loops and multiple branches;

[0113] Figure 12 is the thermohydraulic system topological structure. Specific implementation manners

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

[0115] A method for assembling and storing the coefficient matrix of two-fluid six equations applicable to Modelica includes the following steps:

[0116] Step 1: Automatically generate the thermohydraulic system topological structure, and restore the thermohydraulic system model with hierarchical encapsulation characteristics to a flattened topological connection general diagram completely composed of orderly connected basic node unit models, specifically including the following:

[0117] In the thermohydraulic hierarchical modeling method based on basic node units within the FMU model framework, there are both basic node unit models such as control volume v, connection pipe j, time-dependent control volume tdv, and time-dependent connection pipe tdj, and there are also equipment models such as pipes, containers, and heat exchangers assembled and encapsulated by basic node unit models, as well as models with higher encapsulation levels such as steam generators and reactor pressure vessels. Since the matrix equation of the two-fluid six equations is obtained by converting the conservation control discrete equations based on basic node units and does not have the characteristics of hierarchical encapsulation, and the thermohydraulic system model constructed by the user is formed by connecting thermohydraulic component models with different encapsulation levels, it is necessary to restore the thermohydraulic system model constructed by the user to a system model completely composed of orderly connected basic node unit models through flattening processing, so as to make the topological structure of the thermohydraulic system model completely correspond to the discrete grid of the two-fluid six equations and obtain the upstream and downstream connection relationships between the control volume and the connection pipe.

[0118] The method for automatically generating the thermohydraulic system model topological structure includes the following steps:

[0119] Step 11: Identify the Modelica-based thermal-hydraulic model through model annotation keywords. The model annotation keywords are divided into basic node unit model identification keywords (e.g., tpfBaseComp) and thermal-hydraulic component model identification keywords (e.g., tpfPhysicComp). The model annotation keywords belong to the extended Modelia syntax and semantics, that is, by introducing the FMU model framework, hierarchical modeling of thermal-hydraulics based on basic node units is supported. The basic node unit models include control volume models v, nozzle models j, time-dependent control volume models tdv, and time-dependent nozzle models tdj. The thermal-hydraulic component models refer to thermal-hydraulic models such as pipes, vessels, and other models with discretized nodes that can be defined internally, as well as higher-level thermal-hydraulic models formed by connecting and encapsulating thermal-hydraulic component models.

[0120] Step 12: Flatten the thermal-hydraulic component model. Convert the thermal-hydraulic component model into a component model consisting of control volumes, nozzles, flow interfaces, and directed connection edges according to the number of discretized nodes and hierarchical encapsulation information. For example, for a pipe model, according to the number of discretized nodes N defined in the pipe model, convert the pipe model into a component model with N control volumes connected by connectors, N - 1 nozzles, 4N - 2 basic node unit flow interfaces, 2 component model flow interfaces, and 2N directed connection edges, as Figure 1 shown; for a thermal-hydraulic model encapsulated by a pipe model and a pipe model, decompose the model into a pipe model along the hierarchical encapsulation information, and then decompose the pipe into a model consisting of control volumes, nozzles, component interfaces, and directed connection edges, as Figure 2 shown.

[0121] Step 13: According to the connection relationship of the model interfaces in the thermal-hydraulic system modeling process, the compiler connects the basic node unit model and the flattened thermal-hydraulic component model to form a general model connection diagram. The vertices of the diagram include the basic node unit model and component interfaces, where the component interfaces are the inlet flow interfaces and outlet flow interfaces of the thermal-hydraulic component model. The vertices of the diagram are stored using the data container NumapNode, with each vertex corresponding to a data container. The data container NumapNode corresponding to the basic node unit model stores the data information required for solving the six equations of two fluids, including the model name where the vertex is located, structural parameters, characteristic parameters, initial parameters, state parameters, upstream and downstream node information, generalized source terms, etc. The data container NumapNode corresponding to the component interface only stores the thermal-hydraulic component model information, including the component name and interface type, and does not store the data information required for solving the six equations of two fluids.

[0122] Step 14: Traverse the vertices of the overall model connection diagram, and obtain the types and information of the vertex NumapNode, as well as the quantities and information of the upstream and downstream vertices (K_Node, L_Node). According to the connection rules of the thermo-fluid model interface, when the vertex is a component interface, check the quantities of the upstream and downstream vertices (K_Node, L_Node) to confirm whether they are both 1. If not, generate an error message and jump out of the topology generation logic; if so, delete the vertex and the directed connection edges connecting it to the upstream and downstream vertices, and then generate a directed connection edge from the upstream vertex K_Node to the downstream vertex L_Node. After the traversal and processing of the vertices of the overall connection diagram are completed, only the flattened model topology connection diagram formed by the connection of the basic node unit models remains.

[0123] Step 15: Traverse the directed connection edges of the flattened topology overall diagram. If the types of the two graph vertices NumapNode connected by the directed connection edge are both takeover-type graph vertices, including takeover j or time-related takeover tdj, an error message is output; traverse the takeover-type graph vertices of the flattened topology overall diagram. If the quantity of the directed connection edges at the inlet or outlet of the takeover-type graph vertex, an error message is output; if the types of the two graph vertices K_Node and L_Node connected by the directed connection edge are both control volume v or time-related control volume tdv, delete the directed connection edge between the graph vertices K_Node and L_Node, add a graph vertex of type takeover j between the two graph vertices, and then establish directed connection edges between the graph vertices K_Node and L_Node and the newly added graph vertex according to the upstream and downstream relationships. The flow area and hydraulic equivalent diameter of the added graph vertex of type takeover j are equal to the flow area and hydraulic equivalent diameter of the graph vertex with the smallest flow area among the graph vertices K_Node and L_Node, and the length of takeover j is equal to half of the sum of the lengths of the graph vertices K_Node and L_Node. By traversing the directed connection edges of the flattened topology overall diagram, a flattened topology overall diagram in the form of "V-J-V" can be formed.

[0124] Step 2: Traverse the basic node units of the flattened topology connection diagram corresponding to the thermo-hydraulic system model, and uniformly number the basic node units of the control volume type and the takeover type respectively according to the direction, branches, and loops of the topological structure in an orderly connection.

[0125] Traverse the basic node units of the flattened topology connection diagram corresponding to the thermo-hydraulic system model, and uniformly number the basic node units of the control volume type and the takeover type respectively according to the direction, branches, and loops of the topological structure in an orderly connection. The specific steps for numbering the basic node units are as follows:

[0126] Step 21: Traverse the graph vertices of the flattened topological structure connection general graph, obtain the loop inlet graph vertices and their quantities, and select a certain loop inlet graph vertex as the first numbering starting point; if there is no loop inlet graph vertex, randomly select a certain control volume class graph vertex as the first numbering starting point. The first numbering starting point is denoted as V1. The loop inlet graph vertex refers to the control volume class graph vertex with only outlet directed connection edges, and the numbering starting point refers to the initial starting point when advancing along the topological structure connection general graph numbering, that is, there is no upstream graph vertex in the numbering advancement route, and there may be multiple numbering starting points in the flattened topological structure connection general graph.

[0127] Step 22: Take the first numbering starting point as the path starting point, determine the first path for numbering advancement, and denote the first takeover class graph vertex of the path as J1. If there are branches downstream of the first numbering starting point, select a certain sub-path of the numbering starting point as the first path for numbering advancement; if there are no branches downstream of the first numbering starting point, take the path where the outlet takeover of the numbering starting point is located as the first path for numbering advancement. The path refers to a single line with V-J-V directed connection and no branches in the topological structure connection general graph, which represents the flow route of the fluid in a single pipeline in the thermal-hydraulic system from upstream to downstream.

[0128] Step 23: Advance along the V-J-V directed connection of the path, sequentially number the control volume class graph vertices and the takeover class graph vertices in order until the path end point, and then return to the path starting point. The path end point is a control volume class graph vertex, including branch points, numbered graph vertices, and outlet graph vertices. The branch point refers to the control volume class graph vertex with more than two outlet connection takeovers, and the outlet graph vertex is the control volume class graph vertex with only inlet directed connection edges. The above-mentioned logic of numbering advancement along the path is defined as the path advancement numbering logic S1, which is applicable to the numbering advancement of the first path and all subsequent paths.

[0129] After returning to the path starting point, judge whether the path starting point is a branch point. If not, directly return to the numbering starting point. At this time, the numbering starting point and the path starting point are the same graph vertex; if it is, it is necessary to execute the path advancement numbering along other sub-paths of the path starting point. The execution method uses the path advancement numbering logic S1 until all sub-paths of the path starting point have executed the path advancement numbering, and then return to the first branch point upstream of the path starting point. Take this branch point as the path starting point, execute the path advancement numbering for all sub-paths of this branch point, and then return to a more upstream branch point to perform the above path advancement numbering operation until there is no upstream branch point, and return to the numbering starting point. The above-mentioned numbering advancement logic that incorporates the influence of branch points is defined as the branch coverage logic S2. The logic S2 includes the path advancement numbering logic S1 and is applicable to the numbering advancement of the topological structure graph vertices with the numbering starting point as the upstream initial point.

[0130] Step 25: After returning to the starting point of numbering, select other loop inlet graph vertices as the next starting point for numbering the graph vertices and continue the numbering. The numbering advancement logic adopts the branch coverage logic S2 until all loop inlet graph vertices are traversed as the starting points for numbering. If all loop inlet graph vertices have been selected as the starting points for numbering, or the flattened topological connection general graph does not contain loop inlet graph vertices that have not been selected as the starting points for numbering, then traverse the graph vertices of the flattened topological connection general graph to determine whether all graph vertices have been numbered. If so, the numbering ends; if not, it means there are multiple loops, then execute the branch coverage logic S2, traverse the graph vertices of the flattened topological connection general graph again and execute the branch coverage logic S2 until all graph vertices have been numbered. The above numbering advancement logic that incorporates the influence of loop inlet graph vertices and the influence of multiple loops is defined as the inlet graph vertex and loop coverage logic S3. The logic S2 includes the path advancement numbering logic S1 and the branch coverage logic S2, and is applicable to the numbering advancement of the entire flattened topological structure connection general graph.

[0131] Number the graph vertices of the flattened topological connection general graph, i.e., the basic node units, through steps 21 to 25. One is to ensure the uniqueness and full coverage of the numbers of each graph vertex; the other is to ensure the orderliness and efficiency of the graph vertex traversal, so that the numbering advancement path is consistent with the fluid flow direction, which not only increases the readability of the graph vertex numbers and the diagonal dominance of the equations assembled based on the graph vertex numbers, but also quickly retrieves and traverses the graph vertices of the topological structure diagram according to the flow direction.

[0132] Step 3: According to the flattened topological structure and the basic node unit numbers corresponding to the thermal-hydraulic system model, process and calculate the data information of the control volume type basic node units and the pipe connection type basic node units to obtain the complete data information required for the assembly of the two-fluid six-equation coefficient matrix;

[0133] The assembly of the two-fluid six-equation coefficient components is carried out based on the discrete formats of the mass conservation sum equation, the mass conservation difference equation, the vapor-phase energy conservation equation, the liquid-phase energy conservation equation, the momentum conservation sum equation, and the momentum conservation difference equation. The discrete formats of each equation are as follows:

[0134] (1) Mass conservation sum equation

[0135]

[0136] (2) Mass conservation difference equation

[0137]

[0138] (3) Vapor-phase energy conservation equation

[0139]

[0140] (4) Liquid-phase energy conservation equation

[0141]

[0142] (5) Momentum conservation and equation

[0143]

[0144] (6) Momentum conservation difference equation

[0145]

[0146] Note: Indicates that the variable is a temporary result and needs to be further corrected after the equation is solved. Indicates that the upwind scheme is adopted, i.e.:

[0147] Among them, V is the control volume, A is the cross-sectional area of the control volume, P is the pressure, T is the temperature, U is the specific internal energy, H is the heat transfer coefficient, ρ is the density, α is the void fraction, v is the flow velocity, ∈ represents whether it is wall boiling or wall condensation, Δt is the time step, Δx is the length of the control volume, the subscript g represents the vapor phase, the subscript f represents the liquid phase, the subscript L represents the downstream control volume, the subscript K represents the upstream control volume, the subscript j represents the nozzle number, j - 1 represents the upstream nozzle, j + 1 represents the downstream nozzle, the superscript n represents the nth time moment, the superscript n + 1 is the (n + 1)th moment, the superscript s and the subscript s both represent the saturated state, the subscript w represents the wall, fx represents the intermediate variable for calculating the interfacial friction coefficient, Cp represents the specific heat at constant pressure, the subscript m represents the mass average of the vapor-liquid mixture, HLOSSG represents the local resistance term of the vapor phase at the nozzle, HLOSSF represents the local resistance term of the liquid phase at the nozzle, Q represents the heat transfer amount between phases or at the wall, DISS represents the energy dissipation term, Γ represents the interfacial mass exchange amount, FWG represents the vapor-phase wall friction term, FWF represents the liquid-phase wall friction term, VISG represents the vapor-phase artificial viscosity term, VISF represents the liquid-phase artificial viscosity term, FI represents the interfacial frictional force, C0 and C1 are the correction coefficients related to the interfacial friction, By is the body force in the y direction of the coordinate, and other symbols all belong to the source term relation expression.

[0148] After adopting semi-implicit discretization, the unknowns in equations (1) to (4) include those at the current (n + 1)th moment Among them, the first 5 variables do not belong to the solution variables of the two-fluid six equations and need to be expanded by Taylor series, i.e.:

[0149]

[0150] Substitute equations (7) to (11) into equations (1) to (4), and through arranging the variable coefficients and combining like terms, equations (1) to (4) can all be uniformly expressed in the following form:

[0151]

[0152] where a1 to a4 are respectively the variable coefficients in the equations, and s1 to s4 are respectively the variable coefficients in the equations, and b is the constant term of the equation. Among them, s1 to s4 are the convective transport coefficients of the mass or energy of the pipe connection, which are related to the control volume number, the pipe connection number and the equation number, and can be written as: s n，j，q，k , where n is the control volume number, j is the pipe connection number, q is the equation number, and k is 1 (vapor phase g) or 2 (liquid phase f).

[0153] In the momentum conservation equation (5) and the momentum conservation difference equation (6), the unknown variables only include those at the current n + 1 moment P L n+1 、P K n+1 . Through arranging the variable coefficients and combining like terms, equations (5) and (6) can both be uniformly expressed in the following form,

[0154]

[0155] where g1, g2, w1, w2 are respectively P L n+1 、P K n+1 the variable coefficients in the equations, and f is the constant term of the equation.

[0156] To provide the complete data information required for assembling the coefficient matrix of the two-fluid six-equation, combining equations (12) and (13), the content stored in the data container corresponding to the basic node element includes the model name where the basic node element is located, the basic node element type, the basic node element number, the upstream and downstream node information, the structural parameters, the characteristic parameters, the initial parameters and the generalized source term. The generalized source term refers 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, pipe connection variable calculation based on control volume variables, control volume variable calculation based on pipe connection variables, etc.

[0157] According to the flattened topological structure corresponding to the thermal-hydraulic system model and the basic node unit numbers, calculate and update the data information of the basic node units of the control volume type and the basic node units of the pipe connection type in real time, and obtain the complete data information required for the assembly of the two-fluid six-equation coefficient matrix. The specific calculation steps and content are as follows:

[0158] Step 31: Obtain the upstream and downstream pipe connection numbers of the basic node units of the control volume type according to the overall diagram of the flattened topological structure, so as to establish the constraint relationship between the solution variables such as pressure, void fraction, vapor specific internal energy, and liquid specific internal energy based on the control volume number and the solution variables such as vapor velocity and liquid velocity based on the pipe connection number in the mass difference conservation discrete equation and the energy conservation discrete equation;

[0159] Step 32: Obtain the upstream and downstream control volume numbers of the basic node units of the pipe connection type according to the overall diagram of the flattened topological structure, so as to establish the constraint relationship between the solution variables such as vapor velocity and liquid velocity based on the pipe connection number and the pressure variable based on the control volume number in the momentum difference conservation discrete equation;

[0160] Step 33: The vapor velocity and liquid velocity of the control volume required for the momentum conservation equation are obtained by averaging the vapor-liquid velocities of the upstream and downstream pipe connections and stored in the corresponding data container of the basic node unit of the control volume type;

[0161] The physical property parameters such as density, viscosity coefficient, thermal conductivity, and specific heat at constant pressure of the pipe connection, as well as the scalar quantities such as void fraction, vapor specific internal energy, and vapor specific internal energy, are all the physical property parameters and scalar quantities of the incoming flow control volume and are stored in the corresponding data container of the basic node unit of the pipe connection type. They are directly used in the calculation of the source term closure relation of the momentum conservation discrete equation and do not need to find the incoming flow control volume through the topological structure when assembling the matrix coefficient matrix;

[0162] Step 35: For the basic node units of the control volume type, calculate the generalized source terms of the mass and energy conservation discrete equations, and store the calculation results in the data container of the control volume. The calculation of the generalized source terms includes the calculation of physical properties, interphase mass transfer, interphase heat transfer, wall heat transfer, energy dissipation and other source terms. According to the upstream and downstream pipe connection numbers of the control volume, calculate the coefficients and constant terms of each solution variable of the mass and energy conservation discrete equations sorted by equation (12), including a1 to a4, s1 to s4, and b, which can be directly used when assembling the coefficient matrix. Among them, s1 to s4 need to correspond to the upstream and downstream pipe connection numbers of the control volume, the equation number, and the type of vapor-liquid velocity.

[0163] Step 36: For the basic node unit of the takeover type, calculate the generalized source term of the discrete equation of momentum conservation, and store the calculation result in the data container of the control volume. The calculation of the generalized source term includes the calculations of wall friction, interfacial friction, local resistance, artificial viscosity terms, etc., as well as the processing of special pumps, critical flows, etc. According to the control volume numbers of the upstream and downstream of the takeover, calculate the solution variable coefficients and constant terms of the momentum conservation sum equation and the momentum conservation difference equation sorted by Equation (13), including g1, g2, w1, w2, f, which can be directly used when assembling the coefficient matrix.

[0164] Step 37: A control volume may have multiple upstream and downstream takeovers (i.e., multi-way), but each takeover has only one upstream control volume and one downstream control volume. When dealing with the multi-way situation, since the number of upstream and downstream control volumes of the takeover is always 1, the number of variables and the structure of the discrete equation of momentum conservation remain unchanged, and Equation (13) can still be used for expression; and since the number of takeovers upstream and downstream of the control volume changes, the discrete equations of mass and energy conservation will increase the takeover flow velocity changes. For example, when the control volume has one more upstream takeover, the mass conservation sum equation of the control volume increases two terms and Sorted according to Equation (12), that is:

[0165]

[0166] It can be seen from Equation (14) that the scalar array, the variable coefficient array of the scalar, and the constant term on the left side do not change, while the vector array on the right side increases and The corresponding variable coefficient array of the right-side vector increases and This shows that when there are multiple upstream and downstream takeovers in the control volume, the scalar array, the variable coefficient array of the scalar, and the equation constant term on the left side of the equation remain unchanged, while the vector array on the right side is always included in the takeover flow velocity variable array of the thermal-hydraulic system. The corresponding variable coefficients all adopt the convective transport expressions of mass or energy, and the positive and negative signs of the variable coefficients are related to whether the takeover is upstream or downstream of the control volume. Therefore, the equation form of the discrete equations of mass and energy conservation remains unchanged and no special treatment is required.

[0167] Through data information processing and calculation based on basic node units, first, complete data information required for assembling the two-fluid six-equation coefficient matrix on basic node units is obtained and stored in the corresponding data containers of the basic node units. During the coefficient matrix assembly, it is directly read according to the basic node unit numbers, avoiding repeated loading of topological structure information and obtaining data information of upstream and downstream basic node units based on topological structure information during the coefficient matrix assembly. Second, it supports establishing a parallel computing mode for two-fluid six-equation data information based on basic node units. Since the data information calculation of basic node units is only related to the basic node unit numbers, parallel computing of data information of multiple basic node units is supported, effectively improving the efficiency of coefficient matrix assembly.

[0168] Step 4: Based on the basic node units of the control volume type, according to the topological structure and the set equation sorting and variable sorting, read the data information of the basic node units, and carry out the assembly of the coefficient matrix of the mass and energy conservation discrete equations. The coefficient matrix consists of the coefficient matrix of the solution scalars such as the pressure, void fraction, vapor specific internal energy, and liquid specific internal energy of the control volume, the coefficient matrix of the solution vectors such as the vapor velocity and liquid velocity of the upstream and downstream nozzles of the control volume, and the constant term. The solution scalar coefficient matrix and the solution vector coefficient matrix are sparse matrices and are stored in row-compressed format.

[0169] The coefficient matrix of the two-fluid six equations is obtained by organizing according to equations (12) and (13). The j-th element a of the i-th row of the matrix ij corresponds to the coefficient of the j-th variable in the i-th equation. The j-th variable refers to the variable ranked in the j-th row of the variable array. Therefore, the final form of the matrix is related to both the arrangement order of the equations and the arrangement order of the variables. Since the discrete conservation equations are established according to the control volume and the nozzles, the arrangement order of the equations and variables can be set according to the numbering order of the control volume and the nozzles.

[0170] The row elements of the coefficient matrix correspond to the coefficients of all solution variables. When a variable is not in this equation, its coefficient is 0. According to equations (12) and (13), the mass and energy conservation discrete equations of the control volume only contain the variables and parameters of this control volume and its upstream and downstream nozzles, while the momentum conservation discrete control equation of the nozzle only contains the variables and parameters of this nozzle and its upstream and downstream control volumes. The coefficients of the variables of the nodes not adjacent to this control volume or nozzle in the matrix row elements must be 0. Therefore, there will be a large number of 0 elements in the matrix. When the number of 0 elements in a matrix is much larger than the number of non-0 elements, it is also called a sparse matrix. When assembling a large sparse matrix, only the non-0 elements in the matrix need to be considered, and the row-compressed storage method of the sparse matrix is used to store only the non-zero elements in the sparse matrix, thereby greatly reducing the memory occupation.

[0171] The present invention stores large sparse matrices in the Compressed Sparse Row (CSR) format. The method of the row compression format is briefly introduced below.

[0172] Taking a matrix in the following form as an example,

[0173]

[0174] It can be stored in the following way. First, load all non-zero elements into a one-dimensional array in the order of rows first and columns second as the constitutive data,

[0175] data = [a, b, c, d, e] (16)

[0176] data=[a,b,c,d,e] (16)

[0177] Then record the length of data,

[0178] nnz = 5 (17)

[0179] Next, record the column numbers where all the data in data are located in sequence,

[0180] colIndex = [1, 2, 2, 2, 3] (18)

[0181] Finally, record the serial numbers of the first non-zero elements in each row in data, and supplement the number of all non-zero elements plus 1 at the end,

[0182] rowIndex = [1, 3, 4, 6] (19)

[0183] nnz, collndex, rowIndex and the number of rows N of the matrix are called the structural parameters of the matrix. Combining with the data array, the complete information representing the large sparse matrix can be obtained.

[0184] In the thermal-hydraulic system model, assuming there are N non-boundary control volumes and M non-boundary nozzles, then there are 4N control volume variables (pressure, void fraction, specific internal energy of the vapor phase, specific internal energy of the liquid phase) and 2M nozzle variables (vapor phase flow rate, liquid phase flow rate). The mass conservation equations, mass conservation difference equations, vapor phase energy conservation equations, and liquid phase energy conservation equations of the N control volumes are sorted according to Equation (12) and combined into the matrix form:

[0185] AX = SV + B (20)

[0186] Among them,

[0187]

[0188] A is a 4N×4N coefficient matrix for solving scalar quantities such as control volume pressure, void fraction, vapor specific internal energy, and liquid specific internal energy. X is an array of 4N control volume variables. S is a 4N×2M coefficient matrix for solving vectors such as vapor velocity and liquid velocity at the upstream and downstream nozzles of the control volume. V is an array of 2M nozzle variables. B is an array of 4N constant terms.

[0189] The assembly steps of the coefficient matrix of the mass and energy conservation discrete equations are as follows:

[0190] Step 41: Obtain the number of non-boundary control volumes N, the number of non-boundary nozzles M, the number of flow boundary nozzles M_TDJ, and the number of multi-pass extra nozzles M_ExtraConnectJun according to the topological structure information.

[0191] Step 42: Calculate the number of non-zero elements of matrix A as 16N, and the number of non-zero elements of matrix S as 12N + 6(M_ExtraConnectJun - M_TDJ).

[0192] Step 43: Apply for the memory of the row-compressed format arrays of matrix A and matrix S according to the number of non-zero elements of matrix A and matrix S, including: A_data, A_colindex, A_rowindex, and S_data, S_colindex, S_rowindex, and apply for the memory of the constant term value B.

[0193] Step 44: Set the mass and energy conservation equations and variable sorting as:

[0194]

[0195] Step 45: In the range of non-boundary control volumes, select control volumes in sequence according to the number from i = 1 (i is a natural number from 1 to N), obtain the data container of the i-th control volume, read the structure parameters, characteristic parameters, initial parameters, state parameters, physical property parameters, generalized source term calculation results stored in the data container, and the coefficients and constants of each solution variable of the mass and energy conservation discrete equations sorted according to equation (12), including variable coefficients a1 to a4, s1 to s4, b. The variable coefficients and constants are different in different conservation discrete equations, and s1 to s4 need to correspond to the control volume number, nozzle number, and equation number. s1 to s4 can all be rewritten as: s n，j，q，k , where n is the control volume number, j is the nozzle number, q is the equation number, and k is 1 (vapor phase g) or 2 (liquid phase f). If there is a multi-pass, the number of variable coefficients of the solution vector will increase. Assuming that m upstream nozzles or downstream nozzles are added, then the variable coefficients of the solution vector of each type of equation will increase from s1 to s4 to s1 to s 4+2m .

[0196] Step 46: Calculate the number of rows and columns of the variable coefficient of the scalar in matrix A according to the set equation and variable sorting. Taking the variable coefficient of void fraction in the mass conservation difference equation as an example, the equation sorting of the mass conservation difference equation of the i-th control volume is 4i - 2, and the sorting of the void fraction of the i-th control volume in the variable array X is N + i. Therefore, the number of rows of the variable coefficient of void fraction in the mass conservation difference equation in matrix A is 4i - 2, and the number of columns is N + i, and its value is the variable coefficient a2 corresponding to the void fraction of the mass conservation difference equation stored in the data container of the i-th control volume.

[0197] Step 47: Store the variable coefficient in the row-compressed format array of matrix A according to the number of rows and columns of the variable coefficient of the solved scalar in matrix A.

[0198] Step 48: Obtain the upstream and downstream connection pipe information of the i-th control volume according to the topological structure, clarify the numbers and data information of each upstream and downstream connection pipe. Assume the number of the connection pipe is j, and confirm whether the connection pipe numbered j is a flow boundary connection pipe, which is divided into the following situations:

[0199] Case 1: If the connection pipe numbered j is a flow boundary connection pipe, then read the physical property parameters, structural parameters, vapor-liquid flow velocity and other data information of the data container of the connection pipe numbered j, as well as the variable coefficients of the vapor-liquid flow velocity of the connection pipe numbered j in the mass conservation and equation, mass conservation difference equation, vapor-phase energy conservation equation and liquid-phase mass conservation equation in the data container of the i-th control volume, calculate the expression term of the vapor-liquid flow velocity of the connection pipe numbered j in the equation. This expression term is equal to the product of the vapor-liquid flow velocity of the connection pipe and its corresponding variable coefficient in the equation, and add it to the constant term of the corresponding equation of the i-th control volume, update the constant term of the data container, and jump to Step 410.

[0200] Case 2: If the connection pipe numbered j is a non-boundary connection pipe, then calculate the number of rows and columns of the variable coefficient of the solved vector in matrix S according to the equation sorting of formula (21) and the variable sorting of formula (22). Taking the mass conservation difference equation as an example, the number of rows of the variable coefficient corresponding to the vapor-phase flow velocity of the connection pipe numbered j in the mass conservation difference equation of the i-th control volume in matrix S is 4i - 2, and the number of columns is j - m_tdj (where m_tdj is the number of flow boundary connection pipes among the connection pipes numbered 1 to j), and its value is the variable coefficient s of the vapor-phase flow velocity of the connection pipe numbered j stored in the data container of the i-th control volume in the mass conservation difference equation n，j，2，1 , where: n is the number of the i-th control volume, and j is the number of the connection pipe.

[0201] Step 49: Store the variable coefficient in the row-compressed format array of matrix S according to the number of rows and columns of the variable coefficient of the solved vector in matrix S.

[0202] Step 410: Sort according to the mass and energy conservation equations, and store the constant terms b of the mass conservation and equation, the mass conservation difference equation, the vapor-phase energy conservation equation, and the liquid-phase mass conservation equation of the i-th control volume data container into the corresponding positions of the constant term array B.

[0203] Step 411: Repeat steps 45 to 410 to sequentially complete the storage of the variable coefficients and constant terms for the control volume numbers i from 1 to N, and then the A matrix, S matrix, and constant term array B can be assembled.

[0204] For example, for the following thermal-hydraulic system model, it includes two non-boundary control volumes V1 and V2, two non-boundary nozzles, and one flow boundary nozzle.

[0205] According to the assembly steps of the coefficient matrix of the mass and energy conservation discrete equations, the arrangement order for assembling the coefficient matrix AX = SV + B is:

[0206]

[0207] The arrangement order of the control volume variable array on the left side of the equation is:

[0208]

[0209] The arrangement order of the flow velocity variable array on the right side of the equation is:

[0210]

[0211] According to the assembly steps of the coefficient matrix of the mass and energy conservation discrete equations, it can be known that Figure 12 In the coefficient matrix of the mass and energy conservation equations of the thermal-hydraulic system model, the number of non-zero elements in the A matrix is 16×2 = 32, and the number of non-zero elements in the S matrix is 12×2 + 6(0 - 1) = 18.

[0212] When assembling the A matrix, select the control volume V1, and the solution scalars of the mass and energy conservation discrete equations are P1, α1, U g，1 、U f，1 , and the variable coefficients corresponding to the solution scalars on the mass conservation and equation, the mass conservation difference equation, the vapor-phase energy conservation equation, and the liquid-phase energy conservation equation are expressed in matrix form, specifically:

[0213]

[0214] Among them, a 111 to a 114 are the variable coefficients a1 to a4 of the solution variables regarding the mass conservation and equation stored in the control volume V1 data container, and a 121 to a 124The variable coefficients a1 to a4, a for the solution variables of the mass conservation difference equation stored in the V1 data container 131 to a 134 The variable coefficients a1 to a4, a for the solution variables of the vapor-phase energy conservation equation stored in the V1 data container 141 to a 144 The variable coefficients a1 to a4 for the solution variables of the liquid-phase energy conservation equation stored in the V1 data container.

[0215] According to the equation sorting of equation (23) and the variable sorting of equation (24), a 111 to a 114 are respectively the 1st, 3rd, 5th, and 7th columns of the first row of the A matrix, a 121 to a 124 are respectively the 1st, 3rd, 5th, and 7th columns of the second row of the A matrix, a 131 to a 134 are respectively the 1st, 3rd, 5th, and 7th columns of the third row of the A matrix, a 141 to a 144 are respectively the 1st, 3rd, 5th, and 7th columns of the fourth row of the A matrix.

[0216] Similarly, the variable coefficients corresponding to the scalar solutions of the mass and energy conservation discrete equations stored in the control volume V2 data container are expressed in matrix form as follows:

[0217]

[0218] According to the equation sorting of equation (21) and the variable sorting of equation (22), a 211 to a 214 are respectively the 2nd, 4th, 6th, and 8th columns of the fifth row of the A matrix, a 221 to a 224 are respectively the 2nd, 4th, 6th, and 8th columns of the sixth row of the A matrix, a 231 to a 234 are respectively the 2nd, 4th, 6th, and 8th columns of the seventh row of the A matrix, a 241 to a 244 are respectively the 2nd, 4th, 6th, and 8th columns of the eighth row of the A matrix.

[0219] Therefore, the assembled A matrix is:

[0220]

[0221] In row-compressed format, the A matrix is stored as follows:

[0222] A_data = [a 111 ,a112 , a 113 , a 114 , a 121 , a 122 , a 123 , a 124 , … a 241 , a 242 , a 243 , a 244 (29)

[0223] A_nnz = 32 (30)

[0224] A_colindex = [1, 3, 5, 7, 1, 3, 5, 7, …, 2, 4, 6, 8] (31)

[0225] A_rowindex = [1, 5, 9, 13, 17, 21, 25, 29, 33] (32)

[0226] Assemble the S matrix and select the V1 control volume. Since there is only an outlet with a non-boundary nozzle in the V1 control volume, the solution vector for the mass and energy conservation discrete equations of the V1 control volume is v g,1 , v f，1 . The variable coefficients corresponding to the solution vector in the mass conservation and equation, mass conservation difference equation, vapor phase energy conservation equation, and liquid phase energy conservation equation are expressed in matrix form, specifically as follows:

[0227]

[0228] Among them, s 1111 and s 1112 are the variable coefficients of the vapor velocity v g,1 and the liquid velocity v f,1 stored in the data container of the control volume V1 for the mass conservation and equation, respectively. s 1121 and s 1122 are the variable coefficients of the vapor velocity v g，1 and the liquid velocity v f，1 stored in the data container of the control volume V1 for the mass conservation difference equation, respectively. s 1131 is the variable coefficient of the vapor velocity v g，1 stored in the data container of the control volume V1 for the vapor conservation equation. s 1142 is the variable coefficient of the liquid velocity v f，1 stored in the data container of the control volume V1 for the liquid conservation equation.

[0229] According to the equation sorting in equation (23) and the variable sorting in equation (25), s 1111 and s 1112They are the 1st column and the 3rd column of the first row of the S matrix, s 1121 and s 1122 They are the 1st column and the 3rd column of the second row of the S matrix, s 1131 is the 1st column of the third row of the S matrix, s 1142 s 1131 is the 3rd column of the fourth row of the S matrix.

[0230] Select the V2 control volume. There is one non-boundary nozzle at both the inlet and outlet of the V2 control volume. The solution vectors include v g，1 、v g,2 、v f,1 、v f,2 . The variable coefficients corresponding to the solution vectors in the mass conservation equation, the mass conservation difference equation, the vapor phase energy conservation equation, and the liquid phase energy conservation equation are expressed in matrix form, specifically as follows:

[0231]

[0232] According to the equation sorting in equation (23) and the variable sorting in equation (25), s 2111 、s 2211 、s 2112 and s 2212 are the 1st column, the 2nd column, the 3rd column, and the 4th column of the fifth row of the S matrix respectively, s 2121 、s 2221 、s 2122 and s 2222 are the 1st column, the 2nd column, the 3rd column, and the 4th column of the sixth row of the S matrix respectively, s 2131 and s 2231 are the 1st column and the 2nd column of the seventh row of the S matrix respectively, s 2142 and s 2242 are the 3rd column and the 4th column of the eighth row of the S matrix respectively. Therefore, the assembled S matrix is:

[0233]

[0234] Adopting the row compression format, the S matrix is stored as follows:

[0235] S_data = [s 1111 , s 1112 , s 1121 , s 1122 , s 1131 , s 1142 , … s 2142 , s 2242 (36)

[0236] S_nnz = 18 (37)

[0237] S_colindex = [1, 3, 1, 3, 1, 3, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4] (38)

[0238] S_rowindex = [1, 3, 5, 6, 7, 11, 15, 17, 19] (39)

[0239] The array form of the constant term B is:

[0240] B = [b 11 , b 12 , b 13 , b 14 , b 21 , b 22 , b 23 , b 24 (40)

[0241] Among them, the first subscript of b 11 represents the control volume number, and the second subscript represents the equation number. In addition, b 11 to b 14 consider the influence of the flow boundary connection on the equation constant term.

[0242] After completing the assembly of the coefficient matrix, the matrix form of the AX = SV + B of the mass and energy conservation equations is:

[0243]

[0244] Step 5: Based on the basic node elements of the connection type, according to the topological structure and the set equation sorting and variable sorting, read the data information of the basic node elements, and carry out the assembly of the coefficient matrix of the momentum conservation discrete equation. The coefficient matrix consists of the coefficient matrices of the solution vectors such as the vapor phase flow rate and liquid phase flow rate of the connection, the coefficient matrix of the control volume pressures upstream and downstream of the connection, and the constant term. The solution vector coefficient matrix and the pressure coefficient matrix are sparse matrices and are stored in row compression format.

[0245] The assembly of the coefficient matrix of the momentum conservation discrete equation includes the following:

[0246] In the thermal-hydraulic system model, assume that there are N non-boundary control volumes and M non-boundary connections. Then there are 4N control volume variables (pressure, void fraction, vapor specific internal energy, liquid specific internal energy) and 2M connection variables (vapor phase flow rate, liquid phase flow rate). Rearrange the momentum conservation equations and momentum conservation difference equations of the M connections according to Equation (13) and combine them into matrix form:

[0247] GV = WP + F (42)

[0248] Among them,

[0249] G is a 2M×2M matrix, V is an array of 2M variables of the vapor and liquid flow rates at the inlets and outlets, W is a 2M×N matrix, and P is an array of pressure variables for N control volumes.

[0250] The steps for assembling the coefficient matrix of the discrete equation for momentum conservation are as follows:

[0251] Step 51: Obtain the number N of non-boundary control volumes, the number M of non-boundary inlets and outlets, and the number N_P of pressure boundaries TDV according to the topological structure information.

[0252] Step 52: Calculate the number of non-zero elements 4M of the G matrix and the number of non-zero elements 4M - 2N_P of the W matrix.

[0253] Step 53: Apply for the memory of the row-compressed format arrays of the G matrix and the W matrix according to the number of non-zero elements of the G matrix and the W matrix, including: G_data, G_colindex, G_rowindex and W_data, W_colindex, W_rowindex, and apply for the memory of the constant term value F.

[0254] Step 54: Set the momentum conservation equation and variable sorting as:

[0255]

[0256] Step 55: Among the number M of non-boundary inlets and outlets, sequentially select the inlets and outlets with j = 1 (j is a natural number from 1 to M) according to the inlet and outlet numbers, obtain the data container of the j-th inlet and outlet, and read the structural parameters, characteristic parameters, initial parameters, state parameters, physical property parameters, calculation results of the generalized source term, and the variable coefficients and constant terms of each solution variable of the momentum conservation equation sorted according to Equation (13) stored in the data container, including g1, g2, w1, w2, f. The variable coefficients and constants are different in the momentum conservation equation and the momentum conservation difference equation. It should be noted that since each inlet and outlet has exactly one upstream control volume and one downstream control volume, in the case of multi-pass, the momentum conservation equation is still expressed in the equation form of Equation (13) without adding additional variable coefficients.

[0257] Step 56: According to the sorting of the momentum conservation equation in Equation (43) and the variable sorting in Equation (14), calculate the number of rows and columns of the variable coefficients of the solution vector in the G matrix. Taking the variable coefficient of the liquid flow rate in the momentum conservation equation as an example, the sorting of the momentum conservation equation of the j-th inlet and outlet is 2j - 1 in the equation sorting, and the sorting of the liquid flow rate of the j-th inlet and outlet in the variable array V is M + j. Therefore, the number of rows of the variable coefficient of the liquid flow rate in the momentum conservation equation in the G matrix is 2j - 1, and the number of columns is M + j, and its value is the variable coefficient g2 corresponding to the liquid flow rate of the momentum conservation equation stored in the data container of the j-th inlet and outlet.

[0258] Step 57: Store the variable coefficients in the row-compressed format array of the G matrix according to the number of rows and columns of the G matrix corresponding to the variable coefficients of the solution vector.

[0259] Step 58: According to the topological structure, obtain the upstream control volume and downstream control volume of the j-th takeover, clarify the numbers and data information of the upstream control volume and downstream control volume. Assume the number of the upstream control volume or downstream control volume is n, and confirm whether the control volume numbered n is a pressure boundary control volume. It is divided into the following situations:

[0260] Situation 1: Assume that the control volume numbered n is a pressure boundary control volume. Then read the physical property parameters, structural parameters, pressure and other data information of the data container of the control volume numbered n, as well as the variable coefficients of the pressure of the control volume numbered n in the momentum conservation sum equation and momentum conservation difference equation in the data container of the j-th takeover. Calculate the expression term for the pressure of the control volume numbered n. This expression term is equal to the control volume pressure multiplied by its corresponding variable coefficient in the equation, and add it to the constant term of the corresponding equation of the j-th takeover. Update the constant term of the data container, and jump to Step 510.

[0261] Situation 2: Assume that the control volume numbered n is not a pressure boundary control volume. Then, according to the sorting of the momentum conservation equation in Equation (43) and the variable sorting in Equation (44), calculate the number of rows and columns of the variable coefficients of the control volume pressure in the W matrix. Taking the momentum conservation sum equation as an example, the number of rows of the variable coefficients of the momentum conservation sum equation of the j-th takeover for the pressure of the control volume numbered n in the W matrix is 2j - 1, and the number of columns is n - n_P (n_P is the number of pressure boundary control volumes among the control volumes numbered from 1 to n). If the control volume numbered n is the upstream control volume of the j-th takeover, then its variable coefficient is w1, which is the variable coefficient of the upstream control volume pressure in the momentum conservation sum equation stored in the data container of the j-th takeover.

[0262] Step 59: Store the variable coefficients in the row-compressed format array of the W matrix according to the number of rows and columns of the W matrix corresponding to the variable coefficients of the solution vector.

[0263] Step 510: Store the constant terms of the momentum conservation sum equation and momentum conservation difference equation of the data container of the j-th takeover into the corresponding positions of the constant term array F according to the sorting of the momentum conservation equation in Equation (43).

[0264] Step 511: Repeat Step 55 to Step 510, and successively complete the storage of the variable coefficients and constant terms for the control volume numbers j from 1 to M, then the G matrix, W matrix, and constant term array F can be assembled and completed.

[0265] Take Figure 12Taking the thermal-hydraulic system model as an example, it includes 2 non-boundary control volumes, 2 non-boundary connection numbers M, and 1 pressure boundary TDV.

[0266] According to the assembly steps of the coefficient matrix of the momentum conservation discrete equation, the arrangement order of the coefficient matrix GV = WP + F is:

[0267]

[0268] The arrangement order of the control volume variable array on the left side of the equation is:

[0269]

[0270] The arrangement order of the flow velocity variable array on the right side of the equation is:

[0271]

[0272] According to the assembly steps of the coefficient matrix of the momentum conservation discrete equation, it can be known that Figure 12 In the coefficient matrix of the momentum conservation equation of the thermal-hydraulic system model, the number of non-zero elements in the G matrix is 4×2 = 8, and the number of non-zero elements in the W matrix is 4×2 - 2×1 = 6.

[0273] When assembling the G matrix and selecting the J1 connection, the solution vector of the momentum conservation equation is v g，1 、v f，1 . The variable coefficients corresponding to the solution vector in the momentum conservation and equation, and the momentum conservation difference equation are expressed in matrix form, specifically:

[0274]

[0275] Among them, g 111 、g 112 are the variable coefficients g1 and g2 of the solution variables v g，1 、v f，1 stored in the data container of connection J1 regarding the momentum conservation and equation respectively. g 121 、g 122 are the variable coefficients g1 and g2 of the solution variables stored in the data container of connection J1 regarding the momentum conservation difference equation respectively.

[0276] According to the equation sorting of formula (45) and the variable sorting of formula (46), g 111 、g 112 are the 1st column and the 3rd column of the first row of the G matrix respectively. g 121 、g 122 are the 1st column and the 3rd column of the second row of the A matrix respectively.

[0277] Similarly, the variable coefficients corresponding to the vector for solving the discrete equation of conservation of momentum stored in the takeover J2 data container are expressed in matrix form as follows:

[0278]

[0279] According to the equation sorting in formula (45) and the variable sorting in formula (46), g 111 and g 112 are the 2nd and 4th columns of the third row of the G matrix respectively, and g 121 and g 122 are the 2nd and 4th columns of the fourth row of the G matrix respectively.

[0280] Therefore, the assembled G matrix is:

[0281]

[0282] In row-compressed format, the G matrix is stored as follows:

[0283] G_data = [g 111 , g 112 , g 121 , g 122 , g 211 , g 212 , g 221 , g 222 (51)

[0284] G_nnz = 8 (52)

[0285] G_colindex = [1, 3, 1, 3, 2, 4, 2, 4] (53)

[0286] G_rowindex = [1, 3, 5, 7, 9] (54)

[0287] Assemble the W matrix, select the takeover J1. The right side of the discrete equation of conservation of momentum for the takeover J1 solves for the pressure variable including the upstream control volume pressure P1 and the downstream control volume pressure P2. The variable coefficients corresponding to the solved pressure variable on the conservation of momentum equation and the difference equation of conservation of momentum are expressed in matrix form, specifically:

[0288]

[0289] Among them, w 111 and w 112 are the upstream control volume pressure variable coefficient w1 and the downstream control volume pressure variable coefficient w2 for the conservation of mass equation stored in the takeover J1 data container respectively, and w 121 and w 122They are the upstream control volume pressure variable coefficient w1 and the downstream control volume pressure variable coefficient w2 for the mass conservation difference equation stored in the takeover J1 data container respectively.

[0290] According to the equation sorting of equation (45) and the variable sorting of equation (47), w 111 , w 112 are the 1st column and the 2nd column of the first row of the W matrix respectively, and w 121 , w 122 are the 1st column and the 2nd column of the second row of the W matrix respectively.

[0291] Select the J2 takeover. Since the downstream of the J2 takeover is a pressure boundary control volume, the pressure variable solved on the right side of the momentum conservation discrete equation of the J2 takeover only includes the upstream control volume pressure P2. The variable coefficients corresponding to the solved pressure variable on the momentum conservation and equation and the momentum conservation difference equation are expressed in array form, specifically:

[0292]

[0293] Among them, w 211 is the upstream control volume pressure variable coefficient w1 for the mass conservation and equation stored in the takeover J2 data container, and w 221 is the upstream control volume pressure variable coefficient w1 for the mass conservation difference equation stored in the takeover J2 data container.

[0294] According to the equation sorting of equation (45) and the variable sorting of equation (47), w 211 is the 2nd column of the third row of the W matrix, and w 221 is the 2nd column of the third row of the W matrix. Therefore, the assembled W matrix is:

[0295]

[0296] In row compression format, the W matrix is stored as follows:

[0297] W_data = [w 111 , w 112 , w 121 , w 122 , w 211 , w 221 (58)

[0298] S_nnz = 6 (59)

[0299] W_colindex = [1, 2, 1, 2, 2, 2] (60)

[0300] W_rowindex = [1, 3, 5, 6, 7] (61)

[0301] The array form of the constant term F is as follows:

[0302] F = [f 11 , f 12 , f 21 , f 22 (62)

[0303] where the first subscript of f 11 represents the control volume number, and the second subscript represents the equation number. Additionally, f 21 , f 22 consider the influence of the pressure boundary on the equation constant term.

[0304] After completing the assembly of the coefficient matrix, the matrix form of the momentum conservation equation GV = WP + F is as follows:

[0305]

Claims

1. A method for assembling and storing the coefficient matrix of two-fluid six-equation applicable to Modelica, characterized in that The steps include: Step 1: Automatically generate the topology of the thermal-hydraulic system, and restore the thermal-hydraulic system model with hierarchical encapsulation characteristics to a flattened topological connection diagram consisting entirely of orderly connected basic node unit models; Step 2: Traverse the basic node units of the flattened topological structure connection diagram corresponding to the thermal-hydraulic system model, and uniformly number the control body type basic node units and the pipe type basic node units according to the direction, branches, and loops of the orderly connection of the topological structure; Step 3: Based on the flattened topology and basic node unit numbers corresponding to the thermal-hydraulic system model, the data information of the control body basic node units and the nozzle basic node units are processed and calculated to obtain the complete data information required for assembling the coefficient matrix of the two-fluid six equations; Step 4: Based on the basic node units of the control body, according to the topological structure and the set equation and variable sorting, read the data information of the basic node units, and assemble the coefficient matrix of the discrete equations of mass and energy conservation. The coefficient matrix consists of the coefficient matrix of scalar solutions such as the pressure of the control body, the cavity fraction, the vapor phase internal energy, and the liquid phase internal energy; the coefficient matrix of vector solutions such as the vapor phase flow rate and liquid phase flow rate of the upstream and downstream pipes of the control body; and the constant term. The scalar coefficient matrix and the vector coefficient matrix are sparse matrices and are stored in a row-compressed format. Step 5: Based on the basic node units of the takeover type, according to the topological structure and the set equation sorting and variable sorting, read the data information of the basic node units, and assemble the coefficient matrix of the momentum conservation discrete equation. The coefficient matrix is composed of the coefficient matrix of the solution vectors such as the vapor phase flow rate and liquid phase flow rate of the takeover, the coefficient matrix of the upstream and downstream control body pressures of the takeover, and the constant term. The solution vector coefficient matrix and pressure coefficient matrix are sparse matrices and are stored in row compression format.

2. The method for assembling and storing the coefficient matrix of two-fluid six-equation applicable to Modelica according to claim 1, wherein The method for automatically generating the topological structure of the thermal hydraulic system model in step 1 comprises the following steps: Step 11: Identify the Modelica-based thermal-hydraulic model through model annotation keywords. The model annotation keywords are divided into basic node unit model identification keywords and thermal-hydraulic component model identification keywords. The model annotation keywords belong to the expanded Modelia syntax semantics, that is, by introducing the FMU model framework to support thermal-hydraulic hierarchical modeling based on basic node units; the basic node unit model includes the control body model v, the takeover model j, the time-dependent control body model tdv, and the time-dependent takeover model tdj. The thermal-hydraulic component model refers to the thermal-hydraulic model of pipes, containers and other internally defined discrete nodes, as well as the higher-level thermal-hydraulic model formed by connecting and encapsulating the thermal-hydraulic component models; Step 12: Flatten the thermal hydraulic component model, converting it into a component model consisting of control bodies, nozzles, flow interfaces, and directed connection edges according to the number of discrete nodes and hierarchical packaging information; For the pipeline model, according to the defined number of discrete nodes N of the pipeline model, the pipeline model is converted into a component model consisting of N control volumes, N - 1 takeovers, 4N - 2 basic node unit flow interfaces, 2 component model flow interfaces, and 2N directed connection edges connected by connectors; For the pipeline model and the thermal - hydraulic model after encapsulating the pipeline model, the model is decomposed into a pipeline model along the hierarchical encapsulation information, and then the pipeline is decomposed into a model composed of control volumes, takeovers, component interfaces, and directed connection edges; Step 13: According to the connection relationship of the model interfaces in the thermal - hydraulic system modeling process, the compiler connects the basic node unit model and the flattened thermal - hydraulic component model to form a general model connection diagram. The vertices of the diagram include the basic node unit model and component interfaces. The component interfaces are the inlet flow interfaces and outlet flow interfaces of the thermal - hydraulic component model. The vertices of the diagram are stored using the data container NumapNode. Each vertex corresponds to a data container. The data container NumapNode corresponding to the basic node unit model stores the data information required for solving the two - fluid six - equation, including the model name where the vertex is located, structural parameters, characteristic parameters, initial parameters, state parameters, upstream and downstream node information, and generalized source terms. The data container NumapNode corresponding to the component interface only stores the thermal - hydraulic component model information, including the component name and interface type, and does not store the data information required for solving the two - fluid six - equation; Step 14: Traverse the vertices of the general model connection diagram to obtain the type and information of the vertex NumapNode, as well as the number and information of the upstream and downstream vertices. According to the connection rules of the thermal - fluid model interface, when the vertex is a component interface, check the number of upstream and downstream vertices to confirm whether both are 1. If not, generate an error message and jump out of the topological structure generation logic; if so, delete the vertex and the directed connection edges connecting it to the upstream and downstream vertices, and then generate a directed connection edge from the upstream vertex K_Node to the downstream vertex L_Node; Step 15: Traverse the directed connection edges of the flattened topological general diagram. If the NumapNode types of the two graph vertices connected by the directed connection edge are both takeover - type graph vertices, including takeover j or time - related takeover tdj, then output an error message; The flattened topology graph is traversed through the takeover graph vertices. If the number of directed connection edges of the takeover graph vertex import or export is too large, an error message is output. If the two graph vertices K_Node and L_Node connected by the directed connection edge are both of the control body v or the time-related control body tdv, the directed connection edge between the graph vertices K_Node and L_Node is deleted, and a graph vertex of the type takeover j is added between the two graph vertices. Then, the graph vertices K_Node and L_Node and the directed connection edges between the newly added graph vertices are established based on the upstream and downstream relationships. The flow area and hydraulic equivalent diameter of the added graph vertex of the type takeover j are equal to the flow area and hydraulic equivalent diameter of the graph vertex with the smallest flow area between the graph vertices K_Node and L_Node. The length of takeover j is equal to half of the sum of the lengths of the graph vertices K_Node and L_Node. By traversing the directed connection edges of the flattened topology graph, a flattened topology graph in the form of "VJV" is formed.

3. A method for assembling and storing the coefficient matrix of two-fluid six equations applicable to Modelica according to claim 1, characterized in that, The specific steps of numbering the basic node units in step 2 are as follows: Step 21: Traverse the graph vertices of the flattened topology connection graph, obtain the loop import graph vertices and their number, and select a loop import graph vertex as the first numbering starting point; if there is no loop import graph vertex, randomly select a control body class graph vertex as the first numbering starting point, and the first numbering starting point is recorded as V1. The loop import graph vertex refers to the control body class graph vertex with only outgoing directed connection edges. The numbering starting point refers to the initial starting point when the numbering is advanced along the topology connection graph, that is, there is no upstream graph vertex in the numbering advancement route; Step 22: Using the first numbered starting point as the path starting point, determine the first path for numbered advancement. The first junction class graph vertex of the path is denoted as J1. If there is a branch downstream of the first numbered starting point, select a subpath of the numbered starting point as the first path for numbered advancement. If there is no branch downstream of the first numbered starting point, select the path where the outlet junction of the numbered starting point is located as the first path for numbered advancement. The path refers to a single line without branches in the topological structure connection diagram that uses a VJV directed connection, expressing the flow route of the fluid from upstream to downstream in a single pipeline in the thermal-hydraulic system. Step 23: Advance along the VJV directed connection of the path, sequentially number the control body class graph vertices and the takeover class graph vertices respectively until the end of the path, and return to the starting point of the path. The end of the path is the control body class graph vertex, including branch points, numbered graph vertices and outlet graph vertices. The branch point refers to the control body class graph vertex with two or more outlet connection takeovers, and the outlet graph vertex is the control body class graph vertex with only inlet directed connection edges. The above logic of numbering and advancing along the path is defined as the path advancement numbering logic S1, which is applicable to the numbering advancement of the first path and all subsequent paths; Step 24: After returning to the path starting point, determine whether the path starting point is a branch point. If not, directly return to the numbered starting point. At this time, the numbered starting point and the path starting point are the same graph vertex; If so, it is necessary to perform path advancement numbering along other subpaths of the path starting point. The execution method adopts the path advancement numbering logic S1 until all subpaths of the path starting point have performed path advancement numbering. Then, it returns to the first branch point upstream of the path starting point, takes this branch point as the path starting point, and performs path advancement numbering on all subpaths of this branch point. Then, it returns to the branch point further upstream and performs the above-mentioned path advancement numbering operation until there are no upstream branch points. Then, it returns to the numbering starting point. The above-mentioned numbering advancement logic incorporating the influence of branch points is defined as branch coverage logic S2. Logic S2 includes the path advancement numbering logic S1 and is applicable to the vertex numbering advancement of the topological structure graph with the numbering starting point as the upstream initial point. Step 25: After returning to the numbering starting point, select other loop import graph vertices as the next numbering starting point to continue numbering the graph vertices. The numbering advancement logic adopts branch covering logic S2 until all loop import graph vertices are traversed as numbering starting points. If all loop import graph vertices have been selected as numbering starting points, or the flattened topology connection graph does not contain loop import graph vertices that have not been selected as numbering starting points, then the graph vertices of the flattened topology connection graph are traversed to determine whether all graph vertices have been numbered. If so, the numbering ends; if not, it means that there are multiple loops, then the branch covering logic S2 is executed, and the graph vertices of the flattened topology connection graph are traversed again and the branch covering logic S2 is executed until all graph vertices have been numbered. The above-mentioned numbering advancement logic that incorporates the influence of loop import graph vertices and multiple loops is defined as import graph vertex and loop covering logic S3. Logic S2 includes path advancement numbering logic S1 and branch covering logic S2, which is applicable to the numbering advancement of the entire flattened topology structure connection graph.

4. A method for assembling and storing the coefficient matrix of two-fluid six-equation applicable to Modelica according to claim 1, characterized in that, The assembly of the six equation coefficient components of the two fluids in step 3 is based on the discrete format of the mass conservation equation, mass conservation difference equation, vapor phase energy conservation equation, liquid phase energy conservation equation, momentum conservation equation and momentum conservation difference equation.

5. A method for assembling and storing the coefficient matrix of two-fluid six equations applicable to Modelica according to claim 4, characterized in that, The discrete formats of the equations in step 3 are as follows: (1) Mass conservation and equation (2) Mass conservation difference equation (3) Vapor phase energy conservation equation (4) Liquid phase energy conservation equation (5) Momentum conservation and equations (6) Momentum conservation difference equation Indicates that the variable is a temporary result, Indicates the use of the upwind scheme, i.e.: Where V is the volume of the control volume, A is the cross-sectional area of the control volume, P is the pressure, T is the temperature, U is the specific internal energy, H is the heat transfer coefficient, ρ is the density, α is the cavity fraction, v is the flow rate, ∈ indicates whether the wall is boiling or condensing, Δt is the time step, Δx is the length of the control volume, subscript g indicates the vapor phase, subscript f indicates the liquid phase, subscript L indicates the downstream control volume, subscript K indicates the upstream control volume, subscript j indicates the pipe number, j-1 indicates the upstream pipe, j+1 indicates the downstream pipe, superscript n indicates the nth time moment, and superscript n+1 indicates the n+1th time moment; After semi-implicit discretization, the unknowns in equations (1) to (4) include those at the current time step n+1 Among them, the first five variables are not the variables to be solved in the two-fluid six-equation model and need to be expanded using Taylor series, i.e.: Substituting equations (7) to (11) into equations (1) to (4), and by sorting out the variable coefficients and merging similar terms, equations (1) to (4) are uniformly expressed as follows: where a1 to a4 are respectively the variable coefficients in the equation, and s1 to s4 are respectively the variable coefficients in the equation, and b is the constant term of the equation. Among them, s1 to s4 are the convective transport coefficients of the mass or energy of the nozzle, which are related to the control volume number, nozzle number, and equation number, and can be written as: s n，j，q，k , where n is the control volume number, j is the nozzle number, q is the equation number, and k is 1 for the vapor phase g or 2 for the liquid phase f; In the momentum conservation equation (5) and the momentum conservation difference equation (6), the unknown variables only include those at the current (n + 1)th moment Through the arrangement of variable coefficients and the combination of like terms, equations (5) and (6) are uniformly expressed in the following form wherein, g1, g2, w1, and w2 are respectively the variable coefficients in the equation, and f is the constant term of the equation; to provide the complete data information required for the assembly of the two-fluid six-equation coefficient matrix, combining equations (12) and (13), the content stored in the data container corresponding to the basic node unit includes the model name where the basic node unit is located, the type of the basic node unit, the number of the basic node unit, the upstream and downstream node information, the structural parameters, the characteristic parameters, the initial parameters, and the generalized source term. The generalized source term 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 physical property calculation, nozzle variable calculation based on control volume variables, and control volume variable calculation based on nozzle variables.

6. A method for assembling and storing the coefficient matrix of two-fluid six-equation applicable to Modelica according to claim 5, characterized in that, In step 3, according to the flattened topological structure corresponding to the thermal-hydraulic system model and the basic node unit numbers, calculate and update in real time the data information of the basic node units of the control volume type and the basic node units of the nozzle type, and obtain the complete data information required for assembling the two-fluid six-equation coefficient matrix. The specific calculation steps are as follows: Step 31: Obtain the upstream and downstream nozzle numbers of the basic node units of the control volume type according to the overall diagram of the flattened topological structure, so as to establish the constraint relationship between the solution variables of pressure, void fraction, vapor specific internal energy, and liquid specific internal energy based on the control volume number and the solution variables of vapor velocity and liquid velocity based on the nozzle number in the mass conservation discrete equation and the energy conservation discrete equation; Step 32: Obtain the upstream and downstream control volume numbers of the basic node units of the nozzle type according to the overall diagram of the flattened topological structure, so as to establish the constraint relationship between the solution variables such as vapor velocity and liquid velocity based on the nozzle number and the pressure variable based on the control volume number in the momentum conservation discrete equation; Step 33: The vapor velocity and liquid velocity of the control volume required by the momentum conservation equation are obtained by averaging the vapor and liquid velocities of the upstream and downstream nozzles, and are stored in the corresponding data container of the basic node unit of the control volume type; The physical property parameters such as density, viscosity coefficient, thermal conductivity, and specific heat at constant pressure of the nozzle, as well as the scalar quantities such as void fraction, vapor specific internal energy, and vapor specific internal energy, are all the physical property parameters and scalar quantities of the incoming control volume, and are stored in the corresponding data container of the basic node unit of the nozzle type and are directly used in the calculation of the source term closure relationship of the momentum conservation discrete equation; Step 35: For the basic node units of the control volume type, calculate the generalized source terms of the mass and energy conservation discrete equations, and store the calculation results in the data container of the control volume. The calculation of the generalized source terms includes the calculation of physical properties, interphase mass transfer, interphase heat transfer, wall heat transfer, energy dissipation and other source terms, and according to the upstream and downstream nozzle numbers of the control volume, calculate the coefficients and constant terms of each solution variable of the mass and energy conservation discrete equations sorted by formula (12), including a1 to a4, s1 to s4, b, which can be directly used when assembling the coefficient matrix. Among them, s1 to s4 need to correspond to the upstream and downstream nozzle numbers of the control volume, the equation number, and the type of vapor and liquid velocity; Step 36: For the basic node units of the nozzle type, calculate the generalized source terms of the momentum conservation discrete equation, and store the calculation results in the data container of the control volume. The calculation of the generalized source terms includes the calculation of wall friction, interphase friction, local resistance, artificial viscosity terms, etc., as well as the treatment of special pumps, critical flow, etc., and according to the upstream and downstream control volume numbers of the nozzle, calculate the coefficients and constant terms of each solution variable of the momentum conservation sum equation and the momentum conservation difference sorted by formula (13), including g1, g2, w1, w2, f, which can be directly used when assembling the coefficient matrix; Step 37: A control volume may have multiple upstream and downstream nozzles, but each nozzle has only one upstream control volume and one downstream control volume. When dealing with the multi-pass situation, it is expressed by equation (13). When the control volume has one more upstream nozzle, two more terms are added to the mass conservation and equation of the control volume. and It is sorted out according to equation (12), that is: As can be seen from equation (14), the left scalar array, the variable coefficient array of the scalar, and the constant term remain unchanged, while the right vector array increases and The variable coefficient array of the corresponding right vector increases and 7. A method for assembling and storing the coefficient matrix of two-fluid six-equation applicable to Modelica according to claim 1, characterized in that, In step 4, the coefficient matrix of the two-fluid six equations is obtained by arranging according to equations (12) and (13). The j-th element a of the i-th row of the matrix ij corresponds to the coefficient of the j-th variable in the i-th equation. The j-th variable refers to the variable ranked in the j-th row of the variable array. The row elements of the coefficient matrix correspond to the coefficients of all the variables to be solved. When a certain variable is not in this equation, its coefficient is 0. The large sparse matrix is stored in row compression format. In the thermal-hydraulic system model, it is assumed that there are N non-boundary control volumes and M non-boundary nozzles. There are 4N control volume variables, including pressure, void fraction, specific internal energy of the vapor phase, specific internal energy of the liquid phase, and 2M nozzle variables, including vapor phase flow velocity and liquid phase flow velocity. The mass conservation and equation, mass conservation difference equation, vapor phase energy conservation equation, and liquid phase energy conservation equation of the N control volumes are arranged according to equation (12) and combined into matrix form: AX = SV + B (20) Where, A is a 4N×4N coefficient matrix for solving scalars such as control volume pressure, void fraction, vapor specific internal energy, and liquid specific internal energy. X is an array of 4N control volume variables. S is a 4N×2M coefficient matrix for solving vectors such as vapor velocity and liquid velocity at the upstream and downstream nozzles of the control volume. V is an array of 2M nozzle variables. B is an array of 4N constant terms.

8. A method for assembling and storing the coefficient matrix of two-fluid six equations applicable to Modelica according to claim 7, characterized in that, The assembly steps of the coefficient matrix of the mass and energy conservation discrete equations in step 4 are as follows: Step 41: Obtain the number of non-boundary control volumes N, the number of non-boundary nozzles M, the number of flow boundary nozzles M_TDJ, and the number of multi-pass extra nozzles M_ExtraConnectJun according to the topological structure information. Step 42: Calculate the number of non-zero elements of matrix A as 16N, and the number of non-zero elements of matrix S as 12N + 6. Step 43: Apply for the memory of the row-compressed format arrays of matrices A and S according to the number of non-zero elements of matrices A and S, including: A_data, A_colindex, A_rowindex, and S_data, S_colindex, S_rowindex, and apply for the memory of the constant term value B. Step 44: Set the mass and energy conservation equations and variable sorting as: Step 45: In the non-boundary control volume range, select control volumes in sequence according to the numbering from i = 1, where i is a natural number from 1 to N, obtain the data containers of the i-th control volume, and read the structural parameters, characteristic parameters, initial parameters, state parameters, physical property parameters, generalized source term calculation results stored in the data containers, as well as the solution variable coefficients and constants of the mass and energy conservation discrete equations sorted according to Equation (12), including variable coefficients a1 to a4, s1 to s4, and b. The variable coefficients and constants are different in different conservation discrete equations, and s1 to s4 need to correspond to the control volume number, connection pipe number, and equation number. s1 to s4 can all be rewritten as: s n，j，q，k , where n is the control volume number, j is the connection pipe number, q is the equation number, k is 1 for the vapor phase g or 2 for the liquid phase f. If there are multiple passes, the number of variable coefficients of the solution vector will increase. Assume that m upstream or downstream connection pipes are added. Then, the variable coefficients of the solution vector of each type of equation will increase from s1 to s4 to s1 to s 4+2m ; Step 46: Calculate the number of rows and columns of the variable coefficients of the scalars to be solved in matrix A according to the set equations and variable sorting. Step 47: Store the variable coefficients in the row-compressed format array of matrix A according to the number of rows and columns of the variable coefficients of the scalars to be solved in matrix A. Step 48: Obtain the upstream and downstream nozzle information of the i-th control volume according to the topological structure, clarify the numbers and data information of each upstream and downstream nozzle. Assume the number of the nozzle is j, and confirm whether the nozzle numbered j is a flow boundary nozzle. It is divided into the following situations: Situation 1: If the nozzle numbered j is a flow boundary nozzle, then read the physical property parameters, structural parameters, vapor-liquid velocity and other data information of the data container of the nozzle numbered j, and the variable coefficients of the vapor-liquid velocity of the nozzle numbered j in the mass conservation and equation, mass conservation difference equation, vapor phase energy conservation equation and liquid phase mass conservation equation in the data container of the i-th control volume. Calculate the expression term of the vapor-liquid velocity of the nozzle numbered j in the equation. This expression term is equal to the vapor-liquid velocity of the nozzle multiplied by its corresponding variable coefficient in the equation, and add it to the constant term of the corresponding equation of the i-th control volume, update the constant term of the data container, and jump to step 410. Case 2: If the nozzle numbered j is a non-boundary nozzle, then according to the equation sorting in Equation (21) and the variable sorting in Equation (22), calculate the row number and column number of the variable coefficient of the solution vector in the S matrix. Taking the mass conservation difference equation as an example, the row number of the variable coefficient corresponding to the vapor velocity of the nozzle numbered j in the mass conservation difference equation of the i-th control volume in the S matrix is 4i - 2, and the column number is j - m_tdj, where m_tdj is the number of nozzles from No. 1 to j that belong to the flow boundary nozzles, and its value is the variable coefficient s of the vapor velocity of the nozzle numbered j in the mass conservation difference equation stored in the data container of the i-th control volume n，j，2，1 , where: n is the number of the i-th control volume, and j is the nozzle number. Step 49: Store the variable coefficients in the row-compressed format array of matrix S according to the number of rows and columns of the variable coefficients of the solution vector in matrix S. Step 410: Store the constant terms b of the mass conservation and equation, mass conservation difference equation, vapor phase energy conservation equation and liquid phase mass conservation equation of the data container of the i-th control volume in the corresponding positions of the constant term array B according to the mass and energy conservation equation sorting. Step 411: Repeat steps 45 to 410, and sequentially complete the storage of the variable coefficients and constant terms for the control volume numbers i from 1 to N, and then the matrices A, S and the constant term array B can be assembled.

9. A method for assembling and storing the coefficient matrix of two-fluid six equations applicable to Modelica according to claim 1, characterized in that The assembly of the coefficient matrix of the momentum conservation discrete equation in step 5 includes the following: In the thermal-hydraulic system model, assuming there are N non-boundary control volumes and M non-boundary nozzles, then there are 4N control volume variables, including pressure, void fraction, specific internal energy of the vapor phase, specific internal energy of the liquid phase, and 2M nozzle variables, including vapor-phase flow velocity and liquid-phase flow velocity. The momentum conservation equations and momentum conservation difference equations of the M nozzles are sorted according to Equation (13) and combined into matrix form: GV = WP + F (42) Among them, G is a 2M×2M matrix, V is an array of 2M vapor-phase and liquid-phase flow velocity variables of the nozzles, W is a 2M×N matrix, and P is an array of pressure variables of the N control volumes.

10. A method for assembling and storing the coefficient matrix of two-fluid six-equation applicable to Modelica according to claim 9, characterized in that, The assembly of the coefficient matrix of the momentum conservation discrete equation in Step 5 described above includes the following: Step 51: Obtain the number N of non-boundary control volumes, the number M of non-boundary nozzles, and the number N_P of pressure boundary TDVs according to the topological structure information; Step 52: Calculate the number of non-zero elements 4M of the G matrix and the number of non-zero elements 4M - 2N_P of the W matrix; Step 53: Apply for the memory of the row-compressed format arrays of the G matrix and the W matrix according to the number of non-zero elements of the G matrix and the W matrix, including: G_data, G_colindex, G_rowindex and W_data, W_colindex, W_rowindex, and apply for the memory of the constant term value F; Step 54: Set the momentum conservation equation and variable sorting as: Step 55: Among the number M of non-boundary nozzles, sequentially select the nozzle with j = 1, where j is a natural number from 1 to M according to the nozzle number, obtain the data container of the jth nozzle, and read the structure parameters, characteristic parameters, initial parameters, state parameters, physical property parameters, calculation results of the generalized source term, and the variable coefficients and constant terms of each solution variable of the momentum conservation equation sorted according to Equation (13), including g1, g2, w1,; Step 56: Calculate the number of rows and columns of the variable coefficients of the solution vector in the G matrix according to the sorting of the momentum conservation equation of Equation (43) and the variable sorting of Equation (14); Step 57: Store the variable coefficients in the row-compressed format array of the G matrix according to the number of rows and columns of the variable coefficients of the solution vector in the G matrix; Step 58: According to the topological structure, obtain the upstream control volume and downstream control volume of the jth nozzle, clarify the numbers and data information of the upstream control volume and downstream control volume, assume the number of the upstream control volume or downstream control volume is n, and confirm whether the control volume numbered n is a pressure boundary control volume, which is divided into the following situations: Situation 1: Assume that the control volume numbered n is a pressure boundary control volume. Then read the data information such as physical property parameters, structure parameters, and pressure of the data container of the control volume numbered n, as well as the variable coefficients of the pressure of the control volume numbered n in the momentum conservation equation and the momentum conservation difference equation in the data container of the jth nozzle, calculate the expression term of the pressure of the control volume numbered n, which is equal to the control volume pressure multiplied by its corresponding variable coefficient in the equation, and add it to the constant term of the corresponding equation of the jth nozzle, update the constant term of the data container, and jump to Step 510; Case 2: Assume that the control volume numbered n is not a pressure boundary control volume. According to the sorting of the momentum conservation equation in Equation (43) and the sorting of variables in Equation (44), calculate the number of rows and columns of the variable coefficients of the control volume pressure in the W matrix; Step 59: According to the number of rows and columns of the variable coefficients of the solution vector in the W matrix, store the variable coefficients in the row-compressed format array of the W matrix; Step 510: According to the sorting of the momentum conservation equation in Equation (43), store the constant terms of the momentum conservation sum equation and the momentum conservation difference equation of the j-th nozzle data container in the corresponding positions of the constant term array F; Step 511: Repeat Step 55 to Step 510 to sequentially complete the storage of variable coefficients and constant terms for the control volume numbers j from 1 to M, and complete the assembly of the G matrix, the W matrix, and the constant term array F.

Images