\chapter{8m diameter rotor design}\label{Chapter5}

\noindent This chapter presents a design methodology and analysis for a turbine for which the flow to the blade is fully turbulent. This rotor's development allows for the computationally efficient simulation of tidal rotors as fully turbulent conditions permit the $k-\omega$ SST turbulence model with enhanced wall treatment to be used. Smaller turbines, such as the one used in Chapter \ref{Chapter3}, require laminar-turbulent transitional turbulence closure models to determine the boundary-layer fluid mechanics. This is impractical when analysing turbine fences due to the computational requirements of such models. The 8 m diameter turbines developed in this Chapter allows for detailed blade resolved simulations of turbine fences (Chapter \ref{Chapter6}), improving the understanding of fence thrust and power distribution and mitigation options to reduce turbine fence end-losses. The tidal turbine was designed to operate in Strangford Lough, Northern Ireland.

% \noindent As observed in Section \ref{Chapter3}, the presence of a lateral constraint on a turbine fence can enhance the outboard rotors performance due to fluid, which would otherwise be diverted around the fence, realigning normal to the rotor. Section \ref{Chapter3} examined the use of a turbine fence end wall to provide this constraint. A turbine fence endwall may be impractical for larger rotors due to the increased amount of raw material and buoyancy forces associated with the endwall. An alternative for larger turbines may be to yaw the outboard rotors to align with the skewed flow around the turbine fence. An 8 m diameter rotor was designed for use in a turbine fence to examine this effect. This chapter presents a design methodology and analysis of a single turbine designed to operate in Strangford Lough, Northern Ireland. This turbine will also be used in Section \ref{Chapter6} to investigate yaw of outboard turbines in a 4 turbine fence.

% may need to justify the 8m diameter - Impractical for sidewall and possible for scale testing. Higher $Re$ makes for easier CFD

\section{Design process}

A turbine was designed using the process presented in Figure \ref{fig5:Design_method}. There were four phases of the design methodology, the first was to define the design parameters, and the following three were numerical simulations of the rotor with increasing levels of refinement and model complexity. The subsequent sections highlight the details of each design phase.

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.65\textwidth]{Figures/Chapter_5/Design_process.png}
    \caption{Turbine design methodology.}
    \label{fig5:Design_method}
\end{figure}

\section{Design definition}

Flow conditions used for the rotor design were derived based on the hydrodynamic conditions at Strangford Lough in Northern Ireland. This site provides a sheltered environment from wave action owing to the narrow channel leading into the Lough (Figure \ref{fig6:Strangford_lough}) and has previously been the test location of the SeaGen turbine \citep{MacEnri2011}. The channel's width is approximately 550 m at the test location and has a depth of 30 m according to a multibeam bathymetric survey of the Lough from 2014 \citep{AFBI2015}. The area surrounding the site has a relatively flat terrain (between 0-1.9$^\circ$).

Tidal currents through Strangford Lough's narrows can reach up to 3.5 m/s during spring tides with more typical peak velocities during flood and ebb tides of 2.5 m/s \citep{Kregting2014}. A survey by \cite{MacEnri2013} indicated that the turbulence intensity in the channel ranges between 3\% and 9\%, with the values being inversely proportional to flow speed. The hydrodynamic conditions selected for the rotor design are summarised in Table \ref{tab7:hydro_design}.

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.8\textwidth]{Figures/Chapter_6/strangford_lough_map.jpg}
    \caption{Strangford Lough showing the SeaGen testing location within the Narrows.}
    \label{fig6:Strangford_lough}
\end{figure}

% Do I need to credit Google there?

\begin{table}[]
\centering
\caption{Summary of hydrodynamic design conditions}
\label{tab7:hydro_design}
\begin{tabular}{l|l|l}
\textbf{Parameter}                                                         & \textbf{Value} & \textbf{Notes}                                                                                                                                   \\ \hline
Design velocity                                                            & 2.5 m/s          & \begin{tabular}[c]{@{}l@{}}Based on the study of \\ \cite{Kregting2014} this value is \\ reached on 55\% of tidal cycles\end{tabular} \\ \hline
\begin{tabular}[c]{@{}l@{}}Freestream \\ turbulence intensity\end{tabular} & 6\%            & \begin{tabular}[c]{@{}l@{}}This is the average measurements from\\ \cite{MacEnri2013} for the design \\ velocity\end{tabular}   \\ \hline
Width                                                                      & 550 m          & \begin{tabular}[c]{@{}l@{}}Survey data from \\ \cite{AFBI2015}  \end{tabular} \\ \hline
Depth                                                                      & 30 m           & \begin{tabular}[c]{@{}l@{}}Survey data from \\ \cite{AFBI2015}  \end{tabular}                                                                                                
\end{tabular}
\end{table}


Table \ref{tab7:turb_design} lists rotor design parameters. These parameters were selected based on an overview of existing turbine technology and applying scaled values to an 8 m diameter turbine. The nacelle diameter was selected based on the Atlantis AR1500 \citep{SimecAtlantas2016}, Hammerfest Strøm \citep{AndritzHydro2016}, and the experimental rotor of \cite{Cao2019} which all had a nacelle to rotor diameter ratio between 10\% to 15\%. A design range for $\lambda$ was selected based on the AR1500 and the Hammerfest Strøm turbines, designed for $\lambda$ between 4.15 and 5.75. This range is lower than in the smaller rotor of \cite{Cao2019} which was designed for $\lambda=7$; this is partly due to the larger diameter rotors needing to limit tip speed to avoid cavitation inception. Finally, an operating depth placing the nacelle at 1/3 depth was selected based on comparison with the Orbital Marine O2 device \citep{OrbitalMarinePower2019}. Positioning the turbine at this depth ensures it is located away from the seafloor in the more energetic flow region of the water column and provides additional cavitation protection due to the static pressure head. Additionally, this depth provides some tip clearance if vessels are to pass over the fence. 


\begin{table}[]
\centering
\caption{Summary of turbine design conditions}
\label{tab7:turb_design}
\begin{tabular}{l|r}
\textbf{Parameter}        & \multicolumn{1}{l}{\textbf{Value}}                                                 \\ \hline
Diameter (D)              & 8 m                                                                                \\ \hline
Inter-turbine spacing (s) & 2 m                                                                                \\ \hline
Operating depth of nacelle ($Z_{hub}$) & 10 m                                                                                \\ \hline
Number of blades          & 3                                                                                  \\ \hline
Nacelle diameter          & 1 m                                                                                \\ \hline
Aerofoil                   & \begin{tabular}[c]{@{}r@{}}NACA-63-415 with \\ modified trailing edge\end{tabular} \\ \hline
$\lambda$                 & 4-6                                                                                \\ \hline
$\Omega$                  & 2.5-3.75 rad/s                                                                      \\ \hline
$Re_{chord}$ at 80\% span & $2 \times 10^6$ - $3.2 \times 10^6$                                            \\ \hline
$B_g$                       & 1.2\%                                                                              \\ \hline
$B_l$                       & 16.8\%                                                                            
\end{tabular}
\end{table}

\section{RANS-BE actuator disk design}\label{sec5:RANS_DES}

The RANS-BE actuator disk (AD) design tool developed by \cite{McIntosh2011} was used to determine the twist ($\beta$) and solidity ($\sigma$) distribution along the blade length. This tool applied the RANS-BE formulation outlined in Section \ref{CH2:BE_model} with an additional iterative loop that adjusted the blade twist and solidity to achieve a target $\alpha$ and target local thrust coefficient ($c_x$). The design tool did not consider discrete blade effects, which prevented the use of a tip correction model. Consequently, the initial rotor designs overestimated the performance; however, once the initial design was confirmed, addition of simulations were run with a tip correction model to improve the performance estimates.

% Owing to a limitation of this method to include the tip correction model in the design tool, the initial rotor design did not consider discrete blade effects and was run without a tip correction model. Instead, a design was first generated and a second run was carried out with the tip correction model activated to predict performance.

