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

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

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

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

この例では、Symbolic Math Toolbox™ を使用して非線形 PDE をシンボリックに表現し、Partial Differential Equation Toolbox™ の有限要素解析を使用してこの PDE の問題を解く方法について説明します。この例では、過渡解析を行い、板の温度を時間の関数として求めます。この過渡解析では、定常状態の板が平衡温度に達するまでの時間が示されます。

板の熱伝導方程式

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

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

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

$Q_{c} = h_{c} (T - T_{a})$ ,

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

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

$Q_{r} = ϵ σ (T^{4} - T_{a}^{4})$ ,

ここで、 $ϵ$ は板表面の放射率、 $σ$ はステファン・ボルツマン定数です。放射で伝わる熱量は表面温度の 4 乗に比例するため、この問題は非線形となります。

この薄板の温度を表す PDE は以下のようになります。

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

ここで、 $ρ$ は板の材質密度、 $C_{p}$ は板の比熱、 $t_{z}$ は板の厚さ、 $k$ は板の熱伝導率、そして 2 つの因子は板の各面の熱伝導を表します。

PDE パラメーターの定義

パラメーターの値を定義して PDE の問題を設定します。この板は銅でできており、以下の特性を備えています。

kThermal = 400;             % thermal conductivity of copper, W/(m-K)
rhoCopper = 8960;           % density of copper, kg/m^3
specificHeat = 386;         % specific heat of copper, J/(kg-K)
thick = 0.01;               % plate thickness in meters
stefanBoltz = 5.670373e-8;  % Stefan-Boltzmann constant, W/(m^2-K^4)
hCoeff = 1;                 % convection coefficient, W/(m^2-K)
tAmbient = 300;             % the ambient temperature
emiss = 0.5;                % emissivity of the plate surface

PDE 係数の抽出

板の温度を従属変数 $T (t, x, y)$ としてもつシンボリック形式で PDE を定義します。

syms T(t,x,y)
syms eps sig tz hc Ta rho Cp k
Qc = hc*(T - Ta);
Qr = eps*sig*(T^4 - Ta^4);
pdeeq = (rho*Cp*tz*diff(T,t) - k*tz*laplacian(T,[x,y]) + 2*Qc + 2*Qr)

pdeeq(t, x, y) = 
 $2 eps sig ({T (t, x, y)}^{4} - {Ta}^{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 (Ta - T (t, x, y)) + Cp ρ tz \frac{\partial}{\partial t} T (t, x, y)$

次に、Partial Differential Equation Toolbox で使用できるように、PDE モデルの入力として与える係数を作成します。これを行なうには、まず、関数 pdeCoefficients を使用して、シンボリック PDE の係数をシンボリック式の構造体として抽出します。

symCoeffs = pdeCoefficients(pdeeq,T,'Symbolic',true)

symCoeffs = struct with fields:
    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 はシンボリックオブジェクトであるため、Partial Differential Equation Toolbox への有効な入力となるように、関数 pdeCoefficientsToDouble を使用して係数を double データ型に変換します。

coeffs = pdeCoefficientsToDouble(symCoeffs)

coeffs = struct with fields:
    a: @makeCoefficient/coefficientFunction
    c: [4x1 double]
    m: 0
    d: 3.4586e+04
    f: 1.0593e+03

PDE モデル、幾何形状、および係数の指定

次に、Partial Differential Equation Toolbox を使用し、これらの係数に基づいて有限要素解析で PDE の問題を解きます。

まず、1 つの従属変数をもつ PDE モデルを作成します。

numberOfPDE = 1;
model = createpde(numberOfPDE);

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

width = 1; 
height = 1;

幾何形状を表す行列を定義します。関数decsg (Partial Differential Equation Toolbox)を使用して、正方形の幾何形状を作成します。矩形形状の場合、最初の行に 3 を格納し、2 番目の行に 4 を格納します。次の 4 つの行には辺の始点を表す $x$ 座標を格納し、その次の 4 つの行には辺の始点を表す $y$ 座標を格納します。

gdm = [3 4 0 width width 0 0 0 height height]';
g = decsg(gdm,'S1',('S1')');

DECSG 幾何形状を幾何形状オブジェクトに変換し、このオブジェクトを PDE モデルに組み込みます。

geometryFromEdges(model,g);

幾何形状をプロットし、辺のラベルを表示します。

figure; 
pdegplot(model,'EdgeLabels','on'); 
axis([-0.1 1.1 -0.1 1.1]);
title('Geometry with Edge Labels Displayed');

Figure contains an axes object. The axes object with title Geometry with Edge Labels Displayed contains 5 objects of type line, text.

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

hmax = 0.1; % element size
msh = generateMesh(model,'Hmax',hmax);
figure; 
pdeplot(model); 
axis equal
title('Plate with Triangular Element Mesh');
xlabel('X-coordinate, meters');
ylabel('Y-coordinate, meters');

Figure contains an axes object. The axes object with title Plate with Triangular Element Mesh, xlabel X-coordinate, meters, ylabel Y-coordinate, meters contains 2 objects of type line.

PDE モデルの係数を指定します。

specifyCoefficients(model,'m',coeffs.m,'d',coeffs.d, ...
    'c',coeffs.c,'a',coeffs.a,'f',coeffs.f);

過渡解の算出

境界条件を適用します。板の 3 辺は断熱されています。有限要素法の定式化では、ノイマン境界条件が既定で 0 となっているため、これらの辺に境界条件を設定する必要はありません。板の下辺は 1,000 K に固定されています。これを指定するには、下辺 (辺 E1) のすべてのノードにディリクレ条件を適用します。

applyBoundaryCondition(model,'dirichlet','edge',1,'u',1000);

下辺を除き、すべての場所の初期温度を 0 K に指定します。下辺 E1 の初期温度を、一定の境界条件の値 (1,000 K) に設定します。

setInitialConditions(model,0);
setInitialConditions(model,1000,'edge',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.NodalSolution;
figure; 
plot(tlist,u(3,:));
grid on
title 'Temperature Along the Top Edge of the Plate as a Function of Time'
xlabel 'Time (s)'
ylabel 'Temperature (K)'

Figure contains an axes object. The axes object with title Temperature Along the Top Edge of the Plate as a Function of Time, xlabel Time (s), ylabel Temperature (K) contains an object of type line.

プロットを見ると、6,000 秒後に過渡解が定常状態に達しつつあるのが分かります。6,000 秒後、上辺の平衡温度が 450 K に近づきます。

10,000 秒後の板の温度プロファイルを表示します。

figure;
pdeplot(model,'XYData',u(:,end),'Contour','on','ColorMap','jet');
title(sprintf('Temperature in the Plate, Transient Solution (%d seconds)\n', ...
  tlist(1,end)));
xlabel 'X-coordinate, meters'
ylabel 'Y-coordinate, meters'
axis equal;

Figure contains an axes object. The axes object with xlabel X-coordinate, meters, ylabel Y-coordinate, meters contains 12 objects of type patch, line.

10,000 秒後の上辺の温度を表示します。

u(3,end)

ans = 449.7950

参考

関数

pdeCoefficients | pdeCoefficientsToDouble | subs