\chapter{Simulation of tidal turbines}\label{CH3:Model_turb}

This chapter presents some common methodologies for simulating tidal turbines which are used in this thesis. There are two broad classes of models for simulating tidal turbines, analytical - where assumptions are made about the turbines to simplify the model, and numerical - where at least some part of the analytical problem is solved numerically. This chapter explores analytical and numerical models, discussing the incremental changes in the models as complexity is increased and accuracy is improved.

\section{Analytical models}

Analytical models are typically a low computational cost method for simulating turbines. These models tend to be steady approximations, requiring the resistance from the turbine blades to be distributed over the entire swept area. Computationally they are treated as an actuator disk with a surface of negligible thickness that imposes, at a minimum, a momentum loss term on the fluid domain. The turbine model is usually embedded with a secondary set of equations to simulate the fluid flow. These include the momentum, continuity and energy equations \citep{GARRETT2007,Nishino2013}, the shallow water equations \citep{Bonar2018} and Reynolds-Averaged Navier-Stokes (RANS) equations \citep{Hunter2015}. 

The practical application of these models are to provide insights by simulating the leading order physics of the system. Using this model \cite{GARRETT2007} were able to demonstrate the performance improvement when turbines are placed in a constrained channel (blockage effect).

\subsection{Porous disk models}

The porous disk model is a type of actuator disk model which does not require the specification of turbine geometry or rotational speed, instead simplifying the problem by using a single resistance coefficient ($K$) to determine the thrust ($T$) a turbine imposes on the flow following,

\begin{equation}
    T=K\frac{\rho}{2}u_x|u_x|A_t
\end{equation}

\noindent where $\rho$ is the density, $u_x$ is the velocity through the rotor and $A_t$ is the swept area of the turbine. An estimate for the power production of porous disk is often made by assuming that there are no efficiency losses of the turbine allowing power ($P$) to be calculated directly from the turbine thrust,

\begin{equation}
    P=Tu_x
\end{equation}

When incorporated within the momentum equations, the porous disk model gives rise to linear momentum actuator disk theory (LMADT) which has been used to evaluate the theoretical maximum efficiency of a turbine \citep{GARRETT2007,Houlsby2008,Nishino2012}. This model is popular as it removes turbine geometry as a variable, replacing the unconstrained geometry problem with a single parameter which can be iterated over to determine optimal turbine performance.

\subsection{Blade element models}\label{CH2:BE_model}

Blade element (BE) models are commonly used in the design and performance analysis of tidal turbine rotors \citep{Malki2013,Olczak2016,Belloni2017,Cao2019}. The BE method uses lift and drag data from 2D aerofoil sections at varied angles of attack ($\alpha$) in conjunction with the blade twist angle ($\beta$) and chord length ($c$) to approximate the forces at different radial stations on the blade. This imposes a pressure drop ($\Delta p$) and swirl velocity ($u_{\theta}$) across the actuator disc. The use of 2D data assumes no radial flow component along the span of the rotor which is typically violated by span-wise pressure gradients on rotor blades. Despite this, the flow component can be considered small relative to the primary inflow velocity over the majority of the blade with the exclusion of the hub and the tip region; correction models exist in these regions and are discussed further in this section.

% Fix the scaling on this figure
\begin{figure}[h]
\centering
\includegraphics[width=0.8\textwidth]{Figures/Chapter_3/numerical/BEM_001.png}
\caption{(a): Rotor discretisation; (b): resolution of lift and drag forces at radius $r$.}
\label{Fig3:BEM}
\end{figure}

The BE equations can be derived by considering the velocity triangle in Figure \ref{Fig3:BEM} at a point along the blade and equating the relevant forces. This yields the following expressions for pressure drop and swirl velocity for an element at radial location $r$ within the swept area: 

\begin{equation} 
\label{EqPressureDrop}
\begin{aligned}
\Delta p &=\frac{1}{2} \rho W^{2} \sigma\left(c_{l} \cos \phi+c_{d} \sin \phi\right) ,
\end{aligned}
\end{equation}

\begin{equation}
\label{EqSwirlVel}
\begin{aligned} 
u_{\theta} &= \frac{W^{2} \sigma\left(c_{l} \sin \phi-c_{d} \cos \phi\right)}{4u_{x}}, \end{aligned}
\end{equation}

