学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度協同研究 最終報告書 2012 年 5 月

11-NA24３次元有限要素法による光導波路解析の高速化と最適設計に関する検討
辻

寧英 (室蘭工業大学)

概要 光通信の高速・大容量化への要求を実現するため，そこで使われる光導波路のさらなる小型化・高性能化への要求がますます高まっている．こうした光導波路デバイスの性能を追求していくとデバイス構造が複雑化する傾向にあり，デバイス特性を正確に把握するためには解析的な近似を用いずに直接３次元構造を解析する必要がある．また，そこからさらに進めて，所望の特性を与えることでそれを実現する光導波路デバイスを発現させることができる最適設計への期待も高まっているが，有限要素法は任意形状への適用性に優れた汎用的な解析法として知られている．一方，解析領域を要素に分割する必要があり，さらに最終的に得られる方程式が大規模な連立一次方程式となるため大きな計算資源を必要とするなどの課題がある．
本研究では，有限要素法による光導波路の３次元解析の高速化・高性能化と，有限要素法とトポロジー最適化を利用した光導波路の最適設計についての基礎的な検討を行っている．

1

研究の目的と意義をつけるアダプティブメッシュの開発を行うとともに，連立一次方程式の解法としてメモリの削減が可能な反復法

インターネットをはじめとする情報通信技術の発展にの解の収束性と領域分割法による計算の効率化について伴い，そこで扱われる情報は多種多様なものとなり，光通検討を行う．さらに，ここで開発する 3 次元有限要素法

信システムでは動画像まで含めた大容量のデータを高速とトポロジー最適化を組み合わせたデバイス最適設計法に伝送でき多様なデータに対応できるような柔軟なネッを開発することにより，高性能な光導波路デバイスの設

トワークを構成できることが期待されている．光通信シ計を試みる．ここで開発する自動最適設計法では，必ずステムの性能を改善するためには，光ファイバ自身の特しも初期条件を与えなくてもデバイスの最適設計が可能

性はもちろんのこと，そこで使われる光導波路デバイス

であり，これまで知られていない新たなデバイス構造をの性能向上が必須である．これまで，光導波路デバイス

見出せる可能性がある．

の設計は解析的近似解法や 3 次元構造を 2 次元構造に近似して行われることが多かった. しかし，光導波路デバ

2

イスはより高性能化を目指して小型化・複雑化する傾向にあり，計算機を利用した数値解析法により直接 3 次元

当拠点公募型共同研究として実施した意義
超大規模数値計算系応用分野において，北海道大学情解析することが求められるようになってきている．この報基盤センターと共同研究を行った．

ようなことから，さまざまな数値解析法が開発されてい光導波路デバイスの 3 次元有限要素法解析では，最終

るが，なかでも有限要素法は時間領域，周波数領域両方的に大規模な連立一次方程式を解く必要があり，そのたの解析が可能で，任意形状への適用性に優れた汎用性のめには大容量のメモリを有する計算機が必要となる．大高い解析法である [1]-[3]．また，2 次元有限要素法とト型計算機を利用することにより，通常のパーソナルコン

ポロジー最適化を組み合わせたデバイスの自動最適設計

ピュータでは解析が困難であった大規模な解析が可能に法の開発も行われている [4]-[5]．しかしながら，有限要なり，より実際的なモデルでの光導波路の解析・設計が素法は最終的に大規模な連立一次方程式を解く必要があ

行えるようになる．また，北海道大学情報基盤センター

り，特に 3 次元問題に対しては計算規模がすぐに大きくがこれまでに蓄積してきている大規模計算の効率化に関なってしまうという課題がある．

する多くの知識やセンターで開発している行列計算ライ

本研究では，高性能な光導波路デバイスの最適設計を

ブラリーを利用することで，計算の効率化を図れること目的として，光導波路解析のための汎用性の高い 3 次元が期待される．

有限要素法の高速化，光導波路デバイスの最適設計法に北海道大学情報基盤センターにおいては，計算の高速

