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 2)

Ví dụ 3: Thực hiện các bước tính theo phương pháp phần tử hữu hạn, xác định ứng suất và biến dạng dầm bị kéo. Dầm thẳng, đặt trong hệ tọa độ x0y mà trục 0x trùng với trục dọc của dầm. Với mô hình giản đơn này không xuất hiện các trường hợp chuyển đổi hoặc hoán vị góc xoay hệ tọa độ. Những vấn đề phức tạp khác nẩy sinh trong khi áp dụng PPPTHH sẽ được giải đáp trong phần sau của cùng tài liệu. Theo phương pháp PTHH kết cấu từ mô hình tính được chia làm nhiều phần tử. Tiến hành giải bài toán trong phạm vi mỗi phần tử sau đó tổng hợp kết quả cho bài toán chung. Dưới tác động của lực dọc trục đặt tại đầu bên phải dầm chiều dài L = l1 + l2, dầm bị biến dạng và ứng suất xuất hiện trong dầm. Xác định biến dạng và ứng suất trong dầm theo thứ tự sau 1. l1 l2 Phần tư û 1 Phần tư û 2 Hình 1.6 Dầm chịu kéo, nén. a) Rời rạc hóa không gian của bài toán Dầm được chia thành hai phần tử, mang chỉ số 1 và 2. Mỗi phần tử được đánh dấu bằng hai nút, nút đầu mang chỉ số i, nút sau mang chỉ số j. Trên hình vẽ phần tử đầu được đánh dấu với nút 1 bên trái, nút 2 bên phải. Phần tử kế tiếp có cách đánh dấu tương tự là nút 2 và 3. b) Mô hình chuyển vị Chuyển vị trong mỗi phần tử được miêu tả dạng hàm tuyến tính, ví dụ: u(x) = a + bx (a) trong đó các hệ số a, b là hằng số. 1 Tài liệu sử dụng cho phần này chọn trong các sách: 1) C.S.Desai,”Elementary Finite Element Method”, Prentice-Hall,1979 2) K.H.Huebner and E.A.Thornton, “The Finite Element Method for Engineers”, Wiley, N.Y.,1982 3) J.S.Przemieniecki, “Thery of Matrix Structural Analysis”, McGraw-Hill,N.Y.,1968. 4) O.C.Zienkievicz and Y.K.Cheung,”The Finite Element Method in Structural and Continuum Mechanics”,McGraw-Hill, London, 1967. 5) G. Dhatt, G. Touzot, “Une présentation de la méthode des éléments finis”, Paris,1984 17 Nếu ký hiệu chuyển vị nút 1 là u1, của nút 2 là u2, còn chiều dài phần tử l, chúng ta có thể viết: ⎭⎬ ⎫ += += 22 11 bxau bxau , với x1 = 0; x2 = l. Sau khi giải hệ phương trình sẽ nhận được: a = u1 ; b = ( u2 - u1)/l (b) từ đó: u(x) = u1 + (u2 - u1). l x trong phần tử 1, với l = l1. (c) Tương tự vậy trong phần tử số 2 chúng ta có: u(x) = u2 + (u3 - u2).x/l2. (d) c) Ma trận cứng Trong trường hợp cụ thềû này, ma trận cứng của phần tử xác lập từ nguyên lý năng lượng tối thiểu của hệ cân bằng. Năng lượng này bằng tổng của công biến dạng và công ngoại lực tác động lên dầm: Π = U - W (e) trong đó ∫∫ == LL dxAEdxAU 0 2 0 2 1 2 εσε (f) A - diện tích mặt cắt ngang phần tử, L - chiều dài, E- mođun đàn hồi vật liệu, σ - ứng suất, ε - biến dạng, Công thức trên đây đúng cho mỗi phần tử của dầm đang xét. Nếu coi năng lượng của hệ là tổng năng lượng tất cả phần tử cấu thành: . ∑ = =Π E e e 1 π Biến dạng trong phần tử tính bằng đạo hàm chuyển vị dx du e =ε : L uu 12 −=ε (g) Từ đó: 18 ( 212221 0 2 21 2 1 2 2 2 2 2 2 uuuu l AEdx L uuuuAEU l e −+=−+= ∫ ) (h) Mặt khác biểu thức Ue còn được hiểu là Ue = 2 1 {δ}eT[k]e{δ}e (i) với {δ}e = ⎭⎬ ⎫ ⎩⎨ ⎧ 2 1 u u e - chuyển vị hai nút của phần tử. Từ đó tính được [k]: [ ] ⎥⎦ ⎤⎢⎣ ⎡ − −= 11 11 l AEke - gọi là ma trận cứng phần tử. (j) Công ngoại lực tác động được tính: We = u1P1 + u2P2 cho một phần tử. Trong bài toán trên với hai phần tử, tổng cộng 3 nút, lực trên ba nút được ký hiệu trong vector {P} = [ p1 p2 p3 ]T, còn Wp = u1P1 + u2P2 + u3P3. Trên cơ sở giả thuyết năng lượng của hệ bằng tổng năng lượng tất cả phần tử cấu thành Π = ∑πe, nguyên lý năng lượng tối thiểu của hệ cân bằng cho phép viết: 0=Π iu∂ ∂ ; i=1,2,3 (k) Phương trình cuối được viết đầy đủ như sau: 0 2 =⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ −=Π ∑ pe ii W uu π∂ ∂ ∂ ∂ ; i= 1,2,3 (l) hoặc là: ∑ = 2 1e ( [k]e{δ}e - {p}e ) = 0 (m) d) Tập họp ma trận cứng và vector lực cho toàn hệ Trong giai đoạn này tiến hành tập họp ma trận cứng toàn hệ đồng thời với tập họp vector lực tác động lên hệ theo (m). Để thuận lợi lúc đọc tiếp chúng ta hãy gán những giá trị xác định cho kết cấu trên đây và tính các đặc trưng kết cấu bằng số. Giả sử A1 = 2cm2, A2 = 1 cm2, l1 = l2 = 10 cm, E1 = E2 =2x106 kG/cm2, P3 = 1 kG. Thành phần hai ma trận cứng từ hai phần tử có giá trị sau: 19 u1 u2 u2 u3 [k]1 = 105. và [k]⎥⎦ ⎤⎢⎣ ⎡ − − 44 44 2 = .10⎥⎦ ⎤⎢⎣ ⎡ − − 22 22 5 Sau khi tập họp [k]1 với [k]2, ma trận cứng toàn hệ có dạng: u1 u2 u3 u1 u2 u3 [K] = 105 = 2.10 3 2 1 22 2)24(4 044 u u u ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ − −+− − 5 3 2 1 110 132 022 u u u ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ − −− − (n) Vector lực có dạng: { } ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = 1 0 R P , trong đó R - phản lực tại nút 1. (o) e) Giải hệ phương trình đại số, xác định chuyển vị các nút Hệ phương trình giờ đây có dạng: [K]{Δ} = {P} (p) trong đó {Δ} - chuyển vị các nút 1, 2, 3 của hệ, là ẩn số cần tìm. Từ kết cấu của hệ có thể nhận biết, điều kiện biên của bài toán cho phép nhận chuyển vị nút 1: u1 bằng không. Xử lý điều kiện biên trước khi giải hệ phương trình chúng ta nhận được: 2.105. ⎭⎬ ⎫ ⎩⎨ ⎧= ⎭⎬ ⎫ ⎩⎨ ⎧⎥⎦ ⎤⎢⎣ ⎡ − − 1 0 11 13 3 2 u u Từ đó: u2 = 0,25.10-5 cm, và u3 = 0,75.10-5 cm. Xác định biến dạng và ứng suất trong phần tử Tại mỗi phần tử, về sau này chúng ta sẽ sử dụng đến khái niệm hệ toạ độ cục bộ để chỉ trường hợp liên quan đến mỗi phần tử, biến dạng và ứng suất được tính nhờ vào chuyển vị vừa thu nhận được từ bước (e). trong phần tử số 1: ε = dx du = (u2 - u1)/l = 0,25.10-6 trong phần tử số 2: ε = dx du = (u3 - u2)/l = 0,50.10-6 Ứng suất kéo trong phần tử tính bằng biểu thức σ = E. ε trong phần tử 1: σ = (2. 106)(0,25. 10-6) = 0,5 kG/cm2 20 trong phần tử 2: σ = (2. 106)(0,50. 10-6) = 1,0 kG/cm2. Ví dụ 4: Áp dụng nguyên lý công ảo xác định chuyển vị dầm liên tục, ngàm tại hai đầu. Dầm dài L =2l, độ cứng dầm EJ. Dầm chịu tác động lực tập trung P đặt tại giữa dầm. Mô hình hoá dầm như tại hình 1.7. Công thức tính mang dạng: T i T V S S T V T URdSUPdVUPdVD }{}{}{}{}{}{}]{[}{ 1 δδδεεδ ∑∫∫∫ ∫∫∫∫∫ ++= Hình 1.7 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. 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 hàm nội suy: N1 = H01(1) (x) = 3 1 l (2x3 -3l x2 + l3) N3 = H02(1) (x) = - 3 1 l (2x3 -3lx2) N2 = H11(1) (x) = 2 1 l (x3 -2lx2 +l2x) N4 = H12(1) (x) = l 1 (x3 -lx2) (b) Nếu thay ξ = x/l các hàm trên có dạng: N1 = 2ξ3 - 3ξ2 + 1 N2 = l(ξ3 -2 ξ2 + ξ) N3 = -2ξ3 +3 ξ2 21 N4 = l(ξ3 - ξ2) (b’) Vector chuyển vị hai nút của phần tử thứ nhất {δ} = [ u1 u2 u3 u4 ]T, của phần tử thứ hai {δ} = [ u3 u4 u5 u6 ]T, trong đó từ điều kiện ngàm hai đầu có thể thấy u1 = u1 = u5 = u6 = 0. Như đã trình bày trên, biểu thức U tính bằng công thức: dx dx wdEJdx xEJ xMU L L 2 0 2 22 2 1 )( )( 2 1 ∫∫ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛== (c) Biểu thức công ngọai lực: (d) ∫= L dxxPwW 0 )( [ ]{ }δ2 2 2 2 dx Nd dx wd = = [N’’]{δ}, và do vậy: Với mỗi phần tử πe = 2 1 {δ}T [N’’]∫l0 T EJ [N’’]dx {δ} = 21 {δ}T [k]{δ} (e) và từ đó: [k]e = [N’’]∫l0 T EJ [N’’]dx (f) Có thể tính các thành phần ma trận [k] như sau. ⎥⎦ ⎤⎢⎣ ⎡= ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − − − = 2221 1211 2 22 3 4 612 264 612612 ][ kk kk lDX l lll ll l EJk e (g) Tại phần tử 1, hàm chứa chuyển vị nút 1 và 2 trượt tiêu, chúng ta giữ 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, chỉ còn lại bốn thành phần, ứng với u3 và u4: 3 4 4 3 46 612 ][ 31 ⎥⎦ ⎤⎢⎣ ⎡ − −= l EJk (h) Với phần tử số hai ma trận cứng, giữ lại phần ứng với u3 và u4 có dạng: 3 4 4 3 46 612 ][ 32 ⎥⎦ ⎤⎢⎣ ⎡= l EJk (h’) Tập họp ma trận cứng từ hai ma trận phần tử: [ ] ⎥⎦ ⎤⎢⎣ ⎡ + +−+== ∑ = 44 661212 ][ 3 2 1 DXl EJkK e e (i) 22 Từ biểu thức U = 2 1 {u}T [K]{u }, có thể viết biểu thức δU cho hàm U: δU = δ{u}T [K]{u} (j) Trong trường hợp này vector chuyển vị có dạng: { } ⎭⎬ ⎫ ⎩⎨ ⎧= 4 3 u u u (k) Đại lượng δV của công ngoại lực suy từ biểu thức W = , bằng: ∫L dxxPw0 )( δW = -δ{u}T P/2 (l) Thỏa mãn điều kiện δU - δW = 0 có thể viết: { } 0 02 1 80 024 2 8 42 1 3 3 =⎭⎬ ⎫ ⎩⎨ ⎧− ⎭⎬ ⎫ ⎩⎨ ⎧⎥⎦ ⎤⎢⎣ ⎡ P u u L EJu Tδ (m) Vì rằng δ{u}T ≠ 0, biểu thức trong ngoặc phải bằng 0. ⎭⎬ ⎫ ⎩⎨ ⎧= ⎭⎬ ⎫ ⎩⎨ ⎧⎥⎦ ⎤⎢⎣ ⎡ 080 0248 42 1 3 3 P u u L EJ Từ đây nhận được giá trị cho u3 và u4. 0 ; 192 4 3 3 = = u EJ PLu (n) Bài toán động lực học Áp dụng nguyên lý công ảo xác định tần số dao động riêng dầm liên tục, dài L = 2l, độ cứng dầm EJ, nêu tại hình .. . Khối lượng riêng của dầm ρ, diện tích mặt cắt ngang dầm A. Dầm bị ngàm hai đầu. 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 như tại hình đã nêu trong bài toán tĩnh. Trong bài toán dao động tự do, ngoại lực được gán bằng 0, hàm công ngoại lực W cũng trượt tiêu. Thay vào đó động năng của hệ thống xuất hiện trong quá trình rung, tham gia vào phương trình năng lượng. Nếu ký hiệu U thế năng của hệ, T – động năng cũng của cùng hệ ấy, phương trình Lagrange áp dụng cho hệ dao động có dạng: 0=∂ ∂+⎟⎠ ⎞⎜⎝ ⎛ ∂ ∂ w U w T dt d & , trong đó ∫= L dxwxmT 0 2)( 2 1 & , còn U tính theo phần trên sẽ là dx dx wdEJdx xEJ xM L L 2 0 2 22 2 1 )( )( 2 1 ∫∫ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛= 23 Hàm thế năng và động năng trong mỗi phần tử được viết thành: ∫= 1 0 2 2 ξρ dwAlTe & (a) ξξπ dd wdEJe 2 2 21 02 1 ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛= ∫ ; (b) Nếu ký hiệu theo cách quen thuộc, hàm năng lượng của phần tử được viết thành: { } { } { } { }δδπ δδ ][ ;][ K MT e T e = = && (c) Phương trình cân bằng bài toán động lực học có dạng: }0{}]{[}]{[ =+ wKwM && (d) Ma trận [K] đã được giải trình tại (i) trong phần tĩnh học. Ma trận khối lượng [M] được tính theo cùng cách thức. Biểu thức tính ma trận khối lượng theo cách giải trong phần này có dạng: ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − − − = 2 22 6 33234 5,45,19/6 5,198133234 630 ][ lDX l ll ll Alm e ρ (e) Thay biểu thức tính we vào công thức ∫= 1 0 2 2 ξρ dwAlTe & Trong đó: 44332211 uNuNuNuNw &&&&& +++= Với phần tử thứ nhất 021 == uu && , do vậy 4433 uNuNw &&& += ( ) ( )2443231 0 443321 633223463022 uuuxu x AlduNuNAlT &&&&&& +−=+= ∫− ρξρ (f) Cân bằng biểu thức cuối với { } { };][ δδ && MT Te = có thể viết biểu thức của [M]e. ⎥⎦ ⎤⎢⎣ ⎡ − −=− 633 33234 639221 x AlM ρ (g) ⎥⎦ ⎤⎢⎣ ⎡ + +=− 633 33234 639232 x AlM ρ Sau tập họp ma trận khối lượng toàn hệ sẽ là: 24 ⎥⎦ ⎤⎢⎣ ⎡= 120 0468 6302 ][ x AlM ρ (h) Ma trận [K] như đã biết, mang giá trị [ ] ⎥⎦ ⎤⎢⎣ ⎡= 80 024 3l EJK Nếu xét dao động riêng của hệ dạng wi = Aicosωt, trong đó A – biên độ dao động, ω - tần số dao động riêng, phương trình dao động trở về dạng: -[M]ω2 + [K] = {0} (i) hay [K]{w} = λ[M]{w} với {w} = {u3 u4 }T (j) Lời giải hệ phương trình cuối sẽ là: ⎭⎬ ⎫ ⎩⎨ ⎧= ⎭⎬ ⎫ ⎩⎨ ⎧= 0 1 ;7,22 4 3 21 u u m EJ L ω (k) ⎭⎬ ⎫ ⎩⎨ ⎧= ⎭⎬ ⎫ ⎩⎨ ⎧= 1 0 ;0,82 4 3 22 u u m EJ L ω (k’) trong đó m – khối lượng của dầm dài L đang đề cập. Ổn định dầm Áp dụng nguyên lý công ảo xác định hệ số ảnh hưởng của gối đỡ χ≡ k2 trong bài toán ổn định dầm. Dầm được xét có kích thước hình học như đã nêu trong ví dụ 1. Dầm bị ngàm hai đầu như ví dụ trình bày tại hình . Hệ số ảnh hưởng của gối đỡ χ≡ k2 trong công thức xác định lực giới hạn Pcryt như sau: Pcryt = k2EJ/L2. Biểu thức: dx dx wdEJdx xEJ xMU L L 2 0 2 22 2 1 )( )( 2 1 ∫∫ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛== . Biểu thức công ngọai lực trong quá trình dẫn đến mất ổn định dầm tính bằng công thức: ∫= L dxxwPW 0 2)]('[2 . Vector chuyển vị hai nút của phần tử thứ nhất {δ} = [u1 u2 u3 u4]T, của phần tử thứ hai {δ} = [u3 u4 u5 u6]T, trong đó từ điều kiện ngàm hai đầu có thể thấy u1 = u2 = u5 = u6 = 0. Bằng cách tương tự đã thực hiện cho bài toán tính chuyển vị nêu trước, thay các giá trị tương ứng vào biểu thức U1-2, có thể xác định ma trận cứng [k]1 và [k]2ø: ⎥⎦ ⎤⎢⎣ ⎡ − −= 46 612 ][ 31 l EJk ; ⎥⎦ ⎤⎢⎣ ⎡ + += 46 612 ][ 32 l EJk . Tập họp ma trận cứng cho toàn hệ có dạng: [K] = ⎥⎦ ⎤⎢⎣ ⎡=⎥⎦ ⎤⎢⎣ ⎡ + +−+=∑ = 80 024 )44( 661212 ][ 33 2 1 l EJ DXl EJk e e 25 Công ngoại lực trong mỗi phần tử tính theo cách sau: ∫ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ +−=− 1 0 2 4 4 3 3 21 2 ξξξ dud dNu d dN l PW Nếu viết { } }{][ 2 1 2121 δδ −− = SW T với [S] ma trận cứng hình học của phần tử 1: ⎥⎦ ⎤⎢⎣ ⎡ − −−=− 15 2 10 1 10 1 5 6 21 l PS Tương tự vậy: ∫ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ +−=− 1 0 2 4 2 3 1 32 2 ξξξ dud dNu d dN l PW và ⎥⎦ ⎤⎢⎣ ⎡−=− 15 2 10 1 10 1 5 6 32 l PS Ma trận [S] cho toàn hệ có dạng: [S] = ⎥⎦ ⎤⎢⎣ ⎡=⎥⎦ ⎤⎢⎣ ⎡ + −−= + = ∑ 15 4 5 12 15 2 15 2 10 1 10 1 5 662 1 0 0 )( ][ DXl PS e e Từ quan hệ δ(U -W) = δ{u}T( [K] –[ S]) {u} có thể viết: ( [K] – [S] ) {u} = 0 Hệ phương trình cuối được khai triển dưới dạng: ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − − =⎥⎦ ⎤⎢⎣ ⎡−⎥⎦ ⎤⎢⎣ ⎡ 5 480 0 5 1224 0 0 80 024 2 2 15 4 5 122 EJ Pl EJ Pl EJ Pl Tại vị trí l = L/2: 0 15 18 5 324 22 2 1 =⎟⎠ ⎞⎜⎝ ⎛ −⎟⎠ ⎞⎜⎝ ⎛ − kk Từ phương trình cuối có thể tính được: ;40 2 12 11 === EJ LPkχ ;120 2 22 22 === EJ LPkχ 3. Phương pháp đưa hàm sai số về bằng 0 Phương trình cân bằng dạng tổng quát có thể diễn đạt như sau: 26 ‹ (u) = p trong miền V (a*) và thoả mãn các điều kiện biên: B(u) = q trên biên S (b*) Biến u được tìm ở dạng hàm xấp xỉ: ∑ = = N i ii fau 1 * (c*) trong đó ai - các hằng số, còn fi các hàm được chọn nhằm thỏa mãn tất cả điều kiện biên đã kể. Đại lượng thường ký hiệu bằng R - số dư, hoặc ε - hàm sai số dùng để miêu tả sai số sau khi thay u trong (a*) bằng biểu thức (c*) : R = ‹ ( ) - p (d*) ∑ = N i ii fa 1 Trong các phương pháp tính vì rằng R rất khó đạt hoặc không thể đạt 0 người ta tìm mọi cách để R → min hoặc tồn tại với một gía trị có ý nghĩa thực tế. Để xử lý bài toán dạng này trong trường hợp chung tiến hành giải phương trình sau: ∫ V f(R).w.dV = 0. (e*) trong đó w - hàm trọng lượng, f(R) hàm phụ thuộc vào R. Trong nhiều trường hợp kể dưới đây người ta nhận f(R) = R. Phương pháp chọn các điểm nhằm thỏa mãn biểu thức (a*) tiến hành theo phương thức đưa vào (e*) hai biểu thức sau: f(R) = R và w = δ( x - ξ) (f*) Phương pháp Galerkin Cách làm của Galerkin là đặt các hàm fi của (c*) trực giao với hàm R. ∫ V fi R dV = 0. Phương pháp tổng nhỏ nhất các bình phương V ∫ w. R 2 dV = min. Đây là phương pháp được nghĩ ra sớm nhất trong các phương pháp số nhằm thực hiện hàm hoá, tìm nghiệm các bài toán, xử lý kết quả đo đạc. ∫ V w. R 2 dV = w{‹ ( ) - p }∫ V ∑ = N i ii fa 1 2 dV = min. 27 trong đó ẩn số của hệ phương trình chỉ là ai. Điều kiện cần để tích phân này đạt minimum là: ia∂ ∂ ∫ V w{‹ (∑ ) - p } = N i ii fa 1 2 dV = 0 với i=1,2,....,N. Sau khi lấy đạo hàm theo ai hệ phương trình đại số có dạng: ∫ V w. ‹ (fi) {‹ ( ) - p } dV =0 ∑ = N i ii fa 1 Dưới đây sẽ trình bày ví dụ sử dụng phương pháp tính đã giới thiệu tại ứng dụng trong khuôn khổ phương pháp phần tử hữu hạn. Ví dụ 1: Giải phương trình vi phân (a) ghi dưới đây theo phương pháp PTHH, cách giải Galerkin. 02 2 =++ xu dx ud trong miền 0 ≤ x ≤ 1 (a) Điều kiện biên: u(0) = u(1) =0; (b) Trong trường hợp này phương trình R có dạng: R = xu dx ud ++2 2 (c) Hàm u được hàm hóa dạng hàm tuyến tính: u(x) = Ni(x) ui - Nj(x) uj = [N(x)]{δ} (d) trong đó: Ni(x) = l xx j − và Nj = l xxi − , {δ} = [ u1 u2 ]T. Phương trình Galerkin: 0)( 1 0 2 2 =⎥⎦ ⎤⎢⎣ ⎡ ++∫ dxxNxudxud k , k = i, j. (e) Trong mỗi phần tử tính từ xi, đến xj tích phân Galerkin sẽ là: [ ] 0)(1 2 2 =⎥⎦ ⎤⎢⎣ ⎡ ++∫ dxxudxudxN x x T (f) Vế đầu của tích phân có thể giải như sau: 28 [ ] ∫∫ −= jj x x Tx x x x T dx dx du dx xNd dx duxNdx dx udxN )]([)()]([ 1 2 2 (g) Mặt khác: ⎭⎬ ⎫ ⎩⎨ ⎧−= ⎭⎬ ⎫ ⎩⎨ ⎧ − −= 1 11 /)( /)()]([ llxx lxx dx d dx xNd i j T [ ]{ } [ ] ⎭⎬ ⎫ ⎩⎨ ⎧−== 2 1111 u u l N dx d dx du δ Thay thế (g) vào (f) sẽ nhận được: [N(x)] j i x xdx du - ∫ jx x T dx dx du dx xNd )]([ - [N(x)]T{δ} – [N(x)]T x }dx = 0 (i) Hoặc là [k]{δ} = {p} , tính cho mỗi phần tử. (j) Sau khi thực hiện tích phân trên, biểu thức của [k]e và {p}e sẽ như sau: [ ] ⎥⎦ ⎤⎢⎣ ⎡−⎥⎦ ⎤⎢⎣ ⎡ −−= 21 12 611 111 l l ke (k) { } (l) ⎪⎭ ⎪⎬⎫⎪⎩ ⎪⎨⎧ −− −+= )2( )2( 22 22 ijij ijij e xxxx xxxx p Các phép tính tiếp theo tiến hành như đã thực hiện tại phần bàn về phương pháp biến phân. Ví dụ 2: Giải phương trình Laplace trong miền hạn chế Ω theo công thức Galerkin. Phương trình Laplace: Ω=+=∇ trong y u x uu 02 2 2 2 2 ∂ ∂ ∂ ∂ (a) Tại biên S = S1 + S2 cần thỏa mãn các điều kiệïn Dirichlet và điều kiện Neumann. Điều kiện Dirichlet: u = u0 trên S1, Điều kiện Neumann: 0v=n u ∂ ∂ trên S2. (b) Các bước tiến hành: (1) Chia miền Ω ra làm E phần tử, ví dụ các phần tử tam giác, tứ giác, Ω = e1 + e2 +.... + eN 29 (2) Xác lập {u} = [N]{δ} trong đó, với phần tử tam giác công thức tính [N] có dạng, xem phần hàm nội suy: [ ] ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ++ ++ ++ = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = Aycxba Aycxba Aycxba N N N N kkk jjj iii k j i 2/)( 2/)( 2/)( (c) { } ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = k j i u u u δ (d) (3) Áp dụng công thức Galerkin cho bài toán tổng quát: ∫ Ω Nk(x,y).∇2u dA = 0 trong miền Ω, k = 1,2,3,... (e) (4) Giải tích phân cuối theo công thức thứ nhất của Green: ∫ Ω Nk(x,y).∇2u dA = N∫ S k(x,y) ∂ ∂ u n .dS- gradN∫ Ω k(x,y).gradu dA = 0 (f) Hoặc là: ∫ Ω gradNk(x,y).gradu dA + N∫ S k(x,y).v0d S1 = 0; (f’) (5) Tính ma trận “cứng” và vector “lực” Tích phân thứ nhất trong mỗi phần tử đóng vai trò ma trận cứng, có dạng: ∫∫ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ +=⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ + AA kk dAN y N x dA yy N xx N e }{][][}{}{ 2 2 2 2 δ∂ ∂ ∂ ∂ ∂ δ∂ ∂ ∂ ∂ δ∂ ∂ ∂ (g) Thay {N} và {δ} vào biểu thức cuối, đưa biểu thức về dạng ma trận cứng quen thuộc, biểu thức [k] có dạng: [k] = [B]∫ eA T [D] [B] dA (h) trong đó: [ ] ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ = y N y N y N x N x N x N B kji kji ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ và [ ] ⎥⎦ ⎤⎢⎣ ⎡= 10 01 D từ đó: 30 [ ] ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ + ++ +++ = )( )()( )()()( 4 1 22 22 22 kk kjkjjj kikijijiii cb ccbbcbDX ccbbccbbcb A k (i) Biểu thức {p} = N∫ S k(x,y).v0dS1 tính dọc biên Si - Sj có dạng: ∫ S Nk(x,y).v0dS1 = -v0 ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ −= ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ∫ 0 1 1 2 v 0 0 ijSj Si j i s dsN N (j) Các bước tiếp theo gồm tập họp matrận “cứng”, vector “lực”, giải hệ phương trình đại số vv... thực hiện theo cách quen thuộc. 31