3783＜チュートリアル＞ 汎用マルチフィジックス有限要素解析ソフトウェア COMSOL Multiphysics を取り上げて新し

汎 用 マ ル チ フ ィ ジ ッ ク ス 有 限い計算科学工学の仕事の流れを解説していきます。今回は拡散、固有値、最適化を中心に解
要 素 解 析 ソ フ ト ウ ェ ア C O M S O L M u l t i p h ys i c s を 取 り 上 げ て 新 しチュートリアルい 計 算 科 学 工 学 の 仕 事 の 流 れ説していただきます。なお、チュートリアル記事は1ページ目のみを本誌掲載し、続きは日を 解 説 し て い き ま す 。今 回 は 拡 散 、固 有 値 、最 適 化 を 中 心 に 解 説

HP
し て い た だ き ま す 。 な お 、 前 本計算工学会号よりチュー
ト上で公開していますので、そちらも併せてご参照ください。
リ ア ル 記 事 は 1 ペ ー ジ 目 の み を 本 誌 掲 載 し 、続 き＜チュートリアル＞

は 日 本 計 算 工 学 会 HP 上 で 公 開 し て い ま す の で 、 そ ち ら も 併 せ て ご 参 照 く だ さ い 。

汎 用 マ ル チ フ ィ ジ ッ ク ス 有 限 要 素 解 析 ソ フ ト ウ ェ ア COMSOL Multiphysics を 取 り 上 げ て 新 しい 計 算 科 学 工 学 の 仕 事 の 流 れ を 解 説 し て い き ま す 。今 回 は 拡 散 、固 有 値 、最 適 化 を 中 心 に 解 説
し て い た だ き ま す 。な お 、前 号 よ り チ ュ ー ト リ ア ル 記 事 は 1 ペ ー ジ 目 の み を 本 誌 掲 載 し 、続 きCOMSOL Multipysics による計算科学工学 –拡散、固有値、最適化（６）は 日 本 計 算 工 学 会 HP 上 で 公 開 し て い ま す の で 、 そ ち ら も 併 せ て ご 参 照 く だ さ い 。

COMSOL Multiphysics による計算科学工学

Computational Science and Engineering by using COMSOL Multiphysics–Diffusion, Eigenvalue and Optimization (6)COMSOL Multipysics による計算科学工学 –拡散、固有値、最適化（６）– 拡散、固有値、最適化（6）Computational Science and Engineering by using COMSOL Multiphysics橋口真宜–Diffusion, Eigenvalue and Optimization (6)橋口 真宜

積内に流入してくる流束（フラックス）ベクトルの法線成分を

１．拡散
今まで様々な分野での計算科学工学を見てきました．そこで
1 拡散

橋口真宜
境界で積分することで，拡散をモデル化します．検査体積へのベクトルには、熱に関するフーリエ則、化学種に関する流束の法線方向の出入りはガウスの発散定理でdiv へ変換されは現象の数学モデルを偏微分方程式の初期値・境界値問題によ今まで様々な分野での計算科学工学を見てきました。
フィック則、流体におけるニュートン粘性則などがあり
１．拡散って表現し，解析的あるいは数値的に解を求める手法が使われ積内に流入してくる流束（フラックス）ベクトルの法線成分をます．良く出てくる流束ベクトルには，熱に関するフーリエ則，

