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

合集下载

热传导模型方程matlab

热传导模型方程引言热传导模型方程是描述热传导过程的数学模型，它在物理学、工程学等领域中有着广泛的应用。

该方程可以用来求解物体内部的温度分布以及热量的传输过程。

本文将重点介绍热传导模型方程在Matlab中的实现方法及应用。

热传导模型方程热传导模型方程是一种偏微分方程，用于描述物体内部温度分布随时间的变化。

它的一般形式如下：∂u/∂t = κ∇²u其中，u是温度分布函数，t是时间，κ是热传导率，∇²是拉普拉斯算子。

离散化处理为了在计算机中求解热传导模型方程，我们需要对其进行离散化处理。

通常，我们将物体分成若干个小区域，每个小区域的温度可以看作是一个离散节点。

然后，我们利用数值差分方法将偏微分方程转化为差分方程。

最常用的数值差分方法是有限差分法，其中最简单的方法是二阶中心差分法。

二阶中心差分法二阶中心差分法是一种常用的数值差分方法，它能够较为准确地近似热传导模型方程。

该方法利用了节点的前后两个时间步长的温度值和空间步长的温度梯度。

具体来说，对于一个二维网格，节点(i, j)的二阶中心差分近似可以表示为：(∂u/∂t)(i, j) ≈ (u(i+1, j) - 2u(i, j) + u(i-1, j))/(Δx²) + (u(i, j+1) - 2u(i, j) + u(i, j-1))/(Δy²)其中，Δx和Δy分别是空间步长。

Matlab实现在Matlab中，我们可以通过编写代码来求解热传导模型方程。

下面是实现该方程的一个示例代码：% 定义初始条件和边界条件N = 50; % 网格节点数L = 1; % 系统尺寸T = 1; % 总时间dx = L/N;dt = dx^2/4;x = linspace(0, L, N+1);y = linspace(0, L, N+1);[X, Y] = meshgrid(x, y);u = sin(pi*X).*sin(pi*Y); % 初始温度分布% 迭代求解热传导模型方程for t = 0:dt:Tu_new = u + dt*(del2(u)/dx^2);u = u_new;end% 可视化温度分布surf(X, Y, u)xlabel('x')ylabel('y')zlabel('Temperature')在上面的代码中，我们首先定义了初始条件和边界条件，然后利用循环迭代的方式求解热传导模型方程。

传热学中几种常用的软件及数值解法的介绍

传热学中几种常用软件及数值解法的介绍一、常用软件介绍：1、FLUENT软件简介FLUENT软件是美国FLUENT公司开发的通用CFD流场计算分析软件，囊括了Fluent Dynamic International、比利时Polyflow和Fluent Dynamic International（FDI）的全部技术力量（前者是公认的粘弹性和聚合物流动模拟方面占领先地位的公司，而后者是基于有限元方法CFD软件方面领先的公司）。

FLUENT是用于计算流体流动和传热问题的程序。

由于采用了多种求解方法和多重网格加速收敛技术，因而FLUENT能达到最佳的收敛速度和求解精度。

灵活的非结构化网格和基于解的自适应网格技术及成熟的物理模型，使FLUENT在转捩与湍流、传热与相变、化学反应与燃烧、多相流、旋转机械、动/变形网格、噪声、材料加工、燃料电池等方面有广泛应用。

采用的数值解法有限体积法（Finite Volume Method）程序的结构FLUENT程序软件包由以下几个部分组成：(1)GAMBIT——用于建立几何结构和网格的生成。

(2)FLUENT——用于进行流动模拟计算的求解器。

(3)prePDF——用于模拟PDF燃烧过程。

(4)TGrid——用于从现有的边界网格生成体网格。

(5)Filters(Translators)—转换其他程序生成的网格，用于FLUENT计算。

FLUENT程序可以求解的问题(1)可压缩与不可压缩流动问题。

(2)稳态和瞬态流动问题。

(3)无黏流，层流及湍流问题。

(4)牛顿流体及非牛顿流体。

(5)对流换热问题(包括自然对流和混合对流)。

(6)导热与对流换热耦合问题。

(7)辐射换热。

(8)惯性坐标系和非惯性坐标系下的流动问题模拟。

(9)用Lagrangian轨道模型模拟稀疏相(颗粒，水滴，气泡等)。

(10)一维风扇、热交换器性能计算。

(11)两相流问题。

(12)复杂表面形状下的自由面流动问题。

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

(7)
Ti , j Ti , j 1 gx T0 2
图 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); tmp1=zeros(nx+2,ny+2); tmp2=zeros(nx+2,ny+2); c=zeros(nx+2,ny+2)+0.25; ut=zeros(nx+1,ny+2); vt=zeros(nx+2,ny+1); uu=zeros(nx+1,ny+1); vv=zeros(nx+1,ny+1); T=zeros(nx+1,ny+1); TT=zeros(nx+1,ny+1); c(2,3:ny)=1/3; c(3:nx,2)=1/3; c(2,2)=1/2; c(nx+1,2)=1/2; c(nx+1,3:ny)=1/3; c(3:nx,ny+1)=1/3; c(2,ny+1)=1/2; c(nx+1,ny+1)=1/2; beta=1./((Tw+Te)/2+273); % 边界条件

matlab 热阻热容模型

matlab 热阻热容模型Matlab热阻热容模型热阻热容模型是描述热系统的一种常用方法。

在许多热传导、热辐射和热对流等现象中，热阻和热容是重要的热学参数。

在本文中，我们将讨论如何使用Matlab来建立热阻热容模型，并通过一个实例来深入说明其应用。

1. 引言热阻热容模型是指一种将热系统中的热阻和热容组合在一起描述热传导或热辐射等现象的模型。

热阻表示物体传导热量的阻力，热容则表示物体吸收热量时的变热能力。

热阻和热容的比例决定了物体在热过程中的温度变化速率。

2. 热阻模型热阻模型用于描述物体传导热量的能力。

传统上，热阻可用下列方程表示：R = l / (k * A)其中，R为热阻，l为物体长度，k为热导率，A为物体截面积。

然而，在复杂的热传导现象中，热阻可能受到多种因素的影响，如边界条件、材料特性等。

因此，具体的热阻模型需要根据具体问题进行调整。

在Matlab 中，我们可以使用符号计算工具箱来建立热阻模型，根据具体问题设定不同的参数。

3. 热容模型热容模型用于描述物体变温能力。

传统的热容模型可以用下列方程表示：C = m * Cp其中，C为热容，m为物体质量，Cp为物体比热容。

同样，在复杂的热过程中，热容可能也受到其他因素的影响，如相变或化学反应等。

因此，在具体问题中，我们需要根据实际情况进行模型的调整。

在Matlab中，我们可以使用符号计算工具箱来建立热容模型，实现根据问题设定参数的灵活调整。

4. 热阻热容模型热阻热容模型是将热阻和热容组合在一起描述物体的热过程。

热阻热容模型可用下列方程表示：Q = R * T + C * (dT / dt)其中，Q为热量流，T为物体温度，t为时间，dT/dt为温度变化率。

在Matlab中，我们可以使用ODE（ordinary differential equation）求解器来求解上述方程。

具体而言，我们需要设定初始条件、边界条件和模型参数，然后使用ODE求解器对方程进行数值求解。

matlab传热计算程序

matlab传热计算程序
传热计算在工程学和科学领域中是一个重要的应用。

Matlab是一个功能强大的工程计算软件，可以用于传热计算。

在Matlab中，你可以使用各种方法来进行传热计算，比如有限元法、差分法、有限体积法等。

