２自由度集中定数系熱容量の温度

質問

集中定数系熱容量の過渡温度の問題が解けず困っています。
熱容量C１[J/K]、C2[J/K]の温度をそれぞれθ1(t)[℃]、θ2(T)[℃]とする。 θ1は熱通過率α1[W/K]を介して高温側温度T1[℃]とつながっている。θ1(t)とθ2(t)は熱通過率α2[W/K]を介してつながっている。θ2は熱通過率α3を介して低温側温度T2[℃]とつながっている。θ1(t)、θ2(t)の初期温度をそれぞれθ10（℃）、θ20とする。時間ｔ[ｓｅｃ]に対するθ1(t)、θ2(t)を求めよ。
α1、α2、α3の熱通過量をＱ1(t)、Ｑ2(t)、Ｑ3(t)[Ｗ]とし、C1、C2に対する微分方程式をたて、ラプラス変換すれば解けると思うのですが、C1の微分方程式にθ2(t)が入り、C2の微分方程式にθ1(t)が入り、変数分離の考え方が分かりません。
そのへんも含めて詳しくお願いします。

回答No.3

FT56F001 · ベストアンサー

> どこが間違っているのでしょうか？
> C1θ1'=Q1-Q2=α1(T1-θ1)-α2(θ1-θ2) (4)
> C2θ2'=Q2-Q3=α2(θ1-θ2)-α3(θ2-T2) (5)
質問者さんは，解くべき微分方程式を微分して新しい微分方程式(6)を作ったところで，
無縁解を引き入れてしまったのではないかと思います。
さらに部分分数展開する時の因子は(s+β1)ではなくて(s-β1)のはずです。
この二つにより，発散する解を引き込んでしまったと思います。

もともとが2変数の連立微分方程式ですから，
ラプラス変換を使うなら，(4)(5)をまず変換してしまいます。
C1*(θ1*s-θ10)=α1*T1/s -(α1+α2)*θ1 + α2*θ2
C2*(θ2*s-θ20)=α3*T2/s -(α2+α3)*θ2 + α2*θ1

これをθ1とθ2に関する連立方程式として解きます。
θ1(s)=分子1/分母
θ2(s)=分子2/分母

分子1=C1*C2*θ10*s^2 + (C2*α1*T1+C1*(α2+α3)*θ10+C2*α2*θ20)s
+ α1(α2+α3)T1+α2*α3*T2
分子2=C1*C2*θ20*s^2 + (C1*α3*T2+C2*(α1+α2)*θ20+C1*α2*θ10)s
+ α1*α2*T1+α3(α1+α2)T2
分母=s*{C1C2s^2 + {C1(α2+α3)+C2(α1+α2)}s+(α1α2+α2α3+α3α1)}

分母にあたる特性方程式
C1C2 s^2 + {C1(α2+α3)+C2(α1+α2)}s+(α1α2+α2α3+α3α1)=0
の2根をβ1，β2とします。
β1，β2は負の実根となるはずです。(でないと発散してしまう)

θ1(s)，θ2(s)を部分分数展開します。
θ1=分子1/{C1C2*s(s-β1)(s-β2)}
=M1/s + M2/(s-β1) + M3/(s-β2)
θ2=分子2/{C1C2*s(s-β1)(s-β2)}
=M4/s + M5/(s-β1) + M6/(s-β2)

M1=(α1(α2+α3)T1+α2*α3*T2)/(α1α2+α2α3+α3α1)
M4=(α1*α2*T1+α3(α1+α2)*T2)/(α1α2+α2α3+α3α1)
M1とM4は十分時間がたった後の定常温度を表します。

M2={α2*Δθ2-(C1*β2+α1+α2)Δθ1}/{C1(β1-β2)}
M3={-α2*Δθ2+(C1*β1+α1+α2)Δθ1}/{C1(β1-β2)}
M5={α2*Δθ1-(C2*β2+α2+α3)Δθ2}/{C2(β1-β2)}
M6={-α2*Δθ1+(C2*β1+α2+α3)Δθ2}/{C2(β1-β2)}
ただし，Δθ1 = θ10-M1，Δθ2 = θ20-M4
これらは初期温度から定常温度までの変化分を表します。

部分分数展開されたθ1(s)，θ2(s)を逆ラプラス変換して，
θ1(t)=M1 + M2*exp(β1*t) + M3*exp(β2*t)
θ2(t)=M4 + M5*exp(β1*t) + M6*exp(β2*t)
が答えです。

