様々な応力変数や結果表示オプションへのアクセスを提供する COMSOL Multiphysics® ソフトウェアでは, 異なる応力量を最適に評価する方法についてよくご質問をいただきます. このブログでは, これらの問題について詳しく解説します. 応力解析の詳細に入る前にまず, 一般的な結果評価方法についてご説明します.
目次
- 応力の特徴とは?
- 有限要素の基礎知識
- 物理的な連続性
- 要素間の平均化
- 結果はどこで評価されるのか?
- 応力の計算方法
- 特殊なケース: 熱膨張
- 点ごとの状態
- 最大値の評価
- どの応力結果をプロットすればよいのか？
- 最後に
応力の特徴とは?
物理学の多くの分野では, 主な関心の量は場そのものであり, その勾配ではありません. 例えば, 熱伝導の解析では, 通常, 温度が注目され, 詳細な温度勾配や熱流束は二次的な関心事となります. 一方, 構造力学では, 変位よりも局所的な応力や歪みが重要視されることが多くあります. ほとんどの場合, 高い応力は静的または疲労による破壊の原因となります. したがって, 信頼性の高い応力評価が重要となるのです.
有限要素の基礎知識
有限要素法では, ジオメトリが有限要素と呼ばれる小さなパッチに分割され, メッシュを構成します. 各要素内には, 解くべき場の変化に関する仮定が存在します. これは, 形状関数によって近似されます. 最も一般的な形状関数は, 座標の関数として1次または2次多項式です. 場は, スカラー場またはベクトル場です.
すべてのメッシュ要素は, メッシュノードで互いに連結しています. いくつかの流束の平衡考察を用いることで, すべてのメッシュノードにおける場の値に対する連立方程式を設定することができます.
求解される場について, よくあるケースを下表にまとめます:
| フィジックス | 場 | 勾配 | 流束 |
|---|---|---|---|
| 固体力学 | 変位 (ベクトル) | 歪み (テンソル) | 応力 (テンソル) |
| 伝熱 | 温度 (スカラー) | 温度勾配 (ベクトル) | 熱流束 (ベクトル) |
| 拡散 | 濃度 (スカラー) | 濃度勾配 (ベクトル) | 質量流束 (ベクトル) |
| 静電気 | 電位 (スカラー) | 電場 (ベクトル) | 電気変位 (ベクトル) |
有限要素近似によって生成された方程式を解くと, 従属変数 (自由度またはDOF) はすべてのメッシュノードで既知となります. 場はジオメトリ全体にわたって連続的であり, 要素内の任意の点における値は, 要素ノードにおける値と想定される補間多項式 (すなわち形状関数) によって一意に定義されます.
しかし, 場の勾配は, 一般に要素間で連続的ではありません. 要素内では, 形状関数の空間導関数と個々の要素の節点における節点値によって決定されます. したがって, 勾配も各要素内では一意に定義されますが, 要素の境界では定義されません. 次に, 流束は, 勾配の関数として計算されます. 線形の場合, この関係は勾配に特定の材料パラメーターを乗じるだけです.
物理的な連続性
物理的な理由から, 流束はある程度連続的であることを認識することが重要です. ある要素から流出するものは, 次の要素に流入するはずです. 連続でないのは数値表現だけです. しかし, 有限要素法による解は, メッシュを細分化するにつれて, 真の連続解に向かって収束していきます. しかし, 応力解析に携わる者にとって不幸なことに, 勾配の収束は場の収束よりも遅くなります.
ただし, 勾配と流束の連続性の特性は些細なものではありません. 2つの異なる材料の間に境界がある場合, 流束と勾配の一部の成分のみが連続的になります.
一般に, 流束は境界に垂直な方向には連続的でありますが, 接線方向には連続的ではありません. 勾配については, その逆です.
以下の例では, 材料特性の異なる2つの正方形からなる長方形における定常熱伝導を解析しています.
鋼製の正方形とアルミニウム製の正方形で構成される長方形内の温度分布. 温度は, 図に示す値ですべての境界に沿って規定されています.
2つのドメインを隔てる境界に沿って温度勾配と熱流束をプロットすることで, 連続性の特性を可視化することができます.
材料間の境界に沿った温度勾配 (左)と, 材料間の境界に沿った熱流束 (右).
予想通り, x方向の流束と同様に, y方向の温度勾配も界面の両側で同じになります. もし, 固体力学で同じ実験を行った場合, \sigma_\mathrm{x}と\tau_\mathrm{xy}という2つの応力成分と, \varepsilon_\mathrm{y}という1つの歪み成分に連続性が見られたでしょう. 別の見方をすれば, 自由境界で規定できるような磁束の成分は, 内部境界でも連続でなければならないということです.
等価応力や流束ノルムなど, ほとんどの不変量は, 構成テンソルやベクトルのすべての成分が連続的ではないため, 材料の不連続点でジャンプすることにご注意ください.
要素間の平均化
私たちが解析しているほとんどの量は連続的であることが知られているため, 例えば, 要素間の応力を平均化したくなるものです. これは, ほとんどの場合, デフォルトの設定でもあります. 多くの場合, 見た目を改善するだけでなく, 収束した真の解に近くなります.
平均化 (結果表示の用語では平滑化) が適用されているかどうか, またどのように適用されるかを制御することができます. ここで重要なのは, ほとんどの結果表示ノードで利用可能な品質セクションです.
平滑化のオプション.
さまざまな平滑化 (スムージング) オプションは, 下表にまとめられています:
| 選択 | 効果 |
|---|---|
| なし | 隣接する要素間の平滑化なし. |
| 材料ドメイン内 | 平滑化は, 同じ材料ドメインに属する隣接した要素間で行われます. 材料ドメインとは, 最も単純な形では, 同じ材料の割り当てを持つドメインの集合のことです. しかし, フィジックスインターフェースの中には, 独自の定義を実装しているものがあります. そのような例の1つがシェルインターフェースです. 要素は同じ材料と厚さを持ち, 折れ線によって接続されていない場合にのみ, 同じ材料ドメインに属します. これは, 新しいプロットを追加するときのデフォルトのオプションであることにご注意ください. |
| ジオメトリドメイン内 | 平滑化は, 同じジオメトリドメインに属するがドメイン境界を越えていない要素間で実行されます. |
| 全体 | 平滑化は, 隣接するすべての要素間で実行されま. |
| 式 | ユーザー定義のブール式がゼロ以外の値に評価されると, 平滑化が実行されます. |
平滑化を適用する際, 平滑化の閾値を使用することもできます. これは, 特定のメッシュノードで隣接する要素間の値の差が, 平滑化が無効になるまでにどの程度まで大きくなるかを制限するものです. 重要な不連続性を隠すことなく, プロット全体が美しく滑らかになるという妥協点を得られるということです.
COMSOL®ソフトウェアのバージョン 6.0 以降, 構造力学インターフェースで作成されるデフォルトの応力プロットは, この機能を使用しています. デフォルトでは, 閾値は 0.2 に設定されていますが, 必要に応じてこの値をカスタマイズすることができます.
固体力学インターフェースのデフォルトプロットの品質セクション. デフォルトの閾値0.2では, 最大 20% の差異を受け入れて平滑化が実行されます.
下の図では, 単一材料の 2D 固体力学モデルにおける応力プロットについて, さまざまなタイプのの例を示しています. 効果を強調するために, 1次三角形要素が使用されています. これは低性能な要素で, 各要素内で歪みや応力が一定となるものです. 図を見るときに注意したいのは, 次のことです.
- 平滑化が行われない場合 (左下プロット), メッシュの各要素がはっきりと見えます.
- すべての要素間で平滑化を実行した場合 (中央下のプロット), 平滑化を実行しない場合よりも見た目がずっと良くなっています. 細かいメッシュと2次形状関数を用いて生成された「正しい」高解像度解 (右下のプロット) との全体的な類似性は, それほど悪くはありません. しかし, 応力の高い領域の細部は正確に表現されていません.
- この例では1つの材料しか使用していないため, “全体” (中央下のプロット) と「材料ドメイン内」のオプションに違いはありません.
- >”ジオメトリドメイン内” オプションを使用すると (右上のプロット), 各ドメイン内で平滑化が適用され, ドメイン間で応力がジャンプしているのがわかります. もし, 異なる材料が様々なドメインに割り当てられていた場合, “材料ドメイン内” オプションを使用したプロットは次のようになります.
- >“材料ドメイン内” オプション (左上のプロット) では, 構造力学のデフォルト設定が使用されています. 応力ジャンプは, 応力の解像度が低いところでははっきりと見えますが, 勾配があまり顕著でないところでは等高線が滑らかになっていることにご留意ください. このタイプのプロットでは, 解像度が不十分な領域が目立つ一方で, 全体的に滑らかな概観を得ることができ, 両方の長所を兼ね備えています.
さまざまな種類の平滑化の効果.
これまで, 体積プロット, 曲面プロット, 等高線プロットにおいて, 平均化やジャンプがどうなるかを調査してきました. 折れ線グラフの場合は, さらに複雑な問題があります. 線に隣接する面が1つ以上ある場合は, 2方向での平均化が必要となります:
- 線が1つまたは複数の面に接続されているかどうかに関係なく, 線に沿って. これは上記と同じ方法で行われ, 同じ“品質”設定セクションが利用可能です.
- 線を挟んで, 隣接する面の間 (複数ある場合). この操作は最初に行われ, 無条件に実行されます.
プロットされる量が線上に連続してはならない場合, もはや平均値は必要なく, 境界 (3D) またはドメイン (2D) ごとに1つのグラフが必要です. このような状況に対処するためには，特別な演算子を使用する必要があります.
線が2Dの境界を表す場合, up()と down()演算子を使って, どのドメインから結果を取得するかを指定することができます. “up”と”down”は, 境界上で法線ベクトルがどのような向きになっているかを指しています. この方法は, 前述の熱流束のグラフで使用されていました. 下図は, x 方向の温度勾配を繰り返したグラフですが, デフォルトの平均結果も含まれています. 境界を越える際に物理的な不連続性がある場合, デフォルトの平均化は明らかに良い選択ではありません.
x 方向の温度勾配, 平均値, 各ドメインの値.
同じ手法を3Dの内部境界のサーフェスプロットで使用することもできますが, この場合は当然ながら一度に片側からしか値をプロットすることができません.
選択した線が 3D のエッジである場合, そのエッジを共有する境界線が任意の数だけ存在するため, 状況はやや複雑になります. この場合, side()演算子を使用して, 個々の境界から値を選択する必要があります. 副演算子は, side(12,shell.mises)のような構文を持ち, 第一引数に境界の番号を指定します. したがって, まず, 結果が必要な境界の番号を把握する必要があります. これを行う1つの簡単な方法は, 組み込み変数 domのサーフェスプロットを作成し, 目的の境界をクリックしてその番号を確認することです.
サーフェスプロットを使用して境界の数を求める
境界番号を確認するもう1つの方法としては, “モデルビルダー”のツリーで境界が選択されているノードに移動し, マウスポインターを境界の上に置くこともできます. 境界番号は, “グラフィックス”ウィンドウの左上に動的に表示されます.
境界選択機能を用いて境界番号を確認.
結果はどこで評価されるのか?
カラープロットや折れ線グラフを作成する場合, プロットされる式は通常, 各要素, 要素境界, または要素エッジの内の数点にて評価されます. 評価点の数は, “品質”セクションの“解像度”設定で制御されます.
“解像度”のオプション.
定義済みの解像度の1つを選択すると, 評価されるポイントの数は, モデルサイズや空間次元などのいくつかの要因に依存します. 高い解像度を使用すると, 通常, 要素内により詳細なプロットが生成されます. しかし, これは評価された式が実際に要素内で高度な滑らかさを持つ場合にのみ有効です. そうでない場合は, 偽の変動を引き起こすだけです.
応力の計算方法
ここまでは, 主に一般的な結果の表示方法についてご説明してきました. 応力の可視化に関しては, もう一つ考慮すべきことがあります. それは, 各要素内で実際に応力がどのように計算されているかということです.
ひずみの合計は変位の微分であるため, 常に滑らかであると言えます. しかし, 応力は, 多くの場合, いくつかの異なる効果の組み合わせの結果です. 応力の計算方法の正確な詳細は, トップレベルの材料モデルと幾何学的な非線形性が有効かどうかによって異なります. 一般的な原理を説明するために, トップレベルの材料が線形弾性体であり, 解析が幾何学的に線形であると仮定してみましょう.
そうすると, 応力テンソル\sigmaは, 次のように計算されます.
ここで, Cは弾性テンソル, \varepsilon_\mathrm{el}は弾性ひずみ, \sigma_\mathrm{add}は余分な応力の寄与です. 弾性ひずみは, 形状関数と節点変位を用いて計算された全ひずみ\varepsilon_\mathrm{tot}と非弾性ひずみ場\varepsilon_\mathrm{inel}との差です.
非弾性ひずみの例としては, 以下のようなものがあります:
- 熱歪み
- 吸湿性膨潤歪み
- 塑性歪み
- クリープ歪み
- 初期ひずみ
余分なストレスの寄与の例として, 以下のようなものがあります:
- 粘弾性応力
- ダンピングによる応力
- 初期応力
したがって, 一般に, 全応力はいくつかの異なる寄与の合計となります. 特に, \varepsilon_\mathrm{tot} と \varepsilon_\mathrm{inel}は同じ桁の大きさであることが多いため, 弾性歪テンソルには2つ以上の大きな数値から別の数値を引き算して得られる小さな数値が含まれす.
では, なぜこれが問題になるのでしょうか? 応力と歪みの寄与が異なる場合, 要素上で同じタイプの補間によって表現されないと, 平均的な意味では結果が一致していても, 大きな局所的な変動が生じることがあります.
特別な場合: 熱膨張
ここで少し回り道をして, このようなことが起こり得る一般的なケースを見ていきましょう. 熱応力連成解析では, “波状”の応力パターンが発生することが, 有限要素法の研究者たちによって早くから観察されていました. 熱応力連成解析では, 変位と温度の両方に 2次関数を使用するのが最も一般的です. 熱歪みは温度に比例するため, \varepsilon_\mathrm{inel}は各要素内で2次関数的に変化することになります. しかし, 全歪みは変位形状関数の微分であるため, 線形分布となります. さて, 弾性歪みは, 1次関数と2次関数の差として計算されます. その結果, 各要素内に奇妙な弾性歪み (ひいては応力) パターンが存在する可能性があります. このため, 温度を使用して固体の熱膨張を駆動する場合は, 線形形状関数を使用して熱問題を求解することが推奨されることがあります. COMSOL Multiphysics では, この問題は“熱膨張”マルチフィジックスカップリング (および“吸湿膨張”や“インターカレーション歪み”などの類似機能) で内部的に処理されます. 非弾性歪み場は, 変位場から計算された歪みと一致する多項式次数に再補間されます. これにより, 局所的な応力解におけるこのようなアーチファクトを低減することができます.
点ごとの状態
場を全く扱わず, 積分点 (ガウス点) に点状に格納された局所的な状態を扱う, 別のタイプの非弾性歪みがあります. ガウス点についての詳細は, ブログ““数値積分とガウス点への入門””をご参照ください. 可塑性やクリープなど, ほとんどの非線形材料モデルは, このような仕組みになっています. この場合, 要素内の任意の位置での応力の評価はより複雑になります.
solid.sxなどの応力成分を指定すると, 実際には要素内の各位置で”突発的に”計算されます. 非線形材料法則は, 最も近いガウス点から格納された値を利用して, そこで評価されます. 要素全体の非弾性歪みのばらつきが大きい場合, 大きな誤差が生じることがあります. また, ガウス点で評価できるにもかかわらず, 材料則を評価できない場合もあります.
より良い方法は, ガウス点の値に基づいて滑らかな近似値を作成することです. gpeval()演算子によってこれは可能になります. 例えば, solid.sxではなく, gpeval(4,solid.sx)を指定すると, 要素上に滑らかな応力がプロットされます. この場合, 非線形応力–歪み関係がすでに満たされているガウス点でのみ評価されます. そして, ガウス点の値の最小二乗フィットを用いて, 滑らかな場を定義します.
このバージョンでは, gpeval()演算子は, 最初の引数として適切な積分順序を入力する必要があります. 幸い, ほとんどの場合, これを気にする必要はないでしょう. フィジックスインターフェースに, solid.sGpxやcode>solid.misesGpなどの定義済み変数が多数用意されています. もし, 式に適した組み込み変数が見つからない場合は, フィジックスインターフェース固有の演算子, 例えば, solid.gpeval(expr)を使用することができます. 正しい積分順序が, フィジックスインターフェースに特有の演算子に埋め込まれています.
以下の例では, (myDOFというユーザー定義の従属変数によって) ガウス点で格納された場について, さまざまな評価オプションが調査されています. モデルは1つの2次元要素からなり, xと yの座標はともに0から1までです. 場は (Y+1)/(2*X+1) という式で記述されます. 3×3のガウス点スキームが使用されるため, 9つの独立した数値が格納されます. これらの値は, 各ガウス点で評価された元の場の値と全く同じです.
ガウス点データのさまざまな評価方法
同じデータを高さプロットで表現したも.
左上のプロットでは, 元の関数が可視化されています. その下に, 格納された状態のmyDOFのサーフェスプロットの挙動が示されています. 同じプロットに, ガウス点の位置と格納された値も示されています. ガウス点の状態は, 最も近いガウス点の値で表されるというデフォルトの動作が, はっきりと目に見えます. この場合, より良い方法は, gpeval()演算子を用いて滑らかな場を評価することです. これは, 下段の右端の2つのプロットに示されています. これらの2つのプロットの違いは, 平滑化場の多項式の次数です. デフォルトでは双一次場が使用されます. この例では, 2次場へのフィッティングがより良い選択肢です. 中上と右上のプロットでは, 平滑場と関数の正確な値との差を示しています.
最後に, よくある2つの誤解について説明していきます:
solid.sxのような式をガウス点で評価すると, 確かにガウス点での計算値が得られます. ただし,solid.sGpxやgpeval(4,solid.sx)のような式をガウス点で評価しても, その点に格納された結果は得られません. なぜなら, これらの式が最小二乗法によるフィッティングで得られた滑らかな場を与えるからです. フィットされた多項式は必ずしもそのデータ点を通過するわけではありません.- >
gpeval(4,solid.sx)のような式を境界で評価すると, 境界上のガウス点の値が補間多項式の設定に使用されます. これは, 固体力学モデルの境界でサーフェスプロットを行う場合, おそらく意図したものではありません. 幸いなことに,solid.sGpxやsolid.gpeval(solid.sx)といった組み込みの量は, 期待通りの動作 (ドメイン補間)をしてくれます.
最大値評価の例.
これらのアプローチで得られる値は, 必ずしも同じではないことは明らかです. いくつかのヒントをご紹介します:
- プロットやグラフでピーク値を強調表示したい場合は, 対応する“マーカー”サブノードを使用する. 値が一定でないとプロットを見ている人は混乱してしまいます.
- 最悪の場合の値を取得したい場合は, 平均化を抑制する評価を使用します.
- 平滑化が使用されているか, 応力計算に点ごとの状態が含まれる限り, 予想される最大値が表面にあっても, 同じ値が得られないことがあります.
- 構造力学の問題では, 最大応力が境界で発生することがよくあります. もし, 境界に荷重がかかっていない場合, その境界で正確な結果を得るための良い方法があります. 境界の上に, 固体と同じ材料データを持つ非常に薄いメンブレンを追加することです. このメンブレンは, 一種の仮想歪みゲージとして機能します.
以下は, サーフェスおよびボリュームプロットでマーカーを使用し, 応力成分の最小値を評価するためのさまざまなオプションの影響を示すプロットです. solid.syと solid.sGpyの両方が, 平均化あり, なしでプロットされています. 予想通り，平滑化なしで計算された各値は, 平滑化ありの対応する値よりも高くなっています. また，平均化を行わない場合，サーフェス評価とボリューム評価は同じ結果になることにもご注目ください. これは, ピーク値が表面にあるため, どちらの場合も同じピーク位置がプローブされるため, 当然のことです.
さまざまな評価オプションの効果.
どの応力結果をプロットすればよいのか?
プロットする結果を選択すると, 以下のようなメニューが表示されます:
応力結果の選択.
利用できる正確なオプションは, フィジックスインターフェースと, どの材料モデルや他のオプションを使用しているかによって異なります. オプションの数に圧倒されるかもしれませんが, 他よりも重要な量がいくつかあります.
等方性材料では, フォン・ミーゼス応力やトレスカ相当応力などのスカラー応力を使用するのが一般的です. フォン・ミーゼス応力は, 簡単に評価でき, 感度などの計算が容易であるため, 人気があります. しかし，安全な方を取りたい場合, トレスカ応力の方が適しています. トレスカ相当応力の値は, 特定の応力状態のフォン・ミーゼス応力の値よりも 0 ～ 17% 大きくなります. 圧力容器などの工学分野では, 常にトレスカ応力が使用され, 応力強度という名前で呼ばれることもあります.
フォン・ミーゼス相当応力やトレスカ相当応力は, 応力状態を概観するのに非常に有効ですが, 破壊のリスクを評価したり, 特定のクラスの材料について”最大荷重点”を見つけるためにのみ使用すべきです. これらの応力測定は, もともと金属の降伏を予測するために設計されたもので, 平均応力には敏感ではありません. 大気中にある金属は, マリアナ海溝の底にある金属と同じくらい降伏に近いのですから! 一方, 土やコンクリートなどの材料の破壊は, 平均応力に強く依存します. 圧縮された状態が安定します.
材料が金属でない場合, 赤い点が臨界点となるようなプロットはどのように作成すればよいのでしょうか? 答えには2部分あります. まず, 一般的には, 破壊を記述する材料特性を導入しなければ, それを行うことはできません. これらの値を得ることができたら, 構造力学インターフェースの“安全”ノードを使用します. この機能では, さまざまな種類の材料について, 破壊までのマージンを計算することができます. これらの基準は応力に基づくことが多いため, 応力評価について上記で述べたことは, ここでも適用されます.
“安全”機能で見られる失敗モデル.
土質力学問題のフォン・ミーゼス応力 (上) および安全係数 (下) の解析. 以下の解析では, Drucker–Prager 基準を使用しています.
場合によっては, 個々の応力成分を調べたい場合があります. これは, 材料が異方性であったり, 応力分布をよりよく理解したい場合です. 後者の場合, 主応力もプロットすると便利なことが多いです.
個々の応力成分を扱う場合, グローバル座標系以外の方向で評価したいことがよくあります. “ローカル座標系”という名前の結果には, いくつかの種類があります. これは, 材料ノードの座標系選択 (例えば, “線形弾性材料”)によって方向が決定されることを意味します.
材料モデルでの座標系の選択.
等方性材料の場合, 座標系選択の効果は結果の局所的な方向を定義することのみです. 異方性材料の場合, この選択の主な目的は, 材料データの参照方向を提供することです. 副次的な効果として, 応力と歪みの結果を局所的な方向で得ることができます.
COMSOL Multiphysics では, 3 種類の応力テンソルがあります. Cauchy, 第1種Piola–Kirchhoff, 第2種Piola–Kirchhoff です. このトピックの詳細については, ブログ"なぜこのように多くの応力と歪みが必要なのか?"および Multiphysics Cyclopedia の応力運動方程式のページをご覧ください.
最後に
COMSOL Multiphysics には, 結果の評価と表示を微調整するための多数のオプションが用意されています. シミュレーションを最大限に活用するためには, これらのオプションに精通することが重要であり, 何が最適かは解析内容や解析の目的によって異なります.
ここで解説されていないこととして, コーナーやその他の特異点で発生する高い応力に対処する方法があります. 特異点の影響については, ブログ “有限要素モデルにおける特異点: レッドスポットへの対処”で解説しています. 実用上非常に興味深いこの話題は, また登場することになりそうです.
次のステップ
当ソフトウェアの応力解析機能について詳しく知りたいですか? , 以下のボタンをクリックし, COMSOL へお問合せください!
コメント (0)