Adaptive computation method for differential-algebraic network topology in multiphase thermal fluid simulation

Abstract

The application provides a differential algebraic network topology adaptive calculation method for multiphase thermal fluid simulation, comprising the following steps: establishing a multiphase thermal fluid differential algebraic graph network and a static superset graph network; initializing an activation matrix, a primary-secondary mapping table and an idle node linked list; solving full implicit equations simultaneously; reading a convergent solution and calculating an adaptive comprehensive criterion; updating a topology adaptive state based on the activation matrix; and reallocating a conservation parameter and promoting a source term compensation. The application can realize local topology adaptation and full implicit stable promotion of multiphase thermal fluid simulation without reestablishing a global unknown quantity set, without regenerating a connection relationship and without changing a sparse structure of a Jacobian matrix.

Description

Technical Field

[0001] This invention belongs to the field of system-level simulation calculation technology of aerospace cryogenic propellant thermofluid systems. Specifically, it relates to a calculation method that achieves structure-preserving topology adaptation within a differential algebraic graph network with a fixed maximum topology capacity by pre-setting a static superset graph network, activation matrix, free linked list and conservation mapping. It is particularly suitable for fast and high-precision solution of long-period multiphase thermal stratification and phase change processes in cryogenic propellant tanks. Background Technology

[0002] Cryogenic fluid management involves complex spatiotemporal multiphase thermodynamic processes characterized by cross-scale density ratios and intense thermal stratification. Existing numerical simulation tools face a trade-off between model accuracy and computational efficiency when handling long-term simulations.

[0003] The technical contradiction between computational efficiency and high-fidelity analysis in long-term, multi-scale simulations: While traditional computational fluid dynamics (CFD) methods offer high analytical accuracy, simulations lasting several hours (e.g., on-orbit missions of cryogenic liquid rockets) suffer from long computation cycles and enormous computational costs due to the vast computational domain, limitations in mesh reconstruction, and extremely small step sizes. Even with adaptive meshing techniques (CFD-AMR), data structure updates are typically still required around partial differential equation (PDE) meshes, error estimation, and local reconstruction, making it difficult to directly migrate to the fixed sparse structure solution framework of system-level differential algebraic equation (DAE) networks.

[0004] The simulation accuracy of simplified models of differential algebraic systems is insufficient: traditional 0D / 1D lumped parameter models have fast calculation speed, but they usually assume that physical quantities are uniform inside the nodes, making it difficult to increase local resolution as needed in thermal stratification zones, gas-liquid interfaces and near wall heat transfer boundaries, and also difficult to dynamically adjust the number of nodes without changing the structure of DAE unknowns.

[0005] Algorithm limitations: Traditional adaptive mesh or variable structure modeling often involves changes in the set of unknowns, connectivity relationships, or equation graphs during refinement, coarsening, or equation switching. This leads to changes in the global Jacobian matrix structure, failure of symbolic decomposition, and memory reallocation. For system-level simulation software that relies on a fixed structure for long-period fully implicit time integration, these structural changes can create real-time computational bottlenecks. Summary of the Invention

[0006] This invention provides an adaptive computation method for differential algebraic network topology in multiphase thermofluid simulation, addressing the contradiction in long-period simulations where local thermal stratification and phase transition interfaces require higher spatial resolution, while system-level DAE solvers rely on fixed unknown dimensions and fixed Jacobian matrix structures. This method achieves node refinement and merging without reconstructing the global sparse structure through pre-setting a static superset graph, managing activation matrices and idle lists, redistributing conserved mappings, and injecting external algebraic residual source vectors.

[0007] The method of this invention is as follows: Step 1: Establish the multiphase thermofluid differential algebraic graph network and the static superset graph network; The multiphase thermofluid computational domain is discretized into a basic fluid network consisting of master nodes and connector branches using the finite volume method. For each master node, candidate secondary nodes and candidate connection branches in a dormant state are pre-planned and programmed to construct a static differential algebraic superset graph network with pre-assignment characteristics and containing all potential topological connection relationships. Step 2: Initialize the activation matrix, primary and secondary mapping table, and free node linked list; It provides underlying data structure support for the dynamic evolution of superset graph networks. It initializes the 0 / 1 activation matrix for controlling the computing state of nodes and branches through unified numbering, constructs the primary and secondary mapping table for determining spatial subordinate relationships, and generates a one-dimensional idle node linked list for dynamically scheduling idle computing resources, thereby realizing static and orderly management of the state of nodes across the entire domain. Step 3: Solve the fully implicit equations simultaneously; We perform large-step fully implicit time integration on the static superset graph network. By introducing a continuously differentiable smooth step function into the convection term of the conservation equation, we eliminate non-smooth numerical abrupt changes caused by flow reversal or flux approaching zero, and ensure the continuity of the global Jacobian matrix and the robust convergence of nonlinear iteration under the adjustment of activation matrix parameters. Step 4: Read the converged solution and calculate the adaptive synthesis criterion; After the implicit equation solution converges, transient physical features reflecting the local non-uniformity of the flow field, such as temperature gradient, heat flux density and mass flow rate, are extracted from the current global solution vector. These features are then normalized without dimension and weighted summed using a preset feature scale to generate a dynamic adaptive comprehensive criterion for accurately driving local topology refinement or coarsening. Step 5: Topology-adaptive state update based on activation matrix; Based on the logical comparison between adaptive comprehensive criteria and preset dual thresholds, the state bits in the activation matrix are modified in a targeted manner, so as to achieve adaptive switching of local fluid topology resolution without changing the global equation dimension and matrix sparse structure. Step 6: Conservative parameter reallocation and compensation source term advancement; To ensure the physical conservation of the system before and after the topology state transition, the small numerical truncation error caused by the topology switch is transformed into an external algebraic residual source term and injected into the global equation system by performing strict proportional mapping of conserved quantities and gradient reconstruction of non-conserved quantities on active or dormant nodes. This allows the system to smoothly advance to the next time step without triggering Jacobian reorganization.

[0008] In one embodiment of the present invention, step 1 is specifically as follows: The computational domain size, initial thermodynamic state, and boundary conditions of the multiphase thermal fluid simulation object are read. Based on the phase differences of the medium, the computational domain of the multiphase thermal fluid is divided into three main categories: gas phase region, liquid phase region, and phase interface region. Using the Finite Volume Method (FVM), the computational domains of the same type are discretized into gas phase control volume nodes, liquid phase control volume nodes, and phase interface control volume nodes. Energy conservation, mass conservation, and state parameter equations for in-phase fluids are given within the control volume nodes. Momentum and heat transfer equations for in-phase flow, interphase heat and mass transfer, phase change, and boundary heat exchange are given in the connector branches. Nodes and branches are independent sets of computational equations, exposing variables such as temperature, pressure, and flow rate as interfaces. The connection between nodes and branches refers to mapping the port variables of the nodes at both ends of the connector branch to the same set of branch equations based on node number and connection relationship, and incorporating this branch equation into the global differential-algebraic equation set, enabling coupling solutions between nodes through mass, momentum, and energy relationships. After connection, a multiphase thermal fluid differential-algebraic graph network is formed, which serves as the master node network for subsequent topology adaptation. The fluid domain nodes are divided into two categories, master nodes and secondary nodes, according to adaptive logic. The initially discrete control volume nodes are the master nodes. For each master node, candidate secondary nodes and candidate connector branches are pre-established. During refinement, the two-dimensional network is refined once in the vertical and once in the radial directions, and the number of nodes in each direction doubles. Therefore, one master node in the first-level refinement corresponds to four candidate secondary nodes. In the second-level refinement, each first-level candidate secondary node corresponds to four second-level candidate secondary nodes. The refinement level can be customized by the user as needed. A primary and secondary branch is set between the primary node and the secondary node. The equations of the primary and secondary connection branches are the same as those of the primary node's connection branches, except that there is additional code to determine activation and dormancy. Inactive secondary nodes are called candidate secondary nodes, and inactive connection branches are called candidate connector branches. Activated candidate secondary nodes and candidate connector branches are directly called secondary node and connector branches. All candidate secondary nodes and candidate connector branches are also mapped to the same set of branch equations based on node number and connection relationship, and the branch equations are incorporated into the global differential algebra equation system to achieve connection. The method of pre-compiling candidate nodes and candidate connection branches for adaptation into the global differential algebraic equations before the start of the operation is called pre-allocation. Based on the multiphase thermofluid differential algebraic graph network, a pre-allocated static superset graph network is formed.

[0009] In another embodiment of the present invention, step 2 is specifically as follows: After the static superset graph network is established, the candidate secondary nodes and candidate connector branches that are set in advance are uniformly numbered. For each main node, there is a set of candidate secondary nodes in advance. Taking the first level of refinement as an example, one main node corresponds to four candidate secondary nodes, and the four candidate secondary nodes can be numbered in the spatial position order of lower left, lower right, upper left, and upper right. An activation matrix is ​​established, which is a 0 / 1 matrix used to record whether candidate secondary nodes and candidate connector branches participate in the current solution. For candidate secondary nodes, the elements of the activation matrix represent node state records, used to mark whether the candidate secondary node participates in the calculation of mass, energy and state variables at the current moment. For candidate connector branches, the elements of the activation matrix represent activation coefficients in the branch equations, used to control whether the mass flow rate, heat flow rate or phase change mass flow of the branch is written into the global equation set. The elements in the activation matrix take values ​​of 0 or 1, where 0 indicates that the corresponding candidate sub-node is in a dormant state, and 1 indicates that the corresponding candidate sub-node participates in the current physical flux calculation and the solution of the global differential algebraic equation. When the activation element corresponding to the candidate connector branch is 0, the mass flow rate and heat flow rate of that branch are taken as 0 in the equation. When the activation element is 1, the branch calculates the flux according to the pressure difference, temperature difference, heat transfer or phase change relationship, and solves it together with the global differential algebraic equation. During initialization, all elements in the activation matrix are set to 0, indicating that all candidate secondary nodes and candidate connector branches are in a dormant state and do not participate in the current physical flux calculation. A primary-secondary mapping table is established to record the fixed correspondence between each primary node and its preset candidate secondary nodes and candidate connector branches. This table is written once when the static superset graph is established. The primary-secondary mapping table includes at least: primary node number, phase domain type, refinement level, slot number of four candidate secondary nodes, four primary-secondary connector branch numbers, two horizontal connector branch numbers, two vertical connector branch numbers, and the current status of the group of candidate secondary nodes. An idle node linked list is established, which is a one-dimensional index table for storing the slot numbers of candidate secondary nodes. During initialization, a unique slot number is generated for each candidate secondary node in the order of phase domain number, master node number, refinement level, and candidate secondary node local number, and this slot number is added to the idle node linked list in sequence. Each element in the idle node linked list is an integer slot number, used to point to a candidate secondary node that has been pre-established in the static superset graph network but is not currently participating in the solution. When a candidate secondary node is activated, its slot number is removed from the idle node linked list; when a candidate secondary node is dormant, its slot number is added back to the idle node linked list.

