turbine_area_sigma.m 3.91 KB
Newer Older
1 2 3 4
% Calculate the fraction of the rotor swept area is occupying each sigma layer
%
% Example Usage:
%
5
% sigma_frac = turbine_area_sigma(H, Ht, r, sigLay, plot_fig, subplot_info)
6 7 8 9
%
% Input Parameters:     H  - mean sea level (m)
%                       Ht - height of turbine hub above seabed (m)
%                       r  - turbine rotor radius (m)
10
%                       sigLay - number of sigma layers (not levels) in the model
11 12 13 14 15
%                       plot_fig - optional; flag to plot a figure
%                       subplot_info - optional; if present should be an
%                       array containing the three parameters to subplot
%                       that should be used to put the figure into a
%                       subplot.
16 17
%
% Rory O'Hara Murray, 19-Nov-2014
18
% Simon Waldman, 2016.
19
%
20
function sigma_frac = turbine_area_sigma(H, Ht, r, sigLay, plot_fig, subplot_info)
21

22 23 24 25
assert(nargin >= 4, 'Not enough arguments.');
assert(isnumeric(sigLay) && sigLay - fix(sigLay) < eps, 'sigLay (4th parameter) must be an integer number of sigma layers.');

if nargin<5
26 27
    plot_fig = false;
end
28 29 30 31 32
if nargin<6
    splot = false;
else
    splot = true;
end
33

34
dT = H - Ht; % turbine hub depth
35

36 37
assert(dT>r, 'Turbine will stick out of water');

38 39
dLay = H./sigLay;
zLev = [0:-dLay:-H]';
40 41 42 43 44 45 46

% 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
47 48 49 50 51
    if splot
        subplot( subplot_info(1), subplot_info(2), subplot_info(3) )
    else
        figure
    end
52 53 54
    plot([-r r], zLev*[1 1])
    xlabel('Distance (m)')
    ylabel('Depth (m)')
55
    title([num2str(H, '%2.0f') ' m water depth'])
56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
    
    % 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);