无滑移边界的水冲小球例子，用浸没式边界（速度差法）求曳力--修改自王富海硕士论文

%这个例子采用 MRT-LBM 模拟旋转圆柱绕流
%左边速度边界-泊肃叶流，右边压力边界，上下无滑移壁面（全部用非平衡外推格式）
 %力的计算来自：基于 MRT-LBM 的流场与声场仿真计算 --王富海2017
 %对照例题来自：P33，王露. 一种浸没边界格子Boltzmann方法的算法改进及在CFD中的应用[D].

% 还是老规矩先宣传一下QQ群： 格子玻尔兹曼救星：293267908。
 %计算雷诺数=40时静置小球的阻力和升力系数
clc
clear
close all
 
%  设置仿真参数

uMax=0.1; %中间最大速度
Re=40;%雷诺数
yLen=15*24;%垂直方向格子数
xLen=30*24;%水平方向格子数
cylinderD=24;%圆的半径
nu=uMax*cylinderD/Re;%运动粘性
Cs=sqrt(1/3);%格子声速
tau=1/2+3*nu;%松弛时间
omega=1/tau;%松弛频率
step=8000;%最大迭代次数
rho_out=1;%出口密度
uo2=0 ;%圆柱的旋转速度
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=6.25*cylinderD ;%圆心坐标 x
r=cylinderD/2;
yPos=7.5*cylinderD;%圆心坐标 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