采用 MRT-LBM 模拟旋转圆柱绕流2---MATLAB代码--王富海2017--基于 MRT-LBM 的流场与声场仿真计算
%这段代码之前发过，结束后生成图也都贴出来了，但是很多地方没有写出详细的说明，加上王富海的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
采用 MRT-LBM 模拟旋转圆柱绕流2---MATLAB代码--王富海2017--基于 MRT-LBM 的流场与声场仿真计算相关推荐
- 采用 MRT-LBM 模拟旋转圆柱绕流---MATLAB代码--王富海2017--基于 MRT-LBM 的流场与声场仿真计算
%这个例子采用 MRT-LBM 模拟旋转圆柱绕流 %基于 MRT-LBM 的流场与声场仿真计算 --王富海2017 %左边速度边界-泊肃叶流,右边压力边界,上下无滑移壁面(全部用非平衡外推格式) %还 ...
- 【光学】基于matlab模拟双孔干涉附matlab代码
1 内容介绍 通过Matlab软件编程,实现光学双缝干涉的计算机仿真.仿真结果对学生理解光学原理的基本概念很有帮助,提高了教学效果. 2 部分代码 %REDME!!!该仿真模拟了双孔干事实验中光屏逐 ...
- 采用POD以及DMD方法实现圆柱绕流流动分解（POD篇）
背景 笔者通过对POD及DMD两种模型降阶方法进行学习后,尝试对Brunton大佬书中Re=100圆柱绕流case进行复现,现将学习及复现过程的代码及结果进行分享,希望多提宝贵建议. 本文圆柱绕流案例 ...
- catia圆柱转化为圆台_浅析actran气动噪声仿真技术，以圆柱绕流气动噪声仿真为例...
一.写在前面Actran是fft(Free Field Technologies)公司的旗舰产品,"号称"市场上最先进最完善的声学模拟软件(引用官方语言),覆盖振动声学和流动声学的 ...
- 浅析actran气动噪声仿真技术，以圆柱绕流气动噪声仿真为例
附赠仿真学习包,包含结构.流体.电磁.热仿真等多学科视频教程,点击领取: 仿真秀粉丝专属礼包 作者:小禹治水,仿真秀科普作者 一.写在前面 Actran是fft(Free Field Te ...
- LBM学习记录5 Python实现IB二维圆柱绕流 1.0
浸润边界(IB)下二维圆柱绕流,bug1.0版 import numpy as np from math import sqrt from numpy import pi as PI from num ...
- Python-LBM(格子玻尔兹曼)程序源码实例分析—圆柱绕流篇
初次学习LBM计算方法,找到一个比较优秀的用python语言编写的圆柱绕流的实例,对每段代码详细添加了注释,帮助自己总结,也为初学的朋友们提供一点帮助(全部代码在文章最后). 先放一张结果图像: 1. ...
- fluent瞬态计算终止条件在哪里设置_Fluent案例7【圆柱绕流】
一个瞬态的圆柱绕流案例 知识点: 瞬态圆柱绕流的模拟 一个后处理的方法:将瞬态模型中一个点的速度变化绘成图表并将数值导出excel文件 模型如下图所示,左边界为速度边界进口速度0.5m/s,试模拟出计 ...
- Fluent UDF 实现用Newmark-β方法计算圆柱绕流流固耦合时的位移振动响应
Fluent UDF 实现用Newmark-β方法计算圆柱绕流流固耦合时的位移振动响应 问题描述 代码 尚未解决的问题 问题描述 拟用Fluent模拟圆柱振子在不同流速的风作用下的横向振动.采用二维模 ...
- OpenFOAMGmshCFD圆柱绕流（两个圆柱）
问题: 圆柱绕流问题,模拟仿真有两个圆柱.一个源的流体变化情况. 解决步骤: 1.使用Gmsh画出网格,并保存cylindertwo.msh 2.以Cavity为基础创建新的Case:Cylinder ...
最新文章
- 几种特征选择方法的比较，孰好孰坏？
- HTTP协议的挑战者：RSocket
- EMS server Tibco
- MinGW和MSYS的自动安装 【转】
- linux主机解析虚拟机超时_Linux 内核超时导致虚拟机无法正常启动
- kitti数据集_KITTI数据集激光雷达坐标系下的里程计真值
- 阅读react-redux源码(五) - connectAdvanced中store改变的事件转发、ref的处理和pure模式的处理
- mysql 添加唯一索引_浅谈Mysql索引
- Datawhale 一周年，生日快乐！
- Chilkat -- python又一个强大的库
- python读取word文件内容_[python]读取word文档中的数据，整理成excel表
- cad字体安装_装了1个G的CAD字体后，我的CAD崩了怎么办？
- libaio.so.1 is needed by MySQL-server-5.5.48-1.linux2.6.i386
- get 和 post请求的区别
- matlab变压器温度仿真
- VUE3-Cesium（entity、primitive总结及材质的修改）
- 思考伯努利试验的两种组合思想
- 计算机组成原理期中考,计算机组成原理期中考卷
- vue实现点击不同按钮展示不同内容
- 元宇宙之XR（02）VR概念解读 分类说明