读入文件实现涡识别
读入文件“vortex2.dat”中格式为tecplot输出的瞬时流场数据。本文章将采用涡识别准则对流场进行处理,显示出流场中的涡结构,并对比分析不同涡识别方法的性能差异。
if ~exist("vortex2","var")
load('matlab.mat')
end
vortex3=vortex2.';
vortex3=vortex3(:);
len=length(vortex3);
len9=len/9;
mysj=zeros(len9,9);
for i=1:9
mysj(:,i)=vortex3((1+(i-1)*len9):(i*len9));
end
[x,y,z,u,v,w]=deal(mysj(:,1),mysj(:,2),mysj(:,3),mysj(:,4),mysj(:,5),mysj(:,6));
xmm = sort(unique(x));
ymm = sort(unique(y));
zmm = sort(unique(z));
alldatalength=length(x);
xlength=length(xmm);
ylength=length(ymm);
zlength=length(zmm);
xgg= 1:xlength;
newint_x = round(interp1(xmm, xgg, x, 'linear'));
ygg= 1:ylength;
newint_y = round(interp1(ymm, ygg, y, 'linear'));
zgg= 1:zlength;
newint_z = round(interp1(zmm, zgg, z, 'linear'));
u_3d=zeros(xlength,ylength,zlength);
v_3d=zeros(xlength,ylength,zlength);
w_3d=zeros(xlength,ylength,zlength);
for i=1:alldatalength
u_3d(newint_x(i),newint_y(i),newint_z(i))=u(i);
v_3d(newint_x(i),newint_y(i),newint_z(i))=v(i);
w_3d(newint_x(i),newint_y(i),newint_z(i))=w(i);
end
gradu_x=diff(u_3d,1,1)./repmat(diff(xmm),1,ylength,zlength);
gradu_y=diff(u_3d,1,2)./permute(repmat(diff(ymm),1,xlength,zlength),[2,1,3]);
gradu_z=diff(u_3d,1,3)./permute(repmat(diff(zmm),1,ylength,xlength),[3,2,1]);
gradv_x=diff(v_3d,1,1)./repmat(diff(xmm),1,ylength,zlength);
gradv_y=diff(v_3d,1,2)./permute(repmat(diff(ymm),1,xlength,zlength),[2,1,3]);
gradv_z=diff(v_3d,1,3)./permute(repmat(diff(zmm),1,ylength,xlength),[3,2,1]);
gradw_x=diff(w_3d,1,1)./repmat(diff(xmm),1,ylength,zlength);
gradw_y=diff(w_3d,1,2)./permute(repmat(diff(ymm),1,xlength,zlength),[2,1,3]);
gradw_z=diff(w_3d,1,3)./permute(repmat(diff(zmm),1,ylength,xlength),[3,2,1]);
Q=zeros(xlength,ylength,zlength);
for ii=1:xlength
for jj=1:ylength
for kk=1:zlength
i=min(ii,xlength-1);
j=min(jj,ylength-1);
k=min(kk,zlength-1);
Q(i,j,k)=-0.5*(gradu_x(i,j,k).^2 ...
+gradv_y(i,j,k).^2 ...
+gradw_z(i,j,k).^2 ...
+2*gradu_y(i,j,k)*gradv_x(i,j,k) ...
+2*gradu_z(i,j,k)*gradw_x(i,j,k) ...
+2*gradv_z(i,j,k)*gradw_y(i,j,k));
end
end
end
[xdd, ydd,zdd] = meshgrid(ymm, xmm,zmm);
figure;
isosurface(xdd,ydd,zdd,Q,0.1);
view([1,1,-1]);
title('Q');
axis equal;
camlight;
lighting phong;
PP=zeros(xlength,ylength,zlength);
RR=zeros(xlength,ylength,zlength);
QQ=zeros(xlength,ylength,zlength);
RR2=zeros(xlength,ylength,zlength);
QQ2=zeros(xlength,ylength,zlength);
delta=zeros(xlength,ylength,zlength);
for ii=1:xlength
for jj=1:ylength
for kk=1:zlength
i=min(ii,xlength-1);
j=min(jj,ylength-1);
k=min(kk,zlength-1);
PP(i,j,k)=-(gradu_x(i,j,k)+gradv_y(i,j,k)+gradw_z(i,j,k));
QQ(i,j,k)=(-gradu_y(i,j,k)*gradv_x(i,j,k)-gradu_z(i,j,k)*gradw_x(i,j,k) ...
-gradv_z(i,j,k)*gradw_y(i,j,k)+gradu_x(i,j,k)*gradv_y(i,j,k) ...
+gradw_z(i,j,k)*gradu_x(i,j,k)+gradw_z(i,j,k)*gradv_y(i,j,k));
RR(i,j,k)=gradu_x(i,j,k)*(gradv_z(i,j,k)*gradw_y(i,j,k)-gradv_y(i,j,k)*gradw_z(i,j,k))...
+gradv_y(i,j,k)*(gradv_x(i,j,k)*gradw_z(i,j,k)-gradw_x(i,j,k)*gradu_x(i,j,k))...
+gradu_z(i,j,k)*(gradw_x(i,j,k)*gradv_y(i,j,k)-gradv_x(i,j,k)*gradw_y(i,j,k));
RR2(i,j,k)=(RR(i,j,k)+(2/27)*PP(i,j,k)*3-QQ(i,j,k)*PP(i,j,k)./3);
QQ2(i,j,k)=(QQ(i,j,k)-PP(i,j,k)*2/3);
delta(i,j,k)=(QQ2(i,j,k)./3)*3+(RR2(i,j,k)./2)*2;
end
end
end
[xdd, ydd,zdd] = meshgrid(ymm, xmm,zmm);
figure;
isosurface(xdd,ydd,zdd,delta,0.1);
view([1,1,-1]);
title('delta');
axis equal;
camlight;
lighting phong;