固体力学

目次固体力学境界条件
はじめに付録
概要的な例と解析タイプ用語集
方程式参考文献

目次

はじめに

荷重および制約下における固体の解析および挙動は，力学において基本的な重要性を持つ．固体力学では三次元における固体の力学を扱い，構造力学にはシェル構造や梁等の幅広い物体が含まれる．

このチュートリアルでは，偏微分方程式を使った固体力学のモデル化を紹介する．固体力学解析を行うのに関係する方程式および境界条件を導出し，説明する．

偏微分方程式(PDE)で固体力学をモデル化することは，モデル化の唯一の方法ではない．常微分方程式(ODE)を設定するという別の方法もある．このアプローチはWolfram System Modelerで使われている．簡単に言うと，System Modelerのアプローチは大規模な固体系の相互作用により適しており，偏微分方程式のアプローチは特定の物体の細かい解析により適している．両方のアプローチを組み合せて使うのがよい場合もある．

ここでは，まず最初のセクションで単一物体（本棚の棚受け）を使ってさまざまな固体力学解析の種類および利用できる機能を紹介する．その後で根底にある考えや概念というより理論的な説明をする．理論的背景は，さまざまな解析タイプについての知識があった方が格段に理解しやすい．さらに利用可能な境界条件について解説する．

固体力学では，通常三次元の固体の物体を考える．二次元の簡約化を扱う特殊な場合もある．しかしこの簡約化には，三次元の理解が正しければ避けることのできる危険性もある．このような訳で，初めの例は三次元のものである．終盤で平面（二次元）応力およびひずみ，その限界について述べる．

固体力学解析の目標は，荷重下における物体の変形を求めることである．次のステップで，その変形した物体内部のひずみや応力を求める．これらの物理量の解析および解釈は，検討中の物体の工学的設計の質を高めるのに便利である．例えば，物体で構造的に弱い部分を見付けて改善することができる．

固体力学解析の過程は通常段階を追って行われる．まず，解析する物体の幾何モデルを作成する必要がある．幾何モデルは通常コンピュータ支援設計(CAD)で作成される．CADモデルはインポートすることも，製品内で作成することもできる．幾何モデルをインポートするために，DXF，STL，STEP等の一般的なファイル形式がサポートされている．幾何モデルはImportを使ってインポートすることができる．またはOpenCascadeLink等を使って，製品内で幾何モデルを作成することもできる．一旦幾何モデルが利用可能になると，どのような解析を行いたいかを考える必要がある．サポートされている解析の種類には，静力学的解析，時間依存解析，固有モード解析，周波数応答解析がある．次のステップでは，境界条件と制約条件を設定する．さらに使用される材料でPDEモデルが指定される．PDEモデルが完全に指定されたら，続いて有限要素解析で，調査中の物体の所望の数量が計算される．それからこれらの数量は，可視化されるか導出された数量が計算されるかして，後処理される．このノートブックは，CADモデル生成以外の必要なステップを示す．

モデリング過程はそれ自体，結果としてNDSolve，ParametricNDSolve，NDEigensystemで解くことができる偏微分方程式(PDE)系になる．

固体力学PDEモデルの正確さおよび効率は「固体力学モデル検証テスト」という別のノートブックで検証されている．

このノートブックで示されているシミュレーション結果のアニメーションの多くはRasterizeを呼び出して生成される．これは，このノートブックに必要なディスクスペースを減らすためである．この呼出しを使った場合の欠点は，使わなかったときほどくっきり見えないということである．高品質のグラフィックスが得たい場合は，Rasterizeの呼出しを削除するかコメントアウトするかするとよい．

高忠実度の可視化を得るためには，ラスタライズ処理をコメントアウトする．

このチュートリアルで使われる記号と対応する単位は用語集のセクションにまとめてある．

概要的な例と解析タイプ

固体力学における有限要素法の使い方を説明するために，簡単な例を示し，設定の概要，さまざまな解析タイプ，可能な後処理のステップを解説する．

有限要素パッケージをロードし，$HistoryLengthを0に設定する：

最初の例では，簡単な固体力学PDEモデルを設定するためのワークフローを紹介する．その後のセクションでは，より高度で完全なモデルを使う．ここでは，左端が固定されて下向きの力がかかっている長い梁の以下の設定を考える．

1.gif

固体力学モデルの作成は，常に同じステップで構成される：

幾何モデルを設定する：
変数を設定する：

従属変数 ，， はそれぞれ ，， 方向における変位を表す．

ブロックはチタン製である：
境界条件を設定する：
固体力学PDEの要素を設定する：

出力は大きく読みにくいので短縮されているが，このような訳でPDEModelの要素が便利なのである．

PDEを解く：

このような多くの場合，メッシュの生成は不要であり，幾何モデルをソルバに渡すだけで十分である．

物体の変形の弧長表現を，もとの物体の輪郭とともに可視化する：

以上がワークフローの概要である．次のセクションでは，もっと複雑な幾何モデルについてステップを詳しく説明する．

形状

以下の利用可能な固体力学解析の概要では，定義済みの境界要素メッシュがロードされる．この境界要素メッシュがどのように生成されたかについての情報は，この幾何モデルの生成方法を説明した「本棚の棚受け」チュートリアルをご参照いただきたい．

本棚の棚受けについての定義済み境界メッシュをインポートする：
定義済み境界メッシュを可視化する：

インポートされた境界メッシュに欠陥がないかどうか確かめた方がよい場合がある．欠陥が見つかった場合，欠陥が書かれた凡例が表示される．この場合，欠陥はない．

境界メッシュに欠陥があるかどうかチェックする：

一つ留意しなければならいのは，幾何モデルで使われたスケールである．境界メッシュの長さの単位が例えばメートルだとすると，材料パラメータも同じ単位で指定しなければならない．この例の場合，境界メッシュの単位はメートルなので，棚受けの長さは0.2メートルである．

幾何モデルの詳細は，有限要素メッシュがそのモデルを表すのに必要な要素数に影響する．メッシュの要素数は，特定の問題を解くのに必要なCPU時間とメモリ量に直接影響を及ぼす．要素数を減らすためには，ねじ穴等の幾何的詳細は通常その付近にしか影響しない[11, c.1]ことを考慮するとよい．サンブナンの原理によると，特性長 の部分は，約 の距離の範囲に影響を及ぼす．もし，特定の部分から十分離れている領域の物体の力学的性能に関心があるなら，その部分は幾何モデルから切り離され，計算が節約できる．また，振動解析では，幾何学的断面積の約10%より小さい幾何学的部分は通常無視することができる．

この例では，棚受けにはねじ山がない．これは境界条件を指定するときに扱う．

一般に，できるだけ簡単な幾何モデルから始めて，全体的な解に影響を与えるのか，またどのように与えるのかを見るために詳細を加えていく．

3D幾何モデルの生成およびインポートについての詳細は，「OpenCascadeLinkを使う」チュートリアルに書かれている．

材料パラメータ

次のステップは材料パラメータの割当てである．一般に固体力学モデルのパラメータはすべて，必要なパラメータ値を含むAssociation で集められる．

チタンについて特定の材料特性を設定する：

材料パラメータの設定方法では，Entityから材料を選んだ方が簡単である．Entityを得るためには，WolframAlphaの自由形式入力を使うと便利である．

WolframAlphaを使って材料パラメータを設定する：

Quantityオブジェクトとして値を与えることのできる材料パラメータを使うことができる．材料が利用可能ではない，あるいは異なる単位が使われている場合は，YoungModulusやPoissonRatio等のパラメータを直接指定して材料特性を加えることができる．必要とされる正確な特性名は，SolidMechanicsPDEComponentの参照ページで見付けることができる．

デフォルトの材料モデルは線形等方性弾性材料である．経験則で，線形弾性材料モデルは最大伸張が5%に達するまで適用できる[8, p. 159]．

いくつかの材料でできている幾何モデルも使うことができる．これの例は「複数の材料」のセクションで取り上げる．

単位

幾何モデルの単位が材料単位と異なる場合は，材料単位をスケールすることができる．

材料を設定して，単位をスケールする：

内部的にはすべての材料データの単位は"SIBase"単位に変換される．その結果長さのデフォルトの単位は"Meters"になる．幾何モデルの単位もメートルならば，何も変更されない．幾何モデルの単位がメートルでない場合は，PDEと材料特性は幾何モデルの単位にスケールされるか，幾何モデルが"Meters"にスケールされるかする必要がある．PDEおよび材料パラメータの単位をスケールするためには，パラメータ"ScaleUnits"を与えることができる．明示的に述べられていない場合は，このノートブックの例ではデフォルトの"SIBase"単位を使う．

境界条件

固体力学は荷重と制約条件下における物体の変形について調べるので，境界条件は固体力学モデルの重要な要素である．境界荷重は物体の表面にかかる力または圧力である．荷重が物体に適用されるためには，壁にねじで止められている等，その物体も何らかの方法で制約されていなければならない．そうでないと物体は荷重に対して抵抗しないからである．荷重と制約条件は，境界条件を指定することによって設定する．さまざまな境界条件を使うことができる．これについては境界条件のセクションで詳しく述べる．この概要では，壁とねじによる境界面荷重と制約条件で十分である．

このセクションでは，境界条件が適用される位置を決めることである．境界条件を適用する位置は，実行される解析タイプにかかわらず，同じままである．しかし，境界条件の値の厳密な指定は解析タイプによって異なる．これはそれぞれのセクションで説明する．

境界条件を適用する位置を求める方法として，幾何モデルのセクションで定義した幾何モデルの輪郭と境界条件を一緒に可視化するというものがある．

境界荷重に対する面を確立する：
境界荷重が適用される面を可視化する：
面荷重値の境界条件を設定する：

次のステップでは，棚受けが壁に取り付けられているモデルの制約条件を指定する．棚受けを壁に固定するためには2つの制約条件が必要となる．一つは棚受けの背面が壁を突き抜けないことで，もう一つはねじで棚受けを壁に固定するということである．

ねじを付けるという制約条件から始めると，先も述べたように幾何モデルにはねじ山がない．これは，ねじ穴を作る表面のすべての方向における動きを完全に制約することによって解決できる．つまり，ねじと棚受けが接触する点では全く動きがないということである．

棚受けは壁に突き刺さることがないため，棚受けの背面の動きを負のx方向に制約する．2つのねじが棚受けを壁に押し付け，棚受けは壁に食い込めないため，正の （壁から外側）方向の動きも制限することが妥当である．これは一般的な簡約化である．この簡約化が正当化されない場合は，接触問題が起こる．

背面を構築する：
背面をx0で可視化する：
部分的な変位制約条件を設定する：

次にねじの制約条件を記述する．

ねじ穴と溝に対する面を構築する：

ねじ穴と溝をカバーする十分大きい2つの円筒領域を使う．

ねじが動きを妨げるところを可視化する：
境界条件としてねじを設定する：

もっと複雑な幾何モデルの場合は，境界条件の述語を指定するのに別の方法を使った方が適しているかもしれない．これは付録の境界条件の述語のセクションで説明する．

メッシュ生成

有限要素解析を行うためには，幾何モデルの境界メッシュ表現をメッシュに離散化する必要がある．

メッシュを生成する：

メッシュの生成過程についての詳細は，要素メッシュの生成チュートリアルに記載してある．別の方法としてメッシュをインポートするというものがある．一般的なメッシュファイル形式にはFEMAddOnsでインポートできるものがある．

定常解析

固体力学解析は適用された力および制約条件の結果としての物体の変位を求める．変位の従属変数は ，， と言われ，それぞれ ，， 軸方向の変位を表す．

変数を設定する：
固体力学PDE成分を設定する：

デフォルトの設定では，わずかな変形を想定した線形等方性弾性材料のモデルが生成される．このモデルには物体の自重は含まれない．これはパラメータ"BodyLoadValue"で簡単に加えることができる．これを加えることについてはセクション物体の荷重で説明する．

負の 方向の面荷重を設定する：
背面の制約条件をローラー制約条件として設定する：
ねじの制約条件を設定する：
PDEを解く：

結果は3つのInterpolatingFunctionオブジェクトのリストである．それぞれのInterpolatingFunctionはそれぞれの従属変数についての変位を与える．3つのInterpolatingFunctionオブジェクトで変位ベクトルが構成される．

解く過程とそのオプションについての詳細は，有限要素のためのNDSolveオプションチュートリアルに記載してある．

後処理

固体力学PDEモデルの主要な解は，作用する力によって生じる変位である．この変位は，荷重と制約条件下でどのように物体が変化するかを可視化するために使うことができる．しかし物体内の応力とひずみも面白い．物体のひずみは変位によって回復され，応力はひずみによって回復される．したがって，ひずみと応力はともに二次的な未知数である．回復したひずみと応力は，相当ひずみあるいはミーゼス応力等の総括概念でさらに組み合せることができる．後処理という言葉は，変位が計算された後に行われるすべての計算および可視化に対する総称である．

変形

変位にはよく ，， が使われ，それぞれ ，， 方向において計算された変位である．変位ベクトル が使われることもある．

変形

荷重下における物体の変形は，もとの物体の座標に加えられた計算された変位に過ぎない．物体が荷重下にあるとき，もとはにあった物体の点 がに移動する．

変形した物体のメッシュを作成する：
変形の誇張表現と，棚受けのもとの位置の輪郭線を可視化する：

示した変形は，幾何モデルの伸張の最小サイズおよび最大変位でスケールされている．これはElementMeshDeformationに"ScalingFactor"を指定することで，無効にしたり変更したりすることができる．詳細は変形の可視化をご覧いただきたい．

変形と棚受けのもとの位置の輪郭線を固定倍率10で可視化する：

ElementMeshDeformationはデフォルトでは，変形が非零である限りどれほど小さかろうと可視化する．しかし可視化は真の変形を描写しない．真の変形は"ScalingFactor"->Noneと設定することで表示できる．

倍率を指定せずに変形を可視化する：

ほとんどの場合，変形は小さいので"ScalingFactor"->Noneでは興味深いプロットにはならない．

自動倍率の計算をしたときの結果として，物体の荷重が増加しても，値が再びスケールされてしまうため変形プロットは同じままということがある．変形プロットの値は物体がどのくらい変形するかを見るものである．これによってモデルが想定通りに動作しているかどうかが分かる．

計算される自動倍率の値を知っておくと役に立つことがある．

使用されたAutomatic倍率を計算する：
全変位

全変位を調べると，荷重下で制約のある物体の変位の概要を捉えることができる．全変位は以下で与えられる：

ここで ，， はそれぞれ ，， 方向で計算された変位である．

全変位を可視化する：

ひずみ

ひずみ

ひずみとは物体内の変形の量を示す数量[2]であり，割合であって単位はない．ひずみは変位量と混同してはならない．ひずみは指定された変位の勾配によって変位と関連付けられる．ひずみの概念はひずみについての理論セクションで詳しく説明する．

変位によるひずみを計算する：

SolidMechanicsStrain関数は，独立した各方向についての垂直ひずみとせん断ひずみを計算する．この関数は，さまざまな方向のひずみ成分を表す式を含むSymmetrizedArrayを返す．SymmetrizedArrayが存在しコンパクトな場合，これは対称性を強制する．返されたひずみは以下の順番である：

ここで垂直ひずみは ，せん断ひずみは で表されている．ソフトウェアはデフォルトでは工学ひずみ， を使う．パラメータで"EngineeringStrain"->False と設定すると，工学ひずみの使用を切り替えることができる．これは真ひずみを利用する，線形材料則と非線形材料則を比較する場合に便利である．詳細はSolidMechanicsStrainの関数ページに記載されている．

大雑把に言うと，ひずみは変位の勾配である．ひずみプロットは，大きい変形があるところではなく，変位勾配が大きいところを示す．

と の垂直ひずみを可視化する：

この例では，負の 方向の変位が正の 方向のものよりずっと大きいにもかかわらず， 方向のひずみの絶対値が 方向の絶対値より小さいことに注目する．

ある点におけるひずみ成分をすべて評価する：
すべてのひずみの最小値と最大値を求める：
垂直ひずみ および を同じカラースケールで可視化する：
せん断ひずみ を可視化する：

ひずみが，ある特定の値より大きくなっている領域を見付けたいことがある．

ひずみが >0となっている領域を色付きのセルで表示する：
相当ひずみ

すべてのひずみ成分を調べるのは厄介である．相当ひずみはすべてのひずみ成分をスカラーにまとめる簡約化であるため，理解しやすい．しかし相当ひずみはひずみの完全像は含まない．

相当ひずみは以下を計算する：

相当ひずみの実装の中には（ はポアソン比）の因数を使うものがある．相当ひずみを，相当ひずみの応力版であるミーゼス応力と互換にするために，この因数を含まないことにした．

相当ひずみを計算する：
立体上で相当ひずみの等高線を表示する：

変形した物体上にひずみをプロットした方が便利な場合がある．

相当ひずみの関数を抽出する：
変形したメッシュ上に相当ひずみの補間関数を作成する：
変形した物体上に相当ひずみを可視化する：
主ひずみ値と方向ベクトル

主ひずみ値は主となるひずみ値である．主ひずみ値はひずみ行列の計算された固有値であり，せん断ひずみ が0である局所座標系に対応する．主ひずみ値は局所座標の であり，荷重方向に独立のひずみを与える．

主ひずみの値とベクトルを計算する：

主ひずみ値は大きさによって並べられる．

座標における第1，第2，第3の主固有値を計算する：

それぞれの評価座標には3つの主ひずみ値がある．そのいずれもベクトルとして解釈し可視化することができる．

第1，第2，第3の主ひずみ値をベクトルで示す：
第1，第2，第3の主ひずみ方向ベクトルを局所座標軸として表示する：

セクション境界荷重：張力には，PrincipalEigensystemを使って非軸の荷重がかかった円筒における応力を求める方法を示す例が含まれている．

