Cải tiến phần tử vỏ khối tứ giác 8 nút bằng kỹ thuật khử khóa màng MITC4+ dùng phân tích tĩnh kết cấu tấm/vỏ

Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 45(01/2018) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
43 
CẢI TIẾN PHẦN TỬ VỎ KHỐI TỨ GIÁC 8 NÚT 
BẰNG KỸ THUẬT KHỬ KHÓA MÀNG MITC4+ 
DÙNG PHÂN TÍCH TĨNH KẾT CẤU TẤM/VỎ 
IMPROVEMENT ON 8-NODE QUADRILATERAL SOLID-SHELL 
ELEMENTS BY USING MITC4+ TECHNIQUE TO REMOVE MEMBRANE 
LOCKING FOR STATIC ANALYSES OF PLATE/SHELL STRUCTURES 
Lê Trần Nhật, Châu Đình Thành 
Trường Đại học Sư phạm Kỹ thuật TP.HCM, Việt Nam 
Ngày toà soạn nhận bài 28/6/2017, ngày phản biện đánh giá 2/8/2017, ngày chấp nhận đăng 5/10/2017 
TÓM TẮT 
Phân tích kết cấu tấm/vỏ bằng phần tử vỏ khối tứ giác 8 nút sẽ xảy ra hiện tượng khóa 
cắt, khóa màng và khóa hình thang do sử dụng hàm xấp xỉ chuyển vị dạng C0. Để khắc phục 
các hiện tượng khóa này, các trường biến dạng uốn, biến dạng màng và biến dạng thẳng theo 
phương chiều dày được xấp xỉ lại thông qua giá trị biến dạng tại các điểm buộc được xác 
định trước. Trong bài báo này, ngoài việc khử khóa cắt và khóa hình thang, phần tử vỏ khối tứ 
giác 8 nút, gọi là S8+, còn được khử khóa màng dựa trên kỹ thuật MITC4+ được phát triển 
cho phần tử vỏ suy biến tứ giác 4 nút[1]. Phân tích tĩnh của một số bài toán tấm/vỏ điển hình 
được trình bày. Kết quả số cho thấy, khi sử dụng cùng số phần tử, phần tử S8+ cho kết quả 
chuyển vị tốt hơn so với kết quả cho bởi các loại phần tử vỏ tứ giác và tam giác khác. 
Từ khóa: tấm/vỏ; phần tử vỏ khối tứ giác 8 nút; kỹ thuật MITC4+; khóa màng; phân tích tĩnh. 
ABSTRACT 
Analyses of plate/shell structures by using 8-node quadrilateral solid-shell elements often 
lead to the phenomena of shear, membrane and trapezoidal lockings due to the C
0
-type 
displacement approximation. To overcome these locking phenomena, the bending strains, 
membrane strains and normal strains in the thickness direction are separately interpolated 
from values of these strains at pre-defined typing points. In this paper, beside removing the 
bending and trapezoidal lockings, the present 8-node quadrilateral solid-shell elements, 
namely S8+, are also eliminated the membrane locking based on the MITC4+ technique 
developed for the 4-node quadrilateral degenerated shell elements[1]. Static analyses of some 
benchmark plate/shell structures are presented. Numerical results show that when using the 
same number of elements, the S8+ elements can give better displacements than those provided 
by other quadrilateral and triangular shell elements. 
Keywords: plate/shell; 8-node quadrilateral solid-shell elements; MITC4+ technique; 
membrane locking, static analyses.
1. GIỚI THIỆU 
Tấm/vỏ là kết cấu được sử dụng rất phổ 
biến trong các công trình xây dựng dân dụng, 
xây dựng công nghiệp, máy bay, tàu biển, ô 
tô, do có nhiều ưu điểm vượt trội về khả năng 
chịu lực và thẩm mỹ. Do đó, việc tính toán 
chính xác ứng xử của các kết cấu tấm/vỏ có 
hình dạng, tải trọng và điều kiện biên bất kỳ 
là rất cần thiết và được giải quyết bằng các 
phương pháp số. 
Hiện nay có rất nhiều phương pháp số 
khác nhau được phát triển để giải quyết bài 
toán tấm/vỏ nhưng phương pháp phần tử hữu 
hạn (PPPTHH) vẫn là phương pháp phổ biến 
và hiệu quả nhất. Trong đó, các công thức 
PTHH xây dựng cho lý thuyết tấm/vỏ biến 
44 
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 45(01/2018) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
dạng cắt bậc nhất rồi dùng các kỹ thuật loại 
bỏ các hiện tượng biến dạng vượt trội thường 
được sử dụng do chỉ cần hàm xấp xỉ chuyển 
vị dạng C0 và có thể dùng để phân tích kết 
cấu tấm/vỏ dày hoặc mỏng. Theo hướng phát 
triển này, các công thức PTHH tấm/vỏ có thể 
được chia thành 3 loại: phần tử vỏ phẳng (flat 
shell), phần tử vỏ suy biến (degenerated 
continuum mechanics shell) và phần tử vỏ 
khối (solid shell)[2]. Phần tử vỏ phẳng và 
phần tử vỏ suy biến tính toán các ứng xử trên 
mặt trung bình của phần tử nên ngoài chuyển 
vị thẳng cần có góc xoay. Ứng xử trong phần 
tử vỏ phẳng là sự kết hợp thuần túy giữa ứng 
xử của phần màng và phần uốn nên không 
thể hiện được sự tương tác giữa màng và uốn. 
Để khắc phục vấn đề này, ứng xử trong phần 
tử vỏ suy biến được xây dựng từ lý thuyết 
biến dạng 3 chiều được hiệu chỉnh sao cho 
biến dạng thẳng theo phương vuông góc với 
mặt trung bình của vỏ bằng 0. Trong phần tử 
vỏ khối, ứng xử được xây dựng hoàn toàn từ 
lý thuyết biến dạng 3 chiều nên bậc tự do của 
phần tử vỏ khối chỉ gồm các chuyển vị thẳng, 
không có góc xoay. Vì vậy, phần tử vỏ khối 
rất dễ kết nối với các loại phần tử khác và có 
thể xét đến sự thay đổi của biến dạng thẳng 
theo phương vuông góc với mặt trung bình 
của vỏ. 
Phần tử vỏ khối tứ giác đơn giản nhất là 
phần tử 8 nút vì xấp xỉ trường chuyển vị 
được xây dựng từ các hàm dạng C0. Tuy 
nhiên, biến dạng cắt thuần túy tính từ xấp xỉ 
chuyển vị dạng C0 sẽ không thể bằng 0 và 
dẫn đến hiện tượng khóa cắt khi phân tích 
các kết cấu tấm/vỏ mỏng. Để loại bỏ hiện 
tượng khóa cắt, các phương pháp xấp xỉ lại 
trường biến dạng cắt đã được phát triển thành 
công bởi nhiều tác giả như xấp xỉ biến dạng 
tự nhiên (ANS)[3], xấp xỉ biến dạng nâng 
cao (EAS) [4]hoặc xấp xỉ phối hợp các thành 
phần ten-xơ ứng suất (MITC)[5]. Thêm nữa, 
khi rời rạc các kết cấu vỏ có độ cong lớn thì 
các phần tử vỏ khối có dạng hình thang sẽ 
sinh thêm ứng suất pháp theo phương chiều 
dày và dẫn đến hiện tượng gọi là khóa hình 
thang. Tương tự kỹ thuật khử khóa cắt, 
trường biến dạng thẳng theo phương chiều 
dày của vỏ cũng được xấp xỉ lại để loại bỏ 
hiện tượng khóa hình thang[6]. Phần tử vỏ 
khối tứ giác 8 nút khử khóa cắt và khóa màng 
đã được áp dụng thành công trong việc phân 
tích các kết cấu tấm/vỏ đồng nhất, composite 
đàn hồi tuyến tính, phi tuyến. 
Đối với phần tử vỏ tứ giác, hiện tượng 
khóa màng xảy ra khi các nút của phần tử 
không đồng phẳng. Đặc biệt, khi kết cấu vỏ 
có hình dáng phức tạp chia lưới phần tử 
không đều nhau thì hiện tượng khóa màng 
càng ảnh hưởng nhiều đến độ chính xác của 
phần tử. Nhằm loại bỏ hiện tượng khóa màng, 
một số tác giả xây dựng công thức phần tử vỏ 
dựa trên mô hình toán học biểu diễn chính 
xác mặt trung bình vỏ[7]. Tuy nhiên, cách 
tiếp cận này không tổng quát cho tất cả các 
trường hợp vỏ có hình dáng và điều kiện biên 
bất kỳ. Gần đây, Ko và cộng sự[1] đã phát 
triển công thức xấp xỉ trường biến dạng 
màng có khả năng khử khóa màng dựa trên 
cách tiếp cận MITC cho phần tử vỏ suy biến 
tứ giác 4 nút. Kết quả, phần tử MITC4+ có 
khả năng vượt qua các điều kiện kiểm tra cơ 
bản của phần tử vỏ (patch tests) và có độ hội 
tụ tốt khi phân tích các bài toán vỏ phức tạp, 
thậm chí với lưới không đều. 
Trong nghiên cứu này, ngoài việc khử 
khóa cắt và khóa hình thang, công thức 
PTHH vỏ khối tứ giác 8 nút còn được bổ 
sung thêm kỹ thuật khử khóa màng. Muốn 
vậy, biến dạng màng trong phần tử vỏ khối sẽ 
được xấp xỉ lại theo công thức xấp xỉ biến 
dạng màng của phần tử vỏ suy biến tứ giác 4 
nút MITC4+. Công thức PTHH vỏ khối tứ 
giác 8 nút S8+ này cho kết quả hội tụ tốt hơn 
một số phần tử vỏ khác khi phân tích tĩnh 
một số bài toán tấm/vỏ. 
Trong phần tiếp theo của bài báo, công 
thức PTHH vỏ khối tứ giác 8 nút được trình 
bày. Tính chính xác và hiệu quả của phần tử 
S8+ được kiểm tra và đánh giá ở phần ví dụ 
số. Cuối cùng, một vài kết luận được rút ra. 
2. CÔNG THỨC PHẦN TỬ HỮU HẠN 
VỎ KHỐI TỨ GIÁC 8 NÚT S8+ 
2.1 Xấp xỉ hình học và chuyển vị[8] 
Xét phần tử vỏ khối tứ giác 8 nút như 
Hình 1. Trong hệ trục tọa độ toàn cục OXYZ, 
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 45(01/2018) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
45 
tọa độ của 1 điểm nằm trong phần tử có thể 
được biểu diễn bởi véc-tơ X được xác định 
như sau: 
 X X P (1) 
