*Article* **A Comparison of Physical and Numerical Modeling of Homogenous Isotropic Propeller Blades**

**Luca Savio 1,2,\*,†, Lucia Sileo 1,† and Sigmund Kyrre Ås 1,2**


Received: 10 December 2019; Accepted: 26 December 2019; Published: 3 January 2020

**Abstract:** Results of the fluid-structure co-simulations that were carried out as part of the FleksProp project are presented. The FleksProp project aims to establish better design procedures that take into account the hydroelastic behavior of marine propellers and thrusters. Part of the project is devoted to establishing good validation cases for fluid-structure interaction (FSI) simulations. More specifically, this paper describes the comparison of the numerical computations carried out on three propeller designs that were produced in both a metal and resin variant. The metal version could practically be considered rigid in model scale, while the resin variant would show measurable deformations. Both variants were then tested in open water condition at SINTEF Ocean's towing tank. The tests were carried out at different propeller rotational speeds, advance coefficients, and pitch settings. The computations were carried out using the commercial software STAR-CCM+ and Abaqus. This paper describes briefly the experimental setup and focuses on the numerical setup and the discussion of the results. The simulations agreed well with the experiments; hence, the computational approach has been validated.

**Keywords:** fluid structure interaction; propulsion; CFD

#### **1. Introduction**

The topic of hydroelastic behavior of marine propellers is nowadays often associated with the behavior of composite propellers that aim at using some sort of anisotropy to achieve some set goal. The goal may be reducing vibrations or cavitation (Young [1,2]). Already in 2001, Mouritz et al. [3] published a review article citing a quite extensive list of applications of composite materials in naval applications, including propellers, where the benefits of composite materials could offset the increased costs that are associated with their design and production. However, this paper covers isotropic propellers as a first step towards anisotropic propellers. In fact, as it will be clear further in the paper, the simulation strategy adopted here should in principle work independently whether the material is isotropic or anisotropic. The choice of starting from isotropic materials was mainly made because the production methods needed for making the model scale blades were easier than for anisotropic materials.

The elastic behavior of isotropic blades has received less attention although it may in same cases be noticeable also for full-scale metallic propellers. In the past, Atkinson and Glover [4] have shown numerically how full-scale metallic propellers could bend under the action of hydrodynamic forces to an extent that some particular propeller geometries, namely blades with large skew angles, could have some significance in terms of propeller performance. The numerical method adopted was a combination of lifting line and finite element methods (FEM) in what we would call a weak co-simulation in modern terms. Nowadays, co-simulations are increasingly carried out coupling CFD

(Computational Fluid Dynamics) with FEM (Finite Element Method) code. The coupling between the two codes vary from weak to strong, according to how the two simulations communicate in terms of time increments and variables that are transferred. While this has been a research area for quite some time, modern software implementations have reached a level of maturity that makes these types of analysis attractive also at the level of industrial design engineers. The increased ease of setting up co-simulations should not however be taken as a guarantee for obtaining accurate results. In fact, the accuracy of the computations should be verified and compared with results obtained in controlled scenarios, as, for example, during open water tests in a towing tank or cavitation tunnel. Despite some recent attempts of fill in the gap [5,6], there is still a scarcity of experimental data. The validation of numerical results by means of experiments presents a series of challenges that require careful preparation and interpretation of the results. It is important to limit the scope of what can be validated by means of experiments; it is practically impossible to design experiments that allow for correctly scaling to laboratories size displacements and forces of some full-scale hydroelastic phenomenon. A thorough review of what are the limitations of mode scale experiments can be found in Young [7]. What is possible is to design experiments that allow for their correct representation in the simulation codes with the aim of validating the latter. In this paper, we present the results from experiments carried out on three propeller designs that have been produced with blades both in aluminum and in a casting resin resembling an epoxy resin , hence presenting a rigid and flexible behavior, respectively. The three propellers share the same design parameters apart from the skew distribution. The tests have been performed in a classical open-water configuration at different propeller speeds and at different pitch settings. The open-water configuration was chosen so that the hydroelastic behavior was limited to static deflection, i.e., vibrations of the blades were avoided by performing the tests in a homogenous inflow.