そこでは現象の数学モデルを偏微分方程式の初期値・ます。上述のミクロな粒子の挙動は熱伝導係数、拡散係境界で積分することで，拡散をモデル化します．検査体積への今まで様々な分野での計算科学工学を見てきました．そこで化学種に関するフィック則，流体におけるニュートン粘性則なます．それらの基礎方程式には頻繁に div ( -k grad・) といっ
流束の法線方向の出入りはガウスの発散定理でdiv へ変換されは現象の数学モデルを偏微分方程式の初期値・境界値問題によ
境界値問題によって表現し、解析的あるいは数値的に数、粘性係数といった輸送係数に繰り込まれています。
どがあります．上述のミクロな粒子の挙動は熱伝導係数，拡散た形の拡散項が出てきます。
grad の係数が一定の場合にはラプます．良く出てくる流束ベクトルには，熱に関するフーリエ則，って表現し，解析的あるいは数値的に解を求める手法が使われ解を求める手法が使われます。それらの基礎方程式にさて、マックス・プランクは溶鉱炉で溶けた鉄の温係数，粘性係数といった輸送係数に繰り込まれています．
ラシアン∇2 で表現されます．
実際に拡散を受ける対象は多数の化学種に関するフィック則，流体におけるニュートン粘性則なます．それらの基礎方程式には頻繁にdiv ( -k grad・)
といっ
は頻繁に
div
(
-k
grad・)
といった形の拡散項が出てきま
度が一体何度であるかをその輻射光から知る方法につさて，マックス・プランクは溶鉱炉で溶けた鉄の温度が一体
原子，電子，分子，イオンといった粒子群になり，拡散の具合
どがあります．上述のミクロな粒子の挙動は熱伝導係数，拡散た形の拡散項が出てきます。grad の係数が一定の場合にはラプ
2す。gradで表いて研究をしていました。ウィーンの放射則は知られ2 で表現されます．実際に拡散を受ける対象は多数の何度であるかをその輻射光から知る方法について研究をしていはそれらの運動をラグランジュ的に追跡することで把握できま
係数，粘性係数といった輸送係数に繰り込まれています．
ラシアン∇の係数が一定の場合にはラプラシアン∇
現されます。実際に拡散を受ける対象は多数の原子、ていましたが、可視光域では実験に合うが、赤外域でました．ウィーンの放射則は知られていましたが，可視光域でさて，マックス・プランクは溶鉱炉で溶けた鉄の温度が一体す．良く知られた例として，水に浮かべた花粉の動きがありま
原子，電子，分子，イオンといった粒子群になり，拡散の具合何度であるかをその輻射光から知る方法について研究をしていはそれらの運動をラグランジュ的に追跡することで把握できまは実験に合うが，赤外域では合わない．一方，レイリー・ジー
す．これはブラウン運動と呼ばれており，花粉の見せる動きは
電子、分子、イオンといった粒子群になり、拡散の具合は合わない。一方、レイリー・ジーンズの放射則ではました．ウィーンの放射則は知られていましたが，可視光域です．良く知られた例として，水に浮かべた花粉の動きがありま
ンズの放射則では赤外域では合うが，振動数が高くなると破た周囲の分子運動を反映した酔歩運動です．当時は原子・分子のはそれらの運動をラグランジュ的に追跡することで把握
赤外域では合うが、振動数が高くなると破たんする。
は実験に合うが，赤外域では合わない．一方，レイリー・ジー
す．これはブラウン運動と呼ばれており，花粉の見せる動きはんする．プランクはエネルギーの最小単位を考えつき，３年間
存在の確証がなかったのですが，ブラウン運動によって原子・できます。良く知られた例として、水に浮かべた花粉のプランクはエネルギーの最小単位を考えつき、3年間考ンズの放射則では赤外域では合うが，振動数が高くなると破た周囲の分子運動を反映した酔歩運動です．当時は原子・分子の考え抜いた末に，プランク分布を使えば全振動数領域において分子の存在を知ることができることをアインシュタインが示し動きがあります。これはブラウン運動と呼ばれており、え抜いた末に、プランク分布を使えば全振動数領域にんする．プランクはエネルギーの最小単位を考えつき，３年間
存在の確証がなかったのですが，ブラウン運動によって原子・
実験値ときわめてよく一致することを示しました[1]．これが量ました．酔歩運動では速度が不連続に変化するので微分不可能
花粉の見せる動きは周囲の分子運動を反映した酔歩運
おいて実験値ときわめてよく一致することを示しまし考え抜いた末に，プランク分布を使えば全振動数領域において分子の存在を知ることができることをアインシュタインが示し
子力学の幕開けとなりました．量子力学ではとびとびのエネル
であり，通常の数学は使えないので，確率過程（stochastic
実験値ときわめてよく一致することを示しました[1]．これが量ました．酔歩運動では速度が不連続に変化するので微分不可能
動です。当時は原子・分子の存在の確証がなかったのた [1]。これが量子力学の幕開けとなりました。量子力
ギー準位を考えることの他に，物体の粒子性と波動性の二重性
process）という特別の取り扱いをする必要があります．
子力学の幕開けとなりました．量子力学ではとびとびのエネル
であり，通常の数学は使えないので，確率過程（stochastic
ですが、ブラウン運動によって原子・分子の存在を知る
学ではとびとびのエネルギー準位を考えることの他を考慮する必要があります．アインシュタインは波であると考一方で，それらの粒子群の挙動を偏微分方程式で現象論的にギー準位を考えることの他に，物体の粒子性と波動性の二重性
process）という特別の取り扱いをする必要があります．
ことができることをアインシュタインが示しました。酔に、物体の粒子性と波動性の二重性を考慮する必要がえられていた光の粒子性を示し，ド・ブロイはその逆である物扱う場合には，空間に検査体積を設定し，その境界から検査体を考慮する必要があります．アインシュタインは波であると考一方で，それらの粒子群の挙動を偏微分方程式で現象論的に歩運動では速度が不連続に変化するので微分不可能であります。アインシュタインは波であると考えられて質波（ド・ブロイ波）の概念を提唱，それを受けて，シュレー
えられていた光の粒子性を示し，ド・ブロイはその逆である物扱う場合には，空間に検査体積を設定し，その境界から検査体あり、通常の数学は使えないので、確率過程（stochastic
いた光の粒子性を示し、ド・ブロイはその逆である物
ディンガーがシュレーディンガ－方程式を考案して波動力学を
筆者紹介
質波（ド・ブロイ波）の概念を提唱，それを受けて，シュレー
構築，ハイゼンベルグの行列力学との等価性も証明，といったディンガーがシュレーディンガ－方程式を考案して波動力学を筆者紹介 という特別の取り扱いをする必要があります。
process）
質波（ド・ブロイ波）の概念を提唱、それを受けて、構築，ハイゼンベルグの行列力学との等価性も証明，といったことを経て，現在ではシュレーディンガー方程式が各方面で成一方で、それらの粒子群の挙動を偏微分方程式で現象
シュレーディンガーがシュレーディンガ - 方程式を考案
ことを経て，現在ではシュレーディンガー方程式が各方面で成果を挙げています．シュレーディンガーは彼の著書「生命とは論的に扱う場合には、空間に検査体積を設定し、その境して波動力学を構築、ハイゼンベルグの行列力学との果を挙げています．シュレーディンガーは彼の著書「生命とは何か」で分子生物学への道を開きました．
界から検査体積内に流入してくる流束（フラックス）ベク
等価性も証明、といったことを経て、現在ではシュ
何か」で分子生物学への道を開きました．
シュレーディンガー方程式は次の形をしています．なんと，
トルの法線成分を境界で積分することで、拡散をモデル
レーディンガー方程式が各方面で成果を挙げていま
シュレーディンガー方程式は次の形をしています．なんと，
この式にも拡散項が出てきます．
この式にも拡散項が出てきます．
化します。検査体積への流束の法線方向の出入りはガウ
す。シュレーディンガーは彼の著書「生命とは何か」でスの発散定理で div へ変換されます。良く出てくる流束
分子生物学への道を開きました。
シュレーディンガー方程式は次の形をしています。
はしぐち まさのり Masanori HASHIGUCHI
はしぐち まさのり Masanori HASHIGUCHI
なんと、この式にも拡散項が出てきます。
宮崎大学工学部卒，同大学院機械工学専攻修了，九州大学工学
宮崎大学工学部卒，同大学院機械工学専攻修了，九州大学工学
筆者紹介
部航空工学教室助手，マツダ（株）研究開発本部空力実研課課
部航空工学教室助手，マツダ（株）研究開発本部空力実研課課∂ψ   2∂ψ 2   2 2=  − i  ∇ =+V−ψ ∇ + V ψi長，計算流体力学研究所部長，流体物理研究所部長，（財）濱野
長，計算流体力学研究所部長，流体物理研究所部長，（財）濱野
はしぐち まさのり∂t 2m∂ t  2m

