So-net無料ブログ作成
前の10件 | -

今日のアニソン、「キノの旅」から『Bird』 [今日のアニソン]

今日のアニソンは、劇場版「キノの旅」から『Bird』です。


OP曲も紹介するケロよ。



nice!(1)  コメント(0) 
共通テーマ:音楽

バーガース方程式の数値解法6 中心差分 [数値解析]

バーガース方程式の数値解法6 中心差分

 

バーガス方程式

  

の対流項を

  

と中心差分を用いて離散化すると、(1)は

  

となる。

(2)の対流項を反復計算1つ前のを用いて線形化すると、(2)式は

  

となり、3重対角係数の連立1次方程式に変換でき、TDMAThomasのアルゴリズム)を用いて解くことができる。

 

また、

  

と、クーラン数Cと格子レイノルズ数Rを定義すると、(2)式は

  

 

となる。

C>0、格子レイノルズ数R>2のとき

  

C<0、格子レイノルズ数R<−2のとき、

  

となり、(4)式の右辺の項の係数が負となり、数値的に不安定になる恐れがある。

したがって、中心差分を用いるとき、

  

 

を満たすように空間刻み、分割Δxをとることが望ましい。

 

この計算スキームを用いて、

を、ν=0.01Δx=0.05Δt=0.025で数値的に解くと、次のような振動解が得られる。

 

 

 

数値解が振動するのは、0≦x≦0.5では、u=0.5の前後。このときの格子レイノルズ数Rを計算すると、

になっている!!

 

格子レイノルズ数|R|>2のときに必ず振動解が現れるわけではないけれど、計算条件によっては、このように振動する数値解が得られる場合があるので、対流項に中心差分を用いる時、|R<2となるようにΔxを十分に小さく取る必要があることが理解できるだろう。

 

Δx=0.025にとると、数値解が振動することはないけれど、実線で示される厳密解(もどき)と比較すると、振動の傾向は残っており、x=0.5近傍での誤差が大きくなっている。

 

 

計算に使用したプログラムは次の通り。

 

 

parameter (nx=40, nt=40)
real u(0:nx,0:nt)
real a(nx),b(nx),c(nx),d(nx)
real nyu

pi=acos(-1.0)
a=0. ; b=0. ; c=0.; d=0
u=0.

nyu=0.01
dx=1./nx; dt =1./nt
dx2=dx*dx

! initial condition
do i=0,nx
    u(i,0)=sin(2*pi*i*dx)
end do
!u(nx,0)=0.

! boudary condition

! caliculate coefficients
do j=1,nt
    do iter=1,100
        do i=1, nx-1
            u(i,j)=u(i,j-1)
        end do
        do i=1,nx-1
            c0=u(i,j)
            aw=nyu/dx2+c0/(2*dx)
            ae=nyu/dx2-c0/(2*dx)
            a(i) = -aw
            b(i)=aw+ae+1/dt
            c(i)=-ae
            d(i)=u(i,j-1)/dt
        end do
        d(1)=d(1)-a(1)*u(0,j)
        d(nx-1)=d(nx-1)-c(nx-1)*u(nx,j)
        call tdma(a,b,c,d,nx-1)
        do i=1,nx-1
            u(i,j)=d(i)
        end do
    end do
end do

do j=0,nt,4
    write(*,100) j*dt, (u(i,j),i=0,nx,4)
end do

write(*,*)

100 format(12(f8.5,1x))
end

!  TDMA
subroutine tdma(a,b,c,d,n)
real a(n), b(n), c(n),d(n)

! 前進消去
do i=1,n-1
    ratio=a(i+1)/b(i)
    b(i+1)=b(i+1)-ratio*c(i)
    d(i+1)=d(i+1)-ratio*d(i)
end do

! 後退代入
d(n)=d(n)/b(n)
do i=n-1,1,-1
    d(i)=(d(i)-c(i)*d(i+1))/b(i)
end do

end

 

 


nice!(0)  コメント(0) 

数値計算流体力学のノウハウを学べるのはここだけだと思うにゃ [ひとこと言わねば]

非線形のバーガース方程式の数値的解法やその計算プログラムを、惜しむことなく、ブログに書く。このようなブログは、世界中を探しても、この「ねこ騙し数学」以外にないと思うにゃ。嘘だと思うならば、Burgers Equationで検索をかけてみるといいにゃ。英語で書かれたpdfファイルなどは幾つか存在するけれど、プログラムまでは示されていないにゃ。したがって、このように断言してもいいうさ

