タンクやお風呂への貯水や水抜きにかかる時間を予測する

化学工場の操作の一つにタンクへの貯水や水抜きがあります。
また、液面を所望の高さにするためにどのように流体を流入させたり流出させたりすればいいのか考えたり、制御系を組んでその仕組みを自動化させたりします。

身近な現象ではお風呂に水を貯めるのにどれくらいの時間がかかるのか、お風呂の水抜きにどれくらいの時間がかかるのか考えたことはあると思います。
貯水は単なる掛け算で計算できますが、抜水は微分方程式を解いて求めなければいけない問題になります。
水位が高ければ高いほど流出流量は多く、そしてその水位は時間変化するからです。

本記事ではタンクやお風呂に水を貯める・水抜きをする、そしてその速度をコントロールして液面の高さを所望の高さにすると言ったことを目的に
ある流入流量とバルブ抵抗（≒バルブの開度）を与えたときに、タンクの水位がどのように変化していくのかを計算してみたいと思います。

問題設定

①低面積30m²、高さ10mの空タンクに対して、流量q_in = 100 m³/hで水を貯めたい。高さ8mに達するまでの時間を求めよ。
②上記と同じ空タンクにおいて、流量q_in = 100 m³/h、バルブの抵抗を0.08とした。このタンクの水位の時間変化を求め、定常状態となったときの液面の高さを求めよ。
③前問で達成した定常状態において、流量q_in = 50 m³/hに変更した。このタンクの水位の時間変化を求めよ。

まず、状況を図に起こしたいと思います。それぞれの問題のイメージ図は以下の通りです。

①は単なる貯水、②は排水しながら水を入れ続ける、③は②で水位が安定しているところで流量を変更して、その水位への影響を調べる（ステップ応答問題）という状況設定です。

物質収支

あるプロセスに作業を与えて、そのプロセスの応答がどのようになるのかを考える時（今回は、ある流量、バルブ抵抗でタンクに水を入れたとき、液面の高さがどのようになるのかを考えます）、そのプロセス周りで物質収支をとります。タンクに入れる流体の密度をρ[kg/m³]、タンクに貯まる水の体積をV[m³]としたとき

タンク内流体の蓄積量の時間変化＝流入量-流出量

という物質収支式が立ちますので、タンク内水量のモデル式は以下の通りです。（※両辺の単位はkg/s）

ここで、水面の高さをy[m]とすると、タンク内の水の体積＝タンクの底面積×水面の高さとなります。さらに流出流量q_outはバルブ抵抗R[hr/m²]を用いて、以下の様に定義します。

基本的にこの問題において、時間に依存する関数は液面の高さのみであり、流出流量は液面の高さに比例する。という考え方から、上記のような置き換えを行っています。上記に従い式変形すると

となり、解くべき微分方程式を導出することができました。（※両辺の単位はm/hr）

物質収支を解く

ルンゲクッタ法

これから、先ほど導いた物質収支をルンゲクッタ法を用いて解きます。ルンゲクッタ法に関する詳しい説明はこちら
[blogcard url=”http://shimaphoto03.com/science/rk-method/”]

右辺をf(t,y)とし、以下のアルゴリズムに従って、解曲線y=y(t)を求めます。

ソースコード

ルンゲクッタによる微分方程式の求解にはエクセルのVBAマクロを使用しました。

/* –level formula　rk-method ——– */
Public a, q, R As Double

Sub tori()

Dim y0, init, ed, h, h2, ky(4) As Double
Dim i, j As Integer

‘Calculate properties
init = Cells(1, 3)
ed = Cells(2, 3)
h = Cells(3, 3)
h2 = Cells(4, 3)

‘Constants
a = Cells(5, 3)
q = Cells(6, 3)
R = Cells(7, 3)

‘Initial Condition
y0 = Cells(4, 7)
y = y0
t = init

‘Runge-kutta roop

For i = 0 To ((ed – init) / h)

j = 6 + i / h2

Cells(1, 6) = i
Cells(2, 6) = j – 6

If i Mod h2 = 0 Then

Cells(j, 6) = t
Cells(j, 7) = y

End If

ky(1) = h * F1(t, y)
ky(2) = h * F1(t + h / 2, y + ky(1) / 2)
ky(3) = h * F1(t + h / 2, y + ky(2) / 2)
ky(4) = h * F1(t + h, y + ky(3))
ny = y + (ky(1) + 2 * ky(2) + 2 * ky(3) + ky(4)) / 6
nt = t + h

t = nt
y = ny

Next i

End Sub

Function F1(ByVal t As Double, ByVal y As Double) As Double
F1 = (1 / a) * q – (1 / (a * R)) * y
End Function

/* –level formula　rk-method End——– */

上記のソースコードを以下のエクセルシートで走らせます。

問題を解く

問題①

①低面積30m²、高さ10mの空タンクに対して、流量q_in = 100 m³/hで水を貯めたい。高さ8mに達するまでの時間を求めよ。
貯水の問題です。水を貯めていますから、バルブは締めた状態にします。この時、右辺第二項（1/AR）*y=0となりますから、解くべき方程式が変わってしまうのですが、今回は疑似的にR=10と置くことで、十分にバルブ抵抗があると考えて解きます。
※R=5~10まで変化させて解曲線に変化がないことを確認したうえでR=10としています。

答えは以下の様になります。

当たり前じゃんと言われれば当たり前なのですが・・・。
答えとしては高さ8mに達するまでの時間は2.4時間です。

ただし、タンクから流体を溢れさせたら大惨事ですので、実際には制御系（PI、PID制御）を組んで操作します。
[blogcard url=”http://shimaphoto03.com/science/vba-pid/”]

問題②

②上記と同じ空タンクにおいて、流量q_in = 100 m³/h、バルブの抵抗を0.08とした。このタンクの水位の時間変化を求めよ。
バルブを開けながら水を貯めていきます。バルブの抵抗を0.08に変えて再度ルンゲクッタ法で計算します。

今度は、直線ではなく、カーブを描きながら水面の高さが変化していることが分かります。これは、立てた微分方程式の右辺第二項にyの関数が現れたためです。
そして、バルブを開けながら水を貯めるとある高さで一定になることが分かります。この状態になったプロセスのことを「定常状態になった」と表現します。
このプロセスでは、定常状態における液面の高さは8mです。

問題③

③前問で達成した定常状態において、流量q_in = 50 m³/hに変更した。このタンクの水位の時間変化を求めよ。
②において、流量q_in = 100 m³/hで水を貯めながらバルブ抵抗を0.08としたとき、8mで水面が落ち着く（定常になる）ということがわかりました。この状態で、流量を50 m³/hに変更したらどのようになるのか？という問題です。

先ほどのエクセルシートにおいて、Ｇ４セルのy0を8に変更し、qを50に変更して、ルンゲクッタ法で計算します。

つまり、液面高さの初期条件を8ｍとして再度微分方程式を解くということです。

答えは以下のようになります。

10時間もの時間をかけて、水位が4ｍまで落ちるという計算結果になりました。

プロセス制御

これまで解いた問題は制御という操作を全く行わなかったときにどうなるか？を考えていました。
制御という操作を行わないと、例えば問１のような状況で流出バルブを締めて貯水を始め、流入バルブを開けっぱなしにしていたら、タンクから流体が溢れてしまったという惨事を招きます。特に流体が毒劇物だったり石油精製物だったら危険です。

こういったことを防ぐためにプロセスには自動制御系が組まれています。次回の記事では、この自動制御系の仕組みについてまとめてみたいと思います。
[blogcard url=”http://shimaphoto03.com/science/vba-pid/”]