\noindent where $W$ is the relative speed of the flow to the blade, $\phi$ is the flow angle between the approach flow and rotor plane (which is the sum of $\alpha$ and $\beta$), $c_l$ is the 2-dimensional lift coefficient at $\alpha$, $c_d$ is the 2-dimensional drag coefficient at $\alpha$, $u_x$ is the stream-wise velocity component at the radial position $r$ within the rotor plane and $\sigma$ is the blade solidity defined as 

\begin{equation}
    \label{EqSolidity}
    \begin{aligned}
        \sigma = \frac{N_b c}{2 \pi r}
    \end{aligned}
\end{equation}

\noindent where $N_b$ is the number of blades and $c$ is the chord length at radial position $r$. 
% The $c_l$ and $c_d$ are the non-dimensionalised lift and drag forces on the section defined as

% \begin{equation}
%     \begin{aligned}
%         c_l = \frac{F_l}{\tfrac{1}{2}\rho W^2A}, \; \; \; \; c_d = \frac{F_d}{\tfrac{1}{2}\rho W^2A}
%     \end{aligned}
% \end{equation}

% \noindent where $F_l$ is the lift force, $F_d$ is the drag force and $A$ is the representative area (taken as the chord length ($c$) in the 2D simulations). 
$W$ can be determined by analysing the components of the flow and the rotation of the blades to arrive at the expression,

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

\noindent where $\Omega$ is the rotational speed of the turbine.

The derivation of the BE model implies that fluxes are constrained to each concentric streamtube as flow is only permitted in the stream-wise and azimuthal directions, not allowing an exchange of radial momentum. In reality a radial pressure gradient is introduced between the pressure and suction surfaces of the hydrofoil at the blade tip driving a spanwise flow \citep{Masters2011}. Various tip correction models attempt to account for the reduction in power resulting from this flow, an overview of the performance of some common implementation can be found in \cite{Koh2016}. For this study the commonly used tip correction model of \cite{Glauert} was included in the BE model; this factor reduces the loading on the outboard blade span following:

\begin{equation}\label{eq3:TCF}
    \begin{aligned}
        F_{tip} = \frac{2}{\pi} \cos^{-1}\left[\exp{\frac{N_b\left(1-\tfrac{R}{r}\right)}{2 \sin{\phi}}}\right] 
    \end{aligned}
\end{equation}

\noindent where $R$ is the turbine radius. 

Like the porous disk model, the blade element model can be embedded within the momentum equations. This gives rise to the blade element moment (BEM) model. Various forms of this model exist \citep{Masters2011,Galloway2014,Vogel2018} which serve as quick estimates of turbine power production and allow the effect of basic control strategies such as pitch-to-feather to be investigated. One noteworthy example is proposed by \citet{Vogel2018} who analytically accounted for the effects of blockage in BEM theory. \citet{Vogel2018} verified the model by showing that close agreement with power output (within 3\% of blade resolved simulation) can be achieved using the computationally inexpensive constrained BEM method. 

The implementation of the BE model within this thesis embeds the BE equations into a RANS model (RANS-BE) using a user defined function (UDF) within ANSYS Fluent 15.0 \cite{ANSYSInc.2013}. When called, the UDF applies a pressure drop and swirl velocity across a cell face within the turbine swept area as a boundary condition for the RANS model. This approach has previously been used by \cite{McIntosh2011} and \cite{Cao2019} to simulate a turbine within a highly blocked channel. Further refinement of this model is detailed in Chapter \ref{Chapter5} of this thesis. Details of the RANS model can be found in Section \ref{CH2:RANS}.

Careful consideration of the model limitations have to be made prior to the implementation. As discussed previously in this section the set up of the model restricts radial flow and imposes a permanently present resistance at all radial locations in the rotor's swept area. The later restriction limits the use of the model to uniform inflow conditions normal to the rotor plane, however some implementations of BEM models have been shown to have favourable results in comparison to experimental data \citep{Galloway2014} when these restrictions are relaxed. This suggests a similar implementation may be possible in a RANS-BE model. Nevertheless, the restriction to steady conditions limits the ability to simulate turbines under wave motion or rapidly varying flow conditions. Furthermore the restriction limiting the radial exchange of momentum along the blades also has implications for the flow conditions the RANS-BE model can be used in. Conditions where the flow is not orthogonal to the rotor, such as yawed inflow, induce a radial component to the flow, impacting the accuracy of the model. 


\subsection{Actuator line models}\label{AL_model}


