Bài giảng Phương pháp phần tử hữu hạn trong mô phỏng và tính toán ô tô

0 TRƯỜNG ĐẠI HỌC SƯ PHẠM KỸ THUẬT HƯNG YÊN BÀI GIẢNG PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN TRONG MÔ PHỎNG VÀ TÍNH TOÁN Ô TÔ DÙNG CHO CHƯƠNG TRÌNH CAO HỌC Giáo viên: TS Nguyễn Thanh Quang Hưng Yên 2013 1 MỤC LỤC Trang Đề cương bài giảng 2 Phần 1 KHÁI NIỆM CHUNG VỀ PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN 1.1 Mô hình rời rạc hóa kết cấu 3 1.2 Hàm chuyển vị, Hàm dạng 6 1.2.1 Hàm chuyển vị 6 1.2.2 Hàm dạng 7 1.2.3 Lực nút 9 1.3 Phương trình cơ bản của phương pháp PTHH 10 1.3.1 Các quan hệ chuyển vị, biến dạng, ứng suất trong phần tử 10 1.3.2 Phương trình cơ bản của phương pháp phần tử hữu hạn 11 1.3.3 Ma trận độ cứng tổng thể 12 Phần 2 CÁC PHƯƠNG PHÁP GIẢI BÀI TOÁN PHẦN TỬ HỮU HẠN 2.1 Phương pháp sai phân PTHH 14 2.2 Phương pháp năng lượng và PP Niu tơn-Lagrăng 16 2.2.1 Cách 1: Xây dựng phiến hàm 16 2.2.2 Cách 2: Giải theo hệ PTVP cho trước 17 2.3 Phương pháp giải trên phần mềm 2.3.1 Giới thiệu PTHH trong MatLab 43 2.3.2 Giới thiệu PTHH trong Ansys 44 2 ĐỀ CƯƠNG BÀI GIẢNG MỞ ĐẦU Chương 1. Bổ túc về cơ học vật rắn và các phương pháp tính trong cơ học (TL1) 1.1. Cơ học vật rắn biến dạng 1.2. Các phương pháp giải bài toán cơ vật rắn biến dạng, bài toán đàn hồi 1.3. Các nguyên lý năng lượng (các nguyên lý biến phân) Chương 2. Cơ sở và các bước phân tích chung của phương pháp phần tử hữu hạn (TL2) 2.1. Khái niệm về phương pháp phần tử hữu hạn 2.2. Trình tự phân tích bài toán theo phương pháp phần tử hữu hạn 2.3. Hàm xấp xỉ  đa thức xấp xỉ  phép nội suy 2.4. Các phương trình cơ bản Chương 3. Tính toán hệ thanh (TL3) 3.1. Hệ thanh dàn 3.2 Khung phẳng 3.3 Khung không gian Chương 4. Bài toán phẳng của lý thuyết đàn hồi (TL4) 4.1 Các phương trình cơ bản của bài toán phẳng của lý thuyết đàn hồi 4.2 Bài toán phẳng với phần tử dạng tam giác 4.3 Phần tử chữ nhật Chương 5. Phần tử bậc cao và phần tử đẳng tham số (TL5) 5.1 Hệ tọa độ tự nhiên của các loại phần tử 5.2 Phần tử một chiều bậc cao 5.3 Phần tử tam giác bậc cao trong hệ tọa độ tự nhiên 5.4 Phần tử 3 chiều – Khối tứ diện 5.5 Phần tử 2 chiều dạng tứ giác 5.6 Phần tử đẳng tham số Chương 6. Tấm chịu uốn (TL6) 6.1 Các phương trình cơ bản của tấm chịu uốn 6.2 Phần tử tấm không tương thích dạng tam giác 6.3 Phần tử đẳng tham số dạng tứ giác bốn nút Chương 7. Phần tử 3 chiều (TL7) 7.1 Phần tử tứ diện 7.2 Phần tử lục diện 3 Chương 8. Bài toán động lực học kết cấu (TL8) 8.1 Phương trình động lực học kết cấu 8.2 Ma trận khối lượng tương thích và ma trận khối lượng tập trung 8.3 Ma trận khối lượng tương thích trong hệ tọa độ tổng thể 8.4 Dao động tự do – Bài toán trị riêng xác định tần số dao động tự do của kết cấu Chương 9. Một số ví dụ ứng dụng phương pháp phần tử hữu hạn trong phần mềm Matlab-Simulink 9.1. Khảo sát độ bền trục khuỷu động cơ 9.2. Khảo sát độ bền khung vỏ ô tô 9.3 Khảo sát lực cản khí động lực học trên ô tô TÀI LIỆU THAM KHẢO 1. Trần Ích Thịnh, Trần Đức Trung, Nguyễn Việt Hùng (2000), Phương pháp phần tử hữu hạn trong kỹ thuật, Đại học Bách khoa Hà Nội 2. Trần Ích Thịnh, Ngô Như Khoa (2007), Phương pháp phần tử hữu hạn, Giáo trình Đại học Thái Nguyên 3. Ngô Như Khoa (2011), Phương pháp phần tử hữu hạn, Trường Đại học Thái Nguyên 4. Chu Quốc Thắng (1997), Phương pháp phần tử hữu hạn, Nxb Khoa học Kỹ thuật. 5. Trần Vĩnh Hưng (2012), Ứng dụng phần tử hữu hạn. Bài giảng cao học, trường ĐH Sơ phạm Kỹ thuật Hưng Yên, 4 Phần 1 KHÁI NIỆM CHUNG VỀ PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN 1.1 Mô hình rời rạc hóa kết cấu Để tìm ẩn số: chuyển vị, biến dạng, ứng suất tại mỗi điểm bất kỳ trong kết cấu, chi tiết máy, người ta thường áp dụng phương pháp phần tử hữu hạn (PTHH). tiếng Anh: Finite Element Method (FEM). - Ý tưởng cơ bản của phương pháp phần tử hữu hạn toán kết cấu: Coi vật thể liên tục như là tổ hợp của nhiều phần nhỏ liên kết với nhau bởi một số hữu hạn các điểm gọi là nút. - Các phần nhỏ được hình thành gọi là các phần tử hữu hạn (phần tử) - Hình dạng, kích thước các phần tử khác nhau, tạo thành các mạng lưới khác nhau. - Một số sơ đồ rời rạc hóa kết cấu liên tục thành mạng lưới PTHH (hình 1) Hình 1 Một số sơ đồ rời rạc hóa kết cấu - Quan niệm rời rạc hóa là gần đúng Năng lượng bên trong mô hình thay thế bằng năng lượng trong kết cấu thực - Các loại phân tử: Căn cứ vào hình dạng và tình hình chịu lực. + Thanh: Lấy đoạn dầm làm phần tử hữu hạn + Tấm phẳng: Phần tử tam giác, phần tử chữ nhật, tứ giác + Vỏ: Phần tử tấm phẳng, vỏ. + Khối: Phần tử hình tứ diện, hình lập phương, hình lục diện. + Vật thể đối xứng trục: Phần tử hình vành khăn. 5 - Một số loại phần tử hữu hạn thường dùng: Hình 2. Một số loại phần tử thường dùng - Khi phân tích các kết cấu có thể sử dụng các mô hình tính như sau: 1. Mô hình chuyển vị - Chọn chuyển vị ở các nút là giá trị cần tìm. Các giá trị này được xác định từ hệ phương trình cân bằng được thành lập trên cơ sở nguyên lý thế năng toàn phần dừng. - Nguyên lý phát biểu như sau: Trong tất cả các trường chuyển vị thỏa mãn các điều kiện tương thích và điều kiện biên động học, thì trường chuyển vị tương ứng với sự cân bằng của vật thể sẽ làm cho thế năng toàn phần  đạt giá trị dừng (đạt giá trị cực tiểu).  = v + V = 0 (1.1) Trong đó: : v + V: Thế năng toàn phần, là hàm của các chuyển vị v: Thế năng biến dạng đàn hồi của vật thể V: Công của ngoại lực sinh ra trên dịch chuyển của ngoại lực do vật thể biến dạng 6 Hình 3. Thế năng biến dạng đàn hồi Thế năng biến dạng và thế năng bù của vật thể đàn hồi - Nếu hệ ở trạng thái ổn định, thể năng toàn phần có giá trị cực tiểu. - Như vậy sau khi giả thiết một dạng hàm chuyển vị trong phần tử, từ điều kiện dừng của phiếm hàm  ta sẽ nhận được một hệ phương trình cân bằng trong khi các điều kiện liên tục đã được thỏa mãn. 2. Mô hình cân bằng - Chọn các ứng suất hay nội lực ở các nút làm ẩn số. Các ẩn số này được xác định từ phương trình tương thích thành lập trên cơ sở nguyên lý cực tiểu của thế năng bù toàn phần * đạt giá trị dừng * = v* + V* = 0 (1.2) Trong đó: * = v + V* v*: Thế năng bù của biến dạng V*: Công bù của ngoại lực 3. Mô hình hỗn hợp - Coi các ẩn (ứng suất, chuyển vị) là độc lâpj với nhau trong toàn bộ phần tử. Các ẩn này được xác định từ hệ phương trình thiết lập theo nguyên lý biến phân REISSNER. Thông thường người ta hay sử dụng mô hình chuyển vị. Vì nó thuận lợi cho việc tính toán trên máy tính. 1.2 Hàm chuyển vị, Hàm dạng 1.2.1 Hàm chuyển vị 7 Hàm chuyển vị được thiết lập dưới dạng một đa thức. - Bậc của đa thức và số lượng số hạng trong đa thức phụ thuộc vào bậc tự do của phần tử, tức là số chuyển vị ở tất cả các nút của phần tử. Hàm có dạng Bài toàn một chiều: :21)( xf x   tuyến tính 2 321)( xxxf   : bậc hai Bài toán hai chiều yxf yx 321),(   : tuyến tính 2 65 2 4321),( yxyxyxyxf   : bậc hai 1.2.2 Hàm dạng Biểu diễn hàm chuyển vị qua chuyển vị nút, áp dụng bài toán phằng của lý thuyết đàn hồi. - Xét một tam giác: Phần tử có 3 nút là 3 đỉnh của tam giác là khớp nối của các phần tử khác. Hình 4. Phần tử tam giác - Mỗi nút có 2 bậc tự do chuyển vị theo hai phương x, y. Phần tử có 6 bậc tự do, có 6 chuyển vị ở các nút là ui, vi, uj, vj, um, vm: gọi là chuyển vị nút Hợp thành véc tơ chuyển vị nút của phần tử 8                        m m j j i i v u v u v u  (1.3) - Hàm chuyển vị là hàm của các tọa độ, cho phép xác định chuyển vị tại mộ điểm bất kỳ trong phần tử. - Hàm chuyển vị của các phần tử tam giác có dạng       yxv yu yx yx 654),( 321),(   (1.4) Hay                                       m m j j i i yx yx v u v u v u xy xy v u f . 0001 0001 ),( ),( (1.5) Viết gọn lại {f} = [Q]. {} (1.6) Trong đó: {f}: là véc tơ chuyển vị [Q]: là ma trận các đơn thức {}: là véc tơ các hệ số Chuyển vị tại các nút, theo (1.6) ta có: {} = [C] .{} (1.7) Trong đó [C] là giá trị [Q] tại các nút, tức là ma trận tọa độ nút Xác định{} theo [C] từ (1.7) {} = [C]-1.{} (1.8) Thay (1.8) vào (1.6) ta có {f} = [Q].[C]-1{} (1.9) 9 Hay {f} =[N].{} (1.10) Trong đó [N] = [Q].[C]-1 (1.11) [N] gọi là ma trận các hàm dạng. Hay gọi là ma trận các hàm nội suy. Ta có thể suy ra chuyển vị của điểm bất kỳ 1.2.3 Lực nút - Khi vật thể chịu lực, trong các phần tử sinh ra nội lực, PPPTHH giả thiết các nội lực nằm ở điểm nút gọi là lực nút. Đó là lực tương tác giữa các phần tử liên kết với nhau tại nút do chuyển vị nút sinh ra. Tại các nút có thể có ngoại lực (tải trọng). - Trong mỗi phần tử các lực nút hợp thành véc tơ lực nút{F} có số thành phần bằng số thành phần của véc tơ chuyển vị nút, được sắp xếp tương ứng với véc tơ chuyển vị nút. Hình 5. Véc tơ chuyển vị nút của phần tử {F} =[vi Vi vj Vj vm Vm] T Đối với thanh chịu uốn Véc tơ chuyển vị nút {} = [vi i vj j] Véc tơ lực nút {F}e = [Vi Mi Vj Mj] T 10 1.3 Phương trình cơ bản của phương pháp PTHH 1.3.1 Các quan hệ chuyển vị, biến dạng, ứng suất trong phần tử - Chuyển vị ở các nút: Trường chuyển vị theo chuyển vị nút {f} = [N] .{} (1.12) - Biến dạng theo lý thuyết đàn hồi {} = [] .{f} (1.13) [] toán tử vi phân                                                        zz yy yx z y x 0 0 0 00 00 00 (1.14) Ta có véc tơ biến dạng {} =[] . [N] .{f} (1.15) Hay{} =[B] .{} (1.16) Trong đó [B] = [] .  N gọi là ma trận biến dạng - Ứng suất tại một điểm trong phần tử, xác định theo định luật Hooke: {} = [D] .{} (1.17) [D] gọi là ma trận đàn hồi Thay (1.16) vào (1.17) ta có {} = [D].[B]{} (1.18) Hay {} =[S].{} (1.19) Trong đó: [S] = [D] . [B] (1.20) gọi là ma trận tính ứng suất 11 1.3.2 Phương trình cơ bản của phương pháp phần tử hữu hạn (Ma trận độ cứng phần tử và véc tơ tải phần tử) - Sử dụng nguyên lý cực tiểu thế năng toàn phần (mô hình cân bằng) để thiết lập phương trình cơ bản của phương pháp PTHH. - Giả sử một PTHH có thể tích Ve chịu tác dụng của lực thể tích p và lực bề mặt q trên diện tích Se - Thế năng toàn phần tử là e viết dưới dạng               Ve Se TTT ve e dSqfdVpfdV .... 2 1  (1.21) Thay (1.21), (1.15), (1.19) vào (1.21) ta có:                    Se T Ve TTTT Ve e dSqdVpNdVBDB ..... 2 1  (1.22) Biến đổi ta có                            Ve Se T Ve TTTT e dSqNdVpNdVBDB .... 2 1  (1.23) Đặt        Ve T dVBDBK (1.24) gọi là ma trận độ cứng phần tử            Ve Se TTe dSqNdVpNp ... (1.25) gọi là véc tơ tải phần tử bao gồm các thành phần lực đặt tại nút. Các lực này được quy đổi sau khi dời các tải trọng p và q về nút.  ep còn gọi là lực nút tương đương Viết gọn lại ta có        eTTe pK  . 2 1  (1.26) Trong trường hợp ở các nút có tồn tại lực tập trung 12                        n e E R R R . . . 2 1 Thì phải cộng thêm các lực tập trung vào véc tơ tải  ep - Theo nguyên lý cực tiểu thế năng toàn phần. Điều kiện cân bằng tại các nút của phần tử là:   0    e (1.27) Tức là       0;....;0;0 21          n eee  (1.28) Sau khi lấy cực tiểu từ (1.26) ta được     Kp e  (1.29) - Đây là phương trình cơ bản của phương pháp phần tử hữu hạn tính theo mô hình chuyển vị. - Điều này có nghĩa là tại từng nút, lực nút do chuyển vị gây ra      KF e  phải cân bằng tải trọng đặt ở nút    ee pF  - Trong trường hợp PTHH có biến dạng ban đầu 0 và ứng suất ban đầu 0 thì công thức (1.18) được bổ sung thành             00..   DBD (1.30) Do đó công thức véc tơ tải phần tử (lực nút tương đương (1.25)) có thêm thành phần 0 và 0 gây ra:                      Ve Se Ve Ve TTTTe dVBdVBdSqNdVdVpNp 00 ......  (1.31) 13 1.3.3 Ma trận độ cứng tổng thể (Véc tơ tải tổng thể, phương trình cơ bản của hệ) - Sau khi thiết lập được các ma trận độ cứng phần tử và véc tơ tải phần tử của tất cả các phần tử trong mạng lưới kết cấu. Ta cần phải tổ hợp tất cả chúng lại thành ma trận độ cứng tổng thể [K] và véc tơ tải tổng thể [P] của kết cấu, từ đó xây dựng phương trình cơ bản đối với toàn bộ kết cấu. - Việc tổ hợp này có nghĩa là sắp xếp các thành phần trong các ma trận [k] của các phần tử vào các vị trí thích hợp trong ma trận [K] và các thành phần trong các ma trận{p} của các phần tử vào các vị trí thích hợp trong [p]. - Sự sắp xếp này được mô tả bằng ma trận định vị của các phần tử - Gọi véc tơ chuyển vị nút của phần tử là   , véc tơ chuyển vị nút tổng thể của toàn bộ kết cấu là   ta có quan hệ           . 1 endxnndx L (1.32) Trong đó  eL :là ma trận định vị của phần tử nd: số chuyển vị nút trong mỗi phần tử n: số chuyển vị nút trong toàn bộ kết cấu. 14 Phần 2 CÁC PHƯƠNG PHÁP GIẢI BÀI TOÁN PHẦN TỬ HỮU HẠN 2.1 Phương pháp sai phân PTHH Ví dụ: Giải HPTVP tìm hàm y = f(x) Hình 2.1 Ví dụ hàm y = f(x) Xét trường hợp trong khoảng [0, x], điều kiện đầu vào:       0)0( )(' y xxy (2.1) a) Phương pháp giải tích Ta giải dễ dàng bằng giải tích, ta có: 2 2 1 )( xxy  b) Phương pháp sai phân PTHH Chia thành n đoạn trong khoảng [0, x] ta có: o ooxx o x x yy x yy xy        1 )()1(' )( Y X y = f(x) yo yn xo x1 xn ∆x ∆x .. Các điểm nút Các phần tử α y1 15 1 12)1()2( 1 ' )( x x yy x yy xy xx        ... tgx x yy x yy xy n nnxnxnx n         1)()1(' )( Tổng quát: x yy xy k xxkx k     )()(' )( (2.2) Trở lại bài toán trên ta có: Tại điểm xo : oo yxxyy  ).( ' 1 Tại điểm x1 : 11 ' 2 ).( yxxyy  ... Tại điểm xn-1 : 11 ' ).(   nnn yxxyy Tại điểm xn : nnn yxxyy  ).( ' 1 Nối các điểm y1, y2.......yn+1 với nhau ta xác định được dạng hàm số y = f(x) Khi chia khoảng ∆x nhỏ => đường nối “gần” hơn => kết quả chính xác hơn Khi chia khoảng ∆x lớn => đường nối “xa” hơn => kết quả không chính xác. Tuy nhiên việc chia nhỏ là bao nhiêu khoảng còn phụ thuộc vào dung lượng và khả năng của máy tính điện tử nữa. c) Lập trình trong Matlab file *.m Function y=PTVPbac1(f,f0,x0,x1,n) % n la khoang chia Sym x; f=inline(f); Delta_x=(x1-x0)/n; % xay dung mang tu x0 den xn+1 X=[x0:delta_x:x1]; Y=zeros(side(X)); %dua vao hàm zeros de tang toc van de su ly cua phan mem 16 Y(1)=f0; For i=2:n+1 Y(i)=f(X(i-1))*delta_X+Y(i-1); End y=Y; % ve do thi line (X,Y); 2.2 Phương pháp năng lượng và PP Niu tơn-Lagrăng Có hai cách giải:  Xây dựng được phiến hàm mô tả kết cấu: Ví dụ trong cơ học biến dạng đó là phiến hàm năng lượng. Khi kết cấu ổn định (trạng thái đạt được trong giai đoạn đàn hồi của kết cấu khi biến dạng) thì năng lượng cực tiểu.  Có hệ phương trình cho trước: theo Niu tơn hoặc Lagrăng. 2.2.1 Cách 1: Xây dựng phiến hàm Ví dụ 1 Có hai điểm A và B cho trước, khảo sát khoảng cách ngắn nhất giữa A và B. Hình 2.2 Ví dụ về phiến hàm Rõ ràng khoảng cách ngắn nhất là đoạn thẳng nối hai điểm A và B. Ta gọi phiến hàm là đường thẳng AB . Xây dựng phiến hàm: - Tính diện tích vùng S(aABb), ta có Y X A B y = f(x) a b 17   b a xyxyS 2' ))((1)( - Cực tiểu S Ta có Smin khi y(x) là đường thẳng S => là phiến hàm của hàm y(x) => tìm được hàm y(x) 2.2.2 Cách 2: Giải theo hệ PTVP cho trước Ví dụ 2 Cho HPHVP (1):       0)0( )(' y xxy Trên hình 1, chia khoảng xo, x1, ..., xn Phần tử thứ nhất tương ứng với nút xo, x1 Phần tử thứ hai tương ứng với nút x2, x3 ... Phần tử thứ n tương ứng với nút xn-1, xn Xấp xỉ mỗi phần tử như sau: Giả sử ta sử dụng hàm bậc nhất để xấp xỉ: y = ax +b (2.3) y được gọi là hàm xấp xỉ khi đó là gần đúng, ta ký hiệu là y* Phần tử thứ nhất : y* = a1 x +b1 Phần tử thứ hai : y* = a2 x +b2 .... Phần tử thứ n : y* = an-1 x +bn-1 Các hệ số a1, a2,...,an-1, b1,b2,...,bn-1 là các giá trị của hàm tức là bậc tự do của hàm. a) Xét một phần tử thứ nhất Trên mặt phẳng hệ trục tọa độ, theo trục tung có: yo = a1 xo +b1 y1 = a1 x1 +b1 (2.4) Viết dưới dạng ma trận: 18                   11 1 1 1 1 y y b a x x oo (2.5) Ta cần tìm a1,b1 :                    1 1 11 1 1 1 y y x x b a oo Khi y* là chính xác y *= y, ta thấy: yyy  * (2.6) Với ∆y là nhỏ nhất Lấy đạo hàm phương trình (6) ta có: ''*' )()( yyy  (2.7) Và ))(()( 1 ' xfay  (2.8) Xét trong khoảng xo – x1 của phần tử thứ nhất: 0..)( 0..)( 1 1 ' 1 '     dxby dxxy x ox x ox (2.9) Thay (8) vào (9) ta được: 0..))(( 0..))(( 1 1 1 1 1     dxbxfa dxxxfa x ox x ox (2.10) Từ phương trình (10) ta tìm được a1,b1 của phần tử thứ nhất Tương tự ta cũng giải được đối với các phần tử khác. b) Trường hợp tổng quát - Chọn hàm xấp xỉ )(.....)()( 2211 * xfaxfaxfay nn (2.11) Ta có: )(.....)()()( ''22 ' 11 '* xfaxfaxfay nn (2.12) 19 Thay vào phương trình (1) ta có: 0.)())(( ..... 0.)())(( 0.)())(( 1 ''* 1 2 ''* 1 1 ''*          dxxfyy dxxfyy dxxfyy kx kx n kx kx kx kx (2.13) Và: ''*)( yy  là phần dư )(xfn là trọng số Trong đó )(xfn và y’ đã biết, y *’ là hàm xấp xỉ ta đã biết dạng hàm của nó từ việc chọn ban đầu theo phương trình (2.11). Các kết quả mỗi tích phân (2.13) cho ta một biểu thức. Ghép n biểu thức với nhau ta giải tìm được a1, a2,..., an từ đó ta tìm được hàm xấp xỉ (2.11). Ví dụ 3 Giải HPTVP xyyy  ''' (2.14) 0)0( 0)0( '   y y (2.15) Trong khoảng [0 , 10] Hệ phương trình này đã có lời giải bằng giải tích. Hình 2.3 Khoảng giải [0 , 10] gọi là miền không gian - Chọn hàm xấp xỉ 0 1 2 10 20 Ta chọn hàm xấp xỉ theo công thức (2.11) )(.....)()( 2211 * xfaxfaxfay kk (2.16) Trong đó f1, f2, ... fk đã biết trước Theo định nghĩa, ta chọn (tùy chọn) hàm xấp xỉ theo đa thức bậc 3, phương trình (2.17), ta chọn bậc càng cao càng tốt. Việc chọn bậc cao phụ thuộc vào khả năng của máy tính. 43 2 2 3 1 axaxaxay  (2.17) Trong đó các phần tử a1, a2, a3, a4 là khác nhau. Ta có: 32 2 1 ' 23 axaxay  (2.18) 21 '' 26 axay  (2.19) Từ phương trình (14) ta có: 0'''  xyyy (2.20) Thay (2.18), (2.19 vào (2.20) ta được: 0)2()162()3( 234123 2 12 3 1  aaaxaaaxaaxa (2.21) Ta tìm a1, a2, a3, a4 theo phương pháp Galenkin (phần dư có trọng số). Đặt 11 aA  122 3aaA  162 1233  aaaA (2.22) 2344 2aaaA  Ta có: 043 2 2 3 1  AxAxAxA (2.23) Theo (2.13) ta tính hàm trọng số: 0).( 1 3 43 2 2 3 1    kx kx kkkk dxxAxAxAxA 0).( 1 2 43 2 2 3 1    kx kx kkkk dxxAxAxAxA 21 0).( 1 43 2 2 3 1    kx kx kkkk xdxAxAxAxA (2.24) 0).( 1 43 2 2 3 1    kx kx kkkk dxAxAxAxA Thực hiện tiếp ta có: 0 4 1 5 1 6 1 7 1 1 4 4 1 5 3 1 6 2 1 7 1   kx kx k kx kx k kx kx k kx kx k xAxAxAxA 0 3 1 4 1 5 1 6 1 1 3 4 1 4 3 1 5 2 1 6 1   kx kx k kx kx k kx kx k kx kx k xAxAxAxA 0 2 1 3 1 4 1 5 1 1 2 4 1 3 3 1 2 1 5 1   kx kx k kx kx k kx kx k kx kx k xAxAxAxA (2.25) 0 2 1 3 1 4 1 1 4 1 2 3 1 3 2 1 4 1   kx kx k kx kx k kx kx k kx kx k xAxAxAxA Hệ phương trình (2.25) có 4 ẩn số nhưng là hệ suy biến nên chưa giải được ngay vì chưa có điều kiện đầu. Các hệ số kkkk AAAA 4321 ,,, là không độc lập tuyến tính, cần ghép tất cả các phương trình với nhau cùng điều kiện đầu để giải. Triển khai tiếp hệ phương trình (2.25) đưa về dạng chung ( 1 kk xx ) làm cơ sở để lập trình, ta được hệ phương trình (2.26). 0)( 4 1 )( 5 1 )( 6 1 )( 7 1 4 1 4 4 5 1 5 3 6 1 6 2 7 1 7 1   kk k kk k kk k kk k xxAxxAxxAxxA 0)( 3 1 )( 4 1 )( 5 1 )( 6 1 3 1 3 4 4 1 4 3 5 1 5 2 6 1 6 1   kk k kk k kk k kk k xxAxxAxxAxxA 0)( 2 1 )( 3 1 )( 4 1 )( 5 1 2 1 2 4 3 1 3 3 4 1 4 2 5 1 5 1   kk k kk k kk k kk k xxAxxAxxAxxA (2.26) 0)()( 2 1 )( 3 1 )( 4 1 14 2 1 2 3 3 1 3 2 4 1 4 1   kk k kk k kk k kk k xxAxxAxxAxxA Ký hiệu: P p k p k Xxx  1 (với p = 17) ta viết (2.26) dưới dạng ma trận (2.27) 22                     1234 2345 3456 4567 2 1 3 1 4 1 2 1 3 1 4 1 5 1 3 1 4 1 5 1 6 1 4 1 5 1 6 1 7 1 XXXX XXXX XXXX XXXX                   k k k k A A A A 4 3 2 1 =                   0 0 0 0 (2.27) Đặt D1 =                     1234 2345 3456 4567 2 1 3 1 4 1 2 1 3 1 4 1 5 1 3 1 4 1 5 1 6 1 4 1 5 1 6 1 7 1 XXXX XXXX XXXX XXXX (2.28) là ma trận đối xứng Ta đã chọn hàm xấp xỉ là đa thức bậc 3 theo phương trình (2.17) và đạo hàm bậc nhất theo (2.18), kết hợp với điều kiện đầu theo phương trình (2.15), trường hợp này thì : 0 0 3 4   a a (2.29) Từ (2.22), thay vào ta được (2.30)                   k k k k A A A A 4 3 2 1 =                   1120 0126 0013 0001                   4 3 2 1 a a a a +                    0 1 0 0 (2.30) Đặt D2 =                   1120 0126 0013 0001 (2.31) 23 Đặt D3 =                    0 1 0 0 (2.32) Thay thế (2.32) vào (2.27) và chú ý điều kiện đầu (2.29) ta được (2.33) 0: 3 4 3 2 1 21                                       D a a a a DD (2.33) 031 4 3 2 1 21                    DD a a a a DD (2.34) Phương trình (2.34) là phương trình bậc nhất 4 ẩn trong đó 2 ẩn đã biết theo (2.29) Ký hiệu:  kDD 21 (2.35)  PDD 31 (2.36) Và lập trình ma trận theo tọa độ của phần tử ta được (2.37)      0 4 3 2 1                    P a a a a k (2.37) Hay    P a a a a k                    4 3 2 1 (2.38) 24 Ví dụ 4 Giải hệ PTVP:                   44434241 34333231 24232221 14131211 aaaa aaaa aaaa aaaa                    4 3 2 1 x x x x                   4 3 2 1 b b b b (2.39) Ta biến đổi:                    4 3 2 1 x x x x                   00 00 10 01                            4 32 1 0 0 x xx x (2.40) Trong đó: 211 .01. xxx  Thay (2.40) vào (2.38) ta được:  k                                                                4 32 1 0 0 00 00 10 01 a aa a =  P (2.41)  k                           2 1 00 00 10 01 a a =  P                    4 3 0 0 a a (2.42) Phương trình (2.42) có dạng:  1k         2 1 a a =  P (2.43) Hay         2 1 a a =  1k .  P (2.44) Lập trình Matlab 25 D1=matranD1 (10,0); %X1=(xk-xk1); D2=[1 0 0 0; 3 1 0 0; 6 2 1 0; 0 2 1 1]; k=D1*D2; D3=[0; 0; -1; 0]; P=D1*D3; K1=k*[1 0; 0 1; 0 0; 0 0]; K11=K1([1, 2], [1,2]); P1=P([1,2]); NO=K11^-1*(-P1) %giai ra ta duoc NO= -0.0111 0.2019 Vậy ta có: a1 = -0.0111 a2 = 0.0219 Với a3 = 0 a4 = 0  Xấp xỉ tốt nhất ta được: y = - 0.111x3 + 0.2019 x2 Thử lại 02019.02)0111.0(6 2019.02)0111.0(3 '' 2'   xy xxy *Chú ý: Khi chưa biết hệ PT vi phân ta phải xây dựng phiến hàm Trong cơ học => sử dụng phiến hàm năng lượng (biến dạng) => sử dụng nguyên lý Hamiltơn. 26 Trong trường hợp tổng quát, hàm năng lượng có dạng:   AU (2.45) U là thế năng biến dạng đàn hồi (là công của nội lực: lực thế, ứng suất...) A là công của ngoại lực (lực khối, lực bề mặt, lực tập trung...) + Xác định thế năng biến dạng đàn hồi U:   dVU V  )(2 1 (2.46)   là ten xơ biến dạng – là một véc tơ (ma trận vuông)   là ten xơ ứng suất (V) được xét trên toàn vật thể Trong không gian 3 chiều ta có:    zxyzxyzyx   (2.47)                              zx yz xy z y x        (2.48) Thay (2.47), (2.48) vào (2.46) ta có: dVU V zxzxyzyzxyxyzzyyxx  )( )( 2 1  (2.49) + Xác định công của ngoại lực A: Từ công thức tính công suất : iVFP  => công của lực F là:    2 1 2 1 t t i t t i sdFdtVFA (2.50) dtVF i là tích vô hướng Trong không gian 3 chiều thì: 27                z i y i x i z i y i x i u u u F F F Nên : z i z i y i y i x i x i uFuFuFA ...  (2.51) Trong đó chỉ số K đại diện cho nhiều loại lực như áp suất, tải trọng động...  Mối liên hệ giữa ứng suất biến dạng:      C (2.49)  C là ma trận vật liệu (nói chung) – đàn hồi, đẳng hướng và biết trước đối với từng loại vật liệu. Các chỉ số trong đó của vật liệu cơ bản bằng “0”, chỉ còn lại là hệ số Poát xông “e” và mô dun đàn hồi “ ” phụ thuộc vào từng loại vật liệu.  Mối liên hệ giữa biến dạng và chuyển vị: Biến dạng của phần tử liên quan đến hình dáng . Chuyển vị của nút liên quan đến vị trí. Trường hợp tuyến tính mối liên hệ giữa biến dạng và chuyển vị như sau: u w x w x u xyx           z u y u y u yzy           (2.50) x v z v z u xyz           Vậy ta cần xấp xỉ u, v, w bằng hàm có dạng đã biết trước và thực hiện theo các bước sau: Bước 1: - Tính từng thông số , ... - Tính phiến hàm thế năng - Tính chuyển vị, biến dạng Bước 2: 28 - Xấp xỉ hàm chuyển vị u,v, w theo dạng hàm nào đó (thông thường là đa thức bậc ba) - Tính biến dạng và tính ứng suất - Tính thế năng (ra được hàm có các hệ số là xấp xỉ) - Cực tiểu hàm thế năng  Các hệ số là nghiệm của bài toán Ví dụ 4 Xác định ma trận độ cứng phần tử: Xét thanh chịu uốn trên hình 3 Hình 2.4 Thanh chịu uốn Cần tìm chuyển vị uốn vA và chuyển vị góc αA => có hai bậc tự do (nếu khả năng giải được nhiều ta chọn càng nhiều bậc tự do càng tốt nhưng phụ thuộc vào khả năng giải của máy tính). Từ công thức (2.50), bỏ qua chuyển vị theo phương y (sẽ gây ra sai số), ta có: x v x u y u x u xyx              0 00           y v z v y v yzy  (2.51) 00           x w z w z w zyz  Gọi    ndD (2.52) y x u A’ A v o, w αA 29            1 i n v d (2.53) N là số nút Từ phương trình (2.51), (2.52), (2.53) có:                                      A Av   00 00 10 00 00 00 (2.54) Thay (2.54) vào (2.46)   dVU V  )(2 1 ta được:    dVdDCDdU n TT V n 1 )(2 1   (2.55) DCDT 1 = ke là ma trận độ cứng phần tử Trong đó C là ma trận vật liệu (đã biết) D là hàm dạng xấp xỉ => chọn bậc tự do là gì (chuyển vị, góc xoay...) Nếu các phần tử giống nhau thì ke giống nhua và ngược lại. Vì vậy ta cố gắng chia các phần tử giống nhau để giảm bớt việc tính toán các ma trận độ cứng phần tử ke. Ví dụ 5 Thanh chịu kéo nén đúng tâm Xét thanh chịu kéo nén, hình 2.5 Hình 2.5 Thanh chịu kéo nén x F S l A B 30 + Bằng giải tích, ứng suất và biến dạng trong thanh được xác định bởi: SE F E S F .      Tại điểm bất kỳ trên chiều dài l của thanh thì: X SE F dxU x x .0   (2.56) + Bằng phương pháp PTHH Đề dễ hiểu ta chia thanh thành một phần tử và chỉ tồn tại ứng suất pháp x và chuyển vị xU Theo lý thuyết biến dạng đàn hồi thì: dVU xx  2 1 Theo định lý Cô si về kéo nén: E x x    , ta có Ta có: dx ES dxSEdVEU l l xx l xx   0 0 22 0 2 22 1 2 1  (2.57)  ?x => Ta xấp xỉ xU Chọn hàm bậc nhất cho đơn giản (chọn hàm nào thì bậc tự do đó): UA, UB là bậc tự do của phần tử AB; UA là bậc tự do của nút A; UB là bậc tự do của nút B. Giả sử ta chọn hàm xấp xỉ có dạng : xaaU x 21  xU gọi là hàm dạng Lấy điều kiện đầu ta có: 00. 121  aaaU A laalaaU B .. 221  Từ phương trình: x ux x    , đạo hàm ta được: l U a Bx  2 ' (2.58) 31 + Xác định thế năng biến dạng đàn hồi U: Từ (57) ta có: l UES l l UES dx l UES U BB l B x 2 . . 22 2 2 22 0         (2.59) + Xác định công của ngoại lực A: BUFA . (2.60)  Thế năng biến dạng đàn hồi (theo công thức (45)): B B UF l USU . 2 .. 2  (2.61) Ta thấy min khi 0   BU => 0 F l ESU B => l ES F ES lF U B  . Hay : x ES F U x  (2.62) So sánh kết quả với PP giải tích (2.57) thì hai kết quả giống nhau. Tuy nhiên đây là bài toán đơn giản, thanh chịu kéo nén được chai thành một phần tử. Khi chia thanh thành nhiều phần tử thì giải bằng giải tích sẽ rất phức tạp, để thay thế giải tích thường ta phải tiến hành thí nghiệm, còn tính toán bằng PP PTHH sẽ có nhiều ưu điểm hơn. Trình tự giải bằng PP PTHH 1) Chọn bậc tự do của nút (lưu ý chọn bậc tự do của nút phải phù hợp). Tùy theo yêu cầu của bài toán cần tính thông số nào, thông thường là chuyển vị của nút và góc xoay (đạo hàm của chuyển vị theo phương x, y). Bước này còn gọi là rời rạc hóa kết cấu. 2) Xấp xỉ các biến chuyển vị bằng một hàm biết trước (lưu ý số tham số tự do (số bậc tự do) của dạng hàm phải phù hợp với bậc tự do của nút). Ví dụ 6 32 Hình 2.6 Xấp xỉ hàm v = F(vA, αA, vB, αB) Với 4 thông số ta cần 4 phương trình để giải Giả sử ta chọn hàm xấp xỉ dạng: 4 5 3 4 2 321 xaxaxaxaav  (2.63) Ta có: 3 5 2 432 432 xaxaxaa x v     (2.64)  3 5 2 432 3 5 2 432 4 5 3 4 2 321 4 5 3 4 2 321 432 432 BBBB AAAA BBBBB AAAAA xaxaxaa xaxaxaa xaxaxaxaav xaxaxaxaav       (65) Là hệ 4 phương trình 5 ẩn số a1, a2, a3, a4, a5 không giải được. Như vậy việc chọn hàm xấp xỉ (2.63) là không hợp lý, ta cần bỏ thành phần 45xa đi. 3) Xác định ma trận hàm dạng (D ở ví dụ trên). 4) Xác định ma trận độ cứng Ke của phần tử. 5) Xác định véc tơ lực nút (tính công ngoại lực A) 6) Tổng hợp ma trận độ cứng và tổng hợp véc tơ lực nút cho toàn kết cấu 7) Giải các ma trận tổng hợp dạng:     FdK là phương trình đại số, nên giải được dễ dàng. * Tính thế năng biến dạng đàn hồi của dầm chịu uốn (chỉ đúng với dầm chịu uốn) Xem ví dụ hình 1.13 trang 27 (PP PTHH Chu Quốc Thắng) + Lấy kết quả đã có (giải bằng đại số): A B vB, αB vA, αA 33  l x dx EJ M U 0 2 2 1 Với: '' 2 2 .yEJ dx vd EJM  Vậy:           l dx dx vd EJU 0 22 2 1 (2.66)  l qvdxA 0 (2.67)  Thế năng biến dạng toàn phần:           l l qvdxdx dx vdEJ AU 0 0 22 2 (2.68) + Giải bằng PP PTHH Chia dầm thành 2 phần tử, hình 2.7 Hình 2.7. Dầm chịu uốn có hai phần tử với 3 nút Các bước giải như sau: 1) Rời rạc hóa kết cấu a) Chọn bậc tự do của nút Ta chọn hai bậc tự do cho mỗi nút là độ võng y và góc xoay α Với: x y y x y y x y y          3 33 2 22 1 11 ; ; ;    (2.69) b) Thế năng biến dạng đàn hồi 1 3 (y3, (y1, 2 (y2, l 34           2/ 0 22 1 2 1 l dx dx vd EJU (2.70)           l l dx dx vd EJU 2/ 22 2 2 1 (2.71) 2) Xấp xỉ các biến chuyển vị bằng một hàm biết trước Ta thấy trên mỗi phần tử có 4 bậc tự do, ta chọn hàm xấp xỉ bậc 3. 34 2 321)( xaxaxaaxy  (2.72) Biểu diễn y(x) theo các bậc tự do của nó Trên phần tử 1: Chuyển vị tại nút 1:                      4 3 2 1 32 1 1 a a a a xxxy (2.73) Chuyển vị tại nút 2:                      4 3 2 1 2 0001 a a a a y (2.74) Đạo hàm (2.72) ta được: 2432 32 xaxaa x y    (2.75)                         4 3 2 1 23210 a a a a xx x y (2.76) Thay tọa độ x1 =0 vào ta được: 35                        4 3 2 1 1 0010 a a a a x y (2.77) Thay tọa độ x2 = l/2 vào ta được:                                        4 3 2 1 32 2 222 1 a a a a lll y (2.78)                                    4 3 2 1 2 2 2 310 a a a a l l x y (2.79) Tổng hợp lại ta được:                                                                               3 2 1 2 32 2 2 1 1 2 310 222 1 0010 0001 a a a a l l lll x y y x y y o (2.80) Ký hiệu ma trận tọa độ:                                       2 32 1 2 310 222 1 0010 0001 l l lllX (2.81) 36                                            x y y x y y X a a a ao 2 2 1 1 1 1 3 2 1 (2.82)                            x y y x y y Xxxxy 2 2 1 1 1 1 321 (2.83)                            x y y x y y Xxxy 2 2 1 1 1 1 2' 3200 (2.84)                            x y y x y y Xxy 2 2 1 1 1 1 '' 6200 (2.85) Nghịch đảo y” ta có                                   x X x y y x y yy TT 6 2 0 0 1 1 2 2 1 1 '' (2.86)                                                                x y y x y y Xx x X x y y x y yy TT 2 2 1 1 1 1 1 1 2 2 1 1 2 '' 6200 6 2 0 0 (2.87) Trong đó: 37 + Ma trận tọa độ phần tử  TX 11 (tương đương với  TB của Chu Quốc Thắng). Tại phần tử thứ k: 1 2 32 2 11 3 1 2 11 3210 1 3210 1                       kk kkk kk kkk kk xx xxx xx xxx XB (2.88) Ma trận tọa độ của các phần tử khác nhau là khác nhau. + Ma trận bậc tự do nút của phần tử (trang 43)                          x y y x y y qe 2 2 1 1 (2.89) Và tại phần tử thứ k:                            x y y x y y q k k k k k e 1 1 (2.90) 3) Xác định ma trận hàm dạng  x x D 6200 6 2 0 0                    (2.91) Ma trận hàm dạng của các phần tử như nhau sẽ giống nhau. 4) Xác định ma trận độ cứng Ke của phần tử Dạng tổng quát ma trận độ cứng phần tử: dVBDBK V T e  )( (2.92) 38 Ta có ma trận độ cứng phần tử:             dxXDXK dxXDXK Tl l e Tl e 1 1 2/ 1 1 2 1 1 2/ 0 1 1 1     (2.93) Khi đó ma trận độ cứng phần tử thứ k   dxBDBK K V T K k e  )( (2.94) T KB , KB là hằng số, cần tính ma trận hàm dạng D 5) Tính véc tơ lực nút (công ngoại lực A) Xác định lực thuộc phần tử nào? Vì dạng phần tử giống nhau, hàm xấp xỉ trong mỗi phần tử giống nhau và có dạng (2.72). Viết dưới dạng ma trận, hàm xấp xỉ có dạng (93).    eqBxxxy 321 (2.95) Giả sử điểm đặt lực tại tọa độ l/4 thuộc phần tử thứ nhất, ta có:   11 32 444 1 eqB lll y                      (2.96) (chỉ số 1 biểu thị thuộc phần tử thứ nhất) Khi đó công của lực F là:   11 32 444 1. eqB lll FA                      (2.97) Ta có véc tơ lực nút phần tử:                                       3 2 1 11 4 4 4 1 l l l FBqN T T ee (2.98) Thay (2.98) vào (2.97) ta được công của ngoại lực: 39                                                                     3 2 1 2 2 1 1 2 2 1 1 1 4 4 4 1 l l l FB x y y x y y x y y x y y A Te (2.99) 6) Tổng hợp ma trận độ cứng và tổng hợp véc tơ lực nút cho toàn kết cấu (dạng:     FdK là phương trình đại số) Khi chia phần tử ta thấy cứ hai phần tử có chung một nút, ta không để xuất hiện hai bậc tự do ở đây nên ta sẽ tạo một bậc tự do xuất hiện một lần và duy nhất (unique). Từ phương trình (2.89) và (2.90) ta thấy bậc tự do sẽ xuất hiện hai lần tại nút 2 nối hai phần tử: Trên phần tử thứ nhất:                          x y y x y y qe 2 2 1 1 1 (2.100) Trên phần tử thứ hai:                          x y y x y y qe 3 3 2 2 2 (2.101) Xuất hiện x y vày   2 2 hai lần trên kết cấu thanh => cần tổng hợp lại để có véc tơ  uniqueeq là duy nhất. Cách làm: Tìm một ma trận H nhân với  eq để được  uniqueeq    euniquee qHq . (2.102) + Dạng của ma trận H 40                                      100000 010000 001000 000100 001000 000100 000010 000001 H (2.103) + Mã hóa ma trận được nhân Sau khi tổng hợp hai véc tơ  1eq và  2eq bình thường ta được  eq có dạng:                                                x y y x y y x y y x y y qe 3 3 2 2 2 2 1 1 (2.104) Máy tính không nhận được trực tiếp các ký hiệu toán học này => mã hóa các ký hiệu thành: 41                                        23 13 22 12 22 12 21 11 eq (2.105) Cột thứ nhất là thứ tự vị trí nút, cột thứ hai là thứ tự bậc tự do trên nút đó + Nhân (2.103) với (2.105) được ma trận tổng hợp  uniqueeq duy nhất                                                                              23 13 22 12 22 12 21 11 100000 010000 001000 000100 001000 000100 000010 000001 . euniquee qHq (2.106) Tương tự ta có: Ma trận độ cứng tổng hợp: HKHK e T uniquee **  (2.107) Ma trận véc tơ lực nút tổng hợp:  e T uniquee NHN * (2.108) 7) Giải : các ma trận tổng hợp dạng:     FdK là phương trình đại số, nên giải được dễ dàng. 42        ĐKB NqK uniqueeuniqueeuniquee (2.109) ĐKB được cho bởi các điều kiện biên ban đầu, loại bở những thông số nào đã biết, còn lại ta giải dễ dàng. ********** 43 2.4 Phương pháp giải trên phần mềm Ví dụ: Khảo sát ứng suất và chuyển vị của tấm sàn ô tô khách 2.4.1 Giới thiệu PTHH trong MatLab 44 2.4.2 Giới thiệu PTHH trong MatLab 45 **********