matlab找到非定常涡流的每个时间步的涡的涡核位置和这个涡环量(以及重叠网格扑翼流场的涡动力学参数求解的解决方案)

最新推荐文章于 2024-01-14 10:10:05 发布

小尚同学的科研日志小屋

最新推荐文章于 2024-01-14 10:10:05 发布

阅读量2.4k

点赞数 1

文章标签： matlab 流体力学

本文链接：https://blog.csdn.net/Shangjialin1/article/details/106470014

版权

如果提前已经算好fluent每个时间步的cas和dat文件，那么可以用jou文件导出每个节点的u,v,dudy,dvdx，overset-cell-type。fluent中只能直接导出涡量幅值，想要计算涡量只能用dudy,dvdx做差。

追踪涡核位置：
如果是重叠网格，则需要同时导出前景网格和背景网格。根据overset-cell-type来判断是否在matlab端接受数据。如果overset-cell-type<0则说明这个节点是被挖去的节点。之后在用griddata重新插值到自定义的规则网格上面。
追踪涡的位置需要提供涡的初始时刻大致位置。在附近的网格格点寻找涡量的极值点，如果初始位置涡量为正，则寻找极大值点，否则寻找极小值点。寻找极值就是把这个节点的涡量和上下左右四个点的涡量值比较。如果在周围33的范围内找不到就扩大范围到44直到可以找到极值点。但是真实流动中，涡的极值点可能不在网格格点上，而是在4个网格格点之间。根据高斯拟合将附近的网格格点上拟合出一条连续变化的曲线，得到涡量极值的准确位置。
重复这个步骤，总之，之后时间步的涡核位置就可以根据之前的时间步的位置来确定，从而得到每个时间步的涡核位置。

得到这个涡的环量随时间的变化：
假定涡的边界是由涡量等值线确定的，那么沿着这个等值线对速度进行线积分就可以获得环量。
我们可以假设我们规定等值线值为±5的环量等值线决定了一个正的涡或者负的涡的边界。
contour生成的等值线在matlab中的数据结构请参考matlab的帮助。
根据等值线的数据结构，我们可以得知每一条等值线的上面离散点的坐标还有等值线的值。根据inpolygon把包含指定涡核位置且等值线值为±5的等值线选出来就可以确定指定涡核的边界。从而进行曲线积分得到这个涡的环量

针对涡碰撞到扑翼上导致涡逐渐被扑翼吸收的情况的特殊处理：
如果涡环量产生跃变则说明发生了这种情况，涡的边界判定可以从等值线值为±5变为等值线值为±10总之可以缩小等值线的范围，从而把涡捕捉到。如果不缩小范围会把扑翼上的束缚涡也包括进去从而出现错误。
代码有些部分和上面思路不太一样，但是实在懒得修正了：

function [Gama]  =getGamaFromVortexPosition(xPosition,yPosition,U,V,X,Y,VORTICITY,contourLevel)
%获得指定涡核位置的涡的环量
% xPosition=-0.24;
% yPosition=-0.125;
% xPosition=corePosition(1,length(corePosition(1,:)));
% yPosition=corePosition(2,length(corePosition(1,:)));
maxLevel=80;
minLevel=-maxLevel;
space=1;
[C,~]=contour(X,Y,VORTICITY,minLevel:space:maxLevel);
% SIZE=size(C);
sign=-1;
if(interp2(X,Y,VORTICITY,xPosition,yPosition)>0)
    sign=1;
end
level=sign*contourLevel;%假定涡的边界是由值为contourLevel或-contourLevel的等值线构成的
arrayPoint=1;
flag=0;
while(flag==0)
    if((abs(C(1,arrayPoint)-level)<space/2))
        xVector=C(1,(arrayPoint+1):(arrayPoint+C(2,arrayPoint)));
        yVector=C(2,(arrayPoint+1):(arrayPoint+C(2,arrayPoint)));
        if(inpolygon(xPosition,yPosition,xVector,yVector)==1)
            break;
        end
    end
    arrayPoint=arrayPoint+C(2,arrayPoint)+1;
end
% while(abs(C(1,arrayPoint)-level)>space/2)
%     arrayPoint=arrayPoint+C(2,arrayPoint)+1;
%     arrayPoint=round(arrayPoint);
%     xVector=C(1,(arrayPoint+1):(arrayPoint+C(2,arrayPoint)));
%     yVector=C(2,(arrayPoint+1):(arrayPoint+C(2,arrayPoint)));
%     if(inpolygon(xPosition,yPosition,xVector,yVector)==1)
%         break;
%     end
% end
Gama=0;
u=interp2(X,Y,U,xVector,yVector);
v=interp2(X,Y,V,xVector,yVector);
%需要考虑顺时针和逆时针的问题
for i=2:length(xVector)
    Gama=Gama+u(i)*(xVector(i)-xVector(i-1))+v(i)*(yVector(i)-yVector(i-1));
end
Gama=Gama+u(1)*(xVector(1)-xVector(length(xVector)))+v(1)*(yVector(1)-yVector(length(xVector)));

