010_ThermalSystem.mw

熱システムのモデリング

イントロダクション

化学物質を作り出す場合、容器に含まれた液体の温度のコントロールを必要とします。
バッチ操作では、容器は通常、液体で満たされており定められた温度まで閉じた状態で暖められます。

このワークシートでは、液体の温度上昇をモデル化し必要な温度に達するまでの時間を解析します。

厳密または数値的な結果をもとめます。例えば方程式を解くや関数の変換、およびグラフ化を行います。

システムは以下の様に閉じた絶縁体の容器の中に液体が満たされており、電気ヒーターを通して加熱を行います。
:

[Inserted Image]

> restart;
Digits := 4;

Digits := 4

このモデルでは、

1. 容器と絶縁体の熱抵抗は RV であり、RV = 5.0*10^-3 s.K/J である。

2. 電気ヒーターは熱抵抗 RH の金属カバーに覆われており、RH = 1.0*10^-3 s.K/J である。

3. 液体とヒーターの熱容量はそれぞれ CL と CH となり、CL = 1.0*10^6 J/K、CH = 20.0*10^3 J/K である。

4. 液体とヒーターの周囲の初期温度は theta[a] となり、theta[a] = 300 K である。

5. t=0 で液体とヒーターの温度はそれぞれ theta[L] theta[H] で、液体の温度が一様になるよう混ぜられていると仮定します。

6. 時間 t=0 では、ヒーターは一定の割合でエネルギー qi(t) を供給するためにスイッチがオンになります。

この例の目的は、時間関数として液体、theta[L] の温度をモデル化して、必要な温度に達するのにかかる時間が365Kになるように見込みます。

システムモデルの設定

以下はヒーターと液体における温度の変化率を与えます:

> eqn1:= diff(theta[H](t),t) = (1/C[H])* (q[i](t) - q[H]):
eqn2:= diff(theta[L](t),t) = (1/C[L])* (q[H] - q[V]):

> q[H] := (theta[H](t) - theta[L](t))/R[H]:
q[V] := (theta[L](t) - theta[a])/R[V]:

C[H] C[L] R[H] R[V] theta[a] に数値をパラメータを代入します。

> C[H]:= 20E3: C[L]:=1E6: R[H]:= 1E-3: R[V]:= 5E-3: theta[a]:=300:
eqn1; eqn2;

diff(theta[H](t), t) = 0.5000e-4*q[i](t)-0.5000e-1*theta[H](t)+0.5000e-1*theta[L](t)

diff(theta[L](t), t) = 0.1e-2*theta[H](t)-0.1200e-2*theta[L](t)+0.6000e-1

以下を使用し、温度を絶対値ではなく、相対的価値 (theta[a] と比較できるよう) に変換します。

> eqn3:= theta[Hrel](t)= theta[H](t)-theta[a]:
eqn4:= theta[Lrel](t)= theta[L](t)-theta[a]:
eqn5:= q[irel](t)=q[i](t):

> subs({q[i](t)=q[irel], theta[H](t) = theta[Hrel](t)+300, theta[L](t) = theta[Lrel](t)+300}, eqn1):
eqn6:= simplify(%);

eqn6 := diff(theta[Hrel](t), t) = 0.5000e-4*q[irel]-0.5000e-1*theta[Hrel](t)+0.5000e-1*theta[Lrel](t)

> subs({theta[H](t)=theta[Hrel](t)+300,theta[L](t)=theta[Lrel](t)+300},eqn2):
eqn7:= simplify(%);

eqn7 := diff(theta[Lrel](t), t) = 0.1e-2*theta[Hrel](t)-0.1200e-2*theta[Lrel](t)

初期条件を定義します:

> theta[Hrel](0):=0: theta[Lrel](0):=0:

inttrans パッケージの laplace コマンドを使用し、微分方程式をラプラス変換します。

> with(inttrans):
laplace(eqn6, t, s);

s*laplace(theta[Hrel](t), t, s) = 0.5000e-4*q[irel]/s^1.-0.5000e-1*laplace(theta[Hrel](t), t, s)+0.5000e-1*laplace(theta[Lrel](t), t, s)

