## 3 Simulation Model and Simulation Methods
### Governing Equations and Boundary Condition  
The calculation domain adopts the \(k\)-\(\varepsilon\) turbulence model based on the Navier-Stokes equations. The air around the brake disc is regarded as a Newtonian fluid. The flow follows the laws of conservation of mass, momentum, and energy. The equation describing the control principle is as follows [16]. Continuity equation:  
\[\frac{\partial}{\partial x_{i}}(\rho u_{i})=0 \tag{4}\]  
Figure 3: Finite element model of brake disc. Momentum equation:  
\[\frac{\partial}{\partial x_{i}}\big{(}\rho u_{i}u_{j}\big{)}=\frac{\partial}{ \partial x_{j}}\Bigg{(}\mu\frac{\partial u_{i}}{\partial x_{j}}\Bigg{)}-\frac{ \partial p}{\partial x_{i}} \tag{5}\]  
Energy equation:  
\[\frac{\partial}{\partial x_{i}}(\rho u_{i}T)=\frac{\partial}{\partial x_{i}} \bigg{(}\frac{\lambda}{c_{p}}\frac{\partial T}{\partial x_{i}}\bigg{)} \tag{6}\]  
where \(p\) and \(\rho\) are the density and pressure, respectively; \(u\), \(\mu\), \(c_{p}\), \(\lambda\), and \(T\) are the velocity, dynamic viscosity, specific heat capacity, thermal conductivity, and temperature, respectively. Standard \(k\)-\(\epsilon\) turbulence model equations  
\[\rho\frac{\partial}{\partial x_{i}}(ku_{i})=\frac{\partial}{\partial x_{i}} \bigg{[}\bigg{(}\mu+\frac{\mu_{t}}{\sigma_{k}}\bigg{)}\frac{\partial k}{ \partial x_{i}}\bigg{]}+G_{k}-\rho\varepsilon \tag{7}\]  
\[\rho\frac{\partial}{\partial x_{i}}(\epsilon u_{i})=\frac{\partial}{\partial x _{i}}\Bigg{[}\bigg{(}\mu+\frac{\mu_{t}}{\sigma_{k}}\bigg{)}\frac{\partial \epsilon}{\partial x_{j}}\Bigg{]}+C_{1\epsilon}\frac{\epsilon}{k}G_{k}-C_{2 \epsilon}\rho\frac{\epsilon^{2}}{k} \tag{8}\]  
where \(G_{k}=\mu_{t}\frac{\partial u_{i}}{\partial x_{i}}\bigg{(}\frac{\partial u_{i} }{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\bigg{)}\), and the term \(\mu_{t}\) represents the turbulent kinematic viscosity, that is, \(\mu_{t}=\rho C_{t}\frac{k^{2}}{\varepsilon}\). C\({}_{\rm L}\), C\({}_{2\epsilon}\), and \(C_{\mu}\) are all constants; \(\sigma_{k}\) represents the turbulent Prandtl numbers for \(\varepsilon\). During the braking process, there is a temperature difference between the contacting brake discs, leading to heat transfer. Since the temperature of the brake discs changes over time, the heat conduction during braking is a three-dimensional transient heat conduction problem. The heat conduction differential equation is as follows:  
\[\frac{\partial T}{\partial\tau}=\frac{\lambda}{\rho c}\bigg{(}\frac{\partial ^{2}T}{\partial x^{2}}+\frac{\partial^{2}T}{\partial y^{2}}+\frac{\partial^{2 }T}{\partial z^{2}}\bigg{)}+\frac{\dot{\Phi}}{\rho c} \tag{9}\]  
where \(\dot{\Phi}\) is the internal heat source. The high-speed rotation of the brake disc causes irregular airflow in the surrounding air, leading to heat transfer between the air and the brake disc surface. This process falls under forced convective heat transfer and follows [PERSON]'s law of cooling. \[\Phi=hA\Big{(}T_{w}-T_{f}\Big{)} \tag{10}\]  
where \(\Phi\) is the heat flow rate; \(A\) is the solid surface area; \(T_{w}\) and \(T_{f}\) are wall temperature and fluid temperature, respectively; and \(h\) is the average convective heat transfer coefficient. In this study, an average heat transfer coefficient was employed to streamline the numerical model and cut computational expenses. Although it offers a fair approximation of overall temperature, it may overlook localized heat transfer variations. Since the temperature range studied has little effect on the physical parameters, constant physical properties are used for calculation. The simulated physical parameters are shown in Table 5. The boundary conditions used in the calculation are as follows:  
1. The rotational angular velocity of the brake disc solid domain and the cylindrical fluid calculation area is equal to the rotational angular velocity of the aircraft axis, and the rotation axis is the \(X\) axis. 2.