応力

応力

応力とは物体の内部の力の分布を記述する数量[2]であり，またはで測定される．

ひずみから応力を計算する：

Stress関数は独立した各方向についての垂直応力とせん断応力を計算する．この関数はその方向の応力成分を表す式を含むSymmetrizedArrayを返す．返された応力は以下の順番である：

ここで垂直応力は ，せん断応力は で表される．最初の指標は力の方向を示し，2番目の指標は，力が作用する表面の方向を指定する．

垂直応力 と を等高線で表示する：
せん断応力 を可視化する：
ある点における応力成分をすべて評価する：
ミーゼス応力

すべての応力変数を調べるのは厄介である．すべての応力成分を一つの式にまとめる簡約化がミーゼス応力である．ミーゼス応力はスカラーであるため理解しやすい．しかし，ミーゼス応力には物体内に存在する応力の完全像は含まれていない．応力に対するミーゼス応力はひずみに対する相当ひずみと同じである．

ミーゼス応力は以下を計算する：

ミーゼス応力を計算する：
ミーゼス関数を抽出する：
ミーゼス応力の最小値と最大値を求める：
ミーゼス応力の等高線プロットを示す：
上のプロットでミーゼス応力の極値がどこかを求めて表示する：

変形した物体上に応力をプロットした方が便利な場合がある．

変形したメッシュ上にミーゼス応力の補間関数を作成する：
変形した物体上にミーゼス応力の等高線を表示する：
ミーゼス応力が となる領域を色付きのセルで可視化する：
主応力値と方向ベクトル

主応力値は主となる応力値のことである．主応力値は応力行列の計算された固有値であり，せん断応力 が0となる対角応力値 に対応する．主応力値は荷重方向に関わりない応力を与える．

主応力値とベクトルを得る：

座標で評価すると，は第1，第2，第3の主応力値の3つの値のリストを返す．

座標における第1，第2，第3の主応力値を計算する：
第1，第2，第3の主応力値を可視化する：
座標における主応力ベクトルを計算する：
第1，第2，第3の主応力方向ベクトルを可視化する：

固有ベクトルにスカラーを掛けても固有ベクトルのままなので，ここでは凡例は意味をなさない．

セクション境界荷重：張力には，PrincipalEigensystemを使って非軸の荷重がかかった円筒における応力を求める方法を示す例が含まれている．

PrincipalEigensystemは固有値を最も大きいものから小さいものへと並べ替えるので，主固有値の記号計算には使うことができない．記号式の並べ替えは不可能だからである．

安全率

安全率とは，結果となる最大応力が応力の限界を下回る割合がどの程度大きいかを示す値である．安全率はFOSと略されることもあり，複数の定義がある．扱う材料が延性のある脆くない金属の場合，ミーゼスの破壊基準が妥当である．延性のある材料では，塑性変形の開始で破壊の始まりとみなされる（脆い材料では，破砕で破壊とみなされる）．詳細は，破壊理論のセクションでご覧いただきたい．

応力を荷重方向に独立な応力に変換するために，応力のかかった物体の主固有値を計算する．これは後でミーゼス応力を計算するのに使うことができる．

主応力固有系を計算する：

せん断応力が0に設定されていれば，相当ミーゼス応力はミーゼス関数で計算することができる．

相当ミーゼス応力を計算する：
相当ミーゼス応力関数を抽出する：
相当ミーゼス応力の最大値を求める：

延性のある材料では，応力が降伏強度に達すると塑性変形が始まる．脆い材料では応力が極限強度に達するときである．

材料の降伏強度を求める：
材料が降伏するかどうかをチェックする：
安全率を計算する：

安全率が1未満であると問題である．この場合，非線形応力・ひずみ関係を仮定した．非線形応力・ひずみ関係が使われると，最大相当応力は異なる可能性がある．妥当性の範囲外で線形材料モデルが使われる場合は，計算された応力は通常実際の応力よりも高い．それでも降伏応力より高い応力を求めることは非常に注意が必要である．

安全率が10未満である領域を色付きのセルで可視化する：

反力

制約と同時に荷重下にある物体はどれも，荷重の合計と釣りある反力を持つ．この反力は物体の制約された部分にある．

固体力学PDEモデルから方程式の離散系を作成し，反力を計算する：