% Work by \cite{Cao2019} demonstrated that the target $c_x$ had limited bearing on the optimal design and was set to 2 over the majority of the blade. A reduction in the target $c_x$ was applied at the blade root to reduce the chord length and avoid overlapping geometry in this region. 

A modification of the design tool was made to allow interpolation over multiple 2D aerofoil data sets, including 2D aerofoil geometry, chord-based Reynolds number ($Re$) and turbulence intensity ($TI$). Figure \ref{fig5:typicalRETI} plots the expected variation in $Re$ and $TI$ along the length of a blade. Whilst it appears appropriate to interpolate $TI$ owing to the monotonic variation along the blade's length, the variation in $Re$ is relatively small, especially over the high power producing region of $0.4<r/R<0.95$. The advantage that $Re$ interpolation provides is that no new aerofoil data has to be generated when adjusting the target thrust coefficient ($c_x$) or $\lambda$, which change the reference length and velocity in the $Re$ calculations, respectively. 

Furthermore, an alteration of the RANS-BE design code, similar to \cite{Cao2019}, was included to allow interpolation between 2D profiles of different trailing edge thickness to chord ratio. This approach allowed varied aerofoil sections to be simulated without the explicit 2D modelling of the geometry for the constant trailing edge thickness case examined below.

\begin{figure}
    \centering
    \includegraphics[width=0.49\textwidth]{Figures/Chapter_5/Typical_Re.png}
    \includegraphics[width=0.49\textwidth]{Figures/Chapter_5/Typical_Ti.png}
    \caption{Typical Reynolds number (left) and turbulence intensity variation (right) with radial position on a rotor blade of 8 m diameter at $\lambda=6$. The radius of the hub is 0.125R}
    \label{fig5:typicalRETI}
\end{figure}

\subsection{2D aerofoil simulations}

To provide input into the RANS-BE design tool, 2D aerofoil simulations were run for a range of $TI$ and $Re$. The NACA-63-415 profile was selected for this rotor due to its high lift to drag ratio at high $Re$ and the availability of experimental data to validate this geometry \citep{Bak2000}. This profile belongs to the NACA-63-4XX subfamily of aerofoils which are commonly used in wind turbines \citep{Bertagnolio2001} and have also been used to design tidal rotors \citep{Laurens2016}.


\subsubsection{Trailing edge treatment} \label{CH5:TE_treatment}

An analysis varying trailing edge thickness ($\delta$) and thickness to chord ratio ($\delta/c$) of the NACA-63-415 profile examined the trade-off in performance between:
\begin{enumerate}
    \item Sharp trailing edge ($\delta=0$) - This edge treatment is the unmodified geometry of the NACA-63-415 profile. While this profile is likely to be the best performing due to the smooth edge transition, the sharp geometry is difficult to manufacture and may introduce a hazard when handling the blades. 
    \item Truncated trailing edge at constant trailing edge thickness to chord ratio ($\delta/c=constant$) - By applying the thickening function of \cite{Xu2014}, the trailing edge thickness can be adjusted to a suitable size for manufacture and handling. Using the same section along the blade's length means that trailing edge thickness is dictated by the smallest chord length of the turbine (often the tip), resulting in an efficiency loss over the other portions of the blade due to unnecessarily thick trailing edges.
    \item Constant trailing edge thickness ($\delta=constant$) - Applying the thickening function of \cite{Xu2014} at every section along the blade's length augments the aerofoil profile to maintain a constant thickness trailing edge. This reduces the losses due to an unnecessarily thick trailing edge on sections of the blade with large chord if $\delta/c$ is constant, but introduces additional complexity as there is a unique aerofoil profile at each section along the blade's length.
\end{enumerate}

The final two trailing edge treatments in the above involve an iterative approach using the turbine design tool to generate a geometry before applying the desired trailing edge thickness. A minimum trailing edge thickness of 2 mm was chosen for the blade design as the thickness provides durability for manual handling without being too sharp. Minimum and maximum chord lengths of 250 mm ($\delta/c=0.0080$) and 750 mm ($\delta/c=0.0027$) were estimated based on a similar turbine design allowing the range of potential $\delta/c$ to be defined. Applying the thickening function of \cite{Xu2014} the profiles were thickened following 

\begin{equation}
    \begin{array}{cc}
        x=x_{0}, y=y_{0} & 0 \leq x \leq x_{t} \\
        x=x_{0}, y=y_{0} \pm 0.5 \delta\left(\frac{x-x_{t}}{c-x_{t}}\right)^{n} & x_{t} \leq x \leq c,
    \end{array}
\end{equation}

\noindent where $x_0$ and $y_0$ are the original geometry, $x$ and $y$ are the thickened geometry, $x_t$ is the location of maximum thickness along the chord ($c$), $\delta$ is the thickness of the trailing edge and $n$ is a power exponent which is used to control the smoothness of the profile, taken to be 2 in this study. The profiles examined in this study are plotted in Figure \ref{fig5:Airfoil_profiles}.


\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/airfoil_profiles.png}
    \caption{NACA 63-415 aerofoil with three different trailing edge thickness to chord ratio thickened following \cite{Xu2014}}
    \label{fig5:Airfoil_profiles}
\end{figure}

\subsubsection{2D aerofoil analysis}\label{ch5:2d_aero}

ANSYS Fluent 15.0 \citep{ANSYSInc.2013} was used to perform a 2D analysis of the aerofoil sections. Lift, drag and pressure distribution data for four chord-based Reynolds numbers ($Re =$ [\num{8e5}, \num{1.6e6}, \num{3.2e6}, \num{5e6}]) and six turbulence intensities ($TI = [0.98\%,1.92\%,3.67\%,5.27\%,6.77\%,8.15\%]$) were calculated for each of the trailing edge treatments. 

Calibration results between the 2D numerical model and wind tunnel testing by \cite{Bak2000} are presented in Figure \ref{fig5:cali_exp} for the NACA 63-415 profile without any trailing edge treatment. The experiments were conducted at $Re=$\num{1.6e6} and $TI=1\%$, which are slightly lower turbulence conditions than expected at the turbine. A favourable comparison observed with regards to $c_l$ and $C_{pressure,min}$, capturing the trend and magnitude of these quantities around the optimal angle of attack. There is a slight discrepancy in the $c_d$ at low angles of attack, given the relatively low values of $c_d$, this error is compounded in the lift to drag ratio. Nevertheless these errors are within an acceptable range for use in the RANS-BE model. Lift and drag parameters are known to change dramatically during stall, therefore the calibration and analysis in the RANS-BE model was limited to conditions where $\alpha$ was less than 10$^\circ$.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/2d_cali.png}
    \caption{Calibration of $c_l$ (top left), $c_d$ (top right), lift to drag ratio (bottom left) and $C_{pressure,min}$ (bottom right) for a NACA 63-415 section at $Re=$\num{1.6e6} and $TI=1\%$, experimental data was sourced from \cite{Bak2000}}
    \label{fig5:cali_exp}
\end{figure}

Results for $Re=$\num{3.2e6} are presented in Figure \ref{fig5:Re3200000_polars}; this value corresponds to the approximate chord based $Re$ at 80\% span and $\lambda=6$ with a free-stream velocity of 2.5 m/s. The turbulence intensity has a significant impact on the lift to drag ratio (L/D) which the section can achieve. The primary cause for this is the increased skin friction with higher $TI$, which increases $c_d$. Secondary to this, an increase in $c_l$ is observed at lower turbulence intensity; however, as a proportion, this increase is significantly smaller ($\approx$20 times) than the decrease in $c_d$. 