[0010] In yet another embodiment of the present invention, step 3 is specifically as follows: After initialization in step 2, a fully implicit solution is performed on the static superset graph network, which is a system of differential algebraic equations. This system of equations is processed using a fully implicit time integration method, that is, the unknowns such as mass, internal energy, temperature, pressure, flow rate, and heat flow rate at the current time step are solved simultaneously at the same time level, and the equation residuals are made to meet the preset tolerance through iteration. If the residual convergence criterion is not met at the current time step, a step back mechanism is executed, and the iteration is restarted by reducing the time step size. For the convective enthalpy term in the connector branch, a smoothing upwind treatment is adopted. Specifically, when the mass flow rate direction is positive, the enthalpy of the upstream node is taken as the enthalpy entering the branch; the node with negative mass flow rate is the upstream node, and vice versa. When the mass flow rate direction is negative, the enthalpy of the downstream node on the other side is taken as the enthalpy entering the branch. To avoid discontinuous jumps when the mass flow rate is close to 0 or the flow direction is reversed, a continuously differentiable smooth step function is used to smoothly transition the enthalpy values ​​between the two directions. In the formula, The value of the smooth step function that is continuously differentiable is the flow weight, and its value ranges from 0 to 1. This indicates the mass flow rate in the connector branch; This represents the smoothing parameter, a small value greater than 0, used to control the smoothing range of the function from 0 to 1; The flow enthalpy for the connector branch used in convective energy calculations; and These are the enthalpy of the upstream node and the enthalpy of the downstream node, respectively, with the upstream and downstream components determined by mass flow rates. The positive or negative sign of the variable is determined; this treatment ensures that the connector branch maintains equation continuity even when the flow direction changes; During the solution process at one time step, the system of differential-algebraic equations remains unchanged, and the activation matrix is ​​only used as a parameter to control whether the candidate connector branch transmits flux; only at the adaptive control time is the activation matrix updated according to criteria such as temperature gradient, heat flux density or mass flow rate.

[0011] In another embodiment of the present invention, step 4 is specifically as follows: After solving the fully implicit equations in step 3, if the equation residuals at the current time step meet the preset tolerance, the solution at the current time step is considered to have converged. Only when the solution has converged will the adaptive control time step be entered and the adaptive criterion be calculated. When calculating the adaptive criterion, the state variables of each adaptive master node and its adjacent nodes or adjacent connector branches are read from the current global solution vector. The state variables include the temperature, pressure, mass, enthalpy or internal energy of the master node, as well as the mass flow rate, heat flow rate and phase change mass flow rate in the connector branches connected to the master node. Subsequently, the local non-uniformity index is calculated for each adaptive master node. The local non-uniformity index includes at least the temperature gradient index, the heat flux density index and the mass flow rate index. The aforementioned feature indicators are divided by their respective preset feature scales for dimensionless normalization, and then weighted and summed using the weight coefficients of each indicator to generate a dynamic adaptive comprehensive criterion. ; In the formula, For adaptive comprehensive criteria; The total number of physical field characteristic indicators involved in the comprehensive criterion calculation; These are index subscripts for characteristic indicators, representing different physical quantities; For the first The weight coefficients corresponding to the physical field characteristic indicators satisfy the normalization condition. ; Represents the extracted first Key physical field characteristic quantities; Is with Preset feature reference scales or tolerance thresholds with the same dimensions are used to achieve dimensionless representation of each physical quantity; Adaptive Synthesis Criterion Used to update the activation matrix in subsequent steps; when When the refinement conditions are met, it indicates that there is strong local non-uniformity near the main node, and its corresponding candidate secondary node and candidate connector branch need to be activated; when When the condition is below the dormancy condition, it indicates that the state change in this area is weak, and the corresponding candidate secondary node will be put into dormancy.

[0012] In step 4 of a specific embodiment of the present invention, the preset characteristic scales include the characteristic temperature difference of the storage tank and the design heat flux density.

[0013] In step 4 of another specific embodiment of the present invention, the first Key physical field characteristic quantities Represents spatial gradient, physical flux, or mass / energy source term.

[0014] In another embodiment of the present invention, step 5 is specifically as follows: After obtaining the adaptive synthesis criterion in step 4, determine whether each master node needs refinement or hibernation; let the first... The adaptive comprehensive criterion for each master node is: Refine the threshold to The hibernation threshold is And satisfy: Among them, the refined threshold Used to determine whether a local area needs improved spatial resolution; sleep threshold. Used to determine whether an activated candidate secondary node can exit the current solution; an interval is set between the two to avoid the same primary node repeatedly activating and dormant near the critical state; the threshold is refined. Determined after normalization based on temperature gradient threshold, heat flux density threshold, and mass flow threshold; Adaptive comprehensive criteria for primary nodes and candidate secondary nodes The adaptive comprehensive criterion value of different nodes is related to the thermal state of their neighboring nodes and is not necessarily the same; the refinement threshold of the main node is higher than that of the secondary node to avoid large-area activation and control system scale and computational cost; the dormancy threshold of the secondary node is lower than that of the main node. The activation matrix is ​​updated according to the following rules: If a master node is not currently refined, and its adaptive synthesis criterion is satisfied: If the primary node is deemed to need further refinement, then, based on the primary-secondary mapping table established in step 2, the candidate secondary nodes and primary-secondary connector branches corresponding to the primary node are found, and the corresponding activation matrix elements are changed from 0 to 1. If a candidate secondary node under a certain primary node is already active, and its adaptive synthesis criterion satisfies: If the local changes in the region have been weakened, it is considered that the region can be put into hibernation. Similarly, the candidate secondary nodes and candidate connector branches corresponding to the primary node are found according to the primary-secondary mapping table, and the corresponding activation matrix elements are changed from 1 to 0. If the adaptive synthesis criterion satisfies: Then keep the current activation matrix elements unchanged; The adaptive topology state update involves changing the activation matrix elements established in step 2. When an activation matrix element is changed from 0 to 1, the corresponding candidate connector branch begins to transfer mass flow, momentum flow, and heat flow in the static superset graph network established in step 1. When an activation matrix element is changed from 1 to 0, the flux of the corresponding candidate connector branch becomes 0 again. Thus, the adaptive switching of the local topology state is achieved while keeping the global set of unknowns and the equation structure unchanged.

[0015] In another embodiment of the present invention, step 6 is specifically as follows: When candidate secondary nodes adaptive comprehensive criteria Greater than or equal to the refinement threshold At that time, conserved quantities such as mass and internal energy are conserved according to volume ratios: In the formula, Assign a ratio to the parameters. and Representing the master node respectively i and its corresponding candidate secondary nodes j The physical quantity; For non-conservative quantities such as temperature, pressure, velocity, and enthalpy, reconstruction is performed using the states of the master node and its neighboring nodes; for temperature, local linear reconstruction is used. In the formula, The initial temperature of the candidate secondary node; Master node temperature; temperature gradient The temperature difference between the left, right, top, and bottom adjacent nodes of the master node is calculated. The center position of the master node; master node No. j The center position of each candidate sub-node; non-conserved quantities such as pressure and enthalpy are linearly reconstructed using the same method, or recalculated by the physical property parameter equations after mass and internal energy are written. When candidate secondary nodes adaptive comprehensive criteria Less than or equal to the sleep threshold When this happens, the additivity conservation quantities of the candidate secondary nodes are merged back into the corresponding primary node: In the formula, This represents the total number of candidate child nodes pre-assigned to a single master node, i.e., the number of split levels; the conserved quantity of the master node is equal to the sum of the conserved quantities of all directly candidate child nodes corresponding to that master node; this formula is only applicable to additive conserved quantities such as mass, internal energy, and total energy; for non-conserved quantities such as temperature, pressure, velocity, and enthalpy, they should be recalculated based on the combined mass and internal energy through physical property relations and equations of state; During the activation or dormancy of candidate secondary nodes, a small difference in mass or energy may occur due to variable reconstruction, state reinitialization, or numerical rounding. This difference is defined as the compensation amount. The compensation amount generated before and after the topology switch is divided by the time interval corresponding to this adaptive control, making it the source term per unit time. In the formula, The primary node compensation source item generated during topology state switching; This is the conserved quantity of the master node before the topology state switch; This refers to the conserved quantity after switching and completing the conservation mapping; For time step; Write the calculated physical quantity compensation source terms into the right-hand side of the corresponding conservation equation. For the mass and energy equation: in, Master node quality, The mass flow rate flowing into this master node, For the quality flow rate flowing out of the master node, Describe the internal energy of the master node, Enthalpy value, The power flowing into or out of the node; and These are the mass compensation source term and energy compensation source term generated by the topology state transition, respectively; if no difference is generated before and after the conservation mapping, the compensation source term is zero. After the above processing is completed, the solver continues to advance the time step using the updated activation matrix, candidate secondary node state values, and compensation source terms, and continues the fully implicit solution; the topology state evolution only changes the participation state of candidate secondary nodes and candidate connector branches through the activation matrix, and corrects the difference of conserved quantities through the compensation source terms; no global unknowns are added or deleted, and no connection relationships are re-established; since the candidate nodes and candidate connector branches in the static superset graph network have already entered the global equation during the initialization phase, the candidate sparse structure of the Jacobian matrix can be continuously reused, thereby completing the topology adaptive advancement without triggering the global equation reorganization.

[0016] The advantages of this invention are as follows: 1. Combining high accuracy and high computational efficiency: Addressing the issue of insufficient accuracy in existing system-level models (0D / 1D), this invention utilizes graph theory topology to discretize the computational domain into a network composed of control volume nodes and connectors. While retaining the fast characteristics of system-level simulation, this method can increase the local node density near thermal stratification zones, gas-liquid interfaces, and wall heat transfer boundaries, thereby improving the prediction accuracy of local temperature gradients, phase change source terms, and pressure evolution.

[0017] 2. Highly robust numerical solution and memory stability: Within the preset maximum topology capacity, the method of this invention maintains the total dimension of the global unknowns and the Jacobian matrix structure unchanged. Topology adaptation is only reflected in the updates of activation states, local numerical coefficients, and residual source terms. Therefore, compared with schemes that add or delete nodes or rebuild the mesh at runtime, it can reduce sparse structure reconstruction, computational overhead, and memory fluctuations, making large-step fully implicit time integration extremely stable.

[0018] 3. Overcoming the data structure limitations of system-level languages: The method of this invention is suitable for system-level simulation environments that are embedded in Modelica, FMU co-simulation or general program control. For modeling environments where it is inconvenient to change the equation structure at runtime, local topological state changes can be achieved by pre-allocating unknown slots and activation matrices, thereby bypassing the limitation of reconstructing the equation graph at runtime. Attached Figure Description

[0019] Figure 1 The overall flowchart of topology adaptive calculation for multiphase thermofluid simulation is shown; Figure 2 The diagram shown is a schematic of a multi-node mesh structure for discrete tank computation; Figure 3 The diagram shows the logical flow of the topology adaptive method coupled with the loop iteration.

[0020] Figure labels: 1-Gas phase main node; 2-Liquid phase main node; 3-Solid wall node; 4-Environment node; 5-Gas phase connection branch; 6-Liquid phase connection branch; 7-Gas-liquid interface branch; 8-Gas wall heat exchange branch; 9-Liquid wall heat exchange branch; 10-Solid heat conduction branch; 11-Solid and environment heat exchange branch. Detailed Implementation

[0021] The present invention will now be described in detail with reference to the accompanying drawings.

[0022] To achieve the above objectives, the method provided by this invention includes at least the following mutually cooperating steps, wherein static superset graph network pre-allocation, establishment of a free list, activation matrix management, and conservation mapping jointly define the topology adaptive path of this invention. Specifically: Step 1: Establish the multiphase thermofluid differential algebraic graph network and the static superset graph network.

[0023] This step aims to discretize the multiphase thermofluid computational domain into a basic fluid network consisting of master nodes and connector branches using the finite volume method. Based on this, candidate secondary nodes and candidate connection branches in a dormant state are planned and programmed into each master node in advance, thereby constructing a static differential algebraic superset graph network with pre-assignment characteristics and containing all potential topological connections.

