POD.zip_POD_POD 速度场_pod分析_云图_速度矢量分布

% create matrix will all fluctuating velocity components for each snapshot in a column Uall=[reshape(Uf,ni*nj,ns);reshape(Vf,ni*nj,ns)]; Uall = [U3;V3]; %可以修改U1 U2 U3 V1 V2 V3 %Do POD ANALYSIS R=Uall'*Uall; % 构建速度矩阵Autocovariance matrix [eV,D]=eig(R); % 求特征值和特征向量solve: eV is eigenvectors, D is eigenvalues in diagonal matrix [L,I]=sort(diag(D)); % 特征值排序sort eigenvalues in ascending order - I is sorted index vector for i=1:length(D)%特征向量排序 eValue(length(D)+1-i)=L(i); % Eigenvalues sorted in descending order eVec(:,length(D)+1-i)=eV(:,I(i)); % Eigenvectors sorted in the same order end; eValue(length(eValue))=0; % last eigenvalue should be zero menergy=eValue/sum(eValue); % relative energy associated with mode m % 重建前10mode的速度calculate the first 10 modes for i=1:10 tmp=Uall*eVec(:,i); % find mode phi(:,i)=tmp/norm(tmp); % normalize mode end; a = phi'*Uall; phi = phi*a; %区域一重建 position1 = csvread( 'E:\reynoldsstress\2d2cexcel\2-10cm-Q5\S25_H41.6_Q10.093_P1058_2.4yj04fs6.000000.csv' , 9 ,2 ,[ 9 , 2 , 15922 , 3]); % 读取速度坐标位置 P_x = position1(:,1); %取横坐标为P_x P_y = position1(:,2); %取纵坐标为P_y u11 = phi(1:15914,1); %提取mode1的速度 v11 = phi(15915:31828,1); [X,Y] = meshgrid (15.5:16:2335.5,15.5:16:1743.5); %按坐标位置建立坐标系网格 uu11 = griddata(P_x,P_y,u11,X,Y); %差值出离散的点 vv11 = griddata(P_x,P_y,v11,X,Y); figure(1) mesh(X,Y,vv11) %绘图 figure(2) mesh(X,Y,uu11) view grid on mask=csvread( 'E:\reynoldsstress\Export.4y5uf8ys.000000_2.csv' , 9 ,4 ,[ 9 , 4 , 4158344 , 4]); maskreshape=reshape(mask,2352,1768); for m=1:16:2352 for n=1:16:1760 masked = maskreshape(m:m+15,n:n+15); x=(m-1)/16+1; y=(n-1)/16+1; summasked(x,y) = sum(sum(masked))/256; %maskchange(m) = (mask(n+15)+ mask(n+14)+mask(n+13)+mask(n+12)+mask(n+11)+mask(n+10)+mask(n+9)+mask(n+8)+mask(n+7)+mask(n+6)+mask(n+5)+mask(n+4)+mask(n+3)+mask(n+2)+mask(n+1)+mask(n))/16; end end summasked = -1*(summasked-1); position1 = csvread( 'E:\reynoldsstress\2d2cexcel\2-10cm-Q5\S25_H41.6_Q10.093_P1058_2.4yj04fs6.000000.csv' , 9 ,2 ,[ 9 , 2 , 15922 , 3]); P_x = position1(:,1); %取横坐标为P_x P_y = position1(:,2); %取纵坐标为P_y %u11 = UVreconstruction(1:15914,536); x=P_x(:,1)/43.84; y=P_y(:,1)/43.84; u11 = phi(1:15914,1); u11 = reshape(u11,146,109); u11 = u11.*summasked(1:146,1:109); u11 = reshape(u11,15914,1); %v11 = UVreconstruction(15915:31828,536); v11 = phi(15915:31828,1); v11 = reshape(v11,146,109); v11 = v11.*summasked(1:146,1:109); v11 = reshape(v11,15914,1); [X,Y] = meshgrid (15.5:16:2335.5,15.5:16:1743.5); %按坐标位置建立坐标系网格 uu11 = griddata(P_x,P_y,u11,X,Y); %差值出离散的点 vv11 = griddata(P_x,P_y,v11,X,Y); figure(1) mesh(X,Y,vv11) %绘图 figure(2) mesh(X,Y,uu11) view grid on VV1= phi(15915:31828,536); VVV1 = reshape(VV1,146,109); VVVV1 = VVV1(1:146,1:109).*summasked(1:146,1:109); VVVV1 = reshape(VVVV1,15914,1); UU1= phi(1:15914,536); UUU1 = reshape(UU1,146,109); UUUU1 = UUU1(1:146,1:109).*summasked(1:146,1:109); UUUU1 = reshape(UUUU1,15914,1); figure(1) quiver(x,y,UUUU1,VVVV1); axis equal grid on hold on UVSCALER = sqrt(VVVV1.*VVVV1+UUUU1.*UUUU1); z=UVSCALER(:,1); scatter(x,y,5,z)%散点图 axis equal figure [X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x))',linspace(min(y),max(y)),'v4');%插值 pcolor(X-25,Y,Z);shading interp%伪彩色图 axis equal figure,contourf(X-25,Y,Z) %等高线图 hold on quiver(x-25,y,UUUU1,VVVV1,1);%矢量图 figure,mesh(X,Y,Z)%三维曲面 %区域二 position2 = csvread( 'E:\reynoldsstress\2d2cexcel\3-14cm-Q5\S25_H41.6_Q10.093_P1058_2.4y5uy2xr.000000.csv' , 9 ,2 ,[ 9 , 2 , 15922 , 3]); % 读取速度坐标位置 Q_x = position2(:,1); %取横坐标为P_x Q_y = position2(:,2); %取纵坐标为P_y u22 = phi(1:15914,1); v22 = phi(15915:31828,1); [X,Y] = meshgrid (15.5:16:2335.5,15.5:16:1743.5); %按坐标位置建立坐标系网格 uu22 = griddata(Q_x,Q_y,u22,X,Y); %差值出离散的点 vv22 = griddata(Q_x,Q_y,v22,X,Y); figure(1) mesh(X,Y,vv22) %绘图 figure(2) mesh(X,Y,uu22) view grid on mask=csvread( 'E:\reynoldsstress\Export.4y33nbuh.000000_3.csv' , 9 ,4 ,[ 9 , 4 , 4158344 , 4]); maskreshape=reshape(mask,2352,1768); for m=1:16:2352 for n=1:16:1760 masked = maskreshape(m:m+15,n:n+15); x=(m-1)/16+1; y=(n-1)/16+1; summasked(x,y) = sum(sum(masked))/256; %maskchange(m) = (mask(n+15)+ mask(n+14)+mask(n+13)+mask(n+12)+mask(n+11)+mask(n+10)+mask(n+9)+mask(n+8)+mask(n+7)+mask(n+6)+mask(n+5)+mask(n+4)+mask(n+3)+mask(n+2)+mask(n+1)+mask(n))/16; end end summasked = -1*(summasked-1); position2 = csvread( 'E:\reynoldsstress\2d2cexcel\3-14cm-Q5\S25_H41.6_Q10.093_P1058_2.4y5uy2xr.000000.csv' , 9 ,2 ,[ 9 , 2 , 15922 , 3]); % 读取速度坐标位置 Q_x = position2(:,1); %取横坐标为P_x Q_y = position2(:,2); %取纵坐标为P_y x=Q_x(:,1)/43.84; y=Q_y(:,1)/43.84; u22 = phi(1:15914,5); u22 = reshape(u22,146,109); u22 = u22.*summasked(1:146,1:109); u22 = reshape(u22,15914,1); v22 = phi(15915:31828,5); v22 = reshape(v22,146,109); v22 = v22.*summasked(1:146,1:109); v22 = reshape(v22,15914,1); [X,Y] = meshgrid (15.5:16:2335.5,15.5:16:1743.5); %按坐标位置建立坐标系网格 uu22 = griddata(Q_x,Q_y,u22,X,Y); %差值出离散的点 vv22 = griddata(Q_x,Q_y,v22,X,Y); figure(1) mesh(X,Y,vv22) %绘图 figure(2) mesh(X,Y,uu22) view grid on VV2= phi(15915:31828,5); VVV2 = reshape(VV2,146,109); VVVV2 = VVV2(1:146,1:109).*summasked(1:146,1:109); VVVV2 = reshape(VVVV2,15914,1); UU2= phi(1:15914,5); UUU2 = reshape(UU2,146,109); UUUU2 = UUU2(1:146,1:109).*summasked(1:146,1:109); UUUU2 = reshape(UUUU2,15914,1); figure(1) quiver(x,y,UUUU2,VVVV2); axis equal grid on hold on UVSCALER = sqrt(VVVV2.*VVVV2+UUUU2.*UUUU2); z=UVSCALER(:,1); scatter(x,y,5,z)%散点图 axis equal figure [X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x))',linspace(min(y),max(y)),'v4');%插值 pcolor(X-25,Y,Z);shading interp%伪彩色图 axis equal figure,contourf(X-25,Y,Z) %等高线图 hold on quiver(x-25,y,UUUU2,VVVV2,1);%矢量图 figure,mesh(X,Y,Z)%三维曲面 %区域三 %MASK过程 mask=csvread( 'E:\reynoldsstress\Export.4y7k1izg.000000_4.csv' , 9 ,4 ,[ 9 , 4 , 4158344 , 4]); maskreshape=reshape(mask,2352,1768); for m=1:16:2352 for n=1:16:1760 masked = maskreshape(m:m+15,n:n+15); x=(m-1)/16+1; y=(n-1)/16+1; summasked(x,y) = sum(sum(masked))/256; %maskchange(m) = (mask(n+15)+ mask(n+14)+mask(n+13)+mask(n+12)+mask(n+11)+mask(n+10)+mask(n+9)+mask(n+8)+mask(n+7)+mask(n+6)+mask(n+5)+mask(n+4)+mask(n+3)+mask(n+2)+mask(n+1)+mask(n))/16; end end summasked = -1*(summasked-1); position3 = csvread( 'E:\reynoldsstress\2d2cexcel\3-14cm-Q5\S25_H41.6_Q10.093_P1058_2.4y5uy2xr.000000.csv' , 9 ,2 ,[ 9 , 2 , 15922 , 3]); % 读取速度坐标位置 R_x = position3(:,1); R_y = position3(:,2); x=R_x(:,1)/43.84; y=R_y(:,1)/43.84; u33 = phi(1:15914,1); u33 = reshape(u33,146,109); u33 = u33.*summasked(1:146,1:109); u33 = reshape(u33,15914,1); v33 = phi(15915:31828,1); v33 = reshape(v33,146,109); v33 = v33.*summasked(1:146,1:109); v33 = reshape(v33,15914,1); [X,Y] = meshgrid (15.5:16:2335.5,15.5:16:1743.5); %按坐标位置建立坐标系网格 uu33 = griddata(R_x,R_y,u33,X,Y); %差值出离散的点 vv33 = griddata(R_x,R_y,v33,X,Y); figure(1) mesh(X,Y,vv33) %绘图 figure(2) mesh(X,Y,uu33) view grid on VV3= phi(15915:31828,1); VVV3 = reshape(VV3,146,109); VVVV3 = VVV3(1:146,1:109).*summasked(1:146,1:109); VVVV3 = reshape(VVVV3,15914,1); UU3= phi(1:15914,1); UUU3 = reshape(UU3,146,109); UUUU3 = UUU3(1:146,1:109).*summasked(1:146,1:109); UUUU3 = reshape(UUUU3,15914,1); figure(3) quiver(x,y,UUUU3,VVVV3); axis equal grid on hold on UVSCALER = sqrt(VVVV3.*VVVV3+UUUU3.*UUUU3); z=UVSCALER(:,1); scatter(x,y,5,z)%散点图 axis equal figure [X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x))',linspace(min(y),max(y)),'v4');%