function [corePosition,gama]  = FindVortexCoreAndCaculateGama(path,initialPosition)%寻找每个时间步涡核的位置和这个涡的环量
%参数说明：path：文件加路径。 initialPosition:初始位置
%path='D:\flap_case\C19_2_E';
%initialPosition(1)=-0.05802;
%initialPosition(2)=-0.01839;
timestepAutoSave=2;%每隔多少个时间步记录一次
initialTimestep=416+1*2;
totalTimeSteps=416+30*2;%总共的时间步
numberFiles=round((totalTimeSteps-initialTimestep)/timestepAutoSave);
corePosition=zeros(2,numberFiles);
corePosition(1,1)=initialPosition(1);%初始涡核位置的坐标
corePosition(2,1)=initialPosition(2);
gama=zeros(numberFiles,1);
arrayPoint=1;
contFile=0;
contourLevel=5;%等值线的值，用来确定涡的轮廓
for i=initialTimestep:timestepAutoSave:totalTimeSteps
%     i=initialTimestep+2*timestepAutoSave;
%     arrayPoint=3;
    contFile=contFile+1;
   if (i<10)
        newPath=[path,'\flap1220-a-90-3c-1-0000',num2str(i),'.txt'];
    end
    if (i>=10&&i<100)
        newPath=[path,'\flap1220-a-90-3c-1-000',num2str(i),'.txt'];
    end
    if (i>=100&&i<1000)
        newPath=[path,'\flap1220-a-90-3c-1-00',num2str(i),'.txt'];
    end
    if (i>=1000&&i<10000)
        newPath=[path,'\flap1220-a-90-3c-1-0',num2str(i),'.txt'];
    end
    if (i>10000)
        newPath=[path,'\flap1220-a-90-3c-1-',num2str(i),'.txt'];
    end
    [nodenumber,x,y,v,u,dudy,dvdx]=textread(newPath,'%f %f %f %f %f %f %f','headerlines',1);
    vorticity=zeros(1,length(nodenumber));
    for j=1:length(nodenumber)
        vorticity(j)=dvdx(j)-dudy(j);
    end
    %形成griddaata所需要的网格节点
    [X,Y]=meshgrid(min(x):(max(x)-min(x))/2000:max(x),min(y):(max(y)-min(y))/2000:max(y)); 
    VORTICITY=griddata(x,y,vorticity,X,Y);
    U=griddata(x,y,u,X,Y);
    V=griddata(x,y,v,X,Y);
    if(i==initialTimestep)%对第一个时间步的处理，由于初始位置涡核位置已知，只需要计算其涡量
        gama(1)=getGamaFromVortexPosition(corePosition(1,1),corePosition(2,1),U,V,X,Y,VORTICITY,contourLevel);
        arrayPoint=arrayPoint+1;
        continue;
    end
    I_index=length(X(1,:));
    J_index=length(X(:,1));
    %读取上一步的涡核位置,其坐标在网格中对应的索引记录在positionLastTimeStep
    for I=1:I_index
        if((X(1,I)<=corePosition(1,arrayPoint-1))&&(X(1,I+1)>corePosition(1,arrayPoint-1)))
            I_positionLastTimeStep=I;
        end
    end
    for J=1:J_index
        if((Y(J,1)<=corePosition(2,arrayPoint-1))&&(Y(J+1,1)>corePosition(2,arrayPoint-1)))
            J_positionLastTimeStep=J;
        end
    end
    isFouondMinOrMaxVorticity=0;%是否找到极小值点或极大值点，根据上一步的环量来确定
    foundLength=3;%寻找涡核附近的索引范围
    while(isFouondMinOrMaxVorticity==0)
        [currentPositionIndexI,currentPositionIndexJ,isFouondMinOrMaxVorticity]=findMinOrMaxValue(I_positionLastTimeStep,J_positionLastTimeStep,foundLength,VORTICITY);
%         VorticityShow=VORTICITY((J_positionLastTimeStep-foundLength):1:(J_positionLastTimeStep+foundLength),(I_positionLastTimeStep-foundLength):1:(I_positionLastTimeStep+foundLength));
        foundLength=foundLength+1;%如果找不到就扩大范围
    end
    x_test=[X(currentPositionIndexJ,currentPositionIndexI-1),X(currentPositionIndexJ,currentPositionIndexI),X(currentPositionIndexJ,currentPositionIndexI+1)];
    y_test=[Y(currentPositionIndexJ-1,currentPositionIndexI),Y(currentPositionIndexJ,currentPositionIndexI),Y(currentPositionIndexJ+1,currentPositionIndexI)];
    Vor_x_test=[VORTICITY(currentPositionIndexJ,currentPositionIndexI-1),VORTICITY(currentPositionIndexJ,currentPositionIndexI),VORTICITY(currentPositionIndexJ,currentPositionIndexI+1)];
    Vor_y_test=[VORTICITY(currentPositionIndexJ-1,currentPositionIndexI),VORTICITY(currentPositionIndexJ,currentPositionIndexI),VORTICITY(currentPositionIndexJ+1,currentPositionIndexI)];
    corePosition(1,arrayPoint)=crit_interp_p(Vor_x_test,x_test);
    corePosition(2,arrayPoint)=crit_interp_p(Vor_y_test,y_test);
    gama(contFile)=getGamaFromVortexPosition(corePosition(1,arrayPoint),corePosition(2,arrayPoint),U,V,X,Y,VORTICITY,contourLevel);
    arrayPoint=arrayPoint+1;
    if(arrayPoint>3)%运行两个时间步之后检测环量跃变
        if(abs(gama(contFile)-gama(contFile-1))>1.8*abs(gama(contFile-1)-gama(contFile-2)))%跃变的系数是1.8
            contourLevel=contourLevel+1;
            gama(contFile)=getGamaFromVortexPosition(corePosition(1,arrayPoint-1),corePosition(2,arrayPoint-1),U,V,X,Y,VORTICITY,contourLevel);
        end
    end
    if(abs(interp2(X,Y,VORTICITY,corePosition(1,arrayPoint-1),corePosition(2,arrayPoint-1)))<contourLevel)
        break;%如果预测涡核位置的涡量绝对值小于涡的预测边界的涡量则认为涡消失
    end
