以下の拡散方程式を種々の初期条件、境界条件の下に解

質問

以下の拡散方程式を種々の初期条件、境界条件の下に解く。とのことですが、力を貸していただければ幸いです。

＊ FOR DIFFUSION EQUATION BY KOUSUKE TANAKA
＊ NFUNC:FUNCTION NUMBER
＊ N:NO.OF DIVISIONS
＊ T:TIME VARIABLE
＊ DT:TIME INCREMENT
＊ TMAX:UPPER LIMIT OF TIME
＊ H:SPATIAL MESH SIZE
＊ C1,C2,C3:COEFFICIENT MATRIX
＊ FNU:DIFFUSIVITY
＊ Q:PAINT THICKNESS
PARAMETER(NDIM=100)
DIMENSION A(NDIM),B(NDIM),C(NDIM)
WRITE(＊,＊)’INPUT:FNU=’
READ(＊,＊)FNU
WRITE(＊,＊)’INPUT:DT=’
READ(＊,＊)DT
WRITE(＊,＊)’INPUT:TMAX=’
READ(＊,＊)TMAX
WRITE(＊,＊)’INPUT:Q=’
READ(＊,＊)Q
H=1.0/FLOAT(N)
＊ INITIAL CONDITION
DO I=1,N
ここまでも合っているかわかりませんが、ここからが全く手に負えません。

厚さがq(μm)あるペイントの表面の酸素濃度を周期的に変えますC(t)。

拡散方程式 式(1)
C：気体濃度、Z：厚さ方向、t：時間、Dm(μm2/s)：拡散係数
というのがあり、これを中心差分近似したものが、 式(2)

となります。この拡散方程式に加えて、初期条件、境界条件を以下のように加えます。
・初期条件
t=0;C=0(初期状態は厚さ方向に一様の真空状態)
・境界条件
ペイントの底面付近では厚さ方向にて濃度が同じ。
→式(3)
この条件を離散化したものが、
→式(4)とおそらくなります。
また、境界条件といたしまして、表面付近の濃度変化を正弦関数C(t)=Ksin2πft
(K：0より大きい定数、f：周波数、t：時間)
として与えます。
問題は、ペイント表面の気体濃度を正弦関数的に変化させると、ペイント内において拡散係数に応じた気体の拡散が生じます。ペイント内のある時間での気体の厚さ方向での濃度分布、かつ、その厚さ方向の濃度を合計したペイント全体の濃度の時間変化をfortranでプログラミング化するのですが。。。

その他の条件
ペイントを厚さ方向にN(=100)等分する。
ペイントの厚さpaint thicknessをq(μm)とする。
時間増分time incrementの設定
とりうる時間の最大値upper limit of timeの設定
ペイントの厚さ方向における一つの区切りspatial mesh sizeの設定
拡散係数diffusivityの設定

Fortran上では
NFUNC：FUNCTION NUMBER
N：NO.OF DIVISIONS
T：TIME VARIABLE
DT：TIME INCREMENT
TMAX：UPPER LIMIT OF TIME
H：SPATIAL MESH SIZE
FNU：DIFFUSIVITY

回答No.5

nagaihdhr

絵で与えられている式の境界条件はディリクリ条件で、質問の境界条件を考慮するとノイマン条件で解くことになるのではないでしょうか。ということで、次のようの感じでいかがでしょうか。

! for diffusion equation　by　kousuke tanaka

program diffuse
implicit none

real, parameter:: pai = 4.0 * atan(1.0)

integer n ! no.of divisions
real t ! time variable
real dt ! time increment
real tmax ! tmax:upper limit of time
integer jmax ! no. of time division
real h ! spatial mesh size
! c1,c2,c3:coefficient matrix
real q ! paint thickness
real fnu ! diffusivity
real k ! constant
real f ! frequency of density
real alpha !
real, allocatable, dimension(:, :) :: c
! diffusion density
integer i, j

