分子運動に着目した音波伝搬の数値解析
Numerical Analysis of Sound Propagation Focusing on Kinetic Theory of Gases
1W090198-2 小坂 勇人 指導教員 及川 靖広
KOSAKA Yuto Prof. OIKAWA Yasuhiro
概要： 本研究は，音波伝播を従来の波動方程式に基づいた解析ではなく分子シミュレーションを用いて解析することで，そ の詳細な物理現象を捉える手法を確立しようとするものである．分子レベルで音波伝搬の解析を行うことができれば，音波 伝播の厳密な解析などに有用である．そこで分子動力学法（Molecular Dynamics method: MD法）を用いてシミュレー ションプログラムを作成し，分子レベルでの音波伝播解析の可能性の検討を行った．2次元音場において，単原子分子気体の 音波伝搬シミュレーションを行った結果を記す．
キーワード：気体分子運動論，分子動力学法，非線形音響
Keywords：Kinetic theory of gases, Molecular dynamics method, Nonlinear acoustics
1.
ま え が き一般的に音波伝搬の解析には，音線法，鏡像法に代 表されるような幾何音響解析，波動方程式を有限要 素法（
Finite Element Method: FEM
）や時間領域 差分法（Finite-Difference Time-Domain method:
FDTD method）を用いて計算機で解く手法などが
あり，周波数・時系列解析などに対し用途に応じて 用いられている．しかし前提条件として，これらの 方法は波動現象を線形近似したうえで解析を行うも のであり，粘性などによる非線形性が無視できない ような場合，実際の現象との誤差が生じてくる．一方で音波伝搬の現象を分子レベルでとらえるこ とが出来れば，前述のような非線形性が無視できな い音場の解析が可能になると考えられる
[1]．本研究
では，分子レベルによる音場解析を適用することで，その厳密な挙動の計算を目指す．
2.
分子運動解析2. 1
分子動力学法分子動力学法（
Molecular Dynamics method
，MD
法）とは，分子シミュレーションに用いられる計算 手法の一つであり，2体以上の原子間ポテンシャル の下に各原子に対して古典力学におけるニュートン 方程式を解いて，系の静的，動的安定構造や，動的 過程の解析を行う手法である．粒子数が
N
個の多粒子系を考える．ここでは簡略 化をねらい，各粒子はそれぞれ単原子分子とし，粒 子自体の回転運動などは考慮する必要がないものと する．粒子の質量をm i kg，粒子の速度を v i m/s，
着目する粒子がその他の粒子から受ける力を
f i
とす ると，その粒子に関する運動方程式は，m i
dv i (t)
dt = f i (t) =
∑ N
j=1( =i) |
f ij (t) (1)
となる．ここで
f ij
は，粒子j
が粒子i
に及ぼす力 である．N粒子系の場合は各粒子に対して類似の運 動方程式がN
個存在する．これらを解くことによ り，シミュレーション区間内で原子や分子の挙動の 解析を行う．図―
1
分子動力学法2. 2
原子間ポテンシャル原子間には，その相互作用によるポテンシャルエ ネルギーが発生する．本研究では原子間相互作用 ポテンシャルとしてレナード＝ジョーンズ・ポテン シャルを用いた
[2]．
u(r ij ) = 4 [( σ
r ij
) 12
− ( σ
r ij
) 6 ]
(2)
本式の微分値より粒子間にはたらく力を求め，分子 運動の解析を行う．以上を踏まえた分子動力学法の 概要図を図―
1
に示す．2. 3
計 算 条 件剛体壁に囲まれた
2
次元平面上のシミュレーショ ン領域を考える．各気体分子はそれぞれ単原子分子 として，アルゴンをシミュレーション空間に入れて 計算を行った．初期配置，初期速度を標準状態（SATP：298.15
K，10 5 Pa）に基いて考える．系内の粒子数を N，
粒子
1
つ当たりの体積をV 0
，気体の体積をV
とし たとき，その体積分率φ V
は，φ V = N V 0 /V (3)
となる．粒子間隔
d a
は，d a = √ 3
V 0 /φ V = √ 3
V /N (4)
として求められる．標準状態（
SATP
：298.15 K
，表―
1
シミュレーション条件x
方向シミュレーション長さ1.7664
μm y
方向シミュレーション長さ2.415 × 10 −3
μm
初期温度
298.15 K
初期圧力
10 5 Pa
音源周波数196.16 MHz
時間刻み幅（1ステップあたり）1.1 × 10 15 s
最大シミュレーション時間5.2 × 10 −9 s
図―
2
シミュレーションの様子10 5 Pa）の理想気体を考えると，N = 6.02 × 10 23
，V = 24.8 L
であるから，da ≈ 3.45 × 10 − 9 m
とな る．本シミュレーションでは初期位置として，この 粒子間隔毎に，シミュレーション区間内に単原子分 子を格子状に配置した．初期速度はマクスウェル分布を利用した．これは，
熱力学平衡状態において，気体分子の速度が従う分 布関数である
[3]．通常の静止した気体では，それを
構成する分子の速度ベクトルv
は，f (v) = ( m
2πkT ) 3 2
exp
( − mv T v 2kT
)
(5)
という分布に従う．ここで，mは気体分子の質量，
k
はボルツマン定数，Tは系の温度である．系の温 度が設定された場合，分子の初期速度は式(5)
の確 率密度関数を用いて確率的に設定すればよい．平面波の伝搬を観測することを目標に音波の入力 を考える．本研究ではシミュレーション区間左端で 剛体の振動板を
196.16MHz
で振動させることで音 波の入力を行った．なお，壁面上では気体分子が弾 性衝突するものとしてシミュレーションを行った．これらの計算条件とシミュレーション時間条件を 表―
1
に示す．3.
解 析 結 果シミュレーション区間を長さ
h
で分割した検査領 域を考える．分子数密度や速度分布などの熱力学平 衡状態に対する処量は，長さh
で分割した検査領域 内にある分子を数え上げ，すべての分子の位置とそ の分子速度から得ることができる．密度分布は，
ρ(x 1 , t) = m L x L y h
∑ N n=1
χ(x 1 , x 1 | n ; h) (6)
X [μm]
Time [ps]
0.2208 0.4416 0.6624 0.8832 1.104 1.3248 1.5456 1.7664
0 0.65 1.3 1.95 2.6 3.25 3.9 4.55 5.2
Pressure [Pa]
2 4 6 8 10 12 14 16 18 x 104
図―
3
解析結果の音圧分布として求められる．温度分布
T (x 1 , t)
を，T(x 1 , t) = m k b
∑ N n=1
χ(x 1 , x 1 | n ; h)
∑ ν(x 1 , x 1 | n ; h)
× {( dx i
dt
n
) 2
+ ( dy i
dt
n
) 2
+ ( dz i
dt
n
) 2 }
(7)
とすると，圧力分布
P (x 1 , t)
は，P(x 1 , t) = ρ(x 1 , t)T(x 1 , t)R (8)
として求められる．ただし，χ(x 1 , x 1 | n ; h) =
{ 1 ( x 1 − x 1 | n < = h/2) 0 (otherwise) (9)
であり，k b
はボルツマン定数，R
は気体定数である．なお本シミュレーションの場合，上式を
2
次元空間 に対応して計算を行う．シミュレーション過程における分子の挙動の様子 を図―
2
に，シミュレーション結果より解析した音 圧分布の時系列変化を図―3
に示す．図の縦方向が 時間変化，横方向がシミュレーション区間x
方向と 対応する．図―3
より，初期状態では初期条件によ り与えた標準状態の圧力1.0 × 10 5 Pa
で分布してい るが，振動板の移動に伴い，入力された圧力分布の 時間変化が観察できる．また，シミュレーション時 間が4.55 ns〜5.2 ns
へと移る際に，シミュレーショ ン区間右端において波面が反射している様子も見て とれる．4.
む す び本研究では，分子動力学法を利用したシミュレー ションにより，分子レベルでの音波伝播解析が可能 であることを示した．これにより，非線形音場の解 析も可能になるであろう．
今後は
2
原子分子以上の多原子分子で構成された 実在の空気からなる系を考え，現実に即した厳密な 音場の解析を進める所存である．参 考 文 献
