Bài giảng Phương pháp thể tích hữu hạn giải các bài toán

PHƯƠNG PHÁP THỂ TÍCH HỮU HẠN 
GIẢI CÁC BÀI TOÁN
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Böôùc moät: Taïo löôùi.
B A
W
P Caùc ñieåm nuùtThe åtích kieåm soaùt 
Bieân cuûa theå tích kieåm soaùt
E
N
T
W
B
S
E
P
t
n s
w
b
e
S J-
1
N
E
W ew
s
n
j
J
j-1
J+
1
I+
1
i+
1
IiI-1
Thể tích kiểm 
soát vô hướng 
(phương trình 
liên tục)
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Sai phân hóa
φφφρρφ Sgraddivudivt
+Γ=+
∂
∂ )()()(
dtdVSdtdA)grad.(ndtdA)u.(ndVdt)(
t
tt
t V
tt
t A
tt
t AV
tt
t
∫ ∫∫ ∫∫ ∫∫ ∫
∆+
φ
∆+∆+∆+






+





φΓ=





φρ+





ρφ
∂
∂
( ) ( ) dtV.Sdt)uA(dt)A(dVdt)(
t
tt
t
tt
t
tt
tV
tt
t
∫∫∫∫ ∫
∆+∆+∆+∆+
∆+φρ∇−φΓ∆=





ρφ
∂
∂
Tích phân theo thể tích hữu hạn rời rạc 
V).(dVdt)(
t
0
PP
V
tt
t
∆φ−φρ=





ρφ
∂
∂
∫ ∫
∆+
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Sai phân hóa
( ) ( )
[ ] [ ] [ ]( )∫
∫
∫∫
∆+
∆+
∆+∆+
φρ−φρ+φρ−φρ+φρ−φρ
−









∂
φ∂
Γ−
∂
φ∂
Γ+





∂
φ∂
Γ−
∂
φ∂
Γ+



∂
φ∂
Γ−
∂
φ∂
Γ
=φρ∇−φΓ∆
tt
t
btsnwe
tt
t
btsnwe
tt
t
tt
t
dt)uA()uA()uA()uA()uA()uA(
dt)
z
A()
z
A()
y
A()
y
A()
x
A()
x
A(
dt)uA(dt)A(
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Rời rạc hoá phương trình tích phân





 φ−φ
Γ−




 φ−φ
Γ=





∂
φ∂
Γ−





∂
φ∂
Γ
PW
wP
ww
PE
PE
ee
we x
A
x
A
x
A
x
A





 φ−φ
Γ−




 φ−φ
Γ=





∂
φ∂
Γ−





∂
φ∂
Γ
PN
SP
ss
PN
PN
nn
sn y
A
y
A
y
A
y
A





 φ−φ
Γ−




 φ−φ
Γ=





∂
φ∂
Γ−





∂
φ∂
Γ
PB
BP
bb
PT
PT
tt
bt z
A
z
A
z
A
z
A
Đặt:
F = Aρu; D = AΓ/xi,j
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Rời rạc hoá phương trình tích phân
[ ] [ ] [ ]{ }
( ) ( ) ( ) ( ) ( ) ( ){ }
dtV.S
dt)(D)(D)(D)(D)(D)(D
dtFFFFFFV).(
tt
t
tt
t
BPbPTtSPsPNnwPwPEe
tt
t
bbttssnnwwee
0
PP
∫
∫
∫
∆+
∆+
∆+
∆+
φ−φ−φ−φ+φ−φ−φ−φ+φ−φ−φ−φ
+φ−φ+φ−φ+φ−φ−=∆φ−φρ
(*)
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Rời rạc hoá phương trình tích phân
Để xác định vế phải của phương trình (*), tham số trọng 
lượng θ nằm trong khoảng từ 0 đến 1 sẽ được áp dụng. 
Các tích phân bên vế phải sẽ được viết lại như sau:
t])1(.[dtI 0PP
tt
t
P ∆φθ−+φθ=φ= ∫
∆+
φ (**)
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Rời rạc hoá phương trình tích phân
Sử dụng phương trình (**) for φE, φW, φN, φS, φT, φB vào 
phương trình (*) và chia phương trình này cho ∆t, ta được:
[ ] [ ] [ ]{ }
[ ] [ ] [ ]{ }
( ) ( ) ( ) ( ) ( ) ( ){ }
( ) ( ) ( ) ( ) ( ) ( ){ }
VS
)(D)(D)(D)(D)(D)(D)1(
)(D)(D)(D)(D)(D)(D
FFFFFF)1( 
FFFFFF
t
V).(
0
B
0
Pb
0
P
0
Tt
0
S
0
Ps
0
P
0
Nn
0
W
0
Pw
0
P
0
Ee
BPbPTtSPsPNnwPwPEe
0
bb
0
tt
0
ss
0
nn
0
ww
0
ee
bbttssnnwwee
0
PP
∆+
φ−φ−φ−φ+φ−φ−φ−φ+φ−φ−φ−φθ−
+φ−φ−φ−φ+φ−φ−φ−φ+φ−φ−φ−φθ
+φ−φ+φ−φ+φ−φθ−
−φ−φ+φ−φ+φ−φθ−=
∆
∆φ−φρ
(***)
Khi θ = 0, phương trình (***) trở nên tường minh, nếu 0<θ<1,
phương trình (***) không tường minh, còn nếu θ = 1, thì phương
trình (***) hoàn toàn không tường minh. Khi θ = 1/2, phương
trình (***) được gọi là phương trình Crank-Nicolson. Trong phần
này, phương pháp rời rạc hóa không tường minh hoàn toàn sẽ
được áp dụng để rời rạc hóa các phương trình tổng quát.
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Rời rạc hoá phương trình tích phân
u
0
P
0
PTTBBNNSSEEWWPP Saaaaaaaa +φ+φ+φ+φ+φ+φ+φ=φ
Bởi vì phương pháp này áp dụng cho quá trình thay đổi tức thời
(transient), nên người ta sử dụng các phương trình khuếch tán-đối
lưu và các sơ đồ chuyển đổi qua lại. Do đó, ta có:
Trong đó:
P
0
PTBNSEWP SFaaaaaaaa −∆+++++++=
t
Va
0
P0
P ∆
∆ρ
=
Với:
PPu SSV.S φ+=∆
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Rời rạc hoá phương trình tích phân











 += 0,
2
F
D,Fmaxa wwwW 










 −−= 0,
2
F
D,Fmaxa eeeE











 += 0,
2
F
D,Fmaxa sssS 










 −−= 0,
2
FD,Fmaxa nnnN











 += 0,
2
F
D,Fmaxa bbbB











 −−= 0,
2
F
D,Fmaxa tttT
∆F = Fe – Fw + Fn – Fs + Ft – Fb
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Thuật toán ma trận ba đường chéo TDMA
φ1 = C1 
-β2φ1 + D2φ2 - α2φ3 = C2 
 -β3φ2 + D3φ3 - αP1φP1 = C3 
 .. 
 -βnφn-1 + Dnφn - αnφn+1= Cn 
 φn+1 = Cn+1 
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Thuật toán ma trận ba đường chéo TDMA
Trong các phương trình trên, φ1 và φn+1 được xem là những giá trị
biên. Phương trình dạng tổng quát được viết như sau:
-βjφj-1 + Djφj - αjφj+1 = Cj
2
2
1
2
2
3
2
2
2 D
C
DD
+φ
β
+φ
α
=φ 
3
3
2
3
3
4
3
3
3 D
C
DD
+φ
β
+φ
α
=φ 
 .. 
n
n
1n
n
n
1n
n
n
n D
C
DD
+φ
β
+φ
α
=φ −+ 
CuuDuongThanCong.com https://fb.com/tailieudientucntt
BÀI TẬP 
MÔ HÌNH HÓA VÀ MÔ PHỎNG 
SỬ DỤNG 
PHƯƠNG PHÁP THỂ TÍCH HỮU HẠN
CuuDuongThanCong.com https://fb.com/tailieudientucntt
1. Bài toán truyền dẫn nhiệt
Quá trình truyền nhiệt ổn định qua tấm phẳng có bề dày L = 2cm;
hệ số dẫn nhiệt k = 0,5 (W/m.độ); nguồn cấp nhiệt không đổi
q/ρCp = 1000 (kW/m3). Các bề mặt A và B của tấm phẳng có nhiệt
độ 100oC và 200oC. Giả sử rằng kích thước theo phương y và z đủ
lớn để gradient nhiệt độ chỉ xuất hiện theo phương x.
Xác định:
1) Phương trình mô tả quá trình truyền nhiệt
2) Biết rằng: Nghiệm giải tích của bài toán trên có dạng:
A
AB TxxL
k
q
L
TTT +


 −+
