2013年9月29日日曜日

空間2次精度


空間2次精度についての話です。

河川の流れであれば、
1次精度風上差分では不十分でも、せいぜい2次精度があれば十分なようです。
ダムブレークのような問題を解く必要はあまりないですから。

2次精度だとベースをオリジナルのLaxWendroff法(2STEPではない)を用い、FluxLimiterを施すことで、大局2次精度を維持できます。

CIP法を用いた非保存系のスキームも使われているようですが、既存のプログラムがある場合、書き換えるのが煩雑です。
また、CIP法だとわずかに数値振動が残りますがTVDだと1次精度となって振動を抑えます。
TVDだとMUSCLも使われているようですが、TVD-LaxWendroffは使われることは少ないようです。
精度の高い乱流モデルを後々組み込むことを考えると、MUSCLという選択もわるくないと思います。

以上の話は、ナビエストークス方程式(or レイノルズ方程式)のうちの移流項に対する解法であり、
運動方程式、連続式をシステム系としてとらえた解法ではありません。なぜなら、生成項となる摩擦項や地形勾配項があるためであり、
C-Propertyを満たすことが難しいためです。
アカデミックの世界でも研究されている方はいるようですが、マニアックなまでに研究されている方はいないようです。
すでに解決された問題だと捉えられているのでしょう。
まあ、分っている人は少ないと思いますが・・・。




流体の数値計算法について(2)

藤井先生の本についての続編です。

一つ面白い記述に

「最近の多くの実用的な計算結果を見る限りでは、一般的には2次精度を維持することが必要といえる」

というのがあります。

乱流の解析では2次精度では不十分な場合がありますが、多くは2次精度程度である程度事足りるとの指摘です。

河川の流れの解析のように、平面2次元流れの解析であれば、この指摘は当てはまると思われます。

平面2次元流れは、どの程度の流れの解像度を得ることが必要かで乱流モデルを選択する必要がありますが、多くの場合は、ブシネスク近似を使って、ゼロ方程式を用いたモデルが多いようです。

①格子形状
格子形状は境界適合座標系を使い、スカラー量とベクトル量を互い違いに配置するスタッガード格子を採用します。

②解析手法
有限体積法を用いることが多いですが、モデルによっては有限差分法が多いです。有限要素法は滅多にありませんが、アカデミックな世界では実施されていることもあります。
また、移流項に対しては、1次風上差分が非常に多いはずです。藤井先生の本を読まれる方は「は?」と思われるかもしれませんが、実用的な解析では、計算が飛ばないことが非常に重要です。
はっきり言って、解析についてあまり知らなくても、計算結果を確実に得ることが求められます。2次精度以上を確保しようとすると、不安定になる可能性が多くなりますので、1次風上が多いのです。
一方、移流項以外は、2次精度の中央差分が多いです。各項の大きさについての議論は別にして、移流項で3次精度以上を確保しようとすると、他の項も同様に3次精度以上を確保する必要があります。
したがって、移流項に対しては大局的空間2次精度程度を確保することが適切のようです。

さて、藤井先生の本に戻りますが、前述の指摘がある程度的を得ています。
分野が違えど、面白い指摘だなと思った次第です。


なお、ブシネスク近似のうちk−εやLESを用いる場合はこの指摘は当てはまりません。
乱流モデルは、小生の間隔では3次精度以上が必要のようです。

2013年9月16日月曜日

書籍


「プログラミング言語AWK」
AWKの本です。最近、USP Labから出版されたものです。
もちろん著者はAWKの開発者であり、AWK版のK&Rと考えてもらえれば良いとおもいます。

ただ、出版といっても、絶版しては再販されている書籍であり、
需要がある書籍と言えましょう。

読んだ感想ですが、古くさいです。
ただ、K&Rよりも実用的であると感じます。
スクリプト言語であり、ワンライナーによるフィルタをかけるだけであれば、序盤を読めばかなり書けるようになると思います。

まあ、敢えてAWKを使う人はあまりいないと思います。今であれば、Perl、Python、Rubyを使えば比較的簡単にできます。
それがテキスト処理だとしてもです。

ただし、プログラミングの歴史を学ぶ点においては、知っておいて損はないかと。
加えて、「プログラミングが、なぜ、今の形になったのか?」ということを理解する一つの助けとなる気がします。

やはり、周辺知識は独立して使うことが出来るものではないですが、メインの知識と組み合わさって威力を発揮するものだと思います。

等間隔QUICKEST


# 久々にpython

for i in range(3,imx-3):#0~imx-1
c = u*dt/dx #
if u >0:
fn[i] = f[i] - c/6*(2*f[i+1]+3*f[i]-6*f[i-1]+f[i-2]) \
   + c*c/2*(f[i+1]-2*f[i]+f[i-1]) \
   - c*c*c/6 *(f[i+1]-3*f[i]+3*f[i-1]-f[i-2])
else:
fn[i] = f[i] - c/6*(2*f[i-1]+3*f[i]-6*f[i+1]+f[i+2]) \
   + c*c/2*(f[i+1]-2*f[i]+f[i-1]) \
   - c*c*c/6 *(f[i-1]-3*f[i]+3*f[i+1]-f[i+2])

