Microsoft PowerPoint - 7

Transcription

1 環境システムモデリング基礎 ( 有限要素法 ~ 発展的知識とまとめ ) 環境システム学専攻愛知正温

2 内容 有限要素法続き 有限要素法拡散方程式の振動解とその回避方法 移流項の処理 多次元への展開 有限差分法 有限要素法 非線形問題への拡張 連成問題への拡張 連立一次方程式の解法の紹介 直接法 反復法 モデリングの全体の流れ

3 拡散方程式を SGFEM で解く場合に現 れる振動解とその対策 拡散方程式の数値解法として 標準ガラーキン法による有限要素法 (Standard Galerkin Finite Element Method: SGFEM) もよく使われており 市販コードもある しかし SGFEM によって拡散方程式を離散化すると 完全陰解法であっても 条件によって非物理的な振動解が発生してしまう ~Vermeer and Verruijt (1981) と Zhu and Yin (1999; 2000), Zhu et al. (2004) の一連の研究成果の整理 ~

4 既存研究の流れ 実際には 拡散方程式と同型の Terzaghi 方程式 ( 一次元圧密理論 ) で研究されてきた 建築物等による荷重 静水= = 0 余( 初圧建築前期剰条間件隙水) 圧発生σ 排水 ε σ つぶれていく体積ひずみ発生定式化 ( 鉱物と水は 構造に比べて十分固いと仮定 ) 有効応力 σ = σ 構成関係 = 質量保存式 = 1 = 境界条件 = 0 at = 0 1D 体積弾性率 ( 拡散方程式型 )

5 既存研究の流れ Vermeer and Verruijt (1981) 非物理的な振動解の存在と 発生条件を明らかにした Zhu and Yin (1999; 2000), Zhu et al. (2004) 非物理的な振動解を避けるための手法を提案した (Zhu et al., 2004)

6 SGFEM のおさらい h h = 0 これをSGFEMで離散化する h h = 0 部分積分 h h h + = 0 ( 弱形式 ) 1 2 i-1 i i+1 n n+1 x=0 x=l 1 形状関数 = と近似すると h + h = 0 これが 任意の数列 h h に対して 恒等的に成立するためには 全てのに対して h + = 個の連立方程式

7 SGFEM のおさらい h h + h = i-1 i i+1 n-1 n x=0 x=l h = h と近似し 完全陰解法を適用すると h h h + h = 0 0 の領域だけ計算すれば良い 1 形状関数 h h h + h = 0 これを計算すると 0 6 h h + 6 h = 6 h + 4h + h この連立方程式を解くことで h を得る

8 振動解発生の原理 6 h h + 6 h = 6 h + 4h + h これは 数列 h に関する隣接三項間漸化式と見なすと = 6 = 3 + = 6 h + 4h + h とおくと h + 2 h + h = ここで = Δ 12Δt =, = + + > 0 かつ + = Δ Δt 2Δt 0 とおくと なので 0 で h = h h 1 h h 1, < 0 だと振動してしまう!

9 振動解が発生しない条件 =, = +, > 0 であるためには = = 6 < 0 つまり 3 + > 0 6 < (Vermeer and Verruijt (1981) が導いた結論 ) ちなみに = 0 では 何が起こるか? h = 1 6 h + 4h + h これは 振動しない まとめると 6 6 > ならば振動しない ならば振動する 結論 : 拡散方程式を SGFEM で離散化すると 拡散係数 タイムステップ 格子刻み幅に制約がある

10 SGFEM 無条件で振動を回避する方法 h h x=0 + h = 0 simplefem (Zhu and Yin, 1999; 2000; Zhu et al., 2004) h h これを計算すると h + + h = i-1 i i+1 n-1 n 1 h h = 形状関数 貯留項の積分をシンプソン則 ( 二次関数近似 ) から台形積分に置き換えていることになる h x=l 強形式の差分近似 h + + h h = h と等価 したがって 台形積分だからと言って 精度的な問題が発生するわけではない

