CFD-aided approach of modelling and dynamic characteristic optimization for a highly nonlinear auxiliary braking system

It is a considerable challenge to develop a mathematical method for modelling the coupling relationship between the control signal and braking performance in the design and optimization, especially with a highly nonlinear system. The present study thus establishes an integrated CFD-aided method, elaborately proposed based on a 1D/3D co-simulation to accomplish a goal of collaborative calculation between multi-physics fields, and conducting the braking and thermal properties that changes dynamically with the control signal. In particular, the 3D model is actually a multi-domain coupling CFD transient calculation, and the 1D model provides the dynamic boundary conditions. These results manifest some interesting issues, such as modelling the prediction interaction between the isolated control signal and dynamic braking capability. The synergies of the blade agitation and centrifugal force in the multi-phases development affect the oil-filling control and braking generation. The control pressure (Pair ) and throttle area (Aout ) are predominated in enhancing braking torque, shortening response time, and improving heat dissipation efficiency. The operations of Pair  ≤ 3.2 bar and Aout  ≥ 64π mm2 ensure a safe and reliable design and then optimize braking and control characteristics. The resulting maximum braking torque is 4070Nm, and the outlet oil temperature is reduced to 138°C.


Introduction
Heavy-duty vehicles are an important transport means and are widely used to deliver people and goods over long distances, such as in public transportation (Bravo et al., 2018;Broniszewski & Piechna, 2019;Kumar & Rajagopal, 2020). In China, 31% of the length of highways passes through and over mountains, hills, and plateaus, and therefore, have long downhill sections (Qian et al., 2020). Most of the commercial trucks are heavy-duty vehicles and need high braking power while slowing over a long distance. Traditional frictional braking transforms the kinetic energy into heat energy, which can cause brake fade; this not only is a safety risk but also prevents energy recovery and reuse (Xu & Cai, 2020;Ye et al., 2019;Yevtushenko et al., 2015). The braking system of a heavy-duty vehicle requires a high braking torque and quick response time, and it must be able to handle the frequent starts and stops as the vehicle travels at full speed and along downhill sections. It is difficult for traditionally available braking systems to simultaneously fulfil all of these requirements.
CONTACT Chunbao Liu liuchunbao@jlu.edu.cn Supplemental data for this article can be accessed here. https://doi.org/10. 1080/19942060.2022.2103589 Auxiliary braking systems have been designed and installed in heavy-duty vehicles to ensure their driving safety and economy, as shown in Figure 1(a) (Anwar & Stevenson, 2011;Bae et al., 2015;Kumar et al.,2021a;Tian et al., 2020;Zhang et al., 2014). In particular, hydrodynamic auxiliary braking systems are currently applied in heavy-duty vehicles because of their advantages regarding braking performance and control characteristics (Liu et al., 2017a). Recently, far-ranging studies have focused on conventional hydrodynamic auxiliary braking systems in recent years. Hur et al. used the steady CFD model to numerically study the multi-phase flow of the working chamber to clarify the braking torque generation in different oil-filling stages (Hur et al., 2018). Yan et al. combined CFD with a proper orthogonal decomposition model to reconstruct and predict the evaluation of a 3D flow field. The same group also studied bionic surface structures with the aim of reducing the drag torque under non-braking conditions Mu et al., 2016). Liu et al. developed a new CFD method that coupled the scale-resolving simulation method with a dynamic mod-ified viscosity coefficient, optimizing the cooling system of the heat exchanger and the blade structure of the working chamber to improve the maximum braking efficiency and achieved a maximum braking torque of 2800 Nm Liu et al., 2017b). Nonetheless, the response time, heat exchange capability, and maximum braking torque were unsolved, restricting the further development and application of auxiliary braking systems in heavy-duty vehicles (Zheng et al., 2017).
To solve the above problems, a novel auxiliary braking system integrated with an advanced mechanical-electroliquid system and control strategy was designed and optimized, as shown in Figure 1(b), to adjust the speed and ensure safety . To mechanistically reveal how the machine works, the schematic diagram of the auxiliary braking system is shown in Figure 1(c). The 3D zone is a map of the working chamber, which elucidates the hydrodynamic transfer and the braking torque generation. The 1D zone is map of a mechanical-electroliquid control system, including the information control flow consisting of electric signal, control airway, oilcharging way, and oil-discharging way. However, without a robust mathematical prediction model between the oilfilling process of the working chamber and the control signal of the hybrid control system (Mu et al., 2019;Zheng et al., 2018), the nonlinear relationship between the control information and the details of flow field cannot be quantified effectively (as shown in Figure 1(c)).