Trong đó, X là véc-tơ vị trí của mặt trung 
bình phần tử, P là véc-tơ có giá trị bằng một 
nửa véc-tơ nối từ tọa độ nút mặt dưới với tọa 
độ nút mặt trên của phần tử và -1  1 là tọa 
độ thuộc hệ tọa độ tự nhiên . 
Hình 1. Hình dạng hình học phần tử vỏ khối 
Véc-tơ X và P được xấp xỉ thông qua 
các giá trị tại nút phần tử bởi công thức: 
4
1
i i
i
N
 X X (2) 
4
1
i i
i
N
 P P (3) 
Với Ni(,) là các hàm dạng C
0
 được định 
nghĩa trong hệ tọa độ tự nhiên như sau: 
 1 20,25 1 1 ; 0,25 1 1N N    
 3 40,25 1 1 ; 0,25 1 1N N    
Và 
1
2
T
i i
B
i X X X (4) 
1
2
T B
i i i P X X (5) 
Trong đó, Xi
T
 và Xi
B
 (i = 1, 2, 3, 4) lần lượt là 
tọa độ nút ở mặt trên và mặt dưới của phần tử 
(xem Hình 1).
Thế công thức(2), (3) vào (1), véc-tơ X 
có thể được viết lại như sau: 
4
1
i i
i
N 
 X X Pi (6) 
