*2.1. Definition of Structural Simulations*

As mentioned above, the structural solver was used to determine the deformation of the blades for each timestep. One structural model, whose scheme is presented in Figure 3, was designed to perform a series of simulations. Different eccentricity magnitudes and various angular positions of the eccentricity line with respect to the incoming wind were attained by changing the position of the guide ring center (red X) with respect to the rotor shaft (black X) as indicated by white arrows. That position was fixed in particular simulations.

**Figure 3.** Turbine motion constrains in the structural simulation.

The kinematics of the blades (blue and green arcs) was restricted by several constraints indicated schematically in Figure 3. The turbine shaft rotated around the axis indicated by a black X, where the blades were tangent. Their outer tips marked with green and blue points were constrained by frictionless movement along the guide ring (red circle) as indicated by green and blue arrows. The torque generated by the blades was transferred to the shaft by rode–slider mechanisms. The rode marked with a purple dashed line was fixed to the turbine shaft. The sliders, marked as blue and green cylinders, were connected to the blade tips. The rode–slider mechanisms allowed for frictionless linear movement of the blade outer tips with respect to the axis of rotation, but the mechanisms kept them in the same relative angular position during the turbine rotation. Additionally, the blade outer tips could rotate without friction with respect to the rods as indicated by yellow arrows.

The blade was made of structural steel from the standard ANSYS material library, with the Young modulus and Poisson's ratio equal to 210 GPa and 0.3, respectively. The material was selected to limit deformation to the elastic range to ensure stability of the structural solution for selected eccentricity magnitudes.

#### *2.2. Definition of Flow Simulations*

In comparison to typical simulations of Savonius turbines with blades of fixed shapes (where the internal domain rotates), deformations of the blades enforced an application of an advanced meshing approach in that fluid flow analysis. In the case under study, both flow domains are in the stationary frame of reference.

The computational mesh was generated in ANSYS Meshing. Following the mesh dependence test, described in the next subsection, a mesh composed of 1.16 million control volumes was used in the simulations. The external domain did not change during the computational campaign, thus, it was meshed with a single layer of hexahedral elements through the slice. Refinements were applied around the interface with the internal domain and in the wake downstream of the turbine (Figure 4). The internal domain consisted of two regions marked in green and in red. Due to the requirements of the ANSYS Fluent concerning the dynamic mesh options of deformation and remeshing, tetrahedral elements were used in those regions. In order to ensure a high-quality mesh in the region around the blades (marked in red), a highly refined mesh with 12 layers of prismatic elements at the wall was generated. That number of layers was sufficient to satisfy the condition *y*+ < 1 of the first mesh element at the wall for almost all the simulations. The value was exceeded only for very limited regions at the blade tips. Thus, the mesh was sufficient to solve fully the flow in the boundary layer, which is very important as the flow is characterized by numerous boundary layer separations, especially at convex

sides of the blades. Refinement of the mesh around the blades (especially at their tips) caused the number of elements through its width to be increased to two and even four elements at the blade tips in this region of the domain. In order to avoid numerical errors, the sizes of elements at the interface between the outer and inner domains did not differ significantly.

**Figure 4.** Computational mesh in the external and internal domains with details of its refinement in the blade vicinity presented in the successive magnifications.

The problem of rotation of turbine blades with a constant angular velocity as well as their deformation due to eccentricity was addressed by deforming, remeshing and smoothing algorithms of the dynamic mesh offered by ANSYS Fluent [37]. The instantaneous blade geometry was transferred from the structural analysis every timestep. In the region marked in red, the mesh was following the blade and the elements were deforming as the blade was changing its shape. In the region marked in green, the continuous motion of the blade caused deformation of mesh elements and remeshing was launched wherever they degenerated and the quality measures were not satisfied, namely, the minimal and maximal length scale (1 and 5 mm) with the maximum face skewness of 0.6 and maximum cell skewness of 0.8 were set. The highest skewness of the grid cells was in the region of the deformed and remeshed grid, with the maximal value of 0.84. In the remaining regions, the skewness was below 0.6. The maximal cell aspect ratio at the blades was equal to 68 and was almost constant during the solution as the deformation of the prismatic elements was limited.

The flow simulations were performed in a transient mode with the pressure-based ANSYS Fluent solver. Reynolds-averaged Navier–Stokes (RANS) equations were calculated with the k-ω SST (Shear Stress Transport) turbulence model of Menter [38], which is one of the most frequently used in Savonius turbine simulations [16,39,40], due to its good performance in the adverse pressure gradient and separated flows. Despite the fact that changes in air density are negligible, the fluid was defined as compressible according to the ideal gas law, in order to compensate for numerical instabilities resulting from mesh deformation and remeshing.

