A FEM Structural Analysis of a Francis Turbine Blade Parametrized Using Piecewise Bernstein Polynomials

Abstract

Several methodologies have successfully described the runner blade shape as a set of discrete sections joining the hub and shroud, defined by 3D geometrical forms of considerable complexity. This task requires an appropriate parametric approach for its accurate reconstruction. Among them, piecewise Bernstein polynomials have been used to create parametrizations of twisted runner blades by extracting some cross-sectional hydrofoil profiles from reference CAD data to be approximated by such polynomials. Using the interpolating polynomial coefficients as parameters, more profiles are generated by Lagrangian techniques. The generated profiles are then stacked along the spanwise direction of the blade via transfinite interpolation to obtain a smooth and continuous representation of the reference blade. This versatile approach makes the description of a range of different blade shapes possible within the required accuracy and, furthermore, the design of new blade shapes. However, even though it is possible to redefine new blade shapes using the aforementioned parametrization, a remaining question is whether the parametrized blades are suitable as a replacement for the currently used ones. In order to assess the mechanical feasibility of the new shapes, several stages of analysis are required. In this paper, bearing in mind the standard hydraulic test conditions of the hydrofoil test case of the Norwegian Hydropower Center, we present a structural stress–strain analysis of the reparametrization of a Francis blade, thus showing its adequate computational performance in two model tests.

1. Introduction

One of the most critical aspects of turbomachinery design is to find the optimal shape of the turbine blades to reach their highest possible performance [1,2]. Once an actual optimized blade has been operating for a while, it is necessary to take corrective measures after the natural wear of the blades due to their operating conditions and the consequent reduction in the turbine efficiency. This could mean either replacing the whole blade or else to redesigning it by reconstructing it through the use of parametrization tools. Nevertheless, the high costs of the manufacturing process and materials required in the turbine blade construction make the reconstruction process a more economical solution [3].

An essential step in the reconstruction methodology is the implementation of a parametrization of the blade geometry using a low number of parameters, which defines a suitable design space for finding optimal solutions [4]. Blade parametrization methodologies can be classified into two main groups: those that use a sectional approach and those that use surface patches [5,6,7]. The sectional approaches are useful to describe strong shape variations; however, a high number of sections required to define a twisted blade shape involve a large number of parameters [1]. The 3D radial–axial blades of a Francis turbine are such twisted shapes that they can be considered free-form surfaces. This process increments the complexity of accurately adjusting the blade shape by interpolation using one or various analytic functions [8]. Once a numerical model is developed, it is necessary to assess the mechanical technical feasibility of the blade, employing a structural analysis to determine the state of stresses and strains under standard working conditions [9,10,11,12,13,14]. One of the most used methods for the analysis is the finite element method (FEM), a versatile technique that can handle complex shapes in problem domains [15,16,17]. In the stress–strain analysis, all the effects produced by the acting forces must be considered; however, the primary and most studied effect when the whole runner is considered is the centrifugal force formed by high rotational speeds [18]. Its magnitude depends on the rotor’s size and speed [19]. Suha et al. [19] presented a static structural analysis of a turbine blade considering different alloys with Inconel^© and titanium as the principal materials. The results showed a stress–density dependence, with lower stress values in the titanium alloy.

On the other hand, Kauss et al. [20] developed a similar mechanical structural analysis that included thermal loads and studied a turbine blade model to estimate the stress–strain condition using different alloys. They considered the effect of the temperature by breaking up the strain tensor into mechanical and thermal components based on the temperature changes. As a result, it was observed a clear dependence of the stress–strain values on the physical properties of the materials and a predominance of the thermal component over the mechanical stresses and strains. Fatigue life prediction is one of the most important turbine research fields, and there have been many advances in this area [21,22,23]. Claudio et al. [24] applied FEM to the life prediction of a gas turbine disk and compared their results with experimental data obtained from real test discs, obtaining accurate results. Some of the turbine research is focused on computational fluid dynamics (CFD) analysis and the determination of the turbine (or blade) pressure distribution under different flow conditions to evaluate the acting forces and stress–strain state of the element [25,26,27].

Considering more demanding conditions, another alternative for structural analysis is the use of the extended finite element method (XFEM), in which crack generation can be modeled utilizing discontinuous functions to represent the two-dimensional asymptotic crack-tip displacement fields. The XFEM methodology considers the presence of flaws or inhomogeneities during the mesh generation, employing a mesh enrichment in all the nodes whose nodal shape function intersects the interior of the crack [28]. Lukas et al. [29] implemented an integrated methodology for analyzing turbine blades, in which they incorporated aerodynamic and multiple structural simulations in a surface geometry scanned from actual blades. The XFEM was applied to demonstrate the influence of the cracks on the leading edge (LE) behavior and to predict the influence of geometry variations due to wear.

