Giáo trình Phương pháp phần tử hữu hạn - Chương 1: Phương pháp phần tử hữu hạn (Phần 1)

CHƯƠNG 1 PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN Trong phần bàn về cơ học kết cấu, hai nhóm phương pháp số xây dựng bài toán cơ học kết cấu đã được đề cập: phương pháp thành lập phương trình vi phân và phương pháp biến phân. Nhóm thứ nhất bao gồm ba phương pháp chúng ta đã sử dụng trong phần I giáo trình, phương pháp chuyển vị, phương pháp lực và phương pháp hỗn hợp chuyển vị – lực; nhóm thứ hai thuộc về các phương pháp năng lượng, nền tảng của cơ học kết cấu, gồm các nguyên lý: - Nguyên lý tổng thế năng tối thiểu. - Nguyên lý công bù nhỏ nhất. - Nguyên lý công dừng do Reissner đề xuất giành cho hàm năng lượng. Hai phương pháp xây dựng bài toán cơ học kết cấu 1/ Phương pháp phương trình vi phân 2. Phương pháp biến phân Phương pháp chuyển vị Nguyên lý thế năng tối thiểu Phương pháp lực Nguyên lý công bù nhỏ nhất Phương pháp chuyển vị-lực Nguyên lý Reissner Trong thực tế tính toán có thể thấy rõ rằng, phương pháp chuyển vị trong nhóm phương trình vi phân liên hệ gắn bó đến mức không tách rời nguyên lý thế năng tối thiểu, phương pháp lực có chung gốc nguyên lý công bù nhỏ nhất, trong khi đó phương pháp hỗn hợp chuyển vị – gần gũi với nguyên lý công dừng của Reissner. Các phương trình vi phân xây dựng theo phương pháp lực, phương pháp chuyển vị và phương pháp hỗn hợp trình bày tóm tắt như sau: Từ phương trình cân bằng: có thể viết dạng chung: D ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ + ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ ∂∂∂ ∂∂∂ ∂∂∂ 0 0 0 000 000 000 Bz By Bx yz zx xy z y x yzz zxy xyx f f f τ τ τ σ σ σ L M M M T σ + fB = 0 . Định luật Hooke trình bày theo quan hệ: σ = Cε hoặc ngược lại ε = E.σ, hay là ε = C-1σ. Áp dụng nguyên lý công ảo giải phương trình cân bằng DTσ + pv = 0 trong miền V, cùng các điều kiện biên: p = p* hay là p - p* = 0 tại biên Sp, có thể viết: ( ) ( )∫ ∫ =−−+V Sp *TVTT dSppudVpσDu 0δδ 2 Biên S có thể coi là tập họp từ hai biên: S = Sp + Su, tại Su điều kiện biên thuộc chuyển vị u được áp đặt. Tích phân thứ hai phương trình cuối được viết thành: ∫ ∫∫ −= S S TTS T up dSupdSupdSup δδδ Nguyên lý thế năng tối thiểu Thay thế quan hệ ε = Du và δε = Dδu tại V, nêu tại phần đầu tài liệu, vào phương trình trên đây, cùng điều kiện biên δu = 0 tại Su có thể viết: ∫ ∫ ∫ =−− V V Sp T V TT pdSudVpuσdVε 0* δδδ Phương trình dạng này về mặt hình thức chứa ứng suất và chuyển vị, mặc dầu chuyển vị là ẩn cơ bản. Để chuyển phương trình về hẵn dạng chỉ chứa chuyển vị cần thiết tiến hành các bước biến đổi như chúng ta đã thực hiện trong phần đầu sách. Với vật liệu đàn hồi, σ = Eε, có thể sử dụng quan hệ động học ε = Du khi xác định ứng suất: σ = Eε = E(Du) = EDu.u Vì rằng εT = (Du)T, chúng ta có thể viết: δεT = δ(Du)T = δuTDuT. Thay các biểu thức vừa xác định vào phương trình sau: (a) ∫ ∫ ∫ =−− V V Sp T V TT pdSudVpuσdVε 0* δδδ có thể viết tiếp: ( )[ ]∫ =−− V S Vu T u T dSpdVpuEDDu 0*δ ∫ p (b) Trong đó ma trận E và Du được định nghĩa như sau: D ; E ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ ∂∂ ∂∂ ∂∂ ∂ ∂ ∂ = yz xz xy z y x u 0 0 0 00 00 00 LLL ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ − − − −= − − − 2 21 2 21 2 21 0 1 01 1 21 2 ν ν ν ννν ννν ννν ν M M M LLLLLL M M M G ( )ν+= 12 EG , với E – mô đun đàn hồi. Nguyên lý công bù nhỏ nhất Từ quan hệ ε = Du trong miền V, và các điều kiện biên u = u* tại Su, trong khuôn khổ nguyên lý công bù ảo có thể viết: ( ) ( ) 0* =−−− ∫∫ Su T V T dSuupdVDuεσ δδ (c) hay là 0* =− ∫∫ Su T V T dSupεdVσ δδ (d) 3 Phiếm hàm Hellinger-Reisner Từ phương trình năng lượng áp dụng cho vật thể đàn hồi có thể viết: ( )∫ ∫ ∫ ∫∫ −++=+− V V Sp Su TT V B TT dSuupdSpudVpudVdVU ***0 )( εσσ Trong đó U0*(σ) = σTε - U0(ε), xem hình dưới. σ σ = d== σ u* u Hình 1.2 Từ phương trình cuối có thể viết: =)(*0 σU σTε - 2 1 σTε = 2 1 σTEσ, hay là ∫ ∫ −= V V T σdVEσdVU 1*0 2 1)(σ Hàm năng lượng trên đây còn được biết dưới dạng thường gặp: ( ) ∫ ∫ ∫∫∫ −−−−+−=Π V Sp Su iiiiiiBi V jiij V R dSuupdSpudVpudVudVU )( ** ,, * 0 σσ (e) Hàm ΠR mang tên gọi phiếm hàm Hellinger-Reisner, là phương trình chủ đạo của phương pháp hỗn hợp (mixed). Thoả mãn điều kiện dừng, biến phân δΠR tính theo cách sau: ( ) ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ −−−−+− ∫ ∫ ∫∫∫ V Sp Su iiiiiiBi V jiij V dSuupdSpudVpudVudVU )( **,, * 0 σσδ phải bằng 0. Nếu đưa biểu thức vào công thức cuối, đồng thời để ý đến quan hệ U ⎥⎦ ⎤⎢⎣ ⎡ − ∫∫ V ijij V ijij dVdV σεσεδ 0 *(σ) = σTε - U0(ε) có thể thấy: ( ) ∫∫ ∫ =+− VV V ijij dVUdVdVU )(0 * 0 εεσσ . Biến phân hàm năng lượng dạng pha trộn: ( ) ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ −++=+∫ ∫ ∫ ∫∫ V V Sp Su *T*T V B TT* dSuupdSpudVpuεdVσ(σσ)dU 0δ còn được hiểu theo cách khác như sau: ( ) ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ −−−−−−∫ ∫∫∫ ∫ V Su *T V T V Sp *T B T* )dSu(updVDuεσdSpudVpu(εε)dU 0δ 4 Ba thành phần trong dấu ngoặc lớn biểu thị nguyên lý công ảo, còn hai thành phần còn lại ∫ ∫∫ −−− V Sp *T B T V dSpudVpudVU )(*0 ε ( ) ∫∫ −−− Su *T V T )dSu(updVDuεσ diễn tả điều kiện động học của vấn đề. Từ điều kiện dừng, biến phân đang đề cập bằng 0, chúng ta có thể xác định phiếm hàm của thế năng: ( ) ( )∫ ∫∫∫ ∫ −−−−−−=Π V Su iii V ijji V Sp iiiiB dSuppdVudSupdVupdVU * , * , * 0 )( εε (f) Giả sử rằng chuyển vị u thoả mãn điều kiện động học trong miền V, hay là ε = Du trường hợp này biến phân vừa xác lập mang dạng: = 0 ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ −−−−∫ ∫∫ ∫ V S T V Sp T B T dSuupdSpudVpudVU )()( ***0 εδ u u Thành phần là thành phần bổ sung cho điều kiện biên δp ≠ 0 tại S ∫ − Su T dSuup )( * p. Phiếm hàm (g) trong trường hợp này mang tên gọi phiếm hàm lai trộn (hybrid functional). ∫ ∫∫ ∫ −−−−=Π V S T V Sp T B T H dSuupdSpudVpudVU )()( *** 0 ε Phương pháp phần tử hữu hạn trong nhóm các phương pháp số, sử dụng cả ba cách làm đang nêu khi xây dựng bài toán. Những bài toán thường gặp khi xác định ứng suất, chuyển vị trong cơ học kết cấu thông thường được xử lý theo phương pháp chuyển vị. Trong giáo trình này chủ yếu đề cập phương pháp chuyển vị, hay là phương pháp tính dựa trên nguyên lý thế năng tối thiểu. 1. Phương pháp biến phân Để giải phương trình cân bằng dạng tổng quát: ‹(u) = p trong miền V (a) và thoả mãn các điều kiện biên: Bk(u) = qk trên biên S , k =1,2,... (b) Bước đầu tiên phải tìm phiếm hàm tương ứng bài toán (a) trong miền V, và phiếm hàm đó phải thỏa mãn điều kiện biên (b). Phiếm hàm, ký hiệu I(u), thông thường là tích phân chứa u và các đạo hàm riêng của u trong toàn miền V và biên S. Nếu ký hiệu toán hạng của phiếm hàm là F, công thức chung của phiếm hàm sẽ là: I = (c) ∫ V x dxxuuF ,...),,( Trong đó F - phương trình Euler-Lagrange Biến hàm u được tìm dưới dạng: 5 ∑ = = N i ii fau 1 * (d) trong đó các hàm fi - là hàm tuyến tính, hay còn gọi hàm thử; ai - các hệ số phải tìm từ quan hệ sau: 0*)( = ia u ∂ ∂I ; i=1,2,... (e) Các bước thực hiện khi sử dụng phương pháp biến phân như sau: (1) Áp dụng nguyên lý Euler-Lagrange để tìm phiếm hàm tương ứng cho bài toán dạng (a). Phiếm hàm I(u) xác định trong miền V, thỏa mãn điều kiện biên trên biên S. I = (f) ∫∫ + 1 ,...),,( SV qdSdxxuuF (2) Miền V được chia thành E phần tử, liên kết với nhau qua các nút, biên. (3) Biến hàm u, nghiệm cần tìm của bài toán đưa về dạng hàm chuyển vị các nút phần tử: {u }= [N]{δ} (g) trong đó: [N] = { N1 N2....]; [δ ] = [ u1 u2..... ] (4) Thực hiện phép đạo hàm: ∑ = = E e i e k uu 1 ∂ ∂ ∂ ∂ II = 0, i =1,2,..., N (h) Trong các bài toán cơ học kết cấu phương trình trên tương đương biểu thức: }{δ∂ ∂I = [k]{δ} - {p}, tính trong mỗi phần tử. (i) Ma trận [k] và {p} được tính cho mỗi phần tử, mang ký hiệu [k]e, {p}e. (5) Tập họp ma trận cứng và vector lực: [K]{u} = {P} (j) trong đó: [ ] và { } ∑ = = E e ekK 1 ][ ∑ = = E e epP 1 }{ Ma trận [K] và vector {P} xác định trong hệ toạ độ chung, do vậy trước khi tập họp cần tiến hành tính chuyển [k]∑∑ == E e e E e e pk 11 }{][ e sang hệ tọa độ chung. Ví dụ minh họa có thể xem ở trường hợp giản đơn nhất là phần tử BAR sau đây. 6 Hình 1.3 Chuyển vị trong hệ tọa độ cục bộ, với ký hiệu * trên ký tự, gắn liền với phần tử, từ hình vẽ có thể thấy: ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − −= ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ 2 2 1 1 * 2 * 2 * 1 * 1 cossin00 sincos00 00cossin 00sincos v v v v u u u u αα αα αα αα hoặc: {u*} = [T]. {u}, trong đó {u* } tính trong hệ tọa độ cục bộ. Ma trận chuyển trường hợp này được định nghĩa là: [ ] ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − −= αα αα αα αα cossin00 sincos00 00cossin 00sincos T Ma trận [k]e và vector {p}e tính trong hệ tọa độ chung có dạng như sau: [k]e = [T]T[k]e* [T], trong đó [k]* tính trong hệ tọa độ riêng. {p}e = [T]T.{p}e*. Cách tính ma trận [T} cho các kiểu phần tử được trình bày tại phần ma trận chuyển. Vector {u} sau khi giải hệ phương trình đại số, được tính trong hệ tọa độ chung. (6) Xác định vector {u} và lấy đạo riêng của u nếu cần. Ứng suất, biến dạng trong phần tử tính trong hệ tọa độ gắn liền phần tử, do vậy trước khi tính cần tính chuyển vector chuyển vị nút của phần tử đang miêu tả trong hệ tọa độ chung sang hệ tọa độ riêng. {u} = [T]-1 {u*} Phương pháp biến phân áp dụng cho cơ học kết cấu có ý nghĩa đặc biệt. Phiếm hàm cần tìm cho hệ kết cấu như hệ dầm, tấm, vỏ, kết cấu 3D vv... chính là hàm thế năng của hệ. Điều này sẽ được làm rõ hơn trong các phần tiếp theo. Ví dụ đơn giản nhất khi ứng dụng nguyên lý biến phân xác lập bài toán cơ học vật rắn có thể là mô hình xác định lực kéo, nén thông qua chuyển vị của lò xo trụ 7 sau đây. Dưới tác động của lực P lò xo bị kéo dài ra thêm u, thoả mãn điều kiện ku = P. Công biến dạng được tính bằng diện tích hình tam giác như đã dẫn giải trong chương lý thuyết đàn hồi, xem hình 1.4, bằng ½ Pu = ½ k.u2. Công mà ngoại lực P thực hiện trên quảng đường u là W = Pu. Từ đó PukuWU −=−=Π 221 . Biến phân của hàm Π được tính δΠ = (ku-P) δu. Từ đó có thể thấy, hàm Π đạt cực trị, cụ thể tại đây là đạt minimum khi δΠ = 0. Trường hợp đang xem xét biến phân δu ≠ 0, chúng ta có thể viết: ku - P = 0 (*) Hình 1.4 Sử dụng công thức trên cho hệ kết cấu, theo đó thay vì chuyển vị u sẽ thay thế bằng vector {u}, lực P thay bằng vector lực {R} và độ cứng k thay bằng ma trận [K], sẽ nhận được công thức tổng quát: U = 2 1 {u}T[K] {u} và W = {u}T{R} Π = U - W = 2 1 {u}T[K] {u} - {u}T{R} Áp dụng nguyên lý biến phân vào trường hợp này sẽ được: [K]{u} = {R} (**) Trong giáo trình sử dụng công thức (a) nêu tại nguyên lý thế năng tối thiểu: cho các bước tính tóan tiếp theo. ∫ ∫ ∫ =−− V V Sp T V TT pdSudVpuσdVε 0* δδδ Trường hợp chung của vật thể 3D từ vật liệu đàn hồi, các đại lượng đặc trưng chuyển vị, ứng suất trình bày tại lý thuyết đàn hồi được biểu diễn lại như sau đây. Hàm chuyển vị: { } [ ]{ }δN w u u = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = v (a) Với vật thể 3D, matrận [N] có thể viết dưới dạng: [ ] ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ = NlNkNjNi NlNkNjNi NlNkNjNi N 00000000 00000000 00000000 (b) vector ngoại lực: 8 { } ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = z y x p p p P (c) Vector :{ } }]{[v 0 0 0 00 00 00 v v v δ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂∂ ∂ γ γ γ ε ε ε ε B w u xz yz xy z y x z u x w y w z xy u z w y x u zx yz xy z y x = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ = ⎪⎪ ⎪⎪ ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪⎪ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪⎪ ⎪⎪ ⎪ ⎨ ⎧ + + + = ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ = (d) Matrận [B] được thể hiện cụ thể cho mỗi trường hợp i, i=1,2,...: [ ] ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ = ii ii ii i i i i N x N z N y N z N x N y N z N y N x B ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂∂ ∂ 0 0 0 00 00 00 (e) Vector ứng suất đươc viết dạng đầy đủ: { } = [D][B]{δ}-[D]{ε[ ] { } { }( 0εε τ τ τ σ σ σ σ −= ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ = D zx yz xy z y x ) 0} (f) Năng lượng biến dạng: ∫∫∫∫∫∫ == V T V T dVDdVU }]{[}{ 2 1}{}{ 2 1 εεσε (g) Dưới dạng đầy đủ, biểu thức cuối được triển khai thành: 9 ∫∫∫∫∫∫ −= V T V T dVDdVDU }]{[}{}]{[}{ 2 1 0εεεε (g’) Công ngoại lực tính theo công thức chung: 1}{}{}{}{ 1 dSUPdVUPW T V S S T∫∫∫ ∫∫+= (h) Trong đó {P}T – vector lực thể tích, {PS}T – vecto lực bề mặt. Hàm thế năng có dạng: [ ] ∫∫∫∫∫∫∫∫ −−−=Π 1 10 }{}{}{}{2}{][}{2 1 S T S V T V T dSPdVUPdVD εεε (i) hoặc dưới dạng thường dùng sau: ∫∫∫∫∫ ∫∫∫∫∫∫ − −−=Π 1 1 0 }{][}{}{][}{ }]]{[][}{}]{][[][}{ 2 1 S S TT V TT V TT V TT dSPNdVPN dVDBdVBDB δδ εδδδ (i’) 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 toàn hệ đúng bằng tổng thế năng tất cả phần tử họp thành: ∑ = =Π E e e 1 π (j) Thế năng phần tử tính theo công thức dạng (i) song áp dụng chỉ cho một phần tử [ ] ∫∫∫∫∫∫∫∫ −−−= eS T S Ve e T Ve e T e dSPdVUPdVD e 1 10 }{}{}{}{2}{][}{2 1 εεεπ . Công thức tính hàm Π cho toàn hệ thống (j) có thể viết thành: }{}{}{][}{][ }]]{[][}{}{]][[][}{ 2 1 1 1 0 11 C T eS S T Ve T Ve T E e T Ve T E e T PdSPNdVPN dVDBdVBDB Δ−⎟⎟⎠ ⎞+ ⎜⎜⎝ ⎛ +Δ−Δ⎥⎦ ⎤⎢⎣ ⎡Δ=Π ∫∫∫∫∫ ∫∫∫∑∫∫∫∑ == ε (k) trong đó {PC} – lực tập trung. Sau đạo hàm Π theo {Δ} sẽ nhận được hệ phương trình: 10 }{}{][}{][ }]]{[][}{]][[][ 1 1 0 11 C eS S T Ve T Ve T E eVe T E e PdSPNdVPN dVDBdVBDB −⎟⎟⎠ ⎞+ ⎜⎜⎝ ⎛ +=Δ⎥⎦ ⎤⎢⎣ ⎡ ∫∫∫∫∫ ∫∫∫∑∫∫∫∑ == ε (l) Trong các công thức trên, mỗi biểu thức có tên gọi riêng như sau. ∫∫∫= Ve T e dVBDBk ]][[][][ , ma trận cứng phần tử, ∫∫∫= Ve Ti e dVDBP }]]{[][}{ 0 )( ε , vector lực do biến dạng ban đầu, ∫∫= eS S TS e dSPNP 1 1 )( }{][}{ , vector lực mặt, ∫∫∫= Ve TB e dVPNP }{][}{ )( , vector lực thể tích. Dưới dạng thường gặp phương trình cân bằng toàn hệ được viết: [K]{Δ} = {P} (m) với: e E e kK ][][ 1 ∑ = = ∑∑∑ === +++= E e B e E e S e E e i eC PPPPP 1 )( 1 )( 1 )( }{}{}{}{}{ {Δ} – vector chuyển vị các bậc tự do hệ thống, xác định trong hệ tọa độ chung. Ví dụ 1: Bằng phương pháp biến phân giải phương trình vi phân, ghi tại (a) sau đây, theo cách làm trong khuôn khổ phương pháp PTHH. 02 2 =++ xu dx ud trong miền 0 ≤ x ≤ 1 (a) Điều kiện biên: u(0) = u(1) =0; (b) Lời giải: Phiếm hàm tương ứng phương trình (a) có dạng: dxuxu dx duI ∫ ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ++⎟⎠ ⎞⎜⎝ ⎛−= 1 0 2 2 2 2 1 Tiến hành chia đoạn 0 -1 thành E phần tử, chiều dài phần tử L = 1/E. 11 Hàm u trong mỗi phần tử được tìm dưới dạng quen thuộc từ các ví dụ trước. Thay u = [N]{δ} vào biểu thức của phiếm hàm, tiến hành lấy tích phân của phiếm hàm trong phạm vi phần tử: ∫ =⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ++⎟⎠ ⎞⎜⎝ ⎛−= j i x x e dxuxudx duI 22 2 ∫ ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ++⎥⎦ ⎤⎢⎣ ⎡⎥⎦ ⎤⎢⎣ ⎡− j i x x TT T T dxNNN dx dN dx dN }]{[2}]{[][}{}{}{ 2 1 δδδδδ Còn phiếm hàm của toàn hệ bằng tổng các phiếm hàm phần tử: ∑ = = E e eII 1 Sau khi lấy đạo hàm riêng theo ui: 0 1 == ∑ = E e i e i u I u I ∂ ∂ ∂ ∂ , i=1,2,... hoặc là: 0][}]{[][}{ 1 = ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ++⎥⎦ ⎤⎢⎣ ⎡⎥⎦ ⎤⎢⎣ ⎡−∑∫ = E e x x TT T j i dxNxNN dx dN dx dN δδ sẽ nhận được: [ ]{ } { }∑∑ == = E e ee E e e Pk 11 δ [k]e tính từ hai biểu thức đầu tiên của phương trình, có dạng: [ ] = ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ ⎥⎦ ⎤⎢⎣ ⎡ −−− ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ − −− −⎥⎦ ⎤⎢⎣ ⎡− ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧− = ∫ dxl xxl xx l xx l xx ll l lk j i x x ij i j e 11 1 1 )( = 1 1 1 1 1 6 2 1 1 2l l− − ⎡ ⎣⎢ ⎤ ⎦⎥ − ⎡ ⎣⎢ ⎤ ⎦⎥ và { } ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ −− −+= ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ − − = ∫ )2( )2(61 22 22 )( ijij ijij x x i j e xxxx xxxx dx l xx l xx xp j i Nếu chọn E = 2, chiều dài phần tử l = 0,5; Với e = 1: xi = 0,0 và xj = 0,5. Với e = 2: xi = 0,5 và xj = 1,0. Và như vậy: [ ] [ ] ⎥⎦ ⎤⎢⎣ ⎡ − −=⎥⎦ ⎤⎢⎣ ⎡−⎥⎦ ⎤⎢⎣ ⎡ − −== 2225 2522 12 1 21 12 6 5,0 11 11 5,0 1 )2()1( kk 12 { } { } ⎭⎬ ⎫ ⎩⎨ ⎧= ⎭⎬ ⎫ ⎩⎨ ⎧= 5 4 24 1; 2 1 24 1 )2()1( pp Tập họp ma trận cứng cho toàn hệ tiến hành theo cách sau: [ ] ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ − −− − = ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ − −+− − =∑ = 22250 25)44(25 02522 12 1 22250 25)2222(25 02522 12 1 1 E e ek { } { } ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ == ∑ 3 2 1 u u u E eδδ { } ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ += 5 )6( 1 5 )42( 1 24 1P Điều kiện biên bài toán: u1 = u3 = 0. Phương trình cân bằng giờ có dạng: )6( 24 1)44( 12 1 2 =u Từ đó: u2 = 3/44 = 0,06817. Ví dụ 2: Giải phương trình xoắn dầm, mặt cắt ngang hình dạng bất kỳ. Bài toán cơ bản xoắn dầm luôn được đề cập trong các tài liệu lý thuyết đàn hồi và cơ học kết cấu. Phương trình Poisson áp dụng cho trường hợp xoắn dầm chiều dài l, mặc cắt hình dạng bất kỳ như sau: ∇2u = - C khi viết dạng khai triển có dạng: 2 2 2 2 y u x u ∂ ∂ ∂ ∂ + = - 2Gθ trong đó G- mođun xoắn, θ - góc xoắn. (1) Phiếm hàm tương ứng phương trình Poisson có dạng: ∫ ∫+⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛+⎟⎠ ⎞⎜⎝ ⎛= V V Cudxdydxdy y u x uI 22 2 1 ∂ ∂ ∂ ∂ (2) Nếu biểu diễn: 13 u(x,y) = [N(x,y)]{δ} = [ Ni Nj Nk] ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ k j i u u u với Nm = (am bmx cm y) / 2Ae, m= i, j, k có thể tìm ma trận cứng phần tử và vector lực dưới dạng: dxdy y N y N x N x N k eA jiji ij ∫ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ += ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ Triển khai công thức trên cho m= i, j, k ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ + ++ +++ == )( )()( )()()( 4 1][ 22 22 22 kk kjkjjj kikijijiii e cbDX ccbbcb ccbbccbbcb A k . Vector lực: { } ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = 1 1 1 3 θGP . Tập họp ma trận cứng và vector lực thành hệ: [K]{u} - {P} = 0. trong đó: [K] = [k]∑ = E e 1 e ; {P} = {p}∑ = E e 1 e Sau khi giải hệ phương trình đại số sẽ nhận được: {u} = [K]-1{P} (3) Xác định momen xoắn và góc xoắn: MZ = 2Gθl2 ∫ A udxdy còn góc xoắn: θ = Mz / 2Gl2 ∫ A udxdy Ứng suất do xoắn: ⎭⎬ ⎫ ⎩⎨ ⎧ yz xz τ τ = Gθl2 2. Thứ tự giải bài toán cơ học kết cấu theo phương pháp PTHH 14 Từ cơ sở cách làm trong phương pháp chuyển vị, phương trình cân bằng được xác lập từ phương trình thế năng của đối tượng đang được nghiên cứu. Ẩn bài toán của phương pháp này là chuyển vị các nút và các ẩn này tham gia vào hàm thế năng. Theo cách đặt vấn đề từ các chương trước, điều kiện để thế năng đạt gía trị nhỏ nhất đạo hàm của hàm thế năng tính theo các bậc tự do các nút phải bằng 0. Thỏa mãn điều kiện này chúng ta nhận được phương trình (hệ phương trình) cân bằng. Những bước thực hiện như sau. 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ể. Bước 2: Mô hình chuyển vị trong mỗi phần tử tìm dưới dạng vector: { } [ ]{ }ee Nu δ= trong đó [N] – ma trận các 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 toá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 ee W 1 )(π Bước 4: Xử lý hệ phương trình và giải hệ phương trình đại số tuyến tính ghi tại (m). 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ệ giữa ứng suất – biến dạng. 15 Kết cấu thật Mô hình hóa Tính từng phần tử Tập họp Kết quả cuối cùng Hình 1.5 16