ついての検討を行う．有限要素法では要素分割の粗密を化のためのプログラムチューニングおよび大規模データ変えることで計算の効率化を図ることが可能であり，有解析のための反復法ライブラリの整備とカスタマイズを限要素法により生成される行列が疎行列であることを利担当してもらった．

用することでメモリの削減が可能である．ここでは，光導波路中を伝搬する光波の分布に応じて要素分割に粗密1

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度協同研究 最終報告書 2012 年 5 月

いま図 1(b) に示すような入射面 Γin を考え，領域 Ω0ΩBを部分領域 ΩA ，ΩB に分け，各領域に対して有限要素ΩA法を適用すると最終的に以下の行列方程式を得る．

ΦB,in ΦA,in[Pi ] {Φi } = {ui }

ΩA

(3)

Φscat

ΩB

ここに 下添字 i = A, B は領域 i に関する量であること∑を表し， ei は小領域 i に関する全要素の和を意味する．

Γin(a) 全体図(b) 上面図

また，行列 [Pi ]，ベクトル {ui } は以下のように与えられる．

図 1: 3 次元光導波路[Pi ] = [Ki ] − k02 [Mi ]研究成果の詳細と当初計画の達成状況

3[Ki ] =
本研究では，有限要素法に基づく光導波路デバイスの∑ ∫∫∫路デバイスの最適設計についての検討を行った．

[Mi ] =

完全整合層を用いた有限要素法{ui } = −jk0

材料まで取り扱えるように，導波路を構成する材料の比µxx µxy µxz

{N }Γ · (in × Ψ)|Γi dΓi(7)の方程式にするために，入射面 Γin 上の Φ の値からなµzx µzy µzz

るベクトルを {Φi }Γ ，{Φi } から {Φi }Γ を取り除いたベ
クトルを {Φi }0 とすると，式 (3) は
[
][
] [
]
[Pi ]00 [Pi ]0Γ{Φi }0{ui }0=[Pi ]Γ0 [Pi ]ΓΓ{Φi }Γ
{ui }Γ

と表される場合を考える．このとき，マクスウェルの方
程式から，ベクトル波動方程式は∇ × ([p]∇ × Φ) − k02 [q]Φ = 0Γi

次に，領域 ΩA と ΩB に関する方程式をまとめて一つ[ε] =  εyx εyy εyz  , [µ] =  µyx µyy µyz 
εzx εzy εzz∑ ∫∫Γi に接する要素についての和を表し，境界 Γi は入射面
Γin と外部境界 Γ∞ からなる．

誘電率テンソル [ε]，比透磁率テンソル [µ] がεxx εxy εxz(6)であり， {N }Γi は境界 Γi における形状関数ベクトル，∑∑ei は領域 Ωi 内の全の要素についての和， Γi は境界

入射する場合を考える．ここでは，一般的な任意異方性

{N } · ([q]{N }T )dxdydzeiΓi

図 1 に示すような光導波路を考え，入射面 Γin から光を

(5)

∑ ∫∫∫
ei

3.1

(∇ × {N }) · ([p]∇ × {N }T )dxdydz

eiei

3 次元解析の高速・高精度化と，これを利用した光導波

(4)

(1)

(8)

いま，ΩA ，ΩB それぞれから見た入射面 Γin への入射電で与えられる．ここに k0 は自由空間波数であり，Φ，Ψ，

磁界を Φi,in ，散乱電界を Φi,scat (i = A, B) とすると，p，q は未知変数を電界とするか磁界とするかにより

領域 Ωi から見た Γin における電磁界 Φi は√√
Φ = ε0 E, Ψ = µ0 H, [p] = [µ]−1 [q] = [ε]√√Φ = µ0 H, Ψ = − ε0 E, [p] = [ε]−1 [q] = [µ]Φi = Φi,in + Φi,scat(9)と表され，入射面 Γin における電磁界の連続性を考慮すと表される．

ると，領域 ΩA ，ΩB に対する行列方程式 (8) を結合でき解析領域を四面体エッジ要素を用いて離散化すると，る．いま，入射は領域 ΩA 方向のみとすると，最終的に電磁界ベクトル Φ は各要素内で解くべき行列方程式は
[PA ]00 [0][PA ]0Γ[PB ]0Γ [0] [PB ]00

