MATLAB Examples

Complex geometry

The wind field is now simulated on a circle defined in the vertical plane (yz) normal to the flow. The function windSim2 is therefore used with two input files:

  • 'windData.txt'
  • 'circle.txt';

'windData.txt' is similar to the input file used in Example1.m and Example2.m, but does not contain information about the grid to be built. The main reason, is that there is no more grid ! The geometry is defined in the second input file.

'circle.txt' is a .txt files with two columns, where the first one contains the coordinates of the nodes along axis y, and the second columns contains the coordinates of the nodes along axis z.

Contents

Function call

clearvars;close all;clc;
rng(1)
windData = 'windData.txt';
myGeometry = 'circle.txt';
[u,v,w,t,nodes] = windSim2(windData,myGeometry);

U = u+nodes.U*ones(size(t)); % get fluctuating + mean value
Elapsed time is 30.350374 seconds.

Display the geometry

figure
plot(nodes.Y,nodes.Z,'ko','markerfaceColor','r')
ylabel('z')
xlabel('y')
title('geometry')
set(gcf,'color','w')
axis square

Time series overview

figure
box on; hold on
plot(t,U(1:2:6,:))
xlabel('time (s) ');
zlabel('u (m/s)')
ylabel('z (m)')
set(gcf,'color','w')
legend('node1','node3','node 5')
xlim([0,300])

Vertical mean wind speed profile

zr = 10; % 10 meters
a = 0.15 ; % cf windData.txt
U_target= nodes.U(1).*(nodes.Z./zr).^(a);
figure
hold on; box on;
[sortZ,indZ]=sort(nodes.Z,'ascend');
plot(mean(U(indZ,:),2),sortZ,'k')
plot(nodes.U,nodes.Z,'ko','markerfacecolor','r')
legend('target','simulated','location','NorthWest')
ylabel('z (m) ');
xlabel('U (m/s)')
set(gcf,'color','w')

Standard deviation of wind velocity fluctuations

figure
box on;hold on,grid on
plot(nodes.Y, std(u,0,2),'ko','markerFaceColor','r')
plot(nodes.Y, std(v,0,2),'ko','markerFaceColor','y')
plot(nodes.Y, std(w,0,2),'ko','markerFaceColor','b')
ylim([0,3])
xlabel('y (m)')
ylabel('\sigma_{i} (m/s), i = (u,v,w)')
set(gcf,'color','w')
legend('\sigma_{u}','\sigma_{v}','\sigma_{w}','location','SouthWest')

Turbulence intensity

The present code specifies a turbulent intensity that decreases for increasing mean wind velocity, i.e. for increasing altitudes. This is more realistic than the constant turbulent intensity introduced in a previous version.

Iu = std(u,0,2)./nodes.U*100;
Iv = std(v,0,2)./nodes.U*100;
Iw = std(w,0,2)./nodes.U*100;

figure
box on;hold on,grid on
plot(nodes.U,Iu,'ko','markerFaceColor','r')
plot(nodes.U,Iv,'ko','markerFaceColor','y')
plot(nodes.U,Iw,'ko','markerFaceColor','b')
xlabel('U (m/s)')
ylabel('TI (%)')
set(gcf,'color','w')
legend('Iu','Iv','Iw','location','SouthWest')

Co-coherence

time step & sampling frequency

dt = median(diff(t));
fs = 1/dt;
tmax = t(end);
f0 = 1/tmax;
[nNodes,N] = size(u);
% block duration (each)
tBlock = 240;
% number of blocks
NBlock = round(tmax/tBlock);
% number of data point per block
Ncoh = round(N/NBlock);
clear Su Sv Sw

if mod(Ncoh,2),
    cocoh = zeros(nNodes,nNodes,round(Ncoh/2));
else
    cocoh = zeros(nNodes,nNodes,round(Ncoh/2)+1);
end
% computation of the cocoherence using the function cohere
for ii=1:nNodes,
    for jj=1:nNodes,
        [cocoh(ii,jj,:),~,freq] = coherence(u(ii,:),u(jj,:),Ncoh,round(Ncoh/2),Ncoh,fs);
    end
end

Verification: coherence between nodes 1 and 3

nodeID1 = 1;
nodeID2 = 3;

dz = abs(nodes.Z(nodeID1)-nodes.Z(nodeID2));
dy = abs(nodes.Y(nodeID1)-nodes.Y(nodeID2));

figure
box on;hold on
plot(freq,squeeze(cocoh(nodeID1,nodeID2,:)),'ko','markerfacecolor','r');
DecY_u = dy.*7.*freq;
DecZ_u = dz.*10.*freq;
myU = 1/2*(nodes.U(nodeID1)+nodes.U(nodeID2));
plot(freq,exp(-sqrt(DecY_u.^2+DecZ_u.^2)./myU),'k');
legend(['dy = ',num2str(dy,2),' m, and dz = ',num2str(dz,2),' m'],'target');
xlim([0,0.2])

ylabel('co-coherence')
if ii>2,
    xlabel('f (Hz)')
end
set(gcf,'color','w')

Verification: coherence between nodes 1 and 10

nodeID1 = 1;
nodeID2 = 10;

dz = abs(nodes.Z(nodeID1)-nodes.Z(nodeID2));
dy = abs(nodes.Y(nodeID1)-nodes.Y(nodeID2));

figure
box on;hold on
plot(freq,squeeze(cocoh(nodeID1,nodeID2,:)),'ko','markerfacecolor','r');
DecY_u = dy.*7.*freq;
DecZ_u = dz.*10.*freq;
myU = 1/2*(nodes.U(nodeID1)+nodes.U(nodeID2));
plot(freq,exp(-sqrt(DecY_u.^2+DecZ_u.^2)./myU),'k');
legend(['dy = ',num2str(dy,2),' m, and dz = ',num2str(dz,2),' m'],'target');
xlim([0,0.2])
ylabel('co-coherence')
if ii>2,
    xlabel('f (Hz)')
end
set(gcf,'color','w')

Verification: coherence between nodes 20 and 25

nodeID1 = 20;
nodeID2 = 25;

dz = abs(nodes.Z(nodeID1)-nodes.Z(nodeID2));
dy = abs(nodes.Y(nodeID1)-nodes.Y(nodeID2));

figure
box on;hold on
plot(freq,squeeze(cocoh(nodeID1,nodeID2,:)),'ko','markerfacecolor','r');
DecY_u = dy.*7.*freq;
DecZ_u = dz.*10.*freq;
myU = 1/2*(nodes.U(nodeID1)+nodes.U(nodeID2));
plot(freq,exp(-sqrt(DecY_u.^2+DecZ_u.^2)./myU),'k');
legend(['dy = ',num2str(dy,2),' m, and dz = ',num2str(dz,2),' m'],'target');
xlim([0,0.2])
ylabel('co-coherence')
if ii>2,
    xlabel('f (Hz)')
end
set(gcf,'color','w')

Conclusions

A turbulent wind field is simulated on a circle, masde of 30 nodes in the plane (yz). The grid has been removed, and replaced by a geometry generated by a text file. It seems that the spatial and temporal correlation implemented works pretty well. The removal of the grid improves the application range of windSim2.m