Figure 1.
Overview of the mechanical-electro-liquid hydrodynamic auxiliary braking system, showing (a) application and configuration layout of the auxiliary braking system on heavy-duty vehicles; (b) working principle and main components of the auxiliary braking system; and (c) the key problems in the design including the interaction between the control system and braking systems.
These issues have previously been overlooked, but they are regarded as helping in solving the coupling connection between the isolated control signal and dynamic braking capability, thus improving the auxiliary system design's effectiveness and reliability.
Extensive investigations have illustrated the mathematical calculations between the control information of the hydraulic system and the braking characteristics of the working chamber using a 1D simplified numerical model of hydrodynamic retarder (Mu et al., 2019). This method has developed into the integrated model of an auxiliary braking system, which is favourable for modeling the interaction of both systems. However, the braking torque generation varies according to the unsteady turbulent flow of the working chamber rather than being a linear change (Kong et al., 2020). In this circumstance, considering the mathematical prediction model between the multi-physics coupling field of the working chamber and the control signal of a hybrid control system, an integrated model based on the 1D/3D co-simulation may currently be a proper choice (Singh & Abbassi, 2018;Zhang et al., 2021). A 1D/3D co-simulation method that combines a 1D control system model and a complicated 3D CFD model have been used in various fields to resolve the multi-physics coupling relationships between the control information, thermodynamics, and nonlinear flow dissipation of complex mechanical-electro-liquid systems (Lu et al., 2016). Numerical approaches with the dynamic boundary conditions that can accurately predict the synergy of the complex flow and mechanical-electroliquid control signal deeming to be rather serviceable. Such an approach may quantitatively and qualitatively clarify the dynamic coupling relationship betwwen the control information and brake performance, especially in a nonlinear auxiliary braking system design.
In this study, we present an integrated CFD-aided method. Its novelty could be described as a 1D/3D cosimulation elaborately proposed to accomplish the goal of collaborative calculation between multi-physics fields, and conducting the braking and thermal properties that change dynamically with the control signal. It solved the previous problem of developing a robust mathematical method to master the coupling relationship between the control signal and braking properties in a highly nonlinear auxiliary braking system. Simulations and experiments are performed for verification and manifest how to resolve the connection between the isolated control signal and dynamic braking torque. Moreover, the indepth philosophy of the multi-physics coupling field of the unsteady rotor/Stator domain during braking is analyzed to optimize the response of the design parameter to external braking characteristics. This will guide the further design and development of braking strategies and control methods to ensure the overall braking efficiency and reliability of auxiliary braking system.

3d CFD numerical model
(1) Governing equations. The volume of fluid (VOF) model is mainly used to track the multi-phase flow in the working chamber (Wang et al., 2019;Yang et al., 2022). Each calculation unit contains two phases, gas (air) and liquid (engineer-oil), where the sum of air volume fraction (ϑ 1 ) and oil volume fraction (ϑ 2 ) is 1. The governing equations of the flow field are the conservations of mass, momentum, and energy (Peng et al., 2018). To simplify the calculation, assume that the hydraulic oil medium is considered a continuous, incompressible, and isotropic Newtonian fluid. The governing equations of the CFD method can be expressed as follows: Mass conservation equation: where ∇, ∇ = ∂ ∂x , ∂ ∂y , ∂ ∂z is Hamiltonian, and α k is the volume fraction of the k th phase.
Momentum conservation equation: Where ρ m , μ m and v m is the average density (ρ m = n k=1 α k ρ k ), average dynamic viscosity (μ m = n k=1 α k μ k ), and average mass velocity (v m = n k=1 α k ρ k v k ρ m ), respectively. ρ m g is the volumetric gravitational force and F is the volumetric force vector. α k is the volume fraction of the k th phase.
Energy conservation equation: Heat transfer model: (2) Multi-domain interface coupling model. The sliding mesh method can be set between different computational subdomains. Adjacent subdomains then slide along the interface according to the defined motion (Jung et al., 2011). A new mesh interface between subdomains needs to be determined for each time step to realizing real-time coupling between different domains through real-time flow transfer. When the cell volume is zero, this results in: Then, the following is obtained: In addition, the following simplification is obtained: where n f is the number of faces on the control volume, u g is the velocity of moving mesh, and A j is the face area vector. A coupling wall is an alternative approach for modelling non-moving walls between multi-domains. This method better handles the interface between domains that intersect or have gaps. Nevertheless, this can increase the stability and reliability of simulations between subdomains of the hydrodynamic auxiliary braking system because the interaction of interface data can cause the simulation to fail.