The simulations also demonstrate that increasing the trailing edge thickness form $\delta/c=0.008$ to $\delta/c=0.00267$ results in a lower L/D at the optimal angle of attack, $\alpha^*$. Like the turbulence intensity, this is due to increased drag (increased drag was observed on thickened trailing edges at $0^\circ<\alpha<10^\circ$). The thickened trailing edges had a higher $\alpha^*$ which shifted from $\alpha^*=5^\circ$ for the sharp edge to $\alpha^*=6^\circ$. This results from a slight reduction in the gradient of the drag coefficient with angle of attack as thickness increases. The cause of this appears to be due to the thickened sections delaying separation on the upper surface of the aerofoil as demonstrated by the zero intercept of the wall shear stress ($\tau_{wall}$) in Figure \ref{fig5:wall_shear_stress}.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/re3200000_ldcpl2d.png}
    \caption{NACA 63-415 $c_l$ (top left), $c_d$ (top right), lift to drag ratio (bottom left) and $C_{pressure,min}$ (bottom right) variation with $\alpha$ for different trailing edge thicknesses and turbulence intensities at $Re=$\num{3.2e6}.}
    \label{fig5:Re3200000_polars}
\end{figure}



\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/Wall_shear_stress_2d_upper.png}
    \caption{Wall shear stress on the suction surface with chord position for a NACA-63-415 profile with varied trailing edge thicknesses at $Re=$\num{3.2e6}, $\alpha=[0,10,15]^\circ$ and a turbulence intensity of 0.98\%}
    \label{fig5:wall_shear_stress}
\end{figure}

The minimum pressure coefficient ($C_{pressure,min}$) appears to be relatively insensitive to the trailing edge treatment and turbulence intensity at $\alpha <5^\circ$. Beyond this value, the magnitude of $C_{pressure,min}$ rapidly increases and $TI$ has greater influence; these results match the wind tunnel testing of \cite{Bak2000}. Low $C_{pressure,min}$ indicates a lower pressure over the section, if this pressure drops below the vapour pressure of the liquid small vapour cavities may start to form (cavitation). As these cavities move into high pressure regions they collapse, generating a shock waves which can cause damage in their immediate vicinity. Measures to protect the tidal turbine from cavitation are discussed in Section \ref{TSR_select}.

$Re$ also impacts aerofoil performance, as demonstrated by Figure \ref{fig5:TI2_polars}, which plots the polars for a fixed $TI$ of 1.92\%, corresponding to the blade's expected turbulence conditions. Increased $Re$ leads to a higher lift to drag ratio and magnitude of $C_{pressure,min}$.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/ti2_ldcpl2d.png}
    \caption{Aerofoil $c_l$ (top left), $c_d$ (top right), lift to drag ratio (bottom left) and $C_{pressure,min}$ (bottom right) variation with $\alpha$ for different trailing edge thicknesses and Reynolds numbers at a turbulence intensity of 1.92\%.}
    \label{fig5:TI2_polars}
\end{figure}


\subsection{Model domain}

A cylindrical domain of diameter 72.47 m (9.06$D$) was used for the initial rotor design to match the global blockage of 1.2\% at the Strangford Lough location. Figure \ref{fig5:Turbine_design_model_domain}  plots this domain and highlights the velocity inlet, pressure outlet and symmetry boundary conditions applied in the model. The velocity inlet is 10 $D$ upstream of the turbine with a specified inflow of 2.5 m/s with a turbulence intensity of 8\%, which decays to 6\% at the rotor plane. A pressure outlet boundary condition was applied 20D downstream of the rotor plane.


\begin{figure}[h]
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/Turbine_domain.png}
    \caption{Model domain for RANS-BE turbine design.}
    \label{fig5:Turbine_design_model_domain}
\end{figure}



\subsection{Mesh convergence study}

A mesh convergence study analysed three meshes of varied minimum cell sizing at the rotor's root and tip. A constant mesh growth rate of 1.05 was used in the area immediately surrounding the rotor and nacelle; this was relaxed towards the inlet and outlet, where a value of 1.1 was applied. Mesh convergence was assessed based on turbine $C_p$ and $C_t$. The results are summarised in Table \ref{tab5:MCS_turb}. The difference in $C_p$ and $C_t$ between the coarse and fine mesh was 1.5\% and 1\% respectively. The primary cause for this discrepancy is the mesh's resolution in the root region, where the influence of the nacelle and decreased angular velocity causes a rapid change in the required blade twist and solidity. This effect is highlighted in Figure \ref{fig5:MCS_Beta_sigma} with the steep gradient in the near root region; as this location provides little of the overall power output, the medium resolution mesh was adopted. 


\begin{table}[]
\centering
\caption{Turbine design tool performance across mesh resolutions.}
\label{tab5:MCS_turb}
\begin{tabular}{|l|r|r|r|r|}
\hline
\multicolumn{1}{|c|}{Mesh name} & \multicolumn{1}{c|}{\begin{tabular}[c]{@{}c@{}}Minimum element\\ size ($L_{element}/D$)\end{tabular}} & \multicolumn{1}{c|}{\begin{tabular}[c]{@{}c@{}}Number of \\ elements \\ (Millions)\end{tabular}} & \multicolumn{1}{c|}{$C_p$} & \multicolumn{1}{c|}{$C_t$} \\ \hline
Coarse                          & 0.00500                                                                                               & 2.7                                                                                              & 0.574                      & 0.908                      \\ \hline
Medium                          & 0.00250                                                                                               & 7.3                                                                                              & 0.576                      & 0.911                      \\ \hline
Fine                            & 0.00125                                                                                               & 20.0                                                                                             & 0.582                      & 0.917                      \\ \hline
\end{tabular}
\end{table}

\begin{figure}
    \centering
    \includegraphics[width=0.49\textwidth]{Figures/Chapter_5/MCS_beta.png}
    \includegraphics[width=0.49\textwidth]{Figures/Chapter_5/MCS_sigma.png}
    \caption{Blade twist (left) and solidity (right) distribution from the RANS-BE design tool for varied mesh refinement.}
    \label{fig5:MCS_Beta_sigma}
\end{figure}


\subsection{Tip speed ratio selection}\label{TSR_select}

Determining the appropriate $\lambda$ for design requires balancing rotational speed for optimal performance against the onset of rotor blade cavitation and turbine thrust as it provides more resistance to the flow. As observed in the 2D aerofoil plots in Figure \ref{fig5:TI2_polars}, lower turbulence intensity and higher Reynolds number result in an improved lift to drag ratio for the section. From a design perspective, this would suggest that a higher $\lambda$ should be targeted (if flow diversion effects are ignored due to the higher thrust) as the higher velocity increases the $Re$ and lowers turbulence intensity. However, constraining this problem is the onset of blade cavitation, which can degrade blade surfaces and reduce performance. It is recommended that cavitation be avoided for tidal turbines \citep{DNVGL2015}.

To satisfy the no cavitation inception criteria when defining $\lambda$, an approach modified from \cite{Wimshurst2018} was applied to determine the maximum permissible tip speed ratio ($\lambda_{max}$) at a specified depth. The no cavitation requirement constrains the cavitation number ($C_a$) always to be greater than 0. $C_a$ is defined as 
\begin{equation}
    C_a = \frac{p - p_{v}}{\tfrac{1}{2}\rho W^2}
\end{equation}

\noindent where $p_v$ is the vapour pressure, taken to be 2.4 kPa for seawater, $W$ is the relative speed of the flow to turbine blade and $p$ is the static pressure, inclusive of static pressure over the 2D section ($p_{section}$), atmospheric pressure ($p_{atm}$) and hydrostatic pressure ($p_z$) at depth $z$ such that $p = p_{section} + p_{atm} + p_{z}$.

Using the $C_{pressure,min}$ coefficient for the 2D section, an approximate minimum static pressure over the section ($p_{section}$) can be estimated following

\begin{equation}
    p_{section} = C_{pressure,min} \tfrac{1}{2} \rho W^2.
\end{equation}


The critical location for cavitation is at the minimum $p$ along the blade; this was taken to be the blade's tip when positioned at top-dead-centre (a depth of $Z_{hub} - D/2$) as the hydrostatic pressure is at a minimum. Blade resolved simulations of \cite{Wimshurst2018} demonstrated that this is a conservative assumption as minimum $p$ occurs slightly closer to the nacelle (94\% span) due to tip loss effects reducing the suction peak. This effect was ignored for the initial estimates of $\lambda_{max}$ and a conservative assumption for $p_z$ was used,
\begin{equation}
    p_z = \rho g (Z_{hub} - D/2).