以下是一些常见的传热计算程序的示例：
1. 热传导方程求解，你可以编写一个Matlab程序来求解热传导方程，根据给定的边界条件和初始条件，使用差分法或有限元法来离散方程，并进行时间步进求解，得到温度场的分布。

2. 对流换热计算，对于流体内部的对流换热问题，你可以编写一个Matlab程序来求解Navier-Stokes方程和能量方程，结合有限体积法来进行流场和温度场的耦合求解。

3. 辐射换热计算，针对辐射换热问题，你可以编写一个Matlab程序来计算辐射传热，比如使用辐射传热方程和辐射传热模型，结合离散方法进行求解。

4. 传热系统优化，除了单一的传热计算，你还可以使用Matlab进行传热系统的优化设计，比如通过建立传热模型和耦合其
他工程模型，使用优化算法来寻找最优的传热系统设计参数。

总之，Matlab提供了丰富的工具和函数，可以用于传热计算的各个方面。

通过编写程序，你可以灵活地进行传热计算，并且可以根据具体的问题需求进行定制化的计算和分析。

希望这些信息对你有所帮助。

matlab求解一维带内热源传热问题

matlab 求解一维带内热源传热问题解一维带有内部热源的传热问题通常涉及到热传导方程的求解。

热传导方程描述了温度场随时间和空间的变化。

一维热传导方程通常写作：22()T T Q x t xα∂∂=+∂∂ 其中：• T 是温度，• t 是时间，• x 是空间坐标，• α 是热扩散系数，• Q(x) 是热源。

解这个方程需要适当的边界条件和初始条件。

为了简化问题，我们可以考虑一个稳态（0T t∂=∂）情况。

以下是使用 MATLAB 求解一维带有内部热源的传热问题的简单示例代码：% 参数设置L = 1; % 区域长度alpha = 0.01; % 热扩散系数Q = @(x) 1; % 内部热源% 空间离散化N = 100; % 离散网格数x = linspace(0, L, N);% 热传导方程T = zeros(1, N);T(1) = 0; % 初始条件T(N) = 100; % 边界条件% 离散格式求解dx = x(2) - x(1);dt = 0.01;num_steps = 1000;for step = 1:num_stepsfor i = 2:N-1T(i) = T(i) + alpha * dt / dx^2 * (T(i+1) - 2*T(i) + T(i-1)) + Q(x(i)) * dt;endend% 结果可视化plot(x, T);xlabel('空间坐标');ylabel('温度');title('一维带内部热源传热问题');请注意，这是一个简化的例子，具体的问题可能需要更多的考虑，例如更精确的数值方法、不同的边界条件和初始条件、更复杂的热源分布等。

这个示例主要用于演示MATLAB 中解决这类问题的基本方法。

方腔自然对流fortran编程

方腔自然对流fortran编程自然对流是物理学上的一个重要现象，在许多领域中都有应用。

其中，方腔自然对流是一个常见的现象，其研究涉及到流体力学、热传导等多个方面。

本文将介绍方腔自然对流的数值模拟方法，并使用Fortran语言进行编程实现。

一、数学模型我们考虑一个长方体的固定方腔，假设其一侧被加热，另一侧被冷却，四侧和顶部均为绝热。

我们假设流体是不可压缩、稳态、热传导系数恒定的流体。

此时，方腔内部温度的分布遵循如下的Navier-Stokes方程组和能量守恒方程：$ \frac{\partial u}{\partial x}+\frac{\partialv}{\partial y}=0$$ \frac{\partial u}{\partial t}+u\frac{\partialu}{\partial x}+v\frac{\partial u}{\partial y}=-\frac{1}{\rho}\frac{\partial p}{\partialx}+\nu(\frac{\partial ^2 u}{\partial x^2}+\frac{\partial ^2 u}{\partial y^2})$$ \frac{\partial v}{\partial t}+u\frac{\partialv}{\partial x}+v\frac{\partial v}{\partial y}=-\frac{1}{\rho}\frac{\partial p}{\partialy}+\nu(\frac{\partial ^2 v}{\partial x^2}+\frac{\partial ^2 v}{\partial y^2})$$ \frac{\partial T}{\partial t}+u\frac{\partialT}{\partial x}+v\frac{\partial T}{\partialy}=\frac{\alpha}{\rho c_p}(\frac{\partial ^2 T}{\partialx^2}+\frac{\partial ^2 T}{\partial y^2})$其中，$u$和$v$分别表示$x$和$y$方向上的速度分量，$p$表示压强，$T$表示温度，$\rho$表示密度，$\nu$表示运动粘度，$\alpha$表示热扩散系数，$c_p$表示等压比热容。

MATLAB语言在对流—扩散问题中的简单应用

MATLAB语言在对流—扩散问题中的简单应用
黄家友;陈燎原
【期刊名称】《煤矿机械》
【年(卷),期】2004(0)6
【摘要】介绍了在MATLAB软件的基础上 ,结合其计算、绘图等工具箱 ,重点分析了MATLAB语言在解决对流—扩散问题中的应用。

通过编程计算一个具体的二维稳态导热的例子 ,阐述了在对流—扩散问题中使用MATLAB可以提高工作效率 ,有利于方便快捷地解决问题。

【总页数】3页(P56-58)
【关键词】MATLAB;中心差分法;对流-扩散
【作者】黄家友;陈燎原
【作者单位】安徽理工大学
【正文语种】中文
【中图分类】TP30
【相关文献】
1.MATLAB在函数极限问题中的简单应用 [J], 张恩绮;陈文静;张旭;
2.MATLAB语言在计算流体力学中的简单应用 [J], 李志博;林麒;胡国清;许高攀;陈广文;刘文艳
3.非对称对流扩散问题有限元法的MATLAB实现 [J], 范馨月;杨一都
4.QUICK格式在对流占优扩散问题中的应用分析 [J], 赵文娟; 黄凌
5.浓度对流扩散方程并行计算与MATLAB高效实现方法 [J], 谢悦;林建国;芦静因版权原因，仅展示原文概要，查看原文内容请购买。

等温边界与对流边界换热完整matlab代码

eps=abs(T-temp)/t;%精度 temp=T; k=k+1;%迭代次数 end A=flipud(T);%将矩阵上下颠倒 pcolor(A);shading interp % figure(1),contour(A,1000);%等温线集合 figure,contour(A,[12,18,24]) %作出三条等温线 %计算内外壁导热量，偏差 dx=0.1;dy=0.1;lamda=0.53; temp1=zeros(13,1);temp2=zeros(13,1);temp3=zeros(13,1);temp4=zeros(13,1); for j=3:15 temp1(j-2)=(T(1,j)-T(2,j))*lamda*dx*1/dy; end for i=3:11 temp2(i-2)=(T(i,1)-T(i,2))*lamda*dy*1/dx; end heat_out=sum(temp1)+sum(temp2)+(T(2,1)-T(2,2))*lamda*dy*1/dx+(T(1,2)T(2,2))*lamda*dx*1/dy+(T(1,16)-T(2,16))*lamda*dx/2*1/dy+(T(12,1)T(12,2))*lamda*dy/2*1/dx; for j=6:15 temp3(j-5)=(T(5,j)-T(6,j))*lamda*i-5)=(T(i,5)-T(i,6))*lamda*dy*1/dx; end heat_in=sum(temp3)+sum(temp4)+(T(5,16)-T(6,16))*lamda*dx/2*1/dy+(T(12,5)T(12,6))*lamda*dy/2*1/dx; heat=(heat_in+heat_out)/2; e=abs(heat_in-heat_out)/heat; fprintf('外壁传热量=%2.3f\n',heat_out) fprintf('内壁传热量=%2.3f\n',heat_in) fprintf('平均传热量=%2.3f\n',heat) fprintf('相对偏差=%d\n',e) fprintf('墙角总传热量=%f\n',4*heat)

