新世紀 主導数値解析プ ロ グラ ム の 整備と 公開 基 本解近似解法を 用い た 円 柱に よ る 波の 散乱の 図示 千葉 文浩 牛 島 照夫 電気 通信大学 電気 通信学部 情報工学科 計算科学講座 平成 18 年 03 月 01 日 基 本解近似解法を 用い た 円 柱に よ る 波の 散乱の 図示 1 は じ め に 基 本解近似解法に よ り 障害物に ぶ つ か る 波の 振る 舞い を も と め ， こ れ が 波数が 変わ る と ど う なる か を 見た ． 今回は 円 柱ま わ り の 波の 散乱の 様子を 入射波が 平面波の 場合に つ い て 計算した ． こ こ で は ， 波数パ ラ メ ー タと して 正規 化さ れ た 波数 κ（ 波数 k に 円 柱の 半径 a を 掛け た κ = ka） と して ， κ = 1 か ら κ = 100 ま で を と っ た ． 可視化さ れ た 結果を 見る と ， と く に 波数 κ = 50 あ た り か ら 特徴的な現象が 現わ れ る よ う に 思わ れ る ． こ の 報告書は 情報工学科計算科学講座竹田辰興教授の 支援 の 下に 著者た ち が 実施して い る 研究 の 経過を 紹介す る べ く 作成さ れ た も の で あ る ． 2 問題設定 以 下の ２ 次元の 境界値問題を 考え る ． 今回は こ の 境界値問題の 数値解の 図を 構成す る こ と を 目的と す る ． 次式に お い て ， Γa は 平面 R2 上の 原点を 中心と す る 半径 a の 円 で あ る ． ま た ， k は 波数ベ クトル (l, m) の 長さ で あ り ， Ωe は 円 Γa の 外部領域 で あ る ． −∆ϕ − k 2 ϕ = 0 in Ωe , (F1 ) ∂ϕ = 0 on Γa . ∂n こ の 問題 (F1 ) の 解 ϕ を ϕ = u + ei(lx+my) , k= √ l 2 + m2 と お き ， u は 無限遠で ゾン マ ー フ ェル ドの 外向き 放射条件を 満た す も の と す る ． こ の と き ， u は 次の 非斉次ノ イ マ ン 境界条件付き 帰 着波動問題の 解で あ る ． −∆u − k 2 u = 0 in Ωe , ∂u = ψ on Γa , (Fψ ) ∂n √ ∂u − iku = 0. lim r r→∞ ∂r 上の 式で ２ 番目の 式は u の 法線方向の 導関 数が 境界 Γa 上で 連続関 数 ψ に 等しい こ と を 示す ． こ こ で ， u と ψ は 複素数値関 数で あ り ， ∂/∂n = −∂/∂r は Γa 上の 点か ら 中心に 向か う 法線 方向の 導関 数で あ る ． ３ 番目の 式は ， u が 無限遠で ゾン マ ー フ ェル トの 外向き 放射条件を み た して い る こ と を 示す ． 記 法と して ， r = r(θ) で 平面上の 点を 表す ． こ れ は 複素平面上の 点 reiθ に 対応す る ． こ こ で ， r = |r| で あ り ， |r| は R2 の ユ ー クリ ッドノ ル ム で あ る ． 同様に ， a = a(θ) に aeiθ ， a = |a| を 対応さ せ ， ρ = ρ(θ) に ρeiθ ， ρ = |ρ| を 対応さ せ る ． ψ は 以 下で 与え ら れ る ． ψ=− ∂ i(lx+my) e ∂n 1 on Γa . 以 下で 行う 計算で は k = l, m = 0 と す る ． 注 以 下の 補助問題 (F2 ) を 考え る ． 2 = 0, Z −k Z (F2 ) Z (0) − λZ(0) = 0, Z (−h) = 0. −h < z < 0, 問題 (F1 ) の 解 ϕ は 問題 (F2 ) の 解 Z と 組と なり ， 次の 問題 (HLWW) の 解 Φ = e−itω ϕ(x, y)Z(z) を 構成す る ． こ れ は ， 円 柱の 外部に 無限に 広が る 一 定水深 h の 水域 Ω に お い て ， 入射して く る 平面波が 円 柱に よ っ て 散乱さ れ る と き の 合成波の 振る 舞い を 記 述して い る ． こ こ で ， Γ0 は 水面を 表し， Γ1 は 水面以 外， す なわ ち ， 円 柱の 側面と 水域 の 底を 合わ せ た も の で あ る ． 問題 (F1 ) は ， 言わ ば ， こ れ ら を 真上か ら 見た も の に 相当し， そ の 円 柱の 水平の 断面が Γa で あ る ． −∆Φ = 0 in Ω, ∂Φ (HLWW) Φtt + g ∂n = 0 on Γ0 , ∂Φ = 0 on Γ1 . ∂n こ の Φ = e−itω ϕ(x, y)Z(z) が (HLWW) の 解で あ る た め に は ， λ = k tanh(kh), (ω = gk tanh(kh)) が 成立す る こ と が 必要十分条件で あ る ． 3 基 本解近似解法と そ の 数値計算 波源点と 選点 自然数 N を ひ と つ 選ん で 固定し， θ1 と θj を 次の よ う に 定義 す る ． θ1 = 2π , N θj = jθ1 , j ∈ Z. 正数 ρ と a を ρ < a を みた す よ う に 選ん で 固定し， ρj と aj を 次の よ う に 定義 す る ． ρj = ρ(θj ), aj = a(θj ), 0 ≤ j ≤ N − 1. 点 ρj と aj は ， そ れ ぞ れ 波源点と 選点と 呼ば れ る ． こ の よ う に 波源点と 選点を 等間 隔 ， 同相 に 配置す る こ と を 波源点選点等間 隔 同相配置と 呼ぶ こ と と す る ． 基 本解近似問題 問題 (Fψ ) に 対応す る ， 波源点選点等間 隔 同相配置の 場合の 基 本解近似問 (N) 題 (Fψ ) は 次の よ う に 表現さ れ る ： N −1 (N ) u (r) = qj Gj (r) in Ωe , (N) (Fψ ) j=0 ∂ (N ) u (al ) = ψ(al ), 0 ≤ l ≤ N − 1. ∂n 2 (N) 問題 (Fψ ) に お い て ， 我々 は 以 下の よ う な基 底関 数 Gj (r) を 用い る ． こ れ は ヘ ル ム ホ ル ツ方 (1) 程式の 基 本解か ら 来て い て ， 0 次の 第１種ハ ン ケル 関 数で あ る H0 を 用い て 次の よ う に 表 さ れ る ． (1) Gj (r) = H0 (k|reiθ − ρeiθj |), 0 ≤ j ≤ N − 1. 以 後， 次の 正規 化さ れ た パ ラ メ ー タを 用い る ． こ れ ら は ， 0 < γ < 1 ≤ δ, 0 < κ を みた す ． r ρ γ = , δ = , κ = ka. a a こ れ ら を 用い る と ， Gj (r) は 以 下の よ う に 表現さ れ る ． (1) Gj (r) = H0 (κ|δeiθ − γeiθj |), 0 ≤ j ≤ N − 1. 基 本解近似問題か ら 導か れ る 未知係数を 定め る 連立１次方程式 周期 2π を 持つ 連続関 数 h(θ) を h(θ) = (1) (1) 1 − γ cos θ (1) H (κ|1 − γe−iθ |) |1 − γe−iθ | 1 で 定め る ． H0 (z) = −H1 (z) を 用い る と ， − 1 k ∂ Gj (al ) = h(θl − θj ) ∂r (1) で あ る こ と が わ か る ． こ こ で H1 (z) は １次の 第１種ハ ン ケル 関 数で あ る (N) (Fψ ) に お け る 第１式を 第２ 式に 代入す る と ， こ の 問題に お け る 未知係数 qj に つ い て の 次の 連立１次方程式を 得る ． 1 ψ. k こ こ で H と q と ψ は 以 下の よ う に 定義 さ れ る ． (3.1) Hq = Hjl = h(θj−l ), 0 ≤ j, l ≤ N − 1, q = (qj )0≤j≤N −1 , ψ = (ψ(θj ))0≤j≤N −1 . 連立１次方程式の 解 定義 よ り ， H は 巡回行列で あ る ． よ っ て ， 離散フ ー リ エ変換 に よ り ， 方程式 (3.1) の 解 q = (qj )0≤j≤N −1 を 次の よ う に して 得る ． 1 qj = kN (N ) こ こ で Hn (N ) と Ψn N −1 (N ) Ψn (N ) n=0 Hn eijθn , 0 ≤ j ≤ N − 1, は そ れ ぞ れ 離散フ ー リ エ係数で あ り ， 全て の 整数 n に 対して 次で 与え ら れ る ． Hn(N ) 1 = N N −1 h(θj )e −inθj , Ψn(N ) j=0 1 = N N −1 ψ(aj )e−inθj j=0 for n ∈ Z. (N) こ う して 得ら れ た qj を (Fψ ) の 第１式に 代入す る こ と に よ り ， 円 柱の 外の 点 r で の 関 数の 値を 計算す る こ と が 出来る ． い く つ か の 点 r で 関 数値を 求 め れ ば ， そ れ ら の 点で の 散乱波 の 様子が 分か る ． 3 計算結果の 可視化 4 計算点の 決め 方 計算点を 次の よ う に と る ． 領域 Ωe の 点を 極座標 (r, θ) で 表す ． 原点を 中心 と す る 境界 Γa の 外側に 広が る 同心円 を い く つ か 考え ， そ の j 番目の も の の 半径を rj ， 一 番 外側の も の ， J 番目の も の の 半径を rJ = δJ と す る ． 特に ， 境界 Γa が 半径 r0 = δ0 の 0 番目 の 円 に なる ． ま た ， 境界 Γa の 正規 化さ れ た 半径を δ0 = δ = a/a = 1 と し， 他の 円 も 正規 化 さ れ た 半径で 考え る ． そ れ ぞ れ の 円 周で 角度 θ = 0 か ら 反時計回り に 数え た n 番目の 点の 角 度を θn と す る ． 各円 周毎に 点を N 個と る と す る ． こ れ ら の 同心円 の 間 隔 と 円 周上の 角度の 間 隔 を そ れ ぞ れ 等しく と る と ， j 番目の 円 周上の n 番目の 点の 位 置の 極座標表示は 次の よ う に なる ． (4.1) (rj , θn ) = δ0 + 2π δJ − δ 0 × j, ×n J N for 0 ≤ j ≤ J and 0 ≤ n ≤ N − 1. 以 下で は ， J 番目の 円 に は 計算点を 配置しない が ， 図を 構成す る た め に 用い る 台形の 頂点 が 計算点と 同様に 配置さ れ る ． こ う して ， 計算点に 番号 (j, n) が 付け ら れ る ．（ 図 1 参照．） こ れ に 対応して J 行 N 列の 配列 A を 用意 し， 計算点 (j, n) に 対応す る 配列要素 A(j, n) に 計算 点 (j, n) で の 近似解の 値を セットす る ． 図 1: 計算点の 番号の 付け 方 次に 点の 個数 J, N の 決め 方を 説明す る ． 問題に す る 曲線の 上で は 平面波の １波長相当分 の 長さ あ た り 少なく と も 計算点が ８ 点以 上と れ る こ と を 条件に して 点の 数を 決め る ． 曲線 と して は 上に 決め た 同心円 と 原点か ら 出る 半直線上の も の を 考え る ． 正規 化さ れ た 波数を κ と す る と ， 一 番外側の 円 周で は １波長の 長さ が δJ × κ 個以 上含 ま れ る か ら ， 円 周上で は J は 8 × δJ × κ 個と れ ば 良い ． 半径方向で は １波長の 長さ が (δJ − δ1 ) × κ/(2π) 個以 上含 ま れ る か ら ， 円 周上で は N は 8 × (δJ − δ1 ) × κ/(2π) 個以 上と れ ば 良い ． ま と め る と ， 次の 条件 4 を 満た す よ う に J, N を 選べ ば 良い ． J≥ 8 × (δJ − δ1 ) × κ , 2π N ≥ 8 × δJ × κ. 図の 描き 方 今回報告す る 図は 平面波と 散乱波の 和で あ る 合成波の 実部に 関 す る も の で あ る ． 以 下で は 合成波の 実部の 値を 関 数値と 書く ． 境界 Γa と も っ と も 外側の 半径 δJ 円 の 円 で 囲 ま れ る 領域 に 波の 様子を 図示す る ． 計算点 (j, n) を 一 つ 選び ， こ れ を 起 点と す る 台形の 頂点 A, B, C, D を 反時計回り に 以 下の 式の よ う に 選ぶ ．（ A, B, C, D は 計算点で あ る ．） 表示は 極座 標で あ る ． ま た ， (rj , θn ) 等は (4.1) で 計算さ れ る ．（ 図 2） A = (rj , θn ), B = (rj+1 , θn ), C = (rj+1 , θn+1 ), D = (rj , θn+1 ), 0 ≤ j ≤ J − 1 and 0 ≤ n ≤ N − 1. こ こ で θN = θ0 と した ． そ れ ぞ れ の 台形は 計算点 A で の 関 数の 値に 従っ て 色を つ け る ． こ の 台形を 各計算点で 作る こ と に よ り ， 領域 が 埋め 尽く さ れ ， 波の 図が 得ら れ る ． 図 2: 台形の 頂点の 取り 方 報告書に 収め て い る 図の 色の 決め 方は 次の 通り で あ る ． 各計算点で の 関 数値の 実部が 1/2 よ り 大き い と き は 白， そ れ 以 外の と き は 黒と して そ の 点で の 図形に 色を つ け る ． 上で 決め た 計算点毎に 散乱波の 関 数値を 基 本解近似解法で 計算す る こ と に よ り ， 散乱波 の 様子を 表す デー タが 得ら れ る ． こ の 結果を Mathematica で 読み込み， こ れ と Mathematica で 計算した 平面波と 重ね る と ， 求 め た い 合成波の デー タが 得ら れ る ． 5 作図パ ラ メ ー タと 計算環 境 以 下に 本報告書に 収録した 図の 作図パ ラ メ ー タを 表に して 示す ． ま ず パ ラ メ ー タと して 正 規 化した 波数 κ = ka(a は 円 柱の 断面の 半径) を ， 低 “1 ≤ κ ≤ 20”（ 表 1） ， 中 “1 ≤ κ ≤ 100” （ 表 2） の 2 グル ー プ に 分け た ． 正規 化さ れ た 距離 δ = r/a(r は 原点か ら の 距離) を と り ， κ の 低， 中に 対応して ， そ れ ぞ れ 1 ≤ δ ≤ 10， 1 ≤ δ ≤ 5 の 範囲 で δ を 動か した ． 5 波数 κ 毎に 対応す る 計算点の 数を 表す J と N は 前節の 規 則に 従っ て 与え ら れ る ． 半径方 向の 計算点の 数 J が 不規 則に と ら れ て い る の は ， 他の 用途に 使っ た デー タを 一 部流用して ， 大き い 領域 の も の を 部分表示（ トリ ミ ン グ） した か ら で あ る ． しか し， 計算点の 数 J は 規 則の 条件を 満た す よ う に なっ て い る ． 図の 描き 方で あ る が ， 計算点毎に 台形で 埋め る か た ち で 図を 描い た ． 表 1 は 低い 波数， 表 2 は 中程度の 波数の 場合に 対応して い る ． 色は ， 近似解の 実部が 1/2 よ り 大き い と き は 白， 1/2 以 下の と き は 黒と して 付け た ． なお ， 基 本解近似解法に お い て は γ = 0.9 に 固定した ． 散乱波の 計算は 多倍長計算ラ イ ブ ラ リ MPFR を 用い て 計算桁数 30 桁で 行い ， そ の 結果 を 倍精度に 変換 して ， 作図に 用い た ． 計算に は 計算科学講座の Alpha マ シン (bay02, river, lake) を 用い ， 作図は PC の Mathematica で 行っ た ． 使用した PC の メ モ リ ー は 1GB の 大き さ で あ る ． 表 1: κ = 1, 2, 5, 10, 20; 1 ≤ δ ≤ 10 J N κ = 1 κ = 2 κ = 5 κ = 10 κ = 20 128 128 77 128 230 512 512 1024 1024 2048 表 2: κ = 1, 2, 5, 10, 20, 50, 100; 1 ≤ δ ≤ 5 J N 6 κ = 1 κ = 2 κ = 5 κ = 10 κ = 20 κ = 50 κ = 100 64 64 34 57 103 256 512 512 512 1024 1024 2048 2048 4096 観 察 図化して わ か っ た こ と を 述べ る ． 波数 κ を 大き く して い く と (1) κ = 50 ま た は そ れ 以 上で 円 柱後方の 合成波の 振幅が お だ や か に なる ． (2) κ = 100 ま で は 合成波の 円 柱後方以 外で の 縞は 入射平面波に 対応して 細か く なっ て い く よ う に 見え る ． の ２ 点が 観 察さ れ た ． 6 7 図 今回報告す る 波の 図を 次の 頁以 降に 載せ る ． 低い 波数の 場合に つ い て 3 頁， 中く ら い の 波数に つ い て 4 頁， 計 7 頁で あ る ． 作図パ ラ メ ー タは 節 5 で 述べ た 通り で あ る ． 7 list01B2.nb 1≤k≤20, 1≤d≤10 k=1 k=2 1 list01B2.nb 1≤k≤20, 1≤d≤10 k=5 k=10 2 list01B2.nb 1≤k≤20, 1≤d≤10 k=10 k=20 3 list02B2.nb 1≤k≤100, 1≤d≤5 k=1 k=2 1 list02B2.nb 1≤k≤100, 1≤d≤5 k=5 k=10 2 list02B2.nb 1≤k≤100, 1≤d≤5 k=10 k=20 3 list02B2.nb 1≤k≤100, 1≤d≤5 k=50 k=100 4
© Copyright 2024 ExpyDoc