\end{equation}

Finally, an approximation of the velocity magnitude was made by assuming that the axial and tangential flow induction factors at the tip were 0 and 1, respectively, allowing
\begin{equation}
    W \approx \sqrt{u_{\infty}^2+(\Omega R)^2}.
\end{equation}

\noindent This assumption is reasonable given the rotational component of velocity is much greater than the axial component and is also conservative as the higher velocity results in a more significant suction peak. $\lambda_{max}$ can now be expressed as a function of $Z_{hub}$ for a given $\alpha$ following

\begin{equation}
    \lambda_{max} =  \frac{\sqrt{\frac{2 (p_v - p_{atm} - \rho g (Z_{hub} - D/2)) }{\rho C_{pressure,min}} - u_{\infty}^2}}{u_{\infty}}.
\end{equation}


Figure \ref{fig5:Ca_lambda} plots $\lambda_{max}$ for flow conditions expected at the blade tip for varying depths. It can be observed from this figure that $\alpha$ plays an essential role in defining $\lambda_{max}$ and that reducing the blade twist at the tip may permit a higher $\lambda_{max}$. Secondary to this, increasing trailing edge thickness decreases the $\lambda_{max}$ of the rotor; this is due to the higher suction peaks on the thickened aerofoils. 

\begin{figure}
    \centering
    \includegraphics[width=0.9\textwidth]{Figures/Chapter_5/Cavitation_with_depth_Re3200000_TI1.png}
    \caption{Variation in $\lambda_{max}$ with tip submersion depth for $Re=$\num{3.2e6} and turbulence intensity 0.98\% which are typical conditions expected for the D=8 m rotor design. Depending on the trailing edge thickness, $\alpha^*$ is either 5$^\circ$ ($\delta/c=0$) or 6$^\circ$ ($\delta/c>0$).}
    \label{fig5:Ca_lambda}
\end{figure}

To maximise the available energy passing the turbine, it is advantageous to select a rotor depth as close as practical to the water surface due to the flow's shear profile, this height has to be balanced against the impact of waves on the rotor. A smaller submersion depth is also beneficial for a floating turbine as the smaller lever reduces the moment on the platform. A rotor hub submersion depth of 1/3rd the water depth (10 m) was determined based on a comparison to Orbital Marine's O2 device \citep{OrbitalMarinePower2019}. This provides a tip submersion depth of 6 m which corresponds to a $\lambda_{max} \approx 4.5$ for $\alpha^*=6$.

The analysis suggests several ways to increase the $\lambda_{max}$ if desired. These include modifying the section towards the tip to have a lower suction peak or increasing the blade twist angle to target a lower angle of attack. The latter strategy is indirectly employed in this study as a Glauert tip correction factor is not used in the blade design process. The tip correction accounts for the effect of the span-wise flow and tip vortex by reducing the free-stream flow component; not including the factor will, therefore, reduce the attack angle at the rotor's tip for the design.


\subsection{Blade trailing edge}

Three rotors were designed using the modified RANS-BE design tool with the edge treatments described in Section \ref{CH5:TE_treatment} for $\lambda=4.5$, $c_x=2$ and $\alpha^*=6$. Turbine $C_p$ and $C_t$ of each design case are presented in Table \ref{tab5:trailing_edge_perf}. Despite the sharp trailing edge section being the best performing for the lift to drag ratio, maintaining $\delta=0.002$\:m over the blade results in a slightly higher $C_p$. This is primarily due to the higher lift of the thickened sections at low turbulence intensity and the slight difference in drag between the sections (discussed further in Section \ref{ch5:2d_aero}). Given this result, it may be surprising that the performance of the $\delta/c=0.008$ design case is slightly worse than the $\delta=0$ m case. Figure \ref{fig5:Lift_TI} illustrates that $c_l$ reduces towards the inboard section with increasing $TI$. However, the $\delta=0$ m is slightly less sensitive to this change allowing the section to have a higher performance in the hub region than the $\delta/c=0.008$ case. The tangential thrust coefficient ($c_{\theta}$) and the section torque coefficient ($c_q)$ are plotted in Figure \ref{fig5:des_coefficents} to highlight each radial station's contribution to the rotor performance. These coefficients are defined as

\begin{equation}
\label{eq5:ct}
    c_{\theta} =  c_l \sin(\phi) - c_d \cos(\phi)
\end{equation}

\begin{equation}
\label{eq5:cq}
    c_q = c_{\theta} \tfrac{r}{R}.
\end{equation}

Figure \ref{fig5:des_coefficents} shows that the hub region is the reason the $\delta=0$ m design has superior performance to the $\delta/c=0.008$ design. It is worth noting that both designs have a similar $c_q$ towards the tip due to the slightly higher $c_d$ for the $\delta/c=0.008$ design negating the improvement in $c_l$. The $\delta=0.002$ m case is able to maintain a higher $c_l$ towards the root due to the increased chord length reducing $\delta/c$ closer to $\delta/c=0.002$ which has a higher $c_l$ than the $\delta=0$ m dataset for the $TI$ (Section \ref{ch5:2d_aero}).

\begin{table}[h]
\centering
\caption{Rotor performance for turbines designed with different trailing edge treatment at $\lambda^*=4.5$, $c_x=2$ and $\alpha^*=6^\circ$.\\}
\label{tab5:trailing_edge_perf}
\begin{tabular}{l|l|l}
\hline
\textbf{Trailing edge treatment} & $\mathbf{C_p}$ & $\mathbf{C_t}$ \\ \hline
Sharp ($\delta=0$ m)                           & 0.575          & 0.927          \\ \hline
Varied $\delta$ ($\delta/c=0.008$)        & 0.573          & 0.927          \\ \hline
Constant $\delta$ ($\delta=0.002$ m)               & 0.580          & 0.928          \\ \hline
\end{tabular}
\end{table}


\begin{figure}
    \centering
    \includegraphics[width=0.8\textwidth]{Figures/Chapter_5/Lift_TI.png}
    \caption{$c_l$ and $TI$ variation with radial position for rotors designed with various trailing edge treatments.}
    \label{fig5:Lift_TI}
\end{figure}

Overall, the slight differences in performance between the rotors and the similar blade twist and solidity distributions, as observed in Figure \ref{fig5:angles}, imply the choice of blade trailing edge is of little consequence to the overall design. Therefore, practical design and safety (not having a sharp edge) are likely to be the limiting factor in trailing edge choice. For this reason, the $\delta=0.002$ m case was selected for the final design.

\begin{figure}
    \centering
    \includegraphics[width=0.49\textwidth]{Figures/Chapter_5/DES_ct_bem.png}
    \includegraphics[width=0.49\textwidth]{Figures/Chapter_5/DES_cq_bem.png}
    \caption{$c_{\theta}$ (left) and $c_q$ (right) variation with radial position for rotors designed with varied trailing edge treatments.}
    \label{fig5:des_coefficents}
\end{figure}




\begin{figure}
    \centering
    \includegraphics[width=0.49\textwidth]{Figures/Chapter_5/angles.png}
    \includegraphics[width=0.49\textwidth]{Figures/Chapter_5/DES_chord.png}
    \caption{Span-wise distribution of blade twist angle ($\beta$), angle of attack ($\alpha$) and their sum ($\phi$) (left) and design chord length (right)}
    \label{fig5:angles}
\end{figure}


\subsection{$c_x$ adjustments}

Optimising the target $c_x$ is an important step in the design process as it impacts overall rotor thrust and performance. Physically, increasing $c_x$ largely corresponds to increasing chord length. To develop higher values of $c_x$ rotor solidity and thrust are higher, causing a reduction in the stream-wise velocity and requiring the blade twist to be adjusted to achieve the target angle of attack. Consequentially $c_{\theta}$ (and $c_q)$ is augmented following Equation \ref{eq5:ct} and \ref{eq5:cq}; these values should be maximised over the blade through the optimisation of $c_x$. 