复杂方腔内自然对流换热数值模拟分析

复杂方腔内自然对流换热数值模拟分析近年来，随着科学技术的发展，计算机应用于工程力学领域中也取得了巨大进步。

多个数值模拟技术被广泛应用于各种工程领域，其中复杂方腔内自然对流换热数值模拟分析也得到了广泛关注。

复杂方腔内的自然对流换热是指当一个腔体内部的流体在温度不同的情况下发生自然运动时产生的换热作用，这种换热过程是不可逆的，而且其影响范围会深入到腔体内的每一个方向。

因此，对复杂方腔内的自然对流换热做出准确的数值模拟，不仅有助于研究腔体内的温度分布，而且对有关腔体温度控制系统的设计及改进也非常重要。

复杂方腔内的自然对流换热数值模拟分析也称为三维流体动力学模型，是基于泊松、拉格朗日等物理学方程构建的。

首先，根据复杂方腔内自然对流换热的物理过程，建立温度场和流场数学模型，包括温度-湿度条件和流体的动能方程以及传热方程。

接着，使用基于网格方法的求解工具，解决各种复杂的运算难题。

传统的复杂方腔内自然对流换热数值模拟分析方法，一般采用网格技术，分割区域，将区域网格化，将求解转化为一系列线性系统方程求解的问题，但这种方法的结果不够准确，而且计算耗时，并且存在速度与精度的矛盾。

最新的数值模拟方法，采用了新的“现代网格技术”，除了提高计算的准确性，更重要的是提高计算的速度，可以减少计算耗时，有效改善了复杂方腔内自然对流换热的数值模拟分析效果。

因此，复杂方腔内自然对流换热数值模拟分析也是一门非常有用的技术，在腔体内温度分布的研究以及腔体温度控制系统的设计上都有重要的应用价值，从而可以更好的改进腔体的热物理性能。

在腔体内的温度分布有助于了解物体的散热特性；腔体温度控制系统的设计为进一步改进腔体结构提供了可能性。

作为一项基础性研究，复杂方腔内自然对流换热数值模拟分析在腔体内温度分布研究、腔体温度控制系统设计以及腔体设计改进等方面，都有重要的实际意义。

但是，由于自由度较多、变量多及计算精度的要求，复杂方腔内的自然对流换热数值模拟分析仍然处于一个探索阶段。

数值传热学-CMAC算法-源程序代码

