风电塔架风速时程曲线_风荷载时程曲线模拟_matlab仿真_matlab

%代码"Wind"由University of Waterloo的Chad Van der Woude提供 % %代码由Jiawei修改并扩写，以用于5MW风机混凝土塔架的风荷载模拟 % % 本代码依据IEC 61400-1标准，用于模拟具有一定功率谱密度和自相关性的风速时程曲线 % 方法为基于三角级数叠加的谐波合成法（Hansen，2008） flow = 1/10; % 频率下限 (Hz) fhigh = 40; % 频率上限 (Hz) rad = 12; % 模拟叶片上风荷载的位置数 num = 10; % 每个叶片划分段数 fL = 63; % 包括转子半径的叶片长度 hL = 3; % 转子半径 HH = 100; % 转子中心高度 numT = 10; % 塔架分段数 alpha = 0.2; % 幂率分布的指数 Trev = 3.5; % 风机自转周期 L = fL-hL; % 不计转子半径的叶片长度 dL = L/num; % 每个叶片段长度 ang = 2*pi/rad; % 每两个叶片计算位置间的夹角 T = 1/flow; % 风荷载模拟时间长度（s） Tr =T/5 ; % 风速上升到最大的时间 N = fhigh*2*T; % 风荷载模拟时刻数 nF = T*fhigh; % 风荷载模拟频率数 dt = T/N; % 时间微元 dF = 1/T; % 频率微元 % 产生时间序列向量 time = zeros(N,1); for c1 = 1:N time(c1,1) = c1*dt; end % 规定风场特性 V = 10; % 平均风速 (m/s) I = 0.18; % 湍流强度 sd = I*V; % 脉动风速标准差 (m/s) Ls = 340; % 湍流整体尺度 (m)，按 IEC 61400-1 取值。 npts = rad*num+numT; pts = zeros(npts,2); % 计算风荷载作用点坐标，以确定不同点间互相干函数 display('计算风荷载作用点坐标') for c1 = 1:num; for c2 = 1:rad; numpt = c2+(c1-1)*rad; az = (c2-1)*ang; pts(numpt,1) = (hL + (c1-0.5)*dL)*cos(pi/2-az); pts(numpt,2) = HH+(hL + (c1-0.5)*dL)*sin(pi/2-az); end end for c1 = 1:numT numpt = rad*num+c1; pts(numpt,1) = 0; pts(numpt,2) = (HH/numT)*(c1-0.5); end display('计算谱密度函数矩阵S') % 计算谱密度函数矩阵S S = zeros(npts,npts,nF); Ssum = zeros(npts,npts); for c1 = 1:npts for c2 = 1:npts for c3 = 1:nF sep = ((pts(c1,1)-pts(c2,1))^2+(pts(c1,2)-pts(c2,2))^2)^0.5; freq = c3/T; S(c1,c2,c3) = exp(-12*((freq*sep/V)^2+(0.12*sep/Ls)^2)^0.5)*(4*sd^2*Ls/V)/(1+6*freq*Ls/V)^(5/3); end end end display('将S矩阵分解为下三角矩阵H与其共轭矩阵的转置的积') % 递归算法将S矩阵分解为下三角矩阵H与其共轭矩阵的转置的积 H = zeros(npts,npts,nF); for c1 = 1:nF H(1,1,c1) = sqrt(S(1,1,c1)); H(2,1,c1) = S(2,1,c1)/H(1,1,c1); H(2,2,c1) = (S(2,2,c1)-H(2,1,c1)^2)^0.5; for c2 = 3:npts for c3 = 1:c2 if c3 ~= 1 sum1 = 0; sum2 = 0; for c4 = 1:(c3-1) sum1 = sum1 + H(c3,c4,c1)^2; sum2 = sum2 + H(c2,c4,c1)*H(c3,c4,c1); end end if c3 == 1 H(c2,c3,c1) = S(c2,c3,c1)/H(c3,c3,c1); elseif c2 == c3 H(c2,c3,c1) = (S(c2,c3,c1) - sum1)^0.5; else H(c2,c3,c1) = (S(c2,c3,c1)-sum2)/H(c3,c3,c1); end end end end display('生成随机相位角序列') % 生成随机相位角序列 phi = zeros(npts,nF); for c1 = 1:npts for c2 = 1:nF phi(c1,c2) = rand*2*pi; end end display('生成风场模拟所需向量') % 生成 "Real", "Imag", "Amp", "Pha" 向量 Real = zeros(npts,nF); Imag = zeros(npts,nF); Amp = zeros(npts,nF); Pha = zeros(npts,nF); for c1 = 1:npts for c2 = 1:nF for c3 = 1:c1 Real(c1,c2) = Real(c1,c2) + H(c1,c3,c2)*cos(phi(c3,c2)); Imag(c1,c2) = Imag(c1,c2) + H(c1,c3,c2)*sin(phi(c3,c2)); Amp(c1,c2) = sqrt(Real(c1,c2)^2 + Imag(c1,c2)^2); if Imag(c1,c2) > 0 if Real(c1,c2) > 0 Pha(c1,c2) = atan(Imag(c1,c2)/Real(c1,c2)); else Pha(c1,c2) = pi - atan(abs(Imag(c1,c2)/Real(c1,c2))); end else if Real(c1,c2) > 0 Pha(c1,c2) = 2*pi - atan(abs(Imag(c1,c2)/Real(c1,c2))); else Pha(c1,c2) = pi + atan(abs(Imag(c1,c2)/Real(c1,c2))); end end end end end display('生成风速变化时程') % 生成风速变化时程 u = zeros(npts,N); %u的结构：每行表示一个点，每列表示一个时刻 for c1 = 1:npts for c2 = 1:N for c3 = 1:nF u(c1,c2) = u(c1,c2) + 2*Amp(c1,c3)*cos(2*pi*(c3/T)*(c2*dt)-Pha(c1,c3)); end end end display('归一化风速变化时程') % 归一化风速变化时程，生成总时程 for c1 = 1:npts sdbefore = std(u(c1,:)); u(c1,:) = (sd/sdbefore)*u(c1,:); for c2 = 1:N u(c1,c2) = u(c1,c2) + V; end end % 整理数据 display('静态风速时程') uSb = zeros(3*num*N,3); uSt = zeros(numT*N,3); display('叶片1静态风速时程') for c1 = 1:rad:rad*num-2/3*rad newC = floor(c1/rad)*3+1; uSb(1+(newC-1)*N:newC*N,1) = ones(N,1)*newC; uSb(1+(newC-1)*N:newC*N,2) = time; uSb(1+(newC-1)*N:newC*N,3) = u(c1,:); end display('叶片2静态风速时程') for c1 = 1+rad/3:rad:rad*num-1/3*rad newC = floor(c1/rad)*3+2; uSb(1+(newC-1)*N:newC*N,1) = ones(N,1)*newC; uSb(1+(newC-1)*N:newC*N,2) = time; uSb(1+(newC-1)*N:newC*N,3) = u(c1,:); end display('叶片3静态风速时程') for c1 = 1+2*rad/3:rad:rad*num newC = floor((c1-1)/rad)*3+3; uSb(1+(newC-1)*N:newC*N,1) = ones(N,1)*newC; uSb(1+(newC-1)*N:newC*N,2) = time; uSb(1+(newC-1)*N:newC*N,3) = u(c1,:); end display('塔架静态风速时程') for c1 = num*rad+1:num*rad+numT uSt(1+(c1-num*rad-1)*N:(c1-num*rad)*N,1) = ones(N,1)*(c1-num*rad); uSt(1+(c1-num*rad-1)*N:(c1-num*rad)*N,2) = time; uSt(1+(c1-num*rad-1)*N:(c1-num*rad)*N,3) = u(c1,:); end % 旋转叶片动态风速时程 display('旋转叶片动态风速时程') sector = zeros(rad+1,1); for c1 = 1:rad sector(c1,1) = (c1-1)*Trev/rad; sector(c1,2) = c1; end sector(rad+1,1) = Trev; sector(rad+1,2) = 1; uRb = zeros(3*num*N,3); for c1 = 1:N [mindiff closest] = min(abs(sector(:,1) - mod(time(c1),Trev))); count = 0; for c2 = 1:rad:rad*num-2/3*rad newC = floor(c2/rad)*3+1; cr = floor(c2/rad)*rad+sector(closest,2); uRb(c1+(newC-1)*N,1) = newC; uRb(c1+(newC-1)*N,2) = time(c1); uRb(c1+(newC-1)*N,3) = u(cr,c1); end [mindiff closest] = min(abs(sector(:,1) - mod(time(c1)+1/3*Trev,Trev))); for c2 = 1+rad/3:rad:rad*num-1/3*rad newC = floor(c2/rad)*3+2; cr = floor(c2/rad)*rad+sector(closest,2); uRb(c1+(newC-1)*N,1) = newC; uRb(c1+(newC-1)*N,2) = time(c1); uRb(c1+(newC-1)*N,3) = u(cr,c1); end [mindiff closest] = min(abs(sector(:,1) - mod(time(c1)+2/3*Trev,Trev))); for c2 = 1+2*rad/3:rad:rad*num newC = floor(c2/rad)*3+3; cr = floor(c2/rad)*rad+sector(closest,2); uRb(c1+(newC-1)*N,1) = newC; uRb(c1+(newC-1)*N,2) = time(c1); uRb(c1+(newC-1)*N,3) = u(cr,c1); end end display('塔架风速时程') uRt = uSt; display('考虑风速上升时间') for c1 = 1:3*num*N if uSb(c1,2) <= Tr uSb(c1,3) = uSb(c1,2)/Tr*uSb(c1,3); end if uRb(c1,2) <= Tr uRb(c1,3) = uRb(c1,2)/Tr*uRb(c1,3); end end for c1 = 1:numT*N if uSt(c1,2) <= Tr uSt(c1,3) = uSt(c1,2)/Tr*uSt(c1,3); end if uRt(c1,2) <= Tr uRt(c1,3) = uRt(c1,2)/Tr*uRt(c1,3); end end display('考虑风速沿高度分布') for c1 = 1:numT*N uSt(c1,3) = uSt(c1,3)*(pts(rad*num+uSt(c1,1),2)/HH)^alpha; uRt(c1,3) = uRt(c1,3)*(pts(rad*num+uRt(c1,1),2)/HH)^alpha; end display('计算风荷载') % 用叶素动量法计算风荷载 % 叶片风荷载 c=linspace(6,3,11); % 叶片计算宽度（弦长,m） rho=1.25; % 空气密度(kg/m^3) theta=linspace(0.2,0,11); % 叶片扭转角（rad） r=linspace(6,60,10); % 叶片位置（到转轴距离,m） w=2*pi/Trev; % 转速（rad/s） vr=w*r; % 自转引起的速度 % 叶片1线荷载计算 PN1=zeros(N,10); vref=zeros(N,10); phi=vref; % 合速度与旋转平面夹角（rad），随时间变化 AA=phi; % 攻角（rad），随时间变化 CN=AA; % 沿风速方向风力系数 for sec=1:10 for t=1:N vref(t,sec)=sqrt((vr(sec))^2+(uRb(t+N*(sec-1)*3,3))^2); % 1,4,7... phi(t,sec)=acos(vr(sec)/vref(t,sec)); AA(t,sec)