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

如果提前已经算好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;endendarrayPoint=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'];endif (i>=10&&i<100)newPath=[path,'\flap1220-a-90-3c-1-000',num2str(i),'.txt'];endif (i>=100&&i<1000)newPath=[path,'\flap1220-a-90-3c-1-00',num2str(i),'.txt'];endif (i>=1000&&i<10000)newPath=[path,'\flap1220-a-90-3c-1-0',num2str(i),'.txt'];endif (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;endI_index=length(X(1,:));J_index=length(X(:,1));%读取上一步的涡核位置,其坐标在网格中对应的索引记录在positionLastTimeStepfor I=1:I_indexif((X(1,I)<=corePosition(1,arrayPoint-1))&&(X(1,I+1)>corePosition(1,arrayPoint-1)))I_positionLastTimeStep=I;endendfor J=1:J_indexif((Y(J,1)<=corePosition(2,arrayPoint-1))&&(Y(J+1,1)>corePosition(2,arrayPoint-1)))J_positionLastTimeStep=J;endendisFouondMinOrMaxVorticity=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;%如果找不到就扩大范围endx_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.8contourLevel=contourLevel+1;gama(contFile)=getGamaFromVortexPosition(corePosition(1,arrayPoint-1),corePosition(2,arrayPoint-1),U,V,X,Y,VORTICITY,contourLevel);endendif(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;endendend
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;endendend
end

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

AI For Science— 基于AI求解2D非定常圆柱绕流，真的很流体！！
!In [2] %cd work/ !unzip PaddleScience_CubeDomain.zip AI求解CFD基础案例:圆柱绕流好看吗? 如果你是CFD计算流体力学领域的大牛,看了是不是 ...
NUMECA的数值计算模块Fine/Turbo的5种定常交界面方法
商用软件NUMECA的数值计算模块Fine/Turbo提供了5种定常交界面方法,分别为当地守恒(localcon servative coupling)法､周向守恒(conservative coup ...
人工压缩算法--定常原始变量不可压缩N-S方程
流体从非稳态到稳态的变化过程,这种非定常流动和定常流动的转化,实际上是方程类型的改变.人工压缩算法就是通过改变方程类型来求解不可压缩粘性流动的方法. 人工压缩算法的基本思路人工压缩算法基本思路为:如 ...
利用matlab实现POD分解(在一维信号或二维流场矢量中的应用）
利用matlab实现POD分解(在一维信号或二维流场矢量中的应用) 0 前言 0.1 matlab中特征值计算 0.2 matlab中SVD分解计算 0.3 信号的正交性 1 一维信号POD分解 1. ...
寅辞旧岁，卯定常虹丨ASKO洗碗机“净”护新春团圆时刻
农历新年是一年中最重要的节日,但过去三年的特殊时光阻碍了很多人的归乡之行,如今当阴霾逐渐散去,必然会引来大规模的新年归乡潮,奔赴一个久违的团圆年.美馔佳宴是新春佳节的永恒命题,新年家里少不了亲友的光临 ...
分支定界 matlab,使用MATLAB实现分枝定界法求解整数规划的详细资料说明
分支定界法是一种求解离散最优化问题的计算分析方法.它是由Land Doig和Dakin等人在20世纪60年代初提出的.分支定界法可求纯整数或混合整数线性规划问题,求解方法由分支和定界组成." ...
建立飞机的六自由度运动方程，并对飞机定常直线平飞状态进行配平
是当年大三的专业课,但是当时其实完全不懂要干什么,做的还蛮糟糕的.现在上了研究生至少明白了一些些,所以这两天把从前的作业又做了一下,当然还是有很多不懂的地方,期待进步吧. 一.飞机配平(定常直线平飞 ...
ORA-01858: 在要求输入数字处找到非数字字符 13行
文章目录 1. 现象 2. 分析 3. 解决方案 ORA-01858: 在要求输入数字处找到非数字字符13行 1. 现象 insert /*+append*/ into ASSET_LOAN(sele ...
报错“在要求输入数字处找到非数字字符”
在Java链接Oracle数据库的时候报错"在要求输入数字处找到非数字字符",一开始完全不知道问题所在,就开始在请求字段中找数字类型的变量,还是找不到问题所在就开始百度了: 回合1 ...

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

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

最新文章

热门文章