In summary, due to its importance, several structural studies have been carried out to assess the performance of blades [30,31,32,33,34,35]. This paper aims to analyze the structural performance of a twisted Francis 99 turbine blade parametrized by piecewise Berstein interpolation as described by Pérez et al. [3]. FEM was used to determine its stress–strain state in two numerical experiments similar to those of the hydrofoil test case of the Norwegian Hydropower Center [36,37,38,39,40], which were used to determine the form of the pressure load distribution over the blade surface. Due to the structural performance obtained, the results suggest that the parametrized model could be used for CAD design and, eventually, as a previous step in the reconstruction process of a whole turbine. However, more intensive mechanical structural analyses are required for the design of the runner and will be reported in a future paper.

2. Materials and Methods

2.1. Blade Description

The Francis 99 turbine has been used as a study object by the Norwegian University of Science and Technology (NTNU) and the Luleå University of Technology (LTU) (Figure 1). Both institutions have designed, built, and experimentally investigated an axial–radial-type Francis turbine on a test bed. This turbine is a prototype reduced to a 1:5.1 scale of the turbines that operate at the Tokke hydroelectric power station in Norway. Through a series of workshops called Francis 99, open access to data related to different test cases (including CAD models, meshes for simulations, and experimental data) was provided for the turbine and hydrofoils [21]. According to the laboratory tests at the best efficiency point (BEP), the Francis 99 turbine obtained a runner speed of 335.4 rpm, hydraulic efficiency of 92.61%, and generator input torque of 619.56 Nm (torque was measured in between the thrust block and the generator).

The turbine runner consists of 15 main blades and 15 splitter ones (Figure 2a). Figure 2b offers a cross-sectional view of the runner, its main components, and the fluid inlet and outlet. The runner inlet height is 0.06 m, and the runner inlet and outlet diameters are 0.63 m and 0.349 m, respectively.

Figure 3 shows the main runner blade used in this work. It consists of a pressure side (PS) and a suction side (SS), which play an essential role in energy transformation. In addition, the blade has a leading edge (LE), through which the fluid enters with a higher pressure and a trailing edge (TE), through which it leaves with a lower pressure. The main blade is twisted around 180 degrees from the inlet to the outlet of the runner, and the blade thickness at the trailing edge is 3 mm.

2.2. Reconstruction Methodology

The reconstruction methodology process can be summarized in three main stages: data extraction, blade reconstruction, and numerical blade evaluation (see Figure 4) [1,3,8].

2.2.1. Data Extraction Stage

In the data extraction stage, the reconstruction process begins with the extraction of the coordinate points $(x_i,y_i,z_i)$ from the Francis 99 turbine blade workshop data as proportioned by the Norwegian Hydropower Center [41]. Due to the low amount of scanned data on some blade zones, it was required to redefine the trailing and leading edges at the hub and shroud cross-sections using an algorithm based on the local curvature of the interpolating polynomials. Once the LE and TE were redefined, the mean camber line was calculated as a spatial average of the pressure and suction sides. Finally, ten cutting blade sections were calculated as described in the following subsection.

2.2.2. Blade Reconstruction Stage

Each blade section was reconstructed with a piecewise fourth-degree Bernstein polynomial, a technique which has been previously applied to obtain accurate results, see [3]. This procedure allows us to define the 3D camber surface as a series of orthogonal tensor products on the transfinite interpolation of the known boundaries of the camber surface. One of the advantages of this reconstruction methodology is that the blade adjustment is relatively independent of its design parameters, which is the key to using a rather low number of the latter to generate the blade surface. In order to reduce the complexity given by the twisted shape in the section curves, both faces are divided into three or four polynomial pieces to improve the approximation.

Each fourth-degree Bernstein polynomial piece has the form [1,3,42]

$$P{\left(x\right)}^k=\sum _{i=0}^{n=4}{C}_i^{\left(k\right)}{B}_i^4\left(x\right)$$

(1)

where $C_i^{\left(k\right)},i=0,...,4$ are five control points, which act as the parameters to describe the shape of the blade section. An interpolation condition is defined at each piece’s beginning and ending points, which leaves three control points free; the latter ones are set by minimizing the mean square error to the scanned blade points on the corresponding section, thus determining the computational blade cross-section parametrization (Figure 5). Next, the 3D shape of the turbine blade is created by heaping up ten surfaces throughout the spanwise using an orthogonal 3D profiling technique [3,43]. Since these sections are piled up along the spanwise of the camber skeleton, using a Lagrangian interpolation, it is possible to generate more blade sections if required to render a smooth blade surface (Figure 6). A more detailed discussion of this reconstruction process can be found in [3].

