Moneytizerによる広告を掲載しています。本記事ではPythonでの化学工学計算例を紹介しています。今回はNewton-Raphson法を使用して、対数平均温度差を含む熱収支を計算します。熱交換器の設計をする際に計算する式で、実務でもよく登場します。化学工学シミュレータを使用すればあっという間に解けるような式ですが、Pythonでいちから式を立てると意外に解けなくて苦戦することがあるかもしれません。本記事ではプログラミング系の言語を今まで触ったことがない人をターゲットとし、自分が躓いたところを中心になるべく基礎的なところから解説しています。逆に化学工学の計算式についてはある程度理解のあるものとし、解説は適宜割愛します。する式は以下の記事で解説したものを使用します。https://chemical-engineering-review.com/newton-raphson/上の記事ではExcelを使用してNewton-Raphson法で解いていますが、今回は同じ式をPythonで解いてみます。低温側流体流量F:45,000kg/h低温側流体比熱Cp:1.0kcal/(kg・°C)熱交換器の総括伝熱係数U:550kcal/(m2・h・°C)熱交換器の伝熱面積A:20m2高温側流体の入口温度Th1:78°C低温側流体の入口温度Tc1:31°C低温側流体の出口温度Tc2:38°C$$f(Th2)=FCp(Tc2-Tc1)-UA\frac{(Th1-Tc2)-(Th2-Tc1)}{{\rm{ln}}(\frac{Th1-Tc2}{Th2-Tc1})}・・・(1)$$前提条件と計算式は上に示す通りです。(1)式がf(Th2)=0となるようなTh2をNewton-Raphson法で求めます。計算できていれば、Th2=50.66°Cになるはずです。新しい計算ソフトに手を出すときは、本記事のようにとりあえず答えのわかる問題から解くのが良いと思います。Python記述ライブラリ・モジュールのインポートまず最初に使用するライブラリやモジュールをインポートします。今回はmathモジュールを使用します。対数平均温度差の項を計算するためにlogの値を計算する必要があるからです。Excelではlogとかlnと打てば対数計算できますが、Pythonではモジュールを最初にインポートしないと計算できません。次に問題で与えられる前提条件を入力します。°C念のためコメントで単位を付けておくと、後から見たときに単位変換ミスに気づきやすくてよいです。次に初期条件が必要な場合はその値を入力します。今回はNewton-Raphson法を使用するためTh2の初期値が必要ですね。°CTh20=40°Cと設定します。今回はざっくりとその他計算条件でまとめました。細かく分類できる場合は分けた方がよいかもしれません。f(Th2)の微分を求める際に温度刻みΔTが必要なのでその値を入力します。次にメインの方程式を定義します。定義の記述の仕方も様々あるようですが、今回はdefとreturnで関数を定義します。#(2)式定義2,3行目で(1)式を定義しています。対数については、math.log()と記述することでeが底の対数を計算できます。defでfTh2(Th2)という関数を定義し、3行目のreturn文で式を入力します。することで以降の記述でfTh2(Th2)と記載すると、return文に入力した式を計算して値を返してくれます。fTh2(Th2)の括弧内のTh2は数学で言うところの変数(プログラミング用語では引数と呼ばれます)で、様々な値を入力することができます。今回のようにTh2の初期値が40°Cから始まり、繰り返し計算でTh2の値が変わっていくような計算ではよく使用します。5,6行目でf(Th2)の微分f'(Th2)を定義しています。(1)式が簡単に微分可能であれば微分すればよいですが、今回は難しそうですね。その場合には数値微分の形でf(Th2)を算出します。$$f'(Th2(n+1))=\frac{f(Th2(n)+ΔT)-f(Th2(n))}{ΔT}・・・(2)$$こちらの記事で解説していますが、(2)式が数値微分した式になります。前進Euler法を使用しています。5,6行目に(2)式をdefで定義してメインの方程式定義は終了です。続いてNewton-Raphson法の収束計算式を入力していきます。今回はdefで収束計算式を定義し、while文で繰り返し計算を行ないます。#次のステップ2行目でnewton(Th20)という関数を定義し、4行目から収束計算式の条件を書いていきます。4行目は収束計算回数をnという文字で定義します。初期値は1回目なのでn=1とします。5行目はエラーの初期値を入力します。Th2(n+1)-Th2(n)\right|・・・(3)$$(3)式のようにn回目とn+1回目の計算値の絶対値を取ります。このerror値が繰り返し計算でどんどん小さくなり、ほとんどゼロに近くなれば収束したとみなすことができます。最初の1回目は適当な値を初期値として設定する必要があり、今回はerror=1としています。このerror値の初期値が小さすぎると最初から収束していると判定されかねないので、6行目の収束判定値と比較して十分大きい値を設定しましょう。6行目に収束判定値を入力します。化学工学の分野では10-4くらいの精度があれば十分でしょう。8行目で初期値の受け渡しを行なっています。Th2nをn回収束計算させた後のTh2の値とし、その初期値としてTh20=40°Cとしています。10行目で繰り返し計算のメインとなるwhile文を入力しています。while文は後に入力する条件式がFalseになるまで:以降の実行文を繰り返す処理をします。これはNewton-Raphson法の原理の通り、繰り返し計算のerror値が収束判定値exより大きければ繰り返し計算を続行することを意味しています。なる、つまりerror値が収束判定値ex以下になり収束したと判定できれば計算を終了します。似たような計算はfor文でも行なうことができますが、for文では繰り返し計算回数をあらかじめ指定する必要があります。今回の収束計算のように何回繰り返し計算するかわからない場合は、for文よりwhile文の方がコードを書きやすく向いていると思います。11~17行目にwhile文の実行文を入力しています。実行文であることを示すために、必ずTabキーで字下げして入力位置を揃えるようにしましょう。実行文の記述なのに字下げしていなかったり、文頭の位置が揃っていないとエラーになります。Pythonは特に字下げには厳しいです。11行目にNewton-Raphson法のメインの収束計算式である(4)式を入力します。$$x_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})}・・・(4)$$fTh2やf'Th2の()内をTh2nと入力することで、定義した関数のTh2にTh2nの値を入力して計算させることができます。12行目でn+1回目の計算値Th2n1とn回目の計算値Th2nの絶対値を取り、error値を算出しています。14行目でprint関数を使用して計算回数n、計算値Th2、errorの値を出力しています。この14行目はプログラミング言語特有の表現で、defで定義した関数の計算値を出力したい場合にはreturnやprintでわざわざ出力させる必要があります。returnとprintの使い分けに関してはプログラミング言語特有の関数の考えを理解する必要があり、自分もなかなか理解が追い付いていないです。両方試してみてうまくいったprintの方を採用しています。この箇所をreturnにすると収束計算1回目の値しか出力されずうまくいきませんでした。16行目で、11行目で計算したTh2n1の値を次のTh2nに受け渡す処理を行なっています。これもプログラミング言語特有の表現で、今回のような繰り返し計算で変数を計算値によって更新し新たな計算をする場合にはよく使用される表現です。今回は(4)式で計算させたTh2n1の値を新たなTh2nの値とし、もう一度(4)式に新たなTh2nの値を代入して計算する処理を行なっています。17行目も同様に繰り返し計算回数の更新を行なっています。この式はnに1を足したものをnに代入して、nの値を更新させることを意味しています。17行目の計算後はwhile文の冒頭に戻り、条件式がFalseになるまで同じ計算を繰り返します。最後に計算した結果を出力します。#関数newtonを出力このように記述することで、関数内で出力している計算回数n、計算値Th2、error値が出力されます。入力ミス等がなければ計算が収束し結果が以下のように表示されます。8.045155539093685e-06左から計算回数n、Th2、errorの順番で出力されていますが少し見づらいですね。そんな場合には出力方法を少し工夫しましょう。#計算回数、Th2、error値の出力while文のprintの入力で、各出力変数の間に"\t"を挟みます。これは水平タブと呼ばれ、水平方向の出力間隔を揃えるのに便利です。8.045155539093685e-06上のような表示になります。ブログへの記入の都合上、でこぼこしているところがありますが、実際の出力画面はちゃんと文字同士の間隔が揃って見やすくなるはずです。Excelで計算したこちらの記事の結果と同様に、4回で計算が収束しTh2=50.66°Cとなりました。下に計算コードを最初から最後までまとめて載せておきます。#関数newtonを出力Scipyモジュールを使用した別解PythonのモジュールにはNewton-Raphson法を簡単に計算できるものがあるようです。今回はScipy(サイパイ)というモジュールを使用してみます。全体の記述を上に示します。2行目でScipyのnewton法をインポートしています。前提条件、その他計算条件の記述は変わりません。大きく変わる点はfTh2の微分であるf'Th2を定義しなくてよいwhile文等の収束計算式を入力しなくてよい以上の点です。自分で全て入力したケースと比較するとだいぶ楽ですね。21行目でscipyのnewton計算式を入力しています。括弧内は計算する関数とTh2の初期値をそれぞれ入力します。この記述でコードをRunさせると、50.66047825978493このように収束したTh2の値が表示されます。慣れてくればこのようなモジュールを使用してどんどん楽に計算することができます。まとめPythonで対数平均温度差を含む熱収支計算をNewton-Raphson法で計算する方法を解説しました。非線形方程式の解法としてよく使用するのでぜひ覚えておきましょう。コメントをキャンセルメールアドレスが公開されることはありません。※サイト次回のコメントで使用するためブラウザーに自分の名前、メールアドレス、サイトを保存する。date*日本語が含まれない投稿は無視されますのでご注意ください。本記事では非線形最小二乗法の計算手法の1つである最急降下法について解説します。使用する数値計算ソフトの比較:メリット・デメリットを解説この記事では私が実務で使ったことがある、あるいはこれから使おうと考えている計算ソフトのメリット・デメリットを比較して紹介します。計算ソフトは1つにこだわる必要はなく、様々なソフトを使用することで新しい発見が生まれます。【Direct法(逐次代入法)】を解説:非線型方程式の解法Direct法は非線形方程式を解く手法の1つです。逐次代入法とも呼ばれており、初期値xから得られたyの値をそのまま次のステップのxの値に代入して計算する手法です。【Euler法】をわかりやすく解説:常微分方程式の数値解法Euler法(オイラー法)は常微分方程式を解く手法の1つです。十分に小さい刻み幅で差分を取ることにより、近似的に解を得ることができます。本記事ではEuler法の計算手法について解説し、実際にタンク内の濃度変化に関する例題を解いてみます。【ラグランジュ補間】をわかりやすく解説:データを多項式で補間する手法ラグランジュ(Lagrange)補間はあるデータを多項式で補間する手法の1つです。補間法の中でも基礎的な位置付けにあり、初心者はまずこの手法から勉強するのが良いでしょう。【Wegstein法】をわかりやすく解説:非線型方程式の解法Wegstein法は非線形方程式を解く手法の1つです。Direct法(逐次代入法)をベースとして定数qで補正することで収束性を改善しています。【Euler法】をわかりやすく解説:常微分方程式の数値解法プロフィールこのブログを運営しているルートです。化学メーカーで働くケミカルエンジニアです。