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

这个例子使用了:

打开实时脚本

这个例子展示了如何求解薄板非线性传热的偏微分方程。

板是方形的，其温度沿底部边缘固定。没有热量从其他三个边缘转移，因为边缘是绝缘的。热量通过对流和辐射从板的顶部和底部传递。因为辐射是包含在内的，所以这个问题是非线性的。

本例的目的是展示如何使用符号数学工具箱™以符号方式表示非线性PDE，并使用有限元分析在偏微分方程工具箱™。在本例中，进行瞬态分析并求解板内温度随时间的函数。瞬态分析显示了在稳定状态下，板达到平衡温度所需要的时间。

平板的传热方程

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

假设在一定环境温度下，板的两面与环境之间发生对流和辐射换热。

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

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

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

每单位面积上由于辐射而从每个板面传递的热量定义为

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

在哪里 $ϵ$ 面和的发射率 $σ$ 为斯特凡-玻尔兹曼常数。由于辐射传递的热量与表面温度的四次幂成正比，所以这个问题是非线性的。

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

$ρ 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;铜的密度%，kg/m^3specificHeat = 386;铜比热% J/(kg-K)厚度= 0.01;%板厚，单位为米stefanBoltz = 5.670373e-8;% Stefan-Boltzmann常数，W/(m^2-K^4)hCoeff = 1;%对流系数W/(m^2-K)tAmbient = 300;%环境温度Emiss = 0.5;板表面的%发射率

提取PDE系数

以板温为因变量，用符号形式定义偏微分方程 $T （ t ， x ， y ）$ ．

信谊T (T, x, y)信谊每股收益团体tzhc助教ρCpkQc = hc*(T - Ta);Qr = eps*sig*(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模型输入的系数。方法将符号PDE的系数提取为符号表达式的结构pdeCoefficients函数。

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

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

接下来，将表示PDE参数的符号变量替换为它们的数值。

symVars = [eps sig tz hc Ta rho Cp k];symVals = [emiss stefanBoltz thick hCoeff tAmbient rhoCopper specificHeat kThermal];symCoeffs = subs(symCoeffs,symVars,symVals);

最后，由于字段symCoeffs是符号对象，使用pdeCoefficientsToDouble函数将系数转换为双数据类型，这使它们成为偏微分方程工具箱的有效输入。

coeffs = pdeCoefficientsToDouble(symCoeffs)

多项式系数=带字段的结构:a: @makeCoefficient/ coefficients function c: [4x1 double] m: 0 d: 3.4586e+04 f: 1.0593e+03

指定PDE模型、几何形状和系数

现在，利用偏微分方程工具箱，利用这些系数进行有限元分析，求解偏微分方程问题。

首先，创建带有单个因变量的PDE模型。

numberOfPDE = 1;model = 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]);标题(显示边缘标签的几何图形）;

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

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

Hmax = 0.1;元素尺寸百分比msh = generateMesh(模型，“Hmax”, hmax);图;pdeplot(模型);轴平等的标题(“三角网格板”）;包含(“x坐标,米”）;ylabel (“坐标,米”）;

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

指定PDE模型中的系数。

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

寻找瞬态解

应用边界条件。三个板边是绝缘的。因为诺伊曼边界条件等于零是有限元公式中的默认值，所以不需要在这些边上设置边界条件。板底边固定在1000k。在底边(边E1)上的所有节点上使用Dirichlet条件指定这一点。

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

指定初始温度为0 K，底边除外。设置底边E1的初始温度为常数边界条件的值，1000k。

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

定义时域，求偏微分方程问题的暂态解。

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

设置解算器选项的公差。

model.SolverOptions.RelativeTolerance = 1.0e-3;model.SolverOptions.AbsoluteTolerance = 1.0e-4;

解决问题的方法solvepde．画出沿板上边缘的温度随时间的函数。

R = solvepde(model,tlist);u = r . nodesolution;图;:情节(tlist u (3));网格在标题“沿板上边缘的温度随时间的函数”包含“时间(s)”ylabel“温度(K)”

图中包含一个轴对象。标题为Temperature Along The Top Edge of The Plate as a Function of Time的坐标轴对象包含一个line类型的对象。

根据图，6000秒后暂态解开始达到稳态值。6000秒后，顶部边缘的平衡温度接近450k。

显示10,000秒后板材的温度分布。

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

图中包含一个轴对象。标题为Temperature in The Plate, Transient Solution(10000秒)的轴对象包含12个类型为patch, line的对象。

显示顶部边缘10,000秒的温度。

u (3)

Ans = 449.8031