1d control system mathematical model
(1) electromagnetic control model. The electromagnetic force can be expressed by the electromagnetic energy balance equation (Chen et al., 2020): (2) proportional valve model. The control current (I cr ) and the spool displacement (x) are important parameters that affect the outlet pressure. The force balance equation of the spool is.
For the convenience of analysis and calculation, the flow is usually regarded as the isentropic flow of ideal gas. The mass flow rate (Q) is given by The maximum mass flow rate (Q max ) is: The flow rate of the throttle increases with the increase of the pressure difference between the inlet and outlet, where P A P u = 0.528, and Shrinkage phenomenon and friction loss occur in the valve, which needs to be corrected with shrinkage coefficient correction (C d ): (3) gas-liquid conversion model (oil tank). The oil tank is a gas-liquid conversion module, which determines the dynamic parameters between the working chamber and the hybrid control system, including gas-liquid energy conversion. The flow rate of engine oil in the tank is calculated as follows: (4) heat dissipation model. The heat dissipation model has two parts: passive and active (Kumar et al., 2021b). The passive dissipation part is the engine oil's heat conduction (W s ) through the working chamber shell and the heat convection (W a ) between the shell and the air. The active dissipation part is the heat exchanger (W w ), which exchanges heat between oil and water. The dissipation model is given by (Lu et al., 2016): The heat conduction by the shell of the working chamber is the energy transfer caused by the temperature difference: When the cold air flows through the hot shell, part of the heat is taken away by convective heat transfer: A heat exchanger circulates coolant to take away heat. The energy dissipated by the engine oil is equal to the heat absorbed by the coolant: Where, η 0oil and η 0water each corresponds to the total heat transfer area, so the fin heat transfer is equal to the plate heat transfer in the heat exchanger.

CFD-aided approach framework
A 1D/3D co-simulation that relies on the critical components given in Figure 1(c) is established for the integrated CFD-aided method. The CFD-aided approach framework is illustrated to elucidate the nonlinear coupling relationship between the working chamber and the hybrid control system (as shown in Figure 2). The framework considers the complex physical quantities that not only affect individual components but also interact with each other. The 3D CFD model is mapped to elaborate the thermodynamic characteristics and details of the flow field. The 1D control sub-system model is mapped from a combined electromagnetic control model, proportional valve model, gas-liquid conversion model, heat exchanger, and other components. It is established to provide the dynamic boundary conditions for 3D CFD model. Moreover, the co-simulation control block is proposed to update the exchange parameters between the working chamber and the hybrid control system in simulation; namely the inlet flow rate (Q in ), outlet flow rate (Q out ), inlet pressure (P in ), outlet pressure (P out ), inlet temperature (C in ), and outlet temperature (C out ). We envision that this approach will elaborately solve the interaction between the hybrid control system and the working chamber, conducting the braking and thermal properties that change dynamically with actual conditions in the simulation.
The key issue for the collaborative calculation between 3D CFD models and the 1D control system is the data exchange between the two subdomains . A collaborative calculation control block is established to transform and update boundary conditions of the 1D/3D co-simulation. The inlet flow rate (Q i in ), inlet pressure (P i in ), and oil temperature (T i in ) are obtained after a certain number of iterations, as shown in Figure 2. The 1D control sub-system reads the inlet pressure P i in and Figure 2. Framework of the integrated CFD-aided solution method for the coupling relationship between the working chamber and the hybrid control system, including the 1D model for the dynamic boundary conditions, 3D model for the CFD transient calculation, and the co-simulation control block.
obtains the inlet flow rate Q i in by using the revised Eq.s 10 in the mechanical-electro-liquid control system. It then provides the flow rate Q i+1 in to update the boundary conditions of the 3D CFD model for the next coupled step. Q i+1 in is updated according to.
where a is a weighting factor ranging from 0 to 1 that determines the convergence rate and stability and i is the coupled time steps. Owing to the limited paper space, the updates of other boundary conditions are not given in this paper. Particularly, the setting principles for the boundary conditions of the inlet and outlet of the working chamber are highly consistent. As shown in Figure 3, the co-simulation control block is established for multi-time steps based on synchronous serial coupling (SSC) mode . This is based on a linear interpolation scheme and is simpler than other methods currently available, and it recognizes the boundary data interchange and update. Also, a sub-cycling method is adopted for a smaller time step (t n+1 = t n + Δt) in the 3D CFD calculation of the working chamber. This allows a CFD process to run with smaller time steps in a coupling step. Then, the CFD data are transferred to the 1D control system. The mechanicalelectro-liquid control system is analyzed and simulated in the 1D control sub-system with a large time step (T m+1 = T m + ΔT). The boundary conditions of the 3D CFD sub-system are then updated.

3d CFD transient calculation
A CFD model of the working chamber is modelled to investigate the flow conditions and multi-physics coupling field . The geometric parameters are extracted from the design model, listed in Table 1. Figure 4(a) shows the computational domain and the boundary conditions, including a stator, rotor, inlet, and outlet. Moreover, the inlet and outlet boundaries receive the dynamic boundary information from the 1D control sub-system in real-time. As shown in Figure 4(b), the hexahedral cells are used for the whole passage. To better satisfy the requirements of the enhanced wall function method or the mesh in near-wall regions, a grid refinement method is used to guarantee an average dimensionless wall distance y + < 2 around the blades . The flow field changes periodically in the working chamber. The sliding mesh method is used in the unsteady rotor/Stator domain to deal with the multidomain coupling interface (as shown in Figure 4(c)). The coupling wall is used otherwise.
The Courant-Friedrichs-Lewy (CFL) number can be used to evaluate the convergence for the unsteady Figure 3. Synchronous serial coupling mode for the data exchange and coupling calculation of the co-simulation control block, resolving between the 1D control sub-system and the 3D CFD sub-system. CFD calculation, defined as CFL = U ∞ t/ x = 1, where x is the typical grid size, and the time step, t = 5×10 −4 s, is set in the transient calculation. It ensures sufficient calculation time for the gird sliding along the interface during the calculation. Stressblended eddy simulation (SBES) was chosen to capture the thermal-flow details of the unsteady rotor/Stator domain in CFD. Moreover, SBES is employed to achieve the blending on the stress level between the Reynoldsaveraged Navier-Stokes (RANS) and large eddy simulation (LES) formulations. The flow field structure can be solved by switching more quickly from RANS to LES in the separated shear layer. The solver can be set to the pressure-velocity coupling model in Fluent, and the SIMPLEC algorithm is applied. Table 2 lists the other setting conditions. The outlet pressure and braking torque are monitored at the start of the calculation. The convergence criterion is set for each inner iteration so that the mass, momentum, and energy residuals are below 10E-5. The grid sensitivity is analyzed to reduce the influence of the grid number on the calculation accuracy, as shown in Figure 4(d). The change rate (R(n)) is the numerical result of the braking torque calculated with the adjacent grid density, defined as: Where T B (n) is the braking torque result under the current calculation grid, and T B (n+1) is the braking torque result for the next grid density. The grid density is considered satisfactory in engineering applications (i.e. the calculated results are the independent of the grid) when R(n) is less than 3% and the increase in grid density is neglected . Based on the analysis results, as the number of grids increases, the calculation error decreases while the calculation time increases. The grid number is set to 7.2 million by the precision and computation time trade-offs. Herein, R(n) is less than 3% under the adjacent grid density, and the calculation time required per 100-time steps is 2.6 h. The maximum cell size is 3 mm, the first layer for the nearwall region of the O-block has a height of 0.024 mm, and the near-wall region has a growth rate of 1.2. A convective Courant number (CCN) of 0.8 is obtained for time-explicit simulations of the unsteady Rotor/Stator domains. This ensures that the distance travelled in a single time step is less than through the neighboring grid. Figure 3(d) also shows the changes in grid scale with the dimensionless wall distance, where the dimensionless distances 0 and 1 correspond to adjacent blades. For the near-wall region, the minimum grid scale is 1.634 × 10 −4 mm 3 , and the grid volume has a growth rate of less than 1.5.

1d model for dynamic boundary conditions
Mechanical-electro-liquid control system manipulates the braking and thermal properties in the retarder; herein, the 1D control system model results are regarded as providing the boundary conditions for the 3D CFD transient calculation in the 1D/3D co-simulation. Based on the mathematical model described above, a 1D numerical model of the control system is developed using the Pneumatic Component Design (PCD) and Mechanical libraries of AMESim and the bond-graph method (Chen et al., 2020). As shown in Figure 5(a), the 1D control system numerical model has simulated the energy and control information during the braking process through the electronic control model, oil tank model, heat dissipation model, and throttle valve model given in Section 2.1.1. To ensure the effectiveness and robustness of the 1D numerical model, the physical environment is set according to actual working conditions. The parameters are set as listed in Table 3.
For the numerical analysis, the average current of the pulse width modulation (PWM) is taken as the input signal to control the air pressure of the pneumatic proportional valve. Also, 1D numerical model results are continuously updated as the boundary condition of the inlet and outlet for 3D CFD model. The gas-liquid conversion model of the oil tank determines the pressure and flow at the inlet boundary. Furthermore, the pressure and flow data of the throttle valve are calculated in the simulation and are regarded as a boundary condition for the working chamber outlet. The hydrodynamic auxiliary braking system's heat dissipation is obtained using the active and passive dissipation models in the 1D sub-system (Dong et al., 2020).