生命科学研究財団事務局長を経て，現在，計測エンジニアリン
生命科学研究財団事務局長を経て，現在，計測エンジニアリン
宮崎大学工学部卒、同大学院機械工学専攻修了、グシステム（株）第１技術部部長．元宮崎大学地域共同研究セグシステム（株）第１技術部部長．元宮崎大学地域共同研究セ h / 2π= h / 2π =九州大学工学部航空工学教室助手、マツダ（株）研ンター客員教授，元法政大学非常勤講師．日本機械学会認定計
続きは Web でンター客員教授，元法政大学非常勤講師．日本機械学会認定計
続きは Web で究開発本部空力実研課課長、計算流体力学研究所
続きは Web で算力学技術者上級アナリスト．日本機械学会畠山賞受賞．
日本計算工学会誌「計算工学(Vol.23,
算力学技術者上級アナリスト．日本機械学会畠山賞受賞．
日本計算工学会誌「計算工学(Vol.23,
No.3)」HP: No.3)」HP:
部長、流体物理研究所部長、（財）
濱野生命科学研

日本計算工学会誌「計算工学（Vol.23, No.3）」HP:
http://www.jsces.org/activity/journal/

趣味：トランペット 趣味：トランペット究財団事務局長を経て、現在、計測エンジニアリ

ングシステム（株）第1技術部部長。元宮崎大学地域共同研究センター客員教授、元法政大学非常勤講師。日本機械学会認定計算力学技術者上級アナリスト。日本機械学会畠山賞受賞。
趣味：トランペット

Vol.23, No.3 2018

注．COMSOL, COMSOL Multiphysics, COMSOL Desktopおよび Livelinkは
COMSOL ABの商標または登録商標です。その他の製品名、サービ
ス名、組織名は、各組織の商標または登録商標です。なお、本文中に™、®マークは明記しておりません。本稿の説明は著者独自の考えに基づいており他を代表するものではありません。

（19）

チュートリアルCOMSOL Multiphysics による計算科学工学 – 拡散、固有値、最適化（6）

それでは、定常状態におけるこの方程式の有限要素解を求めてみましょう。（教科書は鈴木 [2] を参照）
問題の設定は図1にあるように、質量 m の粒子が一次元空間に設定された無限井戸型ポテンシャルの中で運動する場合の波動関数を計算します。
数値解析においては Veriﬁcation が重要です。そのために、解析解を導出してみます。井戸の中ではポテン
シャルは0であり、次式を解きます。

−

2 2∇ ψ = Eψ2m

この式は左辺の演算子が未知関数ψに作用した結果が自身の E 倍になっていることを意味しており、固有値問題になります。
さて、解きやすいように式を変形します。

図1∇ 2ψ + λ 2ψ = 0

無限井戸型ポテンシャル中の量子の運動

ここで、

λ2 =

2mE
2

という関係を設定しています。すると、解は次のように記述できます。

ψ ( x) = A sin λ x + B cos λ x
境界条件としては、次のディリクレ条件になります。これは未知関数の連続性を確保するためです。

ψ (0) = ψ (a) = 0
この条件を満たすことから次式を得ます。

B = 0, λ a = nπ係数 A は規格化条件から決めることができます。
よって、具体的な解の形を得ることができます。

ψ n ( x) =
En =

2
nπ x
sin
a
a

 2 n 2π 22m a 2

ここで、n は n=1,2,3,…という離散的な数値しか許されません。
つまり、この問題のように、束縛された状態にある量子の波動関数は、離散的な固有値、それに応じた離散的なエネルギー準位と波形（モード）を持つことが示されました。
では、COMSOL Multiphysics の係数形式 PDE（偏微分方程式）インターフェースを利用して、この問題の有限要素解を求めてみます。
COMSOL Multiphysics を 起 動 し ま す。 モ デ ル ウ ィ
ザードで空間1D を選択、次のフィジックス選択に進
み、 数 学 の と こ ろ へ い き、 図2の リ ス ト で 係 数 形 式
PDE(c) をクリックし、追加ボタンをクリックします。
スタディに進み、固有値を選択後、完了ボタンを押します。これでプロトタイプがモデルビルダーに構成さ
計算工学