Để thuận lợi cho việc thiết lập các quan 
hệ biến dạng – chuyển vị và ứng suất – biến 
dạng, hệ trục tọa độ cục bộ oxyz trong từng 
phần tử được xây dựng với các véc-tơ cơ sở 
theo phương x, y, và z lần lượt là Rx, Ry vàRz 
được định nghĩa như sau: 
, , z x y xz
  
  
R R R
R R R R
R R
R
R
 (7) 
Với 
4
10
i
i
i
N


  

 

X
R X (8a) 
10
4
i
i
i
N


  

 

X
R X (8b) 
Trường chuyển vị trong hệ trục tọa độ 
cục bộ oxyz được xấp xỉ như sau: 
4
1
i i i
i
N 
 u u uρ ρ (9) 
Trong đó, u = [ux, uy, uz]
T
 là véc-tơ chứa các 
thành phần chuyển vị thẳng theo phương x, y, 
z, 
T
x y zu u u u là véc-tơ chuyển vị trên 
mặt trung bình, = [ x y z]
T 
là một nửa véc-tơ 
nối chuyển vị nút dưới và nút trên của phần tử 
và: 
1
2
T
xi yi
T B
i i iziu u u u u u (10) 
1
2
T
xi yi
T B
i i izi u uρ (11) 
Với T
iu và
 B