Boundary data interchange and update
The 1D/3D co-simulation is encapsulated to solve the multi-physics coupling interaction during braking, which will quantitatively and qualitatively clarify the synergistic effect between the working chamber and hybrid control system. The main procedure of 1D/3D co-simulation method can be summarized as follows: (1) Solve the 1D control system subdomain with the bond graph method to obtain 1D data at the time ΔT; (2) Set the boundary conditions of the 3D CFD subdomain based on the 1D data, and solve the dynamic 3D flow field at the time Δt; (3) Calculating the 1D zone based on the interchange parameters from 3D zone and provide the updated boundary conditions for 3D CFD.
Also, a co-simulation control block is built to realize a dynamic information interaction between the 3D CFD model and 1D control system model, following in Figure 3. The block also realizes the whole-system model's real-time monitoring and batch computing. The control block reads and writes a data table to complete the exchange of interaction parameters.
The 1D/3D co-simulation model is established by mapping and integrating the relationship between the physical and mathematical models in Figure 5. The cosimulation model includes an engine-oil tank, mechanical-electro-liquid control system, heat dissipation, and open working chamber, as shown in Figure 5(a) and (b). At the coupling interface of the co-simulation control block shown in Figure 5(c), the 1D data from the boundary of the inlet of the working chamber are transformed into 3D data for the inlet boundary of the 3D CFD model. Meanwhile, the 3D data of the 3D CFD model are converted into 1D data for the boundary of the outlet of the 1D control system. These conversions facilitate the use of different data and combine information for control and braking of the auxiliary braking system. This method can be used to realize the interaction of flow and control information and the dynamic analysis of the thermodynamic, braking, and control systems.

Verification and comparison
A test rig of the hydrodynamic auxiliary braking system is built to compare and verify the braking performance and control characteristics. Part 1 in Figure 6 shows the hybrid mechanical-electro-liquid control system of the test rig. Here the signals of the combined pneumatic and hydraulic control system are acquired. These include C in , Q in , P in , C out , Q out , and P out . Part 2 is the main test part and comprises the speed/torque sensor, drive motor, and working chamber. Part 3 is the data collection system, which obtains the dynamic characteristics. The test results are analyzed to obtain the characteristics of the sub-system and the interaction between the control information and braking generation.

Verification of sub-system properties
The 1D control subdomain's simulation model focuses on the proportional valve's response to specific input signals. Thus, the controller is not introduced to form a closed control loop. When the air source pressure is 8 bar, the hybrid control system sets the input current signal (Part 1). The control pressure (absolute pressure) is output by the simulation model of the solenoid proportional valve. The current of the proportional valve is adjusted by a handle, which controls the outlet pressure of the proportional valve. Figure 7(a) shows that I cr is linearly related to the output pressure; that is, the control current determines the filling rate of the working chamber. Figure 6. Schematic diagram of characteristics of the test rig for the auxiliary brake system, showing the hybrid mechanical-electroliquid control system, braking generation system, and data collection system.
Meanwhile, the output control pressure (P air ) is measured with the test rig. The maximum error is 5.1% at 210 mA. Thus, the simulation results for the 1D control system model meet the requirements for application.
The braking torque is obtained using the 3D CFD model for a constant filling rate at q 0 = 0.6. This is then compared to the test rig results for speeds of 400-800 rpm. Figure 7(b) shows that the CFD simulation can accurately predict the braking torque with an absolute estimation error of less than 3.7%. The results of the CFD simulation agree with the theoretical formula, which indicates that T B tends to change with n 2 B . The torque is expressed as.
The pump torque coefficient (λ B ), calculated from the pump torque, quantifies the braking performance. The variations in the numerical results correlate well with the theoretical calculations at different speeds. The verification with the test rig shows that the CFD method with the SBES model can reproduce the braking characteristics. The section of the chord plane and the nodal plane is taken as the reference surface to obtain the flow mechanism in Figure 7(c). Obviously, find that the flow field of the working chamber is the complex movement and evolution of turbulence of the working medium. Turbulent kinetic energy (TKE) and the rate of dissipation of rotor/Stator domains are qualitatively analyzed to clarify the flow mechanism and demonstrate the reliability of CFD. To validate the accuracy of the CFD-aided method of braking characteristics prediction, we selected the high-power hydrodynamic retarders mentioned in the past literature . The CFD physical model and the mesh model were conducted using the methodology that we mentioned in Li's study (refer to Figure S1 in Supporting Information). The resulting braking torques agree well with the data in reference, and the maximum absolute error is less than 3.6%. Within the research scope of the current study, the CFD-aided method can accurately predict the braking characteristics and captures the details of the flow field.