A set of the boundary conditions applied in the simulations is indicated in Figure 2. At the inlet, the velocity of 4 m/s was specified with 1% of the turbulence intensity and the turbulent viscosity ratio equal to 1. This velocity corresponds to the Reynolds number based on the turbine diameter equal to 2.5 <sup>×</sup> 105. The pressure outlet condition (absolute pressure of 101,325 Pa) was defined at the opposite end of the domain. At the upper and lower surfaces of the domain in Figure 2, a free-slip wall condition was imposed. As the task was solved in a quasi-2D way, the symmetry condition was defined on the side surfaces, perpendicular to the axis of turbine rotation. The motion of the blades was governed in the structural simulation and was transferred to the Fluent solver by system coupling. The blades of the turbine were defined as non-slip and smooth walls. The interface between both domains was applied to exchange data. The initial conditions were determined on the basis of the boundary conditions.

The SIMPLE (semi-implicit method for pressure-linked equations) pressure-velocity coupling method was used. The second-order spatial discretization schemes were applied for mass, momentum and energy equations. A first-order implicit transient scheme was also used. Maximally 20 iteration loops were solved in every timestep of transient simulations. Typically, it was enough to reach the residual target of 1e-3, selected as the convergence criterion for all equations solved in Fluent. This criterion was not reached only for a few timesteps, due to numerical problems resulting from mesh deformation or remeshing.

## *2.3. Validation of the Numerical Procedure*

In order to assess the precision of the numerical method with mesh deformation and remeshing used in this investigation, it was decided to compare its results with the results of a typical simulation of the Savonius turbine, applied, e.g., in [13,16,17]. It was possible to do it for the eccentricity magnitude equal to zero, i.e., when rotor blades were not deformed. In this case, the flow was additionally simulated with a fixed mesh, where continuous changes of the rotor position were obtained by rotation of the circular inner domain and exchange of data at the sliding mesh interface between the stationary and rotational domains. The rest of the simulation conditions was preserved.

Taking advantage of much faster simulations for such a task arrangement, the mesh size dependence was also verified by means of Richardson's extrapolation, similarly as in [41]. The same task was solved on three different grids with a mesh refinement ratio of 2. The numbers of control volumes of the grids and the average power coefficient values obtained for them are presented in Table 1. According to the Richardson's extrapolation procedure described in [35], the extrapolated value (*RE*), the apparent order *p*, the ratio of error *R* and the fine-grid convergence index (*FGCI*) were determined. The negative value of R indicates an oscillatory convergence. Only a slight variation of the average power coefficient was observed (1.2%), however, the *y*+ < 1 condition was fulfilled for the mesh composed of 1.16 million control volumes shown in Figure 4. The uncertainty due to discretization determined for the average *Cp* on the basis of the Richardson's extrapolation procedure described in [42] for this mesh was low, i.e., 1.7%. Thus, it was decided to use it in the further numerical investigations.


**Table 1.** Mesh size dependence study.

In the cases of typical 2D simulations of fixed-shape Savonius rotors within the RANS method [16,17], after simulations of a few or more of rotor revolutions, one can reach the limit cycle of the blade load. In this case, a change in the period-averaged value of the power coefficient was lower than 1% of its

value after 7 revolutions. The same criterion for this coefficient for particular instants of the rotor revolution (limit cycle) was reached after 10 revolutions. This limit cycle is marked in Figure 5 with a dashed line.

**Figure 5.** Comparison of the power coefficient *Cp* for one revolution of the non-deformable rotor from simulations with the rotating inner domain and the deformed/remeshed one. Subsequent positions of the turbine illustrated in the schemes above the graph correspond to the turbine revolution angle.

In the case of the simulations with mesh deformation and remeshing applied in the further investigations, such a limit cycle was not obtained, mostly due to numerical instabilities. They resulted from the continuous deformation of mesh control volumes, which successively deteriorated in quality up to the moment when the specific regions were remeshed. Remeshing introduced some further numerical errors due to the data interpolation onto a new mesh. Thus, slight random fluctuations of the blade load were present. Therefore, the data presented in Figure 5 with a solid line were averaged for the last 5 out of 15 rotor revolutions. This number of periods was sufficient to reach no significant changes of the data averaged for the same angular positions during the rotor revolution. The period-averaged values of the power coefficient changes were lower than 1% after 10 revolutions.

As one can see, no significant differences can be observed for two parts of the period for the fixed-mesh simulations (dashed line in Figure 5). In the case of simulations with mesh deformation and remeshing (solid line), differences for two parts of the period are significant. They are observed especially for the revolution angles where the rotor performance is the lowest (75–120◦ and 255–300◦). In these particular ranges, random fluctuations of the blade loads were the highest. The most probable explanation is a very complex flow pattern with numerous vortex structures at the convex side of the advancing blade. The vortex in particular at its tip is very strong with high velocity gradients and, thus, it can be vulnerable to some fluctuations due to the mesh quality deterioration or the remeshing procedure. The poorer flow prediction in those ranges caused that the period-averaged value of the power coefficient was 4.6% higher with respect to the fixed mesh (*Cp* = 0.206 for remeshed, *Cp* = 0.197 for fixed-mesh simulations). This difference is definitely not negligible, but it is an order of magnitude lower than the pressure coefficient increase predicted for rotors with deformable blades. Thus, despite

its drawbacks, the method with mesh deformation and remeshing can be considered as sufficiently reliable for the needs of these investigations.
