前回のブログでは, 線形押出し演算子 を紹介し, ソースと行先間の変数マッピングにおけるその使用例を示しました. 前述の通り, この手法はソースと行先がアフィン変換で対応付けられる場合に限られます. 今日は, 非線形マッピングや異なる次元のジオメトリエンティティ間の変数マッピングを扱うために設計された, 一般的な押出し演算子について解説します.
押出し演算子の簡単な概要
行先エンティティ上の点 P_d において, ソースエンティティで定義された別の量に依存する量を計算したいとします. このとき, 後者のソース点 P_s からの後者の量を行先エンティティにコピーする必要があります. 押出し演算子は, ソースエンティティ内のどの点が行先エンティティ内のどの点に対応するかを特定するために使用されます. 言い換えれば, 演算子は点から点への写像を定義します.
マッピングが アフィン写像 である場合なら, ソース内のいくつかのポイントが行先エンティティの点にどのように対応するかを把握すれば十分となります. このようなソースと行先のペアから, 重ね合わせによって一般的な写像を推測できるからです. しかし, 一般にマッピングを行うには, 写像の数学的表現を与える必要があります. これには, ソース点 P_s を P_d の関数として明示的に定義することも, あるいは P_d と P_s の間の陰的な関係として記述することもできます.
COMSOL Multiphysics における一般押出し演算子の使用
線形押出し 演算子を使用する場合, 十分な点 (基点) に対する写像を視覚的に指定すると, COMSOL Multiphysics が残りの点の変換方法を自動的に決定します. 一般押出し 演算子の場合は, 行先における任意の点に対する写像の数学的記述を与えます.
まず, 線形押出し演算子を一般押出し演算子で再現する方法を考えてみましょう. その後, 一般押出し演算子の使用が必須となる事例を検討します.
アフィン関係に対しては, 一般押出し演算子は線形押出し演算子の代替として使用することができます. 一般的な非線形マッピングにおいては, 一般押出し演算子が必須となります. 一般押出し演算子を追加するには, 定義 > 非ローカルカップリング > 一般押出し と移動します.
例 1
前回の線形押出し演算子に関するブログ記事 では, ソースドメインの点 1, 4, 2 を行先ドメインの点 1, 5, 3 に対応させるアフィン写像について検討しました. 下図をご覧ください. この図形内の2つの円は原点を中心とし, 半径はそれぞれ 1.0 と 1.5 となっています.
任意のアフィン変換は, 線形変換と平行移動の和として表すことができます. したがって, 次の関係が成り立ちます
次に, 定数 a,b,c,d,e,f を求める必要があります. ソース点 (0, 0), (1.0, 0), (0, 1.0) は, それぞれ行先点 (0, 0), (1.5, 0), (0, 1.5) に対応することから, 次の式が得られます
行先の任意の点 (x, y) に対して, これに対応したソース点座標を求める方法が分かったことになります. この式の右辺 (添字なし) を, 一般押出し設定ウィンドウの行先マップに入力します.
一般押出し演算子を用いて構築された線形マッピング.
例 2
次に, 一般押出し演算子を使用して2D軸対称コンポーネントから3Dコンポーネントへデータをコピーする方法を探ります. なお, ソース点と行先点は空間内の同一の点に対応するものとします. 軸対称な伝熱境界条件と材料特性を用いた 熱膨張 を考えます. もし構造境界条件が軸対称でない場合は, あるコンポーネントで軸対称伝熱解析を実行し, 2D軸対称ドメインの温度を, 別のコンポーネント内の構造解析用の3Dドメインへマッピングすることで時間を節約できます.
マッピングを構築する際には, 次の問いが重要となります. 行先地点の座標が与えられた場合に, どのようにしてソース点へ到達するのか? 今回の場合, 関係は次のように与えられます.
例 1 にもあったように, この式の右辺を行先マップに入力します.
一般押出し演算子を使用して, 2D軸対称ドメインから対応する3Dドメインへデータをコピーします. 軸対称コンポーネントの場合, 結果ノードでは, 変数は回転2Dデータセットを用いて3D表示できます. ただし, 3Dコンポーネントのフィジックスノード (例: 熱膨張) で2D軸対称コンポーネントの変数を使用したい場合は, 一般押出し演算子を利用する必要があります.
演算子 genext1 は, 3Dコンポーネント comp2 内で定義されておらず, T も同様となります. このため, 2D軸対称コンポーネントの温度を3Dコンポーネントの入力として使用する場合は, comp1.genext1(comp1.T) を使用する必要があります. この方法により, 2番目のコンポーネント内で genext1 と呼ばれる押出しや別の演算子, あるいは T という名前の変数が存在していた場合の混乱を回避できます.
ここでは, 線形押出し演算子は使用できないことに注意してください. ソースオブジェクトと行先オブジェクトの次元が異なるため, アフィン変換は不可能となっています.
最初の2つの例では, 設定ウィンドウのソースセクションにある ソースマップ使用 チェックボックスはオフのままでした. COMSOL Multiphysics は, 最初のケースでは x と y を自動的に入力しています. 2番目のケースでは r と z を自動的に入力しています. このチェックボックスがオフの場合, COMSOL Multiphysics はソースの各座標が行先座標の関数として陽的に与えられていると仮定します. しかし実際には, 陽的な式が得られない場合がよくあります.
次に, 一般押出し演算子を使用して陰的な関係を指定する方法を見ていきましょう.
例 3
\frac{y}{d} =(\frac{ x}{d})^2 で示される二次元の放物線は, 辺の長さ d の正方形のドメイン内に存在します. 私たちは, この曲線 (下図で青で表示) 上のデータを正方形の異なる部分へ写像する演算子を構築する必要があるとします. この放物線がソースとなり. 放物線上の点から正方形の点へコピーする演算子が必要ということになります. この際, 行先の原点からの距離が, 原点とソース点間の放物線の長さに等しくなるようにします.
微積分を少々行うと, 原点とソース点 (x,y) の間における放物線の弧長が求められます.
したがって, ソース点と行先点の関係は
もしも次のような形式の陽的なソース・行先マッピングが必要だとすると
まず式 L=\frac{x_s}{2}\sqrt{1+4(\frac{x_s}{d})^2}+\frac{d}{4}\ln(2\frac{x_s}{d}+\sqrt{1+4(\frac{x_s}{d})^2}) を解き, x_s を L の式で表さなければなりません. これはまったく面白くありませんね!
まさにこの理由から, COMSOL Multiphysics では, 行先座標とソース座標の間の陰的な関係を, 行先マップとソースマップという2つのマッピングを用いて指定することが可能です. ここでは, 次を満たす T_d と T_s を与える必要があります.
一般押出し演算子において, ソースマップと行先マップを使用して, ソース座標と行先座標間の陰的な関係を定義します.
COMSOL Multiphysics は, ソース座標を特定するために必要な手順である逆変換 T_s^{-1}(T_d(x_d,y_d)) を自動的に処理します. なお, 逆変換が存在するためには, ソースマップが単射 (一対一) である必要があることに注意してください. 実際には, COMSOL Multiphysics はソースマップの逆関数に対する解析式を構築するのではなく, 代わりに各行先点において, まず T_d(x_d,y_d) を評価し, この評価値が T_s(x_s,y_s) と一致するようなソース上の点を見つけるメッシュ探索を実行します. 1対1のソースマップにより, 特定の行先点に対して探索が返すソース点は最大でも 1 つとなります.
上記の一般押出し設定ウィンドウでは, 行先マップ と ソースマップ の下にあるラベルが, x 式 と y 式 ではなく, x^i 式 と y 式 と表示されています. これは, x^i および y^i が, ソースと行先の関係を陰的に定義するためための最初の式のペアと2番目の式のペアに対するインデックスであるためです. これらは, 必ずしもソースまたは行先の x または y 座標に関係しているわけではありません. これらのインデックスは中間メッシュの座標のようなものであり, 一般押出し演算子は同じ中間座標を持つソース点と行先点を対応付けます. この例では, 正方形領域内の任意の行先点を放物線上のソース点に一意に関連付けるには, 式は一つで十分です. したがって, 2行目の y^i 式 は空白のままとしています.
この一般押出し演算子が変数をどのようにマッピングするかを確認するために, 左右のエッジがそれぞれ 300 K と 400 K である平面定常熱伝導問題を考えてみましょう. 上面と下面は断熱されており, 熱源はありません. 温度は x と共に線形に変化します. 下のグラフから, 右側の arcext(T) のプロットが放射状に変化している理由がわかるでしょうか?
左: 温度は左から右へ線形に変化. 中央: 放物線上の温度分布. 右: 放物線からドメインへマッピングされた温度. ドメイン内の原点からの距離が等しいすべての点は, 放物線上の同じ点から温度をコピーしています.
すべてをつなぎ合わせる
これまでに学んだことを応用するために, COMSOL Multiphysics の 電流 フィジックスインターフェースを使用してダイオードモデルを構築しましょう. 押出し演算子を用いることで, 理想的な p-n 接合の両側に通常の電流密度境界条件を設定できます. 下の図に示すように, 両側をそれぞれ 1 と 2 として識別します. 接合部では, 電流-電圧 (I-V) 関係のショックレーダイオード方程式を用います. パラメーター J_s, q, k, T は, それぞれ飽和電流密度, 電子電荷, ボルツマン定数, 温度を表します.
押出し演算子を使用して, 接合部の反対側にある電位にアクセスできます.
面 1 に法線電流境界条件を実装するには, 面 2 の電位 V_2 へのアクセスが必要であり, 同様に面 2 では, 接合部の反対側の電位 V_1 へのアクセスが必要です. したがって, 2つの押出し演算子が必要となります. 下図に示すように, 接合部の両側をそれぞれの押出し演算子のソースエンティティとします.
両方のケースにおいて, 同じ x 座標を持つ点間のマッピングを含みます. ソースエンティティが異なるため, 2つの演算子が必要です.
次に, フィジックスノードにてこれらの演算子を用いて境界条件を実装します. 上面の境界条件を以下に示します. V は上面の点における電位を表し, genext2(V) は真下の下面電位を表すことに注意してください.
接合部の反対側の点における電位を参照するために, 一般押出し演算子を使用.
接合部の下面側にも同様の境界条件を適用します. エッジ 3 に適用された 法線電流密度 2 ノードの対応する法線電流密度は -Js*(exp((V-genext1(V))/kTbyq)-1) となります. ここでは, V は下側の点における電位を, genext1(V) は真上の上面での電位を表します.
実は, どちらの面の設定かに関わらず, 電圧差を genext2(V)-genext1(V) と表して手間を少々省くことができます. 分かりやすさのため, 今回はこの設定を避けています.
装置の下部に電圧端子を配置し, 上部を接地すると, 以下の結果が得られます.
押出し演算子を使用して, 同一コンポーネント内または異なるコンポーネントを跨いだ点同士のカップリングを作成できます. ここでは, ジオメトリの細いギャップがダイオードの p-n 接合に対応しています. ギャップを横切る電流密度を計算するため, 押出し演算子を用いてギャップの一方の電位をもう一方からアクセスします.
他の非ローカルカップリング
押出し演算子を使用すると, ソース点と行先点間の点ごとの対応を設定できます. 時には, ソースとなる線, 表面, または体積に対する積分, 平均, 最大値, 最小値が必要となる場合があります. そのような場合には, 射影, 積分, 平均, 最大, 最小といった非ローカルカップリングが使用可能となります. 射影演算子の使用方法については, 以前のブログ記事 をご参照ください.
おわりに
今日は, 一般押出し演算子を用いて, シミュレーションドメイン内の変数をある部分から別の部分へコピーするマッピングを作成する方法について触れました. p-n 接合の例で示した通り, これらの演算子は既知の量を単純にコピーするだけでなく, 未知変数間の非局所的なカップリングを作成するためにも使用できます. この手法は構造接触や伝熱における表面間放射など, 他種の解析においても有用です. COMSOL Multiphysics には, こうした物理現象に関連する組み込み機能が用意されています.
もしシミュレーション対象の非局所的なカップリングが COMSOL Multiphysics の標準機能に含まれていないとしても, 本日学習した手法を用いると実装することが可能です. ご質問がございましたら, お気軽に お問い合わせ ください!
関連文献
他のシチュエーションでの一般押出し演算子の使用例については, 以下のブログ記事をご覧ください:
コメント (0)