Comparison of external characteristics
A comparative analysis is performed to verify the dynamic braking process for I cr = 560 mA and from n = 400 rpm to n = 1300 rpm. Figure 8 shows the variations in the braking performance and control characteristics under different conditions, where the subscript T indicates the test value and the subscript 1D/3D indicates the simulation value. Figure 8(a) shows that the 1D/3D cosimulation model is more accurate than the calculation results using the traditional 1D model (Mu et al., 2019) when predicting the braking torque at higher rotation speeds (900-1300 rpm). This is because the traditional model does not consider the dynamic data interaction between the control information and braking generation. The braking torque obtained by the 1D/3D co-simulation model agrees well with the experimental data with a maximum error of less than 5%.
The dynamic results of the numerical simulation for the control current (I cr ), control pressure (P air ), output pressure (P out ), braking torque (T B ), and oil temperature (C) are compared with the results obtained using the test rig when n = 900 rpm, as shown in Figure 8(b), (c), and (d). The simulated inlet pressure, outlet pressure, and braking torque are consistent with the experiments. Furthermore, the results demonstrate that the response time of the control air pressure experiences 0.31 s (i.e. the control air pressure signal is obtained at 2.31 s), and the braking torque reaches a maximum value of 3420 Nm at 4.6s. Figure 8(d) demonstrates the braking generation and oil temperature, and the simulation results closely fit the experimental results. Thus, the braking torque increases the energy dissipation of the medium in the braking process, which is then rapidly converted into internal energy and causes the oil temperature to increase. The computational accuracy is considered acceptable. The proposed method combines the information of the control and braking systems with a multi-physics coupling field for systematic analysis of the coupling among the control systems, energy consumption, and braking generation.

Analysis and optimization during dynamic braking
The interaction mechanism of the multi-physics field is investigated for the dynamic oil-filling in the unsteady Rotor/Stator multi-domain. It reveals an important issue for the coupling relationship between the oil-filling process of the working chamber and the control signal of the hybrid control system. The optimize designs' responses to external properties are investigated to improve braking characteristics and ensure a safe and reliable operation.

Coupling relationship between the control signal and flow conditions
The flow flux and air pressure difference at the inlet and outlet of the working chamber are directly related to the maximum braking torque and the response time. The hydrostatic transmission control is affected by the flow flied of the inlet and outlet of the working chamber. The inlet and outlet flow rates are dictated by a number of elements, such as control air pressure, oilfilling rate, outlet throttling area, etc. (Mu et al., 2016). The inlet air pressure and outlet throttling area are constants in this process, and the oil-filling rate in the working chamber depends directly on the inlet and outlet flow rates. The flow flux in the open working chamber is. Figure 9 shows that Q in and Q out change dynamically with time, and Q in first becomes nonzero when the pressure signals (I cr = 560 mA, n = 1000 rpm) are received. However, when the oil-filling rate in the working chamber quickly increases, Q in reaches a maximum value at 3.02 s, and Q in is 479 L/min. The inlet flows flux decreases and finally reaches a stable region at 4.75 s when Q in is 415 L/min. Owing to the small amount of oil in the working chamber in the initial oil-filling stage, there is a time delay before Q out first becomes nonzero at 2.73 s. Additionally, the actual flow flux (Q act = Q in-Q out ) is applied to analyze dynamic performance during the oilfilling process. The hydrostatic pressure of Q A increases at 2.26 s, which decreases the thrust and resistance at the inlet and outlet of the working chamber when ΔP A = 2.14 bar in the pressure field. In this stage, the working chamber is filled with hydrodynamic oil, which is thrown by centrifugal force to the outer ring of the working chamber away from the inlet. This results in decreasing thrust and resistance at the inlet and outlet of the working chamber. Q B reaches its maximum value Figure 9. Coupling relationship between the air control characteristics and flow conditions for the inlet and outlet flow domains when P air = 2.82 bar and d out = 10 mm (Q A is the initial oil-filling stage, Q B is the maximum flow flux difference, and Q C reaches the oil-filling balance stage). at 2.88 s when ΔP B = −3.15 bar, and the pressure difference increases significantly at this time. Meanwhile, Q act then drops rapidly with time. The pressure difference at the inlet and outlet reaches its maximum value when Q act = 0 L/min at 5.1 s, and oil-filling equilibrium is achieved in the working chamber. In contrast, the traditional calculation method sets the inlet and outlet pressures and flow flux as constant values, which incorrectly describes the boundary conditions for the numerical model of complex flow conditions (Liu et al., 2017a;Mu et al., 2016). With the proposed method, by coupling the information interaction between the control characteristics and flow conditions of multi-domains in this study, the maximum oil-filling rate is accompanied by dynamic changes in the pressure difference between the inlet and the outlet.