[0024] The computational domain size, initial thermodynamic state, and boundary conditions of the multiphase thermal fluid simulation object are read. Based on the phase differences of the medium, the computational domain of the multiphase thermal fluid is divided into three main categories: gas phase region, liquid phase region, and phase interface region. Based on the finite volume method (FVM), the computational domains of the same type are discretized into gas phase control volume nodes, liquid phase control volume nodes, and phase interface control volume nodes (the finite volume method is a common discretization method for numerical calculations in fluid mechanics and is well-known to researchers in this field). Energy conservation, mass conservation, and state parameter equations between in-phase fluids are given within the control volume nodes; momentum and heat transfer equations for in-phase flow, interphase heat and mass transfer, phase change, and boundary heat exchange are given in the connector branches (Majumdar, AK, et al. "Generalized Fluid System Simulation Program (GFSSP) Version 6.0" .NASA TM-2016-218224). This document provides detailed forms of the node mass-energy equations and branch momentum equations; however, the specific equations and mathematical models are well-known to researchers in this field and are not the focus of this invention. Nodes and branches are independent sets of computational equations, each exposing variables such as temperature, pressure, and flow rate as interfaces. The connection between nodes and branches refers to mapping the port variables of the nodes at both ends of a connector branch to the same set of branch equations based on node numbers and connection relationships. These branch equations are then incorporated into the global differential-algebraic equation system, enabling coupled solutions between nodes through mass, momentum, and energy relationships. After connection, a multiphase thermofluid differential-algebraic graph network is formed, which serves as the master node network for subsequent topology adaptation.

[0025] The fluid domain nodes are divided into two categories based on adaptive logic: master nodes and secondary nodes. The initially discrete control volume nodes are the master nodes. For each master node, candidate secondary nodes and candidate connector branches are pre-established. During refinement, the 2D mesh is refined once in the vertical and once in the radial directions, doubling the number of nodes in each direction. Therefore, one master node in the first-level refinement corresponds to four candidate secondary nodes; in the second-level refinement, each first-level candidate secondary node corresponds to four second-level candidate secondary nodes. The number of refinement levels can be customized by the user as needed.

[0026] A primary and secondary branch are set up between the primary node and the secondary nodes. The equations of the primary and secondary connection branches are the same as those of the primary node's connection branch, except for the addition of code to determine activation and dormancy. Inactive secondary nodes are called candidate secondary nodes, and inactive connection branches are called candidate connector branches. In the following description, activated candidate secondary nodes and candidate connector branches will be directly referred to as secondary nodes and connector branches.

[0027] All candidate secondary nodes and candidate connector branches are also mapped to the same set of branch equations based on node number and connection relationship, and the branch equations are incorporated into the global differential algebraic equation system to achieve connection.

[0028] The method of pre-compiling candidate nodes and candidate connection branches for adaptation into the global differential algebraic equations before the start of the operation is called pre-allocation. Based on the multiphase thermofluid differential algebraic graph network, a pre-allocated static superset graph network is formed.

[0029] Step 2: Initialize the activation matrix, primary and secondary mapping table, and free node linked list.

[0030] This step aims to provide underlying data structure support for the dynamic evolution of superset graph networks. It initializes a 0 / 1 activation matrix for controlling the computational state of nodes and branches by uniform numbering, constructs a primary-secondary mapping table for determining spatial subordination, and generates a one-dimensional list of idle nodes for dynamically scheduling idle computing resources, thereby achieving static and orderly management of the state of nodes across the entire domain.

[0031] After the static superset graph network is established, the pre-defined candidate secondary nodes and candidate connector branches are first uniformly numbered. For each primary node, there is a pre-defined set of candidate secondary nodes. Taking first-level refinement as an example, one primary node corresponds to four candidate secondary nodes, and the four candidate secondary nodes can be numbered in the spatial order of lower left, lower right, upper left, and upper right.

[0032] An activation matrix is ​​established, which is a 0 / 1 matrix used to record whether candidate secondary nodes and candidate connector branches participate in the current solution. For candidate secondary nodes, the elements of the activation matrix represent node state records, used to mark whether the candidate secondary node participates in the calculation of mass, energy, and state variables at the current moment; for candidate connector branches, the elements of the activation matrix represent activation coefficients in the branch equations, used to control whether the mass flow rate, heat flow rate, or phase change mass flow of that branch is written into the global equation system.

[0033] The elements in the activation matrix take values ​​of 0 or 1, where 0 indicates that the corresponding candidate subnode is in a dormant state, and 1 indicates that the corresponding candidate subnode participates in the current physical flux calculation and the solution of the global differential algebraic equation. When the activation element corresponding to a candidate connector branch is 0, the mass flow rate and heat flow rate of that branch are taken as 0 in the equation; when the activation element is 1, the branch calculates the flux according to the pressure difference, temperature difference, heat transfer, or phase change relationship, and solves it together with the global differential algebraic equation.

[0034] During initialization, all elements in the activation matrix are set to 0, indicating that all candidate secondary nodes and candidate connector branches are in a dormant state and do not participate in the current physical flux calculation.

[0035] A primary-secondary mapping table is established to record the fixed correspondence between each primary node and its preset candidate secondary nodes and candidate connector branches. This table is not a new connection relationship temporarily generated during operation, but is written once when the static superset graph network is established. The primary-secondary mapping table includes at least: primary node number, phase domain type, refinement level, slot numbers of four candidate secondary nodes, four primary-secondary connector branch numbers, two horizontal connector branch numbers, two vertical connector branch numbers, and the current status of the group of candidate secondary nodes.

[0036] An idle node linked list is established, which is a one-dimensional index table for storing the slot numbers of candidate secondary nodes. During initialization, a unique slot number is generated for each candidate secondary node in the order of phase domain number, master node number, refinement level, and candidate secondary node local number, and this slot number is added to the idle node linked list sequentially. Each element in the idle node linked list is an integer slot number, used to point to a candidate secondary node that has been pre-established in the static superset graph network but is not currently participating in the solution. When a candidate secondary node is activated, its slot number is removed from the idle node linked list; when a candidate secondary node is dormant, its slot number is added back to the idle node linked list.

[0037] Step 3: Solve the fully implicit equations simultaneously.

[0038] This step aims to perform large-step fully implicit time integration on the static superset graph network. By introducing a continuously differentiable smooth step function into the convection term of the conservation equation, non-smooth numerical abrupt changes caused by flow reversal or flux approaching zero are eliminated, ensuring the continuity of the global Jacobian matrix under the adjustment of the activation matrix parameters and the robust convergence of the nonlinear iteration.

[0039] After initialization in step 2, a fully implicit solution is performed on the static superset graph network, i.e., a system of differential-algebraic equations (a common solution method for differential-algebraic equation systems, well-known to those skilled in the art). The system of equations is processed using a fully implicit time integration method, meaning that unknowns such as mass, internal energy, temperature, pressure, flow rate, and heat flow rate at the current time step are solved simultaneously at the same time level, and the equation residuals are iterated to meet a preset tolerance. If the current time step does not meet the residual convergence criterion, a step-back mechanism (a common solution method for differential-algebraic equation systems, well-known to those skilled in the art) is executed, and the iteration is restarted by reducing the time step size.

[0040] For the convective enthalpy term in the connector branch, a smoothing upwind treatment can be used (FRANKE R, CASELLA F, OTTER M et al. Stream Connectors – An Extension of Modelica for Device-Oriented Modeling of Convective Transport Phenomena [C / OL] / / The 7 International Modelica Conference, Como, Italy. 2009: 108-121 [2026-05-25]. DOI:10.3384 / ecp09430078). Section 4 (Numerical Properties) of that paper clearly points out that when the mass flow rate approaches zero, i.e., the flow reverses, directly using the upwind scheme will lead to singular solver matrices or numerical oscillations. Therefore, the convective enthalpy mixing equation at the connector must be regularized and smoothed. The paper uses a piecewise function implementation, which has the highest order continuity, while this invention uses an analytical function implementation, which has the continuity of second-order or higher derivatives. Specifically, when the mass flow rate is positive, the enthalpy of the upstream node (nodes with negative mass flow rates are upstream, and vice versa) is used as the enthalpy for entering the branch; when the mass flow rate is negative, the enthalpy of the downstream node on the other side is used as the enthalpy for entering the branch. To avoid discontinuous jumps when the mass flow rate approaches zero or the flow direction reverses, a continuously differentiable smooth step function can be used to smoothly transition between the enthalpy values ​​in the two directions: In the formula, This represents the value of a continuously differentiable smooth step function, also known as the flow weight, and its value ranges from 0 to 1. This indicates the mass flow rate in the connector branch; This represents the smoothing parameter, a small value greater than 0, used to control the smoothing range of the function from 0 to 1; The flow enthalpy for the connector branch used in convective energy calculations; and These are the enthalpy of the upstream node and the enthalpy of the downstream node, respectively, with the upstream and downstream components determined by mass flow rates. The sign of the variable is determined by the variable. This treatment ensures that the connector branch maintains equation continuity even when the flow direction changes, which is beneficial for the convergence of the fully implicit solver.

[0041] During the solution process at one time step, the system of differential-algebraic equations remains unchanged, and the activation matrix is ​​used only as a parameter to control whether the candidate connector branch transmits flux. Only at the adaptive control time is the activation matrix updated according to criteria such as temperature gradient, heat flux density, or mass flow rate (the update method is described in step 5).

[0042] Step 4: Read the converged solution and calculate the adaptive synthesis criterion.

[0043] This step aims to extract transient physical features reflecting the local non-uniformity of the flow field, such as temperature gradient, heat flux density, and mass flow rate, from the current global solution vector after the implicit equation solution converges. These features are then normalized dimensionlessly and weighted summed using a preset feature scale to generate a dynamic adaptive comprehensive criterion for accurately driving local topology refinement or coarsening.

[0044] After solving the fully implicit equations in step 3, if the equation residuals at the current time step meet the preset tolerance, the solution at the current time step is considered to have converged. Only when the solution has converged will the adaptive control time step begin and the adaptive criterion be calculated.

[0045] When calculating the adaptive criterion, the state variables of each adaptive master node and its adjacent nodes or adjacent connector branches are read from the current global solution vector (this reading is done automatically by retrieving variable values ​​within the code, a method well-known to researchers in the field). The state variables include the master node's temperature, pressure, mass, enthalpy, or internal energy, as well as the mass flow rate, heat flow rate, and phase change mass flow rate in the connector branches connected to that master node. Subsequently, a local non-uniformity index is calculated for each adaptive master node. This local non-uniformity index includes at least a temperature gradient index, a heat flux density index, and a mass flow rate index.

[0046] The aforementioned characteristic indicators are divided by their corresponding preset characteristic scales (such as tank characteristic temperature difference, design heat flux density, etc.) for dimensionless normalization, and then weighted and summed using the weight coefficients of each indicator to generate a dynamic adaptive comprehensive criterion. .

[0047] In the formula, For adaptive comprehensive criteria; The total number of physical field characteristic indicators involved in the comprehensive criterion calculation; These are index subscripts for characteristic indicators, representing different physical quantities; For the first The weight coefficients corresponding to the physical field characteristic indicators usually satisfy the normalization condition. ; Represents the extracted first Key physical field characteristics (which may represent spatial gradient, physical flux, or mass / energy source terms, etc.); Is with Preset feature reference scales or tolerance thresholds with the same dimensions are used to achieve dimensionless and normalized representation of each physical quantity.

[0048] Adaptive Synthesis Criterion Used to update the activation matrix in subsequent steps. When the refinement conditions are met, it indicates that there is strong local non-uniformity near the main node, and its corresponding candidate secondary node and candidate connector branch need to be activated; when When the condition is below the sleep condition, it indicates that the state change in this area is weak, and the corresponding candidate secondary node can be put into sleep mode.

