Giáo trình Cơ học kết cấu - Chương 3, Phần 1: Phương pháp tính giải tấm

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) theo phương pháp sai phân hữu hạn được lập trên cơ sở biểu thức (a): Δx4 ( 4 4 x w ∂ ∂ + 2 22 4 yx w ∂∂ ∂ + 4 4 y 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 (ui,k+2 - 4ui,k+1 + 6uik - 4ui,k-1 + ui,k-2 ) = D xqik 4Δ (c’) trong đó: β = 2 2 y x Δ Δ ; 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ố thep 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 56 2w1 -8w2 - 8w3 + 22w4 + w5 + w9 = C (d’) trong đó C = D qa 256 4 (e’) Điều kiện biên 2 2 x w ∂ ∂ = 0 trên x = a/2 được thỏa mãn trong hệ phương trình: w2 + w6 = 0 w4 + w5 = 0 w7 = 0 (f’) Còn điều kiện y w ∂ ∂ = 0 đưa đến: w3 - w8 = 0 w4 - w9 =0 w10 = 0 (g’) Giải hệ phương trình đại số trên đây cho phép xác định hàm w(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: σx = ∂∂ 2 2y w = 2 ,1,1 )( 2 y www kiikki Δ +− −+ σy = ∂∂ 2 2x w = 2 1,1, )( 2 x www kiikki Δ +− −+ ; τxy = - ∂∂ ∂ 2 x y w = yx wwww kikikiki ΔΔ +−− −−+−−+++ 4 1,11,11,11,1 (h’) Mặt khác phương trình (a’) còn được được hiểu là: )( 22224 www ∇∇=∇∇=∇ (i’) Thay u = ∇2w vào phương trình Karman, có thể viết như sau: D quw =∇=∇∇ 222 . (j’) Phương trình ∇2u = C thường gặp trong lý thuyết trường và cơ học chất rắn, viết dạng đầy đủ sẽ là C y u x u =+ 2 2 2 2 ∂ ∂ ∂ ∂ . Với C ≡ 0 chúng ta nhận được phương trình Laplace. Với C = -2Gθ, trong đó G= E/2(1-ν), θ - góc xoắn, phương trình Poisson miêu tả momen xoắn trục tiết diện bất kỳ, dài hữu hạn. 57 Hệ phương trình suy từ phương pháp sai phân hữu hạn áp dụng cho phương trình Laplace và phương trình Poisson có dạng: 0 )( 2 )( 2 2 1,1, 2 ,1,1 =−Δ +−+Δ +− −+−+ C x uuu y uuu kiikkikiikki (k’) Như vậy chúng ta có thể giải phương trình Karman theo cách tùy chọn, tính trực tiếp hoặc thông qua phương trình Poisson khi dùng phương pháp sai phân. Phương pháp sai phân hữu hạn thuộc nhóm phương pháp số, có thể được chương trình hoá. Dưới đây giới thiệu một hàm viết bằng C, giải phương trình Poisson theo phương pháp sai phân hữu hạn, trích trong thư viện toán của người viết tài liệu này (1992). Giải phương trình Poisson ∇2Φ = Ψ(x,y) Gọi chương trình: Poison (int N, double length, double (*Phi)[MAX], double (*Psi)[MAX], int iters ); Với: N Số đoạn phải chia trong D. length Chiều dài cạnh L. Phi Matrận Φ . Psi Matrận Ψ . iters Số lần lặp theo phương pháp Gauss-Seidel MAX Kích cỡ. Phương trình cơ bản: Φ i-1,j - 2Φ i,j + Φi+1,j Φ i,j-1 - 2Φ i,j + Φ i,j+1 ------------------------------------ + ------------------------------------ = Ψ i,j (Δx)2 ( Δy)2 và Φ i,j ← (1/4) [ Φ i-1,j + Φ i+1,j + Φ i,j-1 + Φ i,j+1 - (Δx)2 Ψ i,j ] Phương trình vi phân miêu tả uốn tấm có dạng chung: ∇4w = ∂∂ 4 4 w x + 2 ∂∂ ∂ 4 2 2 w x y + ∂∂ 4 4 w y = q D (*) hàm w thỏa mãn điều kiện biên: w = 0 và ∂2w/∂η2 = 0 dọc các cạnh với η- pháp tuyến đến cạnh tấm, D = Et 3 12 1( )−ν E - mođun đàn hồi, t - chiều dầy tấm υ - hệ số Poisson. Nếu ký hiệu u = ∇2w, để giải bài toán (*), cần thực hiện hai lần tính như sau: ∇2u = q D với u = 0 dọc 4 cạnh, 58 ∇2w = u với w = 0 dọc 4 cạnh. /* to solve ellipstic equation #include #include #define IMAX 12 #define JMAX 12 int i,j,N, ITMAX; double length; double u[JMAX][IMAX], w[JMAX][IMAX], qdD[JMAX][IMAX]; ************************/ void Poisson( int N, double length, double ( *Phi)[IMAX], double ( *Psi)[IMAX], int ITMAX ) { register int iteration; register int i, j; double xsq, fn; fn = (double) N; xsq = (length/fn) * (length/fn); for (iteration =1; iteration <= ITMAX; iteration++) for (i=2; i <= N; i++) for ( j = 2; j <= N; j++) Phi[i][j] = 0.25*(Phi[i-1][j]+Phi[i+1][j] +Phi[i][j-1] + Phi[i][j+1] - xsq * Psi[i][j] ); return; } /**************************** main( ) { double q, t, Sigma, E, D, qD; int nplus1; printf("\n\tEnter N, ITMAX, q, t, Sigma, length, E:\n"); scanf("%d %d %lf %lf %lf %lf %lf", &N, &ITMAX,&q,&t, &Sigma, &length, &E); printf(" %d %d %9.5lf %9.5lf %9.5lf %9.5lf %9.5lf\n", N, ITMAX, q, t, Sigma, length, E); nplus1 = N +1; D = E*t*t*t / ( 12.0*(1.0 - Sigma*Sigma ) ); qD = q / D; for ( i =1; i <= nplus1; i++) for ( j =1; j <= nplus1; j++) { u[i][j] = w[i][j] = 0.0; qdD[i][j] =qD; } Poisson(N, length, u, qdD, ITMAX); /* solve d2(u)=q/D */ Poisson(N, length, w, u, ITMAX); /* solve d2(ww) = u */ printf("\n\t u are:\n"); for (i = 1; i <= nplus1; i++) { for (j=1; j <= nplus1; j++) printf("%9.5lf ", u[i][j]); printf("\n"); } printf("\n\t w are:\n"); for (i=1; i <= nplus1; i++) { for (j=1; j <= nplus1; j++) printf("%9.5lf ", w[i][j]) ; printf("\n"); } } ************************************/ 59 3.2 Phương pháp Rayleigh – Ritz. Phương pháp Rayleigh-Ritz, viết tắt R-R, là phương pháp tính nằm trong khuôn khổ phép biến phân. Tiền đề của phương pháp là phiếm hàm tương đương bài toán đã được xác định ở dạng I = F( u, u x x 1 2 ∫ x, ,x )dx, theo cách đặt vấn đề của Rayleigh. Phương pháp Rayleigh_Ritz tìm cách thay thế biến u trong phiếm hàm I bằng nghiệm gần đúng dưới dạng hàm xấp xỉ, theo cách làm của Ritz, ví dụ: u = (a) a fi i i N = ∑ 1 Hàm fi, i =1,2,..., N phải thoả mãn các điều kiện biên bài toán. Phiếm hàm I chứa u thường được tìm dưới dạng: I = F( u, u x x 1 2 ∫ x, ..., x )dx (b) trong đó ký hiệu ux = ∂u / ∂x. Các điều kiện biên có dạng: u(x1) = u(x2) = 0; (c) Hàm xấp xỉ u liên tục đến bậc N-1 và thỏa mãn các điều kiện biên. Đưa (a) vào (b), sau khi lấy đạo hàm theo ai cho phiếm hàm I, sẽ nhận được hệ phương trình đại số chứa các ẩn phải tìm ai: ∂I / ∂ai = 0, i=1,2,...,n Sau khi thay gía trị ai giải từ hệ phương trình cuối cùng, sẽ nhận được hàm xấp xỉ lập cho u, tính theo (a). Ví dụ:Tìm độ võng tấm chữ nhật và ứng suất 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: w(x,y) = u(x,y) = (x,y) (a’) a fn n n = ∞∑ 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’) 60 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 đủ cho trường hợp chung của tấm trực hướng: 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”) Trong công thức (c*) để tiện phân biệt ký hiệu người viết đã tạm thời sử dụng cách ghi sau: ν12 ≡ νxy; ν21 ≡ νyx ; E1 ≡ Ex; E2 ≡ Ey; G12 ≡ Gxy 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 = )1(12 2112 3 νν− tEx , D2 = )1(12 2112 3 νν− tEy và DT = G t12 3 12 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: 61 U = 1 2 t∫∫ D wx D wx wy D wy D wx y dxdyT1 2 2 1 21 2 2 2 2 2 2 2 2 2 2 4∂∂ ν ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ⎛ ⎝⎜ ⎞ ⎠⎟ + + ⎛ ⎝⎜⎜ ⎞ ⎠⎟⎟ + ⎛ ⎝⎜ ⎞ ⎠⎟ ⎡ ⎣ ⎢⎢ ⎤ ⎦ ⎥⎥ (e’) Trường hợp tấm đẳng hướng, Ex = Ey = E; ν12 = ν21 = ν; G = E2 1( )− ν Hàm thế năng có dạng: U = ( )D w wx y wx wy dxdy2 2 12 2 2 2 2 2 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: V = q(x,y).w(x,y)dxdy (h’) ∫∫ Á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 tải phân bố đều q = const. Tấm được ngàm cả bốn cạnh, chiều dầy tấm t = const. Tâm hệ tọa độ đặt tại góc trái phía trên của tấm. Điều kiện biên: Tại x = 0; x = a: 0;0 =∂ ∂= x ww Tại y = 0; y = b: 0;0 =∂ ∂= y ww Thế năng của tấm tính bằng quan hệ Π = U – V, tính như sau: ( ) dxdy y w x w yx wwD ∫∫ ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ −⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛−+∇ 2 2 2 222 22 )1(2 2 ∂ ∂ ∂ ∂ ∂∂ ∂ν - q(x,y).w(x,y)dxdy Thoả mãn điều kiện biên vế phải biểu thức, trình bày công biến dạng tấm trở thành: ∫∫ dxdy x w x wDU A 2 2 2 2 2 2 ∫∫ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ∂ ∂+∂ ∂= Hàm chuyển vị tấm trong trường hợp này được biểu diễn dạng: ∑∑∞ = ∞ = ⎟⎠ ⎞⎜⎝ ⎛ −⎟⎠ ⎞⎜⎝ ⎛ −= 1 1 2cos12cos1 m n mn b yn a xmaw ππ Hàm w thỏa mãn các điều kiện biên đang nêu. Đưa biểu thức tính w vào phương trình công biến dạng có thể thấy: 62 dxdy a xm b yn b n b yn a xm a maDU a b m n mn 2 2 2 0 0 1 1 2 2 2 2cos12cos 2cos12cos4 2 ⎭⎬ ⎫ ⎥⎦ ⎤⎟⎠ ⎞⎜⎝ ⎛ −+ ⎩⎨ ⎧ ⎢⎣ ⎡ ⎟⎠ ⎞⎜⎝ ⎛ −= ∫ ∫ ∑∑∞ = ∞ = ππ πππ Từ đây có thể khai triển hàm U: ⎪⎭ ⎪⎬ ⎫⎟⎠ ⎞⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛+ +⎪⎩ ⎪⎨ ⎧ ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ⎟⎠ ⎞⎜⎝ ⎛⎟⎠ ⎞⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛= ∑∑∑∑∑∑ ∑∑ ∞ = ∞ = ∞ = ∞ = ∞ = ∞ = ∞ = ∞ = snrn r s m msmr m r s mn m n aa b naa a m a b n a m b n a mabDU 1 1 1 4 1 1 1 4 2 1 1 2244 4 22 2332π với r ≠ s. Công do lực tập trung tác động: ∑∑∫ ∫ ∑∑ ∞ = ∞ = ∞ = ∞ = =⎥⎦ ⎤⎢⎣ ⎡ ⎟⎠ ⎞⎜⎝ ⎛ −⎟⎠ ⎞⎜⎝ ⎛ −= 1 10 0 1 1 2cos12cos1 m n mn a b m n mn aqabdxdyb yn a xmaqV ππ (d*) Trong công thức trên đây D = )1(12 2 3 ν− Et là độ cứng tấm, t - chiều dầy, E - mođun đàn hồi vật liệu, ν - hệ số Poisson. Tiến hành lấy đạo hàm của (U - V) theo amn: mna VU ∂ ∂ )( − sẽ nhận được hệ phương trình đại số, trong đó amn đóng vai trò ẩn số. 022 2334 1 4 1 4 2244 4 =−⎪⎭ ⎪⎬ ⎫⎟⎠ ⎞⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛+ ⎪⎩ ⎪⎨ ⎧ + ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ⎟⎠ ⎞⎜⎝ ⎛⎟⎠ ⎞⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛ ∑∑ ∞ = ∞ = qaba b na a m a b n a m b n a mabD rn r mr r mnπ với r ≠ n và r ≠ m. Trường hợp chỉ chọn m = n = 1, phương trình cuối cho phép rút ra: 242 4 11 )/(2)/(33 1 4 babaD qaa ++= π Với tấm hình vuông, a = b, công thức cuối mang dạng: a11 = qa4/32π4D. Chuyển vị lớn nhất, tại giữa tấm tính theo a11 vừa trình bày sẽ là: 63 D qaw 4 max 00128,0= Trường hợp tấm tựa bốn cạnh, biểu thức tính thế năng của tấm Π = U – V: ( ) dxdy y w x w yx wwD ∫∫ ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ −⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛−+∇ 2 2 2 222 22 )1(2 2 ∂ ∂ ∂ ∂ ∂∂ ∂ν - q(x,y).w(x,y)dxdy ∫∫ được triển khai như sau: ( ) dxdy b yn a xm b yn a xm b n a ma b yn a xm b n a maDU mn a b m n mn ⎭⎬ ⎫⎟⎠ ⎞⎜⎝ ⎛ −×−− ⎩⎨ ⎧ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ += ∫ ∫ ∑∑∞ = ∞ = ππππππν ππππ 2222 2 22 2 22 2 0 0 1 1 22 2 22 2 22 2 coscossinsin12 sinsin 2 và dxdy b yn a xmqaV a b mn∫ ∫= 0 0 sinsin ππ Sau tích phân biểu thức tính U và V sẽ nhận được: 2 1 1 2 2 2 2 2 4 8 ∑∑ ∞ = ∞ = ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ += m n mn b n a maDabU π và ( )( 1cos1cos 1 1 2 2 −−= ∑∑∞ = ∞ = πππ nmmn aqabV m n mn ) Thực hiện phép lấy đạo hàm riêng theo amn, ( ) 0=∂ −∂ mna VU có thể nhận hệ số a11 dạng sau đây: ( ),...3,1,16 222611 = ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ⎟⎠ ⎞⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛ = nm b n a mmn q D a π 3.3 Phương pháp Galerkin Phương pháp do Galerkin đề 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 64 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, (a) và các điều kiện biên: B (u) -q = 0 trên biên S (b) Để giải bài toán có thể tìm hàm u dưới dạng hàm xấp xỉ: u* = (c) 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 (a) vào bài toán (c), 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 ( a i ∑ ifi ) - p ≠ 0. 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,.... (d) Dưới dạng đầy đủ công thức cuối được viết lại: 0= ⎭⎬ ⎫ ⎩⎨ ⎧ −⎟⎠ ⎞⎜⎝ ⎛∫ ∑ dVfpfa i i iiL , i=1,2,.... (e) Ví dụ: Áp dụng công thức Galerkin giải bài toán uốn tấm chữ nhật 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 (phương trình Karman): 65 D.∇2 ∇2 w - 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 w, trong trường hợp này có thể là: w(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 = ∂w/ ∂x = 0; tại y= 0; y = b: u = ∂w/ ∂x = 0; (d’) Nếu chỉ chọn m = 1 và n =1, hàm w(x,y) sẽ có dạng: w(x,y) = a11. 1 4 (1 - cos 2πx a )(1 - cos 2πy b ) (e’) Sau khi thay thế giá trị w từ (e) vào phương trình Galerkin: ∫[D.∇2 ∇2 w - p(x,y)] δwdxdy = 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ị w tính tại tâm tấm là lớn nhất: wmax = 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 wmax = 0,00126 pa D 4 . 3.4 Phương pháp phần tử hữu hạn Trong bài toán phẳng phương trình chuyển vị u(x,y) và v(x,y) được viết dưới dạng vecto như sau: U= Pa 66 với U = và a = u x y v x y ( , ) ( , ) ⎧⎨⎩ ⎫⎬⎭ a a an 1 2 .... ⎧ ⎨ ⎪⎪ ⎩ ⎪⎪ ⎫ ⎬ ⎪⎪ ⎭ ⎪⎪ Hàm chuyển vị được tính theo cách quen thuộc:U = [N] {δ } với [N] - hàm hình dáng N(x,y), còn chuyển vị các nút {δ } = u v u v n n 1 1 ..... ⎧ ⎨ ⎪⎪⎪ ⎩ ⎪⎪⎪ ⎫ ⎬ ⎪⎪⎪ ⎭ ⎪⎪⎪ Trường hợp biến dạng phẳng, vector biến dạng được biết là: {ε} = = ε ε γ x y xy ⎧ ⎨⎪ ⎩⎪ ⎫ ⎬⎪ ⎭⎪ ∂ ∂∂ ∂ ∂ ∂ ∂ ∂ u x v y u y v x + ⎧ ⎨ ⎪⎪⎪ ⎩ ⎪⎪⎪ ⎫ ⎬ ⎪⎪⎪ ⎭ ⎪⎪⎪ = ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ x y y x 0 0 ⎡ ⎣ ⎢⎢⎢⎢⎢⎢⎢ ⎤ ⎦ ⎥⎥⎥⎥⎥⎥⎥ = [B]{δ} (a) u v ⎧⎨⎩ ⎫⎬⎭ Từ đó B = ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ x N y N y N x N [ ] [ ] [ ] [ ] 0 0 ⎡ ⎣ ⎢⎢⎢⎢⎢⎢⎢ ⎤ ⎦ ⎥⎥⎥⎥⎥⎥⎥ (b) Và vector ứng suất phẳng: {σ} = = [D] (c) σ σ τ x y xy ⎧ ⎨⎪ ⎩⎪ ⎫ ⎬⎪ ⎭⎪ ε ε γ ε x y xy ⎧ ⎨⎪ ⎩⎪ ⎫ ⎬⎪ ⎭⎪ − ⎛ ⎝ ⎜⎜⎜ ⎞ ⎠ ⎟⎟⎟{ }0 trong đó matrận cứng [D] được xác định cho mỗi loại vật liệu. Phiếm hàm thế năng có dạng:Π = 1 2 { } { }ε σT A dA∫ − { } { }p u dsT S ∫ (d) hoặc khai triển:Π = 1 2 dAxyxyyyx A x )( γτεσεσ ++∫ - (e) Ma trận cứng phần tử tính theo biểu thức: dsvpup y S x )( +∫ [k] = [B]∫∫ A T [D] [B] dxdy. (f) như đã trình bày tại tài liệu về phương pháp PTHH, cùng người viết. Cách tính trên đây áp dụng cho một phần tử tấm sẽ đưa lại kết qủa sau. 67 Phần tử hình chữ nhật cạnh axb. Hàm chuyển vị nút được trình bày dưới dạng: u(x,y) = P a (a’) trong đó uT = [ u1 u2 u3 u4 v1 v2 v3 v4], (b’) aT = [ a1 a2 a3 a4 b1 b2 b3 b4], (c’) P = ; với P = [ 1 x y xy] (d’) P P 0 0 ⎡ ⎣⎢ ⎤ ⎦⎥ Nếu coi U = Qa, trong đó Q = và α α 0 0 ⎡ ⎣⎢ ⎤ ⎦⎥ α = 1 1 1 1 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 x y x y x y x y x y x y x y x y ⎡ ⎣ ⎢⎢⎢⎢ ⎤ ⎦ ⎥⎥⎥⎥ a = Q-1 U (e’) Từ đó: U(x,y) = N(x,y)δ, (f’) trong đó: N(x,y) = PQ-1 (g’) Sử dụng nguyên lý biến phân năng lượng như đã thực hiện cho các phần đã nêu, có thể xác định được matrận cứng tấm hình chữ nhật. Phương trình tính [k] đã nêu tại (f) phần trên, [k] = [B] V ∫ T [D] [B] dV (h’) Giá trị kij của matrận 8x8 có dạng sau đây. Trong ma trận các ký hiệu mang ý nghĩa qui ước: k0 = 1 3 [(b/a) + ½ (1 - ν ) (a/b)] k00 = 1 3 [(a/b) + ½ (1 - ν ) (b/a)] k1 = - 1 6 [2(b/a) - ½ (1 - ν ) (a/b)] k11 = - 1 6 [2(a/b) - ½ (1 - ν ) (b/a)] 68 k2 = - 1 6 [(b/a) - ½ (1 - ν ) (a/b)] k22 = - 1 6 [(a/b) - ½ (1 - ν ) (b/a)] k3 = + 1 6 [(b/a) - (1 - ν ) (a/b)] k33 = + 1 6 [(a/b) - (1 - ν ) (b/a)] k4 = 1 8 (1+ν) k5 = 1 8 (3ν -1) ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ − − − − −− −− −= 00 40 33500 4140 22411500 425340 11522433500 53425140 21 ][ k kk kkkDX kkkk kkkkk kkkkkk kkkkkkk kkkkkkkk Etk ν Uốn tấm. Dưới tác động của lực và moment (tập trung hoặc phân bố) tác động theo phương pháp tuyến tấm bị uốn. ∫∫∫ −−− === 2/ 2/ 2/ 2/ 2/ 2/ ;;; t t xyxy t t yy t t xx dzNdzNdzN τσσ 69 t/ 2 Ứng suất t/ 2 Momen và lựcσ σ τ τ τ τ τ σ τ τ τ σ z y x xy xO q yq M M M M yx y x xyxy xy xy xy y yx yx xz xz x x y Hình 3.2 Ký hiệu lực, momen của tấm bị uốn. Ứng suất trong phần tử tấm gồm: ứng suất cắt τyz, τxz, τxy và ứng suất pháp tại các mặt cắt ngang σx, σy. Moment và lực được định nghĩa như sau: Mx = ; Mzdz t t x∫+− 2/2/ σ y = ; Mzdztt y∫+− 2/2/ σ xy = zdztt xy∫+− 2/2/ τ qx = τ− +∫ tt//22 xz dz; qy = τ−+∫ tt//22 yz dz Trong tài liệu còn sử dụng ký hiệu Mx ≡ Mxx và My ≡ Myy. Quan hệ giữa chúng thể hiện trong hệ phương trình cân bằng lực. 0=++ p y q x q yx ∂ ∂ ∂ ∂ x xyx q y M x M =+ ∂ ∂ ∂ ∂ y yxy q y M x M =+ ∂ ∂ ∂ ∂ trong đó p - lực phân bố, tác động theo pương pháp tuyến với mặt tấm. Biến dạng tấm thể hiện bằng công thức: 70 { } ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ −= ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = yx w y w x w z xy y x ∂∂ ∂ ∂ ∂ ∂ ∂ γ ε ε ε 2 2 2 2 2 .2 (a) Nếu biểu diễn chuyển vị tấm theo cách quen thuộc: w = [N]{δ} (b) còn {δi } = ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ ⎟⎠ ⎞⎜⎝ ⎛ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛−= ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ i i i yi xi i x w y w ww ∂ ∂ ∂ ∂ θ θ (c) Vec tor biến dạng: {ε} = = [B]{δ} (d) ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ xy y x γ ε ε trong đó: [Bi] = − − ⎧ ⎨ ⎪⎪⎪ ⎩ ⎪⎪⎪ ⎫ ⎬ ⎪⎪⎪ ⎭ ⎪⎪⎪ ∂ ∂ ∂ ∂ ∂ ∂ ∂ 2 2 2 2 22 [ ] [ ] [ ] Ni x Ni y Ni x y (e) Vector ứng suất tấm được xét trong bài toán này: {σ} ≡ {M} = [D]( {ε} - {ε0 } ) (f) Matrận [D] được ghi tại phần đầu tài liệu. Với vật liệu đẳng hướng, [D] trong bài toán phẳng có dạng sau. [D] = ( )Et 3 212 1− ν 1 0 1 0 0 0 1 2 ν ν ν( ) /− ⎡ ⎣ ⎢⎢⎢ ⎤ ⎦ ⎥⎥⎥ Phương trình chuyển vị của tấm dưới tác động lực phân bố p(x,y) có dạng: w yyxx ).2( 4 4 22 4 4 4 ∂ ∂ ∂∂ ∂ ∂ ∂ ++ = D p 71 trong đó D = )1(12 2 3 ν− Et Ứng suất trong phần tử tấm, tính tại các mặt cắt có dạng: ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪ ⎬ ⎫ = = = 12/ 12/ 12/ 3 3 3 t zM t zM t zM xy xy y y x x τ σ σ Ứng suất chính tính bằng công thức: σ1 = 12 (σx + σy) + 1 2 22 4)( xyxy τσσ ++ σ2 = 12 (σx + σy) - 1 2 22 4)( xyxy τσσ +− Phiếm hàm áp dụng cho tấm bị uốn có dạng chung: Π = 1 2 V ∫ {ε}T {σ}dV - p(x,y)wdV (g) V ∫ Từ đó matrận cứng phần tử tấm được tính theo công thức quen thuộc: [k] = [B] A ∫∫ T [D] [B] dxdy (h) Vector lực tính theo: {F} = - [N] A ∫∫ T p(x,y)dxdy (i) Tấm mỏng hình chữ nhật, cạnh axb, 12 bậc tự do. Phương trình chuyển vị: w(x,y) = [P]{a} = [ 1 x y x2 y2 xy x2y xy2 x3 y3 x3y xy3 ] {a} (a) Phương trình trên đây thỏa mãn điều kiện liên tục của bài toán uốn tấm: ∇2 ∇2 w = 0 (b) Áp dụng phương trình trên cho chuyển vị các nút sẽ nhận được: {δ} = [C]{a} (c) Matrận [C] có kích thước 12x12, phụ thuộc vào toạ độ các nút i, i= 1, 2, 3, 4. {a} = [C]-1 {δ} (d) Từ đó: 72 w = [N] {δ} = [P] [C]-1 {δ} (e) [N] = [P] [C]-1 (f) Biểu thức Ni, i=1, 2, 3, 4 sau khai triển: [Ni] = 1 2 [ ( ξ0 + 1)(η0 + 1) (2 + ξ0 + η0 - ξ2 - η2), aξi(ξ0 +1)2(ξ0 -1)( η0 +1), bηi(ξ0 +1)( η0 +1)2(η0 -1) ], (g) trong đó: ξ = x x a c− ; η = y y b c− ; ξ0 = ξξi ; η0 = ηηi Thay các giá trị trên vào biểu thức tính {ε} sẽ nhận được: {ε} = [Q]{a} = [Q]{C]-1 {δ} (h) và từ đó: [B] = [Q]{C]-1 (i) [Q] = 0 0 0 2 0 0 6 2 0 0 6 0 0 0 0 0 0 2 0 0 2 6 0 6 0 0 0 0 2 0 0 4 4 0 6 62 2 − − − − − − − ⎡ ⎣ ⎢⎢⎢ ⎤ ⎦ ⎥⎥⎥ x y xy x y xy x y x y − Nếu đưa biểu thức (d) vào vị trí [a] tại (h) sẽ nhận được biểu thức tính [B] và tiếp đó tính [k] theo công thức chung [B] A ∫∫ T [D] [B] dxdy: [k] = {[C]-1 }T ( ∫∫ [Q]T [D] [Q] dxdy ) [C]-1 . (j) Bằng cách tương tự biểu thức tính vector lực giờ có thể mang dạng: {F} = { -[C]-1 }T ∫∫ [P]T pdxdy. (k) Công thức chuẩn bị sẵn tính [k] theo (j) tham khảo trong các sách bàn về phương pháp phần tử hữ hạn. ( ) [ ] [ ] [ ] [ ] ][21][1180][ 13212 3 RkkkkR ba Etk ⎭⎬ ⎫ ⎩⎨ ⎧ −+++−×= ννν (l) trong đó: a và b – chiều dài hai cạnh phần tử tấm chữ nhật. Nếu ký hiệu a