3750

汎用マルチフィジックス有限要素解析ソフトウェア COMSOL Multiphysics を取り上げて新しい計算科学工学の仕事の流れを解説していきます。今回は波動系を中心に解説していただき
チュートリアルます。実験的な方法にも触れています。なお、前号よりチュートリアル記事は1ページ目の汎
s i ci cs s を
汎用
用マ
マル
ルチ
チフ
フィ
ィジ
ジッ
ック
クス
ス有
有限
限要
要素
素解
解析
析ソ
ソフ
フト
トウ
ウェ
ェア
ア CCOOMMSSOOLL MMuul tl it pi phhyys
を取取りり上上げげてて新新ししみを本誌掲載し、続きは日本計算工学会 HP 上で公開していますので、そちらも併せてご参
い
い計
計算
算科
科学
学工
工学
学の
の仕
仕事
事の
の流
流れ
れを
を解
解説
説し
して
てい
いき
きま
ます
す。今
。今回
回は
は波
波動
動系
系を
を中
中心
心に
に解
解説
説し
してていいたただだききまま
照ください。
す
す。
。実
実験
験的
的な
な方
方法
法に
にも
も触
触れ
れて
てい
いま
ます
す。
。な
なお
お、
、前
前号
号よ
より
りチ
チュ
ュー
ート
トリ
リア
アル
ル記
記事
事は
は 11ペ
ペーージジ目目ののみみをを
本
本誌
誌掲
掲載
載し
し、
、続
続き
きは
は日
日本
本計
計算
算工
工学
学会
会 HHPP 上
上で
で公
公開
開し
して
てい
いま
ます
すの
ので
で、
、そ
そち
ちら
らも
も併
併せ
せててごご参参照照くくだだ＜チュートリアル＞

さ
さい
い。
。

COMSOL Multiphysics による計算科学工学
COMSOL Multipysics
Multipysics による計算科学工学
による計算科学工学 ––波動系（５）波動系（５）– 波動系（5）Computational
Computational Science
Science andand Engineering
Engineeringby
byusingusingCOMSOLCOMSOLMultiphysics
Multiphysics–– Wave
System(5)
Wave System (5) 橋口 真宜、藤井 知、平野 拓一、坂東 弘之

1 音場への適用

1.1 はじめに
今回は波動系を取り上げて解説します。波動は何ら
１．音場への適用
１．音場への適用
かの物理量の空間パターンが質量輸送を伴わずにエネ
１．１
１．１ はじめにはじめにルギーを伝播する現象です。音声や楽器によるコミュ
今回は波動系を取り上げて解説します。波動は何らかの物理量
今回は波動系を取り上げて解説します。波動は何らかの物理量
ニケーションや芸術活動、携帯端末によるデータ転送、の空間パターンが質量輸送を伴わずにエネルギーを伝播する現の空間パターンが質量輸送を伴わずにエネルギーを伝播する現宇宙のかなたから地球に到達する背景放射の観測、光象です。音声や楽器によるコミュニケーションや芸術活動、携象です。音声や楽器によるコミュニケーションや芸術活動、携通信などに応用されています。これらはいずれも空気
帯端末によるデータ転送、宇宙のかなたから地球に到達する背帯端末によるデータ転送、宇宙のかなたから地球に到達する背中あるいは真空中の波動伝播の応用ですが、固体の波景放射の観測、光通信などに応用されています。これらはいず
景放射の観測、光通信などに応用されています。これらはいず
れも空気中あるいは真空中の波動伝播の応用ですが、固体のれも空気中あるいは真空中の波動伝播の応用ですが、固体の動伝播を利用すると、表面弾性波（SAW）の反射具合を
橋口真宜
橋口真宜,藤井智,藤井智, , 平野拓一
平野拓一,坂東弘之,坂東弘之
観測することでパイプラインの腐食箇所の検出もできます。波は、各応用分野においては、音波、電磁波、波動伝播を利用すると、表面弾性波(SAW)の反射具合を観測す
波動伝播を利用すると、表面弾性波(SAW)の反射具合を観測す
弾性波と呼ばれます。光は電磁波の周波数が非常に高ることでパイプラインの腐食箇所の検出もできます。波は、各ることでパイプラインの腐食箇所の検出もできます。波は、各くなった場合に相当します。アインシュタインの予言
応用分野においては、音波、電磁波、弾性波と呼ばれます。光応用分野においては、音波、電磁波、弾性波と呼ばれます。光から100年経過後に、重力波（gravitational wave）が検出は電磁波の周波数が非常に大きくなった場合に相当します。アは電磁波の周波数が非常に大きくなった場合に相当します。ア
されました。波は一般相対性理論とも関係しています。
イ
イン
ンシ
シュ
ュタ
タイ
イン
ンの
の予
予言
言か
から
ら 100
100 年
年経
経過
過後後にに、、重重力力波波
ここでは、音波、電磁波を取り上げます。
(gravitationalwave)が検出されました。波は一般相対性理論と(gravitational wave)が検出されました。波は一般相対性理論とも関係しています。ここでは、音波、電磁波を取り上げます。
も関係しています。ここでは、音波、電磁波を取り上げます。
1.2 音波と電磁波の類似性と差異
１．２
１．２ 音波と電磁波の類似性と差異
音波と電磁波の類似性と差異

音波における音圧 p（大気圧からの変動分）の支配方
音波における音圧
音波における音圧p（大気圧からの変動分）の支配方程式を図p（大気圧からの変動分）の支配方程式を図
程式を図1に示します。音場全体が同じ周波数（角周波１に示します。音場全体が同じ周波数（角周波数ω）で調和振１に示します。音場全体が同じ周波数（角周波数ω）で調和振

筆者紹介
筆者紹介
筆者紹介
はしぐち まさのり Masanori HASHIGUCHI
はしぐちまさのり
まさのり Masanori HASHIGUCHI
はしぐち
計測エンジニアリング(株)第１技術部部
計測エンジニアリング(株)第１技術部部
計測エンジニアリングシステム（株）第1技術部部長。日本機械学会認定計算力学技術者上級長。日本機械学会認定計算力学技術者上級アナリ長。日本機械学会認定計算力学技術者上級
アナリスト。
スト。
アナリスト。

数ω）で調和振動をする場合にはヘルムホルツ方程式に動をする場合にはヘルムホルツ方程式に変換できます。この場
動をする場合にはヘルムホルツ方程式に変換できます。この場
変換できます。この場合は空間微分を含む代数方程式
合は空間微分を含む代数方程式になり、解きやすくなります。
合は空間微分を含む代数方程式になり、解きやすくなります。

になり、解きやすくなります。応用では時間的に変化
応用では時間的に変化する過渡状態よりは定常振動状態を考察
応用では時間的に変化する過渡状態よりは定常振動状態を考察する過渡状態よりは定常振動状態を考察する場合が多する場合が多いのでこの性質は有用です。COMSOL
は時間領する場合が多いのでこの性質は有用です。COMSOL は時間領いのでこの性質は有用です。COMSOL は時間領域と周域と周波数領域、いずれも解析を行なうことができます。
域と周波数領域、いずれも解析を行なうことができます。
波数領域、いずれも解析を行なうことができます。

第1章
第１章
第１章

第2章
第２章
第２章

ふじい さとし
さとし Satoshi FUJII
ふじい
ふじい さとし Satoshi FUJII
筑波大・基礎工卒、修士課程了。京工繊大学院博筑波大・基礎工卒,修士課程了．京工繊大学
筑波大・基礎工卒,修士課程了．京工繊大学士課程了。住友電工、エプソンを経て平21千葉大院博士課程了．住友電工,エプソンを経て平
院博士課程了．住友電工,エプソンを経て平特任教授、平26年東工大特任教授（現兼務）、平27
21 千葉大特任教授,平 26 年東工大特任教授
21 千葉大特任教授,平26 年東工大特任教授
年沖縄高専（現職）、電子デバイス、マイクロ波工（現兼務）,平 27 年沖縄高専（現職）,電子学を専門とし、ダイヤモンド
SAW やマイクロ波加（現兼務）,平 27 年沖縄高専（現職）,電子デバイス,マイクロ波工学を専門とし、ダイ熱の研究とドローンの社会実装に従事
デバイス,マイクロ波工学を専門とし、ダイ
ヤモンド SAW やマイクロ波加熱の研究とドヤモンド
SAW やマイクロ波加熱の研究とドひらのたくいち
ローンの社会実装に従事
平成10名工大・工・電気情報卒。平成12東工大・ローンの社会実装に従事
ひらのたくいち Takuichi HIRANO
電気電子・修士課程了。平成14平成28東工大 国際
ひらのたくいち Takuichi HIRANO
開発工学専攻助教。平成28-30東工大 環境・社会理
工学院助教。平成30現在 東京都市大学知識工学部准教授。博士（工学）。電磁界解析、ミリ波測定技

第3章
第３章
第３章

術およびオンチップアンテナ等の研究に従事。著図1 音波と電磁波の基礎方程式
図１ 音波と電磁波の基礎方程式
図１ 音波と電磁波の基礎方程式

書（共著）「電磁気学」など。