−
= .)(
2
Đánh giá kết quả gần đúng theo kết quả mô phỏng và
kết quả chính xác theo nghiệm giải tích trên (5 nút)
CuuDuongThanCong.com https://fb.com/tailieudientucntt
φφφρρφ Sgraddivudivt
+Γ=+
∂
∂ )()()(
Xuất phát từ phương trình tổng quát
Phương trình truyền nhiệt dạng tổng quát được viết 
dưới dạng:
Q
z
Tk
y
Tk
x
Tk
z
Tu
y
Tu
x
TuC
t
TC zyxzyxpp +





∂
∂
+
∂
∂
+
∂
∂
=





∂
∂
+
∂
∂
+
∂
∂
+
∂
∂
2
2
2
2
2
2
ρρ
Q
x
Tk
t
TC xp +∂
∂
=
∂
∂
2
2
ρ
Phương trình truyền nhiệt trong tấm phẳng có kích 
thước chiều y>>x và z >>x có dạng:
Trong trường hợp truyền nhiệt ổn định:
02
2
=+
∂
∂ Q
x
Tkx
CuuDuongThanCong.com https://fb.com/tailieudientucntt
uPPEEWWPP STaTaTaTa +++=
00
B
A
W
P Caùc ñieåm nuùt
Bieân cuûa theå tích kieåm soaùt
E
∆x ∆x∆x/2
w e
0
)()(
=+





∫∫
VV
qdVdV
dx
dTk
dx
d
0)( =+




−




 Vq
dx
dTAk
dx
dTAk
we
0=∆+




 −
−




 − xqA
x
TTkA
x
TTkA
PW
wP
ww
PE
PE
ee
Điểm giữa
CuuDuongThanCong.com https://fb.com/tailieudientucntt
w
w
wwW D
FDFa =










 += 0,
2
,max e
e
eeE D
FDFa =










 −−= 0,
2
,max
PEWP Saaa −+=
F = 0 D = Ak/∆x
Sp = 0 SU = qA∆x
CuuDuongThanCong.com https://fb.com/tailieudientucntt
uPPEEWWPP STaTaTaTa +++=
00
B
A
W
P Caùc ñieåm nuùt
Bieân cuûa theå tích kieåm soaùt
E
∆x ∆x∆x/2
w e
0
)()(
=+





∫∫
VV
qdVdV
dx
dTk
dx
d
0)( =+




−




 Vq
dx
dTAk
dx
dTAk
we
0
2/
=∆+




 −
−




 − xqA
x
TTkA
x
TTkA
PW
BP
ww
PE
PE
ee
Điểm 1
CuuDuongThanCong.com https://fb.com/tailieudientucntt
00,
2
,max =










 += wwwW
FDFa e
e
eeE D
FDFa =










 −−= 0,
2
,max
PEWP Saaa −+=
F = 0 D = Ak/∆x
SU = qA∆x+2kATB/∆x
x
kASP ∆
−=
2
Làm tương tự đối với điểm 5
CuuDuongThanCong.com https://fb.com/tailieudientucntt
wW Da = 0=Ea
PEWP Saaa −+=
F = 0 D = Ak/∆x
SU = qA∆x+2kATA/∆x
x
kASP ∆
−=
2
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Lập bảng tính
Điểm aW aE SU SP aP = aW + aE - SP
1 0 k/∆x q∆x+2kTB/∆x 3k/∆x
2 k/∆x k/∆x q∆x 0 2k/∆x
3 k/∆x k/∆x q∆x 0 2k/∆x
4 k/∆x k/∆x q∆x 0 2k/∆x
5 k/∆x 0 q∆x+2kTA/∆x 3k/∆x
x
k
∆
−
2
x
k
∆
−
2
Kết quả tính
Điểm aW aE SU SP aP = aW + aE - SP
1 0 125 29000 -250 375
2 125 125 4000 0 250
3 125 125 4000 0 250
4 125 125 4000 0 250
5 125 0 54000 -250 375
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Lập ma trận:
375 -125 0 0 0
-125 250 -1250 0 0
0 -125 250 -1250 0
0 0 -125 250 -1250
0 0 0 -125 375
T1
T2
T3
T4
T5
=
29000
4000
4000
4000
54000
T1
T2
T3
T4
T5
=
150
218
254
258
230
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Điểm 
nút
1 2 3 4 5
X(m) 0.002 0.006 0.01 0.014 0.018
Kết quả 
gần 
đúng
150 218 254 258 230
Kết quả 
chính 
xác
146 214 250 254 226
Sai số 2.73% 1.86% 1.6% 1.57% 1.76%
CuuDuongThanCong.com https://fb.com/tailieudientucntt
2. Bài toán lan truyền chất
Quá trình lan truyền ổn định một chất khí màu trong không khí,
biết: khối lượng riêng lưu chất tại điều kiện lan truyền là ρ =
1kg/m3; vận tốc đối lưu theo phương x là ux = 0.1m/s; hệ số khuếch
tán D = 0.1 m2/s; Chiều dài quãng đường lan truyền là L = 1m. Giả
sử rằng quá trình lan truyền chỉ xảy ra theo một phương x.
Xác định:
1) Phương trình mô tả quá trình truyền vật chất
2) Biết rằng: Nghiệm giải tích của bài toán trên có dạng:
Đánh giá kết quả gần đúng theo kết quả mô phỏng và
kết quả chính xác theo nghiệm giải tích trên (5 nút)
71828.1
)exp(71828.2)( xx −=Φ
CuuDuongThanCong.com https://fb.com/tailieudientucntt
Nút Khoảng 
cách
Kết quả 
mô phỏng
Kết quả 
giải tích
Sai số
1 0.1 0.9421 0.9387 0.36%
2 0.3 0.8006 0.7963 0.53%
3 0.5 0.6276 0.6224 0.83%
4 0.7 0.4163 0.4100 1.53%
5 0.9 0.1579 0.1505 4.91%
CuuDuongThanCong.com https://fb.com/tailieudientucntt
3. Bài toán 2D
Cho tấm phẳng có bề dày 1cm. Độ dẫn nhiệt của vật liệu tấm
phẳng là k = 1000W/mK. Mặt biên phía Tây tiếp nhận thông lượng
nhiệt ổn định 500kW/m2 và các mặt biên phía Đông và Nam được
cách nhiệt.mặt biên phía Bắc được duy trì ở nhiệt độ 100oC, sử
dụng lưới không đổi ∆x=∆y=0,1m để tính toán sự phân bố tại các
điểm nút.
Xác định:
1) Phương trình mô tả quá trình truyền nhiệt trong tấm
phẳng
2) Mô phỏng gradient nhiệt độ trong tấm phẳng nêu trên với
chiều x = 0,3m; y = 0,4m
CuuDuongThanCong.com https://fb.com/tailieudientucntt