【物理应用】水下浮动风力涡轮机的尾流诱导动态模拟风场附matlab代码

1 内容介绍

风力发电机的空气动力学性能是决定风力机安全与效率的最重要因素之一.但由于影响风力机气动性能参数众多,更加高效精确地模拟风力机的气动特性一直是风力机研究的重要发展方向.本研究采用浸入边界法对风力机的不同翼型,单级风力机和两级风力机的气动力学进行了一系列的研究.

​2 仿真代码

%% WInDS Driver -> Wake Induced Dynamics Simulator% % Driver script to compute wind turbine performance via unsteady lifting% line method.%% Uses FAST input and output files to define wind turbine geometry and % operating conditions. WInDS then predicts wind turbine performance due % to wake evolution via free vortex wake method and lifting-line theory. %%% ****Function(s)****% constants          Load constants used by other functions% elliptical         Generate geometry and variables for elliptical wing% rotor              Generate geometry and variables for rotor% input_import       Import FAST-formatted input files% output_import      Import FAST-formatted output files% input_mod          Modify inputs, remove discontinuities% kinematics         Compute positions of blade stations% velocity           Compute velocity contributions due to kinematics% initials           Set initial conditions and preallocate memory% performance        Compute performance and load values% %% This work is licensed under the Open Source Initiative BSD 2-Clause % License. To view a copy of this license, visit % http://www.opensource.org/licenses/BSD-2-Clause% %% Written by Thomas Sebastian (tommy.sebastian@gmail.com)% Last edited December 16, 2011%​%% Clear command window and workspaceclear allclose allclc​%% !!!User-defined variables!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!user.t=[0 5 5]; %Initial t, final t, and frequency in Hzuser.filename='NRELrotor'; %Test case (elliptical, rotor type, or .fst file)user.tol=1e-8; %Tolerance value for convergence of numerical methodsuser.d='visc1'; %Core model for filaments (numerical values are the squared cutoff radius,                 %'viscX' applied viscous model of index X) user.co=1000; %Distance from wake nodes beyond which influence is negligibleuser.integ='pcc'; %Numerical integration schemeuser.ns=20; %Number of spanwise stationsuser.maxiter=30; %Maximum number of iterations for Kutta-Joukowski theoremuser.roll='true'; %If 'true', will apply induction to all wake nodesuser.anim='true'; %If 'true', will generate animation of wake evolutionuser.time=datestr(now ,'mm-dd-yyyy_HHMM'); %Date and time of code executionuser.kjtype='fixed'; %Use either fixed point or Brent's method for convergence (Brent is                     %still a bit coarse)user.relax=0.25; %Relaxation value for fixed-point iteration​%%Variables for user.ellip.* used only if user.filename='elliptical'user.ellip.b=10; %Elliptical wingspanuser.ellip.AR=6; %Elliptical wing aspect ratio (AR=b^2/S)user.ellip.wind=[1 0 0]; %Wind velocity vectoruser.ellip.pitch=[5 5 0]; %Pitch angle of elliptical wing (in degrees)user.ellip.pitchrate=0; %Pitch rate of elliptical wing (in degrees)user.ellip.yaw=0; %Yaw angle of elliptical wing (in degrees)​%%Variables for user.rotor.* used only if user.filename='rotor'user.rotor.wind=[11.4 0 0]; %Wind velocity vectoruser.rotor.tsr=7; %Tip speed ratiouser.rotor.casetype='static_rated';user.rotor.pitch=0; %Pitch angle of rotor blade (in degrees)user.rotor.yaw=0;user.rotor.modes=[];%{'Surge' 0.72520 0.00740 -1.16256 -0.44205 0.07750 2.60940 13.60156 10};​addpath(genpath(fullfile(cd))); %Add directories to search path​%% Load constants (physical and derived)[const]=constants;​​%% Load test case (elliptical wing, rotor, or FAST-generated)if strcmp(user.filename,'elliptical')    [blade,turbine,platform,fastout,airfoils,wind]=elliptical(user);elseif strcmp(user.filename,'NRELflat')    [blade,turbine,platform,fastout,airfoils,wind]=NRELflat(user);elseif strcmp(user.filename,'NRELrotor')    [blade,turbine,platform,fastout,airfoils,wind]=NRELrotor(user);elseif strcmp(user.filename,'FAST')    [airfoils,blade,turbine,platform,wind]=input_import(user.filename);    [fastout]=output_import(user.filename,user.t);end​%% Compute positions of blade stations in inertial reference frame[pos]=kinematics(blade,turbine,platform,fastout);​%% Compute velocities of blade stations due to external motions[vel,pos]=velocity(pos,blade,turbine,wind,fastout);​%% Define initial values (wake strength, geometry, etc)[wake,vel,perf]=initials(pos,vel,blade,turbine,wind,airfoils,fastout,const,user);​%% !!!PRIMARY LOOP OVER TIMESERIES!!!%Determine size of test vectors/arraysnt=length(fastout.Time); %Number of timestepsnb=turbine.NumBl; %Number of bladesns=length(blade.RNodes); %Number of shed nodes (stations)tm=zeros(nt,1); %Preallocate memory for timer (time for each timestep)​for p=2:nt    tic; %Begin timing this timestep%Update shed and trailing filament strength    %Bound filament for previous timestep becomes new bound filament    wake.gamma.shed{p}(:,:,1,:)=wake.gamma.shed{p-1}(:,:,1,:);    %Compute spanwise change in bound filament to compute first set of trailing filaments    wake.gamma.trail{p}(:,:,1,:)=diff([zeros(1,1,1,nb) ; wake.gamma.shed{p}(:,:,1,:) ; ...        zeros(1,1,1,nb)],1);    %Previous set of trailing filaments becomes new set of trailing filaments    wake.gamma.trail{p}(:,:,2:end,:)=wake.gamma.trail{p-1};    %Shed filaments computed via spanwise summation of trailing filaments (ensure Kelvin's     %theorem is satisfied)    wake.gamma.shed{p}(:,:,2:end,:)=diff(cat(3,cumsum(wake.gamma.trail{p}(1:end-1,:,:,:),1), ...        zeros(ns,1,1,nb)),1,3);     %Modify vortex core size via Ramasamy-Leishman model and include effect of filament stretching%from previous timestep    wake=vcore(wake,const,fastout,user,p);    %Compute induced velocity at all points    %Velocity induced by shed filaments on all nodes in wake    if strcmp(user.roll,'true')        vel.uind_shed=BiotSavart(wake.domain{p}(1:end-1,:,:,:),wake.domain{p}(2:end,:,:,:), ...            wake.domain{p},wake.gamma.shed{p},wake.rc_eff.shed{p},user.d,user.co,'full');        %Velocity induced by trailing filaments on all nodes in wake        vel.uind_trail=BiotSavart(wake.domain{p}(:,:,2:end,:),wake.domain{p}(:,:,1:end-1,:), ...            wake.domain{p},wake.gamma.trail{p},wake.rc_eff.trail{p},user.d,user.co,'full');        %Sum the induced velocity contributions due to shed and trailing filaments        vel.uind{p}=vel.uind_shed+vel.uind_trail;    end    %Add the total induced velocity in the wake to the freestream velocity    vel.domain{p}=vel.domain{p}+vel.uind{p};        %Numerically convect wake nodes to time+1    if strcmp(user.integ,'fe') && p~=nt        wake=fe(wake,vel,user,p); %Foward euler    elseif strcmp(user.integ,'ab2') && p~=nt        wake=ab2(wake,vel,user,p); %2nd-order Adams-Bashforth    elseif strcmp(user.integ,'ab4') && p~=nt        wake=ab4(wake,vel,user,p); %2nd-order Adams-Bashforth    elseif strcmp(user.integ,'pcc') && p~=nt        wake=pcc(wake,vel,const,fastout,user,p); %Predictor-corrector, central-difference    end​%Compute strength of new bound vortex via Kutta-Joukowski theorem    [wake,perf,vel,ctj]=KuttaJoukowski(pos,vel,blade,turbine,wake,airfoils,user,perf,p, ...        user.kjtype);    %Determine time spent on current timeloop and estimate time remaining    tm(p-1)=toc; %Time spent on current loop    if p>2        pt=polyfit([0 ; (2:p)'],cumsum([0 ; tm(1:p-1)]),2);        tr=polyval(pt,nt)-sum(tm(1:p-1)); %Extrapolate to determine time remaining        clc; disp([num2str(ctj) ': ' num2str(p/nt*100) ...            '% complete, estimated time remaining: ' num2str(tr/60) ' minutes'])    endend​%% Compute performance metricsperform;​%% Tidy up the workspaceclear yn j nb nt wb1 vs vt pg nst ns trsave(['savedsims\' user.time '_' user.filename '_' user.rotor.casetype '.mat'])​%% Generate wake figureif strcmp(user.anim,'true')    j=length(fastout.Time);    wakeplot(pos,vel,turbine,blade,wake,fastout,j);end

3 运行结果

4 参考文献

1] T. Sebastian and M. Lackner. Development of a Free Vortex Wake Model Code for Offffshore Floating

Wind Turbines. Renewable Energy, Online:1–15, 2011.

[2] T. Sebastian and M. Lackner. Unsteady Aerodynamics of Offffshore Floating Wind Turbines.

Wind

Energy, Online:1–14, 2011. doi: 10.1002/we.545.

[3] Sheila E. Widnall. The Structure and Dynamics of Vortex Filaments. Annual Review of Fluid Mechanics,

7:141–165, 1975.

[4] Mahendra J. Bhagwat and J. Gordon Leishman. Stability, Consistency and Convergence of Time

Marching Free-Vortex Rotor Wake Algorithms. Journal of the American Helicopter Society, 46(1):59–71,

January 2001.

[5] J. Gordon Leishman. Principles of Helicopter Aerodynamics (Cambridge Aerospace Series). Cambridge

University Press, 2006. ISBN 0521858607.

[6] Jack B. Kuipers. Quaternions and Rotation Sequences. Princeton University Press, 1998.

博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。

部分理论引用网络文献，若有侵权联系博主删除。