using System;using System.Collections.Generic;using ponentModel;using System.Data;using System.Drawing;using System.Text;using System.Windows.Forms;using System.IO;using System.Threading;using System.Collections;using ZedGraph;using ZedGraph.Web;using System.Drawing.Imaging;using Excel = Microsoft.Office.Interop.Excel;namespace CMAC{public partial class Form_CAMC : Form{public Form_CAMC(){InitializeComponent();}private void btn_js_Click(object sender, EventArgs e){int M = 41;//x方向网格数int N = 41;//y方向网格数//定义数组double[,] P = new double[M + 1, N + 1];//压力p值double[,] P1 = new double[M + 1, N + 1];//计算过程中的p值，用于过渡double[,] U = new double[M + 1, N + 1];//速度u,v的值（这里存储值的方式是以定位网格线来存储，与我们通常用的行和列刚好相反，请特别注意！）double[,] V = new double[M + 1, N + 1];double[,] U1 = new double[M + 1, N + 1];//计算过程中的u,v值，用于过渡double[,] V1 = new double[M + 1, N + 1];double[,] U11 = new double[M + 1, N + 1];//记录尖括号的u值double[,] V11 = new double[M + 1, N + 1];//记录尖括号的v值double[] XX = new double[M];//网格x坐标double[] YY = new double[N];//网格y坐标double[] XP = new double[M + 1];//压力p控制点x坐标double[] YU = new double[N + 1];//速度u控制点y坐标double[] XV = new double[M + 1];//速度v控制点x坐标double[] YV = new double[N + 1];//速度v控制点y坐标double[,] DUX1 = new double[M + 1, N + 1];//U的动量方程中对流相的第一项，x方向double[,] DUX2 = new double[M + 1, N + 1];//U的动量方程中对流相的第二项，y方向double[,] KSX1 = new double[M + 1, N + 1];//U动量方程扩散相的全部（两相一起）double[,] KSX2 = new double[M + 1, N + 1];double[,] DVY1 = new double[M + 1, N + 1];//V的动量方程中对流相的第一项，x方向double[,] DVY2 = new double[M + 1, N + 1];//V的动量方程中对流相的第二项，y方向double[,] KSY1 = new double[M + 1, N + 1];//V动量方程扩散项的第一项double[,] KSY2 = new double[M + 1, N + 1];//V动量方程扩散项的第二项double AP, AE, AW, AS, AN;//压力泊松方程系数double[,] B = new double[M + 1, N + 1];//压力方程源项double T, DT, DX, DY;//定义时间、时间步长、x方向步长、y方向步长int i, j, INT_P;int DA;//外部循环次数int K;//求解压力方程的循环次数double RE = 100.0;//雷诺数double X = 1.0;//方腔长度double Y = 1.0;//方腔高度double ERRP = 0.0;//压力P残差double ERRU = 0.0;//速度U残差double ERRV = 0.0;//速度V残差T = 0;//时间DT = 0.0001;//时间步长DX = X / (M - 1);//x方向步长DY = Y / (N - 1);//y方向步长#region//初始化变量for (i = 0; i < M + 1; i++){for (j = 0; j < N + 1; j++){P[i, j] = 0.0;}}for (i = 0; i < M + 1; i++){P1[i, j] = 0.0;}}for (i = 0; i < M + 1; i++){for (j = 0; j < N + 1; j++) {U[i, j] = 0.0;}}for (i = 0; i < M + 1; i++){for (j = 0; j < N + 1; j++) {V[i, j] = 0.0;}}for (i = 0; i < M + 1; i++){for (j = 0; j < N + 1; j++) {U1[i, j] = 0.0;}}for (i = 0; i < M + 1; i++){for (j = 0; j < N + 1; j++) {V1[i, j] = 0.0;}}for (i = 0; i < M + 1; i++){for (j = 0; j < N + 1; j++) {U11[i, j] = 0.0;}}for (i = 0; i < M + 1; i++){for (j = 0; j < N + 1; j++) {V11[i, j] = 0.0;}}for (j = 0; j < N + 1; j++){B[i, j] = 0.0;}}#endregion#region//坐标初始化//网格坐标初始化for (i = 0; i < M; i++){XX[i] = DX * i;//网格x坐标}for (j = 0; j < N; j++){YY[j] = DY * j;//网格y坐标}//压力p控制点坐标初始化XP[0] = 0.0;//压力边界点x坐标XP[M] = 1.0;for (i = 1; i < M; i++){XP[i] = (XX[i - 1] + XX[i]) / 2.0;//压力内点x坐标 }YP[0] = 0.0;//压力边界点y坐标YP[N] = 1.0;for (j = 1; j < N; j++){YP[j] = (YY[j - 1] + YY[j]) / 2.0;//压力内点y坐标 }//速度u控制点坐标初始化for (i = 0; i < M + 1; i++){XU[i] = XP[i];//速度u控制点x坐标}for (j = 0; j < N + 1; j++){YU[j] = YP[j];//速度u控制点y坐标}//速度v控制点坐标初始化XV[i] = XP[i];//速度v控制点x坐标}for (j = 0; j < N + 1; j++){YV[j] = YP[j];//速度v控制点y坐标}#endregion#region//边界条件for (j = 0; j < N + 1; j++){U[0, j] = 0.0;//为速度u的W、E边界赋初值。

4.7 方腔自然对流(8)

4.7方腔自然对流4.7.1物理模型一个边长为1m的正方形箱体，右墙温度为2000K，左墙温度为1000K，上下墙面绝热，重力方向向下，由于热重引起密度梯度，所以发展为浮力流。

箱中的介质具有吸收性和散射性的，，因此墙壁间的辐射交换因存在吸收被减弱，同时也因为介质的散射作用而增强了，所有墙壁被认为是黑体。

计算区域如图4-7-1所示。

图4-7-1 计算区域示意图4.7.2在Gambit中建立模型Step1：启动Gambit并选择求解器为Fluent5/6。

Step2：创建面操作：→→打开对话框，输入长度和宽度100，在Direction中选择XY Centered。

Step3：划分面网格1．操作：→→打开对话框如图4-7-2.（1）在Edges中选择正方形的四个边；（2）在Type中的下拉菜单中选择Successive Ratio；（3）选中Double sided前面的复选框；（4）输入Ratio1和Ratio2的值1.03；（5）点击Apply确认。

图4-7-2 网格划分设置对话框图4-7-3 划分面网格设置对话框2．操作：→→打开对话框如图4-7-3所示，Internal size=1，其它保留默认。

Step4：设置边界类型操作：→●在Name栏输入边界名称wall-1，将Type栏选为Wall，在Entity栏选取Edges，并选中方腔左边边线。

●在Name栏输入边界名称wall-2，将Type栏选为Wall，在Entity栏选取Edges，并选中方腔右边边线。

●在Name栏输入边界名称wall-3，将Type栏选为Wall，在Entity栏选取Edges，并选中方腔其它两条边线。

Step5：输出网格文件选中Export 2-D mesh 前面的复选框，输出网格文件。

4.7.3求解计算Step1：导入并检查网格1．读入网格文件操作：Fil e→Read→Case...找到文件后，单击OK按键确认。

1_数值传热-传热学作业-matlab

L P Rn5. 一厚度为 250mm 无限大平壁，其导热系数λ=43+0.08t w/(m ·k)，平壁一侧温度为 250℃，另一侧温度为 46℃，试用数值方法确定平壁内的温度分布，并确定通过该平壁的热流密度。

解：（1）建立模型本题属于非常物性，无热源的一维稳态导热过程将平壁沿厚度方向（x 方向）划分为 N 个均匀相等的间距。

节点布置如图所示12 3 4 N N+1本题给出 N=10。

（2）通过热流法建立离散方程a) 内节点离散方程对节点 P(i)所代表的微元体，在 x 方向上与节点 P 相邻的节点分别为 L(i-1)和 R(i+1)。

由于节点之间的间距很小，可以认为相邻节点间的温度分布是线性的。

节点 P 所代表的网格单元与它周围各网格单元之间的导热量可根据傅里叶定律直接写为:ΦLP = Z LP其中t i -1 - t i ⨯1∆xΦRP = Z RPt i +1 - t i ⨯1∆xZ RPZ LP= 43 + 0.08⨯ (t i + t i +1 )2 = 43 + 0.08⨯ (t i -1 + t i )2对节点 P(i)所代表的微元写热平衡式，即可得节点 P(i)温度的离散方程ΦLP + ΦRP = 0Z LPt i -1 - t i ⨯1+ Z ∆x RPt i +1 - t i ⨯1 = 0 ∆xt = Z LP t i -1 + Z LP t i +1 Z LP + Z RPb) 边界节点离散方程由于本题的壁面温度属于第一类边界条件，因此t 1 = 250 t n +1 = 46c) 热流密度的计算Z∑Z RP(t i +1 - t i )q = ∆t =i =166i46, t i = 次数 m, t k=k+1k>m不收停机（3）编写流程图计算机程序中使用的变量标识符： i 节点的坐标变量 t(i) 节点 i 的温度tt 前一次算出的节点温度 k 迭代次数计算机程序中输入数据： n 沿 x 方向的网格划分数 e 控制迭代过程终止的误差 m 允许的最大迭代次数开始输入 n, e, k, t 1, t n+1, t iN=10, t 1=250, t n+1= 1, e=0.01, k=10000迭代 k=0tt i =t iNoYesNoYes打印i , q 打印“敛”（4）编写程序本题使用matlab 软件，所编写的程序如下：clear;clc;t=ones(1,11); % 设定各项初始值q=ones(1,11);t(1)=250;t(11)=46;e=0.01;k=1;while 1 % 迭代程序tt=t;for i=2:10a=43+0.08*(t(i-1)+t(i))/2;b=43+0.08*(t(i)+t(i+1))/2;t(i)=(a*t(i-1)+b*t(i+1))/(a+b);endk=k+1;ttt=abs(t(5)-tt(5));if(ttt<e)break;endendfor i=2:11 % 计算热流密度a=43+0.08*(t(i-1)+t(i))/2;q(i)=q(i-1)+a*(t(i)-t(i-1));endktq=q(i)/0.25w=0:25:250;v=w/25+1;y=t(v); %绘制温度分布图plot(w,y),xlabel('位置（mm）'),ylabel('温度（℃）')legend('平壁内的温度分布',0),grid（5）计算结果迭代次数k=75平壁温度分布见下表热流密度q= -4.47E+04 w/m2（6）温度分布图∂⎬ ⎝ ⎭ i i -1i 10. 一砖墙厚 240mm ，内、外表面的表面传热系数分别为 6.0 w/(m 2·k)和 15 w/(m 2·k)，墙体材料的导热系数λ=0.43 w/(m ·k)，密度ρ=1668kg/m 3，比热容 c=0.75kJ/(kg ·K)，室内空气温度保持不变为 20℃，室外空气温度周期性变化，中午 12 点温度最高为 3℃，晚上 12 点温度最低为-9℃，试用数值计算方法确定内、外墙壁面温度在一天内的变化。

热传导问题的MATLAB数值计算