#### **2. Model Tests**

#### *2.1. Propeller Geometries*

Three propeller geometries were generated starting from a master geometry varying the skew distribution. The master geometry is labeled P1374 and was designed using SINTEF Ocean proprietary propeller design code AKPD by Achkinadze et al. as presented in [8,9]. Propeller P1374 has a balanced skew distribution that amounts to 23 degrees. The first variation, P1565, was obtained by removing entirely any skew distribution. In the second variation, P1566, the same total skew as propeller P1374 was kept but in an unbalanced form. The geometries of the blades are available, since the FleksProp consortium decided to allow disclosure of the geometries of the three propellers. They are reported in Appendix A, while the geometrical definition of the propellers are described in Appendix B. A visual impression of how the three propellers look like is given in Figure 1.

**Figure 1.** Outline of the three geometries showing the different skew distributions.

#### *2.2. Production of Flexible Blades*

The flexible blades were manufactured using resin casting under vacuum, which is a commonly used technique in rapid prototyping. There are several variants of the resin casting technique. In the variant that was adopted in this case, the first step is to produce a positive object that is then used to form a silicon mold where the resin is cast under vacuum and left to cure at constant temperature. Three complete controllable pitch propellers, including the hubs, were produced in aluminum according to the production standard used in SINTEF Ocean. One blade per kind was then sent to a company specializing in rapid prototyping that produced 5 copies of the blade.

According to the producer of the resin used for casting, Young's modulus is equal to 2.2 GPa. Poisson's ratio is not reported, but given the nature of the material, it was assumed to be equal to 0.33. The material properties reported by the producer have not been checked by mechanical testing; however, the good agreement between computations and experiments indicates that the reported values can be considered accurate.

#### *2.3. Test Executions*

The model tests have been carried in the SINTEF Ocean large towing tank between the months of February and June 2018. The setup adopted was the classical open-water configuration where the propeller is placed in front of the dynamometer. The dynamometer used in the test was a Kempf & Remmers H29 model. Figure 2 shows that only the surfaces in contact with water located on the left side of the 1-mm gap contribute to the forces measured by the dynamometer.

**Figure 2.** Experimental setup showing which surfaces contribute to the forces measured during the experiments.

According to standard procedures for open-water tests, a dummy hub without blades was tested to obtain forces on the cap and hub at different velocities. The forces measured during the test with the dummy are then subtracted from the forces measured with the actual propeller to obtain what is traditionally referred to as open water curves. The traditional way of presenting the data tries to isolate the blades from the propeller. When the experimental data is to be compared with CFD results, it is more consistent to use the uncorrected (raw) force measured during the experiments and to compute the forces on the abovementioned surface on the measuring side of the gap.

The open-water tests were performed at 3 different propeller rotational speeds: 7, 9, and 11 rps in the range of advance coefficients from 0 to 1.2 in 0.2 steps for the pitch setting *P*/*D* = 1.1 and from 0 to 1.05 in 0.15 steps for the pitch setting 0.9. Tests were also performed at *P*/*D* = 1.2, but the results from this pitch setting have not been simulated yet.

The experimental relative uncertainity at 95% confidence for the thrust and torque coefficients has been evaluated using the methodology described in Reference [10] and found to be weakly dependent on the advance coefficient J as Figure 3 shows for P1374 at n = 7 rps. Slightly smaller values of the uncertainty were found for the other rps. The uncertainty reported here is relative only to propeller P1374, but similar values were found also for the other propellers.

