Chương 5
BÀI TOÁN KHÔNG GIAN
5.1. Sơ đồ tính. Phần tử tứ diện
Trong phần này sẽ giới thiệu các quan hệ và công thức cơ bản nhất khi sử dụng
phương pháp PTHH để giải bài toán không gian của Lý thuyết đàn hồi dựa trên các hiểu
biết đã đề cập ở các phần trước. Đối tượng của bài toán là những vật thể hình khối,
trong đó các trường ứng suất, biến dạng và chuyển vị được biểu diễn bằng những hàm
của ba biến x,
Vật thể khối đàn hồi liên tục (dầm lớn, ống vỏ dày, nền móng, thân máy) được
rời rạc hóa thành một tập hợp các phần tử hình khối (phần tử ba chiều). Các phần tử liên
kết với nhau tại các nút bằng các khớp không gian. Khi theo một phương nào đó không
có chuyển vị thì người ta gắn theo phương đó một gối khớp không gian. Có nhiều loại
phần tử khác nhau, với số nút và bậc tự do khác nhau. Hai loại phần tử thường hay dùng
là phần tử tứ diện có 4 điểm nút và phần tử lục diện có 8 điểm nút (hình 5.1). Ngoài ra
sử dụng phần tử bậc cao và phần tử đẳng tham số cũng rất có hiệu quả khi giải bài toán
loại này.
So với bài toán phẳng thì bài toán không gian có khối lượng tính toán lớn hơn rất
nhiều. Nếu chọn hàm xấp xỉ của chuyển vị có dạng giống như bài toán phẳng thì số nút
tăng lên dẫn đến số bậc tự do tăng và bề rộng dải của ma trận độ cứng cũng tăng. Thí
dụ, Nếu trong bài toán phẳng lấy mạng lưới 10 x10 có 100 nút, số phương trình là 200
và bề rộng dải là 20, thì trong bài toán không gian tương đương, với mạng lưới 10 x 10
x10, số nút là 1000, số phương trình là 3000 và bề rộng dải là 300. Nếu lấy mạng lưới
30 x30 thì trong bài toán phẳng số phương trình là 1800, bề rộng dải 60, còn trong bài
toán không gian là 81000 phương trình, bề rộng dải 2700. Nghĩa là khối lượng tính toán
tăng lên rất nhiều. Vì vậy chọn dạng phần tử hợp lý trong bài toàn không gian là rất
quan trọng.
Sau đây giới thiệu cách sử dụng phần tử tứ diện tuyến tính và phần tử đẳng
tham số.
Hình 5.1
Ta xét một phần tử tứ diện có 4 nút i, j, m, p (hình 5.2). Vectơ chuyển vị nút có
dạng
{δ } = ui
vi
wi
uj
vj
wj
um
vm
wm
up
vp
wp
T
(5.1)
Vectơ tải phần tử có dạng
{P}
e
= U i Vi Wi U j V j W j U m Vm Wm U p V p W p
T
(5.2)
Để tìm ứng suất sau khi xác định được chuyển vị nút ta vẫn sử dụng quan hệ:
{σ } = [ D ][ B ]{δ } = [ S ]{δ }
trong đó vectơ ứng suất gồm 6 thành phần:
{σ } = σ x σ y σ z τ xy τ yz τ zx
T
5.2. Hàm chuyển vị
Hình 5.2
Đặt phần tử tứ diện ijmp trong hệ tọa độ không gian xyz. Trạng thái chuyển vị của
một điểm bất kỳ trong phần tử được xác định bởi 3 thành phần chuyển vị u, v, w tương
ứng với 3 phương x,y,z. Ta có
u
{ f } = v
w
(5.3)
Ta biết rằng vị trí của tứ diện phụ thuộc hoàn toàn vào vị trí 4 điểm nút, tức là phụ
thuộc vào 12 thành phần chuyển vị nút. Do đó có thể giả thiết hàm chuyển vị có dạng
sau đây:
u = α1 + α 2 x + α 3 y + α 4 z
v = α5 + α 6 x + α7 y + α8 z
(5.4)
w = α 9 + α10 x + α11 y + α12 z
Tại các nút i, j, m, p ta có
PPPTHH 115
ui = α1 + α 2 xi + α 3 yi + α 4 zi
u j = α1 + α 2 x j + α 3 y j + α 4 z j
um = α1 + α 2 xm + α 3 ym + α 4 zm
(5.5)
u p = α1 + α 2 x p + α 3 y p + α 4 z p
Giải hệ phương trình trên đối với α1 , α 2 , α 3 , α 4 , sau đó thay lại vào (5.4) ta được
biểu thức liên hệ giữa chuyển vị u và các chuyển vị nút:
1
(ai + bi x + ci y + d i z )ui − (a j + b j x + c j y + d j z )u j
6V
+ (am + bm x + cm y + d m z )um − (a p + bp x + c p y + d p z )u p
u=
trong đó:
(5.6)
V là thể tích của tứ diện ijmp:
1
V=
xi
yi
zi
1 1 xj
6 1 xm
1 xp
yj
ym
yp
zj
zm
zp
(5.7)
Còn các hệ số ai , bi ,..., d p là:
xj
yj
zj
ai = xm
xp
ym
yp
zm
zp
1
yj
zj
bi = − 1 ym
1 yp
zm
zp
xj
1 zj
ci = − xm 1 zm
xp 1 z p
xj
yj
1
d i = − xm
ym 1
xp
yp 1
(5.8)
(i , j , m , p )
Để cho thể tích V của tứ diện không mang giá trị âm thì các số hiệu nút i, j,m,p
phải xác định theo qui tắc vặn nút chai: khi i,j,m quay theo chiều i → j → m thì p
hướng đi lên.
Bằng cách tương tự ta được các biểu thức đối với các chuyển vị còn lại:
1
(ai + bi x + ci y + di z )vi − (a j + b j x + c j y + d j z )v j
6V
+ (am + bm x + cm y + d m z )vm − (a p + bp x + c p y + d p z )v p
(5.9)
1
(ai + bi x + ci y + di z ) wi − (a j + b j x + c j y + d j z ) w j
6V
+ (am + bm x + cm y + d m z ) wm − (a p + bp x + c p y + d p z ) wp
(5.10)
v=
w=
Thay (5.6),(5.9) và (5.10) vào (5.3) được vectơ chuyển vị
{ f } = [ ]{δ } = Ii
trong đó:
I j
I m
I p {δ }
I là ma trận đơn vị cấp 3x3, và
i = (ai + bi x + ci y + di z ) / 6V
j = (a j + b j x + c j y + d j z ) / 6V
m = (am + bm x + cm y + d m z ) / 6V
(5.11)
p = (a p + bp x + c p y + d p z ) / 6V
Dễ dàng thấy rằng, do hàm chuyển vị là tuyến tính, các mặt biên chung của hai
phần tử lân cận vẫn liền nhau, tức là đảm bảo điều kiện liên tục.
5.3. Biến dạng và ứng suất. Ma trận đàn hồi
Theo lý thuyết biến dạng tuyến tính ta có vectơ biến dạng gồm 6
thành phần:
∂u
∂x
∂v
ε x ∂y
ε
y ∂w
ε ∂z
{ε } = z = ∂u ∂v
γ xy +
γ yz ∂y ∂x
γ zx ∂v + ∂w
∂z ∂y
∂w + ∂u
∂x ∂z
(5.12)
Thay (5.6), (5.9), (5.10) vào biểu thức trên ta có
PPPTHH 117
{ε } = [ B ]{δ } = Bi
− Bj
Bm − B p {δ }
(5.13)
trong đó:
∂ i
∂x
0
0
[ Bi ] = ∂
i
∂y
0
∂ i
∂z
0
∂ i
∂y
0
∂ i
∂x
∂ i
∂z
0
0
0
bi
0
∂ i
1 0
∂z
=
6V ci
0
0
∂ i
di
∂y
∂ i
∂x
0
ci
0
bi
di
0
0
0
di
0
ci
bi
( i , j , m, p )
(5.14)
Ta thấy các thành phần trong [ B ] đều là hằng số, nên biến dạng trong phần tử là
hằng số.
Ứng suất trong phần tử cũng có thể biểu diễn qua chuyển vị nút:
{σ } = [ S ]{δ }
trong đó:
[ S ] = Si
−Sj
S m − S p
với
bi
Ab
1i
A1bi
E (1 −ν )
[ Si ] =
6(1 + ν )(1 − 2ν )V A2ci
0
A2 d i
A1 =
ν
1 −ν
,
A2 =
A1ci
ci
A1ci
A2bi
A2 d i
0
A1d i
A1d i
di
0
A2ci
A2bi
1 − 2ν
2(1 −ν )
Rõ ràng là ứng suất trong mỗi phần tử cũng là hằng số.
(i , j , m , p )
(5.15)
(5.16)
Ma trận đàn hồi đối với vật liệu đẳng hướng có dạng:
1
E (1 −ν )
[ D] =
(1 + ν )(1 − 2ν )
ν
ν
(1 −ν )
(1 −ν )
1
0
0
0
0
0
0
1 − 2ν
2(1 −ν )
0
ν
(1 −ν )
1
®èi
xøng
1 − 2ν
2(1 −ν )
0
0
0
0
1 − 2ν
2(1 −ν )
0
(5.17)
5.4. Ma trận độ cứng
Trước đây ta đã thành lập được biểu thức tổng quát của ma trận độ cứng của
phần tử:
[ k ] = ∫ [ B ] [ D ][ B ] dV
T
V
Vì các thành phần trong [ B ] là hằng số nên ta có
[ k ] = [ B ] [ D ][ B ]V
T
(5.18)
Cũng có thể viết
kii
−k
[ k ] = k ji
mi
−
k
pi
− kij
kim
k jj
− kmj
k pj
− k jm
kmm
− k pm
− kip
k jp
− kmp
k pp
(5.19)
Các thành phần krs (r = i, j , m, p ; r = i, j , m, p ) trong ma trận trên được tính
như sau:
Ab
Ab
brbs + A2 (cr cs + dr ds )
1 r cs + A2cr bs
1 r ds + A2dr bs
E(1−ν )
Ac
cr cs + A2 (brbs + dr ds )
Ac
[ krs ] =
1 r bs + A2br cs
1 r ds + A2dr cs
36(1+ν )(1− 2ν )V
Ad
b
+
A
b
d
Ad
c
+
A
c
d
d
d
+
A
(
b
b
+
c
c
)
1 r s
2 r s
1 r s
2 r s
r s
2 r s
r s
(5.20)
5.5. Dời tải trọng về nút
PPPTHH 119
Việc dời các tải trọng tác dụng trong phần tử về bốn nút cũng vẫn dựa trên nguyên
lý công ảo, tuy nhiên việc tính toán tương đối phức tạp, ở đây chỉ giới thiệu kết quả đối
với một số trường hợp.
* Lực thể tích có mật độ là p:
X
{ p} = Y
Z
Các lực nút tương đương là
{P}
e
= ∫[ ]
T
{ p} dV
(5.21)
V
Các thành phần của nó là
Pi = ∫ i { p} dV
(5.22)
V
Thí dụ, trọng lượng bản thân của phần tử là γ V , trong đó γ là trọng lượng riêng,
được phân đều lên 4 nút:
T
1
1
1
1
0 0
0 0
0 0
(5.23)
{P} = −γ V 0 0
4
4
4
4
* Tải trọng phân bố bậc nhất lên một mặt bên của tứ diện và vuông góc với mặt
đó, thí dụ áp lực thủy tĩnh tác dụng lên mặt ijm.
e
Gọi mật độ phân bố tại các nút i, j, m là qi , q j , qm thì lực nút tương đương theo
phương vuông góc với mặt ijm là
1
1
1
Pi = qi + q j + qm ∆ijm
6
2
2
1
1
1
Pj = q j + qi + qm ∆ ijm
6
2
2
(5.24)
1
1
1
Pm = qm + qi + q j ∆ ijm
6
2
2
trong đó:
∆ ijm là diện tích của mặt ijm.
5.6. Ứng suất nhiệt
Cũng giống như trong bài toán phẳng, ở đây ta vẫn sử dụng quan hệ vật lý
{σ } = [ D ] ({ε } − {ε o })
trong đó:
α T
α T
α T
{ε o } = = α T
0
0
0
1
1
1
0
0
0
(5.25)
Từ đó ta có
{σ } = [ D ][ B ]{δ } − [ D ]{ε o }
= [ S ]{δ } − [ D ]{ε o }
(5.26)
Ma trận [ S ] tính theo công thức (5.15). Lực nút tương đương
{P} = [ k ]{δ } − [ B ] [ D ]{ε o }V
e
T
(5.27)
Số hạng thứ hai ở vế phải
{R} = [ B ] [ D ]{ε o }V
e
T
biểu thị lực nút tương đương do ảnh hưởng của thay đổi nhiệt. Thay
[ B ] ở (5.14), [ D ] ở (5.17) và {ε o } ở (5.25) vào biểu thức này ta được
{R}
e
trong đó:
Xi
Y
i
Zi
X j
Yj
Eα (Ti + T j + Tm + Tp )
Eα T
Zj
= =
=
24(1 − 2ν )
X m 6(1 − 2ν )
Ym
Zm
X p
Yp
Z
p
bi
c
i
di
−b j
−c j
−d j
bm
cm
dm
−b p
−c p
−d
p
(5.28)
Ti , T j , Tm , Tp là thay đổi nhiệt độ ở các nút, có thể lấy giá trị trung
bình của chúng T thay cho độ thay đổi nhiệt độ của phần tử.
5.7. Về cách chia phân tử
PPPTHH 121
a)
b)
Hình 5.3
Giả sử ta xét một tường chắn như trên hình 5.3a. Nếu dùng các
mặt cắt song song với các mặt phẳng tọa độ để chia kết cấu thành các
phần tử rời rạc, ta sẽ được hai loại phần tử: phần tử lăng trụ tứ giác
gồm 8 điểm nút và loại phần tử lăng trụ tam giác gồm 6 điểm nút.
Lại có thể dùng một mặt cắt chéo để chia một phần tử lăng trụ tứ giác thành hai
phần tử lăng trụ tam giác. Có thể đưa hai loại phần tử này về loại phần tử tứ diện bằng
cách dùng một số mặt cắt chéo. Thí dụ, lăng trụ tam giác trên hình 5.3b có thể phân
thành 3 hình tứ diện: ABCD, EBFD và FBCD. Như vậy ta có thể sử dụng tất cả các
công thức đã trình bày trên đây.
Tuy nhiên cách chia vật thể khối thành các phần tử tứ diện như trên đây thường
làm cho mạng lưới dễ bị rối, và việc đánh số nút cũng dễ bị lầm lẫn. Trong nhiều trường
hợp, sẽ thuận tiện hơn. Nếu ta chia vật thể thành các khối có dạng hộp (hình 5.4) có 6
mặt với 8 điểm nút.
Hình 5.4
Phần tử loại này có tất cả 24 bậc tự do, tức là cũng có thể lấy 24 thông số chưa
biết trong các biểu thức biểu diễn qui luật gần đúng của chuyển vị. Nhờ đó, trạng thái
ứng suất và biến dạng trong phần tử được miêu tả chính xác hơn. Chẳng hạn, giả thiết
các thành phần chuyển vị được biểu diễn gần đúng bằng đa thức:
u = α1 + α 2 x + α 3 y + α 4 z + α 5 xy + α 6 yz + α 7 zx + α 8 x 2
v = α 9 + α10 x + α11 y + α12 z + α13 xy + α14 yz + α15 zx + α16 y 2
w = α17 + α18 x + α19 y + α 20 z + α 21 xy + α 22 yz + α 23 zx + α 24 z 2
(5.29)
Cách thiết lập ma trận độ cứng cũng tương tự như đối với phần tử tứ diện, ở đây
không trình bày lại.
Ta xem một thí dụ đơn giản về bài toán không gian vẽ trên hình 5.5. Đó là bài
toán Boussinesq quen thuộc đối với nửa không gian đàn hồi chịu lực tập trung theo
phương thẳng đứng. Nửa không gian được thay bằng hình hộp 6 mặt. Do tính đối xứng
trên hình vẽ chỉ thể hiện một góc phần tư. Vùng ở gần điểm đặt lực tập trung, các phần
tử khối lấy kích thước nhỏ hơn.
Điều kiện biên lấy như sau:
u = v = w = 0 trên mặt ABCD tức là ở vùng đủ xa điểm đặt lực;
u = 0 trên mặt AEHD và v = 0 trên mặt AEFB vì lý do đối xứng.
Các biên khác tự do.
Hình 5.5
5.8. Khái niệm về phần tử đẳng tham số
Ta bắt đầu từ một phần tử tứ giác phẳng có hình dạng bất kỳ để tìm hiểu một số
khái niệm cơ bản về phần tử đẳng tham số (nhưng mục đích là để dùng vào phần tử
đẳng tham số trong bài toán không gian).
Như trên đã nói, đối với bài toán phẳng, loại phần tử có 3 điểm nút được dùng
nhiều nhất, sau đó là phần tử chữ nhật có 4 điểm nút. Đối với phần tử chữ nhật do mô
hình chuyển vị là hàm bậc hai của tọa độ nên ứng suất trong phần tử không phải là hằng
số mà là thay đổi tuyến tính, cho nên nó có thể phản ánh tốt hơn tình hình phân bố ứng
suất so với phần tử tam giác. Nhưng phần tử chữ nhật không thích hợp với vật thể có
biên cong hoặc các biên không trực giao cũng như kích thước phần tử không dễ tùy ý
thay đổi. Nếu bây giờ ta sử dụng một loại phần tử tứ giác bất kỳ như hình 5.6 mà vẫn
dùng mô hình chuyển vị của phần tử chữ nhật thì ở trên biên chung của hai phần tử lân
cận chuyển vị sẽ không phải là thay đổi tuyến tính nữa, do đó tính liên tục của chuyển vị
ở biên chung không được bảo đảm. Sử dụng phương pháp biến đổi tọa độ ta có thể giải
quyết được mâu thuẫn này.
PPPTHH 123
Hình 5.6
Hình 5.7
Trên hình 5.6 là một phần tử tứ giác bất kỳ. Ta dùng hai họ đường thẳng phân đều
bốn cạnh. Lấy trung tâm của hai họ này ( ξ = η = 0 ) làm gốc, vẽ các trục ξ và η theo
các phương tăng của ξ và η , đồng thời lấy giá trị ở trên bốn biên là ±1 , ta sẽ được một
hệ tọa độ mới. Hệ tọa độ này là hệ tọa độ cục bộ, chỉ dùng trong phạm vi phần tử.
Sau đây ta xây dựng mô hình chuyển vị cũng như thành lập công thức biến đổi tọa
độ giữa tọa độ cục bộ và tọa độ tổng quát x,y. Để đơn giản ta xem xét một phần tử hình
vuông bốn nút (hình 5.7).
Như ta đã biết ở Chương III, mô hình chuyển vị của phần tử này là:
u = 1u1 + 2u2 + 3u3 + 4u4
v = 1v1 + 2 v2 + 3v3 + 4v4
(5.30)
trong đó các hàm dạng:
1
(1 − ξ )(1 − η )
4
1
2 = (1 + ξ )(1 − η )
4
1
3 = (1 − ξ )(1 + η )
4
1
4 = (1 + ξ )(1 + η )
4
1 =
(5.31)
Cũng có thể đưa vào biến mới:
ξo = ξi ξ ,
η o = ηi η
(i = 1, 2, 3, 4)
Thay vào (5.31) và viết chung vào một biểu thức:
i =
1
(1 + ξo )(1 + ηo )
4
(i = 1, 2, 3, 4)
(5.32)
Đem mô hình chuyển vị (5.30) và hàm dạng (5.31) hoặc (5.32) dùng vào phần tử
tứ giác ở hình 5.7 và quan niệm ξ và η là tọa độ cục bộ của phần tử đó thì có thể thấy
rõ, mô hình chuyển vị (5.30) sẽ cho các chuyển vị tại bốn nút, và trên các cạnh của phần
tử chuyển vị sẽ thay đổi tuyến tính, do đó bảo đảm được tính liên tục của chuyển vị. Do
đó (5.30) là mô hình chuyển vị cần tìm. Đồng thời theo (5.30) công thức biến đổi tọa độ
là:
x = 1 x1 + 2 x2 + 3 x3 + 4 x4
(5.33)
y = 1 y1 + 2 y2 + 3 y3 + 4 y4
nghĩa là công thức này sẽ cho các tọa độ tổng thể tại bốn nút, và trên bốn cạnh của phần
tử thì một tọa độ cục bộ bằng ±1 còn tọa độ cục bộ còn lại là thay đổi tuyến tính. Do đó
(5.33) là công thức biến đổi tọa độ cần tìm.
Như vậy có thể hiểu rằng phần tử hình vuông ở hình 5.7 là phần tử cơ bản, còn
phần tử tứ giác bất kỳ ở hình 5.6 là phần tử thực tế do phần tử cơ bản qua biến đổi mà
thành. Do đó trong biểu thức mô hình chuyển vị và trong biểu thức biến đổi tọa độ ta
dùng chung một hàm dạng, vì lẽ đó phần tử thực tế này được gọi là phần tử đẳng tham
số hoặc phần tử đồng tham số.
Phương pháp biến đổi trên đây có thể suy rộng đối với những phần tử có nhiều nút
hơn, thí dụ đối với phần tử thực tế ở hình 5.8 có tám nút và bốn cạnh cong. Để thiết lập
mô hình chuyển vị có thể lấy phần tử cơ bản là phần tử hình vuông tám nút như hình
5.9.
Hình 5.8
Hình 5.9
Mô hình chuyển vị của phần tử cơ bản là:
8
u = ∑ i ui ,
i =1
8
v = ∑ i vi
(5.34)
i =1
trong đó các hàm dạng là:
PPPTHH 125
1
1
(1 − ξ )(1 − η )(−ξ − η − 1) , 2 = (1 + ξ )(1 − η )(ξ − η − 1)
4
4
1
1
3 = (1 − ξ )(1 + η )(−ξ + η − 1) , 4 = (1 + ξ )(1 + η )(ξ + η − 1)
4
4
1
1
5 = (1 − ξ 2 )(1 − η ) ,
6 = (1 − ξ 2 )(1 + η )
2
2
1
1
7 = (1 − η 2 )(1 − ξ ) ,
8 = (1 − η 2 )(1 + ξ )
2
2
(5.35)
1
(1 + ξo )(1 + ηo )(ξ o + ηo − 1)
4
1
i = (1 − ξ 2 )(1 + ηo ) (i = 5, 6)
2
1
i = (1 − η 2 )(1 + ξ o ) (i = 7, 8)
2
(5.36)
1 =
hoặc
i =
(i = 1, 2, 3, 4)
Mô hình chuyển vị (5.34) sẽ cho các chuyển vị tại 8 nút. Đồng thời ở trên mỗi
cạnh của hình vuông thì có tọa độ bằng ±1 , do đó đều là hàm bậc hai của tọa độ còn lại.
Vì vậy chỉ cần 3 nút trên biên chung của hai phần tử lân cận có cùng chuyển vị thì tính
liên tục của chuyển vị được đảm bảo.
Đem mô hình (5.34) và các hàm dạng (5.35) hoặc (5.36) đưa vào phần tử thực tế
ở hình (5.9), có thể thấy đó chính là mô hình chuyển vị cần tìm. Đồng thời công thức
biến đổi tọa độ là:
8
x = ∑ i xi ,
i =1
8
y = ∑ i yi
(5.37)
i =1
tức là từ đó xác định được tọa độ ở 8 nút. Trên một biên nào đó luôn luôn có một tọa độ
cục bộ bằng ±1 , đồng thời x và y đều là biểu thức bậc hai của tọa độ cục bộ còn lại, do
đó đường cong biên là đường cong bậc hai (khi chia phần tử các tọa độ tổng thể của 8
nút phải xác định chính xác rồi từ đó vẽ các đường cong).
Để nâng cao độ chính xác có thể chọn phần tử cơ bản có nhiều nút hơn, sử dụng
phần tử bậc cao hơn dẫn đến phần tử thực tế cũng có bậc cao hơn. Đương nhiên lúc đó
khối lượng tính toán sẽ tăng nhiều.
Bây giờ ta xét đến phần tử không gian.
Cách tính toán đối với phần tử không gian cũng tương tự như đối với phần tử
phẳng. Thí dụ đối với phần tử cơ bản hình lập phương có 8 nút (hình 5.10) có thể căn cứ
vào (5.30) và (5.33) có thể chọn mô hình chuyển vị và biểu thức biến đổi tọa độ như
sau:
8
u = ∑ i ui ,
i =1
8
x = ∑ i xi ,
i =1
8
v = ∑ i vi ,
i =1
8
y = ∑ i yi ,
i =1
8
w = ∑ i wi
i =1
8
z = ∑ i zi
i =1
trong đó:
1
1
i = (1 + ξiξ )(1 + ηiη )(1 + ζ iζ ) = (1 + ξo )(1 + ηo )(1 + ζ o )
8
8
Đối với phần tử cơ bản hình lập phương có 20 nút (hình 5.11) ta có
20
u = ∑ i ui ,
i =1
20
x = ∑ i xi ,
i =1
20
v = ∑ i vi ,
i =1
20
y = ∑ i yi ,
i =1
20
w = ∑ i wi
i =1
20
z = ∑ i zi
i =1
trong đó:
1
i = (1 + ξo )(1 + ηo )(1 + ζ o )(ξo + ηo + ζ o − 2) (i = 1, 2,..., 8)
8
1
i = (1 − ξ 2 )(1 + ηo )(1 + ζ o )
(i = 9,10,11,12)
4
1
i = (1 − η 2 )(1 + ζ o )(1 + ξ o )
(i = 13,14,15,16)
4
1
i = (1 − ζ 2 )(1 + ξo )(1 + ηo )
(i = 17,18,19, 20)
4
Sau khi biến đổi, từ phần tử hình lập phương ta được phần tử khối 6 mặt cong
(hình 5.11b)
b)
a)
Hình 5.10
Hình 5.11
PPPTHH 127
5.9. Tính toán phần tử không gian đẳng tham số
Như trên đã biết, đối với phần tử không gian đẳng tham số có n nút thì mô hình
chuyển vị và công thức biến đổi tọa độ có dạng:
n
u = ∑ i ui ,
i =1
n
x = ∑ i xi ,
i =1
n
v = ∑ i vi ,
i =1
n
y = ∑ i yi ,
i =1
n
w = ∑ i wi
(5.38)
i =1
n
z = ∑ i zi
(5.39)
i =1
Khi tính toán, vì phải dùng đến đạo hàm của hàm dạng i đối với tọa độ tổng thể,
nên trước hết ta cần thành lập công thức tính đạo hàm này.
Đạo hàm của hàm đa biến i là:
∂ i ∂ i ∂x ∂ i ∂y ∂ i ∂z
=
+
+
∂ξ
∂x ∂ξ ∂y ∂ξ ∂z ∂ξ
. . . . . . . . . . . . . . . .
do đó ta có
∂ i
∂ξ
∂ i
=
∂η
∂ i
∂ζ
∂x
∂ξ
∂x
∂η
∂x
∂ζ
∂y
∂ξ
∂y
∂η
∂y
∂ζ
∂z
∂ξ
∂z
∂η
∂z
∂ζ
∂ i
∂ i
∂x
∂x
∂ i
∂ i
= [J ]
∂y
∂y
∂ i
∂ i
∂z
∂z
(5.40)
từ đó được
∂ i
∂ i
∂ξ
∂x
∂ i
−1 ∂ i
= [J ]
∂y
∂η
∂ i
∂ i
∂z
∂ζ
trong đó
(5.41)
∂x
∂ξ
∂x
[J ] =
∂η
∂x
∂ζ
∂y
∂ξ
∂y
∂η
∂y
∂ζ
∂z
∂ξ
∂z
∂η
∂z
∂ζ
(5.42)
gọi là ma trận Jacobien. Để xác định ma trận này, đưa (5.39) vào (5.42) ta được
n ∂ i
∑ ∂ξ xi
i =1
n ∂
[ J ] = ∑ i xi
i =1 ∂η
n ∂ i
xi
∑
i =1 ∂ζ
∂1
∂ξ
∂
= 1
∂η
∂1
∂ζ
∂ 2
∂ξ
∂ 2
∂η
∂ 2
∂ζ
n
∂ i
∑ ∂ξ
n
yi
i =1
n
∂ i
∑ ∂η
∂ i
∑ ∂ξ
i =1
n
yi
i =1
∂ i
yi
∑
i =1 ∂ζ
n
∂ n
∂ξ
∂ n
L
∂η
∂ n
L
∂ζ
L
∂ i
∑ ∂η
i =1
n
∂ i
∑ ∂ζ
i =1
x1
x
2
M
xn
y1
y2
M
yn
zi
zi
zi
z1
z2
M
zn
(5.43)
Như vậy sau khi tìm được i (là hàm của ξ ,η , ζ ), tiến hành đạo hàm đối với tọa
độ cục bộ rồi thay vào (5.43) tìm được [ J ] , tính tiếp [ J ] và từ (5.41) tìm được các đạo
hàm của các hàm dạng đối với tọa độ tổng quát.
Trong quá trình tính toán còn cần dùng tới công thức tích phân các hàm của tọa độ
cục bộ đối với tọa độ tổng quát, do đó có thể sử dụng quan hệ đã biết trong hình học
vi phân:
−1
dV = dxdydz = J d ξ dη d ζ
trong đó:
(5.44)
J là định thức của ma trận Jacobien. Ở đây J lấy giá trị tuyệt đối.
Sau đây ta tiến hành tính toán đối với phần tử không gian đẳng tham số, tức là
thiết lập vectơ tải, ma trận ứng suất và ma trận độ cứng.
* Khi trong phần tử có lực tập trung đặt tại điểm bất kỳ
{P} = Px
Py Pz
T
(5.45)
cần phải dời về các nút, ta vẫn dùng công thức quen thuộc:
PPPTHH 129
{P}
= [ ]
{P}
= [ X 1 Y1 Z1
e
{P}
T
(5.46)
trong đó
e
và
[ ] = [ I1
X 2 Y2 Z 2 ... X n Yn Z n ]
T
I 2 ... I n ]
(5.47)
(5.48)
với I là ma trận đơn vị cấp 3.
* Trường hợp có lực phân bố thể tích
{ p} = [ X
Y Z]
T
(5.49)
có thể sử dụng quan hệ (5.52), và ta có
{P}
e
= ∫[ ]
T
V
{ p} dV = ∫−1 ∫−1 ∫−1[ ] { p} J
1
1
1
T
d ξ dη d ζ
(5.50)
* Trường hợp lực phân bố bề mặt
{ p} = X
Y Z
T
(5.51)
cũng có thể sử dụng công thức (5.46).
Thí dụ trên mặt ξ = 1 có lực bề mặt { p} thì
{P}
e
= ∫[ ]
T
T
{ p} dA = ∫−1 ∫−1[ ] ξ =1 { p} J ξ =1 dη d ζ
1
1
(5.52)
Các lực bề mặt khác cũng xử lý tương tự.
Bây giờ ta thiết lập ma trận ứng suất và ma trận độ cứng. Trước hết thay (5.38)
vào quan hệ
{ε } = [ B ]{δ } = [ B1
với
{δ } = [u1
B2 ... Bn ]{δ }
v1 w1 u2 v2 w2 ... un vn wn ]
(5.53)
T
(5.54)
∂ i
∂x
0
0
[ Bi ] = ∂
i
∂y
0
∂ i
∂z
0
0
∂ i
∂z
0
∂ i
∂y
∂ i
∂x
0
∂ i
∂y
0
∂ i
∂x
∂ i
∂z
0
(5.55)
Sau đó tính ứng suất trong phần tử theo quan hệ sau đây:
{σ } = [ D ][ B ]{δ } = [ S ]{δ }
(5.56)
trong đó
[ S ] = [ S1
với
S 2 ... S n ]
[ Si ] = [ D ][ Bi ]
(i = 1, 2, ..., n )
Lực nút trên phần tử là:
{P}
e
= [U1 V1 W1 U 2 V2 W2 ... U n Vn Wn
]
T
Như ta đã biết
{P} = [ k ]{δ }
e
trong đó ma trận độ cứng
[ k ] = ∫ [ B ] [ D ][ B ] dV
T
V
=∫
(5.57)
1
∫ ∫ [ B ] [ D ][ B ] J dξ dη dζ
1
1
T
− 1 −1 −1
[ B ] tính theo tọa độ tổng thể, còn [ J ] viết theo tọa độ cục bộ.
Ma trận độ cứng của phần tử:
k11
k
[ k ] = M21
kn1
k12 L k1n
k21 L k21
M
M
kn1 L knn
(5.58)
PPPTHH 131
trong đó các ma trận con
[ krs ] = ∫−1 ∫−1 ∫−1[ Br ] [ D ][ Bs ] J
1
1
1
T
d ξ dη dζ
(5.59)
Qua trình bày trên đây ta thấy phần tử đẳng tham số có ưu điểm là nó có thể phản
ánh tương đối tốt sự thay đổi chuyển vị và ứng suất trong vật thể đàn hồi so với phần tử
thường khi số nút như nhau. Trên hình 5.12 là sơ đồ tính một công xon chịu lực tập
trung theo các phương án: sử dụng các phần tử tam giác (hình a), phần tử đẳng tham số
4 nút (hình b), phần tử đẳng tham số 8 nút (hình c), phần tử đẳng tham số 12 nút (hình
d). Theo [I], cùng số nút như nhau nhưng kết quả tính chuyển vị của điểm đặt tải trọng
theo sơ đồ (a) bằng 53% giá trị tính theo lời giải chính xác, trong khi đó theo sơ đồ (c)
là 65%, theo sơ đồ (d) thì kết quả trùng với lời giải chính xác.
a)
b)
c)
d)
Hình 5.12
5.10. Tích phân Gauss
Khi thiết lập ma trận độ cứng và vectơ tải của phần tử đẳng tham số ta cần tính
tích phân sau đây:
1
1
∫ ∫ ∫
1
−1 −1 −1
f (ξ ,η , ζ ) dξ dη dζ
(5.60)
Hàm f (ξ ,η , ζ ) dưới dấu tích phân nói chung rất phức tạp, ngay cả khi viết được
dưới dạng tường minh, cho nên để được kết quả người ta thường dùng phép tích phân
số. Nội dung của phép tích phân này là người ta chọn một số điểm trong phần tử, gọi là
điểm tích phân, rồi tìm giá trị của hàm f tại các điểm đó, sau đó căn cứ vào các giá trị số
đó để tìm giá trị số của biểu thức tích phân. Với phép tích phân Gauss, có thể dùng
tương đối ít số điểm tích phân mà vẫn cho độ chính xác tương đối cao.
Trước hết giới thiệu công thức tích phân Gauss đối với hàm một biến.
Giả sử có hàm f (ξ ) là một đa thức thì
∫
1
−1
trong đó:
n
f (ξ ) d ξ = ∑ H i f (ξi )
(5.61)
i =1
f (ξi ) là giá trị hàm f tại điểm tích phân ξi , H i là trọng số, n là số
lượng điểm tích phân. Bảng 5.1 giới thiệu các giá trị của ξi và H i với
n = 2 đến n = 7 .
Công thức tích phân Gauss đối với hàm hai biến có dạng như sau:
1
∫ ∫
1
−1 −1
n
n
f (ξ ,η ) d ξ dη = ∑∑ H i H j f (ξi ,η j )
(5.62)
j =1 i =1
Số điểm tích phân theo một phương là n, số điểm tích phân theo hai phương là n 2 .
Công thức tích phân Gauss đối với hàm ba biến (ta sẽ dùng khi tính ma trận độ
cứng của phần tử đẳng tham số) có dạng:
1
1
∫ ∫ ∫
1
−1 −1 −1
n
n
n
f (ξ ,η , ζ ) d ξ dη dζ = ∑∑∑ H i H j H m f (ξi ,η j , ζ m )
m =1 j =1 i =1
(5.63)
Số điểm tích phân theo cả 3 phương là n3 .
Tọa độ điểm tích phân và trọng số của tích phân Gauss
Bảng 5.1
∫
1
−1
f (ξ ) d ξ =
n
∑H
i
f (ξi )
i =1
________________________________________________________
±ξi
Hi
n=2
0,5773502692
1,0000000000
n=3
0,7745966692
0,5555555556
0,0000000000
0,8888888889
n=4
0,8611363116
0,3478548451
0,3399810436
0,6521451549
PPPTHH 133
- Xem thêm -
.PDF 
88 
234 
70 