end
maxLevel=80;
minLevel=-maxLevel;
space=1;
plot(corePosition(1,:),corePosition(2,:),'*-','LineWidth',3);hold on;
back_time_step=0;
plot(corePosition(1,1:(end-back_time_step)),corePosition(2,1:(end-back_time_step)),'*-','LineWidth',3);hold on;
contour(X,Y,VORTICITY(),minLevel:space:maxLevel);
plot(gama);

function  [currentPositionIndexI,currentPositionIndexJ,isFouondMinOrMaxVorticity]=findMinOrMaxValue(I_positionLastTimeStep,J_positionLastTimeStep,foundLength,VORTICITY)
isFouondMinOrMaxVorticity=0;
currentPositionIndexI=-1;
currentPositionIndexJ=-1;
%VORTICITY((J_positionLastTimeStep-foundLength):1:(J_positionLastTimeStep+foundLength),(I_positionLastTimeStep-foundLength):1:(I_positionLastTimeStep+foundLength));
if(VORTICITY(J_positionLastTimeStep,I_positionLastTimeStep)>0)
    for i=(I_positionLastTimeStep-foundLength):1:(I_positionLastTimeStep+foundLength)
        for j=(J_positionLastTimeStep-foundLength):1:(J_positionLastTimeStep+foundLength)
%             j=J_positionLastTimeStep-foundLength+5;
%             i=I_positionLastTimeStep-foundLength+3;
            if(VORTICITY(j,i)==max(max(VORTICITY((J_positionLastTimeStep-foundLength):1:(J_positionLastTimeStep+foundLength),(I_positionLastTimeStep-foundLength):1:(I_positionLastTimeStep+foundLength)))))%>VORTICITY(j,i+1)&&VORTICITY(j,i)>VORTICITY(j,i-1)&&VORTICITY(j,i)>VORTICITY(j+1,i)&&VORTICITY(j,i)>VORTICITY(j-1,i))
                isFouondMinOrMaxVorticity=1;
                currentPositionIndexI=i;
                currentPositionIndexJ=j;
                break;
            end
        end
    end
end
if(VORTICITY(J_positionLastTimeStep,I_positionLastTimeStep)<=0)
    for i=(I_positionLastTimeStep-foundLength):1:(I_positionLastTimeStep+foundLength)
        for j=(J_positionLastTimeStep-foundLength):1:(J_positionLastTimeStep+foundLength)
            if(VORTICITY(j,i)==min(min(VORTICITY((J_positionLastTimeStep-foundLength):1:(J_positionLastTimeStep+foundLength),(I_positionLastTimeStep-foundLength):1:(I_positionLastTimeStep+foundLength)))))%<VORTICITY(j,i+1)&&VORTICITY(j,i)<VORTICITY(j,i-1)&&VORTICITY(j,i)<VORTICITY(j+1,i)&&VORTICITY(j,i)<VORTICITY(j-1,i))
                isFouondMinOrMaxVorticity=1;
                currentPositionIndexI=i;
                currentPositionIndexJ=j;
                break;
            end
        end
    end
end

小尚同学的科研日志小屋

关注

1
点赞
踩
20

收藏

觉得还不错? 一键收藏
0
评论
matlab找到非定常涡流的每个时间步的涡的涡核位置和这个涡环量(以及重叠网格扑翼流场的涡动力学参数求解的解决方案)

如果提前已经算好fluent每个时间步的cas和dat文件，那么可以用jou文件导出每个节点的u,v,dudy,dvdx，overset-cell-type。fluent中只能直接导出涡量幅值，想要计算涡量只能用dudy,dvdx做差。追踪涡核位置：如果是重叠网格，则需要同时导出前景网格和背景网格。根据overset-cell-type来判断是否在matlab端接受数据。如果overset-cell-type<0则说明这个节点是被挖去的节点。之后在用griddata重新插值到自定义的规则网格上面。
复制链接

扫一扫