call init()

do j=0 , jmax-1
do i=-1 , n+1
if (i == -1) then
c(i,j+1) = (1-2*alpha)*c(i,j) + 2.0*alpha*c(i+1,j)
else if (i == n+1) then
t = j * dt
c(i,j+1) = (1.0-2.0*alpha)*c(i,j) + 2.0*alpha*c(i-1,j) + 2*alpha*K*sin(2.0*pai*f*t)*h
else
c(i, j+1) = alpha*c(i-1,j) + (1.0 - 2.0*alpha)*c(i,j) + alpha*c(i+1, j)
end if
end do
end do

call print()

contains
subroutine init()
integer ii

write(*,fmt='(a)',advance='no') 'input:q='
read(*,*) q

write(*,fmt='(a)',advance='no') 'input:n='
read(*,*) n

h = q / float(n)

write(*,fmt='(a)',advance='no') 'input:fnu='
read(*,*) fnu

write(*,fmt='(a)',advance='no') 'input:tmax='
read(*,*) tmax

write(*,fmt='(a)',advance='no') 'input:dt='
read(*,*) dt

write(*,fmt='(a)',advance='no') 'input:k='
read(*,*) k

write(*,fmt='(a)',advance='no') 'input:f='
read(*,*) f

alpha = (fnu * dt) / (h * h)
if (alpha >= 0.5) then
write(*,*) 'Value emits it and cannot calculate.'
stop -1
end if

jmax = int(tmax / dt)

allocate( c(-1:n+1, 0:jmax) )
do ii=-1 , n+1
c(ii, 0) = 0.0
end do
end subroutine init

subroutine print()
integer ii, jj

open(11,file='out.d',status='unknown')
do jj=0 , jmax
write(11,*) jj, (c(ii,jj),ii=-1,n+1)
end do
close(11)
end subroutine print
end program diffuse

i, jの添え字が逆かも知れませんが、あとは適当に考えてください。

お礼

かなり詳しくてわかりやすいです。。。

回答No.4

sat000

参考にどうぞ、私が今まで書いたことをまとめました。
???には適当なものを入れてください。
ちゃちゃっと数分で書いたので間違ってたらごめんなさい。

real*8:: cold(10000), cnew(10000)
data pi/3.14159d0/
data dt/ ???/
data nt/ ???/
data dz/ ???/
data nｚ/ ???/
data dm/ ???/

do ii = 1, nz
cold(ii) = 0.0d0
cnew(ii) = 0.0d0
end do
time = 0.0d0
open(10, file = '???', status = 'unknwon' )
write(10, fmt='(???)') time, ( cnew(ii), ii = 1, n )
do ii = 1, nt
time = dt * dble( nt )
cnew(nz) = K * dsin( 2.0d0 * pi * f * time )
do jj = nz - 1, 1, -1
if ( jj .eq. 1 ) then
cnew(jj) = cold(jj) + dm * dt / dz / dz * 2.0d0 * ( cold(jj + 1) - cold(jj) )
else
cnew(jj) = cold(jj) + dm * dt / dz / dz * ( cold(jj + 1) - 2.0d0 * cold(jj) + cold(jj - 1) )
end
end do
write(10, fmt='(???)') time, ( cnew(ii), ii = 1, n )
do ii = 1, nz
cold(ii) = cnew(ii)
end do
end do
close(10)

お礼

返信ありがとうございます！！うれしいです！！

補足

返信ありがとうございます。
申し訳ないのですが、、、これは私の質問に対するプログラムですか？
知識不足で申し訳ないです。。。

回答No.3

sat000

> 左辺のCの添え字がt+1,i/t,i
> 右辺のCの添え字がt,i+1/t,i/t,i-1となります。
じゃあ、私が書いた式のiをt、ｊをiと読み替えてください。
それ以外は問題無く使えるはずです。

