015_HeatingLiquid.mw

液体の加熱

イントロダクション

このワークシートでは、時間に依存する発熱体によって加熱される水を含む壁の薄い容器について考えます。
ヒーターがスイッチをオフにすると、水温は40℃未満になります。

[Inserted Image]

> restart;

問題の定義

容器には、表面 A [m] と熱伝達率 k [W/(m^2*K) ] があり、容量 M は [kg] で表され、比熱 c は [J/(kg*K) ] となります。

液体の温度は T0、大気の温度は Tu です。

ヒーターのパワーは、[W] の最大加熱パワーの Q と不変のτで与えられます。

パラメータの値は以下の通りです。

> val:=k=20,M=1,c=4182,A=0.5,T0=20,Tu=20:
tmax:=200:

発熱体のモデル化

まず最初に、発熱体をモデル化します:

以下はヒーターのパワーを与えます。

> q[on] := Q*(1-exp(-tau*t));

q[on] := Q*(1-exp(-tau*t))

[W] の最大加熱パワーが Q とτが不変の場合、ヒーターが最大の力に達した後に、それは t_off で離れます。

熱の容量は液体が暖まり続けるのを意味し、オフにしている状態では、以下はパワーを与えます。

> q[off] := Q*exp(-tau*(t-t_off));

q[off] := Q*exp(-tau*(t-t_off))

の時、

> valQ := tau = 1/10, Q=1000;

valQ := tau = 1/10, Q = 1000

発熱体のパワーモデルを区分関数で表します。

> Q_ := piecewise(t <t_off,q[on], q[off]);

Q_ := PIECEWISE([Q*(1-exp(-tau*t)), t < t_off], [Q*exp(-tau*(t-t_off)), otherwise])

> Q_:=unapply(subs(valQ,Q_),(t,t_off)):

グラフから発熱体のパワーは 100 秒以降でオフになると仮定すします。

> plot(Q_(t,100),t = 0 .. 200, title="Power of the Heat Source",
labels=["Time [s]","Power [W] "] );

[Plot]

理論

液体の全ての点での温度が同じように、もしミキサーを使用した状況で時間だけに依存すると仮定すると、容器の熱容量は無視されます。

液体の熱容量が容器の熱容量より高いと想定するとこの仮定は有効になります。

容器の温度を表す微分方程式は下記の通りです。

> deq := diff(T(t),t) = -(k*A*(T(t)-Tu))/(M*c) + q/(M*c);

deq := diff(T(t), t) = -k*A*(T(t)-Tu)/(M*c)+q/(M*c)

右辺の1項目は容器の表面を通る熱の損失を表し、2項目は熱を表します。

解決法

初期条件 T(0) = T0 を使用し加熱中の状態を表す微分方程式を解きます。

そして、解を T_off に割り当てます。

> sol1:=dsolve({subs(q=q[on],deq),T(0)=T0},T(t));

sol1 := T(t) = Tu+Q/(k*A)-Q*exp(-k*A*t/(M*c)+t*(k*A-tau*M*c)/(M*c))/(k*A-tau*M*c)+exp(-k*A*t/(M*c))*(T0-Tu-Q/(k*A)+Q/(k*A-tau*M*c))

> T_off:=unapply(rhs(sol1),t):

次に、初期条件 T(t_off) = T_off(t_off) を使用して、ヒーターが消された後の状況を表す微分方程式を解き、Tの解を求めます。

> sol2:=dsolve({subs(q=q[off],deq),T(t_off)=T_off(t_off)},T(t));

