热传导方程PROJECT

热传导方程PROJECT

问题描述

热传导方程是一类典型的抛物线方程，它描述一个区域内温度随时间的变化规律，是一类重要的偏微分方程。求解热传导方程主要有数值解和解析解两种方法。我们以一维热传导问题为例，研究这两种方法的特点。

求解下面热传导方程：

2

2

2

2,0,0

(,0)sin ,0(0,)1,(,),0

u u u x l t t

x

u x x x l u t u l t t t π??-

=≤≤≥??=<<==≥

(1- 1)

给出上述热传导方程的解析解，并分别用显示格式、隐式格式及C-N 格式求其数值解，比较数值解与解析解之间的差别，分析三种差分格式误差精度。

1热传导方程解析解

将式（1-1）转化成热传导问题Cauchy 标准形式，令2(,)(,)t u x t v x t e =，则有

2

2

2220

(0,),(1,),0(,0)sin ,01

t

t

v v t

x

v t e v t t e

t v x x x π--??-

=??==≥=<< (1- 2)

再将边界条件齐次化，令1(,)(,)(,)v x t v x t w x t =+，()22(,)11t

w x t t x e -??=-+??

，则有

()2

2

2112

111211(0,)0,(1,)0,0(,0)sin 1,01

t

v v t t x e

t

x

v t v t t v x x x x π-????-

=--++-??

??==≥=+-<< (1- 3)

此定解问题的解可分解为112(,)(,)(,)v x t w x t w x t =+，其中1w ，2w 分别满足

2

112

1110

(0,)0,(1,)0,0(,0)sin 1,01

w w t

x

w t w t t w x x x x π??-

=??==≥=+-<<

(1- 4)

和

()2

22222

222211(0,)0,(1,)0,0(,0)0,01

t

w w t t x e

t

x

w t w t t w x x -????-

=--++-??

??==≥=<< (1- 5)

求解第一个定解问题，利用分离变量法得

2

22

12

2

2(,)(1)sin ()sin t

n t

n w x t e

x e

n x n πππππ

π

∞

--==-

+

-

∑

(1- 6)

对第二个定解问题采用齐次化原理，得

22

2

2

22

2

2

22

2

2

201

22

()

10

2(2)(2)10

(2)2

22

(,)(,;)(1)14[()

]sin()sin 4((1)())sin 14[(1)

((21)1(2)

t

t

n n t n t

t

n t

n n n n n t

n

n t

w x t w x t d e

e

n x d n n n x e

e

d e

d n n x e

t e

n n τ

πτππτ

πτ

ππτττ

ττπτ

π

π

πττττπππ

π+∞

---=∞

---=--=

--=

--+-

=--++

=--+-?∑

?∑

??2

2

2

2

2

2

1

(2)2

(2)(2)2

2

2

2

2

2

3

)

1

2(1)

(1)

(1)],01,0

2

2

(2)

n n t n

n t

n

n t

e

t t e

e

x t n n n ππππππ∞

=-----+

-----≤≤≥---∑

(1- 7)

最终，原问题的解即为

212(,)((,)(,)(,)),01,0t

u x t w x t w x t w x t e x t =++≤≤≥

(1- 8)

解析解得出的温度场分布图见 Figure 1。

Figure 1 解析解温度场分布

2热传导方程数值解

选取空间长度1l =，时间长度0.01T =，空间步长/dx l N =，时间步长/dt T M =，

1500

M =，100N =。令i x idx =，m t m dt =，[0,]i N ∈，[0,]m M ∈。

记(,)(,)m i u u idx m dt u i m ==。

显示格式

对所求解的偏微分方程做向前差分， 差分格式如下：

1

2

11

22

1

11

2

222m m

i i

m

m

m

i i i m m

m m m

m

i

i

i i i i

u u u t

dt

u u u u x dx

u u u u u u dt

dx

++-++--?=

?-+?=?--+-

= (2- 1)

最终，可得到该问题向前差分的显示格式：

1

11(122),11,0

m m m m

i

i i i u dt r u ru ru i N m ++-=+-++≤≤-≥ (2- 2)

其中，初、边界条件为：

02

sin(),01,0(),0i m

m

N u i dx i N

u m M u m dt m M

π=≤≤=≤≤=≤≤ (2- 3)

记矩阵,(,)[]i m U i m u =，由显示差分格式不难看出，如果已知第m 个时间m t 上各位置的函数值{}0

N

m m i i U U ==，则由式（2-2）即可相互独立地直接算出1(1,,1)m i U i N +=- 。 温度求解结果：

Figure 2 显示格式温度场分布

隐式格式

隐式差分格式如下：

1

1112

1

1

22

2m m

i

i

m m m i i i u u u

t dt u

u

u

u x

dx

+++++--?=

?-+?=

? (2- 4)

则根据偏微分方程可得：