Φ = ({U } ix + {V } iy + {W } iz ){Φ}eT= {N }T {Φ}e

T

T

(2)




{ΦA }0  {ΦB }0 

[PA ]Γ0 [PB ]Γ0 [PA ]ΓΓ + [PB ]ΓΓ{Φscat }Γ
−[PA ]0Γ {ΦA,in }Γ


=
{0}

{uA,in }Γ − [PA ]ΓΓ {ΦA,in }Γ

と表される．ここに {Φ}e は各要素内辺上での電磁界の辺方向の大きさ，{U }，{V }，{W } はそれぞれ x，y ，z方向成分に対する形状関数であり，T は転置することを意味する．また，ix ，iy ，iz はそれぞれ x, y, z 方向の単位ベクトルである．

となる．

2

(10)

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度協同研究 最終報告書 2012 年 5 月

port 2

光導波路デバイスの伝送特性を有限要素法により解く

port 3ものとすると，最終的に以下のような行列方程式を得る．

[P ] {Φ} = {u}port 1
Design region(11)

いま port 1 から光を入射したとして式 (11) を解いて伝搬界分布 {Φ} が求まると，第 n 番目の導波路の導波路への固有モードの振幅透過係数 Sn1 はport N

Sn1 = {Φ}T {gn }図 2: 最適化の問題設定(12)として求めることができる．ベクトル {g} は

start{gn } =

