時空間の離散化による方程式ベースのモデリング

2020年 9月 24日

COMSOL Multiphysics® ソフトウェアの主な強みの 1 つは, 計算モデルのほぼすべての式を変更できることです. この機能は慎重に使用する必要がありますが, 素晴らしいことができます. 並進のある 2 次元の熱モデルを取り上げ, いくつかの材料特性をゼロに設定して, これが 1 次元の過渡モデルを解くのと似ていることを観察してみましょう. これにより, ある種の最適化問題を簡単かつ迅速に実行できることが確認できます.

2つの空間次元から1つの空間次元と1つの時間次元へ

時空離散化問題を示す1D過渡モデルに類似した2D定常熱モデルを示す概略図.

1次元過渡モデルに類似した2次元定常熱モデルの模式図.

ここでは, 非常にシンプルな長方形の2次元熱伝導モデルを考えてみましょう. 体積加熱はなく, 対流項がある場合の温度, Tに対する定常(時間不変)の伝熱支配方程式, \mathbf{u}を以下のように解きます.

\rho C_p \mathbf{u} \nabla T + \nabla \cdot \mathbf{k} \nabla T = 0

ここで\rhoは材料の密度, \rhoは比熱, そして\mathbf{k}は熱伝導率で, ここでは次のような対角線上の行列となっています.

\mathbf{k}= \begin{bmatrix}k_{xx} & 0\\0 & k_{yy}\end{bmatrix}

すべての項を展開して支配方程式をもう少し詳しく書いてみるとよいでしょう.

\rho C_p \left( u_x \frac{\partial T}{\partial x}+ u_y \frac{\partial T}{\partial y}\right) + \frac{\partial}{\partial x} \left(k_{xx} \frac{\partial T}{\partial x}\right)+ \frac{\partial}{\partial y} \left(k_{yy} \frac{\partial T}{\partial y}\right)= 0

ここで, 面白いことをします. ここでは, 速度ベクトルが純粋に+y方向であると仮定し, u_x = 0とし, 対角線上の熱伝導度テンソルのy成分をゼロにしますk_{yy} = 0. 上の式は次のようになります.

\rho C_p u_y \frac{\partial T}{\partial y} + \frac{\partial}{\partial x} \left(k_{xx} \frac{\partial T}{\partial x}\right) = 0

熱伝導度テンソルの項をゼロにすることは非常に珍しいことのように思われますが, これによって興味深いメリットが得られます. これは, 体積加熱や対流項を含まない1次元の過渡的な熱伝導支配方程式を書き出すことで明らかになります.

\rho C_p \frac{\partial T}{\partial t} + \frac{\partial}{\partial x} \left(k_{xx} \frac{\partial T}{\partial x}\right) = 0

上の2つの式の最初の項はほとんど同じに見え, 速度項を正しく選択すれば, 2つの式は同じになることに注意してください. さて, ここで非常に興味深いことが起こります. 2 次元の定常モデルは, 1次元の過渡モデルと同じように見えます.

実際, 方程式は同じですが, それを数値的に解くときには, いくつかの点に注意しなければなりません. まず, 2Dドメインの下部境界で固定温度条件を適用することは, 過渡的な1Dケースの初期条件と類似しています. 第二に, このような2次元モデルでは, マッピングされたメッシュ(Y軸に沿ったもの)を使用することに意味があり, Y方向の要素数は, 固定の時間ステップと考えることができます. これは, デフォルトのアダプティブタイムステッピングを使用するタイムドメインモデルとは全く異なります. また, 2Dモデルには(デフォルトで)数値安定化項が含まれているため, どちらかのモデルがメッシュや許容範囲を高度に細分化しない限り, 両者の結果は微妙に異なるものとなります.

次の質問は, 当然のことながらこれで何ができるのか?ということです. 下図に示すように, 材料のスラブを加熱する古典的な1次元熱伝達モデルを見てみましょう. 一見シンプルですが, 工学的に興味深い問題の多くは, この古典的なケースに集約されます. 左側の表面から材料を加熱すると, 材料の体積全体が徐々に, しかし不均一に熱くなります. また, 放射冷却と, 一定の熱伝導率に近似した左側からの自由対流が発生します.

