熱効率を計算するプログラム

突然ですが問題

熱効率はいくつ?

あ、言い忘れてましたが、以降この記事では全て単原子分子理想気体とします

 

これはまあ基本問題です

この手の問題は表を作るのが良いでしょう

経路 dU Q W
A→B $\dfrac{3}{2} P_0V_0$ $\dfrac{3}{2} P_0V_0$ $0$
B→C $3 P_0V_0$ $5 P_0V_0$ $-2P_0V_0$
C→D $-3 P_0V_0$ $-3 P_0V_0$ $0$
D→A $-\dfrac{3}{2} P_0V_0$ $-\dfrac{5}{2} P_0V_0$ $-P_0V_0$

 

Qの列を見て、符号がプラスのものとマイナスのものに分ける

A→BとB→Cがプラスで、C→DとD→Aがマイナスである。

で、符号がプラスのものがエネルギーの入力、符号がマイナスのものがエネルギーの出力となる。

で、入力は$\dfrac{13}{2}P_0V_0$で、出力は$\dfrac{11}{2}P_0V_0$です。入力されたエネルギーにおいて、一部分は仕事をしますが残りは出力エネルギーとなる。

エンジンをイメージしてほしい。このようなエンジンがあったとき、最初に$\dfrac{13}{2}P_0V_0$の分だけのエネルギーが投入さるが、そのうちエンジンがする仕事は$\dfrac{13}{2}P_0V_0-\dfrac{11}{2}P_0V_0=P_0V_0$であり、残りの$\dfrac{11}{2}P_0V_0$はただエンジンが熱くなったりするだけで仕事をせずに無駄になるエネルギーになると思っておこう。

すると、熱効率は

$$1-\dfrac{\dfrac{11}{2}P_0V_0}{\dfrac{13}{2}P_0V_0}=\dfrac{2}{13}$$

となる。

 

次の問題にしようと思ってたものがむずいのでもうちょっと良心的なものを置く

断熱変化では$PV^\gamma,\gamma=\dfrac{5}{3}$が一定なので、$8^{5/3}=32$よりポアソンの式は満たしている。

これも表を書く

経路 dU Q W
A→B $\dfrac{93}{2} P_0V_0$ $\dfrac{93}{2} P_0V_0$ $0$
B→C $-36P_0V_0$ $0$ $-36P_0V_0$
C→A $-\dfrac{21}{2}P_0V_0$ $-\dfrac{35}{2} P_0V_0$ $7P_0V_0$

入力は$\dfrac{93}{2}P_0V_0$で、出力は$\dfrac{35}{2}P_0V_0$なので、

熱効率は

$$1-\dfrac{\dfrac{35}{2}P_0V_0}{\dfrac{93}{2}P_0V_0}=\dfrac{58}{93}$$

 

 

ということで次

C→Aは、PV図で斜めの直線となっている。つまりこれは定積変化でも定圧変化でも等温変化でも断熱変化でもない。

仕事の部分は、PV図の面積を求めればOK(つまり積分する)、内部エネルギーの変化は温度の変化だけ見ればOKなので、Qはその差を求めればOKなのでこう書けば良さそう???

経路 dU Q W
A→B $-\dfrac{3}{2} P_0V_0$ $-\dfrac{5}{2} P_0V_0$ $P_0V_0$
B→C $\dfrac{3}{2}P_0V_0$ $\dfrac{3}{2} P_0V_0$ $P_0V_0$
C→A $0$ $\dfrac{3}{2} P_0V_0$ $-\dfrac{3}{2}P_0V_0$

なので熱効率は

$$1-\dfrac{\dfrac{5}{2}P_0V_0}{\dfrac{3}{2}P_0V_0+\dfrac{3}{2}P_0V_0}=\dfrac{1}{6}$$

…としたいところだけど、残念ながらこれは間違っている

結論から言うと、正しい表はこれだ!

経路 dU Q W
A→B $-\dfrac{3}{2} P_0V_0$ $-\dfrac{5}{2} P_0V_0$ $P_0V_0$
B→C $\dfrac{3}{2} P_0V_0$ $\dfrac{3}{2} P_0V_0$ $0$
C→D $\dfrac{21}{128} P_0V_0$ $\dfrac{49}{32} P_0V_0$ $-\dfrac{175}{128}P_0V_0$
D→A $-\dfrac{21}{128} P_0V_0$ $-\dfrac{1}{32} P_0V_0$ $-\dfrac{17}{128}P_0V_0$

