コンピュータ・シュミレーション － Ｅｘｃｅｌを使って差分法で偏微分方程式を解く。
 高校物理の波動分野では定常波の例として必ず気柱の共鳴の話が登場します。そして、閉端での固定端反射の説明まではよいのですが、開口端では開口端補正だけ離れたところで自由端反射をするのはなぜかと質問されても、音波が開口部から広い空間に広がるときの波の乱れで一部反射するからだとしか答えられません。
専門書を読んでも、インピーダンス概念を用いたよく分からない説明が載っているだけです。本当のところは、偏微分方程式を解くしかないのだろうと思いますが、境界条件によっても変わるでしょうから厳密な数式解などは求めようがないのだろうと思います。
そこでコンピュータの出番です。以前に偏微分方程式を差分法で解を求めようとしたことがあるのですが、当時のパソコンのメモリー量の制限、境界条件の設定の煩雑さや、結果をグラフィカルに表示できるように設計する手間やら、おまけに粗く近似した所為で解が発散する結果になったため諦めた経緯があります。
今回はEｘｃｅｌを使って、殆どプログラムなしで、単にセルに式を入れるだけで自動計算させることにしました。方法はすべて自己流で工夫してみました。
その結果、Eｘｃｅｌは２次元の偏微分方程式を解くのに大変便利なツールとして利用できると、実感した次第です。
 使ったEｘｃｅｌファイルをｚｉｐ形式に圧縮してありますので、解凍して開口端反射の様子をご覧ください。（見るにはEｘｃｅｌのマクロを有効にする必要があります。）
左図のように長さ100、幅30の矩形枠（閉管に相当する）の内部にパルス波を発生させ、管口へ向けて進ませます。その後の波の様子をシュミレートするという設定です。図の点線の境界は透過壁で反射が起りにくい条件にしています。
シュミレーションからは、開口端での反射が自由端反射であること、何度も反射波が往復する内に定常波が現れてくることや、定常波状態で外部に音が出るタイミング等が分かります。
（これは、あくまで２次元でのシュミレーションであり、この差分近似でどの程度正しいか分かりません。 ですが結構リアルです。 ある程度これに似た現象が起こっていると思われます。）
波動方程式は２階偏微分方程式なので、過去と現在の時間変化と現在の空間方向の変位のゆがみから現在から未来への変位の時間変化が予測できるわけです。
Eｘｃｅｌシート上に150行×200列の領域を４領域作って、それぞれ、シートの上から作業領域、未来領域、現在領域、過去領域を表しています。
シートの未来領域で濃い青色のついたセルが閉管の壁を表しています。未来領域のみセルに式が入っていますが、セルの色が違うとセル内の式も異なります。
シートの最上部あたりにあるＢＶＡのボタン（ｃｌｏｃｋ）をクリックすると横のステップ数の回数だけ計算され時間が進みます。このステップ数を変えることで１クリックで進むステップ数を変更することができます。
・結果のグラフを見る。
Eｘｃｅｌで偏微分方程式を解く方法の参考として以下の説明をお読み下さい。
ＥＸＣＥＬを使って差分法で偏微分方程式を解く
（波動方程式と境界条件－音波のシュミレーション） ｂｙ Ｗａｔｓｏｎ
はじめに
コンピュータ・シュミレーションの醍醐味は、数学的には単純に解けない問題を、有限回の単純計算の繰り返しで近似的に解くことが出来ることである。物理現象に限らず、社会現象を含めて「数学的モデル」の殆どは微分方程式（特に偏微分方程式）の形で表現される場合が多いが、これらは単純に数学的には解けず、コンピュータの数値計算に頼らざるを得ない。
物理教育においても波動や電界など、現象の基本が偏微分方程式で表される分野が多く存在する。この分野の教育においてコンピュータの利用がされてはいるが、殆どがコンピュータ・アニメーションの域を出ないもので、本格的なシュミレーションとは言えないものが多いように思われます。
この原因の一つは、プログラムで偏微分方程式を解く際の境界条件の設定の煩雑さや、結果のグラフィカルな表示に多大な労力を要する点にあります。この点、ＥＸＣＥＬを使えば、これらの難点を容易にカバーすることができます。境界条件はシート上に描いた境界セルに条件を考慮した式をコピーして貼り付けるだけで行えるし、またＥＸＣＥＬの優れたグラフ機能を利用して結果を労なく３Ｄグラフとして視覚化できるなど、すばらしい点が多い。
先ず、音の波動方程式の説明をします。
差分法による２次元波動方程式（双曲型偏微分方程式）の近似解法
２次元波動方程式は一般に
と書ける。ここで、ｃは波の伝搬速度である。
音波の場合にはｕは媒質（空気）の変位ではなく、密度あるいは圧力の変位を表すと考えるのが適当である。また音速ｃは
で与えられ、ここでＫは空気の体積弾性率、ρは密度である。
また空気の圧縮は断熱変化であるから、 ＝一定より Ｋ＝ となり、
 とも表せる。ともかく、ｃ＝定数 を仮定して上の２次元波動方程式を解けば良い。
・差分の格子
時間、空間を格子に分割して  ＝ 
で表す。ここでｋは時間ｔ方向、ｉはｘ方向、ｊはｙ方向の格子番号を表す。
・時間２階微分
は格子間隔をΔｔとして
で与えられる。（右図参照）
・空間の２階微分の和（ラプラシアン）
は格子間隔をｈとして
で与えられる。（右図参照）
ゆえに、波動差分方程式は
となる。ここで、簡単のため  ＝1/4と置く。すなわち、１格子時間に1/2格子距離進む波を仮定する。
（ｃΔt/h＝１と採っても良いように思うが、差分シュミレーション結果は不安定となった。
ｃΔt/h＝１／２だと安定する。）
この式から
 （＊）