計算違ってるかもしれませんので，確認下さい。

お礼

丁寧な回答を頂き、ありがとうございました。
回答を解き、アンサーに入力して頂くだけでも相当の時間を要したと思います。
連立微分方程式を解いたことがなく、独学で勉強していましたので、非常に助かりました。
この２週間悶々としていたのですが、すっきりしそうです。
最後に
β1とβ2は負の実根ですから部分分数に分解したときの分母はｓ+β1とs+β2になると思います。
なぜなら逆ラプラス変換したときに
1/(s＋β）はexp(-β*t)に、1/(s-β)はexp(+β*t)とβの符号が逆転するからです。
（exp(-β*t)は収束しますが、exp(+β*t)は発散します。）
ともかく、回答の正しいプロセスをご教授頂いたので、丁寧に解き直してみます。
またわからない事があったら教えてください。
どうもありがとうございました。

回答No.2

FT56F001

>各熱容量に対する入力熱量と出力熱量の差が、熱容量*温度変化、となるので、
>C1 dθ1/dt = α1(T1-θ1) - α2(θ2-θ1)
>C2 dθ2/dt = α2(θ1-θ2) - α3(T2-θ2)
>となるの思うのですが、間違いでしょうか？
とはならないです。
通常の熱の流れだと，T1>θ1>θ2>T2の温度分布になりますね。
θ1-θ2(>0)の温度差の分に比例してC1から熱が出て行きます。
また，C2から出て行く熱もθ2-T2(>0)に比例します。
C1 dθ1/dt = α1(T1-θ1) - α2(θ1-θ2)
C2 dθ2/dt = α2(θ1-θ2) - α3(θ2-T2)
が正しいです。

あと訂正です。先回の回答で「連立同次線形微分方程式」と書きましたが，
定数項α1T1，α3T2があるので，同次ではありませんでした。
α1(T1-θ1) - α2(θ1-θ2)=0
α2(θ1-θ2) - α3(θ2-T2)=0
の解θ1，θ2が，定常状態(t->∞のときの温度)を与えます。
定常状態の温度分だけ，θ1，θ2をずらせると，定数項が消えて，
同次になります。

お礼

ご教授、ありがとうございます。
連立微分方程式をたて、ラプラス変換で解いてみました。
が、エクセルで時間tに対し温度θ1、θ2　と　α1、α2、α3を通過する熱量Q1、Q2、Q3のグラフを描いてみたところ、定性的挙動が物理現象と一致せず、うまく解けていないようです。
下記に解法プロセスと結果を示します。
どこが間違っているのでしょうか？

記述を簡単にするため、θ1(t),θ2(t)　は　θ1,θ2と、　Q1(t),Q2(t),Q3(t)　はQ1,Q2,Q3と、
dθ/dt,d2θ/dt2　は　θ',θ''と記述します。

Q1=α1(T1-θ1) (1)
Q2=α2(θ1-θ2) (2)
Q3＝α3(θ2-T2)　　　(3)
C1θ1'=Q1-Q2=α1(T1-θ1)-α2(θ1-θ2) (4)
C2θ2'=Q2-Q3=α2(θ1-θ2)-α3(θ2-T2) (5)

まずθ1について求める

(4)式より
θ2=(α1+α2)/α2*θ1-α1/α2*T1+C1/α2*θ1' (6)
(6)式を微分して
θ2'=(α1+α2)/α2*θ1'+C1/α2θ1'' (7)

(6)を(5)に代入
C2θ2'=α2θ1-(α2+α3){(α1+α2)/α2θ1+C1/α2θ1'-α1/α2T1}+α3T2 (8)
(7)を(8)に代入し、θ1'',θ1',θ1について整理すると
C1C2/α1*θ1''+{(α1+α2)C2/α2+(α2+α3)C1/α2}θ1'+{(α1+α2)(α2+α3)/α2-α2}θ1-(α2+α3)α1/α2*T1-α3T2=0 (9)

簡単のため各項の係数を
A=C1C2/α1 (10)
B=(α1+α2)C2/α2+(α2+α3)C1/α2 (11)
C=(α1+α2)(α2+α3)/α2-α2 (12)
D=(α2+α3)α1/α2*T1+α3T2 (13)
とおくと(9)式は

Aθ1''+Bθ2'+Cθ1-D=0 (14)
θ1''+B/Aθ2'+C/Aθ1-D/A=0 (15)