[0049] Step 5: Topology-adaptive state update based on activation matrix.

[0050] This step aims to directly modify the state bits in the activation matrix based on the logical comparison between the adaptive comprehensive criterion and the preset dual thresholds (refinement threshold and dormancy threshold), thereby achieving adaptive switching of local fluid topology resolution without changing the global equation dimension and the sparse structure of the matrix.

[0051] After obtaining the adaptive synthesis criterion in step 4, each master node is individually assessed to determine whether refinement or hibernation is necessary. Let the first master node be... The adaptive comprehensive criterion for each master node is: Refine the threshold to The hibernation threshold is And satisfy: Among them, the refined threshold Used to determine whether a local area needs improved spatial resolution; sleep threshold. This is used to determine whether an activated candidate secondary node can exit the current solution. An interval is set between the two to prevent the same primary node from repeatedly activating and dormant near a critical state. Refined threshold. The weights can be determined after normalization based on the temperature gradient threshold, heat flux density threshold, and mass flow threshold. For example, under a multiphase thermofluid calculation condition, the dimensionless weights of each physical index can be set as follows: Temperature gradient weight =0.6, heat flux density weight =0.3, phase change mass flow weight =0.1. When the temperature change in the region where a master node is located is extremely drastic, and its local temperature gradient just reaches the set allowable limit (i.e., the temperature term normalization result is 1), and the heat flux density at that location reaches 1.33 times the critical threshold (i.e., the heat flux term normalization result is 1.33), even if no phase change occurs at this time (the mass flow normalization result is 0), the comprehensive criterion calculation result for this node is: At this point, the comprehensive judgment criteria are... Refinement threshold The system immediately determines that the spatial resolution of the currently pre-allocated grid is insufficient to support accurate solutions and triggers the refinement of the candidate secondary nodes under the master node.

[0052] Adaptive comprehensive criteria for primary nodes and candidate secondary nodes The calculation method in step 4 yields the following results. The adaptive synthesis criterion value for different nodes is related to the thermal state of their neighboring nodes and is not necessarily the same. The refinement threshold for the primary node is higher than that for the secondary node to avoid large-scale activation, thus controlling system size and computational cost. Similarly, the dormancy threshold for the secondary node is lower than that for the primary node. (Refinement threshold) With sleep threshold The method for determining the values ​​is consistent with that of CFD-AMR theory and is not the focus of this invention. The relevant principles can be found in the published paper: Beisiegel, N., Castro, CE, & Behrens, J. Metrics for Performance Quantification of AdaptiveMesh Refinement. Journal of Scientific Computing, 87, 36 (2021). https: / / doi.org / 10.1007 / s10915-021-01423-0 The main innovation in step 5 lies in the updating of the activation matrix. The activation matrix is ​​updated according to the following rules: If a master node is not currently refined, and its adaptive synthesis criterion is satisfied: If the primary node is deemed to require further refinement, then, based on the primary-secondary mapping table established in step 2, the candidate secondary nodes and primary-secondary connector branches corresponding to the primary node are found, and the corresponding activation matrix elements are changed from 0 to 1.

[0053] If a candidate secondary node under a certain primary node is already active, and its adaptive synthesis criterion satisfies: If the local changes in the region have weakened, it can be put into dormancy. At this point, the candidate secondary nodes and candidate connector branches corresponding to the primary node are found according to the primary-secondary mapping table, and the corresponding activation matrix elements are changed from 1 to 0.

[0054] If the adaptive synthesis criterion satisfies: Then the current active matrix elements remain unchanged.

[0055] The so-called topology adaptive state update in this step does not involve rebuilding nodes, branches, or equations, but rather changing the activation matrix elements established in step 2. When an activation matrix element is changed from 0 to 1, the corresponding candidate connector branch begins to transfer mass flow, momentum flow, and heat flow in the static superset graph network pre-established in step 1; when an activation matrix element is changed from 1 to 0, the flux of the corresponding candidate connector branch returns to 0. Thus, adaptive switching of the local topology state is achieved while maintaining the global set of unknowns and the equation structure unchanged.

[0056] Step 6: Conservation parameter redistribution and compensation source term advancement.

[0057] This step aims to ensure the physical conservation of the system before and after the topology state transition. By performing strict proportional mapping of conserved quantities and gradient reconstruction of non-conserved quantities on active or dormant nodes, the small numerical truncation error caused by topology switching is transformed into external algebraic residual source terms and injected into the global equation system, so that the system can smoothly advance to the next time step without triggering Jacobian reorganization.

[0058] When candidate secondary nodes adaptive comprehensive criteria Greater than or equal to the refinement threshold At that time, conserved quantities such as mass and internal energy are conserved according to volume ratios: In the formula, Assign a ratio to the parameters. and Representing the master node respectively i and its corresponding candidate secondary nodes j The physical quantity.

[0059] For non-conserved quantities such as temperature, pressure, velocity, and enthalpy, direct summation or simple proportional allocation is not advisable. Instead, reconstruction should be performed using the states of the master node and its neighboring nodes. Taking temperature as an example, local linear reconstruction can be used: In the formula, Initial temperature of candidate secondary nodes; Master node temperature; temperature gradient It can be calculated from the temperature difference between the left, right, top, and bottom adjacent nodes of the master node. , and These represent the unit vectors in the horizontal and vertical directions of the master node, respectively. Represents the local absolute temperature of a fluid. and (These represent the spatial dimensions of the two-dimensional discrete computational domain, respectively, in the horizontal and vertical directions). The center position of the master node; master node No. j The center position of each candidate secondary node. Non-conserved quantities such as pressure and enthalpy can be linearly reconstructed using the same method, or recalculated from the physical property equations after mass and internal energy are written. (This type of linear reconstruction belongs to the conventional local reconstruction method in the finite volume method, which is well known to researchers in this field. The focus of this invention is not to propose this reconstruction format alone, but to use it for writing the initial activation values ​​of candidate secondary nodes in a static superset graph network.) When candidate secondary nodes adaptive comprehensive criteria Less than or equal to the sleep threshold When this happens, the additivity conservation quantities of the candidate secondary nodes are merged back into the corresponding primary node: In the formula, This represents the total number of candidate child nodes pre-assigned to a single master node, i.e., the number of split levels. The conserved quantity of a master node is equal to the sum of the conserved quantities of all directly candidate child nodes corresponding to that master node. This formula only applies to additively conserved quantities such as mass, internal energy, and total energy; for non-conserved quantities such as temperature, pressure, velocity, and enthalpy, they should be recalculated based on the combined mass and internal energy using property relations and equations of state.

[0060] During the activation or dormancy of candidate secondary nodes, a small difference in mass or energy may occur due to variable reconstruction, state reinitialization, or numerical rounding. This difference is defined as the compensation amount. The compensation amount generated before and after the topology switch is divided by the time interval corresponding to this adaptive control, making it the source term per unit time: In the formula, The primary node compensation source item generated during topology state switching; This is the conserved quantity of the master node before the topology state switch; This refers to the conserved quantity after switching and completing the conservation mapping; For time step.

[0061] Furthermore, the calculated physical quantity compensation source terms are written into the right-hand side of the corresponding conservation equation, taking the mass and energy equation as an example: in, Master node quality, The mass flow rate flowing into this master node, For the quality flow rate flowing out of the master node, Describe the internal energy of the master node, Enthalpy value, This refers to the power flowing into or out of the node. and These are the mass compensation source term and energy compensation source term generated by the topology state transition, respectively. If no difference is generated before and after the conservation mapping, the compensation source term is zero.

[0062] After the above processing, the solver continues time-stepping using the updated activation matrix, candidate secondary node state values, and compensation source terms, continuing the fully implicit solution. Topology state evolution only changes the participation states of candidate secondary nodes and candidate connector branches through the activation matrix and corrects the differences in conserved quantities through the compensation source terms; it does not add or delete global unknowns, nor does it re-establish connections. Since the candidate nodes and candidate connector branches in the static superset graph network have already entered the global equations during the initialization phase, the candidate sparse structure of the Jacobian matrix can be continuously reused, thus achieving adaptive topology advancement without triggering global equation reorganization.

[0063] In this invention, "node" primarily refers to the control volume in a system-level DAE graph network, while "branch" or "connector" refers to the connection between nodes via equations that transmit mass, momentum, or energy conservation relationships. Unless otherwise specified, "mesh refinement" in this specification refers to changing the activation state of control volume nodes.

[0064] It should be noted that the terms "node refinement," "mesh refinement," "cell refinement," "activation," and "enablement" used in the following descriptions can be interpreted as equivalent in case of ambiguity. Similarly, "node merging," "mesh recycling," "cell coarsening," "inactive," "frozen," and "dormant" can also be interpreted as equivalent.

[0065] The terms "node," "mesh," and "cell" appearing in the following description should be interpreted broadly, as should "branch" and "connection," which are commonly used expressions in discrete systems. For example, adaptive node and adaptive mesh have the same meaning, and nodes connected by branches can also be described as meshes transmitting conserved information. Those skilled in the art can understand the specific meaning of these terms in this invention according to the specific circumstances.

[0066] All technical details not described in detail in this invention are conventional techniques in the field and will not be repeated here.

[0067] For those skilled in the art, the present invention can be practiced without requiring some of the details described above. The following description of the embodiments is merely intended to provide a better understanding of the invention by illustrating examples.

[0068] Example 1 Figure 1This is a calculation process according to an embodiment of the present invention, applied to the thermodynamic simulation of cryogenic storage tanks. This process can be executed by Julia, Python, C++ or a system-level co-simulation program. Figure 2 This is a schematic diagram of a multi-node mesh structure for spatial discrete modeling according to an embodiment of the present invention. Figure 3 It is a logical flowchart of the coupling between the topology adaptive method and the loop iteration.

[0069] It should be noted that, due to Figure 2 This diagram illustrates the connection structure of various physical characteristic nodes in a cryogenic tank. Nodes of the same type have been merged, and the connection status between individual nodes is used to demonstrate the most intuitive discrete structure. In a real discrete mesh system, the number of each type of node is not unique, and the initial number of nodes is determined by the simulation accuracy. Specifically, the spatial arrangement, node combination, and connection method between nodes are not fixed and are also related to the initial physical domain calculated, and can be customized according to the simulation task.

[0070] like Figure 1 As shown, this embodiment is a complete instantiation of the aforementioned general method, "Adaptive Calculation Method of Differential Algebraic Network Topology for Multiphase Thermofluid Simulation," in a specific aerospace engineering scenario (long-term thermodynamic simulation of cryogenic propellant tanks). This embodiment deeply integrates the aforementioned general core architecture (i.e., superset graph network, activation matrix, and residual injection) with specific physical processes within the cryogenic tank (such as gas-liquid interface phase transitions, thermal stratification in the gas pillow region, and complex solid-liquid-gas heat exchange). It details the mathematical transformation and parameter configuration of the general feature criteria under specific multiphase thermodynamic scenarios, thereby verifying and fully demonstrating the invention's objective of maintaining a static macroscopic Jacobian matrix while balancing system-level computational efficiency and local high-fidelity accuracy. The steps include: The simulation object in this embodiment is a cryogenic propellant tank containing a gas phase region, a liquid phase region, a gas-liquid interface, a solid wall surface, and an external environmental boundary.

[0071] First, read the basic structural parameters of the tank, including diameter, height, filling rate, wall thickness, initial pressure, initial fluid temperature, initial solid domain wall temperature, external heat input, and thermal radiation boundary conditions of the external vacuum environment.

