Giáo trình Cơ học kết cấu - Chương 6: Các phương pháp tính

n n ⎭⎬ ⎫⎥⎦ ⎤⎢⎣ ⎡ ⎟⎠ ⎞⎜⎝ ⎛∑∞ = ππππ coscos 1 = 0. (e) Nhờ tính trực giao của hệ hàm L xkπcos và L xkπsin , k =1,2, 3,... trong đoạn (0, L): =∫L dxLxkLxn0 coscos ππ 0 nếu k≠ n = 2 L nếu k = n; Thay biểu thức cuối vào (e ) kết quả sẽ được: 171 0 2 = ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ −⎟⎠ ⎞⎜⎝ ⎛=Π T L kEIa a kk π ∂ ∂ (f) Phương trình trên, 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: 2 22 L EIkT π= (g) Chúng ta quan tâm đến giá trị nhỏ nhất của T, 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, T đạt nhỏ nhất cho trường hợp k =1, do vậy biểu thức lực Euler có dạng: 2 2 L EITE π= (h) Kết quả trên đây phù hợp với kết quả tính trong các chương trước. Ví dụ 3: Tìm độ võng tấm và ứng lực trong tấm bị uốn. Trường hợp tổng quát, chuyển vị theo hướng 0z của tấm được tìm dưới dạng: (x,y) (a) n n n fayxuyxw ∑∞ = == 1 ),(),( Để áp dụng phương pháp R-R 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 toán cụ thể chính là hàm tổng năng lượng của tấm. Từ lý thuyết đàn hồi chúng ta đã biết, thế năng biến dạng được thể hiện bằng công thức: U = 1 2 V∫∫∫ (σxεx + σyεy + σzεz + τxyγxy + τxzγxz + τyzγyz)dxdydz (b) Với tấm bằng vật liệu trực hướng, γxz = γyz = 0; σz = 0; do vậy phương trình U chỉ còn các thành phần sau: U = 1 2 V∫∫∫ (σxεx + σyεy + τxyγxy )dxdydz (c) Các hàm σx,σy, τxy được trình bày trong quan hệ với biến dạng trong trạng thái ứng suất phẳng, trong lý thuyết đàn hồi. Thay tất cả quan hệ đó vào công thức tính U, sẽ nhận được biểu thức đầy đủ: 172 U = 1 2 V∫∫∫ [ ] [ ]E E Gx x x y y y x y yx1 112 21 0 21 0 0 12 21 0 21 0 0 12 0−⎧⎨⎩ + + − + +ν ν ε ν ε ε ν ν ε ν ε ε γ − 2z E w x w y w x E w x w y w x Gx x x y x x x y1 1 2 1 2 1 1 2 1 2 2 12 21 0 2 2 21 0 2 2 21 0 2 2 12 21 0 2 2 21 0 2 2 21 0 2 2 12 0 − + + ⎛ ⎝⎜⎜ ⎞ ⎠⎟⎟ + − + + ⎛ ⎝⎜⎜ ⎞ ⎠⎟⎟ + ⎡ ⎣ ⎢⎢ ⎤ ⎦ ⎥⎥ν ν ε ∂ ∂ ν ε ∂ ∂ ν ε ∂ ∂ ν ν ε ∂ ∂ ν ε ∂ ∂ ν ε ∂ ∂ γ xy + z2 E w x w x w y x 1 12 21 2 2 2 2 2− + ⎛ ⎝ ⎜⎜ ⎞ ⎠ ⎟⎟ + ⎡ ⎣ ⎢⎢ ν ν ∂ ∂ ν ∂ ∂ ∂ ∂ E w y w x w y G w x y1 4 2 2 2 2 2 12 2 2 − + ⎛ ⎝ ⎜⎜ ⎞ ⎠ ⎟⎟ + ⎛ ⎝⎜ ⎞ ⎠⎟ ⎤ ⎦⎥ ⎫⎬⎭νν ∂ ∂ ν ∂ ∂ ∂ ∂ ∂ ∂ ∂ dxdydz (c*) Tiến hành tích phân vế phải phương trình trên theo hướng trục 0z, từ -0,5t đến +0,5t trong đó các tham số phụ thuộc vào z, còn w không phụ thuộc z, kết quả tính sẽ như sau: ε ε γx y xy0 0 0, , U= 1 2 t∫∫ ( ) ( )E Ex x x y y y x1 112 21 0 21 0 0 12 21 0 12 0 0− +⎧⎨⎩ + − +ν ν ε ν ε ε ν ν ε ν ε ε y + Gxyγxy0 ⎫⎬⎭ dxdy+ + 1 2 t∫∫ E t wx wx wy E t wx wx wyx y 3 12 21 2 2 2 21 2 2 2 2 3 12 21 2 2 2 12 2 2 2 212 1 12 1( ) ( )− ⎛ ⎝⎜ ⎞ ⎠⎟ + ⎡ ⎣ ⎢⎢ ⎤ ⎦ ⎥⎥ ⎧ ⎨⎪ ⎩⎪ + − ⎛ ⎝⎜ ⎞ ⎠⎟ + ⎡ ⎣ ⎢⎢ ⎤ ⎦ ⎥⎥ν ν ∂ ∂ ν ∂ ∂ ∂ ∂ ν ν ∂ ∂ ν ∂ ∂ ∂ ∂ + + G t w x y dxdy12 3 2 2 3 ∂ ∂ ∂ ⎛ ⎝⎜ ⎞ ⎠⎟ ⎫⎬⎭ (d) Tích phân trên tiến hành cho toàn bộ diện tích tấm. Nếu ký hiệu độ cứng tấm theo hai hướng 0x và 0y làD1 và D2: D1 = )211(12 12 3 νν− tEx , D2 = )211(12 12 3 νν− tE y ; và DT = 12 3 12tG phương trình thế năng biến dạng tấm, tính từ vế sau của phương trình trên đây có dạng: ∫∫ ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛+⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛++⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛= dxdy yx wD y wD y w x wD x wDtU T 22 2 2 22 2 2 2 2112 2 1 422 1 ∂∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ν∂ ∂ (e) Trường hợp tấm đẳng hướng, Ex = Ey = E; ν12 = ν21 = ν; )1(2 ν−= EG Hàm thế năng có dạng: 173 U = ( ) dxdy y w x w yx wwD ∫∫ ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ −⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛−+∇ 2 2 2 222 22 )1(2 2 ∂ ∂ ∂ ∂ ∂∂ ∂ν (f) trong đó D = )1(12 2 3 ν− Et (g) Công ngoại lực trong trường hợp này có dạng: W = q(x,y).w(x,y)dxdy (h) ∫∫ Tích phân được tiến hành bao quát bề mặt tấm. Áp dụng cách giải trên để xác định mặt võng tấm hình chữ nhật có các cạnh axb, chịu tác động lực tập trung đặt theo phương pháp tuyến tại điểm giữa tấm x0, y0. Tấm được ngàm cả bốn cạnh, chiều dầy tấm t = const. Hàm chuyển vị mặt qua giữa tấm, tính theo trục z, dạng công thức (a) nêu trên được viết theo cách diễn đạt của Lévy: w(x,y) = u(x,y) = a m x a n y bmnnm = ∞ = ∞ ∑∑ 11 sin sinπ π (a’) Với các tấm kiểu Timoshenko, độ võng w không lớn nếu so với kích thước tấm, các hàm năng lượng U, W tính cho tấm diễn đạt bằng (f) và (h) có thể viết: ( ) dxdy y u x u yx u y u x uDU a b∫ ∫ ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ −⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛−+⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ += 0 0 2 2 2 222 2 2 2 2 2 12 2 ∂ ∂ ∂ ∂ ∂∂ ∂ν∂ ∂ ∂ ∂ (b’) Sau khi thay giá trị các đạo hàm riêng tính từ (a’) vào (b’), biểu thức (b’) trở thành: dxdy b yn a xm b n a maDU a b m n mn∫ ∫ ∑∑ ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ⎟⎠ ⎞⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛= ∞ = ∞ =0 0 1 1 222 sinsin 2 ππππ (c’) Công do lực tập trung tác động: b yn a xmaPuPW m n mn ππ sinsin. 1 1 ∑∑∞ = ∞ = == (d’) Trong công thức trên đây )1(12 2 3 ν−= EtD là độ cứng tấm, t - chiều dầy, E - mođun đàn hồi vật liệu, ν - hệ số Poisson. 174 Phiếm hàm I được tìm dưới dạng I ≡ Π = U-W như đã làm cho ví dụ 1. Tiến hành lấy đạo hàm của (U - W) theo amn: ∂ ∂ (U W amn )− sẽ nhận được hệ phương trình đại số, trong đó amn đóng vai trò ẩn số. amn = Dab 4 [(mπ/a)2 + (nπ/b)2 ] 2 - Psin(mπx0/a) sin(nπy0/b) = 0 với m = 1,2,...; n = 1,2,... Thay amn từ (e’) vào (a’), hàm w = u(x,y) có dạng: u(x,y) = 4P Dab nm = ∞ = ∞ ∑∑ 11 sin sin m x a n y b m a n b π π π π 0 0 2 2⎛ ⎝⎜ ⎞ ⎠⎟ + ⎛ ⎝⎜ ⎞ ⎠⎟ ⎡ ⎣ ⎢⎢ ⎤ ⎦ ⎥⎥ sin m x a n y b π πsin (f’) 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 R-R 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ị tổng quát nhằm xác định lực tới hạn. Với trường hợp tấm cứng hình chữ nhật cạnh axb, mép tấm bị giữ chặt, không được phép dịch chuyển trong mặt x0y, cách giải sẽ như sau. Tìm hàm chuyển vị của tấm theo trục 0z: w(x,y) = a nm = ∞ = ∞ ∑∑ 11 mn ϕm(x) Ψn(y) (a) 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 cũng cùng tấm U = U1 + U2. Trong đó U1 tính theo công thức (f) từ ví dụ trước: ( ) dxdy y w x w yx wwDU ∫∫ ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ −⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛−+∇= 2 2 2 222 22 1 )1(22 ∂ ∂ ∂ ∂ ∂∂ ∂ν (b) U2 tính theo công thức: dxdy y w x wT y wT x wTU xyyx∫∫ ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ +⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛= ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ 2 2 1 22 2 (c) trong đó, từ lý thuyết tấm xác định được: ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ +−= y v x uEtTx ∂ ∂ν∂ ∂ ν 00 21 và ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ +−= x u y vEtTy ∂ ∂ν∂ ∂ ν 00 21 175 ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ +== x v y uGtTT yxxy ∂ ∂ ∂ ∂ 00 Trong trường hợp Tx = -σxt; Ty = Txy = 0, hàm U2 trở thành: dxdy x wTU a b x 2 0 0 2 2 1 ⎟⎠ ⎞⎜⎝ ⎛= ∫ ∫ ∂∂ (d) Quay trở lại điều kiện biên, có thể viết: tại x = 0 và x = a: w = 0; ∂∂ w x ≠ 0; Tại y =: w = 0; ∂∂ w y ≠ 0; Tại y = b: w ≠ 0; ∂∂ w y ≠ 0; (e) Từ điều kiện biên có thể nhận thấy rằng, hàm w(x,y) phù hợp cho các trường hợp vừa nêu nếu được viết dưới dạng: w(x,y) = a nm = ∞ = ∞ ∑∑ 11 mn y b m x a n⎛ ⎝⎜ ⎞ ⎠⎟ sin π (f) Để nhanh chóng có kết quả, chưa quan tâm nhiều đến độ chính xác có thể chỉ chọn 1 hằng số khi tính chuỗi cho w. w(x,y) = a11. y b x a sin π (g) 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: ( ) dxdy y w x w yx wwDU ∫∫ ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ −⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛−+∇= 2 2 2 222 22 1 )1(22 ∂ ∂ ∂ ∂ ∂∂ ∂ν = = ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ⎟⎠ ⎞⎜⎝ ⎛−+⎟⎠ ⎞⎜⎝ ⎛ 2 1)1(2 62 242 11 a ab ab a Da πνπ (h) U2 = - 1 2 11 2a .σxt( πa ) 2 ab 6 (i) Tiến hành tính đạo hàm theo biểu thức: 0 )( 21 =+ mna UU ∂ ∂ (j) sẽ nhận được phương trình: 176 a11 D a ab b a a t a ab x π ν π σ π⎛⎝⎜ ⎞ ⎠⎟ + − ⎛ ⎝⎜ ⎞ ⎠⎟ ⎡ ⎣ ⎢⎢ ⎤ ⎦ ⎥⎥ − ⎛ ⎝⎜ ⎞ ⎠⎟ ⎧ ⎨⎪ ⎩⎪ ⎫ ⎬⎪ ⎭⎪ 4 2 6 2 1 1 2 6 ( ) 2 = 0; (k) Từ tính chất cơ học của bài toán, hằng số a11 ≠ 0, do đó để phương trình trên đạt giá trị 0, biểu thức trong dấu ngoặc lớn phải bằng 0, từ đó xác định lực Euler theo công thức: ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ⎟⎠ ⎞⎜⎝ ⎛+−= ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ −+⎟⎠ ⎞⎜⎝ ⎛== 2 22 2 2 2 )1(66)1( a b tb D aat D Ex π νπνπσσ (l) 2. Phương pháp Galerkin Phương pháp do Galerkin1 đề xướng từ 1915, song từ những năm ba mươi trở đi ứng dụng mới thực sự rộng rãi sau khi tác giả phải trình bày kỹ tại Viện Hàn lâm về ứng dụng của phương pháp kể cả trong trường hợp không gia 3D. Phương pháp Galerkin thích hợp cho những bài toá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ếu ký hiệu L() - toán tử tuyến tính, B() - toán tử vi phân, bài toán kỹ thuật tổng quát có thể viết dưới dạng phương trình: L(u) - p = 0 trong miền V, (6.4) và các điều kiện biên: B(u) -q = 0 trên biên S (6.5) Để giải bài toán có thể tìm hàm u dưới dạng hàm xấp xỉ: u* = (6.1’) a fi i i N = ∑ 1 Hàm u* được tìm phải thỏa mãn điều kiện B(u*) = q(s). Với số lượng i hữu hạn, thông thường khi thay (6.1) vào bài toán (6.4), không thể 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 số, ký hiệu R hoặc ε: ε = L( ) - p ≠ 0. a fi i i N = ∑ 1 Cách làm của Galerkin là đặt ε trực giao với hệ hàm cơ sở fi, từ đó có thể nhận tích vô hướng dạng sau: = 0, i=1,2,.... (6.6) Dưới dạng đầy đủ công thức cuối được viết lại: 1 Viện sĩ hàn lâm khoa học USSR, trung tướng công binh. 177 0= ⎭⎬ ⎫ ⎩⎨ ⎧ −⎟⎠ ⎞⎜⎝ ⎛∫ ∑ dVfpfa i i iiL , i=1,2,.... (6.7) Ví dụ 1:Giải phương trình vi phân: L(u) - p = d u dx 2 2 + u + x = 0 trong miền 0 < x <1 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ố: ε = L(u) - p = x( -2 +x -x2)a1 + (2 - 6x + x2 - x3 )a2 Thay các hàm u và ε vào phương trình Galerkin (6.7), chúng ta nhận được hệ phương trình đại số: 0 1 ∫ ε.x(1-x)dx =0 0 1 ∫ ε.x2(1-x)dx = 0 (c) Giải hai tích phân trên, hệ phương trình đại số sẽ có dạng: 3 10 3 20 3 20 13 105 1 12 1 20 1 2 ⎡ ⎣ ⎢⎢⎢ ⎤ ⎦ ⎥⎥⎥ ⎧⎨⎩ ⎫⎬⎭ = ⎧ ⎨⎪ ⎩⎪ ⎫ ⎬⎪ ⎭⎪ a a (d) Sau khi giải hệ phương trình, các hệ số được xác định như sau: a1 = 71 369 và a2 = 7 41 Từ đó hàm xấp xỉ u sẽ là: u = x(1-x)( 71 369 + 7 41 x ) (e) 178 Phương pháp Galerkin được tổng quát hoá, sử dụng rộng rãi như một phép biến phân. Theo hướng này biểu thức (6.7) có thể viết dưới dạng hệ hàm trực giao hệ thống hàm fi với (6.4) như sau: ∫ { L(u) - p} fidV = 0, i=1,2,3,...,N (6.7a) Mặt khác phân sai của u có thể khai triển dạng chuỗi: δu = δa1f1 + δa2f2 +....... + δaNfN (6.8) δu = ∑δaifi Thay biểu thức cuối vào (6.7a) sẽ nhận được hệ phương trình: ∫ { L(u) - p} δu dV = 0 (6.9) Ví dụ 2: Xác định độ võng dầm như đã nêu tại ví dụ trình bày trong phương pháp R-R, phần trước. Phương trình vi phân miêu tả độ võng dầm có dạng: EJ.w(4) (x) - q = 0 (a) Đ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) phần trước: w(x) = a i= ∞∑ 1 i sin i x L π Áp dụng công thức Galerkin vào bài toán sẽ nhận được hệ phương trình: 0 L ∫ [ EJ.w(4) (x) - q] sin i xLπ dx = 0, i=1,2,...,∞ 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. w(x) = 4 1 4 5 5 1 3 qL EJ kkπ = ∞∑ , ,... sin k x 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ụ 3: Áp dụng công thức tổng quát của Galerkin (6.9) 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: ∂ ∂ ∂ ∂ 2 2 2 2 u x u y C+ = (a) Điều kiện biên: 179 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 (6.9) chúng ta nhận được: ∂ ∂ ∂ ∂ 2 2 2 2 00 u x u y C ba + −⎛⎝⎜⎜ ⎞ ⎠⎟⎟∫∫ δudxdy = 0 (d) Trong đó: δu = δα x(x-a)y(y-b) (e) Sau khi thay thế (c) vào (d) sẽ được: ( )[ ]α90 363 3 2 2 3 3a b a b Ca b+ − = 0. từ đó: α = 5 2 2 2 C a b+ (f) và hàm u: u = 52 2 C a b+ (x 2 - ax) (y2 - by) (g) Ví dụ 4: Áp dụng công thức Galerkin (6.9) giải bài toán uốn tấm có các cạnh axb, ngàm tại các mép tấm. Độ cứng tấm D, chiều dầy tấm t. Tấm chịu tải trọng phân bố đều p(x,y) tác động theo hướng pháp tuyến. 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ỉ: u(x,y) = w(x,y) = c i= ∞∑ 1 ifi(x,y) (a) 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à: u(x,y) = a nm = ∞ = ∞ ∑∑ 11 mn. 1 4 (1-cos 2m x a π ) ( 1- cos 2n y b π ) (c) 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) 180 Nếu chỉ chọn m =1 và n=1, hàm u(x,y) sẽ có dạng: u(x,y) = a11. 1 4 (1-cos 2πx a )(1-cos 2πy b ) (e) Sau khi thay thế giá trị u từ (e) vào phương trình Galerkin (6.9): ∫[D.∇2 ∇2 u - p(x,y)] δudxdy = 0 kết quảsẽ nhận được: 00 ba ∫∫ [D.a11.∇2 ∇2 f1(x,y) - p(x,y)] fi(x,y)dxdy = 0. (f) Thực hiện phép lấy đạo hàm riêng ∇2 ∇2 f1(x,y) = 4 4 2 2 π a b [ (1-cos 2π b y)cos 2π a x +cos 2π a x cos 2π b y +(1-cos 2π a x)cos 2π b y ] (g) và thay thế (g) vào (f), hệ số đầu tiên của chuỗi có dạng: a11 = pa b D 2 2 48 π 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: umax = 0,00128 pa D 4 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 = 0,00126 pa D 4 . 3. Phương pháp tính thỏa mãn trên một số điểm hoặc miền được chọn Để giải phương trình: L(u) - p = 0 trong miền V, (6.4) và các điều kiện biên: B(u) -q = 0 trên biên S (6.5) có thể tìm hàm u dưới dạng hàm xấp xỉ: u* = (6.1) a fi i i N = ∑ 1 và hàm u được tìm phải thỏa mãn điều kiện Bk(u) = q(s), k =1,2,... Trong trường hợp bài toán có lời giải chính xác u = u0, phương trình (6.4) và (6.5) sẽ đúng: L(u0) = p trong V và B(u0) = q trên S. 181 Trường hợp u = u0 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 toá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. Ngoài các điểm được chọn những điều kiện đặt ra trong bài toán không chắc được thỏa mãn. Thứ tự giải bài toán như sau: - Tìm hàm xấp xỉ thỏa mãn điều kiện biên như đã làm theo (6.1): u* = a fi i i N = ∑ 1 - Xác lập hàm sai số: ε = L(u*) - p = a i= ∑ 1 Ν i L(fi) - p. 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 (6.10) trong đó hàm Dirac δ được định nghĩa: ξ ξ − + ∫ c c δ(x-ξ)dx = 1 với c → 0, ngoài phạm vi trên hàm δ(x) = 0. Ví dụ 1: Xác định nghiệm phương trình vi phân nêu tại ví dụ trước: L(u) - p = d u dx 2 2 + u + x = 0 trong miền 0 < x <1 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ố ε: ε = L(u) -p = x + ( -2+x- x2)a1 + ( 2-6x +x2 -x3)a2. Tại x =0,25 và x =0,5 các hệ số a1, a2 phải thỏa mãn: 182 29 16 35 64 7 4 7 8 1 4 1 2 1 2 −⎡ ⎣ ⎢⎢⎢ ⎤ ⎦ ⎥⎥⎥ ⎧⎨⎩ ⎫⎬⎭ = ⎧ ⎨⎪ ⎩⎪ ⎫ ⎬⎪ ⎭⎪ a a Từ đó: a1 = 6 31 và a2 = 40 217 Hàm u(x) tính theo phương pháp chọn điểm: u(x) = 1 217 x(x-1)(42+40x). Ví dụ 2: Xác định độ võng dầm liên tục dài L, EJ = const, chịu tác động tải trọng phân bố đều p = cont. 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 theo độ võng dầm: EJw(4) (x) - p = 0 (a) Đ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: w(x) = a k N = ∑ 0 k x L k⎛ ⎝⎜ ⎞ ⎠⎟ (c) 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 toá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à: 6 6 5 5 4 4 3 3 2 210)( ⎟⎠ ⎞⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛+= L xa L xa L xa L xa L xa L xaaxw (d) 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 (6.10 ), hệ phương trình chứa ai sẽ là: 183 24a4 + 30a5 + 22,5a6 = EJ pL4 24a4 + 60a5 + 90a6 = EJ pL4 24a4 + 90a5 + 202,5a6 = EJ pL4 Từ đó: a0 = 0; a1 = 0; a2 = EJ pL 16 ; a3 = - EJ pL 48 5 4 ; a4 = EJ pL 24 4 ; a5 = 0; a6 = 0. Phương trình độ võng dầm theo cách tính trên: ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ⎟⎠ ⎞⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛−⎟⎠ ⎞⎜⎝ ⎛= 4324 24 1 48 5 16 1)( L x L x L x EJ pLxw 4. 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 chứa đạo hàm riêng. Phương pháp dùng cho bài toán một chiều bằng cách chia đoạn thẳng đang xem xét thành nhiều bước, gọi là sai phân, với tổng số bước hữu hạn. Trong bài toá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ự. Chuỗi Taylor. Với các bài toá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í x0 + Δx dựa vào chuỗi Taylor, dùng trong trường hợp tiến: y(x0 + Δx) = y(x0) + Δx1! y’(x0) + Δx 2 2! y’’(x0) + Δx 3 3! y’’’(x0) +.... (6.11) Khi tính thụt lùi chuỗi Taylor có dạng: y(x0 - Δx) = y(x0) - Δx1! y’(x0) + Δx 2 2! y’’(x0) - Δx 3 3! y’’’(x0) +.... (6.12) Nếu chỉ lấy hai thành phần đầu của chuỗi, khi trừ (6.11) với (6.12) chúng ta nhận được: y(x0 + Δx) - y(x0 - Δx) = 12Δx [y’(x0)Δx - (-y’(x0)Δx) ] từ đó công thức tính đạo hàm bậc 1, tại x0 sẽ là: 184 y’(x0) = 1 2Δx . [y(x0 + Δx) - y(x0 - Δx) ] (6.13) 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(x0 + Δx) - y(x0 - Δx) = 2y(x0) + y’’(x0). (Δx)2 từ đó: y’’(x0) = 1 2( )Δx .[ y(x0 + Δx) - 2y(x0) + y(x0 - Δx) ] (6.14) Công thức đơn giản 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: y’k = 1 2Δx (yk+1 - yk-1) y”k = 1 2 2Δx (yk+1 -2yk + yk-1) y”’k = 1 2 3Δx (yk+2 - 2yk+1 + 2yk-1 - yk-2 ) y””k = 1 2 4Δx ( yk+2 - 4yk+1 +6yk - 4yk-1 + yk-2) (6.15) Aùp dụng phân tích chuỗi Taylor với hàm y = f(x) xác định trong miền 0 -L, có thể tính giá trị hàm tại xk+1, k = 0,1,2,..., dựa vào các biểu thức y’, y’’,... x0 x1 x2 xL Hình 2 Sử dụng các biểu thức ghi trong (6.15) giải bài toán trong không gian một chiều được giới thiệu trong ví dụ 1 dưới đây. Ví dụ1: Xác định độ võng dầm consol dài L, độ cứng EJ=const, chịu tác động lực phân bố đều q. Dầm bị ngàm bên trái, và tự do phía phải. Phương trình vi phân của độ võng dầm: EJw (IV) (x) = q trong miền 0 < x < L Điều kiện biên: tại x = 0: w = 0; w’ = 0; tại x = L: w’’ = 0; w’’’ = 0; Chia chiều dài L dầm làm 5 đoạn bằng nhau với bước h= Δx = L/5. Số thứ tự các nút được đánh từ trái sang phải, bắt đều từ nút đầu phía trái là 0, sau đó 1,2,3,4,5. 185 Ngoài đoạn L, các nút giả được đánh số như sau: trước nút 0 thêm nút -1, cách nút 0 đoạn h; sao nút 5 tiếp tục đánh số nút 6 và 7 về nên phải dầm cách nút thứ 5 đoạn tương ứng h và 2h. Đạo hàm bậc 4 xuất hiện trong phương trình cân bằng được tìm theo phương pháp sai phân hữu hạn có dạng: w(IV) (x) = 14h (wk+2 - 4wk+1 +6wk - 4wk-1 + wk-2) trong đó: tại x =0: w0 = 0; góc uốn: -w-1 + w1 = 0; còn tại x = L điều kiện biên có dạng: w4 -2w5 + w6 =0; -w3 + 2w4 -2w6 + w7 =0; Với các nút k= 1,2,3,4,5 phương trình cân bằng sau khai triển có dạng: wk+2 - 4wk+1 +6wk - 4wk-1 + wk-2 = qh EJ 4 Giải hệ phương trình đại số tuyến tính,trong đó các giá trị wk đóng vai trò ẩn số, sẽ xác định độ võng của dầm. 4.1 Phương pháp lưới cho bài toá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 bằng i, còn dọc Oy ký hiệu k. Theo cách qui ước này, công thức tính đạo hàm riêng ∂∂ ∂ m n m nx y u + sẽ như sau: Tại điểm x = xk, y = yi, đạo hàm riêng bậc từ 1 đến 4 có dạng: ∂ ∂ u x u u x ik i k i k= −+ −, ,1 1 2Δ ; ∂ ∂ u y u u y ik i k i k= −+ −1 1 2 , , Δ ∂ ∂ 2 2x uik = u u u x i k ik i k, ( ) + ,− +1 1 2 2 Δ − ; uik = u u u y i k ik i k+ −− +1 1 2 2, , ( )Δ ∂ ∂ 2 2y ∂ ∂ ∂ 2 x y uik = u u u u x y i k i k i k i k+ + + − − + − −− − +1 1 1 1 1 1 1 1 4 , , , , Δ Δ ∂ ∂ 4 4x uik = u u u u u x i k i k ik i k i k, , , , ( ) + + −− −+ − +2 1 1 4 4 6 4 Δ 2 186 ∂ ∂ 4 4y uik = u u u u u y i k i k ik i k i k+ + − −− + − +2 1 1 4 4 6 4, , , ( )Δ 2, ∂ ∂ ∂ 4 2 2x y uik = [1 2 1 1 1 1 1 1 1 1( ) , , , ,Δ Δx y u u u ui k i k i k i k+ + + − − + − −+ + + + (6.16) ]+ − + + ++ − + −4 2 1 1 1 1u u u u uik i k i k i k i k( , , , , ) Ví dụ 1: 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 2 tựa trên gối, dọc y= ± a 2 cạnh bị ngàm. Tính ứng suất trong tấm. Phương trình chung uốn tấm: 4 4 x w ∂ ∂ + 2 22 4 yx w ∂∂ ∂ + 4 4 y w ∂ ∂ = D q (a) trong đó D - độ cứng tấm. Điều kiện biên: tại x = ± 2 a : w=0; 2 2 x w ∂ ∂ = 0; tại y = ± 2 a :w = 0; y w ∂ ∂ = 0; (b) Tấm được chia thành lưới 4x4 với Δx = Δy = a/4. Từ điều kiện biên (b) có thể thỏa mãn điều kiện w = 0 cho tất cả các nút nằm trên cạnh tấm x = ± a 2 và y = ± 2 a . Hệ phương trình đại số để xác định hàm chuyển vị w(x,y