多孔介质流的算例（附Fortran计算代码及MATLAB后处理代码）
! 我又来啦，做了个小算例，在100*100个格子的空腔中加了几个矩形障碍，然后自己设置了雷诺数=400，水流速度为0.08.这些
!大家都可以随意设置哈，设置方法看下之前我分享的教程。.
!这次后处理遇到点问题！当设置格子数为2500的时候（当时我想算个精细点的），MATLAB提示：Error using repmat: Requested 252601*252601(475.4GB) array exceeds maximum array size preference. 大家有什么其他办法进行大图片后处理请告诉我！QQ1542170984.
!上代码！西边界为水流进入区
parameter (n=100,m=100)
real f(0:8,0:n,0:m)
real feq(0:8,0:n,0:m),rho(0:n,0:m)
real w(0:8), cx(0:8),cy(0:8)
real u(0:n,0:m), v(0:n,0:m)
integer i
open(2,file='uvfield.txt')
open(3,file='uvely.txt')
open(4,file='vvelx.txt')
open(8,file='timeu.txt')
uo=0.08 !有变动，设Re=400，H=100,L=100, v=0.01，密度1000（水）
sumvelo=0.0
rhoo=1000
dx=1.0
dy=dx
dt=1.0
alpha=0.01
Re=uo*m/alpha
print *, "Re=", Re
omega=1.0/(3.*alpha+0.5)
mstep=4000 !注意修改步数
w(0)=4./9.
do i=1,4
w(i)=1./9.
end do
do i=5,8
w(i)=1./36.
end do
cx(0)=0
cx(1)=1
cx(2)=0
cx(3)=-1
cx(4)=0
cx(5)=1
cx(6)=-1
cx(7)=-1
cx(8)=1
cy(0)=0
cy(1)=0
cy(2)=1
cy(3)=0
cy(4)=-1
cy(5)=1
cy(6)=1
cy(7)=-1
cy(8)=-1
do j=0,m
do i=0,n
rho(i,j)=rhoo
u(i,j)=0.0
v(i,j)=0.0
end do
end do
do j=1,m-1
u(0,j)=uo
v(0,j)=0.0
end do
! main loop
1 do kk=1,mstep
call collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
call streaming(f,n,m)
! ——————————–
call sfbound(f,n,m,uo)
call rhouv(f,rho,u,v,cx,cy,n,m)
print *, u(0,m/2),v(0,m/2),rho(0,m/2),u(n,m/2),v(n,m/2),rho(n,m/2)
write(8,*) kk,u(n/2,m/2),v(n/2,m/2)
END DO
! end of the main loop
call result(u,v,rho,uo,n,m)
stop
end
! end of the main program
subroutine collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
real f(0:8,0:n,0:m)
real feq(0:8,0:n,0:m),rho(0:n,0:m)
real w(0:8), cx(0:8),cy(0:8)
real u(0:n,0:m), v(0:n,0:m)
DO i=0,n
DO j=0,m
t1=u(i,j)*u(i,j)+v(i,j)*v(i,j)
DO k=0,8
t2=u(i,j)*cx(k)+v(i,j)*cy(k)
feq(k,i,j)=rho(i,j)*w(k)*(1.0+3.0*t2+4.50*t2*t2-1.50*t1)
f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j)
END DO
END DO
END DO
return
end
subroutine streaming(f,n,m)
real f(0:8,0:n,0:m)
! streaming
DO j=0,m
DO i=n,1,-1 !RIGHT TO LEFT
f(1,i,j)=f(1,i-1,j)
END DO
DO i=0,n-1 !LEFT TO RIGHT
f(3,i,j)=f(3,i+1,j)
END DO
END DO
DO j=m,1,-1 !TOP TO BOTTOM
DO i=0,n
f(2,i,j)=f(2,i,j-1)
END DO
DO i=n,1,-1
f(5,i,j)=f(5,i-1,j-1)
END DO
DO i=0,n-1
f(6,i,j)=f(6,i+1,j-1)
END DO
END DO
DO j=0,m-1 !BOTTOM TO TOP
DO i=0,n
f(4,i,j)=f(4,i,j+1)
END DO
DO i=0,n-1
f(7,i,j)=f(7,i+1,j+1)
END DO
DO i=n,1,-1
f(8,i,j)=f(8,i-1,j+1)
END DO
END DO
return
end
subroutine sfbound(f,n,m,uo)
real f(0:8,0:n,0:m)
! bounce back on north boundary
do i=0,n
f(4,i,m)=f(2,i,m)
f(7,i,m)=f(5,i,m)
f(8,i,m)=f(6,i,m)
end do
! bounce back on south boundary
do i=0,n
f(2,i,0)=f(4,i,0)
f(5,i,0)=f(7,i,0)
f(6,i,0)=f(8,i,0)
end do
! moving lid west boundary
do j=1,m
f(1,0,j)=f(3,0,j)+2*uo*rhow/3.0
rhow=(f(0,0,j)+f(2,0,j)+f(4,0,j)+2.*(f(3,0,j)+f(6,0,j)+f(7,0,j)))/(1-uo)
f(5,0,j)=f(7,0,j)+rhow*uo/6.0 !注意此处与中文P86不同 为什么默认(f(2,0,j)-f(4,0,j))*0.5 这一项等于零？？？？？？？？
f(8,0,j)=f(6,0,j)+rhow*uo/6.0 !注意此处与P86不同
end do
! the east outlet：open boundray condiction
do j=1,m
f(1,n,j)=2*f(1,n-1,j)-f(1,n-2,j)
f(5,n,j)=2*f(5,n-1,j)-f(5,n-2,j)
f(8,n,j)=2*f(8,n-1,j)-f(8,n-2,j)
end do
!!!!!自己随便加--边界 1
!bounce back on north boundary
do i=30,40
f(2,i,70)=f(4,i,70)
f(5,i,70)=f(7,i,70)
f(6,i,70)=f(8,i,70)
end do
!bounce back on south boundary
do i=30,40
f(4,i,50)=f(2,i,50)
f(7,i,50)=f(5,i,50)
f(8,i,50)=f(6,i,50)
end do
!bounce back on west boundary
do j=50,70
f(3,30,j)=f(1,30,j)
f(7,30,j)=f(5,30,j)
f(6,30,j)=f(8,30,j)
end do
!bounce back on east boundary
do j=50,70
f(1,40,j)=f(3,40,j)
f(5,40,j)=f(7,40,j)
f(8,40,j)=f(6,40,j)
end do
!!!!!自己随便加--边界 2
!bounce back on north boundary
do i=25,35
f(2,i,30)=f(4,i,30)
f(5,i,30)=f(7,i,30)
f(6,i,30)=f(8,i,30)
end do
!bounce back on south boundary
do i=25,35
f(4,i,15)=f(2,i,15)
f(7,i,15)=f(5,i,15)
f(8,i,15)=f(6,i,15)
end do
!bounce back on west boundary
do j=15,30
f(3,25,j)=f(1,25,j)
f(7,25,j)=f(5,25,j)
f(6,25,j)=f(8,25,j)
end do
!bounce back on east boundary
do j=15,30
f(1,35,j)=f(3,35,j)
f(5,35,j)=f(7,35,j)
f(8,35,j)=f(6,35,j)
end do
!!!!!自己随便加--边界 3
!bounce back on north boundary
do i=60,70
f(2,i,60)=f(4,i,60)
f(5,i,60)=f(7,i,60)
f(6,i,60)=f(8,i,60)
end do
!bounce back on south boundary
do i=60,70
f(4,i,40)=f(2,i,40)
f(7,i,40)=f(5,i,40)
f(8,i,40)=f(6,i,40)
end do
!bounce back on west boundary
do j=40,60
f(3,60,j)=f(1,60,j)
f(7,60,j)=f(5,60,j)
f(6,60,j)=f(8,60,j)
end do
!bounce back on east boundary
do j=40,60
f(1,70,j)=f(3,70,j)
f(5,70,j)=f(7,70,j)
f(8,70,j)=f(6,70,j)
end do
return
end
subroutine rhouv(f,rho,u,v,cx,cy,n,m)
real f(0:8,0:n,0:m),rho(0:n,0:m),u(0:n,0:m),v(0:n,0:m),cx(0:8),cy(0:8)
do j=0,m
do i=0,n
ssum=0.0
do k=0,8
ssum=ssum+f(k,i,j)
end do
rho(i,j)=ssum
end do
end do
do i=1,n
rho(i,m)=f(0,i,m)+f(1,i,m)+f(3,i,m)+2.*(f(2,i,m)+f(6,i,m)+f(5,i,m))
end do
DO i=1,n
DO j=1,m-1
usum=0.0
vsum=0.0
DO k=0,8
usum=usum+f(k,i,j)*cx(k)
vsum=vsum+f(k,i,j)*cy(k)
END DO
u(i,j)=usum/rho(i,j)
v(i,j)=vsum/rho(i,j)
END DO
END DO
do j=0,m
v(n,j)=0.0 !注意！设出口处没有竖直方向的速度了！只是个合理假设！
end do
do i=30,40 !注意！边界 1
do j=50,70
v(i,j)=0.0
u(i,j)=0.0
end do
end do
do i=25,35 !注意！边界 2
do j=15,30
v(i,j)=0.0
u(i,j)=0.0
end do
end do
do i=60,70 !注意！边界 3
do j=40,60
v(i,j)=0.0
u(i,j)=0.0
end do
end do
return
end
subroutine result(u,v,rho,uo,n,m)
real u(0:n,0:m),v(0:n,0:m)
real rho(0:n,0:m),strf(0:n,0:m)
open(5, file='streamf')
! streamfunction calculations
strf(0,0)=0.0
do i=1,n
rhoav=0.5*(rho(i-1,0)+rho(i,0))
if(i.ne.0) strf(i,0)=strf(i-1,0)-rhoav*0.5*(v(i-1,0)+v(i,0))
do j=1,m
rhom=0.5*(rho(i,j)+rho(i,j-1))
strf(i,j)=strf(i,j-1)+rhom*0.5*(u(i,j-1)+u(i,j))
end do
end do
! ———————————–
!write(2,*)"VARIABLES =X, Y, U, V, S"
!write(2,*)"ZONE","I=",n+1,"’J=",m+1,",","F=BLOCK"
do j=0,m
do i=0,n
write(2,*),i,j,u(i,j),v(i,j)
end do
end do
!do j=0,m
!write(2,*)(i,i=0,n)
!end do
!do j=0,m
!write(2,*)(j,i=0,n)
!end do
!do j=0,m
!write(2,*)(u(i,j),i=0,n)
!end do
!do j=0,m
!write(2,*)(v(i,j),i=0,n)
!end do
!do j=0,m
!write(2,*)(strf(i,j),i=0,n)
!end do
do j=0,m
write(3,*)j/float(m),u(n/2,j)/uo,u(n/4,j)/uo,u(3*n/4,j)/uo
end do
do i=0,n
write(4,*) i/float(n),v(i,m/2)/uo
end do
return
end
!============end of the program
后处理还是MATLAB代码：
%画出障碍流的流线图
close all;clc;clear all;
text= load('uvfield.txt');
x=text(:,1);
r=text(:,2);
dvx=text(:,3);
dvr=text(:,4);
%%
Fx = scatteredInterpolant(x,r,dvx); %对数据集执行插值
Fr = scatteredInterpolant(x,r,dvr);
%%
xx=linspace(min(x),max(x),2000); % xx= linspace(x1,x2,n) 生成 n 个点。这些点的间距为 (x2-x1)/(n-1)。 调节此处可以调整疏密度！
rr=linspace(min(r),max(r),2000);
[xgg,rgg]=meshgrid(xx,rr);
xstream = Fx(xgg,rgg);
ystream = Fr(xgg,rgg);
%%
scrsz = get(0,'ScreenSize'); %得到屏幕参数
figure1 = figure('Position',[0.05*scrsz(3) 0.05*scrsz(4) 0.9*scrsz(3) 0.9*scrsz(4)]); % 改变画图大小位置
% scrsz(1): 屏幕最左坐标；scrsz(2): 屏幕最下坐标
% scrsz(3): 屏幕宽(像素)；% scrsz(4): 屏幕高(像素)
[xs,rs] = meshgrid(x,r);
[dvxs,dvrs] = meshgrid(dvx,dvr);
quiver(x,r,dvx,dvr,'r');
numstream=600;
strx=randi([0,99],numstream,1); %randi([imin,imax],...) 返回一个在[imin,imax]范围内的伪随机整数
stry=randi([0,99],numstream,1); %r = randi(imax,[m,n])，返回一个在[1,imax]范围内的的m*n的伪随机整数矩阵 原始为stry=randi([0,12],numstream,1);
%值strx、y代表流线的起始位置
strx=[strx,strx];
stry=[stry,-stry];
h=streamline(xgg,rgg,xstream,ystream,strx,stry); %streamline(X,Y,Z,U,V,W,startx,starty,startz) 根据三维向量数据 U、V 和 W 绘制流线图。
%数组 X、Y 和 Z 用于定义 U、V 和 W 的坐标，它们必须是单调的，无需间距均匀。X、Y 和 Z 必须具有相同数量的元素，就像由 meshgrid 生成一样。
%startx、starty 和 startz 定义流线图的起始位置。
set(h,'LineWidth',0.5,'Color','k')
axis equal
axis tight
box on
%% 2nd
figure2 = figure('Position',[0.05*scrsz(3) 0.05*scrsz(4) 0.9*scrsz(3) 0.9*scrsz(4)]);
[xs,rs] = meshgrid(x,r);
[dvxs,dvrs] = meshgrid(dvx,dvr);
quiver(x,r,dvx,dvr,'r');
numstream=800;
strx=randi([2,99],numstream,1);
stry=randi([0,40],numstream,1); %此处原为[0,12]
strx=[strx,strx];
stry=[stry,-stry];
h=streamslice(xgg,rgg,xstream,ystream); %流线图
set(h,'LineWidth',0.7,'Color','b')
axis equal
axis tight
box on