Matlab涡识别准则可视化

读入文件实现涡识别

读入文件“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;

运行结果：