11 無条件で振動を回避する方法 h + + h h = h この隣接三項間漸化式の特性方程式の解は = 2 + ± 2 + > 0 つまり 等比数列部分の公比は常に正 simplefem を使えば 無条件で振動しない ( 当然 差分法についても真 )

12 有限要素法の移流項処理 有限差分法のような簡単な方法が あまり開発されていない 現在 よく使われている方法の紹介 移流項を Lagrangian 計算する Eulerian-Lagrangian 法 (Neuman S.P. (1981) A Eulerian-Lagrangian numerical scheme for the dipersion-convection equation using conjugate space-time grids. J.Comput.Phys., 41, ) の系統 数値拡散項を追加する 多次元の場合 要素サイズの決め方は工夫が必要 Streamline Upwind Petrov Galerkin (SUPG) 法は おそらく一番有名で よく使われているが 解がオーバーシュートするのが難点 Finite Increment Calculus (FIC) 法 (Oñate E. (2004) Possibilities of finite calculus in computational mechanics. Int. J. Numer. Meth. Eng., 60, ) は将来性がありそうだが まだ研究段階

13 多次元への展開 ( 有限差分法 ) 質量保存式 拡散方程式 = = = 0

14 多次元への展開 ( 有限差分法 ) 移流方程式 = 0

15 多次元への展開 ( 有限要素法 ) 質量保存式 弱形式への展開 + + = = 0 ガウスの発散定理 + + = 0 これが 任意の数列 + = 0 に対して 恒等的に成立するためには 全てのに対して + = 0

16 多次元への展開 ( 有限要素法 ) simple FEM + = 0 要素積分 と を計算する必 要がある しかし 任意要素形状では 手計算は大変 座標変換と数値積分を使う

17 座標変換 任意の六面体要素 変換 = Z +1 Y X z y x 積分結果を座標変換して返す -1-1 こちらの座標系で数値積分 = =

18 座標変換 shape function = Z +1 Y 7 = X = iso-parametric approximation = = これで以下の積分は計算できそう = =

19 Gauss 積分 線形内挿要素における 8 点の場所と重み 1,1,1 2,1,1 1,2,1 2,2,1 1,1,2 2,1,2 1,2,2 2,2,2 3 1, 1, 1 +1, 1, 1 1, +1, 1 +1, +1, 1 1, 1, +1 +1, 1, +1 1, +1, +1 +1, +1, 積分を 8 点の値の重み付き和で近似してしまう ( この方法は二次関数の積分まで正確 ) つまり を計算すれば 積分可能 (1 点にする減退積分というテクニックもあるが 素人は手を出さない方が無難 )

20 非線形問題への拡張 比熱や熱伝導率が温度によって変化したりするなど 実際には非線形問題も多い たとえば これを陰解法差分展開すると = = 0 これは線形代数演算では解くことができない

21 Newton-Raphson 法 有限差分法 有限要素法問わず 一般形として = 0 を解くことになる = の成分 = = + として修正し すべての点で <, 準値内に収束するまで繰り返す < のようにあらかじめ決めた基 は どう計算するか?,,,, は 桁落ちで精度が悪いこともある,,,, は 複素数演算で計算負荷がかかるが精度は良い :imaginary number

22 連成問題への拡張 連成問題の例 熱 流体 ( 熱が流体の流れに影響 流体の流れが熱の移動に影響 ) 多孔質体変形 間隙流体 ( 間隙流体の圧力が変形に影響 変形が間隙流体の圧力に影響 ) 半透膜 ( 化学ポテンシャルが流れに影響 流れが化学ポテンシャルに影響 ) etc 相互に密接に関わる問題は 微分方程式も連立して扱う 数値計算上 離散化までは それぞれの式を離散化するだけなので やることは同じ 実際に解く段階において 少しだけ違いが出る

