解偏微分方程的非线性传热

这个示例使用:

打开生活的脚本

这个例子展示了如何解决偏微分方程(PDE)薄板非线性传热。

板是广场,其温度是固定的底部边缘。没有热量转移从其他三个边缘边缘是绝缘的。热转移从顶部和底部的面孔板通过对流和辐射。因为辐射包括,问题是非线性的。

这个例子的目的是展示如何表示的非线性PDE象征性地使用符号数学工具箱™和解决使用PDE问题有限元分析在偏微分方程的工具箱™。在这个例子中,进行瞬态分析和解决温度板作为时间的函数。瞬态分析显示时间,直到板温度达到平衡稳定状态。

板的传热方程

板的平面尺寸1 m×1米,厚1厘米。由于板相对薄与平面尺寸相比,温度可以假定为常数在厚度方向上,以及由此产生的问题是二维的。

假设对流和辐射热量转移之间发生的两副面孔板与一个指定的环境温度和环境。

从每个盘子的热量转移面临由于对流被定义为每单位面积

$问_{c} = h_{c} (T - - - - - - T_{一个})$ ,

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

从每个盘子的热量转移面临由于辐射的定义是单位面积上

$问_{r} = ϵ σ (T^{4} - - - - - - T_{一个}^{4})$ ,

在哪里 $ϵ$ 脸的发射率和吗 $σ$ 是斯蒂芬玻尔兹曼常数。因为由于辐射热量正比于表面温度的四次方,问题是非线性的。

PDE描述这个薄板的温度

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

在哪里 $ρ$ 材料的密度板, $C_{p}$ 是它的比热容, $t_{z}$ 是它的板厚, $k$ 是它的导热系数和传热的因素两个帐户的面孔。

PDE参数定义

建立了PDE问题通过定义参数的值。板是由铜,具有以下特性。

kThermal = 400;%导热铜,W / (m k)rhoCopper = 8960;^ %密度的铜,公斤/米3specificHeat = 386;%铜的比热,J / (kg-K)厚= 0.01;%在米板厚度stefanBoltz = 5.670373 e-8;%斯蒂芬玻尔兹曼常数,W / (m ^ 2 k ^ 4)hCoeff = 1;%对流系数,W / (m ^ 2 k)tAmbient = 300;%的环境温度工作= 0.5;%板表面的发射率

提取PDE系数

定义PDE符号形式与板温度作为因变量 $T (t, x, y)$ 。

信谊T (T, x, y)信谊每股收益团体tzhc助教ρCpkQc = hc * (T - Ta);Qr =每股收益*团体* (T ^ 4 - Ta ^ 4);pdeeq =(ρ* Cp * tz *差异(T, T) - k * tz *拉普拉斯算子(T (x, y)) + 2 * Qc + 2 * Qr)

pdeeq (t, x, y) =
                  
                    $2 每股收益 团体 ({T (t, x, y)}^{4} - - - - - - {助教}^{4}) - - - - - - k tz (\frac{\partial^{2}}{\partial x^{2}} T (t, x, y) + \frac{\partial^{2}}{\partial y^{2}} T (t, x, y)) - - - - - - 2 hc (助教 - - - - - - T (t, x, y)) + Cp ρ tz \frac{\partial}{\partial t} T (t, x, y)$

现在,创建PDE模型中的系数作为输入根据偏微分方程的工具箱。要做到这一点,首先提取系数的象征性符号表达式的结构使用PDEpdeCoefficients函数。

symCoeffs = pdeCoefficients (pdeeq T“象征”,真正的)

symCoeffs =结构体字段:m: 0: 2 * hc + 2 * eps *团体* T (T, x, y) ^ 3 c: [2 x2 sym] f: 2 * eps *团体* Ta ^ 4 + 2 * hc * Ta d: Cp *ρ* tz

接下来,替代符号变量代表PDE参数数值。

symVars = (eps团体tz hc TaρCp k);symVals =[工作stefanBoltz厚hCoeff tAmbient rhoCopper specificHeat kThermal);symCoeffs =潜艇(symCoeffs symVars symVals);

最后,因为这一领域symCoeffs是符号对象,使用pdeCoefficientsToDouble转换系数的函数双数据类型,这使他们有效的输入为偏微分方程的工具箱。

多项式系数= pdeCoefficientsToDouble (symCoeffs)

多项式系数=结构体字段:@makeCoefficient / coefficientFunction c: [4 x1双)m: 0 d: 3.4586 e + 04 f: 1.0593 e + 03

指定PDE模型、几何和系数

现在,使用偏微分方程的工具箱,解决PDE问题使用基于这些系数的有限元分析。

首先,创建PDE模型与一个因变量。

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

指定PDE的几何模型在这种情况下,广场的维度。

宽度= 1;身高= 1;

定义一个几何描述矩阵。创建几何广场使用decsg(偏微分方程工具箱)函数。对于一个矩形几何,第一行包含3第二行包含4。接下来的四行包含 $x$ 坐标起始点的边缘,和包含的四行之后 $y$ 坐标起始点的边缘。

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

DECSG几何转换成一个几何对象并将其包括在PDE模型。

geometryFromEdges(模型中,g);

绘制几何图形和显示标签。

图;pdegplot(模型,“EdgeLabels”,“上”);轴([-0.1 1.1 -0.1 1.1]);标题(与边缘的几何标签显示的);

图包含一个坐标轴对象。坐标轴对象与标题几何边缘标签显示包含5线类型的对象,文本。

接下来,创建PDE的三角网格模型的网格大小大约0.1在每个方向。

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

图包含一个坐标轴对象。坐标轴对象与网格三角形元素包含标题板2线类型的对象。

指定PDE模型中的系数。

specifyCoefficients(模型,“米”coeffs.m,' d 'coeffs.d,…“c”coeffs.c,“一个”coeffs.a,“f”,coeffs.f);

找到临时解决方案

应用边界条件。三个板边是绝缘的。因为诺伊曼边界条件等于零是默认的有限元公式,你不需要设置这些边缘的边界条件。底部的边缘板固定在1000 K。指定使用狄利克雷条件在所有节点在底部边缘(边缘E1)。

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

指定初始温度为0 K无处不在,除了在底部边缘。设置初始温度对底部边缘E1的值常数边界条件,1000 K。

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

定义了瞬态时域找到PDE问题的解决方案。

endTime = 10000;tlist = 0:50: endTime;

设置公差的解决方案。

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

解决问题用solvepde。情节的温度在顶部边缘板作为时间的函数。

R = solvepde(模型、tlist);u = R.NodalSolution;图;:情节(tlist u (3));网格在标题板的顶部温度作为时间的函数的包含“时间(s)”ylabel“温度(K)”

图包含一个坐标轴对象。坐标轴对象标题温度在顶部的边缘板作为时间的函数包含一个类型的对象。

根据情节,瞬态解开始6000秒后达到稳态值。顶部边缘方法的平衡温度450 K后6000秒。

显示板的温度曲线后10000秒。

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

图包含一个坐标轴对象。坐标轴对象与标题板温度,瞬态解(10000秒)包含12块类型的对象。

顶部边缘显示温度在10000秒。

u (3)

ans = 449.8031