[0072] Figure 2 As shown, the fluid domain of the cryogenic tank is discretized (discretization method: the finite volume method is a common discretization method for numerical calculations in fluid mechanics research and is well known to researchers in this field) into gas phase principal nodes 1 and liquid phase principal nodes 2. The sidewalls, top wall, and bottom wall are discretized into solid wall nodes 3, and the external environment is represented as environment nodes 4. Each fluid principal node stores state quantities or number of iterations such as mass, internal energy, temperature, pressure, density, and enthalpy.

[0073] Furthermore, basic connector branches are established between the aforementioned nodes. A gas phase connection branch 5 is established between adjacent gas phase main nodes, a liquid phase connection branch 6 is established between adjacent liquid phase main nodes, a gas-liquid interface branch 7 is established at the gas-liquid interface, a gas wall heat transfer branch 8 and a liquid wall heat transfer branch 9 are established between fluid nodes and solid wall nodes respectively, a solid heat conduction branch 10 is established between solid wall nodes, and a solid-environment heat transfer branch 11 is established between solid wall nodes and environmental nodes. Each connector branch is used to calculate mass flow rate, heat flow rate, momentum pressure drop, or phase change mass flow.

[0074] In this embodiment, the equation-oriented modeling approach of the Julia language is used to define node components, connection branch components, and their port variables. Adjacent nodes are coupled through the port relationships of the connection branches. Port variables include pressure, mass flow rate, specific enthalpy, temperature, and heat flow. Based on the port connection relationships, the program writes the mass flow, enthalpy flow, and heat flow of each port into the corresponding node's mass conservation equation and energy conservation equation. Simultaneously, the connection branches provide the pressure drop-flow rate relationship or momentum transfer relationship between adjacent nodes, thus forming a complete set of discrete equations. The branches and nodes thus form a multiphase thermofluid differential algebraic graph network for the cryogenic tank, as shown in the schematic diagram below. Figure 2 As shown.

[0075] Furthermore, after the basic graph network is established, a static superset graph network is created for the gas phase master nodes and liquid phase master nodes. Taking two-dimensional first-level refinement as an example, each master node corresponds to four candidate secondary nodes in advance. The four candidate secondary nodes can be refined according to their spatial positions: lower left, lower right, upper left, and upper right. Four master-secondary connector branches are pre-set between the master node and the four candidate secondary nodes. Two horizontal connector branches are pre-set between adjacent left and right candidate secondary nodes, and two vertical connector branches are pre-set between adjacent up and down candidate secondary nodes. If two-level refinement is allowed, each first-level candidate secondary node continues to serve as a master node in a relative sense, and pre-sets its next-level candidate secondary node and candidate connector branches.

[0076] The aforementioned candidate secondary nodes and candidate connector branches are incorporated into the global unknown set and candidate equation structure before the solution process begins. When inactive, they do not transmit actual physical fluxes; when activated, they are solved along with the global differential-algebraic equations. This network, with its pre-established maximum number of candidate nodes and candidate connector branches, is a static superset graph network.

[0077] Further, initialize the activation matrix, primary and secondary mapping table, and free node linked list.

[0078] After the static superset network of the cryogenic storage tank is established, the main nodes and their pre-assigned candidate secondary nodes and candidate connector branches in the gas phase region, liquid phase region, gas-liquid interface region, and solid wall region are uniformly numbered. Among them, the candidate connector branches include the main-secondary connector branches between the main nodes and the candidate secondary nodes, the lateral connector branches between adjacent candidate secondary nodes, and the longitudinal connector branches.

[0079] Subsequently, initialization is performed according to the activation matrix, primary-secondary mapping table, and idle node list defined in the invention method: each primary node, its candidate secondary node slots, and candidate connector branch numbers are written into the primary-secondary mapping table. The activation matrix elements corresponding to the candidate secondary nodes and candidate connector branches are set to 0, so that only the basic primary node network participates in the solution during the initial simulation of the cryogenic tank. All candidate secondary node slot numbers are written into the idle node list in the phase domain order of gas phase, liquid phase, gas-liquid interface, and solid wall.

[0080] Therefore, in subsequent adaptive control moments, the program can read the corresponding slots and branches through the primary and secondary mapping table based on the local temperature gradient, heat flux density, or mass flow changes in the cryogenic tank, move candidate secondary node slots out or write them back through the idle node linked list, and change the state of whether candidate secondary nodes and candidate connector branches participate in the global differential algebraic equation solution through the activation matrix.

[0081] After initializing the index relationships and activation states of the static superset graph network of the cryogenic storage tank, the system of differential-algebraic equations corresponding to the static superset graph network is solved implicitly.

[0082] The set of equations includes mass conservation equations, energy conservation equations, momentum conservation equations, heat transfer equations, phase transition equations, state parameter equations, heat conduction equations, radiation and convection heat transfer equations, and physical property relation equations for gas phase principal nodes, liquid phase principal nodes, gas-liquid interface nodes, solid wall nodes, and external boundary nodes.

[0083] Within each time step, the solver solves both the nodal equations and the connector branch equations simultaneously. The currently participating nodes and connector branches are determined by the activation matrix. For the basic master node and basic connector branch, their equations are always involved in the solution. For candidate secondary nodes and candidate connector branches, their state variables and fluxes only participate in the solution of the current equation when the corresponding activation matrix element is 1. When the corresponding activation matrix element is 0, the mass flow rate, heat flow rate, or phase change mass flow rate of the candidate connector branch is set to 0 in the equation, and the corresponding candidate secondary node does not exchange mass and energy with surrounding nodes.

[0084] For a candidate linker branch, its activation state can be written in the form of a branch flux multiplier: In the formula, Indicates a candidate connector branch; This is the activation matrix element corresponding to this branch, with a value of 0 or 1; and These are the actual mass flow rate and heat flow rate written into the equation system, respectively; and These are the candidate mass flow rate and candidate heat flow rate calculated based on pressure difference, flow resistance, temperature difference, heat transfer coefficient, or phase change relationship, respectively. When When, the actual flux of this branch is 0; when At that time, this branch calculates the flux according to the corresponding physical relationship.

[0085] For the convective energy term in the connector branch, a continuously differentiable smooth step function can be used for smoothing the upwind treatment. That is, when the mass flow direction is clear, the upstream node enthalpy is taken as the flow enthalpy; when the mass flow is close to 0 or the flow direction reverses, a continuously differentiable function is used to smoothly transition the flow enthalpy between the two node enthalpies. This treatment avoids sudden jumps in convective enthalpy values ​​when the flow direction changes, ensuring the continuity of the branch equations during the fully implicit iteration process, which is beneficial for improving the convergence stability of the coupled solutions for gas-liquid flow, phase change, and heat transfer in cryogenic tanks.

[0086] During the solution process at the current time step, the solver iteratively calculates the global differential-algebraic equation system and determines whether the residuals meet the preset convergence tolerance (the convergence tolerance is typically 1e-5 to 1e-6, a value well-known to researchers in the field of numerical computation). If the residuals do not meet the convergence requirements, the time step is reduced and the iteration is restarted; if the residuals meet the convergence requirements, the converged solution at the current time step is saved as input for subsequent transient physical feature extraction and adaptive criterion calculation.

[0087] After the fully implicit solution converges, the system reads the state variables or fluxes such as temperature, pressure, mass, internal energy, enthalpy, heat flow, mass flow and phase change mass flow of each adaptive master node and its adjacent nodes from the solution vector of the current time step.

[0088] Specifically, for cryogenic tank simulation, the adaptive criterion for master nodes can be determined based on the local non-uniformity of the gas phase region, liquid phase region, gas-liquid interface region, and wall heat transfer region. These non-uniformity indices include, but are not limited to: temperature gradient between adjacent master nodes, heat flux density near the wall or phase interface, mass flow rate of in-phase connector branches, and phase change mass flow at the gas-liquid interface. k Example calculation of temperature gradient of one master node: In the formula, Indicates in At time t, the first time in the spatial discrete network Temperature gradient characteristic index at each master node. Indicates the current physical simulation time; Indicates the first k Temperature of each master node; The temperature of its adjacent master node; For two adjacent control body nodes (the first one) Node and the Spatial characteristic distance between the center points of nodes.

[0089] When considering temperature gradient, heat flux density, and mass flow changes simultaneously, the first k The comprehensive adaptive criterion for each master node can be expressed as: In the formula, For the first k A comprehensive adaptive criterion for each master node; For real-time calculation of the longitudinal temperature gradient; Heat flux density; This is the phase change rate or mass flow rate term; , , These represent the weighting coefficients of temperature gradient, heat flux density, and phase change rate in the adaptive judgment, respectively. , , These are the pre-set trigger thresholds for the temperature gradient, heat flux density, and phase change rate, respectively. These weights and thresholds can be pre-set based on the importance of the gas-liquid interface, wall heat exchange zone, or localized high-flow zone in the cryogenic tank.

[0090] Furthermore, after obtaining the comprehensive adaptive criterion Then, for each adaptive master node, determine whether its corresponding candidate secondary node needs to be activated or suspended. To refine the threshold, The threshold value for sleep mode is defined, and it satisfies the following conditions: Setting an interval between the two can reduce the repeated activation and dormancy of the same master node near the critical state. The specific judgment conditions and corresponding operations are as follows: 1. If the first k Each node satisfies: At this time, the pre-allocation activation mechanism is triggered: the system reads the candidate secondary node slot number, primary and secondary connector branch number, horizontal connector branch number, and vertical connector branch number corresponding to the primary node according to the primary-secondary mapping table; then, the corresponding candidate secondary node slot number is removed from the idle node list, and the corresponding element in the activation matrix is ​​changed from 0 to 1. After that, the corresponding candidate secondary node and candidate connector branch begin to participate in the solution of the global differential-algebraic equation in the next time step.

[0091] 2. If the first k Each node satisfies: When a candidate secondary node is active, a sleep operation is performed. The system also reads the corresponding candidate secondary node slot and candidate connector branch number from the primary-secondary mapping table. After completing the merging of conserved quantities and the calculation of compensation source terms, the corresponding element in the activation matrix is ​​changed from 1 to 0, and the corresponding candidate secondary node slot number is rewritten back to the idle node linked list.

[0092] 3. When the comprehensive adaptive criterion satisfies: If the activation matrix elements corresponding to the master node remain unchanged, continue using the current topology state to proceed to the next time step for solution.

[0093] Specifically, when an activation matrix element changes from 0 to 1, it indicates that the corresponding master node needs to increase its local spatial resolution. For additively conserved quantities such as mass, internal energy, and volume, they are allocated from the master node to the candidate secondary nodes according to the volume ratio of the candidate secondary nodes. For quantities that can be allocated geometrically, such as mass, volume, and adjacent heat exchange areas, they are initialized according to the volume or area ratio of the candidate secondary nodes relative to the master node. The first master node corresponding to the Taking a candidate secondary node as an example, it can be written as: In the formula, Indicates the first The first master node corresponding to the The initial quality of each candidate secondary node; Indicates the first The quality of each master node before activation; Indicates the first The control volume of each candidate secondary node; Indicates the first The volume of the control body of each master node; Indicates the first The heat exchange area between each candidate secondary node and its adjacent nodes or boundaries; Indicates the first The heat exchange area of ​​each master node before activation; Indicates the first The volume allocation ratio of each candidate secondary node relative to the primary node; Indicates the first The heat exchange area allocation ratio of each candidate secondary node relative to the primary node.