收稿日期:2002-05-09.作者简介:李灿(1968-),女,副教授;株州,湖南冶金职业技术学院冶金系(412000).¹Partial different ial equation toolbox user c s guide.T he M ath Works,Inc.,2000.热传导问题的M AT LAB 数值计算李灿湖南冶金职业技术学院冶金系高彦栋黄素逸华中科技大学能源与动力工程学院摘要:分析了应用M AT LAB 中PDE 工具箱解热传导问题的方法和步骤,编制了三个难以用解析方法求解的算例.采用有限元法求解导热偏微分方程,应用PDE 工具箱得到数值解.对适合圆柱坐标描述的问题,通过公式变化将其转换为能用PDE 工具箱求解的形式.算例表明,用M AT LAB 对复杂形状和有内热源的非稳态导热问题进行数值计算和图形处理是方便高效的.关键词:热传导;非稳态导热;M AT L AB;数值计算中图分类号:T K 124 文献标识码:A 文章编号:1671-4512(2002)09-0091-03许多工程问题需要确定物体内部的温度场或确定其内部温度到达某一限定值所需要的时间,因此研究导热问题特别是非稳态导热问题十分重要.目前非稳态导热问题的描述方程为多维非线性的偏微分方程,这些方程只在几何形状与边界条件都较简单的情况下才能求得理论解,而对于几何形状和边界条件复杂的情况多用数值解法,需借助于计算机将时间和空间坐标划分成数量巨大的网格才能得到较精确的数值解.本文应用M ATLAB 中PDE 工具箱,求解复杂边界条件的热传导问题.1 求解方法求解方法是基于数值解法中的有限元法[1],其基本原理是把计算区域划分成一系列的三角形单元,每个单元上取一个节点,选定一个形状函数(抛物线形或双曲线形),并通过单元中节点上的被求解变量值表示该函数.通过对控制方程作积分来获得离散方程.有限元法的最大优点是对不规则区域适应性好,故用MATLAB 方法求解的结果在边界上也较精确.对于适合圆柱坐标和球坐标描述的问题,通过对其热传导方程的变换,也能在MATLAB 中求解.应用MATLAB 的PDE (Partial Differential Equation)工具箱可以解如下四类偏微分方程¹-$#(c $u)+au =f ;d(5u /5t)-$#(c $u)+au =f ;d(52u/5t 2)-$#(c $u)+au =f ;-$#(c $u)+au =E du,(1)式中,u 为域8上的求解变量;E 为特征值;d,c,a,f 为常数或变量;t 为时间变量.前3个方程分别称为椭圆方程、抛物线方程和双曲线方程,第4个方程称为特征值方程.导热问题的通用微分方程可写成[2]Q c p (5u /5t )=$(K $u )+q v ,(2)式中,u 为求解变量,此处表示被求解物体内的温度;K 为导热系数;q v 为热源的发热率密度;Q 为密度;c p 为定压比热容.可以看出,式(1)和式(2)中的抛物线方程有着类似的形式.其中,求解变量为区域的温度,d 与Q c p ,c 与K ,f -au 与q v 可以一一对应.M ATLAB 中的PDE 工具箱定义了两类边界条件hu =r ;n #(c $u)+qu =g ,(3)式中,n 为垂直于边界的单位矢量;h ,r ,q 和g 为常量或与u 有关的变量.方程(3)中的第1个方程称为狄利克雷(Dirichlet)边界条件,第2个方程称为纽曼(Neumann)边界条件.可以看出,导热问题中的第一类边界与狄利克雷边界条件对应,第二类和第三类边界条件与纽曼边界条件对应.这些对应关系可以使用MATLAB 中的PDE 工具箱第30卷第9期华中科技大学学报(自然科学版) V ol.30 No.92002年 9月 J.Huazhong U niv.of Sci.&T ech.(Nature Science Edition)Sep.2002来求解.对一个导热问题的计算可以按图1的步骤进行.图1 M AT L AB 计算流程图2 算例2.1 三维非稳态无内热源的导热问题边长为0.5m,0.7m 和1.0m 的长方形钢锭,置于炉温u f =1200e 的加热炉内,计算5h 后钢锭的温度.已知钢锭的K =40.5W/(m #e ),A =0.722@10-5m 2/s ,u 0=25e ,钢锭与外界的对流换热系数h =348W/(m 2#e ).由对应关系可得d =Q c p =K /A ;c =K ;f =0;a =0,边界条件为纽曼边界条件,且钢锭的6个边界条件均相同,由对应关系有:q =h ; g =hf .求得5h 后钢锭内部的温度分布如图2,温度梯度如图3.两图还显示了有限元求解的网格,图3底平面的箭头方向为热流密度方向.图2 5h时刻钢锭的温度分布云图图3 5h 时刻钢锭的温度梯度云图如果导热体物性系数K 为温度u 的函数,只要写出K (u)的函数关系式,就可以得到解.2.2 有内热源的圆柱体非稳态导热问题有一半径为0.2m,长为3m 的圆柱形核电站用燃烧棒置于u f =100e 的水中,由于链式反应,棒内有恒定的产热率密度q v =20000W/m 3,计算10h 后燃烧棒内的温度分布.已知,燃烧棒的密度Q =7800kg /m 3,c p =500W #s/(kg #e ),K =40W/(m #e ),u 0=0e ,燃烧棒右端恒温t r =100e ,左端有一恒热流q l =5000W/m 2,燃烧棒外表面与外界水的对流换热系数h =50W/(m 2#e ).此问题宜采用圆柱坐标,由于燃烧棒内温度沿半径对称分布,因此可以转换为(r ,z )坐标的二维问题.将圆柱坐标内的热传导方程改写为Q c p r 5u 5t -55r K r 5u 5r -55z K r 5u 5z =q v r ,(4)以使其形式与式(1)拟合.式(4)与式(1)中的抛物线方程对比可以得出:d =Q c p r ;c =K r ;a =0;f =q v r ,式中,z 对应第一个坐标方向(在直角坐标中为x 方向);r 对应第二个坐标方向(在直角坐标中为y ).燃烧棒左端的边界条件为:n #(K $u )=q r ,为纽曼边界条件,由对应关系得:q =0; g =q l r ,燃烧棒右端为狄利克雷边界条件u =100.燃烧棒上(外)边界条件n #(K $u)=h(u f -u)为纽曼边界条件,由对应关系得q =hr ;g =hu f r.解析域的下边界为棒的中心,其边界条件为n #(K $t)=0,也为纽曼边界条件.q =0,则把q 和g 都设为0即可.求得10h 时刻燃烧棒内部的温度分布如图4所示,热流密度分布如图5所示.图4 10h 时刻燃烧棒的温度分布云图92 华中科技大学学报(自然科学版) 第30卷图5 10h 时刻燃烧棒的热流密度云图如果内热源是时间或空间的函数,写出函数关系式,也可以得到解.2.3 复杂边界的热传导问题考虑这样一个问题:一个正方形内嵌一菱形,其中正方形区域的密度为2W/m 2,导热系数为10W/(cm #e ),菱形的密度为1W/m 2,导热系数为5W/(cm #e ),并有内热源的发热率密度为10W/m 2,两个区域的定压比热容均为0.1J/(kg #e ).初始温度为0e .计算0.1s 之后的等温图如图6所示,箭头所指为热流方向,热流密度图如图7所示.由此可见,应用MATLAB 可以方便快捷地解出复杂几何形状和复杂边界条件的非稳态问题.并且其强大的图形可视化功能使计算结果形象、直观而便于理解.图6 0.1s时刻等温图图7 0.1s 时刻热流密度云图参考文献[1]陶文铨.数值传热学(第2版).西安:西安交通大学出版社,2001.[2]程俊国,张洪济.高等传热学.重庆:重庆大学出版社,1990.Numerical simulation of problems in heat conduction using MATLABL i Can Gao Yandong H uang S uyiAbstract:The method and steps for finding the solutions for problems in heat conduction w ith the PDE toolbox in M ATLAB are described.T hree ex amples difficult to resolve w ith the analy tical method are g iven.The partial differential equation (PDE)for heat conduction is solved w ith the finite element method and the PDE toolbox is adopted to obtain the num erical simulation.Problems suitable for description w ith cylindrical coordinates are transformed into forms that are capable of solution w ith the PDE toolbox through formula variation.Examples of calculation show that M ATLAB is convenient and highly efficient for numerical simulation and graphic processing of com plex g eometry and non -steady -state heat conduction problems w ith internal thermal source.Key words:heat conduction;non -steady state heat conduction;MATLAB;numerical simulationLi Can Assoc.Prof.;Dept.of M etallurgy,Hu c nan Metallurg y Professional and Technical College,Zhuzhou 412000,Hu c nan,China.93第9期李灿等:热传导问题的M AT LAB 数值计算。

matlab计算流体力学管道传热

matlab计算流体力学管道传热下载提示：该文档是本店铺精心编制而成的，希望大家下载后，能够帮助大家解决实际问题。