図2

PDE インターフェース

れます。デフォルト設定でほとんど完成しています。
それに追加・修正する形で操作をしていくことになります。
まず、図3に示すグローバル定義の下のパラメタで、
ディラックの定数 hD（プランク定数 hP を2π で割り算したもの）、質量 m を hD^2/2、井戸の幅 a を pi（π を意味する）を入力します。ここで使われている _const の形の変数は COMSOL Multiphysics に用意されている物理定数です。Ctrl+ スペースを押すと、画面に一覧表が表示されるので便利です。

図3

グローバル定義の下のパラメタ設定

ジオメトリ（図4）で間隔を選択し、右端点に a と入力
し、1D 空間を作成します。COMSOL Multiphysics はこのようにパラメタでジオメトリ表現ができるので、大きさを変更するときなどに便利です。
ではシュレーディンガー方程式の設定に入ります。
係数形式 PDE(c) の下の係数形式 PDE1をクリックし、
設定ウィンドウで、方程式セクションを開くと、方程（19-2）

Vol.23, No.3 2018

チュートリアルCOMSOL Multiphysics による計算科学工学 – 拡散、固有値、最適化（6）図4

図5波動関数に相当する変数 u を二乗したもの（conj は共役複素数）を空間で積分（intop1）した結果が1になるように係数 A を決めています。それを使って規格化された数値解が uN です。
続いて、有限要素解析の準備をします。離散化する
ためのメッシュの作成です。メッシュを右クリック
し、エッジを選択します。そのエッジを選択し、分布を追加します。分布はメッシュの粗密を制御できます。ここでは64個の等分割とします。先ほど導出した解析解ともしも一致しない場合にはこのところを調整するといった作業が生じます。
図7に示すように、スタディ1の下にステップ：固有値があるのでそこをクリックして、設定ウィンドウにいきます。

ジオメトリの作成

固有値計算用のステップ(b)

固有値探索部分の設定例

係数形式 PDE の設定画面

式のテンプレートが表示されます。先ほど解きやすい形に変更した式を、係数を図5のように入力します。
これでシュレーディンガー方程式の入力が完成しました。
次は、境界条件です。係数形式 PDE(c) を右クリック
して、コンテクストメニューを表示させ、そこからディリクレ境界条件を選択します。グラフィックウィ
ンドウに表示されている一次元空間の両端点をマウスで選択します。デフォルトで関数値を0に設定する内容になっているのでこれ以上の作業は不要です。
ここまでで、偏微分方程式の境界値問題の定義が終わりました。
固有値問題の解はモード形状のみに意味があり、その絶対値はそのままでは決まらないのは周知のとおりです。そこで規格化をします。それには、コンポーネント1の下の定義を右クリックし、コンポーネントカッ
プリングを選択し、積分演算子（intop1( )）を1次元空間に割り当てます。同じ定義を右クリックし、変数を選択後、図6のように変数を定義します。

図6
計算工学

(a)

図7

スタディ1の内容

スタディ1を右クリックし、計算をします。図8に計算結果を示します。図8(a) は固有値の結果です。ここでは a を π に選んでいるので解析解は n です。図8(b) は波動関数の結果で、解析解は実線で、COMSOL の有限
要素解は〇印で示しています。いずれも両者はよく一致しています。結果として、メッシュの調整は不要としました。
数値解と解析解を比較するための描画を行う際に便利な方法を紹介します。固有値解析では複数個の固有値が一つのデータセット（結果の下に表示）に格納されます。それらを図8のように一つのグラフとしてプロットすることを考えます。そのためにはデータセットに格納されている固有値の番号を指定することになります。そのために、ここでは次のような関数を利用しています。いずれも1D プロットグループのライングラフ
で使用しました。

解の規格化（19-3）

Vol.23, No.3 2018

チュートリアルCOMSOL Multiphysics による計算科学工学 – 拡散、固有値、最適化（6）図9(a)波動関数による期待値の計算式の設定例

図9には、座標、座標の2乗、運動量、運動量の2乗の各期待値を算出する場合の式の例を示しています。
続いて、この定常解を使って、波動関数の時間変化を調べてみます。説明をしやすくするために、解析解を使って話をします。すでに説明した定常状態の波動関数
は固有関数であり、互いに直交します。そして、任意の状態を表す波動関数は定常状態の波動関数の線形結合で表現されます（任意のベクトルは基底ベクトルの線形結合で表現できます）。
いま、初期状態が次式で与えられているとします。

固有値の計算値（解析解 n と一致）

φ0 ( x) =2
 3π x 
πx 
sin 
 cos a
 2a 
 2a 

すると初期時刻 t=0では次式が成り立ちます。

φ0 ( x) = c1ψ 1 ( x) + c2ψ 2 ( x) + c3ψ 3 ( x) +･･･(b)

未定係数 c を決めるには、固有関数同士の直交性を利用します。

波動関数の計算値（解析解：実線）
図8a∫ ψ ( x)ψ ( x)dx = δ

有限要素解と解析解の比較

0

解析解については、

∫

0

m,n

1  πx
2π x  2
mπ x
+ sin
sin
dx = cm
 sin

a
a
a  a
a