23 連成問題への拡張 エクセルシートで解く場合 1 タイムステップあたり 式の数だけ行 ( 列 ) を割り当てる それぞれの式に対応するように 記入すれば良い 連立一次方程式で解く場合 並べ方はいろいろあるが mesh/node ごとに ひとまとまりにしておいたほうが良い 例

24 連立一次方程式の解法 直接法 (direct method) 連立方程式の解を直接計算する 計算量は 方程式数の 3 乗の比例 丸め誤差を除いて解が正確に計算される ほぼ確実に解に到達する 反復法 (iterative method) 連立方程式の近似解から出発し 誤差から解を少しずつ修正していく 計算量は 方程式数の 2 乗 ~3 乗に比例 解の精度は 修正の打ち切りによっても発生する 係数行列によっては 修正に失敗して解に到達できないことがある 多くの場合 収束させるために 前処理 (preconditioning) を行う必要がある

25 代表的な直接法 LU 分解 下三角行列 L 上三角行列 U を用いて A=LU と変形する おそらく最も使われている解法 ライブラリとしては スイス バーゼル大学開発の pardiso などが有名 ( アカデミックは無料で使用可能 ) Cholesky 分解 ( 対称行列用 ) 下三角行列 L 対角行列 D を用いて A=LDL T と変形する QR 分解 直交行列 Q 上三角行列 R を用いて A=QR と変形する

26 LU 分解 (inner-product form/left-looking) do j=1, n do k=1,j-1 a(k+1:n,j)=a(k+1:n,j)-a(k+1:n,k)*a(k,j) enddo a(j,j)=1/ a(j,j) a(j+1:n,j)=a(j+1:n,j)*a(j,j) enddo do j=1, n x(j+1:n)=x(j+1:n)-a(j+1:n,j)*x(j) enddo do j=n,1,-1 x(j)=x(j)*a(j,j) x(1:j-1)=x(1:j-1)-a(1:j-1,j)*x(j) enddo = と 下三角行列と上三角行列の積にする LU 分解 ( ループ回数 n 3 ) 前進代入 ( ループ回数 n 2 ) 後進代入 ( ループ回数 n 2 ) 実際は ゼロ成分が多いため かけ算するとゼロになるので 計算しなくても良い部分がある 専門的には それを考慮して工夫されたものが実装されていることが多い それらの工夫により 普通の数値解析の場合は ずっと速く終わる

27 QR 分解 (Givens Rotation) : 単位行列に対し =, =, =, = の要素だけ置き換えたもの ただし + = 1 のj 行 i 列目をゼロにするには =, =, = + このような を他の列にも適用 を掛けることを繰り返す ( = を目指す ) ついでに = も計算しておく 終わったら = になっているので 後退代入法で解が求められる と 正規直交行列 Qと上三角行列 Rの積にする r(:,:)=a(:,:) do j=1,n do i=j+1,m d=sqrt(r(j,j)*r(j,j)+r(i,j)*r(i,j)) c=r(j,j)/d; s=r(i,j)/d r(j,j)=d; r(i,j)=0d0 do k=j+1,n tmp=c*r(j,k)+s*r(i,k) r(i,k)=-s*r(j,k)+c*r(i,k) r(j,k)=tmp enddo tmp=c*x(j)+s*x(i) x(i)=-s*x(j)+c*x(i) x(j)=tmp enddo 該当列は手計算

28 代表的な反復法 SOR(Successive Over-Relaxation) 法 LSOR(Linear Successive Over-Relaxation) 法 Krylov 部分空間法 CG (Conjugate Gradient: 共役勾配法 ) BiCG-STAB (Bi-Conjugate Gradient STABilized: 安定化双共役勾配法 ) CG 法の非対称行列版 GMRES (Generalized Minimal RESidual: 一般化最小残差法 ) Krylov 部分空間法の前処理 不完全 LU 分解 ILU や 不完全 Cholesky 分解 (IC) が多い メモリや計算量を省略して 精度が悪い LU 分解や Cholesky 分解を行って 近似解と単位行列に近い係数行列に変換し 収束性能を改善する 省力化 = 不完全さ と 収束性 のバランスが必要だが そこは経験的に決める ( 厳密に求めようとすると その計算に時間がかかって本末転倒になる )