**Figure 3.** Example of relative uncertainity for propeller P1374 at *P*/*D* 1.1 and *n* = 7.

#### **3. FSI Simulations**

Fluid-structure interaction (FSI), in the widest sense, is the thermomechanical interaction of one or more solid structures with an internal or surrounding fluid flow. FSI problems play prominent roles in many scientific and engineering fields, such as mechanical, aerospace, and biomedical engineering. Thanks to the recent advances of commercial software, numerical simulations have become a feasible and efficient way to investigate the fundamental physics involved in the complex interaction between fluids and solids. Nevertheless, there is still the need for validation when these tools are used in new applications.

When using a numerical procedure, an FSI problem is usually described in terms of a fluid domain and a structural domain, communicating through an in-between fluid-structure interface. This multi-physics problem with adjacent domains can be simulated in a monolithic or in a partitioned way [11]: the former simultaneously solves the governing equations for fluid and structure within a unified algorithm, and it requires a code developed for this particular combination of physical problems. In the partitioned approach, on the other hand, the fluid and the structure are treated as two computational fields which can be solved separately with their respective mesh discretizations and numerical algorithms, and the boundary conditions at the interface are used explicitly to transfer information between the fluid and structure solutions. The partitioned approach enables the use of existing well-established fluid and structural solvers. This represents a significant motivation for adopting this approach in addition to the fact that, for many problems, the staggered approach works well, is very efficient, and preserves software modularity. Difficulties may arise, however, in terms of robustness and accuracy for specific problems like flutters or parachute problems [12,13].

In the present case, we consider the problem of viscous incompressible fluid flow interacting with an elastic body (the propeller blade) immersed in the flow and being deformed by the fluid action. A staggered approach is used, coupling two different solvers that use different methods and codes: a CFD/FVM (Computational Fluid Dynamics/Finite Volume Method) solver for fluids and FEA/FEM (Finite Element Analysis/Method) for structure mechanics.

In "two way" FSI, the results of the CFD solver, in terms of hydrodynamic forces on the interface, are mapped from the fluid model to the structural model and are provided to the structural equation solver; the results obtained by the structural solver are mapped back to the first model in terms of displacements of the interfacial surface, involving the modification/morphing of the mesh of the fluid model. The procedure is repeated until convergence. How often this mapping needs to be performed depends on the degrees of coupling, and it is related to the response times of the structure compared to the fluid: in case of strong coupling, the mapping may be needed in the internal time-step iterations inside each time step, following an implicit method.

In the present case, the fluid and structure are "weakly" coupled because the response of the structure to a disturbance in the fluid is slow compared to the fluid. As seen from the experiments, a "static" solution is expected, which means that one can search for a steady-state condition, corresponding to the shape the structure takes in the stationary fluid when it reaches the hydroelastic equilibrium, corresponding to the steady-state flow around the deformed blade. All these considerations have important consequences on the flexible-blade simulations:


#### *Numerical Setup*

In the present work, STAR-CCM+ v12.04 by Siemens PLM Software [14] is used as CFD code and Abaqus 2016 by Dassault Systèmes is used as FEM code. An open-water simulation was used as a starting point, with the following characteristics: the computational domain is corresponding to one-blade passage only, the MRF (Multi-Reference Frame) approach is used. Using the MRF approach means that the propeller rotation is not simulated with the actual rotation of the mesh but that it is modeled using a reference system rotating with the blades. The equations of motion are solved in this reference system, once the boundary conditions have been correctly set up on the propeller surfaces with the proper value of velocity.

The computational domain and boundary conditions are illustrated in Figures 4 and 5. It is a 90◦ sector of a cylinder, of which the axis corresponds to the propeller axis of rotation and radius equal to 10 propeller diameters. The propeller is located in the region highlighted in red in Figure 5, where a polyhedral mesh is generated and then extruded further upstream and downstream toward the inlet and outlet boundary respectively.

**Figure 4.** Computational domain used in the simulation: The picture shows the size of the computational domain compared to the blade size. Note the fact that only one of the four blades was simulated.

**Figure 5.** Domain size and boundary conditions: The figure depicts the domain size in terms of propeller diameters and the function of the different blocks used in the meshing procedure.

The geometry of the numerical propeller model is kept as close as possible to the one used in the experiments, including the gap of 1 mm between the hub and the cone; the housing is modeled with infinite length, extended up to the outlet.

A longitudinal section of the volume mesh is depicted in Figure 6. The mesh for all the considered cases consists of about 1.6 M cells, including 10 layers of prismatic layers at the wall. The thickness of the first prismatic layer is set to approximately have *Y*<sup>+</sup> of the order of 1 on the blades. The Gamma–ReTheta transition model is used for turbulence, which performed better than the Shear Stress Turbulent (SST)-K-ω and the Reynolds Stress Turbulent (RST) models; the reason is because, at the model scale, there are likely areas of the blades subjected to laminar flow with transition to turbulent flow.

**Figure 6.** Longitudinal section of the volume mesh of the fluid dynamic solver.

The structural mesh consisted of 2.1 × <sup>10</sup><sup>5</sup> tetrahedral elements and was generated by the Abaqus embedded mesher. To determine the number of cells that was necessary for a correct structural simulation, a mesh refinement study was performed where the mesh was refined until the displacements at the blade tip were no longer changing with the mesh size. Mesh refinements were applied at the leading and trailing edges of the blade. While the mesh target size was set to 0.5 mm, the mesh refinement at the edges dominated the mesh size as it can be seen in Figure 7, where the structural mesh of propeller P1565 is shown as an example.

**Figure 7.** As an example of the mesh used by the strucural solver, the mesh relative to propeller P1565 is presented.

For the structural part of the computations, it is important to model the root of the blade and the fastening system. The blade root is built in the same resin as the rest of the blade and has several contact surfaces that constrain it to the propeller hub. The different contact surfaces determine different boundary conditions. In Figure 8, the different boundary conditions are identified with different colors; the red surfaces are an encastre-type boundary condition, while the green surface identifies a surface that allows only rotation along the direction perpendicular to the surface and translations in the plane identified by the surface.

**Figure 8.** Boundary condition at the blade root: The color codes can be found in the text.

The three propellers designs, P1374, P1566, and P1565, are analyzed. Initially the open-water calculations are carried out for the rigid case, with reference to the model tests with the aluminium propellers. The resulting flow field is then used as initial flow state for the CFD-FEM co-simulations, activating the Abaqus and the morpher solvers.

The mapping is managed by StarCCM+, where the co-simulation setup is accomplished by identifying the coupled model parts and the exported/imported fields (pressure/displacements) and by specifying the external code execution details (commands, input files, etc.).

#### **4. Results**

A converged solution of the co-simulation is reached well before the simulation time, initially fixed to 10 s, with time step of 0.005 s, as shown in Figure 9, where the axial displacements at different locations of the blade are monitored during the simulation for the propeller P1374, *P*/*D* = 1.1, and *n* = 11 rps. The results obtained for the three propellers at *P*/*D* = 1.1 and *P*/*D* = 0.9 are reported in Figures 10 and 11, respectively: *kT*, *kQ*, and *η* obtained by the co-simulation are plotted as a function of the advance ratio *J* both for the aluminum propeller and for the (flexible) resin propeller and compared to the experimental values. The agreement is very satisfactory, also in terms of the relative increase of the coefficients in the flexible case with respect to the rigid one. Except for the bollard condition and for the *J* = 1 case, where the blades can experience flow with "negative" angle of attack and subsequent separation and unsteadiness, the difference between the numerical results and the experimental data is less than 2% in terms of thrust coefficients, as shown for example in Figure 12, where the differences for *kT*, *kQ*, and *η*<sup>0</sup> are shown for the propeller P1374, *P*/*D* = 1.1. Note that, with reference to Figure 3, the deviation seen between CFD and experiments is, in most cases, within the uncertainty bands of

the experiments. The only case where the deviations seem to be large is for the propeller efficiency at *P*/*D* = 0.9 and *J* = 1.0; however, it should be noted that the *KT* coefficient is slightly negative, leading to a negative efficiency, which is strictly speaking not consistent with the definition of propeller efficiency, and hence, the propeller efficiency for that specific condition should be disregarded; further, the thrust and torque coefficients are correctly captured by the simulations.