The optimal $c_x$, where $c_x$ is kept uniform along the blade, for $\delta=0.002$ m and $\lambda$=4.5 was determined through a two-stage iteration process. This process first determined $C_p$ for three $c_x$ before interpolating a quadratic between the data and selecting the optima for the next iteration. The process converged rapidly, avoiding unnecessary simulations. Results are plotted in Figure \ref{fig5:cx_optimisation}, which shows a maximum $C_p$ of 0.584 when $c_x=2.1$ was used.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/cx_perfromance.png}
    \caption{$C_p$ and $C_t$ variation with design $c_x$ at $\lambda$=4.5. $c_x$ is uniform along the blade span.}
    \label{fig5:cx_optimisation}
\end{figure}

While optimising a uniform $c_x$ along the blade's length is essential for improving power production, physical limits exist for the chord length. Specifically, high target $c_x$ can cause intersecting blade geometries, especially in the root region. This was observed for all $c_x$ examined in this study. Consequentially, $c_x$ was curtailed in the root region using a smoothing function to increase $c_x$ from 0.5 at the nacelle (r/R=0.125) to 2.1 at r/R=0.4; this blending function is plotted in Figure \ref{fig5:cx_adj}. A limit for blade solidity is non-trivial to include in the design code due as it requires a consideration for the twist angle and the profile thickness as blades can overlap if the twist allows. The curtailment of $c_x$ resulted in a marginally higher $C_p$ ($C_p$=0.584) than the $c_x$=2.1 case ($C_p$=0.581). This is due to the optimisation of $c_x$ uniformly along the blade and not at every radial station. When this parameter is allowed to vary a slightly higher $C_p$ is possible.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/cx_adjustment.png}
    \caption{Variation in design $c_x$ along the rotor blade.}
    \label{fig5:cx_adj}
\end{figure}

\subsection{Final geometry}

The geometry was adjusted at the blade tip and root to limit rapid changes in blade shape. At the root section, a limit on the rate of reduction in $\sigma$ was applied due to the influence of the boundary layer from the nacelle on the RANS-BE design tool. As this region is close to the nacelle, it has a relatively small influence on the turbine's power production. The adjustment made at the tip is a minor correction to limit the reduction in $\sigma$ and $\beta$ in this region. The need for this adjustment at the tip is due to the RANS-BE design tool not including a tip correction model, resulting in artificially high thrusts in this region, reducing blade $\sigma$ and $\beta$. Applying this tip correction causes a reduction in power of 8.8\%, this is the product of two factors. Firstly, not including a tip correction in the design leads to sub-optimal design parameters over this region. Secondly, artificially high thrusts exist at the tip when the correction is not applied as no load reduction is accounted for, resulting in an overestimate of power.

The geometry corrections were made by extrapolating $\beta$ and $\sigma$ using the gradient of the profile at r/R=0.175 and r/R=0.95 for the root and the tip respectively (producing a smooth geometry to the first derivative). The resulting geometry profiles are presented in Figure \ref{fig5:geom_adjust}, and the impact of this adjustment on the power and thrust can be observed in Table \ref{tab5:design_mod}. The table shows a change in $C_t$ of 0.1\% when the geometry adjustment is applied to the design.


\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/Interpolation_Check.png}
    \caption{Original (orange) and extrapolated (blue) turbine geometry; extrapolation has used the gradient at the outer 5\% of the blade to define the geometry beyond these extents.}
    \label{fig5:geom_adjust}
\end{figure}


\begin{table}[]
\centering
\caption{Impacts of design modifications on performance.}
\label{tab5:design_mod}
\begin{tabular}{|l|l|l|}
\hline
\textbf{Turbine design condition}                   & \textbf{$C_p$} & \textbf{$C_t$} \\ \hline
$c_x=2.1$                                          & 0.581          & 0.945          \\ \hline
Curtailed $c_x$                                    & 0.584          & 0.935          \\ \hline
Curtailed $c_x$ + tip correction                   & 0.530          & 0.888          \\ \hline
Curtailed $c_x$ + tip correction + geometry adjustment & 0.530          & 0.889          \\ \hline
\end{tabular}
\end{table}


\section{Blade resolved MRF simulation}

MRF (Multiple Reference Frames) simulations were run to improve the estimates for power and thrust, investigate blade tip treatment and develop the meshing strategy for sliding mesh simulations. The MRF method modifies the Navier-Stokes equations within a sub-domain surrounding the blade and solves them in a moving reference frame. This approach permits steady solution methods, provided that the inflow conditions are steady, allowing for a computationally efficient simulation that resolves the blade geometry without needing a sliding mesh. Further details on the MRF method can be found in Chapter \ref{Chapter2}.

\subsection{Model domain}

The model domain was identical in size to the RANS-BE simulations except for the flow cross-section, which was reduced to a 120$^\circ$ sector by exploiting the rotor's azimuthal periodicity. The additional faces introduced to the model from splitting the domain were replaced with periodic boundaries to replicate the complete cylindrical domain. The model was divided into two sub-domains, the first being the outer domain and the second being an inner domain. The inner domain used a rotating reference frame to simulate the rotation of the turbine. Figure \ref{fig5:MRF_domain} highlights the two domains and the boundary conditions used in this model.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/MRF_Domain.png}
    \caption{MRF Domain highlighting the inner and outer domains and boundary conditions.}
    \label{fig5:MRF_domain}
\end{figure}

\subsection{Mesh convergence}

Mesh convergence was assessed using structured and semi-structured meshing approaches. Structured meshes offer control to users, requiring them to specify node distribution on the edges within the domain, as a consequence this method requires a significant amount of user input, especially for complex geometries. In contrast semi-structured methods take input parameters (such as a global mesh growth rate and maximum element size) to generate a mesh from the geometry. Consequentially, semi-structured meshes are faster to create but may lack the resolution around the area of interest which can be achieved using structured methods. To address this short coming, improvements to semi-structured meshing software have allowed user specification of node distribution and mesh growth rates normal to surfaces using inflation layers. This leads to a mesh very similar to a structured mesh close to surfaces and a conventional unstructured mesh in the remaining domain; this approach is used for the semi-structured mesh presented in this section. 


% have been used almost exclusively by the Oxford Tidal Research Group to simulate tidal rotors \citep{Cooke2016,Vogel2017,Cao2018}. However, improvements in unstructured mesh controls and solvers call to question the time-intensive structured meshing approach and if satisfactory convergence can be achieved using an unstructured methodology.

\subsubsection{Structured mesh convergence}

A block-structured hexahedral mesh was developed using ANSYS ICEM \citep{Inc.2012}. The mesh parameterisation was similar to that of \cite{Cao2018}; however, the near wall resolution was modified due to the higher turbulence conditions permitting the use of the $k-\omega$ SST turbulence model. This turbulence model used an enhanced wall treatment which models the boundary layer and requires either resolution in the viscous sublayer ($y^+<5$) or in the outer log-law region ($y^+>30$). Model resolution $5<y^+<30$ places cells in the buffer region where this model has been shown to be less accurate \citep{ANSYSInc.2013}. 

Table \ref{TAB5:StrucConverge} summarises the results from a mesh convergence study for the blade area-weighted-average $y^+$ ($\overline{y^+}$) and near-wall growth rate ($GR_{NearWall}$). A $GR_{NearWall}=1$ indicates a buffer layer with elements of uniform thicknesses normal to the surface. Overall the $C_p$ and $C_t$ for the rotor was most sensitive to the $\overline{y^+}$ range tested with a $\overline{y^+}<66$ having $<$1\% error; all $GR_{NearWall}$ tested also resulted in $<$1\% error in $C_p$ and $C_t$.

\begin{table}[h]
\centering
\begin{threeparttable}

