方腔内自然对流MATLAB程序数值传热学.pdf

自然 对流传热 问题描述 ： 一个 二维矩形腔体， 高 H1 ，宽 L1 ， 物理 模型如 图 1 绝热 壁面 绝热 壁面 图 1 物理 模型 上下 壁面为绝热壁面， 左壁面 为热 壁面 100T ，右壁面 为冷壁面 0T 腔体 内为 初始温度 0 0T 的 空气，计算由此引起的自然对流传热 数学 模型： 控制 方程 ： 使用 Boussinesq近似来 考虑温差所引起的自然对流 ，在 N-S方程中添加一个 体积力 项  201 p g T Tt         u u u u(1) 其中， 0T 是 初始温度 ，  是运动粘度 ，  是 流体的体胀系数 能量 守恒 方程 如下： 2T TTt      u(2) 其中，  是热 扩散系数 边界 条件： 由 无滑移边界 条件 得，四个壁面上的速度均为零，即： 0un us uw ue    在 热壁面上 100T ，在冷壁面上 0T ，在 上下 绝热壁面上 处 0Ty  数值 处理： 区域 离散化 如 图 2所示 对于 动量守恒方程 ，在 不考虑压力 的 情况下先计算出一个临时速度  20n n n nuu g T Tdt         u u u (3) 用 SOR（超松弛迭代 ）法 求压力 场   1 1 1, 1 , 1 , , 1 , 1, 1 , , 1 , ,1i j i j i j i j i ji j i j i j i j i jp c p p p ph u u v v pdt                    (4) 根据 求得的临 时 速度，在考虑压力的情况下， 得 修正速度 冷 壁 面 T=0 热 壁 面 T=100 1, , 1 , ,1ni j i j i j i ju u p pdt dx (5) 方程 (3)对应的离散 方程如 下： 22, , , 1 , , 1 , , 1 , , 1 ,, 1 1 , 1 , , 1 + 1 , , 1 ,212 2 2 2222t t t t t t t t t t ti j i j i j i j i j i j i j i j i j i jt t t t t t ti j i j i j i j i j i j i ju u u u u u v v u ud t hv v u u u u uh                                                                , 1 , , 12, 1 ,022t t ti j i j i ji j i ju u uhTTg x T        (6) , , , , 1 , 1 , 1 , 1 , 1 , 1 ,22, , 1 , , 1 + 1 , , 1 ,212 2 2 2222t t t t t t t t t t ti j i j i j i j i j i j i j i j i j i jt t t t t t ti j i j i j i j i j i j i jv v u u v v u u v vdt hv v v v v v vh                                                                 , 1 , , 12, , 1022t t ti j i j i ji j i jv v vhTTgx T        (7) 图 2 交错 网格 数值 结果： 速度场 图 3 速度 场 温度 等值线和温度云图 图 4 温度等值线 图 5 温度 云图 由 图 3-5可知 ， 程序 ： %========================================================================== H=1;L=1;nx=32;ny=32;h=1/nx;dt=0.002; gx=0;gy=-100;T0=0.0;kviscosity=0.01; alpha=0.01;Wallhot=100;Wallcold=0;nstep=2000; maxit=100;beta0=1.2; %========================================================================== u=zeros(nx+1,ny+2); v=zeros(nx+2,ny+1); p=zeros(nx+2,ny+2); ut=zeros(nx+1,ny+2); vt=zeros(nx+2,ny+1); tmp1=zeros(nx+2,ny+2); uu=zeros(nx+1,ny+1); vv=zeros(nx+1,ny+1); tmp2=zeros(nx+2,ny+2); T=zeros(nx+1,ny+1); TT=zeros(nx+1,ny+1); c=zeros(nx+2,ny+2)+0.25; %========================================================================== c(2,3:ny)=1/3; c(nx+1,3:ny)=1/3; c(3:nx,2)=1/3; c(3:nx,ny+1)=1/3; c(2,2)=1/2; c(2,ny+1)=1/2; c(nx+1,2)=1/2; c(nx+1,ny+1)=1/2; un=0.0; us=0.0; vw=0.0; ve=0.0; Tw=100.0; Te=0.0; T(1,1:ny+1)=100; TT(1,1:ny+1)=100; beta=1./((Tw+Te)/2+273); %========================================================================== u(1:nx+1,1)=2*us-u(1:nx+1,2); u(1:nx+1,ny+2)=2*un-u(1:nx+1,ny+1); % 边界 条件 v(1,1:ny+1)=2*vw-v(2,1:ny+1); v(nx+2,1:ny+1)=2*ve-v(nx+1,1:ny+1); %========================================================================== for i=1:nx+1 xh(i)=h*(i-1); end for j=1:ny+1 yh(j)=h*(j-1); end %========================================================================== for is=1:nstep for i=2:nx for j=2:ny+1 % 临时速度 u ut(i,j)=u(i,j)+dt*( -(0.25/h)*( (u(i,j)+u(i+1,j))^2-(u(i,j)+... u(i-1,j))^2+(v(i,j)+v(i+1,j))*(u(i,j+1)+u(i,j))-... (v(i,j-1)+v(i+1,j-1))*(u(i,j)+u(i,j-1)) )+... (kviscosity/h^2)*(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1)-... 4*u(i,j))+beta*gx*(0.5*(T(i,j-1)+T(i,j))-T0) ); end end for i=2:nx+1 for j=2:ny % 临时速度 v vt(i,j)=v(i,j)+dt*( -(0.25/h)*( (u(i,j)+u(i,j+1))*(v(i,j)+... v(i+1,j))-(u(i-1,j)+u(i-1,j+1))*(v(i-1,j)+v(i,j))+... (v(i,j)+v(i,j+1))^2-(v(i,j)+v(i,j-1))^2 )+... (kviscosity/h^2)*(v(i+1,j)+v(i-1,j)+v(i,j+1)+v(i,j-1)-... 4*v(i,j))+beta*gy*(0.5*(T(i-1,j)+T(i,j))-T0) ); end end for it=1:maxit % 求 压力 for i=2:nx+1 for j=2:ny+1 p(i,j)=beta0*c(i,j)*(p(i+1,j)+p(i-1,j)+p(i,j+1)+p(i,j-1)-... (h/dt)*(ut(i,j)-ut(i-1,j)+vt(i,j)-vt(i,j-1)))+... (1-beta0)*p(i,j); end end end %========================================================================== u(2:nx,2:ny+1)=ut(2:nx,2:ny+1)-(dt/h)*(p(3:nx+1,2:ny+1)-p(2:nx,2:ny+1)); % 求 修 正 速 度 v(2:nx+1,2:ny)=vt(2:nx+1,2:ny)-(dt/h)*(p(2:nx+1,3:ny+1)-p(2:nx+1,2:ny)); %========================================================================== u(1:nx+1,1)=2*us-u(1:nx+1,2); u(1:nx+1,ny+2)=2*un-u(1:nx+1,ny+1); v(1,1:ny+1)=2*vw-v(2,1:ny+1); v(nx+2,1:ny+1)=2*ve-v(nx+1,1:ny+1); %========================================================================== uu(1:nx+1,1:ny+1)=0.5*(u(1:nx+1,1:ny+1)+u(1:nx+1,2:ny+2)); vv(1:nx+1,1:ny+1)=0.5*(v(1:nx+1,1:ny+1)+v(2:nx+2,1:ny+1)); %================================。