非線形熱伝導に関する偏微分方程式の求解

この例では次を使用します:

この例では，薄板の非線形熱伝導に関する偏微分微分方程式(pde)を解く方法を説明します。

この板は正方形で，下辺の温度は一定です。他の3辺は断熱されており，これらの辺から熱は伝わりません。熱は，板の表裏両面から対流と放射によって伝わります。放射が含まれているため，この問題は非線形となります。

この例では符号数学工具箱™を使用して非線形PDEをシンボリックに表現し,偏微分方程工具箱™の有限要素解析を使用してこのpdeの問題を解く方法にいて説明します。この例では，過渡解析を行い，板の温度を時間の関数として求めます。この過渡解析では，定常状態の板が平衡温度に達するまでの時間が示されます。

板の熱伝導方程式

板の面積は1米× 1米で，厚さは1厘米です。この板は,面積と比べて厚さが薄いため,厚さ方向の温度は一定とみなすことができ,結果として2次元の問題となります。

対流と放射による熱伝導は，板の両面，および指定された周囲温度にある環境との間で起こるものとします。

対流によって板の各面から伝達される単位面積あたりの熱量は，次の式で与えられます。

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

ここで， $T_{一个}$ は周囲温度， $T$ は板表面の特定の $x$ 座標および $y$ 座標における温度， $h_{c}$ は指定された対流係数です。

放射によって板の各面から伝達される単位面積あたりの熱量は，次の式で与えられます。

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

ここで， $ϵ$ は板表面の放射率， $σ$ はステファン·ボルマン定数です。放射で伝わる熱量は表面温度の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$ は板の熱伝導率，そして2の因子は板の各面の熱伝導を表します。

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 ）$ としてもシンボリック形式でpdeを定義します。

信谊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モデルの入力として与える係数を作成します。これを行なうには，まず，関数pdeCoefficientsを使用して，シンボリックpdeの係数をシンボリック式の構造体として抽出します。

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の問題を解きます。

まず，1の従属変数をもpde。

numberOfPDE = 1;model = createpde(numberOfPDE);

Pdeモデルの幾何形状を指定します。この例では，正方形の寸法を指定します。

宽度= 1;高度= 1;

幾何形状を表す行列を定義します。関数decsg(偏微分方程工具箱)を使用して，正方形の幾何形状を作成します。矩形形状の場合，最初の行に3.を格納し，2番目の行に4を格納します。次の4の行には辺の始点を表す $x$ 座標を格納し，その次の4の行には辺の始点を表す $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的对象。

次に,各方向のメッシュサイズが約0.1となるように,PDEモデル内に三角形のメッシュを作成します。

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);

過渡解の算出

境界条件を適用します。板の3辺は断熱されています。有限要素法の定式化では,ノイマン境界条件が既定で0となっているため,これらの辺に境界条件を設定する必要はありません。板の下辺は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.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秒後，上辺の平衡温度が450 kに近づきます。

1万秒後の板の温度プロファ电子邮箱ルを表示します。

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

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

10000秒後の上辺の温度を表示します。

u (3)

Ans = 449.8031