Actuator line models are an unsteady implementation of the blade element equations which permit an approximation for the temporal and spatial variation of turbine loading and performance to be modelled. The model approximates turbine blades as a number of discrete segments which experience 2-dimensional flow locally. At each timestep, the inflow velocity is determined based on the relative motion between the flow field and blade motion 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 and the angle of attack ($\alpha$) is calculated following Figure \ref{Fig3:BEM}. Individual element forces can then be determined through looking up $c_l$ and $c_d$ from a set of 2-dimensional aerofoil data, allowing a blade force and drag to be determined for a given actuator line section following


\begin{equation}
\label{AL_forces}
    \begin{aligned}
      L_m = \frac{1}{2} \rho W^2  c_l(\alpha) c \Delta r \\
      D_m = \frac{1}{2} \rho W^2  c_d(\alpha) c \Delta r
    \end{aligned}
\end{equation}

\noindent where $L_m$ and $D_m$ is the lift and drag on element $m$, respectively, and $\Delta r$ is the length of the actuator line segment in the radial direction. 


Unlike the BE method these forces are not smeared over the rotor's swept area, instead, forces are applied over a spatially and temporally varying region. Given the model resolves the forces along a discrete line, the application of this line source in a numerical model may lead to some numerical instabilities. For this reason, a technique of smearing of the force over adjacent cells is used. This is achieved by applying a scaled force over the adjacent cells using a Gaussian kernel function, $\eta$, defined as 


\begin{equation}
    \eta_m = \frac{1}{\epsilon^3 \pi^{3/2} \exp{ \left[  -(\frac{\mathbf{p}}{\epsilon})^2 \right]  }   }
\end{equation}

\noindent where $\epsilon$ controls the width of the Gaussian and $\mathbf{p}$ is the distance vector between the grid point and the actuator line element. The application of the Gaussian is such that the integrated force applied over the adjacent cells matches that of the forces defined from Equation \ref{AL_forces}. This process is repeated every timestep within the simulation. 



The actuator line method has been implemented in studies of both wind \citep{Schmitz2013ModelingSolver} and tidal \citep{Baratchi2017} turbines. With both models capable of reproducing observed physical phenomenon, such as the formation of a tip vortex. Despite the actuator line models promising performance, it has not been used in this thesis due to the additional simulation time required to achieve converged results that are periodic in time. 



\section{Numerical models and methods}

As discussed in the above section analytical models are useful for resolving the leading order physics of a tidal turbine. However, additional factors such as turbulence, tip vortices and sheared inflow profiles, which are not readily resolved in these models, are often important when designing and determining the performance of turbines. In these conditions the use of a numerical model is often required. While there are many methods of numerically solving this problem (such as Large Eddy Simulations and Direct Numerical Simulation) only the Reynolds-Averaged Navier-Stokes is discussed in this section. Two approaches for turbine simulations within this framework are discussed in this section, the Multiple Reference Frames (MRF) and the Sliding mesh approach. 

\subsection{Reynolds-Averaged Navier-Stokes equations}\label{CH2:RANS}

The Reynolds-Averaged Navier-Stokes (RANS) equations are a fundamental component of modern computational fluid dynamics (CFD). These equations are a time-averaged version of the Navier-Stokes equations which themselves are a per-volume expression of Newton's second law:

\begin{equation}
    \frac{ma}{V} = \frac{\sum F}{V}
\end{equation}

\noindent where $\sum F$ are the forces acting on a control volume of fluid, $a$ is the acceleration of the fluid, $V$ is the volume of the fluid element and $m$ is the fluid mass. The fluid density, $\rho$ can be used in place of the $\frac{m}{V}$ fraction. Considering the shear and normal stresses drawn on the infinitesimally small control volume in Figure \ref{fig:3-Stresses} and making the assumption of an isotropic incompressible fluid, an expression for fluid motion can be derived:


\begin{equation}
    \label{EQ:NS}
    \rho \frac{\mathrm{D} u }{\mathrm{D} t} = \rho \textbf{g} -\bigtriangledown T
\end{equation}

\noindent where $u$ is the velocity, $g$ is the acceleration due to gravity, $\frac{\mathrm{D} }{\mathrm{D} t}$ is the material derivative consisting of both temporal and convective terms, $\bigtriangledown$ is the grad operator and $T$ is the Cauchy stress tensor. In 3D $T_c$ is defined as