底を物理的に抜けてこないという意味で、マスフラックスが0です。
そして差分化と、底の境界に適用すべき差分方程式は前回書いてありますので、もう一度良く見てください。

> この式(4)の意味は、底面側は表面からの気体の拡散の影響をほとんど受けないという意味(物理的な底だから)
影響をほとんど受けないという意味がちょっと不明ですが、底を気体なりなんなりが通過しないという意味ですかね？
だとしたら、上に書いたマスフラックスが0という記述をします。
「ほとんど」受けないという曖昧な記述を数値解析は受け付けません。
受けないのか、ちょっとでも受けるのならその影響をどう記述するのか、受けても無視できるくらい小さいから考えないことにするのか、そのあたりがモデリングの腕です。
底を透過しないのなら、マスフラックス0と普通言います。
式は前回書いてますが、別の表現をすると、濃度勾配が0と言っても良いです。

一方、上の境界はsinで変化するものの際限なく入ってきますから、内部にはどんどん溜まっていきます(sinの影響で濃度分布にムラみたいなものが生じるはずだが)。
大雑把にはそういう時間変化をします。

お礼

詳しく原理の方を教えていただきありがとうございます！！うれしいです！！

補足

この式(4)の意味は、底面側は表面からの気体の拡散の影響をほとんど受けないという意味(物理的な底だから)
影響をほとんど受けないという意味がちょっと不明ですが、底を気体なりなんなりが通過しないという意味ですかね？
→大変申し訳ありません。濃度勾配が０という意味です。ペイントの底面では厚さ方向による濃度勾配がないのです。逆に底面以外は場所が少しずれるだけで濃度が違ってくる(拡散の影響を受けやすい)のでしょう。
今回扱う問題が大変難しいので許してください。質問への回答本当にありがとうございます。

回答No.2

sat000