**Figure 9.** Propeller P1374, *P*/*D* = 1.1, and *n* = 11 rps: axial displacements at different locations of the blade monitored during the simulation.

**Figure 10.** *Cont*.

**Figure 10.** Results for *P*/*D* = 1.1: The results are presented in a matrix form where the different geometries are arranged in rows and the rps is in columns. The continuous lines are relative to the experimental results, while the dots are relative to the simulations. Further, the dashed lines refer to the aluminum models and the solid lines refer to the resin models; the hollow dots are relative to the computation with the rigid propeller, whereas the the solid dots are relative to the flexible propeller.

**Figure 11.** *Cont*.

**Figure 11.** Results of *P*/*D* = 0.9: The results are presented in a matrix form where the different geometries are arranged in rows and the rps is in columns. The continuous lines are relative to the experimental results, while the dots are relative to the simulations. Further, the dashed lines refer to the aluminum models and the solid lines refer to the resin models; the hollow dots are relative to the computation with the rigid propeller, whereas the the solid dots are relative to the flexible propeller.

**Figure 12.** *Cont*.

**Figure 12.** Propeller P1374, *P*/*D* = 1.1, *n* = 7, and *n* = 11 rps: relative differences between numerical and experimental data in terms of thrust and torque coefficients and propeller efficiency.