最初の図も訂正しておこう。問題文の図もこうした方がいい

ということで正しい答えは

$$1-\dfrac{\dfrac{5}{2}P_0V_0+\dfrac{1}{32}P_0V_0}{\dfrac{3}{2}P_0V_0+\dfrac{49}{32}P_0V_0}=\dfrac{16}{97}$$

そう、新しい点Dを置かないとだめなのである!

なぜこうなるかというと、C→Aの経路をたどる途中で、最初はエネルギーが注入されるけどそれはDまでで、Dを通り過ぎた瞬間エネルギーが注入から放出に切り替わるのである!つまり、熱効率を変化するためにはδQの符号を追跡しなくてはいけないけど、この場合はδQの符号が途中で変わることを考慮しなくてはいけないのである!

δQの符号が途中で変わること、は定圧変化でも定積変化でも等温変化でも断熱変化でも起こり得ないので、このような特殊な変な経路をたどるときに初めて考慮する必要が出てくるのである。

 

で、気になるのがどうして途中で切り替わるのか、切り替わるとしてなぜ$\dfrac{9}{8}P_0,\dfrac{15}{8}V_0$だとわかるのかという問題である。

 

一般論として、このような熱機関において、各パラメータは変数で置くことができる。

内部エネルギー$U$は$(P,V)$の関数で以下のようにかける

$$U=\dfrac{3}{2}PV$$

なので、$U$の微小変化は以下のように書ける

$$dU=\dfrac{3}{2}\{(P+dP)(V+dV)-PV\}$$

$$=\dfrac{3}{2}(PdV+VdP)$$

(微小量×微小量は無視できるので$dPdV=0$である)

で、

仕事$\delta W$は以下のように書ける

$$\delta W=-PdV$$

すると、熱力学第一法則より、

$$\delta Q=\dfrac{3}{2}VdP+\dfrac{5}{2}PdV$$

となる。

なぜ$U$の微小変化は$dU$なのに$Q,W$の微小変化は$\delta Q,\delta W$と書いているかというと、$\dfrac{3}{2}(PdV+VdP)=f(P+dP,V+dV)-f(P,V)$となるような関数$f$は存在する$f(P,V)=\dfrac{3}{2}PV$とおけばいい

けど、$-PdV=g(P+dP,V+dV)-g(P,V)$なる関数$g$や$\dfrac{3}{2}VdP+\dfrac{5}{2}PdV=h(P+dP,V+dV)-h(P,V)$となる関数$h$は存在しないのである!

これは、$U$は状態量であるが、$Q$や$W$は状態量じゃないからこうなってる

状態量とはどういうことなのか?という人のために説明

上の表とかを書くときに、$U$を計算するときは始点と終点の情報だけわかってれば求めることができる。正直言って、始点と終点さえ同じだったら、途中の経路がどのようなものであっても内部エネルギーの変化は同じなのである。

しかし$Q$や$W$はそうではない。$Q$や$W$を計算するときは、最終的な答えが途中の経路に依存するのである。なので、始点と終点の情報だけでは計算できない。途中の経路がどのようなものだったのかという情報が必要である。

で、最終的にQやWの計算については$\delta Q$や$\delta W$を線積分すれば求めることができる。

 

で、線積分だけど、$(P,V)=(f(t),g(t))$というような媒介変数で書ける経路において$a(P,V)dP+b(P,V)dV$を線積分したいなら

$$\oint a(P,V)dP+b(P,V)dV=\int a(f(t),g(t))f'(t)dt+\int b(f(t),g(t))g'(t)dt$$

というようになる。

 

本題に戻る。このようなC→Aの直線経路において、

$$(P,V)=((2-t)P_0,(1+t)V_0)$$

というような媒介変数で書くことができる。ただし$t$の動く範囲は$0\leq t\leq 1$

$$\delta Q=\dfrac{3}{2}VdP+\dfrac{5}{2}PdV$$

$$=\dfrac{3}{2}(1+t)V_0 \cdot (-P_0)dt+\dfrac{5}{2}(2-t)P_0 \cdot V_0 dt$$