文档下载后可定制修改，请根据实际需要进行调整和使用，谢谢!本店铺为大家提供各种类型的实用资料，如教育随笔、日记赏析、句子摘抄、古诗大全、经典美文、话题作文、工作总结、词语解析、文案摘录、其他资料等等，想了解不同资料格式和写法，敬请关注!Download tips: This document is carefully compiled by this editor. I hope that after you download it, it can help you solve practical problems. The document can be customized and modified after downloading, please adjust and use it according to actual needs, thank you! In addition, this shop provides you with various types of practical materials, such as educational essays, diary appreciation, sentence excerpts, ancient poems, classic articles, topic composition, work summary, word parsing, copy excerpts, other materials and so on, want to know different data formats and writing methods, please pay attention!在工程领域中，流体力学管道传热一直是一个备受关注的重要问题。

MATLAB语言在对流_扩散问题中的简单应用

文章编号 :100320794 (2004) 0620056203
MATLAB 语言在对流 — 扩散问题中的简单应用
黄家友 , 陈燎原
( 安徽理工大学 , 安徽淮南 232001)
摘要 : 介绍了在 MAT LAB 软件的基础上 , 结合其计算、绘图等工具箱 , 重点分析了 M一个具体的二维稳态导热的例子 ,阐述了在对流 — 扩散问题中使用 MAT LAB 可以提高工作效率 ,有利于方便快捷地解决问题。关键词 : MAT LAB ; 中心差分法 ; 对流 — 扩散中图号 : TP30
xi yi zi 4 结语
5 <) 5 <) 5 <) 5 <) - (Γ A ] + [ (Γ A - (Γ A ] 5x e 5x w 5x n 5x s ( 3) + SΔV
数即可完成 ,其编程是非常简单和快捷的。
参考文献 :
[1 ] 卢佐潮 ,等 . 计算机辅助机械设计 [ M] . 广州 : 华南理工大学出版
( Dw Aw Fs Fw
( 15)
2
Aw ) <P + ( De A e + Fn
Fe
2
A e ) <P + ( Ds A s
控制方程为 5 5T 5 (Γ 5 T ) ( 19) (Γ ) + =0 5x 5x 5y 5y 依据上面给出的初始条件用数值积分方法可求出精确解为 ∞ n π 200 1 - ( - 1) n T ( x , y ) = 300 + sin x× π n∑ n L =1 π n sh ( y ) L ( 20) π ) sh ( n 对平板进行网格划分后 , 根据上面的分析求出

哈尔滨工程大学传热学大作业数值计算matlab程序内容

传热学作业数值计算数值计算matlab程序内容：>> tw1=10; % 赋初值tw2=20;c=1.5;p2=20;p1=c*p2;L2=40;L1=c*L2;deltaX=L2/p2;a=p2+1;b=p1+1;ti=ones(a,b)*5;m1=ones(a,b);m1(a,2:b-1)=zeros(1,b-2);m1(2:a,1)=zeros(a-1,1);m1(2:a,b)=zeros(a-1,1);m1(1,:)=ones(1,b)*2;k=0;max1=1.0;tn=ti;while(max1>1e-6)max1=0;k=k+1;for i=1:1:afor j=1:1:bm=m1(i,j);n=ti(i,j);switch mcase 0tn(i,j)=tw1;case 1tn(i,j)=0.25*(tn(i,j+1)+tn(i,j-1)+tn(i+1,j)+tn(i-1,j));case 2tn(i,j)=tw1+tw2*sin(pi*(j-1)/(b-1));ender=abs(tn(i,j)-n);if er>max1max1=er;endendendti=tn;endktimax1t2=ones(a,b); %求解析温度场for i=a:-1:1for j=1:1:by=deltaX*(a-i);x=deltaX*(j-1);t2(i,j)=tw1+tw2*sin(pi*x/L1)*(sinh(pi*y/L1))/(sinh(pi*L2/L1));endendt2迭代次数k =706数值解温度场ti数值解每次迭代的最大误差max1 =9.8531e-07解析温度场t2取第11行的解析解和数值解的点曲线为第11行的解析解的直线，散点为其数值解的点第11行的误差=[数值解(11行) –解析解(11行)]/解析解数值温度场图像解析温度场图像数值解与解析解的误差数值计算matlab程序内容：>> tw1=10;tw2=20;c=1.5;p2=20;p1=c*p2;L2=20;deltaX=L2/p2;L1=c*L2;a=p2+1;b=p1+1;ti=ones(a,b)*5;m1=ones(a,b);m1(a,2:b-1)=zeros(1,b-2);m1(2:a,1)=zeros(a-1,1);m1(2:a,b)=zeros(a-1,1);m1(1,:)=ones(1,b)*2;k=0;max1=1.0;tn=ti;while(max1>1e-6)max1=0;k=k+1;for i=1:1:afor j=1:1:bm=m1(i,j);n=ti(i,j);switch mcase 0tn(i,j)=tw1;case 1tn(i,j)=0.25*(tn(i,j+1)+tn(i,j-1)+tn(i+1,j)+tn(i-1,j));case 2tn(i,j)=tw2;ender=abs(tn(i,j)-n);if er>max1max1=er;endendendti=tn;endktimax1tx=ones(a,b);for i=1:1:afor j=1:1:by=(a-i)*deltaX;x=(j-1)*deltaX;m=sym('m');g=(((-1)^(m+1)+1)/m)*sin(m*pi*x/L1)*sinh(m*pi*y/L1)/sinh(m*pi*L2/L1); h=symsum(g,m,1,100);tx(i,j)=2*h*(tw2-tw1)/pi+tw1;endendtx迭代次数k = 695数值解温度场ti数值解每次迭代的最大误差max1 =9.8243e-07解析温度场tx =取第11行的解析解和数值解的点曲线为第11行的解析解的直线，散点为其数值解的点第11行的误差=[数值解(11行) –解析解(11行)]/解析解数值温度场ti图像：解析温度场tx图像：数值解与解析解的误差数值计算matlab程序内容：>> t0=90;tf=10;L=10;c=0.25;p2=20;p1=p2/c;B=c*L;d=0.5*B;h=10;a=p2+1;b=p1+1;deltaX=B/p2;lambda=160;Bi=h*deltaX/lambda;ti=ones(a,b)*10;m1=ones(a,b)*3;m1(2:a-1,1)=zeros(a-2,1);m1(a,2:b-1)=ones(1,b-2);m1(1,2:b-1)=ones(1,b-2)*6;m1(2:a-1,b)=ones(a-2,1)*2;m1(1,b)=ones(1,1)*4;m1(a,b)=ones(1,1)*5;m1(1,1)=7;m1(a,1)=8;tn=ti;max1=1.0;k=0;while ( max1>1e-6)k=k+1;max1=0;for i=1:1:afor j=1:1:bm=m1(i,j);n=tn(i,j);switch mcase 0tn(i,j)=t0;case 1tn(i,j)=(2*tn(i-1,j)+tn(i,j-1)+tn(i,j+1)-4*tf)/(4+2*Bi)+tf;case 2tn(i,j)=(2*tn(i,j-1)+tn(i-1,j)+tn(i+1,j)-4*tf)/(4+2*Bi)+tf;case 3tn(i,j)=0.25*(tn(i,j-1)+tn(i,j+1)+tn(i-1,j)+tn(i+1,j));case 4tn(i,j)=(tn(i,j-1)+tn(i+1,j)-2*tf)/(2*Bi+2)+tf;case 5tn(i,j)=(tn(i,j-1)+tn(i-1,j)-2*tf)/(2*Bi+2)+tf;case 6tn(i,j)=(2*tn(i+1,j)+tn(i,j-1)+tn(i,j+1)-4*tf)/(4+2*Bi)+tf;case 7tn(i,j)=t0;case 8tn(i,j)=t0;ender=abs(tn(i,j)-n);if er>max1max1=er;endendendti=tn;endktita=ones(a,b);Bi1=h*d/lambda;sbi=sqrt(Bi1);for i=1:1:afor j=1:1:bif i>(a+1)/2y=-(i-(a+1)/2)*deltaX;else y=((a+1)/2-i)*deltaX;endx=deltaX*(j-1);ta(i,j)=(cosh(sbi*(L-x)/d)+sbi*sinh(sbi*(L-x)/d))*(t0-tf)/(cosh(sbi*L/d)+sbi*sinh(sbi*L/d))+tf;endendta迭代次数k =1461数值解温度场ti解析温度场ta取第11行的解析解和数值解的点曲线为第11行的解析解的直线，散点为其数值解的点第11行的误差=[数值解(11行) –解析解(11行)]/解析解图像如下数值温度场图像解析温度场图像数值解与解析解的误差数值计算matlab程序内容：>> tw=10;L2=15;c=0.75;L1=L2/c;p2=24 ;p1=p2/c;deltaX=2*L2/p2;a=p2+1;b=p1+1;lambda=16;qv0=24;ti=ones(a,b)*5;m1=ones(a,b);m1(1,:)=zeros(1,b);m1(2:a,b)=zeros(a-1,1);m1(2:a,1)=zeros(a-1,1);m1(a,2:b-1)=zeros(1,b-2);tn=ti;max1=1.0;k=0;while(max1>1e-6)max1=0;k=k+1;for i=1:1:afor j=1:1:bm=m1(i,j);n=tn(i,j);switch mcase 0tn(i,j)=tw;case 1tn(i,j)=0.25*(tn(i-1,j)+tn(i+1,j)+tn(i,j-1)+tn(i,j+1)+qv0*(deltaX^2)/lambda);ender=abs(tn(i,j)-n);if er>max1max1=er;endendendti=tn;endtik;tx=ones(a,b);for i=1:1:afor j=1:1:bif i>(a+1)/2y=-(i-(a+1)/2)*deltaX;elsey=((a+1)/2-i)*deltaX;endif j>(b+1)/2x=(j-(b+1)/2)*deltaX;elsex=-((b+1)/2-j)*deltaX;endm=sym('m');xi=(2*m-1)*pi/2;g=((-1)^m)/(xi^3)*(cosh(xi*y/L1)/cosh(xi*L2/L1))*cos(xi*x/L1); h=symsum(g,m,1,100);tx(i,j)=2*qv0*L1^2/lambda*h+qv0*(L1^2-x^2)/(2*lambda)+tw; endendtx数值温度场ti解析温度场tx取第13行的解析解和数值解的点曲线为第13行的解析解的直线，散点为其数值解的点第13行的误差=[数值解(13行) –解析解(13行)]/解析解数值温度场图像解析温度场图像数值解与解析解的误差。

对流方程各种格式代码matlab

对流方程各种格式代码matlab对流方程——偏微分方程的数值解法用迎风格式解对流方程function u = peYF(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);if a>0for j=1:(n+M)u0(j) = IniU(minx+(j-M-1)*h);endelsefor j=1:(n+M)u0(j) = IniU(minx+(j-1)*h);endendu1 = u0;for k=1:Mif a>0for i=(k+1):n+Mu1(i) = -dt*a*(u0(i)-u0(i-1))/h+u0(i);endelsefor i=1:n+M-ku1(i) = -dt*a*(u0(i+1)-u0(i))/h+u0(i);endendu0 = u1;endif a>0u = u1((M+1):M+n);elseu = u1(1:n);endformat long;用拉克斯-弗里德里希斯格式解对流方程function u = peHypbLax(a,dt,n,minx,maxx,M) format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M)u0(j) = IniU(minx+(j-M-1)*h);endu1 = u0;for k=1:Mfor i=k+1:n+2*M-ku1(i) = -dt*a*(u0(i+1)-u0(i-1))/h/2+(u0(i+1)+u0(i-1))/2;endu0 = u1;endu = u1((M+1):(M+n));format short;用拉克斯-温德洛夫格式解对流方程function u = peLaxW(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M)u0(j) = IniU(minx+(j-M-1)*h);endu1 = u0;for k=1:Mfor i=k+1:n+2*M-ku1(i) = dt*dt*a*a*(u0(i+1)-2*u0(i)+u0(i-1))/2/h/h - ...dt*a*(u0(i+1)-u0(i-1))/h/2+u0(i);endu0 = u1;endu = u1((M+1):(M+n));format short;用比姆-沃明格式解对流方程function u = peBW(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M)u0(j) = IniU(minx+(j-2*M-1)*h);endu1 = u0;for k=1:Mfor i=2*k+1:n+2*Mu1(i) = u0(i)-dt*a*(u0(i)-u0(i-1))/h-a*dt*(1-a*dt/h)* ...(u0(i)-2*u0(i-1)+u0(i-2))/2/h;endu0 = u1;endu = u1((2*M+1):(2*M+n));format short;用Richtmyer多步格式解对流方程function u = peRich(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+4*M)u0(j) = IniU(minx+(j-2*M-1)*h);endu1 = u0;for k=1:Mfor i=2*k+1:n+4*M-2*ktmpU1 = -dt*a*(u0(i+2)-u0(i))/h/4+(u0(i+2)+u0(i))/2;tmpU2 = -dt*a*(u0(i)-u0(i-2))/h/4+(u0(i)+u0(i-2))/2;u1(i) = -dt*a*(tmpU1-tmpU2)/h/2+u0(i);endu0 = u1;endu = u1((2*M+1):(2*M+n));format short;用拉克斯-温德洛夫多步格式解对流方程function u = peMLW(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M)u0(j) = IniU(minx+(j-M-1)*h);endu1 = u0;for k=1:Mfor i=k+1:n+2*M-ktmpU1 = -dt*a*(u0(i+1)-u0(i))/h/2+(u0(i+1)+u0(i))/2;tmpU2 = -dt*a*(u0(i)-u0(i-1))/h/2+(u0(i)+u0(i-1))/2;u1(i) = -dt*a*(tmpU1-tmpU2)/h+u0(i);endu0 = u1;endu = u1((M+1):(M+n));format short;用MacCormack多步格式解对流方程function u = peMC(a,dt,n,minx,maxx,M)format long;h = (maxx-minx)/(n-1);for j=1:(n+2*M)u0(j) = IniU(minx+(j-M-1)*h);endu1 = u0;for k=1:Mfor i=k+1:n+2*M-ktmpU1 = -dt*a*(u0(i+1)-u0(i))/h+u0(i);tmpU2 = -dt*a*(u0(i)-u0(i-1))/h+u0(i-1);u1(i) = -dt*a*(tmpU1-tmpU2)/h/2+(u0(i)+tmpU1)/2;endu0 = u1;endu = u1((M+1):(M+n));format short;用拉克斯-弗里德里希斯格式解二维对流方程的初值问题function u = pe2LF(a,b,dt,nx,minx,maxx,ny,miny,maxy,M)%啦-佛format long;hx = (maxx-minx)/(nx-1);hy = (maxy-miny)/(ny-1);for i=1:nx+2*Mfor j=1:(ny+2*M)u0(i,j) = Ini2U(minx+(i-M-1)*hx,miny+(j-M-1)*hy);endendu1 = u0;for k=1:Mfor i=k+1:nx+2*M-kfor j=k+1:ny+2*M-ku1(i,j) = (u0(i+1,j)+u0(i-1,j)+u0(i,j+1)+u0(i,j-1))/4 ...-a*dt*(u0(i+1,j)-u0(i-1,j))/2/hx ...-b*dt*(u0(i,j+1)-u0(i,j-1))/2/hy;endendu0 = u1;endu = u1((M+1):(M+nx),(M+1):(M+ny));format short;用拉克斯-弗里德里希斯格式解二维对流方程的初值问题function u = pe2FL(a,b,dt,nx,minx,maxx,ny,miny,maxy,M)%近似分裂format long;hx = (maxx-minx)/(nx-1);hy = (maxy-miny)/(ny-1);for i=1:nx+4*Mfor j=1:(ny+4*M)u0(i,j) = Ini2U(minx+(i-2*M-1)*hx,miny+(j-2*M-1)*hy);endendu1 = u0;for k=1:Mfor i=2*k+1:nx+4*M-2*kfor j=2*k-1:ny+4*M-2*k+2tmpU(i,j) = u0(i,j) - a*dt*(u0(i+1,j)-u0(i-1,j))/2/hx + ...(a*dt/hx)^2*(u0(i+1,j)-2*u0(i,j)+u0(i-1,j))/2;endendfor i=2*k+1:nx+4*M-2*kfor j=2*k+1:nx+4*M-2*ku1(i,j) = tmpU(i,j) - b*dt*(tmpU(i,j+1)-tmpU(i,j-1))/2/hy + ...(b*dt/hy)^2*(tmpU(i,j+1)-2*tmpU(i,j)+tmpU(i,j-1))/2;endendu0 = u1;endu = u1((2*M+1):(2*M+nx),(2*M+1):(2*M+ny));format short;。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

数学模型：
控制方程：使用 Boussinesq 近似来考虑温差所引起的自然对流，在 N-S 方程中添加一个体积力项。
u 1 uu p 2u g 是运动粘度，是流体的体胀系数。能量守恒方程如下：
(6)
vit, j vit, j 1 vit, j vit, j 1 2 2
2
2

vit+1, j 2vit, j vit1, j h2
vit, j 1 2vit, j vit, j 1 h2
(7)
Ti , j Ti , j 1 gx T0 2
图 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); tmp1=zeros(nx+2,ny+2); tmp2=zeros(nx+2,ny+2); c=zeros(nx+2,ny+2)+0.25; ut=zeros(nx+1,ny+2); vt=zeros(nx+2,ny+1); uu=zeros(nx+1,ny+1); vv=zeros(nx+1,ny+1); T=zeros(nx+1,ny+1); TT=zeros(nx+1,ny+1); c(2,3:ny)=1/3; c(3:nx,2)=1/3; c(2,2)=1/2; c(nx+1,2)=1/2; c(nx+1,3:ny)=1/3; c(3:nx,ny+1)=1/3; c(2,ny+1)=1/2; c(nx+1,ny+1)=1/2; beta=1./((Tw+Te)/2+273); % 边界条件
T u T 2T t
(2)
其中，是热扩散系数。边界条件：由无滑移边界条件得，四个壁面上的速度均为零，即： un us uw ue 0 。在热壁面上 T 100 ，在冷壁面上 T 0 ，在上下绝热壁面上处
T 0。 y
数值处理：
区域离散化如图 2 所示。对于动量守恒方程，在不考虑压力的情况下先计算出一个临时速度
%==========================================================================
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 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 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); v(1,1:ny+1)=2*vw-v(2,1:ny+1); u(1:nx+1,ny+2)=2*un-u(1:nx+1,ny+1); v(nx+2,1:ny+1)=2*ve-v(nx+1,1:ny+1); % 求修正速度 % 求压力 % 临时速度v % 临时速度u
u u n un un 2un g T T0 dt
(3)
用 SOR（超松弛迭代）法求压力场
1 1 pi, j1 c pi 1, j pi 1, j pi , j 1 pi , j 1
h ui , j ui1, j vi, j vi1, j 1 pi, j dt
%========================================================================== 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));
(4)
根据求得的临时速度，在考虑压力的情况下，得修正速度
1 uin, j ui , j
dt