[0094] For non-conservative quantities such as velocity, pressure, and temperature, a linear or quadratic distribution is reconstructed using the gradient information of the master node and its neighboring units. Then, the initial values ​​of each sub-unit center are interpolated from this reconstructed distribution. Taking temperature as an example, a local linear reconstruction can be used: In the formula, Initial temperature of candidate secondary nodes; Master node temperature; temperature gradient It can be calculated from the temperature difference between the left, right, top, and bottom adjacent nodes of the master node; The center position of the master node; Let j be the center position of the candidate sub-node. Non-conserved quantities such as pressure and enthalpy can be reconstructed using similar linear methods, or recalculated from the physical property equations after mass and internal energy have been written.

[0095] When an element in the activation matrix changes from 1 to 0, it indicates that the corresponding candidate secondary node needs to go to sleep. At this point, the additivity conservation values ​​of this group of candidate secondary nodes are merged back into the corresponding primary node. In the formula, , and These represent the mass, internal energy, and heat transfer area of ​​the n candidate subnodes before hibernation; Indicates the first The quality of each master node before activation; Indicates the first The volume of the control body of each master node; Indicates the first The heat exchange area corresponding to each master node before activation.

[0096] This summation relation applies only to additable conserved quantities such as mass, internal energy, and total energy. Non-conserved quantities such as temperature, pressure, density, and enthalpy should be recalculated based on the combined mass, internal energy, volume, and physical property relationships, rather than being directly added together.

[0097] During the activation or dormancy of candidate secondary nodes, small differences in mass or energy may occur due to variable reconstruction, state reinitialization, or numerical rounding. These differences are defined as compensation amounts. The compensation amounts generated before and after the topology switch are divided by the time interval corresponding to this adaptive control, making it the compensation source term per unit time. Let the mass and internal energy of the region corresponding to the master node before the topology state switch be... , The mass and internal energy after the topological state switch and completion of the conservation mapping are respectively , Then the compensation source term is: In the formula, This is the quality compensation source term for the region corresponding to the k-th master node; For energy compensation source terms; is the time step. If there is no difference before and after the conservation mapping, then the compensation source term is 0.

[0098] Furthermore, the calculated physical quantity compensation source terms are written into the right-hand side of the corresponding conservation equation, taking the mass and energy equation as an example: in, and These are the mass compensation source term and energy compensation source term generated by the topology state transition, respectively. If no difference is generated before and after the conservation mapping, the compensation source term is zero.

[0099] In this way, the small conserved differences caused by topology state changes are factored into the current adaptive control time interval for processing. After completing the activation, dormancy, conserved quantity redistribution, and compensation source term writing of candidate secondary nodes, the solver continues the fully implicit solution for the next time step using the updated activation matrix, node state variables, and compensation source terms. Throughout the process, the global unknown dimension, node number, candidate connector branch number, and sparse structure of the Jacobian matrix remain unchanged; topology adaptation is only manifested as the updating of activation matrix elements, candidate secondary node states, and compensation source terms. The above process is executed cyclically with each time step until the set termination time of the cryogenic tank transient simulation is reached.

[0100] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

[0101] Example 2 This embodiment also provides a multiphase thermofluid topology adaptive simulation system for performing the above method. This system is a functional implementation of the "differential algebraic network topology adaptive calculation method for multiphase thermofluid simulation" described in Embodiment 1. The functional logic of each module of this system corresponds one-to-one with the calculation steps in Embodiment 1. Through the collaborative work of each functional unit, dynamic topology adaptation is achieved on the basis of a static superset graph network. The system's logical architecture includes: 1. Discrete construction module for fluid graph network.

[0102] The discrete modeling unit is used to read the geometry, gas-liquid phase domain division, wall structure, initial temperature, initial pressure, filling state, boundary conditions, and external heat flow conditions of the cryogenic tank, and discretizes the simulation object into gas phase master nodes, liquid phase master nodes, gas-liquid interface nodes, solid wall nodes, and external boundary nodes. The connector branching unit is used to establish in-phase flow branches, interphase heat transfer branches, phase change branches, wall heat conduction branches, and environmental heat transfer branches based on the spatial adjacency and physical interaction relationships between nodes, and to provide the corresponding mass, momentum, and energy equations.

[0103] The output of this module includes the set of basic master nodes, the set of basic connector branches, the initial state variables of the nodes, the boundary conditions, and the corresponding differential-algebraic equation structure. These outputs are passed to the static superset graph pre-assignment module to determine which master nodes are allowed to establish candidate secondary nodes and candidate connector branches; and to the fully implicit solution and smoothing module as the foundational equations and initial conditions for subsequent simultaneous solutions.

[0104] 2. Superset Graph Network Pre-allocation and Initialization Module. This module is used to build a static superset graph network based on the base graph network and initialize the activation matrix, primary-secondary mapping table, and idle node list. This module includes a slot pre-allocation unit and an index state initialization unit.

[0105] The slot pre-allocation unit receives the set of main nodes, connector branches, and adaptive region information output by the fluid network discrete construction module. Based on the set refinement level and the number of candidate secondary nodes, it pre-sets candidate secondary node slots and candidate connector branches for the main nodes in the gas phase region, liquid phase region, gas-liquid interface region, and wall heat exchange region. The index state initialization unit is used to establish a main-secondary mapping table, establishing a fixed correspondence between each main node and its candidate secondary node slots, main-secondary connector branches, horizontal connector branches, and vertical connector branches; at the same time, it establishes an activation matrix and sets the initial elements to 0; and then writes all candidate secondary node slot numbers into the idle node linked list.

[0106] The output of this module is a static superset graph, an activation matrix, a primary-secondary mapping table, and a list of idle nodes. These outputs are passed to the fully implicit solution and smoothing module to determine which candidate secondary nodes and candidate connector branches participate in the equation solving at the current time step; they are also passed to the activation adaptive topology reorganization module to subsequently change the activation state based on criteria.

[0107] 3. Fully Implicit Solver Module. This module performs a fully implicit simultaneous solution of the differential-algebraic equations in the current active state. Specifically, it includes a fully implicit coupled solver unit that receives the fundamental node equations and connector branch equations output by the fluid graph discretization module, and simultaneously receives the activation matrix, primary-secondary mapping table, and candidate connector branch information output by the static superset graph pre-allocation module.

[0108] During the solution process, the basic master node and basic connector branch always participate in the equation solution. Whether the candidate secondary node and candidate connector branch participate in the solution is determined by the 0 or 1 state in the activation matrix. The smoothing upwind treatment unit is used to employ a continuously differentiable smooth step function in the convective energy transfer branch, so that the flow enthalpy transitions continuously when the flow direction changes, reducing non-smooth disturbances in the coupled solution of gas-liquid flow, phase change and heat transfer in the cryogenic tank.

[0109] The output of this module is the converged solution vector at the current time step, the residual convergence status, and necessary time step control information. If the equation system fails to converge, the module returns a step size reduction signal to the time progression process and solves the problem again. If the equation system converges, the converged solution vector is passed to the feature extraction and criterion generation module.

[0110] 4. Feature Extraction and Adaptive Criterion Generation Module. This module extracts transient physical features from the fully implicit solution results and generates adaptive comprehensive criteria. Specifically, it includes a solution vector extraction unit, a physical index calculation unit, and a criterion synthesis calculation unit.

[0111] The solution vector extraction unit receives the converged solution vector output by the fully implicit solution and smoothing module, and reads data such as temperature, pressure, mass, enthalpy, heat flux, mass flow rate, and phase change mass flow of each master node and its adjacent nodes. The physical index calculation unit is used to calculate local non-uniformity indices such as temperature gradient, heat flux density, in-phase flow intensity, and phase change mass flow. The criterion synthesis calculation unit is used to normalize the above indices according to preset thresholds and sum them according to weights to obtain the comprehensive adaptive criterion corresponding to each adaptive master node.

[0112] The output of this module is a comprehensive adaptive criterion for each adaptive master node, as well as local physical indicators used to determine whether to activate or suspend. These outputs are then passed to the activation adaptive topology reorganization module.

[0113] 5. Activate the adaptive topology reorganization module. This module updates the activation matrix based on a comprehensive adaptive criterion and performs the removal or writing back of candidate secondary node slots. It includes a dual-threshold judgment unit, an activation matrix update unit, and a slot status update unit.

[0114] The dual-threshold judgment unit receives the comprehensive adaptive criterion output by the feature extraction and criterion generation module and compares it with the activation threshold and the dormancy threshold. If the comprehensive adaptive criterion is greater than or equal to the activation threshold, it is determined that the corresponding primary node needs to activate the candidate secondary node; if the comprehensive adaptive criterion is less than or equal to the dormancy threshold, it is determined that the corresponding candidate secondary node can be dormant; if the comprehensive adaptive criterion is between the two, the current activation state remains unchanged.

[0115] The activation matrix update unit receives the primary-secondary mapping table provided by the static superset graph network pre-allocation module, finds the corresponding candidate secondary node slot and candidate connector branch number based on the primary node number, and changes the corresponding element in the activation matrix from 0 to 1, or from 1 to 0. The slot status update unit operates on the idle node list according to the activation matrix update result: when a candidate secondary node is activated, its slot number is removed from the idle node list; when a candidate secondary node is dormant, its slot number is written back to the idle node list.

[0116] The output of this module is the updated activation matrix, the updated list of idle nodes, and the candidate secondary nodes and candidate connector branch numbers for which state transitions have occurred. These outputs are passed to the conservation mapping and compensation source term module for performing state initialization, conservation merging, and compensation source term calculation; simultaneously, they are fed back to the fully implicit solution and smoothing module for solving the equations in the next time step.

[0117] 6. Conservative Mapping and Compensation Source Item Module. This module performs conservative quantity reallocation, non-conservative state quantity reconstruction, and compensation source item writing when candidate secondary nodes are activated or asleep. It includes a conservative quantity mapping unit, a non-conservative state reconstruction unit, and a compensation source item generation unit.

[0118] The conservation quantity mapping unit receives the state switching information output by the activated adaptive topology reorganization module and reads the corresponding master node and candidate secondary node slots according to the master-secondary mapping table. When a candidate secondary node is activated, the conservation quantities such as the mass and internal energy of the master node are allocated to the candidate secondary node according to the volume ratio of the candidate secondary node; when a candidate secondary node is dormant, the conservation quantities such as the mass and internal energy of the candidate secondary node are summed and merged back into the corresponding master node. The non-conservative state reconstruction unit is used to process state quantities such as temperature, pressure, density, and enthalpy. Temperature can be reconstructed based on the gradient information of the master node and its adjacent nodes, while pressure, density, and enthalpy can be recalculated based on the merged or allocated conservation quantities and physical property relationships.

[0119] The compensation source term generation unit is used to calculate the mass or energy difference that may occur before and after the topology state switch, and converts the difference into a mass compensation source term or an energy compensation source term according to the adaptive control step size. The compensation source term is written into the right-hand side of the mass conservation equation and energy conservation equation of the corresponding master node or candidate secondary node, rather than changing the dimension of the global unknowns or rebuilding the equation structure.

[0120] The output of this module is the updated node state variables and compensation source terms. These outputs are then returned to the fully implicit solver and smoothing module as the initial state and source term inputs for solving the global differential-algebraic equations at the next time step.

[0121] Through the aforementioned module relationships, the system forms a closed-loop workflow: the fluid graph discretization construction module provides basic nodes, branches, and equations; the static superset graph pre-allocation module establishes candidate secondary nodes, candidate connector branches, and index states based on these; the fully implicit solution and smoothing module solves the equations under the current activation state; the feature extraction and criterion generation module calculates the comprehensive adaptive criterion based on the converged solution; the activation adaptive topology reorganization module updates the activation matrix and the idle node list based on the criterion; the conservation mapping and compensation source term module updates the node state and writes the compensation source terms, and returns the results to the solution module. This process is repeated until the set termination time for the cryogenic tank multiphase thermofluid simulation is reached.

