## 3 Numerical Procedure  
The current computational model and numerical scheme are based on those for hydrodynamic TPJBs considering pivot stiffness in Ref. [28], which is extended further to consider hydrostatic pressurization at the recess region, flow turbulence, and pressure/temperature-dependent thermochemical properties of the lubricant (sCO\({}_{2}\)). The finite element method solves the Reynolds equation for pressure, the finite difference method solves the energy equation for the temperature distribution, and the Newton-Raphson method determines the equilibrium position. Figure 2 shows a flowchart of the numerical procedure. Each step in Figure 2 is briefly explained as follows:  
1. The bearing geometry, lubricant properties at the initial condition, the initial guesses for the pad tilting angle \(\delta\), pivot deflection \(\zeta\), eccentricity ratio \(\epsilon\), the angle between the \(X\)-axis and the line of the centers \(\theta_{c}\), and the pressure ratios \(p_{ratio}\) (\(p_{re}/p_{s}\)), \(p_{ratio}\) (\(p_{de}/p_{s}\)), \(p_{ratio}\) (\(p_{da}/p_{s}\)) are provided as the input. For the initial guess for \(\delta\), Equation (4) assumes zero film thickness at the trailing edge of a pad, which indicates the maximum \(\delta\). Half of the maximum \(\delta\) is then taken as the initial guess for \(\delta\) of each pad. The initial guess for \(\zeta\) is 0.1% of the bearing clearance. Similarly, the initial guess for \(\epsilon\), \(p_{ratio}\), \(p_{ratio}\), and \(p_{ratio}\) must be positive and less than 1, and \(\theta_{c}\) is guessed to be 90\({}^{\circ}\), indicating that the rotor is displaced only in the vertical direction. 2. The finite element method solves the Reynolds equation to compute the pressure distribution over each pad. The generated hydrodynamic pressure exerts a force \(F_{pad}\) on the surface of the pad, which deflects the supporting pivot and gives rise to a restoring force \(F_{p}\) in it [28]. As shown in Equation (13), TPJB with the number of pads \(N_{pad}\) in the static equilibrium requires that \(F_{p}\) balance \(F_{pad}\) all the moments and bearing horizontal force \(F_{X}\) balance to zero, the bearing vertical force \(F_{Y}\) balance the static load \(F_{p0}\), and the orifice flow rate \(Q_{o}\) balance the summation of the flow rates at  
Figure 2: Flowchart of the numerical procedure. the recess boundaries \(Q_{r}\). The circumferential and axial pressure drops at the recess edges due to inertia are also calculated iteratively as follows:  
\[\left\{ \begin{array}{l}\left[F_{pdd}-F_{p}\right]_{n=1\ldots N_{pad}}\\ M_{1\ldots N_{pad}}\\ \sum F_{X}\\ \sum F_{Y}-F_{p0}\\ \left[Q_{r}+Q_{o}\right]_{n=1\ldots N_{pad}}\\ \left[p_{dc_{i}}-p_{dc_{i-1}}\right]_{n=1\ldots N_{pad}}\\ \left[p_{da_{i}}-p_{da_{i-1}}\right]_{n=1\ldots N_{pad}}\\ \end{array} \right\}=\left\{ \begin{array}{l}0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ \end{array} \right\} \tag{13}\]  
where \(i\) and \(i-1\) denote the current and previous iterations, respectively. The Newton-Raphson method calculates the equilibrium position. 3. Once the static equilibrium position converges, the thermal energy mixing model calculates the fluid's inlet temperature at the pad's leading edge. The control-volume finite difference technique with an upwind scheme solves the two-dimensional energy equation (Equation (9)) for the temperature distribution (see Ref. [34] for more details). 4. The converged pressure and temperature are then used to calculate the density and viscosity. 5. Once the density and viscosity are converged iteratively using the successive substitution method, the journal is perturbed to a new position, and static equilibrium is achieved again. 6. The stiffness and damping coefficients are calculated from the ratio of the difference between the forces at the new and old equilibrium positions to the difference in distance and velocity, respectively, at both points.