\caption{Structured mesh convergence checks}
\label{TAB5:StrucConverge}
\begin{tabular}{l|rr|rr}
\textbf{Variable} & \textbf{$\overline{y^+}$} & \textbf{$GR_{Near Wall}$} & \textbf{$\frac{C_{t}^{c}-C_{t}}{C_{t}^{c}}$} & \textbf{$\frac{C_{p}^{c}-C_{p}}{C_{p}^{c}}$} \\ \hline
\multirow{4}{*}{\textbf{$\overline{y^+}$}} & \textbf{217.6} & 1 & -2.46\% & -3.10\% \\
 & \textbf{107.8} & 1 & -1.02\% & -1.36\% \\
 & \textbf{65.9} & 1 & -0.45\% & -0.66\% \\
 & \textbf{32.7} & 1 & 0.00\% & 0.00\% \\ \hline
\multicolumn{1}{r|}{\multirow{4}{*}{\textbf{$GR_{Near Wall}$}}} & 32.6 & \textbf{1.2} & 0.42\% & -0.08\% \\
\multicolumn{1}{r|}{} & 32.7 & \textbf{1.1} & 0.33\% & -0.11\% \\
\multicolumn{1}{r|}{} & 32.4 & \textbf{1.05} & -0.22\% & -0.32\% \\
\multicolumn{1}{r|}{} & 32.7 & \textbf{1} & 0.00\% & 0.00\% \\ \hline
\end{tabular}
\begin{tablenotes}
\small
\item The superscript c in $C_t^c$ and $C_p^c$ indicates that results are taken from the highest resolution for the given test variable.
\end{tablenotes}
\end{threeparttable}
\end{table}


Despite similar performance for integrated quantities for the range of $GR_{NearWall}$ tested, there is a significant difference in the pressure distribution around the blade in the range of $0.2<r/R<0.5$. The stability of the solution in this region is insufficient for the two most coarse meshes ($GR_{NearWall}=[1.10,1.20]$) but appears to be stable for the lower growth rates. This is demonstrated in Figure \ref{fig5:SMESH_GR}.

The mesh selected for this design adopted a $GR_{NearWall}$ of 1.05 and had a $\overline{y^+}$ of 66.



% plot of pressure distributions

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/SMESH_GR.png}
    \caption{Tangential (top) and stream-wise (bottom) force distribution per unit span along the length of the rotor blade for varied $GR_{NearWall}$.} 
    \label{fig5:SMESH_GR}
\end{figure}



\subsubsection{Semi-structured mesh convergence}

A semi-structured meshing approach was used in this study to develop a quick-to-build mesh with the capacity to have controlled refinement in the areas of interest. Primarily this method relied on developing a gridded face meshing on the blade surface (similar to structured approaches) and utilised inflation layers over the surface to generate a mesh capable of resolving a close-to-uniform $y^+$. Outside of the near-wall region, unstructured tetrahedral elements were used to fill the remaining model domain, reducing the user input time requirements. An example mesh demonstrating this parameterisation is presented in Figure \ref{fig5:Unstruc_mesh}. As can be observed in this figure the level of control over the near blade mesh resolution is nearly identical to that of structured methods. 

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/Stuctured_Mesh_Example.png}
    \caption{Mesh on blade surface (top). Cross section of mesh through blade geometry demonstrating the semi-structured mesh (bottom left). Zoomed view of blade tip showing inflation layers and the transition to unstructured mesh (bottom right). } 
    \label{fig5:Unstruc_mesh}
\end{figure}

A mesh convergence check, similar to that presented in the previous section, was undertaken. This study included additional variables to further examine the solution accuracy. These variables included: 
\begin{itemize}
    \item $N_{chord}$ - Number of elements in the 2D cross-section
    \item $GR_{chord}$ - Growth rate in the chord-wise direction
    \item $N_{radial}$ - Number of elements in the radial direction (along the blades length)
    \item $N_{blunt}$ - Number of element on the trailing edge
\end{itemize}


% table of performance

\begin{table}[h]
\centering

\begin{threeparttable}
\caption{Semi-structured mesh convergence checks}
\label{TAB5:unstrucConverge}
\begin{tabular}{l|rrrrrr|rr}
\textbf{Variable} & \textbf{$\overline{y^+}$} & \textbf{$N_{chord}$} & \textbf{$GR_{chord}$} & \textbf{$GR_{outer}$} & \textbf{$N_{radial}$} & \textbf{$N_{blunt}$} & \textbf{$\frac{C_{t}^{c}-C_{t}}{C_{t}^{c}}$} & \textbf{$\frac{C_{p}^{c}-C_{p}}{C_{p}^{c}}$} \\ \hline
\multirow{3}{*}{\textbf{$\overline{y^+}$}} & \textbf{120.7} & 160 & 1.1 & 1.2 & 175 & 4 & -0.95\% & -1.08\% \\
 & \textbf{60.8} & 160 & 1.1 & 1.2 & 175 & 4 & -0.40\% & -0.24\% \\
 & \textbf{30.3} & 160 & 1.1 & 1.2 & 175 & 4 & 0.00\% & 0.00\% \\ \hline
\multirow{3}{*}{\textbf{$N_{chord}$}} & 29.5 & \textbf{80} & \textbf{1.2} & 1.2 & 175 & 4 & 0.42\% & -0.10\% \\
 & 30.3 & \textbf{160} & \textbf{1.1} & 1.2 & 175 & 4 & 0.43\% & 0.46\% \\
 & 30.2 & \textbf{320} & \textbf{1.05} & 1.2 & 175 & 4 & 0.00\% & 0.00\% \\ \hline
\multirow{3}{*}{\textbf{$GR_{outer}$}} & 30.3 & 160 & 1.1 & \textbf{1.2} & 175 & 4 & 0.06\% & -0.17\% \\
 & 30.3 & 160 & 1.1 & \textbf{1.1} & 175 & 4 & 0.00\% & 0.00\% \\
 & N/A* & 160 & 1.1 & \textbf{1.05} & 175 & 4 & N/A* & N/A* \\ \hline
\multirow{2}{*}{\textbf{$N_{radial}$}} & 30.3 & 160 & 1.1 & 1.2 & \textbf{175} & 4 & -0.07\% & 0.04\% \\
 & 30.3 & 160 & 1.1 & 1.2 & \textbf{525} & 4 & 0.00\% & 0.00\% \\ \hline
\multirow{3}{*}{\textbf{$N_{blunt}$}} & 30.3 & 160 & 1.1 & 1.2 & 175 & \textbf{4} & 0.13\% & 0.16\% \\
 & 30.2 & 160 & 1.1 & 1.2 & 175 & \textbf{8} & 0.00\% & 0.00\% \\
 & N/A* & 160 & 1.1 & 1.2 & 175 & \textbf{16} & N/A* & N/A* \\ \hline
\end{tabular}
\begin{tablenotes}
\small
\item * Note that these simulations were unstable and results are not available.
\item The superscript c in $C_t^c$ and $C_p^c$ indicates that results are taken from the converged mesh, assumed to be the highest resolution for the given test variable.
\end{tablenotes}
\end{threeparttable}
\end{table}


Like the structured mesh, the range of parameters tested resulted in small changes in $C_p$ and $C_t$; these results are presented in Table \ref{TAB5:unstrucConverge}. Additionally, pressure distributions at radial stations and force distributions along the blade's length were examined for mesh convergence. Figure \ref{fig5:Unstruc_Pressure_chordN} and Figure \ref{fig5:Unstruc_force_ChordN} show the pressure and resultant force distribution along the blade's length for a varied $N_{chord}$. These plots demonstrate that force and pressure distributions are similar at the rotor's mid-span; however, there are some discrepancies between simulations towards the hub and the tip. Most noticeably, the force distribution in the stream-wise direction at $r/R=0.80$ in the lower resolution mesh predicted a higher force than the others. 

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/Unstruc_pressure_ChordN.png}
    \caption{Pressure distribution around at $r/R=0.35$ (left), $r/R=0.75$ (middle) and $r/R=0.99$ (right) for varied chord resolution using a semi-structured mesh.} 
    \label{fig5:Unstruc_Pressure_chordN}
