矩形障碍算例（附Fortran计算代码及MATLAB后处理代码）

! 我又来啦，做了个小算例，在100*100个格子的空腔中加了几个矩形障碍

!上代码！西边界为水流进入区

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

矩形障碍算例（附Fortran计算代码及MATLAB后处理代码）相关推荐

matlab sift代码解读,MATLAB SIFT 代码
[实例简介] matlab 实现的 sift 变换的代码,包含整个过程的详细步骤. [实例截图] [核心代码] sift-0.9.0 ├── data │ ├── img3.jpg │ ├─ ...
matlab画基尼系数,Dagum基尼系数分解的MATLAB程序代码(更新)
很多同学在初次使用Dagum基尼系数法做研究时苦于找不到计算软件,手动计算又耗时费力.去年写论文时,参考Dagum(1997)论文中提供的算例,自己编写了一个MATLAB版的计算程序,在这里分享给大家 ...
动态规划原理介绍（附7个算例，有代码讲解）
动态规划思想动态规划(Dynamicprogramming)是一种通过把原问题分解为相对简单的子问题的方式求解复杂问题的方法. 动态规划常常适用于有重叠子问题和最优子结构性质的问题,动态规划方法所耗 ...
Fortran 算例
初学Fortran 相信大家都希望有一个简易的算例上手看看整个代码编写.运行的模式以下就是用Fortran编写的计算圆柱体表面积的一个代码实现功能:输入半径和高,自动计算圆柱体表面积 progra ...
有限元方法入门：有限元方法简单的二维算例（矩形剖分）
#有限元方法简单的二维算例(矩形剖分) 算例描述我们对下述椭圆边值问题 \label{eq1}{−Δu=fu∣∂Ω=0\left \{ \begin{aligned} & -\Delta u ...
中文文本纠错算例实现（有算例完整代码）
概述文本纠错又称为拼写错误或者拼写检查,由于纯文本往往来源于手打或者OCR识别,很可能存在一些错误,因此此技术也是一大关键的文本预处理过程,一般存在两大纠错类型. 1拼写错误第一种是Non-wor ...
微电网日前优化调度。算例有代码（1）
个人电气博文传送门学好电气全靠它,个人电气博文目录(持续更新中-) 符号说明问题1 求解经济性评估方案: 若微网中蓄电池不作用,且微网与电网交换功率无约束,在无可再生能源和可再生能源全额利用两 ...
matlab中右三角形方向,《有限元基础教程》_【MATLAB算例】4.7.1(2) 基于3节点三角形单元的矩形薄板分析(Triangle2D3Node)...
[MATLAB 算例]4.7.1(2) 基于3节点三角形单元的矩形薄板分析(T riangle2D3Node) 如图4-20所示为一矩形薄平板,在右端部受集中力100 000F N =作用,材料常数为 ...
python-pandapower电力系统短路电流计算（算例1：短路计算讲解））
本系列讲解电力系统潮流计算和最优潮流等,用pandapower求解,语言为python. 专栏订阅后,可以查看该专栏所有文章.希望学完这个专栏后完全掌握,建议从第一个算例看起. 建议先看前面几章的潮流 ...

矩形障碍算例（附Fortran计算代码及MATLAB后处理代码）

矩形障碍算例（附Fortran计算代码及MATLAB后处理代码）相关推荐

最新文章

热门文章