29 SOR 法 = + + と分解する ( は 下三角 対角 上三角行列 ) + + = より = k 回目の反復で得られた解を として = + 1 の計算を繰り返す を上から順に計算していけば 右辺のの成分は 更新済みになるので大丈夫 (acceleration factor) は 0~2の間の数で 経験的に決められる数 どので収束が速いかは あらかじめ知ることは難しい ( 行列を解くより計算負荷!) 正定値対称行列では 必ず収束する それ以外でも比較的広い範囲で収束する w<1は 概して遅くなるが安定化する

30 SOR 法の実際 たとえば 三次元有限差分法の式をメッシュごとに作成し , 1, +, + 1, = とおく これに対して 適当な を仮定してスタート = 1 + { 1 1, 1, +, + 1, }/ 何度も繰り返して < が全ての点でみたされたら 終了する 計算を開始する方向 (x,y,z) を ある程度の頻度で交換したほうが良い エクセルの反復計算は = 1 のときと等価な計算過程である = 1 のとき の対角要素の絶対値がそれ以外の要素の絶対値よりも大きい対角優位行列であれば 収束する

31 LSOR 法 (SOR 法の改良 ) 三次元有限差分法の式をメッシュごとに作成し , 1, +, + 1, = とおく これに対して 適当な を仮定してスタート =, 1, +, + 1, (*) = + + と分解する ( は 下三角 三重対角 上三角行列 ) k 回目の反復で得られた解を として = + 1 とみなして 三重対角行列を Thomas 法で解く 得られた と解く前の を用いて = + 1 を新しい とする (acceleration factor) は 0~2の間の数で 経験的に決められる数 (*) を全てのj,kに対して行い さらに何度も繰り返して < が 全ての点でみたされたら 終了する 三重対角行列と見なす方向 (x,y,z) を ある程度の頻度で交換したほうが良い

32 CG 法 ( 正定値対称行列 ( 移流なし ) 専用 ) 解が 係数行列に関して共役な n 個のベクトルは n 次元空間の基底を構成するので 解はその重み付き和で表せるはずである =, = = = より = = 0 = 初期推定解 からスタートし = = とする = + = = + は に近い方が良いので = + とおく = 0なので = 理論上は n 回以内に良い近似解 が得られるはずだが

33 前処理付き CG 法 係数行列 A を適当にサボりながら LU 分解する (Incomplete LU) たとえば が非ゼロのところだけ計算する ( フィルインなし ILU(0) とも書く ) フィルイン次数を制限する (ILU(n)) 計算してみて ゼロに近い数はゼロと見なして それ以降計算しないなど 計算手順 初期値を決める =, = を解く k=0,1,.( ~0になるまで ) = = = + = = を解く = = + 行列 ベクトル積 ( ループ回数 n 2 ) 前進 後退代入 ( ループ回数 n 2 ) ベクトル ベクトル積 ( ループ回数 n) 理論通りn 回の反復で 全 n 3 オーダーだが前処理がうまくいっていれば 早く終わる さらに 疎行列では 実効的に 行列 ベクトル積 ( ループ回数 mn) ILU の前進後退代入 ( ループ回数 mn) 全体で mn 2 オーダー以下が期待できる

34 BiCGSTAB 法 BiCG 法 = と変形して 無理矢理対称行列にしてから CG 法を適用する BiCGSTAB 法 さらに 補正時に局所的に残差最小化するように修正を加えて安定化させたもの 同様の安定化は CG 法にも適用でき それはCGS 法と呼ばれている 安定化と銘打ってはいるが 実際には 結構破綻するので 前処理は大事