\begin{equation}
    T_c = \begin{bmatrix}
        \sigma_{x} & \tau_{xy} & \tau_{xz} \\
        \tau_{yx} & \sigma_{y} & \tau_{yz} \\
        \tau_{zx} & \tau_{zy} & \sigma_{z} \\
    \end{bmatrix}
\end{equation} 

\noindent where $\sigma$ are normal forces and $\tau$ are shear forces. The expression for motion prescribed by Equation \ref{EQ:NS} is useful as it describes how fluid will behave when stresses act upon it. However, knowing what these stresses are requires the formation of a constitutive relation between the stress and strain rate of the fluid (transverse velocity) through the viscosity ($\mu$). Details of how this constitutive law is applied can be found in  \cite{Schlichting2000}. This leads to the more common expression for the Naiver-Stokes equation as:

\begin{equation}
    \rho \frac{\mathrm{D} u }{\mathrm{D} t} = \rho \textbf{g} -\bigtriangledown p +\mu \bigtriangledown^2 u 
\end{equation}

\noindent where $p$ is the fluids static pressure.

\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{Figures/Chapter_3/numerical/stress_state.png}
    \caption{Stresses acting on an infinitesimally small control volume in 2D}
    \label{fig:3-Stresses}
\end{figure}

When appropriately discretised, the equation presented above is capable of describing fluid flow across all length and time scales with a simulation methodology known as direct numerical simulation (DNS). This is exceptionally computationally intensive both in terms of memory and CPU requirements. For this reason additional assumptions are made to reduce the complexity of the problem. The approach adopted for the RANS equations is to decompose the velocity into a time-averaged and a fluctuating components such that

\begin{equation}
    \mathbf{u}(\mathbf{x},t) = \overline{\mathbf{u}}(\mathbf{x}) + \mathbf{u}^{\prime}(\mathbf{x},t),\;   \overline{\mathbf{u}^{\prime}} = 0
\end{equation}

\noindent where $\mathbf{u}$ is the velocity vector which has a component direction ($\mathbf{x}$), primes represent fluctuating components and the overbar indicates a time average. Substituting this expression into Equation \ref{EQ:NS} yields

\begin{equation}
    \rho \frac{\mathrm{D} \overline{u} }{\mathrm{D} t} = \rho \textbf{g} - \frac{\partial p}{\partial x_{i}} +  \frac{\partial}{\partial x_{j}}(\mu \frac{\partial \overline{u_i}}{\partial x_j} - \rho \overline{u^{\prime}_i u^{\prime}_j})
\end{equation}

\noindent the $\rho  \overline{u^{\prime}_i u^{\prime}_j}$ is commonly referred to as the Reynolds stress tensor. Similar to $T_c$ this expression describes how fluid moves when stresses are applied, in this case the stresses arrive from fluctuations in the flow. However, unlike $T_c$ there is not a constitutive law which can be applied. This is to say the velocity fluctuations are not directly related to the time average velocity. This leads to a closure problem and the need for additional models to account for the effect of turbulence.

\subsubsection{Turbulence modelling}\label{Ch2:turbulence}

Various models exist for approximating turbulence in flows dependant on the application. One of the most common turbulence models simulates the transport and dissipation of energy using an eddy viscosity. This is achieved by applying the Boussinesq hypothesis to relate the turbulent stresses to the mean velocity gradients similar to how the viscous stresses are related to the complete velocity gradient by

\begin{equation}
-\overline{u_{i}^{\prime} u_{j}^{\prime}}=\nu_{t}\left(\frac{\partial \overline{u_{i}}}{\partial x_{j}}+\frac{\partial \overline{u_{j}}}{\partial x_{i}}\right)-\frac{2}{3} k \delta_{i j}
\end{equation}

\noindent where $\nu_{t}$ is the eddy viscosity, $k$ is the turbulent kinetic energy ($k=\tfrac{1}{2} \overline{u_{i}^{\prime}u_{i}^{\prime}}$) and $\delta_{i j}$ is the Kronecker delta. 

This relation extends into the 2-equation $k-\epsilon$ turbulence model which tracks the transport of $k$ and $\epsilon$ (dissipation rate of $k$). The $k-\epsilon$ model is commonly used in industrial and environmental flow applications without adverse pressure gradients \citep{Mansour1988}. In the presence of an adverse pressure gradient (such as flow over aerofoils) the $k-\omega$ turbulence model ($\omega$ in this model is the \emph{specific} dissipation rate of $k$) has been shown to be more reliable. While better in the near-wall region the $k-\omega$ model has been noted to be sensitive to the free-stream turbulence conditions \citep{Kok2000} and unreliable away from the wall region. 


