光線追跡アルゴリズムの選択が解に与えるの影響

2017年 6月 13日

光線追跡は, 高周波光学シミュレーションに効果的なツールです. COMSOL Multiphysics®ソフトウェアの光線光学モジュールは, 光線追跡にマルチフィジックス対応の波面法を使用します. このブログでは, COMSOL Multiphysics の光線追跡アルゴリズムが, 標準の幾何光学教科書に記載されている従来の光線追跡アルゴリズムと異なる理由を探り, シミュレーション結果を最大限に活用するための一連の設定方法を提案します.

シーケンシャル, ノンシーケンシャル, および厳密光線追跡アルゴリズム

光線追跡アルゴリズムは, シーケンシャル法とノンシーケンシャル法に分類できます. 光線光学モジュールで使用される光線追跡法は, ノンシーケンシャル波面法として分類できます.

光学系を介した光線追跡は, 使用されているアルゴリズムの詳細に関係なく, 主に2つの交互のステップで行われます.

  1. オブジェクト上の点または光学系の面上の初期位置と光線の方向が与えられると, 光線追跡アルゴリズムは, 光線が次の面に当たる場所を決定し, そのポイントまで光線を追跡します.
  2. 次に, 光線が面に当たるときに反射や屈折などの境界条件を適用することによって光線の方向を更新し, 追跡が続行されます. これにより, 反射光線または屈折光線の方向が変わり, 後続の材質を介して追跡できるようになります.

多数の面を持つ複雑な系でも, 通常, 光線追跡プロセス全体を, 以下に示すように, これら2つのステップの連続する反復に分割できます.

光線追跡プロセスの模式図.

いくつかのシーケンシャル光線追跡アルゴリズムは, 幾何光学に関する標準テキストで説明されています.(例えば, Refs. 1, 2, と3を参照). シーケンシャル光線追跡とは, 光線が遭遇するレンズまたはミラーの順序が確定しており, 光線追跡アルゴリズムへの入力として指定されることを意味します. つまり, 光線が最初にレンズ1を通過し, 次にレンズ2を通過し, 等々と続くことは分かっていますが, 光線が各レンズ面のどこを通過するか, あるいは反射および屈折する光線の角度がどうなるかは正確には分かりません.

ノンシーケンシャル光線追跡法では, 光線がレンズやミラーにぶつかる順番をあらかじめ指定しません. その代わり, 光源を指定すると, 光線と系内の光学部品との相互作用が自動的に決定されます.

いくつかの重要なシーケンシャル光線追跡アルゴリズムは, 近軸近似を利用し, 光線と光軸の間の角度が小さいことを前提としています. 対照的に, 特定のシーケンシャルおよびノンシーケンシャル光線追跡法は, 光線とレンズ表面の交点を予測するときにレンズの形状が考慮されるという点で, 厳密な方法 (1, 2, 3)として分類できます. 近似は使用されません. したがって, 光線は平面ではなく, 湾曲したレンズの実際の表面と交差します.

厳密な分析光線追跡法では, 光線とレンズまたはミラーの表面との交点は, 代数方程式を解くことによって計算されます. たとえば, レンズの表面が球形の場合, 光線の初期位置と方向が与えられると, レンズの表面との交点は2次方程式を解くことによって正確に計算できます. より精巧な非球形の表面の場合, 光線と表面の交点の方程式には解析解がない可能性があり, 数値的アプローチが必要になる場合があります.

以下の説明では, COMSOL®ソフトウェアで使用されているアプローチを, 上記の厳密な光線追跡方法と比較します. COMSOL Multiphysics の光線追跡は常にスネルの法則の完全な実装を使用し, 小角度近似を使用しないため, 近軸光学のより単純な方法は無視します.

均一および不均一材質

多くの従来の光線追跡方法では, 各材料が均質である必要があります. つまり, 屈折率は各材質で勾配がゼロで, 境界でのみ不連続に変化します. 上記のように, 次のレンズとの交点を代数的に解くことに基づくアプローチは, 均一な材質に対してのみ有効です. これは, 不均一な材質中では光線が湾曲している可能性があるため, 光線の初期位置と方向だけでは, 面と交差する場所を決定するのに十分でないためです.