iu lần lượt là véc-tơ chuyển vị 
trong hệ trục tọa độ cục bộ của các nút ở mặt 
trên và mặt dưới của phần tử. 
Quan hệ giữa tọa độ nút X và chuyển vị 
nút U trong hệ tọa độ toàn cục OXYZ với tọa 
độ nút x và chuyển vị nút u trong hệ tọa độ cục 
bộ oxyz được cho bởi công thức sau: 
T x H X hoặc X Hx (12a) 
T u H U hoặc U Hu (12b) 
Với ma trận chuyển trục H = [RxRyRz] (13) 
1B 2B
4B4
B
1T 2T
4T4
T
X1 X2
X3
X4
2P4 2P3
2P22P1
Z
Y
X
O

x
y
z
46 
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 45(01/2018) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
2.2 Quan hệ giữa biến dạng-chuyển vị[8] 
Từ xấp xỉ chuyển vị cho bởi (9), các 
thành biến dạng trong hệ tọa độ cục bộ được 
xác định như sau: 
x x x
xx
u u
x x x
 
  
  
 (14a) 
y y y
yy
u u
y y y
 
  
  
 (14b) 
z z z
zz z
u u
z z z z
 
  
   
   
 (14c) 
y y yx x x
xy
u uu u
y x y x y x
 
     
      
 (14d) 
x x xz z z
xz x
u uu u
z x z x z x z
 
  
      
       
 (14e) 
y y yz z z
yz y
u uu u
z y z y z y z
 
  
      
       
 (14f) 
Trường biến dạng cho bởi (14) có thể viết lại:
m b
T
xx yy xy z xz yz zz
q

      
  
  
 
ε ε
ε
ε
 (15) 
Trong đó, biến dạng màng m, biến dạng uốn 
b và biến dạng cắt ngoài mặt phẳng q được 
định nghĩa bởi: 
T
y yx x
m
u uu u
x y y x
   
    
ε (16) 
T
y yx x
b
x y y x
    
    
ε (17) 
x x z z
q
y z
y
z
x
y
u u
z z x z x
u u
z z y z y
 
 
 
 
     
     
    
     
ε (18) 
Thế xấp xỉ chuyển vị cho bởi (9), (10), 
(11) vào các công thức biến dạng, ta có thể 
xác định được quan hệ giữa biến dạng và các 
chuyển vị nút qe của phần tử trong hệ tọa độ 
cục bộ như sau: 
; ; ; m m e b b e q q e zz z e ε B q ε B q ε B q B q (19) 
Trong đó, qe = [qe1
T
qe2
T
qe3
T
qe4
T
]
T
với qei = [uxi
B
uyi
B
uzi
B
uxi
T
uyi
T
uzi
T
]
T
 và Bm, Bb, Bq, 
Bz là các ma trận quan hệ giữa biến dạng và 
chuyển vị. 
2.3 Kỹ thuật khử khóa màng, khóa cắt và 
khóa hình thang 
Nếu sử dụng các trường biến dạng cho 
bởi công thức (19) được xác định thuần túy 
từ xấp xỉ chuyển vị dạng C0 thì công thức 
PTHH của phần tử vỏ khối sẽ xuất hiện các 
ứng xuất tăng thêm dẫn đến sai lệch kết quả 
tính toán, gọi là các hiện tượng khóa. Để 
khắc phục các hiện tượng khóa màng, khóa 
cắt và khóa hình thang, các trường biến màng, 
biến dạng cắt và biến dạng thẳng theo chiều 
dày trong hệ tọa độ tự nhiên sẽ được xấp xỉ 
thông qua các giá trị biến dạng tương ứng 
tính tại các điểm buộc được thiết kế trước. 
Để có thể tính lại được ma trận quan hệ 
giữa biến dạng và chuyển vị trong hệ tọa độ 
cục bộ từ xấp xỉ biến dạng trong hệ trục tọa 
độ tự nhiên, ta cần biểu diễn quan hệ giữa 
biến dạng trong hệ tọa độ cục bộ và hệ tọa độ 
tự nhiên như sau[8]: 
1
T
xx yy zz xy xz yz
T
          
      
 Q
 (20) 