#### *Blade Displacements*

It is rather straightforward to extract the blade cylindrical section displacements from the numerical computations and to compare the effect of the different skew distributions of the three blades. One way that seems quite natural is to decompose the displacements of the sections in terms of bending and twisting. In this approximation, the bend is defined as the average displacement of the blade sections in the YZ plane, where Z is the propeller advance direction, X is the radial direction, and Y is the direction perpendicular to both X and Z. Twist is defined as the rotation of the blade section in the same YZ plane. Some results for this twist definition are shown in Figures 13 and 14. Close to the tip, the twist component is correctly identified, while erroneous results appear closer to the blade root. Still, this presentation will be used here since it offers a compact way of presenting the data. Note also that the reservations on the bend-twist presentation of the data apply only to twist.

It was expected that the elastic behavior of the low aspect ratio propeller blades could only partially be represented by pure bend and twist of the blade sections. It is important to keep in mind this finding when trying to implement design codes that are based on lower order theories, as for example lifting line theory coupled with beam theory. It may be beneficial to consider also the warping in addition to the bending and twisting of the blade sections, as shown in Figure 14.

Even with the reservation that has just been described, it is possible to use the bend-twist approximation to explain the effect of the propeller skew on propeller blades that are otherwise identical. In Figures 15 and 16, the bend and twist of the sections from the radial location 0.4 to 0.975 are plotted.

**Figure 13.** Bend and twist of a typical blade section towards the tip: The top plot shows the same cylindrical section unloaded and loaded. The bottom picture shows the same sections once the bending is removed from the loaded blade so that the blade twist can be more easily seen.

**Figure 14.** Total displacement and twist of a section close to the blade root: The top plot shows the same cylindrical section unloaded and loaded. In the bottom picture, the procedure of removing the section bending is applied showing that, close to the blade root, the sections are warping rather than twisting.

Although it may appear to an extent fictitious, since the skew distribution is often used to control cavitation and little room is left for other uses of it, in a hypothetical case where cavitation is not a concern, it may be in principle possible to think using the skew distribution to achieve an hydroelastic response of some sort. The mechanism by which the bend and twist are coupled in isotropic materials is the so-called geometrical bend twist coupling. In the geometrical bend twist coupling, the skew distribution influences the relative position of the pressure center and the elastic center leading to different elastic behaviors even for the same, for the sake of simplicity, uniform pressure distribution. In reality, the skew distribution influences also the pressure distribution on the blade, but in a first approximation ,we assume the geometrical bend-twist coupling as the main driving factor for the different behavior of the three geometries. In Figures 15 and 16, the bend and twist distribution along the blade radius for the three propellers at two operating conditions are shown. The operating conditions are chosen to be close to the design point and an off-design condition where the propeller is heavy-loaded. Not surprisingly, the different skew distributions of the three geometries have limited effects on the propeller bending, while the effect is clearly seen on the twist of blade sections. The blade section twist is defined as positive if, as a result of the deformation, the blade section ends up having a larger pitch than in the unloaded condition. At a design point, the balanced skew (P1374) and the zero skew (P1565) propellers are subjected to almost no twist while the unbalanced skew propeller (P1566) has a slight negative twist. In the off-design condition, the skew distribution affects to a larger extent the elastic behavior of the propellers; in this condition, all blades show a twist distribution that is increasing from the root towards the tip of the blade. The zero skew blade shows an almost linear, comparatively large twist of the blade sections. The two blade designs with skew show a somewhat similar behavior with the section towards the tip, bending considerably more than those at the root. Even though, as it has already been pointed out, the bend-twist decomposition of the blade deformation may not be accurate especially for the sections close to the blade root, a few conclusion can be drawn. The first conclusion is that it is indeed possible to use the skew distribution to obtain an hydroelastic bend-twist coupling for isotropic materials. The second conclusion, at least for the geometries that have been tested, it is that the bend-twist coupling that is obtained is hardly useful for achieving the goals that are typically expected from self-adaptive blades; in fact, all the skew distributions tend to show a tendency to increasing the twist angle and, hence, the pitch, with increasing loading, a result that is the opposite of what is generally desirable. These result are far from being surprising; however, it may be useful in some special applications to consider geometric bend-twist coupling as well when designing isotropic blades that are expected to work in off-design conditions.