$$=\left(\dfrac{7}{2}-4t\right)P_0V_0dt$$

となっている。

で、tごとに$\delta Q$の微小変化を見ると、

$0\leq t\leq\frac{7}{8}$のときに$\delta Q\geq 0$

$\frac{7}{8}\leq t\leq 1$のときに$\delta Q\leq 0$

となる。なので、$t=7/8$となる地点、つまり$(P,V)=(\dfrac{9}{8}P_0,\dfrac{15}{8}V_0)$で一旦分けないといけないのである。

 

これを踏まえた上で次の問題

いや、さっきの問題でやることはわかったけどただクソめんどいだけじゃん

 

ということなので

 

退屈なことはPythonにやらせよう

 

ということでソースコード、ペタッ!w

gistb858bb113423b471eeaa29a8bd22622c

 

説明

PV図上の点において、$N$個の点からなるサイクルを考える、各点の間の経路は全部直線であると仮定した場合、熱効率はいくらになるか?

 

この記事を読んだ人々が真似しやすいように仕組みも書きます。

 

まずは直線経路上において受け取った熱の量と放出した熱の量を計算する方法から

$a,b,c,d$を正の定数とします。(変数のdと微小量のdでダブってるので以降微小量は$\Delta$で書きます)で、$(aP_0,bV_0)\to (cP_0,dV_0)$に向かう直線の経路を考えたとき、

$$P=((c-a)t+a)P_0,V=((d-b)t+b)V_0$$

となります。

で、このとき、

$$\delta Q=\dfrac{3}{2}V\Delta P+\dfrac{5}{2}P\Delta V$$

$$=4(c-a)(d-b)t+\dfrac{3}{2}b(c-a)+\dfrac{5}{2}a(d-b):=At+B$$

とします。

$0\leq t\leq 1$において、$At+B\geq 0$ならば熱を受け取っていて、$At+B\leq 0$ならば熱を放出しています。なので、$0\leq t\leq 1$の範囲内で$At+B$の符号が変化するかどうかを確かめます。符号が変化するなら符号の変化点で区切ります。符号が変化しないならば、

$$U=\dfrac{3}{2}(cd-ab)P_0V_0$$

と、

$$\Delta W=(d-b)((c-a)t+a)P_0V_0dt$$

で、これをtについて0から1まで積分すればOK

こうすることで、直線経路上において、(受け取った熱の量、放出した熱の量)の組を計算することができます

 

次に、熱効率について

これは全ての線分上において、受け取った熱の量の合計$O_{in}$と放出した熱の量の合計$Q_{out}$を合計して、それについて$$1-\dfrac{Q_{out}}{Q_{in}}$$とすればOK

 

で、このプログラムは、$Q_{out}\gt Q_{in}$のときは熱効率が計算しないようにしている。これを判定するために仕事を計算している。$\oint_C PdV\lt 0$なら熱効率は計算できない。これは$PV$図において時計回りの経路なら熱効率が計算できるけど、反時計回りの経路ならできないということである。

 

で、質問の答え(八角形の経路)としては、プログラムに以下の入力を打ち込むだけ

 

print(thermal_efficiency([(2,1),(3,1),(4,2),(4,3),(3,4),(2,4),(1,3),(1,2)] ) )

 

 

というわけで答えは32/119

 

よくある質問

なんで断熱変化や等温変化に対応してないのですか?

答えが有理数にならないから大変そうだと思った。

とりあえず等温変化について

等温変化だと$P=tP_0,V=\dfrac{V_0}{t}$と書くことができる。$\alpha \leq t\leq \beta,0\lt \alpha \lt \beta$

$$dP=P_0dt,V_0=-\dfrac{V_0}{t^2}dt$$

なので、

$$\dfrac{3}{2}VdP+\dfrac{5}{2}PdV=\dfrac{3}{2}\cdot \dfrac{1}{t}-\dfrac{5}{2}\cdot \dfrac{1}{t}P_0V_0dt=-\dfrac{1}{t}P_0V_0dt$$

なので、これを積分すると自然対数が出てくる。

断熱変化についても、$\delta Q=0$なので対処しやすそうに見えるけど、$PV^{\gamma}=const$なので始点終点によってはパラメータに三乗根が出てくる。

 