\end{figure}

\begin{figure}[]
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/Unstruc_force_ChordN.png}
    \caption{Tangential (top) and stream-wise (bottom) force distribution along the length of the rotor blade for varied $N_{chord}$.} 
    \label{fig5:Unstruc_force_ChordN}
\end{figure}

Additionally, the suction side's pressure distribution at $r/R=0.99$ showed two suction peaks in the higher resolution meshes. The secondary peak in suction pressure is not a result of improved aerofoil performance; instead, it is due to pressure equalisation between the suction and pressure surfaces towards the middle chord due to the tip vortex.

The converged mesh used a $\overline{y^+}=30.3$, $N_{chord}=160$, $GR_{chord}=1.1$, $GR_{outer}=1.2$, $N_{radial}=175$ and $N_{blunt}=4$, for all these parameters there was $<0.5\%$ difference in $C_p$ when using a higher resolution model.




\subsubsection{Comparison between mesh approaches}

Results for $C_p$ and $C_t$ between the two mesh types were similar with the structured mesh having a $C_p=0.494$ and $C_t=0.871$ and the semi-structured mesh producing a $C_p=0.492$ and $C_t=0.871$ at the design $\lambda=4.5$. Examining the force and pressure distributions along the blade's radius (Figure \ref{fig5:Force_compare_mesh} and Figure \ref{fig5:Pressure_compare_mesh}), it is apparent that differences in the meshing approaches are subtle, with only a small difference in the force distribution at 0.85$R$. Given the similarity in the pressure distribution on the rotors, secondary factors were considered when selecting the appropriate mesh. Ultimately the semi-structured mesh was selected due to the simulations being ~25\% faster than the structured despite the semi-structured mesh containing more elements ($3.1 \times 10^6$ for the structured and $3.5 \times 10^6$ for the semi-structured). This result is specific to the geometry tested and is unlikely to be true for all other applications.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/Force_compare.png}
    \caption{Tangential (top) and stream-wise (bottom) force distribution along the length of the rotor blade for structured (orange) and semi-structured (blue) meshes.} 
    \label{fig5:Force_compare_mesh}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/pressure_dist1.png}
    \caption{Pressure distribution around at $r/R=0.35$ (left), $r/R=0.75$ (middle) and $r/R=0.99$ (right) for structured (orange) and semi-structured (blue) meshes.} 
    \label{fig5:Pressure_compare_mesh}
\end{figure}


\subsubsection{Turbine tip}
The effect of turbine tip shape was investigated by comparing the flat tip used in the mesh convergence study with a rounded tip. The rounding treatment used a semi-circle centred on the chord line to connect the suction and pressure surfaces; this resulted in a slight increase in blade radius from 4 m to 4.026 m. Results from the MRF simulations had no observable difference between the two tip treatments outside the tip region. Consequently, the change in $C_p$ and $C_t$ for the curved tip was small, with a 0.37\% reduction in power and a 0.63\% reduction in thrust compared to the flat tip. 

Figure \ref{fig5:Pressure_tip} and Figure \ref{fig5:vel_tip} show the static pressure distribution on the blade tip and the radial velocity near the blade. The pressure plot demonstrates that rounding the tip reduces the pressure on the pressure surface and increases pressure on the suction surface as radial flow along the blades permits pressure equalisation. This is confirmed in Figure \ref{fig5:vel_tip} which shows a slightly more significant impediment to the flow as it passes from the pressure to the blade's suction surface.

Ultimately, due to the flat blade tip's higher power generating capacity, this design was selected despite the marginally higher stream-wise thrust. 


% pressure plot
\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/Tip_pressure.png}
    \caption{Static pressure on the blade tip for flat (left) and curved (right) geometries.} 
    \label{fig5:Pressure_tip}
\end{figure}

% velocity plot
\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/Tip_radial_velocity.png}
    \caption{Radial velocity on a plane intersecting the 33\% chord and aligned with the free-stream flow near the blade tip (blade twist at the tip is 2.1$^\circ$) for flat (left) and curved (right) geometries. Free-stream flow direction is left to right.} 
    \label{fig5:vel_tip}
\end{figure}



\section{Sliding mesh}\label{CH5:sliding}

A sliding mesh was constructed by tessellating the 120-degree MRF mesh about the rotor hub to generate a complete turbine's geometry in a cylindrical channel. The domain was identical to that used in the RANS-BE model (Section \ref{sec5:RANS_DES}) except that the mesh defines the whole rotor geometry. A rotational mesh motion was defined for the inner domain to simulate the turbine operating at $\lambda=4.5$ and a sliding interface was prescribed between the inner and outer domains. This removed the need for azimuthal periodicity imposed from the MRF approach. The model was run using a second-order unsteady implicit solver ANSYS Fluent 19.2 \citep{ANSYSInc.2018}.

\subsection{Time-step convergence}
An analysis of time-step convergence was performed by running a series of simulations with varied time-step size and the number of iterations carried out per time-step. The results of this study are shown in Figure \ref{fig5:timestep_converge}. This figure presents the magnitude of the continuity residual and the real-time taken to complete one rotation of the turbine on a 16 core compute node from the Oxford ARCUS cluster with the marker's colour and size, respectively. This figure highlights how important the choice of time-step is to produce accurate results within an acceptable time frame. Further to this there is a limitation on the time-step size to prevent the sliding domain from jumping past cells any cells on the interface between the static and sliding domains.

Ultimately a time-step size of 0.00621 seconds with 20 iterations per time-step was selected (0.5$^\circ$ rotation per time-step) due to it balancing the performance with the simulations time, taking approximately 1 day to simulate a rotation. Results from this parameterisation appeared to converge within the first 10 rotations; an additional 5 rotations (15 total) were used to confirm this. The average CFL over the domain was 0.498; using the implicit scheme avoided restrictive conditions imposed by the CFL within the explicit solver.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/ContinuityResiduals.png}
    \caption{Magnitude of the continuity residual for varied time-step and number of iterations per time-step. Size of the marker reflects the real time taken to simulate a complete rotation of the turbine on a 16 core node on the Oxford ARCUS cluster.} 
    \label{fig5:timestep_converge}
\end{figure}



\subsection{Single turbine performance}
The sliding mesh results were similar to those generated from the MRF simulations, confirming that the MRF model's assumptions are valid for this model domain. Results for $C_p$ and $C_t$ were constant over the rotation, just marginally higher ($\sim$0.5\%) than those produced in the steady MRF simulations. Some difference is expected between the results as the solution switches from a steady to an unsteady solver. Given the close agreement between the models, the mesh and time-stepping method was deemed fit for purpose and suitable for further fence scale simulations. 


\section{Comparison of turbine models}


Three numerical methods of simulating tidal turbines have been used in this chapter to design an 8 m diameter turbine. As the numerical methods increased in complexity with each iteration (RANS-BE $\rightarrow$ MRF $\rightarrow$ Sliding mesh), the capacity to alter the design was reduced. Ultimately, this resulted in a close to optimal design (and practical) but not the best design for performance. This is due to the lower order models' underlying limitations, carrying design errors through to the higher-order models. Figure \ref{fig5:Compare_performance} plots the $C_p$ and $C_t$ predicted from each of the three model types used in this study. The plot shows that the MRF and sliding meshes have close agreement at the design $\lambda$ and highlights how the accuracy of the RANS-BE model changes (relative to the MRF model) over the range of $\lambda$. $C_t$ appears to be captured reasonably well by the RANS-BE model ($\pm 3\%$) at $\lambda>3$. Below $\lambda=3$ the RANS-BE model accuracy is reduced, most likely due to the increased angle of attack on the blade inducing stall; stall is not accurately captured by the 2D dataset or within a steady simulation. RANS-BE $C_p$ has a more significant error than $C_t$, which is expected given that $C_p$ is proportional to $u^3$ whereas $C_t$ is proportional to $u^2$, compounding any errors in the velocity in this parameter. $C_p$ is also sensitive to the angle of attack and reliant on the torque distribution along the blade, increasing errors further towards the blade's tip due to the higher lever arm. 


