VBAでの拡散計算
エクセルのVBAを使って添付画像のようなグラフを作成しようと考えています。
以下の計算で作成できそうなのですが、⊿tを10-5より小さく設定し、1000sec～の濃度変化が知りたいため表計算ではなくVBAを使ってみました。
表計算では、
t=0のときx0=C0(飽和濃度)、x>x0でC=0とし(初期条件)
x0では、
C(j+1,i) = C(j,i) +D * ( C(j,i-1) - C(j,i) ) / dx / dx * dt
x >x0では、
C(j+1,i) = C(j,i) +D * ( C(j,i-1) - 2 * C(j,i) + C(j,i+1) ) / dx / dx * dt
の計算を行い、セル表記が以下のようになりました。どの時間も物質量は一定です。(たぶん・・・)
t0、t1、t2、t3t、・・・
x1,60、58.8、57.6・・・
x2, 0、 1.2、2.38・・・
x3, 0、 0、0.02・・・
上の計算をVBAで以下のように書いてみました。
Sub diffusion_dry_up2()
Dim n As Integer, nt As Integer
Dim i As Integer, j As Integer
Dim b As Double, te As Double, dt As Double
Dim c0 As Double, cs As Double, d As Double
Dim x As Double, dx As Double, t As Double
Dim a As Double, cjp As Double, cj0 As Double
Dim cjm As Double
b = InputBox("配管の長さb(m)")
n = InputBox("配管の長さの方向分割数n")
te = InputBox("計算する時間長t(sec)")
dt = InputBox("時間増分dt(m)")
c0 = InputBox("配管底部の濃度c0(vol.%)")
cs = InputBox("時刻t=0の時の配管内の濃度cs(vol.%)")
d = InputBox("拡散係数d(m^2/sec)")
Sheet1.Cells(1, 2) = "配管の長さb(m)"
Sheet1.Cells(2, 2) = "配管の長さの方向分割数n"
Sheet1.Cells(3, 2) = "計算する時間長t(sec)"
Sheet1.Cells(4, 2) = "時間増分dt(sec)"
Sheet1.Cells(5, 2) = "配管底部の濃度c0(vol.%)"
Sheet1.Cells(6, 2) = "時刻t=0の時の配管内の濃度cs(vol.%)"
Sheet1.Cells(7, 2) = "拡散係数(m^2/sec)"
Sheet1.Cells(1, 3) = b
Sheet1.Cells(2, 3) = n
Sheet1.Cells(3, 3) = t
Sheet1.Cells(4, 3) = dt
Sheet1.Cells(5, 3) = c0
Sheet1.Cells(6, 3) = cs
Sheet1.Cells(7, 3) = d
nt = te / dt
dx = b / n
Sheet1.Cells(1, 1) = nt
t = 0
a = d * dt / dx ^ 2
Sheet1.Cells(1, 5) = t
Sheet1.Cells(1 + 1, 4) = -dx
Sheet1.Cells(1 + 2, 5) = c0
Sheet1.Cells(1 + n + 2, 4) = b + dx
Sheet1.Cells(1 + n + 3, 5) = cs
Sheet1.Cells(3, 5 + i) = a * (-cj0 + cjp) + cj0
For i = 1 To nt
t = t + dt
Sheet1.Cells(1, 5 + i) = t
Sheet1.Cells(1 + n + 3, 5 + i) = a * (cjp - 2 * cjo + cjm) + cj0
For j = 2 To n + 2
cjp = Sheet1.Cells(1 + j + 1, 5 + i - 1)
cj0 = Sheet1.Cells(1 + j + 0, 5 + i - 1)
cjm = Sheet1.Cells(1 + j - 1, 5 + i - 1)
Next j
Next i
linegraph 2, 4, 4 + n, nt + 2
End Sub
Sub linegraph(sr As Integer, sc As Integer, lr As Integer, lc As Integer)
ActiveSheet.ChartObjects.Add(200, 10, 240, 200).Select
ActiveChart.ChartWizard _
Source:=Range(Cells(sr, sc), Cells(lr, lc)), _
gallery:=xlXYScatter, _
Format:=3, _
PlotBy:=xlColumns, _
categoryLabels:=1, _
SeriesLabels:=0, _
HasLegend:="false", _
Title:="ex2", _
categoryTitle:="t", _
ValueTitle:="y", _
ExtraTitle:=""
End Sub
しかしまったく表計算のようになりませんでした。
a = d * dt / dx^2以降の書き込みが変だと思うのですが、どのようにすればよいのでしょうか。
また、上のような表記ではtを大きくdtを小さくするとエラーになってしまいます。
質問項目が多いですが、よろしくお願いします。
