ryujimiyaの日記

C#を使って数値解析したい

周波数を与えて伝搬定数を求める周期構造導波路のFEM表示式

周期構造導波路の固有値問題を解くFEMには2種類の定式化があります。

 

(1) Φ(x, y) = φ(x, y) exp( - jβx) と置く方法

   FEM表示式としては、伝搬定数を与えて固有周波数を求める式

      [S]  {φ} -  k0^2  [T]  {φ} = {0}

または、動作周波数を与えて伝搬定数を求める式

      [K]  {φ} - jβ  [C]  {φ} - β^2  [M]  {φ} = {0}

が得られます。

 

(2) Φ(x, y)を直接解く方法

   ブロッホの周期境界条件 Φ|2 = exp( - jβd)Φ|1 、dΦ/dx|2 = exp( - jβd) dΦ/dx|1 を適用すると伝搬定数を与えて固有周波数を求める式

      [S]  {Φ} - k0^2  [T]  {Φ} = {0}

が得られます。

※下記記事に引用元を記しています。

 三角形格子のフォトニック結晶欠陥導波路の固有モード計算

  http://ryujimiya.hatenablog.com/entry/2013/01/16/074832

 

さて、(2)の定式化で計算する場合、[S]や[T]が exp( - jβd) の関数になるため、伝搬定数を与えて固有周波数を求める方法をとりました。

これは動作周波数を与えて伝搬定数を求めることができないと思っていたからですが、実はできるようです。

FEM表示式を観察するとexp( - jβd)の項は、exp( - jβd)と、exp(jβd)のみなので、

両辺をexp( - jβd)倍すると、exp( - jβd)に関する2次の固有方程式になっています。

つまり、λ = exp( - jβd) と置くと、

     [K]  {Φ} + λ  [C]  {Φ} + λ^2  [M]  {Φ} = {0}

の形に変形できます。ここで

    [K] = [ [P21]  [P20] ]

            [   [0]    [0] ]

    [M] = [ [P12]    [0] ]

             [ [P02]    [0] ]

    [C] = [ [P11] + [P22]  [P10] ]

            [       [P01]        [P00] ]

    {Φ} = [ {Φ}1   {Φ}0 ] 

    [P] = [S] - k0^2  [T]

です。[K]と[M]が特異行列になっているので[M]-1を用いて標準固有値問題に変換することはできませんが線形化は可能で一般化固有値問題として解くことができます。

 ※2次固有値問題の線形化については、下記の記事を参照してください。

  2次元導波路伝搬定数の周期構造解析用FEMによる計算(2) ー非線形固有値方程式の解法についてー

  http://ryujimiya.hatenablog.com/entry/2012/11/28/015112

 
【計算結果】
エアホール型の三角形格子で計算してみた結果を示します。エアホールの半径r = 0.30a、誘電体基盤の屈折率n = 2.76です。横軸は a/λ、縦軸は波数ベクトル βd / 2π です。
偶モードの分散特性

f:id:ryujimiya:20130130012720p:plain

計算できない周波数はこの方法でもありますが、それ以外は前回計算したβを与えて動作周波数を求める方法の結果と一致しています。