initial refractive
index distribution)cβi ([Mtt ]{φt,n } + [Mtz ]{φ0z,n }
k0

(13)

で与えられる．
finite element
analysis
stopping
criterion
Yes

update refractive
index distributionNo

まず，導波路が等方性媒質からなる場合を考え，比誘電率分布が M 個のパラメータ ai (i = 1, 2, · · · M ) により

sensitivity
analysisεr = εr (a1 , a2 , · · · , aM )

(14)

stop

と表されているとする．このとき，i 番目のパラメータ2

ai が微少変化したときの |Sn1 | の変化は以下の式で表される．

図 3: 最適化の流れ

3.2

トポロジー最適化を用いた 3 次元光導波路デバイ∂ |Sn1 |∂S ∗∂Sn1 ∗= Sn1 n1 +
S
∂ai
∂ai
∂ai n12

スの最適設計(15)

トポロジー最適化では，最適化領域内の屈折率分布をい

ここで，式 (15) 中の ∂Sn1 /∂ai は随伴変数法 (AVM) を

くつかの数値パラメータにより表現し，そのパラメータ用いることで効率的に計算することができる．∂Sn1 /∂aiの変化による特性の変化 (感度) を求め，その感度に基づをパラメータ ai の変化に伴う電磁界分布の変化を用いて{}T∑ ∂Sn1 ∂Φj
∂Sn1
∂Sn1
∂{Φ}
=
=
(16)
∂ai
∂Φj ∂ai∂Φ∂ai

いてパラメータを更新することで特性が改善する方向に
最適化領域内の屈折率分布を更新する．これを繰り返しjていくことで最終的に所望の特性を実現する光導波路デのように表すことができる．一方，入出力および仮想境

バイスを見出すことができる．

界は最適化領域に含まれないとして，有限要素法の最終

図 2 に示すような光導波路デバイスを考え，port 1 から光を入射してそれぞれのポートに所望の出力が得られ

式 (11) を εr で微分し，これを式 (16) に代入するとる設計領域内の屈折率分布を求める最適化問題を考える．

∂Sn1
T ∂[P ]
= − {λn }
{Φ}
∂ai
∂ai

図 3 に最適化の流れをフローチャートとして示す．まず，設計領域内に適当な初期屈折率分布を与え，有限要素法(17)の表現を得る．ここで，{λn } を表す式はによる光導波路解析を行い入出力特性を求める．次に，{感度解析により屈折率分布を微小変化させた場合の特性T{λn } =の変化を調べ，特性が改善する方向に屈折率分布を更新する．この手順を繰り返すことで所望の特性を実現する

∂Sn1
∂Φ

}T

[P ]−1(18)であり，式 (12) を考慮すると式 (18) は

光導波路デバイス構造を得ることができる．屈折率分布の変化に対する感度を知るためには，実際に構造を変化[P ]T {λn } = {gn }させてみて特性を調べるのもひとつの方法であるが，そと表すことができる．したがって {λn } は有限要素法解(19)のためには新たな光導波路解析が必要であり，最適化す

析を 1 回行う程度の計算量で求めることができ，[P ] がるパラメータが増えるとその計算量は膨大になる．ここ対称行列で有限要素法解析のときに分解計算がなされてでは，構造の変化に対する特性の変化を効率的に調べるいる場合には単なる代入計算により高速に求めることが方法として，随伴変数法 (AVM) を用いている．

できる．

3

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度協同研究 最終報告書 2012 年 5 月

また ∂[P ]/∂ai は具体的には未知変数が電界の場合には∑ ∫ ∫ ∂εr∂[P ]= −k02{N } · {N }T dxdydz (20)
∂ai
∂a
i
e
e
未知変数が磁界の場合には∑ ∫ ∫ ∂εr 1∂[P ]=−(∇ × {N })2∂aie ∂ai εre

· (∇ × {N }) dxdydz(21)で与えられる．
光導波路デバイスの構造を表す屈折率分布は様々な方法により定義することが可能であるが，一般的に広く用いられている密度法では，設計で想定する材料の中間的

図 4: 3 分岐光導波路の最適化モデル

な屈折率を持った材料領域 (グレイ領域) が発生する問題



0





 ξ+h
h2
×
h−ξ


1−


h2


0

がある．ここでは文献 [5] で提案している関数展開法を用いている．この方法ではグレイ領域の問題を本質的に取り除くことができる．
標準的な 2 媒質を対象とした関数展開法では最適化領域内の比誘電率分布を適当な解析関数 w(y, z) を用いて以下のように表現する．

(ξ ≤ −h)(−h < ξ < 0)(0 ≤ ξ < h)(25)(ξ ≥ h)と求めることができる．この方法の特徴の一つとして，

εr (x, y, z) = εra + (εrb − εra ) H(w(x, y, z))

なんらかの対称性を有する構造を見出したい場合に，(22)fi (x, y, z) として目的の対称性を持った関数を選ぶこと
ここで，εra ，εrb は使用可能な 2 つの材料の比誘電率，で容易に対称性を導入できることがある．fi (y, z) の選

H(ξ) は ξ の値によって 0 か 1 かの値を取る関数であり，び方は特に決まったものはなく様々なものが考えられる．

εr は w(x, y, z) の値によって εra あるいは εrb のどちらかの比誘電率となる．ただし，実際には εr が微分可能となるように，H(ξ) は以下のように定義される連続関

ここでは関数 w(x, y, z) を以下のようにフーリエ級数の形で表現することにする．

数とする．


0


 (
)2



 1 ξ+h2h()2
H(ξ) =

1 ξ−h


1−


2
h



1

w(x, z) =(ξ ≤ −h)θij =(−h < ξ < 0)
(23)

2πi
2πj
x+
z
Lx
Lz

(27)

ここに Nx ，Nz はそれぞれ x 方向，z 方向の展開項数を(0 ≤ ξ < h)表し，Lx ，Lz はフーリエ級数の周期を表す．フーリエ(ξ ≥ h)級数は本来周期関数を表現するためのものであるが，最適化構造が周期構造となる制約条件が加わることは必ずしも好ましくない．その場合には，最適化領域の x，zであり，最終的に h → 0 とすることで，明確な屈折率境

方向の幅を Wx ，Wz として，Lx > Wx ，Lz > Wz とす

界を表現できる．

る必要がある．感度解析に基づき ai ，bi を更新していく

最適化領域内の屈折率分布を決める関数 w(x, y, z) はことにより最適な光導波路デバイス構造を見出すことが一般的にw(x, y, z) =(aij cos θij + bij sin θij ) (26)i=−Nx i=−Nz

ここに h は H(ξ) が連続関数となるように導入された量∑N∑x −1 N
z −1
∑

できる．

ai fi (x, y, z)

(24)

i

3.3

の形で与えられる．このとき，式 (20)，(21) 中の ∂εr /∂ai

光デバイスの設計例

ここでは光デバイスの設計例として，3 分岐導波路，3 波は

長分離素子，交差導波路についての設計例を示す．

∂εr= fi (x, y, z) (εrb − εra )∂ai4学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度協同研究 最終報告書 2012 年 5 月

初期構造

最適化構造

伝搬波形

(a) 初期構造を与えない場合

初期構造

最適化構造

(a) 最適化構造(b) 伝搬波形 (λ = 1.45 µm)(c) 伝搬波形 (λ = 1.55 µm)(d) 伝搬波形 (λ = 1.65 µm)

伝搬波形(b) 単純な 3 分岐を初期構造とした場合

初期構造

最適化構造

normalized power

0.5

伝搬波形

(c) ギャップ持つ 3 分岐を初期構造とした場合

0.4
0.3

0.1

port 1
port 2,4
port 3

0.0
1.45

1.50

0.2

λ i =1.55 µm
λ i =1.45, 1.50, 1.55
1.60, 1.65 µm

1.55
1.60
wavelength [ µ m]

1.65(e) 出力パワーの波長依存性
初期構造

最適化構造

伝搬波形

(d) 目的波形から初期構造を決めた場合

図 6: 波長依存性を低減した 3 分岐素子

図 5: 初期構造の違いによる最適化構造の違い

出力があらかじめある程度等しくなるように直進ポートとの間にギャップを設けておいた場合である．また，(d)3.3.1

はあらかじめ曲り導波路の最適化を行い得られた界分布

3 分岐光導波路と直線導波路の界分布から 3 分岐導波路の目的とする界

図 4 に示すような光導波路デバイスを考え，port 1 か分布を決め，そこから初期構造を決めた場合である．初ら入射した光を port 2∼4 に等分配させるような設計領期構造により最適化構造は異なるが，いずれの場合にも域内の屈折率分布を求める．コア，クラッドの屈折率は各ポートに出力がほぼ等分配される結果が得られている．

n1 = 3.2，n2 = 1 とし，設計領域内は屈折率 n1 ，n2 の 2
つの材料のみで構成されるものとする．入出力導波路の幅を w = 0.2 µm，最適化領域のサイズを Wx = Wy =これはここで取り上げたデバイスが大きな比屈折率差を持つため，目的の特性を実現する構造が多数存在するためと考えられる．したがって，目的に合った初期構造を

2 µm とし，波長 1.55 µm の TE 基本モードが入射する場合を考える．目的関数を，port 2∼4 へ出力を等分配させるように，以下のように与える．
)2
4 (
∑
1
− |Si1 |2
Minimize C =3

決めることが重要である．

表 1: 各ポートの規格化透過パワー
port 1
port 2port 3
port 4(28)i=2

図 5 に 4 種類の初期構造による最適化結果の違いを示し，
表 1 にそれぞれの最適化構造に対する各ポートの規格化
出力パワーを示す．(a) は初期構造を与えず，(b) は初期構造を単純な 3 分岐とした場合，(c) は 3 つの導波路への

5

(a)

0.000

0.332

0.331

0.332

(b)
(c)
(d)

0.000
0.000
0.000

0.333
0.332
0.327

0.331
0.332
0.326

0.333
0.332
0.327

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度協同研究 最終報告書 2012 年 5 月

air

t1

t2

Si
SiO 2

(a) 最適化構造(b) 伝搬波形 (λ = 1.31 µm)
n3 (air)

n2

n1

n2

n3 (air)

図 8: コア部の高さと等価屈折率(c) 伝搬波形 (λ = 1.49 µm)

3.3.3(d) 伝搬波形 (λ = 1.55 µm)

交差導波路

導波路交差部でのクロストークを低減するための交差

1.0
normalized power

導波路の設計例を示す．ここでは図 8 に示すような高さの異なるコア領域を持つ場合に対して検討する．
port 1
port 2port 3
port 4図 4 の最適化モデルにおいて，波長 λ = 1.50 µm∼
1.60 µm，導波路幅 w = 0.5 µm，最適化領域の幅 Wx =0.0
1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60
wavelength [ µ m]

および t2 = 97.3 nm とすると，等価屈折率は n1 = 2.73,

0.5

Wy = 6 µm とする．屈折率は，コアの屈折率を 3.475，ク
ラッドの屈折率を 1.444 とし，コアの高さを t1 = 194 nm
n2 = 2.1, n3 = 1.0 である．目的関数は(e) 出力パワーの波長依存性

Minimize C =K∑(1 − |S31 (λi )|2 )

(29)

i=1

図 7: 3 波長分離素子の設計例

であり，K = 5，λ1 = 1.50，λ2 = 1.525，λ3 = 1.55，

λ4 = 1.575，λ5 = 1.60 としている．

図 5 に示した最適化構造のうち (b) 以外の構造は比較的波長依存性の大きな構造であった．次に，(b) 以外の

図 9 に最適化に用いた初期構造と，構造表現に用いた構造に対しても波長 λ = 1.45 ∼ 1.65 µm の範囲で波長フーリエ級数の展開項数に対する最適化構造の関係を示

依存性の小さな 3 分岐素子とするため，例として図 2(d)す．展開項数が少ないほど構造は単純であり，展開項数の最適化構造を初期構造として，λ = 1.45, 1.50，1.55,が増すにつれて微細な構造が現れていることがあわ k る．

1.60，1.65 µm の 5 波長において等分配される再度最適

図 10 に，展開項数に対する目的関数の収束の様子を示化を行った．図 6 にその結果を示す．図 6 (e) より，初す．図より N x = N y = 6 以上の展開項数で良好な特性

期構造 (破線) に対して最適化構造 (実線) では波長依存が得られているが，初期構造との関係で Nx = Ny = 7, 8

性が抑圧されていることがわかる．

において特性が劣化していることがわかる．図 11 には最適化の反復に対する各波長の透過パワーの収束の様子

3.3.2

を示す．反復に伴いほぼ一様に特性が改善されていく様

3 波長分離素子子がわかる．図 12 に最適化構造に対する透過パワーの波長選択性を有する光導波路デバイスの設計例として，

波長依存性を示す．図には最適化前の構造での波長特性

図 4 の最適化モデルで波長 λ = 1.31, 1.49, 1.55 µm の 3と導波路高さを一様とした場合の波長特性をあわせて示

つの波長の光をそれぞれ port 2, 3, 4 に分離するデバイしている．ここで得られた構造では非常に広い波長範囲スの最適化を行った．初期構造には図 5 (d) の最適化構

でクロストークを抑圧できていることがわかる．

造を用いた．最適化結果を図 7 に示す．図より，それぞれの波長が目的のポートに分離されていることがわかる．

6

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度協同研究 最終報告書 2012 年 5 月(a) 初期構造

Normalized transmitted power

1.00(b) Nx = Ny = 4

0.95
λ = 1.50 µ m

0.90

λ = 1.525 µ m

0.85

λ = 1.55 µ m

0.80

λ = 1.575 µ m

0.75

λ = 1.60 µ m

0.70
0

200

400

600

800

1000

Number of iteration

図 11: 各波長に対する透過パワーの収束の様子
1.00(d) Nx = Ny = 16

Normalized transmitted power(c) Nx = Ny = 8

図 9: 初期構造と異なる展開項数に対する最適化構造
1

10

0

C

10

10−1
10−2
10−3
0

Three

optimized
optimized (without tiny holes)
non−optimized

0.90 materials

0.80

Two
materials

0.70

0.60
1.45

1.50

1.55
1.60
Wavelength, λ (µm)

1.65

図 12: 交差導波路の透過パワーの波長依存性
4812Nx ( = Ny)16

図 10: 展開項数に対する目的関数の収束の様子(a) λ1 = 1.50 µm

3.4

当初計画の達成状況(b) λ3 = 1.55 µm(c) λ5 = 1.60 µm

図 13: 交差導波路中の伝搬波形

3 次元光導波路解析のための有限要素法の高性能化の検討を行い，有限要素法とトポロジー最適化を用いるこ討を行ってきた関数展開法に基づくトポロジー最適化にとで 3 次元光導波路デバイスの最適設計を行えることをついて初期構造や安定な収束のための反復手順ついてさ示した．また，最適化における反復計算を工夫することらに詳細な検討を行った．これにより，効率的な最適化

で，より安定により単純な構造に解を収束させられるこが行えるようになり，目的の特性を有する光導波路をよとを示した．

り単純な構造で実現できるようになった．

当初計画のうち，計算の効率化の観点から，アダプティ

今回の検討において，当初計画していた 3 次元光導波

ブな有限要素メッシュの生成による計算の効率化と連立路のトポロジー最適設計については，波長依存性を持っ一次方程式の反復解法にによる必要とされる計算機メモ

た光導波路デバイスも含めて効率的な最適設計を行える

リの低減について検討を進めている．これまでに，スーようになった．一方で，大規模な光導波路解析問題へのパーコンピュータ向けの計算プログラムのチューニング

適用という面では，連立一次方程式の反復解法においてと複素反復法ライブラリの開発を通じて，これまでにな

まだ十分な回の収束性が得られない場合があり，引き続

い強力な前処理法を開発し，これまでより多くのモデル

き検討を行っている．

に対して比較的安定して解が収束するようになっている．
目的の特性を実現する構造を自動的に生成するための
光導波路デバイスの最適設計法に関しては，これまで検7

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度協同研究 最終報告書 2012 年 5 月

今後の展望

4

5.2

(1) 木島涼輔，藤本幸太，平山浩一，辻 寧英，佐藤慎悟，“3次元光導波路デバイスのトポロジー最適設計における初期条件の検討”，電子情報通信学会技術研究報告，Vol.111,
No.149，EST2011-36, pp. 73-78，2011 年 7 月.

今回の検討では，連立一次方程式の解法には計算速度の観点から直接解法を用いた．しかしながら，より大規模な問題に対しては直接解法ではメモリの制限が大きく，(2) 清水省吾，瀬田純己，平山浩一，辻 寧英，佐藤慎悟，“波面整合法を利用した 3 次元光導波路デバイスの最適設計の検討”，電子情報通信学会技術研究報告，Vol.111, No.149，
EST2011-37, pp. 79-84，2011 年 7 月.

反復解法を用いた解析が必要となる．現状では，これま
でよりも多くの問題で反復解法が収束するようになっているが，計算速度の面で直接解法よりも劣るという問題がある．直接解法に比べて大規模解析が可能であるとい(3) 仲 祐輔，佐藤 慎悟，平山 浩一，辻 寧英，山崎 慎太郎，西脇 眞二 “レベルセット法による導波管 H 面 T 分岐回路の導体形状最適設計の検討”，電磁界理論研究会，EMT-11-97，
2011 年 7 月.

う利点はあるが，最適設計において多数回の行列計算が必要となるため，この反復解法のさらなる高速化が課題として残されている．

(4) 瀬田純己，辻 寧英，平山 浩一，佐藤 慎悟，“結合係数が連続的に変化する方向性結合型光導波路デバイスの最適設計の検討”，電磁界理論研究会，EMT-11-100，2011 年 7月.

光デバイスの最適化においては，目標とする特性を実現するデバイス構造を比較的単純な構造で実現することが可能になってきている．しかしながら，初期構造の違い(5) 藤本幸太，平山浩一，佐藤慎悟，辻 寧英，“多媒質からなる光回路および非相反光回路に対するトポロジー最適設計の検討”，電磁界理論研究会，EMT-11-101，2011 年 7
月.

により局所的な最適解に陥ってしまうという問題が残っている．局所的な最適解を避けるための目的関数の工夫
や，屈折率表現に用いる関数系の工夫などが検討すべき(6) 江口真史，辻 寧英，“単一偏波単一モード楕円格子コア
フォトニックバンドギャップファイバの設計”，電子情報通信学会技術研究報告，Vol.111, No.298，OFT2011-51，
pp.25-28，2011 年 11 月．

課題として残されている．
今後はより大規模な問題への適用を目標として，連立一次方程式の解法に反復解法を用いることと，最適解に(7) 辻 寧英，木島涼輔，佐藤慎悟，安井 崇，平山浩一，“
光回路デバイス最適設計のための関数展開法に基づくトポロジー最適化に関する検討”，電子情報通信学会総合大会，CS-6-4，2012 年 3 月．

おいて反復解法を用いた場合の効率的な感度解析の方法について検討を進める予定である．また，ここで得られた成果を基に，光導波路デバイスの新たな可能性を探る
予定である．

5

研究成果リスト

5.1

学術論文

国内会議発表(8) 後藤裕之，辻 寧英，平山浩一，“関数展開法に基づくトポロジー最適化による導波モード分散特性の最適化に関する検討”，C-15-10，2012 年 3 月．

参考文献
[1] Y. Tsuji and M. Koshiba, Finite element method using
port truncation by perfectly matched layer boundary
conditions for optical waveguide discontinuity problems, IEEE/OSA Journal of Lightwave Technology,
Vol. 20, No. 3, pp. 463–468, March 2002.

(1) K. Hirayama, Y. Tsuji, S. Yamasaki, and S. Nishiwaki,“Design optimization of H-plane waveguide component
by level set method,” IEICE Transactions on Electronics, Vol. E94-C, No. 5, pp. 874–881, May 2011.
(2) 曽根宏靖，原田康浩，今井正明，辻 寧英，中村真毅，“重水中テーパーファイバーからのスーパーコンティニウム
光パルスのスペクトルと位相の数値解析：テーパーウェス
ト長による影響”，光学, Vol. 40, No. 8, pp. 439–447,
Aug. 2011.

[2] N. Kono and Y. Tsuji, Oriented perfectly matchedlayer with flexible parameters for waveguide discontinuity problems, IEEE Photonics Technology Letters,Vol. 16, No. 4, pp. 1089–1091, April 2004.
[3] N. Kono and Y. Tsuji, A novel finite-element method
for nonreciprocal magneto-photonic crystal waveguides, IEEE/OSA Journal of Lightwave Technology,
Vol. 22, No. 7, pp. 1741–1747, July 2004.

(3) H. Gotoh, M. Koshiba, and Y. Tsuji, “Finite-element
time-domain beam propagation method with perfectly
matched layer for electron waveguide simulations,” IEICE Electronics Express, Vol. 8, No. 16, pp. 1361–
1366, Aug. 2011.

[4] Y. Tsuji, K. Hirayama, T. Nomura, K. Sato, and S.
Nishiwaki, Design of Optical Circuit Devices Basedon Topology Optimization, IEEE Photonics Technology Letters, Vol. 18, No. 7, pp. 850–852, April 2006.

(4) K. Fujimoto, Y. Tsuji, K. Hirayama, T. Yasui, S.
Sato, and R. Kijima, “A study on topology optimization of optical circuits consisting of multi-materials,”
IEEE/OSA Journal of Lightwave Technology, Vol. 30,
No. 13, pp. 2210–2215, July 2012．

[5] Y. Tsuji, and K. Hirayama, Design of Optical Circuit
Devices Using Topology Optimization Method With
Function-Expansion-Based Refractive Index Distribution, IEEE Photonics Technology Letters, Vol. 20, No.
12, pp. 982–984, June 2008.

(5) 仲 祐輔，平山浩一，安井崇，佐藤慎悟，辻 寧英，山崎慎太郎，西脇眞二，“レベルセット法による導波管 H 面回路の導体形状最適設計”，電子情報通信学会論文誌，採録
決定．

8