A model by \cite{Menter1994} introduced a blending function between the $k-\omega$ and the $k-\epsilon$ models to leverage the best of the $k-\omega$ model, which performs well in the near-wall region, and the $k-\epsilon$ in the far field \citep{ANSYSInc.2013}. This model also introduces a limiter on the production of $k$ which was observed to be good in adverse pressure gradient regions \citep{Menter1994}. Menter's $k-\omega$ shear stress transport (SST) model has previously been used to model flow over an airfoil profile with a favourable comparison to observations \citep{Menter1994} and to directly simulate flow through a tidal rotor \citep{Cao2019}.  Due to these benefits and the prior applications the $k-\omega$ SST model was used throughout this study. This model does have some limitations in the application to tidal turbines which have to be carefully considered. Specifically this model does not account for transitional flows and care must be taken to check that the flow is fully turbulent over the entire turbine blade. If it is not, then a four equation model such as the $\overline{Re_{\theta,t}},\gamma$ may be required \citep{Menter2006AFormulation}. This model introduces two additional transport variables $\overline{Re_{\theta,t}}$ (momentum thickness Reynolds number) and $\gamma$ (fraction of time where the flow is fully turbulent) to simulate transitional flows.

\subsection{Numerical Solution Process}


\noindent \underline{Solver} 


The three-dimensional RANS equations were solved using a pressure-based finite volume solver by the commercial software ANSYS Fluent 19.2 \citep{ANSYSInc.2018}. The solution approach first solves the velocity field then uses a pressure correction, derived from the continuity and momentum equations, such that the velocity field is corrected by the pressure. This correction is then applied back onto the model and the process is repeated until the convergence criteria is met. This convergence criteria is based on the scaled residuals of the conserved variables. These residuals are defined as the average magnitude of the error within each cell (for example continuity error) scaled by a term representative of the flux of the variable through the domain. For this study a residual magnitude of $<1\mathrm{e}{-10}$ was targeted. 

\noindent \underline{Spatial discretisation} 

A requirement of the solver is a spatial discretisation of the flow domain. This was performed using a meshing approach with ANSYS ICEM \citep{ANSYSInc.2013} and ANSYS Workbench \citep{ANSYSInc.2018} for structured and unstructured meshes respectively. The finite volume method employed by ANSYS Fluent applies the integral form of the conservation equations to each cell as a control volume. The integral form can be written as


\begin{equation}
    \int_{V} \frac{\partial \rho \phi}{\partial t} d V+\oint \rho \phi \mathbf{u} \cdot d \mathbf{A}=\oint \Gamma_{\phi} \nabla \phi \cdot d \mathbf{A}+\int_{V} S_{\phi} d V
\end{equation}

\noindent for a scalar quantity $\phi$ within a control volume $V$, where $\Gamma_{\phi}$ is a diffusion coefficient for $\phi$, $S_{\phi}$ is the volumetric source term and $\mathbf{A}$ is the surface area vector around the control volume such that $\mathbf{A}=A \cdot \mathbf{n}$, where $A$ is the cell area and $\mathbf{n}$ is the surface normal vector. This equation is applied to the control volume such that 


\begin{equation}
    \frac{\partial \rho \phi}{\partial t} V+\sum_{f}^{N_{\text {faces }}} \rho \phi_{f} \mathbf{u}_{f} \cdot \mathbf{A}_{f}=\sum_{f}^{N_{\text {faces }}} \Gamma_{\phi} \nabla \phi_{f} \cdot \mathbf{A}_{f}+S_{\phi} V
\end{equation}

\noindent where $N_{\text{faces}}$ is the number of cell faces, $f$ is the face and the sub-script $f$ refers to the value passing through the face. Values of $\phi$ are stored at the cell centre however $\phi$ is required at the cell faces. The solver is therefore required to interpolate values to the cell faces. This is achieved using a second-order upwind scheme in ANSYS Fluent.


\noindent \underline{Temporal discretisation} 


Sometimes features of the flow can be considered steady in the time, however in certain cases, such as explicitly resolving the flow around a turbine blade, the unsteady flow features are of interest. This requires a discretisation of the solver in the time domain, both steady and unsteady models are considered in this thesis.