[0122] In the above embodiments, the thermodynamic criteria for driving node and branch refinement or dormancy are mainly the local temperature gradient and heat flux density gradient. For different multiphase fluid physics scenarios, this criterion can also be replaced by physical indicators characterizing space thermodynamic non-uniformity, such as density gradient, specific volume gradient, local pressure jump rate, or drastic change rate of phase change mass source term, all of which can achieve the same dynamic refinement driving effect.

[0123] The differential-algebraic equation (DAE) system constructed in this invention is not limited to time integration using the DASSL solver. Those skilled in the art can replace it with Radau IIA, CVODE, LSODI, or other higher-order implicit solvers suitable for strongly rigid ordinary differential / differential-algebraic equations, depending on the stiffness characteristics and time scale of the fluid system. Regardless of the solver used, the underlying layer still reuses the same pre-allocated activation adaptive network mechanism and zero-memory reallocation logic.

[0124] In terms of engineering implementation and solution methods, this invention can not only run independently in a general-purpose programming language (such as Python, Julia, C++, etc.) with complete features, but also adopt an alternative approach of co-simulation with a system modeling language or system simulation platform (such as Modelica, Simulink). That is, after completing the node modeling and discrete system construction, the fluid nodes and equations of this invention are encapsulated into a functional model prototype (FMU). Through a general-purpose programming language as the main control program, cross-platform co-simulation is performed with a system-level simulation environment (such as Modelica, Simulink, etc.), which can also realize multiphysics coupling calculation and topology adaptation.

[0125] When the three conditions of "the conservation equation can be represented as a DAE graph network, the maximum topological capacity can be predetermined, and the topological changes can be mapped to active states and residual source terms" are met, this method can also be applied to thermofluid networks with lumped and distributed parameters, and can also be extended to simulation modeling tasks of discretized network characteristics such as thermal hydraulic networks, cryogenic transport systems, or building thermal environments.

[0126] The mechanism of this invention lies in transforming the node refinement and dormancy behavior in physical space into activation state updates and external algebraic residual source term updates within a fixed-capacity DAE graph network. Since parent nodes, candidate child nodes, and candidate connectors occupy fixed unknown slots during the initialization phase, the subsequent adaptive process does not change the total dimension of the global unknowns or the Jacobian matrix structure. Simultaneously, conserved projection and residual compensation ensure that mass, volume, and energy remain consistent before and after activation or dormancy. Therefore, the system-level DAE solver can reuse a fixed sparse structure while local spatial resolution changes, reducing runtime structure reconstruction and memory fluctuations, and resolving the technical contradiction between system-level static equation solving and mesh adaptive dynamic reconstruction.

[0127] The core inventive points of this invention are as follows: 1. Construction of variable topology networks based on differential algebraic equations (DAE). The complex fluid domain is discretized in two-dimensional space into a dynamic fluid graph network with coupled conservation of mass, momentum, and energy, consisting of control volumes (nodes) and momentum connectors (branches). This provides the equational basis for subsequent local node activation and dormancy within a fixed maximum topology capacity.

[0128] 2. Adaptive Algorithm for Fixed-Capacity Superset Graph Network and Activation State Management. This algorithm overcomes the bottleneck of expensive memory reallocation caused by the need for dynamic addition and deletion of matrix dimensions in traditional adaptive mesh refinement. It pre-generates slots for parent nodes, candidate child nodes, and candidate connectors based on the maximum splitting level. Activation and sleep are then performed within the allocated slots using the activation matrix, parent-child mapping table, and free list, ensuring that the adaptive process does not change the total dimension of the DAE unknowns or the candidate sparse connection relationships.

[0129] 3. Conservative Mapping and Residual Injection Mechanism. During node activation or dormancy, projective redistribution is performed on conserved quantities such as mass, volume, and energy, while reconstructive interpolation is performed on state quantities such as temperature and pressure. The compensation quantities caused by topological state transitions are represented as external algebraic residual source vectors. After this source vector is injected into the fixed DAE residual equations, only the numerical residuals and local coefficients are updated, rather than dynamically adding or deleting matrix dimensions. Activating residual source terms through continuous algebraic mapping makes the large-step fully implicit time integral extremely stable.

[0130] This invention provides an adaptive computational method for differential algebraic network topology in multiphase thermofluid simulation. The method discretizes the multiphase thermofluid computational domain into a differential algebraic graph network composed of control volume nodes and connector branches. Before computation begins, a static superset graph network containing master nodes, candidate secondary nodes, and candidate connector branches is pre-established. An activation matrix records whether candidate secondary nodes and candidate connector branches participate in the current solution. A master-secondary mapping table records the fixed correspondence between master nodes and candidate secondary nodes and candidate connector branches. An idle node linked list manages the slots of candidate secondary nodes that are not participating in the solution. During simulation, a fully implicit method is used to simultaneously solve the mass, momentum, and energy equations in the current active state. Physical indicators such as temperature gradient, heat flux density, mass flow rate, or phase change mass flow are extracted from the converged solution to generate a comprehensive adaptive criterion. When the comprehensive adaptive criterion reaches the activation criterion, the activation matrix elements of the corresponding candidate secondary node and candidate connector branch are changed from 0 to 1; when the comprehensive adaptive criterion is lower than the dormancy criterion, the corresponding activation matrix elements are changed from 1 to 0. When candidate nodes are activated or dormant, conserved quantities such as mass and internal energy are allocated or merged, and state variables such as temperature, pressure, and enthalpy are reconstructed or recalculated. The mass or energy difference resulting from topology state switching is converted into compensation source terms and written into the corresponding conservation equations. Thus, without re-establishing the global set of unknowns, regenerating connectivity relationships, or altering the sparse structure of the Jacobian matrix, local topology adaptation and fully implicit stable propagation in multiphase thermofluid simulation are achieved.

Claims

1. An adaptive computation method for differential-algebraic network topology for multiphase thermofluid simulation, characterized in that, Specifically as follows: Step 1: Establish the multiphase thermofluid differential algebraic graph network and the static superset graph network; The multiphase thermofluid computational domain is discretized into a basic fluid network consisting of master nodes and connector branches using the finite volume method. For each master node, candidate secondary nodes and candidate connection branches in a dormant state are pre-planned and programmed to construct a static differential algebraic superset graph network with pre-assignment characteristics and containing all potential topological connection relationships. Step 2: Initialize the activation matrix, primary and secondary mapping table, and free node linked list; It provides underlying data structure support for the dynamic evolution of superset graph networks. It initializes the 0 / 1 activation matrix for controlling the computing state of nodes and branches through unified numbering, constructs the primary and secondary mapping table for determining spatial subordinate relationships, and generates a one-dimensional idle node linked list for dynamically scheduling idle computing resources, thereby realizing static and orderly management of the state of nodes across the entire domain. Step 3: Solve the fully implicit equations simultaneously; We perform large-step fully implicit time integration on the static superset graph network. By introducing a continuously differentiable smooth step function into the convection term of the conservation equation, we eliminate non-smooth numerical abrupt changes caused by flow reversal or flux approaching zero, and ensure the continuity of the global Jacobian matrix and the robust convergence of nonlinear iteration under the adjustment of activation matrix parameters. Step 4: Read the converged solution and calculate the adaptive synthesis criterion; After the implicit equation solution converges, transient physical features reflecting the local non-uniformity of the flow field, such as temperature gradient, heat flux density and mass flow rate, are extracted from the current global solution vector. These features are then normalized without dimension and weighted summed using a preset feature scale to generate a dynamic adaptive comprehensive criterion for accurately driving local topology refinement or coarsening. Step 5: Topology-adaptive state update based on activation matrix; Based on the logical comparison between adaptive comprehensive criteria and preset dual thresholds, the state bits in the activation matrix are modified in a targeted manner, so as to achieve adaptive switching of local fluid topology resolution without changing the global equation dimension and matrix sparse structure. Step 6: Conservative parameter reallocation and compensation source term advancement; To ensure the physical conservation of the system before and after the topology state transition, the small numerical truncation error caused by the topology switch is transformed into an external algebraic residual source term and injected into the global equation system by performing strict proportional mapping of conserved quantities and gradient reconstruction of non-conserved quantities on active or dormant nodes. This allows the system to smoothly advance to the next time step without triggering Jacobian reorganization.

2. The adaptive calculation method for differential algebraic network topology for multiphase thermofluid simulation as described in claim 2, characterized in that, Step 1 is as follows: The computational domain size, initial thermodynamic state, and boundary conditions of the multiphase thermal fluid simulation object are read; the computational domain of the multiphase thermal fluid is divided into three categories according to the phase difference of the medium: gas phase region, liquid phase region, and phase interface region. Based on the Finite Volume Method (FVM), the computational domain is discretized into gas phase control volume nodes, liquid phase control volume nodes, and phase interface control volume nodes. Within the control volume nodes, energy conservation, mass conservation, and state parameter equations for in-phase fluids are given. Connector branches provide momentum and heat transfer equations for in-phase flow, interphase heat and mass transfer, phase change, and boundary heat exchange. Nodes and branches are independent sets of computational equations, exposing variables such as temperature, pressure, and flow rate as interfaces. The connection between nodes and branches refers to mapping the port variables of the nodes at both ends of the connector branch to the same set of branch equations based on the node number and connection relationship, and incorporating the branch equations into the global differential-algebraic equation system, so that the nodes can achieve coupled solution through mass, momentum, and energy relationships. After connection, a multiphase thermofluid differential algebraic graph network is formed, which serves as the master node network for subsequent topology adaptation. Fluid domain nodes are divided into two categories, master nodes and secondary nodes, according to adaptive logic; the initially discrete control volume nodes are the master nodes. For each main node, candidate secondary nodes and candidate connector branches are pre-established; during refinement, the two-dimensional graph is refined once in the vertical and radial directions, and the number of nodes in each direction is doubled. Therefore, one main node in the first-level refinement corresponds to four candidate secondary nodes; during the second-level refinement, each first-level candidate secondary node corresponds to four second-level candidate secondary nodes. The refinement level can be customized by the user as needed. A primary and secondary branch is set between the primary node and the secondary node. The equations of the primary and secondary connection branches are the same as the equations of the primary node's connection branches, with the addition of code to determine activation and sleep. Inactive secondary nodes are called candidate secondary nodes, and inactive connection branches are called candidate connector branches; The activated candidate secondary nodes and candidate connector branches are directly referred to as secondary node and connector branches; All candidate secondary nodes and candidate connector branches are also mapped to the same set of branch equations based on node number and connection relationship, and the branch equations are incorporated into the global differential algebra equation system to achieve connection. The method of pre-compiling candidate nodes and candidate connection branches for adaptation into the global differential algebraic equations before the start of the operation is called pre-allocation. Based on the multiphase thermofluid differential algebraic graph network, a pre-allocated static superset graph network is formed.

