Dynamic Modeling and Full-Process Simulation of the Core Engine Test Rig Main Test System

: As a virtual digital model that can reﬂect physical entities or systems, digital twins are revolutionizing industry. The ﬁrst prerequisite for the construction of digital twins is the establishment of high-precision and complex entities or system models. A 47-components numerical system is established for the core engine test rig main test system by using the ﬁnite volume modularization modeling method. A comprehensive solution to the system-level valve-spool/oriﬁce throttling modeling, the key issue of the ﬂuid pipeline system modeling, is presented, and the algorithms of throttling and mixing are deepened and expanded. The full-process simulation study on two tests of normal-temperature 1400 s and low-temperature 1240 s shows that the combined regulation of ﬁve regulator valves and the change of cold source directly decide dynamic change of the system in each stage; the simulation reveals the phenomena such as the gas cylinder cooling with deﬂation, the air cooling when expanding from main pipeline to two branch pipelines, shunting ﬂow by branch pipeline, and the cold and hot gases mixing; the overall variation trends of the simulation curves are consistent with those of all the experimental curves of the test rig normal-temperature/low-temperature air supply lines, exhaust bypass, and engine main line in two operating conditions, and the maximum error between simulation curves and test curves of pressure, total pressure, and total temperature is less than 12%. The numerical system can be used for the construction of virtual models of digital twins, and the modeling method provides a feasible solution to the key technology of digital twins.