数値(計算)流体力学のノウハウを学べるのは、「ねこ騙し数学」だけだと。


しかも、数値流体力学の素人がこういうことをやっているところが味噌だにゃ。


と、例によって、自画自賛する(^^ゞ

誰も褒めてくれないし、たまに、こういった自画自賛をしないと、とてもじゃないけれど、モチベーションを保てない。

これは数値流体力学に限ったことではないけれど、教科書的な知識は意外に役に立たないもので、教科書的な知識だけじゃ〜、数値計算のプログラムは作れない。どうしても、教科書などに書かれていないノウハウや計算テクニックが必要になるものなんだにゃ。そして、こうしたノウハウや計算上のテクニックは、研究者や大学の先生、実際に商用のプログラム開発にあたっているヒトは、飯の種だから、簡単には教えてくれない、本のどこにも書いていない。本当に知りたいことは、本には書かれていないものなんだにゃ。
なのに、ネムネコは、自身の経験で体得した、こうしたあノウハウや計算上のテクニックなどを芸として、タダで披露している。
なんて、ネムネコは寛大で心優しいことか!!
そう、思わないケロか?



nice!(1)  コメント(0) 

今日のアニソン、「ハガレン」から『Period』 [今日のアニソン]

解散したChemistryが再結成し、活動を開始したらしいので、今日のアニソンは『Period』だにゃ。


「ハガレン」の主題歌といえば、やっぱり、この曲↓が一番有名だろうから、こっちも記事に埋め込んでおくにゃ。



nice!(0)  コメント(0) 
共通テーマ:音楽

バーガース方程式の解法5 陰解法 [数値解析]

バーガース方程式の解法5 陰解法

 

線形化されたバーガース方程式ではなく、(純)陰解法を用いた非線形のバーガース方程式

  

の解法について考える。

(1)式の時間微分には前進差分、対流項に風上差分を用いると

  

となるので、(1)式は

  

と離散化することができる。

 

(2)式を解けば、を得ることができるが、線形化されたバーガース方程式

  

の時とは異なり、(2)式は連立2次方程式なので、何らかの方法で(2)を線形化し連立1次方程式に書きなおす必要がある。

そこで、(2)式を

  

と線形化する。

ここで、は反復計算の一つ前のの推測値であり、具体的には反復計算の一つ前のである。

適当な推測値を与え、その推測値に基づき連立1次方程式を何度も繰り返してゆけば、いずれ

  

に収束してくれるであろう。

 ――(3)は優対角の連立方程式だから収束する(はずだ)――

この見込み、予想にしたがって、(2)の代わりに(3)を解き、

  

が所定の精度ε以内に収まったら、それを(2)の解とみなす。

 

(3)式は

  

とおくと、

   

となる。

したがって、

  

となり、TDMA(トマスのアルゴリズム)を用いて解くことができる。

 

この考えにしたがって、

  

を、ν=0.01Δx=0.01Δt=0.01を解いた結果は次の通り。

 

 

このグラフを見ると、x=0.5の近傍でuが大きく変化していることが分かる。したがって、精度よく計算するためには、x=0.5近傍に計算点(計算に使用する空間的な点)を複数個配置する必要がある。(下図参照)
――「x=0.5の近傍における赤い曲線と黒い曲線の不一致の程度は打切誤差の範囲を大きく越えている!!」ことに注意――

 

数値計算では、変化が急なところに計算点を複数個配さないと、このように予期せぬ大きな誤差(?)を生むことがある。

 

計算に使用したプログラムは次の通り。

 

parameter (nx=100, nt=100)
real u(0:nx,0:nt)
real a(nx),b(nx),c(nx),d(nx)
real nyu

pi=acos(-1.0)
a=0. ; b=0. ; c=0.; d=0
u=0.

nyu=0.01
dx=1./nx; dt =1./nt
dx2=dx*dx

! initial condition
do i=0,nx
    u(i,0)=sin(2*pi*i*dx)
end do
!u(nx,0)=0.

! boudary condition

! caliculate coefficients
do j=1,nt
    do iter=1,100
        do i=1, nx-1
            u(i,j)=u(i,j-1)
        end do
        do i=1,nx-1
            c0=u(i,j)
            aw=nyu/dx2+(c0+abs(c0))/(2*dx)
            ae=nyu/dx2-(c0-abs(-c0))/(2*dx)
            a(i) = -aw
            b(i)=aw+ae+1/dt
            c(i)=-ae
            d(i)=u(i,j-1)/dt
        end do
        d(1)=d(1)-a(1)*u(0,j)
        d(nx-1)=d(nx-1)-c(nx-1)*u(nx,j)
        call tdma(a,b,c,d,nx-1)
        do i=1,nx-1
            u(i,j)=d(i)
        end do
    end do
end do

do j=0,nt,10
    write(6,100) j*dt, (u(i,j),i=0,nx,10)
end do

write(*,*)

open (1,file='Burgers1.dat',status='replace')
do i=0,nx
    write(1,110) i*dx, (u(i,j), j=0,nt,10)
end do
close(1)

100 format(12(f8.5,1x))
110 format(11(f8.5,1x),f8.5)

end

!  TDMA
subroutine tdma(a,b,c,d,n)
real a(n), b(n), c(n),d(n)

! 前進消去
do i=1,n-1
    ratio=a(i+1)/b(i)
    b(i+1)=b(i+1)-ratio*c(i)
    d(i+1)=d(i+1)-ratio*d(i)
end do

! 後退代入
d(n)=d(n)/b(n)
do i=n-1,1,-1
    d(i)=(d(i)-c(i)*d(i+1))/b(i)
end do

end

 

なお、このプログラムでは収束の判定をせずに、(3)を100回反復計算させた値を収束値とみなしている。反復回数を1000回にしたときの計算結果とほとんど一致していることは確認済みである。

 

 


nice!(0)  コメント(0) 

今日のアニソン、「うえきの法則」から『Earthship〜宇宙船地球号』 [今日のアニソン]

今日のアニソンは、アニメ「うえきの法則」から『Earthship〜宇宙船地球号』です。


さらに、この曲を♪



nice!(1)  コメント(0) 
共通テーマ:音楽

[放物型偏微分方程式にシンプレクティック法は使えないか?(有限要素法の具体化)] その2 [数値解析]

 ところで(6)左辺1項目の∂p/∂tって何でしょう?。pは一般化運動量でしたが、実質は、

 

 

 

でした。だったら、

nsymp42-001.png

nsymp42-002.png

の事だ!、くらいでいいじゃないですか(^^)(5)(15)を使えば、

nsymp42-003.png 

 

です。陽に書いてある対称行列をM0とすると、

nsymp42-004.png

となります。M=t0Nは構造解析では要素質量マトリックスと言われ対称行列で、これも正定値。

 (6)右辺の熱源密度の積分も、なんか似たように感じでできるのは明らかですよね?。なのでその結果をλで表します(4行縦)。以上まとめれば、一要素に対して、

 

 

 

が得られます。ちなみにこのタイプの質量マトリックスは、Consistent Mass(整合質量マトリックス)と言われます。

 

 後は頑張って、KとMの具体的形を定めるだけです。・・・頑張りますか(^^;)

nsymp42-005.png
nsymp42-006.png

 で、Mは・・・。M0の密度の濃さを見ると、とても計算する気にならない。ここは手を抜こう(^^;)。じつは運動方程式(16)は、さっき言ったように構造解析における梁要素の運動方程式とほとんど同じなんです(Kの前の符号が違うだけ)。梁要素とは、回転バネを持つ線材の事です。к2の符号を変えたものが回転バネのバネ定数に相当します。従って質量マトリックスMの出し方は全く同じ。そのアナロジーで行けば、節点未知数X=(T1T2w1w2)T1T2は、図-3の両端点の鉛直変位,w1w2は同回転角変位を表します。

 ところで(16)の加速度は結局、変位に比例する訳ですから、線形バネでつながれた多質点系とみなす事も可能です。このような場合、要素上の質量を節点へ集めて集中分担させても、経験的に実用上十分な結果が出ます。要素に沿った質量の線密度は1と考えられるので、鉛直変位T1T2については両端の節点に要素全質量1×LLの半分、L/2を分担させればいいでしょう。回転変位w1w2については、質量に対応するのはいわゆる回転慣性なので、節点を中心とした要素長の半分までの慣性モーメントを持たせます。

nsymp42-007.png 

 従って質量マトリックスは、

nsymp42-008.png

 

でもたぶんOK。このタイプの質量マトリックスは、Lamped Mass(離散質量マトリックス)と言われます。エッ、どうして(18)1/2L2/24の入る位置がわかったかって?。(17)の行の次元で判断してます(^^)

 

 残るは要素マトリックスを重ね合わせて、全体系を作る事です。図-4のように解析領域[04]n等分割し、左から要素番号をk12,・・・,nで与えます。同図から明らかにように、要素kの節点は必ずkk+1です。

 各節点jの未知数を並べたXj(Tjwj)tをさらに全部並べて作った未知ベクトルXを、

 

iyada-100.png

 

で定義します。Xの成分数はj1n+1なので、2×(n+1)2n2です。Xの成分番号i,要素番号k,節点番号jの関係は明らかに、k → jkk+1 → Tji2k1wji2kです。

 全体系の係数マトリックスKは、1要素上の(6)の積分結果を全ての要素について足して行ったものなので、要素マトリックスKkの成分を、Kのどこに足し込めば良いかは次の要素・節点対応表があればわかります。

 

 例えば要素マトリックスKkの列番号1234の成分は、上記表より、全体系Kの2k-12k+12k2k+2列に足し込めば良い事がわかります。

 行番号は?。要素系と全体系の行番号の対応は、前回やったように、式(6)右辺にあるはずの境界項が隣の要素と足して消えるという条件から決まります。それは、隣接要素で共有される節点で(6)右辺を足したら境界項がキャンセルされるでした。式(6)左辺は今までやってきたように、離散化すればXの係数マトリックスの1行となるような横ベクトルです。共有節点で右辺を足すという事は、その横ベクトルを同じ節点で足すという事です。実際にさっきの要素・節点対応表で共有節点を考慮すると、列番号と全く同じ対応になるのがわかります。こうしてKも正定値対称行列になります。Mについても同様です。

 

 

 最後に、けっこう無茶な定式化をしてるので、やっぱりお釣りが出ちゃったって話です。図-4の要素1で式(6)を、境界項も無視しないで書きます。

 

 

 ただし要素2と足して消える項と熱源密度は省略してます。(19)u(x)f1(x)にすると、(17)(18)を参照しつつ、

 

 

 

 

になります。d3T1/dx3の値は境界条件に書いてないので、解き終わるまでわかりません。しかしこれをまるごと、境界条件T(0

t)0から、

 

 

におきかえる事は可能です。具体的には初期条件でT1P10を与え、後は強制的にdP1/dt0を与え続ければOKです。

 u(x)f2(x)f4(x)(19)右辺は0なので問題なしです。u(x)f3(x)の場合、

nsymp42-010.png

になります。d2T1/dx2の値も境界条件に書いてませんが、その値はわからないのか?と問われれば、わかります。d2T1/dx2は地点x0の熱収支を表します。境界条件はx0で等温だったので、任意の時刻td2T1/dx20です(一般には、熱源密度も考慮する)。よって、

nsymp42-011.png 

を追加境界条件にできます。同様に要素nでは、

nsymp42-012.png


 

を境界条件にできます。

 

 さすがにExcelで以上の手順を追うのは苦しいので、プログラムを作ります。境界要素法プログラムのテストはどうしたんだ?って話もありますが・・・(^^;)

 

 


nice!(0)  コメント(0) 

今日のアニソン、「妹さえいればいい」から『明日の君さえいればいい』 [今日のアニソン]

今日のアニソンは、アニメ「妹さえいればいい」から『明日の君さえいればいい』です。


さらに、同アニメのED曲を♪



nice!(2)  コメント(0) 
共通テーマ:音楽

放物型偏微分方程式にシンプレクティック法は使えないか?(有限要素法の具体化)] その1 [数値解析]

