有限要素法入門 (pdf) - CoMFoS

有限要素法入門
大塚厚二
広島国際学院大学情報学部
三角形分割
領域  の閉包  が三角形の集合   Jj1K j で，次の性質(T1)～
(T3)をもつとき，集合 Th  K j
j1,,J
(T1)
が領域  の三角形分割
K i , K j  Th に対して， i  jならば K i  K jは，空集合か，
三角形の辺の全体か，または 1 点となる。
(T2)
境界   上の折れ曲がり点は，三角形 K  Th の頂点
(T3)
0  h  max KTh diam  K  。ここで， diam  K  は，三角形
K  Th の辺の最大長。
FreeFem++による生成例
領域  が単位円盤
border C(t=0,2) { x=cos(t); y=sin(t); };
mesh Th = buildmesh(a(50));
とすると，次図の三角形分割が生成される。
三角形分割のデータ構造
border C(t=0,2*pi){ x = cos(t); y = sin(t)};
mesh Tｈ = buildmesh(C(10));
赤色：節点(vertex)番号
青色：要素(element)番号
円周上に 10 個
 x j,yj  , j  0,...,10
x j = cos  0, 2j  ,
yj= sin  0, 2j 
内部へ Delaunay 分割
領域の定義
領域  の境界  Jj1j： jは区分的に滑らかな曲線
j   x, y : x  1j,  t  , y   2,j  t  ,a j  t  bj
各j は反時計回りに向き付けられている。穴を開けるときは，
border a(t=0,pi*2){ x = 3*cos(t); y = 2*sin(t)}; //楕円
border b(t=0,2*pi){ x = cos(t); y = sin(t) }; // 内側円
mesh D2 = buildmesh(a(80) + b(-40)) ; // b の方向逆転
有限要素空間(Finite Element Space)
領域  ，その三角形分割 Th  Tk k 1,,M が与えられている。
P n  Tk  ： Tk 上の n 次多項式の造る集合
FreeFem++は命令


fespace Vh Th , P n ;
で n 次多項式連続要素を生成する。

 

Vh Th , P n  v  C0  h  : v  v11    vDD, v i  
ここで，各 Tk  Th に対し  i
Tk
 P n  Tk  で，添え字 Dは，次数 n
によって変化する自由度(Degree of Freedom)を表す。

Vh Th , P1
線形近似空間(1 次)
qk , k  1,, n v :節点で，nv を節点の総数。D= nv
 
要素関数は，  j qk  ik を満たす。
1.00
1.00
11
1
14
3
1.00
11
6
図 1 Th =buildmesh(C(10)) での  4 の形
要素関数の作成方法
要素関数の作成では，面積座標(重心座標)を使う。
k3
q
P
L2=0
A2
k1
q
A3
L3=0
A1
L1=0
k2
q
三角形要素 Tk 内点P  P  L1,L 2 ,L 3  の面積座標は，
L1  A1 A ,L 2  A 2 A , L 3  A 3 A
A, A1, A 2 , A3 はそれぞれ全体，３つに分割された三角形の面積
L i  x, y  a i  bi x  ci y
a i   x jyk  x k yj  D, bi   yj  yk  D,ci   x j  x k  D
ここで  x  , y  は点qk  ,   1, 2, 3 の座標で，  i, j, k  は 1, 2, 3  の偶置
換  i, j, k   1, 2, 3  ,  i, j, k    2, 3,1 ,  i, j, k    3,1, 2  を示し，
1 x1
D  2A  1 x 2
1 x3
y1
y2
y3
一次連続要素関数 Vh(Th,P1)
Tk
Tk
i
k n i 
q q
qi
1

 i  x, y 
L
x, y otherwise 0
i


k
q
T

n
i


k
2S

b
L
xc
y
 x, y  D a
k n i 
k

k n i 
k n i 
S は節点qi を含む三角形の面積の総和。
k n i 

2 次 Lagrange 要素 Vh(Th,P2)
k3
q
Tk
q
i
k4
qj qk 5
q
Tk
k2
q
qk1
k6
q


1


L
x
y
2L
x, y  1 otherwise 0
 i  x, y 
,
i




k
k
q
T

n
i
n
i
 
 
k
2S
2
  x, y L
  x, y otherwise 0
 j  x, y  qjT , j k,l L
k
l
k

S
S, Sはそれぞれ，節点qi ，中点qjを含む三角形の面積の総和。
不連続要素
Vh  Th ,P0  ：各三角形要素 Tk に対し，
k  x   1 if x  Tk ; k  x   x if x  Tk
Vh  Th ,P1nc  ：各三角形要素 Tk に対し，
m2
q
Tk
m3
q

1


,


x
y
D
L
x, y
i




k
k
m1 i
q
T

n
i
 
k
2S
q
otherwise 0
    ,m 
i m
j
ij
j
k n  j
q
k n  j 1
q
2

有限要素空間への射影
 : 関数 f  x, y の有限要素空間 Vh への射影
P0: f  x, y 
nt
 fkk  x, y, fk 
k 1
nv
P1nc: f  x, y 
P1: f  x, y 
P2: f  x, y 
 fii  x, y, fi
i 1
     
f qk1  f qk 2  f qk 3
3
 f m i 
 
nv
i
f

x
,
y
,
f

f
q


 i i
i
i 1
nv
nv
i 1
j1
 
 
i
i
f

x
,
y

f

x
,
y
,
f

f
q
,
f

f
m




 i i
j j
i
j
Poisson 方程式
弱形式
 u  v   fv
の左辺に未知関数 u h  u11    u n v  n v , u i   を，右辺に
 jを代入すると，連立方程式
nv
 u i  i  j   fj
i 1
 
を得る。 u h ql  0,ql   とするため
  i   j if i  j or i  j,qi  
 
A ji  

1 if i  j,qi  


u l    Fl   A iju j   0, Fl   fl if ql  



lj


この方法の利点は，非斉次問題 u  g on  でも


u l    Fl   A iju j   gl , Fl  Fl  1gl , gl  g ql


lj


 
のように簡単な修正で解くことが出来る。
varf a(u,v) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v))
+ on(,u=0);
matrix A = a(Vh,Vh);
// Stiffness matrix
varf b(v,f)= int2d(Th)(v*f);
matrix B = b(Vh,Vh); Vh F;
F[] = B*f[];
u[] = A^-1*F[];
フィルム(薄膜)のシミュレーション
border ball(t=0,2*pi){ x = 125*cos(t); y = 125*sin(t); };
border bottle(t=0,2*pi) { x=27.5*cos(t); y=27.5*sin(t); }
mesh Th = buildmesh(ball(50)+bottle(50));
fespace Vh(Th,P1);
fespace Ph(Th,P0);
Ph reg=region;
Ph f=(0.3/(pi*27.5^2))*(region==reg(0.,0.));
real G = 0.01;
Vh u,v;
problem antiPlane(u,v) =
int2d(Th)( G*(dx(u)*dx(v) + dy(u)*dy(v)) )
- int2d(Th)( -f*v )
+ on(ball,u=0) ;
antiPlane; // 問題はここで，解かれる。
plot (u,value=true,wait=1);
横せん断係数 G が 10gf/mm2 のとき，10mm の沈み。
境界条件 u=g とした剛性行列による解
Vh u,v,ff, bc,g=1;
varf a(u,v) = int2d(Th)( G*(dx(u)*dx(v) + dy(u)*dy(v)) )
+ on(1,u=1) ;
varf b(v,ff) = int2d(Th)(v*ff);
matrix A = a(Vh,Vh);
matrix B = b(Vh,Vh);
ff = -f;
Vh F; F[] = B*ff[];
bc[] = a(0,Vh);
// bc[]  bcj, bcj  1 if qj   else 0
F[] = F[] + bc[].*g[];
//
Fl  Fl  1gl if ql  
// bc[].*g[]  bcjgj 成分同士の積
u[] = A^-1*F[];
正則な三角形分割列
(1)
lim max diam  Tk  ; Tk  Th   0
h 0
ここで，diam  Tk  は，三角形の辺の最大長。
(2) ある定数   0 が存在し
  Tk 
diam  Tk 

Tk  Th ,h
ここで，   Tk  は三角形に内接する円の直径。
Cea の補題
a を強圧的連続双 1 次形式とする。このとき，
u  uh
 M
  1   inf u  vh
  vh Vh

(田端 p.54，次は p.62 系 2.5)


正則な三角形分割の列 Th に対し，有限要素空間 Vh Th ,P k ,
,
k  1  d2 を考える。このとき，任意の v  W k 12
   に対し， h に依
存しない定数c  0 が存在し，次が成り立つ。
(1.1) v   h v s,2,  ch k 1s v k 12
, ,
Poisson 方程式の 1 次要素近似
u  uh
s  0,1
12
, ,
 ch u
？
2,2,
 が凸領域なら， u  W 2,2    だが，凹んでいるときは駄目。
,
2,2

,
W
   , 0    1を
補完空間  W 12







,
12

v
W
 : v


 W 1,2 ,W 2,2  ,2,




12
  

   t K  t, v  dt 
 0

2


 

で定義する。このとき，
,
2,2
s,2
 W 12


,
W


W
  ,s  11     2







作用素Eh : u    u  u h  は，
Eh

,  ,L2 
L W012
   

 1,
Eh

,   W 2,2  ,L2 
L W012
 
   

 ch
よって，補完定理から
Eh
角 j 
L

W01,2
   W   ,L   
s,1
2
 j
 なら特異項 r
s 1
  ch 
1   j  ,2
W
,s  1
   ,0    1
となるので，Dirichlet 境界なら
u  uh
12
, ,
 ch 
 

f
0,2,
  max  j, 0    1
,
Poisson 方程式で f  u  u h ，弱解を w とすると(Aubin-Nitsche)
u  uh
0,2,
  w    u  u h 

    w  wh     u  u h 

 ch 
 

w  wh

ch
0,2,
 

f
0,2,
Adaptive mesh 法
u  uh
, ,
12
 u  h u
, ,
12


 C  hk u
k
12
 
2
2,2,Tk
三角形要素 Tk について，d k  minxTk x   j とすると
Tk  i  ju
2
2    j 1
d k 2   j  2 
r
rdr  Cd k
0
 C
となる。そこで， h k  d k のとき h k  C0 hd k
u  uh
12
, ,
1   j
とすれば，
 Ch
となる。これは，  k の近くでメッシュの細分になる。

Download Report