微分方程式

 流体の流れ現象は時間と位置の偏微分方程式であるナビエ−ストークスの式で完璧に記述される。これと同様に、時間に関係するあらゆる自然現象は微分方程式で記述される。化学の世界でも反応や移動現象は時間に依存するので微分方程式がモデルと解析の主役である。

 変数分離法による微分方程式の解法

微分方程式が次式のような形式:

で、f (x)とg (y)がそれぞれxとyのみの関数に分解できるとき、解は次式となる。

ここで、C は積分定数であり、初期条件(境界条件)から決定される。

なお、初期条件x =x0y =y0がわかっていれば、次式からでもよい。(この章ではこちらを使う)

微分方程式の境界条件

微分方程式の導関数の最高次数をその微分方程式の階数という。微分方程式の解(一般解)はその微分方程式の階数に等しい個数の任意定数を含む。従って微分方程式の解を得るには階数と同じ数の境界条件が必要である。例えば2階の微分方程式の(任意定数を含まない)解を得るには2個の境界条件が必要である。円管内流れの速度分布に関する基礎式は2階の微分方程式:

であるが、この一般解は 

である。これに2つの境界条件、u =0 at r =R @, du/dr=0 at r =0 A

を考慮することで放物線速度分布式:

が得られる。このような具体的な解(境界条件を満たすような解)を特殊解と呼ぶ。

 反応速度の微分方程式

反応速度とは濃度の時間微分である。

【1次反応】

ある成分の反応速度rが自身の濃度Cの一次式に比例する:

-(dC/dt)=kC

変数分離をおこない、

時間0で濃度C0として、両辺を積分する。

すなわち、

よって、 すなわち、

このように反応は一次である場合は、時間対が直線プロットになる。

【2次反応】

A+A→(生成物)

のような反応では反応物の濃度はその2乗に比例して減少する。

変数分離して積分する。

すなわち、

よって、 すなわち、

 

【実習】

右のプロットをおこない、1次反応と2次反応を比較する。

 

【平衡反応】気相反応など多くの反応では反応物が100%生成物に変化することはなく、反応物と生成物が平衡状態になるところで反応は止まる。最も簡単な平衡反応は次式。

この場合成分Aの濃度変化(反応速度式)は正反応と逆反応の和である。

Aが反応した割合を「転化率x」()であらわすと

 すなわち

ここで、となるxをとおくと、微分方程式は

となる。変数分離して、

左辺は置換積分により、

よって、解が、 となる。(

反応操作の微分方程式

一次反応の連続槽型反応器(CSTR)

液容積V[L]の完全混合槽にともに速度W [L/s]で流入・流出がある。流入流れ中の成分Aの濃度がCA0[mol/L]であり、槽内と流出流れ中の濃度がCA である。槽内では一次反応:

 [mol/(L-s)]

により成分Aが成分Cになる。成分Aについての物質収支式は次式。

これより、CAに関する微分方程式が次式となる。

この式はt =0でCA=CA(0)の初期条件で解くことができる。例えば、CA(0)=0とし、

の場合の解は、 である。

【例】 W= 0.15 L/s, CA0= 10 mol/L, V= 10 L,k=0.005, 初期条件t=0 : CA= 2 mol/L で解く。

【解】

微分方程式:

初期条件 t=0: CA=2 を解く。

変数分離して積分

 流体工学の微分方程式

円柱座標系におけるナビエ−ストークスの式を、円管内流れに適用して、簡略化すると次の微分方程式となる。

uは管軸方向の速度)

  ( または)

境界条件:u =0 at r =R @, du/dr=0 at r =0 A

 

を入れて上式を積分すると、積分定数をC1として

境界条件AからC1=0である。さらにもう一度積分して、境界条件@により

のように解ける。これが円管内流れの放物線速度分布である。

伝熱の微分方程式

【コーヒーの冷める早さ】「カップに80℃のコーヒーを入れた。10分後に50℃になった。30分後の温度を求めよ。外気温は20℃である。」という問題を考える。

コーヒーの温度をT [℃]、時間をt [min]とする。ニュートンの冷却法則により、冷却速度dT/dtは温度差(T -20)に比例する。これを書くと次式の微分方程式になる。

k が伝熱の大きさをあらわす定数である。変数分離して積分する。

これが微分方程式の一般解である。初期条件 t =0; T =80 より、C =60。さらに t =10; T =50より、

 より、k =0.069。よって微分方程式の解は

である。これよりt =30 minでのT = 27.6℃。

 物質収支の微分方程式

【単蒸留のレーリーの式】

蒸発缶(スチル)に原料を仕込み,発生蒸気を凝縮させてそのまま留出液とする回分操作が単蒸留である.蒸留酒の製造などに古くから用いられている.この操作では蒸留の進行とともに留出液中の低沸点成分の濃度も低下するので,留出率と留出液濃度の関係(蒸留曲線)が操作の要である.
 はじめに全量F
0[mol], 低沸点成分(メタノール)濃度 x0[モル分率]の2成分混合液をスチルに仕込む.単蒸留期間中のある時間におけるスチル中の全液量,濃度をF, x,留出液すなわち発生蒸気の低沸点成分濃度をy [モル分率]とする.その後の微小時間に仕込み液がdF[mol]蒸発し,その濃度がdx低下したとすると,低沸点成分の物質収支は,

dFdxの項を省略して,

これを蒸留開始の状態(Fo, xo)から,ある時間の(F,x)まで積分する。左辺は、

であるので、

これが単蒸留のレーリーの式である。仕込み液に対する留出蒸気の割合すなわち留出率θをとすると、

 

2成分系の平衡関係が,一定の相対揮発度α:

により表わせる場合には、右辺が

となる。よって、

 

である。

【実習】解析解によりレーリー式を計算する。dis5.xls

 

物質収支の微分方程式2

 

【濾過の基本 ルースの濾過方程式】

固形分を含む原液を濾過し、濾液を得る。濾過の膜面積をA[m2]とする。単位濾液量に対して得られる(湿潤)ケークの体積をc [m3-ケーク/m3-濾液]とすると、得られた濾液量V [m3]に対してケーク厚みLc [m]は、

である。Darcy則によると、粒子層を通る水の流束JV [m3/(m2・s)]は圧力差ΔPに比例し、水の粘性係数μおよび粒子層の厚さLに反比例する。kが抵抗係数である。

ここで、Lmは濾材(膜)の水の透過抵抗を相当するケーク厚みに換算したものである。(ここがポイント)濾液量Vと流束は、

   の関係なので、

である。濾材の抵抗相当のケーク厚みLmができるための濾液量をV0とすると、

 を考慮して、 である。

(この式の逆数をとると、なので、V対dt/dV のプロットでも濾過定数が求められる。)

この微分方程式を解くために、変数分離して両辺を積分する。

よって、次式となる。

  (

これをRuthの濾過方程式といい、K, V0を濾過定数という。ケークの濾過抵抗であるKと濾材(膜)の濾過抵抗に関係するV0がわかれば、濾過過程をモデル的にあらわせたことになる。ただしKの中に圧力差が含まれているので、この式は圧力が一定すなわち「定圧濾過」に関してのものである。(実際にはRuthは濾過実験のデータ解析により、実験値をあらわせるものとして、先ずこの式を発見し、その後に理論的に導いている。)

 

この式をVで解くと次式となる。

【例題】ある濾過で濾液量と時間との関係が表のようであった。

上式とデータをグラフ上で比較して濾過定数を求めよ。

時間 t [s]

6.8

19.0

34.0

53.4

76.0

102.0

131.2

163.0

濾液量
V[×10-3 m3]

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

 

【解】K= 9.0e-8, V0=1.0e-4でのモデルとの比較が右図。

<mf02.xls>

 

【クロスフローろ過の透過流束経時変化】

粒子濃度Cb[m3-固体/m3-水]の水をクロスフローろ過する。水の透過流束JV[m3-水/(m2・s)]はDarcyの式によりあらわされ、圧力差ΔP[Pa]に比例し、ケーク層厚みL[m]に反比例する。

ここで、kは比例定数、μは水の粘度。Lcが膜面に実際に堆積しているケーク層の厚みであり、Lmはろ過膜の透過抵抗をケーク厚みに換算したものである。ケーク層の厚みの変化速度は粒子の堆積速度とケークのはく離速度Jc [m3-固体/(m2・s)]との差による。

Jvを消去して、Lcの時間変化に関する微分方程式となる。

なお、ろ過の初めではLc=0なので、。ろ過流束の定常値で、

これはというかたちの常微分方程式である。変数分離して積分すると、

左辺の積分は

と、積分公式:

 

とから、解析解は次式となる。(

          <微分方程式 mf01.xls, 解析解 mf03.xls>

物質収支の微分方程式3

タンクとバルブの系の液面高さ

タンクの断面積をA、液面高さをhとすると、V=Ahである。流出ラインにバルブが設置されている。このような場合一般に流出量q’は液面高さhに比例しバルブの抵抗Rに反比例する。

物質収支式を液面高さhで書くと、

 すなわち、

【例】底面積A=2 m2のタンクがある。初めは流入量q=1 m3/sの状態で液面高さがh= 5 mで定常状態にあった。t=0でq=2 m3/sに変えると液面高さはどうなるか。

【解】基礎式はである。初めの状態は定常なので、(dh/dt)=0, h=5, q=1 より、

t>0では、微分方程式:

 初期条件: t=0: h=5

を解く。 

変数分離して、 

 

積分公式

 

より、

 

 ⇒  

 

図のように流量qをステップ状に変化させると、次の定常状態に指数関数的に変化する。流量を入力、液面高さhを出力とみれば、ステップ入力に対する指数関数的応答は多くの現象にみられる。このような系を一般に一次遅れ系とよぶ。

振動系の微分方程式

【単振動の微分方程式】ばねがつるされており、ここにMの質量の物体をつり下げるとsだけ伸びた。フックの法則より、

である。(k: バネ定数、g:重力加速度)この初期状態からxだけ変位したとき、おもりに働く力は下向きMg、上向きk(s+x)だから、(質量)×(加速度)=(力)より次式となる。

 すなわち、 

これを、とおいて、

とあらわしておく。

この2階の常微分方程式を解く。解を試みに

とおいてみる。

を代入すると、

よって、を満たすzをとれば、は解である。の解は

であるから、解として、

となる。これを実数とするために、複素数の世界で指数関数と三角関数の同等性を示すオイラーの公式:

を使う。よって、

これにより、x1,x2の和、差が実数値の関数になる。

これらはもとの微分方程式の解であり、独立である。よって解(一般解)はこのようになる。

 

実際の解(特異解)はt=0における2つの条件により定まる。例えば、

1. (C1だけ引っ張ってはなす場合。)

2. (おもりを押す場合)

など各種の場合が考えられる。

(ph01.xls)

振動を抑制するためにおもりにダッシュポット(油の粘性により振動を抑える機器)をつけた場合は、物体の速度に比例する力が運動と反対方向に働くので、基礎式は次式となる。

すなわち、(

以上は自由振動であったが、さらにおもりに時間変化する外力f(t)が加わる場合(強制振動)の微分方程式は、

となる。