実際の光学系には, 設計上またはマルチフィジックス効果のために, 不均一材質または異種の材質を使用した部品が含まれている場合があります. 設計による段階的な屈折率分布の例には, リューネブルクレンズ, マックスウェルの魚眼レンズ, および電気泳動によって製造されたGRINレンズが含まれます. 不均一な温度分布を持つ環境にさらされるレンズは, 温度とひずみに依存する材料特性のために, 不可避的に段階的な屈折率を持ちます. これらは, 波面法が優れている最も重要なケースのいくつかです.

異種材質の光線経路の例を次のプロットに示します. 画像では, サーマルカラーテーブルはモデリングドメインの温度を示しており, 明るい色は高い温度に対応しています. 屈折率の変化は温度の変化に比例するため, 温度勾配が最大になる光線経路はより顕著に湾曲します. 光線経路に沿った色は, 光線の瞬間的な群速度の大きさを示し, 最大の大きさは赤で示され, 最小の大きさは青と紫で示されます. このシミュレーションは, デモンストレーションの目的で誇張されており, レーザー系のポンピングされたレーザー結晶で何が起こるかを概念的に示しています.

加熱された材料の中を曲がった経路をたどる光線のCOMSOLモデル.

単一の光線は, 加熱された材料を通る湾曲した経路(誇張された)をたどります.

波面光線追跡法

COMSOL Multiphysics の光線追跡法は, 瞬間的な光線位置qと波数ベクトルkの成分について, 連立した1次常微分方程式(ODEs)の組を数値的に解きます. これらの連立方程式は, 古典力学におけるハミルトニアンの定式化に類似しています.

\frac{ \partial {\mathbf q} }{ \partial t} = \frac{ \partial \omega( {\mathbf k} )}{\partial {\mathbf k} }
\frac{ \partial {\mathbf k} }{ \partial t} = -\frac{ \partial \omega( {\mathbf k} )}{\partial {\mathbf q} }

ここで, 角周波数ωは古典力学のハミルトニアンHにおけるが位置変数に相当します(Ref. 4).

屈折率が均一である場合, 上記のハミルトニアンの定式化は, 次のように, 一定の速度と光線の方向を説明する式に還元されます.

\frac{ \partial {\mathbf q} }{ \partial t} = \frac{ \partial (c|{\mathbf k}|/n)}{\partial {\mathbf k} } = \frac{c {\mathbf k} }{n|{\mathbf k}|}
\frac{ \partial {\mathbf k} }{ \partial t} =0

界面で屈折率の不連続性がある場合, COMSOL Multiphysicsは, スネルの法則を使用して屈折光線の方向を数値的に計算します. この定式化により, COMSOL®ソフトウェアは, 均一な屈折率を持つ特殊なケースだけでなく, レーザー工学における熱レンズ(直線光線と曲線光線の両方)などのより一般的なケースも計算できます. これは時間ステッピング法であり, 結果は伝播中のさまざまな時点での光線を示していることに注意してください.

表面の離散化とジオメトリ形状次数

幾何学的形状の例として, 次の解析式で記述される回転対称の非球面レンズ表面のサぐを考えてみます.

(1)

z(r) = \frac{cr^2}{1+\sqrt{1-(1+k)c^2r^2}}+A_1r^2+A_2r^4+ \cdots

ここで, z, \ c,\ k, \ A_i \ (i=1,2,\cdots) はそれぞれ, レンズ表面のサぐ, 曲率, 円錐定数, および偶数多項式の係数です. およびr^2=x^2+y^2, ここで, xおよび yは光軸を横切る2つの方向です.

COMSOL Multiphysicsでは, レンズとミラーの形状を解析関数または本格的な3DCADモデルとして入力できます. 解析的に与えられた形状は, パラメトリック曲線またはサーフェスとして入力され, 高次の不均一有理スプライン(NURBS)で近似されます. 次のステップでは, サーフェス(またはカーブ)を有限要素メッシュで近似します. 有限要素メッシュは, 多項式を使用して形状をさらに近似します.

ユーザーインターフェースでは, これらの多項式の最大次数は, モデルコンポーネントの設定ウィンドウのジオメトリ形状次数によって決定されます. これらの近似の精度は, メッシュを改良したり, ジオメトリ形状次数を増やしたりすることで向上させることができます. 次に, 光線は, 交点を予測し, スネルの法則を適用する目的で, 個々のメッシュ要素と相互作用します.

