以下の拡散方程式を種々の初期条件、境界条件の下に解く。とのことですが、力を貸していただければ幸いです。
* 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
A 回答 (5件)
- 最新から表示
- 回答順に表示
No.5
- 回答日時:
絵で与えられている式の境界条件はディリクリ条件で、質問の境界条件を考慮するとノイマン条件で解くことになるのではないでしょうか。
ということで、次のようの感じでいかがでしょうか。! 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
- 回答日時:
参考にどうぞ、私が今まで書いたことをまとめました。
???には適当なものを入れてください。
ちゃちゃっと数分で書いたので間違ってたらごめんなさい。
real*8:: cold(10000), cnew(10000)
data pi/3.14159d0/
data dt/ ???/
data nt/ ???/
data dz/ ???/
data nz/ ???/
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
- 回答日時:
> 左辺のCの添え字がt+1,i/t,i
> 右辺のCの添え字がt,i+1/t,i/t,i-1となります。
じゃあ、私が書いた式のiをt、jをiと読み替えてください。
それ以外は問題無く使えるはずです。
底を物理的に抜けてこないという意味で、マスフラックスが0です。
そして差分化と、底の境界に適用すべき差分方程式は前回書いてありますので、もう一度良く見てください。
> この式(4)の意味は、底面側は表面からの気体の拡散の影響をほとんど受けないという意味(物理的な底だから)
影響をほとんど受けないという意味がちょっと不明ですが、底を気体なりなんなりが通過しないという意味ですかね?
だとしたら、上に書いたマスフラックスが0という記述をします。
「ほとんど」受けないという曖昧な記述を数値解析は受け付けません。
受けないのか、ちょっとでも受けるのならその影響をどう記述するのか、受けても無視できるくらい小さいから考えないことにするのか、そのあたりがモデリングの腕です。
底を透過しないのなら、マスフラックス0と普通言います。
式は前回書いてますが、別の表現をすると、濃度勾配が0と言っても良いです。
一方、上の境界はsinで変化するものの際限なく入ってきますから、内部にはどんどん溜まっていきます(sinの影響で濃度分布にムラみたいなものが生じるはずだが)。
大雑把にはそういう時間変化をします。
この回答への補足
この式(4)の意味は、底面側は表面からの気体の拡散の影響をほとんど受けないという意味(物理的な底だから)
影響をほとんど受けないという意味がちょっと不明ですが、底を気体なりなんなりが通過しないという意味ですかね?
→大変申し訳ありません。濃度勾配が0という意味です。ペイントの底面では厚さ方向による濃度勾配がないのです。逆に底面以外は場所が少しずれるだけで濃度が違ってくる(拡散の影響を受けやすい)のでしょう。
今回扱う問題が大変難しいので許してください。質問への回答本当にありがとうございます。
No.2
- 回答日時:
図がというか式が解像度が悪くて読み取れないのですが、式(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,j)+(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,j)+(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)の意味は、底面側は表面からの気体の拡散の影響をほとんど受けないという意味(物理的な底だから)、を示すために私個人で考えた条件なので正直正しいかわかりません。。。
式(3)の条件をfortranで示すためにはどうしたらよいのでしょうか?お答えしていただくと幸いです。
この問題は、「ある気体を透過しない母体に、気体透過性のペイントを貼った場合、その透過性のペイントの表面の気体濃度を正弦的に変化させた場合、ペイント内部ではどのように気体濃度が変化していくかということを見るためのプログラム」作りたいのです。
自分の知識不足大変申し訳ないです。
No.1
- 回答日時:
拡散方程式を解くプログラムがあって,それを改造して別の境界条件で解け,という課題でしょうか。
しかし,この質問では答えようがないです。
質問者さんが示したプログラムは,コメントとパラメータ入力の部分だけで,計算をする本質部分のコードではありません。
まずは現在のプログラムを解読して,どこで拡散方程式を解いているのか,どこで境界条件を与えているのかを突き止める必要があります。そこに今回必要となった新しい条件を書き足す形になります。
プログラムリストは数ページ~数十ページ以上のかなり大きなプログラムでしょうから,データ入力部分,出力部分,境界条件を与える部分,方程式を解く部分などに区分けされているはずです。これを手がかりにして,境界条件を与えている部分を突き止めます。境界条件を与える変数の意味などを,プログラム中のコメントとコードを手がかりに,解読します。数学的に境界条件を与える部分はほんの数行で,エラー処理や場合分けなどの本質的でない部分に埋もれているのが通例です。こういった解読作業が一通り終わってから,改造作業にかかります。
他人が書いたプログラムの解読は,けっこう難しい作業です。変数の使い方など,発想が人によって違うからです。(自分で書いたプログラムですら,しばらく放置するとさっぱり分からなくなります)。現在のプログラムの規模とドキュメンテーションの具合によっては,解読をあきらめて,初めから作り直す方が早い場合すらあります。全部作り直すくらいのつもりで,解読作業に当たってください。
一つの方法として,理論解がわかっている一次元の簡単な拡散方程式を中心差分などで解くプログラムを一から書いてみる練習から入る手もあります。一度自分でプログラムを書いてみると,他人のプログラムも読めるようになります。
この回答への補足
プログラミングは努力は惜しめないということですね。。。頑張ります。
拡散方程式(2階偏微分方程式)の離散化は自分で行えます。
その拡散方程式を種々の初期条件、境界条件の下にで解くプログラム
が作りたいのです。
私がわからないのは、拡散方程式や種々の条件をfortranにて示す方法が
まったくわかりません。。。種々の条件は図にて示しております。この条件も自分で考えたものなので正しいかどうかわかりませんが↓↓
周りに相談できる方がいないので本当に困っています↓↓
こういうように質問に反応していただけるだけで嬉しい限りです。
この問題は、「ある気体を透過しない母体に、気体透過性のペイント
を貼った場合、その透過性のペイントの表面の気体濃度を正弦的に変化
させた場合、ペイント内部ではどのように気体濃度が変化していくか
ということを見るためのプログラム」作りたいのです。
言い換えるならば、縦軸:ある時間におけるペイント全体の濃度の合計
横軸:0<=t<=tm(tm→とりうる時間の最大値)
を作りたいのです。
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- 物理学 移流熱拡散方程式の解き方 フーリエ変換 1 2022/08/15 15:25
- Excel(エクセル) Excel2007での条件付き書式について 6 2023/05/02 10:56
- 物理学 量子力学 球面調和関数 導出 方位角成分 微分方程式の解 2 2022/07/02 13:40
- 工学 例えば タンクのなかにVm3という液体がありますし、この液体の濃度はAとなるとき このタンクに 他の 1 2023/03/09 21:46
- 物理学 物理の証明問題についての質問です。 平面内を運動する小球がある。この物体にかかる加速度の方向と大きさ 2 2023/05/16 00:28
- 数学 知ってる式から応用させて 微分式を作りたいですが 例えば タンクのなかにVm3という液体がありますし 6 2023/03/09 21:45
- 高校 対数方程式につきまして 4 2022/05/05 07:55
- Excel(エクセル) 図のような散布図の作り方を教えてください。 1 2022/07/19 11:52
- 物理学 時速 54 km で一直線上を車で走っていると、25 m 先に障害物を見つけた。このとき、どのくら 7 2022/05/29 12:09
- 数学 放物型偏微分方程式 ∂u/∂t =α^2 (∂^2u/∂x^2) + xcost (0<x<1) 境 1 2022/12/29 13:54
おすすめ情報
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
加速度 a=dv/dt = (d^2 x) /dt^2
-
高校・大学の物理(力学)教え...
-
v^2-v0^2=2ax 今日この式を習っ...
-
力学について質問です。 1.棒の...
-
シュヴァルツシルト時空での時...
-
dH/dtとH(t)の関係
-
EXCEL上の数字を自動で振り分け...
-
なぜ力積を時間に関して積分す...
-
質量流量の記号「・ の読み方を...
-
仕事率の表し方について。
-
Id²θ/dt²=-mghsinθの厳密解の...
-
今基礎物理学の問題を解いてい...
-
運動方程式の微分積分の計算
-
物理で微積をつかう。
-
高校物理 授業でこういうのをや...
-
摩擦を考慮したサイクロイド曲線
-
d^2r/dt^2の意味
-
プランクの放射式
-
ルンゲクッタ法で数値解析(C...
-
運動エネルギー
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
v^2-v0^2=2ax 今日この式を習っ...
-
d^2r/dt^2の意味
-
電流の時間微分、電圧の時間微分
-
質量流量の記号「・ の読み方を...
-
EXCEL上の数字を自動で振り分け...
-
dx/dt=√(1-x^2)の一般解の求め...
-
Debug.Printで表示される内容を...
-
力学について質問です。 1.棒の...
-
物理で微積をつかう。
-
微分積分のdの意味
-
最後のdv/dtは何でしょうか。
-
加速度 a=dv/dt = (d^2 x) /dt^2
-
運動方程式の微分積分の計算
-
微分記号“d”について
-
Id²θ/dt²=-mghsinθの厳密解の...
-
雨滴の運動質量が変化する落体...
-
運動方程式を求めてください
-
機械力学の問題です!!!
-
d/dx=dt/dx * d/dt =d/dt * dt/...
-
地動加速度が単位インパルスの...
おすすめ情報