热传导问题的MATLAB数值计算

合集下载

一维介质中的热传导问题 卡尔曼滤波 matlab

一维介质中的热传导问题 卡尔曼滤波 matlab

一维介质中的热传导问题一、概述热传导是物理学中的一个重要问题,特别是对于介质的热传导问题更是如此。

一维介质中的热传导问题是指介质在一维空间内热量的传导过程。

这一问题不仅在物理学中具有重要性,而且在工程领域中也有着广泛的应用。

在实际工程中,我们常常需要对介质中的热传导问题进行分析和研究,以便更好地设计和优化热传导设备,提高能源利用效率。

二、热传导方程介质中的热传导过程可以用热传导方程来描述。

一维情况下,热传导方程可以写为:其中,u(x, t)为介质中的温度分布,k为介质的热导率,c为介质的比热容,ρ为介质的密度,t为时间,x为空间坐标。

三、数值模拟对于介质中的热传导问题,我们常常需要进行数值模拟来解决热传导方程。

数值模拟可以采用有限差分法、有限元法等数值方法来进行。

在进行数值模拟时,我们通常需要借助计算机软件来进行计算,其中Matlab是一种非常实用的数学建模和仿真软件,特别适用于解决热传导问题。

四、卡尔曼滤波卡尔曼滤波是一种最优状态估计算法,可以用于对系统的状态进行预测和估计。

在介质中的热传导问题中,我们可以利用卡尔曼滤波算法来对系统的温度状态进行估计,从而更好地理解和分析热传导过程。

五、Matlab仿真在研究介质中的热传导问题时,我们可以利用Matlab软件进行仿真计算。

通过编写Matlab程序,我们可以对介质中的热传导过程进行模拟,并得到系统的温度分布。

我们也可以借助Matlab提供的工具,如ODE求解器等,对热传导方程进行数值求解,得到系统的温度变化规律。

六、结论介质中的热传导问题是一个具有重要意义的物理问题,对其进行深入的研究不仅有助于提高工程设备的效率,而且可以推动物理学领域的发展。

卡尔曼滤波和Matlab仿真技术的应用为介质中的热传导问题研究提供了新的方法和手段,可以更好地帮助我们理解和解决这一重要问题。

希望未来能够有更多的研究者投入到介质中的热传导问题的研究中,共同推动科学技术的进步。

基于matlab的热力系统计算

基于matlab的热力系统计算

基于Matlab的热力系统计算简介热力系统是指由热能传递和转换组成的系统,常用于工业生产和能源领域。

Matlab作为一种强大的数学计算和模拟工具,提供了丰富的函数和工具包,能够方便地进行热力系统的计算和分析。

本文将介绍如何使用Matlab进行热力系统的计算,包括热能传递、热力转换和效率分析等方面。

热力系统模型热力系统可以用数学模型来描述,常用的模型包括热传导模型、热对流模型和热辐射模型。

下面将分别介绍这几种模型的计算方法。

热传导模型热传导模型用于描述热能在固体或液体中的传递过程。

在Matlab中,可以使用pdepe函数求解定态或非定态的热传导问题。

该函数可以根据给定的边界条件和初始条件,计算出固体或液体中的温度分布。

下面是一个示例代码:function [c, f, s] = pdefun(x, t, u, DuDx)k = 0.5; % 热传导系数c = 1; % 热容量f = k * DuDx;s = 0;该代码定义了一个热传导方程,并通过pdepe函数求解。

其中,c、f和s分别表示热容量、热传导通量和内部热源,x 和t分别表示空间和时间变量,u表示温度分布,DuDx表示温度梯度。

热对流模型热对流模型用于描述热能在流体中的传递过程。

在Matlab 中,可以使用convectiveFlux函数计算热对流传热通量。

该函数可以根据给定的流体流速、温度差和传热系数,计算出热对流传热量。

下面是一个示例代码:function [Q] = convectiveFlux(u, v, T1, T2, h) Q = h * u * v * (T2 - T1);该代码定义了一个热对流模型,并通过convectiveFlux 函数计算传热量。

其中,u和v分别表示流体的速度分量,T1和T2表示流体中的温度,h表示传热系数。

热辐射模型热辐射模型用于描述热能通过辐射的方式传递的过程。

在Matlab中,可以使用radiativeFlux函数计算热辐射传热通量。

热传导方程的求解

热传导方程的求解

应用物理软件训练前言MATLAB 是美国MathWorks公司出品的商业数学软件,用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境,主要包括MATLAB和Simulink两大部分。

