薄板非线性传热

打开生活的脚本

这个例子展示了如何对薄板进行传热分析。

板是方形的，沿底部边缘的温度是固定的。没有热量从其他三个边缘转移(即它们是绝缘的)。热量通过对流和辐射从板的顶部和底部传递。因为辐射也包括在内，所以这个问题是非线性的。本例的目的之一是演示如何处理PDE问题中的非线性。

同时进行了稳态和瞬态分析。在稳态分析中，我们感兴趣的是板中不同点在达到平衡状态后的最终温度。在暂态分析中，我们感兴趣的是板内温度作为时间的函数。这个瞬态分析可以回答的一个问题是，板块达到平衡温度需要多长时间。

平板传热方程

平板的平面尺寸为1米乘1米，厚度为1厘米。由于与平面尺寸相比，平板相对较薄，因此可以假设温度在厚度方向上是恒定的;由此产生的问题是二维的。

对流和辐射换热被假定发生在板的两个面和特定的环境温度之间。

每单位面积上每块板面由于对流而传递的热量定义为

$问_{c} ＝ h_{c} （ T - T_{一个} ）$

在哪里 $T_{一个}$ 是环境温度， $T$ 是在平板表面特定x和y位置的温度，和 $h_{c}$ 是指定的对流系数。

单位面积上每块板面因辐射而传递的热量定义为

$问_{r} ＝ ϵ σ （ T^{4} - T_{一个}^{4} ）$

在哪里 $ϵ$ 是面和的发射率 $σ$ 是Stefan-Boltzmann常数。由于辐射传递的热量与表面温度的四次方成正比，这个问题是非线性的。

描述薄板内温度的偏微分方程为

$ρ C_{p} t_{z} \frac{\partial T}{\partial t} - k t_{z} \nabla^{2} T + 2 问_{c} + 2 问_{r} ＝ 0$

在哪里 $ρ$ 是材料密度， $C_{p}$ 是比热容， $t_{z}$ 是板厚，两个因素的传热从两个板面。

用PDE工具箱所期望的形式重写这个方程是很方便的

$ρ C_{p} t_{z} \frac{\partial T}{\partial t} - k t_{z} \nabla^{2} T + 2 h_{c} T + 2 ϵ σ T^{4} ＝ 2 h_{c} T_{一个} + 2 ϵ σ T_{一个}^{4}$

问题的设置

该板由铜组成，具有以下特性:

k = 400;铜的热导率，W/(m-K)ρ= 8960;铜的%密度，kg/m^3specificHeat = 386;%铜比热，J/(kg-K)厚= . 01;%板厚，单位为米stefanBoltz = 5.670373 e-8;% Stefan-Boltzmann常数，W/(m^2-K^4)hCoeff = 1;%对流系数，W/(m²- k)假设环境温度为300开尔文。。ta = 300;工作= 5;板表面的%发射率

创建具有单个因变量的PDE模型。

numberOfPDE = 1;模型= createpde (numberOfPDE);

对于一个正方形，几何形状和网格很容易定义，如下所示。

宽度= 1;身高= 1;

通过给出4个x坐标和4个角的y坐标来定义正方形。

GDM =[3 4 0宽度宽度0 0 0高度高度]';g = decsg (gdm,“S1 ', (“S1 ') ');

将DECSG几何体转换为一个几何体对象，然后将其附加到PDEModel中

geometryFromEdges(模型中,g);

绘制几何图形并显示边缘标签，以便在边界条件定义中使用。

