Skip to content
Snippets Groups Projects
Commit 115875a6 authored by RaspberryPi's avatar RaspberryPi
Browse files

Merge branch 'master' of https://gitlab.com/kakila/windturbine_HSR

parents 8f5ee667 2d04eb18
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
CFD/DMSM/dmsm.m
+ 103
− 7
View file @ 115875a6
......@@ -45,23 +45,27 @@ function [Torque, Cp, varargin] = dmsm (lambda, param)
rthetaUP = [cos(param.thetaUP(:)) sin(param.thetaUP(:))];
% theta interval for integration
theta_intDN = {0.5 * pi + 1e-6, 1.5 *pi - 1e-6};
theta_intUP = {-0.5 * pi + 1e-6, 0.5 * pi - 1e-6};
theta_intDN = {pi/2 + param.iteration.int_eps, ...
3*pi/2 - param.iteration.int_eps};
theta_intUP = {-pi/2 + param.iteration.int_eps, ...
pi/2 - param.iteration.int_eps};
%interference factors
u_factorUP = u_factorDN = zeros (param.NS,1);
% Inputs for iteration solver.
% Empties are placeholders
args = {[], [], [], ...
v_tang, scaleForce, scaleG, delta, c_over_eta, param.iteration};
args = {[], [], [], []...
v_tang, scaleForce, delta, c_over_eta, param.iteration};
argsUP = args;
argsUP{1} = param.v_inf;
argsUP{2} = theta_intUP;
argsUP{4} = scaleG;
argsDN = args;
argsDN{2} = theta_intDN;
argsDN{4} = -scaleG;
% Pre-allocation of result arrays
FtUP = u_factorUP = zeros (1, param.NS);
......@@ -137,7 +141,7 @@ end%function
%! % With u restart (default)
%! T = Cp = zeros (nl, 1);
%! tic
%! for i = 1:nl;
%! for i = 1:nl
%! [T(i) Cp(i)] = dmsm (l(i), p);
%! end%for
%! toc
......@@ -146,7 +150,7 @@ end%function
%! p.iteration.reset_u = false;
%! Tr = Cpr = zeros (nl, 1);
%! tic
%! for i = 1:nl;
%! for i = 1:nl
%! [Tr(i) Cpr(i)] = dmsm (l(i), p);
%! end%for
%! toc
......@@ -166,3 +170,95 @@ end%function
%! ylabel ('Cp [-]');
%! grid on
%! legend (h(1), {'reseting u'})
%! xlabel ('Tip/Wind speed ratio')
%!demo
%! warning ('This demo takes about 1 minute to run')
%! p = dmsm_parameters ();
%! l = 4.5;
%! % reduce integration interval of divergent integral by this amount
%! n = 10;
%! int_eps = logspace(-10,-1, n);
%!
%! % With u restart (default)
%! T = Cp = zeros (n, 1);
%! tic
%! for i = 1:n
%! p.iteration.int_eps = int_eps(i);
%! [T(i) Cp(i)] = dmsm (l, p);
%! end%for
%! toc
%!
%! % Without u restart
%! p.iteration.reset_u = false;
%! Tr = Cpr = zeros (n, 1);
%! tic
%! for i = 1:n
%! p.iteration.int_eps = int_eps(i);
%! [Tr(i) Cpr(i)] = dmsm (l, p);
%! end%for
%! toc
%!
%! figure (1)
%! subplot (2,1,1);
%! h = semilogx (int_eps,[T Tr], 'o');
%! for i=1:2; set (h(i), 'markerfacecolor', get (h(i), 'color')); end
%! axis tight
%! ylabel ('Torque [Nm]');
%! grid on
%! legend (h(1), {'reseting u'})
%! subplot (2,1,2);
%! h = semilogx (int_eps,[Cp Cpr], 'o');
%! for i=1:2; set (h(i), 'markerfacecolor', get (h(i), 'color')); end
%! axis tight
%! ylabel ('Cp [-]');
%! grid on
%! legend (h(1), {'reseting u'})
%! xlabel ('interval eps')
%!demo
%! warning ('This demo takes about 1 minute to run')
%! p = dmsm_parameters ();
%! p.iteration.int_eps = 0;
%! nl = 30;
%! l = linspace (0.25, 7, nl);
%!
%! % With u restart (default)
%! T = Cp = zeros (nl, 1);
%! tic
%! for i = 1:nl
%! [T(i) Cp(i)] = dmsm (l(i), p);
%! end%for
%! toc
%!
%! % Without u restart
%! p.iteration.reset_u = false;
%! Tr = Cpr = zeros (nl, 1);
%! tic
%! for i = 1:nl
%! [Tr(i) Cpr(i)] = dmsm (l(i), p);
%! end%for
%! toc
%!
%! figure (1)
%! subplot (2,1,1);
%! h = plot (l,[T Tr], 'o-');
%! for i=1:2; set (h(i), 'markerfacecolor', get (h(i), 'color')); end
%! axis tight;
%! ylabel ('Torque [Nm]');
%! grid on
%! legend (h(1), {'reseting u'})
%! subplot (2,1,2);
%! h = plot (l,[Cp Cpr], 'o-');
%! for i=1:2; set (h(i), 'markerfacecolor', get (h(i), 'color')); end
%! axis tight;
%! ylabel ('Cp [-]');
%! grid on
%! legend (h(1), {'reseting u'})
%! xlabel ('Tip/Wind speed ratio')
%!test
%! p = dmsm_parameters();
%! [t c] = dmsm(0.45, p);
%! printf ('Torque: %g\n', t);
%! printf ('Cp : %g\n', c);
CFD/DMSM/dmsm_iteration.m
+ 25
− 16
View file @ 115875a6
......@@ -17,9 +17,9 @@
%% Author: Matthias Pasquon <matthias.pasquon@hsr.ch> (26.06.2016)
%% Adapted by: Juan Pablo Carbajal <ajuanpi+dev@gmail.com>
function [FT fctor varargout] = dmsm_iteration (lambda, v_ref, theta_interval, ...
rtheta, vt, scaleF, scaleG, delta, ...
c_over_eta, iter_params)
function [FT fctor varargout] = ...
dmsm_iteration (lambda, v_ref, theta_interval, rtheta, scaleG, ...
vt, scaleF, delta, c_over_eta, iter_params)
% DMSM_ITERATION
% Function to calculate the iteration for the induction factor fctor and the
% corresponding aerodynamical forces on the blade in the streamtube
......@@ -32,6 +32,7 @@ function [FT fctor varargout] = dmsm_iteration (lambda, v_ref, theta_interval, .
% See also dmsm
% interference factor
fctor = iter_params.u_factor;
tol = iter_params.AbsTol;
......@@ -69,28 +70,36 @@ function [FT fctor varargout] = dmsm_iteration (lambda, v_ref, theta_interval, .
CN = CL * car + CD * sar;
CT = CL * sar - CD * car;
% fdn: function to find interference factor u (eq. 3.13)
g = @(x) abs (sec (x)) .* (CN .* cos (x) - CT .* sin (x));
% Function to find interference factor u (eq. 3.13)
%g = @(x) abs (sec (x)) .* (CN .* cos (x) - CT .* sin (x)) .* ...
% (TSR * (TSR - 2 * sin (x)) + 1);
%G = abs(scaleG) * quadgk (g, theta_interval{:});
% g can be written as
% g = sign (cos(x)) * ( CN - CT * tan (x) )
% But
% sign(cos(x)) == 1 for x in I+ = [0, pi/2) U (3*pi/2,2*pi]
% sign(cos(x)) == -1 for x in I- = (pi/2, 3*pi/2)
% The integral has analytic solution in those intervals
vrel2 = v_relbl^2;
G = scaleG * (vrel2 / v_blade^2) * quadgk (g, theta_interval{:});
% g = sign (cos(x)) * ( CN - CT * tan (x) ) * X^2
% with X^2 = TSR * (TSR - 2 * sin (x)) + 1
% and needs to be integrated between [-pi/2,pi/2]
% sign(cos(x)) == 1 for this interval
% The integral has then an analytic solution
% but it is divergent!
% https://gitlab.com/kakila/windturbine_HSR/wikis/Upwind-Downwind-function-(divergent)
% using that page we separate the contributions from CN and CT
GN = CN * ( pi * ( TSR.^2 + 0.5) + 1);
gT = @(x) (TSR * (TSR - 2 * sin (x)) + 1) .* tan(x);
GT = CT * quadgk (gT, theta_interval{:});
G = scaleG * (GN - GT);
% interference factor for next iteration process (eq. 3.14)
fctor_next = pi / ( G + pi );
if (fctor - fctor_next < tol)
break
else
fctor = fctor_next;
end%if
end%while
% Normal & Tangential force (eq. 2.6-2.7)
vrel2 = v_relbl^2;
FN = scaleF * CN * vrel2;
FT = scaleF * CT * vrel2;
......
CFD/DMSM/dmsm_parameters.m
+ 1
− 1
View file @ 115875a6
......@@ -64,7 +64,7 @@ function param = dmsm_parameters()
param.iteration.u_factor = 1; % Initial vlaue for velocity induction factor
param.iteration.AbsTol = 1e-6; % Stoping tolerance for iterations
param.iteration.reset_u = true;% whether u_factor is used in next iteration
param.iteration.int_eps = 1e-5; % hack to avoid singularity in the integral for the interaction factor
% Function to interpolate Lift and Drag coefficients
% TODO interace with table generator, like xfoil
NACA_profile = 'NACA0012'; % NACA-profile
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment