## 3 Numerical Methods
### Improved Boundary Condition Enforcement  
The numerical enforcement of the boundary condition plays an important role in the CFD solver development and can significantly affect the computational efficiency and robustness of the marching process in a steady simulation. A commonly used characteristic boundary condition (CBC) for perfect gas involves the so-called \"Riemann-invariant\" at the inlet and adopts linearized characteristic equations at the outlet, of which the details and an example can be found in reference [52]. However, in practice, the writers found that difficulties in convergence were encountered when the Riemann-invariant-based CBCs were enforced in the time-marching throughflow simulation of the centrifugal compressor. Especially for the high-pressure-ratio cases, the computation can blow up even near the peak efficiency operating point. To ensure a stable computation, a method of imposing the CBCs based on [PERSON]'s iteration method was developed for the throughflow method withthe help of the local one-dimensional inviscid (LODI) relations. The developed CBCs refer to the methodology in [53], but were adapted for the throughflow equations. Usually, four values need to be specified at the inlet boundary for the subsonic inflow in the throughflow simulation, therefore a residual vector can be defined as  
\[\mathbf{R}= \begin{bmatrix}R_{1}\\ R_{2}\\ R_{3}\\ R_{4}\end{bmatrix} = \begin{bmatrix}h_{t}-h_{t0}\\ s-s_{0}\\ v-u\tan\sigma\\ w-u\tan\alpha\end{bmatrix} = \begin{bmatrix}dh\\ ds\\ v-u\tan\sigma\\ w-u\tan\alpha\end{bmatrix}  \tag{9}\]  
where \(h_{t0}\) and \(s_{0}\) are the specified total enthalpy and entropy, respectively, and \(h_{t}\) and \(s\) are the corresponding calculated values. \(\sigma\) and \(\alpha\) denote the radial and tangential flow angle. The residual vector \(\mathbf{R}\) is dependent on the primitive variables at the boundary and should vanish after the computational convergence. Using [PERSON]'s method:  
\[\mathbf{R}+\frac{\partial\mathbf{R}}{\partial\mathbf{U}_{p}}d\mathbf{U}_{p}=0 \ \Rightarrow\ \frac{\partial\mathbf{R}}{\partial\mathbf{U}_{p}}d\mathbf{U}_{p}=-\mathbf{R} \tag{10}\]  
\(\mathbf{U}_{p}\) represents the primitive variable \(\mathbf{U}_{p}= \begin{bmatrix}p&u&v&w&T\end{bmatrix} ^{T}\). According to the characteristic theory, another characteristic variable propagating outside the inlet is added here as:  
\[dC_{5}=du_{n}-\frac{dp}{\rho c}=0 \tag{11}\]  
where \(u_{n}\) is the local velocity normal to the inlet boundary. Equation (11) represents a left-propagating acoustic wave. Combining Equations (9) and (11), a system can be presented as:  
\[\mathbf{A}_{in}\cdot d\mathbf{U}_{p}=\mathbf{b}_{in} \tag{12}\]  
with  
\[\mathbf{A}_{in}= \begin{bmatrix}0&u&v&w&c_{p}\\ -\frac{R_{g}}{p}&0&0&0&\frac{c_{p}}{T}\\ 0&-\tan\sigma&1&0&0\\ 0&-\tan\alpha&0&1&0\\ 1&-\rho cn_{z}&-\rho cn_{r}&-\rho cn_{\varphi}&0\end{bmatrix} ,\mathbf{b}_{in}= \begin{bmatrix}-R_{1}\\ -R_{2}\\ -R_{3}\\ -R_{4}\\ 0\end{bmatrix} \]  
The analytical solution of Equation (12) can be obtained as  
\[dp=\frac{\rho dH-\frac{pds}{R}-\frac{\rho us_{ln}}{\cos\delta}+\rho\left(u^{2 }+v^{2}+w^{2}\right)}{1+\frac{u_{s}}{c\cos\delta}} \tag{13}\]  
\[d\rho=\frac{\rho}{\gamma}\frac{dp}{p}-\frac{\rho ds}{c_{p}} \tag{14}\]  
\[du_{n}=\frac{dp}{\rho c} \tag{15}\]  
\(\mathbf{s}\left(\frac{1}{A},\frac{\tan\sigma}{A},\frac{\tan\alpha}{A}\right)\) denotes the unit vector of the prescribed flow direction at the inlet boundary. Where \(A=\sqrt{\tan^{2}\sigma+\tan^{2}\alpha+1}\), \(u_{s}=\mathbf{W}\cdot\mathbf{s}\), \(\cos\delta=\mathbf{n}\cdot\mathbf{s}\). For the subsonic outlet boundary, one characteristic variable propagates into the computational domain and the other four propagate outside. In the time-marching throughflow, the static pressure is usually imposed at the outlet boundary, and thus can be chosen as the inward propagating characteristic  
\[R_{5}=p-p_{b} \tag{16}\]where \(p_{b}\) denotes the specified backpressure.