過渡加熱の問題がどのように1Dモデルになるかを示す図.

材料のスラブの過渡的な加熱は, 1Dモデルになります.

この問題は COMSOL Multiphysics® で 2D の矩形領域を設定し, フィジックスインターフェース伝熱 (固体および液体)を使用して解くことができます. このインターフェースでは, 固体機能に並進運動サブ機能を, 下部境界に温度境界条件を, 左側に熱流束および表面から周囲への放射機能を適用して対流冷却および放射冷却をモデル化し, さらに熱流速条件を適用して熱負荷を与えることができます. 時間軸上に一様な熱負荷を与えた場合の結果例を以下に示します.

1D材料スラブモデルの過渡加熱の2Dシミュレーション結果とメッシュ.

1次元スラブの過渡的な加熱を表現した2次元静止モデルの結果例. メッシュも表示されています.

ここで, この加熱プロセスを最適化したいと考えてみましょう. シミュレーション終了時にスラブ全体の温度が目標温度にできるだけ近づくことを目的に, 熱負荷を時間経過とともに変化させ, ピーク温度が規定値を超えないようにするという制約も考えてみましょう. 1次元のスラブを2次元にしたことで, この最適化問題は非常に簡単に設定できることがわかりました. では, その方法を見てみましょう.

最適化問題の定義

定義トポロジ最適化にある密度モデル機能は, 時間経過に伴う熱流束を制御する変数を定義するために使用されます. この機能では, 加熱なしと最大加熱に対応する0と1の間に拘束された境界上の場を定義します. この場の離散化は線形です. つまり, 時間軸に沿った各要素上で熱負荷は線形に変化するため, 熱負荷はy方向のメッシュサイズで決まる時間にわたって変化することになります. これは理想的ではありません. 熱負荷が時間ステップと同じ速さで変化していることになります. これは理想的ではありません. 熱負荷は離散化よりもゆっくりと時間的に変化するようにしたいのです. そこで, もう一歩踏み込んで, ヘルムホルツ・フィルターを導入することにしました. これは, トポロジ最適化の分野では一般的なアイデアで, 設定のフィルタリングセクションで利用できます. フィルター半径は, 単にメッシュサイズよりも大きくする必要があり, スムージング時間を表します. このフィルタリングされた制御変数場の名前は, dtopo1.thetaで, 入射熱流束を乗算します.

トポロジ最適化密度モデル機能の設定ウィンドウのスクリーンショット.

トポロジ最適化密度モデル機能のスクリーンショット. スラブの加熱面に沿って, 時間を表す軸全体で定義されています.

積分される式は, 計算された解答TT_target(到達したい温度)に基づいています. この境界線上の計算された温度と目標温度の二乗差の積分である目的関数を定義することで, 最終温度が目標にできるだけ近いときに最小値を持つ微分可能な関数が得られます. 目的関数の名前 comp1.objで, これは最適化のスタディステップから参照します.

COMSOL Multiphysics の目的関数プローブ機能の設定ウィンドウのスクリーンショット.

目的関数プローブ機能のスクリーンショット. 上部の境界線上に定義され, 最小化する目的関数式を定義しています.

最適化インターフェースの中で定義された, 次のような制約があります.

comp1.constr < 1

この式, comp1.constrは, 上のスクリーンショットに示されているように, グローバル変数probeを介して定義されています. この式は, 最初に以下の平均値を取ります.

\left( \left( \frac{T}{T_{max}}\right)^{p_{exp}} \right)

加熱された境界線上で, p_{exp}-th rootを取ります. これがPノルムです. p_{exp}が無限大に近づくと, この制約は, 加熱された境界に沿って, 温度が最大温度以上にならないように強制することになります. ここで, 数値計算上の理由から, p_{exp}は無限大ではなく, この制約を正確に満たすことはできませんが, p_{exp}が十分に大きければ, この制約はほぼ満たされます. しかし, p_{exp}が十分に大きければ, 制約はほぼ満たされます. p_{exp} を大きくしすぎると, 非常に非線形な制約関数を導入することになり, 収束が遅くなるだけでなく, 数値的なオーバーフローの問題が発生する可能性があるため, この値はモデルのチューニングパラメーターと考えてください.