Development of the unsteady multi-phases in working chamber
The oil-filling rate variation, the air-oil multi-phase flow distribution in the unsteady rotor/stator multi-domains, and the velocity streamline for the dynamic oil-filling process, as described in Figures 10 and 11. As the engine oil fills into the working chamber, the distribution of the air-oil field and unsteady turbulence flow in different flow domains are obviously changed by the combined action of viscous and centrifugal forces. According to the complexity flow condition in the oil-filling process of the working chamber (Mu et al., 2019), the process can be divided into three stages, as shown in Figure 10(a). Stage 1 (0-2.6 s) is the initial oil-filling stage and directly affects the braking response time. The oil is mainly distributed in the inlet domain, and a small amount of oil enters the working chamber, as shown in Figure 10(b). The oil-filling process in the inlet domain is steady and encounters little resistance (see Figure 11(a)). Stage 2 (2.6-3.6 s) is the oil filling of the working chamber stage when the oil is quickly filled into the working chamber. The engine oil fills the working chamber to a certain extent, while a small amount of oil is expelled to the outlet domain, as shown in Figure 11(b). In terms of the volume fraction, the oil is distributed mostly outside the working chamber through the action of the centrifugal force, as shown in Figure 10(c). Stage 3 (3.6-4.54 s) is the oilfilling balance working stage when the outlet domain is quickly filled. In this stage, the flow experiences resistance during the oil-filling process, and it is separated entirely by the action of the blade agitation and centrifugal force (as shown in Figures 10 and 11(c)). This effect determines the maximum filling rate at equilibrium. The interaction of the flow field for the multi-domain that developed during the oil-filling is the main factor that affects the braking generation and control characteristics.  The present study thus analyzes the dynamic evolution and development of gas-liquid multi-phases for unsteady multi-domains, as well as the control signal changes. It will direct the co-design of control systems and braking systems for hydrodynamic auxiliary braking systems.

Interaction of the dynamic oil-filling process
As mentioned in the above analysis, as the oil-filling rate increases, the complex flow field of the working chamber affects the interaction of control characteristics and braking performance. The dynamic changes in the inlet and outlet oil pressures during the oil-filling process are shown in Figures 12 and 13. The oil pressures at the inlet and outlet of the working chamber is discussed, and the pressure field on the nodal plane of the stator is obtained at given times. The oil pressure of the inlet remains constant at 9.1 bar at 2.8 s, and the oil pressure at the outlet reaches a stable value of 18.5 bar at 5.1 s. This shows that Q ct increases first and then decreases and reaches its maximum value at 2.88 s. This is basically consistent with the trend of Q ct shown in Figure 9. Moreover, the qualitative and quantitative results in Figures 12 and 13 manifest that the pressure on the nodal plane increases rapidly with the oil-filling rate in the working chamber, and low-pressure areas are observed near the inlet and outlet. The nodal surfaces with inlet and outlet pressure field at Stage 1, Stage 2, Stage 3, and the balance stage are illustrated in Figure 12. The low-pressure area in the inlet appears at t = 2.1 s and q = 0.05 because transmission oil moves away from the inlet under the centrifugal force. According to the quantification results in Figure 13, the oil-pressure difference is smaller in an initial oil-filling stage. Moreover, with the filling rate increasing from 2.4 s to 4.1 s, the oil-pressure in the working chamber increases remarkably, and a local high pressure is formed between the inlet and outlet. The inlet and outlet pressure became steady at t = 5.1 s and q = 0.79, and high pressure is observed at most distances. The original inlet control air pressure and outlet throttle area remain unchanged, and the oilfilling process is in dynamic balance. Thus, the oil-filling rate can be improved by increasing the control pressure (P air ) and finally increasing the inlet pressure in a stable state. Additionally, the filling rate of the working chamber is increased by reducing the area of the outlet throttle (A out ).

Optimized designs response to external characteristics
P air and A out can be adjusted to optimize the braking generation, heat exchange capability, and response time based on the above analysis. Figure 14 compares the Figure 13. Quantitative analysis of oil pressure for the same streamline at the same position between the inlet and outlet (the dimensionless distances 0 and 1 correspond respectively to the inlet and outlet of the rotor and stator domains).

