Bộ giáo dục và đào tạo bộ quốc phòng
Học viện kỹ thuật quân sự
.......................................... Trần Bá Thnh đề ti luận văn
Bài toán trị riêng trong phơng pháp phần tử hữu hạn
giải cho hệ dầm liên tục
Mã số: 60 58 50 luận văn thạc sỹ kỹ thuật Ngời hớng dẫn khoa học:
PGS.TS. Nguyễn Quốc Bảo
Hà Nội, Năm 2008 Bộ giáo dục và đào tạo bộ quốc phòng
Học viện kỹ thuật quân sự
.......................................... luận văn thạc sỹ kỹ thuật
LỜI CAM ĐOAN
Tôi xin cam đoan luận văn này là công trình nghiên cứu của riêng tôi.
Các số liệu, kết quả nêu trong luận văn là trung thực và chưa từng được
ai công bố trong bất kỳ công trình nghiên cứu nào khác./.
Tác giả luận văn
Trần Bá Thành
MỤC LỤC
PHẦN MỞ ĐẦU................................................................................................... 1
CHƯƠNG 1: TỔNG QUAN................................................................................. 3
CHƯƠNG 2: XÂY DỰNG PHƯƠNG TRÌNH TRỊ RIÊNG HỆ DẦM LIÊN
TỤC THEO PHƯƠNG PHÁP PTHH ........................................... 4
2.1.
Xây dựng tính chất phần tử....................................................................... 4
2.1.1
Ma trận độ cứng phần tử................................................................... 4
2.1.2
Những tính chất chủ yếu của trị riêng, vectơ riêng ................................ 32
3.2.1
Tính chất của véctơ riêng................................................................ 32
3.2.2
Đa thức đặc trưng của bài toán trị riêng ......................................... 34
3.2.3
Trượt trị riêng.................................................................................. 35
3.3.
Chuyển từ bài toán trị riêng tổng quát sang bài toán trị riêng chuẩn ..... 36
3.3.1
Sự cần thiết ..................................................................................... 36
3.3.2
Các bước chuyển từ bài toán tổng quát sang bài toán chuẩn.......... 36
3.4.
Các kỹ thuật giải áp dụng trong giải bài toán trị riêng ........................... 38
Lặp xuôi vectơ ................................................................................ 49
3.6.3
Trượt trị riêng trong lặp vectơ ........................................................ 52
3.6.4
Lặp thương số Rayleigh.................................................................. 52
3.6.5
Tốc độ hội tụ trong phương pháp lặp ............................................. 54
3.7.
Phương pháp biến đổi ma trận hay chéo hóa ma trận............................. 59
3.7.1
Phương pháp xoay Jacobi dùng cho bài toán chuẩn....................... 61
3.7.2
Phương pháp Jacobi dùng cho bài toán tổng quát.......................... 63
3.7.3
Phương pháp lặp ngược Householder – QR................................... 70
Một số chú ý khi chọn vectơ lặp ban đầu ....................................... 79
3.9.4
Sự hội tụ.......................................................................................... 80
CHƯƠNG 4: THỬ NGHIỆM LẬP TRÌNH TRÊN MATLAB ......................... 82
4.1.
Tổ chức số liệu trong chương trình......................................................... 82
4.1.1
Số liệu vào SLV.............................................................................. 83
4.1.2
Số liệu trung gian SLTG................................................................. 83
4.1.3
Số liệu kết quả tính SLR................................................................. 84
4.2.
Tổ chức chương trình và một số hàm cơ bản ......................................... 84
4.2.1
iển hình là các ứng xử liên quan đến tác động động đất,
tác động gió, va xô tàu thuyền...Ví dụ như cầu đường sắt Kevda (Nga)
bị phá hủy năm 1875, cầu Menkhienxtein (Thụy Sỹ) bị phá hủy năm
1891, cầu dàn Quebec (Canada) bị phá hủy năm 1907, cầu dàn Mojur
(Nga) bị phá hủy năm 1925.
Điều 4.7.1.5 của tiêu chuẩn thiết kế cầu 22 TCN 272-05 có ghi: “trừ
khi được chỉ rõ, phải sử dụng các dạng và tần số c
ủa dao động riêng
không giảm rung để đáp ứng yêu cầu thiết kế về ứng xử động học đàn
hồi”. Như vậy, có thể thấy mọi tính toán liên quan đến ứng xử động
lực học đều có liên quan đến tần số và dạng dao động riêng.
Nhằm tìm hiểu và đóng góp một phần vào lĩnh vực này, học viên đã
chọn hướng nghiên cứu là cách tính tần số dao độ
ng riêng của kết cấu,
đặc biệt là các kết cấu có số lượng phần tử lớn và lấy dầm liên tục là
một ví dụ để khảo sát.
Mục đích, đối tượng và phạm vi nghiên cứu
Tìm hiểu phương pháp được sử dụng phổ biến và hiệu quả trong việc
tính toán bài toán trị riêng, vectơ riêng của các hệ kết cấu lớn (trong
trường hợp này là kết cấu dầm liên tụ
c), giúp cho người dùng cũng
Trang 2
như các nhà nghiên cứu có được một công cụ dễ hiểu, trực quan khi
cần phân tích dao động của kết cấu là phương pháp phần tử hữu hạn
(PTHH).
Đi sâu vào việc nghiên cứu thuật toán của một trong phương pháp tính
tần số dao động riêng tương đối hiện đại và phổ biến hiện nay khi tính
tần số và dạng dao động riêng của hệ kết cấu có số lượng phần tử lớ
n
là phương pháp lặp không gian con. Trên cơ sở đó phát triển và giải
như SAP, STAD....Tuy nhiên, mộ
t chương trình chuyên sâu vào việc
phân tích tính toán dao động và tần số dao động riêng để từ đó có thể
nghiên cứu các biện pháp giảm thiểu hoặc loại trừ các tác động gây
hại của Việt Nam hiện vẫn đang trong quá trình hình thành.
Đề tài này đề cập đến những cơ sở lý thuyết dùng trong phương pháp
PTHH để tính toán bài toán trị riêng của các kết cấu có số lượng phần
tử lớn và ứng dụng nó trong ngôn ngữ lập trình Matlab
để tính toán
các trị riêng cho dầm liên tục.
Trang 4
CHƯƠNG 2: XÂY DỰNG PHƯƠNG TRÌNH TRỊ RIÊNG HỆ DẦM
LIÊN TỤC THEO PHƯƠNG PHÁP PTHH
2.1. Xây dựng tính chất phần tử
Tính chất phần tử của dầm liên tục trong bài toán trị riêng bao gồm:
− Ma trận độ cứng phần tử [EK] hay [K]
e
− Ma trận khối lượng phần tử [EM] hay [M]
e
− Ma trận chuyển toạ độ [ET] hay [T]
e
Ngoài ra, nhằm phục vụ cho việc lập trình sau này, ở phần này sẽ trình
bày thêm cách xây dựng thuật toán xếp các ma trận phần tử vào ma
trận tổng quát.
2.1.1 Ma trận độ cứng phần tử
2.1.1.1 Công thức tổng quát xây dựng ma trận độ cứng phần tử
Ma trận [EK] hay [K] chứa các thành phần biểu thị độ cứng được gọi
SXX
Sb
dudVudV
T
S
T
V
T
V
∫∫∫
+=
σε
Thay: {u} = [N]{d}; {e} = [B]{d}; {s} = [C][B]{d} vào công thức
trên ta có phương trình cân bằng phần tử
[K]{d}= {Q}
Với:
[]
∫
=
V
T
CBBK
- là ma trận độ cứng phần tử (2.1)
Công thức (2.1) cũng là công thức tổng quát tính ma trận độ cứng
phần tử.
Ma trận [K] là ma trận thưa có 2 tính chất quan trọng là tính đối xứng
và tính chất băng.
Với các bài toán động lực học khi xấp xỉ gia tốc chuyển vị bằng cùng
hàm dáng [N] với chính chuyển vị người ta thu được công thức
nguyên là thẳ
ng và vuông góc với trục thanh).
y
z
x
O
l
Trang 7
Theo giả thiết trên, về chuyển vị các điểm trên trục thanh chỉ chuyển
vị theo trục y, vectơ chuyển vị {u} = v(x). Nhưng càng xa trục thanh
theo giả thiết tiết diện phẳng xuất hiện chuyển vị dọc trục u. Để tính
chuyển vị này cần đưa thêm thành phần chuyển vị góc xoay
dx
dv
x
z
=)(
θ
Cũng theo giả thiết trên trong thanh chỉ xuất hiện một thành phần biến
dạng
dx
du
x
=
ε
và một thành phần vectơ ứng suất
xx
E
εσ
z
=
)(
θ
nên chỉ cần chọn hàm xấp xỉ cho v. Với thanh dầm ta sử
dụng phương trình bậc 3 để xấp xỉ.
v = a
0
+ a
1
x + a
2
x
2
+ a
3
x
3
= α
1
L
1
3
+ α
2
L
2
3
+ α
3
1
LLLLLLLL
lx
L
L
v
x
L
L
v
dx
dv
z
ααααααθ
+++−−−=
∂
∂
∂
∂
+
∂
∂
∂
∂
==
(2.4)
Sử dụng các điều kiện biên:
Tại nút 1: v = v
1
θ
2
θ
z
= θ
z2
L
1
= 0 L
2
= 1
Trang 8
ta có: v
2
= α
2
)3(
1
242
ααθ
+−=
l
z
Suy ra: α
4
= -(lθ
z1
+ 3v
1
2
212
2
22
2
11
2
1
)23()23(
z
v
z
v
LlLLLLlLLLv
θ
θ
(2.5)
So sánh với biểu thức tổng quát {u} = [N]{d} ta có ma trận hàm dáng
cho thanh uốn:
[]
[ ]
2
212
2
22
2
11
2
1
ε
So sánh với biểu thức tổng quát {ε} = [B]{d} ta có ma trận biến dạng
cho thanh chịu uốn:
[]
)42(126)42(126
1
212121
2
LLlLLLlL
l
yB −−−−−−=
(2.7)
Ma trận vật liệu chỉ chứa một hằng số là moduyn đàn hồi [C]= E
Tính ma trận độ cứng
[] [][][]
∫
=
V
T
dxBCBk
Áp dụng công thức tính tích phân một chiều trong hệ toạ độ tự nhiên.
)!1(
!!
21
++
=
∫
qp
=
l
EI
DX
l
EI
l
EI
l
EI
l
EI
l
EI
l
EI
l
EI
l
EI
l
EI
k
4
612
264
612612
33
2
2323
⎢
⎢
⎣
⎡
−
−
−
−
==
l
EI
DX
l
EI
l
EI
l
EF
l
EI
l
EI
l
EI
l
EI
l
EI
l
EI
T
V
∫∫∫∫
===
ρρρ
Trong đó:
∫
=
S
m
dS
ρρ
- là khối lượng trên một đơn vị chiều dài thanh
Trang 10
Như vậy ta còn phải tính thành phần thứ 2 trong biểu thức tổng quát là
tích phân
[][]
dlNN
T
l
∫
.
Tích phân này phụ thuộc vào hàm dáng [N]
2.1.2.2 Ma trận khối lượng của thanh uốn thuần tuý
Ma trận hàm dáng của thanh dầm có dạng:
[]
[ ]
2
1
)23()23(
)23(
)23(
LlLLLLlLLL
LlL
LL
LlL
LL
NN
T
−−−
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎣
⎡
−
−
−
=
2
3
1
2
1
2
2
3
1
2
4
21
2
2
4
22
3
2
2
121
2
2
2
1
3
2
3
1
2
2
)23()23)(23()23()23(
LLlLLlLLLlLLlL
LLlLLLLLlLLLLL
LLlLLlLLLlLLlL
LLlLLLLLLLlLL
L
NN
T
Áp dụng công thức:
)!1(
!!
21
++
=
∫
qp
lqp
dlLL
qp
(2.11)
Ta có:
[]
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
tính chính trung tâm của mặt cắt ngang của thanh và chiều dương của
các trục x, y, z xác định theo quy tắc tam diện thuận.
Trong một kết cấu hệ thanh (dàn, khung ...) thường các phần tử thanh
có phương khác nhau, nên nói chung hệ tọa độ của từng phần tử
không giống nhau. Hệ toạ độ riêng đó đối với từng phần tử ta gọi là hệ
toạ độ phần tử hay hệ toạ độ địa ph
ương.
Khi tính toán kết cấu gồm nhiều phần tử, để thuận tiện khi thành lập
các phương trình cân bằng, người ta cần sử dụng một hệ tọa độ chung
cho toàn bộ kết cấu, thường gọi là hệ toạ độ kết cấu hoặc hệ toạ độ
chung.
Vì vậy trước khi bắt tay vào việc lập phương trình cân bằng ở các nút
cần phải biến đổi quan hệ
giữa các lực nút và chuyển vị nút trong hệ
tọa độ phần tử thành quan hệ giữa lực nút và chuyển vị nút trong hệ
tọa độ kết cấu, sự biến đổi đó được gọi là biến đổi tọa độ.
Trước hết, ta xét một hệ thanh phẳng như hình dưới. Các hệ tọa độ địa
phương là xyz, hệ tọa độ chung là x'y'z', trục z và z' hướng ra ngoài.
x'
y'
y
−=
100
0cossin
0sincos
3
ϕϕ
ϕϕ
T
Sắp xếp các bậc tự do phần tử theo trật tự
{dz}={u
1
v
1
θ
z1
u
2
v
2
θ
z2
}
Ta có ma trận chuyển cho phần tử thanh phẳng
⎥
⎦
⎤
⎢
⎣
⎡
=
100000
0cossin000
0sincos000
000100
0000cossin
0000sincos
ϕϕ
ϕϕ
ϕϕ
ϕϕ
hay
Nếu chỉ xét riêng thanh dầm hay hệ dầm liên tục, ta có ma trận chuyển
suy ra từ ma trận chuyển thanh phẳng trên bằng cách bỏ đi các thành
phần ứng với bậc tự do u
1
và u
2
, nghĩa là bỏ đi hàng 1 và 4, cột 1 và 4.
Nếu chọn trục x’ trùng với trục dầm x thì ϕ=0.
Ta có ma trận chuyển thanh dầm là ma trận đơn vị, ký hiệu T4.
Trang 13
⎥
⎥
⎥
⎥
⎦
⎤
⎢
ij
= k
ji
.
+ Tính chất băng có thể hiểu là các phần tử có nghĩa của [K] chỉ nằm
xung quanh đường chéo chính. Để xác định chính xác hơn khái niệm
băng của ma trận, ta xét ma trận [K] cấp n, ma trận có n hàng và n cột.
Xét cột i bất kì, đi từ phần tử trên đường chéo chính k
ii
dọc theo cột i
lên trên ta thấy có những phần tử khác 0 và có những phần tử bằng 0
cho tới khi đến một phần tử khác 0 cuối cùng, trên đó chỉ còn các
phần tử bằng 0. Ký hiệu chỉ số hàng của phần tử này là m
i
. Vì đi từ
đường chéo chính nên 1 ≤ m
i
≤ i.
Nếu m
i
luôn luôn bằng 1 với mọi cột
thì ma trận là ma trận đủ (ma trận đầy) và tính chất băng không còn
nữa. Khi có nhiều cột có m
i
>1, gọi h
i
= i - m
i
là chiều cao cột i, kí
, một chiều của mảng có
kích thước bằng hạng của ma trận, chiều thứ hai bằng b+1 với b là bề
rộng băng của ma trận [K]. Cột một của mảng lưu trữ các phần tử nằm
trên đường chéo chính, các cột tiếp theo lưu b giá trị cho mỗi cột. Khi
đó số phần từ cần lưu là n (b+1). Phương pháp này tiết kiệm bộ nhớ
hơn phương pháp đầu khi b+1 < (n+1)/2.
Điều kiện này dễ được thoả
mãn bằng cách đánh số nút hợp lý. Bề rộng nửa dải b càng bé càng tiết
kiệm bộ nhớ. Tuy nhiên cách lưu theo băng là chưa hoàn toàn tiết
Trang 15
kiệm, ngay cả khi bề rộng nửa dải hẹp đến mức tối ưu. Số giá trị
không vô nghĩa buộc phải lưu theo băng là Σ(b-h
i
) i=1÷n.
Ở cách lưu thứ ba, lưu theo chiều cao cột. Theo cách lưu này, người
lập trình sử dụng mảng động một chiều, ký hiệu {SK}, các phần tử
của ma trận độ cứng [K] được xếp lần lượt vào {SK} theo từng cột.
Trong mỗi cột bắt đầu xếp từ phần tử nằm trên đường chéo chính k
ii
lên trên cho đến hết chiều cao cột. Như vậy trật tự các phần tử của [K]
trong mảng một chiều {SK} sẽ là:
{SK} = {k
11
k
22
k
12
k
33
các giá trị thực.
Do đa số các h
i
đều
nhỏ hơn b nên tổng trên nhỏ hơn đáng kể so với
với n x (b+1) là số giá trị cần lưu theo sơ đồ băng.
Dưới đây là bảng so sánh bộ nhớ và thời gian giải cần thiết cho lời
giải trong phương pháp PTHH với ma trận độ cứng M thông thường
cho 2 dạng lưu dữ liệu
Trang 16
Với cách lưu mọi phần tử của ma trận
Hạng ma trận Bộ nhớ
Thời gian
Máy trạm/PC
Thời gian
Supercomputer
10
4
8MB 5giây 0.05giây
10
5
240MB 8phút 5giây
10
6
8GB 15giờ 8phút
Với cách lưu skyline
Hạng ma trận Bộ nhớ
Thời gian
I
II
III
IV
V
+ Bước 1: Xây dựng mảng JF ( mảng bậc tự do nút )
Mỗi nút tổng quát có 6 bậc tự do là 3 chuyển vị theo 3 trục toạ độ và 3
góc xoay quanh 3 trục [ u
i
v
i
w
i
θ
xi
θ
yi
θ
zi
]
Bậc tự do của một nút được chia làm 2 loại:
•
Bậc tự do tích cực: là chuyển vị nút cần tìm
•
Bậc tự do tiêu cực là bậc tự do có sẵn các giá trị xác định do điều kiện
biên hoặc không dùng đến do hạn chế của bài toán
Quản lý các bậc tự do nút là xác định xem một bậc tự do là tích cực
2 1 0 1 1 1 0
3 0 0 1 1 1 0
4 0 0 1 1 1 0
5 1 1 1 1 1 1
TsNut 6 0 1 1 1 1 0
+ Bước 2: Xây dựng mảng NDF (Đánh số bậc tự do nút)
Để đánh số bậc tự do nút, ta dùng mảng 2 chiều ký hiệu NDF(tổng số
nút,6), mảng này được xây dựng từ mảng JF bằng cách đổi các giá trị
1 (bậc tự do tiêu cực) về 0 còn các giá trị 0 (bậc tự do tích cực) thì
được đếm lần lượt từ trái sang phải, từ trên xuống dưới. Nhìn vào
mảng NDF ta nhận thấy, nếu giá tr
ị của phần tử mảng nào bằng 0 thì
ta không cần tìm chuyển vị tương ứng tại nút đó, ngược lại nếu giá trị
của phần tử mảng khác 0 thì ta cần phải tìm chuyển vị tương úng tại
nút đó. Giá trị lớn nhất của phần tử mảng NDF cũng chính là số
chuyển vị phải tìm hay số phương trình cần giải, ký hiệu là NEQ.
Copyright: Tài liệu đại học ©