この式では三角公式を利用してφ0を展開しています。この式で m を1,2,3と順に変えていくと、m が3以上では cm は
全て0になります。m が1と2では cm=1/ √2 に決まります。
以上のことから、t>0での波動関数は、

Ψ ( x,t ) =11
ψ 1 ( x) exp(−iE1 t ) + ψ 2 ( x) exp( −iE2 t )
2
2となります。
では、これをもとに COMSOL で波動関数の時間発展を可視化してみましょう。先ほどまで利用してきたモデル
ビルダの内容に追加をしていきます。
コンポーネント1の下の定義をクリックし、変数を追加します。そこに図10の内容を記述します。

d 2) ψ dxdx

この被積分関数の中には波動関数の空間2階微分が出てきます。もしも有限要素解を1次の形状関数で計算した場合には2階微分値は0になってしまいます。そのような場合に は 離 散 化 次 数 を 変 更 し た う え で 計 算 を 行 い ま す。
COMSOL は離散化セクションで形状関数の次数を変更で
きます。
計算工学

a

0

としました。他の n についてはこれらを複製後、n を各々書き換えれば図8(b) のプロットが完成します。
COMSOL Multiphysics の半導体モジュールはシュレーディンガー方程式のインターフェースが用意されています。本格的な計算を行うにはそちらを利用するのも良いでしょう。
波動関数を利用して、運動量の2乗の期待値を計算することを考えてみます。
理論式は以下の通りです。
an

ここでδm, n はクロネッカーのδです。
初期時刻での関係式の両辺に m 番目の波動関数を掛け算し、空間で積分します（ベクトル間の内積計算と同じ）。