Figure 14.
Factors for the braking performance: trends in the braking torque and temperature according to the various stages (e.g. control air pressure is 2.4bar, 2.8bar, and 3.2bar; and outlet throttle diameter is 8, 10, and 12 mm). braking torque and oil temperature for the different control air pressures and different outlet throttle areas at a speed of 900 rpm. An increase in the control air pressure helps improve the working chamber's oil-filling process, and the response time is shortened in Stage 1, and the braking torque increases substantially in Stage 3. This is made possible by a higher control air pressure to counteract the centrifugal force on the medium-oil, which can extend the equilibration time while continuously raising the oil-filling rate of the working chamber. Also, the increase in control air pressure heightens the flow rate at the inlet and outlet, making the heat transfer capacity of the working chamber significantly higher and lowering the oil temperature. The maximum braking torque at P air = 3.2 bar is 4070 Nm. Increases P air to improve the oil-filling rate and improve the braking torque, and helping the outlet temperature for braking equilibrium is reduced to 138°C. The oil temperature at P air = 3.2 bar is lower than that at P air = 2.4 bar. Increasing the control air pressure ensures both an improvement in braking torque and a reduction in oil temperature. Still, the equilibrium is delayed under the joint action of the control air pressure and the centrifugal force. Reducing the outlet throttle area will increase the working chamber's oil-filling rate, and the braking torque rises substantially in Stage 3. However, equilibrium is reached at 152°C when A out = 16π mm 2 (d out = 8 mm). Compared with the other cases, more substantial braking torque is quickly generated when A out = 16π mm 2 , but the maximum temperature increases slightly. When A out is decreased to reduce Q out , the maximum braking torque eventually increases, and the heat exchange capability is reduced.
To evaluate the effects of P air and A out on the braking performance and control characteristics more reliably, the multi-physics fields are compared qualitatively and quantitatively for different cases in Figure 15. When P air is increased, the area of low pressure on the blades' suction surface increases, as shown in Figures 15(a) and (b). This is consistent with Figure 12, which shows that a higher control air pressure increases the inlet and outlet flow rates of the stable area. The high-temperature area around the blade is higher with increased P air , which has an important effect on blade design. The high-temperature area around the stator and rotor blades expands and intensifies under the combined action of the decreased outlet flow rate and increased braking torque. However, neglecting the thermodynamics at high P air may result in local high temperatures. An analysis of the braking and thermodynamic characteristics shows that adjusting P air and A out can improve the braking characteristics. The safety of blades and other components is affected by an increased temperature in the working chamber. In other words, the control air pressure (P air ≤ 3.2 bar) and throttle area (d out ≥ 8 mm) can be set to ensure a safe and reliable design, which various standards have defined as an oil temperature of less than 155°C . Based on the above analysis, P air and A out can be adjusted to ensure the driving safety and economy of the hydrodynamic auxiliary braking system for heavy-duty vehicles.

Conclusion
We have established an integrated CFD-aided method, which we elaborately proposed based on a 1D/3D cosimulation to accomplish the goal of collaborative calculation between multi-physics fields, and conducting the braking and thermal properties that changes dynamically with the control signal (see Figure 5). The 1D/3D co-simulation model is more accurate than the calculation results achieved the traditional 1D model (Mu et al., 2019). Also, the numerical results agreed with the experimental results, with a maximum error of less than 5% (see Figure 8(a)). The interaction mechanism of the multi-physics field for the dynamic oil-filling in the unsteady rotor/stator multi-domain is investigated. The synergies of the blade agitation and the centrifugal force in the multi-phases development affected the oil-filling control and braking generation. The optimize designs' responses to external properties were investigated to improve braking characteristics and ensure a safe and reliable operation. Increasing the control pressure and reducing outlet throttle were predominated in increasing braking torque, shortening the response time, and improving the heat exchange capability. Thermodynamic and dynamic analyses showed that operations of P air ≤ 3.2 bar and A out ≥ 64π mm 2 can ensure a safe and reliable design to heighten braking and control characteristics. The maximum braking torque was 4070 Nm, and the outlet oil temperature is reduced to 138°C (see Figure 14).

Summary
In this study, we demonstrate an integrated CFDaided method to derive the nonlinear prediction model between the braking propert of the working chamber and the control characteristics of the mechanical-electroliquid control system, thus enhancing the effectiveness and dependability of the auxiliary braking systems in the design process. The significant findings are summarized as follows: (1) An integrated CFD-aided method is established by mapping the newly-developed auxiliary braking system. The method is elaborately proposed based on a 1D/3D co-simulation to accomplish the goal of collaborative calculation between multi-physics fields. The 3D model, interchanging the updated boundary data, is a multi-domain coupling CFD transient calculation, and the 1D model provides the dynamic boundary condition for it according to the driving state of the truck.
(2) Verifications of the braking generation and control characteristics show that the response time and maximum braking torque reached 0.31 s and 3420 Nm. The numerical results agreed with the experimental results, having a maximum error of less than 5%. We found that this method resolves the connection between the isolated control signal and dynamic braking capability. (3) The maximum oil-filling rate of the working chamber is accompanied by a dynamic change in the pressure difference between the inlet and the outlet in the control system. And the blade agitation and centrifugal force synergies in the multi-phases development dictate the oil-filling control and braking generation. (4) Increasing the control pressure and reducing outlet throttle are predominated in increasing braking torque, shortening the response time, and improving the heat exchange capability. Thermodynamic and dynamic analyses show that operations of P air ≤ 3.2 bar and A out ≥ 64π mm 2 can ensure a safe and reliable design to heighten braking and control characteristics. The maximum braking torque is 4070 Nm, and the outlet oil temperature is reduced to 138°C.
The proposed method provides a foundation for designing collaborative control braking systems. However, this method may not yet be applicable to certain fields, so its universality will be not guaranteed. Future investigations will focus on applying and applying an integrated CFD-aided method universally to the different mechanical-electro-liquid systems, shedding light on the optimal design of new products.