[physical application] Wake induced dynamic simulation of underwater floating wind turbine wind field with matlab code

1 Content introduction

The aerodynamic performance of wind turbine is one of the most important factors that determine the safety and efficiency of wind turbine . However, there are many parameters affecting the aerodynamic performance of wind turbines , More efficient and accurate simulation of the aerodynamic characteristics of wind turbines has always been an important development direction of Wind Turbine Research . In this study, the immersed boundary method is used to analyze different airfoils of wind turbines , A series of studies have been carried out on the aerodynamics of single-stage wind turbines and two-stage wind turbines .

​2 Simulation code

%% 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 ([email protected])% 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 Running results

4 reference

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.

About bloggers ： Good at intelligent optimization algorithms 、 Neural networks predict 、 signal processing 、 Cellular automata 、 The image processing 、 Path planning 、 UAV and other fields Matlab Simulation , relevant matlab Code problems can be exchanged by private letter .

Some theories cite network literature , If there is infringement, contact the blogger to delete .