Introduction
With the advent of the industry 4.0 era, the concept of a digital twin is gradually being applied to aerospace, automotive, military, energy, and other industries [1].It is consensual that a digital twin is a virtual representation of a specific physical asset.A virtual representation means a machine learning model and/or a physics-based computational model [2].Modeling methods based on machine learning do not need to consider the physical laws specially and can be used to solve complex engineering problems when the data is sufficient.Min et al. [3] used machine learning-based LightGBM algorithms to build a digital twin model of a petrochemical factory, which provides a solution to control problems in petrochemical production.However, machine learning models are often limited to the training domain and cannot be used to design a new system.In addition, they usually need more data for the training (calibration) to compensate the absence of an initial mathematical structure (materialized as a physical law).Therefore, it is more appropriate to establish a digital twin model using a modeling method based on physical laws for engineering problems with limited amount of experimental data and focusing on physical changes during the entire work process of physical entities.Defraeye et al. [4] used COMSOL to establish a multiphysical simulation model of mangoes, as a digital twin of cold chain transportation of fruits, and studied the thermal behavior of mangoes in the entire cold chain.Vasileiou et al. [5] imported the passive biped robot structure drawn by CAD into the multibody software (MSC Adams) and established the multibody dynamics simulation model as the digital twin of the passive biped robot.Combined with the actual model of the passive biped robot, the design of the passive biped robot was completed.In the aerospace field, more attention is paid to the parameters and physical changes in the real-time working process of complex systems such as satellites and engine pipeline systems, and the establishment of complex simulation systems based on the physical laws is the key technology of digital twins.Shangguan et al. [6] used Modelica based on the object-oriented declarative modeling language to establish a multi-domain integrated satellite system as a digital twin.Combined with the actual satellite system, fault diagnosis and health monitoring were realized.
The modeling methods based on the laws of physics can be roughly divided into two types by application object: one is the high-dimensional modeling and simulation of a single component, such as the nozzle [7] in a liquid rocket engine system and the compressor [8] in an aircraft engine system, and this type of research generally uses computational fluid dynamics (CFD) software; the other is the simulation of the whole or partial system, and this type of research usually uses a one-dimensional (1D) or even zero-dimensional (0D) simplified model.It focuses on the overall characteristic of the system and the interrelationships among various components rather than the physical process within any single component.Common special or commercial software tools include GFSSP [9,10], Amesim [11,12], EcosimPro [13], Flowmaster [14] and so on.To meet the requirements of multipurpose general modeling, the system simulation usually follows the ideology of modularization, that is, to segment the system into several independent functional modules, and then combine them into a numerical model of the whole system according to certain rules.All the above-mentioned software for system simulation follow such modeling idea [15].However, since their components models come from different fields, the biggest problem in the basic simulation theory of the system-level model systems represented by these software is the absence of a system-level unified and comprehensive set of basic equations in integral or differential form, and the absence of a derivation process of component model from such an equation set in combination with the structure of typical component.To solve this problem, [16] proposes a simulation theory and model system for a flow, heat transfer, and combustion multi-field coupling pipe flow system.Its flow-field model can consider many factors, such as multispecies chemical reaction, species diffusion, variable volume, variable properties, fluid axial heat conduction, pipe-wall heat transfer, gravity field, and friction effects.Its temperature-field model can handle two types of boundary conditions, including convective and radiation heat transfer, and can calculate the heat transfer within variable-property pipe wall with multilayer or vacuum interlayer structure.Its chemical reaction kinetic model can determine the rate of each reaction by various reaction mechanisms and then obtain the production rates of various species.This system can also derive two sub-systems, i.e., a thermodynamic-calculation-based multi-field coupling pipe-flow model system and a flow/heat transfer two-field coupling pipe-flow model system without chemical reaction and species diffusion.Nevertheless, if modularization simulation models of the common components (pipe, chamber and branch, valve, etc.) and various complex components (combustion chamber, pressure reducing regulator, steady apparatus, tank, etc.) in the pipeline system need to be derived from the three systems, various reaction mechanisms need to be introduced into the components which consider the chemical kinetic effect, and suitable valve-spool/orifice throttling models need to be used for the components with throttling effect.The previous researches [15,16] commonly provide clear modeling solutions for the throttling problem of the valves with known information (i.e., structure and throttling characteristics) but seldom mentions the throttling modeling method for those valves without such information.This paper takes the modeling of a core engine test rig system with both known and unknown information valves as an opportunity to give a comprehensive solution to the key issue of the system-level spool/orifice throttling modeling.
As propellant feeding, gas pressurization and blowdown, and combustion-gas flow are often involved in various propulsion systems, these systems are typically complex pipe networks consisting of various pipelines.The components which adjust and control such parameters as flow rate and pressure of the liquid or gas in each pipeline are mainly various kinds of regulator valves.Therefore, the modeling accuracy of the regulator valve directly determines the accuracy of the flow-rate and pressure-drop distribution and finally determines the accuracy of the system simulation.The numerical research on regulator valves can also be roughly divided into two types: the first type [17][18][19][20][21][22][23] focuses on certain characteristics of the valve itself or effect of other factors on the valves, and usually uses CFD software to conduct high-dimensional modeling for the target valve and then makes numerical analysis.For example, Shin [17] used Fluent 6.1 to conduced 3D modeling and steady-state and transient-state computation for the pressure control valve and the pipes connected therewith in the natural gas pipeline, analyzed the transient flow characteristics of the control valve at the time of being closed, and found out the failure cause of the pressure control valve system and examined the freezing possibility of the connecting pipes in its closing process.By using Fluent 15.0 to conduct steady-state computation for the triple eccentric butterfly valve, Sun et al. [18] discovered that the surface roughness may result in a big difference between the flow coefficient obtained by simulation without friction and the actual flow coefficient.The second type of research [24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42] focuses on the whole system rather than the valve, and thus generally uses the traditional 0D throttling flow-rate equation instead of the momentum equation to describe the dynamic throttling effect of the valve spool.In order to accurately describe the throttling effect, this simplified modeling method must take flow-rate equation, flow coefficient calculation, and valve-spool opening and closing rule into full consideration.
Firstly, as the core of the throttling model, the valve-spool or orifice throttling flow-rate equation is a semi-theoretical and semi-empirical formula.By the form of relationship between the flow rate and the upstream and downstream pressures, it can be roughly divided into the pressure ratio-based form [15,24,25] and the pressure difference-based form [15,[27][28][29][30][31][32][33][34][35][40][41][42]. The former generally applies to the condition of gas flow medium while the latter usually applies to the condition of liquid flow medium.In addition, to avoid the flow-rate distortion of the foresaid approach in the case of slight differential pressure in front of and behind the valve spool, the local loss model [25] returning to the momentum equation may be used.Secondly, in the case that the pressures and densities at the upstream and downstream of the spool or orifice are known, the calculation value of the flow rate at the throttling point depends on values of the flow cross-sectional area and the flow coefficient, where the flow area is a one-variable function of the spool opening, and the relation between the flow area and the opening can be obtained by mathematical derivation once the shapes of the valve spool and valve seat are determined, whereas the flow coefficient is a function of many factors such as opening, fluid properties, and throttling form, and it is a very complicated challenge to determine the relation between the flow coefficient and the opening which is generally obtained by interpolating or fitting the experimental data of spool characteristics.By treatment method for the flow area and the flow coefficient, the second type of research can be further divided into four cases: (1) both are calculated by their respective relations with the opening; (2) both are integrated into one variable and deemed as the function of opening; (3) the flow coefficient is an empirical constant; (4) both flow coefficient and flow area are constant.For example, in [15], the pressure ratio-based injector orifice model uses the first type of method, while the pressure difference-based orifice model uses the second type of method.Liang et al. [42] investigated the temperature control of a vehicle climate chamber, and the flow coefficient of adjusting valves in the chilled water system also adopted the second type of method.As for the third type of method, Karimi et al. [26] used a nozzle model with variable area for the modeling of the pressure reducer valve when simulating the liquid propellant engine pressurization systems, where the flow area is function of the opening and the flow coefficient is constant; the same treating method also occurred in the modeling of the first and second throttling orifices of a class of three-way pressure reducing valve by Gad [27].The fourth type of method usually applies to the full open valve or the constant-flow-area orifice.When conducting the modeling for the three-way pressure-reducing valve, Gad [27] considered the flow coefficient as constant for the loading orifice with constant flow area.Huang et al. [28] also set the flow coefficient of the pilot metering orifice in the control valve as constant in their numerical research on a two-stage proportional flow control valve with pilot digital pressure compensator.Thirdly, as the opening or closing of the valve spool is not completed instantaneously, it is necessary to take into account the response time and rule of the spool when adjusting its opening.
The second type of simulation research covers both the information-known valves [15,25,[27][28][29][30][31][32][33][34][35][36] and the information-unknown valves [25,37,38,42], as well as the throttling devices, such as some orifices [39] or Venturi tube [40], and can be divided into the following three cases by the valve opening in the experimental process: (1) Zero opening, that is, the valve is in the closed position: for the valves with known or unknown information, the flow rate can be guaranteed at zero by the flow cross-sectional area or flow coefficient being zero, similar to the valves GV2 and GV6 in [25].
(2) Full opening, that is, the valve is in fully open position: for the valves with known information, the corresponding throttling model can be used for calculation; the valves with unknown information can be deemed as pipe or flow resistance component.For example, although a dual-model scheme, which can automatically switch according to the pressure ratio, is set for GV5, a clear-way valve, in [25], the local loss model was actually used for the dynamic simulation of 260 s due to the small pressure loss under the test conditions; for another example, Li et al. [37] used the simple pressure-drop equation for calculation of the valve when it is fully opened in the steady-state simulation of the natural gas pressure-relieving system.
(3) Variable opening, that is, the valve has variable opening in the test.In this case, there are usually three options for valves with known information or known structure: (1) the flow coefficient of the valve is set as a constant.For example, Gad [31] set the flow coefficient of the variable throttling area as a constant when conducting modeling for the pilot operated relief valve; (2) the flow coefficient is calculated by the fitting formula of the relationship between it and opening.For example, for the emergency shutdown valves of the high-pressure carbon dioxide transportation pipelines, Mahgerefteh et al. [29], in view of the structure character of the ball valve and based on the pressure difference-based Wylie flow-rate equation, used the fitting formula based on Wylie ball valve data for the relationship between flow coefficient and opening; (3) in the case of small opening or low Reynolds number, the valve flow coefficient is calculated by the relation formula obtained by experiment between it and the Reynolds number.For example, in the modeling of the disc style check valve conducted by Knutson et al. [34], the flow coefficient is calculated by the relation with the Reynolds number corrected by test data; Gad [27] also used the relation formula between flow coefficient and Reynolds number in the modeling of the damping orifice of the three-way pressure reducing valve.For valves with structure and throttling character unknown, Li et al. [37] provided the Perkins model used for pressure drop calculation of the choke valve in a natural gas pressure-relieving system, but in fact the test data were directly used for determining the conditions behind and in front of the valve under the steady-state condition in the steady-state simulation, and the pipeline flow rate was calculated by the pipe flow equation in combination with pressure distribution data obtained by test rather than by valve flow-rate equation, which is not applicable to the dynamic simulation of the variable-opening valve system.For the water distribution systems/networks, Piller et al. [38] used the constrained optimization method based on pipe-network test data to deal with the valve throttling so as to provide the pressure distribution of the pipe network, but it is not applicable to ordinary variable-opening dynamic adjustment simulation.
The two most commonly used flow characteristics of the regulator valves are linear and equal percentage [41,43], and the variation trends of the flow characteristic curves of the valves with equal percentage flow characteristic are basically the same.This means the flow characteristic curves of most regulator valves are between the equal percentage curve [15,25,35,[41][42][43] and the linear curve [41,43,44].Based on this conclusion, this paper proposes a set of modeling methods for the information-unknown valves by using information-known reference valve model and limited system test data and applies the methods to the modeling of a core engine test rig main test system consisting of a number of regulator valves.Meanwhile, the corresponding valve-spool throttling scheme is used for the modeling of the valves with known information.The gas steady apparatus can be regarded as a valve with many throttling orifices located on the orifice plate, so it is also a throttling component.The model of gas mixing apparatus is often studied on component level and uses Fick's law, and it generally does not consider wall heat transfer in the process of gas mixing [45,46].The modularization simulation models are derived from the flow/heat transfer two-field coupling pipe-flow model subsystem [16,24] in combination with the typical component structure.The effectiveness of the model system is verified by comparing the simulation results of the whole process of the system in two operating conditions with the test results, so it can be used as a building block for a so-called digital twin.

Physical Model of the System
Figure 1 shows the schematic diagram of the core engine test rig system, which can be divided into three parts: the main test system, the fuel supply and ignition system, and the engine system.The main test system is to supply gas to the test engine and can be roughly divided into three lines by test operating conditions: normal-temperature line, low-temperature line, and high-temperature line, which are, respectively, used for environmental simulation tests of different operating conditions; the fuel supply and ignition system is to supply fuel to the engine and ignite the aviation kerosene and air mixture to enable burning.distribution of the pipe network, but it is not applicable to ordinary variable-opening dynamic adjustment simulation.
The two most commonly used flow characteristics of the regulator valves are linear and equal percentage [41,43], and the variation trends of the flow characteristic curves of the valves with equal percentage flow characteristic are basically the same.This means the flow characteristic curves of most regulator valves are between the equal percentage curve [15,25,35,[41][42][43] and the linear curve [41,43,44].Based on this conclusion, this paper proposes a set of modeling methods for the information-unknown valves by using information-known reference valve model and limited system test data and applies the methods to the modeling of a core engine test rig main test system consisting of a number of regulator valves.Meanwhile, the corresponding valve-spool throttling scheme is used for the modeling of the valves with known information.The gas steady apparatus can be regarded as a valve with many throttling orifices located on the orifice plate, so it is also a throttling component.The model of gas mixing apparatus is often studied on component level and uses Fick's law, and it generally does not consider wall heat transfer in the process of gas mixing [45,46].The modularization simulation models are derived from the flow/heat transfer two-field coupling pipe-flow model subsystem [16,24] in combination with the typical component structure.The effectiveness of the model system is verified by comparing the simulation results of the whole process of the system in two operating conditions with the test results, so it can be used as a building block for a so-called digital twin.

Physical Model of the System
Figure 1 shows the schematic diagram of the core engine test rig system, which can be divided into three parts: the main test system, the fuel supply and ignition system, and the engine system.The main test system is to supply gas to the test engine and can be roughly divided into three lines by test operating conditions: normal-temperature line, low-temperature line, and high-temperature line, which are, respectively, used for environmental simulation tests of different operating conditions; the fuel supply and ignition system is to supply fuel to the engine and ignite the aviation kerosene and air mixture to enable burning.The research object of this paper is the main test system.In the normal-temperaturecondition test, the regulator valve F3 is used to regulate the total air flow rate, F9 and F10 are used separately to control the air flow into the engine main line and the exhaust bypass, and F4, F6, F7, and F11 are closed all the time; in the low-temperature condition, F4  The research object of this paper is the main test system.In the normal-temperaturecondition test, the regulator valve F3 is used to regulate the total air flow rate, F9 and F10 are used separately to control the air flow into the engine main line and the exhaust bypass, and F4, F6, F7, and F11 are closed all the time; in the low-temperature condition, F4 is used to regulate the normal-temperature air flow rate, the low-temperature air flow rate is controlled by the low-temperature air source and its corresponding regulator valve (F11 only plays an auxiliary regulation role) which is not involved in this test rig, the functions of F9 and F10 are the same as those in the normal-temperature condition, and F3 and F8 remain closed.The ultimate goal of the two operating conditions is to control the total pressure p tm, total temperature T tm , and the flow rate Q ma of the engine inlet gas by combined adjusting various valves so as to provide the required air working medium for the engine test.
It can be seen from the physical model and working characteristics of the system that the gas flow pipeline of the test rig is complex in constitution, there are numerous regulating devices, and the intake state of the engine is controlled by the combined regulation of multiple valves.As different operating conditions have different requirements on the engine inlet air, there are big differences in the regulating scheme of various valves, the combination of multiple valves makes the time sequence design and regulation more complex, and these valves are regulated manually based on experience and limited technical specifications and features.All these factors together make the current debugging of the test rig an arduous task for engineers.By contrast, the computer numerical simulation method can be used to conduct numerical modeling for the rig-regulating process without starting up the real hardware to study the effect of the opening and closing time sequence of the valves on the system dynamic characteristics.Then, an appropriate regulating scheme can be concluded for a specific operating condition.

Modularization Modeling Approach
The typical components of the main test system include fluid source, gas pipe, gas volume, five types of regulator valves, mixer, filter, stabilizer, and heater.As the structure of the filter is similar to that of the stabilizer, the filter is deemed as stabilizer.Table 1 shows the module information of seven typical components generated by modularization disassembly of the main test system, of which the fluid source, gas pipe, gas volume, gas mixing apparatus, gas valve are five basic modules, based on which gas steady apparatus is obtained by combination development; the heater exchange cold gas (also can be called heater) has an heat-transfer increasing-temperature effect only in the high-temperature operating condition and is equivalent to the gas pipe module component under the normaland low-temperature operating conditions.The derivation source of modularization models of all these typical components is a flow/heat transfer two-field coupling subsystem derived from a multi-field coupling pipe-flow system simulation theory and model system [16], including the finite volume model of quasi-one-dimensional compressible variable-section pipe flow transient flow field and the finite volume model of pipe-wall transient heat transfer.Below gives the model derivation method for three modules, GV, GSA, and GMA.

Finite Volume Model
For the components with small length-to-diameter ratio and sudden change in the cross section in the fluid flow direction, such as stop valve, ball valve, compensator, charge/discharge valve, safety valve, and orifice plate, the key issue to be addressed in the modeling is how to describe the throttling and switching characteristic of the valve spool.
As shown in Figure 2, the gas valve can be regarded as a combined module in which two gas volume units are connected by local flow resistance.

Finite Volume Model
For the components with small length-to-diameter ratio and sudden change in t cross section in the fluid flow direction, such as stop valve, ball valve, compensat charge/discharge valve, safety valve, and orifice plate, the key issue to be addressed in t modeling is how to describe the throttling and switching characteristic of the valve spo As shown in Figure 2, the gas valve can be regarded as a combined module in which tw gas volume units are connected by local flow resistance.Assumptions: (1) The gas is 1D ideal gas flow without considering the effects of axial heat condu tion and gravity field.
(2) The state parameters, such as pressure, density, and total energy per unit volu in the gas volume unit, are instantaneously uniform and homogeneous, and the friction and local losses are ignored.
(3) The wall is a rigid wall, irrespective of the wall elastic deformation.
(4) The valve-spool throttling can be described by the valve mass-flow-rate algebr equation (actually can be used only when (5) The local loss of the spool throttling can be described by introducing a flow co ficient.
As shown in Figure 2, the half elements of the upstream and downstream of the spo are considered as two velocity elements, 1 u and 2 u , and their model equations are follows: Continuity equation: For the wall of the gas valve inlet channel (inlet chamber) and outlet channel (out chamber), one of the radial one-dimensional heat transfer model, the zero-dimension heat transfer model, and the adiabatic model [15,16] can be selected based on the actu condition of the temperature field.

PipeIn PipeOut
Finite control volume grids of gas valve. Assumptions: (1) The gas is 1D ideal gas flow without considering the effects of axial heat conduction and gravity field.
(2) The state parameters, such as pressure, density, and total energy per unit volume in the gas volume unit, are instantaneously uniform and homogeneous, and the frictional and local losses are ignored.
(3) The wall is a rigid wall, irrespective of the wall elastic deformation.
(4) The valve-spool throttling can be described by the valve mass-flow-rate algebraic equation Q m = f (τ, p 1 , ρ 1 , p 2 , ρ 2 ) (actually can be used only when τA vt A v1 , A v2 ).( 5) The local loss of the spool throttling can be described by introducing a flow coefficient.As shown in Figure 2, the half elements of the upstream and downstream of the spool are considered as two velocity elements, u 1 and u 2 , and their model equations are as follows: Continuity equation: Energy equation: For the wall of the gas valve inlet channel (inlet chamber) and outlet channel (outlet chamber), one of the radial one-dimensional heat transfer model, the zero-dimensional heat transfer model, and the adiabatic model [15,16] can be selected based on the actual condition of the temperature field.

Valve-Spool Throttling Models
The injector orifice model of IEC60534-1998 standard-based expansion coefficient method is established for Vanessa 30,000 pneumatic zero leakage butterfly valve.The spool mass-flow-rate equation as follows is obtained by converting to the international unit and converting volume flow rate into mass flow rate: In the case of reverse flow, that is, when the equation meets 0 ≤ p 1 < p 2 , the positions of Var 1 and Var 2 are exchanged and a minus sign is added in front of the mass flow-rate result, where Var ∈ {ρ, p, T}, F k = K/1.4,Y ef = 1 − X/(3F k X t ).The relation formula of X t and τ (θ = 90τ is butterfly valve-spool rotation angle) is: Besides 0~7 eight calculation schemes of flow coefficient [14,24], the ninth scheme can be obtained by fitting the flow-rate characteristic curve of the Vanessa 30,000 series valves provided by the manufacturer: K vPercent = −76.58189τThe error sum of squares between the calculated results of the flow coefficient percentage K vPercent = K v /K vRating and the characteristic data is 5.4356 × 10 −4 , and the standard error of the estimate (SEE) is 2.0141 × 10 −3 , where K vRating is the flow coefficient rating, that is, the value of K v in the case of τ = 1.In case of small opening, the linear interpolation of two points, K v (τ) = τ • K v (τ subcr )/τ subcr , τ subcr = 0.01, is used to supplement the formula near zero-opening zone and ensure Q m = 0 when τ = 0.

The Treatment Scheme for Valves with Unknown Information
The throttling model in Section 3.1.2is applicable to the valves with known structure and throttling characteristics, but there are some valves in the test rig system whose throttling flow-rate equation and flow characteristic curve are unavailable.Some of them directly decide the flow rate of the branches where they are located.For such valves with unknown information (hereinafter referred to as the undetermined valves), a modeling method which uses system test data to correct the spool throttling model of the reference valve is proposed.According to the condition of experimental data and the complexity of modeling, it can be divided into three options.Option A and Option B, applicable to the case with limited test data, can ensure that the flow-rate simulation result is the same as the test result at the calibration opening of the valve.Option A is applicable to the case where the test data on the upstream and downstream of the undetermined valves are available; Option B is applicable to the case where there is lack of upstream and downstream test data of the undetermined valves, but the system test data are available.
Option A: (1) The rated (i.e., when τ = 1) flow coefficient of the reference valve model is calibrated based on the rig-regulating test data: with the test data of the upstream and downstream of the undetermined valve (only the upstream data in the case of supercritical throttling state) as known data, the rated flow coefficient corresponding to the undetermined valve is obtained through inverse calculation of throttling flow-rate equation and flow characteristic fitting formula of the reference valve.Although the calibration should use such test data as flow rate and pressure of the undetermined valve in its fully open position, there is often only a limited opening for the valve in the test process.In this situation, the test data at its maximum opening can be used.(2) Since the experimental data for inverse calculation is not necessarily obtained by the upstream and downstream measurement points near the undetermined valve (they may be far apart), an adjusting step needs to be conducted for the calibration result, i.e., the calibrated rated flow coefficient is finely adjusted according to the difference between the full-system simulation results and the test data.The adjusting criterion is to match the calculated flow rate with the test result.Figure 3 shows the calibration and adjustment cases of undetermined valve F3 when selecting Vanessa 30,000 pneumatic zero-leakage butterfly valve and HCB cage-type double-seat high-capacity sleeve regulator valve as the reference valve, respectively.step needs to be conducted for the calibration result, i.e., the calibrated rated flow coefficient is finely adjusted according to the difference between the full-system simulation results and the test data.The adjusting criterion is to match the calculated flow rate with the test result.Figure 3 shows the calibration and adjustment cases of undetermined valve F3 when selecting Vanessa 30,000 pneumatic zero-leakage butterfly valve and HCB cagetype double-seat high-capacity sleeve regulator valve as the reference valve, respectively.Option B: (1) Full-system simulation trial calculation is conducted with the undetermined valves substituted by the reference-valve model.(2) According to the deviation between the simulation results and the test data of test-rig relevant part, the rated flow coefficient corresponding to the undetermined valve in the reference-valve model can be calibrated through finite number of simulation trials by the dichotomy method.This option is used for the undetermined valve F4 in the test rig.Since there is no flow-rate measurement point in the normal-temperature and the low-temperature air supply lines, the calibrating criterion of the F4 rated flow coefficient is to match the overall amplitude of the total-temperature simulation curves of exhaust bypass and engine inlet with the test curves.
Option C: (1) With the accumulation of test data under multioperating conditions, system experimental data with multiple openings gradually become available for the undetermined valve.These test data can be regarded as known quantities to obtain the flow coefficients corresponding to various openings of the undetermined valve through inverse calculation of throttling flow-rate equation of the reference valve.Option B: (1) Full-system simulation trial calculation is conducted with the undetermined valves substituted by the reference-valve model.(2) According to the deviation between the simulation results and the test data of test-rig relevant part, the rated flow coefficient corresponding to the undetermined valve in the reference-valve model can be calibrated through finite number of simulation trials by the dichotomy method.This option is used for the undetermined valve F4 in the test rig.Since there is no flow-rate measurement point in the normal-temperature and the low-temperature air supply lines, the calibrating criterion of the F4 rated flow coefficient is to match the overall amplitude of the total-temperature simulation curves of exhaust bypass and engine inlet with the test curves.
Option C: (1) With the accumulation of test data under multioperating conditions, system experimental data with multiple openings gradually become available for the undetermined valve.These test data can be regarded as known quantities to obtain the flow coefficients corresponding to various openings of the undetermined valve through inverse calculation of throttling flow-rate equation of the reference valve.(2) The relational expression between flow coefficient and opening is obtained by fitting, i.e., the flow characteristic fitting formula.The combination of it and the reference valve throttling flow-rate equation are the throttling model of the undetermined valve.

Gas Steady Apparatus Model
Figure 4 shows the finite control volume grids of gas steady apparatus.The threechambers and two-orifices structure and the two-chambers and one-orifice structure are considered compatibly in the modeling and programming, where the flow-field grid number variable is used to recognize the specific structure.Take the more complex threechambers and two-orifices structure in Figure 4a as an example, its boundaries are boundary grids of connecting gas pipes; the gas steady apparatus is regarded to be made up of three gas volumes including high-pressure, middle-, and low-pressure chambers, which are connected by the orifice plate with multiple orifices (local flow resistances assembly).Volumes of three chambers are considered as constant.
considered compatibly in the modeling and programming, where the flow-field grid number variable is used to recognize the specific structure.Take the more complex threechambers and two-orifices structure in Figure 4a as an example, its boundaries are boundary grids of connecting gas pipes; the gas steady apparatus is regarded to be made up of three gas volumes including high-pressure, middle-, and low-pressure chambers, which are connected by the orifice plate with multiple orifices (local flow resistances assembly).Volumes of three chambers are considered as constant.It is further assumed, based on the gas valve model assumptions, that without considering the difference in flow characteristics caused by the difference in the radial position of the orifices on the orifice plate, the flow coefficients are the same if the orifice diameters are the same.
Energy equation: Gas pipe inlet

Middle chamber
The last state element of gas pipe The first state element of gas pipe High-pressure chamber Low-pressure chamber Gas pipe inlet

Gas pipe outlet
The last state element of gas pipe The first state element of gas pipe High-pressure chamber Low-pressure chamber It is further assumed, based on the gas valve model assumptions, that without considering the difference in flow characteristics caused by the difference in the radial position of the orifices on the orifice plate, the flow coefficients are the same if the orifice diameters are the same.
Energy equation: Temperature calculation formulas:

The Throttling Points of the First and Second Orifice Plates
Mass flow-rate equation of a single throttling orifice: where the flow sectional area of each orifice on the ith layer orifice plate A oi = πd oi 2 /4, and n oi is the number of throttling orifices on the ith layer orifice plate; C di and C di are flow coefficients which represent orifice throttling characteristic, and the flow coefficient curves need to be established for different shapes of throttling orifices by test and appropriately amended for different fluid media according to viscosity and density.In the case of reverse flow, that is, when the equation meets 0 ≤ p i < p i+1 , the positions of Var i and Var i+1 are exchanged and a minus sign is added in front of the mass flow-rate result, where Var ∈ {ρ, p}.
Local velocity calculation formula at specific cross section:

Wall Heat Transfer Model
One of the radial 1D heat transfer model, the 0D heat transfer model, and the adiabatic model can be selected in line with the actual condition of temperature field.

Gas Mixing Apparatus Model
Gas mixing is a common phenomenon in nature.The mixing calculation needs to consider mass, momentum, and energy transport of various species in mixed gas and calculate the physical properties of various species.Figure 5 gives the finite control volume grids of gas mixing apparatus with three inlets and three outlets.Its boundaries are boundary grids of connecting components, i.e., its boundary is made up of half a velocity element extended from each gas pipe boundary.In the state element of gas mixing apparatus, the state parameters, such as pressure, density, and internal energy per unit mass, are instantaneously consistent and uniform.The frictional and local losses and the effects of axial heat conduction and gravity field are neglected.Assumptions: (1) The various gases flowing into the mixing apparatus are instantaneously and evenly mixed, and the mixing process is adiabatic.
(2) The incoming gases and mixtures are ideal gases.
(3) If the inflow is combustion gas, it is considered as completely frozen flow with State element of Gas Mixing Apparatus The first state element of outlet gas pipe 1 The last state element of inlet gas pipe 1

PipeOut1 PipeIn1
The last state element of inlet gas pipe 2 The first state element of outlet gas pipe 2

Pipe Out3
The last state element of inlet gas pipe 3 The first state element of outlet gas pipe 3 Assumptions: (1) The various gases flowing into the mixing apparatus are instantaneously and evenly mixed, and the mixing process is adiabatic.
(2) The incoming gases and mixtures are ideal gases.
(3) If the inflow is combustion gas, it is considered as completely frozen flow with species and physical properties frozen unchanged.

Flow Model
Continuity equation of species: where fun Energy equation: Mixture fraction equation: Partial density calculation formula of species: Partial mole number calculation formula of species: Total mass calculation formula of mixture: Total mole number calculation formula of mixture: Total density calculation formula of mixture: Mass flow-rate equation: Velocity (mean velocity) calculation formula: Pressure calculation formula of mixture: Temperature calculation formula of mixture: In the equations, ρ M, i and m M, i (i = 1, 2, • • • , n in ) are partial density and mass of various gas species involved in mixing, respectively.When calculating the mass and energy exchange among the state elements of gas mixing apparatus and connected pipes, the upwind scheme is used to ensure compatibility of the state equations, that is The adiabatic model is used for the wall heat transfer.

Thermophysical Properties Mixing Model
(1) Mean molar mass of mixture: (2) Mean gas constant of mixture: (3) Frozen specific heat at constant pressure of mixture: Mass fraction of the ith gas species involved in mixing: Frozen specific heat at constant volume of mixture: (4) Frozen sound velocity of mixture: (5) Viscosity coefficient of mixture µ Method of Wilke [47]: Mole fraction of the ith gas species involved in mixing: (6) Thermal conductivity of mixture:

System Model
The core engine test rig main test system is complex in constitution, there are physical phenomena such as flow, heat transfer, and mixing, and the combined regulation by multiple regulator valves is needed to meet the requirements of air total temperature, total pressure, and flow rate under different operating conditions.Depending on the system operating conditions, there are two modeling options available: (1) Establish a numerical model suitable for a specific operating condition only in consideration of the actually circulating pipelines and ignoring the non-circulating pipelines under such operating condition.
(2) Establish a complete simulation model based on the actual configuration of the system, in which different regulation time sequence, initial, and boundary conditions are set for different operating conditions.
The second option has the most components.Although the near-zero flow-rate noncirculating pipeline section seems to bring unnecessary calculation workload, the second option is most consistent with the actual debugging of the test rig and is more versatile and general.Therefore, this option is chosen for the modeling.Figure 6 shows the numerical simulation model of the core engine test rig main test system established strictly according to the structural size parameters given by the drawing.In the modularization modeling, the system is divided into three fluid sources (FS1-3), twenty-four gas pipes (GP1-24), six gas volumes (GVol1-6), one gas mixing apparatus (GMA1), two gas steady apparatus (GSA1-2), ten gas valves (GV1-10), and one heater exchange cold gas (HECG1) so as to form a 47-components numerical system.option is most consistent with the actual debugging of the test rig and is more versatile and general.Therefore, this option is chosen for the modeling.Figure 6 shows the numerical simulation model of the core engine test rig main test system established strictly according to the structural size parameters given by the drawing.In the modularization modeling, the system is divided into three fluid sources (FS1-3), twenty-four gas pipes (GP1-24), six gas volumes (GVol1-6), one gas mixing apparatus (GMA1), two gas steady apparatus (GSA1-2), ten gas valves (GV1-10), and one heater exchange cold gas (HECG1) so as to form a 47-component numerical system.Numerical simulations are carried out for two test-rig regulating processes in a 0~1400 s normal-temperature operating condition and a 0~1240 s low-temperature operating condition, respectively.In Figure 6, the length (unit: m), the outer diameter, and the wall thickness (unit: mm) of each pipe are given; the flow-field grid of the pipe-type com-   Numerical simulations are carried out for two test-rig regulating processes in a 0~1400 s normal-temperature operating condition and a 0~1240 s low-temperature operating condition, respectively.In Figure 6, the length (unit: m), the outer diameter, and the wall thickness (unit: mm) of each pipe are given; the flow-field grid of the pipe-type component is divided into 100 mm/grid; the length of the lumped parameter component (such as gas valve, gas volume, etc.) in the direction of pipeline is two standard grids, i.e., 200 mm; and the total volume of the air source (composed of many high-pressure cylinders) is set according to specific operating condition, i.e., 10,000 and 500 m 3 for normal-temperature condition and low-temperature condition, respectively.The simulation gives initial flow field and pipe-wall temperature field distribution which are close to the actual situation, and the test data are used for FS1-3 three fluid source boundary conditions.
The pipe-wall radial 1D heat transfer model is used for the GP, GV, and GSA module components; the wall 0D heat transfer model is used for the GVol module components; and the wall adiabatic treatment is given for the GMA module component.The total number of wall radial grid serial number of the components using the 1D heat transfer model is four.The heat exchange between the exterior walls of various components and the outer environment is considered as simple superposition of convection heat transfer and radiation heat transfer, where the convective heat transfer coefficient is set at 6 W/(m 2 K) for the sections of normal-temperature air, at 10 W/(m 2 K) for the section behind GP21, and at 0.1 W/(m 2 K) for the section with coating layer.The outer environment temperature is 291.65 K in normal-temperature operating condition and 289.25 K in low-temperature operating condition.
The ten regulator valves are modeled as gas valve module components, and by a large amount of simulation debugging, the throttling calculation schemes are set, as shown in Table 2.In the table, the flow-rate characteristic fitting formula scheme of Vanessa 30,000 pneumatic zero-leakage butterfly valve is used as the flow coefficient calculation scheme for GV2 (F11), GV9 (F9), and GV10 (F10), and the three valves correspondingly adopt the injector orifice model scheme of IEC60534-1998 standard-based expansion coefficient method; with HCB cage-type double-seat high-capacity sleeve regulator valve [25] as reference valve, GV3 (F3) and GV4 (F4) use the equal percentage flow character fitting formula scheme for which the rated flow coefficients are recalibrated, respectively, and correspondingly adopt the pressure difference-based injector orifice model scheme of Wuzhong Instrument Co., Ltd.[25].GV6 (F6) and GV7 (F7) remain fully closed in the two tests, so it is not important which calculation scheme is used as long as the flow rate calculated by the scheme in the zero-opening position is zero; the relative openings of GV1 (F1) in the two test operating conditions are all constant, and GV5 (F5) and GV8 (F8) are fully open in the normal-temperature condition and fully closed in the low-temperature condition.These five valves do not control the flow rate, and their throttling states do not affect the variation trend of the simulation curves.The HCB model is used for them in accordance with GV3, and their rated flow coefficients, by reference to the calibrationadjusting GV3, are obtained based on the principle that the ratio of rated flow coefficients is equal to the ratio of flow sectional areas.The flow coefficients at GSA1 and GSA2 orifice plate throttling points are set at 0.3.The dynamic process of the system is simulated by using the classical four-stage four-order Runge-Kutta method.

Regulating Time Sequence and Boundary Conditions
Figures 7 and 8 show the test data (time parameter) of the time sequence of key regulator valves and the boundary conditions of fluid source FS1-3 in two operating conditions.The acquisition frequency of the test data is 1 Hz, that is, one data point per second.FS1-3 are the pressure + temperature inlet, the total pressure outlet, and the pressure outlet fluid sources, respectively.FS1 uses the test data of the atmosphere pressure p gs and temperature T gs as boundary conditions in the normal-temperature condition due to closure of the cold air line, and the test data of the cold air inlet pressure p f11q and temperature T f11q (obtained by decreasing the test curve of the exhaust bypass total temperature T tb by 3K through verification of a large number of examples due to failure to install temperature measurement point as designed) as boundary conditions in the low-temperature condition; FS2 uses the test data of the engine upstream total pressure p tm and total temperature T tm as boundary conditions; as the outlet of the exhaust bypass is the atmosphere environment, the test data of the atmosphere pressure p gs and temperature T gs are used as boundary conditions for FS3.Ttm as boundary conditions; as the outlet of the exhaust bypass is the atmosphere environment, the test data of the atmosphere pressure pgs and temperature Tgs are used as boundary conditions for FS3.

Results and Discussion
Figures 9a-f and 10a-f are the comparisons between the simulation results of normaltemperature and low-temperature operating conditions and the pressure, flow rate, total temperature, and total pressure curves actually measured by the main test system, of which the experimental flow rates are obtained indirectly by the total static pressure method with the total pressure and pressure measured by the sensors.The simulation data output frequency in the figures is 10 Hz.The subheading of each figure is named after the corresponding measurement point.See Figure 1 for the location of each measurement point, and the numerical simulation model of the system shown in Figure 6 for the location of each simulation data point.In the figures, the parameter curves not marked by symbol are test data, while the parameter curves marked by triangle, circle, and other symbols are the simulation results at the component grids corresponding to the locations of the sensor measurement points, in which the 0, n/2, and n − 1 along flow direction are

Results and Discussion
Figures 9a-f and 10a-f are the comparisons between the simulation results of normaltemperature and low-temperature operating conditions and the pressure, flow rate, total temperature, and total pressure curves actually measured by the main test system, of which the experimental flow rates are obtained indirectly by the total static pressure method with the total pressure and pressure measured by the sensors.The simulation data output frequency in the figures is 10 Hz.The subheading of each figure is named after the corresponding measurement point.See Figure 1 for the location of each measurement point, and the numerical simulation model of the system shown in Figure 6 for the location of each simulation data point.In the figures, the parameter curves not marked by symbol are test data, while the parameter curves marked by triangle, circle, and other symbols are the simulation results at the component grids corresponding to the locations of the sensor measurement points, in which the 0, n/2, and n − 1 along flow direction are the initial, middle, and last grid of the pipe flow field, respectively.Below is a comparative analysis of the simulation and test curves in the two operating conditions.

Normal-Temperature Operating Condition
As shown in Figure 9a, the variation trend of the pressure simulation curve at GP1 middle grid is consistent with that of the air supply pressure test curve in the whole process of the normal-temperature operating condition; both curves move down with the gradual consumption of air in the high-pressure gas cylinders, except that the test curve moves down in staircase form and by an amplitude smaller than the simulation curve.The staircase characteristic exists because that the measurement accuracy at p a measuring point is only 0.01 MPa, and the pressure change lower than 0.01 MPa cannot be shown.The main reason why the simulation curve decreases slightly greater than the test curve is that the setting (such as volume) of the gas cylinders GVol1 is inconsistent with the actual situation of the air cylinder set of the gas source field which supplies air to this test.As there are a large number of high-pressure gas cylinders in the gas source field which may supply gases to several test rigs at the same time, and there is air charging operation at some certain moments, it is difficult to set the normalized GVol1 exactly the same as the actual situation.However, at 1400 s, the error between simulation and test reaches the maximum, which is still smaller than 4%.
is that the setting (such as volume) of the gas cylinders GVol1 is inconsistent with the actual situation of the air cylinder set of the gas source field which supplies air to this test.As there are a large number of high-pressure gas cylinders in the gas source field which may supply gases to several test rigs at the same time, and there is air charging operation at some certain moments, it is difficult to set the normalized GVol1 exactly the same as the actual situation.However, at 1400 s, the error between simulation and test reaches the maximum, which is still smaller than 4%.As shown in Figure 9b, the variation trend of the flow-rate simulation curve at GP22 outlet grid is consistent with that of the flow-rate test curve at the engine inlet except that the simulation curve is lower than the test curve on the whole.As indicated in Figure 7a, 16~35 s, the relative opening of F3, the main control valve of the normal-temperature air supply line, gradually increases to 0.44 from approximately zero and remains unchanged; As shown in Figure 10a, similar to the normal-temperature operating condition, the variation trend of the pressure simulation curve at GP1 middle grid is consistent with that of the normal-temperature air supply pressure test curve in the 1240 s rig-regulating process of the low-temperature operating condition, and the maximum error is less than 3%, except that the test curve with only 0.01 MPa measuring accuracy moves down in staircase form and by an amplitude a little bigger than the simulation curve.In addition, the test appears to have air replenishing operation around 1050 s.As shown in Figure 10b, the variation trend of the flow-rate simulation curve at GP22 outlet grid is consistent with that of the flow test curve at the engine inlet, and the mean error between simulation and test is less than 10%.In addition, as in the normal-temperature operating condition, the low-temperature operating condition also has the problem As shown in Figure 9b, the variation trend of the flow-rate simulation curve at GP22 outlet grid is consistent with that of the flow-rate test curve at the engine inlet except that the simulation curve is lower than the test curve on the whole.As indicated in Figure 7a, 16~35 s, the relative opening of F3, the main control valve of the normaltemperature air supply line, gradually increases to 0.44 from approximately zero and remains unchanged; approximately in the first 60 s, F10 is fully open while the relative opening of F9 is very small, only 0.0027.Therefore, in the first 60 s, the normal-temperature air supplied by the high-pressure gas cylinders flows through the normal-temperature air supply line mainly controlled by F3 and is discharged from the exhaust bypass mainly controlled by F10, and almost does not enter the engine main line.As a result, as shown in Figure 9b, GP22 flow-rate simulation result is close to zero in the first 60 s while the test result maintained at 0.7 kg/s is obviously inconsistent with the physical reality.During the subsequent period of 60~1400 s, the engine inlet flow rate changes under combined action of the dynamic regulation of F9, F10, and F3 in Figure 7a and the pressure decreasing resulting from the air consumption of gas cylinders in Figure 9a.The relative opening of F3 changes slightly while the relative openings of F9 and F10 change suddenly at 578 s, 758 s, 937 s, and 1060 s, and the simulation curve has obvious peaks at these four moments while the peak of the test curve is less obvious, which, together with the inaccuracy of the test results in the first 60 s, indicates that the flow-rate result measured by the total static pressure method is not necessarily accurate.
As shown in Figure 9c, the variation trend of the simulation curves of temperature, total temperature, and wall temperature at GP22 outlet grid is consistent with that of the total temperature test curve at the engine inlet, and with the opening of F9 at 60 s, these curves keep moving down until the end of this operating condition at 1400 s due to the deflation effect, which can also be seen from the decrease of the GP1 temperature curve in Figure 9a.The maximum error between the simulation result and the test result of GP22 total temperature is less than 3%.
As shown in Figure 9d, the variation trend of the total pressure simulation curves of GP23 initial and middle grids are consistent with that of the total pressure test curve in the front section of the exhaust bypass, and they are only different locally.As F10 is fully open in the first 60 s, the total pressures in the exhaust-bypass front section gradually increase from the initial value close to the atmosphere pressure with gradual opening of F3 during 16~35 s, reach the maximum value around 35 s, and then remain at around 0.18 MPa until 60 s or so.Afterwards, the total pressure curves vary with the change of F10, F9, and F3 opening curves and the pressure reduction of high-pressure gas cylinders.The maximum error between simulation and test is less than 12%.
As shown in Figure 9e, the variation trend of the flow-rate simulation curve at GP23 initial grid is the same as that of the flow-rate test curve of the exhaust bypass.Because of the problem of flow-rate measurement method and the modeling accuracy of F3, F9, and F10 valves at some openings under normal-temperature condition, the simulation curve is mostly higher than the test curve.Due to the limited test data and technical specifications, the modeling accuracy of F3, F9, and F10 valves cannot be further improved.However, the mean error between simulation and test is less than 19%.
As shown in Figure 9f, in the 0~1400 s dynamic change process, the variation trend of the total temperature simulation curves in the front section of the exhaust bypass is consistent with that of the test curve, except that the simulation curves have a peak around 16 s when F3 opens, and the maximum error in the rest of the time period is less than 2%.

Low-Temperature Operating Condition
As shown in Figure 10a, similar to the normal-temperature operating condition, the variation trend of the pressure simulation curve at GP1 middle grid is consistent with that of the normal-temperature air supply pressure test curve in the 1240 s rig-regulating process of the low-temperature operating condition, and the maximum error is less than 3%, except that the test curve with only 0.01 MPa measuring accuracy moves down in staircase form and by an amplitude a little bigger than the simulation curve.In addition, the test appears to have air replenishing operation around 1050 s.
As shown in Figure 10b, the variation trend of the flow-rate simulation curve at GP22 outlet grid is consistent with that of the flow test curve at the engine inlet, and the mean error between simulation and test is less than 10%.In addition, as in the normaltemperature operating condition, the low-temperature operating condition also has the problem of inaccurate test flow rate.In the low-temperature condition, the low-temperature air flow rate is controlled by the valve in front of FS1 and not affected by F11, which remains fully open during 0~1184 s; after the main control valve F4 of the normal-temperature air supply line opens at 125 s, its opening varies around 0.2 and the variation range is less than 0.1, so the changes in the F4 opening and the low-temperature air flow rate only control the total flow rate after mixing of the low-temperature and normal-temperature air, while the variation trend of the flow rate at the engine inlet is similar to that in the normaltemperature operating condition and is mainly affected by the change in the openings of F9 and F10.In low-temperature condition, the modeling accuracy of F9 and F10 valves is higher than that in normal-temperature condition.The Vanessa 30,000 valve throttling model may be associated with fluid temperature.
As shown in Figure 10c, the overall variation trend of the total temperature simulation curve at GP22 outlet grid is consistent with that of the total temperature test curve at engine inlet in the 0~1240 s whole process of the low-temperature operating condition, but the flow-field total temperature simulation curve descends at a speed much higher than the total temperature test curve during 30~400 s after the engine main line valve F9 is open.To be specific, at about 30 s, the valve F9 opens (Figure 7b), the low-temperature air below 250 K (as shown in Figure 10f) begins to flow into the pipeline behind F9, and the engine inlet total temperature T tm begins to drop.At 55 s, when the relative opening of F9 reaches 0.53, the simulation curve has decreased to 256.48 K, close to the incoming flow total temperature 247.85K (Figure 10f), but the test curve slowly drops to only 284.65 K and then to 258.55 K till about 400 s.By analysis on order of magnitude, the total temperature simulation curve drops from the normal temperature to near the inflow temperature on the order of tens of seconds, while the test curve costs the order of several minutes.The total volume of all components from F9 to GP22 outlet approximates to 0.78 m 3 , and it needs to take about 3 s to complete the displacement according to about 0.5 kg/s average flow rate of engine main line during 30~400 s, as shown in Figure 10b.Even considering the heating effect of heat transfer from the pipeline wall to the internal flow field and the mixing effect of the cold incoming flow with the local retained gas, it should not reach the order of several minutes measured in experiment for T tm to drop from normal temperature to a level close to the incoming flow temperature.The main reason why the test curve moves down slowly is that the location of the temperature sensor in the measurement rake is near the pipe wall; the measured temperature is the total temperature of near-wall flow field, while the 1D flow-field simulation result is the average temperature of each grid.Under the condition that F9 relative opening is no more than 0.6 (0.02~0.6) during 30~400 s, the cold jet flow from the F9 valve spool (DN250) rapidly fills and cools the center of the GP21 (inner diameter 261 mm)-GSA2 (inner diameter 834 mm)-GP22 (inner diameter 300 mm) pipeline, but the filling and cooling to the near-wall region is relatively slow, so the temperature stratification occurs along the radial direction of cross section at T tm measuring point.Under the condition that the difference between the temperature of the normal-temperature flow field of engine main line and the temperature of the incoming cold flow is large after F9 opens, the test curve in the near-wall region, which does not reflect the average cross-section temperature, decreases comparatively slowly.This phenomenon also occurs in the normal-temperature operating condition: with the gradual opening of F9 during 60~100 s shown in Figure 7a, the normal-temperature air flow path divides into two paths of air supply line-engine main line and exhaust bypass from one path of air supply line-exhaust bypass, so the temperature of normal-temperature air decreases due to the expansion and cooling effect.Figure 9c shows that the total temperature test curve and the simulation curves of total temperature and temperature all begin to decline after F9 opens, but the falling speed of the test curve is slower than that of the simulation curve; the total temperature curves in the exhaust-bypass front section shown in Figure 9f also present a similar phenomenon.
As shown in Figure 10d, the total pressure simulation curves of GP23 initial and middle grids and the total pressure test curve in the exhaust-bypass front section are close to overlap in the entire time period.
As shown in Figure 10e, the variation trend of the flow-rate simulation curve at GP23 initial grid is consistent with that of the flow-rate test curve of the exhaust bypass, and the mean error between simulation and test is less than 16%.
As shown in Figure 10f, the variation trend of the total temperature simulation curves at GP23 initial and middle grids are consistent with that of the total temperature test curve in the exhaust-bypass front section, and the maximum error is less than 2%.As the main control valve F4 of the normal-temperature air supply line is not open in the first 125 s (Figure 7b), the total temperature in the exhaust-bypass front section keeps decreasing to 243.95 K at 125 s with the temperature decreasing of the cold source air in front of FS1.As F4 opens at 125 s, the normal-temperature air flows into GMA1 and mixes with the low-temperature air, GP9 air total flow rate increases and the temperature also goes up, so the total temperature in the exhaust-bypass front section starts to rise.Afterwards, the total temperature in the exhaust-bypass front section dynamically changes under the combined influence of the F4 dynamic regulation and the cold flow fluid source FS1 till 1184 s when the valve F11 in the low-temperature air supply line is closed.After 1184 s, it slowly returns to the temperature of normal-temperature air flow.

Conclusions
The modeling and simulation research on the core engine test rig main test system can be concluded as follows: (1) The high-precision complex system models reflect the real situation of the physical entities or systems, and its modeling technology is the key technology to realize the digital twins.The modularization modeling approach used in the core engine test rig main test system provides a solution to digital twin modeling.In the dynamic simulation of the whole process in two operating conditions of normal temperature 1400 s and low temperature 1240 s by using the 47-components numerical system, the variation trends of all simulation curves are consistent with those of the test curves.The maximum errors of gas source pressure, exhaust-bypass total temperature, and total pressure under the normaltemperature condition are lower than 4%, 2%, and 12%, respectively; the maximum errors of gas source pressure and exhaust-bypass total temperature under the low-temperature condition are lower than 3% and 2%, respectively; meanwhile, the simulation analysis fully and clearly reveals the dynamic variation process of the system with combining regulation of various valves.These indicate that the simulation models of various components and the established numerical system can be used for virtual test of the test rig in two operating conditions.However, it still needs to further improve accuracy, especially with respect to different operating conditions.
(2) A comprehensive solution is provided to the valve-spool/orifice throttling modeling at the system level, which is a key issue of the fluid pipeline system modeling.Vanessa 30,000 valve-spool throttling model is used for the modeling of F9, F10, and F11, the valves with known information in the test rig system; a modeling approach which utilizes system test data to correct the reference-valve throttling model is proposed for those valves with unknown structure and throttling characteristic information.This approach consists of three options, where for F3, the regulator valve with upstream and downstream system test data, HCB high-capacity sleeve regulator valve is selected as reference valve and the calibration-adjusting scheme (Option A) is used for modeling; for F4, the regulator valve only with indirect system test data, HCB valve is used as reference valve and the full-system simulation trial-calibration scheme (Option B) is used for modeling.The variation trends of the flow-rate simulation curves of engine main line and exhaust bypass in two operating conditions are consistent with those of the test curves, and the maximum error between simulation curves and test curves of the remaining parameters in the two branches is less than 12%, which indicates that the proposed method can solve the modeling problems of various regulator valves in the complex pipeline network system.
(3) In practice, both Vanessa 30,000 and HCB are used as reference-valve models separately for the modeling of F3 and F4 so as to assess their validity.Taking F3 as an example, it is found by assessment that the flow-rate simulation results of the two models are the same at the zero opening point, and the difference between the two models with the opening in the range of 0.42~0.63,including the calibration opening point of 0.44, is less than 20%; outside such range, the difference between the two models with the opening in 0.33~0.42and 0.63~0.94 is less than 20~50%.These indicate that the proposed modeling approach for information-unknown regulator valves has no strict restrictions on the reference valve and has certain versatility.However, it needs to be reassessed based on the corresponding test when used for the opening points outside the permitted error range.At this time, with the accumulation of multi-operating conditions test data, the modeling can use Option C, under which the flow characteristic is obtained by fitting based on the reference-valve throttling flow-rate equation and the test data.
(4) The flow-rate data obtained from the tests of two operating conditions are obtained, instead of from direct measurement, indirectly by calculation of the total static pressure method which uses the pressure and total pressure data measured by the sensor, but the test flow rate is very important to the modeling of the information-unknown valves; if the total temperature measurement point of engine main line can be close to the central flow area of the pipe, it can provide a temperature close to the mean flow for the comparison of the one-dimensional flow field simulation.Therefore, there is room for improvement of flow-rate and temperature measurements of the test.

Figure 1 .
Figure 1.Schematic diagram of the core engine test rig system.

Figure 1 .
Figure 1.Schematic diagram of the core engine test rig system.

Figure 2 .
Figure 2. Finite control volume grids of gas valve.

Figure 3 .
Figure 3.Comparison of mass-flow-rate calibration and adjusting simulation results of F3 selecting two reference-valve models with experimental measurements: (a) Vanessa 30,000 butterfly valve throttling model; (b) HCB cage-type regulator valve throttling model.

Figure 3 .
Figure 3.Comparison of mass-flow-rate calibration and adjusting simulation results of F3 selecting two reference-valve models with experimental measurements: (a) Vanessa 30,000 butterfly valve throttling model; (b) HCB cage-type regulator valve throttling model.

Figure 4 .
Figure 4. Finite control volume grids of gas steady apparatus.(a) Three-chambers and two-orifices structure.(b) Twochambers and one-orifice structure.

Figure 4 .
Figure 4. Finite control volume grids of gas steady apparatus.(a) Three-chambers and two-orifices structure.(b) Twochambers and one-orifice structure.

26 Figure 5 .
Figure 5. Finite control volume grids of gas mixing apparatus.

Figure 5 .
Figure 5. Finite control volume grids of gas mixing apparatus.

Figure 6 .
Figure 6.Numerical simulation model of the core engine test rig main test system.

Figure 7 .Figure 7 .
Figure 7. Regulating time sequence of two operating conditions of the core engine test rig main test system.(a) Normaltemperature condition; (b) low-temperature condition.

Figure 7 .Figure 8 .
Figure 7. Regulating time sequence of two operating conditions of the core engine test rig main test system.(a) Normaltemperature condition; (b) low-temperature condition.

Figure 8 .
Figure 8. Boundary conditions of two operating conditions of the core engine test rig main test system.(a) Normaltemperature condition; (b) low-temperature condition.

Figure 9 .
Figure 9.Comparison of simulation results with experimental measurements on the 1400 s normal-temperature condition.(a) Pressure of normal-temperature air at inflow measurement section.(b) Mass flow rate at engine inlet measurement section.(c) Total temperature at engine inlet measurement section.(d) Total pressure at bypass measurement section.(e) Mass flow rate at bypass measurement section.(f) Total temperature at bypass measurement section.

Figure 9 .
Figure 9.Comparison of simulation results with experimental measurements on the 1400 s normal-temperature condition.(a) Pressure of normal-temperature air at inflow measurement section.(b) Mass flow rate at engine inlet measurement section.(c) Total temperature at engine inlet measurement section.(d) Total pressure at bypass measurement section.(e) Mass flow rate at bypass measurement section.(f) Total temperature at bypass measurement section.

Figure 10 .
Figure 10.Comparison of simulation results with experimental measurements on the 1240 s low-temperature condition.(a) Pressure of normal-temperature air at inflow measurement section.(b) Mass flow rate at engine inlet measurement section.(c) Total temperature at engine inlet measurement section.(d) Total pressure at bypass measurement section.(e) Mass flow rate at bypass measurement section.(f) Total temperature at bypass measurement section.

Figure 10 .
Figure 10.Comparison of simulation results with experimental measurements on the 1240 s low-temperature condition.(a) Pressure of normal-temperature air at inflow measurement section.(b) Mass flow rate at engine inlet measurement section.(c) Total temperature at engine inlet measurement section.(d) Total pressure at bypass measurement section.(e) Mass flow rate at bypass measurement section.(f) Total temperature at bypass measurement section.

Table 1 .
Modularization disassembly of the core engine test rig main test system.

Table 2 .
Types and throttling calculation schemes of ten valves in the 47-components system.
[25]: NTC is the normal-temperature condition; LTC is the low-temperature condition.n_Type_Cd: 3 is the flow character fitting formula scheme of HCB cage-type double-seat high-capacity sleeve regulator valve[25]; 8 is the flow character fitting formula scheme of Vanessa Series 30,000 valve.ppl.Sci.2021, 11, x FOR PEER REVIEW 17 of 26