Trong đó, 
2 2 2
11 21 31 11 21 11 31 21 31
2 2 2
12 22 32 11 22 12 32 22 32
2 2 2
13 23 33 13 23 13 33 23 33
11 12 21 22 31 32 12 21 11 22 12 31 11 32 22 31 21 32
11 13 21 23 31 33 23 11 21 13 13 13 11 33
2 2 2
2 2 2
I I I I I I I I I
I I I I I I I I I
I I I I I I I I I
I I I I I I I I I I I I I I I I I I
I I I I I I I I I I I I I I I
Q =
23 31 21 33
11 13 11 13 11 13 23 12 22 13 13 32 12 33 23 32 22 332 2 2
I I I
I I I I I I I I I I I I I I I I I I
 (21) 
Với Iij là số hạng ở hàng i, cột j của ma trận 
nghịch đảo của ma trận Jacobi J được xác 
định bởi: 
T
  
   
    
x x x
J (22) 
2.3.1 Kỹ thuật khử khóa màng 
Hiện tượng khóa màng xuất hiện khi 
phần tử vỏ khối tứ giác 8 nút có các điểm 
trên mặt trung bình không đồng phẳng. Để 
khắc phục hiện tượng này, kỹ thuật khử khóa 
màng MITC4+ được đề xuất bởi Ko và cộng 
sự [1] cho phần tử vỏ suy biến tứ giác 4 nút 
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 45(01/2018) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
47 
được áp dụng. Theo cách tiếp cận này, các 
thành phần biến dạng màng trong hệ tọa độ 
tự nhiên được xấp xỉ lại dựa vào giá trị biến 
dạng màng tính tại các điểm buộc 
A1(0,-1,0),A2(1,0,0), A3(0,1,0), A4(-1,0,0) 
(xem Hình 2) như sau: 
31 2 4
34 2 1
0,25
 0,5 0,5
AA A A
ij ij ij ij ij
AA A A
ij ij ij ij
    
     
 (23) 
Trong đó, k
A
ij với i = , , j = ,  và k = 1..4 
thể hiện các biến dạng màng trong hệ tọa độ 
tự nhiên tính từ xấp xỉ chuyển vị (9) tại các 
điểm buộc Ak. 
Từ quan hệ (14) giữa chuyển vị và biến 
dạng trong hệ tọa độ cục bộ và công thức 
(20) chuyển biến dạng từ hệ tọa độ cục bộ 
sang hệ tọa độ tự nhiên, xấp xỉ biến dạng 
màng (23) theo kỹ thuật MITC4+ có thể biểu 
diễn bằng quan hệ giữa biến dạng màng 
trong hệ tọa độ cục bộ và chuyển vị nút phần 
tử như sau: 
m m e ε B q (24) 
Hình 2. Các điểm buộc khử khóa màng, cắt 
2.3.2 Kỹ thuật khử khóa cắt 
Theo phương pháp khử khóa cắt 
MITC4[5], biến dạng cắt ngoài mặt phẳng 
trong hệ tọa độ tự nhiên được xấp xỉ lại như 
sau: 
31
2 4
0,5 1 0,5 1
0,5 1 0,5 1
AA
A A
  
  
    
    
 (25) 
Với 31 2 4, , , 
AA A A
       lần lượt là giá trị của 