数値解については、< p 2 >= ∫ ψ * (−im（19-4）

図10

波動関数の時間発展を記述する変数の定義
Vol.23, No.3 2018

チュートリアルCOMSOL Multiphysics による計算科学工学 – 拡散、固有値、最適化（6）intop1(conj(wave)*x*wave)を記述し、解を更新することで結果を参照できます。
解析解は次式で与えられます。

< x >=
図11

時間方向の可視化をするための定数の設定

グローバル定義の下のパラメタに図11を追加します。
モデルビルダの root を右クリックしてスタディ追加を選
択し、時間依存をスタディ2としてモデルビルダに加えます。スタディ2のステップ1：時間依存で図12の内容を設定します。COMSOL は解かないということも設定できます。ここでは係数形式 PDE は解かないので図12のように計算対象のチェックを外しておきます。この機能はマル

チフィジックスを連成したり、連成を切り離したりする際にも大変便利です。

図12

時間依存ステップの設定内容

それではスタディ2を右クリックして計算を実行します。計算結果は結果のデータセットの箇所にスタディ
2/ 解2(sol2) として保存されます。
計算が終了すると、デフォルトで1D プロットグルー
プ2が追加されます。その設定ウィンドウの式に waveを設定します。そのライングラフを右クリックして複製を行います。複製されたライングラフに先ほどの初期値の式を記述します。プロットをクリックすると、
図13のプロットを出力できます。

両者は図14のように一致します。

図14

波動関数の時間変化（初期分布：破線）

最後にこの波動関数によって粒子の位置座標の期待値を計算してみます。
図10の下に変数名を xE として、式に計算工学

粒子位置の期待値の時間変化（実線：解析解）

粒子はポテンシャルの中で往復運動をしています。
ここでは無限井戸型ポテンシャルのある粒子の波動関数を見てきました。この束縛状態ではエネルギーはとびとびの値を持ちますが、ポテンシャルの大きさによってはトンネル効果が生じたり、階段型ポテンシャ
ルでの散乱では連続的なエネルギーを持つことなど興味ある現象を数値的に解析できます。最近注目されているメタマテリアルなどでは周期構造を利用します。
そのような系ではエネルギーの連続的な状態がある幅を も っ て 存 在 す る バ ン ド 構 造 が 重 要 な 事 項 で す。
COMSOL Multiphysics には半導体モジュールがあり、その中にシュレーディンガー方程式を解くためのインターフェースがすでに用意されています。興味をもたれた読者は試してみるとよいでしょう [3]。
シュレーディンガー方程式は、粒子の波動関数に関する空間2階の拡散項をもち、時間に関する1階の発展方程式です。この形は生物や化学反応系でみてきたものに近いのですが、シュレーディンガー方程式は時間微分項の係数に虚数が入っていることが大きな特徴です。つまり、粒子の波動関数は複素数であり、振幅と位相を持っています。位相が激しく変化する箇所では粒子同士の打消しが激しく生じ、位相がゆっくり変化
する箇所では打消しが生じにくいので、位相の変化がゆっくりした箇所を運動軌跡として選択していくように見えることになります [4]。

2

図13

a 16a
 E − E1 
− 2 cos  2
t
2 9π
 


固有値

シュレーディンガー方程式は固有値問題を構成することを見てきました。このチュートリアルの連載の中で見てきた電磁波や光も境界条件を決める際に固有値問題を利用していました。
今度は固体力学を見ていきます。COMSOL Multiphysics
は構造力学モジュールの中に、ソリッドモデルに加えてシェルモデル（2次元ではプレートと称している）を用意しています。板厚方向にメッシュを作成する必要が無いという利点があります。
こ こ で は、 シ ェ ル の 振 動 問 題 を 取 り 上 げ ま す。
COMSOL は単位系は SI 単位系を利用していますが、利用者が SI 以外でも単位に [ ] を付ければ自動的に SI に

（19-5）

Vol.23, No.3 2018

チュートリアルCOMSOL Multiphysics による計算科学工学 – 拡散、固有値、最適化（6）

変換します。単位の整合性もチェックしてくれます。
単位は数式のように記述します。
厚 み1[mm]、 ヤ ン グ 率72[Gpa]、 密 度2700[kg/m^3]、ポアソン比0.35の正方形板に対して固有周波数解析をしてみました。
この問題の理論的な固有角周波数は次式で与えられ
ています。
1/2

 D   m2 n2 
ωmn = π    2 + 2  ρ h   am an 
2

D=

Eh3
12(1 − v 2 )このモデルは読者は COMSOL 社のホームページからダウンロードできます。
モデルビルダ画面の表示が英語・日本語併記になっています。COMSOL に搭載されているアプリケーショ
ンライブラリは英語で記述されており、そこからファイルを開くと英語表記になっていますが読者が日本語環境で使用している場合には、モデルビルダのモデル
ツリーノードテキストボタンをクリックし、タイプを選択することでこのように英語・日本語併記ができます。もちろん、他の設定ウィンドウの内容はタイプを選択しなくても日本語になっています。

この理論解と比較して数値計算の Veriﬁcation を行います。単純支持の結果は図15に示す通りです。

図15

固体力学における単純支持板の固有振動解析

このように COMSOL の有限要素解は理論解と一致しています。長方形板の全周固定、短辺支持の計算も行い、それらも理論解とよく一致することを確認しました。
それでは、このシェルと圧力音響の連成した固有周波数解析を行なってみます。圧力音響の計算の正しさはすでに見てきたように確認しています。
取り上げる問題は、空気の入った半径38[mm]、高さ255[mm] の円筒容器の底面が振動版であるとしてシェルでモデル化し、その連成を解析します。図16にはモ
デル図、それを解析するための COMSOL モデルビルダ
画面を図17に示しています。

図16
計算工学

図17

シェル－圧力音響連成時の固有周波数解析

この解析をするにあたって、スタディで計算対象の選択を切り替えることで、シェルのみ、音響のみ、連成時というスタディを簡単に行うことができます。
結果を図18に示します。半解析式および実験と比較していますが、COMSOL 解はよく一致しています。
続いて、周波数応答解析を行なった結果を図19に示

底面に振動板を有する円筒容器

図18（19-6）COMSOL の数値解と他の比較
Vol.23, No.3 2018

チュートリアルCOMSOL Multiphysics による計算科学工学 – 拡散、固有値、最適化（6）図19

周波数応答解析結果

します。この場合、底面に面荷重を与えることで加振しています。
音響単独の固有周波数がある2箇所においては、音響− シ ェ ル 連 成 に よ っ て（637.1, 707.6）お よ び（2623,
2697）の新しい共振が出現していることがわかりました。このように、固有周波数解析（固有値解析）を行っておき、それらが連成によってどのように影響を受けるかは周波数応答解析を行うことになりますが、共振を見出すには、分析する周波数の刻み幅を小さくしな
いと見落としてしまいます。一方で、小さな刻み幅は
計算時間の長大化を招いてしまいます。
そこで、振動解析と固有値の関係を分析してみます。
時間依存の振動解析が次式で行われるとします。
M
M
2

2

2

+
2 +
2

+

図20

振の例です。周波数応答解析に2分かかっていた計算がモーダル解析では11秒に短縮されました。
固有値に対応する固有モードの直交性を活かすことでこのような解析が可能になっています。

3

2

+

=0
0
=

−1

=0

−1

=0

となります。P は固有ベクトルであるので、次式を満
たします。
−1

= Λ

ここで、Λは対角行列。
左から P-1を掛け、対角行列Λとの関係
−1 −1
Λ
Λ=
=
を用いることで次式を得ます。
2

2

+Λ =0

つまり、新しい解ベクトル v はその成分しか含まない方
程式に分解され、連立方程式を解く必要が無いので、短時間に容易に解を求めることができます。周波数領域では調和振動を仮定して、周波数領域で代数的に解きます。角周波数をωとすると、上式の場合では、−

2+Λ = 0を解くことになります。
COMSOL はモーダル解析を行なえるようにしています。図20に結果を示します。3次元ソリッドモデルの加
計算工学

最適化

数学モデルが完成し、Veriﬁcation を終えたら次は最適 化 に 進 む と よ い で し ょ う。COMSOL は 最 適 化 モ
ジュールの中に図21に示すソルバーを用意しています。

左辺第2項の係数行列の固有ベクトル P を求め、それを使ってその係数行列の対角化を行うと、変数変換 u=Pvによって、
2周波数応答解析におけるモーダル解析の例

図21

最適化モジュールのソルバー群

Veriﬁcation から Validation に進む際に、実験値を利用
して材料特性のパラメタ推定が必要となってきます。
これには定常あるいは時刻列での実験データを CSVファイルとして読み込み、最小二乗法によるパラメタ推定を行うことができます。その際、モデルは数式でも有限要素モデルでも利用できます。この方法は、リチウムイオンバッテリーの電気化学インピーダンス法にも応用できます。
寸法最適化では COMSOL の備えている形状のパラメ
タ表現を利用したり、形状最適化では変形ジオメトリ（ジオメトリの変形に応じてメッシュ再構成）の利用などによって設定操作も簡単に行えます。もちろん、ト
ポロジー最適化も可能です。図22に寸法最適化、形状最適化、トポロジー最適化の例を示します。
寸法最適化という一見単純に見える方法でもうまく利用すれば、電気浸透流というマイクロ流体で重要な問題を簡単に解決できます。図23は流路のベジェ曲線
を最適化することで化学種輸送の精度を大幅改善しました。マイクロ流路では化学種を輸送する際に大きな変形は避けたいものです。iPS 細胞などはわずかな刺激にも反応して予期せぬ分化を生じることは良く知られ（19-7）

Vol.23, No.3 2018

チュートリアル

図22

図23

COMSOL Multiphysics による計算科学工学 – 拡散、固有値、最適化（6）

L を指定して FEM 解析を行なう関数は frequency(L1)といった箇所で利用されています。
これらを含めて、アプリにしたものを図24に示します。ソフトウェアの使い方を知らなくても高度なマルチフィジックス解析や最適化計算をすぐに開始できます。

非一様分布荷重をもつ場合の最適化の例

流路の寸法最適化による化学種輸送の改善例ています。
最適化を行うには目的関数と設計変数を適切に設定する必要があります。自身の設定した設計変数範囲で目的関数が最大あるいは最小となるかが重要です。これは COMSOL ではパラメトリックスイープをすることで概観できます。必ずリハーサルをすることです。
紙数の都合でこのあたりの手順まで十分に記載できませんので資料をぜひ参照してください [5]。
最適化においては、いろいろな手法を試してみたくなります。その際、COMSOL のもつアプリケーション
ビルダーでメソッドエディタによって利用者独自のプ
ログラム (Java) を記述できます。
音叉は腕の長さを変更すると固有周波数が変わります。ここでは、周波数を指定して音叉の腕の長さを決めることを考えてみます [6]。固有周波数 f は腕の長さ Lの関数ですので f(L) と記述できます。指定周波数を
ftarget とすると、F(L)=f(L)-ftarget=0となるような L を決めれば良いということになります。F(L)=0とする L を求めるルーチンはニュートン法がありますが、これはF の微係数も計算する必要が有るので、セカント法を考えます。すると、これは関数値 F(L) のみの計算で済みます。メソッドエディタにおけるポイントとなる箇所のみ考え方を記述します。
関数 f は L を設定して音叉の固有周波数計算を行う部
分を記述します。CAE のシーケンスに沿って作成した
高度な FEM の内容が run というたったの一行で参照されてしまうのは驚きでもあり、寂しさでもあります。

図24

アプリは国内の大学の講義でも利用されるようになってきました [7]。
著者も大学の講義の支援に呼ばれる機会も増えてきました。その際、ソフトウェアを立ち上げて利用手順から説明していくと時間も相当かかり、肝心の利用というところまで十分に到達できません。図25のようにそういった場面に必要な手順を整理し、だれでも5分程度でアプリを作成できる手順とひな形を開発してきています [8]。
図26は著者の所属する会社で営業部員向けに音響の教育を行っている様子です。営業部員からのフィード
バックでアプリの内容もいくつか改善されるという思

オブジェクト指向言語は予期できない境地へ急速に発展しています。設計思想の良さを反映しています。
日本でのソフトウェアづくりにおいても重視すべきところであり、ぜひ参考にすべきところです。
セカント法の適用は概略、次の通りです。
計算工学

アプリの例

（19-8）

図25

簡単かつ応用範囲の広いアプリ作成手順へ
Vol.23, No.3 2018

チュートリアルCOMSOL Multiphysics による計算科学工学 – 拡散、固有値、最適化（6）わぬフィードバックもありました [8]。新人の技術部員
が COMSOL サーバーというソフトをサーバー機にインストールしてくれ、そこでアプリを動かしています。
大勢でアクセスすると運用方式によっては計算負荷も当然増えてきますが、今後、このような方式が研究開
発や業務を根本から革新していくことと思います。

図26

KESCO 営業部員向けの社内教育の様子

専門でない人たちから時間、場所に関わらず、新しいアイデアが提案されてくる可能性があります。100倍の薬品を投入してノーベル賞へつながった例は記憶に新しいと思います。専門の枠を超えて思い切ったチャ
レンジが必要とされていますが、アプリはその強力な道具になり得ます。経営側は自身の考え方を早急に革新する必要があります。社員に適切な道具と環境を積極的に用意し、誰もが解析やフィードバックを行なえる環境をトップダウンで構築すべきです。

4

あとがき

今回で本チュートリアルの連載が終了となります。
関係者の方々には自由な執筆をさせていただき感謝申し上げます。
人生はあっという間に過ぎてしまいます。若い人たちの参考になればと思い、思い出話をします。筆者はロボットをやりたかったのですが大学の先生方が流体関係が多く、流体を専攻することになりました。助手時代は航空分野でマッハ数10の弱電離極超音速流れの熱伝達問題の実験に取り組む過程で、可視化をする重要性に気付き、ランダムチョイス法で燃焼反応を伴う衝撃波の計算を独学開始しました。小口伯郎先生の極超音速強干渉理論やスチュワートソンのトリプルデッキ理論という美しいものに出会ったことも喜びでした。航空分野では人為的なモデル化をきらう傾向があり、これはとても勉強になりました。M 社に入ってから空力開発を開始する中で、私の特命事項であった数値流体の準備をはじめました。PC、エンジンベンチのミニコンや、IBM 大型計算機のメンテナンス時のフリー計算など色々な方々の協力を得てプログラム開発を行い、同僚と離散渦法の結果を国際会議で発表しました。そうしているうちに
DNS を提唱する宇宙科学研究所桑原邦郎先生と出会い、スーパーコンピューティングを導入しました。それまで1か月くらいかかっていた計算が30分もかからず終了していくことを目の前にして、とても驚きました。その後
計算工学

ごく短期間で自動車の空力計算が形状のデジタルデー
タさえあれば可能なところまで到達します。風洞実験との誤差は初めから10% 程度に収まるという素晴らしい成果が出せました。その頃、欧米では誤差はとてつもなく
大きいものでした。一方で、衝突開発が急務となっていた時期でもあり、そちらも支援することになりました。
そのために、夕方まで衝突、夜は空力といった事もやりました。衝突では先輩の多賀宏二さんがその分野の計算のリーダーでした。この方は私が M 社で計算を行な
うための PC を初めて導入してくれた人です。いろんな面で大変お世話になり感謝しています。私は流体の可視化の考えを活かして衝突実験にもいくつかアイデアを出しました。1000万円のカメラを車に溶接して衝突させたときは首を覚悟しましたが、そうはなりませんでした。そのおかげで肝心な部分の可視化ができ、性能が大きく改善されました。先行開発を加速するために横浜に模型風洞を建設し、そこで先行空力開発に携わりました。良い同僚に恵まれ、この風洞では全自動運転・データ自動計測を実現しました。おかげで、先行開発の質は格段に良くなりました。企画設計部やデザイン部のモデラーが風洞を動かして自分で実験できるようにもなりました。さらに、M 社での最後の論文では画期的な発見をしました。その後、桑原邦郎先生と一緒に多方
向風上差分による円柱まわりの計算やソフトの開発などを行いました。桑原先生の研究所には大勢の世界的な研究者や教科書の著者などが多く来所されており、大変勉強になりました。今井功先生が90歳になられても一番シャープな質問をされていたのは印象的でした。神部
勉先生の流体関係のテキストに桑原先生と私でやった円柱の計算結果が引用されたことは大変な喜びでした。
尾川茂さんと共同で神部先生の漸近展開理論に基づいて縦渦の空力騒音の数値解析も実現しました。その後、
濱野生命科学研究財団の事務局長を経験します。ここでは研究者を事務局の立場から観察することができ、勉強になりました。fMRI の BOLD 法の発見者でノーベル賞候補の小川誠二先生のもとで刺激を受けました。
その後、現在の職場で COMSOL Multiphysics を扱うように な り ま し た。 そ れ で BOLD 法 の 数 学 モ デ ル を
COMSOL Multiphysics で計算して、大勢の医者の前で発表することも行いました。
現在の会社に入社したのが50歳を過ぎてからでしたが Maxwell 方程式やらなにやらと幅広い物理の勉強を始めることになり、今も継続中です。人生、勉強に終わりということはなさそうです。
近頃、偶然にも明治大学の萩原一郎先生とお会いすることになりました。萩原先生は自動車開発でとても有名な方で東工大の先生にもなられた先生です。その先生が自動車会社に勤務されていたころの仕事の完成をいま再び目指していることを教えてくださいました。話はいろいろな計画に及びました。とても迫力のあるお話ぷりに私も大きな衝撃を受けました [9]。
いろいろ思うことは、日本が世界を真にリードするとはどういうことだろうかという点です。

（19-9）

Vol.23, No.3 2018

チュートリアルCOMSOL Multiphysics による計算科学工学 – 拡散、固有値、最適化（6）

答えはありませんが、先生方も実践されているように、各人が自分で抱えた懸案事項を解決することに情熱を燃やし続けられるかどうかが重要ではないかと思います。時間が経過すると、学生のころに考えていた
柔軟なロボットの実現例がどんどん出てきていることなど、昔はできにくかったことが出来るようになってくるといったことが起こります。
従って、継続して情熱を持ち続けてとにかく勉強を続けていれば、必ず大きなチャンスが到来するはずです。今は難しくてもそのチャンス到来時にうまく解決
できれば幸せではないかと思います。何かのヒントになればと思い、あとがきに替えます。

■参考文献
[1]
[2]
[3]
[4]
[5]
[6]
[7]

[8]
[9]

計算工学

（19-10）

花 村 克 悟 , マ ッ ク ス・ プ ラ ン ク（1858〜1947）の 功 績 ,J.HTSJ, Vol.48, No.205 (2009).
鈴木克彦 , シュレディンガー方程式 , 共立出版 (2013).
https://www.comsol.jp/model/superlattice-band-gaptool-45281
前野昌弘 , 量子場の理論入門 , http://www.phys.u-ryukyu.
ac.jp/~maeno/ﬁeld.pdf.
橋口真宜 , 米大海 , COMSOL を用いた最適化計算とパラメータ推定 , http://www.kesco.co.jp/ (2018).
https://www.comsol.jp/model/tuning-fork-8499
高野直樹 , COMSOL Multiphysics&Application Builder を用いた計算力学の研究と教育 , COMSOL カンファレンス東
京 (2017).
橋口真宜 , 米大海 , フィジックス教育用のアプリ開発 , 第23回計算工学講演会 (2018).
萩原一郎 , 音 - 構造振動の連成解析の基礎から応用セミ
ナー, http://www.kesco.co.jp/ (2018).

Vol.23, No.3 2018