ジオメトリ形状次数のデフォルトオプションである自動は, ほとんどの場合, 2次多項式を使用してサーフェスを近似しますが, 反転メッシュ要素の作成を回避するために線形形状次数が使用されることもあります. ほとんどのアプリケーションでは, ユーザーは次数を7次多項式に増やすことができます.

COMSOL Multiphysicsのジオメトリ形状次数設定のスクリーンショット.

ジオメトリ形状次数の設定ウィンドウ.

ここで, 前述した非球面レンズについて考えてみましょう. 離散化を見やすくするために, ここではレンズを非常に粗くメッシュ化しています. 次の図では, 線形化されたメッシュ要素は灰色の線で示され, 曲面の数値表現は黒い輪郭で示されています.

左図は, デフォルト(自動, つまりほとんど2次関数)の設定の場合の光線追跡結果で, 右図は1次ジオメトリ形状次数の場合のものです. 光線の色は, 波動ベクトルのy成分(縦軸)を表しています. 光線がどこで屈折するかは, 色が変わる位置を探せばわかります. ジオメトリ形状次数をどのように選択しても, サーフェスの形状は, 上記の非球面の公式では書けなくなります. つまり, ジオメトリ形状次数の設定によって次数が決まる, 区分連続多項式の集合となります.

2つの異なる形状次数に対する非球面レンズ表面のモデル.

非球面レンズの表面と, 2次(左)および1次(右)のジオメトリ形状次数を使用したその近似表面形状.

光線と有限要素メッシュ

COMSOL Multiphysics の光線追跡アルゴリズムが有限要素メッシュを使用して, サーフェスで光線を反射および屈折させる理由はいくつかあります. 前述のように, これは, 屈折率分布が別のフィジックスインターフェースによって計算された場の関数である不均一材質を介して光線を追跡できるようにすることで, マルチフィジックスシミュレーションツールとしての COMSOL®ソフトウェアの強みを強調しています. レンズ表面の解析式が利用できない可能性がある別の状況は, 熱膨張の結果として, またはレンズに対する他の外力のために表面が変形する場合です. 結論として, 光学系と連成されたほとんどすべてのフィジックスは有限要素法を使用してモデル化されているため, 連成は有限要素メッシュを介して行われます.

光線追跡アルゴリズムは有限要素メッシュを使用するため, COMSOL Multiphysics での光線追跡の最小要件は, 光線を反射または屈折するすべての境界に境界メッシュが必要であるということです. 境界メッシュ要素は, 光線が反射または屈折するときに必要な法線方向や主曲率などの界面特性を計算するために使用されます. 境界メッシュは, 光線が境界に到達する瞬間を予測するメッシュ検索アルゴリズムにも情報を提供し, 強度や位相など, 光線経路に沿って計算される量をより正確に計算できるようにします.

場合によっては, 光線が伝播するドメインにもドメインレベルの有限要素メッシュが必要です. たとえば, 不均一材質を含むドメインは常にメッシュ化する必要があります. 屈折率が温度に依存する場合, 材質内の光線の瞬間的な位置での屈折率を計算するには, 温度場から値を照会する必要があります. 温度は, メッシュ要素で定義された有限要素基底関数を使用して, メッシュ要素内の光線の正確な位置で補間できます. 次に, この温度値を屈折率の温度依存式に代入できます.

光線が材質を伝播する場合も,ドメインレベルのメッシュが必要です. 材質に照射された光線パワーは, メッシュ要素で定義された区分的不連続関数を使用して離散化されるため, メッシュは, 各光線の全熱源への寄与を有限範囲に制限する必要があります.

メッシュ解像度が解に与える影響

これまでのところ, COMSOL®ソフトウェアの光線追跡は, 各サーフェスを有限要素メッシュに基づく区分的連続多項式として扱う近似解であることを学びました. したがって, 解の精度はメッシュの解像度に依存します. より高いジオメトリ形状次数とより細かいメッシュにより, 離散化誤差が減少します. 前の図の右側のプロットでは, 一部の光線が同じ境界要素を通過するため, 誤って同じ波数ベクトル(同じ色の線)が生成されます. これは, ジオメトリ形状次数が低く, メッシュが非常に粗い場合に視覚的に確認できる近似誤差の1つのタイプです.

