65
nhiệt sẽ kết hợp cả 3 điều kiện biên

3. áp dụng các phơng pháp xấp xỉ
d
dtije
, sẽ đợc hệ phơng trình
đại số của t
ije
,

k+1, với n ẩn số.
o
r
i
R
r
i
r
i
+
1
r
q
r

Z
l
-
1
Z

ije
k
+
d
dtije

k
,
đặt
ri

= i,
R

= , = và F =
2
a


, ta có hệ phơng trình đại số:
t
ijek+1
= F




)
2
i

) t
ij+1e
+ (1+
2
i

)t
i+1je
k




(ije) V
t
ijek+1
= F










4
i
1

2


-
4

-B)]t
Rje
+t
Rje+1
+
)
4
1(
2
2



t
Rj+1e
k



+
)
4
1(
BF2

5.8. FDM cho bài toán biên phi tuyến:
5.8.1. Điều kiện biên phi tuyến tính:
- Định nghĩa:
1. Phơng trình F(T
n
,
m
x
T
) = 0 có n 1 hoặc m 1 gọi là phơng
trình phi tuyến tính
2. Điều kiện biên đợc mô tả bởi một phơng trình vi phân phi
tuyến gọi là điều kiện biên phi tuyến
- Ví dụ: điều kiện biên loại 3, khi mặt vách tiếp xúc chất khí hoặc
chân không, trao đổi nhiệt với môi trờng chủ yếu bằng bức xạ, xác
định nhờ định luật Stefan-Boltzmann, thì phơng trình cân bằng nhiệt
trên biên có dạng:
-T
x
(W,) =
w

0
[T
4
(W,)-
4
f
T
]

T
x
(0,)=


0
[
4
f
T
-T
4
(0,)]
T
x
(0,) = 0
Nhận xét: bài toán nà
y có cùng
một mô hình TH

với bài toán biên
thuần nhất, khi vách dày 2 và TĐN
O
t
1 2 43
L
x
q = (T - T )
0



tạo ra 5 phần tử có 5 nút là i, (i = 0ữ4)
2. Sai phân theo x: phơng trình cân bằng nhiệt là:
- Trên biên bức xạ x
i
= 0:
cS
2
x
d2
dT
o
=
0
(
4
f
T
-
4
0
T
)S -
x


(T
0
-T
1

)
Nếu chuẩn hoá bằng cách đặt =
Ti
T

X =
L
x

X =
L
x

F =
2
Lc

=
2
L
a

Thì do

T
=

T
.
F

dF
d
4
0
4
f
i
2
00
i
2
0

i
2
40
2
aT)xL(
a
)TT(L2




(t)
hay
S
68
(
F


+2
4
f
3
i0
.
X
LT.




Nếu đặt R=


3
i0
LT.
là đại lợng không thứ nguyên, ta có:

0F
=
2
)X(
1

[-2(1+RX.
3
0

=
2
)X(
1

[
1
- 2
2
+
3
] V

3F
=
2
)X(
1

[
2
- 2
3
+
4
]

4F
=
2

0
,k+2p
1
,k+2pRX
4
f



i
,k+1= p
i-1k
+(1-2p)
ik
+p
1+1,k
với i = 1,2,3

4
,k+1= 2p
3
,k+(1-2p)
4
,k
( Với biên cách nhiệt
X
= 0, do đối xứng, coi
i-1
=
i+1

(1-2p) p


1

0

2
=
p (1-2p) p


2
+
0

3

p (1-2p) p


3

0

4
k+1

2p (1-2p)


...v. v

5.8.3. Bài toán truyền nhiệt qua vách phi tuyến:
Nếu thay ĐKB tại x = L của bài toán tại H.32 bởi điều kiện biên
loại 3, ta có bài toán sau:
T = aTxx
T(x,0) = T
i
T
x
(0,) =


0
[
4
1f
T
-T
4
(0,)]
T
x
(L,) =


[T(L,)-T
f2
]
Bài toán này đợc giải tơng

x
T
4

=
x

(T
3
-T
4
)-[T
4
-T
f2
] hay với B =

x
=

L
X

F
=
2
)X(
1

[2T


]
2p

0

2pR

X
4
f



1
p (1-2p) p

1
0

2
= p (1-2p) p

2
+ 0

3
p (1-2p)
p


0

]} > 0

tức là: 1-2p(1+ RX
3
0

> 0
1-2
2
)X(
F


(1+RX
3
0

) > 0 tức phải chọn;
F <
)](maxXR1[2
)X(
3
0
2
+

, với max
0

0, k
sau mỗi bớc.

71
Chơng 6
phơng pháp phần tử hữu hạn
finite element method (FEM)
6.1. Nội dung và các bớc của phơng pháp phân tử hữu hạn
6.1.1. Nội dung FEM.
T tởng của FEM là thay bài toán giải hệ phơng trình vi phân
(t) bằng một bài toán biến phân tơng ứng, tức là tìm hàm số t làm cực
tiểu một phiếm hàm I tơng ứng bài toán (t).
Bài toán biến phân đợc giải gần đúng nhờ phép xấp xỉ tích phân
bằng cách thay hàm t(x,y,z,) bởi một hệ M hàm thời gian t
n
() tại các
nút (đỉnh) của một số hữu hạn E phần tử tạo ra vật cần xét.
Kết quả cho biết, để cực tiểu biến hàm I, hàm t phải tìm cần thoả
mãn hệ M phơng trình vi phân thờng cấp 1 nh sau:
C
d
d
[t] = -(K+H)[t] + (h+q)
Trong đó C, K, H là các ma trận vuông (MxM) của các hệ số nhiệt
dung C, dẫn nhiệt , toả nhiệt , còn h, q là các ma trận cột (Mx1)
của các giá trị nhiệt độ môi trờng t
f
và dòng nhịêt q qua biên loại 2,
[t] =


I sẽ có dạng tích phân sau:
I[t(x,y,z)] =

v
f(x,y,z,t,t
x
,t
y
,t
z
) dV +

w
(qt +
2
t
2
)dw
Với W là biên của vật V, còn dạng hàm f xác định theo bài toán
(t) và định lý Euler-Lagrange về điều kiện cực tiểu phiếm hàm.
Điều kiện cực tiểu phiếm hàm I là biến phân I = 0.
O
t
e
3
k
i
S
W
t

j
y
i
x
i
x
j
x
k
x

H. 32 Phân bố nhiệt độ t
(e)
trong phần tử e hai chiều
2. Mô tả điều kiện cực tiểu I = 0 ở dạng hệ phơng trình vi phân
thờng cấp 1 của các nhiệt độ nút (hoặc hệ phơng trình đại số khi (t)
là bài toán ổn định). Bớc này có thể chia ra các bớc nhỏ nh sau:
2.1. Chia V (hoặc S, hoặc ) ra một số hữu hạn phần tử có dạng
khối tứ diện (hoặc tam giác, hoặc đoạn x) bởi hệ thống M điểm nút
(coi là đỉnh phần tử). Đánh số thứ tự và ghi địa chỉ nút theo toạ độ (i, j,
k) của mỗi nút.
2.2. Giả thiết rằng: Tại một thời điểm bất kỳ, phân bố nhiệt độ
t
(e)
trong phần tử e là một hàm tuyến tính của các nhiệt độ nút (tức có
dạng mặt phẳng qua 3 điểm t
i
, t
j
, t

M
) là phiếm hàm của M hàm nhiệt độ nút
t
i
() và điều kiện cực tiểu I = 0 trở thành
]t[d
dI
= 0,
với [t] =










M
1
t
..
t
tức
]t[d
dI
=

=

tuỳ ý.
- FEM rất tiện lợi khi hệ vật có biên dạng không quy tắc, vì khi đó
chỉ cần ghi điạ chỉ các nút biên, không cần tính thể tích và diện tích
mặt các phần tử nh trong FDM.
- Khi biên di động (bài toán biên loại 5), cả FDM, FEM sẽ trở nên
phức tạp.
6.2. Cực tiểu của hàm nhiều biến và phép xấp xỉ tích phân
Cơ sở của FEM là điều kiện cực tiểu của hàm số và phiếm hàm
cùng phép xấp xỉ tích phân.
6.2.1. Điều kiện cực tiểu hàm số u = u(x
1
, x
2
,...,x
n
)
- Theo phép tính vi phân, điều kiện cần và đủ để hàm u = u(x) đạt
74
cực tiểu tại x là:
dx
du
= 0 và
2
2
dx
ud
> 0
- Điều kiện cần và đủ để u = u(x,y) đạt cực tiểu tại (x,y) là:




Định nghĩa một ma trận cột [x]













n
..
2
1
x
x
x
, điều kiện trên có thể viết ở
dạng:
0
]x[
u
=








xn
u
...
x
u
x
u
2
1
=
















F(x,y,z) là tuyến tính với các biến tọa
độ (x
i
, y
j
, z
k
) tại 4 đỉnh của tứ diện.
- Nhờ cách chia và giả thiết nói
trên, ta có thể xác định tích phân I nh
tổng các giá trị trung bình của tích
phân trên mỗi phần tử I =

=
E
1e
)e(
I

- Ví dụ: Biểu thức xấp xỉ tích phân trên E = 5 phần tử ở hình H.33.
I =

M
1
x
x
dx)x(F
là I =

=

)/2 là trị trung bình của F(x) trên
phần tử (e) có kích thớc
e
x

= x
j
-x
i
Nếu chọn x
ij
x
j
- x
i
= , M = 4, E = 5 thì ta có:

+

=
M
1
x
x
e
ji
)FF(

[u
1
(x),u
2
(x)], ...v.v
- Phiếm hàm thờng là một tích phân của một hàm F nào đó trên
một miền xác định cho trớc.

=

O
t
(1)
x
1
2
(Fi + fj)
(2) (e) (E)
x
1
x
i
x
j
x
M
x

F(x)
H.33. Để xấp xỉ tích phân I


2
1
s
s
)s(T
ds.
+ Độ dài đờng cong y(x) bất kỳ qua 2 điểm (x
0
,y
0
) (x
1
,y
1
) là
phiếm hàm l = l[y(x)] =
dx'y1
2
1
x
x
2

+

6.3.2. Nội dung của lý thuyết biến phân:
- Lý thuyết biến phân là một ngành của toán học chuyên nghiên
cứu các phơng pháp tìm cực trị của các phiếm hàm. Nó cho phép tìm
đợc một hoặc một số hàm số làm cực tiểu một phiếm hàm đã cho.

A(x
0
,y
0
), B (x
1
,y
1
) cho trớc để chất
điểm lăn không ma sát trên nó từ A đến
B tốn ít thời gian nhất.
Theo định luật bảo toàn năng lợng, tìm thấy quan hệ = [y(x)]
H35. BT đờn
g đoản thời
O
y
x
x
o
x
1
y(x)
ds
dx
g
v
H.34. Bài toán tìm mặt cực tiểu

O
t

1
d
2






+

= [y(x)] =
dx
y
'y1
g2
1
1x
0x
2

+

3. Tìm đờng trắc địa:
Tìm đờng cong f(x,y,z) nối hai
điểm A,B trên mặt cong (x,y,z) = 0
cho trớc, có độ dài ngắn nhất. Theo
ý nghĩa hình học, bài toán (3) theo
biến phân là tìm cực tiểu phiếm hàm
l= l[y(x),z(x)] =

thay đổi tuỳ ý trong lớp hàm
u(x,) nào đó, gần đờng cong u(x).
Tơng tự định nghĩa:
u(x,y) =
),( yxu
- u(x,y)
- ý nghĩa hình học:
)x(u
thay đổi trong lớp đờng cong
qua hai điểm biên (x
0
,y
0
) (x
1
,y
1
) phụ thuộc thông số nhỏ nào đó:
H.36. BT tìm đờng
trắc địa



=
=
)x(zz
)x(yy

O
z

78
)x(u

);x(u
= u(x) + [
)(xu
-u(x)
= u(x) + u.
),( yxu
thay đổi trong lớp mặt cong qua
chu tuyến S, phụ thông số nhỏ nào đó:
),( yxu
u(x,y;) và u(x,y;) = u(x,y) +
[
)y,x(u
-u(x,y)] = u(x,y) + u(x,y).

Biến phân u là sự sai khác nhỏ giữa hai hàm u cùng chung biên
nào đó.
6.3.3.2. Biến phân của phiếm hàm:
- Định nghĩa: Biến phân I của phiếm hàm I[u(x)] là đạo hàm của
phiếm hàm I[u(x) + u(x)] theo tham số , tính tại = 0.

I[u(x)] =
d
d
I[u(x)] + u(x)]

=0


bảo toàn nào đó, so với định luật bảo toàn năng lợng.
6.3.4. Định lý Euler -Lagrange
6.3.4.1. Điều kiện cực tiểu phiếm hàm I[u(x)] =
dxuuxF
x
x

1
0
)',,(

- Cần tìm hàm u = u(x) qua 2 điểm biên cho trớc [x
0
, u
0
= u(x
0
)]
và [x
1
, u
1
= u(x
1
)] làm cực tiểu phiếm hàm có dạng:
I[u(x)] =

1
0
x

u(x) và u(x) là biến phân của hàm
u(x).
- Với đờng cong u(x;) thì phiếm hàm I[u(x,)] là một hàm chỉ
của , do đó điều kiện cần để cực tiểu phiếm hàm của hàm u(x), ứng
với = 0, là:
d
d
I[u(x,)]

=0
= 0 tức là I = 0
I = 0 tức là I =
d
d

1
0
x
x
F(x,u(x;), u'(x;))dx|

=0

=
+










=
0)''.(
1
0
=+

dxuFuuFu
x
x


Phân đoạn tích phân sau, ta có:
I =

1
0
x
x
Fu u dx + Fu' u
|
1
0
x
x
-
udxFu

Vì u (x) 0 x (x
0
, x
1
) nên Fu -
dx
d
Fu' = 0
Suy ra định lý Euler - Lagrange:
- Định lý E - L1:

Để cực tiểu phiếm hàm I[à(x)] =

1
0
x
x
F(x, u, u')dx, hàm u(x) cần
H.39. Để cực tiểu
I[u(x)]

O
x
0
x
1
x
u
(
x

y
-
'y
F
dx
d
= 0 suy ra:
F
y
- F
y'x
- F
y'y
.y'-F
y'y'
.y'' = 0 ,(E)
Phơng trình (E) gọi là phơng trình Euler.
- Phơng pháp tìm đờng cong cực trị là:
1) Tích phân phơng trình Euler thu đợc hàm y(x,c
1
, c
2
)
2) Xác định c
1
, c
2
theo hai điều kiện biên y
0
= y(x

=
dx
d
(F-F
y'
.y')= 0, (dạng đạo hàm toàn phần theo x)
Tích phân đầu theo x có F- F
y'
y' = C
1
là một phơng trình vi phân
không chứa x, giải đợc bằng cách tách biến hoặc đổi biến, ví dụ
* Tìm mặt cực tiểu (bài toán 1) :
Mặt tròn xoay qua (o,y
0
), (x
1
,y
1
)
cực tiểu khi phiếm hàm S[y(x)] cực
tiểu. Do dS = 2yds =
2y
22
dydx +
nên suy ra
S[y(x)] = 2
dxyy
x
x

y
1
y
o
ds
dx
81
F- y'F
y'
= y
2
'y1+
-
2
2
'y1
'yy
+
=
2
'y1
y
+
= C
1

Đổi biến y' = sht thì y = C
1
cht và dx =
==

Cx

là đờng dây xích, quay
quanh ox tạo ra mặt Catenoid, với C
1
, C
2
tìm theo điều kiện biên y
o
=
y(0) , y
1
= y(x
1
).
* Tìm đờng đoản thời: (bài toán 2)
Thời gian chất điểm lăn trên đờng cong y(x) từ AB là:
+
= [y(x)] =
g2
1
dx
y
y
x
x

=
+
1

= C
y(1 + y'
2
) = C
1
là phơng trình vi phân cấp 1 phi tuyến.
Đổi biến y' = cotg t có y =
tg
C
y
C
2
1
2
'
1
cot1
1
+
=
+

= C
1
sin
2
t =
2
C
1

v
d
=


y
d
x
B
y
1
y
0

A
x
1
82
= C
1
(1- cos2t)dt x = C
1
(t-
2
1
sin2t) + C
2

hay x =
2

1
2
1
qua A(0,0) C
2
= 0 đặt C =
2
C
1
,
2t = ta có đờng đoản thời là đờng cycloid, mô tả bởi phơng
trình dạng tham số sau:



=
=
)cos1(Cy
)sin(Cx
hay x = c[arcos(1-
c
y
) - sin(arcos(1-
c
y
))]
với C = bán kính vòng tròn lăn để vẽ ra Cycloid, xác định theo y
1
= y(x
1

u
(x,y) - u(x,y) ]
= u(x,y) + u(x,y)
thì phiếm hàm I[u(x,y,)] chỉ phụ thuộc , đợc cực tiểu khi:
d
d
I[u(x,y,)]

=0
= I = 0, tức là khi:
83
δI =
0
D
y
x
x
u
Fuy
u
Fu
u
Fu
=ε
∫∫
⎥
⎦
⎤
⎢
⎣

∂
∂
δ+δ
∂
∂
=δ
∂
∂
yy
y
y
xx
x
x
uFuu
y
Fu
)uFu(
y
uFuu
x
Fu
)uFu(
x

nªn
)(
yuyx
D
x

x
Fu
(
y
D
x
δ
∂
∂
+
∂
∂
∫∫

= -
udxdy)
y
Fu
x
Fu
(
y
D
x
δ
∂
∂
+
∂
∂

- V× vËy, ®iÒu kiÖn cùc tiÓu phiÕm hµm lµ:
udxdy)
y
Fu
x
Fu
Fu(I
y
x
D
δ
∂
∂
−
∂
∂
−=δ
∫∫
= 0, do δu|
D
≠ 0 suy ra:
→
0)
u
F
(
y
)
u
F

H.42. §Ó cùc tiÓu
phiÕm hµm I[µ(x,y)]

x
y
D
s
U
(
x
,
y
)
U
(
x
,
y
)
U
U
δ
u
u
u(x,y
u(x,y
W
0
84
u(x,y) cần thoả mãn điều kiện:


D
F
(x,y,z,t,t
x
,t
y
,t
z
)dV+

+
S
2
ds)
2
t
qt(

- Khi cực tiểu phiếm hàm I, cần có I = 0, tức là:
I =











dV+
dstttq
S

+ )(

= 0
- Theo t
x
=
)()( t
xx
t



=


, v.v., ta có:
I=









F
)t(
x
.
t
F
t
t
F
d
V+
sdtqt
S

+ )(

= 0
- Sau khi tích phân từng phần (per partes) ta
có:






V
x
dV)t(
x
.







V
y
dV)t(
y
.
t
F
=



S
y
y
lt
t
F
.

dS -





z
lt
t
F
.

dS -






V
z
dV)
t
F
.(
z
t

với l
x
, l
y
, l
z
là cosin chỉ phơng của pháp tuyến
n phía ngoài V tại M S. Do đó, điều kiện I = 0



























)()()(


+

Do t|
v
0 biểu thức trong dấu móc []
v
= 0 suy ra định lý
E-L3:
Định lý Euler -Lagrange N
0
3:
Để cực tiểu phiếm hàm I[t(x,y,z)] cho theo công thức (3), hàm
t(x,y,z) cần thoả mãn điều kiện:







=++


+


+


=



yt
F
xt
F
z
z
y
y
x
x
zyx
),,(,0
),,(,0)()()(


6.4. Ví dụ minh hoạ các bớc áp dụng FEM
6.4.1. Bài toán biên cô lập
Ta sẽ dùng FEM giải bài toán DN
không ổn định một chiều, biên cô lập
W
1
+ W
20
cho bởi mô hình.







(
xt
F
x
=







(phơng trình E-L)
H.44.Bài toán (t)

x
O
t = t
tx = 0
o
t(x,0) = t
i


c


t
=


=

t
c
t
F
và
x
x
t
t
F
=



Tích phân 2 phơng trình trên theo t và t
x
, coi t và t
x
là biến độc
lập, ta có:
F =



=




2
1
)t(g +

So sánh hai tích phân cho thấy hàm tích phân là:
F=








+


2
2
)(
2
1
x
tt
c




Bài toán biến phân tơng ứng với bài toán (t) là tìm hàm t = t(x,)

I = I

+I
c
với I

=
dx)
x
t
(
2
1
2
L
0




, I
c
=
dx
t
c
2
1
2
L

,...,t
i
,t
j
,...,t
M
)
H4
5. FE formulation
x
1
O
t
(1)
x
(2) (e) (E)
x
2
x
i
x
j
x
M
i
t
2
t
i
t

các phần tử hữu hạn, so với hàm t(x) chính xác cha biết thỏa mãn
phơng trình cân bằng nhiệt
2
2
x
tt
c






= 0 .
Tại thời điểm bất kỳ, phiếm hàm I là hàm của M nhiệt độ nút I =
I(t
1
, t
2
, ..., t
M
) nên điều kiện cần để cực tiểu I là sự triệt tiêu của các
đạo hàm I lần lợt theo mỗi nhiệt độ nút t
i.
=


i
t
I

t
...
t
t
(1), thì điều kiện cực tiểu I là:
=
][td
dI
][
td
dI

+
][td
dI
c
= 0, trong đó định nghĩa
][td
dI












e
tại mỗi phần tử e;
I =

=
E
1e
e
I
=

==

+
E
1e
e
c
E
1e
e
II
với
dx)
x
t
(
2
1
I

Chú ý: Các chỉ số e hiểu là "của phần tử e", không phải số mũ.
Định nghĩa ma trận (2x1) nhiệt độ nút của phân tử e: [t]
e







j
i
t
t
(2)
88
và ma trận (Mx2) định vị phần tử e là: D
e













=


























E
e
j






E
1e
.
00
..
10
01
..
00


















=
E
e
e
e
c
e
td
dI
D
1
][

Nh vậy, để cực tiểu I, phải xác định [t] sao cho:
]t[d
dI

=

=

E
1e
e
e
e
]t[d
dI
D

=
e
2
e
1
+
x= [1x]
e
2
1








p
T
[

]
e
, ở đây gọi

T
[1 x] (4),
các hệ số
e

1
j
i
x1
x1














P
e
[

]
e

Với P
e










11
1
ij
ij
xx
xx
(6) là ma trận
nghịch đảo của P
e

6.4.3.3. Tính
]t[d
dI

=
e
E
e
e
e1
dI
D
d[t]

x
t
(
2
1
I
=
2
)][(
2









j
i
x
x
eeT
e
tRp
x

dx
=

)1x2(
e
]t[
) theo ma trận [t]
e
:
ee
e
td
d
td
tAd
][][
)][(
=
([a
1
a
2
]
)
t
t
e
j
i










2
1
a
a
= A
T
Chú ý: Ký hiệu A
T
hiểu là chuyển vị của ma trận A, và quy tắc
chuyển vị tích (AB)
T
= B
T
A
T
(không giao hoán):
- Tính
e
e
]t[d
dI

=
2
x

)pR()]t[Rp(
x
e
x
x
eeT
x
e
T
j
i


dx=

j
i
T
x
x
eeT
xx
ee
tRppR
)][)((

dx, vì R
e
, [t]
e

.
1
0
1
1
1
1
















i
j
ij
ij
e
x
x









=
ee
]t[K

Với K
e










11
11
x
ij
e
(8) (đối xứng) là ma trận dẫn nhiệt của phần







M
j
i
t
t
t
t
1
=
]t[D
T
e
ta có:
e
e
E
1e
e
]t[d
dI
D
]t[d
dI