💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。
⛳️座右铭：行百里者，半于九十。
📋📋📋本文目录如下：🎁🎁🎁
目录
⛳️赠与读者
👨💻做科研，涉及到一个深在的思想系统，需要科研者逻辑缜密，踏实认真，但是不能只是努力，很多时候借力比努力更重要，然后还要有仰望星空的创新点和启发点。当哲学课上老师问你什么是科学，什么是电的时候，不要觉得这些问题搞笑。哲学是科学之母，哲学就是追究终极问题，寻找那些不言自明只有小孩子会问的但是你却回答不出来的问题。建议读者按目录次序逐一浏览，免得骤然跌入幽暗的迷宫找不到来时的路，它不足为你揭示全部问题的答案，但若能让人胸中升起一朵朵疑云，也未尝不会酿成晚霞斑斓的别一番景致，万一它居然给你带来了一场精神世界的苦雨，那就借机洗刷一下原来存放在那儿的“躺平”上的尘埃吧。
或许，雨过云收，神驰的天地更清朗.......🔎🔎🔎
💥1 概述
唤醒诱导动力学模拟器（WInDS）是一个基于提线理论（LLT）的自由涡尾流方法（FVM）代码，由马萨诸塞大学阿默斯特分校风能中心开发，其明确目的是以比动量平衡方法更高的精确度来模拟海上漂浮风力发电机（OFWT）的气动力学。WInDS本身集成了海上漂浮风力发电机中存在的多个自由度，从而实现了对流场更真实的模拟。
唤醒诱导动力学模拟器（WInDS），由Sebastian和Lackner[1]编写，是一个基于提线理论（LLT）的自由涡尾流方法（FVM）代码，该代码在马萨诸塞大学阿默斯特分校风能中心开发，其明确目标是为了模拟海上漂浮风力发电机（OFWT）的气动力学，达到比动量平衡方法所能达到的更高精度。WInDS本身就融入了海上漂浮风力发电机中存在的多自由度，这使得流场的模拟更加逼真。
Sebastian和Lackner[2]证明，与海上漂浮风力发电机平台运动（图1）相关的额外自由度将导致气动不稳定性，这种不稳定性超过了陆上和传统海上系统的经历的不稳定性。
势流理论允许叠加基本流以构建更复杂的流体运动。涡丝是三维基本流的一个例子，是集中涡度的物质线。Helmholtz定理指出，涡丝上的环量强度或涡度沿涡丝是恒定的，涡丝必须形成一个闭合的环或者延伸至无穷远。多个涡丝可以组合起来形成一个随每个时间步长增长的闭合涡格，以此来模拟与尾流相关的复杂且不稳定的流场。由强度为Γ、长度为L、位于节点x1和x2之间的涡丝段（图2）在感兴趣的点P处诱导的速度，可以通过Biot-Savart定律进行计算。详细文章见第4部分。
📚2 运行结果
部分代码：
user.t=[0 5 5]; %Initial t, final t, and frequency in Hz
user.filename='NRELrotor'; %Test case (elliptical, rotor type, or .fst file)
user.tol=1e-8; %Tolerance value for convergence of numerical methods
user.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 negligible
user.integ='pcc'; %Numerical integration scheme
user.ns=20; %Number of spanwise stations
user.maxiter=30; %Maximum number of iterations for Kutta-Joukowski theorem
user.roll='true'; %If 'true', will apply induction to all wake nodes
user.anim='true'; %If 'true', will generate animation of wake evolution
user.time=datestr(now ,'mm-dd-yyyy_HHMM'); %Date and time of code execution
user.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 wingspan
user.ellip.AR=6; %Elliptical wing aspect ratio (AR=b^2/S)
user.ellip.wind=[1 0 0]; %Wind velocity vector
user.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 vector
user.rotor.tsr=7; %Tip speed ratio
user.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/arrays
nt=length(fastout.Time); %Number of timesteps
nb=turbine.NumBl; %Number of blades
ns=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'])
end
end
🎉3 参考文献
文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。
🌈4 Matlab代码、文章
资料获取，更多粉丝福利，MATLAB|Simulink|Python资源获取
385
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
