*3.1. Governing Equations and Numerical Schemes*

The incompressible Navier–Stokes equations were discretised and solved for a computational mesh using the finite volume method [23,24]. The *k*-*ω* SST-SAS eddy viscosity model was employed to take into account the effects of unresolved turbulence. The SAS modifications to the standard *k*-*ω* SST turbulence model allow for a decrease in the turbulent viscosity where the turbulence may be partially or fully resolved [25]. This turbulence model has recently been successfully used in a number of studies concerning transient operations in a hydropower context [14,15,26]. The simulations were carried out on the full 3D computational domain shown in Figure 1b. Each runner domain rotated individually with a solid body rotation, employing an Arbitrary Mesh Interface (AMI) for the transfer of fluxes at the sliding interfaces [27,28].

The convection terms of the momentum equations were discretised using the secondorder accurate Linear-Upwind Stabilised Transport (LUST) scheme [29]. This is a blended scheme, utilising a 75% central differencing scheme (linear), for accuracy, and a 25% second-order upwind scheme (linearUpwind), for stability [30]. The first-order upwind scheme was used to discretise the convection terms in the transport equations for turbulent kinetic energy (*k*) and specific turbulence dissipation rate (*ω*). The linear scheme with the Gauss theorem was employed for the gradient schemes of all variables. The implicit, second-order accurate scheme (backward) was used for temporal discretisation. A constant time step of <sup>Δ</sup>*<sup>t</sup>* = <sup>5</sup> × <sup>10</sup>−<sup>5</sup> <sup>s</sup> was used. At the design point, in pump mode, the time step corresponded to a Courant number of less than 2.2 in 90% of the numerical domain, and the average Courant number was less than 0.065. The maximum runner rotation was 0.44° per time step.

The pressure–velocity coupling was handled using the PIMPLE algorithm [30]. The algorithm allows the Courant number to exceed unity by combining the PISO and SIMPLE algorithms. The solver algorithm in the present work was configured with two inner PISO loops and a maximum of ten outer SIMPLE corrector loops. One non-orthogonal corrector step was used within each inner loop. The PIMPLE algorithm interrupts the outer corrector loops if convergence is reached within each time step. The solution was deemed as converged if the absolute tolerance was below 10−<sup>6</sup> and 10−<sup>5</sup> or the relative tolerance was below 0.1 and 0.01, for velocity and pressure, respectively. The absolute tolerance is the solver residual, and the relative tolerance is in relation to the initial solver residual. The solver algorithm always converged within three or four outer corrector loops.

### *3.2. Boundary Conditions*

A total pressure boundary condition was used for the inlet, while a static pressure boundary condition was imposed on the outlet boundary. The flow rate was thus calculated as a part of the solution. A total-to-static pressure difference of 100 kPa and 134 kPa was used in turbine and pump mode, respectively. The total-to-static pressure difference over the computational domain ensured the net head of the CRPT at the design point, presented in Table 1. Note that the net head in the table was defined by the total pressure difference over just the CRPT, and the total-to-static pressure difference concerned the entire computational domain. In OpenFOAM-v1912, the totalPressure boundary condition was used for pressure at both the inlet and the outlet, together with the velocity boundary condition pressureInletOutletVelocity. This combination of boundary conditions has recently been used to accurately simulate a shutdown transient for a high-head Francis model turbine by Uppström et al. [31]. The totalPressure boundary condition corresponds to a static pressure formulation on the boundary as:

$$p\_{\text{boundary},i} = p\_0 - \frac{\rho |\vec{u}\_i|^2}{2}, \text{for inflow (at face } i\text{)},\tag{4}$$

$$p\_{\text{boundary},i} = p\_{0\prime} \text{ for outflow (at face } i\text{)}.\tag{5}$$

Here, *p*<sup>0</sup> is the user-specified pressure (total for inflow and static for outflow, hence the minus in Equation (4)), *ρ* is the fluid density, and *u* is the velocity vector. Index *i* corresponds to face *i* on the boundary, meaning that the conditions are applied on a face-byface basis. The boundary condition imposes a constant total pressure at the inlet faces and constant static pressure at the outlet faces. The pressureInletOutletVelocity boundary condition for the velocity applies a homogeneous Neumann (zero gradient) condition on the outlet faces. At the inlet faces, the velocity is obtained from the internal field in the face normal direction, as described by Fahlbeck [32]. It was essential to use a pressure–velocity boundary condition combination that allows for reverse flow, since the flow direction changes during the shutdown and startup operations in pump mode.