2.2.3. Numerical Blade Evaluation

The mean square error (Equation (2)) evaluates the fit between the original and parametrized blades:

$$SME=\sqrt{\frac{{\sum }_{i=1}^n{D}_i^2}{n}}$$

(2)

where $D_i$, defined by Equation (3), is the Euclidean distance between a blade point $(x_i$, $y_i$, $z_i)$ to the tangent surface defined by the unit normal vector d (Equation (4)) on the blade surface (see Delgado et al. [8]).

$$D_i=|cos\left(\alpha \right)x_i+cos\left(\beta \right)y_i+cos\left(\gamma \right)z_i+d|$$

(3)

$$d=cos\left(\alpha \right)x+cos\left(\beta \right)y+cos\left(\gamma \right)z$$

(4)

The value of the mean square error between the data and the parametrized blade used in this paper is $1.5496\times {10}^{-3}\%$ on the pressure side, and $9.685\times {10}^{-5}\%$ on the suction side. A thorough discussion on this subject, as well as more reconstruction and evaluation details, can be found in the work of Pérez et al. [3].

2.3. Structural Model

The relation between the divergence of the stress tensor $\mathbf{\sigma }$ and the acting forces $\mathit{f}$ gives the governing equation that defines the structural problem of interest:

$$\nabla ·\mathit{\sigma }+\mathit{f}=0.$$

(5)

The key difference between different models is the constitutive law for the stress–strain relation. In the theory of plasticity, several elastic and plastic stress–strain relations for solids have been studied [44,45]. The most common plasticity models implemented in engineering analysis are those of Tresca, Von Mises, Mohr–Coulomb, and Dracker–Prager [46]. In this paper, a small-strain isotropic Von Mises plasticity constitutive law was considered, and a Von Mises yield and a Von Mises potential criterion was applied. This plasticity model was recently used in the analysis of reinforced concrete structures strengthened with post-tensioned or pre-stressed cable tendons [47], showing interesting results in the plasticity behavior of the material interaction. This model was also employed in the failure analysis of a nuclear power plant reactor subjected to accidental internal pressures, obtaining very satisfactory results [48].

The basic components of the elastoplastic constitutive models could be summarized as follows [46]:

• The elastoplastic strain decomposition (Equation (6));
• An elastic law (Equation (7));
• A yield criterion, and a yield function (Equation (8));
• A plastic flow rule (Equation (9));
• A hardening law (Equation (10)).

In this case, according to the additive decomposition rule [49], the strain tensor is the sum of its elastic part (${\mathit{\epsilon }}^{\mathit{e}}$) and plastic part (${\mathit{\epsilon }}^{\mathit{p}}$) defined by

$$\mathit{\epsilon }={\mathit{\epsilon }}^{\mathit{e}}+{\mathit{\epsilon }}^{\mathit{p}}$$

(6)

and the elastic law as

$$\mathit{\sigma }=\mathit{C}:{\mathit{\epsilon }}^{\mathit{e}}.$$

(7)

The yield function for the Von Mises plasticity model is given by

$$\Phi \left(\sigma \right)=\sqrt{J_2\left(S\left(\sigma \right)\right)}-{\tau }_y,$$

(8)

where $J_2$ is the deviator invariant and ${\tau }_y$ is the shear yield stress. The corresponding flow rule and hardening law are described by

$${\dot{\epsilon }}^p=\gamma \mathit{r}(\mathit{\sigma },\mathit{q}),$$

(9)

$${\dot{\mathit{q}}}^{\mathit{p}}=-\gamma \mathit{h}(\mathit{\sigma },\mathit{q}),$$

(10)

respectively.

All these equations are subject to the Kuhn–Tucker complementary conditions [44,46].

2.3.1. Materials

For the test case, a Francis 99 turbine blade made of an aluminum alloy, whose technical designation is Aluminum 5456-H24, was considered. The material properties taken from Elsherif et al. [50] are described in

Table 1. Aluminum alloy material properties.

Property	Magnitud	Unit
Young’s Modulus	71,000	MPa
Poisson’s Ratio	0.33
Density	2770	kg/m ${}^3$
Yield Strength	280	MPa
Tensile Ultimate Strength	310	MPa
Bulk Modulus	69,608	MPa
Shear Modulus	26,692	MPa
Ultimate Bearing Strength	669	MPa

