基于MATLAB实现的AR法模拟脉动风场风速+使用说明文档.rar

%% Programmed by Duan Zhiyong; %- AR模型模拟脉动风速； v10 = 20; n = 0.01:0.01:10; xn = 1200*n./v10; k = 0.003; ti = 0.1; s1 = 4*k*v10^2*xn.^2./n./(1+xn.^2).^(4/3); %- Davenport谱 %- 定义空间点坐标 for i = 1:20 x(i) = 5+i; z(i) = 5+i; end %- AR模型取4阶，p=4 %- 求R矩阵，分别求出R0,R1,R2,R3,R4 syms f R0 = zeros(20); for i = 1:20 for j = i:20 H0 = inline('(4*k*v10^2*(1200*f/v10).^2)./f./(1+(1200*f/v10).^2).^(4/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2))','f','k','dx','dz','v10'); dx = x(i)-x(j); dz = z(i)-z(j); R0(i,j) = quadl(H0,0.001,10,0.001,0,k,dx,dz,v10); R0(j,i) = R0(i,j); end end R1 = zeros(20); for i = 1:20 for j = i:20 H1 = inline('(4*k*v10^2*(1200*f/v10).^2)./f./(1+(1200*f/v10).^2).^(4/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*1*ti)','f','k','dx','dz','ti','v10'); dx = x(i)-x(j); dz = z(i)-z(j); R1(i,j) = quadl(H1,0.001,10,0.001,0,k,dx,dz,ti,v10); R1(j,i) = R1(i,j); end end R2 = zeros(20); for i = 1:20 for j = i:20 H2 = inline('(4*k*v10^2*(1200*f/v10).^2)./f./(1+(1200*f/v10).^2).^(4/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*2*ti)','f','k','dx','dz','ti','v10'); dx = x(i)-x(j); dz = z(i)-z(j); R2(i,j) = quadl(H2,0.001,10,0.001,0,k,dx,dz,ti,v10); R2(j,i) = R2(i,j); end end R3 = zeros(20); for i = 1:20 for j = i:20 H3 = inline('(4*k*v10^2*(1200*f/v10).^2)./f./(1+(1200*f/v10).^2).^(4/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*3*ti)','f','k','dx','dz','ti','v10'); dx = x(i)-x(j); dz = z(i)-z(j); R3(i,j) = quadl(H3,0.001,10,0.001,0,k,dx,dz,ti,v10); R3(j,i) = R3(i,j); end end R4 = zeros(20); for i = 1:20 for j = i:20 H4 = inline('(4*k*v10^2*(1200*f/v10).^2)./f./(1+(1200*f/v10).^2).^(4/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*4*ti)','f','k','dx','dz','ti','v10'); dx = x(i)-x(j); dz = z(i)-z(j); R4(i,j) = quadl(H4,0.001,10,0.001,0,k,dx,dz,ti,v10); R4(j,i) = R4(i,j); end end A = [R0 R1 R2 R3 ;R1 R2 R3 R0; R2 R3 R0 R1;R3 R0 R1 R2]; B = [R1;R2;R3;R4]; F = A\B; %- 80*20矩阵 q1 = F(1:20,:); %- 取F的第一个20*20矩阵 q2 = F(20+1:2*20,:); %- 取F的第二个20*20矩阵 q3 = F(2*20+1:3*20,:); %- 取F的第三个20*20矩阵 q4 = F(3*20+1:4*20,:); %- 取F的第四个20*20矩阵 Q1 = q1'; Q2 = q2'; Q3 = q3'; Q4 = q4'; RN = R0-(Q1*R1+Q2*R2+Q3*R3+Q4*R4); L = chol(RN); %- cholesky分解 nt = zeros(20,2048); for i = 1:2048 nt(:,i) = normrnd(0,1,20,1); end V(:,1) = L*nt(:,1); %- t1时各点的风速 V(:,2) = Q1*V(1:20,1)+L*nt(:,2); %- t2时各点的风速 V(:,3) = (Q1*V(1:20,2)+Q2*V(1:20,1))+L*nt(:,3); %- t3时各点的风速 V(:,4) = (Q1*V(1:20,3)+Q2*V(1:20,2)+Q3*V(1:20,1))+L*nt(:,4); %- t4时各点的风速 for t = 5:2048 V(:,t) = (Q1*V(1:20,t-1)+Q2*V(1:20,t-2)+Q3*V(1:20,t-3)+Q4*V(1:20,t-4))+L*nt(:,t); %- t时各点的风速 end %- 求某一点的风速 V1 = V(10,:); %- 取第10点的风速 t = (1:2048)*ti; figure; subplot(2,1,1); plot(t,V1,'k-'); %- 第10点的风速时程图 xlabel('t/s'); ylabel('v/(ms-1)'); axis([0 200 -8 8]); %- 与目标谱进行比较 nfft = 2048; window = boxcar(nfft/4); noverlap = nfft/8; fs = 10; [s,f] = pwelch(V1,window,noverlap,nfft,fs,'onesided'); subplot(2,1,2); loglog(f,s,'k-',n,s1,'r--'); %- 第10点的目标谱与模拟谱比较 xlabel('freq/Hz'); ylabel('S/(m2s-1)');