放物型偏微分方程式にシンプレクティック法は使えないか?(有限要素法の具体化)]

 

放物型偏微分方程式(熱伝導方程式)のためのシンプレクティック法の有限要素法を具体的に書きます。最初に形状関数fj(x)を決めます。

 j14に対する形状関数の「形状」は、以下です。

 

 

nsymp41-001.png

 係数ap(j)p03j14は、

nsymp41-007.png
 

 

nsymp41-fig-002.png

 

で決定できます。ejj番目の成分が1の単位ベクトルです。形状関数fj(x)を使うと、区間[0L]上の関数T(x)を、x0Lでの関数値T1T21階微分値w1w23次補間できます。


 じつはこの3次補間のやり方は、構造解析分野のフレーム解析で梁要素という名で頻繁に用いられてまして、昔むかしに面倒くさいよぉ~と、ぼやきながら渋々行った計算結果がノートに残っていたので、それを使います(^^;)

nsymp41-003.png 

 

 ただしξx/Lです。x0→Lに対してξ0→1。こうすると後の計算が楽なんですよ。(4)右辺のマトリックスの列ベクトルをNjj14で表します。そうすると(3)(4)を比較し、

 

 

が得られます。Njは自然に縦ベクトル。また(T1T2w1w2)tをXで表す事にします。

 前回の結果から、

 

 

 

 