これらに対してどうにかするには厳密な値を諦めてfloatで近似値を出す実装にするか、伝家の宝刀sympyを使うかくらいしかなさそう。そしてプログラムのインプットが経路の頂点だけになってるのまずそう。各経路がどのような変化かを明言する必要がありそう。とりあえず直線変化(∋定圧変化、定積変化)、等温変化、断熱変化だけに対応すれば良さそう?実装めんどいので妥協

 

δA=f(P,V)dP+g(P,V)dVみたいなのが状態量であるかどうかってどうやって判定すればいいですか?

これは外微分をすればOK

というわけでこの「微小量みたいな形になっているもの」さらに形式的に微分するみたいなことをする

$$d\delta A=\{\dfrac{\partial f}{\partial P} (P,V)dP+\dfrac{\partial f}{\partial V} (P,V)dV\}\land dP+\{\dfrac{\partial g}{\partial P} (P,V)dP+\dfrac{\partial g}{\partial V} (P,V)dV\}\land dV$$

$$=\dfrac{\partial f}{\partial P}dP\land dP+\dfrac{\partial f}{\partial V}dV\land dP+\dfrac{\partial g}{\partial P}dP\land dV+\dfrac{\partial g}{\partial V}dV\land dV$$

で、微小量と微小量に対して$dx\land dy$になってるのは何やねんとなるのだが、

とりあえず以下のような計算法則が成り立つものだと思っておけばOK

$$dx\land dy=-dy\land dx$$

$$dx\land (bdy+cdz)=b dx\land dy+c dx\land dz$$

$$(a dx+b dy)\land dz=a dx\land dz+b dy\land dz$$

 

ということでかけ算の順序を入れ替えると符号が逆になることに注意。「ベクトルの外積みたいなもの」だと思えば理解できる人がいるかも

というわけで計算練習をすると、先程の計算規則に対して$x=y$とおくと

$$dx\land dx=-dx\land dx$$

なので$$dx\land dx=0$$

となる。

 

詳しく知りたい人へ→

hooktail.sub.jp

 

ということで、

$$d\delta Q=\{\dfrac{\partial f}{\partial V}-\dfrac{\partial g}{\partial P}\}dV\land dP$$

である。これが0にならなければ、確実に存在しない。

つまり、

$$Q=h(P,V)$$ならば

$$\delta Q=\dfrac{\partial h}{\partial P}dP+\dfrac{\partial h}{\partial V}dV$$

なので、

$$d \delta Q=\{\dfrac{\partial^2 h}{\partial P\partial V}-\dfrac{\partial^2 h}{\partial V\partial P} \}dV\land dP$$

であり、お行儀が良い関数ならば偏微分する順番を入れ替えても同じ結果になるのでこれが0になるのである。

つまり、もう一回微分して0になることが状態量が存在することの必要条件である。

一方、もしこれが0になったら、これは状態量である。$f$を$P$で積分したり$g$を$V$で積分したりすると候補が見つかるかも

ja.wikipedia.org

 

 

δQは状態量にはならないけどδQ/Tは状態量になるでは?と思いましたがどうなんですか

おっいい質問ですね

$$\delta Q=\dfrac{3}{2}VdP+\dfrac{5}{2}PdV$$

$$\dfrac{\delta Q}{T}=\dfrac{3}{2}\dfrac{V}{T}dP+\dfrac{5}{2}\dfrac{P}{T}dV$$

で、$PV=nRT$より、

$$\dfrac{\delta Q}{T}=\dfrac{3}{2} nR\dfrac{dP}{P}+\dfrac{5}{2}nR \dfrac{dV}{V}$$

で、

$$d\dfrac{\delta Q}{T}=0$$

なので状態量である。

で、

$$S=\dfrac{3}{2}nR \log\{PV^{\gamma}\}$$

という量を置くと($\gamma$は比熱比で単原子分子理想気体ならば$\gamma=5/3$)

$$dS=\dfrac{\delta Q}{T}$$

となります。

このSには名前がついていて、エントロピーというらしいですね

$\delta Q=0$なら$dS=0$なので、断熱変化のときに変化しない量となります。

 

知らんけど (自分はこのような理解をしているのですが、間違ってたらこの記事を燃やしてください)