熱傳導(dǎo)方程的差分解法市公開課獲獎(jiǎng)?wù)n件

1、第九章 熱傳導(dǎo)方程差分解法第1頁第1頁9.1 熱傳導(dǎo)方程概述 考慮三維空間溫度改變情況, 設(shè) t 時(shí)刻點(diǎn)(x,y,z)處溫度為u(x,y,z,t), 則 t 時(shí)間內(nèi)通過橫截面積 S 傳導(dǎo)熱量為(沿 n 方向): 其中: K(x,y,z,t) 是介質(zhì)熱傳導(dǎo)系數(shù), 為溫度梯度法向量分量. 取空間中一個(gè)小區(qū)域 V, 其邊界面 S 為一封閉曲面. 則 t1 到 t2 時(shí)刻通過包面 S 傳入 V 熱量為: 第2頁第2頁由高斯公式:為哈密頓算子:設(shè)介質(zhì)比熱容為 c, 密度為 , 則 V 內(nèi)溫度改變消耗熱量: 設(shè) V 內(nèi)部熱源密度為 F(x,y,z,t), 則內(nèi)部熱源產(chǎn)生熱量為:第3頁第3頁依據(jù)能量守恒原則

2、: Q2 = Q1 + Q3即:亦即:若 F(x,y,z,t)0, c, , K, 為常數(shù),則:第4頁第4頁其中: 為拉普拉斯算子:因此熱傳導(dǎo)方程為:其中: Kc.第5頁第5頁 9.2 一維熱傳導(dǎo)方程差分解法一維熱傳導(dǎo)方程:初值問題初值條件: 初邊值混合問題初值條件:邊值條件:(關(guān)于邊界點(diǎn)x=0和x=l)第一類.第6頁第6頁第二類:第三類:其中g(shù)1(t), g2(t), 1(t), 2(t) 為給定函數(shù), 要求1(t), 2(t) , 且不同時(shí)為零. 第7頁第7頁設(shè)空間步長為 h, 時(shí)間步長為 .把空間和時(shí)間離散化:近似微分:故可定義: 對(duì)空間一階向前插商:第8頁第8頁對(duì)空間一階向后插商:對(duì)空

3、間二階中心差商:對(duì)時(shí)間一階向前插商:第9頁第9頁代入熱傳導(dǎo)方程:迭代公式:第10頁第10頁 t i-1 i i+1 x k+1 k第11頁第11頁第一類初邊值條件:已知:第12頁第12頁第二類初邊值條件:已知即:第13頁第13頁計(jì)算過程:第14頁第14頁第三類初邊值條件:已知:即:第15頁第15頁例1:差分方程:初邊值條件:第16頁第16頁function u = rcd(lamda,tao,h,H,T)x = 0:h:H;t = 0:tao:T;a = tao*lamda/h2;N = length(x);M = length(t);u(:,1) = (4*x.*(1-x);u(1,2:M)

4、 = 0;u(N,2:M) = 0;for k=1:M-1 for i=2:N-1 u(i,k+1)=a*u(i+1,k)+(1-2*a)*u(i,k)+a*u(i-1,k); endend第17頁第17頁h1=line(Color,1 0 0,Marker,.,MarkerSize,20,EraseMode,xor);for i=1:length(t) set(h1,Xdata,0:0.1:1,Ydata,u(:,i); pause(tao); end第18頁第18頁X,Y=meshgrid(x,0:0.01:0.2);Z = repmat(u(:,1),size(X,1),1);h2 =

5、surface(X,Y,Z);shading interp,axis equal;set(h2,EraseMode,xor);for i=1:length(t) CD = repmat(u(:,i),size(X,1),1); set(h2,Cdata,CD); pause(tao); end第19頁第19頁9.3 二維熱傳導(dǎo)方程差分解法內(nèi)部無熱源均勻介質(zhì)中二維熱傳導(dǎo)方程:初值條件: 邊值條件視詳細(xì)情況而定.設(shè)空間步長為 h, 時(shí)間步長為 . 設(shè)Nh=l, Mh=s, 把時(shí)間和空間離散化:第20頁第20頁即:微分近似:第21頁第21頁代入熱傳導(dǎo)方程得:k+1ki-1 i i+1j+1jj-1第

6、22頁第22頁例:初值條件: 即: 恒溫邊界:絕熱邊界:即:第23頁第23頁差分公式:第24頁第24頁function u = rcd2(lamda,tao,h,T,L,S)% 二維熱傳導(dǎo)方程t = 0:tao:T;x = 0:h:L;y = 0:h:S;a = tao*lamda/h2;if a0.25 error(lamda*tao/h20.25);end第25頁第25頁D = length(t);N = length(x);M = length(y);M1 = ceil(M/2)-3;M2 = ceil(M/2)+3;u = zeros(N,M,D);u(:,:,1) = 0;u(:,1

7、,:) = 0;u(:,M,:) = 0;u(1,M1:M2,2:D) = 1;第26頁第26頁for k=1:D-1 for i=2:N-1 for j=2:M-1 u(i,j,k+1)=(1-4*a)*u(i,j,k)+a*(u(i+1,j,k)+u(i-1,j,k)+u(i,j+1,k)+u(i,j-1,k); u(N,j,k+1) = u(N-1,j,k+1); if (jM2) u(1,j,k+1) = u(2,j,k+1); end end endend第27頁第27頁X,Y=meshgrid(x,y);Z = u(:,:,1);h2 = surface(X,Y,Z);shading interp,axis equal;set(h2,EraseMode,x