传热学上机C程序源答案之一维稳态导热的数值计算
传热学上机C程序源答案之一维稳态导热的数值计算
一维稳态导热的数值计算1.1物理问题一个等截面直肋，处于温度t∞=80的流体中。
肋表面与流体之间的对流换热系数为h=45W/(m2∙℃),肋基处温度tw=300℃，肋端绝热。
肋片由铝合金制成，其导热系数为λ=110W/(m ∙℃)，肋片厚度为δ=0.01m，高度为H=0.1m 。
试计算肋内的温度分布及肋的总换热量。
1.2数学描述及其解析解引入无量纲过余温度θ=t -t∞tw -t∞，则无量纲温度描述的肋片导热微分方程及其边界条件：2220d m dxθθ-=x=0,θ=θw =1 x=H,0xθ∂=∂ 其中 Ahpm =λ上述数学模型的解析解为：[()]()()w ch m x H t t t t ch mH ∞∞--=-⋅()()w hpt t th mH m∞∅=-1.3数值离散1.3.1区域离散计算区域总节点数取N 。
1.3.2微分方程的离散对任一借点i 有：2220i d m dx θθ⎛⎫-= ⎪⎝⎭用θ在节点i 的二阶差分代替θ在节点i 的二阶导数，得：211220i i i i m x θθθθ+--+-=整理成迭代形式：()112212i i i m x θθθ+-=++ (i=2,3……,N-1)1.3.3边界条件离散补充方程为：11w θθ==右边界为第二类边界条件，边界节点N 的向后差分得：10N N xθθ--=，将此式整理为迭代形式，得：N 1N θθ-=1.3.4最终离散格式11w θθ==()112212i i i m xθθθ+-=++ (i=2,3……,N-1) N 1N θθ-=1.3.5代数方程组的求解及其程序假定一个温度场的初始发布，给出各节点的温度初值：01θ，02θ，….，0N θ。
将这些初值代入离散格式方程组进行迭代计算，直至收敛。
假设第K 步迭代完成，则K+1次迭代计算式为：K 11w θθ+=()11112212i i K K K i m xθθθ+-++=++ (i=2,3……,N-1) 111N K K N θθ-++=#include<stdio.h>#include<math.h>#define N 11main(){int i;float cha;/*cha含义下面用到时会提到*/float t[N],a[N],b[N];float h,t1,t0,r,D,H,x,m,A,p; /*r代表λ,x代表Δx，D代表δ*/printf("\t\t\t一维稳态导热问题\t\t");printf("\n\t\t\t\t\t\t----何鹏举\n");printf("\n题目：补充材料练习题一\n");printf("已知：h=45，t1=80, t0=200, r=110, D=0.01, H=0.1 (ISO)\n");/*下面根据题目赋值*/h=45.0; t1=80.0; t0=300.0; r=110.0; D=0.01; H=0.1;x=H/N; A=3.1415926*D*D/4; p=3.1415926*D; m=sqrt((h*p)/(r*A));/*x代表步长，p代表周长，A代表面积*/printf("\n请首先假定一个温度场的初始分布，即给出各节点的温度初值：\n");for(i=0;i<N;i++){scanf("%f",&t[i]);a[i]=(t[i]-t1)/(t0-t1);b[i]=a[i];/*这里b[i]用记录一下a[i]，后面迭代条件及二阶采用温度初场要用到*/ }/*采用一阶精度的向后差分法数值离散*/cha=1;while(cha>0.0001){a[0]=1;for(i=1;i<N;i++)a[i]=(a[i+1]+a[i-1])/(2+m*m*x*x);a[N-1]=a[N-2];cha=0;for(i=0;i<N;i++)cha=cha+a[i]-b[i];cha=cha/N;/*cha代表每次迭代后与上次迭代各点温度差值的平均值*/}for(i=0;i<N;i++)t[i]=a[i]*(t0-t1)+t1;printf("\n\n经数值离散(一阶精度的向后差分法)计算得肋片的温度分布为:\n");for(i=0;i<N;i++)printf("%4.2f\t",t[i]);printf("\n\n");getchar();/*采用二阶精度的元体平衡法数值离散(温度初值还用设定的初场，便于比较)*/ for(i=0;i<N;i++)a[i]=b[i];cha=1;while(cha>0.0001){a[0]=1;for(i=1;i<N;i++)a[i]=(a[i+1]+a[i-1])/(2+m*m*x*x);a[N-1]=a[N-2]/(1+0.5*m*m*x*x);cha=0;for(i=0;i<N;i++)cha=cha+a[i]-b[i];cha=cha/N;}for(i=0;i<N;i++)t[i]=a[i]*(t0-t1)+t1;printf("\n\n经数值离散(二阶精度的元体平衡法)计算得肋片的温度分布为:\n"); for(i=0;i<N;i++)printf("%4.2f\t",t[i]);printf("\n\n");getchar();}-----精心整理，希望对您有所帮助!。
(完整word版)一维非稳态导热的数值计算
int i,j,l;
float cha;
float a,x,y,Fo,Bi;
float t[N][K],b[N][K];
/*打印出题目*/
printf("\t\t\t一维非稳态导热问题\t\t");
printf("\n\t\t\t\t\t\t----何鹏举\n");
printf("\n题目：补充材料练习题三\n");
/*时刻为零时，赋予初场温度*/
for(i=0;i<N;i++)
t[i][0]=1000;
/*循环开始，每次计算一个时刻*/
for(j=0;j<K-1;j++)
{
for(i=0;i<N;i++)
b[i][j]=t[i][j];
/*下面对每一个时刻进行迭代求解对应的温度分布，公式按传热学课本P178页公式*/
y=1;/*y代表Δτ*/
x=0.05/(N-1);
a=34.89/(7800*712);
Fo=(a*y)/(x*x);
Bi=233*x/34.89;
printf("\n显示格式条件:");
printf("\n1、Fo=%3.1f<0.5\t",Fo);
printf("\t2、1-2Fo*Bi-2Fo=%4.2f>0\n\n",1-2*Fo*Bi-2*Fo);
{
printf("\n");
l=0;
}
}ห้องสมุดไป่ตู้
getchar();/*为了是生成的exe文件结果算的后不会立即退出，方便观看*/
传热学上机C程序源答案之一维非稳态导热的数值计算
二维稳态导热的数值计算2.1物理问题一矩形区域，其边长L=W=1，假设区域内无内热源，导热系数为常数，三个边温度为T1=0，一个边温度为T2=1，求该矩形区域内的温度分布。
2.2 数学描述 对上述问题的微分方程及其边界条件为：2222T T 0x y∂∂+=∂∂ x=0，T=T 1=0x=1，T=T 1=0y=0，T=T 1=0y=1，T=T 2=1 该问题的解析解：112121(1)sin n n n sh y T T n L x n T T n L sh W L ππππ∞=⎛⎫⋅ ⎪---⎛⎫⎝⎭=⋅ ⎪-⎛⎫⎝⎭⋅ ⎪⎝⎭∑ 2.3数值离散2.3.1区域离散区域离散x 方向总节点数为N ，y 方向总节点数为M ，区域内任一节点用I,j 表示。
2.3.2方程的离散 对于图中所有的内部节点方程可写为：2222,,0i j i jt t x y ⎛⎫⎛⎫∂∂+= ⎪ ⎪∂∂⎝⎭⎝⎭ 用I,j 节点的二阶中心差分代替上式中的二阶导数，得：+1,,-1,,+1,,-1222+2+0i j i j i ji j i j i j T T T T T T x y --+=上式整理成迭代形式：()()22,1,-1,,1,-12222+2()2()i j i j i j i j i j y x T T T T T x y x y ++=++++ (i=2,3……,N -1)，(j=2,3……,M -1)补充四个边界上的第一类边界条件得：1,1j T T = (j=1,2,3……,M),1N j T T = (j=1,2,3……,M),1i j T T = (i=1,2,3……,N),2i M T T = (i=1,2,3……,N)#include<stdio.h>#include<math.h>#define N 10#define K 11main(){int i,j,l;float cha;float a,x,y,Fo,Bi;float t[N][K],b[N][K];/*打印出题目*/printf("\t\t\t一维非稳态导热问题\t\t");printf("\n\t\t\t\t\t\t----何鹏举\n");printf("\n题目：补充材料练习题三\n");y=1;/*y代表Δτ*/x=0.05/(N-1);a=34.89/(7800*712);Fo=(a*y)/(x*x);Bi=233*x/34.89;printf("\n显示格式条件:");printf("\n1、Fo=%3.1f<0.5\t",Fo);printf("\t2、1-2Fo*Bi-2Fo=%4.2f>0\n\n",1-2*Fo*Bi-2*Fo);/*时刻为零时，赋予初场温度*/for(i=0;i<N;i++)t[i][0]=1000;/*循环开始，每次计算一个时刻*/for(j=0;j<K-1;j++){for(i=0;i<N;i++)b[i][j]=t[i][j];/*下面对每一个时刻进行迭代求解对应的温度分布，公式按传热学课本P178页公式*/ cha=1;while(cha>0.001){for(i=0;i<N-1;i++){if(i==0)t[i][j+1]=Fo*(t[i+1][j]+t[i+1][j])+(1-2*Fo)*t[i][j];/*当计算t[0]时，要用到t[-1],其中t[-1]=t[2]的(对称分布)*/elset[i][j+1]=Fo*(t[i+1][j]+t[i-1][j])+(1-2*Fo)*t[i][j];t[N-1][j+1]=t[N-2][j]*(1-2*Fo*Bi-2*Fo)+2*Fo*t[N-1][j]+2*Fo*Bi*20;/*边界点温度用热平衡法推导出公式*/}cha=0;for(i=0;i<N;i++)cha=cha+abs(t[i][j]-b[i][j]);cha=cha/N;}}/*输出温度分布，其中l控制输出值的排列;这个结果是横轴为x，纵轴为τ的直角坐标下从左上角开始依次的*/printf("\n经数值离散计算的温度分布为：\n");l=0;for(j=K-1;j>=0;j--)for(i=0;i<N;i++){if(t[i][j]>999.99)printf("%6.1f ",t[i][j]);elseprintf("%6.2f ",t[i][j]);l=l+1;if(l==N){printf("\n");l=0;}}getchar();/*为了是生成的exe文件结果算的后不会立即退出，方便观看*/}。
一维稳态热传导方程的数值解法及其
由上两式有：
xe xe xe
e
P
E
（5）
此式即为界面上的当量导热系数调和平均公式，
它可以看成是串联过程中热阻叠加原则的反映。
当网格划分为均匀网格时
e
2P E P E
（4）
6.3 两种方法的比较
1）当λ E 0时，由4式λe 0，说明在一个绝热层
界面）
4）把物性阶跃面设置成一个节点的位置比作为控 制容积分界面，使计算结果会更加精确。（由于此种情
况阶跃面两侧温度梯度不同，如按3处理，相当于用平均值来代替，采用此种方法处理 时，物性阶跃面两侧温度梯度单独计算，提高了计算精度。）
一维稳态导热方程的离散形式可表示成：
aP T PaE T Ea W T W b
的表面上qe=0，合乎实际；但 3 式 λe 0;
2）如 P
分时， e
E
P
2
，按算术平均法，当网格为均匀划
E
P
2
则P，E间的导热阻力为
2
x
e
P
，说明P，E间的导
热热阻由导热系数大的决定 ，这是不对的。
若按调和平均法计算，由5式则导热热阻为
xe xexexe
w
TPxTwW
整理得：
TExT PTPxT WxSCSPTP0
e
w
e
w
TPxee xwwSpxTExeeTWxwwScx
简化成 aP T PaE T Ea W T W b （2）
§6 一维稳态热传导方程的数 值解法及其应用
6.1 一维稳态导热的通用控制方程
一维稳态导热方程离散化、边界条件及源项的处 理及非线性代数方程的求解方法等对对流问题数 值解也适用。 一维稳态导热微分方程的通用形式为
一维稳态导热问题数值计算
一维稳态导热问题数值计算刘强引言❖目前为止，一般稍微复杂的导热问题几乎都依靠数值法求解。
❖导热问题的数值法有三种：有限差分法，有限元法和边界元法。
本教材介绍目前在铸造领域温度场计算中普遍采用的直接差分法，也叫单元热平衡法。
❖基本思想：不用导热微分方程，而是直接通过能量守恒定律，根据相邻单元间的能量交换关系导出差热方程。
❖分析i 单元的热量平衡关系，从t n 到t n+1时间内，由i-1单元流入i 单元的热量为：=1Q x i T i T k n n ∆---)1()(x ∆⋅（1）由i 单元流入i+1单元的热量为：=2Q 由内能计算公式：t x i T i T k n n ∆⋅∆-+-)()1(Tm C Q p ∆=而在该时间内，得出单元的内能增量为：[])()(1i T i T C x Q n n p -∆=+ρ蓄（2）（3）根据能量守恒定律则能得出蓄Q Q Q =-21t x i T i T k n n ∆⋅∆---)1()([])()()()1(1i T i T C x t x i T i T k n n p n n -∆=∆⋅∆-+++ρ或是其中[])1()()2()1(1++-+-i T i T M i T M n n n tx M ∆⋅∆=α/上式即为显式差分格式（4）=+)(1i T n初始条件：边界条件：给定初始温度T （i ），i=1，2,3，…，N由初始和边界条件可计算区域内部各节点随时间t 变化的温度值：代表时间步常数给定边界温度n n N T T nn ,,2,1,0),(),1(⋅⋅⋅=),3,2,1;1,,3,2(),(⋅⋅⋅=-⋅⋅⋅=n N i i T n步骤如下由初始条件和边界条件知图中第0排的温度，知，其中由初始条件提供)1(~)2(T 00-N T 由边界条件提供，与)()1(00N T T 第一排的温度值)1,,3,2)(1(1-⋅⋅⋅=N i T 可由（4）式得到；再利用边界条件，得到),()1(11N T T 与即能得到第一排上的全部节点的温度再由（4）式和边界条件依次算得inT n⋅⋅⋅==⋅⋅⋅i),,),2,1;(,3,2(n显示与隐式差分格式)(1i T n +)(1i T n +)1()()1(+-i T i T i T n n n 、、在4式中，n+1排上的任一节点i 的温度只依赖在n 排上i 节点及相邻节点i-1、i+1的温度值换言之，就是可由明显地来表示出来⇒显示差分格式若用)1()()1(111+-+++i T i T i T n n n 、、时刻的温度去计算1+n t tx i T i Tk Q n n ∆⋅∆---=++)1()(111t x i T i T k Q n n ∆⋅∆-+-=++)()1(112,21Q Q 、则能得到（5）（6）结合（3）式便得到另一种差分格式)()1(1)()21()1(1111i T i T Mi T M i T M n n n n =+-++--+++（7）此式只是表示的时间水平不同，实际上⇒与（4）式形势完全一致式（7）即完全隐式差分格式谢谢。
一维稳态导热数值解法matlab
一维稳态导热数值解法matlab 在导热传输的研究中，解析方法常常难以适用于复杂的边界条件和非均匀材料性质的情况。
因此，数值解法在求解热传导方程的问题上发挥了重要作用。
本文将介绍一维稳态导热数值解法，以及如何使用MATLAB来实现。
稳态导热数值解法通常基于有限差分法(Finite Difference Method, FDM)，它将连续的一维热传导方程离散为一组代数方程。
首先，我们需要将热传导方程转化为差分格式，然后利用MATLAB编写程序来求解。
下面，将具体介绍该方法的步骤。
步骤一：离散化根据一维热传导方程，可以将其离散为一组差分方程。
假设被研究的材料长度为L，将其等分为N个离散节点。
令x为节点位置，T(x)表示节点处的温度。
则可以得到以下差分方程：d²T/dx² ≈ (T(x+Δx) - 2T(x) + T(x-Δx)) / Δx²其中，Δx = L/N是节点之间的间距。
将热传导方程在每个节点处应用上述差分格式后，我们便得到了一组代表节点温度的代数方程。
步骤二：建立矩阵方程将差分方程中各节点的温度代入，我们可以将其表示为一个线性方程组。
这个方程组可以用矩阵的形式表示为Ax = b，其中A是系数矩阵，x是节点温度的向量，b是右侧项的向量。
步骤三：求解方程组使用MATLAB的线性方程求解器可以直接求解上述的线性方程组。
具体而言，通过利用MATLAB中的"\ "操作符，我们可以快速求解未知节点的温度向量x。
步骤四：结果分析与可视化在得到节点温度向量后，我们可以对结果进行可视化和分析。
例如，可以使用MATLAB的plot函数绘制温度随位置的分布曲线，以及温度随节点编号的变化曲线。
这样可以直观地观察到温度的变化情况。
总结：本文介绍了一维稳态导热数值解法以及使用MATLAB实现的步骤。
通过将热传导方程离散化为差分方程，然后建立矩阵方程并利用MATLAB的线性方程求解器求解，我们可以得到节点温度的数值解。
一维稳态导热数值计算
一维稳态导热数值计算引言在工程和科学领域中，热传导是一个重要的问题，它涉及到物体内部的热量传递过程。
一维稳态导热是指物体在一个方向上的热传导过程，且不随时间变化。
为了分析和解决一维稳态导热问题，我们可以使用数值计算方法，如有限差分法。
本文将介绍一维稳态导热数值计算的基本原理和步骤。
基本原理一维稳态导热问题可以描述为以下的热传导方程：$$\\frac{{d}}{{dx}}(k \\frac{{dT}}{{dx}}) = 0$$其中，k是物质的热导率，T是温度。
我们需要根据边界条件和初始条件求解该方程的解析解或数值解。
在数值求解中，我们通常将问题的区域离散化，将连续变量转化为离散变量。
我们可以将区域划分为多个小区间，每个小区间内的温度和导热系数近似为常数。
然后，我们可以使用有限差分法来近似求解。
数值计算步骤为了进行一维稳态导热问题的数值计算，我们需要按照以下步骤进行操作：步骤 1：确定区域和边界条件首先，我们需要确定问题的区域，并确定边界条件。
区域可以是一根导热杆或其他具有一维结构的物体。
边界条件可以是固定温度或热流量。
步骤 2：离散化区域将区域离散化是数值计算的基础。
我们可以将区域划分为多个小区间，每个小区间内的温度和导热系数近似为常数。
确定离散化的步长可以根据问题的要求进行选择。
步骤 3：建立差分方程根据离散化后的区域，我们可以建立差分方程，将热传导方程转化为一个线性方程组。
在一维稳态导热问题中，通常采用中心差分法或其他差分格式进行近似。
步骤 4：求解线性方程组求解差分方程就是求解线性方程组。
我们可以使用常见的数值计算工具或算法，如高斯消元法或迭代法，来求解线性方程组。
根据边界条件的不同，方程组的形式也会有所不同，需要根据具体情况进行选择。
步骤 5：计算结果最后，根据线性方程组的解，我们可以计算出每个小区间内的温度分布。
可以根据具体需求进行进一步计算和分析。
总结本文介绍了一维稳态导热数值计算的基本原理和步骤。
一维稳态热传导方程的数值解法及其
具体步骤如下：（1）先假设一个温度分布初值；
（2）计算相应函数b, a n b 及 a p
（3）求解线性离散方程组； （4）由新的温度再计算函数（改进系数）；
（5）返回2后，再重复计算T，直到 104 为止。
其中
Tn1 Tn
Tn
设初值为T*，迭代后新的温度分布为T，
例如在热传导问题中SP为正值，意味着TP增加，源项热源也增加，如果这时没有有效的散热机构，可能会反 过来导致温度的升高，如此反复下去，造成温度飞升的不稳定现象。
为了保证代数方程迭代求解的收敛。Δν为控制容积的体积， 线性代数方程迭代求解收敛的一个充分条件是对角占优，即
ap anbSPV
ap anb
,这里A是控制体积界面的面积，这里取1，于是ΔV= ΔX
从而有
d dT x e d dT x w xSCSP T P0
对扩散项T 随x 呈分段线性分布得：
dT dx
e
e
TExTeP理得：
TExT ePTPxT w WxSCSPTP0
e
w
TPxee xwwSpxTExeeTWxwwScx
1) S c =4 S p=-5
2) S c =4-5Tp* S p=0
3) S c =4+7Tp* S p=-12
2）中将S作为常数（以上一次迭代计算的T*计算S）处理，使源项相对于T永远有一 个滞后；1）中Tp是迭代计算当前值使S能更快跟上Tp的变化；3）比实际的S~ T 关系更陡的曲线，使迭代收敛速度减慢，相当于欠松弛。
一维稳态导热方程的离散形式可表示成：
aP T PaE T Ea W T W b
1
aE
传热学2.3 典型一维稳态导热问题的分析解
1笛卡尔坐标系中三维非稳态导热微分方程的一般表达式·)()()(Φ+∂∂∂∂+∂∂∂∂+∂∂∂∂=∂∂ztz y t y x t x t c λλλτρ边界条件——导热物体边界上温度或换热情况第一类边界条件()0w t f ττ>=时第二类边界条件20()()w tf nτλτ∂>−=∂时第三类边界条件()()w w f th t t nλ∂−=−∂定解条件初始条件——初始时间温度分布非稳态项扩散项源项物理问题→数学描写→微分方程①导热系数为常数c zt y t x t a tρτ·222222)(Φ+∂∂+∂∂+∂∂=∂∂②导热系数为常数 、无内热源 222222()t t t ta x y zτ∂∂∂∂=++∂∂∂∂③导热系数为常数 、稳态·2222220t t t x y z λ∂∂∂Φ+++=∂∂∂简化④导热系数为常数 、稳态 、无内热源 2222220t t tx y z ∂∂∂++=∂∂∂·)()()(Φ+∂∂∂∂+∂∂∂∂+∂∂∂∂=∂∂ztz y t y x t x t c λλλτρ⑴ 物理问题：大平壁，λ=const.⑵ 数学描写：微分方程边界条件·)()()(Φ+∂∂∂∂+∂∂∂∂+∂∂∂∂=∂∂z tz y t y x t xt c λλλτρ导热微分方程稳态、一维、无内热源、常物性t 1t 2q oδxtdx1. 单层平壁⑶ 解微分方程：⎪⎩⎪⎨⎧=−=⇒12121t c t t c δ112t x t t t +−=δ线性分布带入Fourier 定律δ12d d t t x t −=⇒)(12λδλδδλA ttt t q Δ=ΦΔ=−−=⇒—— 温度分布—— 通过平壁导热的计算公式共同规律可表示为 :2.热阻的含义过程中的转换量 = 过程中的动力 / 过程中的阻力)(λδA tΔ=Φ如：欧姆定律A R R A δδλλ==热阻分析法适用于一维、稳态、无内热源的情况RU I /=平板导热:转移过程的动力转移过程的阻力导热过程的转移量面积热阻热阻tq δλΔ=}多层平壁：由几层不同材料组成}例：房屋的墙壁 — 白灰内层、水泥沙浆层、红砖（青砖）主体层等组成}假设各层之间接触良好，可以近似地认为接合面上各处的温度相等边界条件：⎪⎩⎪⎨⎧====+=∑1110n n i i t t x t t x δ热 阻：nn n r r λδλδ==,,111L 3.多层平壁的导热t 1t 2t 3t 4t 1t 2t 3t 4三层平壁的稳态导热热阻的特点：串联热阻叠加原则：在一个串联的热量传递过程中，若通过各串联环节的热流量相同，则串联过程的总热阻等于各串联环节的分热阻之和。
一维传热问题数值计算
一维传热问题数值计算
一维传热问题是热传导理论中的经典问题，涉及热量在一个维
度上的传递和分布。
数值计算一维传热问题通常涉及使用数值方法
来模拟热量在材料中的传递和分布。
这个问题在工程、物理学和材
料科学等领域都有重要的应用。
首先，我们可以考虑使用有限差分法来数值计算一维传热问题。
有限差分法将材料空间离散化为若干个网格点，然后利用热传导方
程进行离散化，最终转化为一个差分方程。
通过迭代求解这个差分
方程，我们可以得到材料中温度随时间和空间的分布。
这种方法通
常需要考虑边界条件和初始条件，以及选择合适的时间步长和空间
步长。
另外，有限元法也是计算一维传热问题的常用数值方法。
有限
元法将材料分割为有限个小单元，然后利用单元间的热传导关系建
立整个系统的方程。
通过求解这些方程，可以得到材料中温度的分布。
有限元法通常适用于复杂几何形状的材料，并且可以很好地处
理不均匀材料性质的情况。
除了这些基本的数值方法，还可以考虑使用计算流体动力学
（CFD）方法来模拟一维传热问题。
CFD方法可以更全面地考虑流体在传热过程中的影响，适用于液体或气体在管道或其他结构中的传热问题。
在进行数值计算一维传热问题时，需要注意选择合适的数值方法和参数，以确保计算结果的准确性和稳定性。
同时，还需要考虑材料的热物性参数、边界条件、初始条件等因素，以保证数值模拟的真实性和可靠性。
总之，数值计算一维传热问题涉及多种数值方法和复杂的物理过程，需要综合考虑材料性质、边界条件和数值方法的选择，以获得准确而可靠的计算结果。
一维稳态导热数值解法matlab
一维稳态导热数值解法matlab 导热是物体内部热量传递的一种方式，对于一维稳态导热问题，我们可以使用数值解法来求解。
MATLAB是一种强大的数值计算软件，可以方便地实现一维稳态导热数值解法。
首先，我们需要了解一维稳态导热问题的基本原理。
一维稳态导热问题可以用一维热传导方程来描述，即：d²T/dx² = Q/k其中，T是温度，x是空间坐标，Q是热源的热量，k是热导率。
我们需要求解的是温度T在空间上的分布。
为了使用数值解法求解这个方程，我们需要将空间离散化。
假设我们将空间分成N个小区间，每个小区间的长度为Δx。
我们可以将温度T在每个小区间的位置上进行离散化，即T(i)表示第i个小区间的温度。
接下来，我们可以使用有限差分法来近似求解热传导方程。
有限差分法的基本思想是用差分代替微分，将微分方程转化为差分方程。
对于一维热传导方程，我们可以使用中心差分公式来近似求解：(T(i+1) - 2T(i) + T(i-1))/Δx² = Q(i)/k其中，Q(i)是第i个小区间的热源热量。
将上述差分方程整理后，可以得到：T(i+1) - 2T(i) + T(i-1) = (Q(i)/k) * Δx²这是一个线性方程组，我们可以使用MATLAB的矩阵运算功能来求解。
首先，我们需要构建系数矩阵A和常数向量b。
系数矩阵A是一个(N-1)×(N-1)的矩阵，其中A(i,i) = -2，A(i,i+1) = A(i,i-1) = 1。
常数向量b是一个(N-1)维的向量，其中b(i) = (Q(i)/k) * Δx²。
然后，我们可以使用MATLAB的线性方程组求解函数来求解这个方程组。
假设我们将求解得到的温度向量为T_solve，那么T_solve就是我们所求的稳态温度分布。
最后，我们可以使用MATLAB的绘图功能来可视化温度分布。
通过绘制温度随空间坐标的变化曲线，我们可以直观地观察到温度的分布情况。
计算传热学第6节-第2章 一维导热3
TP1 =0.334T =0.334*4.64+2.118=3.67 P2+2.118
TP2 =0.602*5.24+1.49=4.64 TP3 =0.718*5.45+1.328=5.24
P4
P5 Pa Pb
代入可得 aP1=aP10 aP2=aP20 aP3=aP30 aP4=aP40 aP5=aP50
= = = = =
aE1(TP1,TP2)TP2+aW1(TP1)Tpa +6 aE2(TP2,TP3)TP3+aW2(TP2,TP1)TP1+6 aE3(TP3,TP4)TP4+aW3(TP3,TP2)TP2+6 aE4(TP4,TP5)TP5+aW4(TP4,TP3)TP3+6 aE5(TP5) TPb+aW5(TP5,TP4)TP4+6
解方程可得 TP1=TP11 TP2=TP21 TP3=TP31 TP4=TP41 TP5=TP51 代入可得 aP1=aP11 aP2=aP21 aP3=aP31 aP4=aP41 aP5=aP51
Pa Pb
TPa =3 TPb =5 解方程可得 TP1=TP1n TP2=TP2n TP3=TP3n TP4=TP4n TP5=TP5n
一维导热讨论问题
导热系数为x的函数λ
Pa P1 P2
λ=10W/m·k
(x)
P3
λ=50W/m·k P4 P5 Pb T=5℃
T=3℃
Δx=0.6m Δx=0.6m Δx=0.6m Δx=0.6m Δx=0.6m
一维稳态热传导方程
一维稳态热传导方程
一维稳态热传导方程
一维稳态热传导方程是热传导问题的基本方程，其中，热传导系数K是用温度未变量的函数。
热传导方程的求解是热储存问题的基本算法，常用Finite element法和Finite Difference法求解。
一维稳态热传导方程
一维稳态热传导方程用来描述传统的单调方程热导率, 它表示
时间和空间的变化：
/d2T/dx2 + q_e = 0
其中：/d2T/dx2 为温度对x的二阶导数，q_e 为热源或热损失。
一般应用条件
1. 物质性质不变：
流体恒定的密度和热容，固体恒定的力学强度，导热系数和热容。
2. 境界条件：边界条件需要定义，自由表面是定义化学反应率的温度，或热损失。
热传导方程的求解
常用Finite element法和Finite Difference法求解。
Finite element法：采用有限元素方法，把区域分成小的等边形或三角形元素，然后求解每一个元素对应的温度场。
Finite Difference法：将方程化为一组一次方程组，然后由矩阵的求解得到方程组的解。
结论
一维稳态热传导方程是热传导问题的基本方程，其中，热传导系数K是用温度未变量的函数，然后常用Finite element法和Finite Difference法求解。
一维稳态热传导方程
一维稳态热传导方程
* 《二维共轭热传导方程》.doc
### 二、计算方法
1、一维稳态热传导方程求解
稳态热传导方程可以把热温度当成一种流体，速度可以理解为温度的变化率。
求解一维稳态热传导方程，可以根据牛顿热传导方程： $frac{partial T(x,t)}{partial t} = frac{kappa}{
ho C_p} frac{partial ^2 T(x,t)}{partial x^2}$ 其中$T$是温度，$kappa$是传热系数，$
ho$是物质密度，$C_p$是比容热容量。
根据牛顿热传导方程，可以用有限差分法(finite difference method)来求解一维稳态热传导方程。
这种方法基本思路是用一系列
离散点近似地描述出温度分布，用差分公式及多项式估算温度的变化，进而求解出待求温度分布。
2、二维共轭热传导方程求解
二维共轭热传导方程可以表示为：
$frac{partial^2 T(x,y)}{partial x^2} + frac{partial^2
T(x,y)}{partial y^2} = 0$
一般来讲，我们可以将二维共轭热传导方程拆分成一维的，用有限差分法就可以求解出温度分布。
也可以用拉普拉斯变换法求解，也可以用三角形网格法，或者用椭圆分形网格法。
- 1 -。
传热学一维稳态导热
传热学一维稳态导热传热学是物理学和工程学中一个重要的分支，研究热量在物质中的传递过程。
在传热学中，导热是其中一个重要的热传递方式。
导热是指热量通过传导传递，不涉及物质的移动。
在一维稳态导热的条件下，我们将详细介绍导热的基本原理和计算方法。
一维稳态导热的基本理论一维稳态导热是指热量沿一个方向传导，而且在传导过程中温度分布保持不变。
在一维稳态导热中，我们可以使用傅立叶热传导定律来描述热量的传导过程。
傅立叶热传导定律表明，单位时间通过导热展面的热流量与温度变化率成正比。
数学上可以表示为：$$ q = -k\\frac{{dT}}{{dx}} $$其中，q表示单位时间通过导热展面的热流量，k表示导热系数，dT表示温度的变化量，dx表示距离的微小变化量。
导热系数k是物质的属性，用于衡量物质传热的能力。
单位为W/(m·K)。
根据傅立叶热传导定律，可以得到温度随距离变化的微分方程。
在一维稳态导热中，由于温度分布保持不变，微分方程可以简化为：$$ q = -k\\frac{{dT}}{{dx}} = const $$这意味着在一维稳态导热中，热流量在传导过程中保持不变。
这是因为传热过程中能量守恒的原理。
一维稳态导热的计算方法在一维稳态导热的条件下，我们可以通过解微分方程来计算温度分布和热流量。
以下是一维稳态导热计算的基本步骤：1.确定热传导的边界条件：在一维稳态导热中，通常需要给定两个边界条件，例如温度或热流量。
这些边界条件用于确定问题的求解范围和约束条件。
2.确定物质的导热性质：导热系数k是物质传热能力的关键参数，需要根据材料的物性参数进行选择。
通常可以通过查表或实验来获取。
3.设定坐标系和建立微分方程：在一维稳态导热中，需要选择一个坐标系，并根据傅立叶热传导定律建立微分方程。
根据边界条件确定微分方程的边界条件。
4.求解微分方程：通过求解微分方程，可以得到温度随距离变化的数学表达式。
这将给出热流量和温度分布的解析解。
导热问题的数值解法一
T5 = Tr
2 = 0 = λ 1 W / m ⋅K) (
一维稳态常物性导热问题
有限容积积分(内点)
∆x x i-1 i-1/2 i i+1/2 i+1
∫
xi +1/ 2
xi −1/ 2
d dT λ dx dx
dT 0⇒λ dx = dx
对所求变量给出预估值，然后根据离 散方程对其不断改进，直至求得收敛解。
一维稳态常物性导热问题
Jacobi迭代法： 任意点上未知量的更新都用上一轮迭代中 所获得的值来计算
1 (I ) (I ) Ti aiTi + c T = + 1 i i −1 bi
( I +ቤተ መጻሕፍቲ ባይዱ)
(
)
Gauss-Seidel迭代法： 每一步计算取邻点的最新值
界面导热系数的计算方式
λ= i +1/ 2 λ= i −1/ 2
1 ( λi + λi +1 ) 2 1 ( λi + λi −1 ) 2
算术平均
哪种好？
∆x x i-1 i i-1/2 i+1/2 i+1
λi +1/ 2 =
2 1
λi λi −1/ 2 =
1
+
1
λi +1
1
2
调和平均
λi
+
λi −1
Ti +1 − Ti ⇒ = 1 Ti +1 − Ti 2 ⇒λ = i +1/ 2 1 1 1 1 + + 2λi +1 2λi λi +1 λi
一维稳态导热方程推导
一维稳态导热方程推导导热方程简介在理解一维稳态导热方程之前，我们先来了解导热方程的概念。
导热方程描述了热量在介质中的传播规律，通过该方程可以推导出温度在空间中的分布情况。
一维稳态导热方程的定义一维稳态导热方程描述了在一个维度上，介质中的热量分布不随时间变化的情况下的热传导行为。
该方程可以用数学形式表示为：$$\\frac{{d}}{{dx}}(k(x)\\frac{{dT}}{{dx}}) = 0$$其中，k(x)是介质的导热系数，T(x)是温度分布函数，x表示一维空间坐标。
推导过程为了推导一维稳态导热方程，我们需要考虑热传导过程中的热流量和热量守恒。
假设我们有一个长度为L的介质，其中的热传导是沿着x方向进行的。
首先，我们考虑介质内的一个微小段dx，该段温度为T(x)，导热系数为k(x)。
根据热量守恒定律，该微小段内的热量变化率应与进入和离开该段的热流量之和相等。
我们可以将该微小段的热量变化率表示为：$$\\frac{{d}}{{dt}}(Q(x)) = -\\frac{{d}}{{dx}}(F(x))$$其中，Q(x)表示该微小段内的热能，F(x)表示通过该微小段的热流量。
考虑微小段dx与其邻近的微小段之间的热传导，我们可以使用傅里叶定律来表示热流量F(x)：$$F(x) = -k(x) \\frac{{dT}}{{dx}} dx$$代入热量守恒定律的表达式中，我们得到：$$\\frac{{d}}{{dx}}(k(x)\\frac{{dT}}{{dx}}) = 0$$这就是一维稳态导热方程。
方程的物理意义一维稳态导热方程描述了热量在一个维度上的传播规律。
方程表明，在稳态下，介质中的热传导是由温度梯度驱动的。
温度梯度越大，热传导越强。
导热系数k(x)描述了介质内导热的特性，它可以反映介质的性质和结构。
通过求解一维稳态导热方程，我们可以确定介质内不同位置的温度分布。
这对于许多工程问题和科学研究都是非常重要的，例如热传导问题、传热设备设计等。
【石油工程】传热学-第3章-稳态导热的计算与分析
m
tw2 tw1
x
0
t
w1
b 2
t 2w1
这表明，当材料的导热系数随温度呈线性规律变化时
，平壁内的温度分布是二次曲线方程，该二次曲线的凹凸
性主要由温度系数b的正负决定。
利用傅里叶定律分析表明：
——b>0时，温度分布曲线的开口向下；
——b<0时曲线开口向上
17
3.1.3 第一类边界条件下变物性、无内热源的平壁
dx
11
3.1.2 第一类边界条件下的常物性、无内热源的平壁
常物性、无内热源平壁稳态导热的计算公式：
Φ A tw1 tw2
q tw1 tw 2
稳态法测定物质导热系数的基本依据
常物性、无内热源的条件下，平壁一维稳态导热的 热流量或热流密度为常数
由此可以采用另一种方法得到平壁内的温度分布
• 厚度为δ 、侧面积为A的单层平壁 • 没有内热源，导热系数λ为常数 • 两侧表面分别维持均匀稳定的温度 tw1、tw2，且tw1>tw2
d dt Φ 0
dx dx
9
3.1.2 第一类边界条件下的常物性、无内热源的平壁
• 导热问题的数学描述为
d2t dx 2
0
边界条件为：
t x0 tw1
4
3.1.1 平壁一维稳态导热的数学模型
（2） 物理模型 墙壁、玻璃、罐壁等物体具有相
似的几何特征 ——某一方向的尺寸远远小于其他两 个方向的尺寸
将高度和宽度远远大于厚度(8～10倍)的物体称为大
平壁，简称平壁。基本尺寸有平壁厚度δ和面积A。
5
3.1.1 平壁一维稳态导热的数学模型
平壁一维稳态导热简化的基础： ——平壁的几何特征 ——平壁的传热特点： （1）平壁两侧换热均匀（沿高度、 宽度方向），忽略边缘效应 （2）温度变化发生在平壁的厚度方 向上
一维稳态导热
t1 r1
t2 r2 t3 r3
t4
15
推广到n层壁的情况:
qt1n tn1
i
t总 rd总
i1 i
W m2
t1 tn1 t总
n
i
Rd总
A i1 i
W
16
例2－2由三层材料组成的加热炉炉墙。第一层为耐 火砖。第二层为硅藻土绝热层，第三层为红砖，各 层 的 厚 度 及 导 热 系 数 分 别 为 1 ＝ 240mm ， 1=1.04W/(m℃)， 2＝50mm, 2=0.15W/(m℃)， 3＝115mm, 3=0.63W/(m℃)。炉墙内侧耐火砖 的表面温度为1000℃。炉墙外侧红砖的表面温度 为60℃。试计算硅藻土层的平均温度及通过炉墙 的导热热流密度。
的热流密度为
q (t1 t2 ) 0 0 ..8 23 7 5 0 2 5 1.4 W 7/m 2
9
A dt dx
0 1
bt
A
dt dx
；
分离变量
t w 2 (1 bt ) dt dx ；
tw1
0 0A
( t w 2
tw1 )
b 2
t
w
2 2
t
w
2 1
0A
＝3t
ln3 0.1426
23 223 2
两种情况散热量之比为：
q1 0.1426 1.1或 9qL0.84 qL 0.11969 qL
结论：导热系数大的材料在外面，导热系数小 的材料放在里层对保温更有利。
31
三、通过球壁的导热
A r dt 4 r 2 dt
dr
dr
t w 2 dt r2 dr
5
1、单层平壁
数值传热学部分习题答案
习题4－2一维稳态导热问题的控制方程：022=+∂∂S xTλ 依据本题给定条件，对节点2节点3采用第三类边界条件具有二阶精度的差分格式，最后得到各节点的离散方程： 节点1： 1001=T节点2： 1505105321-=+-T T T 节点3：75432=+-T T 求解结果：852=T ，403=T对整个控制容积作能量平衡，有：02150)4020(15)(3=⨯--⨯=∆+-=∆+x S T T h x S q f f B即：计算区域总体守恒要求满足习题4－5在4－2习题中，如果25.03)(10f T T h -⨯=，则各节点离散方程如下：节点1： 1001=T节点2： 1505105321-=+-T T T节点3：25.03325.032)20(4015])20(21[-⨯+=-⨯++-T T T T对于节点3中的相关项作局部线性化处理，然后迭代计算； 求解结果：818.822=T ，635.353=T （迭代精度为10-4）迭代计算的Matlab 程序如下： x=30; x1=20;while abs(x1-x)>0.0001a=[1 0 0;5 -10 5;0 -1 1+2*(x-20)^(0.25)]; b=[100;-150; 15+40*(x-20)^(0.25)]; t=a^(-1)*b; x1=x; x=t(3,1);endtcal=t习题4-12的Matlab程序%代数方程形式A i T i=C i T i+1+B i T i-1+D imdim=10;%计算的节点数x=linspace(1,3,mdim);%生成A、C、B、T数据的基数；A=cos(x);%TDMA的主对角元素B=sin(x);%TDMA的下对角线元素C=cos(x)+exp(x); %TDMA的上对角线元素T=exp(x).*cos(x); %温度数据%由A、B、C构成TDMAcoematrix=eye(mdim,mdim);for n=1:mdimcoematrix(n,n)=A(1,n);if n>=2coematrix(n,n-1)=-1*B(1,n);endif n<mdimcoematrix(n,n+1)=-1*C(1,n);endend%计算D矢量D=(coematrix*T')';%由已知的A、B、C、D用TDMA方法求解T%消元P(1,1)=C(1,1)/A(1,1);Q(1,1)=D(1,1)/A(1,1);for n=2:mdimP(1,n)=C(1,n)/(A(1,n)-B(1,n)*P(1,n-1));Q(1,n)=(D(1,n)+B(1,n)*Q(1,n-1))/(A(1,n)-B(1,n)*P(1,n-1));end%回迭Tcal(1,mdim)=Q(1,mdim);for n=(mdim-1):-1:1Tcal(1,n)=P(1,n)*Tcal(1,n+1)+Q(1,n);endTcom=[T;Tcal];%绘图比较给定T值和计算T值plot(Tcal,'r*')hold onplot(T)结果比较如下，由比较可知两者值非常切合（在小数点后8位之后才有区别）：习题4-14充分发展区的温度控制方程如下：)(1rTr r r x T uc p ∂∂∂∂=∂∂λρ 对于三种无量纲定义w b w T T T T --=Θ、∞∞--=ΘT T T T w 、ww T T T T --=Θ∞进行分析如下1）由wb wT T T T --=Θ得：w w b T T T T +Θ-=)(由T 可得：x T x T x T T T x T w b w w b ∂∂Θ-+∂∂Θ=∂+Θ-∂=∂∂)1(])[(rT r T T r T T T r T w w b w w b ∂∂Θ-+∂Θ∂-=∂+Θ-∂=∂∂)1()(])[( 由b T 与r 无关、Θ与x 无关以及x T ∂∂、rT∂∂的表达式可知，除了w T 均匀的情况外，该无量纲温度定义在一般情况下是不能用分离变量法的； 2）由∞∞--=ΘT T T T w 得： ∞∞+Θ-=T T T T w )(由T 可得：xT x T T T x T w w ∂∂Θ=∂+Θ-∂=∂∂∞∞])[(rT r T T r T T T r T w w w ∂∂Θ+∂Θ∂-=∂+Θ-∂=∂∂∞∞∞)(])[( 由b T 与r 无关、Θ与x 无关以及x T ∂∂、rT∂∂的表达式可知，在常见的四种边界条件中除了轴向及周向均匀热流const q w =的情况外，有0=∂∂rT w，则该无量纲温度定义是可以用分离变量法的； 3）由wwT T T T --=Θ∞得： w w T T T T +Θ-=∞)(由T 可得：xT x T T T x T w w w ∂∂Θ-=∂+Θ-∂=∂∂∞)1(])[(rT r T T r T T T r T w w w w ∂∂Θ-+∂Θ∂-=∂+Θ-∂=∂∂∞∞)1()(])[( 同2）分析可知，除了轴向及周向均匀热流const q w =的情况外，有0=∂∂rT w，该无量纲温度定义是可以用分离变量法的；习题4－181）采用柱坐标分析，写出统一的稳态柱坐标形式动量方程：r r x x w r v r r r u x ∂∂+∂∂∂∂=∂∂+∂∂+∂∂1)()(1)(1)(φλφρθφρφρx 、r 和θ分别是圆柱坐标的3个坐标轴，u 、v 和w 管内的流动方向；对于管内的层流充分发展有：0=v 、0=w ，0=∂∂xu； 并且x 方向的源项：x pS ∂∂-=r 方向的源项：r pS ∂∂-=θ方向的源项：θ∂∂-=pr S 1由以上分析可得到圆柱坐标下的动量方程： x 方向：0)(1)(1=∂∂-∂∂∂∂+∂∂∂∂x pu r r r u r r r θλθλ r 方向：0=∂∂r pθ方向：0=∂∂θp边界条件： R r =，0=u0=r ，0=∂∂r u ；对称线上，0=∂∂θu不考虑液体的轴向导热，并简化分析可以得到充分发展的能量方程为：)(1)(1θλθλρ∂∂∂∂+∂∂∂∂=∂∂Tr r r T r r r x T uc p 边界条件： R r =，w q r T =∂∂λ；0=r ，0=∂∂rTπθ/0=，0=∂∂-θλT2）定义无量纲流速：dxdp R uU 2-=λ并定义无量纲半径：R r /=η；将无量纲流速和无量纲半径代入x 方向的动量方程得：0))1((1))1((122=∂∂-∂-∂∂∂+∂-∂∂∂xp U dx dp R R R R U dx dp R RR R θληλθηηλληηη 上式化简得：01)1(1)(1=+∂∂∂∂+∂∂∂∂θηθηηηηηU U 边界条件：1=η，0=U0=η，0=∂∂ηU；对称线上，0=∂∂θU 定义无量纲温度：λ/0R q T T b -=Θ其中，0q 是折算到管壁表面上的平均热流密度，即：Rq q wπ=0； 由无量纲温度定义可得：b T Rq T +Θ=λ将T 表达式和无量纲半径η代入能量方程得：)(1)(100θληλθηηλληηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂R q R R R R q R R R x T uc b p 化简得：)1(1)(10θηθηηηηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂x T u c q R b p （1）由热平衡条件关系可以得：mm m b m p b p p RU U q R u u R q A u u dx dT A u c x T u c x T uc 020221221)(===∂∂=∂∂ππρρρ 将上式代入式（1）可得：)1(1)(12θηθηηηηη∂Θ∂∂∂+∂Θ∂∂∂=m U U 边界条件：0=η，0=∂Θ∂η；1=η，R q q w πη10==∂Θ∂0=θ，0=∂Θ∂θ；πθ=，0=∂Θ∂θ单值条件： 由定义可知：0/0=-=ΘλR q T T b b b 且： ⎰⎰Θ=ΘAAb U d AU d A 即得单值性条件：0=Θ⎰⎰AA UdAUdA3）由阻力系数f 及Re 定义有：228)(2/Re ⎪⎭⎫ ⎝⎛=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=D D U D u u dx dp D f e m e m me νρ 且：m W b m W b m W R q T T D T T q Nu ,0,,0~2)/(2Θ=-=-=λλ5－21．一维稳态无源项的对流－扩散方程如下所示：xx u 22∂∂Γ=∂∂φφρ （取常物性）边界条件如下：L L x x φφφφ====,;,00上述方程的精确解如下：11)/(00--=--⋅PeL x Pe L e e φφφφ Γ=/uL Pe ρ 2．将L 分成20等份，所以有：∆=P Pe 201 2 3 4 5 6 ………… …………… 17 18 19 20 21 对于中心差分、一阶迎风、混合格式和QUICK 格式分别分析如下： 1） 中心差分中间节点： 2)5.01()5.01(11-∆+∆++-=i i i P P φφφ 20,2 =i2） 一阶迎风中间节点： ∆-∆++++=P P i i i 2)1(11φφφ 20,2 =i3） 混合格式当1=∆P 时，中间节点：2)5.01()5.01(11-∆+∆++-=i i i P P φφφ20,2 =i当10,5=∆P 时，中间节点： 1-=i i φφ 20,2 =i 4） QUICK 格式*12111)35(8122121⎥⎦⎤⎢⎣⎡---++++++=+--∆∆-∆∆+∆i i i i i i i P P P P P φφφφφφφ 2≠i *1111)336(8122121⎥⎦⎤⎢⎣⎡--++++++=+-∆∆-∆∆+∆i i i i i i P P P P P φφφφφφ 2=i数值计算结果与精确解的计算程序如下：%except for HS, any other scheme doesnt take Pe<0 into consideration %expression of exact solutiony=dsolve('a*b*Dy=c*D2y','y(0)=y0,y(L)=yL','x')y=subs(y,'L*a*b/c','t')y=simple(subs(y,'a*b/c*x','t*X'));ysim=simple(sym(strcat('(',char(y),'-y0)','/(yL-y0)')))y=sym(strcat('(',char(ysim),')*(yL-y0)','+y0'))% in the case of Pe=0y1=dsolve('D2y=0','y(0)=y0,y(L)=yL','x')y1=subs(y1,'-(y0-yL)/L*x','(-y0+yL)*X')%grid Pe numbertt=[1 5 10];%dimensionless lengthm=20;%mdim is the number of inner nodemdim=m-1;X=linspace(0,1,m+1);%initial value of variable during calculationy0=1;yL=2;%cal exact solutionfor n=1:size(tt,2)t=m*tt(1,n);if t==0yval1(n,:)=eval(y1);elseyval1(n,:)=eval(y);endend%extra treatment because max number in MATLAB is 10^308if max(isnan(yval1(:)))yval1=yval1';yval1=yval1(:);indexf=find(isnan(yval1));for n=1:size(indexf,1)if rem(indexf(n,1),size(X,2))==0yval1(indexf(n),1)=yL;elseyval1(indexf(n),1)=y0;endendyval1=reshape(yval1,size(X,2),size(yval1,1)/size(X,2));yval1=yval1';end%CD solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);b(n,:)=repmat([0.5*(1-0.5*t)],1,mdim);c(n,:)=repmat([0.5*(1+0.5*t)],1,mdim);d(n,1)=0.5*(1+0.5*tt(1,n))*y0;d(n,mdim)=0.5*(1-0.5*tt(1,n))*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval2=TDMA(a,b,c,d,mdim);yval2=[repmat([1],size(tt,2),1),yval2,repmat([2],size(tt,2),1)]; Fig(1,X,yval1,yval2,tt);title('CD Vs. Exact Solution')% FUS solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);b(n,:)=repmat([1/(2+t)],1,mdim);c(n,:)=repmat([(1+t)/(2+t)],1,mdim);d(n,1)=(1+tt(1,n))/(2+tt(1,n))*y0;d(n,mdim)=1/(2+tt(1,n))*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval3=TDMA(a,b,c,d,mdim);yval3=[repmat([1],size(tt,2),1),yval3,repmat([2],size(tt,2),1)]; Fig(2,X,yval1,yval3,tt);title('FUS Vs. Exact Solution')% HS solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);if t>2b(n,:)=repmat([0],1,mdim);c(n,:)=repmat([1],1,mdim);d(n,1)=y0;elseif t<-2b(n,:)=repmat([1],1,mdim);c(n,:)=repmat([0],1,mdim);d(n,mdim)=yL;elseb(n,:)=repmat([0.5*(1-0.5*t)],1,mdim);c(n,:)=repmat([0.5*(1+0.5*t)],1,mdim);d(n,1)=0.5*(1+0.5*t)*y0;d(n,mdim)=0.5*(1-0.5*t)*yL;endendc(:,1)=0;b(:,mdim)=0;% numerical cal by using TDMA subfuctionyval4=TDMA(a,b,c,d,mdim);yval4=[repmat([1],size(tt,2),1),yval4,repmat([2],size(tt,2),1)]; Fig(3,X,yval1,yval4,tt);title('HS Vs. Exact Solution')%QUICK Solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);b(n,:)=repmat([1/(2+t)],1,mdim);c(n,:)=repmat([(1+t)/(2+t)],1,mdim);d(n,1)=(1+tt(1,n))/(2+tt(1,n))*y0;d(n,mdim)=1/(2+tt(1,n))*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval5=zeros(size(tt,2),mdim);yval5com=yval5+1;counter=1;%iterativewhile max(max(abs(yval5-yval5com)))>10^-10if counter==1yval5com=TDMA(a,b,c,d,mdim);endfor nn=1:size(tt,2)for nnn=1:mdimif nnn==1d(nn,nnn)=((6*yval5com(nn,nnn)-3*y0-3*yval5com(nn,nnn+1))*tt(1,nn))/(8*(2+tt(1, nn)))+((1+tt(1,nn))/(2+tt(1,nn))*y0);elseif nnn==2d(nn,nnn)=((5*yval5com(nn,nnn)-3*yval5com(nn,nnn+1)-yval5com(nn,nnn-1)-y0)*tt (1,nn))/(8*(2+tt(1,nn)));elseif nnn==mdimd(nn,nnn)=((5*yval5com(nn,nnn)-3*yL-yval5com(nn,nnn-1)-yval5com(nn,nnn-2))*tt (1,nn))/(8*(2+tt(1,nn)))+(1/(2+tt(1,nn))*yL);elsed(nn,nnn)=((5*yval5com(nn,nnn)-3*yval5com(nn,nnn+1)-yval5com(nn,nnn-1)-yval5 com(nn,nnn-2))*tt(1,nn))/(8*(2+tt(1,nn)));endendendyval5=TDMA(a,b,c,d,mdim);temp=yval5;yval5=yval5com;yval5com=temp;counter=counter+1;endyval5=yval5com;yval5=[repmat([1],size(tt,2),1),yval5,repmat([2],size(tt,2),1)];Fig(4,X,yval1,yval5,tt);title('QUICK Vs. Exact Solution')%-------------TDMA SubFunction------------------function y=TDMA(a,b,c,d,mdim)%form a b c d resolve yval2 by using TDMA%eliminationp(:,1)=b(:,1)./a(:,1);q(:,1)=d(:,1)./a(:,1);for n=2:mdimp(:,n)=b(:,n)./(a(:,n)-c(:,n).*p(:,n-1));q(:,n)=(d(:,n)+c(:,n).*q(:,n-1))./(a(:,n)-c(:,n).*p(:,n-1));end%iterativey(:,mdim)=q(:,mdim);for n=(mdim-1):-1:1y(:,n)=p(:,n).*y(:,n+1)+q(:,n);end%-------------ResultCom SubFunction------------------ function y=ResultCom (a,b,c)for n=1:max(size(c,2))y(2*n-1,:)=a(n,:);y(2*n,:)=b(n,:);end%-------------Fig SubFunction------------------ function y=Fig(n,a,b,c,d)figure(n);plot(a,b);hold onplot(a,c,'*');str='''legend(';for n=1:size(d,2)if n==size(d,2)str=strcat(str,'''''Pe=',num2str(d(1,n)),''''')''');elsestr=strcat(str,'''''Pe=',num2str(d(1,n)),''''',');endendeval(eval(str));精确解与数值解的对比图，其中边界条件给定10=φ，2=L φ。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
一维稳态导热的数值计算
1.1物理问题
一个等截面直肋，处于温度t ∞=80
的流体中。
肋表面与流体之间的对流换热系数为h =
45W/(m 2∙℃),肋基处温度t w =300℃，肋端绝热。
肋片由铝合金制成，其导热系数为λ=110W/(m ∙℃)，肋片厚度为δ=0.01m ，高度为H=0.1m 。
试计算肋内的温度分布及肋的总换热量。
1.2数学描述及其解析解
引入无量纲过余温度θ
=
t−t ∞t w −t ∞
，则无量纲温度描述的肋片导热微分方程及其边界条
件：22
20d m dx
θθ-=
x=0,θ=θw =1 x=H,
0x
θ∂=∂ 其中 A
hp
m =
λ上述数学模型的解析解为：[()]
()()
w ch m x H t t t t ch mH ∞∞--=-⋅
()()w hp
t t th mH m
∞∅=
-
1.3数值离散
1.3.1区域离散
计算区域总节点数取N 。
1.3.2微分方程的离散
对任一借点i 有：22
20i d m dx θθ⎛⎫-= ⎪⎝⎭
用θ在节点i 的二阶差分代替θ在节点i 的二阶导数，得：211
2
20i i i i m x θθθθ+--+-=
整理成迭代形式：()1122
1
2i i i m x θθθ+-=++ (i=2,3……,N -1)
1.3.3边界条件离散
补充方程为：11w θθ==
右边界为第二类边界条件，边界节点N 的向后差分得：1
0N N x
θθ--=，将此式整理为
迭代形式，得：N 1N θθ-=
1.3.4最终离散格式
11w θθ==
()1122
1
2i i i m x θθθ+-=
++ (i=2,3……,N -1)
N 1N θθ-=
1.3.5代数方程组的求解及其程序
假定一个温度场的初始发布，给出各节点的温度初值：01θ，02θ，….，0
N θ。
将这些初
值代入离散格式方程组进行迭代计算，直至收敛。
假设第K 步迭代完成，则K+1次迭代计算式为：K 1
1
w θθ+=
()
11
11
2212i i K K K i m x
θθθ+-++=
++ (i=2,3……,N -1) 1
11N K K N θθ-++=
#include<stdio.h> #include<math.h> #define N 11 main() { int i;
float cha;/*cha 含义下面用到时会提到*/
float t[N],a[N],b[N];
float h,t1,t0,r,D,H,x,m,A,p; /*r 代表λ,x 代表Δx ，D 代表δ*/ printf("\t\t\t 一维稳态导热问题\t\t"); printf("\n\t\t\t\t\t\t----何鹏举\n");
printf("\n 题目：补充材料练习题一\n");
printf("已知：h=45，t1=80, t0=200, r=110, D=0.01, H=0.1 (ISO)\n"); /*下面根据题目赋值*/
h=45.0; t1=80.0; t0=300.0; r=110.0; D=0.01; H=0.1;
x=H/N; A=3.1415926*D*D/4; p=3.1415926*D; m=sqrt((h*p)/(r*A)); /*x 代表步长，p 代表周长，A 代表面积*/
printf("\n 请首先假定一个温度场的初始分布，即给出各节点的温度初值：\n");
{
scanf("%f",&t[i]);
a[i]=(t[i]-t1)/(t0-t1);
b[i]=a[i];/*这里b[i]用记录一下a[i]，后面迭代条件及二阶采用温度初场要用到*/ }
/*采用一阶精度的向后差分法数值离散*/
cha=1;
while(cha>0.0001)
{
a[0]=1;
for(i=1;i<N;i++)
a[i]=(a[i+1]+a[i-1])/(2+m*m*x*x);
a[N-1]=a[N-2];
cha=0;
for(i=0;i<N;i++)
cha=cha+a[i]-b[i];
cha=cha/N;/*cha代表每次迭代后与上次迭代各点温度差值的平均值*/
}
for(i=0;i<N;i++)
t[i]=a[i]*(t0-t1)+t1;
printf("\n\n经数值离散(一阶精度的向后差分法)计算得肋片的温度分布为:\n");
for(i=0;i<N;i++)
printf("%4.2f\t",t[i]);
printf("\n\n");
getchar();
/*采用二阶精度的元体平衡法数值离散(温度初值还用设定的初场，便于比较)*/ for(i=0;i<N;i++)
a[i]=b[i];
cha=1;
while(cha>0.0001)
{
a[0]=1;
for(i=1;i<N;i++)
a[i]=(a[i+1]+a[i-1])/(2+m*m*x*x);
a[N-1]=a[N-2]/(1+0.5*m*m*x*x);
cha=0;
for(i=0;i<N;i++)
cha=cha+a[i]-b[i];
cha=cha/N;
}
for(i=0;i<N;i++)
t[i]=a[i]*(t0-t1)+t1;
printf("\n\n经数值离散(二阶精度的元体平衡法)计算得肋片的温度分布为:\n");
printf("%4.2f\t",t[i]); printf("\n\n");
getchar();
}。