biến dạng cắt ngoài mặt phẳng tính từ xấp xỉ 
chuyển vị (9) theo các công thức (14) và (20) 
tại các điểm buộc có tọa độ trong hệ tự nhiên 
A1(0,-1,0), A3(0,1,0), A2(1,0,0), A4(-1,0,0) 
như Hình 2. 
Thế các biến dạng cắt ngoài mặt phẳng 
được tính theo chuyển vị nút tại các điểm 
buộc vào xấp xỉ (25),quan hệ giữa biến dạng 
cắt ngoài mặt phẳng trong hệ tọa độ cục bộ 
và chuyển vị nút phần tử được viết lại: 
q q e ε B q (26) 
2.3.3 Kỹ thuật khử khóa hình thang 
Khi vỏ cong hoặc/và chia lưới không 
đều sẽ xuất hiện thêm ứng suất pháp theo 
phương chiều dày vỏ do phần tử có dạng 
hình thang theo phương chiều dày, gọi là 
hiện tượng khóa hình thang. Để khắc phục 
hiện tượng này, biến dạng thẳng  được xấp 
xỉ lại thông qua các giá trị 31 2 4, , , 
BB B B
      
cho bởi xấp xỉ chuyển vị tính tại các điểm 
buộc B1(-1,-1,0), B2(1,-1,0), B3(1,1,0), 
B4(-1,1,0) (xem Hình 3) như sau[6]: 
 1 20,25 1 1,25 11 0B B         
 3 41 10,25 1 0,25 1B B       (27) 
Hình 3. Các điểm buộc khử khóa hình thang 
Thế quan hệ giữa biến dạng  tại các 
điểm buộc với chuyển vị nút phần tử vào xấp 
xỉ (27), ta được: 
zz z e B q (28) 
2.4 Ma trận độ cứng phần tử S8+ 
Từ các công thức (24), (26), (28), sau khi 
sử dụng các kỹ thuật khử khóa cắt, khóa màng 
và khóa hình thang, quan hệ giữa biến dạng và 
chuyển vị nút phần tử được viết lại: 
1 2
34
5 6
78
A1
A2
A3
A4

1 2
34
5 6
78
B1 B2
B3B4

48 
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 45(01/2018) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
m b
zz e
q


  
 
 
ε ε
ε Bq
ε
 (29) 
Đối với vật liệu đồng nhất, đẳng hướng, 
ứng suất được tính từ quan hệ: 
T
xx yy xy zz xz yz      σ Cε (30) 
Trong đó, 
2 0 0 0 0
2 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
E
  
  



C
Với
 2
, 
1 2 1
Ev E
v v
  
và E,  lần lượt là 
mô-đun đàn hồi và hệ số Poisson của vật liệu. 
Thế quan hệ giữa ứng suất – biến dạng và 
biến dạng – chuyển vị nút cho bởi (29) và 
(30)vào nguyên lý công ảo, theo phương pháp 
PTHH, ma trận độ cứng phần tử ke được xác 
định như sau[9]: 
1 1 1
1 1 1
d d dTe   
 k B CB J (31) 
Để tính ke, tích phân số với 2x2x2 điểm 
Gauss được sử dụng, tức là 2 điểm Gauss 
theo mỗi phương , , . 
3. VÍ DỤ SỐ 
Trong phần này, kết quả phân tích một 
số bài toán tấm vỏ điển hình bằng phần tử 
S8+ được trình bày. Kết quả phân tích được 
so sánh với kết quả cho bởi các loại phần tử 
vỏ khối tứ giác 8 nút Xsolid84 dựa trên dạng 
chuyển vị không tương thích[8],Xsolid85 sử 
dụng kỹ thuật ANS khử khóa cắt và khóa 
hình thang[8], RH8s-2 kết hợp kỹ thuật làm 
trơn vào phần tử Xsolid85 [10], và các loại 
phần tử vỏ phẳng tứ giác MITC4, tam giác 
MITC3+[11]. 
Để khảo sát và so sánh độ hội tụ của các 
phần tử, kết cấu ở các ví dụ sau được chia 
lưới đều với nX x nY phần tử trong mặt trung 
bình và 1 phần phần tử theo phương chiều 
dày đối với phần tử vỏ khối. Trong đó, nX, nY 
lần lượt là số phần tử theo phương X và Y. 
3.1 Tấm xiên liên kết tựa đơn bốn cạnh 
chịu tải phân bố đều 
Tấm xiên góc 300có kích thước mỗi cạnh 
L=100 và dày h=1 như Hình 4. Tấm tựa đơn 
4 cạnh và chịu tải trọng phân bố đều q=1. Vật 
liệu làm tấm có mô-đun đàn hồi E = 10x106, 
và hệ số Poisson  = 0.3. 
Kết quả tính chuyển vị tại tâm tấm chia 
cho giá trị chuyển vị tham khảo 0.04455 
[8]ứng các loại lưới và phần tử khác nhau 
được thể hiện trong Bảng 1và Hình 5. Kết 
quả cho thấy, với cùng số phần tử, phần tử 
S8+ cho kết quả gần với lời giải tham khảo 
hơn các loại phần tử khác. Với lưới chia 
20x20x1 phần tử, sai số tương đối của kết 
quả cho bởi phần tử S8+, Xsolid84, MITC3+, 
Xsolid85, MITC4 với kết quả tham khảo lần 
lượt là 0.017%, 0.070%, 0.173%, 0.052%, 
0.130%. 
Hình 4.Tấm xiên liên kết tựa đơn 4 cạnh 
chịu tải phân bố đều 
Bảng 1. Độ võng tương đối tại tâm tấm xiên 
Phần tử 
Chuyển vị theo số lưới chia 
4x4x1 8x8x1 20x20x1 
Xsolid84 0.320 0.676 0.930 
MITC3+ 0.523 0.653 0.827 
Xsolid85 0.879 0.870 0.948 
MITC4 0.805 0.802 0.870 
S8+ 0.884 0.900 0.983 
X
Y
3
0
°
C
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 45(01/2018) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
49 
Hình 5. Chuyển vị lớn nhất của tấm xiên 
theo lưới chia 
3.2 Mái vòm ngàm phẳng chịu tải trọng 
lượng bản thân 
Mái vòm có chiều dài L=50, bán kính 
R=25 và dày h=0.25. Mái vòm chịu trọng 
lượng bản thân p =90, làm bằng vật liệu có E 
= 4.32x10
8
 và = 0. 
Do tính chất đối xứng nên 1/4 mái được 
mô phỏng với điều kiện biên được thể hiện 
trên Hình 6. 
Tỉ số giữa chuyển vị thẳng tại điểm A 
(điểm giữa biên cạnh dài của mái) và kết quả 
tham khảo 0.3024 [8] cho bởi phần tử S8+ và 
các loại phần tử khác nhau được trình bày 
trong Bảng 2. Trong ví dụ này, phần tử S8+ 
cũng cho giá trị chuyển vị gần với lời giải 
tham khảo hơn các loại phần tử khác. 
Hình 6. Mái vòm ngàm phẳng chịu tải 
trọng lượng bản thân 
Bảng 2. Chuyển vị tương đối mái vòm 
tại điểm A 
Phần tử 
Chuyển vị theo số lưới chia 
4x4x1 8x8x1 16x16x1 
Xsolid84 0.549 0.959 1.003 
MITC3+ 0.762 0.916 1.014 
RH8s-2 1.179 1.680 1.028 
MITC4 0.943 0.978 1.008 
S8+ 0.946 0.980 0.995 
3.3 Vỏ bán cầu có lỗ mở 
Vỏ bán cầu có lỗ mở 180có bán kính 
R=10 và dày h = 0.04. Vỏ chịu lực tập trung 
F = 2 như Hình 7 và có E=6.825x107, =0.3. 
Chuyển vị tại vị trí lực tập trung 0.094 [8]. 
Do tính chất đối xứng nên chỉ cần mô 
phỏng 1/4 vỏ bán cầu. Điều kiện biên đối 
xứng xác định như Hình 7, biên trên và biên 
dưới của vỏ tự do. 
Hình 7.Vỏ bán cầu với lỗ mở 180 
Bảng 3. Chuyển vị tương đối vỏ bán cầu lỗ mở 
18
0
 tại điểm B 
Phần tử 
Chuyển vị theo số lưới chia 
8x8x1 16x16x1 20x20x1 
Xsolid85 1.005 - - 
MITC3+ 1.032 1.012 1.012 
RH8s-2 1.003 0.996 0.996 
MITC4 0.971 0.988 - 
S8+ 0.985 0.992 1.000 
0
0.2
0.4
0.6
0.8
1
1.2
4 X 4 X 1 8 X 8 X 1 2 0 X 2 0 X 1 C
h
u
y
ể
n
 v
ị 
tư
ơ
n
g
 đ
ố
i 
tạ
i 
C
Lưới chia 
Xsolid84 MITC3+ Xsolid85
S8+ MITC4 Tham khảo 
40°
u
y = 0
u x
 =
 0
A
ux,uz = 0
Z
X
Y
O
L
18°
R
X YF F
F
F
Z
u y
 =
 0
u
x =
 0
B
50 
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 45(01/2018) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
Kết quả Bảng 3 cho thấy phần tử S8+ 
cho kết quả hội tụ tương đương các phần tử 
khác. Với lưới chia 20x20x1, phần tử S8+ 
cho kết quả tốt hơn các phần tử khác (độ lệch 
so với kết quả tham khảo là 0%). 
4. KẾT LUẬN 
Phần tử vỏ khối tứ giác 8 nút S8+ bổ 
sung khử khóa màng theo kỹ thuật MITC4+ 
cho hiệu quả tính toán tốt hơn các phần tử vỏ 
suy biến tính trên mặt trung bìnhMITC3+, 
MITC4 và các phần tử vỏ khối tứ giác 8 nút 
khác chỉ thực hiện khử khóa cắt và khóa hình 
thang Xsolid84, Xsolid85, thậm chí tốt hơn 
phần tử vỏ khối kết hợp kỹ thuật làm trơn 
RH8s-2. Do mô phỏng và tính toán được ứng 
xử theo phương chiều dày vỏ nên phần tử vỏ 
khối S8+ có thể dùng để phân tích các kết 
cấu tấm/vỏ composite có xét ảnh hưởng của 
ứng suất theo chiều dày, mô phỏng tự tách 
lớp hay sự hình thành và phát triển vết nứt 
theo chiều dày tấm/vỏ.
TÀI LIỆU THAM KHẢO 
[1] Y. Ko, P.-S. Lee, and K.-J. Bathe, “The MITC4+ shell element and its performance,” 
Computers & Structures169, 57–68 (2016). 
[2] H. T. Y. Yang, S. Saigal, A. Masud, and R. K. Kapania, “A survey of recent shell finite 
elements,” Int. J. Numer. Meth. Engng.47, 101–127 (2000). 
[3] K.-J. Bathe and E. N. Dvorkin, “A four-node plate bending element based on 
Mindlin/Reissner plate theory and a mixed interpolation,” Int. J. Numer. Meth. Engng.21, 
367–383 (1985). 
[4] U. Andelfinger and E. Ramm, “EAS-elements for two-dimensional, three-dimensional, 
plate and shell structures and their equivalence to HR-elements,” Int. J. Numer. Meth. 
Engng.36, 1311–1337 (1993). 
[5] K.-J. Bathe and E. N. Dvorkin, “A formulation of general shell elements—the use of 
mixed interpolation of tensorial components,” Int. J. Numer. Meth. Engng.22, 697–722 
(1986). 
[6] K. Y. Sze and L. Q. Yao, “A hybrid stress ANS solid-shell element and its generalization 
for smart structure modelling. Part I—solid-shell element formulation,” Int. J. Numer. 
Meth. Engng.48, 545–564 (2000). 
[7] G. M. Kulikov and S. V. Plotnikova, “A family of ANS four-node exact geometry shell 
elements in general convected curvilinear coordinates,” Int. J. Numer. Meth. Engng.83, 
1376–1406 (2010). 
[8] K. D. Kim, G. Z. Liu, and S. C. Han, “A resultant 8-node solid-shell element for 
geometrically nonlinear analysis,” Comput Mech35, 315–331 (2004). 
[9] K.-J. Bathe, Finite Element Procedures (Prentice Hall International, Inc., 1996). 
[10] X. J.-G. Élie-Dit-Cosaque, A. Gakwaya, and H. Naceur, “Smoothed finite element 
method implemented in a resultant eight-node solid-shell element for geometrical linear 
analysis,” Comput Mech55, 105–126 (2015). 
[11] P. Nguyễn Hoàng, “Phân tích kết cấu vỏ bằng phần tử MITC3+ được làm trơn trên phần 
tử với hàm Bubble (bCS-MITC3+)” (Trường ĐH Sư phạm Kỹ thuật Tp.HCM, 2016). 
Tác giả chịu trách nhiệm bài viết: 
Lê Trần Nhật 
Trường Đại học Sư phạm Kỹ thuật TP.HCM 
Email:trannhat93@gmail.com