Tạp chí Khoa học Công nghệ Xây dựng, NUCE 2020. 14 (4V): 1–15 PHÂN TÍCH PHI TUYẾN ỨNG XỬ UỐN CỦA TẤM BẰNG VẬT LIỆU FGM XỐP ĐẶT TRÊN NỀN ĐÀN HỒI PASTERNAK VỚI CÁC ĐIỀU KIỆN BIÊN KHÁC NHAU CÓ XÉT ĐẾN VỊ TRÍ THỰC CỦA MẶT TRUNG HÒA Nguyễn Văn Longa,∗, Trần Minh Túa , Lê Thanh Hảib , Vũ Thị Thu Trangc a Khoa Xây dựng dân dụng và công nghiệp, Trường Đại học Xây dựng, số 55 đường Giải Phóng, quận Hai Bà Trưng, Hà Nội, Việt Nam b Khoa Xây dựng, Trường Đại học Vinh, số 182 đường Lê Duẩn, thành phố Vinh, tỉnh Nghệ An, Việt Nam c Viện Cơ khí, Trường Đại học Hàng hải Việt Nam, số 484 đường Lạch Tray, quận Lê Chân, Hải Phòng, Việt Nam Nhận ngày 21/06/2020, Sửa xong 22/07/2020, Chấp nhận đăng 31/08/2020 Tóm tắt Bài báo phân tích phi tuyến tĩnh tấm bằng vật liệu FGM xốp đặt trên nền đàn hồi Pasternak, chịu uốn dưới tác dụng của tải trọng phân bố vuông góc với bề mặt tấm. Vật liệu FGM xốp với ba dạng phân bố lỗ rỗng khác nhau: đều, đối xứng, bất đối xứng được khảo sát. Các phương trình chủ đạo được xây dựng trên cơ sở lý thuyết tấm Mindlin có kể đến yếu tố phi tuyến hình học von Kárman và vị trí mặt trung hòa. Bằng việc sử dụng phương pháp Bubnov-Galerkin, lời giải giải tích theo phương pháp ứng suất sử dụng hàm Airy đã được thiết lập với các điều kiện biên khác nhau. Ví dụ kiểm chứng đã được thực hiện qua so sánh với các công bố của các tác giả khác trong trường hợp vật liệu đẳng hướng. Ảnh hưởng của các tham số vật liệu, kích thước hình học, các tham số nền đàn hồi và điều kiện biên đến độ võng và các thành phần nội lực trong tấm được khảo sát cụ thể qua các ví dụ số. Từ khoá: phân tích uốn phi tuyến; tấm vật liệu FGM xốp; mặt trung hòa; lý thuyết biến dạng cắt bậc nhất; điều kiện biên khác nhau. NONLINEAR BENDING ANALYSIS OF FUNCTIONALLY GRADED POROUS PLATES RESTING ON PASTERNAK ELASTIC FOUNDATION UNDER VARIOUS BOUNDARY CONDITIONS BASED ON NEUTRAL SURFACE POSITION Abstract In this paper, the static nonlinear bending analysis of functionally graded porous plates resting on Pasternak elastic foundation is presented. Porous materials with three different types of porosity distribution: uniform, non-uniform symmetric and non-uniform non-symmetric are considered. The governing equations are derived based on Mindlin plate theory and neutral surface position, taking to account von Kárman nonlinearity. The Airy’s stress function and Bubnov-Galerkin method are employed to obtained the analytical solution with different boundary conditions. The verifications are conducted by comparing with the results published in the available literature for the isotropic plates. The effect of material, geometric, elastic foundation parameters, and boundary conditions on deflection, internal force resultants is investigated in detail. Keywords: nonlinear bending analysis; functionally graded porous plate; neutral surface position; first-order shear deformation theory; various boundary conditions. https://doi.org/10.31814/stce.nuce2020-14(4V)-01 © 2020 Trường Đại học Xây dựng (NUCE) ∗ Tác giả đại diện. Địa chỉ e-mail: longnv@nuce.edu.vn (Long, N. V.) 1 Long, N. V., và cs. / Tạp chí Khoa học Công nghệ Xây dựng 1. Giới thiệu Việc tìm kiếm, nghiên cứu ứng dụng các loại vật liệu mới có các tính chất khác biệt dần thay thế các loại vật liệu truyền thống là xu thế của thời đại ngày nay. Vật liệu FGM xốp (porous material) được biết đến như là một loại vật liệu nhẹ, có khả năng hấp thụ năng lượng tốt, thường được sử dụng để chế tạo kết cấu sandwich, tấm tường, sàn cách âm, cách nhiệt. Ở vật liệu FGM xốp, các lỗ rỗng (pore) phân bố theo một phương nhất định trong kết cấu tạo nên sự thay đổi trơn và liên tục các đặc trưng cơ học của vật liệu. Kết cấu sử dụng vật liệu FGM xốp được sử dụng rộng rãi trong nhiều ngành công nghiệp như: hàng không, ô tô, đóng tàu, xây dựng dân dụng, . . . Vì thế việc tìm hiểu ứng xử cơ học của các kết cấu bằng loại vật liệu này luôn là một đề tài hấp dẫn, thu hút sự quan tâm của giới khoa học trong và ngoài nước. Phuong và cs. [1] xây dựng nghiệm giải tích phân tích uốn dầm bằng vật liệu có cơ tính biến thiên (FGM) có lỗ rỗng vi mô. Long và Huong [2] phân tích ổn định dầm FGM có các lỗ rỗng vi mô chịu các điều kiện biên khác nhau. Thang và cs. [3] sử dụng lý thuyết biến dạng cắt bậc nhất nghiên cứu ổn định đàn hồi và dao động riêng của tấm vật liệu FGM xốp với phân bố các lỗ rỗng là đều và không đều. Lời giải chính xác cho tần số dao động riêng của tấm dày làm bằng vật liệu FGM xốp được Rezae và Saidi trình bày trong [4]. Arani và cs. [5] nghiên cứu dao động tự do của tấm chữ nhật làm bằng vật liệu rỗng trên nền Winkler theo lý thuyết biến dạng cắt bậc ba của Reddy, tần số dao động xác định bằng phương pháp DQM (differential quadrature method). Rezaei và cs. [6] phân tích dao động tự do của tấm FGM có bọt rỗng phân bố đều và không đều trên cơ sở lý thuyết tấm với bốn ẩn số chuyển vị. Ảnh hưởng của vi bọt rỗng đến ứng xử uốn và dao động riêng của tấm FGM được Akbas khảo sát trong [7]. Tu và cs. [8] phân tích ổn định và sau ổn định của tấm rỗng không hoàn hảo dựa trên lý thuyết tấm cổ điển. Các nghiên cứu trên đây hầu hết đối tượng là tấm chữ nhật vật liệu FGM xốp liên kết khớp trên chu vi. Tấm bằng vật liệu FGM xốp chịu điều kiện biên bất kỳ đã được một số tác giả nghiên cứu, chẳng hạn Yang và cs. [9] phân tích ổn định và dao động riêng của tấm vật liệu FGM xốp bằng phương pháp Ritz. Zhao và cs. [10] xây dựng lời giải ba chiều chính xác cho tấm dày FGM có vi bọt rỗng với các điều kiện biên đàn hồi bất kỳ. Zhao và cs. [11] sau đó đã sử dụng phương pháp chuỗi Fourier cải tiến để phân tich dao động riêng của tấm vật liệu FGM xốp có các liên kết đàn hồi trên các cạnh. Pradhan và Chakraverty [12] dùng phương pháp Rayleigh–Ritz và lý thuyết tấm mỏng phân tích tĩnh tấm FGM với các điều kiện biên khác nhau. Demirhan và Taskin [13] phân tích uốn và dao động riêng tấm FGM có vi bọt rỗng với điều kiện biên Levy. Hiện tại, theo hiểu biết của các tác giả, các nghiên cứu về phân tích phi tuyến các kết cấu tấm bằng vật liệu FGM xốp không nhiều, các phân tích phi tuyến chủ yếu áp dụng cho vật liệu đẳng hướng [14–16], vật liệu composite [17, 18] và vật liệu FGM [19]. Thêm vào đó, các nghiên cứu về tấm bằng vật liệu FGM xốp với các điều kiện biên khác nhau mới chỉ dừng lại ở các bài toán phân tích dao động và ổn định, và chủ yếu là phân tích tuyến tính. Do vậy, mục đích của bài báo là phân tích phi tuyến ứng xử uốn của tấm bằng vật liệu FGM xốp đặt trên nền đàn hồi Pasternak với các điều kiện biên khác nhau. Hệ phương trình phi tuyến khảo sát đường cong tải - độ võng được thiết lập với tiếp cận theo ứng suất kết hợp phương pháp Bubnov- Galerkin. Ảnh hưởng của ba loại phân bố lỗ rỗng: đều, đối xứng và bất đối xứng cũng như hệ số mật độ lỗ rỗng, điều kiện biên, tham số kích thước tấm, nền đàn hồi đến độ võng và các thành phần nội lực sẽ được khảo sát. 2 Long, N. V., và cs. / Tạp chí Khoa học Công nghệ Xây dựng 2. Mô hình tấm bằng vật liệu FGM xốp Xét tấm chữ nhật bằng vật liệu FGM xốp có chiều dày h, kích thước theo phương các trục x, y là a (chiều dài), b (chiều rộng). Tấm được đặt trên nền đàn hồi Pasternak (Hình 1) với các hệ số nền: Kw hệ số độ cứng uốn (Winkler stiffness), K si (i = x, y) - hệ số độ cứng cắt (shear stiffness). Hình 1. Mô hình tấm chữ nhật xốp trên nền đàn hồi Các hằng số vật liệu FGM xốp biến thiên liên tục theo chiều dày tấm, phụ thuộc vào mật độ phân bố lỗ rỗng [20, 21]: - Phân bố đều: !2 ) 1 1 2p 2 1 − e0 − + 1 − (1) {E, G} = {E1 , G1 } (1 − e0 χ) ; χ = e0 e0 π π - Phân bố đối xứng: - Phân bố bất đối xứng:   πz  {E(z), G(z)} = {E1 , G1 } 1 − e0 cos h (2)   πz π  + {E(z), G(z)} = {E1 , G1 } 1 − e0 cos 2h 4 (3) trong đó E1 , G1 lần lượt là các giá trị lớn nhất của mô đun đàn hồi kéo - nén, mô đun đàn hồi trượt; E2 , G2 là các giá trị nhỏ nhất tương ứng (xem Hình 2). Hệ số Poisson được coi là không thay đổi theo tọa độ chiều dày. Hệ số mật độ lỗ rỗng e0 được tính theo: e0 = 1 − G2 E2 =1− E1 G1 (4) (0 < e0 < 1) Vị trí mặt trung hòa của tấm FGM xốp trong trường hợp phân bố bất đối xứng không trùng mặt trung bình, được xác định từ điều kiện [22]:  h/2   h/2  Zh/2  Z   Z      (z − C) E(z)dz = 0 ⇒ C =  zE(z)dz /  E(z)dz (5)     −h/2 −h/2 3 −h/2 Long, N. V., và cs. / Tạp chí Khoa học Công nghệ Xây dựng (a) Phân bố đều (b) Phân bố đối xứng (c) Phân bố bất đối xứng Hình 2. Tấm bằng vật liệu rỗng với các hàm mật độ phân bố lỗ rỗng khác nhau e ee 3. Lý thuyết biến dạng cắt bậc nhất Sử dụng khái niệm mặt trung hòa, các thành phần chuyển vị u, v, w của điểm bất kỳ có tọa độ (x, y, zns ) trong không gian tấm biểu diễn dưới dạng [23]: u(x, y, zns ) = u0 (x, y) + zns θ x (x, y); v(x, y, zns ) = v0 (x, y) + zns θy (x, y); w(x, y, zns ) = w0 (x, y) (6) trong đó u0 , v0 , w0 là các thành phần chuyển vị của điểm trên mặt trung bình theo các phương x, y, zns ; θ x , θy là các góc xoay của pháp tuyến mặt trung hòa quanh hai trục y, x. Các thành phần biến dạng có kể đến thành phần phi tuyến hình học von Kármán thể hiện như dưới đây [23]:       0   ) ( 0 ) ( ε      κx  ε x      x      γ xz γ xz        0  εy  κy  εy  = ; (7) + zns  =  0          γyz γyz      γ    0    κ , , w γ xy xy xy trong đó , ,w ) ,20,x) ) ε0x w20,y ε0y = u0,x + ; = v0,y + ; γ0xy = u0,y + v0,x + w0,x w0,y ; κ x = θ x,x ; κy = θy,y ; 2 2 ( , y ); = kèm ( x, các thành q qphần 0q q κ xy = θ x,y + θy,x ; γ0xz== w(0,x, y+) θ+xz; γyz = w0,y + θy . Dấu (, ) đi chuyển vị chỉ đạo hàm riêng theo biến tương ứng. Vật liệu FGM xốp được coi là đàn hồi tuyến tính, các thành phần ứng suất được xác định từ định , , w, , w, w luật Hooke:  q; ,;q q,Q,11q Q12 0   , ,,σ,x, ; , ) y, x. #( ) " (    εx            γ Q 0 σ       xz 55 xz σy  =  Q21 Q22 0   εy  (8) = ;         γyz 0 Q44 σyz   0  σ    0 Q γ xy 66 xy νE(zns ) E(zns ) E(zns ) . Tích phân các ; Q12 = Q21 = ; Q44 = Q55 = Q66 = 2 2 2 (1 + v) 1−ν 1−ν thành phần ứng suất theo chiều dày của tấm ta nhận được các thành phần nội lực:  0          εx      Nx  κx  Mx   A11 A12 0    C11 C12 0                           0   εy  Ny  κy  My  =  A12 A11 0   ; ;  =  C12 C11 0                          N       0   0 0 A66 0 0 C66 κ xy M xy γ xy xy (9) ( #( 0 ) ) " s γ xz Q xz A44 0 = s 0 0 A44 Qyz γyz trong đó Q11 = Q22 = 4 Long, N. V., và cs. / Tạp chí Khoa học Công nghệ Xây dựng h/2−C Z trong đó Ai j , Ci j =   −h/2−C   s = ks Qi j 1, z2ns dzns ; i j = 11, 12, 66; A44 h/2−C Z Q44 dzns . Hệ số hiệu chỉnh −h/2−C cắt k s = 5/6 được sử dụng trong nghiên cứu này. Từ (9) có thể thấy rằng, việc sử dụng mặt trung hòa đã giúp loại bỏ được các tương tác màng-uốn trong tấm. Nguyên lý thế năng cực tiểu được sử dụng để thiết lập các phương trình cân bằng của tấm [24], với dạng toán học như sau: 0 = δU P + δU F + δV (10) trong đó δU P , δU F , δV lần lượt là biến phân của thế năng biến dạng đàn hồi của tấm, thế năng biến dạng của nền và thế năng của tải trọng. Hệ phương trình cân bằng thu được có dạng [24]: N x,x + N xy,y = 0; N xy,x + Ny,y = 0; Q xz,x + Qyz,y + N x w0,xx + 2N xy w0,xy + Ny w0,yy − Kw w0 + K sx w0,xx + K sy w0,yy + q = 0; M x,x + M xy,y − Q xz = 0; (11) M xy,x + My,y − Qyz = 0 Các tham số điều kiện biên bao gồm: (un , Nn ) , (u s , Nns ) , (w0 , Qn ) , (θn , Mn ) , (θ s , Mns ) . Các chỉ số dưới n, s thể hiện phương pháp tuyến và tiếp tuyến của biên tấm. Sử dụng hàm ứng suất Airy ϕ(x, y) được định nghĩa: N x = ϕ,yy ; Ny = ϕ,xx ; N xy = −ϕ,xy (12) Khi đó, hai phương trình đầu trong (11) tự thỏa mãn. Sử dụng các quan hệ (9), (7) và (12), ba phương trình còn lại trong (11) được viết lại theo chuyển vị và hàm ứng suất: s s s s A44 w0,yy + A44 w0,xx + A44 θ x,x + A44 θy,y + ϕ,yy w0,xx − 2ϕ,xy w0,xy + ϕ,xx w0,yy −Kw w0 + K sx w0,xx + K sy w0,yy + q = 0; s s C11 θ x,xx + C66 θ x,yy + (C12 + C66 ) θy,xy − A44 θ x − A55 w0,x = 0; (13) s s (C12 + C66 ) θ x,xy + C66 θy,xx + C11 θy,yy − A44 θy − A44 w0,y = 0 Mặt khác, phương trình tương thích biến dạng của tấm chữ nhật nhận được [25]: ε0x,yy + ε0y,xx − γ0xy,xy = w20,xy − w0,xx w0,yy (14) Dựa trên các quan hệ (9) và (12), các thành phần biến dạng màng có thể được xác định thông qua các thành phần lực dọc và hàm ứng suất: A11 A12 A11 A12 Nx − 2 Ny = 2 ϕ,yy − 2 ϕ,xx ; 2 2 2 − A12 A11 − A12 A11 − A12 A11 − A212 A12 A11 A12 A11 Ny − 2 Nx = 2 ϕ,xx − 2 ϕ,yy ; ε0y = 2 2 2 2 A11 − A12 A11 − A12 A11 − A12 A11 − A212 1 1 γ0xy = N xy = − ϕ,xy A66 A66 ε0x = A211 Thay (15) vào phương trình tương thích (14), ta được:   !  ∂2 w0 2 ∂2 w0 ∂2 w0  4   − ∇ ϕ = D  ∂x∂y ∂x2 ∂y2  5 (15) (16) Long, N. V., và cs. / Tạp chí Khoa học Công nghệ Xây dựng   ∂4 ∂4 ∂4 2 . 1 − ν + + 2 ; D = A 11 ∂x4 ∂y4 ∂x2 ∂y2 Hệ gồm ba phương trình trong (13) và phương trình (16) là hệ phương trình chủ đạo để giải bài toán uốn theo phương pháp ứng suất. Đây là hệ phương trình phi tuyến với 4 ẩn số độc lập: w0 , θ x , θy , ϕ. trong đó ∇4 = 4. Lời giải giải tích Theo lý thuyết biến dạng cắt bậc nhất, lời giải giải tích được thiết lập bằng việc sử dụng phương pháp Bubnov-Galerkin cho tấm chữ nhật bằng vật liệu rỗng với các điều kiện biên. Trong bài báo này, các điều kiện biên được xem xét bao gồm: - Liên kết khớp 4 cạnh (SSSS): + Trường hợp 1 (SSSS-1): Tất cả bốn cạnh của tấm tựa bản lề và có thể tự do dịch chuyển (freely movable) trong mặt phẳng tấm. Các điều kiện biên tương ứng là: w0 = θ s = 0, Nns = 0, Mn = 0, Nn = Nn0 = 0 (17) + Trường hợp 2 (SSSS-2): Tất cả bốn cạnh của tấm tựa bản lề và không thể tự do dịch chuyển (immovable) trong mặt phẳng tấm: un = w0 = θ s = 0, Nns = 0, Mn = 0 (18) - Liên kết ngàm 4 cạnh (CCCC): + Trường hợp 1 (CCCC-1): Tất cả bốn cạnh của tấm liên kết ngàm và có thể tự do dịch chuyển trong mặt phẳng tấm: w0 = θn = θ s = 0, Nns = 0, Nn = Nn0 = 0 (19) + Trường hợp 2 (CCCC-2): Tất cả bốn cạnh của tấm liên kết ngàm và không thể tự do dịch chuyển trong mặt phẳng tấm: un = w0 = θn = θ s = 0, Nns = 0 (20) - Liên kết đối xứng ngàm 2 cạnh, khớp 2 cạnh (SCSC): + Trường hợp 1 (SCSC-1): Hai cạnh đối diện của tấm tựa bản lề, hai cạnh còn lại liên kết ngàm và có thể tự do dịch chuyển trong mặt phẳng tấm: Tại x = 0, a : w0 = θy = 0, N xy = 0, M x = 0, N x = N x0 = 0 Tại y = 0, b : w0 = θ x = θy = 0, N xy = 0, Ny = Ny0 = 0 (21) + Trường hợp 2 (SCSC-2): Hai cạnh đối diện của tấm tựa bản lề, hai cạnh còn lại liên kết ngàm và không thể tự do dịch chuyển trong mặt phẳng tấm: Tại x = 0, a : u0 = w0 = θy = 0, N xy = 0, M x = 0 Tại y = 0, b : v0 = w0 = θ x = θy = 0, N xy = 0 (22) trong đó N x0 , Ny0 là các lực dọc màng tác dụng lên các cạnh của tấm chữ nhật theo phương x, y tương ứng trong trường hợp các cạnh đó có thể tự do dịch chuyển; hoặc là phản lực trên các cạnh của tấm trong trường hợp các cạnh không thể dịch chuyển trong mặt phẳng. 6 Long, N. V., và cs. / Tạp chí Khoa học Công nghệ Xây dựng Các điều kiện u0 = 0 (tại x = 0, a) và v0 = 0 (tại y = 0, b) được thỏa mãn theo nghĩa trung bình [26, 27]: Zb Za Zb Za u0,x dxdy = 0; v0,y dxdy = 0 (23) 0 0 0 0 Trong trường hợp tổng quát, với cả ba điều kiện biên được xem xét, hàm ứng suất được chọn dưới dạng: x2 y2 ϕ = ϕ̄(x, y) + N x0 + Ny0 (24) 2 2 Khi điều kiện biên là có thể tự do dịch chuyển trong mặt phẳng: N x0 1 = b Zb N̂ x dy = 0; Ny0 1 = a 0 Za (25) N̂y dx = 0 0 Khi điều kiện biên là không thể dịch chuyển trong mặt phẳng, từ (23) ta xác định được các thành phần phản lực: Zb Za  A12 2  A11 2 1 w0,x + w dxdy; −ϕ̄,yy + N x0 = ab 2 2 0,y 0 0 (26) Zb Za   1 A11 2 A12 2 Ny0 = w0,x + w dxdy −ϕ̄,xx + ab 2 2 0,y 0 0 Với các điều kiện biên SSSS, ta chọn nghiệm chuyển vị dưới dạng khai triển: XX mπ nπ w0mn sin αm x sin βn y; αm = w0 = , βn = ; m, n = 1, 3, 5, . . . ; a b m n XX XX θx = θ xmn cos αm x sin βn y; θy = θymn sin αm x cos βn y m n m n trong đó w0mn , θ xmn , θymn là các hệ số cần xác định. Thay (27) vào (16), ta được:        y x cos β − β K cos α − α q s 1 p r          XXXX   +K2 cos α p + αr  x cos βq + β s  y ϕ̄ = w0pq w0rs    +K3 cos α p − αr x cos βq + β s y   p q r s         +K4 cos α p + αr x cos βq − β s y trong đó các hệ số K1 , K2 K3 , K4 được trình bày trong Phụ lục A. Với các điều kiện biên CCCC đã nêu ta chọn nghiệm chuyển vị dưới dạng: XX w0 = w0mn sin2 αm xsin2 βn y; m, n = 1, 2, 3 · · · m θx = XX m                    n θ xmn sin 2αm xsin2 βn y; θy = n XX m 7 (27) n θymn sin2 αm x sin 2βn y (28) (29) Long, N. V., và cs. / Tạp chí Khoa học Công nghệ Xây dựng Thay (29) vào (16), ta được:       K1 cos 2 α p − αr x cos 2 βq − β s y       +K2 cos 2 α p + αr x cos 2 βq + β s y      XXXX  +K3 cos 2 α p − αr x cos 2 βq + β s y     ϕ̄ = w0pq w0rs   +K4 cos 2 α p + αr x cos 2 βq − β s y + K5 cos 2α p x cos 2β s y p q r s       +K6 cos 2α p x cos 2 βq − β s y + K7 cos 2α p x cos 2 βq + β s y      +K8 cos 2 α p − αr x cos 2β s y + K9 cos 2 α p + αr x cos 2β s y trong đó các hệ số K1 ÷ K8 được trình bày trong Phụ lục B. Với các điều kiện biên SCSC đã nêu ta chọn nghiệm chuyển vị dưới dạng: XX w0 = w0mn sin αm xsin2 βn y; m = 1, 3, 5 · · · ; n = 1, 2, 3 · · · θx = m n X X m θ xmn cos αm xsin2 βn y; θy = n XX m Thay (31) vào (16), ta được:     XXXX  ϕ̄ = w0pq w0rs   p q r s   q r n     K1 cos α p − αr x cos 2 βq − β s y     +K2 cos α p + αr x cos 2 βq + β s y     +K3 cos α p − αr x cos 2 βq + β s y     +K4 cos α p + αr x cos 2 βq − β s y     +K5 cos α p − αr x cos 2β s y + K6 cos α p + αr x cos 2β s y s (31) θymn sin αm x sin 2βn y trong đó: các hệ số K1 ÷ K6 được trình bày trong Phụ lục C. Thay các biểu thức xác định ϕ̄ trong (28), (30) và (32) vào (26), ta được: XXXX XXXX (1) (2) N x0 = w0pq w0rs K pqrs ; Ny0 = w0pq w0rs K pqrs p        (30)      p q r          (32) (33) s Thay ϕ̄ và N x0 , Ny0 vào (24) ta xác định được hàm ứng suất ϕ(x, y); sau đó thay vào (13); theo đó ta được hệ phương trình cân bằng theo w0 , θ x , θy : XX  XXXXXX (33) (34) (35) + w0mn w0pq w0rs g(33) w0mn lmn + θ xmn lmn + θymn lmn mnpqrs + q = 0; m m n XX m n  (45) (43) (44) w0mn lmn + θ xmn lmn + θymn lmn n p q XX = 0; m n r s  (53) (54) (55) w0mn lmn + θ xmn lmn + θymn lmn =0 (34) Nhân các biểu thức trong phương trình (34) với các hàm riêng tương ứng rồi thực hiện tích phân trên toàn bộ miền A của tấm, ta được: XX  XXXXXX (33) (34) (35) w0mn Lmni + θ L + θ L w0mn w0pq w0rsG(33) xmn mni j ymn mni j + j mnpqrsi j + F i j = 0; m n m n XX  (43) (44) (45) w0mn Lmni + θ L + θ L xmn ymn j mni j mni j = 0; m n p q XX m n r s  (53) (54) (55) w0mn Lmni + θ L + θ L xmn ymn j mni j mni j = 0 n o (35) Nghiệm của hệ phương trình đại số phi tuyến (35) là véc tơ chuyển vị w0mn ; θ xmn ; θymn từ đó xác định được các phần chuyển vị, biến dạng, ứng suất, nội lực của bài toán phân tích phi tuyến tĩnh. Trong các phân tích tuyến tính, bỏ qua các thành phần biến dạng phi tuyến trong công thức (7). 8 Long, N. V., và cs. / Tạp chí Khoa học Công nghệ Xây dựng 5. Kết quả số và thảo luận Với nghiệm giải tích đã thiết lập ở phần trên, chương trình tính trên nền Matlab được viết để thực hiện các ví dụ số. Các kết quả phân tích là phi tuyến trừ những trường hợp được nói trước. Các công thức không thứ nguyên được sử dụng [28, 29]: ! K sy b2 K sx a2 1 a b Kw a4 q0 a4 ; J = = ; E = 1.0 GPa; P = w̄ = w0 , ; K0 = 0 0 h 2 2 E0 h3 E0 h3 ν E0 h3 ν E 1 h4 (36) 5.1. Ví dụ kiểm chứng Trong phần này, các tác giả tiến hành kiểm chứng độ tin cậy của chương trình máy tính và lời giải giải tích theo phương pháp ứng suất. Qua nghiên cứu tổng quan, hiện tại chưa có công trình khoa học nào phân tích uốn phi tuyến của tấm chữ nhật vật liệu FGM xốp. Do đó, các tác giả sẽ tiến hành kiểm chứng cho một trường hợp đặc biệt của vật liệu FGM xốp: vật liệu đẳng hướng. a. Ví dụ 1: Kiểm chứng độ võng của tấm đẳng hướng điều kiện biên khớp 4 cạnh Tấm vuông dày đẳng hướng (E = 7,8.106 psi, ν = 0,3) điều kiện biên khớp bốn cạnh (SSSS-1, SSSS-2) với h = 1 inch, a = b = 10h, dưới tác dụng của tải trọng phân bố đều. Độ võng không thứ nguyên tại tâm tấm w̄ được tính toán như trong Bảng 1, và so sánh với các nhóm tác giả: Putcha và Reddy [17] sử dụng phương pháp phần tử hữu hạn dựa trên lý thuyết biến dạng cắt bậc cao 5 ẩn chuyển vị, Kapoor và Kapania [18] sử dụng phương pháp phần tử đẳng hình học dựa trên lý thuyết biến dạng cắt bậc nhất, Yu và cs. [19] sử dụng phương pháp đẳng hình học dựa trên lý thuyết biến dạng cắt bậc nhất đơn giản với 4 ẩn số chuyển vị. Bảng 1. Độ võng không thứ nguyên w̄ của tấm vuông đẳng hướng điều kiện SSSS-1, SSSS-2 dưới tác dụng của tải trọng phân bố đều P Yu và cs. [19] Putcha và Reddy [17] Kapoor và Kapania [18] Bài báo Sai số δ (%)* 0,2812 0,5185 0,8672 1,3147 1,8679 0,2840 0,5244 0,8790 1,3341 1,8918 0,2829 0,5257 0,8848 1,3319 1,8450 0,40 0,24 0,66 0,16 2,48 0,2790 0,4630 0,6911 0,9575 1,2688 0,2784 0,4626 0,6910 0,9579 1,2696 0,2637 0,4455 0,6727 0,9315 1,2166 5,27 3,69 2,65 2,76 4,17 SSSS-1 6,25 12,5 25 50 100 0,2813 0,5186 0,8674 1,3150 1,8688 SSSS-2 6,25 12,5 25 50 100 - *Sai số so với kết quả của Kapoor và Kapania [18]. 9 Long, N. V., và cs. / Tạp chí Khoa học Công nghệ Xây dựng b. Ví dụ 2: Kiểm chứng độ võng của tấm đẳng hướng điều kiện biên SCSC-2 Bảng 2 thể hiện kết quả độ võng không thứ nguyên w̄ tại tâm tấm vuông đẳng hướng điều kiện biên hai cạnh đối diện tựa khớp, hai cạnh còn lại liên kết ngàm (SCSC-2) dưới tác dụng của tải trọng phân bố đều với h/a = 0,05, ν = 0,3, E = 0,3.107 psi. Các kết quả tính toán trong bài báo được so sánh với Lei [15] sử dụng phương pháp phần tử biên (the boundary element method) dựa trên lý thuyết tấm bậc nhất, Azizian và Dawe [16] sử dụng phương pháp dải hữu hạn (the finite strip method) sử dụng lý thuyết tấm Mindlin và phương pháp Rayleigh-Ritz theo lý thuyết tấm cổ điển. Bảng 2. Độ võng không thứ nguyên w̄ của tấm vuông đẳng hướng điều kiện biên SCSC-2 dưới tác dụng của tải trọng phân bố đều P Azizian và Dawe [16] (Lý thuyết tấm cổ điển) Azizian và Dawe [16] (Lý thuyết tấm Mindlin) Lei [15] Bài báo Sai số δ (%) 0,9158 4,5788 6,8681 9,1575 0,0191 0,0951 0,1416 0,1867 0,0199 0,0988 0,1469 0,1936 0,0199 0,0984 0,1455 0,1904 0,0198 0,0982 0,1461 0,1929 3,66 3,26 3,18 3,32 *Sai số so với kết quả của Azizian và Dawe [16] (Lý thuyết tấm cổ điển). Từ các kết quả kiểm chứng ở các ví dụ 1 và ví dụ 2, có thể thấy rằng lời giải giải tích nhận được bằng phương pháp ứng suất và chương trình tính trên nền Matlab tự viết có độ tin cậy. 5.2. Khảo sát ảnh hưởng của các tham số vật liệu, hình học, nền đàn hồi và điều kiện biên Xét tấm chữ nhật bằng bọt kim loại - metal foam (h = 0, 1 m, E1 = 200 GPa, ν = 1/3) đặt trên nền đàn hồi, dưới tác dụng của tải trọng phân bố đều P. Bảng 3. Độ võng không thứ nguyên w̄ của tấm vuông FGM xốp với số số hạng khác nhau trong khai triển chuỗi lượng giác kép Điều kiện biên m, n = 1 m, n = 2 m, n = 3 m, n = 4 SSSS-1 SCSC-1 CCCC-1 0,5977 0,3319 0,2282 0,5758 0,3208 0,2150 0,5782 0,3313 0,2275 0,5776 0,3310 0,2260 SSSS-2 SCSC-2 CCCC-2 0,4996 0,3165 0,2264 0,4774 0,3040 0,2129 0,4799 0,3139 0,2254 0,4793 0,3134 0,2238 Bảng 3 trình bày các kết quả phân tích phi tuyến độ võng không thứ nguyên w̄ của tấm vuông vật liệu FGM xốp (a/h = 10, e0 = 0,5, K0 = 100, J0 = 10) với các loại điều kiện biên khác nhau. Số số hạng trong các khai triển chuỗi lượng giác kép tăng từ m, n = 1 đến m, n = 4. Có thể thấy rằng nghiệm giải tích có sự hội tụ rõ ràng khi tăng m, n; và với chương trình tính bằng Matlab thực hiện trên máy tính cá nhân các kết quả có thể được xem là hội tụ khi lấy m, n = 3 (sai số lớn nhất về độ võng không thứ nguyên w̄ khi lấy m, n = 3 so với khi lấy m, n = 4 là 0,69% trong trường hợp biên CCCC-2). Hình 3 và Hình 4 khảo sát quy luật biến thiên của độ võng lớn nhất không thứ nguyên w̄ theo hệ số lỗ rỗng e0 (P = 10) và đường cong tải-độ võng (e0 = 0,5) của tấm chữ nhật vật liệu FGM xốp liên 10 Long, N. V., và cs. / Tạp chí Khoa học Công nghệ Xây dựng kết khớp 4 cạnh SSSS-1. Các phân tích tuyến tính (LL) và phi tuyến (NL) cho ba quy luật phân bố lỗ rỗng: đều, không đều đối xứng và không đều bất đối xứng được thực hiện. Các kết quả cho thấy: - Khi hệ số mật độ lỗ rỗng tăng, độ võng của tấm tăng; hệ số lỗ rỗng e0 càng tăng thì sự khác biệt về độ võng giữa tính toán tuyến tính và phi tuyến, giữa các dạng phân bố lỗ rỗng càng lớn. Độ võng tính toán theo phi tuyến luôn thấp hơn tính toán theo tuyến tính. Độ võng của tấm với dạng phân bố lỗ rỗng không đều đối xứng luôn thấp nhất, độ võng của hai dạng phân bố còn lại (đều và không đều bất đối xứng) không khác biệt lớn đối với cả hai phương pháp phân tích tuyến tính và phi tuyến. Như vậy tấm có dạng phân bố lỗ rỗng không đều đối xứng có độ cứng lớn nhất. - Khi phân tích tuyến tính: dạng phân bố lỗ rỗng không đều đối xứng có độ võng tăng gần như tuyến tính khi e0 tăng. Hai dạng phân bố lỗ rỗng còn lại (đều và đối xứng) tăng phi tuyến khi e0 tăng, tăng càng nhanh khi e0 càng lớn ww w e0 e wtheo hệ số Hình 3. Biến thiên độ võng w̄ của tấm w xốp rỗng e0 với các quy luật phân bố lỗ rỗng khác nhau Hình 4. Biến thiên độ võng w̄ của tấmwtheo tải trọng phân bố đều P với các quy luật phân bố lỗ rỗng khác nhau Hình 5 và Hình 6 khảo sát quy luật biến thiên của độ võng lớn nhất không thứ nguyên w̄ theo tỷ số kích thước tấm a/h (b/a = 1) và theo tỷ số kích thước cạnh b/a (a/h = 10) của tấm xốp với các điều kiện biên khác nhau. Các kết quả cho thấy: - Về ảnh hưởng của điều kiện biên: rõ ràng là các biên SSSS có độ võng lớn nhất, sau đó đến biên SCSC, biên CCCC có độ võng nhỏ nhất; các biên hạn chế chuyển vị trong mặt phẳng có độ võng bé hơn so với biên không hạn chế chuyển vị trong mặt phẳng tương ứng; ảnh hưởng của việc hạn chế chuyển vị trong mặt phẳng lên độ võng lớn nhất đối với biên SSSS, sau đó là biên SCSC, rồi đến biên CCCC. - Về ảnh hưởng của tỷ số a/h : khi tăng tỷ số a/h, độ võng không thứ nguyên w̄ giảm nhanh khi a/h còn nhỏ (tấm dày, a/h ≤ 10); sau đó độ võng giảm chậm lại và gần như không đổi khi a/h đủ lớn (a/h ≥ 30). - Về ảnh hưởng của tỷ số b/a : khi tỷ số b/a tăng, độ võng lớn nhất không thứ nguyên w̄ tăng, khi b/a đủ lớn (b/a ≥ 2) độ võng hầu như không thay đổi. Hình 7 thể hiện đường cong tải-độ võng của tấm bằng vật liệu FGM xốp, điều kiện biên SSSS-1 với 4 trường hợp hệ số nền khác nhau. Các kết quả cho thấy: khi tăng các hệ số nền K0 , J0 , độ võng của tấm giảm đáng kể; và ảnh hưởng của hệ số nền J0 lớn hơn so với hệ số nền K0 . Hình 8 khảo sát quy luật biến thiên của các thành phần mô men uốn M x , My tại tâm tấm theo tải 11 Long, N. V., và cs. / Tạp chí Khoa học Công nghệ Xây dựng trọng phân bố Pcủa tấm xốp chữ nhật (b/a = 2) với các điều kiện biên khác nhau. Các kết quả cho thấy: với mọi mức tải P, mô men M x lớn hơn khá nhiều so với mô men uốn My . Một điều thú vị là với biên SCSC, khi tăng tải P, mô men uốn M x tăng trong khi mô men My tăng khi P còn nhỏ (P ≤ 12), sau đó lại giảm. w w w Hình 5. Biến thiên độ võng w̄ của tấmw theo w tỷ số kích thước tấm a/h với các điều kiện biênakhác /h nhau Hình 6. Biến thiên độ võng w̄ của tấm w theo tỷ số kích thước cạnh b/a với các điều kiện biên b /akhác nhau a /h w w Mn Hình 7. Biến thiên độ võng w̄ của tấmwtheowtải trọng phân bố đều P với các hệ số nền khác P nhau , của Hình 8. Biến thiên mô men uốn M x , My (Nm/m) tấm theo tải trọng phân bố đều P với các điều kiện biên khác nhau P P 6. Kết luận K , Jmô, hình tính toán phi tuyến độ võng và các thành phần nội lực trong tấm J chữ Bài báo xây dựng , JFGM , xốp đặt trên nền đàn hồi chịu uốn với một số điều kiện biên khác nhau. Với nhật bằng vậtKliệu . hệ quy chiếu là mặt trung hòa, . các tương tác màng uốn được loại bỏ. Nghiệm giải tích thu được bằng phương pháp Bubnov-Galerkin, chương trình tính tự viết trên nền Matlab được kiểm chứng, với tấm đẳng hướng cho thấy đủ tin cậy. Các khảo sát số đã được thực hiện cho phép đánh giá ảnh hưởng của P P / 12 2 P P, . . P 12 Long, N. V., và cs. / Tạp chí Khoa học Công nghệ Xây dựng các tham số hình học, vật liệu, nền đàn hồi và điều kiện biên đến độ võng và nội lực trong tấm. Các kết quả nhận được cho thấy phân tích uốn phi tuyến cho độ an toàn cao hơn phân tích tuyến tính. Có thể lựa chọn dạng phân bố lỗ rỗng và điều chỉnh mật độ phân bố lỗ rỗng để nhận được các kết quả mong muốn trong công tác tính toán, thiết kế kết cấu tấm sử dụng vật liệu FGM xốp. Lời cảm ơn Tác giả chân thành cảm ơn sự hỗ trợ tài chính của Quỹ phát triển Khoa học và Công nghệ Quốc gia (NAFOSTED), mã số 107.02-2018.322. Tài liệu tham khảo [1] Phuong, N. T. B., Tu, T. M., Phuong, H. T., Long, N. V. (2019). Bending analysis of functionally graded beam with porosities resting on elastic foundation based on neutral surface position. Journal of Science and Technology in Civil Engineering (STCE)-NUCE, 13(1):33–45. [2] Long, N. V., Hường, N. T. (2020). Phân tích ổn định kết cấu dầm vật liệu xốp chịu nén dọc trục với các điều kiện biên khác nhau. Tạp chí Khoa học Công nghệ Xây dựng (KHCNXD)-ĐHXD, 14(2V):97–106. [3] Thang, P. T., Nguyen-Thoi, T., Lee, D., Kang, J., Lee, J. (2018). Elastic buckling and free vibration analyses of porous-cellular plates with uniform and non-uniform porosity distributions. Aerospace Science and Technology, 79:278–287. [4] Rezaei, A. S., Saidi, A. R. (2015). Exact solution for free vibration of thick rectangular plates made of porous materials. Composite Structures, 134:1051–1060. [5] Arani, A. G., Khoddami Maraghi, Z., Khani, M., Alinaghian, I. (2017). Free vibration of embedded porous plate using third-order shear deformation and poroelasticity theories. Journal of Engineering, 2017. [6] Rezaei, A. S., Saidi, A. R., Abrishamdari, M., Mohammadi, M. P. (2017). Natural frequencies of functionally graded plates with porosities via a simple four variable plate theory: an analytical approach. Thin-Walled Structures, 120:366–377. [7] Akbaş, Ş. D. (2017). Vibration and static analysis of functionally graded porous plates. Journal of Applied and Computational Mechanics, 3(3):199–207. [8] Tu, T. M., Hoa, L. K., Hung, D. X., Hai, L. T. (2020). Nonlinear buckling and post-buckling analysis of imperfect porous plates under mechanical loads. Journal of Sandwich Structures & Materials, 22(6): 1910–1930. [9] Yang, J., Chen, D., Kitipornchai, S. (2018). Buckling and free vibration analyses of functionally graded graphene reinforced porous nanocomposite plates based on Chebyshev-Ritz method. Composite Structures, 193:281–294. [10] Zhao, J., Choe, K., Xie, F., Wang, A., Shuai, C., Wang, Q. (2018). Three-dimensional exact solution for vibration analysis of thick functionally graded porous (FGP) rectangular plates with arbitrary boundary conditions. Composites Part B: Engineering, 155:369–381. [11] Zhao, J., Wang, Q., Deng, X., Choe, K., Zhong, R., Shuai, C. (2019). Free vibrations of functionally graded porous rectangular plate with uniform elastic boundary conditions. Composites Part B: Engineering, 168:106–120. [12] Pradhan, K. K., Chakraverty, S. (2015). Static analysis of functionally graded thin rectangular plates with various boundary supports. Archives of Civil and Mechanical Engineering, 15(3):721–734. [13] Demirhan, P. A., Taskin, V. (2019). Bending and free vibration analysis of Levy-type porous functionally graded plate using state space approach. Composites Part B: Engineering, 160:661–676. [14] Levy, S. (1942). Square plate with clamped edges under normal pressure producing large deflections. Technical report, National Advisory Committee for Aeronautics. [15] Lei, X.-y., Huang, M.-K., Wang, X. (1990). Geometrically nonlinear analysis of a Reissner type plate by the boundary element method. Computers & Structures, 37(6):911–916. 13 Long, N. V., và cs. / Tạp chí Khoa học Công nghệ Xây dựng [16] Azizian, Z. G., Dawe, D. J. (1985). Geometrically nonlinear analysis of rectangular mindlin plates using the finite strip method. Computers & Structures, 21(3):423–436. [17] Putcha, N. S., Reddy, J. N. (1986). A refined mixed shear flexible finite element for the nonlinear analysis of laminated plates. Computers & Structures, 22(4):529–538. [18] Kapoor, H., Kapania, R. K. (2012). Geometrically nonlinear NURBS isogeometric finite element analysis of laminated composite plates. Composite Structures, 94(12):3434–3447. [19] Yu, T. T., Yin, S., Bui, T. Q., Hirose, S. (2015). A simple FSDT-based isogeometric analysis for geometrically nonlinear analysis of functionally graded plates. Finite Elements in Analysis and Design, 96: 1–10. [20] Chen, D., Yang, J., Kitipornchai, S. (2016). Free and forced vibrations of shear deformable functionally graded porous beams. International Journal of Mechanical Sciences, 108:14–22. [21] Barati, M. R., Zenkour, A. M. (2017). Investigating post-buckling of geometrically imperfect metal foam nanobeams with symmetric and asymmetric porosity distributions. Composite Structures, 182:91–98. [22] Larbi, L. O., Kaci, A., Houari, M. S. A., Tounsi, A. (2013). An efficient shear deformation beam theory based on neutral surface position for bending and free vibration of functionally graded beams#. Mechanics Based Design of Structures and Machines, 41(4):421–433. [23] Reddy, J. N. (2006). Theory and analysis of elastic plates and shells. CRC Press. [24] Reddy, J. N. (2017). Energy principles and variational methods in applied mechanics. John Wiley & Sons. [25] Brush, D. O., Almroth, B. O., Hutchinson, J. W. (1975). Buckling of bars, plates, and shells. McGrawHill, New York. [26] Librescu, L., Stein, M. (1991). A geometrically nonlinear theory of transversely isotropic laminated composite plates and its use in the post-buckling analysis. Thin-Walled Structures, 11(1-2):177–201. [27] Shen, H.-S. (2007). Thermal postbuckling behavior of shear deformable FGM plates with temperaturedependent properties. International Journal of Mechanical Sciences, 49(4):466–478. [28] Thai, H.-T., Choi, D.-H. (2011). A refined plate theory for functionally graded plates resting on elastic foundation. Composites Science and Technology, 71(16):1850–1858. [29] Zenkour, A. M. (2009). The refined sinusoidal theory for FGM plates on elastic foundations. International Journal of Mechanical Sciences, 51(11-12):869–880. Phụ lục A. Các hệ số trong công thức (28) của điều kiện biên SSSS  α p αr βq β s − α2p β2s K 1 =  2 2 ; (Khi p = r và q = s : K1 = 0); 2  α p − αr + β q − β s   D 2 2 4 α p αr βq β s − α p β s K2 =  2 2 ; 2  α p + αr + βq + β s   D 2 2 4 α p αr βq β s + α p β s K3 =  2 2 ; 2  α p − αr + βq + β s   D 2 β2 α α β β + α p r q s p s 4 K4 =  2 2 ; 2  α p + αr + βq − β s D 4  14 Long, N. V., và cs. / Tạp chí Khoa học Công nghệ Xây dựng Phụ lục B. Các hệ số trong công thức (30) của điều kiện biên CCCC   α p αr βq β s − α2p β2s K1 =   2 2 ; 2  α p − αr + βq − β s   D 2 2 16 α p αr βq β s − α p β s K2 =  2 2 ; 2  α p + αr + βq + β s   D − 16 α p αr βq β s + α2p β2s K 4 =  2 2 ; K5 = 2  α p + αr + β q − β s D 16 K7 =  D 2 2 8 α pβs α2p  + βq + β s 2 2 ; K8 =  (Khi p = r và q = s : K1 = 0)   D α p αr βq β s + α2p β2s − 16 K3 =  2 2 ; 2  α p − αr + βq + β s − D4 α2p β2s 2 ;  2 2 α p + βs D 2 2 8 α pβs α p − αr 2 + β2s 2 ; K6 =  D 2 2 8 α pβs α2p K9 =  2 2 ;  + βq − β s D 2 2 8 α pβs α p + αr 2 + β2s 2 ; Phụ lục C. Các hệ số trong công thức (32) của điều kiện biên SCSC  α p αr βq β s − α2p β2s K1 =  2 2 ; (Khi p = r và q = s : K1 = 0) 2  α p − α r + 4 βq − β s i  h  D 2 2 − D4 α p αr βq β s + α2p β2s 4 α p αr βq β s + α p β s K2 =  2 2 ; K3 =   2 2 ; 2 2   α p + αr + 4 βq + β s α p − αr + 4 βq + β s   D 2 2 − D4 α p αr βq β s − α2p β2s − D2 α2p β2s 2 α pβs K4 =  2 ; K 6 =  2 ;  2 2 ; K 5 =  2 2 2  2 2 α p + αr + 4 βq − β s α p − αr + 4β s α p + αr + 4β s D 4  15