图;pdegplot(模型,“EdgeLabels”，“上”）;轴([-。1.1 - 1。1 1.1]);标题(显示边缘标签的几何图形）

图中包含一个axes对象。标题为Geometry with Edge Labels的axes对象包含5个类型为line、text的对象。

指定系数。通过将上面的方程与PDE工具箱文档中的标量抛物方程进行比较，可以很容易地确定PDE工具箱所需系数的表达式。

c =厚* k;

由于辐射边界条件的存在，“a”系数是温度u的函数。它被定义为MATLAB表达式，因此在分析过程中可以对u的不同值求值。

a = @(~，state) 2*hCoeff + 2*emiss*stefanBoltz*state.u.^3;f = 2*hCoeff*ta + 2*emiss*stefanBoltz*ta^4;d = *ρ* specificHeat厚;specifyCoefficients(模型,“m”0,“d”0,“c”c“一个”一个,“f”f);

板的底部边缘设置为1000开氏度。

应用边界条件。板的三个边缘是绝缘的。由于诺伊曼边界条件等于零是有限元公式中的默认值，这些边的边界条件不需要显式设置。对底边1上的所有节点设置狄利克雷条件，

applyBoundaryCondition(模型,“边界条件”，“边缘”, 1“u”, 1000);

指定初始猜测。

setInitialConditions(模型中,0);

在正方形上创建三角形网格，每个方向大约有十个元素。

hmax = 1;%的元素大小msh = generateMesh(模型,“Hmax”, hmax);图;pdeplot(模型);轴平等的标题(三角网格板)包含(“x坐标,米”) ylabel (“坐标,米”）

图中包含一个axes对象。标题为Plate with triangle Element Mesh的axis对象包含两个类型为line的对象。

稳态解

因为a和f系数是温度的函数(由于辐射边界条件)，solvepde自动选取非线性求解器求解。

R = solvepde(模型);u = R.NodalSolution;图;pdeplot(模型,“XYData”u“轮廓”，“上”，“ColorMap”，“喷气机”）;标题(“板内温度，稳态解”)包含(“x坐标,米”) ylabel (“坐标,米”)轴平等的

图中包含一个axes对象。标题为Temperature In The Plate, Steady State Solution的axis对象包含12个类型为patch, line的对象。

p = msh.Nodes;plotAlongY (p, u, 0);标题(“温度与y坐标的函数”)包含(“坐标,米”) ylabel (“温度,开尔文”）

图中包含一个axes对象。标题为Temperature作为y坐标的函数的axes对象包含一个类型为line的对象。

流(['板上边缘的温度='.．.' % 5.1 f degrees-K \ n”), u (4));

板顶边缘的温度= 449.8度- k

临时的解决方案

包括d系数。

specifyCoefficients(模型,“m”0,“d”d“c”c“一个”一个,“f”f);endTime = 5000;tlist = 0:50: endTime;numNodes =大小(p, 2);

将所有节点的初始温度设置为室温，300k。

情况(1:numNodes) = 300;

设置下边缘E1的初始温度为常量BC, 1000k。

1000年setInitialConditions(模型,“边缘”1);

设置以下求解器选项。

model.SolverOptions.RelativeTolerance = 1.0 e - 3;model.SolverOptions.AbsoluteTolerance = 1.0的军医;

通过使用解决问题solvepde．求解器自动选取抛物线求解器求解。

R = solvepde(模型、tlist);u = R.NodalSolution;图;:情节(tlist u (3));网格在标题([“沿顶部边缘的温度”.．.“板块作为时间的函数”])包含(“时间,秒”) ylabel (“温度,开尔文”）

图中包含一个axes对象。标题为“沿板上边缘的温度作为时间的函数”的axis对象包含一个类型为line的对象。

图;pdeplot(模型,“XYData”u(:,结束),“轮廓”，“上”，“ColorMap”，“喷气机”）;标题(sprintf ([“板块温度”，.．.'暂态解(%d秒)\n'), tlist (1)));包含(“x坐标,米”) ylabel (“坐标,米”)轴平等的；

图中包含一个axes对象。标题为Temperature In The Plate,Transient Solution(5000秒)的axis对象包含12个类型为patch, line的对象。

流(['\n温度在顶部边缘(t = %5.1f秒)= '.．.' % 5.1 f degrees-K \ n”), tlist(结束),u (4));

顶部边缘的温度(t = 5000.0秒)= 441.8度- k

总结

稳态解和结束时刻瞬态解得到的板内温度曲线非常接近。也就是说，在大约5000秒后，瞬态解已经达到稳态值。平板上边缘的两种溶液的温度一致在1%以内。