Phân tích tấm reissner-mindlin bằng phần tử es+ns-mitc3

PHÂN TÍCH TẤM REISSNER-MINDLIN BẰNG PHẦN TỬ  
ES+NS-MITC3  
Châu Đình Thành a,*, Hồ Thị Còn b, và Lê Phương Bình a  
a Khoa Xây Dựng, Trường ĐH Sư Phạm Kỹ Thuật TP.HCM,  
Quận Thủ Đức, TP. Hồ Chí Minh  
b Trung Tâm Giáo Dục Nghề Nghiệp – Giáo Dục Thường Xuyên,  
Huyện An Phú, Tỉnh An Giang  
Tóm tắt  
Trong nghiên cứu này, phương pháp phần tử hữu hạn biến dạng trơn (S-FEM) dựa trên cạnh  
và trên nút phần tử (ES-FEM và NS-FEM) được kết hợp để phát triển cho phần tử tấm tam  
giác 3 nút sử dụng kỹ thuật MITC3 để khử hiện tượng khóa cắt. Cách tiếp cận này dựa trên  
trên phương pháp -FEM với tham số là hệ số hiệu chỉnh tỉ lệ miền được làm trơn trên  
cạnh và trên nút phần tử. Phần tử đề xuất ES+NS-MITC3 có các trường biến dạng uốn và cắt  
được làm trơn một phần trên cạnh (ES) và phần còn lại trên nút (NS) phần tử. Phần tử  
ES+NS-MITC3 được áp dụng để phân tích một số loại tấm Reissner-Mindlin. Kết quả số cho  
thấy, phần tử ES+NS-MITC3 cho kết quả tốt hơn các loại phần tử chỉ sử dụng phương pháp  
ES-FEM hoặc NS-FEM, đặc biệt là trong các trường hợp tấm mỏng.  
Từ khóa: Tấm Reissner-Mindlin, MITC3, ES-FEM, NS-FEM.  
Abstract  
In this thesis, the edge and node-based smoothed strain finite element method (ES-FEM and  
NS-FEM) were combined to develop for 3-node using the MITC3 technique to attenuate the  
shear-locking. This approach is based on the -FEM in which the parameter is used to tune  
the contribution ratio of the edge-based and node-based. The ES+NS-MITC3 element in  
which the bending and shear strain fields are smoothed on a part of the edge-based domains  
and the other on the node-based domains. The ES+NS-MITC3 element phenomenon is  
employed some of the similar Reissner-Mindlin plates. Numerical results show that, the  
ES+NS-MITC3 element can give results better than those of the similar elements using the  
ES-FEM or NS-FEM only, particular in the thin plates.  
Keywords: Plates Reissner-Mindlin, MITC3, ES-FEM, NS-FEM.  
1
1. Giới thiệu  
cạnh và trên nút ES-FEM và NS-FEM đã  
được tác giả Chen và LI [18] đề xuất cho  
các vấn đề về tấm Reissner – Mindlin, là  
phương pháp làm trơn dựa trên phiên bản  
Kết cấu tấm được sử dụng phổ biến trong  
xây dựng công trình hiện nay. Từ những  
yêu cầu giải các bài toán phức tạp về lý  
thuyết đàn hồi, phân tích kết cấu trong xây  
dựng. Vì thế mà lý thuyết tấm Kirchoff [1]  
và lý thuyết tấm của Mindlin [2] được sử  
dụng phổ biến. Với ý tưởng thiết lập công  
thức PTHH [3], cho loại phần tử đơn giản  
có thể áp dụng tính toán cho cả tấm dày,  
tấm mỏng với độ chính xác cao, việc xây  
dựng hàm xấp xỉ chuyển vị C0 đơn giản  
hơn so với xây dựng hàm xấp xỉ C1 nhưng  
khi xấp xỉ chuyển vị C0 sẽ dẫn đến hiện  
tượng khóa cắt “shear locking ”. Đến năm  
1985, Bathe và Dvorkin [4 -5] đã tìm ra  
phương pháp nội suy hỗn hợp các thành  
phần ten-xơ (MITC - Mixed Interpolation  
Tensorial Components) và phương pháp  
MITC nhanh chóng giải quyết được vấn đề  
“khoá cắt”. Sau đó, phương pháp MITC đã  
được phát triển rất thành công, cho các loại  
phần tử tứ giác 4 nút như MITC4, 8 nút  
MITC8 [6], 9 nút MITC9 [7], 16 nút  
MITC16 [8],… và MITC4+ [9]. Nối tiếp có  
rất nhiều các nghiên cứu, phương pháp  
PTHH trơn đã được phát triển cho các phần  
tử tấm tam giác 3 nút sử dụng kỹ thuật khử  
khóa cắt DGS3 [10], MIN3 [11], MITC3  
[12] và MITC3+ [13] như các phần tử, NS-  
DSG3 [14], ES-DSG3 [15], ES-MITC3  
[16], NS-MITC3 [17],…  
(
beta PTHH FEM).  
Kỹ thuật làm trơn trên cạnh (ES) và trên  
nút (NS) phần tử của phương pháp PTHH  
làm trơn đã được kết hợp, để có thể tận  
dụng lợi thế về tính chất độ cứng cho bởi cả  
hai kỹ thuật làm trơn này. Cách tiếp cận  
này còn gọi là beta PTHH (-FEM).  
Gần đây, tác giả Li and He [19] phát triển  
và tối ưu hóa vị trí của các điểm tích phân  
Gauss dựa trên phần tử tam giác trên cạnh  
được sử dụng trong ES-FEM kết hợp với  
NS-FEM để xây dựng FEM được làm  
trơn tổng quát cho các bài toán tấm.  
Trong nghiên cứu này, sự kết hợp của làm  
trơn trên cạnh (ES) và trên nút (NS) của  
phần tử sẽ được phát triển cho phần tử tấm  
tam giác 3 nút sử dụng kỹ thuật khử khóa  
cắt MITC3, để phân tích kết cấu tấm đồng  
nhất đẳng hướng, qua đó đánh giá kết quả  
với ý tưởng kết hợp giữa 2 phương pháp  
làm trơn khác nhau trên 1 phần tử sẽ cho  
kết quả tính toán chính xác hơn và tốc độ  
hội tụ nhanh hơn so với chỉ sử dụng  
phương pháp làm trơn đơn lẻ trên 1 phần tử  
như ES-MITC3 [ 16], NS-MITC3 [ 17] qua  
so sánh với các phương pháp hiện tại.  
Ngoài ra phương pháp PTHH trơn dựa trên  
1
chiều dương qui ước được định nghĩa như  
2. Lý thuyết tấm biến dạng cắt  
Hình 2.1.  
bậc nhất  
2.2. Trường biến dạng.  
Trong bài báo này, lý thuyết biến dạng cắt  
bậc nhất [2] được áp dụng cho tấm đồng  
nhất đẳng hướng.  
Trường biến dạng được xác định từ trường  
chuyển vị (2.1)  
2.1. Trường chuyển vị.  
u  
x  
v  
y  
x  
x  
x  
z  
z  
Theo giả thuyết tấm Mindlin-Reissner, khi  
tấm phẳng chịu tải trọng phân bố q thì  
đoạn thẳng vuông góc mặt trung bình tấm  
sau khi biến dạng sẽ chuyển vị võng xuống  
1 đoạn thẳng w theo phương z (trong đó w  
là 1 hàm phụ thuộc x, y) và 2 chuyển vị  
y  
y  
y  
0  
z
(2.2)  
y  
u v  
x  
y  
xy  
z  
y x  
x  
u w  
w  
x   
xz  
z x  
x  
w  
y  
v w  
z y  
góc xoay β , β (là 2 góc xoay của đoạn  
x
y
yz  
y   
thẳng vuông góc với mặt trung bình của  
Trong đó: x, y, z là các biến dạng dài  
theo phương x, y, z; xy, xz, yz lần lượt là  
các biến dạng trượt trong mặt phẳng Oxy,  
Oxz, Oyz.  
tấm lần lượt quay quanh trục y, x) như  
Hình 2.1. Cả 3 thành phần này w, β , β  
x
y
được biểu diễn như sau:  
2.3. Trường ứng suất.  
Xét kết cấu tấm được phân tích theo trạng  
thái ứng suất phẳng Hình 2.2:  
Hình 2.1: Trường chuyển vị trong tấm  
u zx, y  
x   
Hình 2.2: Các thành phần ứng suất phân  
bố theo bề dày tấm  
v zx, y  
y   
(2.1)  
Theo định luật Hooke mối quan hệ giữa  
w w x, y  
ứng suất và biến dạng ở Hình 2.2 được  
biểu diễn như sau:  
Trong đó: u, v, w lần lượt là chuyển vị  
thẳng theo phương x, y, z; β , β lần lượt là  
x y  
góc xoay của đoạn thẳng vuông góc với  
mặt trung bình quanh trục y và trục x với  
3
E
  
