当前位置：首页 > news >正文

MATLAB 二维方腔自然对流 SIMPLE 算法

news 2026/6/30 6:37:43

主脚本 main.m

clear;clc;close all;%% 参数Ra=1e5;Pr=0.71;nx=64;ny=nx;% 网格L=1;dx=L/nx;dy=dx;dt=0.01;alpha=0.1;% 亚松弛maxIt=2000;tol=1e-5;%% 场量（交错网格）u=zeros(ny+2,nx+1);% u(i,j) 位于 (i-0.5,j)v=zeros(ny+1,nx+2);% v(i,j) 位于 (i,j-0.5)p=zeros(ny+2,nx+2);% p(i,j) 位于单元中心T=zeros(ny+2,nx+2)+0.5;% 初始 0.5Th=1;Tc=0;% 左热右冷%% 系数beta=1;g=Ra/(Pr*L^3)*(Th-Tc);% Boussinesq 浮力项%% 主循环fork=1:maxIt% 1) 预测速度（动量方程）[u,v]=momentum(u,v,p,T,dx,dy,dt,Pr,g,alpha);% 2) 压力修正 (SIMPLE)[p,u,v]=pressureCorrection(u,v,p,dx,dy,dt,alpha);% 3) 能量方程T=energy(u,v,T,dx,dy,dt,Pr,alpha);% 4) 收敛监测res=max([max(abs(u(:))),max(abs(v(:))),max(abs(T(:)))]);ifmod(k,100)==0,fprintf('it=%4d res=%.3e\n',k,res);endifres<tol,break;endend%% 后处理figure;contourf(T',20,'LineColor','none');colorbar;axis equal;title(['T 场 Ra='num2str(Ra)]);

动量方程 momentum.m

function[u,v]=momentum(u,v,p,T,dx,dy,dt,Pr,g,alpha)[ny,nx]=size(p);ue=u;ve=v;% 系数（中心差分 + 一阶迎风）forj=2:nx-1fori=2:ny-1% u 方程Fe=0.5*(u(i,j)+u(i+1,j))/dx;Fw=0.5*(u(i,j)+u(i-1,j))/dx;Fn=0.5*(v(i,j)+v(i,j+1))/dy;Fs=0.5*(v(i,j)+v(i,j-1))/dy;De=1/dx^2;Dw=De;Dn=1/dy^2;Ds=Dn;aE=De-max(Fe,0);aW=Dw-max(Fw,0);aN=Dn-max(Fn,0);aS=Ds-max(Fs,0);aP=aE+aW+aN+aS+(Fe-Fw+Fn-Fs)+1/dt;b=(p(i,j)-p(i,j+1))/dx+0.5*(T(i,j)+T(i,j+1))*g+u(i,j)/dt;ue(i,j)=u(i,j)+alpha/aP*b;% v 方程（同理，略）...endendu=ue;v=ve;end

压力修正 pressureCorrection.m

function[p,u,v]=pressureCorrection(u,v,p,dx,dy,dt,alpha)[ny,nx]=size(p);uStar=u;vStar=v;% 构造压力修正方程系数forj=2:nx-1fori=2:ny-1ae=dy/dx;aw=ae;an=dx/dy;as=an;ap=ae+aw+an+as;d=(uStar(i,j-1)-uStar(i,j))*dy+(vStar(i-1,j)-vStar(i,j))*dx;% TDMA 解 p'...endend% 修正速度 & 压力u(2:ny-1,2:nx-1)=uStar(2:ny-1,2:nx-1)-(p(2:ny-1,3:nx)-p(2:ny-1,2:nx-1))/dx*alpha;v(2:ny-1,2:nx-1)=vStar(2:ny-1,2:nx-1)-(p(3:ny,2:nx-1)-p(2:ny-1,2:nx-1))/dy*alpha;p=p+alpha*p;end

能量方程 energy.m

functionT=energy(u,v,T,dx,dy,dt,Pr,alpha)[ny,nx]=size(T);Te=T;forj=2:nx-1fori=2:ny-1% 对流项一阶迎风 + 扩散项中心差分Fe=max(u(i,j),0)*T(i,j)+max(-u(i,j),0)*T(i,j+1);Fw=max(u(i,j-1),0)*T(i,j-1)+max(-u(i,j-1),0)*T(i,j);Fn=max(v(i,j),0)*T(i,j)+max(-v(i,j),0)*T(i+1,j);Fs=max(v(i-1,j),0)*T(i-1,j)+max(-v(i-1,j),0)*T(i,j);De=1/dx^2/Pr;Dw=De;Dn=1/dy^2/Pr;Ds=Dn;aE=De;aW=Dw;aN=Dn;aS=Ds;aP=aE+aW+aN+aS+1/dt;b=(Fe-Fw+Fn-Fs)+T(i,j)/dt;Te(i,j)=T(i,j)+alpha*b/aP;endendT=Te;% 边界：左热右冷，上下绝热T(:,1)=1;T(:,end)=0;T(1,:)=T(2,:);T(end,:)=T(end-1,:);end