sol2 := T(t) = -exp(-k*A*t/(M*c))*(-(k^2*A^2*Tu-k*A*Tu*tau*M*c+Q*k*A-Q*tau*M*c-Q*exp(-t_off*tau)*k*A+exp(-k*A*t_off/(M*c))*k^2*A^2*T0-exp(-k*A*t_off/(M*c))*k*A*T0*tau*M*c-exp(-k*A*t_off/(M*c))*k^2*A^2...sol2 := T(t) = -exp(-k*A*t/(M*c))*(-(k^2*A^2*Tu-k*A*Tu*tau*M*c+Q*k*A-Q*tau*M*c-Q*exp(-t_off*tau)*k*A+exp(-k*A*t_off/(M*c))*k^2*A^2*T0-exp(-k*A*t_off/(M*c))*k*A*T0*tau*M*c-exp(-k*A*t_off/(M*c))*k^2*A^2...sol2 := T(t) = -exp(-k*A*t/(M*c))*(-(k^2*A^2*Tu-k*A*Tu*tau*M*c+Q*k*A-Q*tau*M*c-Q*exp(-t_off*tau)*k*A+exp(-k*A*t_off/(M*c))*k^2*A^2*T0-exp(-k*A*t_off/(M*c))*k*A*T0*tau*M*c-exp(-k*A*t_off/(M*c))*k^2*A^2...sol2 := T(t) = -exp(-k*A*t/(M*c))*(-(k^2*A^2*Tu-k*A*Tu*tau*M*c+Q*k*A-Q*tau*M*c-Q*exp(-t_off*tau)*k*A+exp(-k*A*t_off/(M*c))*k^2*A^2*T0-exp(-k*A*t_off/(M*c))*k*A*T0*tau*M*c-exp(-k*A*t_off/(M*c))*k^2*A^2...sol2 := T(t) = -exp(-k*A*t/(M*c))*(-(k^2*A^2*Tu-k*A*Tu*tau*M*c+Q*k*A-Q*tau*M*c-Q*exp(-t_off*tau)*k*A+exp(-k*A*t_off/(M*c))*k^2*A^2*T0-exp(-k*A*t_off/(M*c))*k*A*T0*tau*M*c-exp(-k*A*t_off/(M*c))*k^2*A^2...

> T:=unapply(rhs(sol2),t,t_off):

最高温度の決定

最高温度を決定するために、時間変数に関して温度関数を微分した解にパラメタの値を挿入します。

> dT := unapply(subs({val,valQ},diff(T(t,t_off),t)),t_off);

dT := proc (t_off) options operator, arrow; 0.2391200383e-2*exp(-0.2391200383e-2*t)*(-102.4497795-2.449779520*exp(-1/10*t_off)+102.4497795*exp(-0.2391200383e-2*t_off))/exp(-0.2391200383e-2*t_off)+.244...dT := proc (t_off) options operator, arrow; 0.2391200383e-2*exp(-0.2391200383e-2*t)*(-102.4497795-2.449779520*exp(-1/10*t_off)+102.4497795*exp(-0.2391200383e-2*t_off))/exp(-0.2391200383e-2*t_off)+.244...

t_off には最高温度の解を使用します:

> MAX:=t_off->fsolve(dT(t_off)=0,t):

さらに使用しやすくするために、線型スプラインへ発展させます。

> l:=[25*i$i=2..8];

l := [50, 75, 100, 125, 150, 175, 200]

> ll:=map(MAX,l);

ll := [72.35139178, 93.50982067, 115.8588490, 138.8653067, 162.2869020, 185.9940729, 209.9094634]

> Max:=unapply(spline(l,ll,t),t,linear);

Max := proc (t, linear) options operator, arrow; piecewise(t < 75, 30.59174482+.835192939200000040*t+0.111022302462515656e-17*(t-50)^2+0.178307462500000002e-4*(t-50)^3, t < 100, 28.36290154+.868625588...Max := proc (t, linear) options operator, arrow; piecewise(t < 75, 30.59174482+.835192939200000040*t+0.111022302462515656e-17*(t-50)^2+0.178307462500000002e-4*(t-50)^3, t < 100, 28.36290154+.868625588...Max := proc (t, linear) options operator, arrow; piecewise(t < 75, 30.59174482+.835192939200000040*t+0.111022302462515656e-17*(t-50)^2+0.178307462500000002e-4*(t-50)^3, t < 100, 28.36290154+.868625588...Max := proc (t, linear) options operator, arrow; piecewise(t < 75, 30.59174482+.835192939200000040*t+0.111022302462515656e-17*(t-50)^2+0.178307462500000002e-4*(t-50)^3, t < 100, 28.36290154+.868625588...Max := proc (t, linear) options operator, arrow; piecewise(t < 75, 30.59174482+.835192939200000040*t+0.111022302462515656e-17*(t-50)^2+0.178307462500000002e-4*(t-50)^3, t < 100, 28.36290154+.868625588...Max := proc (t, linear) options operator, arrow; piecewise(t < 75, 30.59174482+.835192939200000040*t+0.111022302462515656e-17*(t-50)^2+0.178307462500000002e-4*(t-50)^3, t < 100, 28.36290154+.868625588...Max := proc (t, linear) options operator, arrow; piecewise(t < 75, 30.59174482+.835192939200000040*t+0.111022302462515656e-17*(t-50)^2+0.178307462500000002e-4*(t-50)^3, t < 100, 28.36290154+.868625588...

温度のグラフ表示

プロセスを可視化する為に、区分関数に値とモデルの液体の温度を入力します:

> Tp := unapply(subs({val, valQ},piecewise(t<t_off, T_off(t), T(t,t_off))), t,t_off);

Tp := proc (t, t_off) options operator, arrow; piecewise(t < t_off, 120.0000000+2.449779520*exp(-.1000000000*t)-102.4497795*exp(-0.2391200383e-2*t), -exp(-0.2391200383e-2*t)*(-102.4497795-2.449779520*...Tp := proc (t, t_off) options operator, arrow; piecewise(t < t_off, 120.0000000+2.449779520*exp(-.1000000000*t)-102.4497795*exp(-0.2391200383e-2*t), -exp(-0.2391200383e-2*t)*(-102.4497795-2.449779520*...Tp := proc (t, t_off) options operator, arrow; piecewise(t < t_off, 120.0000000+2.449779520*exp(-.1000000000*t)-102.4497795*exp(-0.2391200383e-2*t), -exp(-0.2391200383e-2*t)*(-102.4497795-2.449779520*...

> pl1:=plot3d(Tp(t,t_off), t=0..tmax, t_off=50..tmax, color=red):
pl2:=plot3d(['Max(t_off)', t_off, 'Tp(Max(t_off),t_off)'], t_off=50..tmax, dummy=0..1, thickness=2):
pl3:=plots[textplot3d]([Max(100), 100, Tp(Max(100),100),
      " Maximum Temp. (t_off)"], color=black, font=[TIMES, BOLD, 12]):
plots[display](pl1, pl2, pl3, axes=FRAME, orientation=[-101,62],
      title="Temp. of Liquid (red) & Curve of Maximum Temp. (black)",
      labels=["Time [s]","Turn off Time [s] ","Temp. [°C]"] );

[Plot]

ターンオフ時間の計算

40℃を超えないように加熱をオフにする必要があり、温度の関数に最大値を入力して、それを必要な温度と同一視しグラフからターンオフ時間の近似解を求めます。

> fsolve(Tp(Max(t_off),t_off)=40,t_off=100);

97.44016506

結論

シンボリックな計算は、液体を加熱した際の可視化や必要ターンオフ時間の計算を可能にします。

参考文献

Frank P. Incropera and David P. De Witt, Fundamentals of Heat and Mass Transfer, John Wiley & Sons.