1

1111(12)(12),11m m m m

i

i i i r u ru ru dt u i N +++-++--=+≤≤-

(2- 5)

其中，初、边界格式保持不变。

由式（2-5），对于任意一个时刻m 均有如下矩阵形式：

11

110122133

12

211

1112000

(12)1200(12)01200(12)00012(12)0

12(12)m m m m m

m m

m m N N m m m N N N r

r u dt u ru r r r u dt u r r u dt u r r u dt u r

r u dt u ru +++++--++--+-???++???????-+-+???

???????-++=?

???????????+-+????-+++???????

?

???????

???????

(2- 6)

其中，2

dt r dx

=

为隐式格式的网格比。

求解的温度场分布如下

:

Figure 3 隐式格式温度场分布

C-N 格式

C-N 差分格式如下：

1

111

2

1111222

2212m m

i i

m m m m m m i i i i i i u u u

t dt

u u u u u u u

x dx dx +++++-+--?=

???-+-+?=+?????

(2- 7)

根据偏微分方程可得：

1

111

111122

1

1

11111

221221111(1)(12)2

2

2

2

m m

m m m m m m m i

i

i i i i i i i m m m m

m

m

i

i i i i i u u u u u u u u u dt dx dx r u ru ru dt r u ru ru +++++-+-++++-+-??--+-+-+=????

+-

-

=+-+

+

(2- 8)

根据式（2-8），可得到对于任意一个时刻m 均有如下矩阵形式：

1

1

0,m m

m N AU

BU

C ++=+

(2- 9)

其中，2

dt r dx

=

1

1

1211

3

1211m m m m m N m N u u u U

u u +++++-+-????????=??????

??????

1

2321m

m m m

m N m N u u u U u u --????????=??????

??????

1

1

0,11

200

012

m m N

m N ru C ru +++????????

????=??

????????

?? (2- 10)

110002

11100221010021000121000

12r r

r r r

r

r A r r r

r ?

?

+-???

???

-+-

??????

-+=?

???????+-??????

-+???

?

(2- 11)

112000211120022101200210001221000

122dt r r

r

dt r

r

r

dt r

B dt r

r r

dt

r ?

?

+-??????+-??????+-=?

???????+-??????+-???

?

(2- 12)

求解的温度场分布如下：

Figure 4 C-N格式温度场分布

3误差分析

三种差分格式均在网格比0.0667

r ，网格数一致的情况下进行的误差计算。从表中可以看出，此时隐式格式具有较高的精度，误差较少，C-N格式误差较大。

4 Matlab代码

显示格式代码

clc;

clear all;

N=100;

M=1500;

T=0.01; %时间长度

L=1; %空间长度

dx=L/N; %空间步长

dt=T/M; %时间步长

U=zeros(N+1,M+1);

for i=1:N+1 %初边界条件赋值

U(i,1)=sin((i-1)*pi*dx);

end

U(1,:)=1;

for i=1:M+1

U(N+1,i)=((i-1)*dt)^2;

end

r=dt/(dx^2);

for m=1:M

for i=2:N

U(i,m+1)=(1+2*dt-2*r)*U(i,m)+r*U(i+1,m)+r*U(i-1,m); %显示格式计算式end

end

x=0:dx:N*dx;

t=0:dt:M*dt;

mesh(t,x,U); %显示格式数值解

xlabel('时间长度t');

ylabel('空间长度L');

title('温度场分布T');

for m=1:M+1 %输入解析解

for i=1:N+1

v1(i,m)=(((m-1)^2*dt^2-1)*(i-1)*dx+1)*exp(-2*(m-1)*dt);

end

end

for m=1:M+1

for i=1:N+1

v2(i,m)=0;

for n=2:50

v2(i,m)=v2(i,m)+(-2/(n*pi)*exp(-n^2*pi^2*(m-1)*dt))*sin(n*pi*(i-1)*dx);

end

v2(i,m)=v2(i,m)+(1-2/pi)*exp(-pi^2*(m-1)*dt)*sin(pi*(i-1)*dx);

end

end

for m=1:M+1

for i=1:N+1

v3(i,m)=0;

for n=1:50

a=sin(n*pi*(i-1)*dx)/(n*pi)*exp(-n^2*pi^2*(m-1)*dt);

b=(-1)^n/((n^2*pi^2-2)^2)*((2*(m-1)*dt-1)*exp((n^2*pi^2-2)*(m-1)*dt)+1);

c=(exp((n^2*pi^2-2)*(m-1)*dt)-1)/(n^2*pi^2-2);

d=((-1)^n*(((m-1)*dt)^2-(m-1)*dt)/(n^2*pi^2-2))*exp((n^2*pi^2-2)*(m-1)*dt );

e=2*(-1)^n/((n^2*pi^2-2)^3)*(exp((n^2*pi^2-2)*(m-1)*dt)-1);