1 pi 1, j pi , j dx
2
(5)
方程(3)对应的离散方程如下：
t t t t t t t t 1 ui , j ui 1, j ui , j ui 1, j vi , j vi 1, j ui , j 1 ui , j dt h 2 2 2 2 t t t t t t t t t t ui +1, j 2ui , j ui 1, j ui , j 1 2ui , j ui , j 1 vi , j 1 vi 1, j 1 ui , j ui , j 1 2 2 2 h h2 Ti , j Ti 1, j gx T0 2 t t t t t t t t vi , j vi , j 1 ui , j ui , j 1 vi , j vi 1, j ui 1, j uit1, j 1 vit, j vit1, j dt h 2 2 2 2 t t uit, j ui , j 2
%========================================================================== for i=2:nx for j=2:ny T(i,j)=TT(i,j)+dt*( -(0.25/h)*( (uu(i,j)+uu(i+1,j))*(TT(i,j)+... TT(i+1,j))-(uu(i-1,j)+uu(i,j))*(TT(i-1,j)+TT(i,j))+... (vv(i,j)+vv(i,j+1))*(TT(i,j)+TT(i,j+1))-(vv(i,j)+... vv(i,j-1))*(TT(i,j)+TT(i,j-1)) )+(alpha/h^2)*(TT(i+1,j)+... TT(i-1,j)+TT(i,j+1)+TT(i,j-1)-4*TT(i,j)) ); end end for i=2:nx T(i,1)=T(i,2); T(i,ny+1)=T(i,ny); end TT=T; %========================================================================== end %========================================================================== for i=1:nx+1 for j=1:ny+1 velocity(i,j)=sqrt(uu(i,j)*uu(i,j)+vv(i,j)*vv(i,j)); end end %========================================================================== startx=h:10*h:0.5; starty=zeros(size(startx))+0.5; streamline(xh,yh,flipud(rot90(uu)),flipud(rot90(vv)),startx,starty); clf(figure) quiver(xh,yh,flipud(rot90(uu)),flipud(rot90(vv)),'r'),title('Velocity'),axis equal,axis([0 L 0 H]); clf(figure) contour(xh,yh,flipud(rot90(T))),title('Temperature map'),axis equal,axis([0 L 0 H]); clf(figure) [X,Y,Z]=griddata(xh,yh,flipud(rot90(T)),linspace(0,1,nx+1)',linspace(0,1,ny+1),'v4'); pcolor(X,Y,Z);shading interp;title('Temperature contour'); colorbar,colormap(jet); %========================================================================== %求温度场

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

热传导模型方程matlab

传热学中几种常用的软件及数值解法的介绍

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

matlab 热阻热容模型

matlab传热计算程序

matlab求解一维带内热源传热问题

方腔自然对流fortran编程

最新利用matlab程序解决热传导问题

MATLAB语言在对流—扩散问题中的简单应用

等温边界与对流边界换热完整matlab代码

复杂方腔内自然对流换热数值模拟分析

数值传热学-CMAC算法-源程序代码

4.7 方腔自然对流(8)

1_数值传热-传热学作业-matlab

热传导问题的MATLAB数值计算

matlab计算流体力学管道传热

MATLAB语言在对流_扩散问题中的简单应用

哈尔滨工程大学传热学大作业数值计算matlab程序内容

对流方程各种格式代码matlab