.

2.3.2. Test Case

The test consists of a computer-simulated approximation to the hydrofoil test case developed by the Norwegian Hydropower Center (NHC). That experiment was performed in a long piping loop with a segment with a square section of $0.15$ m $\times 0.15$ m internal dimension. The hydrofoil geometry was $0.012$ m thick and had a cord length of $0.25$ m. The boundary conditions were defined according to the NHC hydrofoil test. The blade was fixed in the internal side walls of the square section in such a way that no rotation was allowed. The flow is axial in each time step. The test was carried out with a flow water velocity ranging from 0 to 25 m/s, with velocity steps of 5 m/s (Figure 7) [41]. This maximum velocity value was used, even though it could be considered critical since it exceeds the standard working conditions for a runner [51].

In the present analysis, to obtain the pressure distribution on the element, the flow rate was also increased linearly from 0 to 25 m/s (Figure 8). The pipeline computer model was designed according to the geometry of the parametrized blade. Here, two cases were considered:

• Model 1. A blade with constant cross section generated by extrusion from the Francis blade profile. The pipeline geometry is $0.80$ m long and has a square section of $0.15$ m × 0.15 m (Figure 9a). This geometry resembles the original test model.
• Model 2. A Francis turbine blade was created using the reconstruction methodology described in the preceding section. For this model, the pipeline geometry was adjusted to wrap the blade Francis geometry, thus having a length of $0.80$ m and a rectangular section of $0.20$ m × 0.22 m (Figure 9b).

Kratos Multiphysics^©, a modular software developed at CIMNE, was used for the CFD setting of this paper. Kratos is a working environment focused on the implementation of numerical methods for the solution of multiphysics problems. For both models, the 3D flow module was employed to generate a numerical solution of the Navier–Stokes equations using the monolithic approach in conjunction with the Bossak method. A large eddy simulation (LES) model was considered for turbulence [52,53,54,55].

GiD v16.0.3, also developed by CIMNE, was used with Kratos as a problem type: Kratos being the solver, and GiD the interpreter, pre-processor, and post-processor of the data [56]. Through GiD, the discretization of the volume enclosed between the pipe and the blade was carried out for both models, generating a mesh of triangular-shaped adaptive elements on the surfaces and one of tetrahedrons for the volumes. In model 1, there were 22,060 triangular elements and 54,032 tetrahedra; in model 2, there were 12,026 triangles and 55,950 tetrahedra.

In both models, the pipe was considered to be initially empty. Figure 10 shows the setting of the boundary conditions for each model. As mentioned, as in the hydrofoil test of the NTNU, the flow velocity at the inlet was increased linearly; a zero pressure was considered at the outlet.

The extremal pressure values are obtained at the blade surface near the LE and TE as expected. It can be seen in model 1 that the particle flow trajectories are modified for the blade geometry, generating some fast velocity changes at the suction side (Figure 11a). Under the same inlet velocity conditions, model 2 has different performance. The maximum pressures occur again at the trailing and leading edges. However, in this case, the particle flow trajectories are strongly altered for the blade geometry, with a remarkable spreading at the middle of the surface face of the turbine blade (Figure 11b).

Table 2. Pressure values on LE.

Time (s)	Hydrofoil Test (kPa)	Model 1 (kPa)	Model 2 (kPa)
0.0	659.5	0.0	0.0
1.0	666.8	13.5	31.7
1.6	675.2	31.1	80.0
2.0	683.7	46.6	123.7
3.0	710.0	97.4	278.5
4.0	747.3	163.9	504.5
5.0	794.2	244.7	773.0

shows pressure values reported by NHC and the corresponding values from the CFD setting. The pressure values on the table were measured on the LE. As the initial condition, the pipeline in the hydrofoil test had a pressure value of $659.5$ kPa in order to avoid cavitation; on the other hand, since the setting of interest here is when the blade is subject to the maximum pressure values, the initial condition for pressure was set to zero. In model 1 the maximum differs from that of the hydrofoil test despite the fact that both models appear to have similar geometries. Furthermore, the maximum pressure in model 2 (773 kPa) differs only $2.7\%$ with respect to hydrofoil one ($794.2$ kPa), which is a good approximation for the maximum pressure load distribution applied in the structural analysis.

The hydrodynamic geometry of model 1 increases the flow velocity upon contact with the LE. Due to the area reduction and the upper blade curvature, the maximum velocities are developed near the top surface. However, at the TE, the velocities tend to equal the input ones (Figure 12). The blade shape of model 2 produces more significant velocity variations. It changes the inlet flow from uniform to a rapidly changing one near the middle surface section and close to the blade TE, producing a fast change gradient zone (Figure 13).

