Luận văn Phân tích động lực học của khung dầm FGM chịu tải trọng động đất bằng phương pháp phần tử hữu hạn

.1) trong đó: M và Q lần lượt là khối lượng và trọng lượng của kết cấu công trình 0,maxx : gia tốc cực đại của nền đất dưới chân công trình. g: gia tốc trọng trường. Ks: hệ số địa chấn. Không lâu sau khi lý thuyết tính toán tĩnh này xuất hiện, việc phân tích tác động của động đất lên kết cấu công trình đã làm nổi lên một số nhược điểm trong phương pháp của Omori. Đầu tiên, có rất ít kết cấu có thể được xem là cứng tuyệt đối. Khi nền đất chuyển động, đa số các công trình đều biến dạng nên chuyển vị và gia tốc tại các vị trí khác nhau trên công trình là khác nhau, thậm chí còn lớn hơn chuyển vị và gia tốc nền đất dưới chân công trình. Tiếp theo, chu kỳ dao động tự nhiên của hệ kết cấu trùng hoặc gần trùng với chu kỳ dao động của nền đất thì có thể xảy ra hiện tượng cộng hưởng làm hiệu ứng tác động của động đất tăng lên nhiều lần. Tiếp đến là chưa xét đến độ cản của kết cấu trong quá trình dao động. Năm 1920, nhà khoa học Nhật Mononobe đã đề nghị đưa các tính chất biến dạng của kết cấu vào trong tính toán tác động động đất. Ông xem kết cấu như một hệ có một bậc tự do dao động không có lực cản và giả thiết trong thời gian xảy ra động đất, nền đất chuyển động theo quy luật điều hòa sau [4]: 0 0,max( ) sin( )x t x t (1.2) Trong phương pháp của mình, Mononobe chỉ xét tới phần dao động cưỡng bức của hệ kết cấu và đã thu được hệ số khuyếch đại động có dạng sau: 5 2 2 0 1 1 T T    (1.3) trong đó: T, T0 lần lượt là chu kỳ dao động của kết cấu và chu kỳ dao động của nền đất (T0 có giá trị từ 0.8s ~ 1s). Trên cơ sở hệ số động đất này, tác động động đất lớn nhất lên hệ kết cấu được xác định theo biểu thức [4]: 0,max s 0,max sMF Q K Q g      x x (1.4) Tuy nhiên Mononobe đã bỏ qua lực cản, chưa xét đến lực động đất sẽ tăng thêm khi có tác dụng của dao động tự do và dao động cưỡng bức, phạm vi áp dụng cũng cho kết cấu có 1 bậc tự do nhưng chưa giải quyết được sự phân bố động đất theo chiều cao công trình (tức kết cấu có nhiều bậc tự do). Năm 1927, nhà khoa học Nga Zavriev đã đưa ra các yếu tố quan trọng trong dao động tự nhiên trong giai đoạn khởi đầu của tác động động đất [4]. Zavriev đã đặt nền móng đầu tiên cho cơ sở lý thuyết động lực học trong tính toán tác động động đất. Năm 1934, nhà khoa học Mỹ Biot đã đề xuất phương pháp tính tải trọng động đất bằng cách dùng các số liệu dao động nền đất thực ghi lại được khi động đất xảy ra. Năm 1949, Housner và Kahn đã đưa ra được cách xác định phổ gia tốc bằng thiết bị tương tự điện. Vào những thập niên 80 của thế kỷ XX, hàng loạt các kết quả nghiên cứu lý thuyết và thực nghiệm được thực hiện, quan điểm thiết kế kháng chấn mới được hình thành. Theo các quan điểm mới này, các công trình được thiết kế sao cho có khả năng chịu được các trận động đất vừa và nhỏ xuất hiện ngẫu nhiên mà công trình không bị hư hỏng, khi gặp các trận động đất mạnh thì công trình không bị sụp đổ. 1.2. Một số phương pháp tính toán 1.2.1. Phương pháp tính toán tĩnh tương đương Phương pháp tính toán tĩnh tương đương (còn gọi là phương pháp lực ngang tương đương) là phương pháp tính toán đơn giản nhất trong số các phương pháp được dùng để xác định phản ứng của kết cấu chịu tác động động đất. Phương pháp này giả định rằng kết cấu làm việc đàn hồi tuyến tính, còn tính phi tuyến hình học được xem xét tới một cách gián tiếp. Các tải trọng ngang tác động lên chiều cao 6 công trình được xem là tương đương với tác động động đất và được tổ hợp với các tải trọng đứng (lực trọng trường). Phương pháp này thường được sử dụng để thiết kế các công trình tương đối đều đặn có chu kỳ cơ bản bằng khoảng 1.5 - 2s. Đối với các công trình có hình dạng không đều đặn hoặc có chu kỳ dài cần sử dụng các phương pháp động chính xác hơn như phân tích dạng hoặc phân tích lịch sử phản ứng không đàn hồi. 1.2.2. Phương pháp tính toán tĩnh phi tuyến Trong phương pháp này, sự phân bố giả định lực quán tính ngang được dựa trên giả thiết cho rằng phản ứng của công trình được kiểm soát bởi một dạng dao động duy nhất và hình dạng của dao động này giữ nguyên không đổi trong suốt thời gian phản ứng. Thông thường, dạng dao động cơ bản được chọn là dạng phản ứng trội của hệ nhiều bậc tự do động, ảnh hưởng của các dạng dao động khác được xem là nhỏ và được bỏ qua. Phương pháp tính toán tĩnh phi tuyến với phân bố tải trọng ngang như vậy được gọi là phương pháp tính toán đẩy dần quy ước và thường được dùng để tính toán phản ứng của các công trình có chiều cao thấp và trung bình. Do tính đơn giản và khả năng xác định với độ chính xác chấp nhận được quá trình biến dạng của hệ kết cấu cũng như các bộ phận thành phần mà không cần phải thực hiện việc mô hình hóa phức tạp và tính toán công phu như tính toán động nên phương pháp tính toán đẩy dần được xem là một phương pháp hữu hiệu và tiện lợi trong tính toán động. 1.2.3. Phương pháp phân tích dạng dao động và phổ phản ứng Phản ứng của kết cấu có nhiều bậc tự do chịu tác động động đất có thể được tính toán bằng cách phân tích hệ kết cấu thành nhiều hệ kết cấu có một bậc tự do tương đương. Tính toán phản ứng mỗi hệ tương đương theo thời gian và sau đó cộng đại số các phản ứng lại để được phản ứng của kết cấu ban đầu. Phương pháp này được gọi là phương pháp phân tích dạng. Nếu việc tính toán chỉ nhằm xác định các đại lượng phản ứng lớn nhất thì tác động động đất sẽ được cho dưới dạng phổ phản ứng và kết quả tính toán theo phương pháp tích phân dạng dao động sẽ là phản ứng lớn nhất của hệ kết cấu. Phương pháp tính toán này có tên gọi là phương pháp phổ phản ứng. Phương pháp tích phân dạng dao động cũng như phương pháp phổ phản ứng có những nhược điểm sau:  Phụ thuộc vào việc tách một cách nhân tạo các dạng dao động.  Phải tổ hợp các kết quả tính toán ở các dạng dao động lại theo nguyên tắc cộng tác dụng nên chỉ giới hạn ở giai đoạn làm việc đàn hồi tuyến tính của vật liệu. 7  Không áp dụng được cho một số hệ kết cấu không sử dụng được kỹ thuật phân tích dạng.  Không cho các chỉ dẫn chính xác về sự hình thành khớp dẻo ở một số cấu kiện. 1.2.4. Phương pháp tích phân trực tiếp phương trình chuyển động Phản ứng của các hệ kết cấu chịu tác động bất kỳ hoặc động đất có thể xác định được bằng cách tích phân trực tiếp phương trình chuyển động theo thời gian. Phương pháp này không cần thay đổi hoặc biến đổi các phương trình chuyển động sang hệ có một hoặc nhiều bậc tự do như ở phương pháp tích phân dạng dao động. Phương pháp tích phân trực tiếp theo thời gian xác định các giá trị gần đúng của nghiệm đối với một tập hợp các giá trị thời gian T được lựa chọn. Có thể tóm tắt nguyên tắc của phương pháp này như sau: (i) giả thiết các hàm mô tả sự biến thiên của chuyển vị, vận tốc và gia tốc trong một khoảng thời gian và (ii) các phương trình chuyển động không phải thỏa mãn ở tất cả ở mọi thời gian T mà chỉ trong khoảng thời gian không đổi Δt. Khoảng thời gian này được gọi là bước thời gian. Điều này cũng có nghĩa rằng điều kiện cân bằng tĩnh của các lực quán tính, lực cản và lực đàn hồi với tải trọng tác động sẽ xảy ra ở nhiều bước thời gian Δt, 2Δt, , nΔt, . Ở mỗi bước thời gian, phương trình chuyển động được giải với các điều kiện ban đầu là chuyển vị, vận tốc được xác định ở bước trước đó. Phương pháp tích phân trực tiếp theo thời gian có thể áp dụng cho các hệ kết cấu tuyến tính lẫn phi tuyến nên có thể xem là phương pháp tổng quát duy nhất tính toán phản ứng động của các hệ kết cấu chịu tải trọng bất kỳ. 1.3. Kết luận chương 1 Chương 1 đã trình bày tóm lược lịch sử phát triển các phương pháp phân tích ứng xử của kết cấu chịu tải trọng động đất. Đồng thời cũng đưa ra một số ưu nhược điểm của các phương pháp đó. 8 Chương 2 XÂY DỰNG MÔ HÌNH PHẦN TỬ HỮU HẠN ĐỂ TÍNH TOÁN ĐÁP ỨNG ĐỘNG LỰC HỌC CỦA KHUNG, DẦM 2D-FGM CHỊU TẢI TRỌNG ĐỘNG ĐẤT 2.1. Dầm 2D-FGM Hình 2.1 minh họa dầm 2D-FGM với chiều dài L, chiều rộng b và chiều cao h trong hệ tọa độ Đề-các (x, z). Hệ tọa độ (x, z) được chọn sao cho trục x trùng với mặt giữa của dầm, trục z vuông góc với mặt giữa và hướng lên trên. Hình 2.1. Mô hình dầm 2D-FGM Dầm được giả định được tạo thành từ bốn vật liệu thành phần, cụ thể là gốm 1, gốm 2, kim loại 1 và kim loại 2. Tỉ lệ thể tích của các vật liệu thành phần phân bố theo quy luật hàm lũy thừa như sau: 1 2 1 1 1 , 2 2 z x z xn n n n c c z x z x V V h L h L                                 1 2 1 1 1 1 , 1 2 2 z x z xn n n n m m z x z x V V h L h L                                                   (2.1) trong đó Vc1, Vc2, Vm1, Vm2 lần lượt là tỉ phần thể tích của vật liệu gốm 1, gốm 2, kim loại 1 và kim loại 2; L và h tương ứng là chiều dài và chiều cao của dầm; nz và nx là các tham số vật liệu. Hình 2.2 minh họa tỉ phần thể tích vật liệu gốm 1 và gốm 2 khi nz và nx thay đổi. Các tính chất hữu hiệu (P) của dầm (mô-đun đàn hồi, mật độ khối, ) có thể được đánh giá theo mô hình Voigt: P=Vc1Pc1 + Vc2Pc2 + Vm1Pm1 + Vm2Pm2 (2.2) trong đó Pc1, Pc2, Pm1, Pm2 biểu thị các tính chất của vật liệu gốm 1, gốm 2, kim loại 1, kim loại 2. Thay phương trình (2.1) vào phương trình (2.2) ta được: 0 9 Hình 2.2. Tỉ phần thể Vc1 và Vc2 khi nz, nx thay đổi 1 1 2 2 1 1 ( , ) ( ) 1 ( ) 2 2 z x z xn n n n c m c m z x z x P x z P P P P h L h L                                                   (2.3) Từ phương trình (2.3), thành phần vật liệu chỉ là gốm 2 và kim loại 2 nếu 0xn  . Khi đó thay phương trình (2.1) vào (2.3) ta thu được: 2 2 2 1 ( , ) ( ) 2 zn c m m z P x z P P P h          (2.4) Còn khi nz = 0, thành phần vật liệu gồm có pha gốm 1 và gốm 2. 2.2. Các phương trình cơ bản Xét một phần tử dầm có chiều dài l, trong hệ tọa độ (x, z). Trục x được chọn trùng mặt giữa của dầm. Dựa trên lý thuyết dầm Timoshenko, chuyển vị dọc trục và chuyển vị ngang của một điểm bất kì trên dầm được cho bởi: 0 0 ( , , ) ( , ) z ( , t) ( ,z, ) ( , ) u x z t u x t x w x t w x t    (2.5) trong đó t là biến thời gian, 0 0( , ), ( , )u x t w x t tương ứng là chuyển vị dọc trục và chuyển vị theo phương ngang của một điểm bất kì trên mặt giữa của dầm; ( , )x z 0 0.5 1 -0.5 0 0.5 0 0.5 1 x/L (a) n z =n x =1/3 z/h V c 1 0 0.5 1 -0.5 0 0.5 0 0.5 1 x/L (b) n z =n x =1/3 z/h V c 2 0 0.5 1 -0.5 0 0.5 0 0.5 1 x/L (c) n z =n x =3 z/h V c 1 0 0.5 1 -0.5 0 0.5 0 0.5 1 x/L (d) n z =n x =3 z/h V c 2 10 là góc quay của thiết diện ngang của dầm. Biến dạng dọc trục ( xx ) và biến dạng trượt ( xz ) thu được từ phương trình (2.5) có dạng: 0 0 ( , ) ( , ) ( , ) , ( , )xx xz u x t x t w x t z x t x x x               (2.6) Theo định luật Hook, ứng suất dọc trục và ứng suất trượt được xác định: ( , ) , ( , )xx xx xz xzE x z G x z      (2.7) trong đó ( , )E x z và ( , )G x z tương ứng là mô-đun đàn hồi và mô-đun trượt hữu hiệu của dầm, ψ là hệ số hiệu chỉnh trượt và chọn bằng 5/6 cho mặt cắt hình chữ nhật. Năng lượng biến dạng đàn hồi (U) và động năng (T) có dạng: 2 2 2 11 , 12 , , 22 , , 0 2 2 2 11 12 22 0 1 2 33( ) 2 1 ( ) 2 2 L x x x x x L U A u A u A A w dx T I u w I u I dx                      (2.8) (2.9) Trong đó Aij là các độ cứng; Iij là các mô-men khối lượng và chúng được định nghĩa trong phương trình theo công thức: 2 11 12 22 33 ( , , ) ( , )(1, , ) , ( , ) A A A A A E x z z z dA A G x z dA     (2.10) và 2 11 12 22( , , ) ( , )(1, , ) A I I I x z z z dA  (2.11) Thay biểu thức các tính chất hiệu dụng ở phương trình (2.3) vào phương trình (2.10), ta có thể viết lại các độ cứng Aij theo dạng 1 1 1 1 2 2 11 11 11 11 1 1 1 1 2 2 12 12 12 12 1 1 1 1 2 2 22 22 22 22 1 1 1 1 2 2 33 33 33 33 ( ) , ( ) , ( ) , ( ) x x x x n c m c m c m n c m c m c m n c m c m c m n c m c m c m x A A A A L x A A A A L x A A A A L x A A A A L                                     (2.12) 11 trong đó 1 1 1 1 1 1 11 12 22, , c m c m c mA A A và 1 133 c mA là các độ cứng sinh ra bởi cặp vật liệu gốm 1 và kim loại 1; 2 2 2 2 2 2 11 12 22, , c m c m c mA A A và 2 233 c mA là các độ cứng được sinh ra bởi cặp vật liệu gốm 2 và kim loại 2. Từ phương (2.12) có thể nhận thấy rằng tham số vật liệu nx chỉ ảnh hướng tới số hạng thứ hai trong biểu thức của Aij và trong trường hợp 1 1 2 2 12 12 0 c m c mA A  thì độ cứng sẽ tăng khi tham số nx tăng. Cũng từ phương trình (2.12), độ cứng phần tử dầm 2D-FGM sẽ biến đổi về dầm FGM thông thường khi ta cho tham số nx = 0 hoặc bằng cách lựa chọn 2 vật liệu gốm và 2 kim loại giống hệt nhau. Với tính chất hiệu dụng của dầm được đưa ra ở phương trình (2.4), chúng ta có thể dễ dàng biểu diễn các độ cứng của dầm FGM một chiều ở dạng hiển. Chẳng hạn như độ cứng 1 1c m ijA (i, j = 1, 2) sẽ có biểu thức dạng hiển như sau: 1 1 1 1 11 2 1 1 1 1 12 3 2 1 1 1 1 22 1 1 1 1 33 ( ) ; ( 1) ( ) ; 2( 1)( 2) ( 2)( ) ; 4( 1)( 2)( 3) ( ) . ( 1) c m c z m z c m z c m z z c m z z c m z z z c m c z m z bh E n E A n bh n E E A n n bh n n E E A n n n bh G n G A n                  (2.13) Tương tự như các độ cứng 1 ij 1c mA và 2ij 2c mA , các mô-men khối lượng cũng có thể được viết lại như sau: 1 1 1 1 2 2 11 11 11 11 1 1 1 1 2 2 12 12 12 12 1 1 1 1 2 2 22 22 22 22 ( ) ( ) ( ) x x x n c m c m c m n c m c m c m n c m c m c m x I I I I L x I I I I L x I I I I L                            (2.14) trong đó 1 ij 1c mI và 2ij 2c mI tương ứng là các mô-men khối lượng sinh ra bởi cặp vật liệu gốm 1, kim loại 1 và gốm 2, kim loại 2. Các biểu thức hiển của 1 ij 1c mI và 2ij 2c mI cũng có dạng tương tự như (2.13). 2.3. Chuyển vị nút và nội suy Phần tử dầm 2 nút, mỗi nút gồm 3 bậc tự do với chiều dài l. Véc-tơ chuyển vị nút cho một phần tử khởi tạo (i, j) bao gồm các thành phần 12  i i i j j ju w u w d T (2.15) trong đó chỉ số trên ‘T ’ được sử dụng để chỉ chuyển vị của một véc-tơ hay một ma trận; ,i iu w và i tương ứng là chuyển vị dọc trục, chuyển vị theo phương ngang và góc xoay tại nút thứ i; ,j ju w và i tương ứng là chuyển vị dọc trục, chuyển vị theo phương ngang và góc xoay tại nút thứ j. Chuyển vị dọc trục u(x), chuyển vị theo phương ngang w(x) và góc xoay θ(x) cho phần tử của dầm được nội suy qua hàm dạng như sau 0 0, ,u wu w   N d N d N d (2.16) trong đó       1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6 u u u u u u u w w w w w w w N N N N N N N N N N N N N N N N N N          N N N (2.17) tương ứng là ma trận các hàm nội suy (hàm dạng) cho các chuyển vị dọc trục 0 ( , )u x t , chuyển vị theo phương ngang 0 ( , )w x t và góc xoay ( , )x t . Các hàm nội suy tuyến tính được dùng để nội suy chuyển vị dọc trục 0 ( , )u x t và hàm nội suy Kosmatka được sử dụng cho chuyển vị theo phương ngang 0 ( , )w x t và ( , )x t . Biểu thức của các hàm nội suy như sau: 1 4 2 3 5 6 ; ; 0 u u u u u u l x x N N l l N N N N        (2.18) 1 4 3 2 2 3 2 3 3 2 5 3 6 0; 1 2 3 1 1 2 1 1 2 2 1 2 3 1 1 1 2 w w w w w w N N x x x N l l l l x x x N l l l x x x N l l l l x N l                                                                                                                 2 2 x x l l                 (2.19) 13 2 1 4 2 2 3 2 5 2 6 6 0; ; (1 ) 1 3 (4 ) 1 ; 1 6 ; (1 ) 1 3 (2 ) . 1 x x N N N l l l x x N l l x x N l l l x x N l l                                                                                             (2.20) trong các phương trình (2.18), (2.19) và (2.20), 2 1 112 /m mE I l G A  là tham số trượt, A và I lần lượt là diện tích, mô-men quán tính của thiết diện ngang. 2.4. Ma trận độ cứng Sử dụng phép nội suy trên ta có thể viết được biểu thức cho năng lượng biến dạng cho một phần tử dầm Ue, dưới dạng sau đây: , , , , , , , , 2 2 2 11 , 12 , , 22 , 33 , 0 11 12 22 0 33 1 2 ( ) 2 1 [ 2 2 ( ) ( )] 1 ( + ) 2 1 = 2 x x x x x x x x l e x x x x x l T T T T u u u T w w T uu u T U A u A u A A w dx A A A A dx                               d N N N N N N N N N N d d k k k k d d kd (2.21) trong đó: +uu u    k k k k k (2.22) là ma trận độ cứng phần tử và , , , , , , , , 11 12 0 0 22 33 0 0 ; ; ; ( ) ( ) x x x x x x x x l l T T uu u u u u l l T T w w A dx A dx A dx A dx                    k N N k N N k N N k N N N N (2.23) tương ứng là ma trận độ cứng phần tử sinh ra từ biến dạng dọc trục, tương hỗ giữa biến dạng dọc trục và uốn, biến dạng uốn và biến dạng trượt. 14 Trong một trường hợp tổng quát, khi phần tử nghiêng so với trục ngang của hệ tọa độ tổng quát một góc α, các ma trận độ cứng phần tử được biểu diễn: T[ ] [ ] [ ][ ]g k S k S (2.24) trong đó [kg] là ma trận độ cứng phần tử trong hệ tọa độ tổng quát, [S] là ma trận quay, có dạng sau cos sin 0 0 0 0 sin cos 0 0 0 0 0 0 1 0 0 0 0 0 0 cos sin 0 0 0 0 sin cos 0 0 0 0 0 0 1                            S (2.25) Ma trận độ cứng tổng thể được xác định thông qua các ma trận độ cứng phần tử như sau: 1 [ ] nELE g e e  K k (2.26) 2.5. Ma trận khối lượng Có nhiều loại ma trận khối lượng khác nhau có thể sử dụng trong phân tích động lực học kết cấu. Tuy nhiên trong nghiên cứu này, tác giả sử dụng ma trận khối lượng nhất quán, là ma trận xây dựng dựa trên các hàm nội suy cho trường chuyển vị (2.15). Với trường chuyển vị này, ta có thể viết biểu thức động năng cho phần tử dầm dưới dạng sau: 2 2 11 12 22 0 11 11 12 22 0 1 [ ( + ) 2 ] 2 1 [ 2 )] 2 1 ( ) 2 1 2 l e l T T T T T u u w w u T uu ww u T T I u w I u I dx I I I I dx                   d N N N N N N N N d d m m m m d d md (2.27) trong đó uu ww u    m m m m m (2.28) là mà trận khối lượng nhất quán của phần tử và 15 11 11 0 0 12 22 0 0 , , l l T T uu u u ww w w l l T T u u I dx I dx I dx I dx              m N N m N N m N N m N N (2.29) tương ứng là các ma trận nhất quán sinh ra từ chuyển vị dọc trục, chuyển vị theo phương ngang, tương hỗ giữa chuyển vị dọc trục và sự quay của thiết diện ngang, sự quay của thiết diện ngang. Khi phần tử nằm nghiêng so với trục ngang của hệ tọa độ tổng quát một góc α, ma trận khối lượng phần tử có dạng T[ ] [ ] [ ][ ]g m S m S (2.30) trong đó [mg] là ma trận khối lượng phần tử trong hệ tọa độ tổng quát, [S] là ma trận quay được định nghĩa trong phương trình (2.25). Từ phương trình (2.30), ma trận khối lượng tổng thể được xác định thông qua việc nối ghép: 1 [ ] nELE g e e  M m (2.31) 2.6. Phương pháp tích phân trực tiếp Phương trình chuyển động cho phân tích kết cấu chịu tải trọng động đất theo phương pháp phần tử hữu hạn có thể được viết dưới dạng [6]: gMD + CD + KD = -MID (2.32) trong đó: D là véc-tơ chuyển vị nút tổng thể, / t  D D là véc-tơ vận tốc nút tổng thể, 2 2/ t  D D là véc-tơ gia tốc nút tổng thể, C là ma trận cản được hình thành từ tổ hợp tuyến tính của ma trận khối lượng tổng thể (2.31) và ma trận độ cứng tổng thể (2.26):   C M K (2.33) trong đó  và  tương ứng là hệ số cản tỉ lệ với khối lượng và độ cứng, chúng được tính từ tỉ số cản tới hạn và tần số tự nhiên của kết cấu như sau: 1 2 1 2 1 2 2 2 ,              (2.34) Một tỉ lệ cản 0.02  được giả định cho kết cấu 2D-FGM. Vế phải của (2.32) xác định ngoại lực do chuyển động nền, trong đó gD là véc-tơ chuyển động nền 16 và I là véc-tơ hệ số ảnh hưởng, có giá trị 1 cho các phần tử tương ứng với bậc tự do theo hướng chuyển động nền và có giá trị 0 cho các bậc tự do khác. Một bảng số liệu của gia tốc nền trong 20 giây đầu tiên của trận động đất El Centro xảy ra tại miền Nam California năm 1940 [6] được minh họa trong Hình 2.3. Phương pháp tích phân trực tiếp Newmark được sử dụng rộng rãi trong việc tính toán đáp ứng động của kết cấu. Trong đó xấp xỉ sai phân hữu hạn được dùng để thay thế cho các đạo hàm riêng ở (2.35), tức thay thế D và D bằng sai phân của chuyển vị nút D tại các thời điểm khác nhau. Ý tưởng trung tâm của phương pháp này là phân chia tổng thời gian ΔT thành các bước thời gian nhỏ Δt. Các đáp ứng động học của kết cấu được tính theo thời gian Δt, 2Δt, 3Δt, ... nΔt ... Các phương trình của chuyển động tại một thời điểm mới (n + 1)Δt là: 1 1 1 ( 1) n n n g n      MD CD KD MID (2.35) Đối với phân tích tuyến tính, ma trận độ cứng K là không thay đổi từ thời gian dừng kế tiếp. Tuy nhiên, để phân tích phi tuyến, K là một hàm của chuyển vị D. Có nhiều cách khác nhau có thể được sử dụng để tính toán các đáp ứng động học tại thời điểm (n + 1) Δt, trong đó họ phương pháp tích phân trực tiếp Newmark là rất phổ biến và được cho bởi [7] 2 1 1 1 1 [(1 2 ) 2 ] 2 [(1 ) n n n n n n n n n t t t                     D D D D D D D D D (2.36) trong đó β và γ là các hằng số do người phân tích lựa chọn để kiểm soát tính hội tụ và độ chính xác của thuật toán số. Bằng việc giải hệ phương trình (2.36), thu được: 1 1 1 12 ( ) 1 1 2 1 1 ( ) 1 2 n n n n n n n n n n t t t t                                            D D D D D D D D D D (2.37) Thay phương trình (2.37) vào phương trình (2.35), thu được: ef 1 g( 1) 2 1 1 1 1 2 1 1 2 n n n n n n n n t t t t                                                K D MID M D D D C D D D (2.38) 17 trong đó : ef 2 1 , t t         K M C K (2.39) Với các giá trị khác nhau của β và γ sẽ thu được các phương pháp tích phân trực tiếp khác nhau dùng trong phân tích động lực học kết cấu như:  Phương pháp vi phân trung tâm: 1 0, 2     Phương pháp gia tốc tuyến tính: 1 1 , 6 2     Phương pháp gia tốc trung bình: 1 1 , 4 2     Phương pháp Fox-Goodwin: 1 1 , 2 2    và trong Luận văn này, tác giả lựa chọn sử dụng phương pháp gia tốc trung bình vì phương pháp này ổn định không điều kiện. Phương trình (2.39) cùng với phương trình (2.37) hoàn toàn xác định chuyển vị nút, vận tốc và gia tốc tại thời gian mới ( 1)n t  . Lưu ý rằng các ma trận độ cứng hiệu dụng không thể là một ma trận đường chéo bởi vì nó có chứa các ma trận độ cứng K. Hình 2.3. Gia tốc nền ghi nhận được của trận động đất El Centro 0 2 4 6 8 10 12 14 16 18 20 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 Time (s) G ro u n d a c c e le ra ti o n ( g ) 18 Hình 2.4. Sơ đồ khối của thuật toán Bắt đầu Số liệu đầu vào i=0 i=i+1 i < nSTEP Hiển thị kết quả Kết thúc Đúng Sai 19 2.7. Kết luận chương 2 Chương 2 đã trình bày mô hình kết cấu khung, dầm 2D-FGM với các tính chất vật liệu thay đổi theo chiều dài và chiều cao của dầm bằng quy luật lũy thừa. Tác giả tiến hành xây dựng các phương trình cơ bản, biểu thức năng lượng cho dầm 2D-FGM. Các biểu thức năng lượng được đưa ra dựa trên trường chuyển vị nhận được từ lý thuyết biến dạng trượt bậc nhất. Các biểu thức cụ thể cho ma trận độ cứng và ma trận khối lượng cho phần tử dầm 2D-FGM được xây dựng từ các biểu thức năng lượng. Tiếp đến, Luận văn xây dựng mô hình phần tử hữu hạn để tính toán đáp ứng động lực học của khung, dầm 2D-FGM chịu tải trọng động đất dưới sự trợ giúp của thuật toán tích phân trực tiếp Newmark. 20 Chương 3 TÍNH TOÁN SỐ VÀ THẢO LUẬN Một chương trình tính dựa trên phần tử được mô tả và phương pháp gia tốc trung bình được phát triển và sử dụng để phân tích một số kết cấu khung, dầm 2D-FGM như trong hình 3.1. Số liệu của các vật liệu thành phần sử dụng trong Luận văn này được cho trong bảng 3.1. Các dầm và khung xem xét trong Luận văn được giả thiết chịu ảnh hưởng của gia tốc nền của trận động đất El Centro như trong hình 2.3. 3.1. Kiểm tra chương trình tính toán Để đảm bảo tính chính xác của các phần tử được phát triển cũng như chương trình tính toán số được xây dựng dựa trên phương pháp tích phân số Newmark, tác giả tiến hành kiểm chứng chương trình tính. Bằng cách cho tham số vật liệu nx = 0 và nz = n, phần tử dầm 2D-FGM quay trở về phần tử FGM thông thường. Hình 3.1. Kết cấu khung, dầm 2D-FGM được nghiên cứu Bảng 3.1. Tính chất vật liệu thành phần cho khung, dầm 2D-FGM Vật liệu Vai trò E (GPa) 3(kg/m ) v Steel Kim loại 1 210 7800 0.3 Aluminum Kim loại 2 70 2702 0.23 Alumina Gốm 1 390 3960 0.3 Zirconia Gốm 2 200 5700 0.3 21 Đầu tiên, một cột FGM, ngàm chặt một đầu và một đầu tự do như thể hiện trong hình 3.1(a) được tác giả phân tích. Chiều cao của cột là 20 m, kích thước của mặt cắt thiết diện ngang của nó là b = h = 0.2 m. Hai mươi phần tử được sử dụng để rời rạc hóa cột FGM. Bảng 3.2. So sánh tần số và phản ứng của cột thép Nguồn f (Hz) max(uL) (m) min(uL) (m) max(vL) (m/s) min(vL) (m/s) max(aL) (m/s2) min(aL) (m/s2) ANSYS 0.4245 0.4503 -0.4519 1.6133 -1.5354 12.7040 -13.0950 Luận văn 0.4245 0.4326 -0.4396 1.6086 -1.5498 13.9556 -14.4007 Hình 3.2. Chuyển vị ngang tương đối theo thời gian tại đỉnh cột Trong bảng 3.2, tần số cơ bản, các giá trị lớn nhất và nhỏ nhất của chuyển vị tương đối, vận tốc và gia tốc tuyệt đối ở điểm trên cùng của cột thép được so sánh với các kết quả thu được bằng cách sử dụng phần mềm ANSYS 15 để mô phỏng. Như đã thấy từ bảng 3.2, các kết quả số thu được trong việc tính toán phù hợp so với việc mô phỏng bằng phần mềm ANSYS. Hình 3.2 và hình 3.3 minh họa chuyển vị tương đối và vận tốc ở đầu cột FGM cho các giá trị khác nhau của các tham số vật liệu n. Trong hình 3.4, không gian pha cho chuyển vị và vận tốc ở đầu cột được biểu thị cho n = 0.5 và n = 5. Như đã thấy rõ từ các kết quả số thu được, tham số vật liệu n cao sẽ ảnh hưởng đến phản ứng địa chấn của cột. Cả chuyển vị theo phương ngang và vận tốc tại 0 2 4 6 8 10 12 14 16 18 20 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 Time (s) u ( m ) n = 0.5 n = 5 fure steel 22 đỉnh cột với tham số n thấp hơn là thấp hơn đáng kể so với cột có tham số vật n cao hơn. Kết quả số thu được cũng chỉ ra rằng, cột FGM có đáp ứng động lực học tốt hơn nhiều so với cột được làm từ vật liệu thép thuần nhất. Hình 3.3. Vận tốc theo thời gian tại đỉnh cột Hình 3.4: Không gian pha của chuyển vị ngang tương đối và vận tốc tại đỉnh cột Tiếp theo, tác giả đi phân tích một khung FGM giản đơn như trong hình 3.1(b). Phần khung được làm bằng ba dầm với chiều dài 5 m, kích thước của mặt cắt ngang là b = h = 0.25 m. Trong hình 3.5, 3.6 và 3.7, chuyển vị ngang tương đối, vận tốc và gia tốc theo thời gian ở góc trên bên trái, điểm A của khung được minh họa cho hai giá trị của tham 0 2 4 6 8 10 12 14 16 18 20 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Time (s) v ( m /s ) n = 0.5 n = 5 fure steel -0.3 -0.15 0 0.15 0.3 -1.5 -1 -0.5 0 0.5 1 1.5 u (m) v ( m /s ) -0.3 -0.15 0 0.15 0.3 -1.5 -1 -0.5 0 0.5 1 1.5 u (m) v ( m /s ) (a) n=0.5 (b) n=5 23 số vật liệu n = 0,5 và n = 5. Các giá trị lớn nhất của các đáp ứng động lực của khung hình 3.1(b) là nhỏ hơn khi tham số vật liệu n nhỏ hơn. Các kết quả thu được có thể được giải thích bằng sự gia tăng của thành phần gốm ứng với tham số vật liệu n nhỏ hơn. Điều này dẫn tới độ cứng của kết cấu tăng, dẫn tới các giá trị đáp ứng động lực học nhỏ hơn khi n nhỏ. Các kết quả thu được khi phân tích cột FGM và khung giản đơn FGM hoàn toàn trùng khớp với các kết quả trong [17]. Hình 3.5. Chuyển vị ngang tương đối theo thời gian tại đỉnh A của khung Hình 3.6. Vận tốc theo thời gian tại đỉnh A của khung 0 2 4 6 8 10 12 14 16 18 20 -0.03 0 0.03 u ( m ) Time (s) n = 0.5 n = 5 0 2 4 6 8 10 12 14 16 18 20 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 v ( m /s ) Time (s) n = 0.5 n = 5 24 Hình 3.7. Gia tốc theo thời gian tại đỉnh A của khung 3.2. Cột 2D-FGM Cột 2D-FGM được ngàm chặt một đầu, một đầu tự do trong hình 3.1(a). Cột có kích thước thiết diện ngang b = h = 0.2 m và chiều dài L=10 m. Mười phần tử đã được sử dụng để rời rạc hóa cột 2D-FGM. Hình 3.8. Chuyển vị ngang tương đối theo thời gian tại đỉnh cột (nz = 0.5) 0 2 4 6 8 10 12 14 16 18 20 -15 0 15 Time (s) a ( m /s 2 ) n = 0.5 n = 5 0 2 4 6 8 10 12 14 16 18 20 -0.1 -0.05 0 0.05 0.1 u ( m ) Time (s) n x = 0.2 n x = 3(nz = 0.5) 25 Hình 3.8, 3.9 và 3.10 thể hiện chuyển vị ngang tương đối, vận tốc và gia tốc của cột 2D-FGM ứng với trường hợp tham số nz cố định. Hình 3.11, 3.12 và 3.13 biểu thị chuyển vị ngang tương đối, vận tốc và gia tốc của cột 2D-FGM ứng với trường hợp tham số nx cố định. Như đã thấy rõ từ các kết quả số thu được, hai tham số vật liệu nz và nx có ảnh hưởng lớn đến phản ứng địa chấn của cột 2D-FGM. Cụ thể, khi nz cố định, các giá trị của các đáp ứng động lực học ở đỉnh cột ứng với tham số nx lớn hơn là thấp hơn. Còn khi tham số vật liệu nx cố định, cả chuyển vị Hình 3.9. Vận tốc theo thời gian tại đỉnh cột (nz=0.5) Hình 3.10. Gia tốc theo thời gian tại đỉnh cột (nz=0.5) 0 2 4 6 8 10 12 14 16 18 20 -1.5 -1 -0.5 0 0.5 1 1.5 v ( m /s ) Time (s) n x = 0.2 n x = 3(nz = 0.5) 0 2 4 6 8 10 12 14 16 18 20 -20 -15 -10 -5 0 5 10 15 20 a ( m /s 2 ) Time (s) n x = 0.2 n x = 3(nz = 0.5) 26 Hình 3.11. Chuyển vị ngang tương đối theo thời gian tại đỉnh cột (nx=0.5). Hình 3.12. Vận tốc theo thời gian tại đỉnh cột (nx=0.5) ngang tương đối, vận tốc và gia tốc ở đỉnh cột có xu hướng tăng khi tham số nz tăng. Các kết quả thu được có thể giải thích bởi sự gia tăng của hàm lượng vật liệu kim loại khi tham số vật liệu nz tăng lên. Điều đó dẫn đến các độ cứng được định nghĩa trong (2.10) giảm đi. Do vậy các kết quả đáp ứng động lực học của cột sẽ giảm khi tham số nz tăng lên. Lập luận tương tự được sử dụng để giải thích ảnh 0 2 4 6 8 10 12 14 16 18 20 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 u ( m ) Time (s) n z = 0.2 n z = 3(nx = 0.5) 0 2 4 6 8 10 12 14 16 18 20 -1.5 -1 -0.5 0 0.5 1 1.5 v ( m /s ) Time (s) n z = 0.2 n z = 3(nx = 0.5) 27 hưởng của tham số vật liệu nx đối với các giá trị đáp ứng động lực học khi tham số nz cố định. Hình 3.13. Gia tốc theo thời gian tại đỉnh cột (nx=0.5) 3.3. Khung giản đơn Một khung giản đơn 2D-FGM trong hình 3.1(b) được phân tích. Khung được cấu tạo từ ba phần tử dầm 2D-FGM với chiều dài 5 m, kích thước của mặt cắt ngang là b = h = 0.25 m. Hình 3.14. Chuyển vị ngang tương đối theo thời gian tại đỉnh A của khung (nz=0.5) 0 2 4 6 8 10 12 14 16 18 20 -1.5 -1 -0.5 0 0.5 1 1.5 v ( m /s ) Time (s) n z = 0.2 n z = 3(nx = 0.5) 0 2 4 6 8 10 12 14 16 18 20 -4 -3 -2 -1 0 1 2 3 4 x 10 -3 u ( m ) Time (s) n x = 0.2 n x = 3(nz = 0.5) 28 Ảnh hưởng của các tham số vật liệu đối với các giá trị đáp ứng động lực học của khung giản đơn được minh họa trong các hình 3.14 đến hình 3.19, trong đó chuyển vị ngang tương đối, vận tốc và gia tốc theo thời gian tại đỉnh A của khung được nghiên cứu với các giá trị khác nhau của tham số vật liệu nz và nx. Với tham số vật liệu nz cố định, hình 3.14, 3.15 và 3.16 cho thấy sự giảm của chuyển vị ngang, vận tốc và gia tốc ở đỉnh A của khung khi giá trị tham số vật liệu nx tăng. Mặt khác, hình 3.17, 3.18 và 3.19 cho thấy sự tăng của chuyển vị ngang tương đối, vận tốc và gia tốc khi tham số nx cho trước và nz tăng. Sự giảm của các giá trị đáp ứng Hình 3.15. Vận tốc theo thời gian tại đỉnh A của khung (nz=0.5) Hình 3.16. Gia tốc theo thời gian tại đỉnh A của khung (nz=0.5) 0 2 4 6 8 10 12 14 16 18 20 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 v ( m /s ) Time (s) n x = 0.2 n x = 3(nz = 0.5) 0 2 4 6 8 10 12 14 16 18 20 -10 -5 0 5 10 a ( m /s 2 ) Time (s) n x = 0.2 n x = 3(nz = 0.5) 29 Hình 3.17. Chuyển vị ngang tương đối theo thời gian tại đỉnh A của khung (nx=0.5) Hình 3.18. Vận tốc theo thời gian tại đỉnh A của khung (nx=0.5) động lực học khi tăng nx có thể được giải thích bằng sự gia tăng hàm lượng gốm 1, như có thể nhìn thấy ở công thức (2.1). Do mô-đun đàn hồi của vật liệu gốm 1 như trong bảng 1, cao hơn nhiều so với mô-đun đàn hồi của vật liệu gốm 2. Kết quả là độ cứng được định nghĩa theo phương trình (2.10) cao hơn trong trường hợp khung kết hợp với tham số vật liệu nx cao. Do vậy đáp ứng động lực học thấp hơn là kết quả của khung có tham số vật liệu nx cao hơn. Lập luận tương tự có thể 0 2 4 6 8 10 12 14 16 18 20 -4 -3 -2 -1 0 1 2 3 4 x 10 -3 u ( m ) Time (s) n z = 0.2 n z = 3(nx = 0.5) 0 2 4 6 8 10 12 14 16 18 20 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 v ( m /s ) Time (s) n z = 0.2 n z = 3(nx = 0.5) 30 được sử dụng để giải thích ảnh hưởng của tham số nz đối với các đáp ứng động lực học của khung khi cố định tham số nx. Hình 3.19. Gia tốc theo thời gian tại đỉnh A của khung (nx=0.5). 3.4. Khung nhiều tầng Khung 2D-FGM nhiều tầng thể hiện trong hình 3.1(c) được xem xét. Phần khung được hình thành từ mười hai dầm có cùng chiều dài và kích thước mặt cắt ngang, L = 5 m, b = h = 0.25 m. Hình 3.20. Chuyển vị ngang tương đối theo thời gian tại đỉnh B của khung (nz=0.5) 0 2 4 6 8 10 12 14 16 18 20 -10 -5 0 5 10 a ( m /s 2 ) Time (s) n z = 0.2 n z = 3(nx = 0.5) 0 2 4 6 8 10 12 14 16 18 20 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 Time (s) u ( m ) n x = 0.2 n x = 3(nz = 0.5) 31 Hình 3.21. Vận tốc theo thời gian tại đỉnh B của khung (nz=0.5) Hình 3.22. Gia tốc theo thời gian tại đỉnh B của khung (nz=0.5) Hình 3.20, 3.21 và 3.22 thể hiện chuyển vị tương đối theo phương ngang, vận tốc và gia tốc theo thời gian tại đỉnh B, góc trên bên trái của khung. Các giá trị đáp ứng động lực học của khung đa tầng là khác so với cột và khung giản đơn 2D- FGM. Với giá trị tham số nz cho trước, vận tốc và gia tốc tại đỉnh B của khung nhiều tầng thậm chí còn tăng khi tham số nx tăng lên. Ảnh hưởng của tham số vật liệu nz đến các đáp ứng động lực học của khung nhiều tầng trong trường hợp tham 0 2 4 6 8 10 12 14 16 18 20 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 v ( m /s ) Time (s) n x = 0.2 n x = 3(nz = 0.5) 0 2 4 6 8 10 12 14 16 18 20 -10 -5 0 5 10 a ( m /s 2 ) Time (s) n x = 0.2 n x = 3(nz = 0.5) 32 số nx cho trước cũng khác so với khung giản đơn. Biên độ lớn nhất của chuyển vị theo phương ngang và vận tốc tại đỉnh B của khung không có sự sai khác nhiều giữa hai trường hợp tham số nz=0.2 và nz=3. Do đó, có thể thấy rằng phản ứng địa chấn của khung dầm 2D-FGM không chỉ phụ thuộc vào tham số vật liệu nz và nx mà còn phụ thuộc vào cấu hình thực của kết cấu. Hình 3.23. Chuyển vị ngang tương đối theo thời gian tại đỉnh B của khung (nx=0.5) Hình 3.24. Vận tốc theo thời gian tại đỉnh B của khung (nx=0.5) 0 2 4 6 8 10 12 14 16 18 20 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 u ( m ) Time (s) n z = 0.2 n z = 3(nx = 0.5) 0 2 4 6 8 10 12 14 16 18 20 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 v ( m /s ) Time (s) n z = 0.2 n z = 3(nx = 0.5) 33 Hình 3.25. Gia tốc theo thời gian tại đỉnh B của khung (nx=0.5) 3.5. Khung bất đối xứng Cuối cùng, một khung bất đối xứng mô tả trong hình 3.1(d) được xem xét. Phần khung được hình thành từ hai cột với chiều dài của cột dọc là L=20 m, kích thước thiết diện ngang b = h = 0.25 m và khoảng cách giữa hai chân cột là L/4 = 5 m. Mười phần tử, năm cho mỗi cột được sử dụng trong phân tích khung. Hình 3.26. Chuyển vị ngang tương đối theo thời gian tại đỉnh C của khung (nz=0.5) 0 2 4 6 8 10 12 14 16 18 20 -10 -5 0 5 10 a ( m /s 2 ) Time (s) n z = 0.2 n z = 3(nx = 0.5) 0 2 4 6 8 10 12 14 16 18 20 -4 -2 0 2 4 6 x 10 -4 u ( m ) Time (s) n x = 0.2 n x = 3(nz = 0.5) 34 Hình 3.26 đến 3.31 mô tả chuyển vị theo phương ngang, vận tốc và gia tốc của khung bất đối xứng. Có thể thấy rằng, trong trường hợp với tham số nz cho trước, sự giảm của các đáp ứng động lực học tại đỉnh C của khung bất đối xứng khi tham số nx tăng lên. Ngược lại, đối với trường hợp tham số nx cho trước, dễ nhận thấy sự tăng của các đáp ứng động lực học tại đỉnh C của khung khi tham số nz tăng lên. Đáng ngạc nhiên, các đáp ứng động lực học của khung bất đối xứng là tốt hơn nhiều so với khung nhiều tầng, trong khi nó chỉ được hình thành từ hai cột. Biên Hình 3.27. Vận tốc theo thời gian tại đỉnh C của khung (nz=0.5) Hình 3.28. Gia tốc theo thời gian tại đỉnh C của khung (nz=0.5) 0 2 4 6 8 10 12 14 16 18 20 -0.01 -0.005 0 0.005 0.01 v ( m /s ) Time (s) n x = 0.2 n x = 3(nz = 0.5) 0 2 4 6 8 10 12 14 16 18 20 -1 -0.5 0 0.5 1 a ( m /s 2 ) Time (s) n x = 0.2 n x = 3(nz = 0.5) 35 độ của sự dịch chuyển, vận tốc và gia tốc ở phía trên, đỉnh C của khung bất đối xứng là thấp hơn nhiều so với khung 2D-FGM nhiều tầng. Các số liệu tính toán một lần nữa cho thấy rằng đáp ứng động lực học của khung 2D-FGM không chỉ phụ thuộc vào các tham số vật liệu nz, nx mà còn phụ thuộc nhiều vào cấu hình thực của khung. Hình 3.29. Chuyển vị ngang tương đối theo thời gian tại đỉnh C của khung (nx=0.5) Hình 3.30. Vận tốc theo thời gian tại đỉnh C của khung (nx=0.5) 0 2 4 6 8 10 12 14 16 18 20 -6 -4 -2 0 2 4 6 x 10 -4 u ( m ) Time (s) n z = 0.2 n z = 3(nx = 0.5) 0 2 4 6 8 10 12 14 16 18 20 -0.01 -0.005 0 0.005 0.01 v ( m /s ) Time (s) n z = 0.2 n z = 3(nx = 0.5) 36 Hình 3.31. Gia tốc theo thời gian tại đỉnh C của khung (nx=0.5) 3.6. Kết luận chương 3 Chương 3 đã tiến hành thiết lập chương trình tính toán số để phân tích một số kết cấu khung, dầm 2D-FGM với các cấu hình khác nhau. Việc phân tích này dựa trên cơ sở các phần tử được xây dựng và thuật toán số trong chương 2. Xem xét ảnh hưởng của hai tham số vật liệu tới ứng xử động lực học của kết cấu dưới tác động của tải trọng động đất đã được khảo sát chi tiết và thảo luận trong chương này. 0 2 4 6 8 10 12 14 16 18 20 -1 -0.5 0 0.5 1 a ( m /s 2 ) Time (s) n z = 0.2 n z = 3(nx = 0.5) 37 KẾT LUẬN Luận văn trình bày phương pháp phần tử hữu hạn trong phân tích địa chấn của kết cấu khung, dầm 2D-FGM. Tính chất của vật liệu được giả định thay đổi theo cả chiều cao và chiều dài của dầm, theo quy luật hàm lũy thừa. Phần tử dầm 2 nút, mỗi nút 3 bậc tự do, sử dụng các hàm dạng Kosmatka để nội suy chuyển vị theo phương ngang và góc xoay của thiết diện ngang dùng trong phân tích được xây dựng trong Luận văn. Ma trận độ cứng và ma trận khối lượng được thiết lập từ các biểu thức năng lượng biến dạng đàn hồi và động năng của phần tử. Đáp ứng động lực học của kết cấu dưới tác động của trận động đất El Centro được tính toán với sự trợ giúp của phương pháp tích phân trực tiếp Newmark. Khung 2D- FGM với các dạng hình học khác nhau đã được phân tích và ảnh hưởng của tham số vật liệu đối với ứng xử động lực học của khung đã được tính toán và thảo luận. Kết quả phân tích số nhận được trong Luận văn có thể tóm lược dưới đây: 1) Phần tử dầm 2D-FGM và thuật toán số xây dựng trong Luận văn đủ tin cậy và hiệu quả trong việc tính toán đáp ứng động lực học của khung, dầm 2D-FGM chịu tải trọng động đất. 2) Hai tham số vật liệu xác định sự phân bố của vật liệu theo chiều cao và chiều dài của dầm có ảnh hưởng khác nhau đến đáp ứng động lực học của khung, dầm 2D-FGM chịu tải trọng động đất. Chuyển vị ngang, vận tốc và gia tốc của kết cấu khung, dầm 2D-FGM không chỉ phụ thuộc vào hai tham số vật liệu nz và nx mà còn phụ thuộc nhiều vào cấu hình thực của kết cấu. Như đã nói trong phần mở đầu, nghiên cứu ứng xử động đất của kết cấu FGM nói chung và khung, dầm 2D-FGM thực hiện trong Luận văn nói riêng mới chỉ là những nghiên cứu ban đầu về ứng xử động đất của kết cấu làm từ vật liệu có cơ tính biến thiên. Để hiểu rõ hơn về ứng xử của kết cấu làm từ loại vật liệu mới này dưới tác động của tải trọng động đất, vì thế rất cần các nghiên cứu tiếp theo. 38 DANH MỤC CÔNG TRÌNH KHOA HỌC CỦA TÁC GIẢ LIÊN QUAN ĐẾN LUẬN VĂN 1. Nguyen Quang Huan, Bui Manh Cuong, and Nguyen Dinh Kien (2016), “Seismic Analysis of Planar Functionally Graded Beams and Frames Using Direct Integration Method”, Proceedings of the 4th International Conference on Engineering Mechanics and Automation (ICEMA4), Hanoi August, pp. 332-339. 2. Nguyen Quang Huan, Nguyen Dinh Kien (2017), “Finite Element Analysis of Planar 2D-FGM Beam and Frame Structures Excited by Earthquake Loads”, Tuyển tập Hội nghị Cơ học toàn quốc lần thứ X, Hà Nội, pp. 504- 511. 3. Dinh Kien Nguyen, Quang Huan Nguyen, Thi Thom Tran, Van Tuyen Bui (2017), “Vibration of bi-dimensional functionally graded Timoshenko beams excited by a moving load”, Acta Mechanica, Vol. 228, pp. 141–155. 39 TÀI LIỆU THAM KHẢO Tiếng Việt 1. Lê Thị Hà (2015), Phân tích dầm FGM có mặt cắt ngang thay đổi chịu tải trọng di động, Luận án Tiến sĩ Khoa học, Viện Cơ học, Viện Hàn lâm Khoa học và Công nghệ Việt Nam, Hà Nội. 2. Nguyễn Đình Kiên (2009), Phương pháp phần tử hữu hạn trong cơ học kết cấu và vật rắn biến dạng. Nhà xuất bản Khoa học tự nhiên và Công nghệ, Hà Nội. 3. Chu Quốc Thắng (1997), Phương pháp phần tử hữu hạn, Nhà xuất bản Khoa học và Kỹ thuật, Hà Nội. 4. Nguyễn Lê Linh (2011), Cơ sở lý thuyết tính toán công trình chịu động đất. Nhà xuất bản Khoa học và Kỹ thuật, Hà Nội. Tiếng Anh 5. ANSYS Inc. (2013), ANSYS 15.0 release. Pittsburgh, PA, USA. 6. A.K. Chopra (1995), Dynamics of Structures. Theories and Applications to Earthquake Engineering. Prentice-Hall Inc, Englewood. 7. R.D. Cook, D.S. Malkus, M.E., Plesha, R.J. Witt (2002), Concepts and applications of finite element analysis. John Wiley & Sons, Inc. USA, fourth edition. 8. B.S Gan, T.H. Trinh, T.H. Le, D.K. Nguyen (2015), “Dynamic response of non-uniform Timoshenko beams made of axially FGM subjected to multiple moving point loads”, Structural Engineering and Mechanics, 53(5), pp. 981- 995. 9. Jha, D.K., T. Kant, R.K. Singh (2013), “A critical review of recent research on functionally graded plates”, Composite Structures, 96, pp. 833–849. 10. Karamanli (2017), “Bending behaviour of two directional functionally graded sandwich beams by using a quasi-3d shear deformation theory”, Composite Structures. 174, pp. 70–86. 11. J.B. Kosmatka, “An improved two-node finite element for stability and natural frequencies of axial-loaded Timoshenko beams”, Computer & Structures, 57, pp. 141–149. 40 12. M. Nemat-Alla and N. Noda (2000), “Edge crack problem in a semi-infinite fgm plate with a bi-directional coefficient of thermal expansion under two- dimensional thermal loading”, Acta Mechanica, 144, pp. 211–229. 13. M. Şimşek (2010), “Vibration analysis of a functionally graded beam under a moving mass by using different beam theories”, Composite Structures, 92(4), pp. 904-917. 14. M. Şimşek (2015), “Bi-directional functionally graded materials (BDFGMs) for free and forced vibration of timoshenko beams with various boundary conditions”, Composite Structures, 133, pp. 968–978. 15. M. Şimşek (2016), “Buckling of timoshenko beams composed of two- dimensional functionally graded material (2D-FGM) having different boundary conditions”, Composite Structures, 149, pp. 304–313. 16. Wang Z., Wang X., Xu G., Cheng S., Zeng T. (2016), “Free vibration of two- directional functionally graded beams”, Composite Structures, 135, pp. 191– 198. 17. Nguyen Quang Huan, Bui Manh Cuong, and Nguyen Dinh Kien (2016), “Seismic Analysis of Planar Functionally Graded Beams and Frames Using Direct Integration Method”, Proceedings of the 4th International Conference on Engineering Mechanics and Automation (ICEMA4), pp. 334-341. 18. Dinh Kien Nguyen, Quang Huan Nguyen, Thi Thom Tran, Van Tuyen Bui (2017), “Vibration of bi-dimensional functionally graded Timoshenko beams excited by a moving load”, Acta Mechanica, 228, pp. 141–155. 19. Nguyen T.K., Sab K., Bonnet G. (2008), “First-order shear deformation plate models for functionally graded materials”, Composite Structures, 83, pp. 25- 36. 20. Thai H.T., Vo T.P., Nguyen T.K. (2014), “Static and vibration analysis of functionally graded beams using refined shear deformation theory”, Meccanica, 49, pp.155–168. DOI 10.1007/s11012-013-9780-1. 41 PHỤ LỤC Chương trình tính toán số % Maincode.m % BELOW GEOMETRIC AND MATERIAL DATA clear; n1=3; n2=0; LT=10; % total beam length h=0.2; % beam height b=0.2; psi=5/6; xE=0; Em1=210*10^9; % Elastic modulus of metal 1 Em2=70*10^9; % Elastic modulus of metal 2 Ec1=390*10^9; % Elastic modulus of ceramic 2 Ec2=200*10^9; % Elastic modulus of ceramic 2 nuM1=0.3; % Poisson ratio of metal 1 nuM2=0.23; % Poisson ratio of metal 2 nuC1=0.3; % Poisson ratio of ceramic 1 nuC2=0.3; % Poisson ratio of ceramic 2 Gm1=Em1/(2*(1+nuM1)); % Shear modulus of metal 1 Gc1=Ec1/(2*(1+nuC1)); % Shear modulus of metal 2 Gm2=Em2/(2*(1+nuM2)); % Shear modulus of ceramic 1 Gc2=Ec2/(2*(1+nuC2)); % Shear modulus of ceramic 2 Rm1=7800; % Mass density of metal 1 Rm2=2702; % Mass density of metal 2 Rc1=3960; % Mass density of ceramic 1 Rc2=5700; % Mass density of ceramic 2 --------------------------------------------------------------------------------------------------------------------------- % Column % nELE=10; % number of elements L=LT/nELE; nNOD=nELE+1; % number of nodes cNOD=zeros(nNOD,2); zNOD=linspace(0,LT,nNOD); cNOD(:,2)=zNOD'; % coordinates of nodes pre=[1 2 3]; act=setdiff([1:nDOF]',pre); --------------------------------------------------------------------------------------------------------------------------- % Portal_frame % % nELE=3; % number of elements 42 % nNOD=nELE+1; % number of nodes % cNOD=[0 0 % coordinates of nodes % 0 5 % 5 5 % 5 0]; % ELE(1,1)=1; ELE(1,2)=2; % nodes of elements % ELE(2,1)=2; ELE(2,2)=3; % ELE(3,1)=3; ELE(3,2)=4; % nDOF=(nELE+1)*3; % number of D.O.F % pre=[1 2 3]; % act=setdiff([1:nDOF]',pre); --------------------------------------------------------------------------------------------------------------------------- % Multi_frame % % nELE=12; % number of elements % nNOD=10; % number of nodes % cNOD=[0 0 % coordinates of nodes % 0 5 % 0 10 % 0 15 % 0 20 % 5 20 % 5 15 % 5 10 % 5 5 % 5 0]; % ELE(1,1)=1; ELE(1,2)=2; % nodes of elements % ELE(2,1)=2; ELE(2,2)=3; % ELE(3,1)=3; ELE(3,2)=4; % ELE(4,1)=4; ELE(4,2)=5; % ELE(5,1)=5; ELE(5,2)=6; % ELE(6,1)=4; ELE(6,2)=7; % ELE(7,1)=3; ELE(7,2)=8; % ELE(8,1)=2; ELE(8,2)=9; % ELE(9,1)=6; ELE(9,2)=7; % ELE(10,1)=7; ELE(10,2)=8; % ELE(11,1)=8; ELE(11,2)=9; % ELE(12,1)=9; ELE(12,2)=10; % nDOF=(nNOD)*3; % number of D.O.F. % pre=[1 2 3 28 29 30]; % act=setdiff([1:nDOF]',pre); --------------------------------------------------------------------------------------------------------------------------- % Assym_frame % % nELE=10; % number of elements % nNOD=nELE+1; % number of nodes 43 % cNOD=[0 0 % coordinates of nodes % 0 4 % 0 8 % 0 12 % 0 16 % 0 20 % 1 16 % 2 12 % 3 8 % 4 4 % 5 0]; % ELE(1,1)=1; ELE(1,2)=2; % nodes of elements % ELE(2,1)=2; ELE(2,2)=3; % ELE(3,1)=3; ELE(3,2)=4; % ELE(4,1)=4; ELE(4,2)=5; % ELE(5,1)=5; ELE(5,2)=6; % ELE(6,1)=6; ELE(6,2)=7; % ELE(7,1)=7; ELE(7,2)=8; % ELE(8,1)=8; ELE(8,2)=9; % ELE(9,1)=9; ELE(9,2)=10; % ELE(10,1)=10; ELE(10,2)=11; % nDOF=nNOD*3; % total number of D.O.F. % pre=[1 2 3 31 32 33]; % act=setdiff([1:nDOF]',pre); % Aij.m Function [Ac1m1,Ac2m2,Bc1m1,Bc2m2,Dc1m1,Dc2m2,Gc1m1,Gc2m2] = Aij(Em1,Em2,Ec1,Ec2,Gm1,Gm2,Gc1,Gc2,n1,n2,h,b) syms z V=z/h+1/2; E1=(Ec1-Em1)*V^n1 + Em1; E2=(Ec2-Em2)*V^n1 + Em2; G1=(Gc1-Gm1)*V^n1 + Gm1; G2=(Gc2-Gm2)*V^n1 + Gm2; Ac1m1=b*int(E1,z,-h/2,h/2); % axial rigidity Ac2m2=b*int(E2,z,-h/2,h/2); Bc1m1=b*int(E1*z,z,-h/2,h/2); % bending coupling rigidity Bc2m2=b*int(E2*z,z,-h/2,h/2); Dc1m1=b*int(E1*z^2,z,-h/2,h/2); % bending rigidity Dc2m2=b*int(E2*z^2,z,-h/2,h/2); Gc1m1=b*int(G1,z,-h/2,h/2); % shear rigidity Gc2m2=b*int(G2,z,-h/2,h/2); 44 % Iij.m function [Ic1m1,Ic2m2,Jc1m1,Jc2m2,Kc1m1,Kc2m2] = Iij(Rm1,Rm2,Rc1,Rc2,n1,n2,h,b) syms z V=z/h+1/2; R1=(Rc1-Rm1)*V^n1 + Rm1; R2=(Rc2-Rm2)*V^n1 + Rm2; % Mass moments Ic1m1=b*int(R1,z,-h/2,h/2); Ic2m2=b*int(R2,z,-h/2,h/2); Jc1m1=b*int(R1*z,z,-h/2,h/2); Jc2m2=b*int(R2*z,z,-h/2,h/2); Kc1m1=b*int(R1*z^2,z,-h/2,h/2); Kc2m2=b*int(R2*z^2,z,-h/2,h/2);