続きは Web で
日本計算工学会誌「計算工学（Vol.23, No.2）
HP:で
続きは」
Web
早稲田大・博士課程了。博士（工学）。ポスドクとばんどう ひろゆき Hiroyuki Bando
続きは Web でhttp://www.jsces.org/activity/journal/
ばんどう ひろゆき Hiroyuki Bando （CRESTO） 日本計算工学会誌「計算工学(Vol.23, No.2)」HP:
して大阪電気通信大学、理化学研究所
ばんどう ひろゆき

を経て、平13年千葉大工学部助手。千葉大大学院
工学研究院助教（現職）。専門は、結晶成長、半導

日本計算工学会誌「計算工学(Vol.23, No.2)」HP:
注．COMSOL, COMSOL Multiphysics, COMSOL Desktopおよび Livelinkは
http://www.jsces.org/Issue/Journal/
COMSOL ABの商標または登録商標です。その他の製品名、サービ
http://www.jsces.org/Issue/Journal/

ス名、組織名は、各組織の商標または登録商標です。なお、本文中に™、®マークは明記しておりません。本稿の説明は著者独自の考えに基づいており他を代表するものではありません。

体物性。

第4章
第４章
第４章

（20）

計算工学

チュートリアルチュートリアルCOMSOL Multiphysics による計算科学工学 – 波動系（5）図１に示すように音圧の方程式は時間に関する 2 階微るもの、二重極は 2 個の単極が近接しているようなも

図1に示すように音圧の方程式は時間に関する2階微の、4重極は空気の噴流中の渦から放射されるようなもの、4 重極は空気の噴流中の渦から放射されるような分を含みます。このことは波動にとって重要な事項で分を含みます。このことは波動にとって重要な事項でものです。
す。質量－ダンパ－バネ系という機械振動系を思い出 のです。
す。質量−ダンパ−バネ系という機械振動系を思い出
すと、質量による慣性項が無視できる場合、系はダンパによる時間１階微分をもつ時間遅れ系となり、指数
すと、質量による慣性項が無視できる場合、系はダン関数解をもつだけで振動解を持ち得ません。また、時パによる時間1階微分をもつ時間遅れ系となり、指数関間に関する 2 階微分によって音波は遅延時間の形で因数解をもつだけで振動解を持ち得ません。また、時間
果律を満たします。そこには有限な伝播速度が存在し、に関する2階微分によって音波は遅延時間の形で因果律
音波は有限な音速で伝わります。
を満たします。そこには有限な伝播速度が存在し、音今度は電磁気における有名なマックスウェル方程式
波は有限な音速で伝わります。
を見てみましょう。音圧 p はスカラー量ですが、マックスウェル方程式から導出される電磁波を支配する波
今度は電磁気における有名なマックスウェル方程式
動方程式は電場 pEはスカラー量ですが、マッと磁場 H に関する各々のベクトルを見てみましょう。音圧量に関する時間 2 階微分項を含む方程式となります。
クスウェル方程式から導出される電磁波を支配する波
図１には真空中の波動方程式を示しています。この場
動方程式は電場
E と磁場 H に関する各々のベクトル量
合にも有限な大きさの波の伝播速度が存在し、それはに関する時間2階微分項を含む方程式となります。図1
光速になります。
このように式の形は全く同じです。一方で、大きなには真空中の波動方程式を示しています。この場合に違いもあります。音圧 p は縦波です。つまり、波の進も有限な伝播速度が存在し、それは光速になります。
図2 係数形式 PDE インターフェースによる音計算
図２ 係数形式 PDE インターフェースによる音計算
行方向に振動します。電磁波は横波です。進行方向に
このように式の形は全く同じです。一方で、大きな垂直な面内で振動しながら進んでいきます。電磁波で英国の数理科学者であるライトヒルは当時問題で違いもあります。音圧
p は縦波です。つまり、波の進
それでは、実用的な系として、複数のギアで構成されは Transverse という用語が良く出てきます。これはあったジェット機の騒音発生機構の検討を依頼され、行方向に振動します。電磁波は横波です。進行方向にる機械振動系から放射される音場の計算を見てみます。
横波の特徴を表しています。
電車の中で解決法を思いつき、封筒に式を書いたよう垂直な面内で振動しながら進んでいきます。電磁波でこれは、COMSOL のもつ応力変位も計算するマルチ
１．３ 音場計算の実際
ボディダイナミックスモジュールと圧力音響モジュー
では、COMSOL
Multiphysics
の係数形式
PDE(偏
です。ここでも封筒が活躍しました。この式がライト
は Transverse という用語が良く出てきます。これは横波ルを連成させています。皆さんにとっては、音響と構
微分方程式)インターフェースを利用して、
周波数領域
ヒル方程式です。圧縮性流体力学の式を変形して波動の特徴を表しています。
造の連成メカニズムであると思いますので、そこの部での解析法を使って平面波の伝播を計算してみましょ
方程式を導出し、ヘルムホルツ方程式とのアナロジー
分を説明します。
う。本チュートリアルの第１回目では一般形式 PDE
1.3 音場計算の実際
ギアの運動は時間依存解析で行います。その部分でを紹介しましたがそこを係数形式 PDE に読み替える からライトヒルテンソルが音源になる、つまり、流れでは、COMSOL
Multiphysicsの係数形式 PDE
（偏微分
機械振動系の音響放射面での時間変位を算出できます。
だけです。
空間次元２Ｄで係数形式
PDE
をスタディ： の中に音源があるという類推をしました。ライトヒル
続いて、
その時間変位を定常で、モデルビルダーに読み込んだ後、係数形式 が 偉
大なのは さ ら に そFFTの 式解析によって周波数領域
から放射音のパワーが方程式）インターフェースを利用して、周波数領域でのの振幅表現に変換します。音場は周波数領域で解析さPDE の設定ウィンドウで、拡散係数 c に 1、吸収係数 ジェットノズルの吹き出し速度の8乗に比例するという音波の式を使って平面波の伝播を計算してみましょう。
れるので、放射面の加速度振幅を音源として利用でき
a に－k^2 を入力します。波数 k はグローバル定義の8乗則を結論し、ジェットエンジンの騒音低減に関する
本チュートリアルの第1回目では一般形式 PDE を紹介します。あとは空気中の音場を解くだけです。結果を図パラメタで音速、周波数による算出式を記述しておき
設計指針まで導出した点です。このように、騒音にはましたがそこを係数形式
PDE に読み替えるだけです。
３に示します １）。
ます。ジオメトリは矩形とし、幅 L、高さ H で作成し空間次元2Dます。ここでの従属変数はで係数形式 PDE をスタディ：定常で、モデu という表式であり、音圧 機械振動による騒音の他に、流体から出る流体騒音があります。吹き出し速度の大きさをマッハ数で表す
ルビルダーに読み込んだ後、係数形式
PDE の設定ウィ
p に相当します。境界条件は左右の境界を流束境界(ノイマン条件に相当)、左側境界に-i*k*(u-2)、右側境界と、マッハ数が0.2を超えるあたりから放射音のパワー
ンドウで、拡散係数
c に1、吸収係数 a に− k^2を入力しに-i*k*uを設定します。有限要素メッシュは
3角形メが急に強くなるという実測と一致します。これは4重極ます。波数 k はグローバル定義のパラメタで音速、周波ッシュとし、最大要素サイズを波長 lambda0 の 6 分
音源によるものです。
数による算出式を記述しておきます。ジオメトリは矩の１であるように設定します。波動を計算する場合のそれでは、実用的な系として、複数のギアで構成さ形とし、幅 L、高さH
で作成します。ここでの従属変数
重要事項は波長の解像です。真っ先にこの点を注意しれる機械振動系から放射される音場の計算を見てみまは u という表式であり、音圧pに相当します。境界条件てください。仮に、波長の解像度は十分確保されているメッシュ系だとしても、もしも領域内に物体があるす。これは、COMSOL のもつ応力変位も計算するマル
は左右の境界を流束境界（ノイマン条件に相当）とし場合にはその形状の曲率などを解像できるメッシュでチボディダイナミックスモジュールと圧力音響モて、左側境界に-i*k*(u-2)、右側境界に -i*k*u を設定しあるかを検討することが重要です。メッシュの役目はジュールを連成させています。皆さんの興味は、音響ます。有限要素メッシュは3角形メッシュとし、最大要波動の解像と、物体形状の解像という二つの役目を果と構造の連成メカニズムにあると思いますので、そこ
素サイズを波長lambda0の6分の1であるように設定しま
たさなければなりません。
す。波動を計算する場合の重要事項は波長の解像です。
いままで見てきた音圧の計算は音源が無いというも の部分を説明します。
図３ 振動系と放射音場の連成解析
ギアの運動は時間依存解析で行います。その部分で真っ先にこの点を注意してください。仮に、波長の解のでした。実際には音源があります。ここでは、音源の取り扱い方を説明します。
音場の計算は広い空間を計算する場合があります
機械振動系の音響放射面での時間変位を算出できま
像度は十分確保されているメッシュ系である場合には、音源には単極音源、二重極音源、四重極音源といっが、波長の 6 分の１（2 次の形状関数利用時）を満たす。続いて、その時間変位をFFT 解析によって周波数領域内に物体がある場合にはその形状の曲率などを解たものがあります。単極音源は伸縮する物体が誘起すす解像度を確保するには膨大な計算になる可能性があ
領域の振幅表現に変換します。音場は周波数領域で解像できるメッシュであるかを検討することが重要です。
析されるので、放射面の加速度振幅を音源として利用
メッシュの役目は波動の解像と、物体形状の解像とい計算工学(??-2)
Vol.22, No.2 2017
できます。あとは空気中の音場を解くだけです。結果
う二つの役目を果たさなければなりません。
を図3に示します [1]。
いままで見てきた音圧の計算は音源が無いというも音場の計算は広い空間を計算する場合がありますのでした。実際には音源があります。ここでは、音源が、波長の6分の1（2次の形状関数利用時）を満たす解の取り扱い方を説明します。
像度を確保するには膨大な計算になる可能性がありま
音源には単極音源、二重極音源、四重極音源といっす。計算量を減らし、いままで困難であった計算を実たものがあります。単極音源は伸縮する物体が誘起す現するために COMSOL では境界要素法（BEM）を導入
るもの、二重極は2個の単極が近接しているようなも計算工学（20-2）