v3(i,m)=v3(i,m)+a*(b+c-d-e);

end

end

end

for m=1:M+1

for i=1:N+1

v(i,m)=(v1(i,m)+v2(i,m)+v3(i,m))*exp(2*(m-1)*dt);

end

end%解析解输入结束

figure

x=0:dx:N*dx;

t=0:dt:M*dt;

mesh(t,x,v); %解析解图像

xlabel('时间长度t');

ylabel('空间长度L');

title('解析解——温度场分布T');

error_max=max(max(v-U)) %最大误差

error_mean=abs(mean(mean(v-U))) %平均误差

ordre_error=dt+dx^2 %局部截断误差

隐式格式代码

clc;

clear all;

N=100;

M=1500;

T=0.01; %时间长度

L=1; %空间长度

dx=L/N; %空间步长

dt=T/M; %时间步长

U=zeros(N+1,M+1);

for i=1:N+1 %初边界条件赋值

U(i,1)=sin((i-1)*pi*dx);

end

U(1,:)=1;

for i=1:M+1

U(N+1,i)=((i-1)*dt)^2;

end

V=zeros(N-1,1);

a=-1/(dx^2);

b=1/dt+2/(dx^2);

c=1/dt+2;

A=b*eye(N-1)+diag(a*ones(1,N-2),1)+diag(a*ones(1,N-2),-1); %隐式格式矩阵

for m=1:1:M %计算隐式格式数值解

V(1)=(1/(dx^2))*U(1,m+1)+U(2,m)*c;

V(N-1)=(1/(dx^2))*U(N+1,m+1)+U(N,m)*c;

V(2:N-2)=c*U(3:N-1,m);

U(2:N,m+1)=inv(A)*V;

end

x=0:dx:N*dx;

t=0:dt:M*dt;

mesh(t,x,U); %隐式格式数值解图像

xlabel('时间长度t');

ylabel('空间长度L');

title('温度场分布T');

for m=1:M+1 %输入解析解

for i=1:N+1

v1(i,m)=(((m-1)^2*dt^2-1)*(i-1)*dx+1)*exp(-2*(m-1)*dt);

end

end

for m=1:M+1

for i=1:N+1

v2(i,m)=0;

for n=2:50

v2(i,m)=v2(i,m)+(-2/(n*pi)*exp(-n^2*pi^2*(m-1)*dt))*sin(n*pi*(i-1)*dx);

end

v2(i,m)=v2(i,m)+(1-2/pi)*exp(-pi^2*(m-1)*dt)*sin(pi*(i-1)*dx);

end

end

for m=1:M+1

for i=1:N+1

v3(i,m)=0;

for n=1:50

a=sin(n*pi*(i-1)*dx)/(n*pi)*exp(-n^2*pi^2*(m-1)*dt);

b=(-1)^n/((n^2*pi^2-2)^2)*((2*(m-1)*dt-1)*exp((n^2*pi^2-2)*(m-1)*dt)+1);

c=(exp((n^2*pi^2-2)*(m-1)*dt)-1)/(n^2*pi^2-2);

d=((-1)^n*(((m-1)*dt)^2-(m-1)*dt)/(n^2*pi^2-2))*exp((n^2*pi^2-2)*(m-1)*dt );

e=2*(-1)^n/((n^2*pi^2-2)^3)*(exp((n^2*pi^2-2)*(m-1)*dt)-1);

v3(i,m)=v3(i,m)+a*(b+c-d-e);

end

end

end

for m=1:M+1

for i=1:N+1

v(i,m)=(v1(i,m)+v2(i,m)+v3(i,m))*exp(2*(m-1)*dt);

end

end%解析解输入结束

error_max=max(max(v-U)) %计算最大误差

error_mean=abs(mean(mean(v-U))) %计算平均误差

ordre_error=dt+dx^2 %计算局部截断误差

C-N格式代码

clc;

clear all;

N=100;

M=1500;

T=0.01; %时间长度

L=1; %空间长度

dx=L/N; %空间步长

dt=T/M; %时间步长

U=zeros(N+1,M+1);

for i=1:N+1 %初边界条件赋值

U(i,1)=sin((i-1)*pi*dx);

end

U(1,:)=1;

for i=1:M+1

U(N+1,i)=((i-1)*dt)^2;

end

V=zeros(N-1,1);

r=dt/(dx^2);

b=1+r;

a=(-1/2)*r;

c=1+2*dt-r;

d=(1/2)*r;

A=b*eye(N-1)+diag(a*ones(1,N-2),1)+diag(a*ones(1,N-2),-1); %C-N格式矩阵A B=c*eye(N-1)+diag(d*ones(1,N-2),1)+diag(d*ones(1,N-2),-1); %C-N格式矩阵B for m=1:M %计算C-N格式数值解