関数PDEModels`ReactionForceは，NDSolveの通常の呼出しのように変位を返し， 反力であるInterpolatingFunctionの第2の集合を返す．

各方向における合計反力を計算する：

反力の合計は境界荷重の値に一致する．

境界荷重を調べる：
反力をベクトルプロットとして可視化する：

示された反力は，棚受けに作用する下向きの力（表示されていない）に対抗する．

現在のバージョンでは，線形定常固体力学モデルについてのみ反力が計算できる．

時間依存解析

PDEモデルの成分の中に時間依存の挙動をするものがある場合は，時間依存解析を使う必要がある．一般的な時間依存成分は時間依存境界条件である．

時間依存の固体力学PDEを設定する：
負のz方向における時間依存面荷重を設定する：

物体に動的に荷重を掛ける際の重要な数値的側面は，荷重は即時に適用されないということである．場合によっては，そのようなモデルは反物理学的であり結果的にシミュレーションの時間が長くなることがある．

ねじの制約条件を設定する：
背面の制約条件を設定する：

時間依存の固体力学方程式は二次時間依存なので，初期条件は初期変位と初期速度場から構成される．初期条件は物体全体に影響を及ぼす．例えば，非零の初期変位を指定するということは，物体全体が初期時間に指定された量だけ変位するということである．

初期条件と初期速度条件をゼロで設定する：
時間依存PDEを解く：
さまざまな時間における変形を示すフレームを作成する：
アニメーションを可視化する：

時間経過に伴う物体内の特定の点の動きを監視するために，クエリポイントを設定する．クエリポイントは時間を通して追跡することができる．

幾何モデルのメッシュの中に赤いクエリポイントを可視化する：
クエリポイントにおける時間経過に伴う垂直変位 w を可視化する：

系は減衰せず，変位の振幅は周期を通して同じままである．

ミーゼス応力を計算する：
ミーゼス応力の時間依存の可視化を作成する：
アニメーションを可視化する：

時間依存解析は時間がかかりメモリも多量に必要とする．このため，解析が本当に過渡的である必要があるかどうかを推定できることが重要である．せん断波が特性長 で弾性立体を通過するのにかかる時間はおおよそ以下になる[11, c.1]．

ここで は質量密度， は剛性率である．応力値は約で停留値に減衰する．物体の回転における加速で引き起こされる応力についても同様に考えられる[11, c.1]．

固有モード解析

固有モード解析は物体の共振周波数を計算する．この周波数は固有振動数とも呼ばれる．これらの固有振動において，調査する物体は固有モードという異なる形状に変形する．固有モード解析の結果は周波数とそれに対応するモードのリストである．

固有振動数と固有モードを計算する理由には以下をはじめとしていくつかある：

線形弾性材料に対する固体力学の一般化運動方程式は以下で与えられる：

ここで は質量行列， は減衰行列， は剛性行列， は荷重ベクトルである． は従属変数 ，， の変位ベクトルである． は変位ベクトルの導関数，速度ベクトルと は変位の二次導関数と加速度である．振動解析では減衰 は通常無視されるため[7]以下のようになる：

固有振動数解析では，制約条件は考慮されるが荷重は考慮されないため である．

固有振動数解析では，以下の形式の方程式の解が関心の対象となる．

ここで は と を提供する固体力学PDE成分演算子である． は固有値 である．

固有モード演算子 を設定する：

他の解析タイプと対照的に，固有モード解析はパラメータとして指定する必要がある．他の場合，SolidMechanicsPDEComponentは変数の指定から解析タイプを推定することができる．固有モード解析は時間依存の設定と非常に似ており，変数の入力だけからでは推測されない．

まず制約を受ける固有モード解析を見てみる．つまり，物体は境界条件によって制約を受けるのである．制約のある固有モード解析に関連のある条件は，従属変数が0に設定されるDirichletConditionになる境界条件だけである．境界力および非零のSolidDisplacementConditionは関係ない．

ねじの穴と背面が制約された固有値と固有ベクトルを計算する：

固有値 は以下によって角周波数 と関連付けられる：

角周波数は以下によってで測定される周波数 と関連付けられる：

固有振動数は固有値から計算することができる：

固有振動数を計算する．
そのその周波数で変形モードを可視化する．

変形した形状における変形の量は実際の変形ではない．表示される変形の大きさは任意である．一定の時間において固有モードはまだ同じ固有モードである．固有モード解析はそれぞれの周波数における固有振動数とモードの質的形状を計算する．

物体に制約がないという特別な状況も起こる．そのとき最初のモードは零であり剛体モードと呼ばれる．剛体モードは変形なしで移動したり回転したりする物体を表す．固有値と固有モードは物体が取り得る変形の基本的な形状を捉える．固有モードがないということは，物体が平行移動または回転できたということである．物体が移動したり回転したりできるということは新しいことではなく，このようなモードは興味の対象にはならない．三次元では，各方向の移動について3つ，各軸を中心とする回転について3つの合計6つの剛体モードがある．

最初の12個の固有値と制約条件のない本棚の関数について解く：
最初の固有振動数を調べる：

最初の6個の固有値は比較的0に近く剛体モードを表す．また剛体モードでは小さい虚部を含む固有値を持つことができる．

変形剛体モードを可視化する：

これらは剛体モードである．平行移動あるいは回転される剛体の固有値はゼロである．このすべてで対応する ，， 軸方向で可能な3つの平行移動と同じ軸を中心とする3つの回転がある．固有ソルバによって見付けられた固有モードは純粋な剛体モードではないが，基本的な剛体運動の重ね合せである．

モーダル解析で剛体モードがあることが分かったら，その物体は十分制約されていない．物体が十分に制約されていることを検証するには固有解析がよい．

パラメトリック解析

PDEモデルのパラメータを変化させ，さまざまなパラメータについて同じPDEを繰り返し解きたい場合があるだろう．これを行う便利な方法にパラメトリック解析がある．本棚の例では，いろいろな面荷重の下でどのように本棚が変形するかを見るのは面白いかもしれない．この場合境界荷重の力がパラメータにされる．シミュレーションはパラメトリック解析ではない場合と全く同じ方法で設定され，関数のParametricNDSolveファミリーを使ってモデルのパラメータ名を指定するだけである．

パラメータの荷重を設定する：
パラメトリック関数を作成する：
さまざまな荷重についてパラメトリック関数を評価する：
総変位を可視化する：

結果は質的には同じであるが，総変位の大きさは適用される荷重に応じて線形に変化する．

荷重変位プロット

パラメトリック解析は荷重・変位プロットの作成に使うことができる．特定のクエリポイントについて，荷重が増加するときの変位が追跡される．

さまざまな荷重に対する変位を計算する：
特定の座標における荷重と総変位を可視化する：

この場合，力と変位の間に線形挙動が観察できる．これはデフォルトで使用される線形材料モデルで想定されることである．

パラメトリック材料法則

パラメトリック材料法則がある可能性がある．例えば，ポワソン比や熱移動係数がパラメトリックである可能性がある．この場合は他のパラメトリックの場合と同様に動作する．一つ気を付けなければならないのは，ひずみと応力を計算するときに，材料法則のパラメータ値を指定しなければならないということである．

記号的ポワソン比 でパラメトリック材料法則を作成する：
記号的材料法則からパラメトリック関数を作成する：
いくつかのポワソン比についてのパラメトリック関数を評価する：

この場合，ひずみの計算はポワソン比パラメータ 依存せず，通常の方法で行われる．

ひずみを計算する：

しかし応力の計算はポワソン比パラメータ に依存する．この場合，パラメトリックパラメータは計算に対する実際のパラメータ値を持つ必要がある．これは記号的パラメータ値を実際の値に置き換えることで行える．

パラメトリックのポワソン比の値 をパラメータに代入する：

周波数応答解析

周波数応答解析は，物体の特定の点が周波数範囲を通してスイープにどのように反応するかに関する情報を与える．周波数応答解析の結果は周波数の範囲における実際の変位間の関係を示す周波数応答プロットである．周波数応答解析は調和解析とも呼ばれる．

周波数応答解析を実行する前に，固有モード解析と静的解析を行うと便利なことがある．固有モード解析は，周波数応答解析に関与する臨界周波数とそれに関連するモードを求める．静的解析は周波数成分がないときの最大変位を与える．

固体力学の領域では，周波数応答解析は実質的に以下を解く：

ここで は角周波数， は虚数の初期値， は結果の変位である．，， は固体力学PDEの質量，減衰，剛性である．これらはSolidMechanicsPDEComponentで与えられる．

この場合，上と同様に となるような減衰しない場合を考える．

変数とPDEモデルを設定する：

周波数応答解析を実行するためには，周波数依存の荷重または定数を設定する必要がある．これは境界荷重が関係しない固有モード解析とは対照的である．

周波数依存の荷重を設定する：

残った境界条件は定常境界条件である．

調和PDEを のパラメトリック関数として設定する：

最小の固有振動数から計算される最大の固有振動数までの間の周波数範囲を選ぶ．

使用する周波数のリストを設定する：
周波数を角周波数に変換する：
要求された周波数における調和解を計算する：
幾何モデル上にクエリポイントを赤で可視化する：
計算された周波数におけるクエリポイントの総変位を求める：
総変位と周波数をプロットする：

結果をよりよく理解するために，周波数応答と定常解を比較する．周波数応答解析では，負の 方向で の周波数依存荷重

を使った．定常解析では負の 方向で の同等な境界荷重
を使った．定常解からのクエリポイントの変位を求めると面白い．これによって定常解のクエリポイントに関して周波数応答の変位を求めることができる．

定常の場合のクエリポイントの変位を計算する：
総変位と，前に計算した静的な総変位の周波数をオレンジ色でプロットする：

共振周波数近くでは，クエリポイントの最大変位は定常解析における同じ点の変位より約10倍大きいことが分かる．

クエリポイントの最大総変位の周波数を求める：

前のセクションでは，制約のある本棚の固有値を計算した．

前に計算した最初の数個の固有振動数と比較する：

これは物体が共振状態にある場合，変形の程度は同等の定常力で起こるものよりずっと大きくなる可能性があるということである．

固有モード解析中に計算された共振周波数におけるパラメトリック調和関数を評価することは，可能ではない場合がある．

固有モード解析から調和共振を設定する：
共振周波数では調和PDEは解けない可能性がある：

方程式

概要

固体力学は荷重と制約条件下の物体の解析である．固体力学のトピックはすべて次の4つの成分に関連している：

これらの成分は以下のように関連している：

134.gif

詳細を説明する前に，さまざまな固体力学成分の概要から始めよう．

平衡方程式

平衡方程式は，力密度 と応力 を関連付ける．大雑把に言うと，応力 は適用された荷重に対する内部の点の抵抗である．時間非依存の平衡方程式は以下で与えられる：

ここで は質量密度， は変位ベクトルである．荷重 は物体全体に作用することもあれば境界に作用することもある．簡単にするためにここでは時間依存項を無視する．応力成分 および物体の力 に対する三次元空間の時間非依存平衡方程式 の成分形式は以下で与えられる：

運動学方程式

運動学方程式は変位 をひずみ を関連付ける．ひずみ は物体内の点と点の間の相対的変位を記述する．ひずみ測定にはいろいろある．ここでは2つのアプローチを区別する．微小変形理論では変位とひずみは小さいものと仮定される．これはSolidMechanicsPDEComponentのデフォルトメソッドである．微小ひずみ測定は以下で与えられる：

有限変形理論では仮定はない．ひずみ変位の関係が非線形であるなら，幾何学的非線形性という．有限変形理論はそのような幾何学的非線形性である．幾何学的非線形性は非線形運動学方程式とも呼ばれる．幾何学的非線形性は荷重のかかる間幾何が発展するということの原因である．

構成方程式

平衡方程式（力）と運動学方程式（変形の記述）の関係は構成方程式で構築される．構成方程式は材料モデルとしても知られており，応力とひずみを関連付ける．一般にこれは以下の形式の関係である．

ここで関数 はひずみや他の場の量と応力を関連付ける．最も重要な関係は，応力とひずみの関係は線形であるというフックの法則であり，以下で与えられる：

ここで比例定数 はヤング率である．ひずみが小さいままならば，ほとんどの金属や合金は弾性材料としてモデル化できる．

応力とひずみの間の線形関係が仮定できない場合は，非線形構成方程式になる．これは非線形材料の法則とも言われる．

線形性の他，材料の法則には面白い特性がある．力を取り除くと物体がもとの設定に戻る場合，この物体の変形は弾性と言われる．荷重が除かれても物体がもとの設定に戻らない場合は，塑性変形と言う．物体が塑性変形するときは常にエネルギーが失われる．塑性変形は永久である．

現在のバージョンでは線形または非線形の弾性材料モデルが可能である．塑性は今後のバージョンで考慮される．

平衡方程式

平衡方程式の導出には材料特性（構成方程式）は使われない．したがって，平衡方程式はすべての材料に適用できる．平衡方程式は以下で与えられる：

提示された固体力学公式は，主未知数としての変位に基づいている．有限要素ソルバが変位を計算するので，これは一般的なアプローチである[4, p.480]．ひずみや応力等の第2の未知数は変位から回復される．方程式はすべての未知数について解かれるように設定されるということは考えられるが，結果の方程式系は妥当な時間と利用可能な限られたコンピュータメモリで解くにはとてつもなく大きくなる．

固体力学モデルの項の単位は力密度()である．

物体の荷重

物体の荷重は物体の体積全体に作用し，外部力場から生じる力である．これらの力は体積力とも言われ，力密度（単位は）で指定される．物体の荷重は物体力と呼ばれることもよくあるが，実際は単位体積あたりの力である．物体の荷重はパラメータ"BodyLoad"で指定することができ，ベクトル場として指定される．重力は一定の物体の荷重の例である．また回転する物体の位置依存遠心力は別の例である．

例えば物体の自重をモデル化するためには，物体の荷重が であるような重力加速度 と質量密度 の積を指定しなければならない．これはよくあることなので，追加のパラメータ"BodyLoadValue"を指定することができる．"BodyLoadValue"が指定されると，質量密度 は自動的に掛けられ物体の荷重を出す．

壁に固定された梁の自重による変形を計算する：
変形した梁を可視化する：
の断面における， 方向の梁の変位をプロットする：

運動学方程式

運動学とは，物体の動きや変形を数学的に定式化するものである．見付けられた定式は，それらの変形を引き起こす力とは独立である．変位 とひずみ は次のセクションで紹介し，定義する．

変位と変形

固体力学方程式の解は，それぞれ ，， 方向の変位である3つの変位関数 ，， の集合を与える．もとはにあった物体の任意の点 は，荷重下ではに移動する．この変位は変位ベクトル で集められることが多い．

物体のすべての点が同じ変位を受けるときは変形はない．つまり物体が変形するのは，一様ではない変位があるときだけである．変位だけあって変形がない場合は，物体の剛体運動と言う．物体の回転や移動は剛体運動である．

変位に関係する数量は， で表されることが多い変位勾配である：

剛体の移動では，変位勾配は0である．

変形は隣合せの粒子の相対運動を記述する一方，変位は絶対運動を記述する．変形は大きさや形状の変化をとらえるが，変位は捉えない．

ひずみ

ひずみは，物体内の変形または変位の量を記述し[2]，割合であり単位はない．

回転や平行移動等の剛体運動を受ける物体は，変形しない．剛体運動という名前がこれを表している．剛体運動は純粋な変位である．しかし，ひずみは大きさまたは形の変化の結果である．正しく測定されたひずみでは，ひずみは変形だけを測定するので，剛体運動に対してはゼロであることが分かる．

ひずみは温度のような物理量ではない．ひずみを測定する方法はいろいろある．今日使用されているひずみ測定方法 [13, c. 4.1, p. 90]はほとんどすべて，変形が小さいときは大体同じ結果を与えるように構築されている．小さい変形の場合は，どのひずみ測定法を使っても関係ない．

変位がゼロのときでもひずみが非零になることはある．左で固定され，右で正の 方向に引っ張られる棒を考えてみよう．固定された端では変位はゼロになるがひずみは非零である．

微小ひずみと工学ひずみ

微小ひずみは工学ひずみ，または無限小ひずみとも言われる．

一次元の棒についての簡単なひずみの定義は以下で与えられる．

ここで は変形していない物体の長さ（単位は）， は変形した物体の長さ（単位は），は適用された荷重による長さの変化（単位は）である．

190.gif

上のひずみの定義はWolfram言語で実際に使われているものではなく，直感を養うのに便利である．実際に使われている定義を導くためには[4, p.475]に従う．物体の変形を記述するために，もとの設定の点および最終の変形した設定における同じ点を考える．大きさと形状の変化の特性を示す無限小の立方体を考える． 大きさの変化は立方体の長さの変化で判定する．形状の変化は，辺と辺の間のもとの角度 からの角度の変化で判定する．

まず，長さの変化を考える．もとの形状では線分は 軸に変更であるとする．変位によりこの線分は最終形状のに移動する．

195.gif

点 がにあるとしよう．線分はx軸に平行なので点 はにある．点 は で動かされ，点 はとなる．点 は で動かされ，最終形状での点の位置はとなる．これらの座標からもとの線分と変形した線分の長さは以下で計算できる：

は 軸に平行なので，増分は のみの関数であり，次のように書ける．

代入し直すと以下が得られる．

方向の垂直ひずみは次のように定義される．

ここまでの計算は厳密である．ここで特に微小ひずみを導く決定的な簡約化を行う．変位は小さく，1よりずっと小さいと仮定する．つまり以下となる．

したがって，垂直ひずみは以下で近似される．

方向と 方向の残りの垂直ひずみにも同じことが言える．

あとしなければならないのは，もとはだった角度の変化を数量化することである．これはせん断ひずみの方程式になる．せん断ひずみは角度の変化を数量化する．もとの形状で線分が線分に対してであると考える．最終形状ではこれらはそれぞれ，に移動している．

もとにあった点 はに移動した．変位は小さいという仮定を立てていたし，角度の変化だけを記述して長さの変化を無視したいので，点 は点 に対してだけ移動する．点はに対してだけ右に動く．角度の変化を記述したいときは長さの変化は考慮しない．

では ，では なので，以下のようになる．

これを使うと，次のせん断応力は以下になる．

同様に進めると．

ここで以下のように成分を行列でグループ化することができる．

しかし，これには不利なこともある．上記の物体はテンソルのように動作せず，物体を表す他の簡潔な方法を除外してしまう．工学者たちはテンソルが発明されるずっと前にせん断ひずみを測定していた．これを解決するために，テンソルせん断ひずみ を工学せん断ひずみ のと定義する．

ひずみテンソル成分は対称なので，完全な微小ひずみテンソルは以下のようになる．

または次のようにも書ける．

これが微小ひずみ測定と言われるものである．このひずみの定義は変位と回転が小さいときだけ使えるということを理解するのは非常に重要である．それ以外の場合は別のひずみ測定を使う必要がある．

パラメータで"EngineeringStrain"->False と設定すると，工学ひずみの使用を切り替えることができる．これは真ひずみを利用する，線形材料則と非線形材料則を比較する場合に便利である．詳細はSolidMechanicsStrainの関数ページに記載されている．

テンソルひずみを使うと，テンソルを次のように表すことができる．

ここでベクトル は変位勾配である．

変位ベクトル からひずみテンソルを生成する：

Wolfram言語のデフォルトの材料モデルは線形弾性材料モデルなので，微小ひずみ測定はWolfram言語で使用されるデフォルトのひずみ測定である．モデル化される場合には変形が小さいことが一般的なのでこれが選ばれた．さらに，微小ひずみ測定は線形であり，解くのに長くかかる非線形方程式系には結果的にならない．微小ひずみ測定はコンクリート，剛性プラスチック，金属，ポリマー材料等の線形粘弾性材料，程度な荷重での土や粘土等の多孔質材料の小さい変形をモデル化するのに便利である．実際，荷重が高すぎなければほとんどどのような材料でもモデル化することができる．微小ひずみ測定は一般にゴム材料，軟組織，大きい変形には向かない[13, c. 4.1, p. 95]．

しかし，微小ひずみ測定には限界がある．微小ひずみは剛体回転についてひずみゼロにはならない．これは例を使って説明する．

最初のメッシュを生成する：
座標を抽出する：
回転変換を適用する：
変位ベクトルを作成する：
最初と最後の設定を可視化する：
ひずみを計算する：
における ひずみ値を検査する：

垂直ひずみがゼロではないことが分かる．微小ひずみ測定は，小さい回転のときだけに有効である．回転が大きいときはグリーン・ラグランジュ法のような他のひずみ測定法を使った方がよい．

残念ながらグリーン・ラグランジュひずみはコーシー応力と互換ではないため，微小ひずみを例えばグリーン・ラグランジュひずみにただ置き換えるほど簡単ではない．これについては超弾性のセクションで詳しく説明する．

変形勾配

大きい変形のひずみ測定を考える前に，変形勾配の概念を確立する．小さい変形では，変形前と変形後のオブジェクトの違いは，オブジェクトのサイズと比較して小さいものであり無視される．大きい変形の場合は，そうはいかず変形の説明が必要になる．

まず変形していないオブジェクトと変形したオブジェクトの違いから始める．このセクションでは，文献で使われる表記法と一貫性を持たせるために，変位ベクトル を太字にはしていない．

261.gif

変形していないオブジェクトは，基底ベクトル で座標系に置かれる．この領域は基準配置と呼ばれる．この領域の実体を言及するのに大文字を使うという慣例に従う．このオブジェクトが変形を受けるとき，すべての材料点 が変形したオブジェクトの材料点 に移動する．変形した領域の実体にはすべて小文字を使う．簡単にするために，基準配置のオブジェクトおよび変形配置のオブジェクトには同じ座標系 を使う．基準配置は材料配置，初期配置，ラグランジュ記述と呼ばれることもある．変形配置は空間配置，現／最終配置，オイラー記述と呼ばれることもある．

変形関数 は を にマップする：

無限小の線分 ：

線形化する：

再び挿入すると以下が得られる：

ここで は単位行列， は以下の成分を持つ変形勾配である：

ここで はクロネッカーのデルタである． は行列形式では以下のように書ける：

簡潔に書くと以下になる：

変形勾配 は空間座標および材料座標で表すこともできる．次のように書く：

これは以下に従う：

変形勾配行列 は変形マップ のヤコビ行列である．変形勾配テンソルでは，変形配置の2つの隣接する粒子の相対位置を，基準配置のその相対的粒子位置で記述することができる[17, p. 81]．

変形勾配は非線形固体力学の中心である．

グリーン・ラグランジュ

グリーン・ラグランジュひずみ測定は，微小ひずみ測定の制約となっている小さい変位や回転とは関係のない非線形ひずみ測定である．概念としてグリーン・ラグランジュひずみ測定は次をモデル化する．

ここで は変形しない物体の長さ（単位は）で， は変形した物体の長さ（単位は）である．

グリーン・ラグランジュひずみ測定が小さい変形の制約を受けないことを示すために，微小ひずみのセクションと同じ例をグリーン・ラグランジュひずみ測定を使って考える．

グリーン・ラグランジュひずみを計算する：
における ひずみ値を検査する：

グリーン・ラグランジュひずみの場合，剛体運動で想定されるように，垂直ひずみとせん断ひずみはゼロである．

グリーン・ラグランジュひずみ測定は，以下で与えられる変形勾配 を計算することによって実装される．

ここで は単位行列である．グリーン・ラグランジュのテンソルひずみ は以下で与えられる：

グリーン・ラグランジュひずみ測定は，基準の設定におけるひずみを記述する材料ひずみテンソルである．これは，例えば変形の設定を記述するオイラー・アルマンジひずみ測定とは対照的である．

オイラー・アルマンジ

包括的にするために，オイラー・アルマンジひずみ も説明する．オイラー・アルマンジひずみ測定は，変形の設定におけるひずみを記述する空間ひずみテンソルである．概念的に，オイラー・アルマンジひずみは以下をモデル化する：

ここで は変形しない物体の長さ（単位は）， は変形した物体の長さ（単位は）である．

オイラー・アルマンジのテンソルひずみは次で与えられる：

ひずみ測定

実際に使われるひずみ測定は調べることができる．

ひずみ測定を調べる：

ポアソン比

ポアソン比 は，物体が縦方向に引っ張られたときの横縮みと縦膨張の間の関係を与える．

これまで見てきたことから，弾性域で等方材料が使われると仮定するなら，以下が言える：

等方性材料はすべての方向で同様に動作する材料である． 方向が縦方向， と 方向が横方向とすると，弾性領域では以下が言える：

材料が等方性 であり， to および に対する の関係はポアソン比なので以下が成り立つ：

等方性線形弾性材料では，ポアソン比 はの間の値になる．ほとんどの材料では，ポアソン比は約1/3である．多くの材料のポアソン比は0.2-0.3である．極端な例を挙げると，コルクのポアソン比は0である．ゴムのポアソン比は1/2に近い．ポアソン比の値 は材料が非圧縮性であり，数値的シミュレーションでは問題を引き起こすということである．オーセティックス材料等の人工材料は負のポアソン比を持つこともある．オーセティックス材料が一方向に伸張されると，別の方向では厚くなる．

構成方程式

構成方程式は，材料モデルを介して平衡方程式と運動学方程式を関連付ける．このことについて話す前に応力の概念を説明する．

応力

応力は物体内部の内力の分布を記述する[2]数量であり，あるいはで測定される．

ここで は力（単位は）であり，は力が作用する初期面積（単位は）である．

327.gif

物体の応力を求めると，いつその物体が破損するかを予想することができるため，その計算は重要である．ある材料が最大応力に耐えられるとすると，となったときにその材料は破損する．材料が破損する前に適用することのできる最大力は以下である．

詳細は破壊理論セクションをご参照いただきたい．

応力・ひずみ関係

応力とひずみの関係は材料特性であり，応力・ひずみ図で可視化することができる．この図のデータは引張試験のものである．この試験では，2つの締め具のある機械に材料の標本が固定され，引っ張られる．破損するまで，適用される力と生み出されるひずみが記録される．

冷延鋼板の引張試験から未加工の測定データをインポートする：

少なくともデータは，機械によって適用された荷重と測定されたひずみを含む必要がある．

荷重とひずみのデータを抽出し，スケールする：

応力を計算するためには，標本の直径が必要である．

荷重と標本の断面積から応力を計算する：

最後のデータ点は通常，材料の破裂の記録であるため削除する必要がある．

材料データの応力ひずみ曲線を作成する：

上は冷延鋼板の材料データであり，これは延性材料である．脆性材料（アルミニウム等）ならば上部の平らな部分はもっと小さい．

331.gif

曲線の下の面積は，吸収されたエネルギーに相当する．

応力ひずみ曲線は，測定されたひずみと力，つまり適用された応力を関連付ける．典型的な応力ひずみ曲線には議論に値する部分がいくつかある．以下は延性材料についての誇張された応力ひずみ曲線である．いくつかの特徴のある点に印が付いている．特に点A，B，Cは互いに非常に近いことがあり，AとCの結合が直線になることもある．応力ひずみ曲線のこのような点で物理的に起きていることは，内部結晶境界で材料が滑るということである．これは非常に速く起こるため検出が難しい．

以下の点に印が付いている：

最初の領域は線形領域である．ここでは，応力とひずみの関係は線形でありフックの法則として知られている．この領域に適用される荷重は完全に可逆的である．この領域は弾性域とも呼ばれる．この領域では以下のようになる．

ここで はヤング率であり材料特性である．ヤング率は弾性率とも呼ばれる．ヤング率は応力と同じ単位を持つ．ヤング率の材料特性は単軸張力試験で測定される．

真応力と真ひずみ

応力・ひずみ曲線は引張試験から得られる．この試験では標本が力によって引き離され，このことが応力を適用し，ひずみが記録される．試験の間に標本は変形し，結果として断面が変化する．しかし断面は応力の方程式式に以下のように入る：

試験中に瞬間的な断面を測定するのは難しい．そのため標本の最初の形だけが記録される．しかし真応力と真ひずみはこの形の変化を考慮に入れる．

2Dの横断面としてここに示す長さ ，面積 の変形していない立方体を考える．この標本は力 を受け，その力によって変形し，現在の長さは ，面積は である．

341.gif

真応力 は以下で与えられる：

工学応力 は以下で与えられる：

試験標本の応力の測定で違いが現れる．通常標本のもとの面積だけが考慮される．ひずみについても同じである．真ひずみは以下で与えられる．

工学ひずみは以下で与えられる：

次に，標本の体積は一定のままであると仮定する．弾性域の体積変化は小さいので，この仮定は弾性域では妥当である．塑性変形の最初の部分では材料は圧縮できないと考えられるため，この過程は塑性領域の最初の部分でも妥当である．しかし，くびれ現象が始まった後は，塑性領域の2番目の部分においてはこの仮定はもはや妥当ではない．

一定の体積を想定する：

体積が一定であるという仮定のもとでは，以下の変換が構築できる：

また，小さいひずみ値では，弾性域のときは塑性領域のときに比べて真応力・ひずみと工学応力・ひずみの差は小さい．顕著な塑性変形がある場合は，真応力・ひずみ曲線を使わなければならない．

上の応力・ひずみ関係のセクションからのデータを使って，工学応力と工学ひずみを真応力と真ひずみに変換する．

測定データから真ひずみを計算する：
測定データから真応力を計算する：
同じ材料データに対する工学応力・ひずみ曲線と真応力・ひずみ曲線を作成する：

Wolfram言語は工学ひずみを使う．

線形弾性材料モデル

構造方程式は材料モデルとも呼ばれる．構造方程式は応力とひずみがどのように関係しているかを記述する．線形の場合，応力 とひずみ は弾性行列 を使ったフックの法則の一般化で関連付けられる：

弾性行列は構成行列，または材料剛性行列とも言われる．

最も一般的な場合，線形弾性構造方程式の弾性行列 は以下で与えられる：

はせん断応力である．これを使って，工学せん断ひずみ が使われたことを示す．

弾性定数の数 は減らすことができる．この簡約化はいくつかの仮定を作成することで行える[4, p.478]．弾性行列はほぼ例外なく対称であると想定される．この仮定により36個の弾性定数は21個に減る．

行列がある平面について対称性を示す場合は，さらに減らすことができる．弾性行列 の厳密な形式はこれらの対称性に依存する．通常等方性材料と異方性材料に区別する．等方性材料では，材料特性はすべての方向で同じである．ほとんどの金属は等方性材料である．等方性ではない材料はすべて異方性材料と言われる．多くの場合異方性とは材料が可能なすべての方向で異なった挙動を示すことを意味する．

一般的な異方性の場合は実に一般的であるため，名前の付いた下位の場合が存在する．例えば直交異方性材料では，材料特性は物体のそれぞれの軸方向で異なる．木材は，材料特性が木目の方向に依存する直交異方性材料と考えられる．

一般的に異方性材料は，材料特性が一部またはすべての方向で異なる複合材料または生物組織が多い．また，織りが全方向の材料特性に影響を与える繊維教科材料は異方性と考えられる．

幾何モデルの方向が材料特性の方向と合致していなかったら，方向行列を指定する必要がある．この手順は材料の方向のセクションで説明する．

線形弾性材料はほとんどの工学設計計算に適しており，成分は降伏を超えてはならない[11, c.1.1.5]．

以下のセクションでは，これらの線形材料モデルを使う．

構成方程式と運動学方程式の両方が線形である場合は，方程式は以下のように簡約することができる：

ここで および は弾性行列 を完全なテンソル形式 に変換する．完全なテンソル形式は，SolidMechanicsPDEComponentの出力で見られる形式である．

これについては「SolidMechanicsPDEComponentの出力」のセクションで詳しく説明する．

等方性線形弾性材料

等方性線形弾性材料モデルはWolfram言語で使われるデフォルトの材料モデルである．材料がすべての方向で同様にふるまう場合に，その材料を等方性という．換言するとその材料特性は方向依存なのである．等方性線形弾性材料モデルはわずかな変形と低い荷重を受ける多結晶金属，セラミック，ガラス，ポリマーに適している[11, c.1]．

等方性材料の弾性行列は以下で与えられる：

弾性行列を指定するのに必要な弾性係数は，ヤング率 とポアソン比 だけである．

等方性材料を記号的に設定する：
記号的固体力学PDE成分を生成する：

これは項の線形の形である．

その他の材料パラメータ

等方性線形弾性の場合，広く使われるヤング率とポワソン比は他の弾性率で表すことができる．

体積弾性係数をヤング率とポワソン比で表す：

現在，さまざまな弾性係数を他のもので表すという機能は3Dでのみ利用できる．

利用可能な弾性係数を求める：

どの弾性係数も他の2つの係数で表すことができる．このことで，一般的な弾性率のどれでも使うことができ，SolidMechanicsPDEComponentはその操作に必要な弾性率を見付ける．

ラメパラメータとせん断弾性係数として与えられた材料パラメータで，等方性材料を記号的に設定する：
選ばれた弾性係数から記号的固体力学PDE成分を作成する：

直交異方性線形弾性材料

直交異方性材料は3つの次元で異なる特性を持つ．直交異方性材料としてモデル化できる例として木材がある．丸太を考えてみよう．木は木目に沿って最も硬く，円周方向に沿ってやや硬く，半径方向で最も硬さが弱い．

直交異方性材料の一般的な形式は以下で与えられる：

簡単にするために，通常弾性行列 の係数は明示的に与えられず，コンプライアンス行列 を介して定義される．ここで がある：

弾性行列 は となるようにコンプライアンス行列 を転置することで求められる．

直交異方性材料を記号的に設定する：
記号的な固体力学PDE成分を作成する：

これは項の線形の形である．

横等方性線形弾性材料

横等方性は直交異方性の特殊な形である．直交異方性材料は3つの直交面について対称であるが，横等方性材料は1つの直交面について対称である．

横等方性材料は微小の薄層の物理的特性と同じものを持つ．堆積岩がその例である．堆積岩はさまざまな層，つまり 方向と 方向の層からなる．硬度のような材料特性は 方向と 方向で同じである．しかし 方向では物理的特性は変化し，材料はその方向について直交異方性となる．

還元すると，横等方性材料は，材料が対称となる中心の対称軸を1つ持っているということである．つまり，その軸について材料を回転させると，軸が垂直となる平面における等方性材料のようにふるまうということである．横等方性材料は等方性材料と直交異方性材料の間に位置する材料区分である．

横等方性物体の別の例として生体膜がある．この場合，膜の材料特性は膜の平面と同じであるが，垂直方向の膜とは異なる．

385.gif

堆積岩の材料特性は 方向と 方向で同じであり，これが等方性平面を形成する． 方向では，材料特性は異なる．

横等方性材料の部類に入る2つ目のタイプの材料は，同じ方向に流れる繊維の族によって強化される材料である．この場合，繊維が材料の対称性の軸の役割を果たす．繊維強化された材料は，通常横方向よりも繊維の方向の方がずっと強い．

389.gif

材料行列に1種類の繊維が埋め込まれた，繊維強化された複合材料．

繊維強化された材料と堆積岩との違いは横断面の構造である．堆積岩は等方性面において均質であるが，繊維強化された材料は異なる．

390.gif

左側は繊維強化された材料を前から見たところで，右側は2Dの横断面を示す．この横断面は等方性面であるが，厳密にいうと等方性面が均質ではないことを示している．

では，材料パラメータはどのようにして得ればよいのだろうか．ある複合材料の材料パラメータはそれに関連した文献で見付けることができる．例えばファイバーグラスでは[25]を見るとよい．材料のパラメータがどこにも見付からない場合は実験を行ってみるのが最善の方法である．しかし，これは常に可能であるとは限らない．行列および繊維の材料パラメータが既知ならば，すべてが失われる訳ではない．一つの方法として，均質化シミュレーションを使って，複合材料の材料パラメータを得るというものがある．そのプロセスについては[26]を，均質化理論およびそれに関連したテクニックの概説は[27]をご覧いただきたい．別の方法として，材料と繊維の材料パラメータを使うというものもある．行列の材料パラメータは横方向で，繊維のパラメータは軸方向で使うことができる．

横等方性の場合は一方向にのみ走る繊維にのみ適しており，織り繊維には適していない．

横等方性材料の仕組みの概要は[16]で見ることができる．繊維強化材料に焦点を当てた概要は[24]に記載されている．

横等方性は直交異方性の特殊タイプなので，9個の独立成分を持つ直交異方性材料の一般形から始めることができる：

対称軸が 方向である横等方性材料を考える：

独立成分の数が5個に減った．コンプライアンス行列は通常以下のように書ける：

ここで はヤング係数， はポワソン比， は条件 の剛性率である． 成分が2回現れていることに注目のこと．この関係は等方性材料について言えるものと似ている．

および は軸の方向の材料パラメータなので，下付き文字 を使って表すことができる．材料パラメータ および は横方向に関連付けられているので，下付き文字 を使って表すことができる．これによって以下の式が得られる．この式では材料パラメータが横方向と軸方向の2つのグループに分けられている：

ここでも条件は である．

次に例を示す．

横等方性材料を記号的に設定する：
記号的な固体力学PDE成分を生成する：

条件 は自動的に強制される．しかし， 方向にカスタムの剛性率を規定することもまだ可能である．

カスタムの剛性率 を持つ横等方性材料を記号的に設定する：
記号的固体力学PDE成分を生成する：

異方性線形弾性材料

異方性材料では，材料特性はすべての材料方向において異なる．異方性線形弾性材料モデルは，強化複合材，木材，金属単結晶およびセラミック単結晶に向いている[11, c.1]．

異方性材料を設定するためには，パラメータ"ElasticityMatrix"を持つ完全な弾性行列 を提供する．

異方性弾性行列を指定する：

完全な異方性弾性行列は，順番がのフォークト記法で指定しなければならない．

記号的な固体力学PDE成分を作成する：

弾性行列の逆行列（コンプライアンス行列）しか使えない場合は，パラメータ"ComplianceMatrix"で指定することができる．それで弾性行列得るための転置が自動的に行われる．弾性行列を指定すると，コンプライアンス行列は上書きされる．

または，完全な弾性テンソルを指定することもできる．

弾性テンソルを指定する：
記号的な固体力学PDE成分を作成する：

これは項の線形の形である．

弾性行列の例は平面ひずみのセクションに記載されている．複合材料のコンプライアンス行列の設定方法の例はSolidMechanicsPDEComponentに記載されている．複合材料弾性テンソルにも同じ手順を使うことができる．

座標軸に平行でない材料

ここまでは，対称軸が ，， 軸に沿っている直交異方性材料だけを見てきた．対称軸が 方向であると仮定した横等方性材料もそうである．しかしこの仮定を常に満たすことはできない．例えば木目やコンクリートの鉄骨はさまざまな向きであるためモデルの設定がより複雑になる．また，内部構造がある点から別の点までで異なるという理由で軸に平行ではない材料もあり得る．ここで示す方法を使うと，例えば曲線の鉄骨を持つコンクリートを考えることができる．このような構造はどの軸にも平行にはならない．立体を 方向に合わせる代りに材料モデルを木目の方向に合わせた方がずっと簡単であることを示す．

この問題を解決する方法として2種類のアプローチがある．一つは弾性（コンプライアンス）行列を変更して，材料の向きを考慮に入れるようにするというものである．このアプローチは[1]で示すものであり，カスタムの弾性（コンプライアンス）行列を使って"MaterialSymmetry" 設定を"Anisotropic"にする．

2つ目のオプションはもとの弾性行列を持つもとの方程式が，適切な向きの座標系で局所的に使えるという考えに基づいたものである．ここではこちらのアプローチについて詳しく話す．簡単にするために例として横等方性材料を考えるが，このアプローチは直交異方性材料や完全異方性材料にも使うことができる．

線形弾性方程式は以下のように書くことができる

ここで は応力， はひずみ， は弾性行列である．各点で，材料モデルの局所的な向きが 方向になるように座標系を選ぶことができる． 方向は適切な回転行列を使って，座標系を固定された座標系に変換される．次にモデルは固定された大域的座標系で設定される．材料は局所的に強制されるので，材料の横等方性は満足される．

軸と整列しない材料対称性の場合，回転行列 は整列しない物体の各点において弾性行列の向きが 方向になるように、固定基本座標系から弾性行列 を回転させる．モデルの向きは横等方性材料の場合対称軸によって与えられる．

固定座標系を取り，対称軸の向きを の参照設定で表す．回転行列 があり，対称軸を 方向に合わせるように固定座標系を回転させるので以下が得られる：

この回転された座標系は特権を持つ座標系，実験室座標系，材料座標系と呼ばれることもあり，この座標系の数量は上付き文字 で表す．以下の関係を使うと，微小ひずみテンソル と応力テンソル を大域的な固定基底から局所的な特権を持つ基底に変換しあうことができる．

ここで は局所的な特権を持つ基底の微小ひずみテンソル， は局所的な特権を持つ基底の応力テンソルである．回転行列 は一定である必要はないので，空間の関数になることができる．これによって向きがさまざまな材料に対処できるようになる．

と を使うと，線形弾性構成方程式(1)を使うことができる．特権を持つ基底では材料は局所的に適切な向きを向いているからである．これで以下の関係が得られる：

横等方性材料の場合は，これは以下のように変換できる：

以下の例は回転行列の使い方を示すものである．ここでは一方で固定され他方が引っ張られる立方体を考える．材料は横等方性であり2つの実験を行う．最初の実験では，材料の対称軸が一般的なデフォルトの向きである 方向に平行になるように向いている．2つ目の実験では 軸に平行になるように向いており，回転行列を使ってそのような材料の向きを得る．どちらの場合も，対称軸は簡便さのため得られた結果の比較と検証が簡単な定数とする．

まず材料がデフォルトの向きと考える．

立方体領域を定義し，材料のデフォルトの向きをプロットする：
変数を定義する：
材料パラメータを定義する：
パラメータを定義する：
方程式を設定する：
方程式を解き，変位を得る：
変位をプロットする：

次に同じ材料を使うが，今回は向きを変える．このために回転行列を作成する．

回転行列を定義し，方向のある材料をプロットする：
方向のある材料の材料パラメータを定義する：

重量な点は，材料パラメータがデフォルトの向きの材料と全く同じように定義されているという点である．材料パラメータを定義するすべての行列が， 方向の対称軸を持つ材料で同じである．これは材料パラメータが大域的な固定基底ではなく局所的な特権基底について定義されるからである．

回転行列を設定する他に，工学ひずみがオフになっている点に注目する．これは回転行列が使われるときは常に必要なことである．この設定はWolfram言語では自動的に行われるが工学ひずみがオフであるという警告メッセージが出力される．手動でオフにすると警告メッセージが出ない．

方向付きの材料についての方程式を設定する：
方程式を解いて変位を求める：
変位をプロットする：

どちらの変位プロットも予想通り，同じ変形で回転したものを示している．例を簡単にするため，両方の材料が同じように振舞うことをチェックすることができる．

計算領域の回転を定義する：
向きのある材料の変位を回転させてプロットする：

2つ目の実験のプロットを回転すると，最初の実験の向きのない行列の材料と同じ変位になる．これは変位を引いて結果を最大化することで検証できる．

変位の差のノルムを最大化する計算を行うことで，変位が同じであることを検証する：

回転行列に基づく変換を利用する場合に，考えなければならないことがいくつかある：

材料特性は，固定された大域的座標系ではなく局所的な特権座標系で指定する必要がある．所望の結果が大域的な固定基底なのでこれは不自然に思えるかもしれないが，大域的な固定基底で向きが変化する材料の材料特性を指定することは非常に難しいので，必要なことである．

次の例は非定数の回転行列を使う．ここでも横等方性材料を考えるが，今回は対称軸がさまざまである．単位立方体を一端で固定し，他方には圧力をかける．

方向変換行列 を設定する．これは材料の向きを記述する：

方向変換行列は空間座標の関数である．

領域を定義し，材料の向きを可視化する：
材料パラメータを設定する：
線形弾性に対する変数とパラメータを設定する：

"EngineeringStrain"をFalseに設定もした．このバージョンでは方向変換行列は工学ひずみと一緒には使えない．工学ひずみは線形弾性材料ではデフォルトのひずみであるがこの場合はそうではない．

モデルを設定する：
変位を計算する：
変形した物体を可視化する：

線形弾性材料

実際に使われる線形材料モデルは調べることができる．

材料モデルを調べる：

線形弾性構成方程式の一般化

線形弾性方程式はいくつかの方法で一般化することができる．初期ひずみ または熱誘導されたひずみ を考慮することや初期応力 をモデル化することができる．これらを加えると，構成方程式は次のようになる：

次のセクションでは，さまざまな成分を説明する．

初期ひずみ

初期ひずみは，温度等の物理現象によって引き起こされるひずみにおけるオフセットをモデル化するために使うことができる．初期ひずみは計算されたひずみから引いたものなので以下のようになる：

初期ひずみを設定する：

右にスクロールすると，ひずみの効果が分かる．

複数材料モデルの場合，初期ひずみはIf分等の分岐コンストラクトを使って設定することができる．

複数材料の初期ひずみを設定する：

ElementMarkerの生成と使用法は要素メッシュの生成チュートリアルで説明してある．

初期ひずみは別々に計算することができる：
初期ひずみを見る：

材料の方向について指定したように，指定可能な方向変換行列は初期ひずみに自動的に適用される．

熱弾性

熱誘導による弾性は追加のひずみとしてモデル化される．熱ひずみ は初期ひずみの一般的な形式であり以下で与えられる：

熱ひずみ

温度変化は熱ひずみを引き起こす．熱ひずみは，次の場合熱応力になることがある．

物体の温度変化には体積の変化，その反応として長さの変化がある．温度で引き起こされる長さの変化は，ひずみの変化，熱ひずみ でモデル化することができる．熱ひずみは， が温度のないひずみとなるようにひずみ から引かれる．熱ひずみのデフォルトモデルは温度とひずみの間の線形関係として与えられる．

ここで は，で測定される線形熱膨張の定数係数(CLTE)であり， はで測定される物体温度，は物体内に熱ひずみのない，Kで測定される基準温度である．

熱ひずみがシミュレーションに適しているかどうかを推定する[11, c.1] ためには，以下を使うことができる．

ここでとはそれぞれ物体の想定される積 の最大値と最小値である． が物体内のかなりの割合のひずみ ならば，熱ひずみを考慮する必要がある．過渡熱伝導解析が必要かどうかを決める際にも同様に推定することができる[11, c.1]．特性長 の立体は時間 で定常状態温度に達する．

ここで は比熱容量， は質量密度， は熱伝導率である．一定熱流束および より大きい時間スケールのときは，定常状態解析で十分である．

自体は温度依存であることを認識する必要がある．が成り立ち，熱膨張 の線形係数が温度範囲上であまり大きく変化しない限り，デフォルトモデルが適している．「温度範囲上であまり大きく変化しない」の詳細は，特定の応用で必要な確度によって異なる．高い確度が必要なときは，温度依存の が必要となる．これは後で説明する．最初の例を分かりやすくするために，熱膨張係数が温度に依存しない設定を選ぶ．この場合熱解析は必要ない．

例として， 軸上で長さ ，半径 の円筒が両端で制約を受けるとする．円筒が熱荷重を受けたとき 方向には移動できず および 方向にしか膨張できない．この設定の場合，エンドキャップに作用する反力についての解析解が利用でき，結果の検証に使うことができる．

長さ ，半径 の円筒のメッシュを作成する：
熱ひずみの温度 と熱ひずみの基準温度 を設定し，円筒が80度熱せられるようにする：

温度差がで指定される場合のみ，設定に"ThermalStrainTemperature"を使わなければならないということを覚えておかなければならない．そうしないと"ThermalStrainReferenceTemperature"が0に設定されるからである．

方向の円筒の運動を制約する：

方向の動きを固定することは円筒を完全に制約するには不十分である．それは従属変数 によって表される変位を制約するだけだからである．変位従属変数 と も制約するために，どのような運動も不可能になるように点とを固定する．

，で円筒を固定する：
PDEモデルを設定する：
方程式を解く：
総変位を可視化する：

結果を検証するために，制約された面の一つで軸反力を計算することができる．これは以下で求められる．

ここで はヤング率， は線形熱膨張の係数， は物体の温度，は物体内に熱ひずみのない基準温度である． は当該表面積である．この形状は離散化され結果の検証では離散化誤差は考慮しないので，離散表面の面積が使われる．

形状の一つの面の面積を計算する：
円筒の一方のキャップにおける予想される反力を計算する：

次に円筒の一方のキャップにおける反力を計算する．

反力を計算する：
予想される反力と計算された反力の差の割合を計算する：

直交異方性材料では，熱膨張係数はそれぞれの方向で指定しなければならない：

直交異方熱膨張係数で直交異方性材料を記号的に設定する：
記号的固体力学PDE成分を生成する：

右にスクロールすると，熱ひずみの効果が分かる．

異方性材料では，熱膨張係数はそれぞれの方向で指定しなければならない：

異方性熱膨張係数で異方性弾性行列を指定する：
記号的固体力学PDE成分を生成する：

右にスクロールすると，熱ひずみの効果が分かる．

熱ひずみは別々に計算することができる：
熱ひずみを見る：
結合熱移動方程式

前の例では，熱分布は領域全体を通して一定であった．次の例では一定ではない温度場を考える．"ThermalStrainTemperature"を前の熱シミュレーションの補間関数として指定すること，または完全結合固体力学熱モデルを生成することによって，温度場を固体力学PDEモデルと結合させる．変形した物体が温度分布に影響を与えず，結合が熱場から固体力学への方向でのみ起こるのであれば，補間関数を使うことが推奨されるアプローチである．"ThermalStrainTemperature"ソースとしてInterpolatingFunctionオブジェクトを使うことは，順次に場を解く最大のメモリ必要量が完全結合PDEよりも少ないという利点もある．完全結合PDEはより一般的で強力である．必ずしも必要な訳ではないがここで示す．

次のモデルの例は，熱電対として使われるバイメタルストリップをモデル化する[9, p. 296]．バイメタルストリップは上下に結合された2つの金属からできている．金属は線形熱膨張係数が異なる．バイメタルストリップは，左側で100に保たれる．ストリップの上面と下面はHeatTransferValueにさらされる．左側で，構造体は壁に固定される．線形膨張係数が異なるため，バイメタルストリップは曲がると想定される．最大のたわみを求める．

複数材料領域を設定することについての詳細は複数材料のセクションに記載されている．

長さ ，高さ ，幅 の複数材料の幾何モデルを作成して可視化する：
完全メッシュを生成する：
熱移動の変数名を設定する：
複数材料のストリップを設定する：

これで"ThermalStrainTemperature"は従属熱変数 となった．値 は結合熱移動モデルで計算される．熱移動モデリングについての詳細は，熱移動モノグラフをご参照いただきたい．"ThermalStrainReferenceTemperature"は熱ひずみが想定されない温度である．

PDEを設定する：
境界条件を設定する：
方程式を解き，監視時間とメモリ使用量を求める：
断面 における温度分布を可視化する：
変形を可視化する：
方向の最大変位を求める：

ここで最大変位は負の 方向であるため，負の符号を持つ．

熱移動と固体力学を組み合せた同様の例として，ディスクブレーキの熱分析と構造解析と鉄骨の熱負荷がある．

熱膨張の温度依存係数

ここまで熱膨張係数は線形で定数であった．線形のアプローチは材料依存の温度範囲でのみ利用できる．範囲がもっと大きかったり，より正確な結果が求められる場合は，線形係数を想定することでは不十分である．通常，材料供給者がデータシートにこの情報を集める．この時点では，指定された温度における熱ひずみを計算するのに，異なる方法を使わなければならない．次のような方法がある：

熱ひずみ関数自体を使うことは，モデルパラメータセクションで"InitialStrain"を指定することで行える．"InitialStrain"の設定は初期ひずみのセクションで説明するので，ここでの説明は避ける．

熱膨張係数の定義には，正接係数と正割係数の2種類ある．Wolfram言語のPDEモデルでは熱膨張の正割係数を使う．このセクションでは熱膨張係数のさまざまな定義，それらの間の変換方法，固体力学PDEモデルでそれらを使う方法について述べる．

以下の図は，熱膨張の正接係数 と正割係数 の違いを理解するのに役立つ．

523.gif

熱膨張の正接係数 は次で定義される：

これは熱ひずみが以下であることに従う．

熱膨張の平均係数 を次のように定義するとする．

すると，熱ひずみは以下のようになる．

Wolfram言語のPDEモデルでは，ではなく を与えることが必要である．

指定された領域上のひずみ温度が線形ならば，正割係数と正接係数は等しい．

方向変換行列を使った熱ひずみ

等方性材料の向きが幾何学モデルの軸に従わない場合，材料モデルに正しく投影するように方向変換行列を指定することができる．熱膨張がある場合も同じことが自動的に起こる．

等方性材料では，熱膨張係数はそれぞれの方向で指定しなければならない：

特権を持つ座標に回転させるために，熱ひずみテンソル は固定基底から特権を持つ基底に変換される：

ここで は特権を持つ基底の熱ひずみテンソルである．回転行列 は定数である必要はなく，空間の関数でもよい．これによって方向がさまざまな材料が扱えるようになる．以下の関係が導かれる：

異方性熱膨張係数を持つ異方性弾性行列を指定する：

熱膨張と方向変換行列を持つ記号的な固体力学PDE成分を作成することにより，大きい記号式が生成されるがここでは示さない．

記号式を熱ひずみだけについて調べる．

熱ひずみは別に計算することができる：
熱ひずみを見てみる：

平面ひずみのセクションに例が挙げてある．

初期応力

初期応力は構成方程式に初期応力項 を加えることでモデル化される．

初期応力を設定する：

右にスクロールすると，出力における初期応力の効果が分かる．

初期応力（残留応力とも呼ばれる）は，考慮する物体が製造過程で大きなひずみや高温にさらされたときに出現することがある．物体はこの初期応力を指定することによってモデル化できる圧縮応力状態にすることができる[10, p. 77]．

複数材料モデルの場合，初期応力はIf文等の分岐コンストラクトを使って設定することができる．

複数材料の初期応力を設定する：

ElementMarkerの生成と使用法は，要素メッシュの生成チュートリアルで説明してある．

初期応力は別々に計算することができる：
初期応力を見る：

平面ひずみ，平面応力，軸対称モデル

このセクションでは，結果的に二次元空間モデルになる固体力学PDEモデルの簡約化に関心を向ける．つまり，領域次元と組込み次元の両方が二次元である，理想化された二次元の立体を扱うのである．これは，領域と組込みの次元が両方3である三次元立体，または領域次元が3であるが組込み次元が2であるシェル要素とは対照的である．

2Dの簡約化を使う一番の理由は時間である．通常2Dモデルの方が設定が速く，解くのに少ないリソースで済む．しかし，2Dモデルは簡約化であるため，完全な3Dモデルで得られるほどの詳細をすべて提供することはできない．

このセクションでは，平面ひずみと平面応力のPDEモデルを紹介する．これらは完全な固体力学PDEモデルを2Dで使うことを可能にする簡約化である．注：これらは2Dモデルであるが，三次元の成分を含んでいる．2Dモデリングでは，3Dオブジェクトの境界についての仮定を強制される．これらのモデルの出どころおよびその限界を理解することが重要である．平面ひずみモデルの典型的な例は，非常に厚い物体の横断面であり，平面応力モデルの典型的な例は薄い板である．

2DのPDEモデルの導出は3D等方性弾性モデルから始める．直交異方性および異方性の2Dモデルは現在サポートされていない．

平面ひずみ

壁の間に固定された棒のように，3Dオブジェクトの運動が 方向に制約され，すべての力が 平面に作用する場合，平面ひずみの状態が利用できる．これは3Dオブジェクトの 方向の端点も 平面に制約されているという意味ではない．実際そうではないのである．これらの点は， 方向だけに制約されるのである．もしそうではなくて， 方向の端点も 平面に完全に制約されるのであれば，そのときに限り，オブジェクトは，平面ひずみモデルが有効となるための物体の 平面の膨張と比較して， 方向に非常に長くなければならない．ダムやシャフトの解析はこのカテゴリーに入る．

平面ひずみモデルの最初の仮定は， 方向には変位がないということである．つまり以下を意味する．

この仮定をひずみ測定値に挿入すると， 方向のひずみは以下のように0となるという結果になる．

構成則から始めて，応力 とひずみ を で関連付ける． 方向に変位がないという仮定をひずみ測定に挿入する． 方向のひずみ は以下のように0となるという結果になる．

簡約化すると以下になる．

これは次のように書くことができる．

ここで理解すべき非常に重要な点は， 方向のひずみがゼロに設定されていても， 方向の垂直応力は0ではないということである．現行バージョンのSolidMechanicsStressは平面ひずみの場合の 方向の応力を回復しない．これは手動で行う必要があるが，以下に示すように関数SolidMechanicsExtendedStressを使うこともできる．

で熱誘起ひずみがある場合は，ひずみベクトルは以下になる．

上記と同じ解析が行われる．以下の

および熱ひずみ

が導かれる．次の例は"PlaneStrain"モデル形式を使う方法を示す．モデル領域はパイプの4分の1断面である．同じ単位で内半径 ，外半径 ，厚さ である．左境界にはパイプが上下に動ける対称制約があり，右底にはパイプが左右に動ける第2の対称制約がある．圧力 はパイプの中で外向きに作用する．残りの境界は自由に動く．ヤング率は で，ポアソン比は で与えられる．

578.gif

モデルの変数とパラメータを定義する：
2Dの定常状態固体力学モデルを設定する：

4分の1のパイプ構造は，左側で 方向に対称条件を使う．

方向の左側の対称条件：
方向の底の対称条件：

内部には外向きの20単位の圧力がある．

内側で外向きに作用する境界荷重を設定する：

残りの側面は自由に動く．

幾何モデルを設定する：
PDEを解いて，ひずみと応力を計算する：
変形の誇張表現と4分の1のパイプのもとの位置の辺のフレームを可視化する：

上で述べたように平面ひずみの場合， 方向のひずみ成分 はないが， 方向の応力成分 はある．これらは手計算で求めることができるが，簡便性のためにこれらの成分を計算する関数がある．この関数は，次元が2×2の構造化配列である，計算された変位，ひずみ，応力を取り，それらを次元3×3の構造化配列に変換する． 成分は{3,3}の位置に見付かる．

拡張されたひずみを計算する：

平面ひずみの設定では， 方向のひずみ成分は0である．

拡張されたひずみを抽出する：
拡張された応力を計算する：

平面ひずみの設定では， 方向の応力成分は以下で与えられる：

拡張された応力を抽出する：

次の例は，熱膨張と方向変換行列が適用された異方性の設定における"PlaneStrain"モデル形式を使った例である．モデル領域は管の四分の一の横断面である．内半径 と同じ単位で外半径は ，厚さは である．右下では，管が左右に動くように対称制約条件が加えられている．および では物体は固定されている．残りの境界は自由に動く．ポワソン比は である．

597.gif

異方性弾性行列を指定する：
熱膨張行列を指定する：
モデルの変数とパラメータを定義する：
2Dの定常状態固体力学モデルを設定する：

四分の一の管の構造では，左側の 方向の対称条件を使う．

方向の下の対称条件：
および における固定条件：

残りの辺は自由に動く．

形状を設定する：
PDEを解き，ひずみと応力を計算する：
変形した物体を可視化する：
拡張したひずみを計算する：

平面ひずみの設定では， 方向のひずみ要素は0である．

拡張したひずみを抽出する：
平面応力

平面応力の導出は平面ひずみの場合と似ているが，後で説明するように非常に異なる場合をモデル化する．平面応力の場合， 方向の応力はすべて0と仮定する．これは近似であり，平面応力モデルは物体の厚さがゼロに近付くほど正確になる． 方向の応力をすべて0に設定することは，これらの表面の境界荷重が0であることに等しい．この境界が近ければ近いほど，平面内の条件はより正確になる．

平面応力の導出は，逆の応力・ひずみ関係から始める．

これは以下のように書ける．

応力・ひずみ関係として表すと，以下のようになる．

方向の垂直ひずみは0ではない．現行バージョンのSolidMechanicsStrainは，平面応力の場合に 方向のひずみを回復しない．

で熱誘起ひずみがある場合，ひずみベクトルは以下のようになる．

次の例は複数材料で熱膨張がある"PlaneStress"モデルの形の使い方を説明する．

境界メッシュを生成する：
赤で微調整した境界メッシュと領域を可視化する：
メッシュ微調整関数を作成する：
材料領域に領域マーカーを使ったメッシュを作成する：
メッシュを可視化する：
変数と複数材料パラメータを作成する：

ElementMarkerの生成と使用法は，要素メッシュの生成チュートリアルで説明してある．

左側と底に対称条件を持つPDEモデルを設定する：
方程式を解く：
変形した幾何モデルを赤で，もとのモデルを黒で可視化する：
ひずみを計算する：
垂直ひずみとせん断ひずみを可視化する：
Compute the stress:

平面応力の場合， 方向のひずみ成分 はあるが， 方向の応力成分 はない．これらの値は手計算で求められるが，簡便性のためにこれらを計算する関数がある．この関数は，次元が2×2の構造化配列である，計算された変位，ひずみ，応力を取り，3×3の構造化配列に変換する． 成分は{3,3}の位置に見付かる．

拡張されたひずみを計算する：

平面応力の場合， 方向のひずみ成分は以下で与えられる：

拡張されたひずみを抽出する：
拡張された応力を計算する：

平面応力の設定では， 方向の応力成分は0である．

拡張された応力を抽出する：
軸対称モデル

3D立体が回転軸を持つ場合，軸対称シミュレーションを行うことができる．軸対称モデルは切頭円筒座標系を使う．直交座標系と比較すると，放射軸方向 は 方向の場所を取り，角度方向 は の場所を取り， は対称軸になる．独立変数の順序は固定され変更できない．

625.gif

この図は，破線の 軸を中心に回転する軸を持つ物体を表す．この物体は， 軸を中心に回転する2D領域を利用して角方向 に従うことで，軸対称で表すことができる．材料特性と負荷も 軸について対称である場合は，2Dの軸対称モデルを使うことができる．これは3Dに埋め込まれた2D領域として下で説明する．

2Dの軸対称固体力学方程式の導出は，3Dの直交座標バージョンから開始する．変位はまだ ，， と呼ばれる．ここで は放射方向 の変位， は角方向 の変位， は軸方向 の変位である．円筒座標では微小ひずみは以下のように定義される：

変位ベクトルから円筒座標のひずみテンソルを作成する：

テンソルせん断ひずみ は工学せん断ひずみ のとして定義されることを思い出す：

以下のようになる：

軸対称モデルの簡単な仮定は， 方向には変位がないということである．つまり ということである．また，放射変位 と軸変位 は とは独立である：

次に 方向には変位がないということを使う．簡約された軸対称ひずみテンソルは以下で与えられる：

方向には変位がないということは，ひずみ ，であることも意味する．ひずみ はゼロではなく応力もゼロではない．

切頭円筒座標ひずみテンソルを求める：

上のひずみテンソルを直接計算するために，従属変数指定を直接使ったことに注意する．

応力・ひずみ関係は座標系には関係なく変化しない．線形等方性弾性関係は以下で与えられる：

軸対称せん断ひずみと応力指定を使うと，軸対称応力・ひずみ関係は以下で与えられる：

上で見てきたように， 方向の変位がないということはひずみ ， でもあることを意味する．これより，応力 ， となる．

簡約された軸対称ひずみテンソルは，同等に以下で与えられる：

平面ひずみと平面応力の場合と比べると，軸対称では平衡方程式を円筒座標系に適応させる必要もある．円筒座標系の平衡方程式は以下で与えられる：

直交座標からの円筒平衡方程式の発散項を記号的に計算して，せん断成分の対称性を考慮に入れて上の式を検証する：
直交座標からの簡約された軸対称平衡方程式の発散項を記号的に計算し，せん断成分の対称性を考慮に入れる：

上では最初と3つ目の方程式が平衡方程式に関係している．

例を見てみよう．

幾何モデルを指定する：
メッシュを作成し可視化する：
変数とパラメータを指定する：

あるいは"ModelForm""Axisymmetric"で"ModelForm""RegionSymmetry"を指定することも可能である．

方向の ， 方向の で物体を固定したままで， 方向の負荷の変位を計算する：

境界条件に関する重要な点の一つとして，形状が から始まる場合3Dの回転体は 軸を中心とする穴を持たないため，についての境界条件は物理的な意味をなさないということがある．

変形しない物体の辺の枠の上に変形した物体を可視化する：
軸対称ひずみを計算する：

軸対称の場合，応力の計算には，計算された変位が引数として必要である．

平面ひずみと平面応力の場合と同様に， 方向にひずみ成分と応力成分がある．拡張されたひずみ計算では，関数は次元2×2の構造か配列である，計算されたひずみと応力を取りそれを3×3の構造化配列に変換する． 成分は{2,2}の位置に見付かる．

拡張されたひずみを計算する：

軸対称の設定では， 方向のひずみ成分は以下で与えられる：

拡張された応力を計算する：

軸対称の設定では， 方向の応力成分は以下で与えられる：

軸対称の場合， 方向のひずみと応力をすぐに得る別の方法がある．これを使うためには，従属変数の設定で従属変数 を0に設定する必要がある．

軸対称従属変数指定で を設定する：

設定を とすると，軸対称は3D設定のように扱われる．しかし実際の計算は2Dで行われる．これにはいくつかの利点がある．1つ目は返される変位も であること，2つ目は応力とひずみの両方が 成分であることである．ベクトル値のパラメータは3成分ベクトルとして指定されるようになる．

拡張された方程式として軸対称の場合を解く：

変位の2番目の成分は である．

変形のプロットでは 成分と 成分を使う：
完全な軸対称ひずみを計算する：
完全な軸対称応力を計算する：

時間依存の軸対称方程式は，3Dの場合と全く同じように設定することができる．

時間依存の軸対称固体力学PDEを指定する：

線形弾性の限界

線形弾性の限界を示す数点を以下に挙げる[10, p. 78]．

いつ弾性モデルが適しているかという問題は応力で定義される：

関数 は降伏関数として参照される．材料は で降伏すると言われる．一般的な2つの降伏関数には，トレスカとミーゼスの基準がある．

トレスカ基準は以下で与えられる：

ここで は の主応力（固有値）， は材料パラメータである．

ミーゼス基準は以下で与えられる：

は材料パラメータである．

材料パラメータ と は，線形弾性モデルがもはや当てはまらなくなり，非弾性変形が起こる条件を定義する．

詳細は破壊理論のセクションに記載されている．

SolidMechanicsPDEComponentの方程式形式

このセクションではSolidMechanicsPDEComponentの出力と平衡方程式の関係を説明する．外部荷重が与えられない定常の場合の平衡方程式は次で与えられる．

SolidMechanicsPDEComponent演算子を見ると，以下で与えられる．

しかしSolidMechanicsPDEComponentの出力を見ると，次のようになる．

ここで は4テンソル， は変位ベクトルである．問題を解説していこう．

SolidMechanicsPDEComponentの出力の形式は である：

デフォルトの出力方程式は，方程式で与えられる平衡形式から逸脱する．問題はなぜそうなるかである．その問題に答える前に，システム関数を上書きするひずみ関数または応力関数を使いたいと仮定する．これは"StrainFunction"と"StressFunction"を指定することで行うことができる．これらの関数を指定すると，SolidMechanicsPDEComponentではこれらの関数が非線形であると想定される．簡単にするためにSolidMechanicsStrainとSolidMechanicsStressだけを使う．

指定されたひずみ関数と応力関数を持つSolidMechanicsPDEComponentの出力の形式はである：

ここでDiv形式が想定される．ソルバがこれを見ると，Div演算子はDiffusionPDETerm，ConvectionPDETerm，DerivativePDETermのいずれかと合致する必要がある．利用できる項（必要とされる項）がこれらだけだからである． の関数はないが， の導関数はあるので，DerivativePDETermとしてパースするという選択肢しかない．これはうまくいく上，可能な限りの最も一般的な形式である．しかし，非線形ソルバを発動してしまう可能性があるという欠点もある． の導関数ではなく について解く．それ自体は問題ではないが，非線形ソルバは少なくとも2回関数を評価して収束しなければならないので，この一般的なアプローチでは想定よりもやや遅くなる．結局線形方程式は1ステップで解かれなければならない．方程式が線形の場合，どのようにして線形定式に到達すればよいのだろうか．そのために方程式をDiffusionPDETermと を使う同等の形式に書き換え，以下のような形式を得る．

線形弾性材料の法則を使う：

この材料法則の形式は簡略形である．フォークト表記を完全な四次テンソル表記に変換すると， の係数が得られる．

弾性行列をフォークト表記で設定する：
ヘルプ関数を書いて，フォークト行列を四次テンソルに変換する：
フォークト形式の弾性行列を弾性四次テンソルに変換する：
拡散項のアプローチがSolidMechanicsPDEComponentの出力と同じであることを確認する：
デフォルトの線形弾性材料モデルについてのSolidMechanicsPDEComponentの出力が一般形式および最適化形式についても同じであることを確認する：

固体力学のまとめ

このセクションでは，固体力学方程式を簡潔に定式化する．これによって，さまざまなコンポーネントがどのように採用しているかを理解することが理解しやすくなる．まず平衡方程式から始める：

応力 は非弾性部分と弾性部分からなる：

非弾性応力は内部応力または外部応力等の応力で構成される．SolidMechanicsPDEComponentを使うと，初期応力 が"InitialStress"パラメータを使うように設定することができる．実質的にすべての非弾性応力 は初期応力を介して設定される．

応力の弾性部分は， と弾性ひずみで与えられる構成式の縮約(:)で定義される：

全ひずみ は非弾性ひずみと弾性ひずみからなる：

これは以下のように書き換えられる：

全ひずみ は微小ひずみ等のひずみ測定で与えられる：

非弾性ひずみは次で与えられる：

初期ひずみ は"InitialStrain"パラメータを使って設定することができる．熱ひずみは"ThermalExpansion"と"ThermalStrainTemperature"の設定である．

固体力学をまとめる：

非線形弾性材料モデル - 亜弾性モデル

変形が小さい非線形弾性モデルの特別なクラスは亜弾性材料モデルと言われる．これは変形が大きくなる可能性のある超弾性材料モデルと混同してはならない．亜弾性材料モデルは，変形が小さい弾性材料モデルであるが，構成方程式はもはや線形ではない．ここで提示する亜弾性材料モデルは[11, c.3.3 and c.8.3]に基づいている．亜弾性材料モデルは実際の材料の正確なモデルではないため，それほど広くは使われない．しかし，ここでそれを扱うには理由がある．亜弾性材料モデルははカスタムの材料モデル（ここでは非線形の応力ひずみ関係）の書き方においてよい例になるため，ここで紹介する．亜弾性材料モデルは，線形応力ひずみ関係の限界の外側にある，応力ひずみ測定データをモデル化するために使われる場合がある．

亜弾性材料モデルは動作的現象をモデル化するために使われる．

亜弾性材料は，ひずみが小さい場合でさえ非線形応力ひずみ関係を示すが，同時に完全に可逆的でもある．完全に可逆的とは，適用された荷重が取り除かれると，オブジェクトがもとの状態に戻るというものである．材料が完全に弾性であれば，塑性は関係ない．このモデルでは，ひずみと回転は小さいものとされ，微小ひずみ測定が使われる．応力モデルとして，コーシー応力が使われる．

概念的に言えば，これらのモデルはひずみエネルギー密度関数 を定義し，ひずみについての導関数を求めることによって導出される：

応力を表すためにエネルギー密度の導関数を求めるということについての詳細は，超弾性材料についてのセクションで説明する．

モデル特定のひずみエネルギー密度関数は[11]では与えられていないが，亜弾性構成方程式は以下で与えられている：

パラメータは以下である：

ここで は における一軸応力ひずみ曲線の傾きである．モデルパラメータ ，，は実験データから想定するパラメータである．これについては後で説明する．

この材料モデルを指定する関数を書く．この関数の引数は変数 vars，パラメータ pars，デフォルトのひずみ測定等のデータを含んだデータ data である．以下の材料モデル関数には，モデルの本質が曖昧にならないように，ハードコードされたモデルパラメータがいくつか含まれている．これらのパラメータは，パラメータ pars に置くことができ，関数内でパースされる．

亜弾性材料モデルを実装する関数：

このカスタム材料モデルを試すために，簡単な設定を使う．矩形領域の左側と下側が囲まれており，圧力 が正の 方向に適用される．

747.gif

形状とメッシュを作成する：
亜弾性材料モデルで変数とモデルパラメータを設定する：
境界条件を含むPDEモデルを作成する：
圧力についてPDEモデルを解く：
変形した領域を赤で可視化する：

非線形性が非常に強いため，直接解が得られないこともある．このような場合，非線形ソルバでは収束できないというメッセージが出力される．そうであってもステップを踏めば最終的な解を求めることができる．まず，最終圧力よりもずっと小さい圧力から始めるのである．次のステップで圧力をわずかに上げる．このとき，非線形ソルバが開始するためのInitialSeedingとして，最後に使った低い圧力のときの解を使うのである．後でこの例を示す．

モデルの挙動が非線形応力ひずみ関係であることを検証するために，力と変位の関係のプロットを作成する．プロットを作成するために，指定された圧力のときにある方向における点でどの程度の変位があるかを計算する．

さまざまな圧力入力と更新されるInitialSeedingについて同じモデルを評価するよい方法に，ParametricNDSolveValueを利用するというものがある．

適用された圧力 のパラメトリックモデルを使う：
最大圧力，ステップ数，初期変位を設定する：

後続のサンプリングでは，直前の解を使って，現在の解についての非線形ソルバを開始する：

圧力が上昇する間の変位のサンプリングを監視する：
力と変位の関係を可視化する：

上記の出力は，適用された力（圧力設定による）とオブジェクトの変位の間に非常に非線形な関係があることを示す．

亜弾性材料モデルのデータフィット

残りはモデルパラメータ ，，を求めるだけである．このためには，測定された応力ひずみ曲線の例をロードし，データのフィットを行ってモデルパラメータを推定する．と は測定された応力ひずみデータの線形部分から推定され， は亜弾性モデルにフィットされる．

スチールの応力ひずみデータをインポートする：
引張試験データから適用された荷重の値を取得する：
ひずみ値を取得し，パーセントに変換する：
標本の領域を設定する：
さまざまな荷重での応力を計算する：

データを可視化する．しかし最後のデータ点は標本が破裂する点なので，これは削除する．

破裂のデータ点を削除する：
測定された応力と力が適用されたひずみの関係を可視化する：

ここで，応力ひずみ曲線の線形部分を抽出する．このために，データを切断する．応力ひずみ曲線の線形部分がどこで終りどこで止まるかが明確ではない可能性がある．これは解析者の裁量による．測定のノイズによって，これが困難になる．

応力ひずみ測定の線形部分を抽出し可視化する：
応力ひずみ曲線の線形部分への線形フィットの係数を求める：
応力ひずみ曲線の線形部分にフィットしたデータを作成し可視化する：

この目的は，最適フィットを作成するモデルパラメータを求めることである．試行錯誤で見付け，曲線の線形部分をちょうど超えてフィットするモデルを評価する．この場合，測定された応力ひずみ曲線への最適フィットになる．

フィットしたモデルから と を推定する：
推定データのヤング率を調べる：

次のステップではモデルパラメータ を推定する．

と を亜弾性モデルに入れる：
応力ひずみデータの大きい部分を取り可視化する：
モデルパラメータ を求める：

これでモデルをデータに合わせるためのモデルパラメータすべてが揃った．最後にこれらのパラメータが測定データのモデルを作成することを検証する．このために，さまざまなステップでモデルを評価し，評価されたモデルと応力ひずみ曲線の測定データを視覚的に比較する．

モデルパラメータ を代入する：
モデルを評価し，測定データとともに可視化する：

超弾性

ゴムや発泡体のような材料は大きい変形にさらされてもまだ完全に弾性のままであり得る．つまり荷重が取り除かれると，変形は塑性変形なしで完全に可逆なのである．これらは超弾性材料と言われる．ゴムや発泡体以外にも，ゴムのような形態の生物組織やポリマーが超弾性材料のカテゴリに入る．亜弾性材料とは対照的に，超弾性材料は大きい変形や回転を受ける．

詳細は超弾性および超弾性材料の法則に関するノートブックを参照されたい．

破壊理論

材料のタイプによって異なる破壊理論が利用できる．例えば延性材料の場合はミーゼスおよびトレスカの破壊理論が，脆性材料の場合はモール・クーロン理論および修正モール理論等が使える．

ほとんどの破壊理論は主応力と材料特性を比較する．延性材料では主応力と降伏強度を比較する．脆性材料では主応力と極限強度または破壊開始点を比較する．主応力は，応力値が，単軸引張試験から導かれた降伏強度および極限強度のような材料データと比較可能になるように計算される．境界荷重：張力のセクションには，PrincipalEigenvalueを使って非軸の荷重円筒における応力を求める例が記載されている．

破壊理論はさまざまな観察挙動を説明する必要がある．例えば，脆性材料は引張強度より大きい圧縮強度を持つ傾向がある．また，静水圧応力は延性材料では降伏を引き起こさない．静水圧応力は，形状の変化を引き起こす偏差応力とは異なり，体積の変化を引き起こす応力である．

次ではさまざまな破壊理論を説明する[12]．

ランキン

ランキンの破壊理論では，主応力の最大または最小値の絶対値が，延性材料では降伏強度に達し，脆性材料では極限強度に達したら物体は破壊されたと考える．ランキンの破壊理論は最大主応力説とも呼ばれる．ランキンの理論は簡単なものであるが，特に延性材料の場合は特に正確な理論とは言えない．延性材料については，ランキンの理論は静水圧応力が延性材料では降伏を引き起こさないということを配慮しない．したがってランキンの破壊理論を使うとしたら脆性材料に対してである．

トレスカ

トレスカの破壊理論は，最大せん断応力説とも言われる．

最大せん断応力が単軸引張試験の降伏時にせん断応力と等しいくなると，降伏が起こる[12]．

単軸試験において は適用された応力に等しい．応力が適用された単軸方向だからである．それに応じて と の応力は0である．したがって ，つまり以下が成り立つ．

トレスカの安全率は以下で計算できる：

安全率が1未満であると問題である．通常物体はどのような使い方をしても安全率が常に1を上回るように作られる．降伏応力よりも高い応力を見付けたら，注意しなければならない．

トレスカ理論は，静水圧応力は延性材料の降伏には関係ないということに一致している．

トレスカ理論の方がミーゼスの破壊理論よりも適用しやすく，保守的である．この2つの理論の最大差を計算すると15.5%になる．ミーゼス理論の方が実験データによりよく一致する．

トレスカ理論はモール・クーロン理論の特殊形である．

以下の例では，応力テンソルに対する安全率を計算する．

応力テンソルを作成する：
主固有値を求める：
最大せん断応力 を求める：
降伏強度 の材料の安全率は1を上回っている：

概要の安全率セクションの例もご覧いただきたい．

ミーゼス

ミーゼスの破壊理論は最大ひずみエネルギー説とも言われる．

最大ひずみエネルギーが単軸引張試験の降伏時のひずみエネルギーと等しくなると降伏が起こる．

トレスカ理論の方がミーゼスの破壊理論よりも適用しやすく，保守的である．この2つの理論の最大差を計算すると15.5%になる．ミーゼス理論の方が実験データによりよく一致する．

モール・クーロン

モール・クーロン理論は脆性材料に使われる．トレスカ理論はモール・クーロン理論の特殊形である．

修正モール

修正モール理論は，モール理論を拡張して実験データによりよくフィットするようにしたものである．

複数材料

複合材とは複数材料で作られた物体である．複数材料の設定を示すために，簡単な2つの棒が両端で制約を受けると同時に表面力を受ける．

厳密に必要な訳ではないが，幾何学モデルで内部境界による材料境界を示した方がよい．

複合梁の境界メッシュを生成する：

梁の中心に分け目があることに注目する．分け目の左側に一つの材料，右側に別の材料がある．

内部境界を持つ形状を作成するための，手動作業が少なくより一般的な方法は，部分領域を持つ要素メッシュのセクションに記載されている．3Dの場合はOpenCascadeLinkを使うとよい．

OpenCascadeLinkを使って材料領域を持つ形状を構築する：
メッシュを生成し可視化する：

内部の材料境界は完全なメッシュの作成中維持されている．2つの材料に働く力の影響をよりよく示すために，2つの極端な要素を選ぶ．左側はチタン，右側は鉛でできているものとする．

2つの材料領域のパラメータを設定する：

使用される材料がEntityとして利用できない場合，材料パラメータはIf文のような条件文として与えることができる．

複数領域に対する別の材料パラメータ設定：
左と右で固定制約され，z=0.1である上部に負の 方向の力があるPDEモデルを設定する：
方程式を解く：
総変位を可視化する：

最大の総変位は材料境界の右にある．鉛はチタンよりもずっと柔らかいのでこの移動は想定通りである．

これに続く解析は，単一材料の場合と全く同じように行われる．

より複雑な形状の場合，ElementMarkerを使って部分領域，境界部分，境界ノードを指定した方が簡単なことがある．ElementMarkerの生成と使用方法は，要素メッシュの生成チュートリアルで説明してある．このセクションには，メッシュ上にElementMarkerを含む関数の評価方法についての記述もある．

減衰

ここまで見てきた時間依存の固体力学モデルは非減衰系であった．つまり，振動している物体は永遠に振動を続けるということである．これは実生活では起こらない．例えば，音叉を叩くとその振幅は時間とともに減衰していき，弱まっていく．減衰するにつれてエネルギーは散逸していく．減衰にはいくつかモデルがある．ここでは，振動している物体が振動を弱める粘性媒体に囲まれているため，エネルギーを散逸する粘性減衰を考える．

減衰のモデル化には時間依存PDEモデルが必要である．時間積分の入門は時間依存解析セクションをご一読いただきたい．

レイリー減衰

線形弾性材料の固体力学についての一般化された運動方程式の行列形式は以下で与えられる：

レイリー減衰は減衰行列 を構築するメソッドを導入する．項 は粘性減衰を表し，物体の速度に影響を及ぼす．レイリー減衰は， は質量行列 と剛性行列 の線形結合であると仮定する．

質量減衰パラメータ の単位はであり，剛性減衰パラメータ の単位はである．特筆すべきは，この線形結合仮定には十分な物理的根拠はないが，広く使われているほどかなりうまくいくという点である．レイリー減衰は現象論的モデルであり物理モデルではないので，パラメータ および の値を特定することは難しい場合がある．以下でこの点について説明する．

レイリー減衰パラメータを含む平衡方程式の一般形式は以下で与えられる：

そして方程式を一般化された波動方程式系に変換する．さらに二階時間微分に変換することで一階時間微分も持つ．

初歩的な例として，平面応力モデル形式で長さ メートル，高さ ， 厚さ の矩形領域を設定する．梁は左端で固定され，右側では下向きの荷重がかかっている．この例では，下向きの力は一定の時間非依存の値である．これは実質的に，荷重は時間積分の最初から単位ステップ関数として適用されることを意味する．非減衰モデルと減衰モデルを作成する．減衰モデルでは質量パラメータと剛性パラメータが定義されている．これらの値は文献のものであるか，共振周波数のために計算されたものであるかのどちらかである．これについては以下で説明する．

領域と時間依存変数を設定する：
一定圧力を含む非減衰の平面応力モデルパラメータを作成する：
質量減衰 ，剛性減衰 という試行値で減衰モデルパラメータを作成する：
シミュレーションの終了時間を設定する：

次に，モデルパラメータ一式を取り変位を計算する補助関数を設定する．これによって，確実にシミュレーションしたい両方の場合に同じ設定が使われるようにする．補助関数にはPDEモデル，境界条件と初期条件，領域，時間範囲が含まれる．

時間積分の間の進行状況を監視するために，EvaluationMonitorが加えられている．また，解の特徴がステップオーバーされないようにMaxStepSizeも加えられている．

PDEモデルのモデルパラメータ依存シミュレーション関数を設定する：
非減衰方程式を解き，時間とメモリ使用量を測定する：
減衰方程式を解く：

完全にするために，定常解も計算する．

定常変位を計算する：
方向のにおける定常変位を評価する：
方向のにおけるクエリポイントの変位に従って揺れている梁の減衰を可視化する：

減衰パラメータを指定すると，減衰PDEモデルはクエリポイントの変位の振幅の減衰を示すが，非減衰モデルでは最小減衰しか示さない．この最小減衰は数値誤差によるものであり，MaxStepSizeを小さく設定することで減少させることができる．どちらも定常解の値付近で振動する．

さまざまな時間における減衰した梁の変形を示す枠を作成する：

レイリー減衰パラメータ と のよい値を得るのは難しい．さまざまな自然の繰返し周波数 の減衰比 は次で計算できる：

以下のプロットは，減衰比 に対する繰返し周波数 をプロットすることでさまざまな減衰タイプを可視化している．

805.gif

のとき，で剛性減衰となる．剛性減衰は線形関係である．プロットから剛性減衰が分かり，剛性減衰パラメータ が高周波数を減衰させる．のとき， で質量減衰となり，質量減衰パラメータ は低周波数を減衰させる．質量減衰は「水中」減衰のようなものである．レイリー減衰は質量減衰と剛性減衰の線形結合である．

前の方程式を操作すると，直接 と を定義することもできるようになる．これは2つの共振周波数 と をで指定し，対応する2つの減衰係数 と を指定することで行うことができる．これにより以下の減衰パラメータが得られる：

制約条件を加える：

ほとんどの工業用金属材料では < 5%という値が想定でき，その多くが2%を下回る．ゴムのようなエラストマー材料では，値は大きくなる．ここで の値を知っているものと仮定すると，次に問題となるのは周波数 と の求め方である．このために，固有モード解析のセクションで示したように，非減衰系の固有値解析を行う．これにより固有振動数とそれに対応するモードが分かる．その後，減衰させたい運動に最も関係のある2つの最低の周波数モードを選ぶ．関連する固有振動数が周波数 と である．

PDEモデルについて最初の6個の固有値とモードを計算する：
固有振動数を計算する：
その振動数における変形モードを可視化する：

下向きの運動を減衰させたいので，最初の2つのモードを選ぶ．3つ目のモードは 方向なので関係がない．高いモードは省略する．ポイントは，最低の最も関係のあるモードを選ぶことである．

最初と2つ目の周波数に対する減衰係数はそれぞれ， ()， ()と仮定する．減衰係数の値が未知の場合，後ほど測定手順を紹介する．

質量減衰の値を計算する：
剛性減衰の値を計算する：
モード減衰モデルパラメータを設定する：
減衰方程式を解く：
非減衰系，減衰系，固有モードベースの減衰系それぞれの減衰挙動を可視化する：

モードベースの減衰は，周波数 における1%の質量減衰の減衰係数 と周波数 における2%の剛性減衰の減衰係数 の結合である．

減衰係数 の値が利用できない場合，衝撃試験で測定することができる．このためには，次のように定義された対数減衰率を使う．

ここで は周期の数，はピークの振幅，は 先の 周期のピークの振幅である．

1つのモードだけを刺激することは重要であるが難しい．特定の一つのモードの減衰係数は以下で計算することができる：

次の例では手順を説明する．このために系の特定の一つのモードを刺激する．変位の振幅と時間の関係を示すダイアグラムはシステムのリングダウンと言われることもあり，これが記録される．手順を示すために，まず架空のリングダウンデータセットを生成する．通常これは実験で得られる．

シミュレーションが行われたリングダウンデータを生成しプロットする：

目的はシステムのリングダウンに対する減衰係数を求めることである．まず，減衰解の2つのピークにおける時間を求める：

0.3に近いピークの時間を求める：
1に近いピークの時間を求める：
ピークを可視化する：
対数減衰率 を計算する：
減衰係数 を計算する：

この手順を使うと，1つの減衰係数が見付かり，2つ目の減衰係数を見付けるために実験が繰り返される．

計算された減衰係数が正しいことを検証するために，シミュレーションが行われたリングダウンのエンベロープをプロットする．これができるのは，架空の実験データで固有振動数 も指定したからである．

非減衰振動のエンベロープを破線で可視化する：

を低く指定すると，振幅と応力が大きくなるため，減衰係数を低く設定する方が保守的である．疑わしいならば，低い減衰係数を指定するべきである．減衰係数は周波数の関数である．このことと2つのモードしか使わないということがレイリー減衰の短所となる．周波数の範囲が大きいと，必要な減衰にほとんど適合しないのである．周波数範囲が大きすぎる場合，中間の周波数領域の減衰がほとんどなく，低周波数領域と高周波数領域での減衰が大きすぎるということが起こりがちである．

実際のデータを入手し利用することは非常に複雑になるので，ここでは架空のリングダウンデータ集合を使った．

有限要素プログラミングチュートリアルには，レイリー減衰をシステム行列レベルでどのように実装できるかを示したセクションがある．

固体力学境界条件

固体力学での境界条件は2つのカテゴリのどちらかに属する．一つは境界制約，も一つは境界荷重である．制約は物体の変位を一方向または複数方向に制限，強制する．例として壁に固定された物体がある．固定点では物体の変位は不可能である．このタイプの境界条件はディリクレ条件で実現される．したがってその名前には「条件」という言葉が含まれる．微分方程式が解けるようにするためには，通常少なくとも一つの条件タイプの境界条件を指定しなければならない．もう一方で境界荷重も表面力と呼ばれる．これは表面に作用する力，圧力，運動である．例として，本箱での本の重みのような，表面に働く外部力がある．このような境界荷重はノイマン値境界条件で実現され，その名前には「値」という言葉が含まれる．境界のある部分で境界条件が指定されていない場合は，その部分への表面力はなく，物体は自由に動く．

さまざまな境界条件を説明するために，デフォルトの単位を持つ固体力学PDE成分を設定する．

固体力学PDE成分を設定する：

境界荷重では，境界条件の述語としてElementMarkerを使った方が便利である．制約条件の場合は幾何述語，あるいはSelectPointMarkerFromBoundaryMarkerでPointMarkers → BoundaryDeducedを指定して生成したメッシュを使う．

変位制約

さまざまな制約を説明するために，簡単な梁をモデル化する．梁は左側で固定されている．

粗いメッシュを生成する：

強制変位

梁を左で固定する．
右側の変位を強制する．
PDEを解く．
変形した梁を可視化する．

右側の面が負の 方向に下側に1単位移動した． 方向と 方向の変位は指定された通り0である．

強制変位から方向を除外するためには，Noneを使うことができる．例えば，これを使うと負の 方向の変位を強制する一方で，物体は他の方向には自由に動くことができる．

負の 方向の変位を強制し他の方向は強制しない．
PDEを解く．
変形した梁を可視化する．

負の 方向だけに1単位の強制変位を指定したとき，右側の 方向にも変位があることに注意する．

ローラー拘束

ローラー拘束は，制約が適用さえる面に垂直な物体の変位を制約するために使われる．つまり，物体は面に垂直な方向には動けないが，面は他の方向に自由に動けるのである．これは面がローラー上に載っているかのようなものである．面もまた，負の垂直方向に制限されている．言い換えると，面は垂直方向に動くことはできないが，法線に垂直な接線方向には動くことができる．

一般的にローラー拘束は法線に沿った従属変数を以下のように制約する：

ここで ，， は従属変数であり， はさまざまな方向の法線である．現在，トップレベルでは，法線 が軸方向で法線成分のうちの2つがゼロである場合のみローラー拘束がサポートされる．その場合，ローラー条件はDirichletConditionで実現することができる．より一般的なアプローチでは，低レベルの有限要素関数を使う必要がある．

以下の画像は，ローラー拘束の働きを示している．梁は黒で示したように左下の端で固定されている．上部では，青い矢印で示したように下向きの圧力がかかっている．赤でハイライトした領域では，ローラー拘束が作用している．この領域では梁は表面に垂直に動くことはできない．

862.gif

左下の端を固定する：
上部に下向きの圧力がかかっている：
方向における変形の部分制約：
PDEを解く：
変形した梁を可視化する：

必要な制約条件の数

DirichletConditionタイプの境界条件を全く適用しないで方程式系を解くと，NDSolveはDirichletConditionまたはロビンタイプのNeumannValueが指定されていないという警告メッセージを発する．DirichletConditionを指定しないことで方程式系が特異になり，解が見付かったとしても定数までに過ぎないからである．この挙動は固体力学の応用とは無関係であるが，一般に真実でありDirichletConditionおよびNeumannValueのリファレンスページで示されている．物体を十分に制約するためにはいくつの制約条件が必要なのだろうか．目標は剛体運動を排除することである．その理由は，物体と荷重が並進または回転する場合，物体とその荷重の全体的な性能に剛体運動は関係ないからである．3Dではどの物体も6自由度である．物体は3方向に並進し3つの軸すべてについて回転することができる．2Dでは物体は2方向に並進し平面で回転することができるので，2Dでは3自由度である．

解くことのできる完全に制約されたPDEモデルを作成するためには，少なくとも剛体運動を含む自由度を制約する必要がある．Wolfram言語で使われる連続体有限要素は回転自由度を提供しないので，それを直接設定することはできない．この場合，並進自由度を制約するには，回転運動を制約するようなものにする必要がある．

2Dの例を見てみよう：

864.gif

上では，のノードは 方向と 方向で制約されており， 変位および 変位における0 DirichletConditionに当たる．回転制約では，右下に別の制約条件が加えられる．ここでは 方向に一つの制約条件がある．つまり，物体は 方向には自由に動けるのである．この2つ目の制約条件は，2つの制約条件の間の方向に垂直な運動を禁止する．

例えば，他を選ぶこともできる：

872.gif

3Dで類似の方法を使う．ここでは ， ， の並進を制約する制約条件， 方向と 方向等の2つの方向を制約する2つ目の制約条件， 方向等，一つの方向を制約する3つ目の制約条件を加える．

879.gif

最小の制約条件の実際の選択が任意であることを示すために，内側に非対称に置かれた穴を持つ矩形の形状を考える．

非対称の形状を生成し可視化する：

この例は平面応力モデルの形式を使う．外側の境界にはすべて圧力がかかっている．

平面応力モデルを設定する：

このモデルを2回解く．どちらも異なる位置で最少数の制約条件を課す．

{-1,1}でuとvに，{1,-1}でvに制約条件を課して変位について解く：
変位を可視化する：

次のステップで再びこの問題を解く．今度は別の位置に制約条件を課す．

{-1,1}でuとvに，{-1,1}でuに制約条件を課して，変位について解く：
変位を比較する：

明らかに変位は異なる．どちらの場合も{0,0}におけるノードにはuとvへの制約条件がある．最初の場合（左側），vに対する2つ目の制約条件は右下で設定される．2回目の場合，uに対する2つ目の制約条件は左上に設定される．異なる制約条件を選ぶことで変位は異なるが，ひずみと応力は異ならない．

両方の変位に対するひずみと応力を計算する：
両方の場合の応力を可視化する：

どちらの場合もミーゼス応力は同じであることが分かる．

境界荷重

境界荷重を考えるとき，表面力という言葉がよく使われる．単位[]の表面力 は以下で定義されるベクトルである：

言い換えると，表面力は応力のように面積あたりの力である．成分形式で書くと以下になる．

応力テンソル が対称の場合は以下のようにも言える：

体積力がない，物体が平衡な場合は，境界に沿った表面力上での積分はゼロになる：

Gaußの定理はテンソルについても当てはまる：

左辺の境界積分はNeumannValueの関数ページまたは「偏微分方程式と境界条件」のノイマン境界値の導出で説明した通りNeumannValueである．右辺は平衡方程式 である．したがってNeumannValueを使って表面力ベクトル成分を指定することができる．

荷重は表面に作用する圧力または力である．このどちらもSolidBoundaryLoadValueで指定することができる．表面に作用する力の場合，その力 は，それが作用している表面積 を考慮することによって，自動的に圧力 に変換される．つまり，表面力は自動的に表面圧力に変換されるのである．

圧力または力は常に面積あたりの力の単位で，また二次元で与える必要がある．二次元の場合，圧力は"Thickness"モデルパラメータの値で内部的に乗算され，長さ単位あたりの力が得られる．長さ単位がメートルの場合，圧力単位から単位に変換される．

SolidBoundaryLoadValueを評価するとNeumannValueのリストになる．

圧力境界荷重を設定する：
力境界荷重を設定する：

力成分は，境界荷重がアクティブとなっている面積で分割される．圧力または力の成分はいずれも0になり得る．しかし，SolidDisplacementConditionでは可能なNoneの指定はない．

境界荷重のさまざまな概念を示すために，標準的な荷重形式を表示する．標準的な荷重タイプを以下の図で示し，説明する．

895.gif

さまざまな荷重を例示するために，簡単な円筒を作成する．

幾何モデルを設定してメッシュを生成する．
面に色とラベルを付けてメッシュを可視化する．

以下の例ではすべて，棒は左で壁に固定される．

棒を左で固定する：
さまざまな場面で作用する圧力を指定する：

境界荷重：圧縮

896.gif

軸方向圧縮
円筒の右面で負の 方向の圧縮圧力を指定する：
PDEを解く：
変形した物体を可視化する：
非軸方向圧縮

境界荷重は軸方向である必要はない．例えば，円筒の物体に垂直に作用する圧力場を考える．

圧力場を可視化する：
明示的な方向ベクトルを設定する：
PDEを解く：
変形した物体を可視化する：

もっと一般的なアプローチとして，BoundaryUnitNormal式を使って必要な法線成分を抽出するというものがある．BoundaryUnitNormalを計算することは，方向ベクトルを明示的に指定するよりも計算的に高価である．方向ベクトルがすぐに利用できる場合は，そちらを利用した方がよい．他の場合，例えば幾何モデルが複雑なときはBoundaryUnitNormalの方がよい．

BoundaryUnitNormalを使って方向ベクトルを設定する：

境界荷重：引張

898.gif

円筒の右の面に正の 方向の引張圧力を指定する：
PDEを解く：
変形した物体を可視化する：

この検証の例もご参照いただきたい．

簡単なテストとして，ひずみと応力を計算することもできる．それを使うとその領域上の力は 方向の応力である：．回復可能な面圧力 を指定した．

ひずみと応力を計算する：
ある点における応力をすべて評価する：

応力は，適用された表面圧力に対応する．

円筒の十分内部の点における圧力は指定された圧力であることを検証する：
比較としてミーゼス応力を可視化する：

固定された端である左に応力集中があるが，これは想定通りである．

背面の応力集中を可視化する：

この応力集中はポアソン効果と呼ばれる．これは左で円筒を固定する境界条件を変更することで防ぐことができる．それについてここで説明する．その前に， 方向の応力が，指定された境界圧力の1%より大きくなる体積を見てみる．

指定された圧力より1%大きい 応力のある体積を可視化する：

上で述べたように，左境界の応力集中は境界条件によって引き起こされる．これは想定内のことで容認できる．境界条件を設定するのに少し異なる方法がある．しかしそれは壁境界をモデル化する少し異なる方法でもある．このためには， 方向を制約し， 方向と 方向も制約される1点を除いて，他の方向はすべて自由に動くようにする．

別の壁制約：
変位を計算する：
変形した物体を可視化する：
残りの相違は数値的ノイズであるミーゼス応力を計算し，可視化する：

全体的な応力分布は，指定された表面圧力に近くなった．残りの差は数値的雑音である．このアプローチは，で固定条件を課すこの例のように，使用できる物体に対称性がある場合にうまくいく．非対称の立体では，その条件をどこに置くかが常に明確ではない可能性がある．

非軸方向引張

前の例では， 軸と平行な円筒と，同じ軸に平行な境界荷重を使った．このため， 応力成分は表面圧力を反映した．物体に軸と並行の荷重がかかっていない場合，応力（またはひずみ）成分は軸と平行にならない．

簡単にするために，円筒を回転させるために任意の回転行列を使い，境界荷重が円筒の右端に垂直になるように回転させる．

任意の回転行列を作成する：
円筒を回転させる：
適用した圧力を青で描いて回転された円筒を可視化する：
メッシュを生成する：
回転した荷重を設定する：
PDEを解く：
ひずみと応力を計算する：
クエリポイントを回転させる：
回転したクエリポイントで応力を評価する：

これで 垂直応力は表面荷重とは同じではなくなった．荷重がもはや軸と並行ではないので，これは想定通りである．ミーゼス応力はさまざまな垂直応力とせん断応力を組み合せるので，これはまだ表面圧力を与えるという点に注意する．

ミーゼス応力を計算する：

円筒と荷重が軸と並行であったかのように，まだ最大垂直応力を求めることができる．これは，主応力値と方向ベクトルのセクションで示した通りPrincipalEigensystemを計算することで行える．

回転したクエリポイントで主固有値を計算し評価する：

円筒と荷重が軸と並んでいたら，最大の固有値である第1主応力は応力値に対応することが分かる．

境界荷重：せん断

913.gif

円筒の右面で正の 方向にせん断圧力を指定する：
PDEを解く：
変形した物体を可視化する：

この検証例もご参照いただきたい．

境界荷重：ねじり

最終的に境界荷重はすべて表面に作用する圧力に変換される必要がある．境界力を圧力に変換することは，表面力が作用している面積を計算し，計算された面積で力を割ることである．ねじりにも同じ過程が必要である．

この例では棒は左で固定され，のトルクが右端に適用される．

916.gif

トルク は表面圧力に変換する必要がある．以下から始める．

ここで はせん断応力（圧力）， は半径，は面積の2次モーメントである．並べ替えると以下になる．

右に作用する境界荷重を設定する：
PDEを解く：
変形プロットを表示する：

このシャフトは右端で放射状に大きくなっているように見える．これは物理的に正しくない．この場合，2つの問題によって引き起こされている．一つはメッシュの変形プロットによって見え方が誇張されているという問題である．要素メッシュ変形プロットは，変形の効果をよりよく可視化するために固定量の拡大を計算する．したがって，解に小さい誤った半径成長がある一方，それは変形プロットによってさらに誇張されるのである．小さい半径方向の変形があるということは，この場合使われる線形理論によるものである．

Automaticスケール因子を使って計算する：

線形弾性モデルとそれに関連したひずみ量は小さいひずみと小さい回転を仮定する．せん断角度は小さいままであることに気を付けなければならない．回転角度 は三角関数の簡約化 およびが成り立つほど小さくなければならない．この近似が成り立たない場合は，大きい変形を考えなければならない．次のグラフィックスを考える．

927.gif

で回転した点は最終的に位置になる．ここで小さい角の仮定をすると，点はになる．これは要素メッシュ変形プロットで選ばれるものである．

基準点を{0.1,0,radius}で指定する：
基準点の変位位置を求める：
基準点と変位点の間のせん断角度を計算する：

せん断角度は，トルクのとき1度よりわずかに小さい．

境界上の点{0,0.01}を青で，同じ変位点を赤で可視化する：
小さい角の仮定が満足されることを検証する：

せん断角度は十分小さいので，小さいひずみと小さい回転の仮定は有効である．小さいせん断角度でも，半径方向にわずかな成長があり，それが変形プロットに取り上げられ誇張される．

ひずみと応力を計算する：
ミーゼス応力を可視化する：

空洞の梁の方が効率的にねじり荷重をもつ．梁の中央部はねじり荷重のわずかな部分に抵抗するだけなので，その材料を取り除くことはねじり荷重の下の梁の性能にたいした影響を及ぼさないからである．

最大のミーゼス応力を計算する：
チタンの降伏強度を抽出する：
安全率：

この大きさのトルクでは，チタンの降伏強度よりわずかに大きい．

検証例もご参照いただきたい．

大きいねじり

このセクションではやや学術的な問題を見てみる．降伏強度を無視して，大きい回転となるさらに大きいモーメントを適用する．ひずみ測定のセクションから，無限小のひずみ測定は大きい回転には向いていないことが分かっている．

右側に作用する境界荷重を誇張したモーメントで設定し，方程式を解く：
基準点の変位位置を求める：
基準点と変位点の間のせん断角度を計算する：

小さい角に必要な条件はもはや満足されない．

境界上の点{0,0.01}を青で，変位点を赤で可視化する：

線形弾性材料モデルを使っているため，デフォルトのひずみ測定は無限小のひずみ測定である．同じPDEモデルを再び計算することができるが，微小ひずみの代りに，小さい回転制限を要求する非線形グリーン・ラグランジュひずみを使う．

固体力学PDE成分を設定する：
グリーン・ラグランジュひずみ測定で非線形PDEを解く：
基準点の変位位置を求める：
基準点と変位点の間のせん断角度を計算する：
境界上の点{0,0.01}を青で，変位点を緑で可視化する：

次に変形をプロットする．

変形プロットを表示する：

ここで，要素メッシュの変形が実際には圧縮を表していることが分かる．要素メッシュ変形プロットは変形が見えるように強調しようとする．この場合，変形プロットで示される効果は純粋に可視化によるものであり，微小ひずみの小さい回転の制限によるものではない．

Automaticスケール因子を使って計算する：

変形は小さいので，変形をスケール因子1でプロットすることができる．これは"ScalingFactor"1または"ScalingFactor"Noneで指定する．

"ScalingFactor"を1として変形プロットを表示する：

ねじりの量は材料の降伏強度を大きく超えていることを強調したい．大きいねじりは，微小ひずみとグリーン・ラグランジュひずみの違いを強調するために使われた．

境界荷重：曲げ

この種の境界荷重の例は，ローラー拘束のセクションに記載されている．

対称条件

次に深海のバチスカーフに付いているような球形シェルへの一様圧縮荷重を調べる．このモデルの立体は，境界荷重がすべての面に垂直に作用する中空の球である．換言すると，表面に作用するSolidBoundaryLoadValueの値しかないということである．このことは，NeumannValueとSolidBoundaryLoadValueを評価して得られるものが一意にPDEを解くのに十分ではない．なぜそうなのかについての詳細はNeumannValueページを参照されたい．解法は2D平面ひずみの例のように領域の対称性を使う幾何学的制約を使うことである．中空の球全体をモデル化する代りに，その1/8だけを使う．これによって評価すると必要なDirichletConditionになるSolidDisplacementConditionを使ってモデル化される対称制約条件を導入することができる．

中空の球の1/8を作成する：
境界マーカーでメッシュを可視化する：
外面に垂直の圧力を設定する：

対称境界条件のそれぞれの切断面が切断の平面上の物体を修正し，同時に別の方向の運動を可能にする．

対称制約条件を作成する：
PDEを解く：
もとの領域をエッジフレームとして，変形した領域を可視化する：

付録

このセクションにはいろいろな役立つ情報があり，固体力学PDEモデルがどのように動作するのかをよりよく理解するのに役立つユーティリティ関数が含まれている．また，境界条件の述語指定についての別のワークフローや，応力特異性について考えられる問題のセクションもある．

モデルパラメータ

与えられるモデルパラメータは，追加情報で強化される．例えば，「Material」を指定すると，その材料から適切な"YoungModulus"と"PoissonRatio"が抽出され，モデルパラメータに保存される．強化されたこのようなモデルパラメータを調べることができる．

モデルパラメータを設定する：
処理されたモデルパラメータを調べる：

境界条件の述語

有限要素解析を行うためには，幾何モデルの境界メッシュ表現はメッシュでなければならない．境界条件の設定を簡単にするために，境界条件が指定される前に完全メッシュが生成される．実質的に，境界条件には2つのタイプがある．一つは表面に作用する境界条件でNeumannValue型のものである．「Value（値）」で終る固体力学境界条件の名前はすべてこの型である．二つ目の境界条件はDirichletCondition型のものであり，メッシュの表面ノードに作用する．「Condition（条件）」で終る固体力学境界条件の名前はすべてこの型である．

境界条件がアクティブな位置を述語と言う．このような述語を指定する形式はいくつもある．その一つとして，のような方程式として述語を指定することができる．別の方法には，ElementMarkerで境界表面と境界点を選ぶというものもある．ElementMarkerの生成と使い方はElementMeshの生成チュートリアルで説明してある．NeumannValueの境界値に対して使われる要素マーカーと，DirichletConditionの境界条件は異なる．境界値は境界要素のマーカーを使う一方，条件は点要素のマーカーを使う．このことで柔軟性が生まれる．複雑な幾何モデルの場合，境界面の要素に基づいて境界値と境界条件の両方を設定した方がよい可能性がある．こうすることで点マーカーと境界マーカーが関連付けられるが，通常は必要だとは限らない．これにはオプション"PointMarkers""BoundaryDeduced"を指定する必要がある．

点のマーカーが境界マーカーから推測される境界要素メッシュからメッシュを生成する：

固体力学は，荷重と制約下の物体の変形の計算に関するものである．荷重は物体の表面に適用される力または圧力である．荷重が物体に適用されるためには，その物体は壁に固定される等何らかの制約を受けていなければならない．そうでないと物体は荷重に対する抵抗をもたらさないからである．荷重と制約は境界条件を指定することで設定することができる．さまざまな境界条件が利用でき，これらについては境界条件のセクションで詳しく述べる．この例では，境界面の荷重および壁とねじによる制約で十分である．

これらの条件はDirichletConditionおよびNeumannValueでモデル化することができる．境界条件が表面に適用される場所を求める簡単な方法は，幾何モデルの境界に目を通すことである．

本棚の物体の境界に目を通す：

表面荷重は本棚の棚受けの上面に下向きに加えられる．上+を押して表面をブラウズすることによって，上面に対応する境界マーカー識別番号(ID)が見付かる．

境界荷重面の位置を保存する：

見付かったElementMarker IDを使って，表面荷重値を指定することができる．

表面荷重値を設定する：

次のステップでは，壁に付けられた棚受けをモデル化する制約を指定する．棚受けを壁に固定するためには，2つの制約が必要となる．一つは棚受けの背が壁に突き刺さらないことである．つまりその境界は 方向には全く動けないという制約が付く．棚受けが壁を突き抜けることができないのは明らかなので，負の 方向の制約は保証される．

最初の制約は棚受けが壁にねじ止めされているということである．つまり，ねじと棚受けの接触点では，運動は全く不可能なのである．2つ目の制約は，棚受けは壁に突き刺さらないので，負の 方向の制約は保証される．2つのねじが棚受けを壁に押し付けているので，適切なアプローチは正の 方向の運動も制限することである．これは一般的な簡約化である．この簡約化が使われなければ，接触問題が発生する．棚受けは実際に正の 方向では壁から離れているからである．

背面要素マーカー番号から導出された背面制約点の位置を選び保存する：
棚受けの背面を構成するノードの位置を可視化する：
ローラー拘束として制約条件を設定する：

境界条件指定の詳細は解析タイプによって異なるため，それぞれ適したセクションで示されている．

次にねじの点マーカーを構築する．このために，際立った穴を構成する関係のある境界要素マーカーを選ぶ．背面からのノードは排除する．

境界要素マーカー番号から導出したねじの制約点の位置を選んで保存し，背面は排除する：
棚受けのねじ穴を構成するノードの位置を可視化する：
ねじをモデル化するために制約条件を設定する：

境界条件指定の詳細は解析タイプによって異なるため，それぞれ適したセクションで示されている．

応力特異性

このセクションでは，内向きの角では応力が特異であることを説明する．[6, p. 21]に従うと，簡単な棚受けは上部が固定されており，下の青い矢印で示した通り，力は右面で下向きに作用する．赤い線は応力特異性が想定される部分を示している．材料パラメータがミリメートルでスケールされるよう長さの尺度はミリメートルなので，厳密な次元や材料パラメータは関係ない．領域が作成され材料パラメータが設定されると，一連の細かいメッシュを使って，応力特異性が内向きの角でどのように発達するかを示す．

938.gif

領域，ミリメートルにスケールされる材料パラメータ，PDEを指定する：
領域からデフォルトのメッシュを生成する：

メッシュの密度をよく理解するためには，生成されたメッシュのワイヤフレームを見るとよい．

ワイヤフレームを可視化する：

ミーゼス応力と総変位の計算を隠す小さいヘルパー関数がある．この関数はヘルパー関数の付録セクションで定義される．左のプロットは最大応力の場所を示したミーゼス応力を可視化する．同時に右のプロットは総変形を可視化する．応力は特異であり，変形はシミュレーションを通して一定である．

mesh1 についてのミーゼス応力と総変形を可視化する．

応力は内向きの角で最大である．応力をよりよく解決するために，次のステップではメッシュを調整する．最初の調整アプローチではメッシュ全体を調整する．これにより要素数が増加する．

メッシュを調整する．
調整したメッシュのワイヤフレームを可視化する．
mesh2 についてのミーゼス応力と総変形を可視化する．

計算されたミーゼス応力の値がどのくらい増加したかに注目のこと．同時に棚受けの総変形は同じままである．さらに調整すると応力値の収束につながるという結論に達するかもしれない．一つの方法として，特に関心のある部分でメッシュを調整するというものがある．

内向きの角でメッシュを調整する．
調整したメッシュのワイヤフレームを可視化する．
mesh3 についてのミーゼス応力と総変形を可視化する．

応力値はさらに向上し，応力が内向きの角で最大になることがより明らかになった．これらの応力を避けるためには，通常角にフィレットを入れる．

内向きの角で半径のフィレットを作成する．
デフォルトのメッシュを生成する．
mesh4 についてのミーゼス応力と総変形を可視化する．
フィレット周辺のメッシュを調整する．
デフォルトのメッシュを生成する．
mesh5 についてのミーゼス応力と総変形を可視化する．

これで，調整後でも応力値がほぼ同じであることが分かる．フィレットを使うことで，非常にとがった角よりも物理的に現実的にもなる．

最大のミーゼス応力と総変形を可視化する関数．

検証

固体力学検証テストをご参照いただきたい．検証モデルは数学モードを検証するということを理解することが重要である．数学モードが実際の状態を表すかどうかは別の問題である．使用するモデルの限界を理解することが必要である．例えば，平面応力モデルの低誤差は，そのモデルが特定の場面に使えなければたいした意味はない．

大規模な有限要素モデル

固体力学PDEモデルを設定し解くことで，最終的に解くのに時間と大量のメモリを必要とする大規模なPDEモデルになる場合がある．これは大規模なPDEモデルすべてにとっての挑戦となるため，有限要素法のドキュメントにはさまざまな解法を紹介するセクションがある．

用語集

参考文献

1.  Backstrom, G; Simple Deformation and Vibration; GB Publishing; 2006

2.  The Efficient Engineer, https://www.youtube.com/watch?v=aQf6Q8t1FQE; retrieved 2021/4/27

3.  Cook, R; et. al.; Concepts and Applications of Finite Element Analysis; Wiley 2002; 4th Edition

4.  Bhatti, M.; Fundamental Finite Element Analysis and Application; Wiley 2005

5.  Bhatti, M.; Advanced Topics in Finite Element Analysis of Structures; Wiley 2006

6.  Brand, M.; Grundlagen FEM mit Solidworks; Vieweg+Teuber 2011; ISBN: 978-3-8348-1306-0

7.  https://en.wikipedia.org/wiki/Modal_analysis_using_FEM, retrieved 01/06/2021

8.  Nasdala, L.; FEM-Formelsammlung Statik und Dynamik; Vieweg+Teuber 2010

9.  Alawadhi, E. M.; Finite Element Simulations Using ANSYS; CRC Press 2010

10.  Constantinescu A., Korsunsky A.; Elasticity with Mathematica; Cambridge University Press 2007

11.  Bower, A.; Applied Mechanics of Solids; CRC Press 2010; ISBN-13: 978-1439802472; (http://solidmechanics.org/contents.php)

12.  The Efficient Engineer, https://www.youtube.com/watch?v=xkbQnBAOFEg; retrieved 2021/9/13

13.  Kelly, P.A.; Mechanics Lecture Notes: An introduction to Solid Mechanics. (http://homepages.engineering.auckland.ac.nz/~pkel015/SolidMechanicsBooks/index.html)

14.  Moosbrugger, C. (Editor); Atlas of Stress-strain Curves; ASM International 2002; ISBN 0-87170-739-X (https://books.google.com/books?id=up5KS9fd_pkC)

15.  Sifakis, E., Barbič, J.; Finite Element Method Simulation of 3D Deformable Solids; M&C Publishers 2016; ISBN 9781627054423

16.  Holzapfel, A.; Nonlinear Solid Mechanics; Wiley 2020; ISBN 978-0-471-82319-3

17.  Bonet, J., Wood, R.; Nonlinear Continuum Mechanics for Finite Element Analysis; Cambridge University Press 1997; ISBN 0-521-57272-X

18.  Mooney, M. Large elastic deformations of isotropic materials. Journal of Applied Physics 1940. 11: 582–592. https://doi.org/10.1063/1.1712836

19.  Rivlin, R. S. Large elastic deformations of isotropic materials. Philosophical Transactions of the Royal Society of London. Mathematical and Physical Sciences 1948. 241: pg 379–397. https://doi.org/10.1098/rsta.1948.0024

20.  Arruda, M. Boyce, M. A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials. Journal of the Mechanics and Physics of Solids. 1993. Volume 41, Issue 2. Pages 389-412,. https://doi.org/10.1016/0022-5096(93)90013-6

21.  Gent, A. N. A New Constitutive Relation for Rubber. Rubber Chemistry and Technology. 1996. 69 (1): 59–61. https://doi.org/10.5254/1.3538357

22.  Treolar, L.R.G.; Stress-strain data for vulcanised rubber under various types of deformation; Transactions of the Faraday Society, 1943

23.  Ning, X., Zhu, Q., Lanir, Y., Margulies, S. S.: A Transversely Isotropic Viscoelastic Constitutive Equation for Brainstem Undergoing Finite Deformation. Journal of Biomechanical Engineering, 2006. 128: pg 925-933. DOI: 10.1115/1.2354208

24.  Spencer, A. J. M.: Deformations of fibre-reinforced materials, 1972

25.  Nekliudova, E. A., Semenov, A. S., Melnikov, B. E., Semenov, S. G.: Experimental research and finite element analysis of elastic and strength properties of fiberglass composite material. Magazine of Civil Engineering 3 (47), 2014

26.  Otero F., Oller S., Martinez X., Salomón O.: Numerical homogenization for composite materials analysis. Comparison with other micro mechanical formulations. Composite Structures. 2015 Apr 1;122:405-16.

27.  Charalambakis, N.: Homogenization techniques and micromechanics: A survey and perspectives. 2010: 030803.

28.  Pascon, J. P. (2019). Large deformation analysis of plane-stress hyperelastic problems via triangular membrane finite elements. International Journal of Advanced Structural Engineering, vol 3, pp. 331-350, https://doi.org/10.1007/s40091-019-00234-w.

29.  Yosibash, Z. Weiss, D. Hartmann, S. (2014). High-order FEMs for thermo-hyperelasticity at finite strains. Computers and Mathematics with Applications, vol 67, pp. 477-496, https://doi.org/10.1016/j.camwa.2013.11.003

30.  Wriggers, P. Nonlinear finite element methods; Springer 2008

31.  Gasser, T. C. Vascular Biomechanics: Concepts, Models and applications; Springer 2021 https://doi.org/10.1007/978-3-030-70966-2

32.  Ting, T.C.: Anisotropic elasticity: theory and applications. Vol. 45. Oxford university press, 1996.