Vol.23, No.2 2018

E をスタディ：
後、係数形式に 1、吸収係数
ローバル定義のを記述しておきさ H で作成し式であり、音圧を流束境界(ノ-2)、右側境界
シュは 3 角形メ
mbda0 の 6 分計算する場合のこの点を注意し分確保されてい内に物体がある
きるメッシュでッシュの役目は二つの役目を果

機械振動系の音響放射面での時間変位を算出できます。
続いて、その時間変位を FFT 解析によって周波数領域の振幅表現に変換します。音場は周波数領域で解析されるので、放射面の加速度振幅を音源として利用でき
チュートリアルます。あとは空気中の音場を解くだけです。結果を図３に示します １）。

COMSOL Multiphysics による計算科学工学 – 波動系（5）を考慮する必要もあり、実際の系へ適用する際にはいろいろと注意を払う必要があります。

2

チュートリアル

図3 振動系と放射音場の連成解析ります。計算量を減らし、いままで困難であった計算
図３ 振動系と放射音場の連成解析を実現するために COMSOL では境界要素法(BEM)をし、有限要素法と境界要素法のハイブリッド解析がで導入し、有限要素法と境界要素法のハイブリッド解析
音 場[2]の計算は広い空間を計算する場合があります
きます。一般に計算境界での波の不必要な反射をふができます ２）。一般に計算境界での波の不必要な反射
重極音源といっが、波長の6 分の１（2次の形状関数利用時）を満た
せぐために完全整合層（PML）を利用しますが、BEM を
チュートリアルをふ せ ぐ た め に 完 全 整 合 層 (PML)を 利 用 し ま す が 、る物体が誘起す
す解像度を確保するには膨大な計算になる可能性があ
利用できる場合であれば代替えができます。さらに周BEM を利用できる場合であれば代替えができます。
波数の高い音を計算する場合には幾何音響を利用しまります。計算量を減らし、いままで困難であった計算
さらに周波数の高い音を計算する場合には幾何音響を(??-2)す。COMSOL では幾何音響は粒子のイメージで音線をVol.22, No.2 2017を実現するために
COMSOLでは境界要素法(BEM)を利用します。COMSOL
では幾何音響は粒子のイメー
追跡する方法を利用できます。
導入し、有限要素法と境界要素法のハイブリッド解析
ジで音線を追跡する方法を利用できます。
[2]２）
２）
自動車に適用した例を図4に示します。。
ができます。一般に計算境界での波の不必要な反射自動車に適用した例を図４に示します
を ふ せ ぐ た め に 完 全 整 合 層 (PML)を 利 用 し ま す が 、BEM を利用できる場合であれば代替えができます。
さらに周波数の高い音を計算する場合には幾何音響を利用します。COMSOL では幾何音響は粒子のイメー
ジで音線を追跡する方法を利用できます。
自動車に適用した例を図４に示します ２）。

が無いというも
ここでは、音源

図4

FEM、BEM および粒子的手法による音場解析

図４FEM, BEM および粒子的手法による音場解析
１．４ まとめ
周波数領域で音波の計算をする上での必須事項は
音場は圧力音響を扱う場合にはスカラーであるのでメッシュ寸法の調整による波長の空間解像度の確保に
電磁波に比べて取扱いは容易です。一方で、狭い流路なります。波長は周波数が高くなると短縮するので周を通過する音場では壁境界層における粘性や熱の影響波数が低い方がメッシュ最大寸法を緩和でき、有限要を考慮する必要もあり、実際の系へ適用する際にはいろいろと注意を払うが必要であります。
素法を適用しても計算機のメモリを多く使わずに済み
図４FEM, BEM および粒子的手法による音場解析
ます。周波数が上がってくるとメッシュ最大寸法を小１．４ まとめ
参考文献
さくしてゆく必要があり、有限要素法では計算しづら
音場は圧力音響を扱う場合にはスカラーであるので
１）
くなってきます。不連続ガラーキン法を利用すればメ
電磁波に比べて取扱いは容易です。一方で、狭い流路
https://www.comsol.jp/blogs/how-to-model-gearbox-vibr
モリの使用量が減ります。BEMで計算できる場合にはを通過する音場では壁境界層における粘性や熱の影響
ation-and-noise-in-comsol-multiphysics/
計算法を FEM から BEM に変更します。場を離散化すを考慮する必要もあり、実際の系へ適用する際にはいる方法で扱えない場合には幾何音響の考え方に基づいろいろと注意を払うが必要であります。
２）て粒子運動による音線追跡をします。
http://www.aes.org/conferences/2017/automotive/downlo
COMSOL Multiphysics は図4にあるようにこれらの手
参考文献
ads/AES_Modeling_Car_Cabin_Acoustics_2017.pdf
法を選択できます。
１）
https://www.comsol.jp/blogs/how-to-model-gearbox-vibr

1.4 まとめ
ation-and-noise-in-comsol-multiphysics/
２．マイクロ波化学へのシミュレーションの適用
音場は圧力音響を扱う場合にはスカラーであるので
２．１ はじめに電磁波に比べて取扱いは容易です。一方で、狭い流路
２）マイクロ波化学では、マイクロ波を照射することで、を通過する音場では壁境界層における粘性や熱の影響
http://www.aes.org/conferences/2017/automotive/downlo
迅速加熱、高効率、高速熱応答、局所もしくは選択加
ads/AES_Modeling_Car_Cabin_Acoustics_2017.pdf
熱、非熱的効果などの現象が報告されています。これ

マイクロ波化学へのシミュレーションの適用

2.1 はじめにマイクロ波化学では、マイクロ波を照射すること
で、迅速加熱、高効率、高速熱応答、局所もしくは選択加熱、非熱的効果などの現象が報告されています。
これらの現象を解明し、劇的に化学反応に必要なエネ
ルギーを削減出来る技術を完成することが求められています。特に、鉄鋼や化学などの素材製造する産業では、大量のエネルギー消費に加え、殆ど、省エネル
返し、マイクロ波化学を汎用的な技術として完成させギーが進んでない現状があります。[3]ることを目指しています。これまで、マイクロ波によ
そこで、我々の東工大・和田研究グループでは、マる加熱現象について、新しい知見を幾つか発見しましイクロ波化学における物理・化学現象のモデル化を行た。 2) その中で、複雑な物理現象を解析できるシミュ
い、シミュレーションと実験により検証や実証を繰りレーションの役割は大きいが、シミュレーション上で返し、マイクロ波化学を汎用的な技術として完成させは、実験室で実施する実験をすべての物理を取り込め
返し、マイクロ波化学を汎用的な技術として完成させることを目指しています。これまで、マイクロ波によているわけでないことやナノからメータサイズのマル
ることを目指しています。これまで、マイクロ波による加熱現象について、新しい知見を幾つか発見しましチスケールまで考慮しなければならない課題がありまる加熱現象について、新しい知見を幾つか発見しましす。現実現象がシミュレーションできない場合、全てた。[4] その中で、複雑な物理現象を解析できるシミュ
2)
た。
その中で、複雑な物理現象を解析できるシミュ計算手法により達成できていないものか、マイクロ波
レーションの役割は大きいが、シミュレーション上でレーションの役割は大きいが、シミュレーション上でによる反応促進のような新しい物理現象なのかを区別
は、実験室で実施する実験をすべての物理を取り込め
は、実験室で実施する実験をすべての物理を取り込めすることが極めて重要です。正確な実験データの取得ているわけでないことやナノからメータサイズのマルや、シミュレーションで得られた解析結果を十分に比ているわけでないことやナノからメータサイズのマル
チスケールまで考慮しなければならない課題がありま較しながら、検証を行うことが必要です。本解説では、チスケールまで考慮しなければならない課題があります。現実現象がシミュレーションできない場合、全て幾つかのモデルについて解析した例を紹介し、今後のす。現実現象がシミュレーションできない場合、全て計算手法により達成できていないものか、マイクロ波
シミュレーションとしての課題について述べます。
計算手法により達成できていないものか、マイクロ波による反応促進のような新しい物理現象なのかを区別による反応促進のような新しい物理現象なのかを区別することが極めて重要です。正確な実験データの取得
２．２ マイクロ波化学における計算手法
することが極めて重要です。正確な実験データの取得や、シミュレーションで得られた解析結果を十分に比解析手法については他の章の先生方にお任せし、今や、シミュレーションで得られた解析結果を十分に比較しながら、検証を行うことが必要です。本解説で回、COMSOL
マルチフィジックスをマイクロ波化学に適較しながら、検証を行うことが必要です。本解説では、は、幾つかのモデルについて解析した例を紹介し、今用した例を述べます。マイクロ波化学とは、その名前幾つかのモデルについて解析した例を紹介し、今後の後のシミュレーションとしての課題について述べます。
が示すとおり、マイクロ波照射による化学反応プロセ
シミュレーションとしての課題について述べます。
スです。