3. The adaptive calculation method for differential algebraic network topology for multiphase thermofluid simulation as described in claim 3, characterized in that, Step 2 is as follows: After the static superset graph network is established, the candidate secondary nodes and candidate connector branches that are set in advance are uniformly numbered. For each main node, there is a set of candidate secondary nodes in advance. Taking the first level of refinement as an example, one main node corresponds to four candidate secondary nodes, and the four candidate secondary nodes can be numbered in the spatial position order of lower left, lower right, upper left, and upper right. Establish an activation matrix, which is a 0 / 1 matrix used to record whether candidate secondary nodes and candidate connector branches participate in the current solution; For candidate secondary nodes, the activation matrix elements are represented as node state records, used to mark whether the candidate secondary node participates in the calculation of mass, energy, and state quantities at the current moment; for candidate connector branches, the activation matrix elements are represented as activation coefficients in the branch equations, used to control whether the mass flow rate, heat flow rate, or phase change mass flow rate of the branch is written into the global equation set. The elements in the activation matrix take values ​​of 0 or 1, where 0 indicates that the corresponding candidate node is in a dormant state, and 1 indicates that the corresponding candidate node participates in the current physical flux calculation and the solution of the global differential algebraic equation; When the activation element corresponding to the candidate connector branch is 0, the mass flow rate and heat flow rate of that branch are taken as 0 in the equation; When the activation element is 1, this branch calculates the flux according to the pressure difference, temperature difference, heat transfer or phase change relationship, and solves it together with the global differential algebraic equation; During initialization, all elements in the activation matrix are set to 0, indicating that all candidate secondary nodes and candidate connector branches are in a dormant state and do not participate in the current physical flux calculation. Establish a primary-secondary mapping table to record the fixed correspondence between each primary node and its preset candidate secondary nodes and candidate connector branches; This table is written once when the static superset graph network is created; The primary and secondary mapping table includes at least: primary node number, phase domain type, refinement level, slot number of four candidate secondary nodes, four primary and secondary connector branch numbers, two horizontal connector branch numbers, two vertical connector branch numbers, and the current status of the group of candidate secondary nodes. An idle node linked list is established, which is a one-dimensional index table for storing the slot numbers of candidate secondary nodes. During initialization, a unique slot number is generated for each candidate secondary node in the order of phase domain number, master node number, refinement level, and candidate secondary node local number, and this slot number is added to the idle node linked list in sequence. Each element in the idle node linked list is an integer slot number, used to point to a candidate secondary node that has been pre-established in the static superset graph network but is not currently participating in the solution. When a candidate secondary node is activated, its slot number is removed from the idle node linked list; when a candidate secondary node is dormant, its slot number is added back to the idle node linked list.

4. The adaptive calculation method for differential algebraic network topology for multiphase thermofluid simulation as described in claim 1, characterized in that, Step 3 is as follows: After initialization is completed in step 2, a fully implicit solution is performed on the static superset graph network, which is a system of differential algebraic equations. The system of equations is processed using a fully implicit time integration method, that is, the unknowns such as mass, internal energy, temperature, pressure, flow rate and heat flow rate in the current time step are solved simultaneously at the same time level, and the equation residuals are made to meet the preset tolerance through iteration. If the residual convergence criterion is not met at the current time step, the step size back-off mechanism is executed, and the time step is reduced and the iteration is restarted. For the convective enthalpy term in the connector branch, a smoothing upwind treatment is adopted. Specifically, when the mass flow rate direction is positive, the enthalpy of the upstream node is taken as the enthalpy entering the branch; the node with negative mass flow rate is the upstream node, and vice versa. When the mass flow rate direction is negative, the enthalpy of the downstream node on the other side is taken as the enthalpy entering the branch. To avoid discontinuous jumps when the mass flow rate is close to 0 or the flow direction is reversed, a continuously differentiable smooth step function is used to smoothly transition the enthalpy values ​​between the two directions. In the formula, The value of the smooth step function that is continuously differentiable is the flow weight, and its value ranges from 0 to 1. This indicates the mass flow rate in the connector branch; This represents the smoothing parameter, a small value greater than 0, used to control the smoothing range of the function from 0 to 1; The flow enthalpy for the connector branch used in convective energy calculations; and These are the enthalpy of the upstream node and the enthalpy of the downstream node, respectively, with the upstream and downstream components determined by mass flow rates. The positive or negative sign of the variable is determined; this treatment ensures that the connector branch maintains equation continuity even when the flow direction changes; During the solution process at one time step, the system of differential-algebraic equations remains unchanged, and the activation matrix is ​​only used as a parameter to control whether the candidate connector branch transmits flux; only at the adaptive control time is the activation matrix updated according to criteria such as temperature gradient, heat flux density or mass flow rate.

5. The adaptive calculation method for differential algebraic network topology for multiphase thermofluid simulation as described in claim 5, characterized in that, Step 4 is as follows: After solving the fully implicit equations in step 3, if the equation residuals at the current time step meet the preset tolerance, the solution at the current time step is considered to have converged. Only when the solution has converged will the adaptive control time step be entered and the adaptive criterion be calculated. When calculating the adaptive criterion, the state variables of each adaptive master node and its adjacent nodes or adjacent connector branches are read from the current global solution vector. The state variables include the temperature, pressure, mass, enthalpy or internal energy of the master node, as well as the mass flow rate, heat flow rate and phase change mass flow rate in the connector branches connected to the master node. Subsequently, the local non-uniformity index is calculated for each adaptive master node. The local non-uniformity index includes at least the temperature gradient index, the heat flux density index and the mass flow rate index. The aforementioned feature indicators are divided by their respective preset feature scales for dimensionless normalization, and then weighted and summed using the weight coefficients of each indicator to generate a dynamic adaptive comprehensive criterion. ; In the formula, For adaptive comprehensive criteria; The total number of physical field characteristic indicators involved in the comprehensive criterion calculation; These are index subscripts for characteristic indicators, representing different physical quantities; For the first The weight coefficients corresponding to the physical field characteristic indicators satisfy the normalization condition. ; Represents the extracted first Key physical field characteristic quantities; Is with Preset feature reference scales or tolerance thresholds with the same dimensions are used to achieve dimensionless representation of each physical quantity; Adaptive Synthesis Criterion Used to update the activation matrix in subsequent steps; when When the refinement conditions are met, it indicates that there is strong local non-uniformity near the main node, and its corresponding candidate secondary node and candidate connector branch need to be activated; when When the condition is below the dormancy condition, it indicates that the state change in this area is weak, and the corresponding candidate secondary node will be put into dormancy.

6. The adaptive calculation method for differential algebraic network topology for multiphase thermofluid simulation as described in claim 1, characterized in that, In step 4, the preset characteristic dimensions include the characteristic temperature difference of the storage tank and the design heat flux density.

7. The adaptive calculation method for differential algebraic network topology for multiphase thermofluid simulation as described in claim 1, characterized in that, In step 4, the first Key physical field characteristic quantities Represents spatial gradient, physical flux, or mass / energy source term.

8. The adaptive calculation method for differential algebraic network topology for multiphase thermofluid simulation as described in claim 1, characterized in that, Step 5 is as follows: After obtaining the adaptive synthesis criterion in step 4, determine whether each master node needs refinement or hibernation; let the first... The adaptive comprehensive criterion for each master node is: Refine the threshold to The hibernation threshold is And satisfy: Among them, the refined threshold Used to determine whether a local area needs improved spatial resolution; sleep threshold. Used to determine whether an activated candidate secondary node can exit the current solution; an interval is set between the two to avoid the same primary node repeatedly activating and dormant near the critical state; the threshold is refined. Determined after normalization based on temperature gradient threshold, heat flux density threshold, and mass flow threshold; Adaptive comprehensive criteria for primary nodes and candidate secondary nodes The adaptive comprehensive criterion value of different nodes is related to the thermal state of their neighboring nodes and is not necessarily the same; the refinement threshold of the main node is higher than that of the secondary node to avoid large-area activation and control system scale and computational cost; the dormancy threshold of the secondary node is lower than that of the main node. The activation matrix is ​​updated according to the following rules: If a master node is not currently refined, and its adaptive synthesis criterion is satisfied: If the primary node is deemed to need further refinement, then, based on the primary-secondary mapping table established in step 2, the candidate secondary nodes and primary-secondary connector branches corresponding to the primary node are found, and the corresponding activation matrix elements are changed from 0 to 1. If a candidate secondary node under a certain primary node is already active, and its adaptive synthesis criterion satisfies: If the local changes in the region have been weakened, it is considered that the region can be put into hibernation. Similarly, the candidate secondary nodes and candidate connector branches corresponding to the primary node are found according to the primary-secondary mapping table, and the corresponding activation matrix elements are changed from 1 to 0. If the adaptive synthesis criterion satisfies: Then keep the current activation matrix elements unchanged; The adaptive topology state update involves changing the activation matrix elements established in step 2. When an activation matrix element is changed from 0 to 1, the corresponding candidate connector branch begins to transfer mass flow, momentum flow, and heat flow in the static superset graph network established in step 1. When an activation matrix element is changed from 1 to 0, the flux of the corresponding candidate connector branch becomes 0 again. Thus, the adaptive switching of the local topology state is achieved while keeping the global set of unknowns and the equation structure unchanged.

9. The adaptive calculation method for differential algebraic network topology for multiphase thermofluid simulation as described in claim 1, characterized in that, Step 6 is as follows: When candidate secondary nodes adaptive comprehensive criteria Greater than or equal to the refinement threshold At that time, conserved quantities such as mass and internal energy are conserved according to volume ratios: In the formula, Assign a ratio to the parameters. and Representing the master node respectively i and its corresponding candidate secondary nodes j The physical quantity; For non-conservative quantities such as temperature, pressure, velocity, and enthalpy, reconstruction is performed using the states of the master node and its neighboring nodes; for temperature, local linear reconstruction is used. In the formula, The initial temperature of the candidate secondary node; Master node temperature; temperature gradient The temperature difference between the left, right, top, and bottom adjacent nodes of the master node is calculated. The center position of the master node; master node No. j The center position of each candidate secondary node; Non-conserved quantities such as pressure and enthalpy are linearly reconstructed using the same method, or recalculated from the physical property equations after mass and internal energy are written in. When candidate secondary nodes adaptive comprehensive criteria Less than or equal to the sleep threshold When this happens, the additivity conservation quantities of the candidate secondary nodes are merged back into the corresponding primary node: In the formula, This represents the total number of candidate child nodes pre-allocated to a single master node, i.e., the number of split levels; the conserved quantity of a master node is equal to the sum of the conserved quantities of all directly candidate child nodes corresponding to that master node. This formula is only applicable to additively conserved quantities such as mass, internal energy, and total energy. For non-conserved quantities such as temperature, pressure, velocity, and enthalpy, they should be recalculated based on the combined mass and internal energy using property relations and equations of state. During the activation or dormancy of candidate secondary nodes, a small difference in mass or energy may occur due to variable reconstruction, state reinitialization, or numerical rounding. This difference is defined as the compensation amount. The compensation amount generated before and after the topology switch is divided by the time interval corresponding to this adaptive control, making it the source term per unit time. In the formula, The primary node compensation source item generated during topology state switching; This is the conserved quantity of the master node before the topology state switch; This refers to the conserved quantity after switching and completing the conservation mapping; For time step; Write the calculated physical quantity compensation source terms into the right-hand side of the corresponding conservation equation. For the mass and energy equation: in, Master node quality, The mass flow rate flowing into this master node, For the quality flow rate flowing out of the master node, Describe the internal energy of the master node, Enthalpy value, The power flowing into or out of the node; and These are the mass compensation source term and energy compensation source term generated by the topology state transition, respectively; if no difference is generated before and after the conservation mapping, the compensation source term is zero. After the above processing is completed, the solver continues to advance the time step using the updated activation matrix, candidate secondary node state values, and compensation source terms, and continues the fully implicit solution; the topology state evolution only changes the participation state of candidate secondary nodes and candidate connector branches through the activation matrix, and corrects the difference of conserved quantities through the compensation source terms; no global unknowns are added or deleted, and no connection relationships are re-established; since the candidate nodes and candidate connector branches in the static superset graph network have already entered the global equation during the initialization phase, the candidate sparse structure of the Jacobian matrix can be continuously reused, thereby completing the topology adaptive advancement without triggering the global equation reorganization.

Images