%这段代码之前发过，结束后生成图也都贴出来了，但是很多地方没有写出详细的说明，加上王富海的3.2图做的一塌糊涂，力的计算引用自王星，但是王星的学位论文画图也是字母全标错了，当时看到这里也是欲哭无泪。拜托你们毕业的能不能认真一点。所以在这再发一次，对一些内容进行补充。
%还是老规矩先宣传一下QQ群： 格子玻尔兹曼救星：293267908。不收费的哦，就是为了早点毕业建的群。
%这个例子采用 MRT-LBM 模拟旋转圆柱绕流 %左边速度边界-泊肃叶流，右边压力边界，上下无滑移壁面（全部用非平衡外推格式） %基于 MRT-LBM 的流场与声场仿真计算 --王富海2017 clc clear close all % 设置仿真参数
uMax=0.05; %中间最大速度 Re=40;%雷诺数 yLen=81;%垂直方向格子数 xLen=161;%水平方向格子数 cylinderD=(yLen-1)/5;%圆的半径 nu=uMax*cylinderD/Re;%运动粘性 Cs=sqrt(1/3);%格子声速 tau=1/2+3*nu;%松弛时间 omega=1/tau;%松弛频率 step=120000;%最大迭代次数 rho_out=1;%出口密度 uo2=0.1;%圆柱的旋转速度 rhoo=1;%初始化密度 checkStep=100;%收敛计算间隔 saveStep=20;%保存结果间隔 % file Path=uigetdir(cd)%仿真中间过程图片的保存路径 VSSum=[];%所有节点格子速度总和-每 saveStep 步记录一次 VSSum2=[];%所有节点格子速度总和-每 checkStep 步记录一次，监视收敛曲线 dx=1;%相邻格子节点水平间隔 dy=1;%相邻格子节点垂直间隔 dt=1; y1=1;%下边界垂直坐标 y2=yLen;%上边界垂直坐标 GG=uMax/((y1-y2)^2/4); for j=1:yLen U_in(1,j)=uMax;%左边的速度剖面入口 V_in(1,j)=0; end %为上下壁面边界点速度赋值 ub(1:xLen,1)=0; vb(1:xLen,1)=0; ub(1:xLen,yLen)=0; vb(1:xLen,yLen)=0; f_u(1,1)=0; f_u(xLen,yLen)=0; f_v(1,1)=0; f_v(xLen,yLen)=0; %导入拉格朗日点坐标--圆的方程式（x-xPos）^2+(y-yPos)^2=r^2; xPos=2*cylinderD+cylinderD/2;%圆心坐标 x r=cylinderD/2; yPos=(yLen+1)/2;%圆心坐标 y % 计算与垂直网格的交点 lagPosChuiZhi=[];%存储拉格朗日点-位置索引 x_Start=1; x_Stop=xLen; y_Start=1; y_Stop=yLen; for i=x_Start:dx:x_Stop delta=r^2-(i-xPos)^2; if delta>0 y1=yPos+sqrt(delta);%垂直交点的真正 y1 坐标 y2=yPos-sqrt(delta);%垂直交点的真正 y2 坐标 y1_index=(y1-y_Start)/dy+1;%垂直交点的 y 网格索引 y2_index=(y2-y_Start)/dy+1;%垂直交点的 y 网格索引 x_index=(i-x_Start)/dx+1;%x 的网格索引 lagPosChuiZhi=[lagPosChuiZhi;x_index y1_index;x_index y2_index]; elseif delta==0 y1=yPos+sqrt(delta);%垂直交点的真正 y 坐标 y1_index=(y1-y_Start)/dy+1; x_index=(i-x_Start)/dx+1; lagPosChuiZhi=[lagPosChuiZhi;x_index y1_index]; end end %计算与水平网格的交点 lagPosShuiPing=[];%存储拉格朗日点-位置索引 for i2=y_Start:dy:y_Stop delta=r^2-(i2-yPos)^2; if delta>0 x1=xPos+sqrt(delta);%水平交点的真正 x1 坐标 x2=xPos-sqrt(delta);%水平交点的真正 x2 坐标 x1_index=(x1-x_Start)/dx+1;%水平交点的 x 网格索引 x2_index=(x2-x_Start)/dx+1;%水平交点的 x 网格索引 y_index=(i2-y_Start)/dy+1;%y 的网格索引 lagPosShuiPing=[lagPosShuiPing;x1_index y_index;x2_index y_index]; elseif delta==0 x1=xPos+sqrt(delta);%水平交点的真正 x 坐标 x1_index=(x1-x_Start)/dx+1; y_index=(i2-y_Start)/dy+1; lagPosShuiPing=[lagPosShuiPing;x1_index y_index]; end end % 为拉格朗日点附上速度 lenLagShuiPing=length(lagPosShuiPing); for i=1:lenLagShuiPing xTemp=lagPosShuiPing(i,1); yTemp=lagPosShuiPing(i,2); if yTemp-yPos>=0 && xTemp-xPos>=0 %逆时针第一象限 thetaTemp=asin((yTemp-yPos)/r); end if yTemp-yPos>=0 && xTemp-xPos<0 %第二象限 thetaTemp=pi-asin((yTemp-yPos)/r); end if yTemp-yPos<0 && xTemp-xPos>=0 %第四象限 thetaTemp=2*pi+asin((yTemp-yPos)/r); end if yTemp-yPos<0 && xTemp-xPos<0 %第三象限 thetaTemp=pi-asin((yTemp-yPos)/r); end %为拉格朗日点速度赋值,uo2 >0 顺时针。 lagPosShuiPing(i,3)=uo2*sin(thetaTemp);%u 将uo分解为水平速度 lagPosShuiPing(i,4)=-uo2*cos(thetaTemp);%v。1象限，uo顺时针方向分解后的v为负值，然而cos当 0~90度却为正值，所以加负号 end lenLagChuiZhi=length(lagPosChuiZhi); for i=1:lenLagChuiZhi xTemp=lagPosChuiZhi(i,1); yTemp=lagPosChuiZhi(i,2); if yTemp-yPos>=0 && xTemp-xPos>=0 thetaTemp=asin((yTemp-yPos)/r); end if yTemp-yPos>=0 && xTemp-xPos<0 thetaTemp=pi-asin((yTemp-yPos)/r); end if yTemp-yPos<0 && xTemp-xPos>=0 thetaTemp=2*pi+asin((yTemp-yPos)/r); end if yTemp-yPos<0 && xTemp-xPos<0 thetaTemp=pi-asin((yTemp-yPos)/r); end %为拉格朗日点速度赋值,uo2 >0 顺时针。 lagPosChuiZhi(i,3)=uo2*sin(thetaTemp);%u lagPosChuiZhi(i,4)=-uo2*cos(thetaTemp);%v end %将索引位置转化为实际坐标 real_lagPosShuiPing=[(lagPosShuiPing(:,1)-1)*dx+x_Start,(lagPosShuiPing(:,2)-1)*dy+y_Start]; %？？？这其实不就是lagPosShuiPing吗？*dy+y_Start干嘛？ % plot(real_lagPosShuiPing(:,1),real_lagPosShuiPing(:,2),'r*')%绘制水平交点-拉格朗日点 %将索引位置转化为实际坐标 real_lagPosChuiZhi=[(lagPosChuiZhi(:,1)-1)*dx+x_Start,(lagPosChuiZhi(:,2)-1)*dy+y_Start]; % plot(real_lagPosChuiZhi(:,1),real_lagPosChuiZhi(:,2),'bo')%绘制垂直交点-拉格朗日点 % axis equal % hold off
% 计算朗格朗日附近格点-索引 lenShuiPing=length(lagPosShuiPing); XY_left_index=[]; XY_right_index=[]; for i=1:lenShuiPing w1=floor(lagPosShuiPing(i,1)); w2=ceil(lagPosShuiPing(i,1)); if w1== w2 %如果对于lagPosShuiPing中的x值相等，则-1、+1之后成为两个相关点 XY_left_index=[XY_left_index;w1-1 lagPosShuiPing(i,2)]; XY_right_index=[XY_right_index; w1+1 lagPosShuiPing(i,2)]; else %如果对于lagPosShuiPing中的x值不相等，则floor 和ceil分别成为两个相关点 XY_left_index=[XY_left_index;w1 lagPosShuiPing(i,2)]; XY_right_index=[XY_right_index; w2 lagPosShuiPing(i,2)]; end end lenChuiZhi=length(lagPosChuiZhi); XY_up_index=[]; XY_down_index=[]; for i=1:lenChuiZhi w1=floor(lagPosChuiZhi(i,2)); w2=ceil(lagPosChuiZhi(i,2)); if w1== w2 XY_up_index=[XY_up_index; lagPosChuiZhi(i,1) w1+1]; XY_down_index=[XY_down_index; lagPosChuiZhi(i,1) w1-1]; else XY_up_index=[XY_up_index;lagPosChuiZhi(i,1) w2]; XY_down_index=[XY_down_index;lagPosChuiZhi(i,1) w1]; end end %round：四舍五入取整数 XY_left_index=round(XY_left_index); %round：四舍五入取整数 XY_right_index=round(XY_right_index); XY_up_index=round(XY_up_index); XY_down_index=round(XY_down_index);
%将 A 和 B 的每一行都视为单个实体，并返回 A 和 B 的共有行，但是不包含重复项。必须指定 A 和 B，setOrder 是可选的。'rows' 选项不支持元胞数组 [~,~,c1]=intersect(XY_left_index,XY_up_index,'rows'); %c1是XY_up_index中重复的部分的编号 [~,~,c2]=intersect(XY_right_index,XY_up_index,'rows'); same_Up_LR=[c1;c2];%重复索引的行索引。 若一个点的速度被垂直和水平都修改了，就在后面取均值 [~,~,c3]=intersect(XY_left_index,XY_down_index,'rows'); [~,~,c4]=intersect(XY_right_index,XY_down_index,'rows'); same_Down_LR=[c3;c4];%重复索引的行索引 % %D2Q9 模型参数 w=[1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36 4/9 ]%各个方向的权重 cx=[ 1 0 -1 0 1 -1 -1 1 0];%各方向 x 速度分量 cy=[ 0 1 0 -1 1 1 -1 -1 0];%各方向 y 速度分量 M=[1 1 1 1 1 1 1 1 1;-4 -1 -1 -1 -1 2 2 2 2;4 -2 -2 -2 -2 1 1 1 1; ... 0 1 0 -1 0 1 -1 -1 1;0 -2 0 2 0 1 -1 -1 1;0 0 1 0 -1 1 1 -1 -1; ... 0 0 -2 0 2 1 1 -1 -1;0 1 -1 1 -1 0 0 0 0;0 0 0 0 0 1 -1 1 -1;] % 由于分布函数 0 索引被放置在 matlab 索引 9 的位置，所以将 M 第一列放到最后。 M=circshift(M,[-1 -1]); s1(9)=0; s1(3)=s1(9); s1(5)=s1(9); s1(7)=omega; s1(8)=s1(7); s1(4)=1.2; s1(6)=s1(4); s1(1)=omega; s1(2)=s1(1)-0.1; s=diag(s1);%对角矩阵-松弛因子 Minv=inv(M); Minv_s=Minv*s; %网格节点初始化，初始化各节点密度=1，初始化各节点速度为0，初始化各节点分布函数=平衡分布函数 for i=1:xLen for j=1:yLen rho(i,j)=rhoo; u(i,j)=uMax; v(i,j)=0; t1=u(i,j)^2+v(i,j)^2;%速度 U*V 内积,t1 写到外面，节省计算量 for k=1:9 t2=u(i,j)*cx(k)+v(i,j)*cy(k);%U*C 内积 f(k,i,j)=rho(i,j)*w(k)*(1+3*t2+4.5*t2^2-1.5*t1);% 计算初始态下的 分布函数=平衡分布函数 Fi(k,i,j)=0;%初始化离散的边界作用力分布函数,初始值赋为 0 end end end % %主循环------------------------------------------------------------------------------ tstart = tic;%计算时钟 %产生网格数据：[meshX,meshY] [meshX,meshY]=meshgrid(1:xLen,yLen:-1:1);% ！由于MATLAB画的图，y的值自下而上是y值减小的，所以在这里对每个x对应的y值进行自上而下的减小 C_zhuli=[0]; %阻力（有兴趣的可以从拼音猜一下作者是哪里人） C_shengli=[0];%升力 figure %召唤图片，阻力升力图！ fplot=plot(C_zhuli) hold on fplot2=plot(C_shengli,'r') % ylim([-3 3]) figure %召唤图片， pc=pcolor(meshX,meshY,ones(yLen,xLen));%绘制伪彩图 shading interp%伪彩色图---shading interp 会区分每个线形区域的颜色,并且插入与其相近的颜色 axis equal ylim([1 yLen]) %y轴上下限设定ylim([a,b]) ，要不然图很奇怪的 for kk=1:step step2=kk%显示程序的步进 mm=max(max(rho))%观察最大密度 if mm>100%判断仿真是否发散 break; end
%1 碰撞过程 for i=1:xLen for j=1:yLen % %速度空间到动量空间 m=M*f(1:9,i,j); %计算平衡动量 m_eq(9,1)=rho(i,j); m_eq(1,1)=rho(i,j)*(-2+3*(u(i,j)^2+v(i,j)^2)); m_eq(2,1)=rho(i,j)*(1-3*(u(i,j)^2+v(i,j)^2)); m_eq(3,1)=rho(i,j)*u(i,j); m_eq(4,1)=-rho(i,j)*u(i,j); m_eq(5,1)=rho(i,j)*v(i,j); m_eq(6,1)=-rho(i,j)*v(i,j); m_eq(7,1)=rho(i,j)*(u(i,j)^2-v(i,j)^2); m_eq(8,1)=rho(i,j)*u(i,j)*v(i,j); %在动量空间松弛 m_temp=Minv_s*(m-m_eq)-Minv*Fi(:,i,j); for k=1:9 f(k,i,j)=f(k,i,j)-m_temp(k); end end end % 2 迁移过程 f2=f;%将碰撞后分布函数保存给变量 f2 for i=1:9 f(i,:,: ) = circshift(f2(i,:,: ), [0,cx(i),cy(i)]); end %非平衡外推格式-上下无滑移壁面 %下边界 for i=2:xLen-1% ？？？这里我也是很奇怪，为何总是不算边界点？加上又怎么样？ %计算相邻点的速度 j=2; rho(i,j)=sum(f(:,i,j));%计算密度 uSum=0; vSum=0; for k=1:9 uSum=uSum+f(k,i,j)*cx(k); vSum=vSum+f(k,i,j)*cy(k); end u(i,j)=uSum/rho(i,j); v(i,j)=vSum/rho(i,j); %计算边界上的平衡分布函数（用相邻点密度的代替边界密度） t1=u(i,j)^2+v(i,j)^2;%速度 U*V 内积 t1b=ub(i,1)^2+vb(i,1)^2; for k=1:9%[2 5 6]%每个节点 9 个分量 t2=u(i,j)*cx(k)+v(i,j)*cy(k);%U*C 内积 t2b=ub(i,1)*cx(k)+vb(i,1)*cy(k);%U*C 内积 % 计算边界平衡分布函数-用相邻点的密度代替边界密度 feq(k,i,1)=rho(i,j)*w(k)*(1+3*t2b+4.5*t2b^2-1.5*t1b); % 计算相邻点非平衡分布函数 fneq(k,i,j)=f(k,i,j)-rho(i,j)*w(k)*(1+3*t2+4.5*t2^2-1.5*t1); %计算边界分布函数 f(k,i,1)=feq(k,i,1)+fneq(k,i,j); end end for i=2:xLen-1%上边界 %计算相邻点的速度 j=yLen-1; rho(i,j)=sum(f(:,i,j));%计算密度 uSum=0; vSum=0; for k=1:9 uSum=uSum+f(k,i,j)*cx(k); vSum=vSum+f(k,i,j)*cy(k); end u(i,j)=uSum/rho(i,j); v(i,j)=vSum/rho(i,j); %计算边界上的平衡分布函数（用相邻点的密度代替边界密度） t1=u(i,j)^2+v(i,j)^2;%速度 U*V 内积 t1b=ub(i,yLen)^2+vb(i,yLen)^2; for k=1:9%[4 7 8]%每个节点 9 个分量 t2=u(i,j)*cx(k)+v(i,j)*cy(k);%U*C 内积 t2b=ub(i,yLen)*cx(k)+vb(i,yLen)*cy(k);%U*C 内积 % 计算边界平衡分布函数-用相邻点的密度代替边界密度 feq(k,i,yLen)=rho(i,j)*w(k)*(1+3*t2b+4.5*t2b^2-1.5*t1b); % 计算相邻点非平衡分布函数 fneq(k,i,j)=f(k,i,j)-rho(i,j)*w(k)*(1+3*t2+4.5*t2^2-1.5*t1); %计算边界分布函 f(k,i,yLen)=feq(k,i,yLen)+fneq(k,i,j); end end %---------------------------------------------左边界！ %非平衡外推模式-施加左右边界条件,修正所有分布函数变量 %左速度边界-非平衡外推格式 u(1,:)=U_in; v(1,:)=V_in; for j=1:yLen %计算相邻点的速度 i=2; rho(i,j)=sum(f(:,i,j));%计算密度 uSum=0; vSum=0; for k=1:9 uSum=uSum+f(k,i,j)*cx(k); vSum=vSum+f(k,i,j)*cy(k); end u(i,j)=uSum/rho(i,j); v(i,j)=vSum/rho(i,j); for k=1:9%[1 5 8]%每个节点 9 个分量 % 计算边界平衡分布函数-用相邻点的密度代替边界密度 i=1; t1=u(i,j)^2+v(i,j)^2;%速度 U*V 内积 t2=u(i,j)*cx(k)+v(i,j)*cy(k);%U*C 内积 feq(k,1,j)=rho(2,j)*w(k)*(1+3*t2+4.5*t2^2-1.5*t1); % 计算相邻点非平衡分布函数 i=2; t1=u(i,j)^2+v(i,j)^2;%速度 U*V 内积 t2=u(i,j)*cx(k)+v(i,j)*cy(k);%U*C 内积 fneq(k,i,j)=f(k,i,j)-rho(i,j)*w(k)*(1+3*t2+4.5*t2^2-1.5*t1); %计算边界分布函数 f(k,1,j)=feq(k,1,j)+fneq(k,i,j); end end % 右边界-压力出口- for j=1:yLen %计算相邻点的速度 i=xLen-1; rho(i,j)=sum(f(:,i,j));%计算密度 uSum=0; vSum=0; for k=1:9 uSum=uSum+f(k,i,j)*cx(k); vSum=vSum+f(k,i,j)*cy(k); end u(i,j)=uSum/rho(i,j); v(i,j)=vSum/rho(i,j); %计算边界上的平衡分布函数（边界上的密度，用相邻点的速度代替边界速度） t1=u(i,j)^2+v(i,j)^2;%速度 U*V 内积 for k=1:9%[3 6 7]%每个节点 9 个分量 t2=u(i,j)*cx(k)+v(i,j)*cy(k);%U*C 内积 % 计算边界平衡分布函数-用相邻点的速度代替边界速度 feq(k,xLen,j)=rho_out*w(k)*(1+3*t2+4.5*t2^2-1.5*t1); % 计算相邻点非平衡分布函数 fneq(k,i,j)=f(k,i,j)-rho(i,j)*w(k)*(1+3*t2+4.5*t2^2-1.5*t1); %计算边界分布函数 f(k,xLen,j)=feq(k,xLen,j)+fneq(k,i,j); end end
%计算宏观量 for i=1:xLen for j=1:yLen rho(i,j)=sum(f(:,i,j));%计算每一点的密度 end end for i=1:xLen for j=1:yLen uSum=0; vSum=0; for k=1:9 uSum=uSum+f(k,i,j)*cx(k); vSum=vSum+f(k,i,j)*cy(k); end u(i,j)=uSum/rho(i,j);%计算速度场 v(i,j)=vSum/rho(i,j); end end
%3.5 圆柱边界 %已经得到了拉格朗日点的坐标点位置，分别为 lagPosChuiZhi 和 lagPosShuiPing %3.5.1 进行水平拉格朗日交点的左右点速度修正 u_Before=u; v_Before=v; len_ShuiPing=length(lagPosShuiPing); for i=1:len_ShuiPing %对左边的点修正采用二阶拉格朗日插值修正 xb=lagPosShuiPing(i,1); x1=XY_left_index(i,1); x2=XY_left_index(i,1)-1; x3=XY_right_index(i,1)+1; y=XY_left_index(i,2); %用 xb,x2,x3 修正 x1 %从左到右依次是 x3 x2 x1? xb u_xb=lagPosShuiPing(i,3);%固定边界 u_x1=u_Before(x1,y); u_x2=u_Before(x2,y); u_x3=u_Before(x3,y); %二阶拉格朗日插值 u_x1_new=u_xb*(x1-x2)*(x1-x3)/((xb-x2)*(xb-x3))+u_x2*(x1-xb)*(x1-x3)/((x2-xb)*(x2-x3))+u_x3*(... x1-xb)*(x1-x2)/((x3-xb)*(x3-x2)); v_xb=lagPosShuiPing(i,4);%固定边界 v_x1=v_Before(x1,y); v_x2=v_Before(x2,y); v_x3=v_Before(x3,y); v_x1_new=v_xb*(x1-x2)*(x1-x3)/((xb-x2)*(xb-x3))+v_x2*(x1-xb)*(x1-x3)/((x2-xb)*(x2-x3))+v_x3*(... x1-xb)*(x1-x2)/((x3-xb)*(x3-x2)); %更新速度值 u(x1,y)=u_x1_new; v(x1,y)=v_x1_new; %计算反馈力 rho_temp=sum(f(:,x1,y)); f_u(x1,y)=2*rho_temp*(u_x1_new-u_x1)/dt; f_v(x1,y)=2*rho_temp*(v_x1_new-v_x1)/dt; %计算离散的边界作用力分布函数 Fi(1:9,x1,y)=getFi(u(x1,y),v(x1,y),f_u(x1,y),f_v(x1,y),s1);
%对右边的点修正 xb=lagPosShuiPing(i,1); x1=XY_right_index(i,1); x2=XY_left_index(i,1)-1; x3=XY_right_index(i,1)+1; y=XY_right_index(i,2); %用 xb,x2,x3 修正 x1 %从左到右依次是 xb x1-? x2 x3 u_xb=lagPosShuiPing(i,3);%固定边界 u_x1=u_Before(x1,y); u_x2=u_Before(x2,y); u_x3=u_Before(x3,y); %二阶拉格朗日插值
u_x1_new=u_xb*(x1-x2)*(x1-x3)/((xb-x2)*(xb-x3))+u_x2*(x1-xb)*(x1-x3)/((x2-xb)*(x2-x3))+u_x3*(... x1-xb)*(x1-x2)/((x3-xb)*(x3-x2)); v_xb=lagPosShuiPing(i,4);%固定边界 v_x1=v_Before(x1,y); v_x2=v_Before(x2,y); v_x3=v_Before(x3,y); v_x1_new=v_xb*(x1-x2)*(x1-x3)/((xb-x2)*(xb-x3))+v_x2*(x1-xb)*(x1-x3)/((x2-xb)*(x2-x3))+v_x3*(... x1-xb)*(x1-x2)/((x3-xb)*(x3-x2)); %更新速度值 u(x1,y)=u_x1_new; v(x1,y)=v_x1_new; %计算反馈力 rho_temp=sum(f(:,x1,y)); f_u(x1,y)=2*rho_temp*(u_x1_new-u_x1)/dt; f_v(x1,y)=2*rho_temp*(v_x1_new-v_x1)/dt; %计算离散的边界作用力分布函数 Fi(1:9,x1,y)=getFi(u(x1,y),v(x1,y),f_u(x1,y),f_v(x1,y),s1); %注意这个位置的x1已经不是第一个Fi用的left点，而是right点了 end
len_ChuiZhi=length(lagPosChuiZhi); for i=1:len_ChuiZhi % 对上边的点修正 yb=lagPosChuiZhi(i,2); y1=XY_up_index(i,2); y2=XY_up_index(i,2)+1; y3=XY_down_index(i,2)-1; x=XY_up_index(i,1); %用 yb,y2,y3 修正 y1 %从上到下依次是 y3 y2 y1-? yb u_yb=lagPosChuiZhi(i,3);%固定边界 u_y1=u_Before(x,y1); u_y2=u_Before(x,y2); u_y3=u_Before(x,y3); %二阶拉格朗日插值 u_y1_new=u_yb*(y1-y2)*(y1-y3)/((yb-y2)*(yb-y3))+u_y2*(y1-yb)*(y1-y3)/((y2-yb)*(y2-y3))+u_y3*(... y1-yb)*(y1-y2)/((y3-yb)*(y3-y2)); v_yb=lagPosChuiZhi(i,4);%固定边界 v_y1=v_Before(x,y1); v_y2=v_Before(x,y2); v_y3=v_Before(x,y3);
v_y1_new=v_yb*(y1-y2)*(y1-y3)/((yb-y2)*(yb-y3))+v_y2*(y1-yb)*(y1-y3)/((y2-yb)*(y2-y3))+v_y3*(... y1-yb)*(y1-y2)/((y3-yb)*(y3-y2)); %更新新的节点速度值 %判断是不是需要 2 次修正 if ismember(i,same_Up_LR) %C = ismember( A,B )--> C 为 A 中元素能否在 B 中找到。例如：A = [1,2,3]; B = [2,1,1,4]; 则 C = [1,1,0]; u(x,y1)=0.5*(u(x,y1)+u_y1_new); v(x,y1)=0.5*(v(x,y1)+v_y1_new); else u(x,y1)=u_y1_new; v(x,y1)=v_y1_new; end %计算 f 反馈力 rho_temp=sum(f(:,x,y1)); f_u(x,y1)=2*rho_temp*(u(x,y1)-u_y1)/dt; f_v(x,y1)=2*rho_temp*(v(x,y1)-v_y1)/dt; %计算离散的边界作用力分布函数 Fi(1:9,x,y1)=getFi(u(x,y1),v(x,y1),f_u(x,y1),f_v(x,y1),s1); % 对下边的点修正，取欲修正点的下边两个点和边界点，采用二阶拉格朗日插值修正 yb=lagPosChuiZhi(i,2); y1=XY_down_index(i,2); y2=XY_up_index(i,2)+1; y3=XY_down_index(i,2)-1; x=XY_down_index(i,1); %用 yb,y2,y3 修正 y1 %从上到下依次是 yb y1? y3 y2 u_yb=lagPosChuiZhi(i,3);%固定边界 u_y1=u_Before(x,y1); u_y2=u_Before(x,y2); u_y3=u_Before(x,y3); %二阶拉格朗日插值 u_y1_new=u_yb*(y1-y2)*(y1-y3)/((yb-y2)*(yb-y3))+u_y2*(y1-yb)*(y1-y3)/((y2-yb)*(y2-y3))+u_y3*(... y1-yb)*(y1-y2)/((y3-yb)*(y3-y2)); v_yb=lagPosChuiZhi(i,4);%固定边界 v_y1=v_Before(x,y1); v_y2=v_Before(x,y2); v_y3=v_Before(x,y3);
v_y1_new=v_yb*(y1-y2)*(y1-y3)/((yb-y2)*(yb-y3))+v_y2*(y1-yb)*(y1-y3)/((y2-yb)*(y2-y3))+v_y3*(... y1-yb)*(y1-y2)/((y3-yb)*(y3-y2)); %判断是不是需要 2 次修正 if ismember(i,same_Down_LR) u(x,y1)=0.5*(u(x,y1)+u_y1_new); v(x,y1)=0.5*(v(x,y1)+v_y1_new); else u(x,y1)=u_y1_new; v(x,y1)=v_y1_new; end %计算 f 反馈力 rho_temp=sum(f(:,x,y1)); f_u(x,y1)=2*rho_temp*(u(x,y1)-u_y1)/dt; f_v(x,y1)=2*rho_temp*(v(x,y1)-v_y1)/dt; %计算离散的边界作用力分布函数 Fi(1:9,x,y1)=getFi(u(x,y1),v(x,y1),f_u(x,y1),f_v(x,y1),s1); end %每 saveStep 步保存下误差，检查收敛性 if mod(kk,saveStep)==0 temp=0; for i=1:xLen for j=1:yLen temp=temp+sqrt(u(i,j)^2+v(i,j)^2);%节点的速度和 end end VSSum2=[VSSum2;temp]; u2=rot90(u); v2=rot90(v); for i=1:xLen for j=1:yLen speed(j,i)=sqrt(u2(j,i)^2+v2(j,i)^2);%节点的速度 end end set(pc,'cdata',speed); drawnow % saveas(gcf,[file Path,'\',num2str(kk/saveStep),'.tif']) % pic Length=kk/saveStep;%保存的图片数量 % ff=0; ff2=0; for i=1:xLen for j=1:yLen %if (i-xPos)^2+(j-yPos)^2-r^2>0 ff=ff+f_u(i,j); ff2=ff2+f_v(i,j); %end end end %阻力系数与升力系数 C_zhuli(kk/saveStep)=-2*ff/(rhoo*uMax^2*cylinderD); C_shengli(kk/saveStep)=-2*ff2/(rhoo*uMax^2*cylinderD); set(fplot,'ydata',C_zhuli); set(fplot2,'ydata',C_shengli); drawnow %每 checkStep 次保存下误差 if mod(kk,checkStep)==0 VSSum=[VSSum;temp]; save; %计算节点速度和变化（网格点总数*误差），判断是否跳出循环 if length(VSSum)>=2 && abs(VSSum(end)-VSSum(end-1))/VSSum(end-1)<=(1e-6) fprintf('迭代步数:%d \n',kk); fprintf('迭代误差:%d\n',abs(VSSum(end)-VSSum(end-1))/(xLen*yLen)); break; end end end end runtime = toc(tstart);%仿真用的时间 %结果后处理 u2=rot90(u); v2=rot90(v); for i=1:xLen for j=1:yLen speed(j,i)=sqrt(u2(j,i)^2+v2(j,i)^2); end end save( ['re' num2str(Re) '-' num2str(uMax) '-' num2str(uo2) '.mat']) %速度场彩图 figure pcolor(meshX,meshY,speed); shading interp axis equal tight ylim([1 yLen]) %绘制流线图 1 figure startY=linspace(1,xLen,100); startX=zeros(size(startY))+1; streamline(meshX,meshY,u2,v2,startX,startY,2) axis equal tight %绘制流线图 2 figure streamslice(meshX,meshY,u2,v2,8,'nearest ' ) ylim([1 yLen]) axis equal tight hold off % 动态观看保存的计算图片 % figure % for i=1:pic Length % i % picdata=imread([file Path,'\',num2str(i),'.tif']); % pause(0.05) % imshow(picdata); % end