2.2 マイクロ波化学における計算手法
２．２ マイクロ波化学における計算手法
真空引き
解析手法については他の章の先生方にお任せし、今解析手法については他の章の先生方にお任せし、今導波管キャビティ
回、COMSOL マルチフィジックスをマイクロ波化学に回、COMSOL マルチフィジックスをマイクロ波化学に適TE103モード
るつぼ
適用した例を述べます。マイクロ波化学とは、その名用した例を述べます。マイクロ波化学とは、その名前前が示すとおり、マイクロ波照射による化学反応プロパワーが示すとおり、マイクロ波照射による化学反応プロセ
マイクロ
試料
セスです。
メータ
波源
スです。
真空引き
赤外線温度計導波管キャビティ

TE103モード
図 1 シングルモードキャビティの実験装置
るつぼすなわち、支配方程式はマックスウエルの方程式と伝熱方程式、材料の温度が決定されると化学反応はアレ
パワー
マイクロ
試料
メータ
波源
ニウス則の速度方程式に支配されます。電磁場はマックスウエルの方程式で記述するものの、周波数領域で赤外線温度計
計算することから、次のように簡略化できます。

図 1 図5
シングルモードキャビティの実験装置
シングルモードキャビティの実験装置 (1)すなわち、支配方程式はマックスウエルの方程式と伝ここで、E は電場、µは透磁率、εは誘電率、は角速度、すなわち、支配方程式はマックスウエルの方程式と熱方程式、材料の温度が決定されると化学反応はアレ
ｃは光速です。次に、伝熱の方程式は、次のように記
伝熱方程式、材料の温度が決定されると化学反応はア
ニウス則の速度方程式に支配されます。電磁場はマッ述されます。
レニウス則の速度方程式に支配されます。電磁場はクスウエルの方程式で記述するものの、周波数領域で(2)計算することから、次のように簡略化できます。

計算工学
Vol.23, No.2
2018
T は温
ここで、ρは密度、c o は比熱、ｋは熱伝達係数、らの現象を解明し、劇的に化学反応に必要なエネルギ（20-3）ーを削減出来る技術を完成することが求められていま
２．マイクロ波化学へのシミュレーションの適用す。特に、鉄鋼や化学などの素材製造する産業では、
２．１
はじめに大量のエネルギー消費に加え、殆ど、省エネルギーが
1)
進んでない現状があります。
マイクロ波化学では、マイクロ波を照射することで、(1)度、 Q はマイクロ波加熱を扱うことから、
ここで、E は電場、µは透磁率、εは誘電率、は角速度、(3)ｃは光速です。次に、伝熱の方程式は、次のように記
発熱量 Q は材料中の電磁波の損失、 H は磁場です。こ述されます。

(3)χ ιjκl-gearbox-vibr
arbox-vibr

赤外線温度計
赤外線温度計図 1 シングルモードキャビティの実験装置
図 1ε シングルモードキャビティの実験装置

NL
(3)
0
P
=
(ω ; −ω , ω , ω ) ⋅ E j (ω , k ) Ek (ω , k ) Ei (ω , k )χ ιjκl∑i すなわち、支配方程式はマックスウエルの方程式と伝すなわち、支配方程式はマックスウエルの方程式と伝チュートリアルCOMSOL Multiphysics による計算科学工学 – 波動系（5）4 jki
熱方程式、材料の温度が決定されると化学反応はアレ