V=B*U(2:N,m);

V(1)=V(1)+d*U(1,m+1);

V(N-1)=V(N-1)+d*U(N+1,m+1);

U(2:N,m+1)=inv(A)*V;

end

x=0:dx:N*dx;

t=0:dt:M*dt;

mesh(t,x,U); %C-N格式数值解图像

xlabel('时间长度t');

ylabel('空间长度L');

title('温度场分布T');

for m=1:M+1 %输入解析解

for i=1:N+1

v1(i,m)=(((m-1)^2*dt^2-1)*(i-1)*dx+1)*exp(-2*(m-1)*dt);

end

end

for m=1:M+1

for i=1:N+1

v2(i,m)=0;

for n=2:50

v2(i,m)=v2(i,m)+(-2/(n*pi)*exp(-n^2*pi^2*(m-1)*dt))*sin(n*pi*(i-1)*dx);

end

v2(i,m)=v2(i,m)+(1-2/pi)*exp(-pi^2*(m-1)*dt)*sin(pi*(i-1)*dx);

end

end

for m=1:M+1

for i=1:N+1

v3(i,m)=0;

for n=1:50

a=sin(n*pi*(i-1)*dx)/(n*pi)*exp(-n^2*pi^2*(m-1)*dt);

b=(-1)^n/((n^2*pi^2-2)^2)*((2*(m-1)*dt-1)*exp((n^2*pi^2-2)*(m-1)*dt)+1);

c=(exp((n^2*pi^2-2)*(m-1)*dt)-1)/(n^2*pi^2-2);

d=((-1)^n*(((m-1)*dt)^2-(m-1)*dt)/(n^2*pi^2-2))*exp((n^2*pi^2-2)*(m-1)*dt );

e=2*(-1)^n/((n^2*pi^2-2)^3)*(exp((n^2*pi^2-2)*(m-1)*dt)-1);

v3(i,m)=v3(i,m)+a*(b+c-d-e);

end

end

end

for m=1:M+1

for i=1:N+1

v(i,m)=(v1(i,m)+v2(i,m)+v3(i,m))*exp(2*(m-1)*dt);

end

end%解析解输入结束

error_max=max(max(v-U)) %最大误差

error_mean=abs(mean(mean(v-U))) %平均误差

ordre_error=dt^2+dx^2 %局部截断误差

热传导方程的初值问题

