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

最新推荐文章于 2025-12-19 15:31:04 发布

原创最新推荐文章于 2025-12-19 15:31:04 发布 · 247 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#算法 #matlab #人工智能

主脚本 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.5
Th  = 1; Tc = 0;          % 左热右冷

%% 系数
beta = 1; g = Ra/(Pr*L^3)*(Th-Tc); % Boussinesq 浮力项

%% 主循环
for k = 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(:)))]);
    if mod(k,100)==0, fprintf('it=%4d  res=%.3e\n',k,res); end
    if res<tol, break; end
end

%% 后处理
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;
% 系数（中心差分 + 一阶迎风）
for j=2:nx-1
    for i=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 方程（同理，略）
        ...
    end
end
u = 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;
% 构造压力修正方程系数
for j=2:nx-1
    for i=2:ny-1
        ae = 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'
        ...
    end
end
% 修正速度 & 压力
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

function T = energy(u,v,T,dx,dy,dt,Pr,alpha)
[ny,nx] = size(T); Te = T;
for j=2:nx-1
    for i=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;
    end
end
T = Te;
% 边界：左热右冷，上下绝热
T(:,1)   = 1; T(:,end) = 0;
T(1,:)   = T(2,:); T(end,:) = T(end-1,:);
end