> subs({q[irel]=Q(s)*s, laplace(theta[Hrel](t),t,s)=Theta[H](s),
laplace(theta[Lrel](t),t,s)=Theta[L](s)}, %);

s*Theta[H](s) = 0.5000e-4*Q(s)-0.5000e-1*Theta[H](s)+0.5000e-1*Theta[L](s)

> Theta[H](s):=solve(%, Theta[H](s));

Theta[H](s) := 0.1000e-2*(Q(s)+1000.*Theta[L](s))/(20.*s+1.)

> laplace(eqn7, t, s):
subs({laplace(theta[Hrel](t),t,s)=Theta[H](s),
laplace(theta[Lrel](t),t,s)=Theta[L](s)}, %):
Theta[L](s):=solve(%, Theta[L](s));

Theta[L](s) := 0.5000e-2*Q(s)/(0.1000e6*s^2+5120.*s+1.)

上記の結果から伝達関数 H(s) がわかります。

> H(s):= Theta[L](s)/Q(s);

H(s) := 0.5000e-2/(0.1000e6*s^2+5120.*s+1.)

システム応答

今、伝達関数 H(s) がわかったので、一定の加熱をステップ関数で表します。

t > 0 において qirel(t) = 2*10^4 W を使用し、逆ラプラス変換を行い、時間による液体の温度が必要な温度に達するまでに掛かる時間がどのように依存するかを確認します。

> Theta[L](s):= H(s) * 2E4 *(1/s);

Theta[L](s) := 100.0/((0.1000e6*s^2+5120.*s+1.)*s)

> Theta[L](s) := convert(Theta[L](s) ,parfrac, s);

Theta[L](s) := .4000/(s+0.5100e-1)-100.4/(s+0.1961e-3)+100.0/s

> theta[Lrel](t):=invlaplace(%, s, t);

theta[Lrel](t) := .4000*exp(-0.5100e-1*t)-100.4*exp(-0.1961e-3*t)+100.

同等に:

> theta[L](t):=solve(eqn4,theta[L](t));

theta[L](t) := .4000*exp(-0.5100e-1*t)-100.4*exp(-0.1961e-3*t)+400.

液体の温度を時間の関数としてプロットします。

> with(plots):

Warning, the name changecoords has been redefined

> aplot:=plot([theta[L](t), op(3,theta[L](t)), 365], t=0..1.5e4,
 temp=280..430, labels=[`Time (secs)`,`Temp (K)`],
 linestyle=[1,2,3], thickness=2, axes=boxed,
 title=`Liquid Temperature vs. Time`):

> bplot:=textplot({[1000, 400,`Equilibrium Temperature of Liquid`],
 [1000,365,`Desired Temperature`]},align={ABOVE,RIGHT}):

> display([aplot,bplot]);

[Plot]

グラフよりヒーターの出力によって、必要な液体の温度 365 K に達するまでのおおよその時間がわかります。続いて、その時の時間 t を求めます。

しかし、厳密に二重指数を定めることは出来ないので、時間に係数の tau1 = 1/0.0510 = 19.6s と tau2 = 1/0.000196 = 5102s の 2 つの指数の相対的なサイズについて考えます。

> evalf(subs(t=10, op(1,theta[L](t)) ));
evalf(subs(t=10, op(2,theta[L](t)) ));

.2402

-100.2

上記から桁数が 3 桁異なるのがわかります。したがって、より小さい用語を無視します、近似をそれにして:

> theta[L](t) := -101.*exp(-.196e-3*t)+400.;

theta[L](t) := -101.*exp(-0.196e-3*t)+400.

液体の温度が 365 K に達する時間を求めます。

> time_secs:=solve(theta[L](t)=365, t);
time_mins:= time_secs/60;

time_secs := 5407.

time_mins := 90.12

結論

このワークシートでは、一般的な化学工学問題をモデル化するのにおいて Maple を使用して解を求めるプロセスを示しました。

Maple の数式処理、数値解析、グラフィック機能は、多くのシステムの理想的なモデル化を行うためのツールとして利用できます