を計算すればOKでした。(6)右辺の熱源密度は具体的に与えていないので後で考える事にし、(6)左辺2項目をまず処理します。

 

 

に注意し、(4)から、

 

 

 

 

同様に(5)から、

 

 

 

 

になります。従って(6)左辺2項目の被積分関数は、

 

 

 

 

ですが、(10)右辺の最初の2つのベクトルの積はスカラーです(Njは縦ベクトル)。スカラーなら転置をとっても同じなので、

 

 

 

 

 

 

 

 そうすると全体で、4列横×4行縦×4列横×(4×4行列)×4行縦となり、全体で一個の行列積を定義できます。行列の結合法則より順序を変えない限り、どこから計算開始してもOKです。という訳で、

nsymp41-004.png

なので、

nsymp41-005.png

です。陽に書いてある対称行列をB,N=(N1N2N3N4)とすると、

 

 

 

 

 (13)右辺の()に注目します。Njtは横ベクトルです。BNは4×4マトリックスの積。よって()全体で4列横の横ベクトルになります。それをj14に従って縦に並べた姿を想像すると、

 

 

 

 

ですよね?(^^)。NtBNは対称行列です。Bは対称だったから、(tBN)t=NttN=NtBN。しかも作り方から正定値対称です。K=1/L3×tBNとします。構造解析の用語で言うと、これは要素剛性マトリックスです。しかし何故Kは正定値対称になったんでしょう?。ここまでの一連の流れは、重み付き残差法ならどんな方法でもいっしょです。つまり最後の結果は常に、K=1/L3×tBNになります。なので特別なのは、Bが(12)で与えられて正定値対称行列になる点です。どうしてこうなったかというと、重み関数とテスト関数に同じものを使ったからです。ガラーキン法だったからです。そして重み関数とテスト関数の微分階数を同じくする、弱形式に持ち込んだからです。

 数値計算において、支配方程式の係数マトリックスが正定値対称になってくれるのは、とてもとても価値があります。正定値対称行列はものすごく数値的性質が良いからです。そのような理由から、たとえ変分原理の存在しない系であっても、ガラーキン法の弱形式は好んで使われます。

 

 (6)左辺1項目を処理します。今まで省略してましたが、

 nsymp41-006.png

って書きますね。

エェと、X(t)(T1(t)T2(t)w1(t)w2(t))tではなくて、正しくはこれも、X(xt)(T1(xt)T2(xt)w1(xt)w2(xt))tでないの?、という意見はあると思うんですが、T(xt)が要素内任意点の解の値だったのに対し、Tjwjなどの場所は要素端に固定されてるので、これぐらいは省略させて下さいよ(^^)

 


nice!(0)  コメント(0) 

ねこ騙し数学、ブログテーマ学問で4位に!! [ひとこと言わねば]

ねこ騙し数学、ブログテーマ学問で4位に!!
ウソじゃないケロ、本当に4位になったにゃ。これがその証拠にゃ。


まずは、この喜びを踊りで表現するにゃ。


過去にも4位になったことはあるんだけれど、「ねこ騙し数学」の上にあるブログは


は、どれも「歴史」のブログで、一般に人気のない「数学」のブログじゃ〜「歴史」のブログに勝てないにゃ。


この3つのブログに勝つ方法は、何かないケロか?



nice!(2)  コメント(0) 
前の10件 | -