Bài giảng Phương pháp số (phương pháp phần tử hữu hạn)

91 trang vanle 3170

Download

Bạn đang xem 20 trang mẫu của tài liệu "Bài giảng Phương pháp số (phương pháp phần tử hữu hạn)", để 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:

bai_giang_phuong_phap_so_phuong_phap_phan_tu_huu_han.pdf

Nội dung text: Bài giảng Phương pháp số (phương pháp phần tử hữu hạn)

Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . LỜI NÓI ĐẦU Để giải và tính toán các bài toán về kêt cấu cơ học, ngoài các phương pháp giải tích ta còn có các phương pháp số. Do các bài toán cơ học thường dẫn đến việc giải các phương trình vi phân với các điều kiện biên xác định nào đó. Vì vậy thời kỳ đầu của các phương pháp số là : các phương pháp tích phân số và phương pháp sai phân hữu hạn. Cùng với sự phát triển của máy tính điện tử, phương pháp phần tử hữu hạn ra đời và phát triển rất mạnh mẽ và là một phương pháp được dùng rất phổ biến hiện nay khi tính toán các bài toán cơ học. Nó cũng đã được áp dụng để có được nhiều chương trình tính cho các dạng bài toán cơ học khác nhau: Tính cho dàn thanh, khung không gian, các kết cấu dạng tấm , vỏ , Phương pháp phần tử hữu hạn là môn học cơ sở của các ngành kỹ thuật liên quan đến tính toán các kết cấu và hiện cũng là một môn học của ngành Xây dựng và Kỹ thuật công trình thuộc trường ĐHLN. Trong năm trước đây chúng tôi có biên soạn nội dung bài giảng : Phương pháp phần tử hữu hạn để phục vụ cho công tác giảng dạy môn học : Phương pháp số. Vẫn biết rằng tài liệu viết về môn học này đã có rất nhiều trên các dạng : sách , bài giảng và trên mạng, song thiết nghĩ thì việc biên soạn một tài liệu dạng bài giảng về phương pháp phần tử hữu hạn với thời lượng 2 tín chỉ cũng là điều cần thiết để các em sinh viên ( và cả các độc giả lần đầu biết về phương pháp này) tiếp cận với môn học này thuận lợi hơn. Tài liệu mới chỉ tiếp cận đến một số nội dung và khái niệm cơ bản của phương pháp phần tử hữu hạn. Các vấn đề trình bày mới dừng đến việc tính toán cho dàn, khung không gian. Tài liệu cũng đã đưa ra một số thủ tục cơ bản trong lập trình tính toán, các thủ tục này được viết trong Visual Basic, độc giả có thể chuyển đổi dễ dàng sang các môi trường lập trình khác. Mong rằng với ý muốn như thế sẽ giúp ích được phần nào cho quá trình học tập môn học này của các em sinh viên, và tất nhiên rất mong được các đóng góp của độc giả về các vấn đề trình bài trong tài liệu. Tác giả 1 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Chương I MỘT SỐ VẤN ĐỀ CƠ BẢN TRONG CƠ HỌC VẬT RẮN BIẾN DẠNG & CÁC PHƯƠNG PHÁP TÍNH TRONG CƠ HỌC I.1 CÁC VẤN ĐỀ CƠ BẢN TRONG CƠ HỌC VẬT RẮN BIẾN DẠNG I.1.1 Ten xơ ứng suất. Dưới tác dụng của lực ngoài, vật thể chịu lực bị biến dạng và bên trong nó sẽ xuất hiện ứng suất. Ứng suất tại mỗi điểm khác nhau là khác nhau, véc tơ ứng suất không những phụ thuộc vào điểm mà còn phụ thuộc vào hướng của thiết diện qua nó mà được xác định bởi pháp tuyến có hướng n . Như vậy tập hợp cặp véc tơ ứng suất Tn và véc tơ n tại điểm P sẽ xác định trạng thái ứng suất tại điểm đó. Trạng thái ứng suất tại điểm hoàn toàn được xác định qua ten-xơ ứng suất – là một ten xơ đối xứng hạng hai, nên nó có 6 thành phần độc lập: σ11 σ 12 σ 13 σ σ σ σ với σ σ ij 21 22 23 ij ji σ31 σ 32 σ 33 Trong hệ tọa độ De-cac các thành phân của ten xơ ứng suất được ký hiệu là : σx ;σ;σ;τ y z xy ;τ xz ;τ yz I.1.2 Phương trình cân bằng Tách phần thể tích V tùy ý giới hạn bởi mặt S của môi trường liên tục ở hình thái biến dạng, xét sự cân bằng các lực tác dụng lên thể tích đó ( không kể lực quán tính) ta được : hay là Tn dS K dV 0 Ti ni dS K dV 0 S V S V T Do i ( công thức Gaoxơ - Ôtrôgratxki) nên ta có : Ti ni dS dV S V xi Ti K dV , vì V là thể tích tùy ý nên biểu thức dưới dấu tích phân bằng không V xi Ti  ij => ta được : K 0 hay lµ K j 0 (I.1) xi xi 2 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Phương trình (I.1) được gọi là phương trình cân bằng, ρ là mật độ khối lượng , K lực khối. Viết (I.1) dưới dạng tường minh ta được : σ x τ xy τ xz ρK x 0 x  y  z τyx  σ y  τ yz ρK y 0 (I.2) x  y  z τzxτzy  σ z ρK z 0 x  y  z Với bài toán hai chiều ( Tấm , vỏ ), phương trình cân bằng có dạng : σ x τ xy ρK x 0 x  y (I.3) τyx  σ y ρK y 0 x  y Còn trong bài toán một chiều phương trình cân bằng sẽ là : σ x ρK 0 (I.4) x x I.1.3 Quan hệ giữa biến dạng – chuyển vị ( hệ thức Cô-si) Khi xây dựng hệ thức quan hệ giữa biến dạng và chuyển vị, xuất phát từ sự thể hiện thay đổi kích thước dài đoạn vô cùng nhỏ bất kỳ dX lấy từ điểm X và thay đổi góc giữa hai đoạn vô cùng nhỏ bất kỳ lấy từ điểm đó người ta đã dẫn đến ten-xơ biến dạng hữu hạn viết trong hệ tọa độ Đề các : Grin và Anmăngxi 1  UU  U  U γ i j k k (Grin) (I.5) ij 2  Xj  X i  X i  X j 1  uu  u  u γ i j k k (Anmăngxi) (I.6) ij 2  xj  x i  x i  x j Trong đó Xn - biến Lagrăng, xk – biến Ơ le , Um và un là các thành phần chuyển vị theo biến Lagrăng và Ơle. 3 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Trường hợp biến dạng nhỏ ( t.là khi bỏ qua VCB bậc cao – là các thành phần phi tuyến trong (1.5) và (1.6) )) khi đó hai ten xơ này xấp xỉ bằng nhau và các thành phần của ten- xơ biến dạng nhỏ sẽ có dạng:    11 12 13 u1u3  u 2  ij  21  22  23 với các thành phần : ε11 , ε 33 , ε 22 , x1  x 3  x 2  31  32  33 1  u1  u 2 1  u 2u3 1  u 1  u 3 ε12 , ε 23 , ε 13 2  x2  x 1 2  x 3  x 2 2  x 3  x 1 Để thuận lợi cho các công thức sau này trong tính toán theo phương pháp PTHH, người ta ký hiệu : - Các thành phần chuyển vị : u , v , w - Các thành phần của ten-xơ biến dạng ε x ; ε y ; ε z ; γ xy ; γ xz ; γ yz với : u  v  w ε ; ε ; ε ; xx y  y z  z u  v  u  w  v  w γ ; γ ; γ (I.7) xyy  x xz  z  x yz  z  y Hay viết dưới dạng ma trận :  0 0 x εx 0 0 y ε y 0 0  u ε z z (I.8) . v γxy   0 y  x w γyz 0   γxz z  y 0  z  x I.1.4 Phương trình liên tục Hệ thức Cô-si (I.8) cho sự liên hệ giữa 6 thành phần biến dạng xác định duy nhất qua 3 thành phần chuyển vị cho trước. Như vậy với 6 thành phần biến dạng cho trước chỉ từ quan hệ (I.8) sẽ không cho duy nhất 3 thành phần chuyển vị, do đó giữa các thành phần 4 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . biến dạng sẽ có 6 phương trình tương thích biến dạng ( còn gọi là các phương trình tương thích biến dạng Xanhvơnăng), các phương trình này đảm bảo cho sự biến dạng liên tục trong môi trường. 2 2 2 γ  γ 2  ε x εy  γ xy  γzx yx zy  ε x 2 2 2 y  x x  y x  y  z  x  y  z 22 2 2 εy εz  γ yz  γzy  γ yxγzx  ε y 2 2 2 (I.9) z  y z  y y  x  z  y  x  z 2 2 2 2  ε x εz  γ xz  γzxγyz  γ xy  ε z 2 2 2 z  x x  z z  y  x  z  y  x 2 2 2ε ε  γ  Trong bài toán 2 chiều : (I.9) còn 1 phương trình x y xy y2  x 2 x  y  Trong bài toán 1 chiều : các phương trình trên đều thỏa mãn. I.1.5 Điều kiện biên Điều kiện biên liên quan đến chuyển vị hoặc ứng suất  Điều kiện biên liên quan đến chuyển vị ( gọi là điều kiện biên động học) thường được cho trước chuyển vị của một điểm hoặc một phần mặt biên nào đó  Điều kiện biên liên quan đến ứng suất ( gọi là điều kiện biên tĩnh học) đòi hỏi sự cân bằng giữa ứng suất trên mặt biên với ngoại lực đặt lên đó. Ví dụ Một thanh chiều dài  , chiều dầy h, bị ngàm chặt một đầu, một đầu tự do, chịu tác dụng của lực phân bố đều có cường độ q như hình vẽ. Chọn hệ tọa độ : 0x theo chiều dài, sát mặt dưới, 0y hướng lên trên Hình 1.1 Đây là bài toán hai chiều và các điều kiện biên được đưa về các hệ thức sau : v(0, y) - Tại x = 0 : u(0,y) = v(0,y) = 0 ; 0 x 5 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . - Tại mặt trên ( y = h) : τxy (x,h) 0; σ(x,h) y q - Tại mặt dưới ( y = 0) : τxy (x,0) 0; σ(x,0) y 0 - Tại đầu B ( x =  ) : τxy (,y) 0;σ(,y) x  0 I.1.6 Phương trình vật lý (phương trình trạng thái)–Quan hệ ứng suất và biến dạng. Trong giáo trình này chỉ xét giai đoạn làm việc của vật liệu ở giai đoạn đàn hồi, biến dạng là nhỏ và đàn hồi là tuyến tính. Như vậy quan hệ giữa ứng suất và biến dạng ở đây được áp dụng bởi định luật Hooke. Xét với vật liệu đẳng hướng : I.1.6.1 Bài toán 3 chiều : định luật Hooke có dạng : 1 1 2(1 ν) ε σ ν(σ σ ) , γ τ τ xEGE x y z xy xy xy 1 1 2(1 ν) ε σ ν(σ σ ) , γ τ τ (I.10) yEGE y x z yz yz yz 1 1 2(1 ν) ε σ ν(σ σ ) , γ τ τ zEGE z x y xz xz xz Chú ý rằng các thành phần của ten-xơ biến dạng được tính theo chuyển vị qua (I.7) Nếu viết dưới dạng ma trận và có kể đến biến dạng ban đầu thì (I.10) có dạng: ε  C . σ ε0 (I.11) trong đó: T ε ε,ε,ε,γx y z xy ,γ yz ,γ zx  - là véc tơ biến dạng T σ σ,σ,σ,τx y z xy ,τ yz ,τ zx - là véc tơ ứng suất T ε0 ε 0x ,ε 0y ,ε 0z ,γ 0xy ,γ 0yz ,γ 0zx - là véc tơ biến dạng ban đầu ( Chữ T – ký hiệu chuyển vị ma trận) [C] – ma trận các hệ số đàn hồi , 6 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 1 ν ν 0 0 0 ν 1 ν 0 0 0 1 ν ν 1 0 0 0 C (I.12) E 0 0 0 2(1 ν) 0 0 0 0 0 0 2(1 ν) 0 0 0 0 0 0 2(1 ν) với E – mô đun đàn hồi Young , G – mô đun trượt , ν - hệ số Poát-xông của vật liệu. 0 T Trường hợp biến dạng ban đầu do nhiệt độ thì ε0 αT 1,1,1,0,0,0 , trong đó α - hệ số dãn nở vì nhiệt, T0 – độ biến thiên của nhiệt độ. Biểu diễn ứng suất qua các thành phần biến dạng ta sẽ có : σ  D ε ε0 hay là : 1 1 E αT0 1 σ  D ε 1 2ν 0 0 0 với ma trận hệ số D là : 1 ν ν ν 0 0 0 ν 1 ν ν 0 0 0 ν ν 1 ν 0 0 0 E 1 2ν D 0 0 0 0 0 (I.13) (1 ν).(1 2ν) 2 1 2ν 0 0 0 0 0 2 1 2ν 0 0 0 0 0 2 I.1.6.2 Bài toán 2 chiều :  Bài toán ứng suất phẳng : ví dụ như xét các bài toán về tấm, vỏ với tải trọng nằm trong mặt phẳng giữa tấm, phân bố đều theo bề dầy của tấm khi đó chọn trục z vuông góc với mặt phẳng tấm, sẽ dẫn đến thể giả thiết rằng : 7 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . σz τ xz τ yz 0 ứng suất không đổi theo chiều dầy của tấm. với giả thiết này các biểu thức của định luật Hooke có dạng : ε  C σ ε0 với : εx  σ x  ε 0x  1 1 ν 0 0 1 ε ε;σσ;ε     ε  αT1;C  ν 1 0 y y 0 0y E γxy  τ xy  γ 0xy  0 0 0 2(1 ν 0 1 EαT Hay biểu diễn ngược lại: σ  D ε ε0  D ε 1 1 ν 0 1 ν 0 E với D 2 ν 1 0 1 ν 1 ν 0 0 2 ν 0 biến dạng theo phương z vẫn tồn tại và εz σ x σ y αT E  Bài toán biến dạng phẳng : Khi xét vật thể hình lăng trụ dài có mặt cắt ngang không đổi theo chiều dài ( theo chiều trục 0z) , chịu tải trọng đều w vuông góc với 0z , khi đó ta có : w = 0 ; ε 0 ; các đại lượng z z ứng suất và biến dạng chỉ phụ thuộc vào các biến x và y. ε  C . σ ε0 ε0x  1 1 ν ν 0 0 1 ν ε ε  1ναT1;C   ν 1ν 0 0 0y E γ0xy  0 0 0 2 0 1 EαT Hay biểu diễn ngược lại: σ  D ε ε0  D ε 1 1 2ν 0 8 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 1 ν ν 0 E với D ν 1 ν 0 (1 ν).(1 2ν) 1 2ν 0 0 2 0 ứng suất theo phương z vẫn tồn tại và σz νσ x σ y EαT;τ xz τ yz 0 I.1.6.3 Bài toán 1 chiều : 1 ε σ αT0 hoặc σ Eε EεT0 ( => D = E – mô đun đàn hồi) xE x x x I.1.7 Đặt bài toán đàn hồi : Thiết lập bài toán đàn hồi bao gồm việc thiết lập các phương trình và các điều kiện biên, chúng phải lập thành một hệ kín để có thể giải ra được các ẩn cần tìm đó là các giá trị của các thành phần ten – xơ biến dạng, ứng suất , véc tơ chuyển vị. Các phương trình gồm có : Phương trình cân bằng Hệ thức Cô-si ( liên hệ chuyển vị và biến dạng) Phương trình trạng thái ( liên hệ ứng suất và biến dạng : định luật Hooke) và các điều kiện biên động học hoặc tĩnh học. Người ta đã chứng minh được sự tồn tại và duy nhất nghiệm của bài toán đàn hồi. I.2 CÁC PHƯƠNG PHÁP GIẢI BÀI TOÁN ĐÀN HỒI Như phần trên đã trình bày, ta thấy rằng để giải bài toán đàn hồi tuyến tính : 3 chiều ta có tới 15 phương trình cùng các điều kiện biên để tìm ra giá trị của 15 ẩn : 3 thành phần chuyển vị, 6 thành phần biến dạng và 6 thành phần ứng suất. Cho đến nay cũng đã có được các phương pháp giải đúng và gần đúng. Có những phương pháp có thể áp dụng tốt cho một lớp các dạng bài toán cơ học biến dạng nhưng lại khó khăn áp dụng nó cho dạng khác. Có thể tổng kết ra đây theo sơ đồ sau 9 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Các phương pháp giải Các phương pháp giải tích Các phương pháp số Các phương Các phương pháp Các phương pháp số Phương pháp pháp chính xác gần đúng ( biến giải các phương phần tử hữu hạn phân) trình vi phân ( PTHH ) Các phương Phương Mô hình Mô hình Mô hình pháp tích pháp sai tương cân bằng hỗn hợp phân số phân hữu thích hạn I.2.1 Phương pháp chính xác Phương pháp giải tích giải các bài toán đàn hồi : Có thể giải theo chuyển vị , hoặc giải theo ứng suất. ví dụ : xét bài toán uốn dầm như hình 1.2 : dầm chịu tải trọng đều với cường độ q , tựa gối khớp tại hai đầu A , B Lời giải bài toán trên đã có trong sức bền 0  x vật liệu : d4 w EJ q 0 ; ( 0 x  ) (1.14) dx4 z với điều kiện biên : hình 1.2 w(0) w( ) 0 Tích phân liên tiếp (1.14) ta được chuyển vị theo w (0) w ( ) 0 4 q x 3 2 phương z : w(x) Cx1 Cx 2 CxC 3 4 , từ các điều kiện biên ta xác định EJ 24 10 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .  3 được các hằng số tích phân : C2= C4 = 0 ; C;C . dẫn đến nghiệm bài toán 112 3 24 q theo chuyển vị w(x) x4 2 x 3  3 x => tìm được mô men uốn: 24EJ 2 4 d w q 2 5q M(x) EJ2 x  x ; độ võng lớn nhất tại x=  / 2 là wmax = dx 2 384EJ I.2.2 Các phương pháp biến phân : Các bài toán cơ học nói chung có thể biểu diễn dưới dạng tổng quát : L(u) g 0 trong miền V C(u) p trên biên S (1.15) trong đó L , C là các toán tử vi phân , u = u (x,y,z) đại lượng cần tìm , g = g(x,y,z) và p=p(x,y,z) là các hàm cho trước. Phương pháp biến phân là phương pháp gần đúng nhằm tìm một nghiệm xấp xỉ dựa trên một tiêu chuẩn nào đó. Các phương pháp biến phân thường gặp đó là phương pháp Ritz, phương pháp Galerkin, 1. Phương pháp Ritz : Trong một số các bài toán đàn hồi thường tồn tại một phiếm hàm I dạng : I F(x, y,u,ux ,u xx ,u y ,u yy )dxdy và từ điều kiện dừng của phiếm hàm này sẽ dẫn ra V các phương trình vi phân của bài toán. Phương pháp Ritz tìm nghiệm xấp xỉ gần đúng trong một dạng tổ hợp tuyến tính các hàm biết trước, chẳng hạn như trong bài toán phẳng , theo Ritz thì chuyển vị có dạng : N u(x, y) φ(x,y)0  Cφ(x,y) i i (1.16) i 1 trong đó :  φ 0(x,y) là hàm được chọn sao cho thỏa mãn các điều kiện biên không thuần nhất.  φ i(x,y) là các hàm thỏa mãn :  Khả vi đến cấp cần thiết 11 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .  Thỏa mãn các điều kiện biên thuần nhất  N Hệ hàm φii là độc lập tuyến tính và đầy đủ Với các điều kiện trên sẽ làm cho nghiệm của Ritz hội tụ đến nghiệm chính xác ( khi N ) Ci là các hằng số cần tìm, chúng được xác định khi thay u(x,y) vào phiếm hàm I, lấy tích phân trên V, khi đó I sẽ là hàm của các Ci, cho thỏa mãn điều kiện dừng của I (tức là các đạo hàm riêng của I theo các Ci bằng 0) sẽ dẫn đến hệ phương trình đại số xác định các hằng số Ci. 2. Phương pháp phần dư có trọng (Weighted Residual Method) Gọi R(u) = L(u) + g được gọi là hàm phần dư, như vậy nghiệm của bài toán L(u) + g = 0 sẽ được chuyển thành tìm u là nghiệm của bài toán R(u) = 0. Người ta cũng N tìm nghiệm dưới dạng uN φ 0  C i φ i , trong đó : i 1 φ 0 : được chọn sao cho thỏa mãn các điều khiện biên không thuần nhất φ i : khả vi đến số lần cần thiết , thoả mãn các điều kiện biên thuần nhất và hệ hàm { φ i} là độc lập tuyến tính. khi đó R(uN) = R(x,y,z,C1, C2, ,CN). Các hằng số Ci được tìm bằng phương pháp : lấy tích phân trên V của tích các hàm ψ k với R(uN) , tức là : ψk R(u N )dv 0 ; k 1,2, , N (1.17) V các hàm ψ k được gọi là các hàm trọng số (Weighted function) Chú ý rằng toán tử vi phân L trong (1.15) là tuyến tính nên thay R(uN) vào (1.17) ta sẽ dẫn đến : N ψR(u)dvk N  C i ψL(φ)dv k i ψ k L(φ) 0 gdv 0 ; k =1,2, ,N (1.18) VVVi 1 Nếu đặt ψk L(φ i )dv = Aki ; ψk L(φ 0 ) g dv = Bk khi đo ta sẽ có hệ V V phương trình đại số tuyến tính : {A}{C} = {B} xác định N hằng số Ci , thay vào uN ta nhận được nghiệm gần đúng. 12 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 3. Phương pháp Galerkin : Phương pháp Galerkin là phương pháp dư có trọng lấy ψ k  φ k . Như vậy nếu phương pháp Galerkin chọn hệ hàm { φ k} là hệ trực giao thì sẽ rất thuận lợi dẫn đến (1.18) . Ví dụ: Trong ví dụ trên phần 1.2.1 ta co phương trình vi phân của độ võng d4 w EJ q 0 ; ( 0 x  ) dx4 nghiệm gần đúng theo phương pháp Galerkin sẽ được chọn : NN πi wN  C i w i (x)  C i sin x , i 1 i 1  w(0) w( ) 0 như vậy wN thỏa mãn các điều kiện biên : w (0) w ( ) 0 N π4 i 4 πi q Thay vào phương trình vi phân của độ võng ta được :  4 Ci sin x 0 , i 1   EJ πk Nhân cả hai vế với sin x và lấy tích phân hai vế ( k = 1, 2,3, , N) dẫn đến :  N  π4 i 4 πi πk q  πk C sin xsin x dx sin xdx  4 i i 1 0    EJ 0  πk  chú ý rằng hệ sin  là một hệ trực giao trong không gian có tích vô hướng là   tích phân xác định trên [0 ,  ] nên ta luôn có : 0 khi i k  πi πk sin xsin x dx  ;   khi i k 0 2  0 khi k 2m πk  k Bk = sin xdx 1 ( 1) 2 khi k 2m 1 0  kπ (2m 1)π π4 k 4 như vậy ta có ma trận (A) là ma trận chéo với A kk = ; 23 2 q B2m = 0 ; B2m -1 = (2m 1)πEJ 13 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 2 q 2 3 4  4 q C2m = 0 ; C2m - 1 = (2m 1)πEJ π(2m4 1) 4 (2m 1)πEJ 5 5 NM 44 q π(2m 1) vậy nghiệm gần đúng wN  C i w i (x)  5 5 sin x i 1 m 1 (2m 1) π EJ  Độ võng lớn nhất tại x =  / 2 => Nếu lấy M =10 ta được M 44 q π(2m 1) 0,013054  4 q 5  4 q 0,013021  4 q wmax  5 5 sin m 1 (2m 1) π EJ 2 EJ 384EJ EJ sai số 0,25% 4. Phương pháp sai phân : Phương pháp sai phân là phương pháp biểu diễn gần đúng các giá trị của đạo hàm theo các giá trị của hàm số tại các điểm lân cận ( xuất phát từ khai triển Tay-lo của hàm số). Chẳng hạn đối với hàm một biến f(x) ta có được : df f (x Δx) f (x Δx) dxx 2Δx d2 f f (x Δx) f (x Δx) 2f (x) 2 2 dx x Δx d3 f f (x 2Δx) 3f(x Δx) 3f(x) f(x Δx) dx3Δx 3 x d4 f f (x 2Δx) 4f (x Δx) 6f (x) 4f (x Δx) f (x 2Δx) 4 4 dx x Δx khi đó việc giải các phương trình vi phân được đưa về việc tìm các giá trị hàm số tại các điểm nút lưới ( mà khoảng cách giữa các điểm chính là Δx ). ví dụ : xét bài toán uốn dầm như hình 1.2 : dầm chịu tải trọng đều với cường độ q , tựa gối khớp tại hai đầu A , B Có phương trình vi phân của độ võng 0  x d4 w EJ q 0 ; ( 0 x  ) dx4 w(0) w( ) 0 z với điều kiện biên : w (0) w ( ) 0 14 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Chia lưới sai phân như hình vẽ : 4 đoạn với 5 nút (từ 0 đến 4, với nút 0 tại x = 0 ; nút 4  tại x =  ), có Δx 4 2 tại nút 0 : w0 = 0 w 1 2w 0 w 1 0.Δx 0 => w-1 = - w1 pΔx4 tại nút 1: w 4w +6w 4w w 1 0 1 2 3 EJ pΔx4 tại nút 2: w 4w +6w 4w w 0 1 2 3 4 EJ pΔx4 tại nút 3: w 4w +6w 4w w 1 2 3 4 5 EJ 2 tại nút 4: w3 2w 4 w 5 0.Δx 0 w4 = 0 pΔx4 0,013672p4 => w2 =3,5. sai số tương đối : 0.049987 EJ EJ I.3 CÁC NGUYÊN LÝ BIẾN PHÂN I.3.1 Nguyên lý thế năng toàn phần dừng ( nguyên lý biến phân về chuyển vị) Thế năng toàn phần của một hệ đàn hồi  được xác định là :  = U - A trong đó : U - thế năng biến dạng của vật thể đàn hồi được tích lũy trong quá trình biến dạng A – công của ngoại lực sinh ra trên các chuyển dời của ngoại lực do vật thể bị biến dạng. Nguyên lý phát biểu rằng : Trong tất cả các trường chuyển vị ( các trạng thái chuyển vị ) khả dĩ động ( tức là thỏa mãn điều kiện tương thích và điều kiên biên động học) thì trường chuyển vị thực ( tương ứng với sự cân bằng của vật thể) sẽ làm cho thế năng toàn phần  đạt giá trị dừng (nhỏ nhất). Tức là khi đó: δ ({u}) = δ U({u}) - δ A({u}) = 0. Thế năng biến dạng U được tính theo công thức : 1 1 U{ ε}TT {σ}dv hay U {ε} [D]{ε}dv 2VV 2 15 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . nếu có biến dạng đầu thì : 1 TT U{ ε} [D]{ε}dv {ε} [D]{ε0 }dv 2 VV Còn công của ngoại lực ( gồm lực khối {g} , lực mặt {p}) trên chyển dời {u} sẽ là : A {u}TT {g}dv {u} {p}ds VSt khi đó thế năng toàn phần sẽ là : 1 ({u}) U A {ε}TTT [D] {ε} 2{ε } dv {u} {g}dv {u} {p}ds (I.19) 2 0 VVSt Ví dụ : Xét dầm có chiều dài  , tựa khớp hai đầu và chịu tải phân bố đều q, Thế năng biến dạng 0  x 1 M2 (x) d2 w U = dx mà M(x) EJ 2 20 EJ dx 2 z EJ d2 w  => U dx ; A = q.w dx 2 2 0 dx 0 2 EJ d2 w  khi đó thế năng toàn phần :  dx q.w dx . 2 2 0 dx 0 Theo nguyên lý thế năng toàn phần dừng : w(x) là thực nếu : δ δU δA 0 , tức là:  d2 w d 2 w  d 2 w d 2δw d 2 w dδw  dδw d 3 w δU EJ δ dx EJ dx EJ EJ . dx dx2 dx 2 dx 2 dx 2 dx 2dx dx dx 3 0 0 0 0  2 3 4  d w dδw d w d w => δU EJ δw δw. dx  dx2dx dx 3 dx 4 0 0     d4 w d 2 w dδw d 3 w còn δA q.δw dx ; vậy δ EJ. q δw dx EJ δw dx4 dx 2dx dx 3 0 0 0 w(0) w( ) 0 do điều kiện biên : w (0) w ( ) 0 16 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . d2 w  d4 w do đó tại x = 0 và x =  ta có : = 0 ; δw =0 => δ EJ. q δw dx 2 4 dx 0 dx d4 w để thỏa mãn với mọi δw thì cần phải có EJ. q = 0 , đây chính là phương trình dx4 cân bằng của dầm chịu uốn viết cho chuyển vị. I.3.2 Nguyên lý cực tiểu của năng lượng bù toàn phần ( nguyên lý biến phân theo ứng suất) Năng lượng bù toàn phần của vật thể đàn hồi * được định nghĩa là : * = U* - A* , trong đó U* là năng lượng bù của biến dạng. 1 Trong giai đoạn đàn hồi tuyến tính thì U* = U = {σ}T {ε}dv , nếu có biến 2 V *T1 dạng đầu {ε0 } thì U{ σ} [C]{σ} 2{ε0 } dv 2 V * * T A là công bù của ngoại lực : A = {p} {u}ds trong đó st là phần biên đã biết chuyển St 1 vị {u}, vì vậy *TT {σ} [C]{σ} 2{ε } dv {p} {u}ds 2 0 VSt Nguyên lý cực trị năng lượng bù toàn phần được phát biểu như sau: Trong tất cả các trường ứng suất khả dĩ tĩnh ( tức là thỏa mãn các điều kiện cân bằng và điều kiện biên tĩnh học trên St ) thì trường ứng suất thực ( tương ứng thỏa mãn điều kiện tương thích) sẽ làm cho năng lượng bù toàn phần * đạt giá trị dừng. δU δU ({σ}) δA ({σ}) 0 17 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Chương II CƠ SỞ VÀ CÁC NỘI DUNG CƠ BẢN CỦA PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN II.1 Khái niệm về phương pháp phần tử hữu hạn (PTHH) Phương pháp phần tử hữu hạn (PP PTHH) là một phương pháp số để tìm nghiệm gần đúng của một hàm chưa biết trong miền xác định V. Tuy nhiên PP PTHH không tìm dạng xấp xỉ của hàm cần tìm trên toàn miền V mà chỉ trong từng miền con Ve ( phần tử) thuộc miền xác định V. Chính vì lẽ đó nên phương pháp này rất thích hợp để tìm nghiệm gần đúng cho các bài toán vật lý, kỹ thuật khi mà hàm cần tìm được xác định trên những miền phức tạp là những vùng nhỏ có các đặc trưng hình học, vật lý khác nhau, chịu các điều kiện biên khác nhau. Phương pháp được phát biểu một cách tổng quát chặt chẽ như một phương pháp biến phân hay phương pháp dư có trọng số trên mỗi phần tử. Trong PP PTHH , miền V được chia thành một số hữu hạn các miền con được gọi là các phần tử. Các phần tử này được kết nối với nhau tại các điểm trên biên được gọi là các nút. Trong phạm vi mỗi phần tử, đại lượng cần tìm ( chẳng hạn đó là các biến dạng, dịch chuyển, ứng suất , ) được lấy xấp xỉ trong một dạng hàm đơn giản – được gọi là các hàm xấp xỉ ( approximation function). Các hàm xấp xỉ này được được tính thông qua các giá trị của nó ( đôi khi qua các giá trị đạo hàm) tại các điểm nút trên phần tử và các giá trị này được gọi các bậc tự do của phần tử mà ta xem như là các ẩn cần tìm của bài toán. Trong bài toán cơ học vật rắn biến dạng và cơ kết cấu tùy theo ý nghĩa vật lý của các hàm xấp xỉ ta có thể áp dụng bài toán theo ba loại mô hình sau: 1. Mô hình tương thích : Xem chuyển vị là đại lượng cần tìm trước và hàm xấp xỉ biểu diễn gần đúng dạng phân bố của chuyển vị trong phần tử. Các ẩn số được xác định từ hệ phương trình thiết lập trên cơ sở nguyên lý thế năng toàn phần ( hay nguyên lý biến phân Lagrange). 2. Mô hình cân bằng : Hàm xấp xỉ biểu diễn gần đúng dạng phân bố của ứng suất hay nội lực trong phần tử. Các ẩn số được xác định từ hệ phương trình thiết lập 18 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . trên cơ sở nguyên lý năng lượng hệ toàn phần dừng ( hay nguyên lý biến phân về ứng suất – nguyên lý Castigliano) 3. Mô hình hỗn hợp : Xem các đại lượng chuyển vị và ứng suất là hai hai yếu tố độc lập. Các hàm xấp xỉ biểu diễn gần đúng dạng phân bố của chuyển vị và ứng suất trong phần tử. Các ẩn cần tìm được xác định từ hệ phương trình thiết lập trên cơ sở nguyên lý biến phân Reisner. Sau khi tìm được giá trị các ẩn số ( bằng việc giải một hệ phương trình đại số), như vậy ta đã tìm được xấp xỉ các đại lượng cần tìm, từ đó tìm được giá trị của các đại lượng còn lại. Mô hình tương thích được áp dụng rộng rãi. Trong giáo trình này chủ yếu các bài toán được giải theo mô hình tương thích. II.2 Trình tự các bước phân tích bài toán theo phương pháp PTHH  Bước 1. Rời rạc hóa miền khảo sát Miền khảo sát V được chia thành các miền con Ve ( phần tử) có dạng hình học thích hợp. Với bài toán cụ thể thì số phần tử, hình dạng hình học của phần tử và kích thước các phần tử phải được xác định cụ thể. Số điểm nút mỗi phần tử không được lấy tùy tiện mà phải phụ thuộc vào dạng hàm xấp xỉ định chọn. Các phần tử có các dạng hình học đơn giản : Phần tử một chiều : Phần tử bậc 1 Phần tử bậc 2 Phần tử bậc 3 Phần tử hai chiều: Phần tử bậc 1 Phần tử bậc 2 Phần tử bậc 3 19 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Phần tử ba chiều : Dạng tứ diện Phần tử bậc 1 Phần tử bậc 2 Phần tử bậc 3 Dạng lăng trụ Phần tử bậc 1 Phần tử bậc 2 Phần tử bậc 3  Bước 2 Chọn hàm xấp xỉ thích hợp : Chọn dạng hàm xấp xỉ sao cho đơn giản đối với tính toán, nhưng vẫn đảm bảo các tiêu chuẩn hội tụ . Thường chọn các hàm này có dạng đa thức. sau khi chọn dạng hàm xấp xỉ ta biểu diễn các hàm này (kể cả đạo hàm của nó) theo tập hợp các giá trị tại các nút của phần tử {q}e  Bước 3 Xây dựng phương trình phần tử, tức là thiết lập ma trận độ cứng phần tử [K]e và véc tơ tải phần tử {P}e . Kết quả nhận được phương trình có dạng : [K]e {q}e = {P}e  Bước 4 Ghép nối các phần tử trên cơ sở mô hình mà kết quả là hệ thống phương trình : [K] {q} {P} trong đó : [K] là ma trận độ cứng tổng thể ( toàn miền V) {q} là véc tơ tập hợp các giá trị đại lượng cần tìm tại tất cả các nút ( tức là véc tơ chuyển vị nút tổng thể) {P} là véc tơ số hạng tự do tổng thể ( tức là véc tơ tải tổng thể ) sau đó sử dụng điều kiện biên của bài toán sẽ nhận được hệ phương trình : [K ] {q } {P } - là hệ phương trình hệ thống hay còn gọi là hệ phương trình để giải 20 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .  Bước 5 Giải hệ phương trình đại số : [K ] {q } {P } , tìm được chuyển vị của các nút. Việc giải hệ phương trình [K ] {q } {P } đối với bài toán tuyến tính không gặp khó khăn, nhưng với bài toán phi tuyến thì sẽ dùng phương pháp lặp ( mà được tuyến tính hóa , chẳng hạn như phương pháp Newton – Raphson) mà ở mỗi bước lặp ma trận độ cứng [K* ] và {P* } sẽ thay đổi.  Bước 6 Hoàn thiện: tính giá trị các đại lượng còn lại: ứng suất, biến dạng II.3 Hàm xấp xỉ - đa thức xấp xỉ . Phép nội suy 1. Hàm xấp xỉ : Tư tưởng chính của PP PTHH là xấp xỉ hóa đại lượng cần tìm trong mỗi miền con – phần tử Ve . Do đó đầu tiên phải chọn hàm số đơn giản mô tả gần đúng đại lượng cần tìm trong mỗi phần tử. Hàm số đơn giản hay được chọn có dạng đa thức, vì: Đa thức được xem như tổ hợp tuyến tính các đơn thức , các đơn thức này thỏa mãn yêu cầu của Ritz , Galerkin. Hàm xấp xỉ dạng đa thức thường dễ tính toán, dễ thiết lập công thức khi xây dựng các phương trình PP PTHH, dễ đạo hàm, dễ lấy tích phân. Có khả năng tăng độ chính xác ( bằng cách tăng số bậc của đa thức), tuy nhiên trong thực tế thường chỉ lấy bậc thấp mà thôi. 2. Phép nội suy Trong PP PTHH , các hệ số trong các hàm đa thức xấp xỉ được biểu diễn qua các giá trị của nó ( cả những giá trị đạo hàm) tại các điểm nút được định trước trên mỗi phần tử. Ví dụ : u(x) u(x) u(x) x x x a a b b a b a b 2 21 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Nội suy xấp xỉ hằng số Nội suy tuyến tính Nội suy bậc hai a b Nội suy hằng số (u)x u u( ) 0 2 b.u(a) a.u(b) u(b) u(a) Nội suy tuyến tính : (u)x u x 0 b a b a 2 Nội suy bậc hai : (u)x u0 = a1 + a2x + a3 x 3. Dạng đa thức xấp xỉ : Hàm xấp xỉ được chọn dưới dạng đa thức đơn giản : Bài toán 1 –D : u(x) = a1 + a2x ( xấp xỉ tuyến tính) 2 u(x) = a1 + a2x + a3 x ( xấp xỉ bậc 2) 2 3 u(x) = a1 + a2x + a3 x + a4x ( xấp xỉ bậc 3) n 1 i 1 như vậy nếu u(x) xấp xỉ bậc n thì u(x) = ai x hay là: i 1 a1  2 3 n a2 u(x) = [1 x x x x ] .  hay u(x) = [P(x)] . {a}  an 1  trong đó [P(x)] – ma trận các đơn thức {a} - véc tơ tọa độ tổng quát (hay véc tơ các tham số) Bài toán 2-D 2 2 Nếu chọn xấp xỉ bậc hai thì khi đó u(x,y) = a1 + a2x + a3y + a4x + a5y + a6xy a1  2 2 a2 hay u(x,y) = [ 1 x y x y xy] .  => u(x,y) = [P(x,y)] . {a}  a6  Bài toán 3 – D có u(x,y,z) = [P(x,y,z)] . {a} (II.1) Trường hợp xấp xỉ tuyến tính có [P(x,y,z)] = [ 1 x y z ] khi đó a1  u(x, y,z)  1 x y z 0 0 0 0 0 0 0 0 a 2 u v(x,y,z)  00001xyz0000.  e  w(x, y,z) 0 0 0 0 0 0 0 0 1 x y z e a12  22 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 4. Chọn bậc của đa thức xấp xỉ. Các đa thức xấp xỉ cần thỏa mãn được các yêu cầu sau:  Các đa thức xấp xỉ phải thỏa mãn điều kiện hội tụ: Do PP PTHH là một phương pháp số nên phải đảm bảo được rằng khi kích thước các phân tử giảm đi thì kết quả tính phải hội tụ đến giá trị chính xác. Để có được điều này thì các đa thức xấp xỉ ue phải thỏa mãn 3 điều kiện sau: - Liên tục trong phần tử Ve . Điều này được thỏa mãn vì xấp xỉ là đa thức. - Bảo đảm trong phần tử có trạng thái đơn vị ( hằng số) và có các đạo hàm riêng của nó đến bậc cao nhất mà phiếm hàm I(u) đòi hỏi. I(u) F(x,u,u ',u", ,u(m) ) dv - x là tọa độ điểm, u là hàm xấp xỉ V - Trên biên phần tử , u và các đạo hàm của nó có đến cấp m -1 liên tục. Chẳng hạn như khi u là chuyển vị thì muốn đảm bảo trạng thái đơn vị và dịch chuyển cứng thì trong đa thức xấp xỉ không được bỏ qua số hạng a1 , hay không được bỏ qua thành phần 1 trong [P(x,y,z)]. Khi làm mịn lưới các phần tử cần tuân theo các quy tắc sau : + Lưới sau mịn hơn trên cơ sở lưới trước, các điểm nút lưới trước cũng có trong tập các nút lưới sau. + Các phần tử có kích thước nhỏ hơn trước, nhưng dạng hình học vẫn phải được giữ như trước. + Dạng đa thức không đổi trong quá trình mịn hóa lưới phần tử.  Các đa thức xấp xỉ được chọn sao cho không làm mất tính đẳng hướng hình học. Dạng đa thức được chọn từ tam giác Passcal ( cho bài toán 2 chiều), tháp Passcal cho bài toán 3 chiều. 23 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Tam giác Passcal Bậc đa thức Số tham số 1 Hằng số 1 x y tuyến tính 3 2 2 x xy y bậc 2 6 x3 x2y xy2 y3 bậc 3 10  Số các phần tử của {a} – là các tham số của đa thức xấp xỉ phải bằng số bậc tự do của phần tử {q}e , khi đó có thể nội suy đa thức xấp xỉ thep giá trị đại lượng cần tìm tại nút phần tử. 5. Biểu diễn đa thức xấp xỉ theo véc tơ các bậc tự do của phần tử. Ma trận hàm dạng Bậc tự do của một nút ( nodal degree of freedom ) là giá trị ( có thể có cả giá trị đạo hàm ) của hàm ( đa thức ) xấp xỉ tại nút. Tập hợp tất cả các bậc tự do của các nút trên phần tử được gọi là véc tơ bậc tự do của phần tử, ký hiệu là { q}e - còn được gọi là véc tơ chuyển vị nút phần tử . Các bậc tự do này là ẩn số của bài toán PP PTHH Ví dụ : Trong bài toán phẳng đàn hồi , khi dùng phần tử tam giác có 3 điểm nút, mỗi nút có 2 bậc tự do : đó là các chuyển vị theo phương x và y , do đó tập hợp chuyển vị ở 3 nút là véc tơ chuyển vị của phần tử : T {q}e = { ui , vi , uj , vj , uk ,vk }e T  { q1 , q2 , q3 , q4 , q5 , q6}e Như vậy nếu phần tử e có r nút , mỗi nút có s bậc tự do thì véc tơ chuyển vị phần tử {q}e có số thành phần ne = s × r Các đa thức xấp xỉ được biểu diễn theo các bậc tự do của phần tử {q}e , và phải thỏa mãn điều kiện : Các giá trị của đa thức xấp xỉ ( có thể cả đạo hàm của nó) tại các điểm nút thuộc phần tử phải đồng nhất bằng giá trị các bậc tự do của phần tử 24 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Tức là : u(nút1)  u(nút1)  u(x1 , y 1 ,z 1 )  [P(x 1 , y 1 ,z 1 )] u(nút 2) u(nút 2) u(x,y,z) [P(x,y,z)] 2 2 2 2 2 2   {q}e mà   . {a}     u(nút r)  u(nút r)  u(xr , y r ,z r )  [P(x r , y r ,z r )] [P(x1 , y 1 ,z 1 )] [P(x , y ,z )] do vậy , ta sẽ có : 2 2 2 . {a} {q} hay là [A]. {a} {q} (II.2)  e e [P(xr , y r ,z r )] Chú ý rằng ma trận [A] là ma trận vuông có ne × ne phần tử và chỉ chứa tọa độ các điểm -1 nút của phần tử => {a} = A .{q}e . Sau khi tính được {a} ta thay vào (II.1) tức là : u(x,y,z) = [P(x,y,z)].{a} -1 => u(x,y,z) = [P(x,y,z)]. A .{q}e hay u(x,y,z) = [N].{q}e (II.3) với [N] = [P(x,y,z)]. A-1 (II.4) được gọi là ma trận các hàm nội suy hay là ma trận hàm dạng Ví dụ 1: Tìm ma trận hàm dạng của phần tử thanh lăng trụ chịu kéo, nén dọc trục ( Biết thanh chiều dài L , diện tích mặt cắt ngang : F , mô đun đàn hồi E) du Giải Tại mọi điểm chỉ tồn tại chuyển vị và biến dạng dọc trục : u(x) và ε x dx do đó phiếm hàm thế năng đàn hồi : 1T F F U{ ε} σdv εx σ x dx ε x E.ε x dx 2VLL 2 2 EF 2 EF du 2 => U( εx ) dx = ( ) dx 2 L 2L dx thấy rằng đa thức xấp xỉ u(x) chỉ đòi hỏi tồn tại đạo hàm bậc nhất, do đó có thể lấy u(x) là hàm xấp xỉ tuyến tính : u(x) = a1 + a2x ( 0 ≤ x ≤ L ) hay : a1 a1 u(x) [1 x ] . = [P(x)] . {a} => [P(x)] = [ 1 x ] , còn {a} = a 2 a2 Như vậy {a} có 2 tham số, nên véc tơ chuyển vị nút {q}e của phần tử cũng chỉ có 2 bậc tự do - là 2 chuyển vị dọc trục x tại hai đầu nút : 25 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . T T {q}e = {q1 q2}e ≡ {u1 u2 }e u(nút1)  u(x0)  a1  10 a 1 q 1  ta có :     u(nút2)  u(xL)  a1 a 2 L  1L a 2 q 2  1 0 P(x ) 1 0  1  -1 vậy [A] = => A = 1 1 => P(x2 ) 1 L LL 1 0 -1 x x hàm dạng : [N]e = [P(x)]. [A] = [1 x ] . 1 1 1 [N1 (x) N 2 (x)] LL LL x x [N]e = 1 - hàm dạng của phần tử thanh chịu kéo , nén dọc trục LL Thay vào (II.3) ta được biểu diễn xấp xỉ chuyển vị u(x) theo các chuyển vị tại các nút: x x u(x) = [N]e . {q}e = 1 q1 + q2 L L Các hàm ui(x) còn được gọi là các hàm nội suy Lagrange bậc 1. Ví dụ 2. Chọn đa thức xấp xỉ và tìm ma trận hàm dạng phần tử dầm 2 điểm nút chịu uốn. Giải Bỏ qua biến dạng dọc trục, khi đó trạng thái chuyển vị của dầm được đặc trưng bởi chuyển vị v(x) theo phương vuông góc với trục dầm ( trục 0x) Do sự tương thích biến dạng, nên tại biên giữa các phần tử cần đảm bảo sự liên tục về dv chuyển vị v và góc xoay θ => số bậc của phần tử là 4 : đó là các chuyển vị và dx góc xoay tại hai đầu: T {q}e = {v1 θ 1 v2 θ 2 }e T ≡ {q1 q2 q3 q4 }e (a) Như vậy để nội suy hàm xấp xỉ v(x) thì véc tơ tham số {q}e cũng phải có 4 tham số => đa thức xấp xỉ của v(x) là đa thức bậc 3 ( vì đa thức bậc 3 có 4 hệ số) 26 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . a1 2 3 2 3 a2 v(x) = a1 + a2x + a3x + a4x = [ 1 x x x ] . = [Pv(x)]. {a} a3 a4 a1 dv 2 2 a2 Góc xoay θ (x) = = a2 + 2a3x + 3a4x = [ 0 1 2x 3x ] . = [Pө(x)]. {a} dx a3 a4 Từ (a) ta có : hay ở dạng ma trận : q1 ≡ v1 = vx 0 a1 1 0 0 0 dv q1  a 1  q2 ≡ θ 1 = a2 0 1 0 0 dx q2 a 2 x 0  2 3  2 3 q3 1 L L L a 3 q3 ≡ v2 = v = a1 + a2L + a3L + a4L x L q 2 a 4  0 1 2L 3L 4  dv 2 q4 ≡ θ 2 = a2 2a 3 L 3a 4 L  {q}e = [A]. {a} dx x L Chú ý rằng u(nút k ) và P(x,y,z) trong (II.2) ở đây sẽ là : v1 v2 u(nút 1) = ; u(nút 2) = θ1 θ2 P (x) 1 x x2 x 3 => [P(x,y,z)] = v P (x) 2 θ 0 1 2x 3x từ đó ta cũng có thể thấy ngay : [P (0)] 1 0 0 0 1 0 0 0 v [P (0)] 0 1 0 0 0 1 0 0 [A] θ => A-1 = 2 3 2 2 3 [Pv (L)] 1 L L L 3/L 2/L 3/L 1/L [P (L)] 2 3 2 3 2 θ 0 1 2L 3L 2/L 1/L 2/L 1/L 27 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . v -1 Do đó ta có ma trận hàm dạng của chuyển vị v(x) sẽ là : [N ] = [Pv(x)] . [A ], tức là : 1 0 0 0 0 1 0 0 v 2 3 v v v v [N (x)] = [1 x x x ] . = [N 1 N 2 N 3 N 4 ] 3/L2 2/L 3/L 2 1/L 3 3 2 3 2 2/L 1/L 2/L 1/L với : 3x2 2x 3 2x 2 x 3 Nv 1 ; N v x 12 3 2 L 2 LLL 3x2 2x 3 x 2 x 3 N;Nv v 3LLL2 3 4 L 2 Khi đó hàm chuyển vị v(x) của phần tử dầm chịu uốn là : 4 v v v(x) = [N (x)] . {q}e =  Ni .q i i 1 ө -1 Còn ma trận hàm dạng của góc xoay Ө(x) sẽ là : [N ] = [Pө(x)] . [A ] 1 0 0 0 0 1 0 0 ө 2 ө ө ө ө [N ] = [0 1 2x 3x ] = [N 1 N 2 N 3 N 4 ] 3/L2 2/L 3/L 2 1/L 3 3 2 3 2 2/L 1/L 2/L 1/L II.4 Các phương trình cơ bản 1. Chuyển vị , biến dạng và ứng suất trong phần tử. Ma trận cứng phần tử , véc tơ tải phần tử. Do ở đây ta giải bài toán theo mô hình tương thích ( phương pháp chuyển vị) , nên đại lượng cơ bản cần xác định đầu tiên là chuyển vị. Chuyển vị được xấp xỉ hóa và nội suy theo véc tơ chuyển vị nút các phần tử {q}e.  Sau khi tìm được ma trận các hàm dạng [N], ta sẽ biểu diễn được trường chuyển vị theo các bậc tự do ở các nút phần tử {q}e : {u}e = [N]. {q}e (II.5) 28 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .  Ta có các phương trình Cô-si cho sự liên hệ giữa trạng thái biến dạng và chuyển vị, do đó trạng thái biến dạng của các điểm trong phần tử sẽ là : εe   u e   [N]{q}e [B] {q} e (II.6) trong đó [B] = [N] (II.7)  0 0 x 0 0 y 0 0  z Chú ý :  , [N] – hàm dạng của chuyển vị   0 y  x 0   z  y 0  z  x và [B] được gọi là ma trận tính biến dạng  Trạng thái ứng suất tại một điểm thuộc phần tử trong trường hợp vật liệu tuân theo Hooke sẽ là : σ [D]. ε ε0 σ 0 (II.8) e e e  e với ε0 , σ0 là biến dạng và ứng suất ban đầu của phần tử. e e Sử dụng (II.6) ta được : σ [D].[B] q D ε0 σ 0 e e e  e hay là : σ [T] q D ε0 σ 0 (II.9) e e e  e trong đó : [T] = [D]. [B] là ma trận tính ứng suất (II.10) Như vậy nhờ có (II.5) , (II.6) , (II.9) ta có thể biểu diễn chuyển vị, biến dạng, ứng suất trong phần tử theo véc tơ chuyển vị nút phần tử {q}e Thế năng toàn phần của phần tử Ve theo (I.19) là : 1 ({u}) U A {ε}TTT σ dv {u} {g} dv {u} {p} ds (II.11) e e e2 ee e e e e VVSe e e Sử dụng (II.5) , (II.6) và (II.9) vào (II.11) ta được : 29 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 1 ({q} ) {q}TT [B] [D][B] q dv e e2 e e Ve T T1 T T 0 1 T T 0 {u}e {g} e dv {u} e {p} e ds {q} e [B] [D] ε dv {q} e [B] σ dv 2 e 2 e VSVVe e e e T T T Chú ý {u}e = [N]. {q}e => {u}e = {q}e . [N] , thay vào ta được : 1 ({q} ) {q}TT [K] q {q} {P} (II.12) e e2 e ee e e trong đó :[K]e - Ma trận cứng phần tử T [K]e [B] [D][B]dv (II.13) V và {P}e - véc tơ tải phần tử : T T1 T 0 1 T 0 {P}e {N} {g} e dv {N} {p} e ds [B] [D] ε dv [B] σ dv (II.14) 2 e 2 e VSVVe e e e Nhận thấy rằng do [D] là ma trận đối xứng nên tích [B]T [D] . [B] là ma trận đối xứng do đó [K]e là ma trận đối xứng. Ví dụ 1: Tìm ma trận cứng và véc tơ tải của phần tử thanh lăng trụ chịu kéo, nén dọc trục ( thanh chiều dài L , diện tích mặt cắt ngang : F , mô đun đàn hồi E) du Giải Tại mọi điểm chỉ tồn tại chuyển vị và biến dạng dọc trục : u(x) và ε x dx thấy rằng đa thức xấp xỉ u(x) chỉ đòi hỏi tồn tại đạo hàm bậc nhất, do đó có thể lấy u(x) là hàm xấp xỉ tuyến tính : u(x) = a1 + a2x ( 0 ≤ x ≤ L ) hay : a1 u(x) [1 x ] . = [P(x)] . {a} , theo kết quả của ví dụ 1 phần (II.2) ta có : a 2 biểu diễn xấp xỉ chuyển vị u(x) theo các chuyển vị tại các nút: 30 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . x x u(x) = [N]e . {q}e = 1 q1 + q2 L L x x như vậy hàm dạng của chuyển vị là [N]e = 1 (II.15) LL ứng suất chỉ còn một thành phần σx E.ε x [D].{ε} => [D] = E d x x d x d x Ma trận tính biến dạng : [B] = [  ] [N] = 1 = 1 dx LL dx L dx L 1 1 1 [B] =  1 1 , do vậy ma trận cứng phần tử theo (II.13): LLL L 1 1 1 E.F 1 1 [K] [B]T [D][B]dv E [ 1 1] .F . dx = e LLL 1 1 1 Ve 0 E.F 1 1 ma trận cứng của phần tử thanh chịu kéo nén dọc trục: [K]e = (II.16) L 1 1 ( E – mô đul đàn hồi , F – thiết diện thanh , L – chiều dài thanh ) Véc tơ tải phần tử {P}e được tính theo (II.14) là : * Do lực phân bố dọc trục : x LL 1 P TT L {P}e = [N] {p} ds [N] {p(x)} dx p(x)dx , nếu p(x) = p0 hằng số e e x Se 0 0 L P p0 L 1 thì {P}e = (II.17) 2 1 * Do nhiệt độ thì : L 0 i 1T 0 1 1 1 0 E α.T 1 {P}e = [B] [D] ε dv E α.T dx (II.18) 2 e 2 L1 2 1 Ve 0 T0 độ biến thiên nhiệt độ 2. Ghép nối các phần tử - Ma trận cứng – Véc tơ tải tổng thể. a) Ghép nối các phần tử: 31 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Vật thể V được chia thành Ne miền con hay là phần tử ( Ve) bởi R điểm nút. Nếu mỗi nút có s bậc tự do thì số bậc tự do của cả hệ sẽ là : n = R × s . Gọi q là véc tơ chuyển vị nút tổng thể ( còn gọi là véc tơ chuyển vị nút kết cấu) , q có n thành phần và là tập hợp tất cả các bậc tự do của tất cả các nút của hệ. Nếu mỗi phần tử có r nút khi đó số bậc tự do của mỗi phần tử sẽ là : ne = r × s. Véc tơ chuyển vị nút phần tử {q}e có ne thành phần gồm tất cả các bậc tự do của r nút của phần tử. Ta có các thành phần của véc tơ {q}e sẽ nằm trong các thành phần của véc tơ q, do vậy ta có : {q}e = [L]e . q (II.19) (ne × 1) (ne × n) (n× 1) trong đó [L]e được gọi là ma trận định vị của phần tử có kích thước (ne × n) Ví dụ 1: Dầm với 4 nút , mỗi nút có 2 bậc tự do như hình vẽ : Véc tơ chuyển vị nút tổng thể q sẽ là : T q = { q1 q2 q3 q4 q5 q6 q7 q8} khi đó ta sẽ có : q1  1 0 0 0 0 0 0 0 q1  q2 0 1 0 0 0 0 0 0 q2 {q}  [L] .{q}  1 1 q3 0 0 1 0 0 0 0 0  q4  0 0 0 1 0 0 0 0 q8  q3  0 0 1 0 0 0 0 0 q 1  q4 0 0 0 1 0 0 0 0 q 2 {q}  [L] .{q}  2 2 q5 0 0 0 0 1 0 0 0  q6  0 0 0 0 0 1 0 0 q 8  q5  0 0 0 0 1 0 0 0 q 1  q6 0 0 0 0 0 1 0 0 q 2 {q}  [L] .{q}  3 3 q7 0 0 0 0 0 0 1 0  q8  0 0 0 0 0 0 0 1 q 8  32 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . b) Ma trận cứng - véc tơ tải tổng thể. Sử dụng (II.12) ta sẽ có thế năng toàn phần : NNe e 1 TT   e({q} e )  {q} e [K] e qe {q} e {P} e e 1 e 1 2 Ne 1 TTTT => (II.20)   {q} [L]e [K] e Le {q} {q} [L] e {P} e e 1 2 Áp dụng nguyên lý Lagrange ( nguyên lý thế năng toàn phần dừng) về điều kiện cân bằng của toàn hệ tại các điểm nút :   0 q 1  0   δ 0 q2 hay dạng ma trận : 0  q    0 qn NNe e   TT Ta có [L]e [K] e Le {q}  [L] e {P} e 0  q e 1 e 1 dẫn đến phương trình K q P (II.21) N e với : K là ma trận cứng tổng thể K [L]T [K] L (II.22)  e e e e 1 Ne T - và P là véc tơ tải tổng thể P  [L]e {P} e (II.23) e 1 Nhận xét 1. Hệ phương trình K q P {0} biểu diễn sự cân bằng của vật thể tại các điểm nút. Phần tử Kij của ma trận cứng tổng thể K biểu thị lực sinh ra ở nút i do chuyển dịch đơn vị ở nút j khi tất cả các nút bị gắn cứng. 33 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Thành phần Pi của véc tơ tải tổng thể P là ngoại lực tác động lên các phần tử (tính cả biến dạng và ứng suất ban đầu) được quy đổi về tương ứng với bậc tự do thứ i. Trường hợp hệ thanh còn phải kể thêm vào P các ngoại lực tập trung tại tác động lên các điểm nút theo các bậc tự do tương ứng. 2. Khi thiết lập phương trình (II.21) ta chưa đưa vào các điều kiện biên động học ( liên quan đến dịch chuyển), vật thể lúc này là tự do, do đó ma trận K là suy biến ( tức là 1 định thức của K bằng 0) do vậy không tồn tại K . Nhưng sau khi đưa vào các điều kiên biên động học ta sẽ dẫn về phương trình : K q P  . 3. Việc sử dụng ma trận định vị [L]e để tính ma trận cứng K và P thực chất là sắp xếp các phần tử K và P vào vị trí của nó trong K và P . Tuy nhiên trong e e  khi thực hành người ta sử dụng ma trận chỉ số tiện lợi hơn trong quá trình ghép nối. 4. * Do K e đối xứng nên K đối xứng. * K có dạng băng * Bề rộng của băng tùy thuộc theo cách đánh thứ tự số nút Ví dụ Xét khung phẳng 10 tầng ( như hình vẽ), 4 nhịp. Nếu không kể đến 5 nút ngàm cứng thì hệ có 50 nút, mỗi nút có chuyển vị ( 3 bậc tự do) chưa biết => số ẩn của hệ là : 2 50 × 3 = 150 => ma trận K có số phần tử = 150 = 22500 thành phần. a) Cách đánh nút theo chiều ngắn : B = 18 b) Cách đánh nút theo chiều dài : B = 33 34 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Nhưng do ý nghĩa cơ học thấy rằng : Kij 0 khi i và j là các bậc tự do thuộc cùng một nút hoặc hai nút kề nhau, còn lại là bằng 0, chẳng hạn : K1,6 0;K 1,16 0;K 1,18 0;K 1,7 K 1,19 0 => chiều rộng của băng B = (S+1) × s , với S sự sai khác lớn nhất của mã số nút kề nhau. Với cách đánh số nút như hình a) có : B = ( 5 + 1). 3 = 18 b) có B = (10 +1) .3 = 33 Như vậy với ví dụ trên cách đánh số nút theo chiều ngắn như hình a) cho bề rộng băng B là nhỏ nhất. Như thể sẽ có lợi cho việc lưu trữ số liệu và việc giải hệ phương trình đại số . 3. Phép chuyển trục tọa độ Trong các phần trên ta đều tính giá trị các đại lượng : chuyển vị, biến dạng , ứng suất, ma trận hàm dạng [N]e , ma trận cứng phần tử [K]e , véc tơ tải phần tử {P}e đều trong hệ tọa độ thích hợp cho mỗi phần tử - là hệ tọa độ địa phương ( local coordinate system), do vậy phương của các bậc tự do cũng lấy theo hệ tọa độ này. Thực tế các kết cấu mà các phần tử khác nhau thì có các hệ tọa độ địa phương khác nhau, do đó các bậc tự do của các phần tử cũng khác nhau về phương. Chính vì lẽ đó cần có một hệ tọa độ chung cho toàn hệ - gọi là hệ tọa độ tổng thể ( global coordinate system). Nếu gọi : * hệ tọa độ địa phương là xyz và hệ tọa độ tổng thể là x’y’z’ * {q}e , {P}e , [K]e là véc tơ chuyển vị, véc tơ tải , ma trận cứng phần tử trong hệ tọa độ địa phương xyz * {q’}e , {P’}e , [K’]e là véc tơ chuyển vị, véc tơ tải , ma trận cứng phần tử trong hệ tọa độ tổng thể x’y’z’ khi đó ta có mối liên hệ giữa chúng : {q}e = [T]e . {q’}e và {P}e = [T]e {P’}e (II.24) 35 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . ma trận {T}e là ma trận biến đổi ( transformation matrix) các thành phần chuyển vị nút từ hệ tọa độ tổng thể x’y’z’ về hệ tọa độ địa phương xyz Chú ý : T Luôn có [T]e .[T]e = [I] ma trận đơn vị, ( trường hợp nếu [T]e vuông thì nó sẽ là ma trận trực giao, khi đó [T]T = [T]-1) Số thành phần của {q’}e có thể khác số thành phần của {q}e 1 Do đó thế năng toàn phần của phần tử  sẽ là :  {q}TT [K] q {q} {P} e e2 e ee e e 1 =>  {q'}TTTT [T] [K] [T] q' {q'} [T] .[T] {P'} e2 e e e ee e e e e 1 =>  {q'}TTT [T] [K] [T] q' {q'} {P'} , e2 e e e ee e e 1 hay là  {q'}TT [K'] q' {q'} {P'} . với [K '] [T]T [K] [T] (II.25) e2 e ee e e e e e e Ne Do   e , áp dụng nguyên lý Lagrange => : [K'] q ' {P '} (II.26) e 1 Ne Trong đó [K']  [K']e - ma trận cứng tổng thể trong hệ tọa độ tổng thể. e 1 q ' - véc tơ chuyển vị nút tổng thể trong hệ tọa độ tổng thể. Ne {P'} = [P']e {P } n véc tơ tải tổng thể trong hệ tọa độ tổng thể. e 1 với {P }n - véc tơ tải tập trung tại các nút tác dụng theo các phương tương ứng của các thành phần véc tơ chuyển vị nút kết cấu q ' - gọi là véc tơ tải trọng nút ( nodal load vector). Ví dụ: Xét phần tử thanh chịu kéo nén dọc trục : 36 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . nếu lấy hệ tọa độ địa phương 0x dọc trục, 0y vuông góc trục thanh, 0x lập với 0’x’ góc α Gọi các véc tơ cơ sở chính tắc của hệ 0xy là E = {e1 , e2} còn cơ sở chính tắc của 0’x’y’ là E’= { e’1 , e’2 } khi đó ma trận chuyển cơ sở từ E’ sang E là P với cosα sin α P = ( ma trận P có các cột là tọa độ của E biểu diễn qua E’) sin α cosα do đó P-1 là ma trận chuyển cơ sở từ E sang E’, do đó với véc tơ x bất kỳ ta sẽ được: 1 1 cosα sinα [x]EE' P [x] , có P , chuyển vị q1 trong hệ tọa độ 0xy có sin α cosα q1 cosα sin α q1 tọa độ là [q1 ] E 0 sinα cosα q 2  q1 = q’1 . cosα + q’2 . sinα , tương tự chuyển vị q2 biểu diễn qua chuyển vị của hệ tọa độ tổng thể là  q2 = q’3 . cosα + q’4 . sinα Do vậy ta có : q1 qcosα sin α 0 0 q 1 2 q2 0 0 cosα sinα q 3 q 4 do đó : cosα sin α 0 0 [T]e = (II.27) 0 0 cosα sin α Nếu gọi  = cosα , m = sin α khi đó ta có 37 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN .  m 0 0 [T]e = (II.27)’ 0 0 m Theo (II.16) ta có ma trận cứng trong hệ tọa độ địa phương là : E.F 1 1 [K]e = L 1 1 Áp dụng (II.25) thì ma trận cứng của phần tử trong hệ tọa độ tổng thể là  0 m 0 1 1 m 0 0 T E.F [K ']e [T] e [K] e [T] e = L 0 1 1 0 0  m 0 m 2 m  2  m E.F m m2  m m 2 => [K’]e = (II.28) 2 2 L  m   m 2 2 m m  m m 4. Xây dựng ma trận cứng tổng thể - Véc tơ tải tổng thể. Sử dụng ở đây hai hệ thống chỉ số để đánh số cho các bậc tự do của các nút : 1. Hệ thống chỉ số tổng thể: Đánh số các bậc tự do của toàn hệ thống tức là thứ tự của các bậc tự do đang xét trong q hoặc q  , nó được đánh thứ tự là 1 , 2 , 3 , , n = R× s 2. Hệ thống chỉ số phần tử : Để chỉ thứ tự các bậc tự do trong phần tử , hay là thứ tự các bậc tự do trong {q}e hoặc {q’}e , nó được đánh số 1, 2 , 3 , , ne = r × s trong đó R – số nút của cả hệ , r – số nút của phần tử , s – bậc tự do của một nút. Để xác định tương ứng mỗi thành phần của {q}e trong q ( hoặc của {q’}e trong q  ) người ta đưa ra khái niệm ma trận chỉ số [b] ( còn gọi là ma trận liên hệ Boolean) mà mỗi giá trị của thành phần bij chính là chỉ số bậc tự do trong hệ thống chỉ số tổng thể tương ứng với bậc tự do thứ j trong hệ thống chỉ số phần tử của phần tử thứ i. 38 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Ma trận chỉ số [b] có số hàng bằng số phần tử của hệ : Ne , có số cột bằng số bậc tự do của một phần tử. Khi sử dụng ma trận chỉ số [b] để xây dựng ma trận cứng tổng thể [K] và véc tơ tổng thể P ( hoặc [K ] và P ) thì cần nhớ rằng : e mỗi thành phần Kij của ma trận cứng phần tử [K]e ( của phần tử thứ e ) sẽ được cộng thêm vào thành phần Kmn của ma trận cứng tổng thể [K] , Điều này cũng thực hiện cho : e mỗi thành phần K ij của ma trận cứng phần tử [K’]e ( của phần tử thứ e ) sẽ được cộng thêm vào thành phần K mn của ma trận cứng tổng thể [K ], với m = bei ; n = bej ( chú ý bei , bej là giá trị của thành phần hàng e cột i và cột j trong ma trận [b] 5. Áp đặt điều kiện biên Hệ phương trình tổng thể (II.26) : [K'] q ' {P '} sắp xếp lại dạng như sau: K' K ' q '  b  1 p ' 11 12 1 b   ( II.29) K ' K' q ' p ' 21 22 2  2  b q ' trong đó : 2 là véc tơ chứa tất cả các bậc tự do ( chuyển vị nút) đã biết q ' q ' 1 véc tơ chứa các bậc tự do chưa biết ( còn lại trong ) b p ' 1 véc tơ tải với các thành phần đã biết p ' p ' 2 véc tơ tải chưa biết ( còn lại trong ) khi đó hệ phương trình (II.29) được viết thành hai hệ như sau : K' q' K' q'b p' b (a) 1  2  1 11 12 K' q' K' q'b p' (b) 21 1 22  2  2 b b từ hệ (a) ta nhận được : K' q ' p ' K' q ' (II.30) 11 1  1 12  2 q ' Giải hệ (II.30) ta tìm được 1 => véc tơ chuyển vị tổng thể được xác định 39 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . q ' p ' Thay 1 vừa tìm được vào (b) ta tính được 2 => véc tơ tải tổng thể được xác định b Chú ý : Trường hợp các chuyển vị nút bằng 0 thì (II.30) có dạng : K' q ' p ' hệ 11 1  1 này nhận được từ hệ [K'] q ' {P '} bằng cách bỏ đi các hàng , các cột ứng với các bậc tự do trong q ' bằng 0. Ví dụ: Xét hệ thanh phẳng cùng thiết diện F và mô đul đàn hồi E , với các kích thước và chịu lực tập trung P như hình vẽ. Hãy thiết lập ma trận cứng tổng thể và véc tơ tải tổng thể trong hệ tọa độ tổng thể. Giải Mô hình hóa hệ thanh bởi 4 phần tử hữu hạn, mỗi phần tử có 2 bậc tự do ( là 2 chuyển vị theo các phương trục 0x , 0y . Áp dụng công thức (II.28) cho từng phân tử: 3 1 Ở phần tử (1) : α = 300 =>  = cosα = , m = sin α = 2 2 1 2 3 4 <= (Vị trí cột trong ma trận cứng tổng thể) 40 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 3 3 3 3 1 1 E.F 3 1 3 1 2 => [K ']e  = cosα = , m = sinα = 2 2 3 4 5 6 3 3 3 3 3 2 E.F 3 1 3 1 4 [K ']e 4.L 3 3 3 3 5 6 3 1 3 1 3 1 3 1 Ở phần tử (3) : α = 150 =>  = cosα = , m = sin α = 2 2 2 2 5 6 7 8 2 3 1 2 3 1 5 2 2 2 2 1 2 3 1 2 3 3 E.F 2 2 2 2 6 [K ']e 4.L 2 3 1 2 3 1 2 2 2 2 7 1 2 3 1 2 3 2 2 2 2 8 1 3 Ở phần tử (4) : α = 600 =>  = cosα = , m = sin α = 2 2 7 8 3 4 1 3 1 3 7 4 E.F 3 3 3 3 8 [K ']e 4.L 1 3 1 3 3 4 3 3 3 3 Xây dựng ma trận chỉ số [b] : 41 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Chỉ số cục bộ Nút i Nút j Phần tử 1 2 3 4 (1) 1 2 3 4 (2) 3 4 5 6 (3) 5 6 7 8 (4) 7 8 3 4 1 2 3 4 5 6 7 8 3 3 3 3 0 0 0 0 1 3 1 0 0 0 0 6 0 3 3 1 3 2 3 1 3 3 2 3 1 2 3 1 E.F 3 3 K' 2 2 2 2 4.L 2 3 1 2 3 Đối xứng 1 2 2 2 2 3 1 1 3 2 2 2 3 3 2 Véc tơ tải tổng thể chỉ do tải trọng nút tạo nên. Thay tại nút 1 và nút 3 bằng các 1 1 3 3 phản lực liên kết : px , py , px , py , và do tại các nút 1 và nút 3 : chuyển vị bằng 0 ta được : 42 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 1 px 0 1 p1 0 2 y 0 q3 3 0 q4 4 véc tơ tải tổng thể P' , véc tơ chuyển vị nút tổng thể : q ' p3 0 5 x 3 0 6 py q 7 0 7 q 8 P 8 => K' . q ' P' Do véc tơ chuyển vị nút tổng thể có các thành phần : 1 , 2 , 5 , 6 bằng 0, nên trong ma trận cứng tổng thể ta bỏ đi các hàng và cột thứ 1 , 2 , 5 , 6 và cũng bỏ đi các thành phần 1 , 2 , 5 , 6 trong véc tơ tải tổng thể. Khi đó sẽ dẫn về : 6 0 1 3 q 0 q 0 2 3 3 3 3 0,17 q 0 q E.F 2 3 1 4 4 4LP 0,89 1 3 1 3 = => = 0,12 4.L 2 2 q7 0 q7 EF 0,52 1 2 3 q8 P q8 3 3 3 3 2 2 Sử dụng ma trận chỉ số [b] và q ' ta xác định được các véc tơ chuyển vị nút của các phần tử : 0 0,17 0 0,12 EF0 EF 0,89 EF 0 EF 0,52 {q'} ;{q'} ;{q'} ;{q'} 14.L 0,17 2 4.L 0 3 4.L 0,12 4 4.L 0,17 0,89 0 0,52 0,89 43 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Chương III TÍNH TOÁN HỆ THANH BẰNG PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN III.1 Hệ thanh dàn Dàn là một hệ gồm các thanh chỉ chịu kéo nén đúng tâm ( tức là chỉ chịu biến dạng dọc trục) 1. Phần tử thanh chịu biến dạng dọc trục (bài toán một chiều) : Xét phần tử thanh có hai điểm nút chịu biến dạng dọc trục, chịu tải trọng phân bố dọc trụ p(x). Như vậy phần tử có 2 bậc tự do là 2 chuyển vị của nút tại hai đầu . Do có 2 bậc tự do nên chuyển vị dọc trục u(x) của phần tử chỉ có thể là xấp xỉ tuyến tính: x x u(x) = a1 + a2x , hay là u(x) = [N]e . {q}e = 1 q1 + q2 , L độ dài thanh. L L x x như vậy hàm dạng của chuyển vị là [N]e = 1 (III.1) LL ứng suất chỉ còn một thành phần σx E.ε x [D].{ε} => [D] = E d x x d x d x Ma trận tính biến dạng : [B] = [  ] [N] = 1 = 1 dx LL dx L dx L 1 1 1 [B] =  1 1 (III.2) LLL ma trận cứng phần tử theo (II.13): L 1 1 1 E.F 1 1 [K] [B]T [D][B]dv E [ 1 1] .F . dx = e LLL 1 1 1 Ve 0 E.F 1 1 ma trận cứng của phần tử thanh chịu kéo nén dọc trục: [K]e = (III.3) L 1 1 ( E – mô đul đàn hồi , F – thiết diện thanh , L – chiều dài thanh ) Véc tơ tải phần tử {P}e được tính theo (II.14) là : 44 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . * Nếu do lực phân bố dọc trục : x LL 1 P TT L {P}e = [N] {p} ds [N] {p(x)} dx p(x)dx , e e x Se 0 0 L P p0 L 1 nếu p(x) = p0 hằng số thì {P}e = (III.4) 2 1 * Nếu do nhiệt độ thì : L 0 i 1T 0 1 1 1 0 E α.T 1 {P}e = [B] [D] ε dv E α.T dx (III.5) 2 e 2 L1 2 1 Ve 0 T0 độ biến thiên nhiệt độ Ví dụ Giải bài toán thanh dưới đây theo PPPTHH với sơ đồ 2 phần tử. Biết chiều dài thanh là 2a. Diện tích thiết diện F và mô đul đàn hồi E không đổi. Thanh chịu tải phân bố đều dọc trục và có cường độ p0 , chịu tại tập trung đầu dưới là Pt . Các bước giải: 1. Thực hiện việc rời rạc hóa vật thể khảo sát bởi việc xác định các nút, các phần tử và đánh số thứ tự các nút, các phần tử. Trong ví dụ này ta chia thanh thành Ne = 2 phần tử với 3 nút ( hình (b)), 1 2 Thiết lập ma trận chỉ số [b] . Ma trận [b] = 2 3 Chỉ số cục bộ Nút đầu Nút cuối Phần tử (1) 1 2 (2) 2 3 45 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 2. Thiết lập ma trận cứng phần tử [K]e rồi thực hiện ghép nối các phần tử để xây dựng ma trận cứng tổng thể K , ( chú ý rằng với bài toán một chiều thì x’  x ). Trong ví dụ này ta có [K]1 = [K]2 và sẽ có: 1 2 EF 1 1 1 K 1 a 1 1 2 1 2 3 2 3 1 1 0 1 EF 1 1 2 => EF K K 1 1 1 1 2  2 a 1 1 3 a 0 1 1 3 2. Thiết lập véc tơ tải phần tử {p}e và thực hiện ghép nối véc tơ tải tổng thể P Do chịu lực phân bố đều có cường độ không đổi p0 và hai phần tử thanh giống nhau nên theo công thức (III.4) ta có : p .a 1 1 p .a 1 2 p 0 ; p 0 1 2 2 1 2 2 1 3 Nút 1 có phản lực liên kết R, nút 3 có tải trọng tập trung p nên véc tơ tải trọng nút p t n R  sẽ là : pn 0  , khi đó véc tơ tải tổng thể sẽ là : pt  p a 0 R 1 R 2 p0 a P 1 1 0 p0 a 2 1 p p a t 0 p 2 t => có phương trình K q P như sau : 46 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . p0 a R 1 1 0 q1 2 EF K 1 2 1 q p a (III.6) a 2 0 0 1 1 q p a 3 0 p 2 t 4. Áp đặt điều kiên biên : Nhận thấy do kết cấu của bài toán nên q1 = 0, do vậy trong hệ thống phương trình (III.6) ta “tạm bỏ đi” hàng và cột ứng với q1 trong ma trận [K], tức là hàng 1 và cột 1 => dẫn về hệ hai phương trình xác định q2 và q3 : p0 a EF 2 1 q2 p a 0 a 1 1 q3 pt 2 3a2 a 2a 2 2a 5. Giải hệ phương trình trên dẫn đến q p p ; q p p 22EF 0 EF t 3 EF 0 EF t  0 3a2 a do đó véc tơ chuyển vị nút tổng thể sẽ là : q p0 p t  (III.7) 2EF EF 2a2 2a p0 p t EF EF  p a2 p a 2 a Thay vào (III.6) xác định được phản lực liên kết : R = - q - 0 = 20 p 2 2EF EF EF t và các véc tơ chuyển vị nút phần tử sẽ là : 2  0 3a a  p0 p t q1  q2  2EF EF q  2  ; q   1 q 3a a 2 q 2 2  p0 p t 3  2a 2a 2EF EF  p p EF0 EF t  x x Có hàm dạng của chuyển vị là [N]1 = [N]2 = 1 do vậy : a a chuyển vị của mỗi phần tử : 0  x x 2 3a 1 u1(x) = [N]1 {q}1 = 1 3a a  = p0 p t x a a 2EF EF p0 p t 2EF EF  47 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 3a2 a  p0 p t x x 2EF EF u2(x) = [N]2 {q}2 = 1  = (III.8) a a 2a2 2a p p EF0 EF t  3a2 a x a2 a => u2(x) = p0 p t + p0 p t 2EF EF a 2EF EF Chú ý : biến x trong biểu thức của u1(x) và u2(x) là biến địa phương : 0 x a 1 1 1 Ma trận tính biến dạng [B] =  1 1 a a a Lực dọc trên mỗi phần tử Ne = EF. εx = EF.[B].{q}e 0  1 3a 2  N1 = EF.  1 1  = p p a 3a a 2 0 t p0 p t 2EF EF  3a2 a  p0 p t 1 2EF EF a  N2 = EF.  1 1  = p0 p t (III.9) a 2a2 2a 2 p p EF0 EF t  Để nghiên cứu kỹ hơn tính chất của PPPTHH, ta xét lời giải chính xác của bài toán này: Ta có phương trình vi phân : d du(x) EF p(x) 0 , trong đó p(x) là cường độ tải trọng phân dx dx bố dọc trục thanh, ở bài toán này p(x) = p0 - hằng số , do đó ta nhận được bài toán : d2 u(x) EF p 0 ; 0 x 2a (III.10) dx2 0 và các điều kiện biên : u(x)x 0 0 du(x) (III.11) Nx 2a EF pt dx x 2a 48 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Giải phương trình (III.10) cùng các điều kiện biên (III.11) ta nhận được nghiệm chính xác của bài toán : 1 u(x) = 2p 4p a x p x2 và N(x) = p 2p a p x (III.12) 2EF t 0 0 t 0 0 (a) – Biểu đồ chuyển vị u(x) (b) - Biểu đồ lực dọc trục N(x) Trên biểu đồ chuyển vị (a): nghiệm theo PPPTHH là đường gấp khúc ABC - đường liền nghiệm chính xác là đường cong ABC – đường nét đứt Trên biểu đồ lực dọc (b): nghiệm theo PPPTHH là đường gấp khúc ABCD - đường liền nghiệm chính xác là đường EF – đường nét đứt Như vậy để tiệm cận đến nghiệm chính xác thì ta cần tăng số nút lên , chẳng hạn số nú là 4 thì ta có biểu đồ chuyển vị và lực dọc như sau : 2. Phần tử thanh (chịu biến dạng dọc trục) trong dàn phẳng : Trong dàn phẳng, xem mỗi mắt dàn là một đỉnh nút và thanh chịu biến dạng dọc trục, do đó ở mỗi nút có 2 bậc tự do đó là chuyển vị theo các phương x’ và y’ 49 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . x x hàm chuyển vị u(x) = [N]e . {q}e = 1 q1 + q2 (III.13) L L x x như vậy hàm dạng của chuyển vị là [N]e = 1 (III.14) LL Nếu trục phần tử thanh lập với 0x’ góc α thì gọi  cosα ; m sin α ta sẽ được ma trận chuyển tọa độ như của (II.27)’ là :  m 0 0 [T]e = (III.15) 0 0 m Chú ý : Nếu phần tử thanh có nút đầu thứ i , nút cuối thứ j và có tọa độ các nút trong hệ tông thể là (xi , yi) và ( xj , yj) thì ta cùng có : xj x i y j y i 2 2  ; m với L xj x i y j y i - độ dài thanh. LL Khi đó chuyển vị nút phần tử cũng như véc tơ tải trong tọa độ địa phương và trong tọa độ tổng thể liên hệ theo các công thức sau : q'1 q m 0 0 q' 1  2 {q}e = [T]e . {q’}e =>  (III.16) q2  0 0 m q' 3 q'4 p'1 p m 0 0 p' 1  2 và {P}e = [T]e {P’}e =>  (III.17) p2  0 0 m p' 3 p'4 E.F 1 1 ma trận cứng phần tử trong tọa độ địa phương sẽ là : [K]e = L 1 1 T ma trận cứng phần tử trong tọa độ tổng thể sẽ là [K ']e [T] e [K] e [T] e hay như (II.28) : 50 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Ma trận cứng của phần tử thanh chịu kéo, nén dọc trục : 2 m  2  m E.F m m2  m m 2 [K’]e = (III.19) 2 2 L  m   m 2 2 m m  m m 1 1 1 Ma trận tính biến dạng : [B] =  1 1 (III.20) LLL Có σ [D].[B] q D ε0 σ 0 , trong trường hợp này ứng suất chỉ còn một e e e  e thành phần σx E.ε x [D].{ε} => [D] = E Lực dọc thanh : Ne = F. σx E.Fε x EF.[B]{q} e do đó lực dọc thanh Ne trong hệ tọa độ tổng thể : Ne = E.F.[B]. [T]e . {q’}e = [S’]e. {q’}e (III.21) EF  m 0 0 EF với [S’]e =  1 1 =   m  m (III.22) L 0 0 m L Công thức (III.21) xác định lực dọc trong phần tử do chuyển vị nút gây ra. Trường hợp chiều dài phần tử có các lực tác dụng dọc trục hay chịu nhiệt độ thì phải kể thêm vào (III.21) lực dọc N0 do tác dụng của lực dọc trục hay nhiệt độ gây lên trên phần tử khi xem các nút là cố định. Ta có biểu đồ của N0 như sau: 51 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Ví dụ . Tính toán hệ thanh cùng thiết diện F và mô đun đàn hồi E cho trên hình vẽ sau : Các bước thực hiện :  Lập bảng thứ tự các thanh ( phần tử) với các nút đầu , cuối :  Lập ma trận chỉ số [b]  Tính độ dài các thanh và các góc lập bởi trục các thanh ( phần tử) ( véc tơ nút đầu – nút cuối thanh) với 0x’.  Tính ma trận cứng phần tử [K’]e ( Trong hệ tọa độ tổng thể)  Ghép nối ma trận cứng tổng thể => [ K ’]  Viết véc tơ tải tổng thể  Áp các điều kiện biên  Loại các hàng, cột trong K' q ' p'’ ( liên quan đến chuyển vị đã biết : ở bài này là các chuyển vị = 0) => nhận được hệ phương trình xác định các chuyển vị chưa biết còn lại : K''' q p   '  Giải hệ phương trình : Tính được các chuyển vị qk  Tính các chuyển vị trên từng phần tử : các u(x)  Tính các lực dọc trục trên từng phần tử : các N(x) 52 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 1). Thứ tự các thanh và các nút đầu , cuối : Bảng thứ tự thanh – nút Bảng thứ tự thanh - bậc tự do ( ma trận chỉ số [b] ) TT Nút đầu Nút cuối TT Nút đầu - i Nút cuối - j Thanh ( i ) ( j ) Thanh q2i-1 q2i q2j-1 q2j (1) 1 2 (1) 1 2 3 4 (2) 2 3 (2) 3 4 5 6 (3) 3 4 (3) 5 6 7 8 (4) 4 5 (4) 7 8 9 10 (5) 5 6 (5) 9 10 11 12 (6) 6 2 (6) 11 12 3 4 (7) 2 4 (7) 3 4 7 8 (8) 2 5 (8) 3 4 9 10 2). Tính độ dài và góc α của thanh thứ i : ký hiệu L(i) và α (i) : 2a a L(1) = L(2) = L(3) = L(7) = a ; L(5) = L(6) = L(8) = ; L(4) = 3 3 α (1) = 300 ; α (2) = 1500 ; α (3) = 300 ; α (4) = 00 ; α (5) = - 600 ; α (6) = 1800 α (7) = 900 ; α (8) = 600 ; 3). Tính ma trận cứng phần tử trong tọa độ tổng thể : áp dụng (III.19)  Ta có thể viết một đoạn thủ tục ( procedure ) để tính [K’]e - là mảng 2 chiều Ke : Khai báo cosine , sine : real { biến toàn cục) ; i , j , e : integer ; Ke : array[1 4; 1 4] of real ; Procedure MaTranCungPhanTu Begin for i := 1 to 4 do for j:=1 to 4 do Ke[i,j] := 0; Ke[1,1]:= cosine^2 ; Ke[1,2]:= cosine * sine ; Ke[1,3]:= - cosine^2 ; Ke[1,4]:= - cosine * sine ; Ke[2,2]:= sine ^2 ; Ke[2,3]:= - cosine* sine ; Ke[2,4]:= - sine ^2 ; Ke[3,3]:= cosine^2 ; Ke[3,4]:= cosine * sine ; Ke[2,1]:= Ke[1,2] ; Ke[3,1]:= Ke[1,3] ;Ke[3,2]:= Ke[2,3] ; Ke[4,1]:= Ke[1,4] ; Ke[4,2]:= Ke[2,4]; Ke[4,3]:= Ke[3,4]; Ke[4,4]:= sine ^2 ; for i := 1 to 4 do for j:=1 to 4 do Ke[i,j] := Ke[i,j] * E*F/L(e) ; end ; 53 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 4). Ghép nối ma trận cứng tổng thể :  Ta có thể viết một thủ tục ghép nối ma trận cứng tổng thể như sau: Khai báo R = 6 { số nút của toàn hệ } m , n , Ne : integer {Ne – số phần tử } K : array[1 2*R; 1 2*R] of real; Procedure GhepMaTranCung Begin for m :=1 to 2*R do for n:=1 to 2*R do K[i,j] : = 0 ; { rửa mảng K } for e:= 1 to Ne do Begin cosin := cos( α (e)) ; m := sin( α (e)); {chú ý : phải đổi góc α (e) sang radiăng } MaTranCungPhanTu for i:=1 to 4 do for j:=1 to 4 do K[b[e,i], b[e,j]] := K[b[e,i], b[e,j]] + Ke[i,j] ; end ; end ; 5). Véc tơ tải tổng thể : Trong ví dụ này chỉ có véc tơ tải trọng nút : tại nút (6) py = p0; các phản lực liên kết tại nút (1) và (3) , do đó có : R(1)   Đoạn chương trình nhập các giá trị thnàh phần của véc tơ tài 1 (1) tổng thể : R 2 0 Phần khai báo: 0 P : array [1 2*R] of real { biến toàn cục- R số nút của hệ } (3) R1 Procedure VecToTai (3) R p ' 2  Begin for i:=1 to 2*R do p[i]:= 0 ; 0 p[12]:= -p0 ; 0 { Ở đây ta tạm gán cho các phản lực liên kết giá trí 0 , vì các giá 0 trị tại vị trí này sẽ bị loại khỏi hệ phương trình khi ta áp các điều 0 0 kiện biên để từ hệ K' . q ' p' ta dẫn về hệ K''' q p   p 0  end; 54 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . ' 6). Giải hệ phương trình : Tính được các chuyển vị qk 7). Tính các chuyển vị trên từng phần tử : các u(x) 8). Tính các lực dọc trục trên từng phần tử : các N(x) 3. Phần tử thanh (chịu biến dạng dọc trục) trong dàn không gian Xem mỗi mắt dàn như là một nút thì mỗi nút có 3 bậc tự do : đó là 3 chuyển vị thành phần theo 3 phương x’ , y’ , z’ của hệ tọa độ tổng thể. Do đó véc tơ chuyển vị nút phần tử với nút đầu là i nút cuối là j trong hệ tọa độ tổng thể sẽ là : u'i  q' 3i 2  Còn trong hệ tọa độ địa phương xyz : với trục thanh là trục v' q' i 3i 1 x thì véc tơ chuyển vị nút của phần tử thanh (với nút đầu i w 'i q' 3i cuối là j ) vẫn là : q'   e u' q' j 3 j 2 u  i v'j q' 3 j 1 qe  u j  w 'j  q' 3j  trong đó ui và uj là chuyển vị nút theo phương trục thanh ( trục x). Quan hệ giữa {q’}e và {q}e vẫn là {q}e = [T]e. {q’}e  m n 0 0 0 với T (III.23) e 0 0 0 m n trong đó  , m , n là các côsin chỉ phương của các trục x , y , z đối với các trục x’ , y’ , z’ trong hệ tọa độ tổng thể và chúng được tính theo công thức : x'x' y'y' z'z'  j i; m j i ; n j i (III.24) LLL 2 2 2 với L x'x'j i y'y' j i z'z' j i E.F 1 1 ma trận cứng phần tử trong tọa độ địa phương sẽ là : [K]e = L 1 1 T ma trận cứng phần tử trong tọa độ tổng thể sẽ là [K ']e [T] e [K] e [T] e tức là : 55 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 2 m  n  2  m  n m2 mn  m m 2 mn 2 2 EF n  n nm n (III.25) K' e L 2 m  n m2 mn Đối xứng 2 n Ví dụ : Tính toán chuyển vị tại các nút và lực dọc của các thanh (chịu kéo nén dọc trục) trong hệ thanh được liên kết khớp có dạng là cạnh của một chóp cụt đều : đáy và đỉnh chóp là hình vuông có cạnh bằng 4a và 2a, chiều dài cạnh bên bằng 3a. Thiết diện các thanh bằng nhau và bằng F, và có mô - dul đàn hồi E. Bốn góc đỉnh chịu lực tập trung hướng xuống dưới và đều bằng P0 . Đáy tựa trên nền cứng tuyệt đối, không ma sát. Giải: Chú ý : Khi thực hiện việc đánh số thứ tự các thanh ( phần tử)và các nút nên chú ý một điểm là : nếu có các thanh đồng dạng ( về kích thước và vị trị trong không gian) thì nên chọn thứ tự nút đầu , cuối như nhau để có các cosin chỉ hướng của trục thanh như nhau => được ma trận cứng phần tử trong tọa độ tổng thể như nhau. 1). Thứ tự các thanh và các nút đầu , cuối : Bảng thứ tự thanh – nút Bảng thứ tự thanh - bậc tự do ( ma trận chỉ số [b] ) 56 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . TT Nút TT Nút đầu – i Nút cuối - j thanh đầu i cuối j thanh q3i-2 q3i-1 q3i q3j-2 q3j-1 q3j (1) 1 2 (1) 1 2 3 4 5 6 (2) 2 3 (2) 4 5 6 7 8 9 (3) 4 3 (3) 10 11 12 7 8 9 (4) 1 4 (4) 1 2 3 10 11 12 (5) 2 8 (5) 4 5 6 22 23 24 (6) 3 7 (6) 7 8 9 19 20 21 (7) 1 5 (7) 1 2 3 13 14 15 (8) 4 6 (8) 10 11 12 16 17 18 (9) 5 8 (9) 13 14 15 22 23 24 (10) 8 7 (10) 22 23 24 19 20 21 (11) 6 7 (11) 16 17 18 19 20 21 (12) 5 6 (12) 13 14 15 16 17 18 2). Tính độ dài và góc α , β , γ của thanh thứ i : ký hiệu L(i) và α (i) ; β (i) ; γ (i) Độ dài cosin góc lập 0x cosin góc lập 0y cosin góc lập 0z TT thanh thanh L với trục thanh:  với trục thanh : m với trục thanh : n (1) 4a 1 0 0 (2) 4a 0 1 0 (3) 4a 1 0 0 (4) 4a 0 1 0 (5) 3a -1/3 1/3 7 / 3 (6) 3a -1/3 -1/3 7 / 3 (7) 3a 1/3 -1/3 7 / 3 (8) 3a 1/3 1/3 7 / 3 (9) 2a 1 0 0 (10) 2a 0 1 0 (11) 2a 1 0 0 (12) 2a 0 1 0 3). Ma trận cứng phần tử trong tọa độ tổng thể: Áp dụng công thức (III.25) (10 11 12 7 8 9) 1 2 3 4 5 6 1 0 0 1 0 0 1 10 0 0 0 0 0 0 2 11 EF 0 0 0 0 0 0 3 12 K';K'3   1 4a 1 0 0 1 0 0 4 7 0 0 0 0 0 0 5 8 0 0 0 0 0 0 6 9 57 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . (1 2 3 10 11 12) 4 5 6 7 8 9 0 0 0 0 0 0 4 1 0 1 0 0 1 0 5 2 EF 0 0 0 0 0 0 6 3 K';K'4   2 4a 0 0 0 0 0 0 7 10 0 1 0 0 1 0 8 11 0 0 0 0 0 0 9 12 4 5 6 22 23 24 1 1 7 1 1 7 4 1 1 7 1 1 7 5 EF 7 7 7 7 7 7 6 K' 5 27a 22 1 1 7 1 1 7 1 1 7 1 1 7 23 24 7 7 7 7 7 7 7 8 9 19 20 21 1 1 7 1 1 7 7 1 1 7 1 1 7 8 EF 7 7 7 7 7 7 9 K' 6 27a 19 1 1 7 1 1 7 1 1 7 1 1 7 20 21 7 7 7 7 7 7 1 2 3 13 14 15 1 1 7 1 1 7 1 1 1 7 1 1 7 2 EF 7 7 7 7 7 7 3 K' 7 27a 13 1 1 7 1 1 7 1 1 7 1 1 7 14 15 7 7 7 7 7 7 58 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 10 11 12 16 17 18 1 1 7 1 1 7 10 1 1 7 1 1 7 11 EF 7 7 7 7 7 7 12 K' 8 27a 16 1 1 7 1 1 7 1 1 7 1 1 7 17 18 7 7 7 7 7 7 (16 17 18 19 20 21) 13 14 15 22 23 24 1 0 0 1 0 0 13 16 0 0 0 0 0 0 14 17 EF 0 0 0 0 0 0 15 18 K';K'9   11 2a 1 0 0 1 0 0 22 19 0 0 0 0 0 0 23 20 0 0 0 0 0 0 24 21 ( 13 14 15 16 17 18 ) 22 23 24 19 20 21 0 0 0 0 0 0 22 13 0 1 0 0 1 0 23 14 EF 0 0 0 0 0 0 24 15 K';K'10   12 2a 0 0 0 0 0 0 19 16 0 1 0 0 1 0 20 17 0 0 0 0 0 0 21 18 Thực hiện ghép nối ma trận cứng tổng thể K' 4) Véc tơ tải tổng thể : Do các thanh đáy trượt trên nền cứng tuyệt đối và có trọng tải tập trung tại các đỉnh đều bằng P0 nên ta có véc tơ tải tổng thể P' : 59 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 0  0 Chú ý Vì có các chuyển vị : q ' 3 = q ' 6 = q ' 9 = q ' 12 = 0 R 1,3 do đó trong hệ phương trình K' .q ' = P' ta tạm loại 0 0 đi các hàng và cột : 3 , 6 , 9 , 12 như vậy sẽ dẫn R 2 ,3 về hệ phương trình K' . q' = P' 0 * * * 0 ( còn lại 20 phương trình , xác định 20 ẩn q ' k ) R 3,3 0 0 R 4 ,3 P'  0 0 p 0 0 0 p 0 0 0 p 0 0 0 p 0  60 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . III.2 Hệ khung Do liên kết giữa các thanh là các liên kết cứng nên trong các thanh có thể xuất hiện biến dạng dọc , biến dạng xoắn tức là khi đó trong thanh có thể có lực dọc trục , lực trượt và các mô men uốn 1. Phần tử dầm chịu uốn – dầm liên tục Xét phần tử dầm chịu uốn , dầm có chiều dài L , thiết diện không đổi F. Do dầm chịu uốn nên có thể bỏ qua biến dạng dọc trục và khi đó chuyển vị tại x bất kỳ được đặc trưng bởi chuyển vị v(x) vuông góc với trục dầm và v(x) là một đa thức xấp xỉ bậc 3 của x : 2 3 v(x) = a1 + a2x + a3x + a4x a1  2 3 a2 hay với [P(x)] = [ 1 x x x ] và {a} =  a3 a4  vậy ta có v(x) = [P(x)].{a} (III.27) Do sự tương thích biến dạng giữa phần tử đang xét và phần tử kề bên nên tại các nút cần có sự liên tục của chuyển vị và góc xoay θ dv(x) 2 2 θ = a2 + 2a3x + 2a4x = [ 0 1 2x 3x ] . {a} (III.28) dx Như vậy số bậc tự do tại mỗi nút là 2 dẫn đến số bậc tựa do cho mỗi phần tử dầm ( 2 nút ) là 4 , tức là ta có : v1  q 1  θ1 q 2 qe   , do có (III .27) và (III.28) tính tại x = 0 và x = L ta có được : v2 q 3 θ2  q 4  61 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 1 0 0 0 v1  q 1  a 1  0 1 0 0 θ1 q 2 a 2 q    hay {q}e = [A] . {a} e v q 1 L L2 L 3 a 2 3 3 θ q 2 a 2  4  0 1 2L 3L 4  1 0 0 0 0 1 0 0 có [A]-1 = 3/L2 2/L 3/L 2 1/L 3 3 2 2/L 1/L 2/L 1/L -1 và { a } = [A] . {q}e -1 Thay vào (III.27) ta có được v(x) = [P(x)]. [A] . {q}e với [N] là ma trận hàm dạng của chuyển vị : -1 [N] = [ P(x)] .[A] = [ N1 N2 N3 N4 ] với 3x2 2x 3 2x 2 x 3 Nv 1 ; N v x 12 3 2 L 2 LLL 3x2 2x 3 x 2 x 3 N;Nv v 3LLL2 3 4 L 2 Khi đó hàm chuyển vị v(x) của phần tử dầm chịu uốn là : 4 v v v(x) = [N (x)] . {q}e =  Ni .q i (III.29) i 1 ө -1 Còn ma trận hàm dạng của góc xoay Ө(x) sẽ là : [N ] = [Pө(x)] . [A ] 1 0 0 0 0 1 0 0 ө 2 ө ө ө ө [N ] = [0 1 2x 3x ] . = [N 1 N 2 N 3 N 4 ] 3/L2 2/L 3/L 2 1/L 3 3 2 3 2 2/L 1/L 2/L 1/L 4 ө θ Ө(x) = [N (x)] . {q}e = Ni .q i (III.30) i 1 2 2 2 2 ө 6x 6x ө 4x 3x ө 6x 6x ө 2x 3x với N 1 = ; N 2 = 1 ; N 3 = ; N 4 = LL2 3 L L2 LL2 3 LL3 2 62 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Gọi y là khoảng cách từ điểm đang xét đến đường trung hòa thì chuyển vị dọc trục u(x) dv(x) và độ võng v(x) có quan hệ : u(x) y , do đó biến dạng dọc trục dx du d2 v(x) d 2 [N] d2 [N] ε y y {q} [B].{q} với [B] y xdx dx2 dx 2 e e dx2 như vậy [B] là ma trận hàng 6 12x 4 6x 6 12x 2 6x [B] y 2 3 2 2 3 2 (III.31) LLLLLL LL Ứng suất tại mọi điểm của dầm chịu uốn : σx  [D]. ε với [ D] = E Ma trận cứng của phần dầm chịu uốn : K [B]TTT [D].[B] dv E [B] .[B]dF .dx = E.F. [B] [B]dx =  e VLFLe do đó : 12 6L 12 6L 2 2 EJ 6L 4L 6L 2L K (III.32) e L3 12 6L 12 6L 2 2 6L 2L 6L 4L với J y2 dF - mô men quán tính của mặt cắt ngang đối với đường trung hòa. F 63 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . bh3 Nếu thiết diện dầm là chữ nhật : rộng - b ; cao - h thì J = 12 Véc tơ tải {P}e : Nếu trên chiều dài phần tử thanh có tải trọng phân bố p(x) , có lực tập trung Qi và các mô men tập trung Mi tác dụng thì trên cơ sở công thức (II.14) ( hoặc từ sự cân bằng về công sinh ra trên các chuyển vị nút với công của ngoại lực) ta có được công thức tính véc tơ tải phần tử sau : T nq nM T T dN {P}e {N} p(x)dx  N(x Qi ) Q i  (x Mi ) M i L i 1 i 1 dx 1) p(x) là cường độ lực ngang trên dầm 2) mô men uốn Mi tập trung tại xMi 3) Lực tập trung Qi tại xQi  Trường hợp tải phân bố đều với cường độ p0 : 2 3 3x 2x  Lp  1 0 LL2 3 2 v N  2 3 2 p1  1 2x x L p0 v x 2 p2 N 2 L p .p dx = L .p dx = 12 e  v 0 2 3  0  p3 LLN 3x 2x Lp0 3 v 2 3 (III.33) p4  LL 2 N4  2 3 2 x x L p0 L L2  12  Trường hợp lực tập trung có giá trị Q đặt cách nút đầu phần tử khoảng cách a : 3a2 2a 3  1 2 3 LL 2 3 p1  2a a a 2 p2 v T L p N (a) .Q Q L e  2 3  p3 3a 2a 2 3 p4  LL a2 a 3 (III. 34) L L2  đặc biệt nếu a = L/2 ( tức là Q đặt giữa phần tử dầm) , khi đó có 64 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 1  2 p1  L p2 T 8 pe   N(a) .Q Q  (III.35) p3 1 p 2 4  L 8  Mô men uốn nội lực trong phần tử dầm : nếu dầm có thiết diện không đổi thì: d2 v d 2 M(x) EJ EJ [N]{q} EJ [N ] {q} (III.36) dx2 dx 2 e e v v v v với [N ] N1 N 2 N 3 N 4 v v Để ý rằng các hàm Ni (x) là các đa thức bậc 3, nên Ni (x) là các hàm bậc nhất dẫn đến M(x) là hàm bậc nhất , do vậy ta chỉ cần tính được các M(nút 1) , M(nút 2) => M(x) M(nút1)  Gọi {M}e =  là véc tơ mô men uốn tại các nút 1 và nút 2 của phần tử dầm, M(nút 2)  M  N (x 0) => M 1 EJ. q S q (III.37) e   e e  e M2  N (x L) M  hay là M 1 S q e   e  e M2  từ đây ta có được biểu đồ mô men Mq(x) – là đoạn thẳng đi qua điểm (0, M1) và (L , M2) EJ 6LLLL 42 6 2 2 Với S (III.38)  e 3 2 2 L 6LLLL 2 6 4 Chú ý trong trường hợp chiều dài phần tử có lực tác dụng cần cộng thêm vào Mq(x) mô men M0(x) do tải trọng gây ra trên phần tử khi xem tất cả các nút được ngàm 65 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Các kết quả cần nhớ : Phần tử dầm độ dài L, chịu tải trọng đều với cường độ q thì ta có : q qL qL2 M (x) x2 x (III.39) 0 2 2 12 Có biểu đồ mô men M0(x) như hình bên Phần tử dầm độ dài L , chịu tải tập trung tại điểm cách nút đầu một khoảng là a có độ lớn p , khi đó có : a a M0 (x) p 1 x với 0 x a L 2 (III. 40) pa L a M0 (x) x với a x L L 2 Ví dụ : Tính toán chuyển vị và nội lực trong dầm liên tục được cho trên hình vẽ : Đoạn AB chịu tải trong phân bố đều với cường độ q ; đoạn BC chịu tải tập trung cường độ P0 = tại điểm giữa đoạn BC . Vẽ biểu đồ mô men uốn khi p0 = q.a 66 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Giải . 1) Rời rạc hóa kết cấu : Chia dầm thành 2 phần tử : Đoạn AB và BC, đánh số phần tử , số nút dẫn đến ma trận chỉ số [b] , chú ý ở mỗi nút có 2 bậc tự do : chuyển vị phương vuông góc trục dầm, góc xoay. Chỉ số đphương Nút đầu (i) Nút cuối (j) Phần tử (1) 1 2 3 4 (2) 3 4 5 6 2) Thiết lập ma trận cứng phần tử Sử dụng công thức (III.32) 1 2 3 4 3 4 5 6 12 6a 12 6a 1 12 6a 12 6a 3 2 2 2 2 EJ 6a 4a 6a 2a 2 EJ 6a 4a 6a 2a 4 K K 1 a3 12 6a 12 6a 3 2 a3 12 6a 12 6a 5 2 2 2 2 6a 2a 6a 4a 4 6a 2a 6a 4a 6 Chú ý : * Hệ tọa độ tổng thể và địa phương cùng hướng (trục x) nên không có sự đổi tọa độ. 67 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . * Để tiện lợi nên viết véc tơ tải, và các điều kiện biên để từ đó loại ra các phương trình dẫn đến ma trận cứng phần tử [ K* ] Xét điều kiện biên: các chuyển vị nút 1 và nút 2 và 3 bằng 0 , góc xoay tại nút 1 bằng 0 => q1 = q2 = q3 = q5 = 0 do vậy trong ma trận [ K ] ta bỏ đi các hàng, cột 1 , 2 , 3 , 5 dẫn đến ma trận [ K* ] còn lại : 4 6 EJ 8a2 2a 2 4 K* 3 2 2 a 2a 4a 6 3) véc tơ tải phần tử aq  2 2 aq  1 1  3 a q 2 2 12 a2 q a aq p 2 4 0 12 8 2 2 p  ; p p0  => tải tổng thể p  1 aq 2 1 a2 q ap 3 5 0 2 2 12 8 2 a a q p0 4 6 12  8  2 ap 0 8  * q4  chú ý rằng trong véc tơ chyển vị nút chỉ còn thành phần q4 , q6 => q   q6  a2 q ap  0 * 12 8 và véc tơ tải p   ap0 8  4) Giải hệ phương trình : 2 a q ap0  2 2 EJ 8a 2a q4  12 8 Phương trình K q p   =    3 2 2 q a 2a 4a 6  ap0 8  68 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . aq p0  2 4 1 q4  a 3 2 =>   1 2 q6  8EJ p0 2  2 2 a 2aq 3p0 a aq 5p0 giải ra ta được : q4 và q6 56EJ 3 2 56EJ 3 2 5) Tính các đại lượng còn lại : M(x) ; v(x) Sử dụng (III.37) và (III.38) ta có 0  0  0  0  a2 2 aq 3p 0 0 0 q4 56EJ 3 2 có q  0  ; q   ; 10 2 0 0 2 a 2 aq 3p q 0 q 2 4  6  a aq 5p0 56EJ 3 2  56EJ 3 2  M S q 1 e  1 0  3ap 2a2 q  0 0 2EJ 3 2a 3 a 1 0 2 3 M 2   1 a 3a 3 2 a 28 6ap 4a2 q 2 0 a 2 aq 3p0 2 3  56EJ 3 2  0  2 a 2 aq 3p0 ap 2  0 a q 2EJ 3 2a 3 a 56EJ 3 2 1 2 M 2 2   a 3a 3 2 a 0 28 7 ap0 2 2  a aq 5p0 56EJ 3 2  Nếu p0 = qa khi đó ta có 69 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 2  3ap0 2a q ap0 2  2 a q 2 1 2 3 a q 1  1 2 a q 1  M ; M 1 2   2   286ap 4a q 84 5  287ap0 56 7  0 2 3  2  q qa qa2 M (x) x2 x ; 01 2 2 12 qa a a M02 (x) x với 0 x ; 2 4 2 qa 3 a a M02 (x) x với x a ( x tọa độ địa phương trên phần tử) 2 4 2 Biểu đồ mô men uốn : Giá trị của các M(x) được nhân với qa2 Kết quả cuối cùng ta được các mô men uốn : qx23 qax qa 2 M (x) ;0 x a 1 2 7 14 70 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 9qa qa2 a x khi0 x 14 7 2 M2 (x) 5qa 5 qa2 a x khi x a 28 28 2 2. Phần tử dầm chịu uốn trong hệ khung phẳng Xét phần tử dầm chịu uốn trong hệ khung phẳng, độ dài phần tử dầm là a , có diện tích thiết diện và mô đun đàn hồi là F và E không đổi, trong phần tử dầm chỉ có chuyển vị v(x) theo phương vuông góc với trục dầm, do đó véc tơ chuyển vị nút phần tử: v1  q 1   q q 1 2 e   v2 q 3 2  q 4  Mỗi chuyển vị vi sẽ có 2 thành phần theo phương x' và y' trong tọa độ tổng thể , do vậy véc tơ chuyển vị nút phần tử trong tọa độ tổng thể sẽ có 6 thành phần : và {q}e = [T]e . {q'}e với u1  q1   ym y 0 0 0 0 v1 q 2 0 0 1 0 0 0 T (III.41)  e 1 q 3 0 0 0 m 0 q    y y e u q 2 4 0 0 0 0 0 1 v q 2 5  2  q 6  các  y và my là các cosin chỉ phương của 0y với các trục 0x' , 0y'. nếu gọi là góc lập bởi trục phần tử dầm với 0x' thì  y = - sin ; my = cos , 71 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . do vậy ta vẫn gọi  = cos ; m = sin và sẽ dẫn đến: m  0 0 0 0 0 0 1 0 0 0 T (III.42) e 0 0 0 m  0 0 0 0 0 0 1 T sử dụng (III.32) và [K']e = [T]e [K]e . [T]e ta được : Ma trận cứng của phần tử dầm chịu uốn trong khung phẳng: 12m2 12 m 6 am 12 m 2 12  m 6 am 122 6a  12  m 12  2 6 a  EJ 4a2 6 am 6 a 2 a 2 K  3 (III.43) e a 12m2 12 m 6 am đối xứng 122 6a  2 4a Mô men uốn tính theo trong tọa độ tổng thể : Theo (III.37) và (III.38) ta có 2 2 M1  EJ 6a 4 a 6 a 2 a M   S q với S eM e e e 3 2 2 2  a 6a 2 a 6 a 4 a Thay {q}e = [T]e . {q'}e dẫn đến : Me  Se q e  S e T e q e  S e q  e với SSTe   e  e EJ 6am 6 a 4 a2 6 am 6 a  2 a 2 S (III.44)  e 3 2 2 a 6am 6 a 2 a 6 am 6 a  4 a 3. Phần tử dầm chịu uốn và kéo , nén trong khung phẳng Chuyển vị tại mỗi nút của phần tử dầm sẽ theo 2 phương : trục dầm và vuông góc với trục dầm . Véc tơ chuyển vị nút của phần tử dầm sẽ là : 72 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . u1  u1  q1  v v q 1 1 2 1 1 q 3 q  và trong tọa độ tổng thể sẽ là : q    (III.45) e e u2 u2 q4 v v q 2 2 5 2   2  q 6  Do tính chất cộng tác dụng nên ta sử dụng ở đây các công thức (III.19) và (III.43) ta được: Ma trận cứng phần tử dầm chịu uốn và kéo nén dọc trục trong tọa độ tổng thể : (III.46) F.2 12 m 2 H F  m 12  mH 6 amH F.  2 12 m 2 H F  m 12  mH 6 amH Fm2 12 2 H 6 aH  Fm  12  mH Fm 2 12  2 H 6 aH  E 4a2 H 6 amH 6 a H 2 a 2 H K  đối xứng e a F.2 12 mH 2 Fm  12  mH 6 amH Fm2 12 2 H 6 a  H 2 4a H trong đó :  = cos ; m = sin ; - góc lập bởi trục phần tử dầm và 0x' J a – độ dài thanh dầm F – diện tích thiết diện ; H = a2 Ví dụ Tính chuyển vị và các nội lực của các thanh trong khung phẳng chịu kéo nén dọc trục và chịu uốn. Các thanh đều có độ dài a, diện tích thiết diện F, mô đun đàn hồi E, chịu tác dụng lực ngoài và các điều kiện biên như hình vẽ. 0 ( lực p0 lập với các thanh góc 45 ) 73 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Giải: 1) Rời rạc hóa các kết cấu: đánh số phần tử và số nút => tính ma trận chỉ số [b] Phần tử Nút đầu Nút cuối (1) 1 2 (2) 2 3 => ma trận chỉ số [b] Chỉ số đ.ph Nút đầu i Nút cuối j Phần tử 1 2 3 4 5 6 (1) 1 2 3 4 5 6 (2) 4 5 6 7 8 9 2)Thiết lập ma trận cứng tổng thể: Chú ý rằng tại các nút 1 và 3 do bị ngàm chặt nên các chuyển vị và góc xoay bằng 0 tức là ta có q'1 = q'2 = q'3 = q'7 = q'8 = q'9 = 0 nên trong các ma trận cứng phần tử các hàng hoặc các cột là các số 1 ; 2 ; 3 ; 7 ; 8 ; 9 thì được "loại ra" do vậy áp dụng (III.46) ta được : Với phần tử (1) :  0; m 1 2 4 5 6 12H 0 6 aH 4 E KF * 0 0 5 e a 2 6aH 0 4 a H 6 Với phần tử (2) : 0  1; m 0 4 5 6 F 0 0 4 E K * 0 12 H 6 aH 5 e a 2 0 6aH 4 a H 6 74 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . * Ma trận cứng tổng thể K ( đã tính đến các điều kiện biên) 12H F 0 6 aH E J K * 0 12 H F 6 aH ; H = a a2 2 6aH 6 aH 8a H 3) Véc tơ tải tổng thể : RR1  1  1 y y RR1 1 2 x x 0 0 3 Có p    1 4 Q2 qa Q qa 5 1 6 0  0  Véc tơ tải trên phần tử dầm (2) do tải trọng đều q và trọng tải nút P0 tác động. aq  2 a2 q q 12 Theo (III.13) ta có p  2 aq 2 2 a q 12  q q mà p [T]T p 21  2 áp dụng (III.42) với : 0 nên  1; m 0 , do đó ta được 75 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . aq  aq  2 Q  1 0 0 0 2 0 2 2 2 Q 0 0 0 0 a q a q 1 0 q 0 1 0 0 12 12 nut p    ; véc tơ tải trọng nút p' ()2  2 0 0 1 0 aq aq 2 R x 0 0 0 0 2 2 R()2 y 0 0 0 1 a2 q 0 0  12  a2 q 12  Do vậy véc tơ tải phần tử (2) trong tọa độ tổng thể là : aq  3 aq  4 aq  Q 22 2 2 Q  0 2 Q1 aq 5 2 Q 2 2 a q 1 a q a q 0 nut q 12 12 12 6 p p p 2  2  2  ()2    aq Rx aq()()2 aq 2 7 RRx x 2 R()2 2 2 y ()()2 2 0 RRy y 8 0  a2 q 2 2 a q a q 12  12  12  9 Vậy có véc tơ tải tổng thể : R()1  1 y ()1 2 Rx 3 0 5aq 2 4 5aq  2qa 5 2 * p ' a2 q  => p' 2 qa  2 12 6 a q aq ()2 12  Rx 2 7 R()2 8 y a2 q 12  9 76 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . và ta có hệ phương trình : 5aq  12H F 0 6 aH q4  2 E J K q 0 12 H F 6 aH q5  2 qa  với H = 2 a a 2 2 6aH 6 aH 8a H q6  a q 12  Giải hệ phương trình trên ta tính được các giá trị q4 , q 5 ,q 6 Chú ý : Để giải hệ phương trình trên ta nên đưa các ẩn cần tìm về cùng một thứ nguyên: 12H F 0 6 H q4  30  a2 q 0 12H F 6 H q 24 5   12E 6H 6 H 8H aq6  1  4) Tính các nội lực: mô men uốn Chú ý: khi tính mô men uốn trên phần tử thanh (2) cần cộng thêm mô men M0(x) theo (III.39) III.3 Hệ khung không gian 1. Ma trận cứng phần tử : Ta xét phần tử khung không gian là dầm thẳng có thiết diện không đổi, trên mặt cắt ngang tồn tại lực dọc, mô men uồn tròn cả hai mặt phẳng quán tính chính và mô men xoắn. Các bậc tự do chuyển vị đặc trưng cho trạng thái chuyển vị - biến dạng của phần tử dầm 2 nút được biểu diễn trên hình. Hệ tọa độ địa phương xyz với trục x là trục dầm, y và z là hai trục chính của mặt phẳng cắt ngang Véc tơ chuyển vị phần tử hai nút là : T {q}(e) = { q1 q2 q3 q4 q5 q6 q7 q8 q9 q10 q11 q12} 77 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . Trong đó : q1 và q7 chuyển vị dọc trục dầm gây ra biến dạng dọc trục dầm q2 và q8 chuyển vị thẳng theo trục y gây ra biến dạng uốn trong mặt phẳng xy q6 và q12 : góc xoay trong mặt phẳng xy q3 và q9 chuyển vị thẳng theo trục z gây ra biến dạng uốn trong mặt phẳng xz q5 và q11 : góc xoay trong mặt phẳng xz gây ra biến dạng xoắn của thanh q4 và q10 góc xoắn quanh trục x Như vậy trong 12 bậc tự do chuyển vị sẽ gây ra 4 nhóm biến dạng độc lập nhau. Vì vậy ma trận cứng phần tử [ K ](e) có kích thước 12 × 12 được thành lập từ 4 ma trận con như sau : a) Biến dạng dọc trục ( q1 và q7) q1 q 7 E.F 1 1 q1 (III.47) [K](e) L 1 1 q7 78 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . b) Biến dạng xoắn theo trục x ( q4 và q10) Thấy rằng đa thức xấp xỉ của hàm góc xoắn θ (x) là hàm bậc nhất : θ (x) = a1 + a2x ( 0 ≤ x ≤ L ) hay : a1 θ(x) [1 x ] . = [P(x)] . {a} , theo kết quả của ví dụ 1 phần (II.2) ta có : a2 biểu diễn xấp xỉ hàm góc xoắn θ (x) theo các góc xoắn tại các nút: x x θ (x) = [N]e . { q}xoắn = 1 q4 + q10 (III.48) L L x x ở đây [N]e = 1 LL d Nếu là thanh tròn thì biến dạng xoắn sẽ là  r x và ứng suất tiếp  G.  yz dx xy xy ( theo Hooke và G là mô đun trượt , r là khoảng cách từ tâm đến điểm tính toán) áp dụng ε xoắn =  θxoắn = [N]{q} xoắn = [B] {q} xoắn σ xoắn = { τyz }= [T] q xoắn (III.49) trong đó : [T] = [D]. [B] là ma trận tính ứng suất gọi [Kxoắn ]e - Ma trận cứng phần tử T [Kxoắn ]e [B] [D][B]dv (III.50) Ve ở đây ta có { } = {  yz } và { } = [ B ] . {q }(xoắn) r r  với B =  ; [D ] = [ G ] LL  Do đó ta được : 79 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . 1 L 2 L 1 1 [Kxoắn ]e = G dx r dF. . (III.51) 1 LL 0 F L q4 q 10 => [K ] = q (III.52) xoắn e GJx 1 1 4 L 1 1 q10 R 4 ở đây Jx là mô men quán tính độc cực của mặt cắt ngang: J ; R – bán kính dầm x 2 Trong trường hợp nếu mặt cắt ngang dầm là hình chữ nhật có các cạnh a × b thì : 3 Jx = c ab với c là hằng số được lấy theo bảng sau ( phụ thuộc vào tỷ số a/b) a/b 1 1,5 2 3 5 10 c 0,141 0,196 0,229 0,263 0,291 0,312 T c) Biến dạng uốn trong mặt phẳng xz ( do { q}xz = { q3 , q5 , q9 , q11} ) q3 q 5 q 9 q 11 12 6L 12 6L q 3 áp dụng EJ 6L 4L2 6L 2L 2 q (III.53) K  y 5 xz e L3 12 6L 12 6L q 9 2 2 q 6L 2L 6L 4L 11 2 với Jy z dF - mô men quán tính của mặt cắt ngang đối với đường trung hòa. F 80 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012
Bài giảng : Phương pháp phần tử hữu hạn – Bộ môn Toán ĐHLN . ab3 Nếu thiết diện dầm là chữ nhật có kích thước các cạnh như hình vẽ thì Jy = 12 R 4 Nếu thiết diện tròn bán kính R thì : J y 4 T d) Biến dạng uốn trong mặt phẳng xy ( do { q}xz = { q2 , q6 , q8 , q12} ) q2 q 6 q 8 q 12 12 6L 12 6L q2 2 2 áp dụng EJ 6L 4L 6L 2L q (III.54) K z 6 xy e L3 12 6L 12 6L q 8 2 2 q 6L 2L 6L 4L 12 2 với Jz y dF - mô men quán tính của mặt cắt ngang đối với đường trung hòa. F a3 b Nếu thiết diện dầm là chữ nhật có kích thước các cạnh như hình vẽ thì Jz = 12 R 4 Nếu thiết diện tròn bán kính R thì : J z 4 81 Biên soạn Vũ Khắc Bảy – Tháng 08 năm 2012