Commit 37954e52 authored by Pierre Cazenave's avatar Pierre Cazenave

Merge branch 'dev' of gitlab.ecosystem-modelling.pml.ac.uk:fvcom/fvcom-toolbox into dev

parents c9a4dbbf 60402a31
%
% Example code showing how to make a tidal turbine parameter input netCDF
% file for FVCOM.
%
% Rory O'Hara Murray
% 6 Oct 2015
%
clear, close all
% load in PFOW model grid and bathymetry
Mobj = read_fvcom_mesh(['\\isilonml\Shelf_Model\PFOW\1_HD Model\4_Climatology Results\input\PFOW_SMS_1_grd.dat']);
Mobj.h = read_fvcom_bath('\\isilonml\Shelf_Model\PFOW\1_HD Model\4_Climatology Results\input\PFOW_SMS_1_dep.dat');
lonc = nodes2elems(Mobj.x,Mobj); % calculate approximate x-coordinate at the centre of each element
latc = nodes2elems(Mobj.y,Mobj); % calculate approximate y-coordinate at the centre of each element
[xc yc] = ll2utm(lonc, latc, 30);
% calculate depth at the centre of each element
Mobj.hc = mean(Mobj.h(Mobj.tri),2);
% setup empty turbine array
turbine.numbers = single(zeros(Mobj.nElems,1));
turbine.sigma_layer = single(zeros(Mobj.nElems,10));
turbine.area = single(zeros(Mobj.nElems,1));
% turbine.thrust = single(zeros(Mobj.nElems,1));
% turbine.drag = single(zeros(Mobj.nElems,1));
% Example turbine positions in the Inner Sound of Stroma
TurbinePositionsll = [-3.137,58.661; ...
-3.12945337127092,58.6596303572796; ...
-3.12960283663035,58.6600269324843; ...
-3.13193517522843,58.658758560404; ...
-3.1320846649985,58.659155132772];
[TurbinePositions(:,1) TurbinePositions(:,2)] = ll2utm(TurbinePositionsll(:,1), TurbinePositionsll(:,2), 30);
% For each turbine find which element it is in and count the number of
% turbines in each element
I = fun_nearest2D(TurbinePositions(:,1), TurbinePositions(:,2), xc, yc);
for ii=1:length(I)
turbine.numbers(I(ii)) = turbine.numbers(I(ii)) + 1;
end
turbinesI = turbine.numbers>0;
numTurbineElems = sum(turbinesI); % total number of elements with turbines in
numTurbines = sum(turbine.numbers); % total number of turbines
% loop through each turbine and work out the fraction spread over each sigma layer
for ii=find(turbinesI)'
turbine.sigma_layer(ii,:) = turbine_area_sigma(Mobj.hc(ii), 15, 10);
end
% an alternative (simpler) way of specifying the fractional split across
% sigma layers
% turbine.sigma_layer(turbinesI,:) = ones(numTurbineElems,1)*[0 0 0 0 0 0.25 0.25 0.25 0.25 0];
% turbine rotor swept area - 10 m radius.
turbine.area(turbinesI) = 10*10*pi;
% The thrust curve is currently not inputted into FVCOM, this could be
% added in the future though.
% turbine.thrust_curve = [0.99 1 1.25 1.5 1.75 2 2.25 2.5 2.75 3 3.25 3.5 3.75 4 4.01; 0 0.85 0.85 0.85 0.85 0.85 0.85 0.85 0.635 0.490 0.385 0.308 0.250 0.205 0];
% write the netcdf input file
write_FVCOM_TT(turbine,['Tidal_Turbines_Example.nc'],'Example Scenario with 5 tidal turbines in the inner sound');
%% plot the location of the turbines
%plot_fvcom_field(Mobj,elems2nodes(turbine.numbers,Mobj), 'grd', 'w')
plot_fvcom_field(Mobj,single(turbinesI), 'grd', 'w')
colormap( jet(2) );
daspect([1 1 1])
function write_FVCOM_TT3(turbine,filename,mytitle)
% Dump tidal turbine parameters to FVCOM forcing file.
%
% function write_FVCOM_TT(turbine,filename,mytitle)
%
% DESCRIPTION:
% Generate a NetCDF file containing tidal turbine parameters for FVCOM.
% netCDF variable names:
% turbine_sigma_layer, turbine_number, swept_area
%
% INPUT:
% turbine.numbers = user defined array of integers specify the number of turbines in each element
% defined on the elements
% turbine.sigma_layer = user defined array specifying the fractional
% split across sigma layers that the turbine rotor swept area occupies
% turbine.area = rotor swept area (m^2)
% filename = filename to dump to
% mytitle = title of the case (set as global attribute)
%
% EXAMPLE USAGE:
% write_FVCOM_TT(turbine, 'tst_file.nc', 'tst tidal turbine parameters')
%
% Author(s):
% Rory O'Hara Murray (MSS)
%
% Revision history
%
%==============================================================================
%warning off
subname = 'write_FVCOM_TT';
global ftbverbose;
if(ftbverbose);
fprintf('\n'); fprintf(['begin : ' subname '\n']);
end;
% check dimensions
nElems = numel(turbine.numbers);
nSiglay = size(turbine.sigma_layer,2);
if(nElems == 0)
error('dimension of turbine_numbers is 0, something is wrong ')
end;
%------------------------------------------------------------------------------
% Dump to turbine NetCDF file
%------------------------------------------------------------------------------
if(ftbverbose);
fprintf('Dumping to turbine parameters NetCDF file: %s\n',filename);
fprintf('Size of turbine_numbers array: %i\n',nElems);
end;
nc = netcdf.create(filename,'clobber');
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'title',mytitle)
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'history', 'File created using write_FVCOM_TT.m from the MATLAB fvcom-toolbox')
% dimensions
nele_dimid=netcdf.defDim(nc,'nele',nElems);
nsiglay_dimid=netcdf.defDim(nc,'siglay', nSiglay);
% variables and attributes
num_tt_varid=netcdf.defVar(nc,'turbine_number','NC_FLOAT',nele_dimid);
netcdf.putAtt(nc,num_tt_varid,'long_name','Number of turbines in each element');
sigma_layer_varid=netcdf.defVar(nc,'turbine_sigma_layer','NC_FLOAT',[nele_dimid nsiglay_dimid]);
netcdf.putAtt(nc,sigma_layer_varid,'long_name','the fraction of each sigma layer that the turbines are located in');
swept_area_varid=netcdf.defVar(nc,'swept_area','NC_FLOAT',nele_dimid);
netcdf.putAtt(nc,swept_area_varid,'long_name','Area swept by turbine rotor blades');
% variables that could be added in the future
% thrust_coeff_varid=netcdf.defVar(nc,'thrust_coeff','NC_FLOAT',nele_dimid);
% netcdf.putAtt(nc,thrust_coeff_varid,'long_name','Turbine thrust ceofficient');
% blade_coeff_varid=netcdf.defVar(nc,'blade_coeff','NC_FLOAT',nele_dimid);
% netcdf.putAtt(nc,blade_coeff_varid,'long_name','Drag coefficient for the turbine blades');
% end definitions
netcdf.endDef(nc);
% write data
% netCDF variable names:
% turbine_sigma_layer, turbine_number, swept_area
netcdf.putVar(nc,num_tt_varid,turbine.numbers);
netcdf.putVar(nc,sigma_layer_varid,turbine.sigma_layer);
netcdf.putVar(nc,swept_area_varid,turbine.area);
% variables that could be added in the future
% netcdf.putVar(nc,thrust_coeff_varid,turbine.thrust);
% netcdf.putVar(nc,blade_coeff_varid,turbine.drag);
% close file
netcdf.close(nc);
if(ftbverbose);
fprintf(['end : ' subname '\n'])
end;
% Calculate the fraction of the rotor swept area is occupying each sigma layer
%
% Example Usage:
%
% sigma_frac = turbine_area_sigma(H, Ht, r, plot_fig)
%
% Input Parameters: H - mean sea level (m)
% Ht - height of turbine hub above seabed (m)
% r - turbine rotor radius (m)
% plot_fig - flag to plot a figure (optional)
%
% Rory O'Hara Murray, 19-Nov-2014
%
function sigma_frac = turbine_area_sigma(H, Ht, r, plot_fig)
if nargin<4
plot_fig = false;
end
% Turbine and sigma layer parameters
elev = 0; % water elevation above/below MSL - change this to see how the sigma layer occupation fraction changes with the tide
depth = H + elev;
dT = depth - Ht; % turbine hub depth
sigLev = 11;
sigLay = sigLev-1;
dLay = depth./sigLay;
zLev = [0:-dLay:-depth]';
% what sigma layer is the hub in?
drsl = zLev+dT; % depth of hub relative to each sigma level
hub_sigma = sum(drsl>=0);
%% draw sigma levels / layers
if plot_fig
figure
plot([-r r], zLev*[1 1])
xlabel('Distance (m)')
ylabel('Depth (m)')
title([num2str(depth, '%2.0f') ' m water depth'])
% draw rotor area
a=0;
b=-dT;
ang = 0:pi/20:2*pi;
x=r*cos(ang);
y=r*sin(ang);
hold on
plot(a+x,b+y, a, b, 'o')
end
%% what fraction of the rotor area is in each sigma layer?
% loop trough all segments below the hub
dBot=-drsl(-drsl>=0);% the minimum of this array is hub height above a sigma level
numBot=sum(dBot<=r);
segmentsBotCum = [];
for ii=1:numBot
phi = acos(dBot(ii)/r);
sector = phi*r*r;
triBot(ii) = r*sin(phi)*dBot(ii);
segmentsBotCum(ii) = sector-triBot(ii);
end
% loop through all the segments above the hub
dTop=flip(drsl(drsl>=0));
numTop=sum(dTop<=r);
segmentsTopCum = [];
for ii=1:numTop
phi = acos(dTop(ii)/r);
sector = phi*r*r;
triTop(ii) = r*sin(phi)*dTop(ii);
segmentsTopCum(ii) = sector-triTop(ii);
end
segmentsTopCum2 = flip(segmentsTopCum);
segmentsTop = segmentsTopCum2-[0 segmentsTopCum2(1:end-1)];
segmentsBot = segmentsBotCum-[segmentsBotCum(2:end) 0];
% check that there are segments below/above hub, i.e. whether the rotors
% actually span multiple sigma layers
% sig_cent is area of the rotors in the layer the hub is in
if numTop>0 & numBot>0
sigCent = pi*r*r-segmentsTopCum(1)-segmentsBotCum(1);
elseif numTop==0 & numBot>0
sigCent = pi*r*r-segmentsBotCum(1);
elseif numTop>0 & numBot==0
sigCent = pi*r*r-segmentsTopCum(1);
elseif numTop==0 & numBot==0 % entire rotor is confined to one sigma layer
sigCent = pi*r*r;
end
% if sigCent is zero then the hub must exactely co-inside with a sigma
% level (unusual...) check it's larger than a very small area (or zero)
if sigCent>0
% if hub co-insides exactely with a sigma level
segments = [segmentsTop sigCent segmentsBot];
segments_frac = segments./(pi*r*r);
sig_span = hub_sigma + [-length(segmentsTop) length(segmentsBot)];
numSigma = numTop+numBot+1; % total number of occupied sigma layers
else
segments = [segmentsTop segmentsBot];
segments_frac = segments./(pi*r*r);
sig_span = hub_sigma + [-length(segmentsTop) length(segmentsBot)-1];
numSigma = numTop+numBot; % total number of occupied sigma layers
end
% work out the fraction of the turbine area in each sigma layer
sigma_frac = zeros(1,sigLay);
sigma_frac(sig_span(1):sig_span(2)) = segments_frac;
total_frac = sum(sigma_frac);
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment