Phương pháp phần tử hữu hạn tự thích ứng và ứng dụng trong phân tích đập bê tông trọng lực

 PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN TỰ THÍCH ỨNG VÀ ỨNG DỤNG 
TRONG PHÂN TÍCH ĐẬP BÊ TÔNG TRỌNG LỰC 
ThS. VŨ HOÀNG HƯNG(1,2), TS. NGUYỄN QUANG HÙNG(1) 
(1) Trường Đại học Thủy lợi 
(2) NCS Học viện Công trình Thủy lợi thủy điện, Đại học Hà hải, Nam kinh, Trung quốc 
Tóm tắt: Phương pháp phần tử hữu hạn (FEM) là một phương pháp số có hiệu quả cao trong giải các 
bài toán cơ học công trình. Tuy nhiên đây là một phương pháp gần đúng nên vấn đề khống chế sai số 
nâng cao độ chính xác lời giải bài toán cũng như đánh giá độ hội tụ bài toán phân tích ứng suất biến 
dạng công trình là một vấn đề mang tính khoa học cao. Bài báo này tập trung đi sâu nghiên cứu vào kỹ 
thuật phi tuyến tính hình học trong việc áp dụng FEM vào phân tích ứng suất đập bê tông trọng lực 
trên nền đá để từ đó hoàn thiện hơn công nghệ tính toán đập bê tông trọng lực. Những kết quả nghiên 
cứu này có ý nghĩa khoa học và thực tế nhằm giải quyết vấn đề ứng suất cục bộ tại những vùng đặc 
biệt của đập bê tông trọng lực đang là vấn đề còn tồn tại của công tác thiết kế đập bê tông trọng lực 
của Việt Nam. 
Từ khóa : Phần tử hữu hạn, ứng suất, chuyển vị, tự thích ứng, sai số 
1 ĐẶT VẤN ĐỀ 
Phương pháp phần tử hữu hạn (Finite Element Method - FEM) là một phương pháp gần đúng để 
giải các bài toán cơ học, nghiệm của bài toán thu được dựa trên giải phương trình cân bằng của toàn 
miền tính toán theo nguyên lý cân bằng cực hạn. Chính vì vậy các đặc trưng cơ học cũng như hình học 
ảnh hưởng rất nhiều đến độ chính xác của nghiệm cũng như sự hội tụ của lời giải. Trong những tính 
toán giải bài toán FEM thông thường, ảnh hưởng chủ quan của việc chia lưới phần tử sẽ dẫn đến sự 
khác biệt rất lớn trong lời giải. Vì vậy theo những kinh nghiệm thông thường khi giải bài toán có miền 
tính toán phức tạp người ta thường chia nhỏ lưới phần tử tại những vùng có khả năng phát sinh ứng 
suất cục bộ cũng như những vùng có các đặc trưng đặc biệt. Tuy nhiên sự thay đổi lưới phần tử mang 
tính chủ quan này có ảnh hưởng như thế nào đến kết quả bài toán? Lời giải toán học có đạt được ý 
nghĩa cơ học hay không? Để trả lời cho câu hỏi này, bài báo tiến hành nghiên cứu phân tích một số 
nhân tố ảnh hưởng của việc chia lưới phần tử tới trường ứng suất biến dạng của đập bê tông trọng lực 
trên nền đá. Từ đó có những kết quả so sánh bước đầu trong trường hợp sử dụng thuật toán tự thích 
ứng trong FEM với những tính toán FEM thông thường nhằm nâng cao độ chính xác của lời giải bài 
toán cơ học cũng như đánh giá mức độ hội tụ của nó. Việc tính toán so sánh phân tích được tiến hành 
dựa trên 5 trường hợp khác nhau như hình 1 và hình 5. 
 Lưới PT dạng 1 (a) Lưới PT dạng 2 (b) Lưới PT dạng 3 (c) Lưới PT dạng 4(d) 
Hình 1 Bốn trường hợp phân chia kích thước mạng lưới phần tử hữu hạn không giống nhau 
 Bảng 1 Kết quả phân tích ứng suất chính nhỏ nhất tại gót đập 
Số trường hợp Mạng lưới tính toán Tải trọng tính toán Ứng suất chính min 
 tại gót đập (MPa) 
1 Lưới PT dạng 1 TLBT + ALNTL 10.96 
2 Lưới PT dạng 2 TLBT + ALNTL 8.22 
3 Lưới PT dạng 3 TLBT + ALNTL 5.79 
4 Lưới PT dạng 4 TLBT + ALNTL 5.20 
(Giá trị lớn nhất – Giá trị nhỏ nhất)/Giá trị lớn nhất 52.6% 
 Trong đó: TLBT : Trọng lượng bản thân. 
 ALNTL: Mực nước thượng lưu là mực nước dâng bình thường. 
Từ bảng 1 có thể thấy rằng cách chia lưới phần tử cũng như kích thước của nó có ảnh hưởng rất 
lớn đến giá trị ứng suất tập trung tại gót đập, giá trị ứng suất tại gót đập trong các trường hợp chia lưới 
khác nhau có thể sai khác nhau đến 52.6%. Vì vậy khi sử dụng FEM phân tích ứng suất – biến dạng 
nếu như chỉ dựa vào một lần lựa chọn kích thước lưới phần tử để tính toán thì khó có thể thu được kết 
quả ứng suất như mong muốn. Chính vì vậy, để có thể giải bài toán sát với thực tế nhất cần phải có 
một quá trình phân tích hiệu chỉnh để đánh giá tính chính xác của lời giải bài toán FEM. Trong quá 
trình tính toán hiệu chỉnh này, sử dụng thuật toán lưới phần tử tự thích ứng nhằm khống chế và phân 
phối sai số trong quá trình giải. 
2 PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN TỰ THÍCH ỨNG 
2.1 Đánh giá sai số phần tử hữu hạn tự thích ứng 
Lời giải của FEM là lời giải gần đúng, giá trị chênh lệch của nó với giá trị giải chính xác được gọi 
là sai số. Sai số chuyển vị được định nghĩa như sau: 
uˆueu (1) 
 Trong đó: u là giá trị chuyển vị giải chính xác; uˆ là giá trị chuyển vị giải gần đúng 
Sai số ứng suất được định nghĩa là: 
   ˆe (2) 
Trong đó:  là giá trị ứng suất giải chính xác; ˆ là giá trị ứng suất giải gần đúng 
Tuy nhiên trong thực tế thường có hiện tượng tập trung ứng suất tại một số vùng không đồng nhất 
về đặc tính cơ học cũng như hình học. Do vậy nên khi sử dụng các sai số được xác định theo (1) và (2) 
sẽ có sự chênh lệch rất lớn với thực tế. Nói cách khác, lời giải của bài toán đối với những vùng cục bộ 
này sẽ không còn ý nghĩa nữa. Điều này hoàn toàn phù hợp với những kết quả nghiên cứu hiện nay 
tương tự như các bài toán thiết kế, khi sử dụng các trị số ứng suất biến dạng cục bộ, người thiết kế 
thường phải lấy các giá trị này tại vùng lân cận để đưa vào thiết kế. Lưới phần tử tự thích ứng sẽ khống 
chế sự sai khác này và tiến hành phân phối sai số năng lượng toàn phần của toàn miền đảm bảo sai số 
tại các phần tử là như nhau. Điều này hết sức có ý nghĩa về mặt cơ học cũng như toán học. 
Định nghĩa sai số phạm vi năng lượng là: 
2/1
1T dˆDˆe
   

 (3) 
Trong đó: D là ma trận hằng số đàn hồi;  là miền tích phân. 
Thông qua giải phương trình cơ bản của FEM có thể thu được miền nghiệm gần đúng { uˆ } và 
{ ˆ }. Tuy nhiên việc xác định miền nghiệm đúng {u} và {} đối với các bài toán cơ học có biên phức 
 tạp là hết sức khó khăn. Do vậy có thể dùng trường ứng suất cải tiến {*} và trường chuyển vị cải tiến 
{u*} thay thế cho miền nghiệm đúng {}, {u} để phân tích 
[2,3]. Khi sử dụng miền nghiệm cải tiến, sai 
số năng lượng toàn miền được xác định như sau : 
2/1
*1T** dˆDˆe
   

 (4) 
Sai số năng lượng đối với phần tử thứ i được xác định theo công thức (5) 
2/1
i
*1T*
)i(
* dˆDˆe
   

 (5) 
Mức độ hiệu quả  dùng để đánh giá sai số được xác định theo công thức (6) 
e
ek *
  (6) 
Trong đó: k là hệ số hiệu chỉnh phụ thuộc vào loại hình phần tử 
Sai số toàn phần của toàn miền tính toán η được xác định theo công thức (7) 
2/1N
1e
2*T
*
adapt
edˆˆ
e
 
 
 
 
 (7) 
Trong đó: Nadapt là số phần tử cần gia tăng cục bộ trong lưới phần tử tự thích ứng. 
Việc đánh giá tính chính xác của lời giải bài toán FEM tự thích ứng được khống chế bằng sai số định 
trước f. Trong trường hợp sai số toàn phần của toàn miền  xác định theo (8) không thỏa mãn điều 
kiện <f thì lưới phần tử cần thiết phải được hiệu chỉnh bằng cách gia tăng số lượng phần tử. Số lượng 
phần tử cần phải bổ sung được xác định theo điều kiện (8): 
adapt
*
)i(
*
N
e
e (8) 
Nếu điều kiện (8) được thỏa mãn phần tử i cần tiến hành nâng cấp. Sự nâng cấp này được tiến 
hành theo các mô hình khác nhau. Tuy nhiên mỗi lần nâng cấp này, lưới phần tử có sự thay đổi và thiết 
lập lại hệ phương trình cơ bản của FEM. Việc giải bài toán FEM được tiến hành theo chu trình giải lặp 
đến khi thỏa mãn điều kiện khống chế sai số  f, lưới phần tử thu được cuối cùng gọi là lưới phần tử 
tự thích ứng. 
2.2 Phương pháp nâng cấp mạng lưới tự thích ứng 
Hiện nay phương pháp FEM tự thích ứng được tiến hành với bai loại mô hình cơ bản: Mô hình h 
(hình 2a) dựa trên nguyên lý chia nhỏ lưới phần tử nhằm nâng cao độ chính xác. Mô hình h phù hợp 
với những miền tính toán tương đối đồng nhất và hàm xấp xỉ gần đúng nhỏ có bậc nhỏ hơn 4. Mô hình 
p (hình 2b) dựa trên cơ sở tăng bậc của hàm gần đúng . Mô hình này tương đối phức tạp đối với những 
hàm gần đúng bậc cao, ma trận cứng phần tử sẽ biến đổi phức tạp và độ lớn của bài toán tăng lên rất 
nhiều. Mô hình r (hình 2c) dựa trên cơ sở phân chia lại kích thước phần tử, kích thước của các phần tử 
được thay đổi không đồng đều nhau nhằm phân tán sai số, nâng cao độ chính xác giải. Ngoài ra có thể 
kết hợp các mô hình để thành tổ hợp các mô hình khác nhau trong bài toán FEM như mô hình (h-p) , 
mô hình (p-r) [4,5]. 
a/ Phương pháp h b/ Phương pháp p c/ Phương pháp r 
Hình 2: Ba phương pháp thay đổi mạng lưới phần tử liên tục 
Trong các mô hình trên, lý luận phương pháp tự thích ứng mô hình h là khá phổ biến, ưu điểm của 
nó là phần mềm phân tích phần tử hữu hạn cơ bản không thay đổi nhằm duy trì tính độc lập, tính tự 
thích ứng của phương pháp khá mạnh, tính linh hoạt của chương trình khá tốt, nhưng nhược điểm là 
dựa vào một bộ phận sinh thành mạng lưới lớn mạnh, bộ phận sinh thành mạng lưới lớn mạnh này yêu 
cầu có khả năng thích ứng các phạm vi phức tạp, tính năng trạng thái của mạng lưới cần tốt, đối với 
ứng dụng của phương pháp mô hình h này yêu cầu khá cao. Đồng thời tốc độ hội tụ của nó khá chậm, 
độ chính xác tương đối thấp, đối với vấn đề điểm kỳ dị hiệu quả khá kém. Phương pháp mô hình p rất 
thích hợp xử lý như là vấn đề tập trung ứng suất đoạn cuối vết nứt, nhưng kết cấu phần mềm phân tích 
phần tử hữu hạn phức tạp. Phương pháp tự thích ứng h – p là loại phương pháp tiên tiến nhất, nó có 
khả năng cung cấp tốc độ hội tụ theo luỹ thừa, phỏng thật máy tính và phân tích vấn đề công trình thực 
tế này rất có ý nghĩa, song phương pháp này cũng là phương pháp khó nhất trong ứng dụng, có rất 
nhiều vấn đề đang trong quá trình nghiên cứu[6]. Hình 3 biểu thị một quá trình chỉnh thể tính toán phần 
tử hữu hạn tự thích ứng [7]. Phần mềm phân tích phần tử hữu hạn thông dụng ANSYS có khả năng 
phân tích tự thích ứng hiệu quả cao. Vấn đề phân tích tự thích ứng có thể tham khảo trong các tài liệu 
hướng dẫn sử dụng phần mềm. 
Hình 3: Quá trình tính toán phần tử hữu hạn tự thích ứng 
3 ỨNG DỤNG FEM TỰ THÍCH ỨNG PHÂN TÍCH ỨNG SUẤT BIẾN DẠNG ĐẬP BÊ TÔNG 
TRỌNG LỰC 
3.1 Mặt cắt tính toán 
Mặt cắt đập được dùng điển hình với chiều cao đập 137m, chiều rộng đáy đập 80m, bề rộng đỉnh 
đập 7m, phạm vi tính toán nền đập về phía thượng lưu 150m, về phía hạ lưu 250m, về phía dưới 250m 
(xem hình vẽ 4). Tải trọng tính toán gồm tải trọng do trọng lượng bản thân đập và áp lực nước thượng 
lưu gây ra. Trục X theo phương ngang hướng từ thượng lưu về hạ lưu, trục Y theo phương thẳng đứng 
hướng lên trên. Hai biên nền đập được khống chế chuyển vị theo phương ngang (phương trục X), biên 
đáy nền đập khống chế chuyển vị theo phương đứng (phương trục Y). Chỉ tiêu cơ lý của bê tông đập 
và nền được cho ở bảng 2. 
Cấu kiện công trình Dừng tính toán 
Lựa chọn mạng 
lưới cơ bản 
Tính toán PTHH Phân tích sai số 
Thực thi ưu hóa mạng lưới 
Thỏa mãn 
độ chính xác 
dự định? 
Có 
Không 
 Bảng 2 Chỉ tiêu cơ lý vật liệu đập và nền dùng trong tính toán 
Loại vật liệu Lực dính 
c (Mpa) 
Góc ma sát 
trong (độ) 
Modun đàn 
hồi E (Mpa) 
Hệ số 
Poisson  
Khối lượng 
riêng (kg/m3) 
Bê tông đập - - 21000 0.167 2400 
Đá nền - - 20000 0.20 2500 
E=21GPa
 
  kg/m3
E=20GPa
kg/m3
 
 
Hình 4 Điều kiện biên dùng trong tính toán 
3.2 Lưới phần tử FEM tự thích ứng 
Sử dụng phương pháp FEM tự thích ứng được sử dụng trong phân tích thông qua áp dụng mô hình 
h có sẵn trong phần mềm ANSYS để phân tích ứng suất biến dạng cho đập bê tông và nền. Giới hạn số 
bước giải là 10 bước (Number of Solutions = 10), sai số cho phép lựa chọn tính toán sử dụng trong 
trường hợp này là 5% (Percent error in energy norm = 5), hệ số nhỏ nhất cho thay đổi kích thước phần 
tử là 0.5 (Minimum element size factor = 0.5), hệ số lớn nhất cho thay đổi kích thước phần tử là 2 
(Maximum element size factor = 2). 
Mệnh lệnh chia lưới phần tử tự thích ứng và giải một mô hình: 
ADAPT,10,,5,0.5,2 
Hình 5 thể hiện lưới phần tử tự thích ứng và sai số toàn phần tương ứng ở các bước giải 2, 4, 6, 8. 
Bước tự thích ứng 2 
Số phần tử: 164 
Sai số:  = 16.1% 
Bước tự thích ứng 4 
Số phần tử: 748 
Sai số:  = 8.4% 
Bước tự thích ứng 6 
Số phần tử: 3408 
Sai số:  = 8.1% 
Bước tự thích ứng 8 
Số phần tử: 18219 
Sai số:  = 4.3% 
Hình 5 Lưới phần tử tự thích ứng và sai số toàn phần tương ứng ở các bước giải 2, 4, 6, 8 
3.3 Phân tích kết quả tính toán 
 Bảng 3 Kết quả phân tích ứng suất chính nhỏ nhất 1 tại gót đập (Mpa) 
 Bước 
giải 
Lưới PT 
2 
( = 16.1%) 
4 
( = 8.4%) 
6 
( = 8.1%) 
8 (cuối cùng) 
( = 4.3%) 
Lưới phần tử 1 10.96 
Lưới phần tử 2 8.22 
Lưới phần tử 3 5.79 
Lưới phần tử 4 5.20 
Lưới tự thích ứng 2.41 6.10 13.06 26.25 
Kết quả tính toán cho thấy giá trị ứng suất kéo tại gót thượng lưu đập theo phương thẳng đứng 
tăng rất lớn đạt tới 26.25 Mpa, tăng 139.5% so với ứng suất kéo lớn nhất trong bốn trường hợp chia 
lưới thông thường theo chủ quan của người tính toán. Hình 6 thể hiện hình ảnh của phân bố ứng suất 
chính nhỏ nhất 1 tại đáy đập ứng với các lưới phần tử khác nhau. 
0 10 20 30 40 50 60 70 80
-4
0
4
8
12
16
20
24
28
1
s
t 
P
ri
n
c
ip
a
l s
tr
e
ss
 (
M
P
a
)
Bottom of dam (m)
 element grid 1
 element grid 2
 element grid 3
 element grid 4
 adaptive element grid
Hình 6 Phân bố ứng suất chính nhỏ nhất 1 tại đáy đập ứng với các lưới phần tử khác nhau 
Từ hình (6) thấy rõ ràng rằng kết quả ứng suất 1 tại chân thượng lưu đập là nơi phát sinh ứng suất cục 
bộ lớn nhất trong trường hợp hồ đầy nước thể hiện rõ sự khác biệt lớn giữa các cách chia lưới phần tử 
khác nhau. Càng về hạ lưu đập, nơi ít xuất hiện ứng suất cục bộ thì sự sai khác này giảm đi đáng kể. 
Chính kết quả này đã thể hiện rõ được hiệu quả của việc sử dụng lưới phần tử tự thích ứng trong bài 
toán phân tích ứng suất biến dạng đập bê tông trọng lực trên nền đá khi có sử dụng phương pháp phần 
tử hữu hạn. 
4 KẾT LUẬN 
Từ những kết quả nghiên cứu đã cho thấy rõ rằng kích thước lưới phần tử hữu hạn ảnh hưởng rất 
lớn đến các giá trị ứng suất tập trung tại gót thượng lưu đập. Khi chia lưới phần tử không giống nhau, 
các giá trị ứng suất tại gót đập lớn nhất và nhỏ nhất sai khác nhau rất nhiều , đạt tới 54.4%. Do vậy nên 
khi sử dụng phương pháp phần tử hữu hạn để phân tích ứng suất – biến dạng đập bê tông trọng lực trên 
nền đá mà chỉ dựa vào một lần lựa chọn kích thước lưới phần tử để tính toán thì không thể thu được 
kết quả ứng suất như mong muốn mà cần phải có một quá trình điều chỉnh lưới phần tử để nâng cao độ 
chính xác của lời giải bài toán. Việc sử dụng kỹ thuật phi tuyến hình học của lới phần tử mà bài báo đã 
đề cập bằng thuật toán lưới phần tử tự thích ứng đã phần nào mô phỏng quá trình làm việc của đập sát 
với thực tế hơn nhằm thu được kết quả phân tích chính xác nhất. 
 Tài liệu tham khảo 
1 Liu Yi. Research on some factors on working performance of high arch dam. Post-doctor 
Dessertation. China Institute of Water Resources and Hydropower Research, 2007. 
2 Zhu Bofang. The Finite Element Method Theory and Applications. China WaterPower Press, 
2004, p590~591. 
3 Yang Qiang, Wu Hao, Zhou Weiyuan. Study on determination of dam stress based on 
h-version adaptive FEM. Journal of Hydraulic Engineering, 2005, Vol.36 No.3, p321~327. 
4 Yang Ling-qiang, MA Jing, CHEN Zu-ping. Study on stress control standards of high arch 
dam with self-adaptive finite element method. Water Resources and Hydropower Engineering, 
2008, 39(4): 53-55. 
5 Guo Shu-xiang. Adaptive finite element method and its applications in engineering. Advances 
in mechanics, 1997, 27(4): 480-488. 
6 Chen Sheng-hong. Phân tích ổn định nền đá phức tạp đập cao với mái dốc đá cao. Nhà xuất 
bản thủy lợi thủy điện Trung quốc, 2001. 
7 Yu Xue-gin, D.Redekop. Applications and Development of Adaptive Finite Element. Journal 
of Nanchang Institute of Aeronautical Technology, 1998, No4: 78-83. 
ADAPTIVE FINITE ELEMENT METHOD AND THE APPLICATION IN 
ANALYSIS OF GRAVITY CONCRETE DAM 
MEng. VU HOANG HUNG(1,2), Dr. NGUYEN QUANG HUNG(1) 
（1） Water Resources University 
（2） College of Water Conservancy and Hydropower Engineering, Hohai Univ., Nanjing 210098, China 
Abstract:With the structural mechanic problems, Finite Element Method is one of the best solutions. It 
is, however, an approximate method; thus there exists the undesired error in solving process. Control 
the error to get the best results and to assess the convergence in analysis the stress-strain of structure is 
posed. This is a problem get the attention of many scientists. The paper focuses on the non-linear 
geometry approach in using Finite element Method to analysis the stress of Concrete gravity Dam on 
rock foundation; through by, the computational technique of gravity dam is amended. The study can 
help the engineers solving the local stress problem in the structure, a disadvantage in design concrete 
gravity dam in Vietnam. 
Keyword: FEM, stress, displacement, adaptive, norm error