図がというか式が解像度が悪くて読み取れないのですが、式(2)の右辺のCの添字がi,j+1、i,j、i,j-1に見えますけど、それで良いんですか？
普通、こういう書き方をすると、二次元を表すんですが、対象は一次元ですよね？
式(2)の左辺のCの添字はi+1,jとi,jですか？
ふーむ、だとするとiは時間、jが場所を表していると見るしかないですね。
あまり見かけない書き方ですが。
仮にiを時間、jを場所とすると、それはそのまま解けば良くて、
C(i+1,j)=C(i,ｊ）＋(D dt)/(dz)^2*(C(i,j+1)-2 C(i,j)+C(i,j-1)) --- (a)
の形にすれば分かりやすいですかね。
陽解法になるので、(D dt)/(dz)^2が1/2よりも大きいと発散する可能性があります(普通は計算時間との兼ね合いの中で、できるだけ小さくする)。
で、式(3)ですが、底面側がマスフラックス0というのは物理的な底だから良いと思いますが、なぜ式(4)では式が2個あるんですか？
しかも式(3)のdC/dzの係数が式(2)のDmと違うように見えますけど？おまけに左辺がzに見える気がしますが、これは図の解像度のせいかもしれませんね。
マスフラックスgamma = -D dC/dz とすると、差分化して、gamma = -D (C(i,j+1)-C(i,j-1))/(2 dz) = 0
より、C(i,j+1)=C(i,j-1)です。
これを式(a)に代入すると、
C(i+1,j)=C(i,ｊ）＋(D dt)/(dz)^2*(2 C(i,j+1)-2 C(i,j))
となり、これが底の境界に適用すべき差分方程式になります。
差分法は、境界と内部で使う式が変わることが多いので、その点がちょっと不便ですかね。
一方、上側の境界は、C(t)=Ksin2πft で与えられている固定境界の一種なので、それをそのまま使えば良いです。
後は逐次解いていけば良いでしょう。
ただ、中心差分だと、拡散のような場合にちょっとまずいことも起こり得ます。
が、まあそういうことも含めて、まずはやってみて色々経験を積むことが重要です。

お礼

返信ありがとうございます！！！

補足

図では
左辺のCの添え字がt+1,i/t,i
右辺のCの添え字がt,i+1/t,i/t,i-1となります。
添え字左は時間、右が場所となります。

私としては、Cの添え字は左から場所、そして時間としたかったので。。。申し訳ありません。

式(3)ですが、底面側がマスフラックス0というのは式(3)ですが、底面側がマスフラックス0というのは物理的な底だから良いと思いますが、なぜ式(4)では式が2個あるんですか？
→この式(4)の意味は、底面側は表面からの気体の拡散の影響をほとんど受けないという意味(物理的な底だから)、を示すために私個人で考えた条件なので正直正しいかわかりません。。。
式(３)の条件をfortranで示すためにはどうしたらよいのでしょうか？お答えしていただくと幸いです。

この問題は、「ある気体を透過しない母体に、気体透過性のペイントを貼った場合、その透過性のペイントの表面の気体濃度を正弦的に変化させた場合、ペイント内部ではどのように気体濃度が変化していくかということを見るためのプログラム」作りたいのです。
自分の知識不足大変申し訳ないです。

回答No.1

FT56F001

拡散方程式を解くプログラムがあって，それを改造して別の境界条件で解け，という課題でしょうか。
しかし，この質問では答えようがないです。
質問者さんが示したプログラムは，コメントとパラメータ入力の部分だけで，計算をする本質部分のコードではありません。

まずは現在のプログラムを解読して，どこで拡散方程式を解いているのか，どこで境界条件を与えているのかを突き止める必要があります。そこに今回必要となった新しい条件を書き足す形になります。
プログラムリストは数ページ～数十ページ以上のかなり大きなプログラムでしょうから，データ入力部分，出力部分，境界条件を与える部分，方程式を解く部分などに区分けされているはずです。これを手がかりにして，境界条件を与えている部分を突き止めます。境界条件を与える変数の意味などを，プログラム中のコメントとコードを手がかりに，解読します。数学的に境界条件を与える部分はほんの数行で，エラー処理や場合分けなどの本質的でない部分に埋もれているのが通例です。こういった解読作業が一通り終わってから，改造作業にかかります。

他人が書いたプログラムの解読は，けっこう難しい作業です。変数の使い方など，発想が人によって違うからです。(自分で書いたプログラムですら，しばらく放置するとさっぱり分からなくなります)。現在のプログラムの規模とドキュメンテーションの具合によっては，解読をあきらめて，初めから作り直す方が早い場合すらあります。全部作り直すくらいのつもりで，解読作業に当たってください。

一つの方法として，理論解がわかっている一次元の簡単な拡散方程式を中心差分などで解くプログラムを一から書いてみる練習から入る手もあります。一度自分でプログラムを書いてみると，他人のプログラムも読めるようになります。

補足

プログラミングは努力は惜しめないということですね。。。頑張ります。

拡散方程式(２階偏微分方程式)の離散化は自分で行えます。
その拡散方程式を種々の初期条件、境界条件の下にで解くプログラム
が作りたいのです。
私がわからないのは、拡散方程式や種々の条件をfortranにて示す方法が
まったくわかりません。。。種々の条件は図にて示しております。この条件も自分で考えたものなので正しいかどうかわかりませんが↓↓
周りに相談できる方がいないので本当に困っています↓↓
こういうように質問に反応していただけるだけで嬉しい限りです。
この問題は、「ある気体を透過しない母体に、気体透過性のペイント
を貼った場合、その透過性のペイントの表面の気体濃度を正弦的に変化
させた場合、ペイント内部ではどのように気体濃度が変化していくか
ということを見るためのプログラム」作りたいのです。
言い換えるならば、縦軸：ある時間におけるペイント全体の濃度の合計
横軸：0<=t<=tm(tm→とりうる時間の最大値)
を作りたいのです。