制約式の定義に使用されるプローブのグローバル変数プローブ設定のスクリーンショット.

制約条件式を定義するプローブを示すスクリーンショット.

最適化モジュールの一部として提供されるスタディブランチ内の最適化スタディステップでは, 目的関数, 制約, および制御パラメーターを導入することができます. 目的関数と制約は, 上述のプローブのようなグローバル変数にすることができます. このケースの設定を見てみましょう.

最適化スタディステップの設定ウィンドウのスクリーンショット.

最適化ソルバーのタイプ, 目的関数, 設計変数, 制約条件を定義した, スタディブランチ内の最適化スタディステップのスクリーンショット

上のスクリーンショットは,スタディブランチ内の最適化 スタディステップを示しています. この機能では, 最適化ソルバー, 目的関数, 設計変数, および制約条件を定義しています. 設定の一番上では, SNOPTソルバーが使用されていることがわかります. このソルバーは, 目的関数と制約条件が, 設計変数に対して微分可能であることを想定しています. 最適性の許容範囲モデル評価の最大数は, ソルバーが最適解を見つけるために何回反復するかを決定します. 次の3つのセクションでは, 目的関数, 制御変数とパラメーター, および制約条件を定義します.

目的関数, 設計変数, 制約条件が定義されれば,あとはモデルを解くだけです.

設計変数, フィルター処理された設計変数, および温度の結果をプロットしたグラフ.

設計変数,フィルタリングされた設計変数, および温度を, 時間を表す軸に沿ってプロットしたもの.

シミュレーションの最終時点での温度と材料スラブの目標温度を比較するプロット.

最終時刻の温度(青)と目標温度(黒)を比較したもので, 1次元領域上の温度場が目標温度に非常に近いことを示している.

このモデルには,熱伝導率の材料非線形性が含まれていますが,上の図に示すような解に最適化するのに約 1 分かかりました.この図では,時間経過に伴うスラブの加熱面の温度,フィルタリングされた変数とフィルタリングされていない設計変数,およびスラブの断面を通る最終温度をプロットしています.特に興味深いのは, 最適化ソルバーが最初にピーク加熱を行い, ピーク温度の制約を超えないように熱負荷を徐々に減少させ, 温度場がスラブ全体で最適に向かって平衡化するとゼロになり, 最後に短時間の加熱を行って周囲への冷却を打ち消すという解決策を見つけていることです. このように, 解決策を定性的に説明することは可能ですが, 具体的な状況に合わせて最適化しないと難しいと思われます.

この時点で, このアプローチの利点は何なのかと疑問に思うかもしれません. 実際のところ, 最適化問題全体を純粋に時間領域で求解することも可能であり, 何が得られるのでしょうか?スピードとシンプルさです.

時空間に再定式化することで, 定常的な最適化問題を解くことができ, 過渡的な最適化問題よりも早く解くことができます. また, ヘルムホルツフィルターを含む内蔵のトポロジ最適化機能を使うことができるので, 任意の, 制約のある, 平滑化された, 時間軸上の強制関数を非常に簡単に設定することができます. では, 概念的な複雑さを除いた欠点は何でしょうか? この手法は, より多くのメモリを使用します. 空間と時間を同時に解くことで, 剛性行列は時間軸に沿ったメッシュに比例して大きくなり, このメッシュはシミュレーションの前に設定する必要があります. しかし, このような2Dモデルでは, 非常に細かいメッシュであっても, 必要なメモリ量はごくわずかです.

2つの空間次元を持つモデルでは, 必要なメモリ容量は大きくなりますが, それでも法外な値ではありません. 2次元のモデルを3次元の時間で最適化する場合の例は, 以下のリンクでも紹介されており, ジュール熱をモデル化するための熱伝導に加えて, 電流のフィジックスも時間経過とともに最適化されています. それ以外は, ここで紹介した技術と同じです. モデリングをお楽しみください.

試してみましょう

時空間モデルを使って, 熱負荷の時間的な最適化を試してみてください. 下のボタンをクリックすると, モデルファイルにアクセスできます.

コメント (0)

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