MATLAB是矩阵实验室(Matrix Laboratory)的简称,和Mathematica、Maple并称为三大数学软件。

它在数学类科技应用软件中在数值计算方面首屈一指。

MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等,主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。

本部分主要介绍如何根据所学热传导方程的理论知识进行MATLAB数值实现可视化。

本部分主要介绍如何根据所学热传导方程的理论知识进行MATLAB数值实现可视化。

本部分主要介绍如何根据所学热传导方程的理论知识进行MATLAB数值实现可视化。

本部分主要介绍如何根据所学热传导方程的理论知识进行MATLAB数值实现可视化。

题目:热传导方程的求解目录一、参数说明 (1)二、基本原理 (1)三、MATLAB程序流程图 (3)四、源程序 (3)五、程序调试情况 (6)六、仿真中遇到的问题 (9)七、结束语 (9)八、参考文献 (10)一、参数说明U=zeros(21,101) 返回一个21*101的零矩阵x=linspace(0,1,100);将变量设成列向量meshz(u)绘制矩阵打的三维图axis([0 21 0 1]);横坐标从0到21,纵坐标从0到1eps是MATLAB默认的最小浮点数精度[X,Y]=pol2cart(R,TH);效果和上一句相同waterfall(RR,TT,wn)瀑布图二、基本原理1、一维热传导问题(1)无限长细杆的热传导定解问题利用傅里叶变换求得问题的解是:取得初始温度分布如下这是在区间0到1之间的高度为1的一个矩形脉冲,于是得(2)有限长细杆的热传导定解问题其中20x 0≤≤,即L=20,取a=10且得的解是(3)非齐次方程定解问题是解析解是其中2、二维热传导问题 定解问题Ut=k^2(Uxx+Uyy) (b y a ≤≤≤≤0,x 0) U(x=0,y,t)=0, u(x=a,y,t)=b3sinyπμ U (x,y=0,t )=0, u(x,y=b,t)=axx ππμcos a 3sinU (x,y,t=0)=03、三维热传导问题球体内的热传导令u=w+Uo,则w 的定解问题是 Wt=w ∆w W (r=ro )=0 W(t=to)=uo-Uo解为rorn enruo Uo w or t a n nnπππsin)1()(22222/1-∞=∑--=r 为空间变量,并用x ,y 表示。

matlab ode解一维热传导偏微分方程

matlab ode解一维热传导偏微分方程

matlab ode解一维热传导偏微分方程一维热传导偏微分方程是在众多领域中经常出现的一个方程,如何用数值方法求解这个方程一直是数学科学家们研究的一个方向。

在这篇文章中,我们将围绕Matlab的Ode求解器,介绍如何使用Matlab 来解决一维热传导偏微分方程。

首先,我们要了解一维热传导方程的形式。

一维热传导方程如下所示:ut = kuxx其中,u表示温度,t表示时间,k是热传导系数,x是空间坐标。

该方程描述了温度随时间和空间的变化情况。

接下来,我们将使用Matlab Ode求解器来解决这个方程。

一个很重要的问题是,我们需要将一维热传导方程转换为一个ODE系统。

这可以通过离散化方法来实现。

我们可以将空间x离散为N个点,用差分来近似求解uxx,进而得到一个差分方程组。

例如,我们可以使用中心差分来近似求解uxx,得到如下方程组:u0 = uN = 0ui,j+1 –ui,j = (kΔt/Δx^2)*(ui+1,j –2ui,j + ui-1,j)其中,ui,j 表示在时间j和位置 i 处的温度,Δx是网格宽度,Δt是时间步长。

现在,我们已经将一维热传导方程转换为一个差分方程组,可以使用Matlab的Ode求解器来解决。

首先,我们需要将差分方程组转换为ODE向量形式。

将所有的ui,j都展开成一个向量u,然后将等式转化为一个向量形式。

我们可以将每一个方程表示为:ui,j+1 – ui,j = F(ui,j)其中,F(ui,j) 表示u的时间导数在i, j的位置。

接下来,我们需要将这个ODE系统输入到Matlab Ode求解器中。

可以使用ODE45或ODE23等求解器解决。

首先,需要定义一个包含所有ODE的函数,该函数接受一个向量u和时间t作为输入,并返回u 的时间导数。

然后,需要指定初始条件 u0 和时间范围。

最后,调用ode45或ode23等求解器,将ODE函数传递给求解器,并得到解。

在得到解之后,可以将解绘制成一维热传导的温度分布图。

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

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

哈佛大学能源与环境学院课程作业报告作业名称:传热学大作业——利用matlab程序解决热传导问题院系:能源与环境学院专业:建筑环境与设备工程学号:姓名:盖茨比2015年6月8日一、题目及要求1.原始题目及要求2.各节点的离散化的代数方程3.源程序4.不同初值时的收敛快慢5.上下边界的热流量(λ=1W/(m℃))6.计算结果的等温线图7.计算小结题目:已知条件如下图所示:二、各节点的离散化的代数方程各温度节点的代数方程ta=(300+b+e)/4 ; tb=(200+a+c+f)/4; tc=(200+b+d+g)/4; td=(2*c+200+h)/4 te=(100+a+f+i)/4; tf=(b+e+g+j)/4; tg=(c+f+h+k)/4 ; th=(2*g+d+l)/4ti=(100+e+m+j)/4; tj=(f+i+k+n)/4; tk=(g+j+l+o)/4; tl=(2*k+h+q)/4tm=(2*i+300+n)/24; tn=(2*j+m+p+200)/24; to=(2*k+p+n+200)/24; tp=(l+o+100)/12 三、源程序【G-S迭代程序】【方法一】函数文件为:function [y,n]=gauseidel(A,b,x0,eps)D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);G=(D-L)\U;f=(D-L)\b;y=G*x0+f;n=1;while norm(y-x0)>=epsx0=y;y=G*x0+f;n=n+1;end命令文件为:A=[4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0;-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0;0,-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0;0,0,-2,4,0,0,0,-1,0,0,0,0,0,0,0,0;-1,0,0,0,4,-1,0,0,-1,0,0,0,0,0,0,0;0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0,0;0,0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0;0,0,0,-1,0,0,-2,4,0,0,0,-1,0,0,0,0;0,0,0,0,-1,0,-1,0,4,0,0,0,-1,0,0,0;0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0,0;0,0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0;0,0,0,0,0,0,0,-1,0,0,-2,4,0,0,0,-1;0,0,0,0,0,0,0,0,-2,0,0,0,24,-1,0,0;0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1,0;0,0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1;0,0,0,0,0,0,0,0,0,0,0,-1,0,0,-1,12];b=[300,200,200,200,100,0,0,0,100,0,0,0,300,200,200,100]';[x,n]=gauseidel(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',1.0e-6) xx=1:1:4;yy=xx;[X,Y]=meshgrid(xx,yy);Z=reshape(x,4,4);Z=Z'contour(X,Y,Z,30)Z =139.6088 150.3312 153.0517 153.5639108.1040 108.6641 108.3119 108.1523 84.1429 67.9096 63.3793 62.4214 20.1557 15.4521 14.8744 14.7746 【方法2】>> t=zeros(5,5);t(1,1)=100;t(1,2)=100;t(1,3)=100;t(1,4)=100;t(1,5)=100;t(2,1)=200;t(3,1)=200;t(4,1)=200;t(5,1)=200;for i=1:10t(2,2)=(300+t(3,2)+t(2,3))/4 ;t(3,2)=(200+t(2,2)+t(4,2)+t(3,3))/4;t(4,2)=(200+t(3,2)+t(5,2)+t(4,3))/4;t(5,2)=(2*t(4,2)+200+t(5,3))/4;t(2,3)=(100+t(2,2)+t(3,3)+t(2,4))/4;t(3,3)=(t(3,2)+t(2,3)+t(4,3)+t(3,4))/4; t(4,3)=(t(4,2)+t(3,3)+t(5,3)+t(4,4))/4; t(5,3)=(2*t(4,3)+t(5,2)+t(5,4))/4;t(2,4)=(100+t(2,3)+t(2,5)+t(3,4))/4;t(3,4)=(t(3,3)+t(2,4)+t(4,4)+t(3,5))/4;t(4,4)=(t(4,3)+t(4,5)+t(3,4)+t(5,4))/4;t(5,4)=(2*t(4,4)+t(5,3)+t(5,5))/4;t(2,5)=(2*t(2,4)+300+t(3,5))/24;t(3,5)=(2*t(3,4)+t(2,5)+t(4,5)+200)/24;t(4,5)=(2*t(4,4)+t(3,5)+t(5,5)+200)/24;t(5,5)=(t(5,4)+t(4,5)+100)/12;t'endcontour(t',50);ans =100.0000 200.0000 200.0000 200.0000 200.0000 100.0000 136.8905 146.9674 149.8587 150.7444 100.0000 102.3012 103.2880 103.8632 104.3496 100.0000 70.6264 61.9465 59.8018 59.6008 100.0000 19.0033 14.8903 14.5393 14.5117【Jacobi迭代程序】函数文件为:function [y,n]=jacobi(A,b,x0,eps)D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=D\(L+U);f=D\b;y=B*x0+f;n=1;while norm(y-x0)>=epsx0=y;y=B*x0+f;n=n+1;end命令文件为:A=[4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0;-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0; 0,-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0; 0,0,-2,4,0,0,0,-1,0,0,0,0,0,0,0,0;-1,0,0,0,4,-1,0,0,-1,0,0,0,0,0,0,0; 0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0,0; 0,0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0;0,0,0,-1,0,0,-2,4,0,0,0,-1,0,0,0,0;0,0,0,0,-1,0,-1,0,4,0,0,0,-1,0,0,0;0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0,0;0,0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0;0,0,0,0,0,0,0,-1,0,0,-2,4,0,0,0,-1;0,0,0,0,0,0,0,0,-2,0,0,0,24,-1,0,0;0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1,0;0,0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1;0,0,0,0,0,0,0,0,0,0,0,-1,0,0,-1,12];b=[300,200,200,200,100,0,0,0,100,0,0,0,300,200,200,100]'; [x,n]=jacobi(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',1.0e-6); xx=1:1:4;yy=xx;[X,Y]=meshgrid(xx,yy);Z=reshape(x,4,4);Z=Z'contour(X,Y,Z,30)n =97Z =139.6088 150.3312 153.0517 153.5639108.1040 108.6641 108.3119 108.152384.1429 67.9096 63.3793 62.421420.1557 15.4521 14.8744 14.7746四、不同初值时的收敛快慢1、[方法1]在Gauss 迭代和Jacobi 迭代中,本程序应用的收敛条件均为norm(y-x0)>=eps ,即使前后所求误差达到e 的-6次方时,跳出循环得出结果。

Matlab在热传递学课程中的应用

Matlab在热传递学课程中的应用

Matlab在热传递学课程中的应用热传递学是研究热能传递和传导的学科,广泛应用于工程、物理和环境领域。

在热传递学课程中,Matlab是一个常用的工具,可以帮助学生理解和分析热传递过程。

下面将介绍在热传递学课程中使用Matlab的步骤和应用。

第一步是建立热传递模型。

在研究热传递过程时,我们需要建立相应的数学模型。

可以使用Matlab来编写这些模型,并通过求解数学方程来分析热传递现象。

例如,我们可以使用Matlab编写热传导方程,并求解得到温度分布。

第二步是处理边界条件。

在热传递过程中,边界条件对结果有着重要的影响。

例如,我们可以设置材料的初始温度、表面的热通量或边界温度等。

Matlab提供了丰富的边界条件处理函数和图形界面,使得处理边界条件变得更加简便。

第三步是求解热传递问题。

在建立了合适的模型和边界条件后,我们可以使用Matlab的数值求解方法来求解热传递问题。

Matlab提供了许多数值求解算法,如有限差分法和有限元法,可以帮助我们得到准确的结果。

通过对求解结果的分析和可视化,我们可以更好地理解热传递过程。

第四步是进行参数敏感性分析。

在研究热传递过程时,我们通常需要考虑不同的参数对结果的影响。

Matlab提供了参数敏感性分析的工具,可以帮助我们理解不同参数对热传递问题的影响程度。

通过参数敏感性分析,我们可以选择最优的参数组合,并优化热传递系统的设计。

第五步是进行热传递实验和数据处理。

除了数值分析,实验也是研究热传递的重要手段。

Matlab 可以辅助我们进行热传递实验的数据处理和分析。

通过编写Matlab程序,我们可以快速地进行数据处理、绘图和拟合曲线,从而更好地理解实验数据和验证理论模型。

综上所述,Matlab在热传递学课程中具有广泛的应用。

它可以帮助学生建立热传递模型,处理边界条件,求解热传递问题,进行参数敏感性分析,并辅助实验数据处理。

通过使用Matlab,学生可以更好地理解和分析热传递过程,提高问题解决能力。

matlab 热传导方程

matlab 热传导方程

matlab 热传导方程
热传导方程是描述物质内部传热过程的数学模型,是研究热传导现象的基础。

在工程领域,热传导方程的求解可以应用于热力学、材料科学、建筑设计、能源开发等方面。

Matlab是一种常用的数学软件,它可以用来求解热传导方程。

在Matlab中,可以通过输入初始条件、边界条件和方程本身来求解热传导方程。

Matlab提供了多种求解方法,包括显式、隐式、Crank-Nicolson等方法。

对于一些复杂的热传导问题,Matlab还可以通过有限元方法来求解。

有限元方法将问题离散化,将连续的问题转化为离散的问题,通过数值计算方法求解得到近似的解。

总之,Matlab是一个非常强大的数学软件,可以应用于各种热传导问题的求解,为工程领域的研究和实践提供了很大的帮助。

- 1 -。

传热学实验—-墙角matlab导热问题

传热学实验—-墙角matlab导热问题

二维导热物体温度场的数值模拟姓名小明学号 111111班级能动学院能动一、问题描述有一墙角模型,尺寸如图1所示,导热系数0.53W/(m·K),墙角内外壁为第一类边界条件。

求解该模型的温度分布及导热量。

图1q=0二、计算原理根据热平衡法列出节点方程,各方向导入单元体的热量之和为零。

内节点和绝热边界点(图1点划线上的点)的方程形式不同。

图2 Array图2所示的内节点和绝热边界节点方程如下:内节点:)()()()(1,,1,,1,1,,1,=⎥⎦⎤⎢⎣⎡-+-+-+-••=+++-+-+x y t t x y t t y x t t y x t t j i j i j i j i j i j i j i j i W E S N ∆∆∆∆∆∆∆∆ΦΦΦΦλ绝热边界点:)(02)(2)(1,,1,1,,1,=⎥⎦⎤⎢⎣⎡-++-+-••=+++--+x y t t y xt t y x t t j i j i j i j i j i j i W E S N ∆∆∆∆∆∆ΦΦΦΦλ三、计算过程用Matlab7.1语言编写计算程序,初取网格步长m y x 1.0=∆=∆运行结果:图1:各个点的温度数值图2:分层设色等温线分布图3:等温线分布(每两条线间隔为三度)四、小结本次数值模拟是运用matlab程序用于数值计算。

小组成员共同讨论并复习了热传导问题的数学描述和热平衡法;从模拟过程中练习了不同节点迭代方程的建立;并简单学习了matlab语言的使用。

这次大作业对于我们以后的学习和可能的研究来说是一个很好的锻炼机会。

matlab移动高斯热源照射物体温度场计算

matlab移动高斯热源照射物体温度场计算

MATLAB是一种流行的数学建模和仿真软件,用于处理各种工程和科学问题。

在热传导领域,MATLAB可以用来计算热源照射物体的温度场。

本文将介绍如何使用MATLAB进行移动高斯热源照射物体的温度场计算。

一、题目背景热传导是工程和科学中的重要问题之一。

当一个物体表面受到热源的照射时,其温度场会发生变化。

通过计算热源照射后物体表面的温度分布,可以更好地理解和预测热传导过程,为工程设计和科学研究提供重要参考。

移动高斯热源照射物体是一种常见的热传导问题,利用MATLAB进行计算可以快速且准确地得到温度场的分布。

二、研究方法1. 建立热传导模型。

需要建立热传导方程和边界条件,以描述热源照射物体的温度场变化。

假设热源为高斯分布,物体表面具有一定的导热性和散热条件,可以建立相应的数学模型。

2. 离散化计算域。

将物体表面离散成网格,利用有限差分或有限元等方法对热传导方程进行离散化处理,以便在计算机上进行数值计算。

3. 计算高斯热源照射。

利用MATLAB编程,实现高斯热源的移动和照射过程。

通过数值方法,对物体表面的温度场进行时间步进计算,得到各个时刻的温度分布。

4. 可视化结果。

将计算得到的温度场数据以图形的形式进行可视化展示,以便更直观地观察高斯热源照射物体的温度场变化。

三、计算示例下面通过一个简单的计算示例,演示如何使用MATLAB进行移动高斯热源照射物体的温度场计算。

1. 建立热传导模型。

假设热源照射的物体为圆形,热源移动轨迹为直线运动,物体表面边界条件为第一类边界条件(即给定表面温度)。

热传导方程可用二维热传导方程表示:$\frac{\partial u}{\partial t} = \alpha \nabla^2 u$其中,u为温度,t为时间,α为热扩散系数。

2. 离散化计算域。

将圆形物体表面离散成网格,采用有限差分方法对二维热传导方程进行离散化处理。

3. 计算高斯热源照射。

编写MATLAB程序,实现高斯热源的移动轨迹和照射过程。

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

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

传热学作业数值计算1 程序内容: 数值计算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; =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= ; while(max1>1e-6) max1=0; k=k+1; for i=1:1:a for j=1:1:b m=m1(i,j); n= (i,j); switch m case 0 tn(i,j)=tw1; case 1 tn(i,j)=0.25*(tn(i,j+1)+tn(i,j-1)+tn(i+1,j)+tn(i-1,j)); case 2 tn(i,j)=tw1+tw2*sin(pi*(j-1)/(b-1)); end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k max1 t2=ones(a,b); %求解析温度场求解析温度场for i=a:-1:1 for j=1:1:b y=deltaX*(a-i); x=deltaX*(j-1); t2(i,j)=tw1+tw2*sin(pi*x/L1)*(sinh(pi*y/L1))/(sinh(pi*L2/L1)); end end t2 迭代次数k =706 数值解温度场 数值解每次迭代的最大误差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; =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= ; while(max1>1e-6) max1=0; k=k+1; for i=1:1:a for j=1:1:b m=m1(i,j); n= (i,j); switch m case 0 tn(i,j)=tw1; case 1 tn(i,j)=0.25*(tn(i,j+1)+tn(i,j-1)+tn(i+1,j)+tn(i-1,j)); case 2 tn(i,j)=tw2; end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k max1 tx=ones(a,b); for i=1:1:a for j=1:1:b y=(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; end end tx 迭代次数k = 695 数值解温度场 数值解每次迭代的最大误差max1 =9.8243e-07 解析温度场 tx = 解析温度场行的解析解和数值解的点取第11行的解析解和数值解的点行的解析解的直线,散点为其数值解的点 曲线为第11行的解析解的直线,散点为其数值解的点解析解 第11行的误差=[数值解(11行) –解析解(11行)]/解析解图像: 数值温度场 图像图像 :解析温度场tx图像:数值解与解析解的误差数值解与解析解的误差程序内容: 数值计算matlab程序内容:>> t0=90; =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; =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= ; max1=1.0; k=0; while ( max1>1e-6) k=k+1; max1=0; for i=1:1:a for j=1:1:b m=m1(i,j); n=tn(i,j); switch m case 0 tn(i,j)=t0; case 1 tn(i,j)=(2*tn(i-1,j)+tn(i,j-1)+tn(i,j+1)-4* )/(4+2*Bi)+ ; case 2 tn(i,j)=(2*tn(i,j-1)+tn(i-1,j)+tn(i+1,j)-4* )/(4+2*Bi)+ ; case 3 tn(i,j)=0.25*(tn(i,j-1)+tn(i,j+1)+tn(i-1,j)+tn(i+1,j)); case 4 tn(i,j)=(tn(i,j-1)+tn(i+1,j)-2* )/(2*Bi+2)+ ; case 5 tn(i,j)=(tn(i,j-1)+tn(i-1,j)-2* )/(2*Bi+2)+ ; case 6 tn(i,j)=(2*tn(i+1,j)+tn(i,j-1)+tn(i,j+1)-4* )/(4+2*Bi)+ ; case 7 tn(i,j)=t0; case 8 tn(i,j)=t0; end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k ta=ones(a,b); Bi1=h*d/lambda; sbi=sqrt(Bi1); for i=1:1:a for j=1:1:b if i>(a+1)/2 y=-(i-(a+1)/2)*deltaX; else y=((a+1)/2-i)*deltaX; end x=deltaX*(j-1); ta(i,j)=(cosh(sbi*(L-x)/d)+sbi*sinh(sbi*(L-x)/d))*(t0- )/(cosh(sbi*L/d)+sbi*sinh(sbi*L/d))+ ; end end ta 迭代次数k =1461 数值解温度场 解析温度场 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; =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= ; max1=1.0; k=0; while(max1>1e-6) max1=0; k=k+1; for i=1:1:a for j=1:1:b m=m1(i,j); n=tn(i,j); switch m case 0 tn(i,j)=tw; case 1 tn(i,j)=0.25*(tn(i-1,j)+tn(i+1,j)+tn(i,j-1)+tn(i,j+1)+qv0*(deltaX^2)/lambda); end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k; tx=ones(a,b); for i=1:1:a for j=1:1:b if i>(a+1)/2 y=-(i-(a+1)/2)*deltaX; else y=((a+1)/2-i)*deltaX; end if j>(b+1)/2 x=(j-(b+1)/2)*deltaX; else x=-((b+1)/2-j)*deltaX; end m=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; end end tx 数值温度场 解析温度场tx 取第13行的解析解和数值解的点行的解析解和数值解的点行的解析解的直线,散点为其数值解的点曲线为第13行的解析解的直线,散点为其数值解的点解析解 第13行的误差=[数值解(13行) –解析解(13行)]/解析解数值温度场图像数值温度场图像 解析温度场图像解析温度场图像数值解与解析解的误差数值解与解析解的误差。

二维热传导方程 matlab

二维热传导方程 matlab

二维热传导方程是描述二维热传导过程的数学模型,它在工程、物理、地球科学等领域都有重要应用。

在实际工程问题中,我们经常需要求解二维热传导方程,以预测物体表面的温度分布、热量传递速率等参数。

Matlab是一个强大的数学软件,通过Matlab我们可以很方便地求解二维热传导方程,并得到预期的结果。

一、二维热传导方程的基本形式二维热传导方程可以用偏微分方程的形式表示为:∂u/∂t = k(∂²u/∂x² + ∂²u/∂y²)其中,u(x, y, t)是温度分布随时间和空间的变化,k是热传导系数。

二、Matlab中求解二维热传导方程的方法在Matlab中,我们可以采用有限差分法(finite difference method)求解二维热传导方程。

有限差分法将偏微分方程离散化,转化为代数方程组,然后通过迭代求解得到数值解。

具体步骤如下:1. 离散化空间和时间变量,将连续的空间区域和时间区间分割成若干个小区间。

2. 利用二阶中心差分格式对二维热传导方程进行离散化,得到代数方程组。

3. 利用Matlab中的矩阵运算和迭代方法,求解代数方程组,得到数值解。

三、Matlab代码示例下面是一个简单的Matlab代码示例,用于求解二维热传导方程:```matlab定义参数和初始条件Lx = 1; Ly = 1; 区域大小Nx = 100; Ny = 100; 离散化网格数T = 1; 总时间Nt = 100; 时间步数k = 1; 热传导系数dx = Lx/Nx; dy = Ly/Ny;dt = T/Nt;x = 0:dx:Lx; y = 0:dy:Ly;[X, Y] = meshgrid(x, y);u = sin(pi*X).*sin(pi*Y); 初始温度分布迭代求解for n = 1:Ntun = u;for i = 2:Nx-1for j = 2:Ny-1u(i, j) = un(i, j) + k*dt/dx^2*(un(i+1, j)-2*un(i, j)+un(i-1, j)) + k*dt/dy^2*(un(i, j+1)-2*un(i, j)+un(i, j-1));endendend可视化结果figure;surf(X, Y, u);xlabel('x'); ylabel('y'); zlabel('Temperature');```以上代码首先定义了区域大小、离散化网格数、总时间、热传导系数等参数,然后利用有限差分法进行迭代求解,最后利用Matlab绘制了温度分布的三维图像。

二维 热传导方程 matlab

二维 热传导方程 matlab

二维热传导方程在实际工程问题中有着广泛的应用。

热传导方程描述了热量在空间和时间中的传播规律,而二维热传导方程则特指热量在二维平面上的传播情况。

在工程领域,我们经常需要利用数值方法来求解二维热传导方程,而MATLAB作为一款强大的工程计算软件,提供了丰富的工具和函数来进行热传导方程的数值求解。

让我们回顾一下二维热传导方程的基本形式,它通常可以表示为:\[ \frac{\partial u}{\partial t} = \alpha \left( \frac{\partial^2u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \]其中,\( u(x, y, t) \) 是温度场的分布,\( \alpha \) 是热扩散系数。

这个偏微分方程描述了温度场随时间和空间的演化规律,是研究热传导问题的基本方程之一。

在实际工程中,我们经常需要利用数值方法求解二维热传导方程,特别是在复杂几何结构和边界条件下。

MATLAB提供了丰富的数值求解工具,最常用的是有限差分法。

有限差分法将空间离散化为网格,用差分格式逼近偏微分方程,然后利用迭代方法求解离散化的代数方程组,得到温度场的数值解。

为了使用MATLAB进行二维热传导方程的数值求解,我们需要进行以下基本步骤:1. 网格划分:将求解区域进行网格划分,确定空间离散化的步长。

2. 边界条件处理:根据实际问题确定边界条件,例如固定温度、热通量等。

3. 数值格式选取:选择合适的差分格式逼近偏微分方程,通常采用显式或隐式格式。

4. 迭代求解:利用MATLAB提供的迭代方法,求解离散化后的代数方程组,得到温度场的数值解。

在实际工程中,二维热传导方程的数值求解往往涉及到复杂的边界条件、材料参数、热源分布等。

我们需要对MATLAB的数值求解工具有深入的理解和熟练的应用,才能准确地模拟和分析实际问题。

在我看来,二维热传导方程的数值求解是工程领域中非常重要的课题之一。

导热的反问题matlab

导热的反问题matlab

导热的反问题matlab
在MATLAB中,导热问题通常涉及热传导方程的数值求解。

热传导方程描述了物体内部温度分布随时间的变化,通常采用偏微分方程来描述。

解决导热问题的一种常见方法是使用有限差分法。

在MATLAB中,可以通过编写代码来离散化热传导方程,并使用迭代方法求解离散化后的方程。

另一种常见的方法是使用MATLAB的偏微分方程工具箱(Partial Differential Equation Toolbox)。

该工具箱提供了一系列函数和工具,可以帮助用户建立和求解偏微分方程,包括热传导方程。

用户可以通过定义边界条件、初始条件和热传导方程的参数来建立模型,并使用工具箱中的函数进行数值求解。

此外,MATLAB还提供了用于可视化和分析结果的丰富工具,例如绘制温度分布图、计算热通量等。

通过这些工具,用户可以全面分析导热问题的结果,并对模型进行验证和优化。

总之,在MATLAB中,可以通过编写代码、使用偏微分方程工具箱以及可视化分析工具来解决导热问题,从而全面深入地研究热传导现象。

热传导问题的MATLAB数值计算

热传导问题的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在导热问题中的应用

MATLAB在导热问题中的应用导热问题简介导热是指物质内部不同温度区域之间的热量传递现象。

在不同的热力学系统中,由于温度差异,导致热量从高温区域流向低温区域,以减少温度差异,直到两个区域相等为止,这个过程叫做导热。

在工业生产和科学研究中,导热问题是一个非常重要的问题,例如,建筑物的两面温度差、内部电子器件的散热等等都涉及到导热问题。

对于一些研究者而言,如何利用数学模型和计算机软件来解决导热问题,就成为了一个非常重要的课题。

MATLAB在导热问题中的应用MATLAB是一个非常强大的工具箱,因其拥有强大的计算功能,可以用于解决一些复杂的导热问题,例如:热传导方程热传导方程是描述物质中热量传递的基本方程,可以用MATLAB进行求解。

假设离散化的计算域中存在一系列温度节点,我们可以用以下公式表示热传导方程。

$$ \\dfrac{\\partial T}{\\partial t} = \ abla \\cdot (k \ abla T) $$其中,T为温度场变量,t为时间变量,k为热导率,abla表示热传导方程的梯度算子。

我们可以用MATLAB中的数值计算工具箱进行矩阵运算、微分运算等维度相关的计算,以求解这个方程。

边值问题在一些实际的导热问题中,会涉及到一些带边界的热传导问题,例如,房屋内的热传导问题,需要考虑外界空气温度对房屋内温度的影响。

这时,我们可以使用MATLAB中的偏微分方程工具箱,以求解带边值条件的问题。

辐射换热问题在一些高温应用场合,例如火车内部电力设备的散热问题,会涉及到辐射换热问题。

与传导换热不同,辐射换热是指物体表面和空间中其他物体表面之间的热量传递现象。

在这种情况下,我们可以使用MATLAB中的图像处理工具箱,通过计算辐射通量的分布来解决辐射换热问题。

结论综上所述,MATLAB可以用于解决一些复杂的导热问题,并且可以通过不同的工具箱进行平面模型、三维模型、带边值条件和辐射换热等不同类型的求解。

matlab偏微分方程

matlab偏微分方程

matlab偏微分方程Matlab可以用于求解偏微分方程(PDE)。

以下是一些示例:1. 热传导方程热传导方程描述了温度随时间和空间的变化,由以下方程给出:$\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partialx^2}$在Matlab中,可以使用“pdepe”函数来求解这个问题。

具体来说,需要指定初始条件和边界条件,并设置物理参数。

2. 波动方程波动方程描述了波的传播,由以下方程给出:$\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partialx^2}$在Matlab中,可以使用“pdepe”函数来求解这个问题。

需要指定初始条件和边界条件,并设置物理参数。

3. Navier-Stokes方程Navier-Stokes方程描述了流体的运动,由以下方程给出:$\frac{\partial u}{\partial t} + u \cdot \nabla u = -\frac{1}{\rho}\nabla p + \nu \nabla^2 u$在Matlab中,可以使用PDE工具箱进行求解。

需要指定初始条件、边界条件和物理参数。

4. Schrödinger方程Schrödinger方程描述了量子力学中的波函数演化,由以下方程给出:$i \hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m}\nabla^2 \psi + V(x) \psi$在Matlab中,可以使用PDE工具箱或ODE工具箱进行求解。

需要指定初始条件、边界条件和物理参数。

以上仅是部分示例,Matlab还可以用于求解其他类型的偏微分方程。

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

收稿日期: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 数值计算。

相关文档
最新文档