\noindent The time evolution of $\phi$ can be expressed as

\begin{equation}
    \frac{\partial \phi}{\partial t}=F(\phi)
\end{equation}

\noindent where the function $F$ incorporates any spatial discretisation. This can be approximated to the first and second-order using


\begin{align}
    \label{EQ-time_dis}
    \frac{\phi^{n+1}-\phi^{n}}{\Delta t}&=F(\phi)   &\mathrm{and}&   & \frac{3\phi^{n+1}-4\phi^{n}+\phi^{n-1}}{2\Delta t}&=F(\phi) 
\end{align}

\noindent where $\Delta t$ is the time-step size, $n$ is the value at the current time-step and $n+1$ and $n-1$ are the values at the following and previous time-steps, respectively. For a forward marching explicit time integration scheme, the solver is restricted to a maximum time-step size defined by the Courant-Friedrich-Lewy condition (CFL). This condition states that 


\begin{equation}
    C = \Delta t \sum_{i=1}^{N} \frac{u_{x_i}}{\Delta x_i} \leq 1
\end{equation}

\noindent where C is the Courant number, $\Delta x_i$ is the length across the cell in the $x_i$ direction and $u_{x_i}$ is the velocity in the $x_i$ direction. This condition is not restrictive when using an iterative approach to solve for time (implicit scheme). This can be achieved using a backward differencing appraoch such that Equation \ref{EQ-time_dis} becomes


\begin{equation}
    \frac{\phi^{n+1}-\phi^{n}}{\Delta t}=F(\phi^{n+1}).
\end{equation}

\noindent A second-order implicit approach was used for the unsteady simulation in this thesis. Further information on the numerical scheme can be found in the Fluent science manual \citep{ANSYSInc.2018}.



\subsection{Multiple reference frames}

A common way to simulate a turbine within the RANS framework in a computationally efficient way is to solve the equations over multiple reference frames (MRF). This method reduces the computational complexity of the model and allows a steady solver to be used by separating the computational domain into stationary and a moving reference frames. For the application to turbines the inner region is solved in a rotating reference frame equal to the turbine's angular speed, and the outer region (encasing the inner region) is solved in the absolute reference frame; a transformation between the reference frames is applied at zone boundaries allowing fluxes between them. Using this method ensures that the rotor remains still in the simulation allowing a single mesh to be used. If the outer domain is also radially symmetric, the mesh can be further simplified by reducing it to a wedge containing one blade with periodic boundary conditions at the sides, significantly reducing the model size.


Within the moving domain a relative velocity formulation is used, transforming the following
\begin{equation}
    \mathbf{u_{rel}} = \mathbf{u} - \omega \times \mathbf{r}
\end{equation}

\noindent where $\mathbf{u_{rel}}$ is the relative velocity vector, $\mathbf{u}$ is the velocity vector in the stationary reference frame, $\omega$ is the angular velocity of the internal reference frame and $\mathbf{r}$ is the position vector from the origin of the moving reference frame.


While this method is computationally efficient, the application is limited to radially symmetric inner domains that are steady in time. Consequently this method is unable to fully resolve the helical vortex which forms at the blade tip, vibration loading \citep{Liu2016} and is not applicable in anisotropic blockage conditions or in sheared profiles. Nevertheless this method is a useful tool for sensitivity analysis of designs \citep{Siddiqui2016} and in the initialisation of sliding mesh simulations due to the speed of the model \citep{ANSYSInc.2013}. 



\subsection{Sliding mesh}

The sliding mesh approach is the most computationally expensive method and accurate method for the simulation of a turbine presented in this thesis. This transient model explicitly resolves the rotation of a turbine through the decomposition of the domain into rotating and stationary sub-domains, similar to the MRF model. Unlike the MRF model the inner domain rotates, sliding within the outer domain. Consequently, the sliding mesh interface is nonconformal between the domains requiring fluxes to be interpolated across the interface. Care is required in the mesh development at the interface zones to ensure that the resolution on either side of the interface is approximately the same to assist in the interpolation. The sliding mesh is capable of resolving the higher detailed physical flow processes around a tidal turbine such as the helical vortex at the blade tip and flow separation. Unlike the MRF approach, this methodology is able to capture non-radially symmetric flow features, such as a rectangular domain and spatially varying inflow conditions. 