x
y
x
x
y
2
z
z
1  
x  
E
1
v
0
0
x
y
  
  
  
  
E
y
x
y
2
v
1
1  
2
y  
(2.3)  
0  
1v  
1v  
2   
Z
0
0
xy  
y   
E
x
z
xy  
xy  
 y  
x   
2(1  
)
)
E
  
xz  
yz  
xz  
yz  
x
2(1  
E
2(1  
x  
  
1
v
0
0
Ez  
)
y
v
1
2
y  
1v  
1v  
Trong đó: x, y, z là các ứng suất pháp  
theo phương x, y, z; xy, xz, yz là các ứng  
suất tiếp trên mặt có vectơ pháp tuyến là  
trục x, y, z; E là mô-đun đàn hồi khi kéo  
nén của vật liệu; ν: là hệ số Poisson.  
0
0
  
  
y
2   
x
y  
x  
1 v  
0
Ez  
1v2  
(2.6)  
σ   
v 1  
0
κ
1v  
Theo giả thuyết, bỏ ứng suất pháp theo  
phương z (z = 0). Khi đó trường ứng suất  
(2.3) có thể được viết lại dưới dạng ma  
trận như sau  
0 0  
2   
w  
x   
   
E
x  
w  
y  
xz  
   
yz  
2(1v)  
   
y   
(2.7)  
Ứng suất do biến dạng uốn:  
E
τ   
γ
   
x  
1 v  
v 1  
0
0
   
x   
2(1v)  
E
(2.4)  
(2.5)  
y  
y  
Trong đó: = [x y xy]T; κ là vectơ biến  
1v2  
   
   
xy  
1v xy  
   
0 0  
dạng uốn, được xác định bởi  
2   
x  
x  
y  
xz   
E
xz   
yz  
yz  
2(1v)  
(2.8)  
κ   
Thế quan hệ giữa biến dạng và chuyển vị  
cho bởi (2.2) vào ứng suất (2.4) và  
y  
  
x  
y  
y   
x  
(2.5), ta có:  
Ứng suất do biến dạng cắt  
= [xz yz]T  
và γ là vectơ biến dạng cắt, được xác định  
như sau  
4
T  
w  
x  
w  
y  
M
γ x  
y   
(2.9)  
1 v  
0
x
(2.11)  
MM D v 1 0 κ  
y
2.4. Nội lực trong tấm.  
1v  
M
0 0  
xy  
2   
Tách 1 đơn vị dài trên chiều dày h của  
tấm, khi đó hợp lực các ứng suất phân bố  
trên đơn vị này được gọi là nội lực, hay  
còn gọi là các thành phần ứng lực của tấm.  
Trong đó, M = [Mx My Mxy]T và  
D
là độ  
cứng trụ của tấm, đặc trưng về hình học và  
vật liệu của tấm khi chịu uốn:  
Eh3  
12(1v2 )  
(2.12)  
(2.13)  
D   
Lực cắt trên 1 đơn vị chiều dài:  
h/2  
Q   
  
xzdA   
yzdA   
1dz  
1dz  
xz   
x
A
h/2  
h/2  
Hình 2.3: Các thành phần nội lực của tấm  
Để đưa bài toán ba chiều về bài toán hai  
Q   
  
yz   
h/2  
y
A
chiều, lấy tích phân biểu thức (2.4) theo  
chiều dày tấm, ta nhận được các thành  
phần nội lực như sau:  
Trong đó, Qx, Qy là lực cắt trên 1 đơn vị  
dài lần lượt nằm trong mặt phẳng Oxz và  
Oyz.  
h/2  
z1dz  
M zdA   
x   
Thế ứng suất cắt (2.7) vào lực cắt (2.13), ta  
x
x
A
h/2  
có véc-tơ lực cắt  
h/2  
z1dz  
(2.10)  
M zdA   
y   
y
y
A
h/2  
Q
   
x
1 0  
0 1  
Q   
kGh  
γ
(2.14)  
h/2  
   
Qy  
   
z1dz  
M zdA   
xy   
xy  
xy  
A
h/2  
E
Với  
là mô-đun đàn hồi trượt;  
G   
Trong đó, Mx, My lần lượt là moment uốn  
trên 1 đơn vị dài mặt cắt quay quanh trục  
y, x có chiều dương như Hình 2.3; Mxy là  
moment xoắn trên 1 đơn vị dài mặt cắt  
quay quanh trục trục x.  
2(1v)  
k
= 5/6 là hệ số hiệu chỉnh cắt.  
3. Công thức phần tử ES+NS-  
MITC3  
Thế ma trận ứng suất do biến dạng uốn  
(2.7) vào trường moment (2.10), ta có:  
5
3.1. Công thức PTHH tam giác 3 nút  
cho tấm Reissner-Mindlin.  
N1 1; N2 ; N3   
wi là chuyển vị thẳng theo phương  
nút thứ i.  
z
của  
θxi, θyi lần lượt là góc xoay của đoạn thẳng  
vuông góc với mặt trung bình quay quanh  
trục x và trục y có chiều như Hình 3.1 của  
nút thứ i.  
Công thức (3.1) biểu diễn dưới dạng ma  
trận:  
Hình 3.1: Phần tử tam giác 3 nút.  
   
w
x  
y  
N
0
0
0
w
   
Xét trên mặt trung bình của tấm, khi ta  
thực hiện rời rạc bề mặt này thành các  
phần tử tam giác 3 nút để tính toán thì mỗi  
nút phần tử sẽ được xác định gồm 3 thành  
phần: chuyển vị w và 2 góc xoay (βx, βy).  
Mỗi phần tử tam giác tính toán sẽ gồm 3  
nút. Mỗi nút sẽ có 3 bậc  
i
3
(3.3)  
0
0
Ni x  
  
   
i1   
   
y  
Ni  
0
   
i
N
i
Từ công thức (2.8), vectơ biến dạng uốn  
được viết lại dưới dạng ma trận:  
x  
x  
y  
tự do. Giả thiết với vật liệu đồng nhất và  
đẳng hướng có module đàn hồi E và hệ số  
poisson .  
  
y  
  
x  
y  
y   
x  
(3.4)  
Xấp xỉ trường chuyển vị:  
i   
N  
x  
0
0
   
w
3
Ni  
y  
Ni N  
   
0   
0   
0
x  
y  
  
   
i1  
   
 i  
i   
x  
y  
Bbi  
Hình 3.2: Hệ tọa độ tổng thể.  
3
Trong đó, Bbi là ma trận quan hệ giữa biến  
dạng uốn và chuyển vị và được xác định  
như sau:  
w N w  
i
i
i1  
3
(3.1)  
N   
x
i
yi  
i1  
B B Bb2 Bb3  
b
b1  
3
  N   
y
i
xi  
0
b c  
0
0
0
0
c
0
0
d  
c
0
0
0
b  
0
0
a
i1  
(3.5)  
1
0
0
0
d a  
2A  
Trong đó: Ni là các hàm dạng trong hệ tọa  
độ tự nhiên (,).  
e
d a b c  
d  
a
b  
7
Với Ae là diện tích của phần tử tam giác có  
tọa độ nút (x1,y1), (x2,y2), (x3,y3) như Hình  
3.3 và  
Với  
1
0
0
Eh3  
12(12  
Cb   
1
(3.11)  
)
a x2 x1; b y2 y1;  
(3.6)  
1  
0 0  
c y3 y1; d x3 x1  
2   
1 0  
0 1  
Eh  
Cs   
(3.12)  
2(1  
)
3.2. Kỹ thuật MITC3.  
Ý tưởng cơ bản của kỹ thuật MITC3 dùng  
để khử hiện tượng khóa cắt cho phần tử  
tấm tam giác 3 nút là xấp xỉ lại các biến  
dạng cắt ngoài mặt phẳng sao cho các biến  
dạng cắt này là hằng số trên cạnh phần tử.  
Lee và Bathe [20] đã xây dựng các hàm  
xấp xỉ lại biến dạng cắt ngoài mặt phẳng  
thông qua giá trị biến dạng cắt tính từ xấp  
xỉ chuyển vị (3.1) tính tại các điểm buộc.  
Từ đó, biến dạng cắt ngoài mặt phẳng xấp  
xỉ theo kỹ thuật MITC3 có thể được biểu  
diễn theo chuyển vị nút phần tử.  
Hình 3.3: Tọa độ phần tử tam giác 3 nút  
Tương tự, từ công thức (2.9), vectơ biến  
dạng cắt được viết lại dưới dạng ma trận:  
w  
x  
w  
i
x  
  
y   
y  
(3.7)  
N  
x  
N  
y  
   
w
0
Ni  
3
   
x  
   
i1   
i
Ni  
0
   
y  
 i  
B
si  
Theo phương pháp PTHH truyền thống, từ  
các quan hệ giữa biến dạng và chuyển vị,  
hệ phương trình rời rạc của phần tử tấm  
tam giác 3 nút dùng phân tích tĩnh tấm  
Reissner–Mindlin được viết dưới dạng:  
Ku = F (3.8)  
Để có thể dễ dàng áp dụng các kỹ thuật  
làm trơn cho biến dạng cắt ngoài mặt  
phẳng xấp xỉ theo kỹ thuật MITC3, Chau-  
Dinh và cộng sự [21] đã tính toán ma trận  
quan hệ giữa biến dạng cắt ngoài mặt  
phẳng theo kỹ thuật MITC3 và chuyển vị  
nút phần tử theo tọa độ nút phần tử bằng  
cách sử dụng 1 điểm Gauss.  
Trong đó, u là véc-tơ chuyển vị nút kết  
cấu, F là véc-tơ tải trọng tại các nút của kết  
cấu và được lắp ghép từ các véc-tơ tải  
trọng tại nút phần tử  
Do đó, công thức PTHH tấm tam giác 3  
nút MITC3 có biến dạng cắt (3.7) được  
T
fe N 0 0 N 0 0 N 0 0 qdA (3.9)  
1
2
3
   
w
A
e
3
   
thay bằng MITC3  
BMITC3 x (3.13)  
ke BTb CbBbdABTs CsBsdA  
si  
   
i1  
   
y  
A
A
 i  
e
e
8
Với:  
(b-c)(bc)  
6
(d-a)(bc)  
(b-c)  
A   
e
1
2A  
MITC3  
6
B
B
s1  
(b-c)(ad)  
(d-a)(ad)  
e
(d-a) -A -  
-
e
6
6
bc c (bc) ac  
d
6
c
(bc)  
1
2
2
3
2
MITC3  
s2  
2A  
bd  
2
c
6
ad  
d
6
e
d  
(ad)  
(ad)  
Hình 3.4: Chia lưới phần tử tam giác làm  
trơn trên cạnh và trên nút.  
Đối với miền làm trơn kết hợp giữa miền  
2
bc b  
bd  
2
a
6
a
6
b  
(bc)  
(bc)  
1
2A  
MITC3  
s3  
2
6
B
ac  
2
b
6
ad  
2
e
a
(ad)  
(ad)  
trơn trên nút và miền làm trơn trên cạnh,  
các phần tử tam giác 3 nút trong miền  
chuẩn PTHH được chia thành miền bao  
quanh một nút để tạo thành miền làm trơn  
(3.14)  
Và công thức ma trận độ cứng phần tử  
(3.10) được viết lại cho phần tử MITC3:  
T
MITC3  
e
T
b
MITC3  
s
MITC3  
s
k
ˆ
k
B
C B dA  
B
C
B
dA  
(3.15)  
trên nút NS (miền chấm chấm trên Hình  
b
b
s
A
A
e
e
3.5) và phần còn lại chung cạnh giữa 2  
Trong đó, BMs ITC3 Bs1  
B
B
s3  
MITC3  
MITC3  
s2  
MITC3  
ES  
phần tử là miền làm trơn trên cạnh  
3.3. Công thức phần tử ES+NS-MITC3.  
(miền gạch chéo trên Hình 3.5).  
Xét 1 tấm được chia thành các phần tử tam  
giác 3 nút. Theo phương pháp PTHH trơn,  
các trường biến dạng được trung bình trên  
miền làm trơn giữa các phần tử chung nút  
k
kES  
hoặc giữa 2 phần tử chung cạnh  
NS  
như Hình 3.4. Trong đó, miền làm trơn  
trên nút được xác định bằng cách nối  
Hình 3.5: Miền làm trơn dựa trên lưới tam  
giác cho mô hình ES+NS-MITC3  
kNS  
Để điều chỉnh tỉ lệ giữa 2 miền làm trơn,  
trọng tâm với trung điểm các cạnh của các  
phần tử chung nút (Hình 3.4a). Miền làm  
hệ số hiệu chỉnh được định nghĩa nhằm  
xác định tỉ lệ chiều dài trên cạnh phần tử  
tham gia vào miền làm trơn trên nút (NS)  
và trên cạnh (ES) như sau  
kES  
trơn trên cạnh  
được định nghĩa bằng  
cách nối trọng tâm của 2 phần tử chung  
cạnh với 2 đỉnh của cạnh chung (Hình  
L
AEkS  
ANkS  
lần lượt là diện tích  
L L 1  
; L2   
L
(3.16)  
3.4b). Gọi  
và  
1
3
2
k
của miền làm trơn trên cạnh  
và trên  
ES  
Trong đó, L là chiều dài cạnh phần tử, L1,  
L3 là chiều dài của cạnh phần tử tham gia  
vào miền làm trơn trên nút, và L2 là chiều  
k
nút  
.
NS  
9
dài của cạnh phần tử tham gia vào miền  
làm trơn trên cạnh.  
Và kỹ thuật làm trơn trên nút  
e
N
k
1
1 i  
A B  
e bi  
ˆ
B
i1  
;
Khi đó, miền làm trơn trên cạnh ES có  
bi  
2
k
3
1  
A
NS  
(3.20)  
diện tích  
e
k
N
1
2
MITC3  
1 i MITC3  
A B  
ˆ
B
1
si  
e si  
k
AES   
AEkS  
3
1  
A
i1  
(3.17)  
NS  
3
Ai  
Trong đó,  
là diện tích của phần tử thứ i  
k
e
ˆ
Và miền làm trơn trên nút NS có diện  
Nke  
tham gia vào miền làm trơn và  
là số  
tích  
phần tử trong miền làm trơn k.  
1
k
ANS 12 Ak  
ˆ
(3.18)  
NS  
Kết quả, ma trận độ cứng của phần tử  
3
ES+NS-MITC3 được xác định như sau  
Làm trơn các biến dạng uốn và cắt ngoài  
keESNS kee  
(3.21)  
MITC3  
mặt phẳng của phần tử tấm tam giác 3 nút  
MITC3 trên các miền làm trơn  
và  
ES  
Với:  
T
k
MITC3  
T
b
MITC3  
MITC3  
ˆ
k
C
(3.22)  
(3.23)  
NS bằng cách trung bình biến dạng giữa  
e
b
b
s
s
s
s
A
A
e
e
các  
T
MITC3  
e
T
b
MITC3  
s
MITC3  
s
ˆ
ˆ
ˆ
ˆ
ˆ
k
B
C B dA  
B
B
dA  
b
b
A
A
phần tử nằm trong miền làm trơn với diện  
tích miền làm trơn, các ma trận quan hệ  
giữa biến dạng và chuyển vị nút phần tử  
được cho bởi kỹ thuật làm trơn trên cạnh  
e
e
Như vậy, từ công thức (3.9) ma trận độ  
cứng phần tử ES+NS-MITC3 được lắp  
ghép vào ma trận độ cứng kết cấu thông  
qua các bậc tự do tại nút đỉnh phần tử.  
Thế (3.21) vào phương trình cân bằng  
(3.8), sử dụng ngôn ngữ lập trình Matlab  
để lập trình tính toán chuyển vị và nội lực  
của tấm.  
e
N
k
1
1 i  
A B  
e bi  
B
i1  
;
bi  
2 k  
3
A  
ES  
(3.19)  
e
k
N
1
MITC3  
si  
1 i MITC3  
A B  
B
2 k  
A  
ES  
e si  
3
i1  
ES-MITC3 [16], NS-MITC3 [17], ES-  
DSG3 [15], QHS [22]. Tham số hiệu chỉnh  
tỉ lệ tham gia của miền làm trơn trên nút và  
trên cạnh được chọn = 0.6.  
4. Ví dụ số  
Bài báo thực hiện việc phân tích tĩnh một  
số bài toán tấm đồng nhất đẳng hướng, sau  
đó so sánh kết quả chuyển vị và mô-men  
cho bởi phần tử ES+NS-MITC3 với kết  
quả tính toán của các phần tử tấm tam giác  
3 nút tham khảo gồm phần tử MITC3 [12],  
10  
tỉ số h/L = 0.1 hoặc 0.001, và chịu tải trọng  
phân bố đều q = 1. Tấm được làm bằng vật  
liệu đồng nhất đẳng hướng, có các thông  
số vật liệu: mô-đun đàn hồi E = 1092000  
và hệ số Poisson = 0.3.  
4.1. Bài toán kiểm tra.  
1 (0.2; 0.12)  
4 (0.44; 0.12)  
5 (0.3; 0.6)  
Tấm được mô hình thành bằng các phần tử  
tam giác 3 nút với lưới chia đều N×N phần  
tử như Hình 4.3. Trong đó, N là số phần tử  
chia trên mỗi cạnh của tấm. Để khảo sát sự  
hội tụ của phần tử đề xuất, các lưới phần  
tử khác nhau với N = 4, 8, 10, 12, 16 được  
khảo sát. Nhằm mục đích so sánh với lời  
giải giải tích và lời giải tham khảo [23],  
kết quả độ võng và mô-men không thứ  
nguyên tại tâm tấm được tính theo công  
thức sau:  
2 (0.2; 0)  
3 (0.44; 0)  
Hình 4.1: Bài toán kiểm tra  
Bảng 4.1: Kết quả bài toán kiểm tra  
Phần tử  
ES+NS-  
MITC3  
Exact  
w5  
θx5  
θy5  
mx5  
my5  
mxy5  
0,6422  
0,6422  
1,13 – 0,64 – 0,0111 – 0,0111 – 0,0033  
1,13 – 0,64 – 0,0111 – 0,0111 – 0,0033  
Trong bài toán này, phần tử ES+NS-  
MITC3 được kiểm tra khả năng tìm lại kết  
quả chuyển vị và nội lực của tấm chữ nhật  
có chiều dày h = 0.01 và có kích thước  
như Hình 4.1. Tấm làm bằng vật liệu đồng  
nhất đẳng hướng có mô-đun đàn hồi E =  
105 và hệ số Poisson = 0.3.  
4
;
w w / qL /100D  
c
c
2
Tấm được rời rạc thành các phần tử tam  
giác 3 nút có tọa độ như Hình 4.1. Kết quả  
độ võng, góc xoay và mô-men tại nút 5 tọa  
độ x = 0.3 và y = 0.6 cho bởi phần tử  
ES+NS-MITC3 và lời giải chính xác được  
cho trong Bảng 4.1 cho thấy phần tử  
NS+ES-MITC3 có khả năng tìm lại được  
chính xác giá trị chuyển vị và nội lực của  
bài toán kiểm tra.  
M M / qL /10  
c
c
Trong đó, wc và Mc lần lượt là độ võng và  
Eh3  
12(1v2 )  
mô-men tại tâm tấm, và  
D   
là độ cứng trụ của tấm.  
4.2. Tấm hình vuông ngàm và tấm tựa  
đơn 4 cạnh chịu tải phân bố điều.  
Xét ví dụ tấm hình vuông ở 2 trường hợp  
ngàm 4 cạnh như Hình 4.2(a) hoặc tựa đơn  
4 cạnh như Hình 4.2(b). Cả 2 trường hợp  
tấm vuông cạnh chiều dài L, độ dày h với  
11  
Hình 4.5: Kết quả hội tụ chuẩn hóa độ  
võng tại tâm tấm dày vuông ngàm với tỉ lệ  
h/L = 0.1  
Hình 4.2: (a) Tấm vuông đồng nhất đẳng  
hướng biên ngàm 4 cạnh; (b) Tấm vuông  
đồng nhất đẳng hướng biên tựa đơn 4  
cạnh.  
Phần tử ES+NS-MITC3 cho kết quả độ  
võng ở tấm mỏng tốt hơn tấm dày. Tuy  
nhiên ở trường họp này chỉ vượt trôi hơn  
phần tử NS-MITC3 nhưng lại kém hơn  
phần tử ES-MITC3, ES–DSG3, MITC3.  
Hình 4.3: Chia lưới phần tử hữu hạn cho  
tấm.  
Kết quả độ võng và mô-men không thứ  
nguyên tại tâm tấm ứng với các loại lưới  
khác nhau cho bởi phần tử đề xuất và các  
phần tử ES-MITC3, NS-MITC3, ES-  
DSG3, MITC3  
Hình 4.6: Kết quả hội tụ chuẩn hóa  
moment tại tâm tấm mỏng vuông ngàm  
với tỉ lệ h/L = 0.001  
Hình 4.7: Kết quả hội tụ chuẩn hóa  
moment tại tâm tấm dày vuông ngàm với tỉ  
lệ h/L = 0.1  
Hình 4.4: Kết quả hội tụ chuẩn hóa độ  
võng tại tâm tấm mỏng vuông ngàm với tỉ  
lệ h/L = 0.001  
Đối với cả tấm dày và mỏng từ lưới 8×8  
đến 16×16 phần tử ES+NS-MITC3 cho kết  
quả mô-men tại tâm tấm khá tốt và hội tụ  
nhanh đến lời giải chính xác khi chia lưới  
12  
càng mịn. Ở tấm dày Hình 4.7, phần tử  
ES+NS-MITC3 cho kết quả mô-men hội  
tụ tương đương phần tử NS-MITC3 và  
vượt qua phần tử ES-MITC3, ES-DSG3,  
MITC3. Khi tấm mỏng Hình 4.6, phần tử  
ES+NS-MITC3 cho kết quả mô men hội tụ  
tương đương phần tử ES-MITC3 và vượt  
qua phần tử NS-MITC3, ES-DSG3,  
MITC3. Tóm lại, phần tử ES+NS-MITC3  
cho kết quả mô-men rất tốt khi giải tấm  
dày hoặc tấm mỏng.  
Hình 4.10: Kết quả hội tụ chuẩn hóa  
moment tại tâm tấm dày vuông tựa đơn  
với tỉ lệ h /L = 0.001  
Khi sử dụng phần tử ES+NS-MITC3 giải  
bài toán tấm vuông đồng nhất đẳng hướng  
có biên tựa đơn 4 cạnh chịu tải trọng phần  
bố đều trong trường hợp tấm mỏng lẫn tấm  
dày nhìn chung đều cho ra kết quả độ võng  
khá tốt, cho ra kết quả chính xác hơn các  
phần tử ES-MITC3, NS-MITC3, ES-  
DSG3, và MITC3.  
Hình 4.11: Kết quả hội tụ chuẩn hóa  
moment tại tâm tấm mỏng vuông tựa đơn  
với tỉ lệ h /L = 0.1  
Hình 4.8: Kết quả hội tụ chuẩn hóa độ  
võng tại tâm tấm mỏng vuông tựa đơn  
với tỉ lệ h/L = 0.001  
Bên cạnh đó ở Hình 4.10- 4.11 phần tử  
ES+NS-MITC3 cho kết quả mô-ment tại  
tâm tấm với 2 loại lưới thô và lưới mịn hội  
tụ như nhau và gần trùng với lời giải chính  
xác và là phần tử cho kết quả tốt nhất so  
với các phần tử ES-MITC3, NS-MITC3,  
ES-DSG3, và MITC3.  
Hình 4.9: Kết quả hội tụ chuẩn hóa độ  
võng tại tâm tấm dày vuông tựa đơnvới tỉ  
lệ h /L = 0.1  
4.3. Tấm hình thoi tựa đơn ở cạnh trên  
và dưới, 2 cạnh bên tự do và chịu tải  
phân bố đều (tấm Razzaque).  
Ở ví dụ này, ta xét tấm hình thoi Razzaque  
[24] đồng nhất đẳng hướng nghiên một  
góc = 600 so với trục x. Tấm chịu tác  
dụng của tải phân bố đều q = 1 như Hình  
4.12. Tấm có cạnh kích thước L = 100,  
13  
chiều dày h = 0.1, tựa đơn cạnh trên và  
cạnh dưới, và 2 cạnh bên tự do. Vật liệu  
làm tấm có mô-đun đàn hồi E = 10.92 và  
hệ số Poisson = 0.3.  
phần tử ES+NS-MITC3 cho kết quả tốt  
hơn phần tử ES-MITC3, NS-MITC3,  
MITC3, và kết quả độ hội tụ gần như trùng  
với lời giải Razzaque.  
Tấm cũng được mô hình thành các phần tử  
tam giác, với các lưới 4×4, 8×8, 12×12,  
16×16, và 24×24 như Hình 4.13. Kết quả  
4
độ võng chuẩn hóa  
w w 100D / qL  
c
  
c
và  
mô-men  
chuẩn  
hóa  
Hình 4.15: Kết quả hội tụ chuẩn hóa mô-  
men tại tâm tấm hình thoi Razzaque với tỉ  
lệ h /L = 0.001  
2
tại tâm tấm cho bởi  
M M / 10/ qL  
y
y
phần tử ES+NS-MITC3 và tài liệu tham  
Qua Hình 4.14 và Hình 4.15 ta thấy kết  
khảo [25]  
quả hội tụ chuẩn hóa độ võng và mô-men  
tại tâm tấm hình thoi phần tử ES+NS-  
MITC3 cho kết quả chuẩn hóa độ võng và  
mô-men gần trùng với lời giải Razzaque.  
Phương pháp ES+NS-MITC3 cải thiện  
đáng kể độ chính xác của độ võng và  
moment so với các phương pháp ES-  
MITC3, NS-MITC3, và MITC3.  
Hình 4.12: Mô hình tấm hình thoi  
Razzaque  
4.4. Tấm hình thoi tựa đơn 4 cạnh chịu  
tải phân bố đều (tấm Morley).  
Hình 4.13: Các cách chia lưới N×N của  
tấm hình thoi Razzaque thành các phần tử  
tam giác 3 nút (với N = 4, 8, 12, 16, 24).  
Tiếp theo, ta xét tấm đồng nhất đẳng  
hướng hình thoi Morley tựa đơn chịu tải  
phân bố đều như Hình 4.16. Lý thuyết tấm  
này được nghiên cứu bởi Morley [26].  
Tấm hình thoi nghiêng một góc = 300 so  
với trục x, cạnh L = 100, chiều dày h = 0.1,  
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 = 10.92 và hệ số Poisson = 0.3.  
Hình 4.14: Kết quả hội tụ chuẩn hóa độ  
võng tại tâm tấm hình thoi Razzaque  
với tỉ lệ h/L = 0.001  
Từ Hình 4.14 cho thấy giá trị độ võng đạt  
được của tấm hình thoi Razzaque bằng  
14  
Tương tự như ví dụ 4.2, tấm được chia đều  
thành các phần tử tam giác 3 nút, với các  
loại lưới 4×4, 8×8, 10×10, 12×12, và  
16×16.  
MITC3. Điều này cho thấy, phần tử  
ES+NS-MITC3 có hiệu quả rất tốt khi giải  
bài toán tấm Morley.  
4.5. Tấm hình tròn ngàm chịu tải phân bố  
đều.  
Cho tấm tròn bán kính R = 5, chiều dày h  
với h/R = 0.2 hoặc 0.02, ngàm trên biên và  
chịu tải phân bố đều q = 1 như Hình 4.18.  
Tấm được làm bằng vật liệu đồng nhất  
đẳng hướng có mô-đun đàn hồi E = 10.92  
và hệ số Poisson = 0.3.  
Hình 4.16: Mô hình tấm hình thoi Morley  
Sử dụng tính chất đối xứng, 1/4 hình tấm  
được mô hình bằng các phần tử tam giác 3  
nút với các loại lưới như Hình 4.20.  
Hình 4.17: Kết quả hội tụ chuẩn hóa độ  
võng tại tâm tấm hình thoi Morley tựa đơn  
với tỉ lệ h/L = 0.001  
Hình 4.19: Tấm tròn với điều kiện biên  
ngàm  
Hình 4.18: Kết quả hội tụ chuẩn hóa  
moment tại tâm tấm hình thoi Morley tựa  
đơn với tỉ lệ h/L = 0.001  
Từ Hình 4.17 và Hình 4.18 ta thấy độ võng  
và moment tấm hình thoi Morley tựa đơn  
sử dụng phần tử ES+NS-MITC3 cho kết  
quả độ võng hội tụ nhanh hơn phần tử ES-  
DSG3, NS-MITC3, MITC3. Tuy nhiên,  
kết quả của phần tử đề xuất kém hơn phần  
tử ES-MITC3. Mô-men chuẩn hóa tại tâm  
tấm của phần tử ES+NS-MITC3 lại hội tụ  
nhanh hơn giá trị mô-men cho bởi phần tử  
ES-MITC3, MITC3, ES-DSG3, và NS-  
Hình 4.20: Mô hình chia lưới 1/4 tấm tròn  
có sử dụng phần tử tam giác 3 nút  
Kết quả khảo sát độ võng và mô-men tại  
tâm tấm khi tỉ số h/R thay đổi với các loại  
lưới khác nhau, kết quả so sánh độ võng và  
mô-men cho bởi phần tử ES+NS-MITC3  
và các kết quả tham khảo [27] được thể  
hiện qua Hình 4.21-4.24.  
15  
Hình 4.23: Kết quả hội tụ mô-men  
của tấm tròn ngàm (h/R = 0.02),  
chịu tải phân bố đều.  
Hình 4.21: Kết quả hội tụ của độ võng của  
tấm tròn ngàm (h/R = 0.02), chịu tải phân  
bố đều.  
Hình 4.24: Kết quả hội tụ mô-men của tấm  
tròn ngàm (h/R = 0.2),  
chịu tải phân bố đều.  
Từ Hình 4.21-4.24 cho thấy tấm tròn đồng  
nhất đẳng hướng có biên ngàm giải bằng  
phần tử ES+NS-MITC3 cho ra kết quả rất  
tốt đối với tấm mỏng và tấm dày. Trong  
trường hợp lưới thô hoặc lưới mịn, độ  
võng và mô-men đều hội tụ gần với lời  
giải chính xác và hội tụ nhanh hơn các  
phần tử ES-DSG3, MITC3, NS-MITC3,  
và ES-MITC3.  
Hình 4.22: Kết quả hội tụ của độ võng của  
tấm tròn ngàm (h/R = 0.2),  
chịu tải phân bố đều.  
5. Kết luận  
Từ những kết quả tính toán độ võng và  
mô-men các bài toán tấm đồng nhất đẳng  
hướng với hình dạng và điều kiện biên  
khác nhau và so sánh với các phần tử tấm  
tam giác 3 nút cùng loại, có thể rút ra một  
số kết luận như sau: Phần tử ES+NS-  
MITC3 cho kết quả tương đối tốt, khá gần  
với lời giải chính xác đối với một số loại  
tấm khác nhau như tấm hình vuông, tấm  
hình thoi Razzaque, Morley và tấm tròn có  
điều kiện biên tựa đơn hoặc ngàm. Đặc  
biệt, kết quả độ võng và mô-men cho bởi  
phần tử ES+NS-MITC3 rất tốt đối với tấm  
tròn, tấm hình thoi Morley, Razzaque, và  
tấm vuông tựa đơn chịu tải phân bố đều.  
Ngoại trừ trường hợp tấm vuông ngàm cho  
kết quả hội tụ chưa tốt nhưng so với phần  
tử NS-MITC3, ES-DSG3, MITC3, QHS  
thì phần tử ES+NS-MITC3 cải thiện đáng  
kể độ hội tụ và chính xác về độ võng và  
16  
mô-men tại điểm khảo sát. Tuy nhiên, so  
với các loại phần tử chỉ sử dụng 1 kỹ thuật  
làm trơn thì phần tử ES+NS-MITC3 tốn  
nhiều thời gian hơn để tính toán ma trận độ  
cứng phần tử do phải xác định miền làm  
trơn trên cạnh (ES) và trên nút (NS) phần  
tử.  
Tài liệu tham khảo  
[1] K. T. S. Oldham, The doctrine of  
description: Gustav Kirchhoff, classical  
physics, and the “purpose of all science”  
in 19th-century Germany. University of  
California, Berkeley, 2008.  
[6]  
K.-J. Bathe and E. N. Dvorkin, “A  
formulation of general shell elements—  
the use of mixed interpolation of  
tensorial components,” International  
Journal for Numerical Methods in  
Engineering, vol. 22, no. 3, pp. 697–722,  
Mar. 1986.  
[2] R. D. Mindlin, ‘Influence of Rotary  
Inertia and Shear on Flexural Motions of  
Isotropic, Elastic Plates,’ Journal of  
Applied Mechanics, Vol. 18, No. 1, 1951,  
[7] K.-J. Bathe, P.-S. Lee, and J.-F. Hiller,  
“Towards improving the MITC9  
pp.  
31-  
shell element,” Computers  
&
38. - References - Scientific Research  
Publish.”  
Structures, vol. 81, no. 8, pp. 477–  
489, May 2003.  
[3] G. R. Liu and N. T. Trung, Smoothed  
Finite Element Methods, 1 edition.  
Boca Raton: CRC Press, 2010.  
[8] M. L. Bucalem and K.-J. Bathe,  
“Higher-order MITC general shell  
elements,” International Journal for  
Numerical Methods in Engineering,  
vol. 36, no. 21, pp. 3729–3754, Nov.  
1993.  
[4] E. Dvorkin and K.-J. Bathe, “A  
Continuum mechanics based four-  
node shell element for general  
nonlinear analysis,” Engineering  
Computations, vol. 1, pp. 77–88,  
Dec. 1984.  
[9] Y. Ko, P.-S. Lee, and K.-J. Bathe, “A  
new MITC4+ shell element,”  
Computers & Structures, vol. 182,  
pp. 404–418, Apr. 2017.  
[5] K.-J. Bathe and E. N. Dvorkin, “A  
four-node plate bending element  
based on Mindlin/Reissner plate  
theory and a mixed interpolation,”  
International Journal for Numerical  
Methods in Engineering, vol. 21, no.  
2, pp. 367–383. 1984.  
[10] K.-U. Bletzinger, M. Bischoff, and E.  
Ramm, “A unified approach for  
shear-locking-free triangular and  
rectangular shell finite elements,”  
Computers & Structures, vol. 75, no.  
3, pp. 321–334, Apr. 2000.  
17  
[11] A. Tessler and T. J. R. Hughes, “A  
three-node mindlin plate element  
with improved transverse shear,”  
Computer Methods in Applied  
Mechanics and Engineering, vol. 50,  
no. 1, pp. 71–101, Jul. 1985.  
[16] Q. Nguyen-Duy, ''Phân tích kết cấu  
tấm bằng phần tử biến dạng trơn ES-  
MITC3,'' luận văn thạc sĩ, Trường  
ĐH Sư phạm Kỹ thuật Tp.HCM.  
[17] D. Nguyen-Van ''Phân tích kết cấu  
tấm bằng phần tử biến dạng trơn NS-  
MITC3’’,luận văn thạc sĩ, Trường  
ĐH Sư phạm Kỹ thuật Tp.HCM,  
2016.  
[12] H.-M. Jeon, P.-S. Lee, and K.-J.  
Bathe, “The MITC3 shell finite  
element enriched by interpolation  
covers,” Computers & Structures,  
vol. 134, pp. 128–142, Apr. 2014.  
[18] F. Wu, W. Zeng, L. Y. Yao, M. Hu,  
Y. J. Chen, and M. S. Li, “Smoothing  
Technique Based Beta FEM (βFEM)  
for Static and Free Vibration  
[13] Y. Lee, P.-S. Lee, and K.-J. Bathe,  
"The MITC3+ shellelement and its  
performance  
Computers  
&
Analyses  
of  
Reissner–Mindlin  
Structures'', vol. 138, pp. 12-23, Jul.  
2014.  
Plates,” Int. J. Comput. Methods, p.  
1845006, Mar. 2018.  
[14] Nguyen – Thoi T, Liu GR, Nguyen –  
Xuan H. Additional properties of the  
node – based smoothed finite element  
method (NS – FEM) for solid  
mechanics problems. International  
Journal of Computational Methods,  
6(4): 633 – 666, 2009.  
[19] E. Li and Z. He ''Development of a  
perfect match system in the  
improvement of eigenfrequencies of  
free  
vibration,”  
Applied  
Mathematical Modelling, vol. 44, pp.  
614–639, Apr. 2017.  
[20] P.-S. Lee and K.-J. Bathe,  
“Development of MITC isotropic  
triangular shell finite elements,”  
Computers & Structures, vol. 82, no.  
11, pp. 945–962, May 2004.  
[15] P. Phung-Van, Chien H. Thai, T.  
Nguyen-Thoi, H. Nguyen-Xuan,  
''Static and free vibration analyses of  
composite and sandwich plates by an  
edge-based smoothed discrete shear  
gap method (ES-DSG3) using  
[21] T. Chau-Dinh, Q. Nguyen-Duy, and  
H. Nguyen-Xuan, “Improvement on  
MITC3 plate finite element using  
triangular elements  
based on  
layerwise theory'', Composites Part  
B: Engineering, Volume 60, Pages  
227–238, April 2014.  
edge-based  
strain  
smoothing  
enhancement for plate analysis,” Acta  
Mechanica, vol. 228, Feb. 2017.  
2
[22] Chen Wanji and Y. K. Cheung.  
Refined 9-Dof triangular Mindlin  
plate elements. International journal  
for Static and Free Vibration  
Analyses of Reissner–Mindlin  
Plates,'' Int. J. Comput. Methods, p.  
1845006, Mar. 2018.  
for  
engineering; 51:1259–1281. 2001.  
[23] S. P. Timoshenko and  
Woinowsky-Krieger,  
numerical  
methods  
in  
[26] L. S. D. Morley, Skew plates and  
structures. Oxford; New York:  
Pergamon Press; distributed in the  
Western Hemisphere by Macmillan,  
New York, 1963.  
S.  
Theory of  
Plates and Shells, Second Edition.  
McGraw-Hill, 1959.  
[24] A.  
Razzaque,  
“Program  
for  
[27] R. Ayad, A. Rigolot, and N. Talbi,  
“An improved three-node hybrid-  
mixed element for Mindlin/Reissner  
plates,” International Journal for  
Numerical Methods in Engineering,  
vol. 51, no. 8, pp. 919–942, 2001.  
triangular bending elements with  
derivative smoothing,” International  
Journal for Numerical Methods in  
Engineering, vol. 6, no. 3, pp. 333–  
343, Jan. 1973.  
[25] J. Chen, and M. S. Li, ''Smoothing  
Technique Based Beta FEM (βFEM)  
Tác giả chịu trách nhiệm bài viết:  
Họ tên: Hồ Thị Còn  
Đơn vị: Trung Tâm Giáo Dục Nghề Nghiệp – Giáo Dục Thường Xuyên,  
Huyện An Phú, Tỉnh An Giang  
Điện thoại: 0939550094  
3
pdf 18 trang Mãnh Khiết 12/01/2024 19300 Free
Bạn đang xem tài liệu "Phân tích tấm reissner-mindlin bằng phần tử es+ns-mitc3", để tải tài liệu gốc về máy hãy click vào nút Download ở trên

File đính kèm:

  • pdfphan_tich_tam_reissner_mindlin_bang_phan_tu_esns_mitc3.pdf