次のプロットは, 前の例の同じレンズによる最初の屈折後の屈折光線の単位波数ベクトルのy成分を示しています. 横軸は, y-軸の初期光線位置を表します. 求解されていないメッシュによって発生するエラーは, ジオメトリ形状字数が1次である場合に最も重要です.

2次および1次形状次数のy成分の2つのプロット.

さまざまなメッシュ要素サイズ(さまざまな色で表示)を使用した, 2次ジオメトリ形状次数(左)と1次形状字数(右)の単位波数ベクトルの y成分のプロット.

スポットダイアグラムに対するメッシュ細分化の影響

メッシュの細分化のレベルが光線軌道の精度に影響することを確認したので, スポットダイアグラムの精度にも影響するのは当然のことです. 初期光線の位置と方向の分布が対称で規則的である場合, 正確な光線追跡法の結果は通常, ほぼ完全な対称性と連続性を示します. 非球面レンズのスポットダイアグラムは常に光軸に対して対称であり, 光線の収差は常に滑らかで連続的です. 有限要素ベースの波面法の場合, メッシュが元のレンズジオメトリの対称性を継承している場合にのみ, 結果は対称になります.

光線の位置と方向の初期分布が対称である対称光学系の場合でも, COMSOL Multiphysics での光線追跡の結果は, メッシュが非対称である場合, わずかに非対称な結果をもたらす可能性があります. ただし, メッシュを細分化すると, 対称面の両側にある小さなメッシュ要素は次第に同じ形状に近づくので, これらの対称性からの逸脱は徐々に減少します. これは, メッシュの離散化誤差のもう1つの兆候です. 構造化メッシュを使用することで, 解の対称性を向上させることができます. ただし, 対称解は必ずしも解がより正確であることを意味するわけではないことに注意してください.

次の図は, スポットダイアグラムと平凸レンズのRMSスポット半径のプロットに対する最大メッシュ要素サイズの影響を示しています. この集束レンズは, 回折限界のスポットを持つように設計されています. メッシュサイズが予想どおりに減少すると, RMSスポット半径が特定の値にどのように収束するかを確認できます. よく見てみましょう. 最も粗いメッシュのスポットは, 回折限界よりも大きくなっています. また, それは非常に非対称であり, 交点は非常に広い領域に分布しています. このメッシュサイズは, 熱や機械的な問題などの他のフィジックスにとっては特に粗いというわけではありませんが, デフォルトの自動ジオメトリ形状次数を使用した光線追跡には粗すぎて, 正確なスポットダイアグラムを作成できません. 最後の2つのメッシュサイズだけが, 集中して対称に見えるスポットを提供します.

2次形状次数と3次形状次数のスポットダイアグラムを視覚的に比較.

最大メッシュサイズに応じた, 2次および3次ジオメトリ形状次数間のフォーカスでのスポットダイアグラムの比較. 赤い円は,0.66umの波長でのレンズの回折限界サイズを示しています.

COMSOL Multiphysicsで作成したRMSスポット半径の両対数プロット.

回折限界サイズと比較した, 2次および3次ジオメトリ形状次数の最大メッシュサイズの関数としてのRMSスポット半径の両対数プロット.

COMSOL Multiphysics®の光線追跡法によって生成された解データの解釈

ここで, 解データの形式が, 上記の従来の(定常)解析的光線追跡アルゴリズムの形式とどのように異なるかを見てみましょう.

以下は, COMSOL Multiphysicsで得られた時間依存の解を, 前の例の同じ集束レンズの標準的なシーケンシャル光線追跡平面間伝搬法からの一般的な出力と比較する例です. 平凸レンズは光線をスポットに集束させます. 従来の平面間伝搬法でおなじみの結果では, すべての光線が面で終端します. 一方, COMSOL Multiphysics では, 各光線は, 指定された時間ステップにあるはずの位置に描画されます. したがって, 一般的に, 光線は同じ平面に配置されません.

これが, COMSOL Multiphysics での光線の軌跡の結果がほとんど常にある程度の曲率を示す理由です. この曲率は, 幾何光学領域の観点からは波面に他なりません. この例では, 波面はレンズが光線に導入する2次位相関数によって湾曲しています. すべての光線を最終面に描画するには, すべての光線が最終面に到達するまで待つ必要があります.

COMSOL Multiphysics における時間依存型および平面間伝搬法光線追跡の例.