(15)式をラプラス変換する
θ10、θ10’は初期値

s2Θ1 + B/AsΘ1 + C/AΘ1　- D/A*1/s - sθ10 - θ10’ - B/Aθ10 =0　　(16)
(s2+B/A*s+C/A)Θ1=D/A*1/s + θ10(s-B/A)+θ10' (17)
θ10'は時間t<=0においてθ10=constでθ10'=0ととし、両辺をs2+B/A*S+C/Aで割ると

Θ1=D/A*1/s*1/(s2+B/A*s+C/A) + θ10(s-B/A)*1/(s2+B/A*s+C/A) (18)

θ1は安定系であるから特性方程式は異なる2実解を持つ
解をβ1,β2(β2>β1)とする
β1,β2は解の公式より
β1={-B-√(B^B-4AC)}/2A (19)
β2={-B+√(B^B-4AC)}/2A (20)

(18)式を部分分数に分解する
右辺係数D/A項について、分子をK0,K1,K3とし
K0/s + K1/(s+β1) + K2/(S+β2) = 1/{s(s+β1)(s+β2)} (21)
とおき、これを求めると
K0=0, K1=-β2/(β2-β1), K2=β1/(β2-β1)　　(22)

θ10の項について、分子をK3,K4とし
K3/(s+β1) + K4/(s+β2) = (s+B/A)/{(s+β1)(s+β2)}　　　(23)
とおき、これを求めると
K3=(B/A-β1)/(β2-β1), K4=-(B/A-β2)/(β2-β1),

よって
Θ1=D/A{1/s - β2/(β2-β1)*1/(s+β1) + β1/(β2-β1)*1/(s+β2)}
+θ10/(β2-β1){(B/A-β1)*1/(s+β1)-(B/A-β2)*1/(s+β2)} (24)

Θ1を逆ラプラス変換して
θ1=L-1(Θ1)=D/A｛1-1/(β2-β1){β2exp(-β1*t)-β1exp(β2*t)}｝
+θ10/(β2-β1)｛(B/A-β1)exp(-β1*t)-(B/A-β2)exp(-β2*t｝ (25)

θ2もθ1と同様にして特性方程式の解をβ1',β2'として解くと
θ2=D'/A'｛1-1/(β2'-β1'){β2exp(-β1'*t)-β1'exp(β2'*t)}｝
+θ10/(β2'-β1')｛(B'/A'-β1')exp(-β1'*t)-(B/A-β2')exp(-β2'*t｝(26)

但し
A'=C1C2/α1 (27)
B'=(α1+α2)C2/α2+(α2+α3)C1/α2 (28)
C'=(α1+α2)(α2+α3)/α2-α2 (29)
D'=(α1+α2)α3/α2T2+α1T1 (30)

β1',β2'は解の公式より
β1'={-B'-√(B'^B'-4A'C')}/2A' (31)
β2'={-B'+√(B'^B'-4A'C')}/2A' (32)

と計算しました。
何度計算しても係数等に間違いはないようで、
解法が間違っているところがあるようです。
どこがおかしいのでしょうか？

回答No.1

FT56F001

電気回路に置き換えれば，
一定電圧源T1,T2の間に，
コンダクタンスα1，α2，α3が横につながって，
コンデンサC1,C2が縦にぶら下がった等価回路で表せます。
(コンデンサの電圧がθ1，θ2)
/* 電気回路に置き換えると
見やすくなるか，かえって混乱するか，
は人によりますが */

微分方程式は，
C1 dθ1/dt = α1(T1-θ1) + α2(θ2-θ1)
C2 dθ2/dt = α2(θ1-θ2) + α3(T2-θ2)
となります。

定係数の連立同次線形微分方程式ですから，
係数行列の固有値を求めて，指数関数で表す，
という標準的な解き方があります。

例えば
http://next1.cc.it-hiroshima.ac.jp/MULTIMEDIA/diffpub/node33.html

お礼

御回答頂まして、有難うございます。
行列は苦手なので、勉強しています。
微分方程式ですが、
各熱容量に対する入力熱量と出力熱量の差が、熱容量*温度変化、となるので、
C1 dθ1/dt = α1(T1-θ1) - α2(θ2-θ1)
C2 dθ2/dt = α2(θ1-θ2) - α3(T2-θ2)
となるの思うのですが、間違いでしょうか？