Contents

classdef JensenWakeModel < handle
    properties
        c
    end

    properties (SetAccess = protected)
        A
        k
    end

    methods
        function O = JensenWakeModel()
            O.setConfig(JensenWakeModelConfig());
        end

        function O = setConfig(O, config)
            O.c     = config;
            O.k     = O.getKFromTI(O.c.TI);
        end

        function k = getKFromTI(O, TI)
            TIs = [.08  .10  .13  .15  .16  .18  .21  .24  .29];
            ks  = [.040 .052 .063 .075 .083 .920 .100 .108 .117];

            k = interp1(TIs, ks, TI);
        end

        function O = calculateFlowField(O, M)
            for ii = 1:M.layout.N
                A(:,:,ii) = O.calculateSingleWake( ...
                    M, ...
                    M.layout.z(ii) , ...
                    M.layout.y(ii), ...
                    1/3 ...
                );
            end

            O.A = O.mixWakes(A);
        end

        function A = calculateSingleWake(O, M, z0, y0, a)

Fetch Information

            k  = O.k;
            R  = M.turbineType.R;
            Y  = M.mesh.Y;
            Z  = M.mesh.Z;
            nY = M.mesh.nY;

Determine Domain

            theta       = atan2((Y-y0),(Z-z0));

            distTotal   = sqrt((Z - z0).^2 + (Y-y0).^2);
            distCentral = distTotal .* sin(theta - M.c.alpha);
            distCalc    = distTotal .* cos(theta - M.c.alpha);

            r    = R + k * distCalc;

            mask = abs(distCentral) < r & distCalc > 0;

Calculate Wake

            K  = mask .* R ./ r;

            A = 1 - (1-(1-2*a)) .* K.^2;

        end

        function A = mixWakes(O, As)

Multiply

            sizes = size(As);
            A     = ones(sizes(1), sizes(2));

            for ii = 1:sizes(3)
                A = A .* As(:,:,ii);
            end

Maximum Wake

           if(sizes(3) > 1)
               for ii = 1:sizes(3)
                   A = min(A, As(:,:,ii));
               end
           else
               A = As(:,:,1);
           end

        end

        function h = plotFlowField(O, M)
            h = surf(M.mesh.Z/M.turbineType.D,M.mesh.Y/M.turbineType.D, O.A);
            shading interp;
            xlabel('z [#D]');
            ylabel('y [#D]');
            zlabel('u [m/s]');
            view(0,90);
            axis tight equal;
            colorbar;
            set(gca, 'CLim', [0, 1]);
            zlim([0 1]);
        end

        function A = getFlowField(O)
            A = O.A;
        end
    end

end

Published with MATLAB® R2014b