35 GMRES 法 の次の補正ベクトルを探すときに これまで使った全ての補正ベクトルと直交する方向を探す (Arnoldi 直交化過程 ) その方向で 誤差が最小となるような補正ベクトルの大きさを最小二乗法で求める Krylov 部分空間法系では一番安定な部類 とはいえ 前処理した方が速くなってよい 反復するごとに補正ベクトルを記憶していくのと 過去のどの補正ベクトルとも直交する方向を探す必要性から 反復が進むごとに計算負荷が増大

36 モデリングの全体 モデル化すべき現象の選定 どういうプロセスが起こっているか? 記述法の選定 統計モデル モデル関数の選定 物理モデル 支配方程式 エージェントモデル 行動アルゴリズム 解法の選定 厳密解 準厳密解 数値解 問題の性質とモデルに要求される仕様の設計 時間方向 定常 / 周期的定常 / 非定常? 対象期間 必要な時間解像度 初期条件 空間方向 1 次元 /2 次元 /3 次元? 直交 / 極座標 / 球面座標? 対象領域 必要な空間解像度 境界条件 空間や材料の性質 均質 / 不均質 等方性 / 異方性 入力データ整備 データ収集 計測データ 既存資料 情報整理 テーブル GIS 予備的処理 間引き 平滑化 Fourier 解析 etc フィッティング法の選定 Fourier 解析 最小二乗法 遺伝的アルゴリズム etc フィッティング結果の評価 モデルの検定 誤差解析 計算結果の可視化 1D/2D/3D グラフ 動画 GIS etc 結果の利用 現状評価 将来予測 etc

37 数値モデリングの全体 空間離散化手法の選定 有限差分法 有限要素法 境界要素法 etc 時間積分法の選定 陰解法 / 陽解法 / クランク ニコルソン法 ルンゲクッタ法 etc 移流項解法の選定 ( 必要なら ) 風上法 Eulerian-Lagrangian 法 CIP 法 etc 離散化方程式の解法の選定 非線形解法 Newton-Raphson 法 行列解法 直接法 Gauss 消去 LU 分解 etc 反復法 SOR GMRES etc 高速化手法 並列化 領域分割 etc 安定化手法 仮想粘性 etc ソフトウェアの選定 選択肢 自作 ( プログラミング ) 市販 エクセル /R など 専用ソフト 選択のポイント 支配方程式 各種手法 解法 安定性 速度 入力性能 ( プリ ) 扱う問題の大きさ 出力性能 ( ポスト ) ソフトウェアの検証 解析解との比較 etc データ入力 時間 空間の範囲 解像度 初期条件 境界条件 生産 / 消滅項 各種物性 計算実行 計算結果の確認 残差は収束しているか? 物理的直感と合っているか? etc フィッティング 評価 可視化

38 実践上のコツ 最初は 仮想的な簡単な問題で練習して慣れる 最初から大きくて複雑なモデルを作らない 大きくて複雑なモデルは作るのも大変 計算するのも大変で挫折しやすい 単純なモデルですめば それがベスト 単純なモデルで見当を付けてから 徐々に高度なモデルにしていけば良い メッシュ数を増やしすぎない 経験上 パソコンで扱う場合は 5~10 万程度が我慢の限界 最初は色々と試行錯誤する必要があり 速く答えが欲しいので 1 万以下に抑えておくのが得策

39 実践上のコツ 範囲と境界条件に注意する 解析期間と初期条件 解析領域と境界条件 定常解析では Neumann 条件 ( フラックス条件 ) のみでは解が不定になるので 必ず Dirichlet( 指定値 ) 条件を入れる 行列解法 最初は直接法がオススメ ( 待っていれば確実に答えが出るし 収束特性に関わるパラメータの調整などの専門的な知識 技能 勘がいらない ) 非定常解析で安定性が悪い場合は タイムステップを小さくすると 対角優位行列になって 解が得やすくなる 非線形問題では 変化量が小さくなると線形問題に近くなり 安定する効果もある モデルと観測データが合わない場合は まずモデルから疑う