total Lagrange と updated Lagrange 2
2. 局所作用の原理 (principle of local action) : 物質点 X の応力の決定において は , X の近傍の運動だけが関与し , その外の物質点の運動は無視できる
3. 物質客観性の原理 (principle of material frame indiffernce or principle of
mate-rial objectivity) : 構成式は基準枠 (referance frame)
1の変化の下で不変でなけ ればならない . すなわち , 2 つの基準枠 O
∗と O が
x
∗= c( t ) + Q( t ) · x (81) の関係にあるとき
T
∗= Q( t ) · T · Q( t )
T(82) の関係があり , また構成式は両枠に対して同一の形をとらねばならない .
1事象(event)と呼ばれる空間内の位置 xと時間 tの組{x, t} を観測する観測者(observer)のこと. 基準枠が異なると事象は {x, t},{x∗, t∗} とベクトル自体が異 なって観測される.
36
構成式 4
• 2つの基準枠 O∗ と O が
x∗ = c(t) + Q(t)· x (83)
• 2 点を結ぶベクトルは b∗ = x∗2 −x∗1, b = x2 − x1 とすると,
x∗2 −x∗1 = Q(t)(x2 − x1) (84) b∗ = Q(t) · b (85) このように両基準枠から観測されるベクトルが,限時刻での Q(t) のみを介した形で関係づ けられるとき,そのベクトルには客観性があるという.
• テンソルはベクトルの線形変換なので,以下の関係を満たすとき客観性があるという
K∗ = Q(t) · K · Q(t)T (86)
• 客観性のあるテンソルの例
V ,B,A,D,T (87)
• 客観性のないテンソルの例
R,L,W (88)
• 観測不変テンソルの例
U,C,E,S (89)
37
構成式 4
• 不合理な例:線形弾性体
εij(u) = 1 2
∂ui
∂Xj + ∂uj
∂Xi
(90) Tij = κ(divu)δij + 2GεDij(u) (91)
• x1 −x2 平面内で 90◦ 回転,応力は0のはず
Fij =
⎡
⎢⎣
0 −1 0 1 0 0 0 0 1
⎤
⎥⎦ (92)
Eij = 1 2
⎛
⎜⎝
⎡
⎢⎣
−1 −1 0 1 −1 0 0 0 0
⎤
⎥⎦+
⎡
⎢⎣
−1 1 0
−1 −1 0 0 0 0
⎤
⎥⎦
⎞
⎟⎠
=
⎡
⎢⎣
−1 0 0 0 −1 0 0 0 0
⎤
⎥⎦ (93)
• 微小変形の場合は,対角項は cosθ ≈ 1 非対角項はsinθ ≈ θ なので,問題ないが,有限変形 を考えると,不合理を生じる.
38
構成式 4
• 以下のような構成式を考えたとする.
T = a E (94)
• 観測者が変わると T
∗= a
∗E
∗となるが,
– T
∗= QT Q
T(客観性がある)
– a
∗= a (スカラー)
– E
∗= E (観測不変量)
• 従って
T
∗= a
∗E
∗→ QT Q
T= a E = T (95) となり,矛盾する.
• たとえば,以下のものは矛盾しない.
T = a A (96)
S = a E (97)
39
長方形断面の梁要素
• 一般的な偏微分方程式を素直に離散化するなら、ソリッド要素が自然
• しかし「構造解析」をする場合には、別な考え方がある
• 梁要素は代表的な構造要素
• 断面の寸法が長手方向の寸法より十分に小さいスレンダーな部材からなる構 造物の解析に使用される .
• 同構造をソリッド要素でモデル化するのに対して , 大幅に総自由度やモデル化 の手間が軽減される .
• 低次のソリッド要素で梁をモデル化すると,ロッキングする
• 時刻 0 において,梁の中立軸に垂直な平面は変形の間も平面を保つが必ずし も変形した中立軸に垂直である必要はない
• 梁の断面は変形しない.
4
座標及び変位の補間関数 1
• 時刻tで節点nの断面を張るベクトルを tV(n)2 , tV(n)3 , また時刻 t での節点 n の位置ベクトルを tx とする
全体座標系
局所座標系
a b
e1 e2
e3
X1 X2
X3 r1
r2 r3
tV2 tV3
1 1
2 2
3 3
4 4
図 1: 座標系の定義
• 時刻 t における要素内の任意の点における位置ベクトルは,
tx(r1, r2, r3) = N(n)(r1)
txn+ a
2r2tV(n)2 + b
2r3tV(n)3
(1)
と, 表すことができる.
• なお, 形状関数 N(n) には1次元要素のものを用いる.
5
座標及び変位の補間関数 2
• 局所座標ベクトル tV (2n), tV (3n) は梁の断面に固定されており, 要素の変形とともに回転するこ とを考慮すれば, 時刻 t から t(= t +∆t) までの間の変位増分ベクトル u は
u = tx−tx
= N(n)
tx(n)−tx(n)+ a 2r2
tV(n)2 −tV(n)2 + b
2r3
tV(n)3 −tV(n)3
となる.
• ここで∆t間の回転増分は微小であると仮定すれば,
tV (in) tV (in) + θ(n) × tV (in) (2)
と表すことができる.
• よって,
tV(n)i −tV(n)i = θ(n)×tV(n)i
=
⎡
⎢⎣
0 −θ3(n) θ2(n) θ3(n) 0 −θ(n)1
−θ(n)2 θ(n)1 0
⎤
⎥⎦
⎧⎪
⎨
⎪⎩
tVi1(n)
tVi2(n)
tVi3(n)
⎫⎪
⎬
⎪⎭ =
⎧⎪
⎨
⎪⎩
θ2(n)tVi3(n) −θ3(n)tVi2(n) θ3(n)tVi1(n) −θ1(n)tVi3(n) θ1(n)tVi2(n) −θ2(n)tVi1(n)
⎫⎪
⎬
⎪⎭
=
⎡
⎢⎣
0 tVi3(n) −tVi2(n)
−tVi3(n) 0 tVi1(n)
tVi2(n) −tVi1(n) 0
⎤
⎥⎦
⎧⎪
⎨
⎪⎩ θ(n)1 θ(n)2 θ(n)3
⎫⎪
⎬
⎪⎭ (3)
6
座標及び変位の補間関数 3
• ここで tx,u の全体座標系での成分表示をしておく.
⎧⎪
⎨
⎪⎩
tx1
tx2
tx3
⎫⎪
⎬
⎪⎭ = N(n)(r1)
⎧⎪
⎨
⎪⎩
tx(n)1
tx(n)2
tx(n)3
⎫⎪
⎬
⎪⎭+ a
2r2N(n)(r1)
⎧⎪
⎨
⎪⎩
tV21(n)
tV22(n)
tV23(n)
⎫⎪
⎬
⎪⎭+ b
2r3N(n)(r1)
⎧⎪
⎨
⎪⎩
tV31(n)
tV32(n)
tV33(n)
⎫⎪
⎬
⎪⎭ (4)
⎧⎪
⎨
⎪⎩ u1 u2 u3
⎫⎪
⎬
⎪⎭ = N(n)(r1)
⎧⎪
⎨
⎪⎩ u(n)1 u(n)2 u(n)3
⎫⎪
⎬
⎪⎭+ a
2r2N(n)(r1)
⎡
⎢⎣
0 tV23(n) −tV22(n)
−tV23(n) 0 tV21(n)
tV22(n) −tV21(n) 0
⎤
⎥⎦
⎧⎪
⎨
⎪⎩ θ1(n) θ2(n) θ3(n)
⎫⎪
⎬
⎪⎭
+ b
2r3N(n)(r1)
⎡
⎢⎣
0 tV33(n) −tV32(n)
−tV33(n) 0 tV31(n)
tV32(n) −tV31(n) 0
⎤
⎥⎦
⎧⎪
⎨
⎪⎩ θ1(n) θ2(n) θ3(n)
⎫⎪
⎬
⎪⎭ (5)
7
シェル要素の定式化
• シェル要素の開発は過去多くの研究者により行われてきたが, 現在も高精度, 高信頼性かつ計 算効率のよい要素を目指して開発が続けられている. 優れたシェル要素の条件は,
– 薄肉, 厚肉双方のシェルに適用できる (薄肉シェル要素においてもロッキングが生じない).
– 任意の形状のシェルを取り扱える.
– 剛体モード以外の虚偽のゼロ固有値を持たない. – 要素のゆがみに対して解の精度が損なわれない.
– 要素内の任意の点においてゼロを含む一定ひずみを表現できる.
などが挙げられる. これらを考慮すると, Batheら[?]が 1984 年以降に発表した MITC (Mixed Interpolation of Tensorial Components) 要素が理論の明解さ, 定式化の容易さの点で最も優れて いる.
50
幾何学的非線形シェル要素の基礎式
• Bathe らの MITC 要素の概要を示す [?] [?]. この要素の特徴は, 1. アイソパラメトリック degenerated シェル要素である.
2. 面外せん断ひずみに関しては, あるサンプリング点の面外せん断ひずみの値から内挿する ように再定義する.
3. ひずみ成分および応力成分には, 自然座標系での共変成分, 反変成分を使用する.
などが挙げられる. 上記の 2, 3 がMixed Interpolation of Tensorial Components 要素と呼ばれる 理由であり, Assumed-Strain 要素と参照されている場合もある.
51
形状および変位補間式
• 図に時刻 0 におけるシェル要素の形状を示す. 以下, 特に記載のない場合には図中の記号で 説明をおこなう. 図 5 より時刻 0 における要素内の任意の点の位置ベクトル X は形状関数 N(r1, r2) を用いて,
X = 4
n=1
Nn
r1, r2
Xn + a 2r3
4 n=1
Nn
r1, r20
V n3 (160)
で表される. 式(160)で 0V k3 はシェルの厚み方向を示すベクトル(以下, ディレクター. 左肩 の添字は時刻を表す)であり, 解析面に対する法線として一節点にただひとつ定義され, その 節点を有する要素間で共有される.
ここでシェルの変形に関して次の 2 つの仮定を設ける.
1. 各節点におけるディレクター(時刻 0 において直線を仮定した)は, 変形の間も直線を保つ が, 必ずしも中立面に垂直である必要はない(面外せん断変形を許す).
2. シェルの肉厚 a は変形の間も変化しない(大ひずみ領域は取り扱わない).
このとき時刻 t における要素内の任意の点 tx の位置ベクトルは,
tx = 4
n=1
Nn
r1, r2t
xn + a 2r3
4 n=1
Nn
r1, r2t
V n3 (161)
52
となる. したがって, 時刻 t における変位 tu , 時刻 t から t(= t + ∆t) までの変位増分 u は それぞれ tu = tx − X,u = tx −tx であるから,
tu = 4
n=1
Nn(r1, r2)tUn + a 2r3
4 n=1
Nn(r1, r2)(tV n3 − 0V n3) (162) u =
4 n=1
Nn(r1, r2)Un + a 2r3
4 n=1
Nn(r1, r2)(tV n3 − tV n3) (163) となる. ここで, 各節点のシェルディレクターの時刻 t から t までの有限回転を表すテンソル を ttR とすると, tV k3 = ttR tV k3 であるから式(163)はつぎのように書き直される.
u = 4
n=1
Nn(r1, r2)Un + a 2r3
4 n=1
Nn(r1, r2)(ttR− I) tV n3 (164) (I は単位テンソル)
53
a
A B
C
D
1 2
3
4
e1
e2
e3
G1
G2
G3
x1
x2 x3
r1 r2
r3
α4
β4
u11 u12 u13
0V41
0V42
0V43
a shell thickness
xi Cartesian coordinates
ei orthonormal base vector in global coordinate system ri natural coordinates ( -1≤ ri ≤1)
Gi covariant base vector of natural coordinate ri
Vki orthonormal base vector in local coordinate system at node k Vk3 shell director vector
uk displacement vector at node k
αk, βk rotations of the director vector aboutVk1 and Vk2 図 5: 4 節点アイソパラメトリックシェル要素
54
面外せん断ひずみ成分補間式
• 通常のアイソパラメトリック要素では, 要素内任意点のひずみ, ひずみ増分は式(162), (163)で 表せる変位, 変位増分の空間微分をとることにより直接求められるが, MITC 要素では, 面外 せん断ひずみの成分についてのみ, あるサンプリング点における面外せん断ひずみ (通常の方 法から求める) から内挿関数を再定義して求める. 図 6 にサンプリング点位置 (点A〜D) を,
式(165)に面外せん断ひずみ成分の補間式を示す.
E13 = 1
2(1 + r2)E13A + 1
2(1 − r2)E13C (165) E23 = 1
2(1 + r1)E23D + 1
2(1 − r1)E23B (166) ここで,
E13, E23 : 要素内の任意点のせん断ひずみ
EA∼D : サンプリング点におけるせん断ひずみ である.
55
2 1
3
4
A
A
B B
C C
D
D
r1
r1 r1
r1 r1
r2 r2
r2 r2
r2
EA13
EC13
EB23 ED23
図 6: 面外せん断ひずみのサンプリング点 56
ALE 流体解析
• Lagrange 表示で表される座標系を物質座標系
• Euler 表示で表される座標系を空間座標系
• これらの座標系とは独立した任意の座標系,すなわ ち参照座標系を設定し,その座標系上で物体の運
動を記述することを参照表示 (ALE 表示 ) と定義する.
• 今, Lagrange , Euler , ALE の各表示により解析領
域および領域内のある点を次のように表す
時間導関数の関係
• 実質時間導関数と物質時間導関数の関係式
• 実質時間導関数 ( 物質時間導関数 ) と空間時間導関数の 関係式
• 実質時間導関数物 ( 物質時間導関数 ) と参照時間導関数
の関係式
各種の速度の物理的解釈
• Euler 座標系に対する物質点の速度 v .
• 参照座標系に対する物質点の速度 w .
• Euler 座標系に対する参照座標系の速度
Euler 領域での時間導関数の式
• 位置ベクトルを x とすれば,速度関係式は
• 参照座標系に対する物質点の相対速度 を導入 すると
• これより
参考までに
Eulerだと
ALE 流体構造連成解析(1)
• Navier-Stokes 方程式と連続の式の ALE 表示
• 非圧縮性 Newton 流体(構成則)
• 弾性体の運動方程式
ALE 流体構造連成解析(2)
• 連成面での力学的平衡条件式と幾何学的連 続条件式
• 連成系の弱形式
有限要素離散化したマトリックス方程式
• 流体
• 構造
連成系のマトリックス方程式
• ただし肩符号 c (common) は連成面上の自由
度で、肩符号 i (independent) はそれ以外の
自由度とする
マトリックスの成分
補足:大まかなプ ログラムの流れ
• メッシュの制御
• 時間積分
補足:メッシュ制御
• 流体部分を構造が取り囲んでいるような系、total
submerged system の場合、連成系のマトリックスを解く と構造の変位が得られる。
• 流体領域を弾性体に見立てて、流体と固体の境界上の 節点は、固体の解析で得られた変位を、変位境界条件と して与えて、弾性体の剛性方程式を解き、内部の節点の 移動量を計算する。
• 流体の節点のメッシュの速度は、単純に変位をΔt で割 るだけでもよい。(経験的にはそれほど差は無い)
• なお、弾性体でなく、Laplace 方程式でも良いし、線形に 配分しても良いし、もっと複雑な構成式(超弾性体)などを 使っても良い。
• ただし、逆に言うと、 ALEはそれで対応できる範囲に限ら れると考えたほうがよい。構造の変形にともなって。流体 メッシュが大きくゆがんだり、極めて間隔が狭くなるような 場合は対応できない。メッシュを切りなおして、再補間す るなどの処理が必要。
?
補足:時間積分は?
• 固体を含まない流体の有限要素法では、速 度と加速度の項だけ。
• 固体があると、速度、加速度に加えて、変位 が未知数になるが、どうやって積分するのか。
• 一般的には陰解法で Newmark- β法をつか
う。
動的解析手法
• 静的な仮想仕事の式
V
SijδEijdV =
V
giρδuidV +
S
tiδuidS (1)
加速度の項がある場合
V
ρaiδuidV +
V
SijδEijdV =
V
giρδuidV +
S
tiδuidS (2)
• 時刻t+ Δtでの運動方程式
M ·t+Δtu¨ +C ·t+Δtu˙ +t+ΔtQ= t+ΔtF (3)
• M は質量マトリックス,C は減衰マトリックス. 実はC は仮想仕事の式からは導かれない。現実問題では、固 定部の摩擦などによって減衰が生じる。また材料も粘性を持つのでこれからも減衰を生じる。これをあわせて C としているので実際には条件を変えて実験を行って C を妥当な値に調整する、ということが行われている。
• ここでは,この方程式をそのままの形で時間方向に積分する直接時間積分法についてのべる.
• 線形問題では,固有モードの重ね合わせで解析する場合もあり,モード重ね合わせ法という.
• 時刻 t での解が得られているとき,この解を基にして t+ Δt の解を厳密に求める方法を陰解法(線形加速度法,
Newmark-β 法など),t の解を基に予測する方法を陽解法(中央差分)という
• 通常は M, C は時刻によらず一定.しかし Q は変化する
• 内力ベクトルについて線形化し,Newton-Raphson 法を用いる.
M ·t+Δtu¨(k) +C ·t+Δtu˙(k)+t+ΔtK(k−1)Δu(k) = t+ΔtF −t+ΔtQ(k−1) (4)
t+Δtu(k) = tu(k−1) + Δu(k) (5)
t+Δtu˙(k) = tu˙(k−1) + Δ ˙u(k) (6)
t+Δtu¨(k) = tu¨(k−1) + Δ¨u(k) (7)
4
• このままでは,未知量として Δu(k), Δ ˙u(k), Δ¨u(k) が含まれるため,通常は,この三者を何らかの仮定の下に関 係づけ,どれか1つの変数にまとめて解析する.
5