›
›

Giáo trình Phương pháp tính cơ học kết cấu tàu thủy

pdf 194 trang ngocly 30
Download
Bạn đang xem 20 trang mẫu của tài liệu "Giáo trình Phương pháp tính cơ học kết cấu tàu thủy", để tải tài liệu gốc về máy bạn click vào nút DOWNLOAD ở trên

Tài liệu đính kèm:

  • pdfgiao_trinh_phuong_phap_tinh_co_hoc_ket_cau_tau_thuy.pdf

Nội dung text: Giáo trình Phương pháp tính cơ học kết cấu tàu thủy

  1. TRẦN CƠNG NGHỊ , ĐỖ HÙNG CHIẾN PHƯƠNG PHÁP TÍNH CƠ HỌC KẾT CẤU TÀU THỦY ĐẠI HỌC GIAO THƠNG VẬN TẢI TP. HỒ CHÍ MINH 6-2009
  2. (trang này để trống)
  3. 2 TRẦN CƠNG NGHỊ, ĐỖ HÙNG CHIẾN PHƯƠNG PHÁP TÍNH CƠ HỌC KẾT CẤU TÀU THỦY ĐẠI HỌC GIAO THƠNG VẬN TẢI TP HỒ CHÍ MINH TP HỒ CHÍ MINH 6-2009
  4. 3 Mục lục Mở đầu 5 Chương 1 Phương pháp biến phân và trọng hàm dư 6 1. Phép biến phân 6 2. Các phương pháp nhĩm trọng hàm dư 23 Chương 2 Phương pháp sai phân hữu hạn 32 1. Hàm một biến 32 2. Phương pháp lưới cho bài tốn 2 chiều 37 3. Xoắn dầm 41 4. Bài tốn trường 2D với biên cong 42 5. Phương pháp sai phân hữu hạn trên cơ sở phép biến phân 44 6. Dao động dầm 51 7. Dao động tấm 52 8. Ổn định tấm 54 Chương 3 Phương pháp phần tử hữu hạn 57 1. Phương pháp phần tử hữu hạn 57 2. Thứ tự giải bài tốn cơ học kết cấu theo phương pháp phần tử hữu hạn 59 3. Ma trận cứng phần tử. Ma trận cứng hệ thống 61 4. Áp đặt tải 63 5. Xử lý điều kiện biên 65 6. Giải hệ phương trình đaị số tuyến tính 66 7. Những phần tử thơng dụng trong cơ học kết cấu 70 8. Trạng thái ứng suất phẳng. Trạng thái biến dạng phẳng 82 9. Tấm chịuu uốn 88 10. Vật thể 3D 93 11. Nén ma trận. Khối kết cấu 95 12. Sử dụng phần mềm SAP và ANSYS tính tốn kết cấu 103 13. Phân tích kết cấu bằng ngơn ngữ MATLAB 112 14. Dao động kỹ thuật 142 Chương 4 Tính tốn độ tin cậy 163 1. Độ tin cậy 164 2. Tính tốn độ tin cậy 164 3. Xác định chỉ số an tồn, xác suất hư hoại 165 4. Phép tính thống kê và biến ngẫu nhiên 170 5. Các phương pháp tính 171 6. Phân tích những điều khơng chắc chắn từ tải và độ bền 181 7. Chọn hàm phân bố 182 8. Phân tích độ tin cậy hệ thống 182 9. Xác định các hệ số sử dụng 183 10. Thủ tục phân tích độ tin cậy kết cấu 189 11. Độ bền thân tàu 190 Tài liệu tham khảo 193
  5. 4 Ký hiệu chính A diện tích, area b chiều rộng, beam, width c, C hệ số, coefficient d đường kính, diameter D độ cứng tấm, flexural rigidity of plate E mộ đun đàn hồi, modulus of elasticity f hàm, function F lực, lực cắt, force, shear force G mộ đun đàn hồi (cắt ), shear modulus g gia tốc trọng trường, gravity constant h chiều cao, depth, heigh I, П, F phiếm hàm, functional I, J momen quán tính mặt cắt, moment of inertia Jp momen quán tính trong hệ độc cực, polar moment of inertia K, k hệ số, coefficient K, k độ cứng L, l chiều dài, length M momen, moment m khối lượng, mass N lực dọc trục, axial force P tải, load P cơng suất, power p áp suất, pressure Q tải, load q tải phân bố, distributed load R hàm sai số, residual function R, r bán kính, radius S diện tích, area T, MT momen xoắn, torque, couple t chiều dày, thickness t thời gian, time U thế năng, potential energy u0 thế năng đơn vị, strain energy per unit volume V lực cắt, shear force W trọng lượng, weight W,w cơng ngoại lực, work α gĩc nĩi chung, angle generally β gĩc nĩi chung, angle generally δ, Δ, w chuyển dịch vị trí, deflection δ tốn tử biến phân, variational operator γ biến dạng gĩc, shear strain θ gĩc, chuyển vị gĩc, angle, angle deflection Π thế năng, potential energy ε biến dạng , strain σ ứng suất nĩi chung, stress, generally η hệ số nĩi chung, coefficient generally ν hệ số Poisson, Poisson’s coefficient φ, ψ vector riêng, eigenvector ρ mậ t độ, density γ trọng lượng riêng, specific weight τ, T chu kỳ, perio ω tần số gĩc, circular frequency, generally ωn tần số riêng , natural frequency, generally
  6. 5 Mở đầu Cuốn sách “PHƯƠNG PHÁP TÍNH CƠ HỌC KẾT CẤU TÀU THỦY” trình bày các phương pháp tính cần cho việc xử lý những vấn đề thuộc cơ học kết cấu. Đây là phần khơng tách rời của bộ sách giành cho cơ học kết cấu tàu thủy, cần cho những người quan tâm cơ học kết cấu tàu thủy và cơng trình ngồi khơi. Ba cuốn sách đã được phát hành: “Cơ học kết cấu tàu thủy”, “Sức bền tàu”, “Dao động tàu thủy” cần đến các phương pháp tính trình bày trong sách này lúc xử lý các đề tài. Các chương trong sách sẽ trình bày những đề tài được quan tâm nhiều hiện nay. Chương đầu bàn về ứng dụng phương pháp biến phân kinh điển, dùng hiệu quả hàng trăm năm trong tốn và cơ học, giải những bài tốn uốn dầm, uốn tấm. Phương pháp Ritz cĩ sử dụng hàm thử và phép biến phân cùng các ứng dụng để giải bài tốn cơ học chất rắn nĩi chung, dầm và tấm nĩi riêng, là phần cần để ý của chương. Các phương pháp cĩ sử dụng hàm thử song khơng qua giai đoạn tính biến phân giới thiệu cùng chương mang tên gọi chung là phương pháp trọng hàm dư. Chương tiếp theo giới thiệu phương pháp sai phân hữu hạn hiện là phương pháp hữu hiệu trong tốn tính và trong cơ học. Tại chương này người đọc gặp cách xây dựng bài tốn và giải bài tốn cơ học kết cấu theo cách làm quen thuộc trước nay. Phương pháp sai phân hữu hạn đang phát huy tác dụng lớn ngày nay và chắc cịn tác dụng dài lâu. Bên cạnh đĩ những cách làm theo hướng đổi mới thủ tục tính tốn cho phương pháp truyền thống trình bày trong chương này giúp bạn đọc xem xét vấn đề đầy đủ, cĩ tính thời sự. Cĩ thể phát biểu rằng những cách làm mới khơng thay đổi nội dung phương pháp sai phân hữu hạn song làm cho nĩ bắt kịp tiến bộ trong lĩnh vực tốn tính. Những cơ sở của phương pháp tính phần tử hữu hạn và ứng dụng của nĩ xử lý những bài tốn cơ học kết cấu giới thiệu trong sách giúp bạn đọc làm quen và cĩ điều kiện nâng cao khả năng tính tốn theo phương pháp rất hữu hiệu này. Chương bốn trình bày các phương pháp tính đang dùng phổ biến trong mơn học “Độ tin cậy kết cấu”. Các thủ tục tính trình bày tại đây giúp người đọc xác định đúng và nhanh trong điều kiện cĩ thể các thơng số liên quan độ tin cậy kết cấu dân dụng nĩi chung và của tàu thủy nĩi riêng. Mỗi chương của sách ngồi phần lý thuyết và hướng dẫn tính tốn đều cĩ những ví dụ minh họa. Những người chuẩn bị sách cố ý trình bày những ví dụ độ phức tạp khơng cao, người đọc dễ dàng kiểm tra bằng các phép tính thủ cơng. Tuy nhiên, với các bài tốn động lực học, khối lương tính tốn thường lớn, đề nghị bạn đọc sử dụng cơng cụ tính thích hợp khi tìm trị riêng và vecto riêng. Những người viết
  7. 6 Chương 1 PHƯƠNG PHÁP BIẾN PHÂN VÀ TRỌNG HÀM DƯ Các bài tốn cơ học kết cấu giải theo nhiều phương pháp khác nhau. Trong chương này của sách đề cập những cách giải dựa trên các phương pháp dùng hàm thử theo nghĩa kinh điển. Các phương pháp trực tiếp tìm lời giải bao gồm: phương pháp biến phân kinh điển (Direct Variational Method) và phương pháp trọng hàm dư (Weighted Residual Method). Phương trình vi phân chính yếu trình bày trạng thái cân bằng vật rắn, xem xét trong chương: L(u) - p = 0 trong miền V, (a) và các điều kiện biên: B(u) - q = 0 trên biên S = Su + Sp (b) trong đĩ u – là hàm chuyển vị, nếu khơng giải thích khác, p – tải Biểu thức (b) được hiểu cụ thể theo cách diễn giải tại hình 1.1: * Điều kiện động học u = u tại Su; * Điều kiện động lực học q = q tại Sp. L và B là những tốn tử vi phân. Tốn tử thường gặp cĩ thể là ∇, ∇2, ∇4 = ∇2∇2 1 PHÉP TÍNH BIẾN PHÂN Hình 1.1 Điều kiện biên Phép tính biến phân liên quan vấn đề xác định cực trị, tức maximum hoặc minimum của các phiếm hàm. Phiếm hàm (functional) hiểu là hàm của các hàm. Trong chừng mức nhất định phiếm hàm cĩ nét tương đồng với hàm số chúng ta vẫn quen, điểm khác nhau cần nhắc đến, hàm số theo nghĩa thơng thường là hàm của các biến, cịn hàm đĩng vai trị phiếm hàm của các hàm. Việc chính của tính tốn biến phân là tìm hàm, ví dụ hàm u(x), với x – biến độc lập, để phiếm hàm dưới dạng tích phân giới hạn x1, x2: x2 I= ∫ F( u , u ', , x ) dx (1.1) x1 đạt cực trị. Trong tích phân này u’ = du/dx, I và F cùng được gọi phiếm hàm. Với những vấn đề thuộc cơ học kết cấu: I ≡ Π = U – W U – cơng biến dạng, W – cơng của ngoại lực. Nghiệm gần đúng tìm từ biểu thức: u~()()() x= u x + δ u x (1.2) Hình 1.2 Hàm u và biến phân δu
  8. 7 trong đĩ u(x) – nghiệm chính xác, nếu tồn tại, δu(x) cĩ tên gọi biến phân. δ là tốn tử biến phân. Phép tính biến phân hàm I: δ (∫Fdx)= ∫ ()δ F dx (1.3) ⎛ du ⎞ d Và δ ⎜ ⎟ = ()δu (1.4) ⎝ dx ⎠ dx ∂F ∂F δF = δu + δu' (1.5) ∂u ∂u' Điều kiện cần để I đạt cực trị: x x 2 ⎛ ∂F ∂F ⎞ 2 δI = ⎜ δu + δu'⎟ dx= δ Fdx = 0 (1.6) ∫ ∂u ∂u' ∫ x1 ⎝ ⎠ x1 1.1 PHƯƠNG PHÁP RITZ Phương pháp Ritz xây dựng trên cơ sở phép biến phân. Phương pháp Ritz1 tìm cách thay thế biến u, cĩ thể chọn ví dụ bài tốn một chiều u(x), trong x2 phiếm hàm (1.1): I= ∫ F( u , u ', , x ) dx , bằng nghiệm gần đúng dưới dạng hàm xấp xỉ: x1 N ~ u= ∑ ai f i (1.7) i=1 Hàm u, chúng ta đã gặp trong các bài tốn cơ học khác nhau, thể hiện tại (a), để tiện xem xét cĩ thể coi là hàm chuyển vị trong ví dụ tiếp theo. Hàm cơ sở hay cịn gọi hàm thử fi, i =1,2, , N phải thoả mãn các điều kiện biên (b) S = Sp + ~ Su, tức là điều kiện động lực học trên Sp, và điều kiện động học tại biên Su. Hàm xấp xỉ u liên tục đến bậc r-1, trong đĩ r – bậc đạo hàm cao nhất trong I. ~ Thay u vào I, cơng thức (1.1), tích phân I trở thành hàm của các ẩn ai. Phiếm hàm I tương đương hàm tổng thế năng Π gặp trong những bài tốn cơ học kết cấu Π = U – W, trong đĩ U – cơng biến dạng2, W – cơng ngoại lực3. Điều kiện cần để I đạt cực trị là: ∂I() u =0i = 1,2,L , n (1.8) ∂ai Xác định hàm I trong các bài tốn cơ học kết cấu cĩ thể tiến hành theo cách gán I bằng tổng năng lượng hệ thống Π = U – W. 1 Cơng biến dạng vật thể làm từ vật liệu đàn hồi: U = ∫{}{}εT σ dV , (1.9) 2 V Cơng do ngoại lực tác động lên vật thể: W= ∫{}{} pT u dS , (1.10) S trong đĩ {ε}= [C]{σ} – vector biến dạng, {σ} = [D]{ε} – vector ứng suất, 1 Ritz W., “Über eine neue Methode zur Lưsung gewissen Variations-Problem der mathematischen Physik”, J. Rein Angew. Math. (1909). 2 strain energy 3 external work due to applied loads
  9. 8 {p} – vector ngoại lực {u} – vecto chuyển vị. 1 Π = {}{}{}{}εT σ dV− p T u dS )( 1 1 . 1 2 ∫ ∫ V S p Thay hàm Π của hàm u vào vị trí phiếm hàm I, xác định hàm u đảm bảo tổng thế năng đạt minimum. Từ phép biến phân xác định biểu thức δΠ: ⎧ ⎫ ⎪1 T T ⎪ δΠ = δ⎨ {}{}{}{} ε σ dV− p u dS⎬ ) ( 2 1 . 1 2 ∫ ∫ ⎩⎪ V S p ⎭⎪ Giải bài tốn cơ học vật rắn chúng ta nhận phương trình: ∫δ{ ε }T [D ]{ ε } dV= ∫{ p }T δ{} udS ) 3 1 . 1 ( V S p Trong đĩ {p} – tải bên ngồi tác động lên biên Sp vật thể đang xem xét. Từ đây cĩ thể viết: δU = δW hoặc δ(U –W) = 0 ).( 4 1 . 1 Trường hợp bài tốn ba chiều hàm chuyển vị cĩ thể thể hiện: {u}= [ u v w]T ⎫ u= ∑ aiϕ i (,,) x y z ⎪ i ⎪ v= ∑biψ i ( x , y , z )⎬ ) ( 5 1 . 1 i ⎪ w= ∑ciθ i ( x , y , z )⎪ i ⎭ trong đĩ ai, bi, ci - các hệ số cần xác định, đĩng vai trị tọa độ suy rộng, ϕi, ψi, θi – các hàm cơ sở hay cịn gọi hàm thử. Hàm cơ sở thoả mãn các điều kiện biên tại S = Sp + Su. Biến phân hàm chuyển vị xác định như sau: ⎫ δu= ∑ δ ai ϕ i (,,) x y z ⎪ i ⎪ δv = ∑ δbi ψ i (,,) x y z ⎬ ) 6( 1 . 1 i ⎪ δw = ∑ δci θ i (,,) x y z ⎪ i ⎭ 2 2 2 2 1 ⎪⎧⎛ ∂u ⎞ ⎛ ∂v ⎞ ⎛ ∂w ⎞ ν ⎛ ∂u ∂v ∂w ⎞ Biết rằng U= u0 dV , trong đĩ u0 = ⎨⎜ ⎟ + ⎜ ⎟ + ⎜ ⎟ + ⎜ + + ⎟ + 2 ∫ ∂x⎜ ∂ y ⎟ ∂z 1− 2ν ⎜ ∂x∂ y ∂z ⎟ V ⎩⎪⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ 2 2 2 1 ⎡⎛ ∂u ∂v⎞ ⎛ ∂ v ∂w ⎞ ⎛ ∂w ∂u ⎞ ⎤⎪⎫ + ⎢⎜ + ⎟ +⎜ + ⎟ +⎜ + ⎟ ⎥⎬ 2 ⎜ ∂y∂ x⎟ ⎜ ∂ z ∂yx ⎟ ∂x ∂z ⎣⎢⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎦⎥⎭⎪ và:
  10. 9 ∂u ∂ ⎫ δεx = δ= δu ⎪ ∂x∂ x ⎪ L ⎬ ⎛ ∂u ∂w ⎞ ∂ ∂ ⎪ δγ= δ + = δu + δw zx ⎜ ⎟ ⎪ ⎝ ∂z ∂x⎠ ∂ z ∂x ⎭ cĩ thể viết: ⎡ ∂ ∂ ∂ ⎛ ∂ ∂ ⎞ ⎛ ∂ ∂ ⎞ ⎛ ∂ ∂ ⎞⎤ δU = σδσδσδτδδτδδτδδu+ v+ w+ u+ w + ⎜ v+ w⎟+ ⎜ u+ v⎟ dxdydz ∫∫∫⎢ x y z zx ⎜ ⎟ yz ⎜ ⎟ xy⎜ ⎟⎥ ⎣ ∂x ∂y∂ z ⎝∂z ∂x ⎠ ⎝∂z∂ y ⎠ ⎝∂y ∂x ⎠⎦ Thay giá trị δu, δv, δw từ biểu thức (1.16) vào phương trình xác định δU và δW tiếp tục xác định biến phân δΠ, nhận được cơng thức sau của phương pháp Ritz. ⎛ ∂Π ∂Π ∂Π ⎞ ⎜ ⎟ δΠ = ∑⎜ δai + δbi + δci ⎟ , i=1,2, ,n (1.17) i ⎝ ∂ai ∂bi ∂ci ⎠ Từ biểu thức (1.17), với δai, δbi, δci khác 0, cĩ thể viết: ∂Π ∂Π ∂Π = 0;= 0;= 0i = 1,2, , n (1.18) ∂ai ∂bi ∂ci Từ đây đưa đến lập hệ phương trình đại số tuyến tính chứa các ẩn ai, bi, ci. Cơng biến dạng dầm Thế năng dầm bị tác động bởi các lực kéo, nén, cắt, momen uốn, momen xoắn: 2 2 2 1 ⎛ M2 dx M dx M2 dx N dx F dx F2 dx ⎞ U = ⎜ T + y +z + + k y + k z ⎟ (1.19) 2 ⎜ ∫GI ∫EI ∫EI ∫ AE z ∫∫GA y GA ⎟ ⎝ t x z ⎠ trong đĩ MT - momen xoắn, My, Mz - momen uốn N – lực kéo, nén Fy, Fy – lực cắt Cơng biến dạng tấm chữ nhật axb, dày t. ⎧ ⎡ 2 2 2 2 ⎤⎫ D ⎪ 2 2 ⎛ ∂ w ⎞ ∂ w ∂ w ⎪ U =⎨() ∇w +2(1 −ν )⎢⎜ ⎟ − ⎥⎬dxdy (1.20) 2 ∫∫ ⎜ ∂x ∂ y ⎟ ∂x 2 ∂y 2 ⎩⎪ ⎣⎢⎝ ⎠ ⎦⎥⎭⎪ E t 3 trong đĩ D = 12(1−ν 2 ) Cơng biến dạng tấm trịn bán kính R, dày t 2 2 D⎪⎧⎛ ∂ 2 w 1∂w 1 ∂ 2 w ⎞ ⎡ ∂ ⎛ 1 ∂w⎞ ∂ 2 w ⎛ 1∂w 1 ∂ 2 w ⎞⎤⎪⎫ U = ⎨⎜ + + ⎟ +2(1 −ν )⎢ ⎜ ⎟ − ⎜ + ⎟⎥⎬rdrdϕ 2 ∫∫ ⎜ ∂r 2 r ∂r r 2 ∂ϕ 2 ⎟ ∂r⎜ r ∂ϕ ⎟ ∂r 2 ⎜ r ∂r r 2 ∂ϕ 2 ⎟ ⎩⎪⎝ ⎠ ⎣⎢ ⎝ ⎠ ⎝ ⎠⎦⎥⎭⎪ Cơng biến dạng vật thể 3D: 1 U= u dV và W= F.T u dS ∫ 0 ∫{}{} 2 V S
  11. 10 trong đĩ ε 1 u = εT σd ε =σεσεσε + + +2 τγ + 2 τγ + 2 τγ 0 ∫{}{ } ()x x y y x x xy xy yz yz zx zx 0 2 Thủ tục giải bài tốn cơ học kết cấu bằng phương pháp Ritz N ~ 1. Xây dựng hàm hoặc hệ hàm fi cho biến u: u= ∑ ai f i i=1 2. Xây dựng phiếm hàm I (ký hiệu tương đương Π) của u, ux, . . . ∂Π ∂I 3. Lập hệ phương trình đại số tuyến tính ≡ = 0 , i=1,2, ,n và xác định ai. ∂ai∂a i 4. Tìm nghiệm u. Phương pháp Ritz giải dầm Từ phương trình cân bằng dầm uốn, hình 1.3, cĩ thể xây dựng quan hệ: d 2 d2 w Phương trình chính: EJ = p dx 2 dx 2 d d2 w Lực cắt: V = EJ dx dx 2 d2 w Momen uốn: M= − EJ dx 2 dw Gĩc xoay: θ = − dx Hình 1.3 Uốn dầm Ví dụ 1: Áp dụng phương pháp Ritz xác định độ võng dầm liên tục dài L, EI = const, tựa trên hai gối tại đầu nút, dưới tác động tải trọng phân bố đều q = const, hình 1.4. Tiến hành xử lý bài tốn theo thủ tục đang nêu. Độ võng dầm:
  12. 11 ∞ ~ w() x= ∑ an f n n=1 trong đĩ fn (x) - hàm thử, an - hệ số cần xác định. Điều kiện biên: tại x =0: w(0) = 0; w’(0) ≠ 0. tại x = L: w(L) = 0; w’(L) ≠ 0. Để thoả mãn điều kiện biên hàm fn (x) cĩ thể mang dạng: nπ x fn ( x )= sin , n =1,2, L Hình 1.4 Dầm thẳng, tựa hai đầu, chịu tải Độ võng tính theo cơng thức: q(x) = const ∞ ~ nπ x w() x= ∑ an sin n=1 L Từ tính đối xứng của phương trình độ võng, chỉ cần giữ lại các hệ số cĩ chỉ số lẻ 1, 3, 5 Hàm đúng giờ đây chỉ cịn là: ~ nπ x w() x= ∑ an sin n=1,3, L Phiếm hàm I cho dầm bị uốn đơn thuần, bằng tổng thế năng biến dạng và cơng ngoại lực tác động lên dầm:I ≡ Π = U – W. Thế năng dầm uốn: 1 L U= ∫ EJ[] w"( x ) 2 dx 2 0 2 2 1 L ⎡ ∞ ⎛ nπ⎞ n π x⎤ U= EJ⎢ − a ⎜ ⎟ sin ⎥ dx 2 ∫ ∑ n L L 0 ⎣⎢ n=1 ⎝ ⎠ ⎦⎥ Cơng ngoại lực: LL nπ x W= q. w~( x ) dx= q . a sin dx ∫∫∑ n 00n=1,3, L 2 2 1 LL⎡ ∞ ⎛ nπ⎞ n π x⎤ ∞ nπ x Π ≡I = EJ⎢ − a ⎜ ⎟ sin ⎥ dx− q a sin dx 2 ∫∫∑ n L L ∑ n L 00⎣⎢ n=1,3, ⎝ ⎠ ⎦⎥ n=1,3, Hệ số ai xác định sau khi lấy đạo hàm của I theo ai. 2 2 ∂ ()UW− L ⎡ ∞ ⎛ nπ⎞ n π x⎤⎛ k π⎞ n π x L nπ x = EJ⎢ a ⎜ ⎟ sin ⎥⎜ ⎟ sin dx− qsin dx = 0 ∂ a ∫ ∑ n L L L L ∫ L {} 0 ⎣⎢ n=1 ⎝ ⎠ ⎦⎥⎝ ⎠ 0 Trong khoảng khơng gian 0 - L họ hàm (sin nxπ ) cĩ tính trực giao: L L mπ x n π x ⎧L / 2 m= n ∫ sin sin dx = ⎨ 0 L L ⎩ 0 m≠ n
  13. 12 4 ∂ ()UW− ⎛ kπ ⎞ L 2L = EJ⎜ ⎟ ak − q =0 k = 1,3, ∂ak ⎝ L ⎠ 2 kπ Từ đĩ: 4qL4 1 a = × k EJ k 5π 5 Biểu thức độ võng dầm w(x): 4qL4 ∞ 1 iπ x w~() x = sin i=1,3,5, 5∑ 5 EJπ i=1 i L Tại vị trí giữa dầm x = L/2 giá trị của hàm xấp xỉ như sau: 4 4 qL 4 w(L/2) = ≈ 0,013071 qL / EJ, khi sử dụng chỉ một hệ số a1. π 5 EJ 4 4 ⎛ 1 ⎞ qL 4 w(L/2) = ⎜1 − ⎟ ≈ 0,0130172 qL / EJ, nếu sử dụng a1 và a2. π 5 ⎝ 35 ⎠ EJ Lời giải “chính xác” theo phương pháp giải tích đưa ra kết quả: 5 qL4 w(L/2) = ≈ 0,0130208 qL4 /EJ. 384 EJ Moment uốn tính theo cơng thức: M(x) = EJ.w’’(x) Tại x = L/2 giá trị momen này được tính như sau: 4qL4 ∞ 1 iπ M(L/2) = - sin , i=1,3,5, 3 ∑ 3 π i=1 i 2 Ổn định dầm Ví dụ 2: Ổn định dầm dài L, độ cứng EI, chịu tác động lực nén N. Hàm chuyển vị được tìm dưới dạng: ∞ nπ x w() x= ∑ an sin n=1 L Với dầm chịu lực nén N, hàm W mang dạng: 1 L W= ∫ Nw'2 dx 2 0 2 1 L ⎡ ∞ ⎛ nπ⎞ n π L⎤ W= N a cos dx ∫ ⎢∑ n ⎜ ⎟ ⎥ 2 0 ⎣ n=1 ⎝ L ⎠ L ⎦ Tổng thế năng của dầm: 2 L ⎧ 2 2 ⎫ 1 ⎪ ⎡ ∞ ⎛ nπ⎞ n π x⎤ ⎡ ∞ ⎛ nπ⎞ n π x⎤ ⎪ Π =U − W = EI − a sin − N a cos dx ∫ ⎨ ⎢ ∑ n ⎜ ⎟ ⎥ ⎢∑ n ⎜ ⎟ ⎥ ⎬ 2 n=1 ⎝ L ⎠ L n=1 ⎝ L ⎠ L 0 ⎩⎪ ⎣⎢ ⎦⎥ ⎣ ⎦ ⎭⎪ Tiến hành đạo hàm Π theo ak, sẽ nhận được hệ phương trình:
  14. 13 2 2 ∂Π 1 L ⎪⎧ ⎡ ∞ ⎛ nπ⎞ n π x⎤⎛ k π⎞ k π x =⎨EI⎢ − an ⎜ ⎟ sin ⎥⎜ ⎟ sin − ∂ a 2 ∫ ∑ L L L L {} 0 ⎩⎪ ⎣⎢ n=1 ⎝ ⎠ ⎦⎥⎝ ⎠ ⎡ ∞ ⎛ nπ⎞ n π x⎤ k πk π x⎫ − N⎢∑ an ⎜ ⎟cos ⎥ cos ⎬dx = 0. ⎣ n=1 ⎝ L ⎠ L ⎦ L L ⎭ kπ x kπ x Nhờ tính trực giao của cos và sin , k =1,2, 3, trong đoạn (0, L): L L L nπ x k π x ∫ cos cos dx = 0 nếu k≠ n 0 L L L = nếu k = n; 2 sẽ nhận được: 2 ∂Π ⎡ ⎛ kπ ⎞ ⎤ = ak ⎢ EI⎜ ⎟ − N ⎥ = 0 ∂ak ⎣⎢ ⎝ L ⎠ ⎦⎥ Với ak ≠ 0, biểu thức trong dấu ngoặc vuơng phải bằng khơng, do vậy cĩ thể viết: k2π 2 EI N = L2 Chúng ta quan tâm đến giá trị nhỏ nhất của N, khi vượt qua giá trị đĩ dầm chuyển sang giai đoạn mất ổn định. Trong cơng thức cuối cĩ thể thấy, N đạt nhỏ nhất cho trường hợp k =1, biểu thức lực Euler cĩ dạng: π 2 EI N = E L2 Xây dựng phương trình ma trận Sử dụng cơng thức biến phân tổng năng lượng như tổng cọng biến phân thế năng và cơng ngoại lực của bài tốn uốn dầm cĩ thể viết: 1 L L Từ U= ∫ EJ[] w"( x ) 2 dx xác định: δU= ∫ δ w". EJ . w " dx 2 0 0 L L Từ W= ∫ q.() w x dx viết δW= ∫ q. δ wdx 0 0 L L δΠ = ∫ δw". EJ . w " dx− ∫ q .δ wdx = 0 (1.21) 0 0 Hàm chuyển vị w(x) theo cách làm ngày nay nên viết dưới dạng vecto như sau: w= Nu Các đại lượng liên quan {w} xác định theo cách sau: δθ = −δw'';"'' = − Nδ uδ w= Nδ u trong đĩ N – hàm hình dáng theo cách gọi ngày nay, u – chuyển vị tại các nút tính tốn.
  15. 14 L L L L TT TT T ⎡ T T ⎤ δΠ = ∫ δu N"EJ N"u dx −∫δ u N" pdx =δu⎢EJ ∫ N" N"dx u− ∫ N" pdx⎥ = 0 (1.22) 0 0 ⎣ 0 0 ⎦ L L Ký hiệu: K = EJ ∫ N"T N"dx và P= ∫ N"T pdx , cĩ thể viết phương trình cuối dạng: 0 0 Ku= P ) ( 3 2 . 1 Ví dụ 3: Xác định độ võng dầm, momen uốn và lực cắt dầm nêu tại hình 1.5. Hình 1.5 Dầm thẳng chịu tải phân bố tuyến tính Điều kiện biên: w(0) = 0; -θ(0) = w’(0) = 0 và w(L) = 0. Ký hiệu ξ = x/L trong các phép tính tiếp theo. Hàm hình dáng, bàn chi tiết tại mục “hàm nội suy”, tại đây chúng ta nhận như sau: []N =[(ξ2 − ξ 3) ( ξ 2−2 ξ 3 + ξ 4 )] Đạo hàm bậc hai của [N]: 1 []N" =[]()2 − 6ξ() 2 − 12 ξ + 12 ξ 2 L2 Thay hai biểu thức này vào cơng thức xác định các thành phần K và P : 2 2 3 EJ L ⎡ 4− 24ξ + 36 ξ 4− 36ξ + 96 ξ − 72 ξ ⎤ EJ ⎡4 0 ⎤ []K = dξ = 3 ∫0 ⎢ 2 3 3 4 ⎥ 3 ⎢ ⎥ L ⎣4− 36ξ + 96 ξ − 72 ξ 4 − 48 ξ + 192 ξ + 144 ξ ⎦ L ⎣0 0,8⎦ 2 3 4 L ⎧ ξ−2 ξ + ξ ⎫ ⎧1/ 30⎫ {}P= p0 L dξ = p0 L ∫0 ⎨ 2 3 4 5 ⎬ ⎨ ⎬ ⎩ξ−3 ξ + 3 ξ − ξ ⎭ ⎩1/ 60⎭ Vector {u} xác định từ quan hệ: 4 ⎧u1 ⎫ ⎧1/120⎫ p L ⎨ ⎬ = ⎨ ⎬ 0 ⎩u2 ⎭ ⎩1/ 48 ⎭ EJ Lời giải: 4 4 p0 L ⎡ 1 2 3 1 2 3 4 ⎤ p0 L 2 3 4 w() x = ⎢ ()ξ− ξ +() ξ −2 ξ + ξ ⎥ = ()7ξ− 12 ξ + 5 ξ EJ ⎣120 48 ⎦ 240EJ p L2 M()" x= − EJw = −0 ()7 − 36ξ + 30 ξ 2 120 p L V()''' x= − EJw =0 ()8 − 15ξ 30
  16. 15 Kết quả vừa trình bày chưa đáp ứng điều kiện momen tĩnh tại gối trái phải triệt tiêu. Cần thiết hiệu chỉnh hàm thử, trong trường hợp này là hàm hình dáng. Cĩ thể chọn hàm bậc cao hơn cho [N], ví dụ cĩ thể chọn: []N =[(ξ2 − ξ 3) ( ξ 3 − ξ 4 ) ( ξ4− ξ 5 )] Các phép tính tiếp theo: ⎡ 4− 24ξ + 36 ξ 2 12ξ− 60 ξ2 + 72 ξ3 24 ξ2 − 112 ξ3 + 120 ξ 4 ⎤ EJ 1 ⎢ ⎥ []K = 12ξ− 60 ξ2 + 72 ξ3 36 ξ2 − 144 ξ3 + 144 ξ4 72 ξ3 − 264 ξ4 + 240 ξ 5 dξ L4 ∫0 ⎢ ⎥ ⎢ 2 3 4 3 4 5 4 5 6 ⎥ ⎣24ξ− 112 ξ + 120 ξ 72 ξ− 264 ξ + 240 ξ 144 ξ−480 ξ + 400 ξ ⎦ ⎧ξ2− ξ 3 ⎫ ⎧1/ 30 ⎫ 1 ⎪ 3 4 ⎪ ⎪ ⎪ {}P= p0 (1 − ξ )ξ− ξ dξ = p0 L 1/ 60 ∫0 ⎨ ⎬ ⎨ ⎬ ⎪ 4 5 ⎪ ⎪ ⎪ ⎩ξ− ξ ⎭ ⎩1/105⎭ Hệ phương trình đại số [K]{u}={P} cĩ dạng: ⎡4 4 4 ⎤⎧u1 ⎫ ⎧1/ 30 ⎫ EJ ⎪ ⎪ ⎪ ⎪ ⎢4 24 / 5 26 / 5 ⎥ u = p L 1/ 60 3 ⎢ ⎥⎨ 2 ⎬ 0 ⎨ ⎬ L ⎪ ⎪ ⎪ ⎪ ⎣⎢4 26 / 5 208/ 35⎦⎥⎩u3 ⎭ ⎩1/105⎭ Từ đĩ xác định: p L {}u =0 []4 − 4 1 T 120 p L4 Và {}w = 0 ()4ξ2− 8 ξ 3 + 5 ξ 4 − ξ 5 120EJ Kết quả tính trính bày tại hình 1.6. Các hình từ trên xuống giới thiệu chuyển vị w(x), momen uốn M(x) Hình 1.6 Kết quả tính theo phương pháp Ritz và lựcc cắt F(x), Trong cùng hình kết quả tính theo phương pháp giải tích ghi lại tại đường cong đánh dấu A, tính theo phương án đầu trong ví dụ này đánh dấu bằng B. Kết quả tính theo phương án cải tiến trùng với đường A. Dao động ngang dầm Từ phương trình xác định momen uốn dầm: d2 w EI = −M() x ( a) dx 2 tiến hành lấy đạo hàm hai vế phương trình, nhận được các biểu thức sau: d d2 w EI = −F () x ( b) dx dx 2 d 2 d2 w EI = q () x ( c) dx 2 dx 2 Áp dụng nguyên lý D’Alembert vào phương trình (c), cĩ nghĩa thay lực quán tính vào vị trí của tải tĩnh q(x):
  17. 16 d2 w q() x thay bằng − m , dt 2 trong đĩ m – khối lượng đơn vị dài của dầm. Cĩ thể viết phương trình chuyển động ngang dầm: d4 w d2 w EI +m = 0 ( d) dx 4 dt 2 Điều kiện biên xác định cho ví dụ đang nêu: w = 0 ⎫ ⎬ tại x = 0 và x = L Hình 1.7 Dao động d ầm (e) dw/ dx = 0⎭ Áp dụng cơng thức Rayleigh để xác định tần số dao động riêng. Dao động ngang dầm thể hiện dạng điều hịa: w( x , t )= u ( x )cos ω t ( f) Thế năng và động năng hệ thống: 2 2 2 2 1 L ⎛ d w ⎞ 1 L ⎛ d w ⎞ U= EI⎜ ⎟ dx T= m() x ⎜ ⎟ dx ∫0 ⎜ 2 ⎟ ∫0 ⎜ 2 ⎟ 2 ⎝ dx ⎠ 2 ⎝ dt ⎠ Giá trị maximum của U và T: 2 2 2 L L 1 ⎛ d u ⎞ ω 2 U max = EI⎜ ⎟ dx T= m() x u dx ∫0 ⎜ 2 ⎟ ∫0 2 ⎝ dx ⎠ 2 Cơng thức Rayleigh áp dụng vào đây mang dạng: 2 2 1 L ⎛ d u ⎞ EI⎜ ⎟ dx 2 ∫0 ⎜ dx 2 ⎟ ω 2 = ⎝ ⎠ ( g) 1 L m() x u2 dx 2 ∫0 Hàm u(x) cĩ thể là: ⎛ 2πx ⎞ u( x )=⎜ 1 − cos ⎟ ( h) ⎝ L ⎠ Thay u(x) vào cơng thức Rayleigh, xác định tần số thứ nhất: 22,792 EI ω = ; ( i) 1 L2 m Thủ tục thực hiện theo phương pháp Ritz tiến hành như sau: ~ ux()()()= afxafx1 1 +2 2 +L+ afxn n ( ) ( j) 2 Thay biểu thức (j) vào cơng thức Rayleigh cho phép nhận ω là hàm của a1, a2, . . . Cần thiết chọn a1, a2, . . . để ω đạt minimum. ∂()ω2 ∂() ω 2 ∂()ω 2 = =L = = 0 ( k) ∂a1 ∂ a2 ∂an Điều này dẫn đến hình thành hệ phương trình đại số gồm n phương trình, chứa n ẩn a1, a2, . . . an
  18. 17 2 2 L L ∂ ⎛ d w ⎞ 2 ∂ 2 EI⎜ ⎟ dx − ω j m( x ) u dx = 0 j = 1, 2, . . . , n ( l) ∫0 ⎜ 2 ⎟ ∫0 ∂a j ⎝ dx ⎠ ∂a j Hãy chọn hàm u chứa 2 thành phần: ⎛ 2πx ⎞ ⎛ 4πx ⎞ u( x )= a1 ⎜ 1 − cos ⎟ +a2 ⎜1 − cos ⎟ (m) ⎝ L ⎠ ⎝ L ⎠ Cơng thức (g) giờ cĩ dạng: 4 8π EI 2 2 3 ()a1 +16 a2 ω 2 = L mL ()3a2 + 3 a2 + 4 a a 2 1 2 1 2 ∂()ω2 ∂() ω 2 Thỏa mãn điều kiện: = = 0 sẽ nhận hệ phương trình đại số: ∂a1 ∂ a2 2 16π EI ⎡1 0 ⎤⎧a1 ⎫ 2 ⎡3 2⎤⎧a1 ⎫ 3 ⎢ ⎥⎨ ⎬ = ω mL⎢ ⎥⎨ ⎬ L ⎣0 16⎦⎩a2 ⎭ ⎣2 3⎦⎩a2 ⎭ Lời giải của hệ phương trình: 22,35 EI ⎧a1 ⎫ ⎧ 1 ⎫ ω1 = 2 ⎨ ⎬ = ⎨ ⎬ L m ⎩a2 ⎭ ⎩0,575⎭ 124 EI ⎧a1 ⎫ ⎧ 1 ⎫ ω2 = 2 ⎨ ⎬ = ⎨ ⎬ L m ⎩a2 ⎭ ⎩−1,4488⎭ 1.2 PHƯƠNG PHÁP RITZ GIẢI TẤM MỎNG Phương pháp Ritz áp dụng tính độ võng, momen uốn, lực cắt tấm mỏng dựa vào phương pháp xác định biến phân hàm năng lượng tấm δΠ. Độ võng màng tấm trình bày dạng chuỗi: n ~ w(,)(,) x y= ∑ ai f i x y )( 4 2 . 1 i=1 Đưa hàm chuyển vị tấm vào phiếm hàm, trường hợp này là hàm năng lượng tồn phần của tấm Π. Các hằng số ai xác định sau khi thực hiện phép tính biến phân hàm năng lượng δΠ. Từ lý thuyết tấm cĩ thể ghi lại những biểu thức của momen uốn tấm: ⎛ ∂ 2 w ∂ 2 w ⎞ ⎛ ∂ 2 w ∂ 2 w ⎞ m= − D⎜ +υ ⎟; m= − D⎜ +υ ⎟ (1.25a) x ⎜ 2 2 ⎟ y ⎜ 2 2 ⎟ ⎝ ∂x ∂y ⎠ ⎝ ∂y ∂x ⎠ ∂ 2 w m= − D()1 −υ (1.25b) xy ∂x ∂ y Năng lượng biến dạng của tấm uốn tính từ biểu thức: 2 2 D⎪⎧⎛ ∂ 2 w ∂ 2 w ⎞ ⎡⎛ ∂ 2 w ⎞ ∂ 2 w ∂ 2 w⎤⎪⎫ U b = ⎨⎜ + ⎟ +2(1 −ν )⎢⎜ ⎟ − ⎥⎬dxdy (1.26) 2 ∫∫ ⎜ ∂x 2 ∂y 2 ⎟ ⎜ ∂x ∂ y ⎟ ∂x 2 ∂y 2 ⎩⎪⎝ ⎠ ⎣⎢⎝ ⎠ ⎦⎥⎭⎪
  19. 18 Hãy xét điều kiện biên các tấm khi xử lý vấn đề. Nếu tất cả các cạnh tấm bị ngàm, thành phần thứ hai biểu thức này trở thành 0. Điều kiện biên trong trường hợp này trở thành: w = 0 tại cạnh. Với tấm bị ngàm cả 4 cạnh, chịu tác động tải phân bố p(x,y), năng lượng tồn phần tấm tính như sau: 2 2 D⎪⎧⎛ ∂ 2 w ∂ 2 w ⎞ ⎡⎛ ∂ 2 w ⎞ ∂ 2 w ∂ 2 w⎤⎪⎫ Π = ⎨⎜ + ⎟ +2(1 −ν )⎢⎜ ⎟ − ⎥⎬dxdy − pw(,) x y dxdy 2 ∫∫ ⎜ ∂x 2 ∂y 2 ⎟ ⎜ ∂x ∂ y ⎟ ∂x 2 ∂y 2 ∫∫ ⎩⎪⎝ ⎠ ⎣⎢⎝ ⎠ ⎦⎥⎭⎪ Tùy thuộc điều kiện biên, hàm chuyển vị nên chọn hợp hồn cảnh. Những cơng thức sau đây cĩ thể giúp bạn đọc chọn lựa hàm hình dáng cho các phép tính gần đúng. w(,) x y= ∑∑ amn Xm( x ). Y n ( y ) ) ( 7 2 . 1 mn Bảng 1.1 Điều kiện biên Hàm hình dáng, ∑Xm(x) mπ x ∑sin m a 1 ⎛ 2mπ x ⎞ ∑ ⎜1− cos ⎟ m = 1,3,5,L m 2 ⎝ a ⎠ πx mπ x ∑sin sin m a a 2 2 x ⎛ x ⎞ m x ⎛ x ⎞ 1 mπ x −1+ − 1 −1 − sin ⎜ ⎟ ∑() 2 ⎜ ⎟ ∑ a ⎝ a ⎠ m a ⎝ a ⎠ m mπ a m 2 2 2 ∑(a− 4 x) x m = 0,1,2,3,L m x ⎛ x ⎞⎛ x ⎞ 1 mπ x ⎜ −1⎟⎜ −1⎟ − ∑ sin a ⎝ a ⎠⎝ 2a ⎠ m mπ a x ⎛ x 2 ⎞ 1 mπ x ⎜ −1⎟ − 1 m − sin ∑ ⎜ 2 ⎟() ∑ m 2a ⎝ a ⎠ m mπ a Ví dụ 3: Áp dụng phương pháp Ritz xác định mặt võng tấm thép hình chữ nhật cĩ các cạnh axb, chịu tác động lực tập trung theo phương pháp tuyến tại điểm (x1, y1), hình 1.8. Tấm tựa cả bốn cạnh, chiều dầy tấm t = const. Mơ đun đàn hồi vật liệu E. Để áp dụng phương pháp Ritz cần thiết tìm phiếm hàm cho phương trình chuyển vị dưới tác động lực. Phiếm hàm của bài tốn này chính là hàm năng lượng của tấm. Hàm chuyển vị mặt qua giữa tấm: mπ x nπ y Hàm f(x,y) hãy là: f( x , y )= sin sin mn a b
  20. 19 ∞ ∞ mπ x n π y Hàm chuyển vị: w(,) x y= ∑∑ amn sin sin mn=11= a b Điều kiện biên: Tại x = 0, x = a: w = 0; ∂w/∂x ≠ 0. Tại y = 0, y = b: w = 0; ∂w/∂y ≠ 0. Cơng biến dạng: 1 U= u dV ∫ 0 2 V ε với thế năng đơn vị: u = σT εd ε 0 ∫{}{} 0 Hình 1.8 Tấm tựa trên 4 cạnh 1 hay là u =()σεσε + + σε +2 τγ + 2 τγ + 2 τγ 0 2 x x y y x x xy xy yz yz zx zx ⎧ ⎡ 2 2 2 2 ⎤⎫ D ⎪ 2 2 ⎛ ∂ w ⎞ ∂ w ∂ w ⎪ U =⎨() ∇w +2(1 −ν )⎢⎜ ⎟ − ⎥⎬dxdy 2 ∫∫ ⎜ ∂x ∂ y ⎟ ∂x 2 ∂y 2 ⎩⎪ ⎣⎢⎝ ⎠ ⎦⎥⎭⎪ E t 3 trong đĩ D = 12(1−ν 2 ) Cơng ngoại lực: W= ∫∫ q( x,)(,) y w x y dxdy Thay w và các đạo hàm của w theo x, y vào biểu thức tính U: 2 Dab⎪⎧⎛ ∂ 2 w ∂ 2 w ⎞ ⎡⎛ ∂ 2 w ⎞ ∂ 2 w ∂ 2 w⎤⎪⎫ U = ⎨⎜ + ⎟ +2() 1 −ν ⎢⎜ ⎟ − ⎥⎬dxdy 2 ∫∫ ⎜ ∂x 2 ∂y 2 ⎟ ⎜ ∂x ∂ y ⎟ ∂x 2 ∂y 2 00⎩⎪⎝ ⎠ ⎣⎢⎝ ⎠ ⎦⎥⎭⎪ 2 2 ab⎧ 2 2 ⎫ D ⎪ ∞ ∞ ⎡⎛ mπ⎞ ⎛ n π ⎞ ⎤ mπ x n π y ⎪ U = a + sin sin dxdy ∫∫⎨∑∑ mn ⎢⎜ ⎟ ⎜ ⎟ ⎥ ⎬ 2 mn=11= ⎝ a ⎠ ⎝ b ⎠ a b 00⎩⎪ ⎣⎢ ⎦⎥ ⎭⎪ Cơng do lực tập trung tác động: ∞ ∞ mπ x n π y W= P∑∑ amn sin sin mn=11= a b Phiếm hàm I: I ≡ Π = U-W. ∂()UW− Tiến hành lấy đạo hàm của (U - W) theo amn: sẽ nhận được hệ phương trình đại số, ∂amn trong đĩ amn đĩng vai trị ẩn số. Sử dụng tính trực giao hàm lượng giác trong phạm vi 0 ≤ x ≤ a và 0 ≤ y ≤ b: a mπ x n π x ⎧ 0 m≠ n ∫ sin sin dx = ⎨ 0 a a ⎩a/ 2 m= n
  21. 20 a mπ x n π x ⎧ 0 m≠ n ∫ sin sin dx = ⎨ 0 b b ⎩b/ 2 m= n ⎡ 2 2 ⎤ Dab⎛ mπ⎞ ⎛ n π ⎞ mπ x0 n π y0 amn × ⎢⎜ ⎟ + ⎜ ⎟ ⎥ − Psin sin = 0 4 ⎣⎢⎝ a ⎠ ⎝ b ⎠ ⎦⎥ a b với m = 1,2, ; n = 1,2, Từ phương trình này xác định amn: 4P mπ x n π y a = sin0 sin 0 mn 2 2 ⎡⎛ mπ⎞ ⎛ n π ⎞ ⎤ a b abD⎢⎜ ⎟ + ⎜ ⎟ ⎥ ⎣⎢⎝ a ⎠ ⎝ b ⎠ ⎦⎥ Nghiệm bài tốn: mπ x nπ y sin0 sin 0 4P ∞ ∞ mπ x n π y w(,) x y = a b × sin sin ∑∑ 2 2 abD mn=11= ⎡⎛ mπ⎞ ⎛ n π ⎞ ⎤ a b ⎢⎜ ⎟ + ⎜ ⎟ ⎥ ⎣⎢⎝ a ⎠ ⎝ b ⎠ ⎦⎥ Ví dụ 3b: Giải bài tốn uốn tấm, tỷ lệ cạnh a/b = 1,5. Tấm bị ngàm 4 cạnh, chịu tác động tải phân bố đều theo phương pháp tuyến, hình 1.9. Hình 1.9 Tấm chữ nhật, ngàm 4 canh Sử dụng tính đối xứng cấu hình tấm, chọn gốc hệ tọa độ nằm chính giữa tấm, hai trục đi qua đường đối xứng tấm. Hàm thử trình bày dưới dạng: ∞∞ amn ⎡ m mπ x⎤⎡ n mπ y ⎤ w(, x,) y = ∑∑ ⎢1− (1) cos⎥⎢ 1− (1) cos ⎥ m,n = 1,3,5 . . . mn 4 ⎣ a ⎦⎣ a ⎦ Điều kiện biên:
  22. 21 ⎛ ∂w ⎞ ()w x=± a = 0, ⎜ ⎟ = 0 ⎝ ∂x ⎠ x=± a ⎛ ∂w ⎞ ()w x=± b = 0, ⎜ ⎟ = 0 ⎝ ∂y ⎠ x=± b Để nhanh chĩng nhận kết quả, cĩ thể độ chính xác cịn bàn thêm, hãy nhận m = n = 1: a ⎛ πx ⎞⎛ πy ⎞ w =11 ⎜1 + cos⎟⎜ 1+ cos ⎟ 4 ⎝ a ⎠⎝ b ⎠ Thay biểu thức này vào cơng thức tính cơng biến dạng tấm sẽ nhận được: +a+b 4 2 D 2 Dπ a ⎛ 3b 3a 2 ⎞ U = ∇ 2 w dxdy = 11 + + ∫∫() ⎜ 3 3 ⎟ 2 −a−b 32 ⎝ a b ab ⎠ Cơng ngoại lực thực hiện: +a+b W= p w(,) x y dxdy= p a ab 0 ∫∫ 0 11 −a−b Từ phép tính tìm minimum tổng năng lượng của tấm: ∂()UW − = 0 ∂a11 Xác định được: 4 16 p0 a 1 a11 = Dπ 4 3+ 3()a / b4 + 2 () a / b 2 Với a/b = 1,5 và υ = 0,3, giá trị độ võng tại x = y = 0 sẽ là: p a 4 w = 0,0791 0 max Et 3 Ví dụ 4: Ổn định tấm. Thứ tự thực hiện phép tính khi tính ổn định tấm trong khuơn khổ phương pháp Ritz bao gồm xác định hàm năng lượng của tấm trong trường hợp chịu nén trong mặt phẳng và các phép tính đạo hàm phương trình năng lượng theo chuyển vị suy rộng nhằm xác định lực tới hạn. Xác định lực giới hạn tác động tấm chữ nhật axb, tựa tự do tại x = 0, x = a và y = 0. Cạnh y = b tự do. Hàm chuyển vị của tấm: ∞ ∞ w(x,y) = ∑∑ amn ϕm(x) Ψn(y) m=1 n= 1 Hàm năng lượng của tấm gồm năng lượng tấm chịu uốn và năng lượng biến dạng mặt trung hịa tấm U = U1 + U2. Trong đĩ U1 tính theo cơng thức: ⎧ ⎡ 2 2 2 2 ⎤⎫ D ⎪ 2 2 ⎛ ∂ w ⎞ ∂ w ∂ w ⎪ U1 =⎨() ∇w +2(1 −ν )⎢⎜ ⎟ − ⎥⎬dxdy 2 ∫∫ ⎜ ∂x ∂ y ⎟ ∂x2 ∂y 2 ⎩⎪ ⎣⎢⎝ ⎠ ⎦⎥⎭⎪ U2 tính theo cơng thức:
  23. 22 2 2 1 ⎪⎧ ⎛ ∂w ⎞ ⎛ ∂w ⎞ ∂w ∂w⎪⎫ U 2 = ⎨N x ⎜ ⎟ + N y ⎜ ⎟ + 2N xy ⎬dxdy 2 ∫∫ ∂x ⎜ ∂y ⎟ ∂x ∂y ⎩⎪ ⎝ ⎠ ⎝ ⎠ ⎭⎪ trong đĩ: Et ⎛ ∂u ∂v ⎞ Et ⎛ ∂v ∂u ⎞ N = ⎜ 0 +ν 0 ⎟ và N = ⎜ 0 +ν 0 ⎟ x 2 ⎜ ⎟ y 2 ⎜ ⎟ 1−ν ⎝ ∂x ∂y ⎠ 1−ν ⎝ ∂y ∂x ⎠ ⎛ ∂u0∂v 0 ⎞ Nxy = Nyx = Gt⎜ + ⎟ ⎝ ∂y ∂x ⎠ Trường hợp Nx = -σxt; Ny = Nxy = 0, hàm U2 trở thành: 2 1 ab ⎛ ∂w ⎞ U = N dxdy 2 ∫∫ x ⎜ ⎟ 2 00 ⎝ ∂x ⎠ Điều kiện biên: ∂w tại x = 0 và x = a: w = 0; ≠ 0; ∂x ∂w Tại y = 0: w = 0; ≠ 0; ∂y ∂w Tại y = b: w ≠ 0; ; ≠ 0 ∂y Hàm w(x,y): ∞ ∞ n ⎛ y⎞ mxπ w(x,y) = ∑∑ amn ⎜ ⎟ sin m=1 n= 1 ⎝ b ⎠ a Để nhanh chĩng cĩ kết quả, chỉ chọn 1 hằng số khi tính chuỗi cho w. y πx w( x , y )= a sin 11 b a Thay giá trị vừa viết vào phương trình của U1, U2 hàm thế năng này trở thành: ⎧ ⎡ 2 2 2 2 ⎤⎫ D ⎪ 2 2 ⎛ ∂ w ⎞ ∂ w ∂ w ⎪ U1 =⎨() ∇w +2(1 −ν )⎢⎜ ⎟ − ⎥⎬dxdy = 2 ∫∫ ⎜ ∂x ∂ y ⎟ ∂x2 ∂y 2 ⎩⎪ ⎣⎢⎝ ⎠ ⎦⎥⎭⎪ 4 2 a2 D ⎡⎛ π ⎞ ab 1 ⎛ π ⎞ a⎤ = 11 ⎢⎜ ⎟ +2(1 −ν ) ⎜ ⎟ ⎥ 2⎣⎢⎝ a ⎠ 6 b⎝ a ⎠ 2⎦⎥ 2 1 2 ⎛ π ⎞ ab U2 + a11σ x t⎜ ⎟ 2 ⎝ a ⎠ 6 Tiến hành tính đạo hàm theo biểu thức: ∂ ()UU+ 1 2 = 0 ∂a11 sẽ nhận được phương trình:
  24. 23 4 2 2 ⎪⎧ ⎡⎛ π ⎞ ab 1 ⎛ π ⎞ a ⎤ ⎛ π ⎞ ab⎪⎫ a11 ⎨D⎢⎜ ⎟ +2(1 −ν ) ⎜ ⎟ ⎥ − σ xt⎜ ⎟ ⎬ ; 0 = ⎩⎪ ⎣⎢⎝ a ⎠ 6 b⎝ a ⎠ 2⎦⎥ ⎝ a ⎠ 6 ⎭⎪ Từ tính chất cơ học của bài tốn, hằng số a11 ≠ 0 biểu thức trong dấu ngoặc lớn phải bằng 0, từ đĩ xác định lực Euler: 2 2 D ⎡⎛ π ⎞ 6 ⎤ π 2 D ⎡6(1−ν ) ⎛ b ⎞ ⎤ σx= σ E = ⎢⎜ ⎟ +(1 −ν ) 2 ⎥ = 2 ⎢ 2 + ⎜ ⎟ ⎥ t⎣⎢⎝ a ⎠ a ⎦⎥ b t ⎣⎢ π ⎝ a ⎠ ⎦⎥ Phương pháp Rayleigh-Ritz giải bài tốn dao động tấm Ví dụ a: Áp dụng phương pháp Rayleigh – Ritz xác định tần số riêng tấm hình chữ nhật cạnh axb, chiều dày t = const, bốn mép bị ngàm. Thế năng tấm chịu uốn được tính như sau: 2 2 1 ⎪⎧⎛ ∂ 2 w ∂ 2 w ⎞ ⎡⎛ ∂ 2 w ⎞ ∂ 2 w ∂ 2 w⎤⎪⎫ UD0 = ⎨⎜ + ⎟ +2() 1 −υ ⎢⎜ ⎟ − ⎥⎬dxdy , 2 ∫∫ ⎜ ∂x 2 ∂y 2 ⎟ ⎜ ∂x ∂ y ⎟ ∂x 2 ∂y 2 ⎩⎪⎝ ⎠ ⎣⎢⎝ ⎠ ⎦⎥⎭⎪ Et 3 trong đĩ D = - độ cứng tấm, với t – chiều dày tấm. 12() 1−υ 2 Khối lượng qui đổi tấm tính bằng biểu thức: m = ∫∫ ρtw2 dxdy , trong đĩ ρ - khối lượng của đơn vị diện tích tấm. ⎛ 2πx ⎞⎛ 2πy ⎞ Hàm chuyển vị tìm theo cách làm của Ritz: w= C⎜1 + cos ⎟⎜1+ cos ⎟ ⎝ a ⎠⎝ b ⎠ Thế năng và động năng tấm tham gia dao động: 2 4 ⎛ 3 2 3 ⎞ 9 2 U0 = C. D .2π ab⎜ + + ⎟; M = Cρ tab ⎝ a4 a 2 b 2 b 4 ⎠ 4 Điều kiện để thế năng và động năng tấm đạt cực trị: ∂ 1 2 ∂ U 0 − ω M =0;k = 1,2, ∂Ck 2 ∂Ck 2U 2U 4 a 2 a 4 D ω 2 = 0 hay là ω =0 =π 2 3 + 2 + 3 × M M 3 b 2 b 4 ρta 4 D Trường hợp a = b cĩ thể thấy rằng: ω = 37,2 ρta 4 Ví dụ b: Xác định hai tần số thấp nhất dầm tựa trên trên hai gối đơn, tại hai đầu theo phương pháp Rayleigh-Ritz. Hàm chuyển vị thành lập dạng chuỗi hữu hạn sau, nhằm thỏa mãn điều kiện biên: 2 3 x ⎛ x ⎞ ⎛ x ⎞ w() x= a0 + a 1 + a2 ⎜ ⎟ + a3 ⎜ ⎟ ( ) a L ⎝ L ⎠ ⎝ L ⎠
  25. 24 Điều kiện biên: w(0) = w(L) = 0 (b) Cĩ thể viết lại phương trình (a) dưới dạng: 2 3 ⎡⎛ x ⎞ x ⎤ ⎡⎛ x ⎞ x ⎤ w( x) = A1 ⎢⎜ ⎟ − ⎥ + A2 ⎢⎜ ⎟ − ⎥ (a’) ⎣⎢⎝ L ⎠ L⎦⎥ ⎣⎢⎝ L ⎠ L⎦⎥ Thành phần ma trận khối lượng mij và ma trận cứng kij, i, j = 1, 2: 2 L ⎛ x 2 x ⎞ mL L ⎛ x 2 x ⎞⎛ x 3 x ⎞ mL m= m⎜ − ⎟ dx = ; m= m⎜ − ⎟⎜ − ⎟dx = ; 11 ∫ ⎜ 2 ⎟ 12 ∫ ⎜ 2 ⎟⎜ 3 ⎟ 0 ⎝ L L ⎠ 30 0 ⎝ L L ⎠⎝ L L ⎠ 20 2 2 L ⎛ x3 x ⎞ 8mL L ⎡ d 2 ⎛ x 2 x ⎞⎤ 4EJ m= m⎜ − ⎟ dx = ; m = m k= EJ ⎜ − ⎟ dx = 22 ∫ ⎜ 3 ⎟ 21 12; 11 ∫ ⎢ 2 ⎜ 2 ⎟⎥ 3 0 ⎝ L L ⎠ 105 0 ⎣dx ⎝ L L ⎠⎦ L 2 2 L d 2 ⎛ x 2 x ⎞ d 2 ⎛ x3 x ⎞ 6EJ L ⎡ d 2 ⎛ x 3 x ⎞⎤ 12EJ k= EJ ⎜ − ⎟ ⎜ − ⎟dx = ; k= EJ ⎜ − ⎟ dx = 12 ∫ 2 ⎜ 2 ⎟ 2 ⎜ 3 ⎟ 3 22 ∫ ⎢ 2 ⎜ 3 ⎟⎥ 3 0 dx ⎝ L L ⎠ dx ⎝ L L ⎠ L 0 ⎣dx ⎝ L L ⎠⎦ L k21 = k12 EJ ⎡4 6 ⎤ ⎡ 1 1 ⎤ −ω 2 mL 30 20 = 0 3 ⎢ ⎥ ⎢ 1 8 ⎥ L ⎣6 12⎦ ⎣ 20 105 ⎦ Phương trình đặc trưng từ hệ phương trình trên đây cĩ dạng: 2 ⎛ β⎞⎛ 8 β⎞ ⎛ β ⎞ ω 2mL 4 ⎜4 − ⎟⎜12 − ⎟ −⎜6 − ⎟ = 0; β = ⎝ 30 ⎠⎝ 105 ⎠ ⎝ 20 ⎠ EJ Từ đây cĩ thể nhận được: EJ ω 2 = 120 ; 1 mL4 EJ ω 2 = 2520 ; 2 mL4 2 CÁC PHƯƠNG PHÁP NHĨM TRỌNG HÀM DƯ Phương pháp hàm trọng lượng dư, gọi cách khác trọng hàm dư (weighted residual method) thích hợp giải gần đúng các bài tốn phương trình vi phân tuyến tính và/hoặc khơng tuyến tính. Sử dụng phương pháp này khơng qua giai đoạn tìm phiếm hàm. Bài tốn kỹ thuật suy rộng cĩ thể viết dưới dạng phương trình nêu tại (a) phần mở đầu: L(u) - p = 0 trong miền V, và các điều kiện biên: B(u) -q = 0 trên biên S Hàm u dưới dạng hàm xấp xỉ: N ~ u= ∑ ai f i i=1 ~ ~ Hàm u được tìm phải thỏa mãn điều kiện B(u ) = q(s).
  26. 25 ~ Với số lượng i hữu hạn, khơng chắc thỏa mãn L (u ) = p trong miền V. Sai số trong trường hợp này cĩ thể biểu diễn dưới dạng hàm các sai số, ký hiệu R, mở đầu từ Residual function, hoặc ký tự mở đầu từ error ε: N R ≡ ε = L ( ∑ai f i ) - p ≠ 0. (1.28) i=1 Yêu cầu tính tốn là phải xây dựng trọng hàm w, để tích w.f(R), là hàm của R, thỏa mãn những điều kiện minimum hoặc một tiêu chuẩn “nhỏ nhất cĩ thể” theo nghĩa cụ thể, hồn cảnh cụ thể. Hàm f(R) được chọn nhằm đạt f(R) = 0 khi R = 0, điều này xảy ra khi u~ tiến đến u. Hàm u~ trong mọi trường hợp phải thỏa mãn điều kiện biên, và tiêu chuẩn “nhỏ nhất cĩ thể”: wf( R ) dV = 0 (1.29) ∫V 2.1 PHƯƠNG PHÁP GALERKIN Phương pháp Galerkin thích hợp cho những bài tốn cơ học thuộc lý thuyết đàn hồi, biến dạng dẻo, lý thuyết trường và cho các phương trình vi phân. N ~ ~ Hàm u dưới dạng hàm xấp xỉ: u= ∑ ai f i , hàm sai số hay hàm dư ε = L(u ) - p i=1 Cách làm của Galerkin là đặt ε trực giao với hệ hàm thử fi, từ đĩ cĩ thể nhận tích vơ hướng sau: = 0, i =1,2, (1.30) Dưới dạng đầy đủ cơng thức cuối được viết lại: ⎧ ⎛ ⎞ ⎫ ⎜ a f⎟ − p f dV = 0 , i =1,2, (1.31) ∫ ⎨L ∑ i i ⎬ i ⎩ ⎝ i ⎠ ⎭ Ví dụ 5: Giải phương trình vi phân: d2 u L ()u− p = +u + x = 0 trong miền 0 < x <1 dx 2 với điều kiện biên: u = 0 tại x =0 u = 0 tại x =1. ~ ~ Hàm xấp xỉ u được tìm dưới dạng: u = x(1-x)( a1 + a2x + ) (a) Hàm u~ thoả mãn các điều kiện biên đặt ra trước đĩ. Nếu chỉ hạn chế chuỗi (a) với 2 thành phần của a1 và a2, hàm u sẽ là: ~ u = x( 1-x) (a1 + a2x) (b) Hàm sai số: ~ 2 2 3 ε = L (u ) - p = x( -2 +x -x )a1 + (2 - 6x + x - x )a2 Thay các hàm u~ và ε vào phương trình Galerkin (1.31), chúng ta nhận được hệ phương trình đại số: 1 ∫εx()1− x dx = 0 0 1 ∫εx2 (1− x ) dx = 0 (c) 0
  27. 26 Giải hai tích phân trên, hệ phương trình đại số sẽ cĩ dạng: ⎡ 3 3 ⎤ ⎧ 1 ⎫ ⎢10 20 ⎥⎧a1 ⎫ ⎪12 ⎪ ⎢ ⎥⎨ ⎬ = ⎨ ⎬ (d) 3 13 a 1 ⎢ ⎥⎩ 2 ⎭ ⎪ ⎪ ⎣20 105⎦ ⎩20⎭ Sau khi giải hệ phương trình, các hệ số được xác định như sau: 71 7 a = và a = 1 369 2 41 Từ đĩ hàm xấp xỉ u sẽ là: 71 7 u~ = x(1- x)( + x ) (e) 369 41 Phương pháp Galerkin suy rộng, sử dụng rộng rãi như một phép biến phân. Theo hướng này cĩ thể viết hệ hàm trực giao hệ thống hàm fi như sau: ∫ { L (u) - p} fidV = 0, i=1,2,3, ,N (1.32) Mặt khác biến phân của u cĩ thể khai triển dạng chuỗi: δu = δa1f1 + δa2f2 + + δaNfN (1.33) δu = ∑δaifi Thay biểu thức cuối vào (1.31) sẽ nhận được hệ phương trình: ∫ { L (u) - p} δu dV = 0 (1.34) Ví dụ 6: Xác định độ võng dầm như đã nêu tại ví dụ trình bày tại ví dụ 1, phần trước. Phương trình vi phân miêu tả độ võng dầm cĩ dạng: d 2 ⎛ d2 w ⎞ ⎜ EI ⎟ −q = 0 0 ≤ x ≤ L 2 ⎜ 2 ⎟ dx ⎝ dx ⎠ Điều kiện biên: tại x = 0 và x= L: w =0; w’’ = 0; Độ võng được tìm theo biểu thức (d) như đã trình bày trong phương pháp Ritz: Hãy chọn lời giải thử nghiệm hai thành phần sau: πx 3πx w~( x )= a sin + a sin 1 L 2 L 4 4 ⎛ π⎞ πx ⎛ 3π⎞ 3 πx R= EIa1 ⎜ ⎟ sin + EIa2 ⎜ ⎟ sin − q ⎝ L ⎠ L ⎝ L ⎠ L Thủ tục của phương pháp Galerkin đưa đến hệ phương trình đại số: L πx L 2L ⎫ sin Rdx= EIa1 −q = 0 ∫0 L 2 π ⎪ ⎬ L 3πx L 2L sin Rdx= EIa2 −q = 0⎪ ∫0 L 2 3π ⎭ Từ hệ phương trình cĩ thể xác định: 4qL4 4qL4 a = a = 1 π 5 EI 2 243π 5 EI
  28. 27 Trường hợp các hệ số vơ cùng lớn, hàm w(x) viết dưới dạng chung: ∞ (2i− 1)π x w( x )= ∑ ai sin i=1 L Áp dụng cơng thức Galerkin vào bài tốn sẽ nhận được hệ phương trình: L ⎡ d4 w() x ⎤ (2i− 1)π x EI − q sin dx = 0 , i =1,2, ,∞ ∫ ⎢ 4 ⎥ 0 ⎣ dx ⎦ L Sau khi giải hệ phương trình sẽ xác định các giá trị của hệ số ai, từ đĩ cĩ thể viết biểu thức cho hàm độ võng dầm. 4qL4 ∞ 1 kπ x w() x = sin 5 ∑ 5 π EIk =1,3, k L Kết quả tính theo phương pháp Galerkin trong trường hợp này trùng với kết quả tính theo phương pháp Ritz. Ví dụ 7: Áp dụng cơng thức suy rộng của Galerkin (1.34) giải phương trình Poisson trong miền hạn chế bằng hình chữ nhật cạnh axb. Phương trình Poisson: ∂ 2u ∂ 2u + = C (a) ∂x 2 ∂y 2 Điều kiện biên: tại x = 0 và x = a: u =0; tại y = 0 và y = b: u =0. (b) Hàm xấp xỉ cĩ thể viết dưới dạng: u(x,y) = αx( x-a)(y-b)y (c) Thay (c) vào (1.9) chúng ta nhận được: ab⎛ ∂ 2u ∂ 2u ⎞ ⎜ + − C ⎟ δudxdy = 0 (d) ∫∫⎜ 2 2 ⎟ 00⎝ ∂x ∂y ⎠ trong đĩ: δu = δα x(x-a)y(y-b) (e) Sau khi thay thế (c) vào (d) sẽ được: α Ca3 b 3 []a3 b 3() a 2+ b 2 − = 0. 90 36 5 C Từ đĩ: α = (f) 2 a2+ b 2 và hàm u: 5C u = (x2 - ax) (y2 - by) (g) a2+ b 2 Ví dụ 8: Áp dụng cơng thức Galerkin (1.34) giải bài tốn uốn tấm thép cĩ các cạnh axb, ngàm tại các mép tấm, chịu tải phân bố đều p(x,y) = const tác động pháp tuyến. Độ cứng tấm D, chiều dày tấm t. Phương trình độ võng w(x,y) trong miền hạn chế nêu trên được tìm dưới dạng hàm xấp xỉ:
  29. 28 ∞ u(x,y) = w(x,y) = ∑ cifi(x,y) (a) i=1 Mỗi hàm fi, i=1, 2, , phải thỏa mãn điều kiện biên và phương trình vi phân: D.∇2 ∇2 u - p(x,y) = 0 (b) trong đĩ p(x,y) tải trọng tác động lên tấm. Biểu thức (a) chọn cho u, trong trường hợp này cĩ thể là: ∞ ∞ 1 2mxπ 2nyπ u(x,y) = amn. (1-cos ) ( 1- cos ) (c) ∑∑ 4 a b m=1 n= 1 Chuỗi (c) thoả mãn điều kiện biên: tại x = 0; x = a: u = ∂u/ ∂x = 0; tại y= 0; y = b: u = ∂u/ ∂x = 0; (d) Nếu chỉ chọn m =1 và n =1, hàm u(x,y) sẽ cĩ dạng: 1 2πx 2πy u(x,y) = a11. (1-cos )(1-cos ) (e) 4 a b Sau khi thay thế u từ (e) vào phương trình Galerkin (1.34): ∫[D.∇2 ∇2 u - p(x,y)] δudxdy = 0 Kết quả sẽ nhận được: a b 2 2 ∫∫ [D.a11.∇ ∇ fi(x,y) - p(x,y)] fi(x,y)dxdy = 0. (f) 0 0 Thực hiện phép lấy đạo hàm riêng 4 2 2 4π 2π 2π 2π 2π ∇ ∇ f1(x,y) = [ (1-cos y)cos x +cos x cos y +(1-cos a2 b 2 b a a b 2π 2π x)cos y ] (g) a b và thay thế (g) vào (f), hệ số đầu tiên của chuỗi cĩ dạng: pa2 b 2 a = 11 8Dπ 4 Với trường hợp tấm hình vuơng cạnh a, giá trị u tính tại tâm tấm là lớn nhất: pa 4 umax = 0,00128 D Kết quả tính cĩ thể so với nghiệm “chính xác” tính bằng phương pháp giải tích umax = pa 4 0,00126 . D 2.2 PHƯƠNG PHÁP TÍNH THỎA MÃN TRÊN MỘT SỐ ĐIỂM HOẶC MIỀN CHỌN LỰA N ~ ~ Hàm u dưới dạng hàm xấp xỉ: u= ∑ ai f i , hàm sai số hay hàm dư ε = L (u ) - p i=1 ~ ~ Hàm u được tìm phải thỏa mãn điều kiện Bk(u )= q(s), k =1,2,
  30. 29 Trong trường hợp bài tốn cĩ lời giải chính xác u~ = u, phương trình (a) và (b) tại phần mở đầu sẽ đúng: L (u) = p trong V và B(u) = q trên S. Trường hợp u~ = u giá trị ε→ 0, cịn trong các trường hợp khác biểu thức của ε được hiểu theo cách ε ≠ 0. Phương pháp chọn điểm tính tốn cố gắng đưa ε đến gần 0 theo nghĩa bình quân, hiểu theo ý sau: ~ = ∫ { L (u ) - p} wi dV = 0, i =1,2, tại một số điểm hoặc miền, được chọn trước. Ngồi các điểm được chọn những điều kiện đặt ra trong bài tốn khơng chắc được thỏa mãn. Thứ tự giải bài tốn như sau: 1. Tìm hàm xấp xỉ thỏa mãn điều kiện biên: N ~ u= ∑ ai f i i=1 2. Xác lập hàm sai số: ~ Ν ε = L ( u ) - p = ∑ ai L (fi) - p. i=1 3. Tại những điểm được chọn tiến hành ép buộc hàm ε = 0. ~ ∫{ L ( u ) - p }δi dV = 0, i =1,2, M (1.32) trong đĩ hàm Dirac δ được định nghĩa: ξ+c ∫ δ(x - ξ)dx = 1 với c → 0, ngồi phạm vi trên hàm δ(x) = 0. (1.33) ξ−c Ví dụ 9: Xác định nghiệm phương trình vi phân nêu tại ví dụ 5: d2 u L ()u− p = +u + x = 0 trong miền 0 < x <1 dx 2 với điều kiện biên: u = 0 tại x = 0 và u = 0 tại x =1. Hàm u được thử nghiệm như sau: ~ u (x) = x(1-x)(a1 + a2x ) Chọn điểm tính tại x = 0,25 và x = 0,5. Hàm sai số ε: ~ 2 2 3 ε = L ( u ) - p = x + ( -2 + x - x )a1 + ( 2- 6x + x - x )a2. Tại x =0,25 và x =0,5 các hệ số a1, a2 phải thỏa mãn: ⎡29 − 35⎤ ⎧1⎫ ⎢16 64 ⎥⎧a1 ⎫ ⎪4⎪ ⎢ ⎥⎨ ⎬ = ⎨ ⎬ 7 7 a 1 ⎢ ⎥⎩ 2 ⎭ ⎪ ⎪ ⎣ 4 8 ⎦ ⎩2⎭ 6 40 Từ đĩ: a1 = và a2 = 31 217 Hàm u~ (x) tính theo phương pháp chọn điểm:
  31. 30 1 u~ (x) = x(x-1)(42+40x). 217 Ví dụ 10: Xác định độ võng dầm liên tục dài L, EI = const, chịu tác động tải trọng phân bố đều q = const. Dầm ngàm bên trái, gối tự do bên đầu phía phải. Phương trình vi phân xác định độ võng dầm: d 2 ⎛ d2 w ⎞ ⎜ EI ⎟ −q = 0 (a) 2 ⎜ 2 ⎟ dx ⎝ dx ⎠ Điều kiện biên: tại x = 0: w = w’ = 0; tại x = L: w = w’’ =0; (b) Các điểm tính được chọn tại x = 0,25L; 0,5L; 0,75L. Lần thử đầu tiên sẽ chọn hàm w(x) theo chuỗi: N k ⎛ x ⎞ w(x) = ak ⎜ ⎟ (c) ∑ ⎝ L⎠ k =0 Hàm (c) trên đây chưa thỏa mãn điều kiện cân bằng phương trình vi phân (a) cùng các điều kiện biên. Do vậy từ (c) chỉ chọn những thành phần nhằm thoả mãn 4 biểu thức của các điều kiện biên đã nêu, cịn tại ba điểm chọn cho tính tốn, phả thỏa mãn điều kiện cân bằng phương trình vi phân bậc bốn. Hàm w(x) phải là: 2 3 4 5 6 ⎛ x ⎞ ⎛ x ⎞ ⎛ x ⎞ ⎛ x ⎞ ⎛ x ⎞ ⎛ x ⎞ w() x= a0 + a 1 ⎜ ⎟ + a2 ⎜ ⎟ + a3 ⎜ ⎟ + a4 ⎜ ⎟ + a5 ⎜ ⎟ + a6 ⎜ ⎟ (d) ⎝ L ⎠ ⎝ L ⎠ ⎝ L ⎠ ⎝ L ⎠ ⎝ L ⎠ ⎝ L ⎠ Sau khi thoả mãn điều kiện biên, các hệ số ai được xác định từ hệ phương trình: a0 = 0 a1 = 0 a0 + a1 + a2 + a3 + a4 + a5 + a6 = 0 2a2 + 6a3 + 12a4 + 20a5 + 30a6 = 0 Tại các điểm được chọn x = 0,25L; 0,5L; 0,75L, sau khi thỏa mãn phương trình (8.23), hệ phương trình chứa ai sẽ là: pL4 24a4 + 30a5 + 22,5a6 = EJ pL4 24a4 + 60a5 + 90a6 = EJ pL4 24a4 + 90a5 + 202,5a6 = EJ pL 5pL4 pL4 Từ đĩ: a0 = 0; a1 = 0; a2 = ; a3 = - ; a4 = ; a5 = 0; a6 = 0. 16EJ 48EJ 24EJ Ký hiệu ξ = x/L phương trình độ võng dầm theo cách tính trên viết lại như sau: 4 pL ⎡ 1 2 5 3 1 4 ⎤ w() x = ⎢ ξ− ξ + ξ EJ ⎣16 48 24 ⎦⎥
  32. 31 2.3 PHƯƠNG PHÁP TỔNG NHỎ NHẤT CÁC BÌNH PHƯƠNG CỦA SAI SỐ Phương pháp tốn cĩ tên gọi Least Squares sử dụng rất rộng rãi trong các phương pháp tính. Phương pháp trọng hàm coi đây cũng là cách làm hữu hiệu khi giải bài tốn phương trình vi phân. Trong khuơn khổ phương pháp này coi rằng tích phân của bình phương hàm sai số với trọng hàm dư trong miền khảo sát đạt minimum. wR2 dV = min ) 4 3 . 1 ( ∫V 2 ⎧ ⎛ ⎞ ⎫ wL ⎜ ai f i ⎟ − p dV = min ) 5 3 . 1 ( ∫V ⎨ ∑ ⎬ ⎩ ⎝ i ⎠ ⎭ Điều kiện cần để hàm mục tiêu đạt minimum: 2 ∂ ⎡ ⎧ ⎛ ⎞ ⎫ ⎤ ⎢ w⎨L ⎜ ai f i ⎟ − p⎬ ⎥ dV=0 i = 1,2,L , n (1.36) ∂a ∫V ∑ i ⎣⎢ ⎩ ⎝ i ⎠ ⎭ ⎦⎥ Hay là ⎧ ⎛ ⎞ ⎫ wLL() fi ⎨ ⎜ ai f i ⎟ − p dV=0 i = 1,2, , n (1.37) ∫V ∑ ⎬ L ⎩ ⎝ i ⎠ ⎭ Trong phương pháp này trọng hàm nhận bằng đơn vị. Ví dụ : Giải bài tốn dầm tựa hai đầu chịu tải trọng phân bố đều đang nêu bằng phương pháp tổng nhỏ nhất các bình phương. Sử dụng hàm thử, hàm sai số như đã nêu tại phương pháp Galerkin: πx 3πx w~( x )= a sin + a sin 1 L 2 L 4 4 ⎛ π⎞ πx ⎛ 3π⎞ 3 πx R= EIa1 ⎜ ⎟ sin + EIa2 ⎜ ⎟ sin − q ⎝ L ⎠ L ⎝ L ⎠ L Hệ phương trình đại số hình thành từ điều kiện cần: 4 4 L ⎧ L ∂ ⎛ 2 ⎞ ∂ ⎪ 2 2 ⎛ π⎞ 2 πx 2 2 ⎛ 3π⎞ 2 3 πx 2 ⎜ R dx⎟ = ⎨ [()EI a1 ⎜ ⎟ sin dx+ () EI a2 ⎜ ⎟ sin dx+ q + ⎝ ∫0 ⎠ ∫0 ∂a1 ∂a1 ⎩⎪ ⎝ L ⎠ L ⎝ L ⎠ L 2 4 4 ⎫ 2 ⎛ π⎞ πx 3 πx ⎛ π⎞ πx ⎛ 3π⎞ 3 πx⎤ ⎪ 2()EI a1 a 2 ⎜ ⎟ sin sin − 2EIqa1 ⎜ ⎟ sin− 2EIqa2 ⎜ ⎟ sin ⎥dx⎬ = 0 ⎝ L ⎠ L L ⎝ L ⎠ L ⎝ L ⎠ L ⎦⎥ ⎭⎪ 4 4 L ⎧ L ∂ ⎛ 2 ⎞ ∂ ⎪ 2 2 ⎛ π⎞ 2 πx 2 2 ⎛ 3π⎞ 2 3 πx 2 ⎜ R dx⎟ = ⎨ [()EI a1 ⎜ ⎟ sin dx+ () EI a2 ⎜ ⎟ sin dx+ q + ⎝ ∫0 ⎠ ∫0 ∂a1 ∂a1 ⎩⎪ ⎝ L ⎠ L ⎝ L ⎠ L 2 4 4 ⎫ 2 ⎛ π⎞ πx 3 πx ⎛ π⎞ πx ⎛ 3π⎞ 3 πx⎤ ⎪ 2()EI a1 a 2 ⎜ ⎟ sin sin − 2EIqa1 ⎜ ⎟ sin− 2EIqa2 ⎜ ⎟ sin ⎥dx⎬ = 0 ⎝ L ⎠ L L ⎝ L ⎠ L ⎝ L ⎠ L ⎦⎥ ⎭⎪ Hoặc 8 4 2 ⎛ π ⎞ L ⎛ π ⎞ ()EI a1 ⎜ ⎟ L− 4 EIq ⎜ ⎟ = 0 ⎝ L ⎠ π ⎝ L ⎠
  33. 32 8 8 L ⎧ L ∂ ⎛ 2 ⎞ ∂ ⎪ 2 2 ⎛ π⎞ 2 πx 2 2 ⎛ 3π⎞ 2 3 πx 2 ⎜ R dx⎟ = ⎨ [()EI a1 ⎜ ⎟ sin dx+ () EI a2 ⎜ ⎟ sin dx+ q + ⎝∫0 ⎠ ∫0 ∂a2 ∂a2 ⎩⎪ ⎝ L ⎠ L ⎝ L ⎠ L 8 4 4 ⎫ 2 ⎛ π⎞ πx 3 πx ⎛ π⎞ πx ⎛ 3π⎞ 3 πx⎤ ⎪ 2()EI a1 a 2 ⎜ ⎟ sin sin − 2EIqa1 ⎜ ⎟ sin− 2EIqa2 ⎜ ⎟ sin ⎥dx⎬ = 0 ⎝ L ⎠ L L ⎝ L ⎠ L ⎝ L ⎠ L ⎦⎥ ⎭⎪ Hoặc 8 4 2 ⎛ π ⎞ L ⎛ 3π ⎞ ()EI a2 ⎜ ⎟ L− 4 EIq ⎜ ⎟ = 0 ⎝ L ⎠ 3π ⎝ L ⎠ Từ hệ phương trình cĩ thể xác định: 4qL4 4qL4 a = a = 1 π 5 EI 2 243π 5 EI
  34. 33 Chương 2 PHƯƠNG PHÁP SAI PHÂN HỮU HẠN Phương pháp này cịn cĩ tên gọi thơng dụng, dễ nhớ là phương pháp lưới. Phương pháp được chọn dùng khi giải phương trình vi phân thơng thường và sau đĩ dùng cho các phương trình vi phân đạo hàm riêng. Phương pháp dùng cho bài tốn một chiều bằng cách chia đoạn thẳng đang xem xét thành nhiều phân đoạn hay là bước, gọi là sai phân, với tổng số bước hữu hạn. Trong bài tốn hai chiều cần tiến hành chia miền đang xét thành lưới với mắt lưới nằm trên các nút chọn lựa. Khơng gian ba chiều cũng được chia làm lưới theo phương pháp tương tự. 1 HÀM MỘT BIẾN Với các bài tốn một chiều, hàm liên tục y = y(x), với giá trị đầu y(x0) = y0, cĩ thể xác định giá trị của hàm tại vị trí x + Δx dựa vào chuỗi Taylor, dùng trong trường hợp tiến: 2 3 y(x + Δx) = y(x) + Δx y’(x) + Δx y’’(x) + Δx y’’’(x) + (2.1) 1! 2! 3! Khi tính thụt lùi chuỗi Taylor cĩ dạng: 2 3 y(x - Δx) = y(x) - Δx y’(x) + Δx y’’(x) - Δx y’’’(x) + (2.2) 1! 2! 3! Nếu chỉ lấy hai thành phần đầu của chuỗi, khi trừ (2.2) với (2.1) chúng ta nhận được: y(x + Δx) - y(x - Δx) = 1 [y’(x)Δx - (-y’(x)Δx) ] 2Δx Từ đĩ cơng thức tính đạo hàm bậc 1, tại x sẽ là: y’(x) = 1 . [y(x + Δx) - y(x - Δx) ] (2.3) 2Δx Nếu lấy ba thành phần đầu của chuỗi để tính, phép tính đạo hàm bậc hai sẽ như sau: y(x + Δx) - y(x - Δx) = 2y(x) + y’’(x). (Δx)2 Từ đĩ: y’’(x) = 1 .[ y(x + Δx) - 2y(x) + y(x - Δx) ] (2.4) ()Δx 2 Cơng thức rút ra từ cách làm trên, áp dụng khi tính đạo hàm bậc từ 1 đến 4 cĩ dạng sau: 1 y ' = ()y− y k 2Δx k+1 k − 1 '' 1 yk = ()yk +1 −2 yk + y k −1 ()Δx 2 ''' 1 yk = (yk+2 - 2yk+1 - 2yk-1 - yk-2 ) 2()Δx 3 (IV )' 1 yk = ( yk+2 - 4yk+1 +6yk - 4yk-1 + yk-2) (2.5) ()Δx 4 Hình 2.1 Hàm u(x) và sai phân
  35. 34 Ứng dung phương pháp sai phân hữu hạn giải phương trình vi phân Sử dụng các biểu thức ghi trong (2.5) giải bài tốn phương trình vi phân thường gặp trong cơ học kết cấu. Dầm thẳng, momen quán tính mặt cắt thay đổi theo luật I(x) = I0 (1 + x/L), với I0 là momen quán tính mặt cắt đầu dầm phía trái. Hãy tính độ võng dầm tại vị trí ¼ chiều dài dầm, tính từ trái. Phương trình vi phân uốn dầm: d2 y EI( x) = −M (*) dx 2 trong đĩ M – momen uốn dầm. Momen uốn dầm tính như sau: p M =0 x() L − x 2 Điều kiện biên: Hình 2.2 Dầmthẳng chịutácđộng tải phân bố đều M(0) = M(L) = 0 Phương trình vi phân (*) cĩ thể viết lại dưới dang: d2 y p L2 ⎡ x ⎛ x ⎞ ⎛ x ⎞⎤ p L2 ξ(1− ξ ) 0 0 2 = − ⎢ ⎜1− ⎟ /⎜ 1+ ⎟⎥ ≡ − dx 2EI 0 ⎣ L ⎝ L ⎠ ⎝ L ⎠⎦ 2EI 0 ()1+ ξ trong đĩ ξ = x/L Điều kiện biên: y(0) = y(L) = 0 Cĩ thể dùng biến khơng thứ nguyên u thay vào vị trí của y: 4 p0 L u= y /(2EI ) 0 Phương trình vi phân giờ cĩ dạng: d u''= −ξ()() 1 − ξ / 1 + ξ u'= dξ cùng điều kiện biên: u(0) = u(1) = 0 Sử dụng các cơng thức tính sai phân phần trên, trong đĩ h = Δx: '' 1 y = (yk+1 -2yk + yk-1) k h 2 cho các phép tính tiếp theo. Trường hợp 1: chia dầm làm 2 đoạn bằng nhau ( n = 2) h =1/ 2 = 0,5 ; 1 u'' =()() u −2 u + u =5 0 − 2u − 0 = − 8u 1 5,0 2 0 1 2 1 1 '' Với ξ = 0,5 và u1 = −8 u1 −8u1 = − 0,5(1 − 0,5) /(1 + 0,5) từ đây tính được u1 = 0,020833 So với lời giải chính xác 0,017601, kết quả đang nêu sai số 18,4%. Trường hợp 2: chia dầm làm 4 đoạn, n = 4; h = 0,25. '' 2 Tại k =1: u1 = (1/ 0,25)()u0− 2 u 1 + u 2 = −0,25( 1 − 0,25) /( 1 + 0,25) = − 0,15
  36. 35 '' 2 Tại k =2: u2 = (1/ 0,25)()u1− 2 u 2 + u 3 = −0,5( 1 − 0,5) /( 1 + 0,5) = − 0,166667 '' 2 Tại k =3: u3 = (1/ 0,25)()u2− 2 u 3 + u 4 = −0,75( 1 − 0,75) /( 1 + 0,75) = − 0,107143 Lập hệ phương trình đại số từ kết quả vừa cĩ: ⎡− 2 1 0 ⎤⎧u1 ⎫ ⎧ 9,375 ⎫ ⎪ ⎪ ⎪ ⎪ ⎢ 1− 2 1 ⎥ u = −10−3 10,417 ⎢ ⎥⎨ 2 ⎬ ⎨ ⎬ ⎪ ⎪ ⎪ ⎪ ⎣⎢ 0 1− 2⎦⎥⎩u3 ⎭ ⎩ 6,697 ⎭ Kết quả tính: {}u = []0,013917 0,018452 0,012567 T Sai số của u2 so với lời giải chính xác 4,8%. Trường hợp n = 8; h = 1/8 Sau khi tính nhận được u4 = 0,17817, sai số 1,2% Trường hợp chung, dầm chia làm N đoạn, h = 1/N " 1 uk = ()uk 1 −2 uk + u k 1 ()1/ N 2 − + Ký hiệu ξk = 1/N, k = 1, 2, . . . , N cĩ thể viết lại phương trình đang nêu: " uk= −ξ k()()1 − ξ k / 1 + ξ k Từ hai cơng thức này hình thành quan hệ: −1()k / N( 1− k / N ) −1 k( N− k) u−2 u + u = = k −1 k k +1 N 2 1+ k/ N N 3 N+ k Dưới dạng ma trận phương trình này cĩ thể viết như sau: ⎡− 2 1 ⎤⎧ u1 ⎫ ⎧ 1(NN− 1)( / + 1 ) ⎫ ⎪ ⎢ 1− 2 1 0 ⎥⎪ u ⎪ 2NN− 2 / + 2 ⎪ ⎢ ⎥⎪ 2 ⎪ ⎪ ()()⎪ ⎢ 1− 2 1 ⎥⎪ u3 ⎪ ⎪ 3()()NN− 3 / + 3 ⎪ ⎢ ⎥⎪ ⎪ ⎪ ⎪ ⎢ O ⎥⎪ M ⎪ −1 ⎪ M ⎪ ⎨ ⎬ = 3 ⎨ ⎬ ⎢ O ⎥ M N M ⎢ ⎥⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎢ 1− 2 1 ⎥ u N −3 ()NNN−3 3/() + () − 3 ⎪ ⎪ ⎪ ⎪ ⎢ 0 1− 2 1 ⎥ u ()NNN−2 2 /() + () − 2 ⎢ ⎥⎪ N −2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎣⎢ 1− 2⎦⎥⎩u N −1 ⎭ ⎩ ()NNN−1 /() + () − 1 ⎭ Sau giải hệ phương trình sẽ xác định các thành phần vecto {u}. Để ý đến ma trận vuơng đang đĩng vai trị ma trận cứng cĩ thể nhận xét, đây là ma trận ba đường chéo và là ma trận đối xứng. Tính chất rất quí này của ma trận cho phép chúng ta sử dụng các thuật tốn xử lý ma trận ba đường chéo chính khi giải hệ phương trình. Điều này cịn dẫn đến khả năng chương trình hĩa cách giải các bài tốn cơ học kết cấu trong khuơn khổ phương pháp sai phân hữu hạn.
  37. 36 Hình 2.3 Ba trường hợp tính tốn: (a) n =2, (b) n = 4, (c) n = 8 Ví dụ 2.1: Xác định độ võng dầm dài L, độ cứng EJ = const, chịu tác động lực phân bố q(x) hình tam giác, p0 tại đầu bên trái, p = 0 tại đầu phía phải, hình 2.4. Hình 2.4 Dầm ngàm trái, tựa phải chịu tải phân bố tam giác Chia dầm thành 5 đoạn thẳng, chiều dài mỗi đoạn h = L/5. Phương trình vi phân chính yếu cĩ dạng: d 2 ⎛ d2 w ⎞ ⎛ x ⎞ EJ ⎜ ⎟ =q( x ) = p 1 − trong miền 0 < x < L (a) 2 ⎜ 2 ⎟ 0 ⎜ ⎟ dx ⎝ dx ⎠ ⎝ L ⎠ Điều kiện biên:
  38. 37 w(0) = 0; w’(0) = 0; w(L) = 0 và M(L) = 0. (b) Tại những điễm chọn làm nút tính tốn, biểu thức (a) trở thành: IV EJwj = p j ( ) c Cơng thức xác định đạo hàm bậc 4 ghi tại (b) phần trên. (IV) 1 w (x) = (wk+2 - 4wk+1 +6wk - 4w k-1 + wk-2 )( ) d h 4 Giá trị pj đọc tại nút thứ j. Từ (c) và (d) tiến hành xây dựng hệ phương trình đại số: 4 j = 1: w-1 - 4w0 + 6w1 - 4w2 + w3 = (4/5)(p0h )/(EJ) 4 j = 2: w0 - 4w1 + 6w2 - 4w3 + w4 = (3/5)(p0h )/(EJ) 4 j = 3: w1 - 4w2 + 6w3 - 4w4 + w5 = (2/5)(p0h )/(EJ) 4 j = 4: w2 - 4w3 + 6w4 + 4w5 = w6 = (1/5)(p0h )/(EI) (e) Từ điều kiện biên cĩ thể viết w0 = w5 = 0. Chuyển vị nút ảo w-1 xác định từ w’0 = 0. Biểu thức tính như sau: 1 w' =() −w + w 0 2h −1 1 Từ đĩ w-1 = w1. " 1 Từ w =0 =()w −2 w + w xác định w6 = -w4. 5 h 2 4 5 6 Điều kiện trên đây miêu tả tại hình 2.5 Hệ phương trình đại số viết tại (e) sau thay điều kiện biên: ⎡ 7− 4 1 0 ⎤⎧w1 ⎫ ⎧4⎫ ⎢ ⎥⎪ ⎪ ⎪ ⎪ 4 −4 6 − 4 1 ⎪w2 ⎪ ⎪3⎪ p h ⎢ ⎥⎨ ⎬ = ⎨ ⎬ 0 ⎢ ⎥ 1− 4 6 − 4 ⎪w3 ⎪ ⎪2⎪ 5 EI ⎢ ⎥⎪ ⎪ ⎪ ⎪ ⎣ 0 1− 4 5 ⎦⎩w4 ⎭ ⎩1⎭ Lời giải phương trình: ⎧w1 ⎫ ⎧0,001243⎫ ⎪ ⎪ ⎪ ⎪ 4 ⎪w2 ⎪ ⎪0,002553⎪ p0 L Hình 2.5 ⎨ ⎬ = ⎨ ⎬ ⎪w3 ⎪ ⎪0,002793⎪ EI ⎪ ⎪ ⎪ ⎪ ⎩w4 ⎭ ⎩0,001788⎭ Momen uốn dầm và lực cắt tính theo cơng thức: M() x= − EIw "() x F () x= − EIw "'() x Diễn đạt theo cách làm trong phương pháp sai phân hữu hạn, hai cơng thức vừa nêu được biết là: 1 M= − EI ()w−2 w + w i h 2 i−1 i i+1 1 F= − EI () −w +2 w − 2 w + w i 2h 2 i−2 i−1 i+1 i + 2 2 Giá trị tính tốn M/(p0L ) và F/(p0L) như sau:
  39. 38 B ảng 2.1 x/L 0 0.4 0,6 1,0 2 M/(p0L ) -0,06215 0,02675 0,3113 0 F/(p0L) 0,40238 0,082 -0,0179 -0,0979 Tương tự cách làm này tính độ võng dầm thẳng dài L, độ cứng EJ = const, chịu tác động lực phân bố q(x) = const sẽ thực hiện theo cách sau. Với các nút k = 1,2,3,4,5 phương trình cân bằng cĩ dạng: qh 4 wk+2 - 4wk+1 +6wk - 4wk-1 + wk-2 = EJ Phương pháp xử lý w0 và w5 theo cách đang nêu. 2 PHƯƠNG PHÁP LƯỚI CHO BÀI TỐN HAI CHIỀU Động tác cần tiến hành đầu tiên là tạo lưới cho mặt hai chiều. Lưới được chia ở dạng đơn giản nhất gồm hai hệ đường thẳng trực giao, hệ thứ nhất song song với trục Ox, hệ kia song song với Oy. Bước của lưới gồm: Δx - dọc trục Ox và Δy - dọc trục Oy. Chỉ số biến thiên vị trí dọc Ox ký hiệu ∂ m+ n bằng i, cịn dọc Oy ký hiệu k. Theo qui ước này, cơng thức tính đạo hàm riêng u sẽ như sau: ∂xm ∂ y n Hình 2.6 Phương pháp lưới Tại điểm x = xk, y = yi, đạo hàm riêng bậc từ 1 đến 4 cĩ dạng: ∂u u− u ∂u u− u ik = i, k+ 1 i , k − 1 ; ik = i+1, k i− 1, k )( 6 . 2 ∂x 2Δx ∂y 2Δy 2 2 ∂ uuuik,,+−11− 2 ki+ ik ∂ ui+1, k− 2 u ik+ u i−1, k uik = ; uik = (2.7) ∂x 2 ()Δx 2 ∂y 2 ()Δy 2 2 ∂ ui+1, k + 1− u i + 1, k − 1− u i − 1, k + 1+ u i − 1, k − 1 uik = )( 8 . 2 ∂∂xy 4ΔΔx y 4 ∂ ui, k + 2− 4 ui , k + 1 + 6 uik− 4 u i, k− 1+ u i , k − 2 uik = )( 9 . 2 ∂x 4 ()Δx 4 4 ∂ ui+2, k− 4 u i+1, k+ 6 u ik− 4 u i−1, k+ u i−2, k uik = 0() 1 . 2 ∂y 4 ()Δy 4
  40. 39 ∂ 4 1 uik = [ui+1, k + 1+ u i + 1, k − 1 + u i − 1, k + 1 + u i − 1, k − 1 + ∂∂xy22 ()ΔΔx y 2 + 4uik − 2( u i, k+ 1+ u i , k− 1+ ) u i + 1, k+ u i−1, k 1)] 1( . 2 Các biểu thức vi phân của momen uốn và lực cắt: ⎡ Δ2 w Δ2 w ⎤ ()mx i, k ≈ − D⎢ 2 +υ 2 ⎥ ⎣()Δx ()Δy ⎦ i, k ⎡ Δ2 w Δ2 w ⎤ m≈ − D +υ ()y i, k ⎢ 2 2 ⎥ ⎣()Δy ()Δx ⎦ i, k ⎡ Δ ⎛ Δw ⎞⎤ ()mxy = () mxy ≈ −()1 −υ D⎢ ⎜ ⎟⎥ 2() 1 . 2 i, k i, k 2Δx ⎜ 2Δy ⎟ ⎣ ()⎝ ⎠⎦ i, k và ⎪⎧ Δ ⎡ Δ2 w Δ2 w ⎤⎪⎫ ()fx ≈ − D⎨ ⎢ + ⎥⎬ i, k ⎪2()Δx Δx 2 Δy 2 ⎪ ⎩ ⎣() ()⎦⎭i, k ⎪⎧ Δ ⎡ Δ2 w Δ2 w ⎤⎪⎫ ()fy ≈ − D⎨ ⎢ + ⎥⎬ ) ( 3 1 . 2 i, k ⎪2()Δy Δx 2 Δy 2 ⎪ ⎩ ⎣() ()⎦⎭i, k Thứ tự giải bài tốn uốn tấm thể hiện qua bài tốn xác định độ võng tấm hình vuơng cạnh a, chịu tải trọng phân bố đều q, tác động theo phương pháp tuyến. Cạnh tấm x = ± a tựa trên gối, dọc y = 2 ± a cạnh bị ngàm. 2 Phương trình chung uốn tấm: ∂ 4 w ∂ 4 w ∂ 4 w q + 2 + = ∂x 4 ∂x2 ∂ y 2 ∂y 4 D trong đĩ D - độ cứng tấm. Điều kiện biên: a ∂ 2 w tại x = ± : w = 0; = 0; 2 ∂x 2 a ∂w tại y = ± : ; w = 0; 0 = 2 ∂y Tấm được chia thành lưới 4x4 với Δx = Δy = a/4. a a Hàm w cần thỏa mãn w = 0 cho tất cả các nút nằm trên cạnh tấm x = ± và y = ± . 2 2 Hệ phương trình đại số để xác định hàm chuyển vị w(x,y) theo phương pháp sai phân hữu hạn: ⎛ ∂ 4 w ∂ 4 w ∂ 4 w ⎞ Δx 4 ⎜ + 2 + ⎟ = (u - 4u + 6u - 4u , + u ) ⎜ 4 2 2 4 ⎟ i+2,k i+1, k ik i-1 k i-2,k ⎝ ∂x ∂x ∂ y ∂y ⎠ik +2β(u i+1, k+1 + u i+1,k-1 - 2ui+1,k + 4uik -2ui,k-1 - 2ui-1,k + ui-1,k+1 + ui-1,k-1) +
  41. 40 4 2 qik Δ x + β (ui,k+2 - 4ui,k+1 + 6uik - 4ui,k-1 + ui,k-2 ) = D Δx 2 trong đĩ: β = ; Δy 2 Trong trường hợp này β =1; qik = q = const. Thay thế các giá trị Δx = Δ; y = a/4 vào phương trình cuối, thực hiện phép tính cho tất cả 4 nút, hệ phương trình đại số theo phương pháp lưới cĩ dạng: 10w1 - 16w2 - 16w3 + 8 w4 = C -8w1 + 21w2 + 4w3 - 16w4 + w6 = C -8w1 + 4w2 + 21w3 - 16w4 + w8 = C 2w1 -8w2 - 8w3 + 22w4 + w5 + w9 = C qa 4 trong đĩ C = 256D ∂ 2 w Điều kiện biên = 0 trên x = a/2 được thỏa mãn trong hệ phương trình: ∂x 2 w2 + w6 = 0 w4 + w5 = 0 w7 = 0 ∂w Cịn điều kiện = 0 đưa đến: ∂y w3 - w8 = 0 w4 - w9 =0 w10 = 0 Giải hệ phương trình đại số trên đây cho phép xác định hàm u(x,y) tại tất cả các nút trên lưới. Ứng suất trong tấm tính theo cơng thức: ∂ 2 u−2 u + u σ =u = i+1, k ik i−1, k x ∂y 2 ()Δy 2 ∂ 2 u−2 u + u σ =u = i, k + 1 ik i, k − 1 ; y ∂x 2 ()Δx 2 ∂ 2 u− u− u+ u τ = − u = i+1, k + 1 i + 1, k − 1 i − 1, k + 1 i − 1, k − 1 xy ∂x∂ y 4ΔΔx y Ví dụ 2.2: Tấm thép hình chữ nhật cạnh dài a, cạnh ngắn ¾a. Tấm bị ngàm một mép chiều dài a, hai mép cận kề tựa gối cứng, mép cịn lại tự do. Tải tác động pháp tuyến cĩ giá trị khơng đổi. Sử dụng lưới bước ¼a xác định chuyển vị w(x,y) tấm. Chia tấm thành 12 hình vuơng cạnh ¼a. Các nút ảo được đánh số trên cùng hình.
  42. 41 Hình 2.7 Phương trình chung uốn tấm: ∂ 4 w ∂ 4 w ∂ 4 w q + 2 + = (a) ∂x 4 ∂x2 ∂ y 2 ∂y 4 D trong đĩ D - độ cứng tấm. Điều kiện biên: Tại mép ngàm y = 0: w = 0; w5 = w6 = w7 = w8 = w9 = 0 (b) ∂w/∂y = 0: w13 = w1; w14 = w2; w15 = w3. (c) 2 2 ∂ w/∂x = 0 dọc trục x: w8 = w10 = 0 (d) Tại mép tựa x = 0 và y = ¾a: w = 0: w12 = w20 = w27 = w28 = w29 = w30 = w31 = 0 2 2 ∂ w/∂x = 0: w13 = -w11 ; w21 = -w19 2 2 và ∂ w/∂y = 0: w21 = w33; w22 = -w34; w23 = -w35, w24 = -w36 (e) 2 2 ∂ w/∂x = 0 dọc mép y =3a/4 : w30 = w32 = 0 Tại mép tự do x = 0: Nz = 0 tại nút 16 và 24. w18 – w14 + (ν-2)(w8 – w25 + w23 – w10) + (6-ν)(w15 – w17) = 0 w26 – w22 + (ν-2)(w15 – w30 – w17 – w32) + (6-ν)(w23 – w25) = 0 Mz = 0 tại nút 16, 24: w15 + w17 + ν(w9 + w24) – 2(1+ν)w16 = 0 w23 + w25 + ν(w16 + w31) – 2(1+ν)w24 = 0 Sau xử lý điều kiện biên, hệ phương trình đại số cịn lại 12 ẩn, tính từ w13 đến w26. Phương trình ma trận cĩ dạng:
  43. 42 ⎡20− 8 1 0 0 0− 8 2 0 0 0 0⎤⎧w13 ⎫ ⎧1⎫ ⎢− 8 21− 8 1 0 0 2− 8 2 0 0 0⎥⎪w ⎪ ⎪1⎪ ⎢ ⎥⎪ 14 ⎪ ⎪ ⎪ ⎢ 1− 8 21− 8 1 0 0 2− 8 2 0 0⎥⎪w15 ⎪ ⎪1⎪ ⎢ ⎥⎪ ⎪ ⎪ ⎪ ⎢ 0 1− 8 21− 8 1 0 0 2− 8 2 0⎥⎪w16 ⎪ ⎪1⎪ ⎢− 8 2 0 0 0 0 18− 8 1 0 0 0⎥⎪w ⎪ ⎪1⎪ ⎪ 17 ⎪ ⎢ ⎥ 4 ⎪ ⎪ ⎢ 2− 8 2 0 0 0− 8 19− 8 1 0 0⎥⎪w18 ⎪ ph ⎪1⎪ ⎨ ⎬ = ⎨ ⎬ ⎢ 0 2− 8 2 0 0 1− 8 19− 8 1 0⎥ w D 1 ⎢ ⎥⎪ 21 ⎪ ⎪ ⎪ ⎢ 0 0 2− 8 2 0 0 1− 8 19− 8 0⎥⎪w ⎪ ⎪1⎪ ⎪ 22 ⎪ ⎪ ⎪ ⎢ 0− 15,4 0− 5,41 0 0− 1,7 0 1,7 0⎥ w 0 ⎢ ⎥⎪ 23 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎢ 0 0− 1,7 0 1,7 0 0− 15,4 0− 5,41⎥ w24 0 ⎢ ⎥⎪ ⎪ ⎪ ⎪ w ⎢ 00 1− 2,6 1 000 0 0,3 0 0⎥⎪ 25 ⎪ ⎪0⎪ ⎢ ⎥⎪ ⎪ ⎪ ⎪ ⎣ 00 0 0,3 0 000 1− 2,6 1 0⎦⎩w26 ⎭ ⎩0⎭ Nghiệm phương trình mang dạng sau, đơn vị đo ph4/D. w13 = 0,25619 w14 = 0,38943 w15 = 0,45037 w16 = 0,51951 w17 = 0,70839 w18 = 1,07433 w21 = 0,3038 w22 = 0,46598 w23 = 0,54558 w24 = 0,63989 w25 = 0,96226 w25 = 2,27756 3 XOẮN DẦM Bài tốn xoắn dầm trụ mặt cắt hình chữ nhật, trong phần “Lý thuyết đàn hồi” xác định dạng: ∂ 2ψ∂ 2 ψ + = −2ψ = 0 tại biên ∂x2 ∂ y 2 Từ cơng thức này chúng ta sẽ xác định hàm ứng suất Prandtl ψ . Ví dụ 2.3 Xác định ứng suất do xoắn dầm mặt cắt hình vuơng, cạnh a. (a) (b) Hình 2.8 Mặt cắt hình vuơng và lưới Lưới mặt cắt thể hiện tại hình 2.8, theo đĩ Δx = Δy = 1/3. Vị trí các nút ψj,k j, k = 1,2,3,4 như thể hiện tại cùng hình. Điều kiện biên của hàm Prandtl như sau: ψ j 1, = 0 ψ ,1 k = 0 và ψ j 4, = 0 ψ ,4 k = 0 j = 1,2,3,4; k = 1,2,3,4 Từ các cơng thứctính đạo hàm riêng: 2 2 ∂ uuuik,,+−11−+2 ik ik ∂ ui+ ,1 k− 2 u ik+ u i− ,1 k uik = ; uik = ∂x 2 ()Δx 2 ∂y 2 ()Δy 2
  44. 43 và điều kiện biên cĩ thể xây dựng hệ phương trình đại số: −4ψ2,2 + ψ 2,3 + ψ 3,2 = −0,22222⎫ ⎪ ψ2,2−4 ψ 2,3 + ψ 3,3 = −0,22222 ⎪ ⎬ ψ2,2−4 ψ 3,2 + ψ 3,3 = −0,22222 ⎪ ⎪ ψ2,3+ ψ 3,2 −4 ψ 3,3 = − 0,22222 ⎭ Lời giải hệ phương trình: {}ψ = []0,11111 0,11111 0,11111 0,11111 T Trường hợp dùng lưới mịn hơn, hình phía phải, hệ phương trình đại số cĩ dạng sau: ⎡− 4 1 ⎤⎧ψ 2,2 ⎫ ⎪ ⎪ ⎢ 1 4 1 ⎥ ψ ⎢ − ⎥⎪ 2,3 ⎪ ⎢ 1− 4 1 0 ⎥⎪ψ 2,4 ⎪ ⎢ ⎥⎪ ⎪ ψ ⎢ 1− 4 1 ⎥⎪ 3,2 ⎪ ⎢ ⎥⎪ ⎪ 1− 4 1 ⎨ψ 3,3 ⎬ ={} − 0,125 ⎢ ⎥ 1− 4 1 ⎪ψ ⎪ ⎢ ⎥⎪ 3,4 ⎪ ⎢ 0 1− 4 1 ⎥⎪ψ ⎪ ⎢ ⎥ 4,2 ⎪ψ ⎪ ⎢ 1− 4 1 ⎥⎪ 4,3 ⎪ ⎢ ⎥ ⎪ ⎣ 1− 4⎦⎩⎪ψ 4,4 ⎭ ⎧ψ 2,2 ⎫ ⎧ 0,0859 ⎫ ⎪ ⎪ ψ ⎪ 0,10937 ⎪ ⎪ 2,3 ⎪ ⎪ ⎪ ⎪ψ 2,4 ⎪ ⎪ 0,0859 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ψ 3,2 ⎪ ⎪ 0,10937 ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ψ 3,3 ⎬ = ⎨0,140625⎬ ⎪ψ ⎪ ⎪ 0,10937 ⎪ ⎪ 3,4 ⎪ ⎪ ⎪ ⎪ψ 4,2 ⎪ ⎪ 0,0859 ⎪ ⎪ψ ⎪ ⎪ ⎪ ⎪ 4,3 ⎪ ⎪ 0,10937 ⎪ ⎪ ⎪ ⎩⎪ψ 4,4 ⎭⎪ ⎩ 0,0859 ⎭ Ứng suất cắt do xoắn xác định từ biểu thức: ∂ψ Gφ ∂ψ Gφ τ= G φ = ()ψ− ψ τ= −G φ = − ()ψ− ψ zx ∂y 2Δ 3,4 3,2 zy ∂x 2Δ 4,3 2,3 4 BÀI TỐN TRƯỜNG 2D VỚI BIÊN CONG Trong trường hợp biên cong cần thiết hàm hĩa biên dạng đa thức: 2 2 wxy(,) = w0 + axayax 1 + 2 + 3 + ay4 ) ( 5 1 . 2 Từ hình 2.9 cĩ thể tính các hằng: 2 2 2 2 h()() w4− w 0 + h 4 w0 − w 2 h( w1− w 0) + h 1 ( w0 − w 3 ) a1 = ; a2 = (2.16) hh4 () h+ h4 hh1() h+ h 1
  45. 44 h()() w4− w 0 + h 4 w 0 − w 2 h( w1− w 0) + h 1( w 0 − w 3 ) a3 = ; a4 = hh4 () h+ h4 hh1() h+ h 1 Hình 2.9 Biên cong ⎛ ∂ 2w ⎞ ⎛ ∂ 2 w ⎞ Tại điểm 0 (x =0; y = 0): ⎜ ⎟ = 2a ; ⎜ ⎟ = 2a ) 7 1( . 2 ⎜ 2 ⎟ 3 ⎜ 2 ⎟ 4 ⎝ ∂x ⎠0 ⎝ ∂y ⎠0 Từ đĩ cĩ thể viết: ⎛ ∂ 2 w ∂ 2w ⎞ 2w 2w 2w 2w ⎛ 2 2 ⎞ h2 ⎜ + ⎟ ≈ 1 + 2 + 3 + 4 −⎜ + ⎟w (2.18) ⎜ 2 2 ⎟ ⎜ ⎟ 0 ⎝ ∂x ∂y ⎠ α1()1+ α 1 ()1+ α4 ()1+ α1 α4 ()1+ α4 ⎝ α1 α 4 ⎠ h trong đĩ α =i ,i = 1,4 i h Trường hợp hi của cả 4 nhánh đều nhỏ so với h, cơng thức trên cĩ dạng: 2w 2w 2w 2w h2()∇ 2 w ≈ 1 + 2 + 3 + 4 ααα1() 1+ 3 ααα2() 2+ 4 ααα3() 1+ 3 ααα4() 2+ 4 (2.19) ⎛ 2 2 ⎞ ⎜ ⎟ −⎜ + ⎟w0 ⎝α1 α 3 α 2 α 4 ⎠ Ví dụ 2.4: Xác định độ võng và momen uốn tấm tại tâm tấm thép hình ellip, tựa tại mép tấm, chịu tác động tải phân bố đều p. Kích thước cho trước: a = 0,15 m, b = 0,1 m, t = 0,05 m. Sử dụng tính chất đối xứng qua cả hai trục của cấu hình tấm, chỉ đưa ¼ tấm vào tính tốn.Theo cách ký hiệu trên hình 2.10. khoảng cách h1 = 0,044 m cịn h2 = 0,0245 m, h3 =0,03 m. Tại nút 1, 2, 3 và 4 sẽ áp dụng cơng thức sai phân chuẩn, cịn nút 5 và 6 áp dụng mơ hình trình bày phía phải hình 2.9. Để ý rằng tại biên tấm w = 0 và M = 0. Từ biểu thức ∇4 w = p/ D , sau biến đổi cĩ thể viết: ∂ 2 M ∂ 2 M ∂ 2 w ∂ 2 w M + = − p;+ = − ∂x 2 ∂y 2 ∂x 2 ∂y 2 D Dưới dạng ma trận biểu thức cuối viết như sau: Hình 2.10
  46. 45 ⎡− 4 2 0 2 0 0 ⎤⎧M 1 ⎫ ⎢ 1− 4 1 0 2 0 ⎥⎪M ⎪ ⎢ ⎥⎪ 2 ⎪ ⎢ 0 1− 4 0 0 2 ⎥⎪M 3 ⎪ 2 ⎢ ⎥⎨ ⎬ = − ph ⎢ 1 0 0− 4 2 0 ⎥⎪M 4 ⎪ ⎢ 0 1,06 0 1− 4,26 1 ⎥⎪M ⎪ ⎢ ⎥⎪ 5 ⎪ ⎣⎢ 0 0 1,34 0 1,25− 7,41⎦⎥⎩⎪M 6 ⎭⎪ Giải hệ phương trình sẽ nhận được các giá trị sau cho vecto M: {}M= ph 2 []1,384 1,23 0,769 1,038 0,884 0,423 T Và ph 4 {}w = []1,52 1,282 0,647 1,066 0,853 0,325 T D Momen uốn tính tại vị trí giữa tấm: D M =[]2w − 2 w +υ() 2 w − 2 w = 0,187 pb 2 x h 2 1 2 1 4 D M =[]2w − 2 w +υ() 2 w − 2 w = 0,263pb2 = M y h 2 1 4 1 2 max 5 PHƯƠNG PHÁP SAI PHÂN HỮU HẠN TRÊN CƠ SỞ PHÉP BIẾN PHÂN Phương pháp sai phân hữu hạn đang xem xét áp dụng hữu hiệu cho những kết cấu cĩ cấu hình chuẩn, hay nĩi rõ là cấu hình đơn giản. Ma trận cứng lập ra trong những trường hợp này thuộc dạng ba đường chéo chính hoặc năm đường chéo chính, và là ma trận đối xứng. Với những kết cấu cĩ cấu hình phức tạp, tính đối xứng của ma trận khơng cịn, những ưu việt vốn cĩ của phương pháp khơng phát huy tác dụng như mong đợi. Ngày nay các nhà nghiên cứu đề xuất thủ tục xử lý phương pháp sai phân hữu hạn trên cơ sở phép biến phân kinh điển. Tên gọi phương pháp cĩ nguồn gốc từ các cơng trình liên quan, đăng tải trên các tạp chí tốn và cơ học, trong những năm cuối sáu mươi, đầu bảy mươi thế kỷ XX. Đây là phương pháp sai phân xây dựng trên nền của nguyên lý năng lượng (finite difference energy method). Từ nguyên lý cơng ảo đề cập tại sách “Cơ học kết cấu tàu thủy” thấy rõ rằng xây dựng bài tốn theo nguyên lý này đảm bảo ma trận trong hệ phương trình đại số luơn đối xứng và luơn dương. Trong chừng mức cĩ thể nhận xét, đây là sự chuyển hĩa những ưu việt tiềm tàng của phương pháp phần tử hữu hạn vào phương pháp sai phân hữu hạn. Nĩi cách khác, giữa hai phương pháp này cĩ sự tương đồng. Tranh thủ tính tương đồng này phương pháp sai phân hữu hạn đ ang tìm cách đổi mới. Những ví dụ tiếp theo trình bày cách làm này. Ví dụ 2.5: Giải bài tốn uốn dầm tại ví dụ 2.1, hinh 2.2, trên cơ sở của nguyên lý năng lượng. Từ nguyên lý cơng ảo cĩ thể viết cơng ảo hệ thống như sau đây: LL δΠ = ∫∫EIw""δ w dx− pδ wdx = 0 (a) 00 Tiến hành chia chiều dài dầm tại hình 2.2 thành 6 phân đoạn. Phương trình năng lượng của dầm xét như tổng năng lượng sáu thành phần thành viên:
  47. 46 5 L δΠ = EIw""δ w− p δ w dx = 0 ( ) b ∑∫ ()i i i i=0 0 Đại lượng w” và EI xem là khơng đổi trong phạm vi phân đoạn, L0 và L5 hãy là L/10, cịn L1 đến L4 là L/5. Lấy tích phân theo cơng thức cuối sẽ nhận được: ⎛ L ⎞ ⎛ L ⎞ ⎛ L ⎞ δΠ = ⎜ ⎟ δW0 + ⎜ ⎟()δWWWW1+ δ 2 + δ3 + δ 4 + ⎜ ⎟δW5 = 0 (c) ⎝10 ⎠ ⎝ 5 ⎠ ⎝10 ⎠ " " trong đĩ δWi= δ w i () EI wi− p iδ w i " 1 Thay cơng thức tính đạo hàm bậc hai: wi = (wi−1 −2 wi + w i+1 ) vào δΠ sẽ nhận được quan hệ: ()Δx 2 EI δΠi =() δw i 1 −2 δ wi + δ w i 1 (wi 1 −2 wi + w i 1 ) − pδ wi (d) − + ()L / 5 2 − + Điều kiện biên: ' w0= w 0 = w5 =0 w−1 = w 1 Từ đĩ: 3 δΠ = [EI/() L / 5]{δ w1[]() 7 w 1− 4 w 2 + w 3 − (4 / 5)p0 ( L / 5) +δw2 []() −4 w1 + 6 w 2 − 4 w 3 + w 4 − (3/ 5)p0 ( L / 5) +δw[]() w −4 w + 6 w − 4 w − ( 2 / 5)p( L / 5) 3 1 2 3 4 0 (e) +δw4[]() w 2 −4 w 3 + 6 w 4 + w 6 − (1/ 5) p0 ( L / 5) +δw6[]() w 4 + w 6 = 0 Dùng ký hiệu ma trận cĩ thể viết phương trình (e) dạng sau: δΠ = δwT () Kw− P = 0 ( ) f Để ý đến điều kiện w4 + w6 = 0 hay là w6 = -w4 cĩ thể viết: ⎡ 7− 4 1 0 0⎤⎧w1 ⎫ ⎧4 / 5⎫ ⎢−4 6 − 4 1 0⎥⎪w ⎪ ⎪3/ 5⎪ ⎢ ⎥⎪ 2 ⎪ ⎪ ⎪ 4 ⎪ ⎪ ⎪ ⎪ p0 () L / 5 ⎢ 1− 4 6 − 4 0⎥⎨w3 ⎬ = ⎨2 / 5⎬ ( ) g ⎢ ⎥ EI 0 1− 4 6 1 ⎪w ⎪ ⎪1/ 5⎪ ⎢ ⎥⎪ 4 ⎪ ⎪ ⎪ ⎢ ⎥ ⎣ 0 0 0 1 1⎦⎩⎪w5 ⎭⎪ ⎩⎪ 0 ⎭⎪ Ma trận cứng của hệ thống là ma trận năm đường chéo, cĩ tính đối xứng và dương. Giải hệ phương trình cho phép xác định chuyển vị các các điểm tính. Kết quả tính trùng với kết quả nhận từ ví dụ trước. Dầm đang đề cập xét dạng chung như tập hợp của nhiều phân đoạn, giới hạn giữa các phân đoạn là nút. Phương trình cân bằng xây dựng cho mỗi phân đoạn, như đã gọi là sai phân, được viết lại như sau:
  48. 47 ⎧ 1 ⎫ ⎧wi−1 ⎫ ⎪ ⎪ EI ⎪ ⎪ δΠi =() δw i−1 −2 δ wi + δ w i+1 ⎨ − 2⎬ 4 []1− 2 1 ⎨ wi ⎬ ⎪ ⎪ h ⎪ ⎪ ⎩ 1 ⎭ ⎩wi+1 ⎭ ⎡ 1− 2 1 ⎤ EI = δw iT ⎢−2 4 − 2⎥wi= δ w i k i w i ( ) h h 4 ⎢ ⎥ ⎣⎢ 1− 2 1 ⎦⎥ trong đĩ wi là vec to chuyển vị các nút, ki là ma trận cứng lập cho phân đoạn. Ma trận cứng vừa lập cĩ tính đối xứng. Ví dụ 2.6: Xác định ứng suất cắt do xoắn dầm mặt cắt hình vuơng nêu tại ví dụ 2.3 Hình 2.11 Chia lưới cho ví dụ 2.6 Phương trình Poisson viết cho trường hợp này: ∂ 2ψ∂ 2 ψ + = −2 ψ = 0 tại biên ∂x2 ∂ y 2 Sử dụng nguyên lý cơng bù, nguyên lý này trình bày tại sách “Cơ học kết cấu tàu thủy”, cĩ thể viết biểu thức xác định cơng bù ảo: * ⎛ ∂ψ ∂ ∂ψ ∂ ⎞ δΠ = ∫⎜ δψ + δψ− 2 δψ ⎟dxdy = 0 ( ) a ⎝ ∂y∂ y∂ x∂ x ⎠ Hãy tạo lưới cho bài tốn 2 chiều này như đề nghị tại hình 2.11. Mỗi tam giác vuơng cân trên mặt cắt dầm qui ước mang tên gọi “phân tử”, tên gọi mươn từ phương pháp phần tử hữu hạn. Thứ tự các nút của mỗi “phần tử” tam giác giới thiệu tại cùng hình 2.11. Tích phân () thực hiện cho tồn mặt cắt. Đại lượng δψ cho mỗi tam giác nhận bằng giá trị trung bình cọng của δψ tính cho 3 nút tam giác đĩ. Các cơng thức tính đạo hàm riêng. Tính cho tam giác bên trái: ∂ψ ψ− ψ ∂ψ ψ− ψ = 2 1 = 3 2 ∂x h∂ y h ∂ δψ− δψ ∂ δψ− δψ δψ = 2 1 δψ = 3 2 ∂x h∂ y h Tính cho tam giác phía phải
  49. 48 ∂ψ ψ− ψ ∂ψ ψ− ψ = 3 2 = 2 1 ∂x h∂ y h ∂ δψ− δψ ∂ δψ− δψ δψ = 3 2 δψ = 2 1 ∂x h ∂y h Cơng bù ảo của dầm bằng tổng cơng bù ảo các phần tử tạo thành: * ⎛ ∂ψ ∂ ∂ψ ∂ ⎞ δΠ = ∑ ∫⎜ δψ + δψ− 2 δψ ⎟dxdy = 0 ( ) b k ⎝ ∂y∂ y∂ x∂ x ⎠ Hay là ⎡⎛ψk− ψ k⎞⎛ ψ k− ψ k ⎞ ⎛ψk− ψ k⎞⎛ ψ k− ψ k ⎞ 2 ⎤ δ * A ⎜ 2 1⎟⎜ 2 1 ⎟ ⎜ 3 2⎟⎜ 3 2 ⎟ +δψk + δψk + δψ k Π = ∑ ⎢⎜ ⎟⎜ ⎟ + ⎜ ⎟⎜ ⎟ ()1 2 3 ⎥ (b’) k ⎣⎢⎝ h ⎠⎝ h ⎠ ⎝ h ⎠⎝ h ⎠ 3 ⎦⎥ k trong đĩ A – diện tích tam giác thứ k, ψ i ,i = 1,2, 3 là giá trị ψ tại ba nút tại ba gĩc tam giác. Biểu thức này nên viết dưới dạng ma trận: δΠ* = ∑ δψ kT(k kψ k−p k ) =0 ( ) c k với k ⎧ψ 1 ⎫ ⎡−1 1 0 ⎤ ⎧1⎫ k ⎪ k ⎪ k ⎢ ⎥ k 2A ⎪ ⎪ ψ = ⎨ψ 2 ⎬ k = 1− 2 − 1 p = ⎨1⎬ ( ) d ⎢ ⎥ 3 ⎪ k ⎪ ⎢ ⎥ ⎪ ⎪ ⎩ψ 3 ⎭ ⎣ 0− 1 − 1⎦ ⎩1⎭ Điều cĩ thể nhận xét là cơng thức (c) vừa lập cĩ dạng tương tự cơng thức chúng ta gặp trong phương pháp phần tử hữu hạn. Mọi thủ tục tính tốn cho các ma trận trong biểu thức này trùng với thủ tục của phương pháp phần tử hữu hạn. Sau xử lý điều kiện biên cĩ thể viết biểu thức (c) dưới dạng: δψT ()KP ψ − = 0 T trong đĩ ψ = []ψ2,2 ψ 2,3 ψ 3,2 ψ 3,3 Kết quả tính theo sơ đồ đang nêu trùng với kết quả thu được từ ví dụ trên. Ví dụ 2.7: Xác định chuyển vị các nút kết cấu 2D trình bày tại hình 2. 12. Chiều dày tấm t = 0,2 m, cạnh a của tấm 2 m, bước h = 1. Mơ đun đàn hồi vật liệu E = 30 GN/m2. Từ nguyên lý cơng ảo áp dụng chung cho bài tốn 3D: δW = δεT σdV − δuT pdS ∫V ∫Sp cĩ thể viết biểu thức tính δW cho baì tốn phẳng đang xem xét: δW = δuTT D EDudA − δuT pdS = δ (Du)T EDudA − δuT pdS ∫∫ASp∫∫ASp trong đĩ
  50. 49 ⎡∂ x 0 ⎤ ⎡1υ 0 ⎤ ⎢ ⎥ Et ⎢ ⎥ ⎧u x ⎫ D = 0 ∂ y E = υ 1 0 u = ⎨ ⎬ ⎢ ⎥ 1−υ 2 ⎢ ⎥ u ⎢ ⎥ ⎢ 1−υ ⎥ ⎩ y ⎭ ⎣∂y ∂ x ⎦ ⎣0 0 2 ⎦ Hình 2.12 Thực hiện phép tich phân theo diện tích A, nhận kết quả dạng: ⎡k k ⎤ δuTT D EDudA= δvT E * 11 12 v= δ vT kv ∫A ⎢ ⎥ ⎣k21 k 22 ⎦ trong đĩ A là diện tích tam giác, E* = tEA/(1 - υ2)h2 T v = []ux1 u x 2 u x 3 u y 1 u y 2 u y 3 ⎡⎧−1⎫ ⎧ 0 ⎫ ⎤ ⎡ 1− 1 0 ⎤ ⎢⎪ ⎪ 1−υ ⎪ ⎪ ⎥ k = 1 −1 1 0 + −1 0− 1 1 = ⎢−1 γ − β ⎥ 11 ⎢⎨ ⎬[]⎨ ⎬[]⎥ ⎢ ⎥ ⎪ ⎪ 2 ⎪ ⎪ ⎣⎢⎩ 0 ⎭ ⎩ 1 ⎭ ⎦⎥ ⎣⎢ 0 − β β ⎦⎥ ⎡ ⎧−1⎫ ⎧ 0 ⎫ ⎤ ⎡ 1 υ− υ⎤ ⎢ ⎪ ⎪ 1−υ ⎪ ⎪ ⎥ k = υ 1 0− 1 1 + −1 −1 1 0 = ⎢ β− β υ ⎥ 12 ⎢ ⎨ ⎬[]⎨ ⎬[]⎥ ⎢ ⎥ ⎪ ⎪ 2 ⎪ ⎪ ⎣⎢ ⎩ 0 ⎭ ⎩ 1 ⎭ ⎦⎥ ⎣⎢− β β 0 ⎦⎥
  51. 50 ⎡ ⎧ 0 ⎫ ⎧−1⎫ ⎤ ⎡ 0 β− β ⎤ ⎢ ⎪ ⎪ 1−υ ⎪ ⎪ ⎥ k =υ −1 −1 1 0 + 1 0− 1 1 = ⎢ υ− β β ⎥ 21 ⎢ ⎨ ⎬[]⎨ ⎬[]⎥ ⎢ ⎥ ⎪ ⎪ 2 ⎪ ⎪ ⎣⎢ ⎩ 1 ⎭ ⎩ 0 ⎭ ⎦⎥ ⎣⎢−υ υ 0 ⎦⎥ ⎡⎧ 0 ⎫ ⎧−1⎫ ⎤ ⎡ β− β 0 ⎤ ⎢⎪ ⎪ 1−υ ⎪ ⎪ ⎥ k = −1 0− 1 1 + 1 −1 1 0 = ⎢− β γ −1⎥ 22 ⎢⎨ ⎬[]⎨ ⎬[]⎥ ⎢ ⎥ ⎪ ⎪ 2 ⎪ ⎪ ⎣⎢⎩ 1 ⎭ ⎩ 0 ⎭ ⎦⎥ ⎣⎢ 0− 1 − 1⎦⎥ ()1−υ ()3 −υ với β = γ = 2 2 Ma trận k tính cho tam giác phía trái, hình 2.11: ⎡1− 1 0 0 υ− υ⎤ ⎢ ⎥ ⎢ γ− β β− β υ ⎥ * ⎢ β− β β 0 ⎥ k = E ⎢ ⎥ ⎢ β− β 0 ⎥ ⎢ ĐX γ −1⎥ ⎢ ⎥ ⎣⎢ 0 ⎦⎥ Ma trận k tính cho tam giác phía phải, hình 2.11: ⎡β β 0 0 β− β ⎤ ⎢ ⎥ ⎢ γ−1 υ − β β ⎥ * ⎢ 1 −υ υ 0 ⎥ k = E ⎢ ⎥ ⎢ 1− 1 0 ⎥ ⎢ ĐX γ− β ⎥ ⎢ ⎥ ⎣⎢ β ⎦⎥ Hình 2.13
  52. 51 Vec to tải tại phần tử 5 và 7: h h 4h 5h p 5 = p p 5 = p p 7 = p p 7 = p 7 12 0 8 6 0 8 12 0 9 12 0 Điều kiện biên: uý= u x3 = u x 6 = u x 9 = 0 Hệ phương trình đại số hình thành từ thủ tục đang nêu: KV = P Trong đĩ ma trận K mang dạng sau: ⎡ 2γ ⎤ ⎢ ⎥ ⎢ −1 2γ ⎥ ⎢− β 0 2γ ⎥ ⎢ ⎥ ⎢ 0− 2β − 2 γ ⎥ ⎢ 0 0 − β 0 γ ⎥ ⎢ ⎥ ⎢ 0 0 0− 2β − 1 2γ ⎥ ⎢ υ− β 0 0 0 0 2γ ⎥ KE= * ⎢ ⎥ ⎢ 0 υ 0 0 0 0 − β γ ⎥ ⎢ β 0 − β α υ− α 0 0 2γ ⎥ ⎢ ⎥ ⎢− α α α− 2 β 0 α −2 0 − 2β 4 γ ⎥ ⎢ ⎥ ⎢ 0 − α 0 α 0 0 0− 1 0− 2β 2 γ ⎥ ⎢ 0 0 β 0 − β υ 0 0− 1 0 0 γ ⎥ ⎢ ⎥ ⎢ 0 0 − α α β− β 0 0 0− 2 0 − β α ⎥ ⎢ ⎣ 0 0 0 − α 0 α 0 0 0 0− 1 0 − β γ ⎦⎥ với α =()1 + 2 / 2 Kết quả tính: Nút ux uy 1 -2,5.10-4 0 2 -2,02.10-4 -0,44.10-4 3 0 -0,55.10-4 4 -0,059.10-4 -0,183.10-4 5 -0,0648.10-4 -0,48.10-4 6 0 0,062.10-4 7 0,15.10-4 -0,25.10-4 8 0,11.10-4 -0,55.10-4 9 0 -0,74.10-4 Đồ thị trình bày kết quả tính giới thiệu tại hình 2.14 Hình 2.14 6 DAO ĐỘNG DẦM Phương trình dao động ngang dầm
  53. 52 d4 w d2 w EI +m = 0 ) ( 4 4 . 2 dx 4 dt 2 Nếu viết w( x , t )= X ( x ). eiω t phương trình đang nêu sẽ mang dạng phương trình vi phân bậc 4: d4 X −β 4 X = 0 ) ( 5 4 . 2 dx 4 ω ω trong đĩ β =n ρA = n m (2.46) EI EI Tần số riêng bài tốn, giải bằng đường giải tích sẽ là: β 4 EI ω 2 = n ) ( 7 4 . 2 n m trong đĩ β1L = 4,73; β2L = 7,85 Hình 2.15 Dầm thẳng và phân sai Để xác định tần số riêng dầm dài L, độ cứng EI, diện tích mặt cắt ngang A, khối lượng riêng vật liệu ρ, tựa 2 đầu cần chia chiều dài dầm thành những đoạn bằng nhau. Trường hợp này hãy chia L làm 3 phần bằng nhau. Số nút ghi lần lượt từ trái sang phải: 0, 1, 2, 3. Nút ảo phía trái, cách nút 0 khoảng cách L/3 về bên trái mang số thứ tự -1, nút ảo phía phải cách nút 3 khoảng cách L/3 mang số thứ tự 4, hình 2.15 Hệ phương trình lập tại nút 1 và 2 cho phương trình vi phân bậc 4 này sẽ là: 4 XXXXXX1 −40 + 61 − 4 2 + 3 −β 1 1 = 0⎫ 4 ⎬ ( ) a XXXXXX0 −41 + 62 − 4 3 + 4 −β 1 2 = 0⎭ 4 4 ⎛ L ⎞ 4 trong đĩ β1 = ⎜ ⎟ β ( ) b ⎝ 3 ⎠ Điều kiện biên: XX= = 0 ⎫ 0 3 ⎪ dX hay là XX= = 0; XXXX= = (c) = 0 nút 0,3⎬ 0 3 −1 1 2 4 dx ⎭⎪ Hệ phương trình đại số: 4 7XXX1 − 4 2 = β 1 1 ⎫ 4 ⎬ ( ) d 4XXX1 + 7 2 = β 1 2 ⎭ ⎡ 7− 4⎤ 4 ⎡1 0⎤ Hay là ⎢ ⎥{}X − β 1 ⎢ ⎥ {}X = {}0 ( ) e ⎣− 4 7 ⎦ ⎣0 1⎦ Tần số riêng xác định từ phương trình: 15,59 EI 29,85 EI ω = ω = ( ) f 1 L2 m 2 L2 m
  54. 53 7 DAO ĐỘNG TẤM Từ phương trình tĩnh học xác định chuyển vị ngang tấm: ∂ 4 w ∂ 4 w ∂ 4 w q + 2 + = , ∂x 4 ∂x2 ∂ y 2 ∂y 4 D trong đĩ D - độ cứng tấm. Thay lực quán tính vào vị trí tải bên ngồi, phân bố theo luật xác định q(x,y) sẽ nhận được phương trình chuyển động tấm. ω 2 ∇2 ∇ 2 w −mw = 0 (2.48) D Trong khuơn khổ phương pháp lưới cơng thức () sẽ triển khai dạng: 4 2 2 Δx( ∇ ∇ w)ik = (u i+2,k - 4u i+1, k + 6uik - 4u i-1,k + ui-2,k) +2β(u i+1, k+1 + u i+1,k-1 - 2ui+1,k + 4uik -2ui,k-1 - 2ui-1,k + ui-1,k+1 + ui-1,k-1) + 2 2 ω + β (ui,k+2 - 4ui,k+1 + 6uik - 4ui,k-1 + ui,k-2 ) = mw (2.49) D i, k Δx 2 trong đĩ: β = ; Δy 2 Điều kiện biên bài tốn phù hợp với hồn cảnh cụ thể. Trường hợp tấm ngàm 2 cạnh đối diện tại y = ± a , hai cạnh cịn lại x = ± a tựa gối cứng như nêu tại ví dụ , điều kiện biên ghi lại như sau: 2 2 a ∂ 2 w tại x = ± : w = 0; = 0; 2 ∂x 2 a ∂w tại y = ± : w = 0; = 0; 2 ∂y Biểu thức (2.21) cĩ thể viết lại dưới dạng chung, quen thuộc trong các bài tốn tìm trị riêng, vecto riêng. Dưới dạng ma trận cĩ thể viết: AX= λ BX (2.50) Trong cơng thức này A = [aij] ma trận vuơng chứa các thành phần xác định theo phương pháp sai 2 2 phân hữu hạn liên quan đến tốn tử ∇ ∇ (.), B = [bij] là ma trận đường chéo chính, liên quan đến m, và λ = ω2. Xử lý phương trình (2.50) sẽ nhận được tần số dao động kết cấu và các vector trình bày dạng dao động kết cấu đĩ, ứng với tần số riêng xác định. Ví dụ 2.5: Xác định hai tần số thấp nhất dao động tấm hình vuơng trình bày tại hình 2.16. Phương trình dao động tự do, khơng cản của tấm: ω 2 ∇2 ∇ 2 w −mw = 0 D Chia tấm nhờ lưới đều, bước h = Δy = Δx = a/4.
  55. 54 Hình 2.16b Hình 2.16c Hình 2.16a Tấm vuơng và lưới Lực quán tính xác định tại các nút: h 2 F = m()100 w+ 10 w + 10 w + w 1 144 1 2 4 5 h 2 F = m()10 w+ 100 w + 10 w + w +10 w + w 2 144 1 2 3 4 5 6 h 2 F = m()10 w+ 100 w + w +10 w 3 144 2 3 5 6 h 2 F = m()10 w+ w +100 w + 10 w + 10 w + w 4 144 1 2 4 5 7 8 h 2 F = m() w+10 w + w +10 w + 100 w + 10 w + w +10 w + w 5 144 1 2 3 4 5 6 7 8 9 h 2 F = m() w+10 w + 10 w + 100 w + w +10 w 6 144 2 3 5 6 8 9 h 2 F = m()10 w+ w +100 w + 10 w 7 144 4 5 7 8 h 2 F = m() w+10 w + w +10 w + 100 w + 10 w 8 144 4 5 6 7 8 9 h 2 F = m() w+10 w + 10 w + 100 w 9 144 5 6 8 9 mω2 h 4 ω 2ma 4 Ký hiệu: λ* = = phương trình sai phân viết cho mỗi nút sẽ như sau đây. 144D 144× 256D Tại nút 1: 3816w−765 w + 120 w −765 w − 180 w + 51 w + 120 w + 51 w + w = 1 2 3 4 5 6 7 8 9 * λ ()13,122w1 +2916 w2 +162 w3 +2916 w4 +648 w5 + 36 w6 + 162 w7 + 36 w8 + 2 w 9 Tại nút 2:
  56. 55 −768w +3888 w −768 w − 132 w −1134 w −1321 w + 64 w + 144 w + 64 w = 1 2 3 4 5 6 7 8 9 * λ (13240w1 +16,524w2 +2340 w3 +720 w4 +3672 w5 +720 w6 + 40 w7 + 204 w8 + 40 w9 ) Tiếp tục tính cho đến nút thứ 9: Tại nút 9: 8w− 51 w + 120 w − 51 w − 180 w −765 w + 120 w −765 w +3816 w = 1 2 3 4 5 6 7 8 9 * λ ()2w1 + 36 w2 + 162 w3 + 36 w4 + 648 w5 +2916 w6 +162 w7 +2916 w8 + 13,122w9 Sau tập hợp nhận được phương trình ma trận: []A{} w= λ* [] B{} w ) ( 1 5 . 2 Hay là ([]A− λ* [] B){} w = 0 ) ( 2 5 . 2 Trị riêng của phương trình xác định theo các phương pháp thơng dụng, như đã bàn trong khuơn khổ phương pháp tính dao động. Kết quả tính như sau, xem hình 2.16b và 2.16c: 1 D Với 1/λ* = 28,999 (tần số gĩc ω = 35,65 ) vecto riêng cĩ dạng: 1 a 2 m []0,275 0,533 0,276 0,533 1 0,533 0,276 0,533 0,276 1 D Với 1/λ* = 7,3336, (tần số gĩc ω = 70,90 ) vecto riêng cĩ dạng: 2 a 2 m − []0,523 0 0,523− 1 0 1− 0,523 0 0,523 8 ỔN ĐỊNH TẤM Tấm mỏng dưới tác động ứng suất nén, một chiều hoặc cả hai chiều, hình 2.17, tùy thuộc tải nén, tấm cĩ nguy cơ chuyển sang trạng thái mất ổn định. Hình 2.17 Phương trình chính:
  57. 56 ⎛ ∂ 2 w ∂ 2 w ⎞ D∇4 w + λ⎜ n + n ⎟ (2.53) ⎜ x0 2 y0 2 ⎟ ⎝ ∂x ∂y ⎠ Chia lưới theo cách giản đơn nhất, h = Δx = Δy, cĩ thể triển khai cơng thức (2.53) dạng sau: D [20w− 8() w + w + w h 4 i, j i+1, j i , j+ 1 i , j− 1 2()wi+1, j + 1+ w i − 1, j + 1 + w i + 1, j − 1 + w i − 1, j − 1 + wi+ ,2 j + w i−2, j + w i , j+ 2 + w i , j− 2 ] (2.54) λ + []n() w−2 w + w + n() w−2 w + w = 0 h 2 x0 i+ 1, j i, j i− ,1 j y0 i , j+ 1i , j i , j− 1 Áp dụng cơng thức (2.54) cho các nút tính tốn, tiến hành áp đặt điều kiện biên, kết quả nhận được hệ phương trình đại số , dạng ma trận sẽ là: []A{} w−λ* [] B{ w} = 0 (2.55) Xác định λ* từ hệ phương trình (2.27) theo các phương pháp tốn đã giới thiệu. Điều lưu ý là, trị riêng nhỏ nhất (thấp nhất) sẽ đĩng vai trị hệ số tải giới hạn, cĩ nghĩa tải gây mất ổn định tấm. Ví dụ 2.6: Xác định tải giới hạn cho tấm vuơng, tựa 4 cạnh, chịu tải nén nx = ny = -λnx0. Cấu hình tấm trình bày tại hình 2.18. Mỗi cạnh tấm chia 3, đều nhau h = Δ = a/3. Tấm chia làm 9 phần tử, các nút trong tấm đánh số 1, 2, 3, 4. Phương trình vi phân chính yếu xét mất ổn định tấm: 4 2 D∇ w +λ nx0 ∇ w = 0 (2.56) Để ý đến tính đối xứng cấu hình và tải bên ngồi, cĩ thể viết: w1 = w2 = w3 = w4 = w, phương trình (2.28) cĩ thể viết dưới dạng gọn sau đây: D λn ()20− 8 − 1 − 8 + 2 − 1 4w +x0 () −4 + 1 + 1 4w = 0 (2.57) h 4 h 2 Hình 2.18 Tấm chịu ứng suất nén hai chiều
  58. 57 Đưa ký hiệu sau vào phép tính: π2 D π 2 D nx0 = = a 2 ()3h 2 Cĩ thể nhận được: 18 1,83π 2 D λ = ≈ 1,83 và n=λ n ≈ cr π 2 cr cr x0 a 2 2π 2 D Các phép tính giải tích đưa lại kết quả và n = , so với lời giải bằng phương pháp sai phân cr a 2 hữu hạn, sai số đạt 9%. Hình 2.19a Mode 1 của mất ổn định tấm Hình 2.19a Mode 2 của mất ổn định tấm
  59. 58 Chương 3 PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN 1 PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN Phương pháp phần tử hữu hạn trong phần này tài liệu thuộc phương pháp chuyển vị của cơ học kết cấu. 1 Tổng thế năng hệ thống đàn hồi tìm ở dạng: Π = U – W, trong đĩ: U = ∫{}{}εT σ dV - cơng 2 V biến dạng, W= ∫{}{} PT u dS - cơng do ngoại lực thực hiện, {P} – vector ngoại lực, {u} – vec to S chuyển vị. Biến dạng và ứng suất trong hệ thống kết cấu xác định từ lý thuyết đàn hồi: {ε}= [C]{σ}, {σ} = [D]{ε}. Ma trận [C] và [D] trình bày tính chất vật liệu. Tổng năng lượng hệ thống: 1 Π = ∫{}{}{}{}εT σ dV− ∫ PT u dS (3.1) 2 V S Điều kiện biên Hình 3.1 Cơng biến dạng, cơng ngoại lực và điều kiện biên Biến phân hàm năng lượng δΠ: ⎧ ⎫ ⎪1 T T ⎪ δΠ = δ⎨ {}{}{}{} ε σ dV− P u dS⎬ = δ{ ε }T [D ]{ ε } dV− { P }T δ { u } dS (3.2) ⎪2 ∫ ∫ ⎪ ∫ ∫ ⎩ V Su ⎭ V S Điều kiện δΠ = 0 cho phép viết: ∫δ{ ε }T [D ]{ ε } dV= ∫ { P }T δ { u } dS (3.3) V S * * Điều kiện động lực học p = p phải được thỏa mãn trên biên Sp thuộc S và điều kiện động học u = u thỏa mãn trên biên Su thuộc S. Điều này cho phép viết (3.3) dạng sau: ∫δ{ ε }T [D ]{ ε } dV− ∫ δ { u }T pdS = 0 (3.3b) V SP hay δUWUW−δ =δ () − = 0 (3.4)
  60. 59 Trong phương pháp tính phần tử hữu hạn (PTHH) tiến hành xử lý phương trình (3.4), đưa bài tốn về hệ phương trình đại số tuyến tính, chuyển vị tại nút {u} đĩng vai trị ẩn số. Trường hợp chung của vật thể 3D từ vật liệu đàn hồi, chuyển vị, biến dạng thể hiện qua quan hệ: Hàm chuyển vị: {}u= [] u w T = [N ]{δ} (v ) 5 . 3 trong đĩ {δ} – vecto chuyển vị nút, [N] – hàm hình dáng. T Vector ngoại lực: {}P= [] px p y p z )( 6 . 3 Vector biến dạng: ⎧ ∂u ⎫ ⎡ ∂ ⎤ ⎪ ⎪ ⎢ 0 0 ⎥ ⎡ ∂ ⎤ ∂x ∂x N 0 0 ⎪ ∂v ⎪ ⎢ ∂ ⎥ ⎢ ∂x i ⎥ ⎧ ε ⎫ ⎪ ⎪ ⎢ 0 0 ⎥ ⎢ ∂ ⎥ x ⎢ ⎥ ⎪ ⎪ ⎪ ∂y ⎪ ⎢ ∂y ⎥ 0 N i 0 ε ⎢ ∂y ⎥ ⎪ y ⎪ ⎪ ∂w ⎪ ⎢ ∂ ⎥⎧ u ⎫ ⎪ ⎪ ⎢ 0 0 ⎥ ⎢ ∂ ⎥ ⎪ ε z ⎪ ⎪ ⎪ ⎪ ⎪ 0 0 N ε = = ∂z = ⎢ ∂z ⎥ v = [B ]{δ } ⎢ i ⎥ {} ⎨ ⎬ ⎨ ∂u ∂v ⎬ ∂ ∂ ⎨ ⎬ []B = ⎢ ∂z ⎥ γ xy ⎢ ⎥ i ∂ ∂ ⎪ ⎪ ⎪ + ⎪ 0 ⎪w ⎪ ⎢ N N 0 ⎥ ⎪ ⎪ ⎪ ∂y∂ x ⎪ ⎢∂y∂ x ⎥⎩ ⎭ i i γ yz ⎢ ∂y ∂x ⎥ ⎪ ⎪ ⎪∂v ∂w ⎪ ⎢ ∂ ∂ ⎥ ⎢ ∂ ∂ ⎥ ⎪γ ⎪ + ⎢ 0 ⎥ 0 N N ⎩ zx ⎭ ⎪ ∂z ∂y ⎪ ∂z∂ y ⎢ i i ⎥ ⎪ ⎪ ⎢ ⎥ ⎢ ∂z ∂y ⎥ ∂w ∂u ⎢ ∂ ∂ ⎥ ∂ ∂ ⎪ + ⎪ 0 ⎢ N 0 N ⎥ ⎩⎪ ∂x ∂z ⎭⎪ ⎣⎢∂z∂ x ⎦⎥ ⎣⎢ ∂z i ∂x i ⎦⎥ Cơng thức (3.3) trở thành: 1 Π = {δ }TT [B ] [ D ][ B ]{δ } dV − { δ }TT [N ] { p } dS = 0 (3.7) 2 ∫V ∫S Trong phương pháp phần tử hữu hạn, miền V được mơ hình hĩa thành E phần tử Ve, e = 1, 2, , E với giả thiết trong mỗi phần tử thỏa mãn các điều kiện cân bằng như vừa nêu. Thế năng mỗi phần tử ghi lại dạng πe, cịn thế năng tồn hệ đúng bằng tổng thế năng tất cả phần tử họp thành: E Π = ∑π e ) ( 8 . 3 e=1 Thế năng phần tử: 1 π= { ε }T [D ]{ ε } dV− uT p dS . e ∫ e ∫{}{ } 2 Ve S Tổng thế năng Π cho tồn hệ thống: E E 1 T ⎡ T ⎤ T T Π ={ Δ }⎢∑∫∫∫ [B ] [ D ][ B ] dV ⎥ {Δ } − { Δ }∑∫∫ [N ] {} p dS ) ( 9 . 3 2 e=1 e=1 ⎣ Ve ⎦ S p Trong các cơng thức trên, mỗi biểu thức cĩ tên gọi sau đây. T [k ]e = [ B ] [ D ][ B ] dV , ma trận cứng phần tử, ∫V e T {}[]{}Pe = N pS dS , vector lực mặt, ∫S Phương trình cân bằng tồn hệ: ) 0 1[K]{ .Δ 3 ( } P { = } E E với: [][]K= ∑ k e , {}PP= ∑{}e (3.11) e=1 e=1 {Δ} – vector chuyển vị các bậc tự do hệ thống, xác định trong hệ tọa độ chung.
  61. 60 2. THỨ TỨ GIẢI BÀI TỐN CƠ HỌC KẾT CẤU THEO PHƯƠNG PHÁP PTHH Bước 1: Phân chia vật thể đang xem xét thành số lượng hữu hạn các phần tử. Quá trình này cịn được gọi là “lý tưởng hĩa” hay “rời rạc hĩa”. Thực tế đây là quá trình mơ hình hĩa kết cấu, chuyển từ kết cấu thực tế thành tập hợp của nhiều cơ cấu vừa tách từ chủ thể, xem hình .2. Bước 2: Xây dựng mơ hình chuyển vị trong mỗi phần tử: {}ue = [] N {}δ e trong đĩ [N] –hàm hình dáng, {δ}e- vector các bậc tự do chuyển vị nút của phần tử. Bước 3: Xác lập ma trận đặc trưng gọi là ma trận cứng và vector lực cho mỗi phần tử trên cơ sở nguyên lý thế năng tối thiểu. Trong những bài tốn thuộc cơ học vật rắn phiếm hàm thế năng của hệ thống được hiểu như tổng thế năng các phần tử cấu thành. E Π =∑()π e −W e e=1 Từ điều kiện δΠ = 0 xây dựng hệ phương trình đại số miêu tả quan hệ lực – chuyển vị. Bước 4: Xử lý hệ phương trình và giải hệ phương trình đại số tuyến tính. Kết quả giải phương trình sẽ là chuyển vị nút trong hệ tọa độ chung. Cần thiết chuyển đổi chuyển vị từ hệ tọa độ chung sang hệ tọa độ cục bộ, gắn liền phần tử. Bước 5: Thực hiện các phép tính lực căn cứ quan hệ ứng suất – biến dạng. Hình 3.3 Các bước tính trong phương pháp PTHH Hình 3.2 Mơ hình hĩa Minh họa thủ tục đang trình bày qua ví dụ xác định chuyển vị dầm liên tục, ngàm hai đầu. Dầm dài L =2l, độ cứng dầm EI, chịu tác động lực tập trung P tại giữa nhịp. 1) Mơ hình hố dầm như tại hình 3.4 Chia dầm làm hai phần tử, đánh số PT 1 cho phần tử nằm bên trái, PT 2 cho phần tử cịn lại. Với hai phần tử trên dầm cĩ 3 nút, đánh theo thứ tự: nút 1 tại ngàm trái, nút 2 giữa dầm, và nút 3 tại ngàm phải. Theo cách làm này, phần tử đầu cĩ hai nút 1 và 2; nút 2 và 3 thuộc phần tử bên phải. Hình 3.4 Dầmngàm
  62. 61 2) Bỏ qua trường hợp xoắn, kéo, nén dọc trục, dưới tác động lực P dầm bị chuyển vị theo hướng thẳng đứng, hướng trục Oy và xoay theo trục Oz. Phương trình chuyển vị trong mỗi phần tử được trình bày dưới dạng quen thuộc: {w}e = [N] {δ} (a) Hàm hình dáng [N] cĩ thể nhận bằng hệ hàm Hermite, trình bày tại phân bàn phần tử BEAM 2D. T Vector chuyển vị hai nút của phần tử thứ nhất {δ} = [ u1 u2 u3 u4 ] , của phần tử thứ hai T {δ} = [ u3 u4 u5 u6 ] , trong đĩ từ điều kiện ngàm hai đầu cĩ thể thấy u1 = u1 = u5 = u6 = 0. 3) Sử dụng phương pháp năng lượng xây dựng ma trận cứng và vector lực tồn hệ thống: Các thành phần ma trận [k]: ⎡ 12 6l − 12 6l ⎤ ⎢ 4l2 − 6 l 2 l 2 ⎥ k k EI ⎢ ⎥ ⎡ 11 12 ⎤ []k e = 3 = ⎢ ⎥ (b) l ⎢ 12− 6l⎥ ⎣k21 k 22 ⎦ ⎢ 2 ⎥ ⎣ĐX 4l ⎦ Tại phần tử 1, chuyển vị nút 1 và 2 trịêt tiêu, cịn lại phần w1-2 = N3 u3 + N4 u4. Ma trận [k]1 tính cho phần tử 1 chỉ cần giữ lại phần k22 gồm bốn thành phần, ứng với u3 và u4: 3 4 EI ⎡12− 6⎤ 3 []k )1( = 3 ⎢ ⎥ (c) l ⎣− 6 4 ⎦ 4 Với phần tử số hai ma trận cứng, giữ lại phần ứng với u3 và u4 : 3 4 EI ⎡12 6⎤ 3 []k )2( = 3 ⎢ ⎥ (d) l ⎣ 6 4⎦ 4 Tập hợp ma trận cứng từ hai ma trận phần tử: 2 EI ⎡12+ 12 − 6 + 6⎤ K=[] k = (e) [] ∑ e 3 ⎢ ⎥ e=1 l ⎣ ĐX 4+ 4 ⎦ 1 Từ U = {u}T [K]{u }, cĩ thể viết biểu thức tính δU: δU = δ{u}T [K]{u} (f) 2 Trường hợp này vector chuyển vị cĩ dạng: ⎧u3 ⎫ {}u = ⎨ ⎬ (g) ⎩u4 ⎭ Biến phân δW của cơng ngoại lực tính từ biểu thức: δW = -δ{u}T P/2. Thỏa mãn điều kiện δU - δW = 0 cĩ thể viết: 8EI ⎡24 0⎤⎧ u ⎫ 1 ⎧P⎫ δ u T 3 − = 0 (h) {} 3 ⎢ ⎥⎨ 1 ⎬ ⎨ ⎬ 2L ⎣ 0 8⎦⎩ 2 u4 ⎭ 2 ⎩0⎭ 4) Xây dựng hệ phương trình đại số [K]{u} = {P}: Vì rằng δ{u}T ≠ 0, biểu thức trong ngoặc phải bằng 0. 8EI ⎡24 0⎤⎧ u ⎫ ⎧P⎫ 3 = (i) 3 ⎢ ⎥⎨ 1 ⎬ ⎨ ⎬ L ⎣ 0 8⎦⎩ 2 u4 ⎭ ⎩0⎭
  63. 62 Từ đây nhận được giá trị cho u3 và u4. PL3 u = ; 3 192EI (j) u4 = 0 3. MA TRẬN CỨNG PHẦN TỬ. MA TRẬN CỨNG HỆ THỐNG Ma trận cứng phần tử Ma trận cứng phần tử tính trong hệ tọa độ cục bộ, gắn liền phần tử, xác định bằng biểu thức: T [k ]e = [ B ] [ D ][ B ] dV ∫V e Tập hợp ma trận cứng của hệ thống Cơng thức chung trong phương pháp phần tử hữu hạn cĩ dạng: [K}{Δ} = {P} E E E trong đĩ []K= ∑ ke ; {}Δ = ∑δ e ; {}P= ∑ pe ; e=1 e=1 e=1 Ví dụ 1: Lập ma trận cứng, vector tải cho giàn hình chữ L, cạnh đều, ngàm tại hai đầu. Tải trọng phân bố đều q. Độ cứng chịu uốn dầm EJ, độ cứng chịu xoắn GJT, diện tích mặt cắt ngang dầm A. Phần tử 1 ngàm tại 3, gặp phần tử thứ hai, vuơng gĩc với phần tử đầu tại 1. Phần tử thứ hai ngàm tại 2, xem hình. Hệ toạ độ chung bố trí để trục Oy trùng với chiều dọc của phần tử thứ hai, trục Oz hướng lên trên. Hệ tọa độ cục bộ gắn liền với phần tử, trong đĩ trục O’x’ trùng tâm dầm. T Chuyển vị mỗi nút thể hiện bằng vector {δ} = [ wj θXj θYj] , j = 1, 2, 3. 3 1 q z y x 2 Hình 3.6 Giàn Ma trận cứng phần tử 3-1, và vector lực tại phần tử tính trong hệ tọa độ cục bộ trùng với hệ tọa độ chung: ⎡12EJ 6EJ ⎤ 0 ⎢ 3 2 ⎥ L L 4 ⎧1⎫4 )1( )1( ⎢ ⎥ ⎡k11 k12 ⎤ )1( GJ T qL ⎪ ⎪ []k 1 = ⎢ ⎥ với k22 = ⎢ 0 0 ⎥ 5 {}p2 = ⎨0⎬5 )1( )1( L 2 ⎣k21 k22 ⎦ ⎢ ⎥ ⎪ L ⎪ 6EJ 4EJ 6 ⎩ 6 ⎭6 ⎢ 0 ⎥ ⎣⎢ L2 L ⎦⎥ Với phần tử 1-2, các phép tính được thực hiện đầu tiên trong hệ tọa độ cục bộ: ⎡12EJ 6EJ ⎤ 0 ⎢ 3 2 ⎥ L L ⎧1⎫ )2( )2( ⎢ ⎥ ⎡k11 k12 ⎤ )2( GJ T )2( qL ⎪ ⎪ []k 2 = ⎢ ⎥ với k11 = ⎢ 0 0 ⎥ {}p2 = ⎨0⎬) )2( )2( L 2 ⎣k21 k22 ⎦ ⎢ ⎥ ⎪ L ⎪ 6EJ 4EJ ⎩ 6 ⎭ ⎢ 0 ⎥ ⎣⎢ L2 L ⎦⎥
  64. 63 Ma trận chuyển, tính cho trường hợp dầm trực giao, gĩc chuyển 90°. ⎡1 0 0⎤ ⎢ ⎥ []R = ⎢0 0 1⎥ ⎣⎢0− 1 0⎦⎥ Tính chuyển sang hệ tọa độ chung, ma trận cứng và vector tải phần tử thứ hai tính theo cơng thức T T [k]chung = [R]{k]cục_bộ[R] cịn {p}chung = [R] {p}cục_bộ : ⎡ 12EJ 6EJ ⎤ − 0 ⎢ 3 2 ⎥ L L 4 ⎧ 1 ⎫4 ⎡k)2( k )2( ⎤ ⎢ 6EJ 4EJ ⎥ qL ⎪ ⎪ 11 12 )2( ⎢ ⎥ )2( L []k 2 = ⎢ )2( )2( ⎥ với k11 = − 2 0 5 {}p2 =⎨ − 6 ⎬5 ⎣k21 k22 ⎦ ⎢ L L ⎥ 2 ⎪ ⎪ GJ 6 ⎩ 0 ⎭6 ⎢ 0 0 T ⎥ ⎣⎢ L ⎦⎥ Sau tập hợp ma trận cứng và vector tải sẽ nhận được: 4 5 6 )1( )1( ⎧ ⎫ ⎡k k 0 ⎤ 4 1 ⎧ 1 ⎫ 11 12 ⎪ ⎪ ⎢ )1( )1( )2( )2( ⎥ pL ⎪ ⎪ pL ⎪ 1⎪ [K ]ch = [] k1 + [] k 2 = k21 k22 + k11 k12 5 và {}p ch = ⎨0⎬ +⎨ − ⎬ ⎢ ⎥ 2 2 6 ⎢ )2( )2( ⎥ ⎪1⎪ ⎪ ⎪ 0 k21 k22 6 0 ⎣ ⎦ ⎩⎪6⎭⎪ ⎩ ⎭ Từ đĩ cĩ thể thấy: 4 5 6 ⎡ 24EJ 6EJ 6EJ ⎤ − ⎧ ⎫ ⎢ 3 2 2 ⎥ ⎪ ⎪ L L L ⎧ w1 ⎫ 2 ⎢ 6EJ GJ 4EJ ⎥⎪ ⎪ pL⎪ L⎪ ⎢ T ⎥ − 2 + 0 ⎨θ X1 ⎬ =⎨ − ⎬ ⎢ L L L ⎥⎪ ⎪ 2 ⎪ 6 ⎪ ⎢ 6EJ GJ 4EJ ⎥⎩θ Y1 ⎭ L 0 T + ⎪ ⎪ ⎣⎢ L2 L L ⎦⎥ ⎩⎪ 6 ⎭⎪ Để ý đến quan hệ θX1 = θY1 cĩ thể viết lại hệ phương trình trên dưới dạng: ⎡ 24EJ 12EJ ⎤ 3 − 2 ⎧ 2 ⎫ ⎢ L L ⎥⎧ w1 ⎫ pL ⎪ ⎪ ⎢ ⎥⎨ ⎬ = ⎨ L⎬ 12EJ 2GJ T 8EJ θ 2 − ⎢− + ⎥⎩ X 1 ⎭ ⎩⎪ 3 ⎭⎪ ⎣ L2 L L ⎦ Cách làm vừa nêu đưa lại kết quả, ma trận cứng tồn hệ [K] được sắp xếp thành ma trận vuơng, đối xứng qua đường chéo chính. Ưu điểm hàng đầu của phương pháp chuyển vị này chính ở chỗ đĩ. Các thành phần mang giá trị 0 trong ma trận này thường rất nhiều, tùy thuộc vào cách tổ chức đánh dấu số thứ tự bậc tự do bài tốn. Thơng lệ đây là ma trận vuơng bậc bằng số bậc tự do. Tùy thuộc mật độ chiếm chỗ trong ma trận, điều này phụ thuộc vào đánh số thứ tự nút, cĩ thể gặp ma trận cĩ rất ít ơ mang giá trị khác 0. Trường hợp này chúng ta gặp ma trận lỗng. Trường hợp đánh số nút theo sơ đồ thỏa đáng người ta thu được ma trận [K] dạng ma trận giải hẹp. Ma trận dạng này cĩ cấu trúc đặc biệt, đường chéo chính điền đầy, giá trị các thành phần trên đường chéo chính luơn lớn. Một số đường chéo nằm gần đường chéo chính, trên và dưới đường này điền gần đầy dữ liệu thu được từ phép cọng ma trận phần tử. Chiều rộng giải ma trận cứng tính bằng số đường chéo chứa dữ liệu khác 0. Từ ma trận giải hẹp, sử dụng tính đối xứng của ma trận cứng tồn hệ người ta chỉ sử dụng một nửa giải cộng
  65. 64 với đường chéo chính khi tính tốn. Cách sắp xếp ma trận cứng dưới dạng ma trận giải hẹp thuận lợi khi giải quyết những bài tốn cỡ bé. Thịnh hành ở USA và các nước châu Âu cách sắp xếp theo đường chân trời (sky line). Các thành phần cĩ nghĩa của nửa ma trận đối xứng cùng các thành phần trên đường chéo chính ma trận [K] được chuyển địa chỉ và sắp vào vector {s}, độ lớn bằng tổng các thành phần đang nêu. Hình 3.7 Sắp xếp ma trận cứng trong skyline 4. ÁP ĐẶT TẢI Tải bên ngồi tác động lên kết cấu tính theo các biểu thức nêu tại phần bàn về phần tử tiêu biểu. Những vật thể về mặt hình học đối xứng qua trục mang tính chất đặc trưng, kết cấu hai phần, hoặc bốn phần, tính qua trục hoặc qua hai trục, giống nhau song hướng ngược nhau. Sử dụng tính đối xứng hình học và đối xứng hoặc phản đối xứng (cịn gọi á đối xứng) tải trọng chúng ta cĩ thể tháo dỡ phần đối xứng, giữ lại một phần kết cấu khi tính. Hình 3.8 Cấu hình đối xứng, tải đối xứng Hình 3.9 Cấu hình đối xứng, tải á đối xứng (phản đối xứng) Trường hợp tiếp theo người tính thường gặp, tải đặt tại vị trí khơng phải tại tâm đối xứng, tác động lên kết cấu cấu hình đối xứng. Cách xử lý thơng dụng là sử dụng nguyên lý cọng tác dụng khi tính, cịn mơ hình kết cấu dung khi tính là phần đối xứng. Hình 3.10 nêu cách làm này.
  66. 65 Hình 3.10a Sử dụng tính á đối xứng kết cấu Hình 3.10b Sử dụng cấu hình đối xứng Trường hợp kết cấu 2D đối xứng qua hai trục, chúng ta cĩ quyền giữ lại chỉ ¼ kết cấu cho việc tính tốn. Mơ hình này chỉ phù hợp khi tải trọng tác động trực tiếp mang tính đối xứng hoặc chí ít phản đối xứng. Hình 3.11 Hình 3.12a
  67. 66 Hình 3.12b Trường hợp chỉ cĩ cấu hình đối xứng, tải á đối xứng, nên tiến hành xử lý điều kiện biên tại trục đối xứng để đưa bài tốn về dạng “phản đối xứng”. 5. XỬ LÝ ĐIỀU KIỆN BIÊN Xử lý điều kiện biên là điều bắt buộc trước khi giải hệ phương trình tìm chuyển vị. Những cách làm thực tế được xem xét dưới đây. 1) Thay thế điều kiện biên vào vị trí của nĩ trong vector chuyển vị Ma trận [K] xác lập cho hệ thống, về nguyên tắc sẽ là bậc đúng bằng bậc của độ tự do các nút. Thực tế trong tổng số bậc tự do đĩ, khơng phải tất cả đều là ẩn số. Giả sử rằng vector chuyển vị {U}, tức vector chứa ẩn theo nghĩa ban đầu gồm hai nhĩm, nhĩm thỏa mãn điều kiện biên, cĩ nghĩa khơng cịn là ẩn và nhĩm thứ hai chỉ chứa ẩn của bài tốn: ⎧U1 ⎫ {}U = ⎨ ⎬ , trong đĩ U2 - giá trị biết trước, U1 - ẩn cần tìm. Bài tốn {K}{U} = {P} giờ cĩ ⎩U 2 ⎭ thể viết lại: ⎡[][]KK11 12 ⎤⎧U1 ⎫ ⎧P1 ⎫ ⎢ ⎥⎨ ⎬ = ⎨ ⎬ (3.12) ⎣[][]KK21 22 ⎦⎩U 2 ⎭ ⎩P2 ⎭ Tiến hành giải phương trình (70) theo cách sau: []KUKUP11{}{} 1 + []12 2 = {}1 (3.13) []KUPKU11{} 1={} 1 − []12{} 2 (3.14) T và []KUKUP12{}{} 1 + []22 2 = { 2 } (3.16) Trường hợp {K11] khơng suy biến, nghiệm U1 sẽ được tìm theo cách thức: −1 {}UKPKU1 = []11 (){ 1}− []12{} 2 (3.17) Trường hợp {U1} được biết trước, hay cịn nĩi theo cách thỏa mãn điều kiện biên, nghiệm {P2} xác định từ biểu thức (3.16 ). Trường hợp {U2} chỉ chứa các giá trị 0, cần tiến hành loại bỏ tất cả dịng và cột của ma trận [K] tương ứng với dịng của {U2} . Cách tính tiếp tục chỉ cịn là: [K11] {U1} = {P1} (3.18) 2) Phương pháp tăng thành phần ma trận cứng nằm trên đường chéo chính, tương ứng với chuyển vị biết trước, làm thay đổi đặc trưng ma trận, thỏa mãn điều kiện biên. Giả sử {K} đã đượcc lập theo cách vừa trình bày trên, vector chuyển vị nút cùng bậc với hệ thống phương trình, trong đĩ cĩ thành * phần thứ j đã xác định Uj = U . Thay Kjj bằng Kjj + α, với α rất lớn. Thay giá trị tương ứng tại vector tải bằng giá trị αU*.
  68. 67 n ⎛ ⎞ * Cĩ thể viết rằng:αUKUUj + ⎜∑ jk k ⎟ = α ⎝ k =1 ⎠ n * * Nghiệm Uj ≈ U nếu αUKU>> ∑ jk k k =1 Ví dụ 5: Xử lý điều kiện biên bài tốn tiếp theo đây. Hệ phương trình đại số tuyến tính thành lập từ phép tính cĩ dạng: ⎡K11 0 0 K14 ⎤⎧U1 ⎫ ⎧P1 ⎫ ⎢ ⎥⎪ ⎪ ⎪ ⎪ 0 KKK22 23 24 ⎪U 2 ⎪ ⎪P2 ⎪ ⎢ ⎥⎨ ⎬ = ⎨ ⎬ ⎢ ⎥ 0 KK23 33 0 ⎪U 3 ⎪ ⎪P3 ⎪ ⎢ ⎥⎪ ⎪ ⎪ ⎪ ⎣KK14 24 0 K 44 ⎦⎩U 4 ⎭ ⎩P4 ⎭ * Điều kiện biên đặt ra U1 = U . Các bước tiếp theo trong việc xử lý điều kiện biên: - Thay đổi thành phần K11 của ma trận cứng 15 15 * ⎡K11 +10 0 0 K14 ⎤⎧U1 ⎫ ⎧10 ×U ⎫ ⎢ ⎥⎪ ⎪ ⎪ ⎪ 0 KKK ⎪U 2 ⎪ ⎪ P ⎪ ⎢ 22 23 24 ⎥⎨ ⎬ = ⎨ 2 ⎬ ⎢ ⎥ U 0 KK23 33 0 ⎪ 3 ⎪ ⎪ P3 ⎪ ⎢ ⎥⎪ ⎪ ⎪ ⎪ ⎣⎢ KK14 24 0 K 44 ⎦⎥⎩U 4 ⎭ ⎩ P4 ⎭ - Chuyển thành phần thứ 1-1 tại đường chéo chính về đơn vị, cĩ nghĩa chia tất cả thành phần trong dịng và cột chứa đại lượng vơ cùng lớn cho chính đại lượng đĩ. * ⎡1 0 0 0 ⎤⎧U1 ⎫ ⎧ U ⎫ ⎢ ⎥⎪ ⎪ ⎪ ⎪ 0 KKK22 23 24 ⎪U 2 ⎪ ⎪ P ⎪ ⎢ ⎥⎨ ⎬ = ⎨ 2 ⎬ ⎢0 KK 0 ⎥ U 23 33 ⎪ 3 ⎪ ⎪ P3 ⎪ ⎢ ⎥⎪ ⎪ ⎪ * ⎪ ⎣0K 24 0 K 44 ⎦⎩U 4 ⎭ ⎩PKU4− 14 ⎭ - Hệ phương trình giờ mang dạng: ⎡KKK22 23 24 ⎤⎧U 2 ⎫ ⎧ P2 ⎫ ⎪ ⎪ ⎪ ⎪ ⎢KK 0 ⎥ U = P ; UU= * ⎢ 23 33 ⎥⎨ 3 ⎬ ⎨ 3 ⎬ 1 ⎪ ⎪ ⎪ * ⎪ ⎣⎢K 24 0 K 44 ⎦⎥⎩U 4 ⎭ ⎩PKU4+ 14 ⎭ 6 GIẢI HỆ PHƯƠNG TRÌNH ĐẠI SỐ TUYẾN TÍNH Phần việc quan trọng hàng đầu trong giải những bài tốn cơ học theo phương pháp PTHH là giải hệ phương trình đại số tuyến tính vừa lập dưới dạng [K]{D} = {P}. Vì rằng ma trận [K] cĩ những đặc tính đáng để ý nên các cách xử lý ma trận được chú ý riêng. Ma trận [K] thành lập trong khuơn khổ phương pháp PTHH cĩ tính đối xứng, và nếu cĩ cách tổ chức hợp lý lúc đánh số thứ tự cho các bậc tự do các nút ma trận này cĩ giải hẹp. Xử lý nh ững ma trận cỡ rất lớn song các thành phần của ma trận lấp khơng đầy cĩ thể sử dụng các thuật tốn lập cho ma trận lỗng. Các phần mềm đang dùng ngày nay đưa cách xử lý này trong nội dung.
  69. 68 Phương pháp loại trừ của Gauss Dùng phương pháp loại trừ phải qua hai giai đoạn: a) chuyển hĩa ma trận [K] về ma trận tam giác, gọi là quá trình “tam giác hĩa”, b) giải hệ phương trình với ma trận tam giác trên vừa lập hay là quá trình lùi. Tam giác hĩa ma trận K Ma trận K trong thành phần phương trình cân bằng cĩ dạng như sau: ⎡ KK11 12 L K 1 n ⎤ ⎧ D 1 ⎫ ⎧ P1 ⎫ ⎢ ⎥ ⎪ ⎪ ⎪ ⎪ KK21 22 L K 2 n ⎪ D 2 ⎪ ⎪ P2 ⎪ ⎢ ⎥ ⎨ ⎬ = ⎨ ⎬ ⎢ MMLM ⎥ ⎪ MM⎪ ⎪ ⎪ ⎢ ⎥ ⎪ ⎪ ⎪ ⎪ ⎣ KKn1 n 2 L K nn ⎦ ⎩ D n ⎭ ⎩ Pn ⎭ Tam giác hĩa tiến hành loại trừ Dj, j = 1, 2, , n –1, tính từ phương trình j + 1 đến n. Trong quá trình này các cột, lần lượt từ j = 1 đến n – 1 bị thay bằng 0. Lần thứ nhất loại trừ tiến hành theo cách cụ thể: 1 D1 =()PKDKD1 − 12 2 L− 1n n K11 Thay D1 vào các phương trình sau 1 của (a) chúng ta cĩ dạng sau đây của ma trận cứng sau lần loại trừ đầu tiên. ⎡ K 11 K 12 L K 1n ⎤ ⎧ P1 ⎫ ⎢ ⎥ ⎧ D1 ⎫ ⎪ ⎪ K 21 K 21 K 21 ⎢ 0 K 22 − K 12 L K 2 n − K 1n ⎥ ⎪ ⎪ ⎪ P2 − P1 ⎪ K K ⎪ D 2 ⎪ ⎪ K ⎪ ⎢ 11 11 ⎥ ⎨ ⎬ = ⎨ 11 ⎬ ⎢ M M L M ⎥ ⎪ M ⎪ ⎪ M ⎪ ⎢ K K ⎥ K 0 K − n1 K K − n1 K ⎪ D ⎪ ⎪ P − 21 P ⎪ ⎢ n 2 12 L nn 1n ⎥ ⎩ n ⎭ n 1 ⎣ K 11 K 11 ⎦ ⎩⎪ K 11 ⎭⎪ Cơng thức chung xác định các thành phần ma trận K và vector P cho lần tính kế tiếp được viết theo cách sau: (0) ⎫ (1) (0) (0) K 1 j KKKij =ij − i1 (0) ⎪ K 11 ⎪ (0) ⎬i, j = 2,3,L , n P PPK(1)= (0) − (0) 1 ⎪ i i i1 (0) ⎪ K 11 ⎭ Sau khi loại trừ D1 D2 . . . Dn-1 sẽ nhận được cơng thứcsau: ⎡ K 11 LLLLL K 1n ⎤ ⎧ D 1 ⎫ ⎧ P1 ⎫ ⎢ 0 K 1 K 1 ⎥ ⎪ D ⎪ ⎪ P 1 ⎪ ⎢ 22 2 n ⎥ ⎪ 2 ⎪ ⎪ 2 ⎪ 2 2 1 ⎢ 0 0 K 33 LLL K 3 n ⎥ ⎪ D 3 ⎪ ⎪ P3 ⎪ ⎢ ⎥ ⎪ ⎪ ⎪ ⎪ ⎢ 0 M 0 O M ⎥ ⎪ M ⎪ ⎪ M ⎪ ⎨ ⎬ = ⎨ ⎬ ⎢ 0 K s −1 K s −1 ⎥ D P s −1 ⎢ M ss L sn ⎥ ⎪ s ⎪ ⎪ s ⎪ ⎢ 0 MM 0 OM ⎥ ⎪ M ⎪ ⎪ M ⎪ ⎪ ⎪ ⎪ ⎪ ⎢ 0 ⎥ ⎢ MM M ⎥ ⎪ ⎪ ⎪ ⎪ n −1 ⎪ ⎪ ⎪ n −1 ⎪ ⎣⎢ 0 0 0 0 0 L K nn ⎦⎥ ⎩ D n ⎭ ⎩ Pn ⎭ Quá trình ngược bắt đầu từ phương trình cuối: n−1 n−1 n−1 Pn KDPnn . n= n từ đây Dn = n−1 K nn Các nghiệm tìm lần lượt: n PKDi − ∑ ij j j= i +1 Di = K ii
  70. 69 Phân tích LU và LDU Ma trận [K] trong phương pháp này phân thành tích hai ma trận [L] – dưới, và [U] – trên, viết là [K] = [L][U]. Cơng thức được minh họa tiếp theo: ⎡KK11 12 L K1n ⎤ ⎡L11 0 0 0 0 ⎤⎡1 UU12 13 L U1n ⎤ ⎢KK K ⎥ ⎢LL 0 0 ⎥⎢0 1 UU ⎥ []A = ⎢ 21 22 L 2n ⎥ = ⎢ 21 22 L ⎥⎢ 23 L 2n ⎥ ⎢ M ⎥ ⎢ M 0 ⎥⎢M 1 M ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎣K n1 LL K nn ⎦ ⎣LLLLn1 n 2 n 3 L nn ⎦⎣0 0 0L 1 ⎦ Các hệ số ma trận tam giác tính như sau: j −1 Lij = Kij − ∑ Lik U kj i≥ j k =1 j −1 KLUij − ∑ ik kj k =1 U ij = ; i< j Lii U ii = 1 Giải hệ phương trình theo cách sau. Thay biểu thức [K] = [L][U] vào phương trình cân bằng cĩ thể thấy: [L][U]{D} = {P}. Nếu đặt {U}{D} = {Z}, phương trình cân bằng trở thành: [L]{Z} ={P} L11Z1 = P1 L11Z2 + L21Z2 = P2 . . . L11Z2 + L21Z2 + . . . Ln1Zn = Pn Nghiệm đầu tiên Z1 giải được sẽ thay vào vị trí vế phải phương trình tiếp theo: D1 + U12 D2 + . . . +U1n Dn = Z1 U2 + U12 D2 + + U2n Dn = Z2 U3 + + U3n Dn = Z3 . . . Dn = Zn Ma trận [U] xác lập theo cách phân chia lại ma trận K (decomposition) như cách làm của Choleski. Ma trận [K] phân thành [U]T.U, trong đĩ [U] chiếm trọn tồn bộ nửa phần trên đường chéo chính của [K]: ⎡u11 u 11 L u11 ⎤ ⎢ 0 u u ⎥ []U = ⎢ 11 L 11 ⎥ ⎢ MMM ⎥ ⎢ ⎥ ⎣ 0 0 0 u11 ⎦ trong đĩ
  71. 70 ⎫ u11 = K11 a ⎪ 1 j ⎪ u1 j =j = 2,3, , n u11 ⎪ ⎪ i−1 ⎪ 2 uii = Kii −∑ uki ; i = 2,3, , n ⎬ k =1 ⎪ 1 i−1 ⎪ uij =Kii − ∑ uki u kj i= 2,3, , n ; j= i + 1, i + 2, , n⎪ uii k =1 ⎪ uij =0; i > j ⎭⎪ Ví dụ 2: Xác định chuyển vị dầm thẳng, tựa trên hai gối đặt hai đầu. Dầm dài L, độ cứng EI, chịu tải phân bố đều q = const. Đây là bài tốn quen thuộc từ sức bền vật liệu, lời giải giải tích đã rõ ràng. Dầm sẽ mơ hình hĩa theo 2 cách: a) xét dầm như hệ thống liên tục, b) lợi dụng tính chất đối xứng cấu hình, tải, chỉ giữ lại phần đối xứng, tức một nửa chiều dài dầm để xem xét và xử lý. a) Cách giải khơng để ý đến tính đối xứng: chia dầm làm 2 phần tử, đánh số 1 và 2 dài bằng nhau, bằng L/2, ba nút 1, 2, 3 ghi lại như tại hình 3.13. Sử dụng những cơng thức xác định độ cứng phần tử, vector tải cho phần tử BEAM 2D nêu tiếp tại điểm 8 chương này để tính tiếp sau đây. Hình 3.13 Dầm thẳng và mơ hình phân tử hữu hạn Ma trận cứng phần tử 1 và 2: ⎡ 12 3L − 12 3L ⎤ ⎢ 2 2 ⎥ 8EI 3LLLL− 3 0,5 []k= [] k = ⎢ ⎥ ( ) a (1) (2) L3 ⎢−12 − 3L 12− 3L ⎥ ⎢ 2 2 ⎥ ⎣ 3LLLL 0,5− 3 ⎦ Ma trận cứng tồn hệ thống: ⎡12 3L − 12 3L 0 0 ⎤ ⎢ 2 2 ⎥ ⎢ LLL− 3 / 2 0 0 ⎥ 8EI ⎢ 24 0− 12 3L ⎥ []K = 3 ⎢ 2 2 ⎥ ( ) b L ⎢ 2LLL− 3 0,5 ⎥ ⎢ ĐX 12− 3L ⎥ ⎢ 2 ⎥ ⎣⎢ L ⎦⎥ Vector tải hệ thống:
  72. 71 T ⎡ qL qL2 qL qL qL2 ⎤ {}P =⎢ = − − −0 − ⎥ (c) ⎣ 48 48 2 4 48 ⎦ Điều kiện biên: v1 = v3 = 0. (d) Sau xử lý điều kiện biên, hệ phương trình chứa 4 ẩn cĩ dạng sau: 2 2 2 ⎡ LLL− 3 0,5 0 ⎤⎧θ1 ⎫ ⎧− qL / 48⎫ ⎢ ⎥⎪ ⎪ ⎪ ⎪ 8EI 24 0 3L ⎪v 2 ⎪ ⎪ − qL 2/ ⎪ ⎢ ⎥⎨ ⎬ = ⎨ ⎬ (e) L3 ⎢ 2LL2 0,5 2 ⎥ θ 0 ⎢ ⎥⎪ 2 ⎪ ⎪ ⎪ 2 ⎪ ⎪ ⎪ 2 ⎪ ⎣⎢ĐX L ⎦⎥⎩θ3 ⎭ ⎩ qL / 48 ⎭ 3 4 Sau giải xác định được: θ1 = -θ3 = -qL /(24EI); v2 = -5qL /(384EI) và θ2 = 0. b) Sử dụng tính đối xứng hình học và tải, chỉ chọn ½ dầm cho mơ hình tính. Điều kiện biên tại tâm đối xứng thể hiện như sau: chuyển vị v2 ≠ 0, đạt maximum; gĩc xoay dầm tại nút 2, θ2 tính bằng đạo hàm chuyển vị v’x=L/2 = 0. Sử dụng 1 phần tử dầm 2D vào tính tốn. Ma trận phần tử (a) tiếp tục được dùng, vector tải tính cho một phần tử cĩ dạng: T ⎡ qL qL2 qL qL2 ⎤ {}P =⎢ − − − ⎥ (a) ⎣ 4 48 4 48 ⎦ Xử lý điều kiện biên v1 = 0, θ2 = 0 sẽ đưa hệ phương trình đại số tới dạng: 2 2 8EI ⎡ LL− 3 ⎤⎧θ1 ⎫ ⎧− qL / 48⎫ 3 ⎢ ⎥⎨ ⎬ = ⎨ ⎬ (b) L ⎣− 3L 12 ⎦⎩v 2 ⎭ ⎩ − qL 4/ ⎭ Sau giải hệ phương trình sẽ nhận kết quả: 3 4 θ1 = -qL /(24EI); v2 = -5qL /(384EI) Cách làm b) cho phép giảm bậc tự do bài tốn, theo đĩ tiết kiệm thời gian tính tốn trên cơ sở đảm bảo độ chính xác phép tính. 7. NHỮNG PHẦN TỬ THƠNG DỤNG DÙNG TRONG CƠ HỌC KẾT CẤU Phần tử 1 chiều (1D) từ hai đến nhiều nút, hình 3.14a. Trong thực tế sử dụng chúng ta cĩ thể dùng phần tử tuyến tính hoặc phần tử cong. Phần tử phẳng (2D) thường dạng hình tam giác hoặc tứ giác, cạnh thẳng hoặc cong. Nút phần tử dạng này nằm tại đỉnh và cả tại các cạnh, xem hình 3.14b. Nhĩm phần tử 3D cĩ thể chia làm ba phân nhĩm, phần tử hình thù giống kim tự tháp gồm bốn mặt, 4 nút hoặc nhiều hơn 4 nút thuộc phân nhĩm đầu tiên, từ kỹ thuật các nước đang dùng đều mang âm hưởng tiếng La tinh tétrahédrique chỉ phần tử này. Phân nhĩm thứ hai trong cụm này là khối 6 mặt (hexahédrique) cĩ từ 8 nút trở lên. Trong tài liệu này chúng ta làm quen với phần tử cĩ số nút đến 20. Phân nhĩm thứ ba đề cập phần tử cĩ mặt cắt ngang khơng đổi, gọi là khối lăng trụ, hình 3.14c.
  73. 72 a) 2 nút 3 nút 4 nút b) 3 nút 9 nút 6 nút c) 12 nút 4 nút 8 nút d) 10 nút 8 nút 4 nút 20 nút Hình 3.14 Phần tử thường dùng Những phần tử nhĩm a) từ hình 3.14 gồm phần tử chỉ chịu lực dọc trục cĩ tên gọi phổ thơng trong các sách chuyên ngành là ROD, BAR, TRUSS, phần tử chịu xoắn TOR, phần tử chịu uốn BEAM. 7.1 Phần tử BAR hay cịn gọi phần tử TRUSS, ROD Xây dựng ma trận cứng phần tử cĩ thể tiến hành theo cách trực tiếp. Phần tử 1-2 dài L, diện tích mặt cắt ngang A, làm từ vật liệu mơ đun đàn hồi E, chịu tác động lực dọc trục F. Chuyển vị nút tính theo cơng thức δ = FL/(AE). Hãy ký hiệu chuyển vị nút 1 là u1, nút 2 là u2. Lực tác động lên phần tử F tính từ biểu thức F = (AE/L)δ. Trường hợp cụ thể này: AE AE FF= = u F= F = u và FFFFFF= − = − + 11 21 L 1 12 22 L 2 1 11 12 2 21 22 AE ⎡ 1− 1⎤⎧u1 ⎫ ⎧F1 ⎫ Hay là ⎢ ⎥⎨ ⎬ = ⎨ ⎬ L ⎣−1 1 ⎦⎩u2 ⎭ ⎩F2 ⎭ Ma trận cứng ke cĩ dạng: EA ⎡ 11− ⎤ [k] = ⎢ ⎥ L ⎣−11⎦ Vector lực: L T ⎧pL / 2⎫ {F} = ∫ [N] {p} dx = ⎨ ⎬ (l) 0 ⎩pL / 2⎭ Dưới tác động lực tập trung P: ⎧ PL()− ξ ⎫ L ⎪ ⎪ {F} = P δ(ξ-x) [N]T dx = L Hình 3.15 Phần tử BAR (TRUSS) ∫ ⎨ Pξ ⎬ 0 ⎪ ⎪ ⎩ L ⎭ Cách tính mang tính hình thức hiểu theo cách sau. Để xác định ma trận [B] trong cơng thức chung: T [k ]e = [ B ] [ D ][ B ] dV ∫V e cần thiết xây dựng quan hệ {u} = {N]{δ}. Trường hợp này nên nhận hàm hình dáng tuyến tính:
  74. 73 ⎡ L− x x ⎤⎧u1 ⎫ {}u = ⎢ ⎥⎨ ⎬ ⎣ L L⎦⎩u2 ⎭ Hàm biến dạng: d ⎡ d ⎤ {}ε ={}u = ⎢ NB⎥{}δ= []{} δ dx ⎣dx ⎦ ⎡ 1 1 ⎤ trong đĩ []B =⎢ − ⎣ LL⎦⎥ Thay [B] vào cơng thức tính [k]: L ⎧−1/ L⎫ ⎡ 1 1 ⎤ AE ⎡ 1− 1⎤ []k = ∫ ⎨ ⎬E⎢− ⎥ Adx = ⎢ ⎥ 0 ⎩ 1/ L ⎭ ⎣ LL⎦ L ⎣−1 1 ⎦ 7.2 Phần tử TRUSS trong khơng gian 3D Chuyển vị nút phần tử hai nút, xét trong hệ tọa độ cục bộ gắn liền phần tử: T {δ}cục_bộ = [ u1 v1 w1 | u2 v2 w2 ] Trong hệ tọa độ chung vector này cĩ dạng: T {δ}chung = [ U1 V1 W1 | U2 V2 W2 ] Quan hệ giữa hai vector này thể hiện bằng biểu thức: ⎧u ⎫ ⎧U ⎫ ⎪ ⎪ ⎪ ⎪ ⎨v⎬ = []R ⎨V ⎬ ⎪ ⎪ ⎪ ⎪ ⎩w⎭ ⎩W ⎭ trong đĩ ma trận chuyển [R] được hiểu như sau: ⎡cos()x , X cos () x , Y cos( x , Z )⎤ ⎢ ⎥ []R = ⎢cos()y , X cos () y , Y cos () y , Z ⎥ ⎣⎢cos()z , X cos () z , Y cos () z , Z ⎦⎥ Ma trận cứng tính trong hệ tọa độ cục bộ, cho phần tử 3D BAR cĩ dạng: ⎡1 0 0− 1 0 0⎤ ⎢ ⎥ ⎢ 0 0 0 0 0⎥ EA ⎢ 0 0 0 0⎥ []k cb = ⎢ ⎥ L ⎢ 1 0 0⎥ ⎢ ĐX 0 0⎥ ⎢ ⎥ ⎣⎢ 1⎦⎥ T Trong hệ tọa độ chung ma trận cứng phần tử tính bằng {R] [k]cb[R]. Vector lực tác động phần tử T tính chuyển sang hệ tọa độ chung bằng biểu thức [R] {p}cb. Ví dụ 3: Sử dụng phần tử dầm chịu kéo, nén (BAR element) giải khung sau: Chiều dài các thanh, qui đổi ra đơn vị l = l23 = 3m: l23 2l l23 l12 = =;l24 = = 2l cosα 3 cos β Ma trận các phần tử: Phần tử số 1 Hình 3.16