Comments
Description
Transcript
RS2行列
固体力学
大学院工学研究科機械工学専攻
固体力学（２単位）
００４D
有限要素法（要素マトリックス計算プログラム）
構成
シミュレーション
計算機としてのＰＣ
シミュレーションのための計算機環境
ハードウエア
ソフトウエア
シミュレーションのための計算法
有限要素法
メッシュ
連立一次方程式の解法
可視化
参考書
有限要素法概説
有限要素法
O.C.Zienkiewicz and K.Morgan, Dover
Finite Element Implementation
矢川元基・吉村忍著，培風館
Finite Elements & Approximation
菊地文雄著，サイエンス社
Y.K.Cheung, S.H.Lo, A.Y.T.Leung, Blackwell
Numerical Computation
C.W.Ueberhuber, Springer
シミュレーションのための計算法
有限要素法
１．２次元定常熱伝導問題
２．２次元応力解析
３．数値積分法
有限要素法（入力データ１）
(2.0, 2.0, 2.0)
25
26
27
y
x
23
22
z
24
19
20
5
21
18
6
15
10
11
1
12
2
9
6
1
(0.0, 0.0, 0.0)
2
3
要素数：８
節点数：２７
有限要素法（入力データ２）
13
10
14
1
11
1
14
11
10
20
14
11
11
21
14
13
23
26
23
17
14
9
6
26
7
15
5
25
15
18
4
8
5
22
12
17
14
4
24
6
17
3
6
23
20
16
13
3
23
5
12
2
22
19
2
5
2
15
14
27
8
24
18
15
有限要素法（入力データ３）
要素データ
1
2
3
4
5
6
7
8
8
8
8
8
8
8
8
8
1
1
1
1
1
1
1
1
1 2 5
2 3 6
4 5 8
5 6 9
10 11
11 12
13 14
14 15
4 10 11 14 13
5 11 12 15 14
7 13 14 17 16
8 14 15 18 17
14 13 19 20 23
15 14 20 21 24
17 16 22 23 26
18 17 23 24 27
22
23
25
26
各行の構成
要素番号，要素タイプ，材料番号，節点番号ｘN
有限要素法（入力データ４）
節点データ
1 0.0 0.0 0.0
2 1.0 0.0 0.0
3 2.0 0.0 0.0
4 0.0 1.0 0.0
5 1.0 1.0 0.0
6 2.0 1.0 0.0
・・・・・・・
25 0.0 2.0 2.0
26 1.0 2.0 2.0
27 2.0 2.0 2.0
各行の構成
節点番号，ｘ座標，ｙ座標，ｚ座標
有限要素法（入力データ５）
材料データ
1 1 2.0e11 0.3
2 1 7.0e10 0.35
各行の構成
材料物性番号，材料タイプ，縦弾性係数，ポアソン比
要素マトリックス計算コード
8
６面体１次要素
7
5
1
N , , 1 i 1 i 1 i
8
e
i
6
1
節点数：８
要素あたり自由度：２４
要素マトリックス：２４ｘ２４
3
2
要素マトリックス計算（全体フロー）
入力データ読み込み，変数(配列）初期化
形状（基底）関数値（微分）計算
積
分
点
ル
ー
プ
Bマトリックス算出
DBマトリックス算出
BDBマトリックス算出
要素マトリックスにBDBを加算
終了
要素マトリックス計算（入出力データ）
void elmlib_03_Dynamic(
element_data element,
node_data *node,
mate_data mate,
double **esm,
int ngauss,
double *gc,
double *gw,
int nfpn)
{
１節点あたりの自由度
（入力データ）
要素を構成する節点の番号を格納した
１次元配列（入力データ）
全節点の座標を格納した１次元配列
（入力データ）
要素の物性値を格納した１次元配列
(入力データ）
要素マトリックス（演算結果）を格納
する２次元配列（出力データ）
ガウス積分の１軸あたり積分点数
（入力データ）
ガウス積分の座標値（ｇｃ）および重
み（ｇｗ）を格納した１次元配列（入力
データ）
要素マトリックス計算（変数・配列）
int
i,j,k,ii,jj,kk,counter,necm=6,
ndeg=8,
nnpe=8,
kdim=24;
double e,v,ro,th,ee,coord[8][3],
ja[3][3],
invja[3][3],
d[6][6],
s,t,u,ra,
rs,sa,ss,ua,us,n[8][7],
b[6][24],
db[6][24],
dtmp01,
ra2,rs2,
ssus,ssua,saus,saua,
det,invdet;
Dマトリックスのサイズ
要素あたり節点数
要素あたり自由度
節点座標
ヤコビアン行列および
その逆行列
Dマトリックス
基底関数（微分値）
Bマトリックス
DBマトリックス
ヤコビアン行列の行列式
およびその逆数
要素マトリックス計算（配列初期化）
for(i=0;i<nnpe;i++){
ii = element.node[i];
for(j=0;j<nfpn;j++)
coord[i][j] = node[ii].f[j];
}
e
v
= mate.d[0];
= mate.d[1];
for(i=0;i<6;i++)
for(j=0;j<24;j++) b[i][j] = 0.0;
for(i=0;i<6;i++)
for(j=0;j<6;j++) d[i][j] = 0.0;
for(ii=0;ii<kdim;ii++){
for(jj=0;jj<kdim;jj++)
esm[ii][jj] = 0.0;
節点番号
全節点の座標値配列node
から当該要素８節点の座標
値をcoordにコピー
縦弾性係数
ポアソン比
Bマトリックス初期化
（ゼロクリア）
Dマトリックス初期化
（ゼロクリア）
要素マトリックス初期化
（ゼロクリア）
要素マトリックス計算（Dマトリックス）
ee = e*(1.0 - v)/(1.0 + v)/(1.0 - 2.0*v);
d[0][0]
d[1][1]
d[2][2]
d[3][3]
d[4][4]
d[5][5]
=
=
=
=
=
=
ee;
ee;
ee;
ee*(1.0 - 2.0*v)/2.0/(1.0 - v);
d[3][3];
d[3][3];
1 v
d[0][1] = ee*v/(1.0 - v);
v
d[0][2] = d[0][1];
v
d[1][2] = d[0][1];
E
0
d[1][0] = d[0][1];
1 v 1 2v
d[2][0] = d[0][2];
0
d[2][1] = d[1][2];
0
Dマトリックス設定
v
v
0
0
v
0
v
1 v
0
0
0
0
0
1 2v
2
0
0
0
1 2v
2
0
0
0
0
1 v
0
0
0
0
0
0
1 2v
2
要素マトリックス計算（積分点ループ）
for(i=0;i<ngauss;i++){
s = gc[i];
for(j=0;j<ngauss;j++){
t = gc[j];
for(k=0;k<ngauss;k++){
u = gc[k];
ra
rs
sa
ss
ua
us
=
=
=
=
=
=
ssus
saus
ssua
saua
(1.0
(1.0
(1.0
(1.0
(1.0
(1.0
=
=
=
=
+
+
+
-
s)*0.5;
s)*0.5;
t)*0.5;
t)*0.5;
u)*0.5;
u)*0.5;
ss*us;
sa*us;
ss*ua;
sa*ua;
積分点ループ
s, t , u , ,
1
ra 1
2
1
us 1
2
1
saus 1 1
4
要素マトリックス計算（基底関数の微分値
１）
n[0][0] = -0.5*ssus;
n[1][0] = 0.5*ssus;
n[2][0] = 0.5*saus;
n[3][0] = -0.5*saus;
n[4][0] = -0.5*ssua;
n[5][0] = 0.5*ssua;
n[6][0] = 0.5*saua;
n[7][0] = -0.5*saua;
rs2 = 0.5*rs;
ra2 = 0.5*ra;
n[0][1] = -rs2*us;
n[1][1] = -ra2*us;
n[2][1] = ra2*us;
n[3][1] = rs2*us;
n[4][1] = -rs2*ua;
n[5][1] = -ra2*ua;
n[6][1] = ra2*ua;
n[7][1] = rs2*ua;
N i
N i
N1 , ,
1
1 1 1
8
N1 , ,
1
1 1
8
N1 , ,
1
1 1
8
N1 , ,
1
1 1
8
要素マトリックス計算（基底関数の微分値
２）
n[0][2]
n[1][2]
n[2][2]
n[3][2]
n[4][2]
n[5][2]
n[6][2]
n[7][2]
=
=
=
=
=
=
=
=
-rs2*ss;
-ra2*ss;
-ra2*sa;
-rs2*sa;
rs2*ss;
ra2*ss;
ra2*sa;
rs2*sa;
n[0][6]
n[1][6]
n[2][6]
n[3][6]
n[4][6]
n[5][6]
n[6][6]
n[7][6]
=
=
=
=
=
=
=
=
rs*ssus;
ra*ssus;
ra*saus;
rs*saus;
rs*ssua;
ra*ssua;
ra*saua;
rs*saua;
N i
N1 , ,
1
1 1
8
1
N1 , , 1 1 1
8
Ni
要素マトリックス計算（ヤコビアン行列）
for(ii=0;ii<nfpn;ii++){
for(jj=0;jj<nfpn;jj++){
ja[ii][jj] = 0.0;
for(kk=0;kk<ndeg;kk++)
ja[ii][jj] += n[kk][ii]*coord[kk][jj];
}
}
x
x
J
x
y
y
y
z N i xi
z N i xi
z N i xi
N i yi
N i yi
N i yi
N i zi N
i
xi
N i zi
N
i xi
N i zi N i
xi
N i
yi
N i
yi
N i
yi
N i
zi
N i
zi
N i
zi
要素マトリックス計算（ヤコビアン行列式）
det =
+
+
-
ja[0][0]*ja[1][1]*ja[2][2]
ja[0][1]*ja[1][2]*ja[2][0]
ja[0][2]*ja[1][0]*ja[2][1]
ja[0][0]*ja[1][2]*ja[2][1]
ja[0][1]*ja[1][0]*ja[2][2]
ja[0][2]*ja[1][1]*ja[2][0] ;
invdet = 1.0/det ;
a b
c
d
e
f a
g
h
i
e
f
h
i
d
b c
h i
g
b
c
e
f
要素マトリックス計算（行列式）
a11 a1n
sgn P a1 p1 a2 p2 a3 p3 anpn
an1 ann
P
1 2 n
P
pn
p1 p2
1
sgn P
1
P : even permutatio n
P : odd permutatio n
要素マトリックス計算（ヤコビアン逆行列）
invja[0][0]
invja[0][1]
invja[0][2]
invja[1][0]
invja[1][1]
invja[1][2]
invja[2][0]
invja[2][1]
invja[2][2]
=
=
=
=
=
=
=
=
=
(ja[1][1]*ja[2][2]
(ja[0][2]*ja[2][1]
(ja[0][1]*ja[1][2]
(ja[1][2]*ja[2][0]
(ja[0][0]*ja[2][2]
(ja[1][0]*ja[0][2]
(ja[1][0]*ja[2][1]
(ja[0][1]*ja[2][0]
(ja[0][0]*ja[1][1]
A
1
i, j
-
1
ja[1][2]*ja[2][1])*invdet;
ja[0][1]*ja[2][2])*invdet;
ja[1][1]*ja[0][2])*invdet;
ja[1][0]*ja[2][2])*invdet;
ja[0][2]*ja[2][0])*invdet;
ja[0][0]*ja[1][2])*invdet;
ja[1][1]*ja[2][0])*invdet;
ja[0][0]*ja[2][1])*invdet;
ja[0][1]*ja[1][0])*invdet;
i j
A
~
A j ,i
要素マトリックス計算（逆行列）
a11 a1n
A
a a
nn
n1
余因子
i j
a~ij 1 Dij
余因子行列
a11
a
i 1,1
Dij
ai 1,1
a
n ,1
a~11 a~21 a~n1
~
~
~
an 2
~ a12 a22
A
~
~
a~
a
a
2n
nn
1n
a1, j 1
a1, j 1
ai 1, j 1
ai 1, j 1
ai 1, j 1
ai 1, j 1
an , j 1
an , j 1
a1n
ai 1,n
ai 1,n
an ,n
1 ~
A
A
A
1
要素マトリックス計算（基底関数の微分値）
for(ii=0;ii<ndeg;ii++){
n[ii][3] = 0.0;
n[ii][4] = 0.0;
n[ii][5] = 0.0;
for(jj=0;jj<3;jj++){
n[ii][3] += invja[0][jj]*n[ii][jj];
n[ii][4] += invja[1][jj]*n[ii][jj];
n[ii][5] += invja[2][jj]*n[ii][jj];
}
}
N i x
N i x
N
i x
N i y N i z N i x
x y z
N i y N i z N i x
x y z
N i y N i z N i x
x y z
y
y
y
N ii
n[ii][3]
x
N ii
n[ii][4]
y
N ii
n[ii][5]
z
z N i
N i
x
x
z N i
N i
J
y
y
N
z N i
i
z
z
要素マトリックス計算（Bマトリックス）
for(ii=0;ii<ndeg;ii++)
jj = ii*nfpn;
b[0][jj]
= n[ii][3];
b[1][1+jj] = n[ii][4];
b[2][2+jj] = n[ii][5];
b[3][jj]
b[3][1+jj]
b[4][1+jj]
b[4][2+jj]
b[5][jj]
b[5][2+jj]
}
=
=
=
=
=
=
n[ii][4];
n[ii][3];
n[ii][5];
n[ii][4];
n[ii][5];
n[ii][3];
x
0
0
B
y
0
z
0
y
0
x
z
0
0
0
N
0
N8 0
0
1 0
z 0 N
0
0
N
0
1
8
0 0
0 N1
0
0 N 8
N 8
y N1
0
0
0
x
x
N 8
N1
x 0
0 0
y
y
N1
0
0
0
0
z
N N
N 8 N 8
1
1
0
y
x
y
x
N 8
N1 N1
0
0
z
y
z
N1
N 8
N1
0
0
x
z
z
0
0
N 8
z
0
N 8
y
N 8
x
要素マトリックス計算（DBマトリックス）
for(ii=0;ii<necm;ii++){
for(jj=0;jj<kdim;jj++){
db[ii][jj] = 0.0;
for(kk=0;kk<necm;kk++)
db[ii][jj] += d[ii][kk]*b[kk][jj];
}
}
1 v
v
v
E
0
D
1 v 1 2v
0
0
v
v
0
0
v
0
v
1 v
0
0
0
0
0
1 2v
2
0
0
0
1 2v
2
0
0
0
0
1 v
0
0
0
0
0
0
1 2v
2
N1
x
0
0
B N
1
y
0
N1
z
0
0
N1
y
0
0
N1
x
N1
z
0
N1
z
0
N 8
x
0
0
N 8
y
0
0
N 8
y
N 8
x
N 8
z
N1
0
y
N 8
N1
x
z
0
0
0
N 8
z
0
N 8
y
N 8
x
要素マトリックス計算（ガウス積分の重み）
dtmp01 = gw[i]*gw[j]*gw[k]*det;
X
(1,1)
f x, y dv
1 1
f , J dd
1 1
f , , ddd f , , H
1
1
1
n
m
l
i
j
k
i
i 1 j 1 k 1
e
e
1 1 1
Y
(1,1)
H j Hk
要素マトリックス計算（要素マトリックス）
for(ii=0;ii<kdim;ii++){
for(jj=0;jj<kdim;jj++){
for(kk=0;kk<necm;kk++)
esm[ii][jj]+=b[kk][ii]*db[kk][jj]*dtmp01;
}
}
e
e
e
N1
x
0
t
BDBd e
N pe
x
0
0
N
y
e
1
0
N pe
y
N1
y e
N
N1e 1
x
x
D 0
e
N p
N1e
y
N pe y
x
N pe
0
N
y
N1e
x
e
1
x
0
N pe
y
0
e
N p
d
y
N pe
x
全体剛性マトリックス構築
for(ii=0;ii<nnpe;ii++){
ig = element.node[ii] ;
for(jj=0;jj<nnpe;jj++){
jg = element.node[jj] ;
for(i=0;i<nfpn;i++){
for(j=0;j<nfpn;j++){
gsm[ig*nfpn+i][jg*nfpn+j] += esm[ii*nfpn+i][jj*nfpn+j];
}
}
}
}