**Figure 15.** Bending of the three different blade geometries at *P*/*D* = 1.1 and 11 rps.

**Figure 16.** Twist of the three different blade geometries at *P*/*D* = 1.1 and 11 rps.

#### **5. Discussion**

In this paper, the validation of fluid-structure interaction computations for flexible propeller blades carried out using the co-simulation approach has been described. The validation data come from a series of ad hoc tests performed with the specific aim of serving as a reference for numerical computations. For the experiments, three propellers were manufactured both in a metallic and plastic material in order to have rigid and flexible models. The propeller geometries and the experimental data are public domain and can be obtained upon request. The simulated values of thrust and torque coefficients lie mostly within the bands of experimental uncertainty, showing how CFD-FEM co-simulations can be used to accurately represent the behavior of homogeneous isotropic blades tested in open water. Furthermore, the thrust coefficients seem to compare better than the torque coefficients. It is often reported that, in model-scale CFD simulation, the torque coefficient suffers from deviation that are larger than those observed on the thrust coefficient, a phenomenon that is often attributed to the laminar-to-turbulent transition of the flow on the model-scale blades. Since the Reynolds numbers of full-scale propellers are higher, it is to be expected that CFD can better predict the full-scale torque coefficient, since the flow is mainly turbulent, a fact that makes the approach described here even more valuable. Finally, in the paper, it is shown that warping of the blade sections towards the root of the blade may be a better representation of the actual deformation of the blade than twist.

**Author Contributions:** Conceptualization, L.S. (Luca Savio) and L.S. (Lucia Sileo); methodology, L.S. (Lucia Sileo); software, L.S. (Lucia Sileo) and S.K.Å.; validation, L.S. (Luca Savio); resources, L.S. (Luca Savio); writing—review and editing, L.S. (Lucia Sileo), L.S. (Luca Savio), and S.K.Å.; project administration, L.S. (Luca Savio); funding acquisition, L.S. (Luca Savio). All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Research Council of Norway and Kongsberg Maritime grant number 267495/O80.

**Acknowledgments:** The present work has been fully supported by the FleksProp project. The FleksProp project is a cooperation between SINTEF Ocean, the Norwegian University of Science and Technology NTNU, and Kongsberg Maritime with the economic support of the Research Council of Norway (RCN) and Kongsberg Maritime. The economic support of the Research Council of Norway and Kongsberg Maritime is greatly appreciated.

**Conflicts of Interest:** The authors declare no conflict of interest.

*J. Mar. Sci. Eng.* **2020**, *8*, 21

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **Appendix A. Propeller Geometry**

*Appendix A.1. P1374*

Main propeller elements


#### **Table A1.** Main geometrical parameters of propeller P1374.


#### *Appendix A.2. P1565*

Main propeller elements



**Table A2.** Main geometrical parameters of propeller P1565.

#### *Appendix A.3. P1566*

Main propeller elements




*Appendix A.4. Blade Section Profiles for P1374/P1565/P1566*

Blade section profiles

x/b—chordwise coordinate, [0;1]; Yu,Yl—ordinates of the upper/lower sides of the profile, Yc—ordinates of the profile mean line, given in mm

**Table A4.** Geometry of the blade sections—P1374/P1565/P1566.




#### **Appendix B. AKPD/AKPA Geometrical Definitions and Sign Conventions**

**Figure A2.** AKPD-AKPA geometrical definitions.

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