§2热传导方程的初值问题 一维热传导方程的初值问题(或Cauchy 问题) ?? ???+∞<<∞-=>+∞<<∞-=??-??x x x u t x t x f x u a t u ),()0,(0 ,),,(2 2 2? （） 偏导数的多种记号xx x t u x u u x u u t u =??=??=??22,,. 问题也可记为 ?? ?+∞ <<∞-=>+∞<<∞-=-x x x u t x t x f u a u xx t ),()0,(0 ,,),(2?. Fourier 变换 我们将用Fourier 变换法求解热传导方程的柯西问题.为此我们将着重介绍Fourier 变换的基本知识.Fourier 变换在许多学科中是重要使用工具. 可积函数,设)(x f f =是定义在),(+∞-∞上的函数, 且对任意A B <，()f x 在[,]A B 上 可积，若积分 ? +∞ ∞ -dx x f )(收敛,则称)(x f 在),(+∞-∞上绝对可积。 将),(+∞-∞上绝对可积函数形成的集合记为),(1 +∞-∞L 或),(+∞-∞L ， 即{ } ∞<=+∞-∞=+∞-∞? +∞ ∞ -dx x f f L L )(| ),(),(1 ，称为可积函数空间. 连续函数空间: ),(+∞-∞上全体连续函数构成的集合,记为),(+∞-∞C , {}上连续在),(|),(+∞-∞=+∞-∞f f C ， {}上连续在),(,|),(1+∞-∞'=+∞-∞f f f C 。 定义 若),(+∞-∞∈L f ,那么积分 ),(?)(21 λπ λf dx e x f x i =? +∞ ∞ -- 有意义,称为Fourier 变换, )(? λf 称为)(x f 的Fourier 变式(或Fourier 变换的象). ? +∞ ∞ --= =dx e x f f Ff x i λπ λλ)(21)(?)( 定理 (Fourier 积分定理)若),(),(1 +∞-∞?+∞-∞∈C L f ,那么我们有

一维热传导方程

一维热传导方程 一． 问题介绍 考虑一维热传导方程： （1） ,0),(22T t x f x u a t u ≤<+??=?? 其中a 是正常数，)(x f 是给定的连续函数。按照定解条件的不同给法，可将方程（1）的定解问题分为两类： 第一类、初值问题（也称Cauthy 问题）：求具有所需次数偏微商的函数),(t x u ，满足方程（1）（∞<<∞-x ）和初始条件： （2） ),()0,(x x u ?= ∞<<∞-x 第二类、初边值问题（也称混合问题）：求具有所需次数偏微商的函数),(t x u ，满足方程（1）（l x <<0）和初始条件： （3） ),()0,(x x u ?= l x <<0 及边值条件 (4) .0),(),0(==t l u t u T t ≤≤0 假定)(x ?在相应区域光滑，并且在l x ,0=满足相容条件，使上述问题有唯一充分光滑的解。 二． 区域剖分 考虑边值问题（1），（4）的差分逼近。去空间步长N l h /=和时间步长M T /=τ，其中N,M 都是正整数。用两族平行直线： 将矩形域}0;0{T t l x G ≤≤≤≤=分割成矩形网格，网格节点为),(k j t x 。以h G 表示网格内点集合，即位于开矩形G 的网点集合；h G 表示所有位于闭矩形G 的网点集合；h Γ=h G --h G 是网格界点集合。 三． 离散格式 第k+1层值通过第k 层值明显表示出来，无需求解线性代数方程组，这样的格式称为显格式。 第k+1层值不能通过第k 层值明显表示出来，而由线性代数方程组确定，这样的格式称为隐格式。 1. 向前差分格式 （5） ,221 11j k j k j k j k j k j f h u u u a u u ++-=--++τ

导热方程求解matlab

使用差分方法求解下面的热传导方程 2 (,)4(,) 0100.2t xx T x t T x t x t =<<<< 初值条件：2(,0)44T x x x =- 边值条件：(0,)0(1,)0 T t T t == 使用差分公式 1,,1,2 2 2 (,)2(,)(,) 2(,)()i j i j i j i j i j i j xx i j T x h t T x t T x h t T T T T x t O h h h -+--++-+= +≈ ,1,(,)(,) (,)()i j i j i j i j t i j T x t k T x t T T T x t O k k k ++--= +≈ 上面两式带入原热传导方程 ,1,1,,1,2 2i j i j i j i j i j T T T T T k h +-+--+= 令2 24k r h =，化简上式的 ,1,1,1,(12)()i j i j i j i j T r T r T T +-+=-++ i x j t j

编程MA TLAB 程序，运行结果如下 1 x t T function mypdesolution c=1; xspan=[0 1]; tspan=[0 0.2]; ngrid=[100 10]; f=@(x)4*x-4*x.^2; g1=@(t)0; g2=@(t)0; [T,x,t]=rechuandao(c,f,g1,g2,xspan,tspan,ngrid); [x,t]=meshgrid(x,t); mesh(x,t,T); xlabel('x') ylabel('t') zlabel('T') function [U,x,t]=rechuandao(c,f,g1,g2,xspan,tspan,ngrid) % 热传导方程：

热传导方程

前言 本文只是针对小白而写，可以使新手对热传导理论由很浅到不浅的认识，如想更深学习热传导知识，请转其它文档。 一、概念与常量 1、温度场： 指某一时刻下，物体内各点的温度分布状态。 在直角坐标系中：； 在柱坐标系中：； 在球坐标系中：。 补充：根据温度场表达式，可分析出导热过程是几维、稳态或非稳态的现象，温度场是几维的、稳态的或非稳态的。 2、等温面与等温线： 三维物体内同一时刻所有温度相同的点的集合称为等温面； 一个平面与三维物体等温面相交所得的的曲线线条即为平面温度场中的等温线。 3、温度梯度： 在具有连续温度场的物体内，过任意一点P温度变化率最大的方向位于等温线的法线方向上。称过点P的最大温度变化率为温度梯度(temperature gradient)。用grad t表示。 定义为： 补充：温度梯度表明了温度在空间上的最大变化率及其方向，是向量，其正向与热流方向恰好相反。对于连续可导的温度场同样存在连续的温度梯度场。

在直角坐标系中： 3、导热系数 定义式：单位 导热系数在数值上等于单位温度降度(即1)下，在垂直于热流密度的单位面积上所传导的热流量。导热系数是表征物质导热能力强弱的一个物性参数。 补充：由物质的种类、性质、温度、压力、密度以及湿度影响。 二、热量传递的三种基本方式 热量传递共有三种基本方式：热传导；热对流；热辐射 三、导热微分方程式（统一形式：） 直角坐标系： 圆柱坐标系： 球坐标系： 其中，称为热扩散系数，单位，为物质密度，为物体比热容，为物体导热系数，为热源的发热率密度，为物体与外界的对流交换系数。 补充： 1处研究的对象为各向同性的、连续的、有内热源、物性参数已知的导热物体。 2稳态温度场，即则有：，此式称为泊松方程。 3无内热源的稳态温度场，则有：，此式称为拉普拉斯方程。 四、单值条件 导热问题的单值条件通常包括以下四项： 1几何条件：表示导热物体的几何形状与大小（一维、二维或三维）

热传导方程及其定解问题的导出

第一章 热传导方程 本章介绍最典型的抛物型方程—热传导方程,在研究热传导,扩散等物理现象时都会遇 到这类方程. §1 热传导方程及其定解问题的导出 1.1热传导方程的导出 物理模型 在三维空间中,考虑一均匀,各向同性的物体Ω,假定它内部有热源,并且与周围介质有热交换,需要来研究物体内部温度的分布和变化. 以函数),,,(t z y x u 表示物体Ω在位置),,(z y x 及时刻t 的温度.物体内部由于各部分温度不同,产生热量的传递,它们遵循能量守恒定律. 能量守恒定律 物体内部的热量的增加等于通过物体的边界流入的热量与由物体内部的热源所生成的热量的总和 . 在物体Ω内任意截取一块D .现在时段],[21t t 上对D 使用能量守恒定律. 设),,,(t z y x u u =是温度(度),c 是比热(焦耳∕度·千克),ρ是密度(千克/米3), q 是热流密度(焦耳/秒·米2),0f 是热源强度(焦耳/千克·秒). 注意到在dt 时段内通过D 的边界D ?上小块dS 进入区域D 的热量为dSdt n q ?-(n 是 D ?的外法向),从而由能量守恒律,我们有 ,)||(21 21 120??????????+?-=-?==t t D t t D D t t t t dxdydz f dt ds n q dt dxdydz u u c ρρ （1.1） 大家知道,热量流动的原因是因为在物体内部存在温差.依据传热学中的傅立叶实验定律,在一定条件下,热流向量与温度梯度成正比 ,u k q ?-= （梯度? ?? ? ????????==?z u y u x u gradu u ,,） （1.2） 这里负号表明热量是由高温向低温流动,k 是物体的导热系数.

热传导方程的求解

应用物理软件训练 前言 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到1 eps是MATLAB默认的最小浮点数精度 [X,Y]=pol2cart(R,TH);效果和上一句相同 waterfall(RR,TT,wn)瀑布图 二、基本原理 1、一维热传导问题 （1）无限长细杆的热传导定解问题 利用傅里叶变换求得问题的解是： 取得初始温度分布如下 这是在区间0到1之间的高度为1的一个矩形脉冲，于是得 （2）有限长细杆的热传导定解问题

热传导方程傅里解

热传导方程傅里解

————————————————————————————————作者：————————————————————————————————日期：

热传导在三维的等方向均匀介质里的传播可用以下方程表达： 其中： ?u =u(t, x, y, z) 表温度，它是时间变量t 与空间变量(x,y,z) 的函数。 ?/是空间中一点的温度对时间的变化率。 ?, 与温度对三个空间座标轴的二次导数。 ?k决定于材料的热传导率、密度与热容。 热方程是傅里叶冷却律的一个推论（详见条目热传导）。 如果考虑的介质不是整个空间，则为了得到方程的唯一解，必须指定u 的边界条件。如果介质是整个空间，为了得到唯一性，必须假定解的增长速度有个指数型的上界，此假定吻合实验结果。 热方程的解具有将初始温度平滑化的特质，这代表热从高温处向低温处传播。一般而言，许多不同的初始状态会趋向同一个稳态（热平衡）。因此我们很难从现存的热分布反解初始状态，即使对极短的时间间隔也一样。 热方程也是抛物线偏微分方程最简单的例子。 利用拉普拉斯算子，热方程可推广为下述形式

其中的是对空间变量的拉普拉斯算子。 热方程支配热传导及其它扩散过程，诸如粒子扩散或神经细胞的动作电位。热方程也可以作为某些金融现象的模型，诸如布莱克-斯科尔斯模型与 Ornstein-Uhlenbeck 过程。热方程及其非线性的推广型式也被应用于影像分析。量子力学中的薛定谔方程虽然有类似热方程的数学式（但时间参数为纯虚数），本质却不是扩散问题，解的定性行为也完全不同。 就技术上来说，热方程违背狭义相对论，因为它的解表达了一个扰动可以在瞬间传播至空间各处。扰动在前方光锥外的影响通常可忽略不计，但是若要为热传导推出一个合理的速度，则须转而考虑一个双曲线型偏微分方程。 以傅里叶级数解热方程[编辑] 以下解法首先由约瑟夫·傅里叶在他于1822年出版的著作Théorie analytique de la chaleur（中译：解析热学）给出。先考虑只有一个空间变量的热方程，这可以当作棍子的热传导之模型。方程如下： 其中u = u(t, x) 是t和x的双变量函数。 ?x是空间变量，所以x∈[0,L]，其中L表示棍子长度。

3热传导方程的初边值问题

例4 周期初始温度分布 求解热传导方程t xx u u =,(,0)x t -∞<<+∞>给定初始温度分布 (,0)1cos 2,()u x x x =+-∞<<+∞。 解 4(,)1cos2t u x t e x -=+. 初始高斯温度分布 例 5求解定解问题22 22 0,(,0) (,0),()kx u u a x t t x u x e x -???-=-∞<<+∞>?????=-∞<<+∞? ， 其中常数0k >. 解 22()4(,)()x s a t u x t s e ds ?-- +∞ -∞ = ? 22 2()4x s ks a t e e ds -- +∞ --∞ = ? 222 2(41)24ka t s xs x a t e ds +-+- +∞ -∞ = ? 222 22224(41)()41414x ka t ka t s x ka t ka t a t e ds +- +++- +∞ -∞ = ? 22 2 222(41)()41 441 k ka t x x s ka t a t ka t e e ds +---+∞ ++-∞ = ? 2241 k x ka t e - += 2241 k x ka t - += . §3初边值问题 设长度为l ,侧表面绝热的均匀细杆,初始温度与细杆两端的温度已知,则杆上的温度分布 ),(t x u 满足以下初边值问题 ?? ? ??≤<==≤≤=<<<<=-T t t g t l u t g t u l x x x u T t l x t x f u a u xx t 0),(),(),(),0(,0), ()0,(0,0),,(212? 对于这样的问题,可以用分离变量法来求解. 将边值齐次化

第四章导热问题的数值解法

第四章导热问题的数值解法 1 、重点内容：①掌握导热问题数值解法的基本思路； ②利用热平衡法和泰勒级数展开法建立节点的离散方程。 2 、掌握内容：数值解法的实质。 3 、了解内容：了解非稳态导热问题的两种差分格式及其稳定性。 §4—1导热问题数值求解的基本思想及内节点方程的建立由前述 3 可知，求解导热问题实际上就是对导热微分方程在定解条件下的积分求解，从而获得分析解。但是，对于工程中几何形状及定解条件比较复杂的导热问题，从数学上目前无法得出其分析解。随着计算机技术的迅速发展，对物理问题进行离散求解的数值方法发展得十分迅速，并得到广泛应用，并形成为传热学的一个分支——计算传热学（数值传热学），这些数值解法主要有以下几种： （１）有限差分法（ 2 ）有限元方法（ 3 ）边界元方法 数值解法能解决的问题原则上是一切导热问题，特别是分析解方法无法解决的问题。如：几何形状、边界条件复杂、物性不均、多维导热问题。 一．分析解法与数值解法的异同点： ?相同点：根本目的是相同的，即确定① t=f(x ， y ， z) ；② 。 ?不同点：数值解法求解的是区域或时间空间坐标系中离散点的温度分布代替连续的温度场；分析解法求解的是连续的温度场的分布特征，而不是分散点的数值。 数值求解的基本思路及稳态导热内节点离散方程的建立 二．解法的基本概念 ?实质 对物理问题进行数值解法的基本思路可以概括为：把原来在时间、空间坐标系中连续的物理量的场，如导热物体的温度场等，用有限个离散点上的值的集合来代替，通过求解按一定方法建立起来的关于这些值的代数方程，来获得离散点上被求物理量的值。该方法称为数值解法。 这些离散点上被求物理量值的集合称为该物理量的数值解。 2 、基本思路：数值解法的求解过程可用框图 4-1 表示。 由此可见： 1 ）物理模型简化成数学模型是基础； 2 ）建立节点离散方程是关键； 3 ）一般情况微分方程中，某一变量在某一坐标方向所需边界条件的个数等于该变量在该坐标方向最高阶导数的阶数。 ?数值求解的步骤 如图 4-2 （ a ），二维矩形域内无内热源、稳态、常物性的导热问题采用数值解法的步骤如下：（1）建立控制方程及定解条件 针对图示的导热问题，它的控制方程（即导热微分方程）为：（ a ） 边界条件： x=0 时， x=H 时， 当 y=0 时，

热传导方程抛物型偏微分方程和基本知识

1. 热传导的基本概念 1.1温度场 一物体或系统内部，只要各点存在温度差，热就可以从高温点向低温点传导， 即产生热流。因此物体或系统内的温度分布情况决定着由热传导方式引起的传热速率（导热速率）。 温度场：在任一瞬间，物体或系统内各点的温度分布总和。 因此，温度场内任一点的温度为该点位置和时间的函数。 〖说明〗 若温度场内各点的温度随时间变化，此温度场为非稳态温度场，对应于非稳 态的导热状态。 若温度场内各点的温度不随时间变化，此温度场为稳态温度场，对应于稳态 的导热状态。 若物体内的温度仅沿一个坐标方向发生变化，且不随时间变化，此温度场为 一维稳态温度场。 1.2 等温面 在同一时刻，具有相同温度的各点组成的面称为等温面。因为在空间同一点不可能同时有两个不同的温度，所以温度不同的等温面不会相交。 1.3 温度梯度 从任一点起沿等温面移动，温度无变化，故无热量传递；而沿和等温面相交 的任一方向移动，温度发生变化，即有热量传递。温度随距离的变化程度沿法向最大。 温度梯度：相邻两等温面间温差△t与其距离△n之比的极限。 〖说明〗 温度梯度为向量，其正方向为温度增加的方向，与传热方向相反。 稳定的一维温度场，温度梯度可表示为：grad t = dt/dx

2. 热传导的基本定律——傅立叶定律 物体或系统内导热速率的产生，是由于存在温度梯度的结果，且热流方向和 温度降低的方向一致，即与负的温度梯度方向一致，后者称为温度降度。 傅立叶定律是用以确定在物体各点存在温度差时，因热传导而产生的导热速率大小的定律。 定义：通过等温面导热速率，与其等温面的面积及温度梯度成正比： q = dQ/ds = -λ·dT/dX 式中：q 是热通量（热流密度），W/m2 dQ是导热速率，W dS是等温表面的面积，m2 λ是比例系数，称为导热系数，W/m·℃ dT / dX 为垂直与等温面方向的温度梯度 “－”表示热流方向与温度梯度方向相反 3. 导热系数 将傅立叶定律整理，得导热系数定义式： λ= q/(dT/dX) 物理意义：导热系数在数值上等于单位温度梯度下的热通量。因此，导热系 数表征物体导热能力的大小，是物质的物性常数之一。其大小取决于物质的组成结构、状态、温度和压强等。 导热系数大小由实验测定，其数值随状态变化很大。 3.1 固体的导热系数 金属：35～420W/(m·℃)，非金属：0.2～3.0W/ (m·℃) 〖说明〗

第八章 热传导和扩散问题的傅里叶解

第八章 热传导方程的傅里叶解 第一节 热传导方程和扩散方程的建立 8.1.1 热传导方程的建立 推导热传导方程和前面弦振动所用的数学方法完全相用，不同之处在于具体的物理规律不同。这里用到的是热学方面的两个基本规律，即能量守恒和热传导的傅里叶实验定律。 热传导的傅里叶实验定律：设有一块连续的介质，选定一定的坐标系，并用(,,,)u x y z t 表示介质内空间坐标为的一点在t 时刻的温度。若沿x 方向有一定的温度差，在x 方向也就一定有热量的传递。从宏观上看，单位时间内通过垂直x 方向的单位面积的热量q 与温度的沿x 方向的空间变化率成正比，即 x u q k x ?=-? （8-1.1） q 称为热流密度，k 称为导热系数。公式中的负号表示热流的方向和温度变化的方向正好相 反，即热量由高温流向低温。 研究三维各向同性介质中的热传导，在介质中三个方向上存在温度差，则有 x u q k x ?=-?，y u q k y ?=-?，z u q k z ?=-? 或 q k u =-?r 即热流密度矢量q r 与温度梯度u ?成正比。 下面以一维均匀细杆为例，根据傅里叶实验定律和能量守恒定律推导介质中的热传导方程。 第一步，定变量。研究介质x 位置处在t 时刻的温度(,)u x t 。 第二步，取局部。在介质内部隔离出从x 到x x +?一段微元长度，在t 到t t +?时间内温度的变化(,)(,)u u x t t u x t ?=+?-。 第三步，立假设。假设均匀介质的横截面积为A ，质量密度为ρ，比热为c ，热传导系数为k 。

第四步，找规律。隔离出来的微元长度在t 到t t +?时间内吸收的热量为： Q c m u c A x u ρ=????=???? （8-1.2） 在t 到t t +?时间内，同过x 位置处的横截面的热量为： 1x x x Q q A t k u A t =???=-?? （8-1.3） 在t 到t t +?时间内，同过x x +?位置处的横截面的热量为： 2x x x x x Q q A t k u A t +?+?=???=-?? （8-1.4） 如果在微元段内有其他的热源，假设在单位时间单位体积内产生的热量为(,)F x t ，则该热源在微元内产生的热量为： (,)3Q F x t t A x =???? （8-1.5） 第五步，列方程。根据能量守恒定律，净流入的热量应该等于介质在此时间内温度升高所需要的热量。 123Q Q Q Q =-+ 即 (,)x x x x x c A x u k u A t k u A t F x t t A x ρ+?????=-??+??+???? (,)x x x x x u u u c k F x t t x ρ+?-??? =+?? 得到： (,)t xx k F x t u u c c ρρ = + 令 a = (,)(,)F x t f x t c ρ= 则得到热传导方程为 (,)2t xx u a u f x t =+ （8-1.6） 当介质内部无其他热源时，热传导方程是齐次的，为 2t xx u a u = （8-1.7） 8.1.2 扩散方程的建立