The results from the fluid dynamic simulation show a concentration of maximum pressure values on the LE for both models. Model 1 exhibits a gradual decrement in pressure, from positive values on the LE to negatives ones on the TE (Figure 14), while model 2 displays a clear pressure gradient, with positive pressure values mainly at the upper face and negative ones at the lower side (Figure 15).

3. Results

Figure 16a,b show different views of the adaptive blade mesh of model 1, and Figure 16c,d present the blade mesh views of model 2. Both meshes are composed of approximately $70,000$ tetrahedral elements. Again, the structural analysis was carried out in Kratos Multiphysics under the assumptions of Section 2.3, using the pressure values described in the previous section.

The structural analysis was carried out independently using Kratos and taking the CFD pressure values as input. Figure 17, Figure 18, Figure 19, Figure 20, Figure 21, Figure 22, Figure 23, Figure 24, Figure 25 and Figure 26 show the results for both models. Figure 17 and Figure 18 present the Cauchy principal stresses, and Figure 19 and Figure 20 display the Cauchy stress tensor for models 1 and 2, respectively. A stress concentration on the TE is observed in all cases, with maximum values of $92.66$ MPa for model 1 and $340.60$ MPa for model 2, and minima of $-181.36$ MPa and $-656.07$ MPa, respectively. It can be noticed that maximum and minimum stresses for model 1 are lower than the yield stress of the aluminum alloy implemented. This is attributed to the fact that maximum pressures developed in the CFD simulation of model 1 are far from the maximum pressure of the hydrofoil test, but, under the complex geometry of model 2, the maximum values almost exceed the ultimate bearing strength.

The Green–Lagrange principal strains (Figure 21 and Figure 22) and the Green–Lagrange strain tensor (Figure 23 and Figure 24) show similar distributions with an absolute maximum value of $2.2734\times {10}^{-3}$ for model 1, and $1.7484\times {10}^{-2}$ for model 2. Both maximum values are attained at the TE.

Figure 25a,b present the maximum Von Mises stress values and deformations for model 1, respectively. The maximum stress value is $169.80$ MPa, and the maximum deformation is $4.36\times {10}^{-4}$ m. Figure 26a,b are the corresponding Von Mises and deformation representations for model 2. In this case, the maximum Von Mises stress equals the yield limit 280 MPa; the maximum deformation is $1.317\times {10}^{-3}$ m. The figures show that the Von Mises stress values are less than the aluminum alloy yield stress only in model 1.

4. Discussion

The use of Bernstein polynomials for piecewise interpolation as the cornerstone for the surface methodology approximation of the Francis turbine blade reconstruction has shown its capability to generate accurate results and be adapted to parametrize complex blade forms. The remaining task is to assert the mechanical feasibility of the generated blade forms under standard test conditions. A first step in that direction was taken in this paper. The computational fluid dynamics simulation addressed in this paper showed that the velocity gradient change caused by the twisted geometry of model 2 (which represents the shape of an actual Francis blade) increased the pressure, stresses, and strain on its surface compared to those of model 1. As expected, the structural analysis exhibits a stress concentration on the TE for both models. In addition, it must be noted that these stress values might exceed the yield stress of the aluminum alloy considered for the simulations. However, such values lie within the limit given by the ultimate bearing strength. It was also observed that maximum strain values are attained at the TE, which is consistent with the damage observed in real blades under working conditions.

The general overview of the performance of both models describes a suitable reconstruction of the geometry using the proposed reparametrization, which, given the versatility of the interpolation technique applied, becomes a promising methodology, even for more challenging tasks, such as the potentiation of the whole turbine.

5. Conclusions

This research was focused on the structural stress–strain evaluation of a reparametrized Francis 99 blade. Based on the accurate geometrical reconstruction of the blade, as well as the proposed operation parameters, which resemble those of the hydrofoil blade test of the Francis 99 workshop, the outcome suggests the feasibility of applying the new blade shape for an actual setting: under similar test conditions, the stress values attain its maximum values at the trailing edge, and the pressure values, its maximum values at the leading edge. However, in both cases, such values lie within the current material limits of the aluminum alloy considered in the tests.

The results suggest that the analyzed mechanical behavior describes a blade well suited for the operation conditions established in this test. Nevertheless, due to the complex dynamics of turbomachines, several studies for the whole runner for its eventual manufacturing are still required to be carried out. That will be reported in future papers.