\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/Combined_perf_curves.png}
    \includegraphics[width=\textwidth]{Figures/Chapter_5/RANS_BE_error.png}
    \caption{$C_t$ (top left) and $C_p$ (top right) variation with $\lambda$ for the three model types used in this study. The bottom plot shows the error in the RANS-BE results relative to the MRF results.} 
    \label{fig5:Compare_performance}
\end{figure}



Figure \ref{fig5:Force_compare} plots the tangential and stream-wise force distribution along the blade's length for each of the model types used in this study. The MRF and the Sliding mesh models produce close to equivalent results; however, significant errors are observed in the RANS-BE model, most noticeable at the hub and the tip regions. The limitation of the RANS-BE model in capturing vortex shedding (both at the tip and the hub of a rotor) has been well studied with correction models being developed to try and account for these effects \citep{Hansen2004,Moriarty2005}. While these models are not perfect, they bring more physics into the model and improve the solution when used correctly. In this case, the Glauert tip correction appears to underestimate the losses from the rotor blade, leading to higher overall $C_p$ and $C_t$ for the RANS-BE case. 

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{Figures/Chapter_5/F_Compare.png}
    \caption{Force distribution along the blade length from RANS-BE (green), MRF (blue) and Sliding mesh (orange) simulations.} 
    \label{fig5:Force_compare}
\end{figure}


The RANS-BE results are only marginally higher across the mid-span of the rotor blade than the other methods. Given that the mid-span of the blade has the weakest span-wise flow components and, therefore, most closely resembles the 2D cases which drive the RANS-BE tool, this result is expected.

The RANS-BE model performs poorly around the hub. While this may be of low consequence for power generation (as torque production is low over this region), it further highlights the differences between the models. The RANS-BE model used in this study did not include a hub correction model meaning the drop in force observed at the hub is not a result of a vortex shedding. Instead, this observation may be credited to the velocity definition used in the RANS-BE model and the change in pressure distribution over the blade due to span-wise flows. Figure \ref{fig5:Vel_profile_compare} shows a comparison of the stream-wise flow component ($u_x$) for the RANS-BE model and the Sliding mesh model at 30$^\circ$, 60$^\circ$, 90$^\circ$ and 120$^\circ$. Note that the Sliding mesh model is periodic, so the selected values represent points at every 30$^\circ$ in a complete rotation, and the blade passes through the plane at 120$^\circ$.

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.85\textwidth]{Figures/Chapter_5/Profiles__x_velocity.png}
    \caption{Stream-wise velocity profiles from RANS-BE and Sliding mesh simulations around the turbine. The blade positions in the Sliding mesh results are at 0$^\circ$, 120$^\circ$ and 240$^\circ$, due to the radial symmetry in the problem data is plotted at 30$^\circ$ increments between 0-120$^\circ$ which is equivalent to 30$^\circ$ increments around the entire rotor. } 
    \label{fig5:Vel_profile_compare}
\end{figure}


Figure \ref{fig5:Vel_profile_compare} shows that the RANS-BE model does an excellent job of capturing the average stream-wise velocity; this is expected given the BE formulation performs azimuthal averaging. What is shown in this plot by the Sliding mesh results is that after the blade passes (30$^\circ$) to just before the subsequent blade passing (90$^\circ$), the stream-wise velocity recovers from the previous blade passing. This exposes the incoming blade to a higher stream-wise velocity than the average (used in the RANS-BE). As the velocity at the blade is determined from stream-wise and rotational contributions following  

\begin{equation}
    \begin{aligned}
       W = \left[ \left(\Omega r + u_{\theta} \right)^{2} +u_{x}^{2} \right]^{1/2}
    \end{aligned}
\end{equation}

\noindent this effect is more prominent towards the hub as the angular velocity component dominates outside this region. As a result of the underestimated stream-wise velocity, the angle of attack is lower in the RANS-BE model than the MRF and sliding meshes.

Furthermore, due to the span-wise flow, the pressure distribution around the 2D aerofoil is changed; this is clear when examining the $C_{pressure}$ variation at different radial stations from the sliding mesh model (Figure \ref{fig5:Pressure_dist_compare}). Towards the tip and the turbine's hub, it is clear that the pressure distribution does not match that for the 2D section at the design angle of attack ($\alpha=6^{\circ}$). At the tip, $\alpha$ decreases to 4$^\circ$ and at the hub, the shape of the pressure distribution is completely changed. Presumably, this is due to the span-wise flow component being significant when compared to the relative velocity of the blade ($W$), modifying the pressure distribution.

% this means that streamlines passing over the blade do not pass over the 2D geometry which has been projected in annular rings, instead streamlines experience a distorted 3D geometry as they move in both the radial and tangential directions.

\begin{figure}
    \centering
    \includegraphics[width=0.9\textwidth]{Figures/Chapter_5/C_pressure_RANS_SLIDING.png}
    \caption{$C_{pressure}$ variation towards the hub ($r/R=0.135$), mid-blade ($r/R=0.65$) and tip ($r/R=0.95$) from the Sliding mesh simulations plotted alongside the $C_{pressure}$ from the 2D simulations at the turbine design target $\alpha$.} 
    \label{fig5:Pressure_dist_compare}
\end{figure}




\section{Conclusions}

A design for an 8 m diameter rotor has been developed to operate in a global blockage of 1.2\%. The design approach adopted gradually increasing model complexity with each design iteration considering more practical constraints to manufacture and aiming to fix design parameters in place before using the more complex models to refine the design. 

2D airfoil data was generated for a NACA-63-415 aerofoil for varied trailing edge treatments operating in a range of turbulence intensities and Reynolds numbers. This data was input to a modified RANS-BE design tool that interpolated for the inflow conditions and between aerofoil profiles to maintain a constant trailing edge thickness of 2 mm. Analysis of the pressure distribution on the aerofoils determined that cavitation initiation may be a limiting factor in the turbine's design operating speed, and an upper limit for $\lambda$ was defined through a simplification of \cite{Wimshurst2018} which took conservative estimates for flow speed relative to the turbine blades and the critical location for the inception of cavitation.

The RANS-BE design was used and further refined in MRF simulations, where both structured and semi-structured meshes were compared. The difference in blade loading between the structured and semi-structured meshes was negligible when the appropriate inflation layers parameterisation was used. This result validates the use of the semi-structured meshing with near-wall inflation layers for the application to tidal turbine rotors.

Sliding mesh simulations were completed to validate the MRF model and determine the appropriate time-stepping for the transient model. A time-step of 0.00621 s (0.5$^\circ$ rotation) with 20 iterations per time-step was sufficient for this application matching the MRF simulation results closely. 

Comparing the three models used in this study highlighted the RANS-BE model's limitations in capturing the physics of the flow near the blade root and tip within a tidal turbine rotor. Specifically, the RANS-BE model's failure to simulate span-wise flows (which violate the underlying assumptions of BE theory), the unsteady fluid mechanics and skewed inflows relative to the turbine plane were all contributors to the error in the RANS-BE model.

The rotor operates at a design $\lambda$ of 4.5 and achieves a $C_p$ of 0.494 and $C_t$ of 0.871 when in a cylindrical channel with a global blockage of 1.2\%. This design meets the requirement of fully turbulent conditions over the blades and is practical for further blade resolved analysis at the fence scale using the $k-\omega$ SST turbulence closure scheme. The design is derived from a first pass through the design procedure. Further optimisation of the rotor is possible, especially at the root and the tip; however, this would come at a high computational cost as lower-order RANS-BE methods do not adequately capture the physical processes in these regions and a more computationally expensive method, such as the MRF technique, would be required.  