等間隔にしか使えないですが・・・。
一応、保存則を使っています。

明日中には不等間隔+TVD化されたQUICKを晒したいです。


2013年9月13日金曜日

書籍

藤井先生著「流体力学の数値計算法」について書いておこうと思います。

双曲型の数値計算法(保存性を意識した差分法)がメインです。基礎方程式は圧縮性流体を想定したオイラー方程式です。解析手法はTVD系です。したがって、非圧縮性流体の解析で使われるスタッガード格子やMAC系の解法に関する手法は述べられてません。

ハッキリ言って、これを読んでも実際に解析が出来るようにはなりません。
実用的な解析をするならもっと簡単な本を読んだ方が良いでしょうし、初心者は読まない方が良いかもしれません。

ただし、味わい深い本だと思います。
小生が購入したのは2004年頃ですが、最近ようやく「あ、これのこと言ってたんだ!」と気づかされることが多いです。FDSやFVSの基本的な考え方が書いてある本はこれくらいだと思います。
浅水流方程式への応用はちゃんと鉛筆動かさないとダメですけどね。
正直、まだ、すべてを理解してません(笑)
ただ、「わかんねぇなー」と思いながら、日々少しずつ理解が深まるところに、楽しさを感じてます。


ところで、

「わかりやすいほにゃらら」
「初心者のためのほにゃらら」
「10日でほにゃらら」

等々。分りやすいのでスイスイ読めますが、長くは読めません(否定している訳ではありません)。

結局、簡単に理解できることは、それだけの価値しかありませんし、他人もすぐ理解できます。
苦しんで覚えたことや、時間をかけて学んだことは忘れませんし、他人もすぐには理解できないことです。

良書を見分ける一つの指標として、「古いけど絶版にならない技術書は良い」というのがあります。
ただし、理解に時間がかかるものが多いです。

技術書に限らず、死ぬまでに多くの良書に出会いたいですね。








2013年9月4日水曜日

蛙跳び2


こんな感じでした。鈍るんです。
LeapFrogは空間積分の右辺をf^{n+1}とf^{n-1}の平均値である0.5*(f^{n+1}とf^{n-1})を使うような良い加減な方法なんであまり好きではないのですが、これ見ると、RungeKuttaを使って計算コスト書けるなら、LeapFrogの方が楽だし、鈍らないし、と判断してしまいますね。何か間違ってることを信じたいです(笑)


2013年9月3日火曜日

蛙跳び


' Leap Frog Method
  dim n as Integer
  dim i as Integer
  dim f(100) as double
  dim fo(100) as double
  dim fn(100) as double
  dim mojiretsu as String

  'Initial Condition
  For i = 1 to ixmx
    if i > 10 And i<20 then
      fo(i) = 1.0
      f(i) = 1.0
      fn(i) = 1.0
    else
      fo(i) = 0.0
      f(i) = 0.0
      fn(i) = 0.0
    end if
  next i
  'Time Marching
  for n =1 to itmx
    'Leap Frog
    for i = 2 to ixmx
      fn(i) = fo(i) - 2*u*dt/dx*(f(i) - f(i-1))
    next i
 
    'Shift Variables
    for i = 2 to ixmx
      f(i) = 0.5*(fn(i)+fo(i))
      fo(i) = fn(i)
    next i
  next n

  'Output
  for i = 1 to ixmx
    mojiretsu = Str(i) +"," + Str(f(i))
    ListBox2.AddRow(mojiretsu)
  next i

4次精度ルンゲクッタと比較すると、2次精度の蛙飛びの方が、鈍りが少ないです。

なぜなんだろう・・・。

2013年9月2日月曜日

ListBox

//ListBoxのクリア

Listbox1.deleteAllRows

Winで作ったバイナリファイルをMacで初めて開きました。
OS依存の機能を使ってないことはありますが、特に問題なく開きました。

ただし、r2でも日本語入力がおかしい箇所がありました。
まあ、日本語使わなければいいんですけどね(笑)

ホビー用途なら、すごくいい環境ですが、
日本語の解説書が充実しないと広まらない気がします。
勿体ない・・・。






2013年9月1日日曜日

値への変換とか(Xojo)

FortranやCを使ってると、値の読み込み時点で型指定をするが、
その他の言語だと、すべて文字列で読み込んで、キャストすることが多いみたいです。

Xojoの型変換は下記のとおり。

Str(数値等):値をstringに変換
Val(文字列):文字列を値に変換

使わないとすぐ忘れます・・・。

ところで、下記に値のやり取りの記録を書きます。

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1.ListBox1を配置して、表示の場合。

//ListBoxの表示をすべて消す。

ListBox1.deleteAllRows

//ListBox1に文字列Program START!を追加
ListBox1.AddRow("Program START!")

2.TextField1に書き込まれた値を変数に代入

//TextField1に書き込まれた文字列を、別途定義した文字列「strings」に代入
//その後、数値型に変換

Dim strings as string
Dim i as Integer
strings = TextField1.Text
i = Val(strings)
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

早くいろいろ晒したいです。