熱方程式、材料の温度が決定されると化学反応はアレニウス則の速度方程式に支配されます。電磁場はマッ∂T
ニウス則の速度方程式に支配されます。電磁場はマッマックスウエルの方程式で記述するものの、周波数領いう画期的なものです。[5] まず、図5のシングルモード
ρc0クスウエルの方程式で記述するものの、周波数領域で− ∇ k∇ T = Q
motive/downlo∂tクスウエルの方程式で記述するものの、周波数領域でive/downlo域で計算することから、次のように簡略化できます。
の実験設備に示す通り、反応温度はアルミナ製るつぼ
計算することから、次のように簡略化できます。
7.pdf
計算することから、次のように簡略化できます。
df2の表面温度から推測するしかなく、次の疑念が生じまωε(3)(3) µ −1∇×Ε ) =(1)∇×(χχ ijκl(1)ιjκl(1) す。真空雰囲気下で断熱状態にあるものの、試料の内c2
部温度は違うのではないかという点があります。この
ここで、Eは電場、μは透磁率、εは誘電率、 は角速度、
は角
ここで、E は電場、µは透磁率、εは誘電率、
ンの適用
NL 1 ε*
ここで、εは誘電率、
は角速度、の適用
ε00E は電場、
2χ (3)
2⋅ E
2) E (ω , k ) 試料の内部温度を確からしい理由で推定するために、
NL(3)(ωµは透磁率、
P
=
;
−
,
,
)
(
,
k
)
E
(
,
k
ω
ω
ω
ω
ω
i ｃは光速です。次に、伝熱の方程式は、次のように記j f速度、cは光速です。次に、伝熱の方程式は、次のようP=(ω , k )l Ei (ω , k )Qσ E χ+ijκlιjκlπ(ωf ε; −0εωr,′′ωE, ω )+⋅ Eπµ,0kµ)rkE′′ k Hi =j (ωｃは光速です。次に、伝熱の方程式は、次のように記
マルチフィジックスシミュレーションを行いました。
2 44 jkljki
述されます。
射することで、に記述されます。
述されます。
することで、
その結果を図6に示します。図6(a) に示すとおり、試料
しくは選択加 ρc ∂
∂TT
くは選択加
(2)
(2)
ρc00 ∂t −−∇∇kk∇∇TT ==QQ
(2) はキャビティ中の磁場が最も強いところに置かれていいます。これ∂tます。これること、また、図6(b) からは、マイクロ波照射中の温2
ここで、ρは密度、c
要なエネルギω ε2 oo は比熱、ko は比熱、ｋは熱伝達係数、T は温−1は比熱、ｋは熱伝達係数、
T は温 度変化について、実測とシミュレーション結果についここで、ρは密度、c
なエネルギ
ここで、
は密度、c
は熱伝達係数、T
は
ω
ε
µ
∇×(
∇
=
×
Ε
)
1Q はマイクロ波加熱を扱うことから、
度、められていま 度、µQ −はマイクロ波加熱を扱うことから、∇×(∇×Ε ) =
c 22
はマイクロ波加熱を扱うことから、られていま 温度、Qて、立上りについては、やや不一致が見られるもの
c
る産業では、
産業では、
(3)
1
2
2
2
(3) (3) の、マイクロ波照射後、約17分前後で試料温度が飽和
エネルギーが Q =1 σ E 2 + π f ε 0ε r′′ E 2+ π f µ0 µr′′ H 2
ネルギーが2 σ E + π f ε 0ε r′′ E + π f µ0 µr′′ H
Q=
していること、また、試料内部とるつぼ温度を計算しQ は材料中の電磁波の損失、 H は磁場です。こ発熱量
2 QQ は材料中の電磁波の損失、Hは材料中の電磁波の損失、 Hは磁場です。
は磁場です。こ た結果を示しているものの、差がないためグラフが重発熱量
発熱量
こで、幸いなことに、マイクロ波は光速度であること
ープでは、マ こで、幸いなことに、マイクロ波は光速度であること
プでは、マ ここで、幸いなことに、マイクロ波は光速度であるこなっており、シミュレーションの結果ではるつぼ側壁から、周囲の境界条件をマイクロ波は一瞬にて感じ取モデル化を行 から、周囲の境界条件をマイクロ波は一瞬にて感じ取
デル化を行 とから、周囲の境界条件をマイクロ波は一瞬にて感じと試料内部の温度の差は殆どないことが示されましり、反応場が形成されます。一方、伝熱に関しては、や実証を繰り り、反応場が形成されます。一方、伝熱に関しては、実証を繰り 取り、反応場が形成されます。一方、伝熱に関してた。つまり、従来法では、1600℃で還元をおこなっていたものが、マイクロ波照射では安い還元材のカーボは、大きさにもよるものの、今、扱う大きさとして(??-3)
Vol.22, No.2 2017
ンを使って約半分に近い850℃で還元が出来ることが示(??-3)
Vol.22,No.22017は、秒単位以上の遷移状態を知りたい場合が殆どでされました。
す。以上のことから、材料の温度変化による誘電率や2つ目の事例として、マルチモードキャビティにて、導電率などの物性定数の変化がない場合、電磁波及び電磁場照射中に電磁波を拡散させる目的で導入されて伝熱の支配方程式を別々に解き、解を求めます。それいるスタラファンの回転を組み込んだ計算について示
らが変化する場合には前述したように支配方程式の時します。前述の例では、マイクロ波による還元反応が間スケールが異なることから、一旦、伝熱計算まで求非常にうまく行った例を示しました。しかしながら、めたのち、再度、電磁波の方程式を解くことで、マル還元した材料は1g 未満です。キャビティも小さいこと
チフィジックスシミュレータとして収束させることがから、実用化には程遠い。そこで、図7に示すマルチ可能です。これ以外にも、例えば、電子レンジのよう
チュートリアルモードキャビティを使って少し大きい実験装置の場合に、マイクロ波を拡散・均一照射させるために、ターについて示します。ここでの解析の目的は、実際に試ンテーブルやスタラファンがある場合、有限要素法で∑∑850℃で還元が出来ることが示されました。
大きさにもよるものの、今、扱う大きさとしては、秒はそもそもメッシュが変わってしまう問題があり、小単位以上の遷移状態を知りたい場合が殆どです。以上
さい変化だと変形メッシュで対応できるものの、大きのことから、材料の温度変化による誘電率や導電率なく形状が変化すると変形メッシュでは対応できなくな
どの物性定数の変化がない場合、電磁波及び伝熱の支ります。しかし、この場合も光速で電磁場の反応場が配方程式を別々に解き、解を求めます。それらが変化形成される現象と機械が動きの速度差から、統計力学の考え方であるエルゴールドの仮説を適用し、時間にする場合には前述したように支配方程式の時間スケー
よる形状変化は場合の数（状態数）で置き換え、各状態
ルが異なることから、
一旦、伝熱計算まで求めたのち、をそれぞれ計算し、それらを平均すれば良いことにな
再度、電磁波の方程式を解くことで、マルチフィジッ
ります。この手法は、パソコンの計算速度が劇的に向
クスシミュレータとして収束させることが可能です。
上したことで可能となったものであり、数十年前ではこれ以外にも、例えば、電子レンジのように、マイク考えられなかったものです。
ロ波を拡散・均一照射させるために、ターンテーブル(a) キャビティ内の磁場分布やスタラファンがある場合、有限要素法ではそもそも
2.3 解析の事例紹介
メッシュが変わってしまう問題があり、小さい変化だ
マイクロ波照射の事例について2つ紹介します。1つと変形メッシュで対応できるものの、大きく形状が変
目はシングルモードキャビティを用い、電磁場と伝熱化すると変形メッシュでは対応できない。しかし、このマルチフィジックスシミュレーションを行った例での場合も光速で電磁場の反応場が形成される現象と機す。この実験では、マイクロ波照射下によりフッ化ス
カンジウムをカーボンにより還元する実験を行いまし 1
械が動きの速度差から、統計力学の考え方であるエル
た。従来法では、フッ化スカンジウムは、カルシウム
ゴールドの仮説を適用し、時間による形状変化は場合
金属を使って凡そ1600℃の反応温度により還元が実施
1の数（状態数）で置き換え、各状態をそれぞれ計算し、されています。我々の試みた還元技術では、還元材をそれらを平均すれば良いことになります。
この手法は、（b）試料内部の温度変化
カルシウムから扱いやすいカーボンに変更し、マイク
パソコンの計算速度が劇的に向上したことで可能とな1
図 図6
2 シミュレーション結果
ロ波照射下において反応温度850℃にて還元が出来るとシミュレーション結果ったものであり、数十年前では考えられなかったものです。
計算工学（20-4）2 つ目の事例として、マルチモードキャビティにて、
Vol.23, No.2 2018
２．３ 解析の事例紹介
電磁場照射中に電磁波を拡散させる目的で導入されてマイクロ波照射の事例について２つ紹介します。1いるスタラファンの回転を組み込んだ計算について示
つ目はシングルモードキャビティを用い、電磁場と伝します。前述の例では、マイクロ波による還元反応が熱のマルチフィジックスシミュレーションを行った例

合について示します。ここでの解析の目的は、実際にる問題について解析を行った例について紹介しました
試料は均一に電磁波が照射されているかどうかです。
ICT 技術の劇的な進展より、過去モデルが複雑すぎて
結局、平均的に試料に偏りなく照射されていることや、解析できないものが簡単に出来るようになってきましチュートリアルCOMSOL Multiphysics による計算科学工学 – 波動系（5）た。適切にツールを使い、現象との比較から新しい発見が見つかる可能性が十分にあると考えています。

、、(??-4)

図8 試料近傍の電場強度とスタラファン状態数
図 4 試料近傍の電場強度とスタラファン状態数

ています。

参考文献
本稿では、電磁界シミュレータを上手く使いこなす[1] ために必要な知識と技術について解説します。電磁界
経済産業省資源エネルギー庁、平成 25 年度エネルギー

シミュレータで使われる電磁界解析手法にもいくつかあり、それらの違いについて説明することは紙面の都図7 マルチモードキャビティの実験装置[2] Y.
Wada, S. Fujii, E. Suzuki, M. M. Maitani, S. Tsubaki,
合上不可能なので、電磁界解析手法として有限要素法
図料は均一に電磁波が照射されているかどうかです。結3 マルチモードキャビティの実験装置を用いたシミュレータの使い方に絞って説明します。
S. Chonan,
M. Fukui, & N. Inazu, “ Smelting Magnesium
局、平均的に試料に偏りなく照射されていることを確
Metal using a Microwave Pidgen Method,”
Scientific
認し、さらには、反応温度を知りたいことが、解析の3.2 電磁界シミュレータの発展と有限要素法解析の特徴
7, 46512; doi:10.1038/srep46512(2017)
ポイントです。そこで、前節で述べたように、スタラ Reports,
電磁界解析法の代表的なものとして、本稿で取り上
Vol.22,
No.2
2017
S. Fujii, S.
Tsubaki,E. Method）
Suzuki,S. Chonan M.
ファンという構造体が回転するため、スタラファンの [3] げる有限要素法（FEM;
Finite Elementの他に、回転状態を細かく分け、それぞれについて計算し、そ Fukui,
モーメント法（MoM;
MethodofMoments）
およびFDTD
and Y. Wada, ”Study on Metal Refining Process
れらの平均を取れば解が得られると考えました。どの法（Finite Element Time Domain Method）が あ り ま す。
of Sc metal using by Microwave Irradiation,” Progress in
くらい細かく分けるのかは、例えば、細分化の指標は
FEM、MoM は周波数領域の解析法、FDTD 法は時間領
Research
Symposium (PIERS), pp.
平均した電場強度が飽和したところになります。その Electromagnetic
域の解析法です。FDTD
法ではメッシュを自在に切るこ
結果を図8に示します。結果のとおり、状態数が9くら 4203-4203とができませんが、FEM、MoM
は自在に切れます。ま(2016)
いでよいことになります。状態数9は回転体を10°ずつた、FEM、MoM ともに行列方程式を解く問題に帰着さ計算してそれを平均化しました。
れます（陰解法と呼びます。対して FDTD は陽解法で
３．電磁界シミュレーターを上手く使いこなすす）。MoM では媒質境界のみに未知数を配置すればよ2.4 まとめいですが密行列になり、FEM
は空間全体に未知数を配３．１
はじめに簡単ではあるものの、COMSOL マルチフィジックス コンピュータの性能が劇的に向上し、現在シミュレー
置する必要がありますが、疎行列が得られるという違によるシミュレーションを使い、マイクロ波化学においがあります。与えられた問題により、FEM、MoM、
ションはあらゆる分野で実用的になっています。支配
ける問題について解析を行った例について紹介しまし
FDTD のどれが適しているか違いが生じる場合もありまた。ICT 技術の劇的な進展より、過去モデルが複雑す 方程式としてマクスウェルの方程式を扱う電磁界分野
すが、一般的にどの解析方法が優れていると判断する
ぎて解析できないものが簡単に出来るようになってき においても同様であり、電磁界シミュレータは高周波ことは難しいです。FEM では導波路ポート（導波路モーました。適切にツールを使い、現象との比較から新しドを入出力とする励振方法）の解析が先行して研究され回路やアンテナ等の素子設計に必要不可欠となっていい発見が見つかる可能性が十分にあると考えています。
ていたことがあり、FEM シミュレータでは導波路ポーます。
トのサポートが充実しています。

に関する年次報告（エネルギー白書 2014）3

電磁界シミュレーターを上手く使いこなす

3.1 はじめに
3.3 電磁界シミュレータを使いこなすための知識につコンピュータの性能が劇的に向上し、現在シミュ 計算工学 いてレーションはあらゆる分野で実用的になっています。
電磁界シミュレータを使いこなすために、電磁界シ支配方程式としてマクスウェルの方程式を扱う電磁界ミュレータの原理を理解している必要があるかという分野においても同様であり、電磁界シミュレータは高質問に対し、少し前であれば基礎をしっかり勉強する
周波回路やアンテナ等の素子設計に必要不可欠となっ必要があるという意見が大半でありました。しかし、計算工学（20-5）

Vol.23, No.2 2018

チュートリアル

電卓を使うのにブール代数と論理回路のロジックを理解している人はどれだけいるでしょうか。この問題を考えても、内部の原理を知らず、ツール（シミュレー
タ）の内部はブラックボックスであっても上手く使える可能性があるのではないでしょうか。ただし、上手く使うためには最低限の知識も必要となります。電卓の例を見るならば、電卓というツール自体はブラックボックスであっても、時間をかければ筆算をして結果を検証する能力をユーザーは持っていなければうまく使えないでしょう。電磁界シミュレータを電卓と対比させるならば、扱う問題が異なるために多少考慮しなければならないことも増えます。この点について、3.4
節で述べます。また、電卓と同様に、結果に対して検証が行える必要もあり、そのためには解析的に解ける、あるいは実験の精度が保証されている問題（規範問題）の結果と比較することは有効な手段です [6, 7]。
3.4 電磁界シミュレーションにおける準備とモデル化
次のことに注意してシミュレーションの準備およびモデル化を行うことをお勧めします。
3.4.1 時間領域と周波数領域
シミュレータの種類（解析手法）によって、時間領域
あるいは周波数領域のどちらの解析を行うか選択することになります。時間波形と周波数はフーリエ変換で関係づけられているので、理論的にはどちらの情報も得られますが、現実問題としては計算時間に差が出てきます。例えば、広帯域な信号を扱いたい場合は、周波数領域の解析で多くの周波数での解析を行うよりも、パルス波形を入力して出力波形をフーリエ変換する時間領域の解析法を用いた方が計算時間は短くて済みます。狭帯域で良い場合は逆に周波数領域の方が計
算時間が短くてすみます。広帯域信号を扱う場合かそうでないかによって使い分けるのが良いと考えます。

COMSOL Multiphysics による計算科学工学 – 波動系（5）析に帰着せざるを得ない電磁界シミュレータでどのようにモデル化すればよいか、実際にモデル化を始めると悩むと思います。また、実際のアンテナの励振部はケーブルで給電していますが、どの部分で打ち切って、どのようにモデル化すべきか悩むと思います。これに関しては、実際の測定と比較することを考えますと、S パラメータの測定ではネットワークアナライザの校正について理解し、対応させたモデル化を行う必要があります。導波路ポートで励振すると、導波路の途中で校正した測定値と直接比較できる S パラメータが得られます。
3.4.4 規範問題を用いた解析例
図9に 導 波 管 ス ロ ッ ト ア ン テ ナ [7] の COMSOLMultiphysics での解析例を示します。方形導波管の上壁からグランド板まで長方形スロットが開いており、外部空間に放射する構造となっています。この問題を例
として解析の手順について説明します。
まずはモデル化です。FEM では解析空間を3次元物
体の組み合わせ（結合、切抜きなどの論理演算も可能）
で作成します。デフォルトでは物体の周囲境界は電気壁（PEC; Perfect Electric Conductor）となっています。物
体の内部の媒質を設定します。物体は完全に他の物体の内部に入れることができ、その場合、内部の物体の媒質定数が優先されます。そして、内部物体の周囲境
界の PEC は除去され、外部の媒質との境界となります。物体同士を接するように配置した場合は、接した面の PEC は除去され、接した2つの媒質境界となります。このようにして、図9の導波管スロットアンテナは3つの直方体（Box1〜Box3）で組み立てられます。この問題では Box1〜Box3の媒質は空気（Air）に設定しています。そして、次に境界条件の設定です。デフォルト
では周囲境界は PEC なので、そうでない境界を上書き指定する必要があります。外部領域（Box3）の下面以外の境界を吸収境界条件（COMSOL では散乱境界条件と呼ばれる）に設定します。完全整合層（PML; Perfectly
Matched Layer）という、より高性能な吸収境界条件も(dB)

3.4.2 解析の種類
チュートリアル電磁界シミュレーションというと、どこかで励振してどのように電磁波が広がっていくかをシミュレー
３．４．１ 時間領域と周波数領域z
ションすることを想像すると思います。しかし、解析
wgx =232 mm
y
シミュレータの種類（解析手法）によって、時間領
wgy =232 mmにも次のようにいくつかの種類があります。
x
域あるいは周波数領域のどちらの解析を行うか選択す● 励振問題ることになります。時間波形と周波数はフーリエ変換●
非励振問題（励振波源がない問題。固有値問題に帰外部領域(Box3)
で関係づけられているので、理論的にはどちらの情報着される）も得られるが、現実問題としては計算時間に差が出て➢ 導波路モードの解析（2次元問題）
きます。例えば、広帯域な信号を扱いたい場合は、周➢ 共振器の共振モードの解析（3次元問題）
ha =116 mm
波数領域の解析で多くの周波数での解析を行うよりも、どの解析を行っているのか、明確に把握しておく必要
パルス波形を入力して出力波形をフーリエ変換する時があります。多くのシミュレータでは特に指定しない
l2 =116 mm
間領域の解析法を用いた方が計算時間は短くて済みま場合、よく使われる「励振問題」をデフォルトとしていPort 2
す。狭帯域で良い場合は逆に周波数領域の方が計算時グランド板ることが多いです。
l1 =116 mm
間が短くてすみます。広帯域信号を扱う場合かそうでスロット(Box2)
長さ 37 mm, 幅 5mm,
ないかによって使い分けるのが良いと考えます。
Port 1管軸からのオフセット 20mm
3.4.3 モデル化（境界条件、励振モデル）
方形導波管(Box1)
３．４．２ 解析の種類
WRJ-4標準導波管 (58.1 mm x 29.1 mm, 壁厚 1.6 mm)
アンテナからの放射を考えた場合、電磁波は広い外電磁界シミュレーションというと、どこかで励振し部空間まで飛んでいきます。FEM などの有限領域の解
導波管スロットアンテナの解析例（4 GHz）
図１図9 導波管スロットアンテナの解析例(4 GHz)てどのように電磁波が広がっていくかをシミュレーシ
ョンすることを想像するのではないだろうか。しかし、計算工学（20-6）
Vol.23, No.2 20180
解析にも次のようにいくつかの種類があります。

励振問題
-5

非励振問題（励振波源がない問題。固有値問題に-10帰着される）
-15

、時間領
か選択す
リエ変換
らの情報
差が出て
合は、周
うよりも、
換する時
て済みま
が計算時
かそうで
す。

で励振し
ュレーシ
。しかし、
す。

チュートリアルありますが、後で目的を述べますように、今回はどちらを用いても結果に大きな違いはありません（大きな違いにならないようにモデル化を工夫します）。次に、励振 波 源 の 設 定 と し て 方 形 導 波 管（Box1）の 前 後2面 をポート（Port）に設定します。各 Port からは基本モードの TE10 モードで励振します。ポートは無限に長い導波路を模擬するのでポートから外に進行する波は無反射となるように動作します。これで、媒質および解析空間の境界条件全体の境界条件の設定が完了しました。
図9中の Box3のサイズは解析の都合上有限サイズに打ち切ったものであり、サイズは大きいほど実験環境に近いですが、大きすぎると計算に時間がかかります。
目的は広い空間を模擬し、スロットから放射された電磁波を吸収境界壁で吸収することなので wgx、wgy、ha をスロット放射部からある程度（指標として、散乱部から1/2波長以上）大きくすれば放射電磁界は吸収境界条件に垂直に近い角度で入射するようになって吸収特性は良くなり、S パラメータは十分収束します。同様に、ポートも無限に長い線路を模擬して吸収する必要があるので、スロットのある散乱部（不連続部）で生じた高次モードが十分減衰するように l1、l2 を大きくしてあるz程度離す（指標として、基本モードの1/2管内波長以上）
wgx =232 mmywgy =232 mm
必要があります。最後に、解析周波数範囲を設定してx解析します。FEM シミュレータは内部で細かいメッ
シュを切り、解析を行います。得られたS パラメータ
外部領域(Box3)をプロットすると図10のようになります。MoM は著者のモーメント法解析の結果であり、両者はよく一致しています。実験との一致を確認する前に、前述したよha =116 mm
うに wgx、wgy、ha、l1、l2 を大きくし、かつメッシュ分
割も変えて収束が十分であることを確認すると安心です。また、図9のように電磁界の瞬時値分布を見たり、
l2 =116 mm
アニメーションを描いて吸収していることを確認する
Port 2
グランド板
ことも検証のために有効な手段です。種々の規範問題
l1 =116 mm
スロット(Box2)のモデル化の詳細およびシミュレータのファイルの例
長さ 37 mm, 幅 5mm,
Port 1管軸からのオフセット 20mm
は [8] でダウンロードすることができます。
方形導波管(Box1)

）
広い外部
域の解析
どのよう
始めると
部はケー
て、どの
に関して
パラメ

ついて、境界条件と励振モデルが重要であり、有限空間による打ち切りの影響について導波管スロットア
レーアンテナを例として説明しました。また、結果を検証するための規範問題について紹介しました。

4

光領域、非線形光学のシミュレーションへの適用

4.1 はじめに
情報通信量の増大に伴い、通信速度のより高速化・広帯域化が求められており、応答速度が1ps 以下の超高速な光デバイスの研究が行われています。我々、千葉大・坂東研究グループでは、電子デバイスの応答速度の限界を超えるべく、電気信号を使用せず光信号を光信号のままスイッチングする全光スイッチの研究を行っており、この原理として2光子吸収の利用を検討しています [9, 10]。2光子吸収は、3次の非線形光学現象であり、光強度が強くなると物質の光吸収係数が増大し透過光強度が低下する物理現象です。そのため、2光子吸収効率のよい物質やデバイス構造を利用することにより、光信号の遮断・通過を光強度の強い制御光にて制御できます。当研究室では、半導体を結晶成長しそれらを半導体プロセスにて導波路形状などに加工した試料について、2光子吸収特性を測定し、3次非線形感
受率テンソル χ(3) を求めています。しかしこれらの測定には、フェムト秒レーザーをレンズで集光した光を使用するため、光強度が空間的にも時間的にも変化し、さらに2光子吸収を生じながらの光伝搬するため、厚膜試料に適用できる解析式がありません。そのため我々
は、実験データから χ(3) を求めるために、いくつかの仮定をいれた伝搬モデルを考案し定式化することにより求めています。そこで、求めた χ(3) の値を使ってシミュレーションし実験結果と比較検討することで、提唱モデルや得られた χ(3) 値、実験方法などの妥当性の検討を行う必要があります。 以下に、我々が検討した例を紹介します。

WRJ-4標準導波管 (58.1 mm x 29.1 mm, 壁厚 1.6 mm)

3.5 まとめ
導波管スロットアンテナの解析例 (4 GHz)
電磁界シミュレーションを使いこなすための知識に

図１

0
-5

値問題に

-10

Amplitude (dB)

問題）
おく必要
定しない
トとして

COMSOL Multiphysics による計算科学工学 – 波動系（5）

-15
-20
-25
S11 (MoM)-30

S21 (MoM)
S11 (COMSOL)

-35
-40

S21 (COMSOL)
3

3.2

3.4

3.6

3.8
4
4.2
Frequency (GHz)

4.4

4.6

4.8

5

4.2 2光子吸収の計算手法
光と物質の相互作用や物質中の光伝搬はマクセル方程式にて記述できるが、今回の場合、主に2つの注意すべき点があります。1つは光伝搬距離の長さであり、もう1つは非線形光学の記述です。
光領域のシミュレーション、特に光伝搬に関するシ
ミュレーションで問題となるのが、屈折率 n の物質中での 波 長 λ= λ0 / n (~0.1−0.5 µm) と モ デ ル 系 サ イ ズ L（>~100 µm）が大きく異なることです。光の位相情報を考慮して計算するには、有限要素法では各場所での波長の1/6以下にメッシュサイズを設定しなければ正しく計算出来ないが、それに比べて L が大きすぎるため、必要となるメモリ量が膨大になります。しかしながら、今回対象としているような系では、光強度の変化は伝搬方向に対しては非常に緩やかです。この様な場合には、
COMSOL Multiphysics の波動光学モジュールに用意されている Beam envelope インターフェイスの使用が有効で

S21
図10 導波管スロットアンテナの反射係数
導波管スロットアンテナの反射係数 S11 と透過係数
図２
S 11 と透過の振幅の周波数特性
係数 S 21 の振幅の周波数特性
媒質定数が優先されます。そして、内部物体の周囲境
計算工学（20-7）
界の PEC は除去され、外部の媒質との境界となります。
物体同士を接するように配置した場合は、接した面のPEC は除去され、接した２つの媒質境界となります。
このようにして、図１の導波管スロットアンテナは 3

Vol.23, No.2 2018

たもので
、大きす
模擬し、
で吸収す

からある

）大きく
い角度で
S パラメ
限に長い
ロットの
が十分減

す（指標

要があり
します。
切り、解
すると図
ト法解析
います。
ように

ュ分割も

全です。
たり、ア
認するこ
範問題の
ルの例は

の知識に
、有限空
ットアレ
結果を検
。

電 子 情報非線形光学現象は、分極 P に非線形分極 PNL 成分が生じることであり、2 光子吸収などの 3 次非線形分極は、行っており、この原理として 2 光子吸収の利用を検討しています [1] 。2 光子吸収は、3 次の非線形光学現象でチュートリアルCOMSOL Multiphysics による計算科学工学 – 波動系（5）あり、光強度が強くなると物質の光吸収係数が増大し透過光強度が低下する物理現象です。そのため、2 光す。これは、x 方向に伝搬している電界 E(x) を、エンベ（レーリー長11.972 µm）であり、出射面（右端）は PML
子吸収効率のよい物質やデバイス構造を利用すること
E1(x) と フ ェ イ ザ ーexp(-jk1x) に て、 例 え ば
ロ
ー
プ
関
数
にて終端しています。計算結果は、電解強度やビームにより、光信号の遮断・通過を光強度の強い制御光にE(x) = (0, 0, E1(x)・exp(-jk1x)) と表し、E1(x) について計算
径がガウスビームの理論式通りに変化していることをて制御できます。当研究室では、半導体を結晶成長しを行うものです。この手法により、伝搬方向のメッ
確認しました。図11より分かるように、面内のメッ
それらを半導体プロセスにて導波路形状などに加工し、と表されます。ここで、は 3 次非線形感受率テンE1(x) の変化に合わせた程度でよくなるた
シュサイズは、λ0 /n に対応して細かく設定しなければ
2シュサイズは
光子吸収特性を測定し、3次非線形感受率テンソル
ならないが、Beam envelope インターフェイスを使用し(3)
χめ、計算時間と使用メモリの大幅な低減ができます。
を求めています。しかしこれらの測定には、フェム
次に、2光子吸収は3次の非線形光学現象であるたているため、伝搬方向は5分割（30 µm ごと）でもガウス
ソル、i,j, k, l = a, b, c（結晶軸方向）です。COMSOL
ト秒レーザーをレンズで集光した光を使用するため、め、通常のマクセル方程式に非線形項を追記しなけれ
ビームの挙動が計算出来ていることが分かります。
光強度が空間的にも時間的にも変化し、さらに 2 光子
では、PNL をユーザーが独自定義し、通常の電束密度(3)ばならないが、COMSOL Multiphysics では、そのよう次に、2光子吸収係数 β = 25 cm/GW となる χ ijκl を設定
吸収を生じながらの光伝搬するため、厚膜試料に適用
L
L
NL
D
=
ε
E
+
P
を
D
=
ε
E
+
P
+
P
となるよう追記しま
なユーザー独自による式の変更などが簡単にできま（集光位
0し、入射ガウスビームのビームウェスト位置0
できる解析式がない。そのため我々は、実験データか
ε0
NL
P
=
∑ χijκl(3) (ω; −ω , ω , ω ) ⋅ E *j (ωi
P に非線形分極 PNL 成分がx
す。非線形光学現象は、分極置）に対する光透過率を計算した結果を、図12に示し1(3)す。なお、試料の座標系として結晶方位に合わせた第4 jkl
ら χ を求めるために、いくつかの仮定をいれた伝搬(3)生じることであり、2光子吸収などの3次非線形分極ます。これは、実際の Z-scan 測定に対応するものであχ ijκl∂T
モデルを考案しそれを定式化することにより求めてい2 座標軸を定義し、電界や分極などの各成分や計算に
ρc0− ∇ k∇ T = Q
は、1光束のみを使用する場合には、り、この結果を我々の解析式にてフィッティングする(3)∂tます。そこで、求めた χ の値を使ってシミュレーシ
ついてはグローバル座標系(x,
y,
z)と第
2
座標系(a,ε0
ことで β を求めています。今回、シミュレーション結 ω 2ε b, c)
NL(3)
ョンし実験結果と比較検討することで、提唱モデルやP=χ ijκl(ω ; −ω , ω , ω ) ⋅ E *j (ω , k ) Ek (ω , k ) El (ω , k )i∇×(µ −1∇×Εβ) =果をフィッティングすることで β_ﬁt を求め、設定値
4 χ(3) 値、実験方法などの妥当性の検討を行う
c2の変換を行う必要があります。
得られた jkl(3)と比較しました。図12から分かるように、シミュレー
1
2
2
χ ijκl は3次非線形感受率テンソと表されます。ここで、∂T必要があります。
以下に、我々が検討した例を紹介し４．３ 解析の事例紹介
Q=
σ E + π f ε 0ε r′′ E + π f µ0
ρc0 j,−k,∇l k=∇a,T b,=c（
Q 結晶軸方向）です。COMSOL では、
2
ション結果は我々の解析式にて非常によくフィッティ
ル、i,
ます。
∂t
ε0
NL
P
=
Ek (ωングされており、β, k ) El (ω , k )と β_ﬁt の値も相対誤差3.98 % にて
PNL をユーザーが独自定義し、通常の電束密度∑ χijκl(3) (ω; −ω , ω , ωD) ⋅ E=*j ε(ω0E, k )最初に、光吸収のない場合のガウスビームの伝搬を
i
2
4 jkl
ω
ε
1
−
L
L
NL
一致しています。この誤差の原因は解析モデルの仮定
P
D
ε
E
P
P
+
を
=
+
+
となるよう追記します。なお、シミュレートした。ガウスビームは、伝搬方向に垂直∇×(µ ∇×2Ε光子吸収の計算手法
0) =
４．２
∂T
c2
ρc0
− ∇k∇T = Qに含まれる誤差と考えられるが、今回の結果は、現実試料の座標系として結晶方位に合わせた第2座標軸を定光 と 物 質 の 相 互 作 用 や 物∂質ルt 中 の 光 伝 搬 は マ ク セな面内の電界強度（光強度）分布がガウシアンになっ12
22
的な β 値に対して、厚膜試料における光伝搬モデルの義し、電界や分極などの各成分や計算についてはグ
2Q=
σ E + π f ε 0ε r′′ E + π−1f µ0 µr′′ωH
方程式にて記述できるが、今回の場合、主に2 つの注
εている光であり、レーリー長と呼ばれる長さを指標と∇×(µ ∇×Ε ) =2
仮定および導出された解析式は概ね正しいことを示し
2 c) の変換を行
y,z)b,
ローバル座標系(x,と第2座標系(a,c
意すべき点があります。1 つは光伝搬距離の長さであていると考えられます。
した伝搬距離に伴いそのビーム径やピーク強度が変化
う必要があります。
1
2
2
2
り、もう 1 つは非線形光学の記述です。
Q=
σ E + π f ε 0ε r′′ E + π f µ0 µr′′ H
最後に、InP 基板において (a)[100] 方向、(b)[110] 方向
2
光領域のシミュレーション、特に光伝搬に関するシします。実験で使用する凸レンズによる集光レーザーにガウスビームを伝搬させ、入射光の直線偏光方向を04.3 解析の事例紹介
ミュレーションで問題となるのが、屈折率n の物質中はこのような光になっています。図1 に、計算に使用
度から90度まで変えたときの、β の偏光方向依存性の結最初に、光吸収のない場合のガウスビームの伝搬をでの波長 λ = λ0 / n (~0.1－0.5 µm)とモデル系サイズ L
果を、図13に示します。測定結果は実際の試料で測定した(a)モデル（メッシュ）と(b)計算結果を示します。
シミュレートしました。ガウスビームは、伝搬方向に（>~100 µm）が大きく異なることです。光の位相情報
垂直な面内の電界強度（光強度）分布がガウシアンに
波長 λされた結果であり、シミュレーションは別実験にて研
屈折率 n = 3.14、試料長 L = 150 µm、を考慮して計算するには、有限要素法では各場所での
0 = 1.64 µm、究室で測定された χ(3) の値 [11] を使用して計算した結果なっている光であり、レーリー長と呼ばれる長さを指波長の1/6 以下にメッシュサイズを設定しなければ正試料径44 µmϕ、スポット径[12]
0 = 5 µmϕ（レーリー長β の理論式に今回の実
です。β
による2w
D は、Dvorak ら標とした伝搬距離に応じてそのビーム径やピーク強度
しく計算出来ないが、それに比べてL が大きすぎるた(3)[11]
χ の値 を適用した結果です。いずれ
11.972験条件およびµm）であり、出射面（右端）はPML にて終端が変化します。実験で使用する凸レンズによる集光
め、必要となるメモリ量が膨大になります。しかしなの結果も、調整パラメータは何もなく結果をそのまま
レーザーはこのような光になっています。図11に、計しています。計算結果は、電解強度やビーム径がガウ
がら、今回対象としているような系では、光強度の変プロットしたものです。どちらの伝搬方向でも、シ算に使用した(a)モデル（メッシュ）と(b)
計算結果を示化は伝搬方向に対しては非常に緩やかです。この様なスビームの理論式通りに変化していることを確認しま
ミュレーション結果と理論式による βD、実際の測定結
λ0 = 1.64 µm、屈折率 n = 3.14、試料長 L =します。波長場合には、COMSOL Multiphysics の波動光学モジューした。図1 より分かるように、面内のメッシュサイズ
果が非常によく一致していることが分かります。
150 µm、 試 料 径44 µmφ、 ス ポ ッ ト 径2w0 = 5 µmφ∑は、λ0 /n に対応して細かく設定しなければならないが、術の 進 展 ,

lems/

への適用
(a)シミュレーションモデル（メッシュ）

高速化・
以下の超
我々、千
の応答速
光信号を
の研究を
(b)シミュレーション結果Vol.22, No.2 2017(??-7)

図1

図11

ガウスビームの伝搬シミュレーション例

ガウスビームの伝搬シミュレーション例

計算工学

(??-8
図12
1

（20-8）

1

透過率のビームウェスト位置依存線（Z-scan シミュ

図 2 レーション）
透過率のビームウェスト位置依存線(Z-scan シミュレーション)
Vol.23, No.2 2018）です。COMSOLにガウスビームを伝搬させ、入射光の直線偏光方向を、通常の電束密度0 度から 90 度まで変えたときの、β の偏光方向依存性
なるよう追記しまの結果を、図 3 に示します。測定結果は実際の試料で方位に合わせた第測定された結果であり、シミュレーションは別実験にの各成分や計算に チュートリアルて研究室で測定された χ(3) の値 [2] を使用して計算した
第 2 座標系(a, b, c)
結果です。βD は、Dvorak ら [3] による β の理論式に今回

■参考文献
[1]

スビームの伝搬を、伝搬方向に垂直
ガウシアンになっれる長さを指標とピーク強度が変化
よる集光レーザー
1 に、計算に使用結果を示します。
料長 L = 150 µm、µmϕ（レーリー長は PML にて終端やビーム径がガウることを確認しまのメッシュサイズればならないが、

位置依存線
ン)

COMSOL Multiphysics による計算科学工学 – 波動系（5）

(a)

(b)

(??-8)

図 3図13
2 光子吸収係数
β の偏光方向依存性。
2光子吸収係数βの偏光方向依存性。
- Vol.22, No.2 2017(a)[100] 方向伝搬、(b)[110] 方向伝搬(a)[100]方向伝搬，(b)[10]方向伝搬4.4 まとめ
簡単ではあるが、COMSOL Multiphysics による2光子吸収を生じる場合の光伝搬シミュレーションの方法とそれらの結果を紹介しました。これらの結果は、実際の測定結果や理論式と非常に良く一致することを紹介しました。これらは、シミュレーション方法や実験方法、得られた χ(3) の値、仮定した伝搬状態や導出した解析式など、いずれもが矛盾ないことを意味していると考えられます。今後は、導波路構造などデバイス形状
にてシミュレーションを行い、デバイス設計などへの利用も期待されます。

計算工学

https://www.comsol.jp/blogs/how-to-model-gearboxvibration-and-noise-in-comsol-multiphysics/[2] http://www.aes.org/conferences/2017/automotive/downloads/
AES_Modeling_Car_Cabin_Acoustics_2017.pdf[3] 経済産業省資源エネルギー庁、平成25年度エネルギーに関する年次報告（エネルギー白書2014）[4] Y. Wada, S. Fujii, E. Suzuki, M. M. Maitani, S. Tsubaki, S.
Chonan, M. Fukui, & N. Inazu, Smelting Magnesium Metal
using a Microwave Pidgen Method, Scientific Reports, 7,46512; doi:10.1038/srep46512(2017)[5] S. Fujii, S. Tsubaki, E. Suzuki, S. Chonan M. Fukui, and Y.
Wada, "Study on Metal Reﬁning Process of Sc metal using by
Microwave Irradiation", Progress in Electromagnetic Research
Symposium (PIERS), pp. 4203-4203 (2016)[6] 平野拓一、 電磁界シミュレータ利用の勘所 、電子情報通信学会誌、小特集 電磁界シミュレーション技術の進展、Vol.100, No.5, pp.342-348, 2017年5月[7] http://www.ieice.org/es/est/activities/canonical_problems/[8] http://www.takuichi.net/em_analysis/canonical/[9] T. Matsusue, H. Bando, S. Fujita and Y. Takayama,
Polarization dependence of two-photon absorption coefﬁcientand nonlinear susceptibility tensor in InP , Phys. Status Solidi
C, 7, pp. 387-389 (2011).
[10] M. Oishi, T. Shinozaki, H. Hara, K. Yamamoto, T. Matsusueand H. Bando, ''Measurement of polarization dependence of
two-photon absorption coefficient in InP using extended
Z-scan technique for thick materials'', Jpn. J. Appl. Phys., 57,
030306, (2018); doi: 10.7567/JJAP.57.030306[11] 原皓、坂東弘之、篠崎智文、大石真樹、松末俊夫、「1光束楕円偏光を用いた Z-scan 法による3次非線形感受率テ
ンソルの虚部の測定」、第75回応用物理学会春季学術講
演会、19a-A26-10、(2014).
[12] M. D. Dvorak et al., "Measurement of the anisotropy of twophoton absorption coefﬁcients in zincblende semiconductors",
IEEE J. Quantum Electron, 30, pp. 256-268 (1994).

（20-9）

Vol.23, No.2 2018