となる。つまり、未来の状態が現在の空間変形状態と過去の状態により決まる。
ＥＸＣＥＬによる実現法
 (1)前処理
(a)ワークシート上に３領域をとり、それぞれ、過去、現在、未来を表す。
(b)過去、現在の領域のセルに初期値を設定する。
(c)未来の各セルに（＊）式に対応する計算式を、境界領域のセルには後述の境界条件の式を設定する。
(2)繰り返し
(a)計算結果（未来領域の式は現在領域と過去領域のセル値を参照して計算する式になっています。）として得られた未来の値を現在領域に、現在の値を過去の領域に順次移動する。（具体的には作業領域を使って、図の①、②、③のコピー操作［値のみ］を行う）
(b)これを繰り返すことで、時間を１ステップづつ進行させる。
境界条件について
 ○自由領域のセル
境界でない各セルには（＊）の計算式を与える。
○媒質変位で固定端（＝密度変位では自由端に相当する）のセル
壁を背にする向きで考えると、前面のセルから受ける圧力と等しい反作用を壁から受けることになるので、前面のセルと等しい状態のセルが背面に存在すると仮定する。
○自由端でも固定端でもない単に領域の端に位置するセル
 これは非常に難しい問題である。シュミレーションでは領域は有限なため、端が必ず存在する。これを仮に自由端あるいは固定端にするとそこで波の反射が生じる。反射せずに波は通り抜けて欲しい。波が通り抜けている状態を想像すると、通り抜けた先のセルには端のセルの過去の状態が反映されているはずである。
そこで、欠けたセルには端のセルの過去の値を入れてみた。このシュミレーションでは波が１格子距離進むのは２格子時間であるから、対応する過去の格子はすでにないので、残っている１格子時間前の値で代用した。実際にこの境界条件で試してみても、殆ど反射なしに波が境界を透過していくことが確かめられた。
気柱の固有振動（閉管）のシュミレーション
閉菅での気柱の振動は、閉じた端で固定端反射するのは納得できるが、開口部での自由端反射はなかなか理解し難い。そこで、開口部ではどのような現象が起こっているのか、シュミレーションで調べてみた。
(1) "Sheet1"上に150行×200列の領域を４領域作る。それぞれ作業領域、未来領域、現在領域、過去領域と呼ぶ。
(2) 各領域に閉菅のパターンをセルの色付けで描く。
(3) ガウス分布に似た１パルス波形を数式で作成して１行を作る。この行の値を現在領域の閉管内部のセルにコピーする。
(4) 同じパルス波形を１／２セル左にずらした波形を数式で作成して別の１行を作る。この行の値を過去領域の閉管内部にコピーする。これで初期条件の作成が完了。
（これで動かせば、波は右へ１／２格子距離ずつ進行するはず！）
・パルス波は開口部（右）へ向かって順調に進む。（図ａ）
・開口部へ到達後パルス波が円形に広がるとともに波の後面が負へ変位し始める。（変位は音波の密度を表すことを考えると、密（＝正）のパルス波から疎（＝負）が生まれたことを意味する）（図ｂ）
 
・この負の波は左へ向かって進む。（密の反射波が疎なので”自由端反射”だ！）この波はかなり不規則な形をしている。（図ｃ）
・左の閉端に到着した波の一部が大きく負へ変位する。（疎が疎で反射＝”固定端反射”）（図ｄ）
 
・これが戻り、開口部で再び反射することを繰り返す内に次第に波形が緩やかになっていく。定常波の出現だ！（図ｅ）（図ｆ）
 
・ 変位が０をまたぐ付近で管の外部への波の送り出しが見られる。（図ｇ）
ＶＢＡマクロ
（プログラムはたったこれだけ）
Dim rg1 As Range
Dim rg2 As Range
Dim rg3 As Range
Dim rg4 As Range
Dim t As Integer
Private Sub CommandButton1_Click()
Dim i As Integer
Dim c As Integer
Set rg1 = Worksheets("Sheet1").Range(Cells(51,
1), Cells(200, 200))
Set rg2 = Worksheets("Sheet1").Range(Cells(201,
1), Cells(350, 200))
Set rg3 = Worksheets("Sheet1").Range(Cells(351,
1), Cells(500, 200))
Set rg4 = Worksheets("Sheet1").Range(Cells(501,
1), Cells(650, 200))
c = Cells(7, 8).Value
If c = 0 Then
c = 10
End If
t = Cells(7, 11).Value
For i = 1 To c
rg2.Copy
rg1.PasteSpecial Paste:=xlPasteValues
rg3.Copy
rg4.PasteSpecial
rg1.Copy
rg3.PasteSpecial
Next i
t = t + c
Cells(7, 11).Value = t
Worksheets("sheet1").Cells(1,
1).Activate
End Sub
このページのトップへ
|
ＰＤＦファイル
このＰＤＦファイル
物理教育
・物理教育序論
・慣性力について（慣性力序説）
・慣性力はみかけの力ではない（慣性力続論）
慣性力についての対話(1)
ＮＥＷ！
慣性力とエネルギー保存則（慣性力対話３）
Ｎｅｗ！
・エネルギーの発見と 筋肉
・力はベクトルか？（力の正体をさぐる。）
ＮＥＷ！
力学系としての電磁気学、Ｍ．プランクの「理論電磁気学」を読む
・物理”診断テスト”と”運動についての常識概念”
物理実験
・新方式の「運動法則の実験」
e-mail
ご感想等
|