COMSOL Multiphysics® での時間依存光線追跡(左)と標準のシーケンシャル平面間伝搬法光線追跡(右)のシミュレーション結果.

生の波面を見るのは便利ですが, シミュレートするすべての光線が, 量を評価する平面を通過していることを確認する必要があります. 次の連続したスポットダイアグラムは, さまざまな計算時間での上記の例のスポットダイアグラムを示しています. COMSOL Multiphysics では, スポットダイアグラムはポアンカレマップタイプを使用して生成されます. 波面が収束しているため, 外側の光線が最初に平面に到達し, 内側の光線が後に到着します. したがって, 面内評価は, 選択した時間によって異なる場合があります.

異なる時間での5つのスポットダイアグラム.

さまざまな時間で生成されたスポットダイアグラム.

関心のある平面内のすべての光線を確実に評価するために, COMSOL Multiphysics にはattimemin attimemaxと呼ばれる演算子が装備されています. たとえば, 次の式は, x = 0.1 m:でのyz平面のRMSスポット半径を示します.

sqrt( gop.gopaveop1( attimemin(0,0.4[ns],(qx-0.1[m])^2, qy^2+qz^2) ) )

演算子gop.gopaveop1()は, すべての光線の括弧内の式の平均を取ります. この演算子への引数,

attimemin(0,0.4[ns],(qx-0.1[m])^2, qy^2+qz^2)

attimemin演算子を呼び出します. これは4つの引数を取ります.

attimeminの最初の2つの引数は, 式の最小値を計算する期間の開始と終了を定義します. 光線が平面を横切る時間がこの範囲内にあることを確認する必要があります. 3番目の引数は, 最小化する式です. この3番目の引数の最小値が光線に対して検出されると, 4番目の引数の値が返されます. attimemin 演算子は, 保存された解の時間の間の補間を使用するため, この交点が解に保存された時間ステップと一致しなくても, 各光線と平面との交点を正確に取得できます.

つまり, 前の式は次のようになります.

平面x = 0.1からの光線の距離が最小になる0〜0.4 nsの時間で, 半径座標の値のすべての光線の平均の平方根を計算します.

次の例は, 光線が境界で反射または屈折したときに, 離散時間ステップでの解の保存が軌道プロットにどのように影響するかを示しています. 通常, 各光線が各サーフェスと相互作用する正確な時間または光路長は事前にわからないため, これらの交差時間は通常, 保存されている解時間と正確に一致しません. その結果, COMSOL Multiphysics の光線追跡の結果は, 次の図に示すように, 境界面または反射壁の近くで不完全な軌道を示しています. この図では, 光線は境界上ではなく境界から少し離れた方向に方向を変えているように見えます. ただし, これは見かけだけで, 計算は正確に行われており, 精度には影響しません. 結果で光線が界面または壁に当たっているように見えない場合でも, COMSOL®ソフトウェアは, 通常, 2次テイラー法を使用して以前に保存された解から光線の位置を推定して, 各光線の正確な相互作用時間を計算します.

壁で反射する光線のCOMSOLモデル.

反射壁での光線の反射. プロット内のすべての光線が壁と交差していませんが, すべての光線の軌跡が正確に計算されます.

おわりに

このブログでは, 光線光学シミュレーションのために COMSOL Multiphysics で使用される独自の時間依存およびメッシュベースの方法について説明しました. また, 特に時間ステップが大きく, メッシュが粗い場合に, 光線追跡の結果で発生する可能性のある異常の種類についても説明しました. 他のフィジックスインターフェースと同様に, メッシュを改良したり, 時間解像度を向上させたりすると, これらの異常な動作が解消され, 正確で収束した解が得られます.

その他のリソース

COMSOLブログで, COMSOL Multiphysicsの光線光学および光線追跡機能の詳細をご覧ください.

参考文献

  1. R.R. Shannon, The Art and Science of Optical Design, Cambridge, 1997.
  2. D.C. O’Shea, Elements of Modern Optical Design, Wiley, 1985.
  3. R. Ditteon, Modern Geometrical Optics, Wiley, 1998.
  4. L.D. Landau and E.M. Lifshitz, The Classical Theory of Fields, 4th ed., Butterworth-Heinemann, Oxford, 1975.

コメント (0)

コメントを残す
ログイン | 登録
Loading...
COMSOL ブログを探索