三角形要素の有限要素法計算
2013年9月17日
後 保範
1
対象方程式
2次元拡散方程式
u
u
k 2 k 2 f in
x
y
2
2
境界条件
u g on 1 ,
u
on 2
n
2
2
三角形要素と要素行列
(x3,y3)
(x1,y1)
(x2,y2)
要素行列(対称)
a11 a12 a13
a22 a23
対 称 a33
右辺
ベクトル
b1
b2
b3
3
要素行列の計算１
1
S det 1
x1
x2
y1
y 2 : 面積の2倍(三角形要素)
1
x3
y3
x1 y 2 x 2 y3 x3 y1 x 2 y1 x3 y 2 x1 y3
( x3 x 2) 2 ( y 2 y 3) 2
a11 k
2S
( x1 x3) 2 ( y3 y1) 2
a 22 k
2S
( x 2 x1) 2 ( y1 y 2) 2
a33 k
2S
4
要素行列の計算２
( x3 x 2)(x1 x3) ( y 2 y3)( y3 y1)
a12 k
2S
( x3 x 2)(x 2 x1) ( y 2 y3)( y1 y 2)
a13 k
2S
( x1 x3)(x 2 x1) ( y3 y1)( y1 y 2)
a 23 k
2S
S
b1 b2 b3 f ,
Sは3角形の2倍の面積
6
5
境界要素と右辺ベクトル
1境界 2境界
u g on 1
u
on 2
n
l
P1
P2
節点P1,P2の値を
gとして,対応する
右辺ベクトへ
(x1,y1)
(x1,y1)
b1
を計算
b2
6
境界要素右辺ベクトル計算
l
l
b1 , b 2
2
2
l
x1 x 2
2
y1 y 2
境界条件
u
on 2
n
2
l
(x2,y2)
(x1,y1)
7
入力データ１（三角形要素）
Element[m][5]
要素
番号
0
1
…
m-1
構成 節点 番号
NO1 NO2 NO3
0
1
2
1
2
3
…
…
…
n-3
n-2
n-1
配列には
持たない
Node[n][3]
物性 値
k
f
0
0
0
1
…
…
2
3
Ktbl[]
Ftbl[]
8
入力データ２（節点）
XY[n][2]
節点
番号
0
1
…
n-1
配列には
持たない
節点 番号
X座標 Y座標
0.0
0.0
0.1
0.0
…
…
1.0
1.0
0は未知数となる節点
Node[n]
節点
拘束
0
0
…
1
Fix[]
9
出力データ１（要素行列）
ElMat[m][3][3]: double
a11
a12
a13
a22
対 称
a23
a33
m
ElB[m]: double --- 右辺ベクトル用
10
出力データ2（密行列）
=
A
Aの形は各種
今回は密行列
x
b
作成
11
共通領域の確保例(Java)
public int Element[][] = new int[1000][5];
public int Node[] = new int[1000];
public double XY[][] = new double[1000][2];
public double ElMat[][][] = new double[1000][3][3];
public double ElB[] = new double[1000];
public double Ktbl[] = new double[100];
public double Ftbl[] = new double[100];
public double Fix[] = new double[100];
public double A[][] = new double[1000][1000];
public double B[] = new double[1000];
12
要素行列の計算1（全体）
public void SetElM(int m)
{ int k, nd1, nd2, nd3;
double x1,y1, x2,y2, x3, y3, S;
double kv, fv;
for (k=0; k<m; k++)
{ 準備計算（節点番号、座標、面積、物性値(k,f)）
要素行列の計算（ElMat[k][][], ElB[k]）
}
}
13
要素行列の計算２（準備計算）
nd1 = Element[k][0]; nd2 = Element[k][1];
nd3 = Element[k][2];
kv = Ktbl[Element[k][3]];
fv = Ftbl[Element[k][4]];
x1 = XY[nd1][0]; y1 = XY[nd1][1];
x2 = XY[nd2][0]; y2 = XY[nd2][1];
x3 = XY[nd3][0]; y3 = XY[nd3][1];
S = Math.abs(x1*y2 + x2*y3 + x3*y1
- x2*y1 - x3*y2 - x1*y3);
//面積*2
14
要素行列の計算３（要素行列）
kv = kv / (2.0*S);
ElMat[k][0][0] = kv*( (x3-x2)*(x3-x2) + (y2-y3)*(y2-y3) );
ElMat[k][1][1] = kv*( (x1-x3)*(x1-x3) + (y3-y1)*(y3-y1) );
ElMat[k][2][2] = kv*( (x2-x1)*(x2-x1) + (y1-y2)*(y1-y2) );
ElMat[k][0][1] = kv*( (x3-x2)*(x1-x3) + (y2-y3)*(y3-y1) );
ElMat[k][0][2] = kv*( (x3-x2)*(x2-x1) + (y2-y3)*(y1-y2) );
ElMat[k][1][2] = kv*( (x1-x3)*(x2-x1) + (y3-y1)*(y1-y2) );
ElB[k] = fv*S/6.0;
ElMat[k][1][0] = ElMat[k][0][1];
//対称なの以下は
ElMat[k][2][0] = ElMat[k][0][2];
//省略可能
ElMat[k][2][1] = ElMat[k][1][2];
15
密行列作成１（全体）
public void SetAB(int m, int n)
{ int i, j, k, IP, JP, IFX, JFX;
行列A,右辺ベクトルのクリア
for (k=0; k<m; k++)
{ for (i=0; i<3; i++)
{ IP = Element[k][i]; IFX = Node[IP];
if (IFX == 0) { i節点は未知数処理 }
}
}
}
16
密行列作成２（i節点は未知数）
B[IP] += ElB[k]; // Set B[IP] by f
for (j=0; j<3; j++)
{ JP = Element[k][j];
JFX = Node[JP];
if (JFX == 0)
//共に未知数なら行列へ
{ A[IP][JP] += ElMat[k][i][j]; }
else
//一方が既知数ならbへ
{ B[IP] -= Fix[JFX]*ElMat[k][i][j]; }
}
17
対称帯行列の作成
M+1
j
A[i][j]: 密行列
A[i][j-i]: 帯行列
i
public double A[][] =
new double[1000][1000];
100: (m+1)以上
18
行列作成２（対称帯行列用）
B[IP] += ElB[k]; // Set B[IP] by f
for (j=i; j<3; j++) j=0をj=iに変更
{ JP = Element[k][j];
JFX = Node[JP]; A[IP][JP]をA[IP][JP-IP]に変更
if (JFX == 0)
//共に未知数なら行列へ
{ A[IP][JP-IP] += ElMat[k][i][j]; }
else
//他方が既知数ならｂへ
{ B[IP] -= Fix[JFX]*ElMat[k][i][j]; }
}
19
入力データ３（境界要素）
Bound[l][3]
境界要素
番号
0
1
…
l-1
配列には
持たない
構成節 点番号
NO1
NO2
0
1
1
2
…
…
l-1
l
物性値
α
0
1
…
3
Node[n][3]
Alph[]
20
© Copyright 2024 ExpyDoc
