›
›

Luận văn Nghiên cứu dao động ngẫu nhiên phi tuyến bằng phương pháp tuyến tính hóa tương đương

pdf 67 trang phuongnguyen 2330
Download
Bạn đang xem 20 trang mẫu của tài liệu "Luận văn Nghiên cứu dao động ngẫu nhiên phi tuyến bằng phương pháp tuyến tính hóa tương đương", để 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:

  • pdfluan_van_nghien_cuu_dao_dong_ngau_nhien_phi_tuyen_bang_phuon.pdf

Nội dung text: Luận văn Nghiên cứu dao động ngẫu nhiên phi tuyến bằng phương pháp tuyến tính hóa tương đương

  1. ®¹i häc quèc gia hµ néi viÖn khoa häc vµ c«ng nghÖ viÖt nam Tr−êng ®¹i häc c«ng nghÖ viÖn c¬ häc NGUYÔN NH¦ HIÕU NGHI£N CøU DAO §éng ngÉu nhiªn phi tuyÕn b»ng ph−¬ng ph¸p tuyÕn tÝnh hãa t−¬ng ®−¬ng LuËn v¨n th¹c sÜ Hµ Néi - 2011
  2. ®¹i häc quèc gia hµ néi viÖn khoa häc vµ c«ng nghÖ viÖt nam Tr−êng ®¹i häc c«ng nghÖ viÖn c¬ häc NGUYÔN NH¦ HIÕU NGHI£N CøU DAO §éng ngÉu nhiªn phi tuyÕn b»ng ph−¬ng ph¸p tuyÕn tÝnh hãa t−¬ng ®−¬ng Ngµnh: C¬ häc Chuyªn ngµnh: C¬ häc VËt thÓ r¾n MÉ sè: 60 44 21 LuËn v¨n th¹c sÜ Ng−êi h−íng dÉn khoa häc: ts. TrÇn D−¬ng TrÝ Hµ Néi - 2011
  3. Lời cam đoan Tôi xin cam đoan đây là công trình nghiên cứu của tôi và chưa được công bố trong bất kỳ công trình nào khác. Các kết quả trong luận văn là trung thực. Người cam đoan Nguyễn Như Hiếu
  4. Lời cảm ơn Tôi chân thành cám ơn các thầy, cô của trường Đại học Công nghệ, Đại học Quốc gia Hà Nội, đặc biệt là các thầy, cô trong Khoa Cơ học kỹ thuật và Tự động hóa đã giúp đỡ và chỉ bảo tận tình trong suốt thời gian tôi học tập tại Khoa. Tôi rất cám ơn Phòng Cơ học Công trình, Viện Cơ học đã tạo điều kiện cho tôi học tập và nghiên cứu tại đây. Tôi gửi lời cám ơn chân thành tới TS. Trần Dương Trí, người đã quan tâm chỉ bảo trong thời gian tôi thực hiện luận văn này. Đặc biệt tôi gửi lời cám ơn chân thành tới GS. Nguyễn Đông Anh và GS. Issac Elishakoff vì những kiến thức bổ ích trong nhiều năm học của tôi.
  5. 1 Mục lục MỞ ĐẦU 3 Chương 1. Tổng quan về phương pháp tuyến tính hóa tương đương 5 1.1. Phương pháp tuyến tính hóa tương đương cho hệ nhiều bậc tự do 5 1.1.1. Phương trình vi phân chuyển động 5 1.1.2. Hệ tuyến tính hóa tương đương 6 1.1.3. Ma trận mật độ phổ 11 1.2. Phương pháp tuyến tính hóa tương cho hệ một bậc tự do 13 1.3. Áp dụng phương pháp tuyến tính hóa tương đương cho các hệ Duffing và Van der Pol 15 1.3.1. Hệ Duffing 15 1.3.2. Hệ Van der Pol 17 Chương 2. Một số mở rộng của phương pháp tuyến tính hóa tương đương 20 2.1. Phương pháp tuyến tính hóa điều chỉnh 20 2.2. Áp dụng phương pháp tuyến tính hóa điều chỉnh cho các hệ Atalik-Utku, Lutes-Sarkani và Van der Pol 23 2.2.1. Hệ Atalik-Utku 23 2.2.2. Hệ Lutes-Sarkani 26 2.2.3. Hệ Van der Pol 30 2.3. Mở rộng của phương pháp tuyến tính hóa điều chỉnh: điều chỉnh hai bước 31 2.3.1. Phương pháp điều chỉnh hai bước cho hệ Atalik-Utku 32 2.3.2. Phương pháp tuyến tính hóa điều chỉnh hai bước cho hệ Lutes-Sarkani 33 Chương 3. Áp dụng phương pháp tuyến tính hóa tương đương cho bài toán dao động của dầm 35 3.1. Phương trình dao động của dầm 35 3.2. Phương pháp tuyến tính hóa điều chỉnh cho bài toán dao động của dầm 39 KẾT LUẬN 48 TÀI LIỆU THAM KHẢO 49 PHỤ LỤC 50
  6. 2 DANH MỤC CÁC BẢNG VÀ HÌNH VẼ CÁC BẢNG Trang Bảng 1. Đáp ứng bình phương trung bình của hệ Duffing với các tham số 2 h 0.5, 0 1,  2 17 Bảng 2. Đáp ứng bình phương trung bình của hệ Van der Pol với các tham số 0.2, 0 1,  2 18 Bảng 3. Phương sai của hệ Lutes-Sarkani với các giá trị khác nhau của a 29 Bảng 4. Đáp ứng bình phương trung bình của hệ Van der Pol với các phương pháp khác nhau ( 0.2, 0 1,  2 ) 31 CÁC HÌNH Hình 1. 2 Đáp ứng bình phương trung bình E w1 theo tham số R với S0 1, 0,  0.1 45 Hình 2. 2 Đáp ứng bình phương trung bình E w1 theo tham số R với S0 1, 1,  0.1 45 Hình 3. 2 Đáp ứng bình phương trung bình E w1 theo tham số R với S0 5, 1,  0.1 46 Hình 4. 2 Đáp ứng bình phương trung bình E w1 theo tham số với S0 5,  0.1,R 1 46
  7. 3 MỞ ĐẦU Dao động là một trong hiện tượng xảy ra phổ biến trong tự nhiên. Nó xuất hiện trong các lĩnh vực như vật lý, hóa học, sinh học, kinh tế và kỹ thuật. Từ các lĩnh vực này, có lớp các bài toán quan trọng là dao động ngẫu nhiên phi tuyến của các hệ động lực. Ta thường bắt gặp những hệ ngẫu nhiên phi tuyến trong thực tế như dao động của các kết cấu xây dựng, nhà cao tầng hay những cây cầu dây văng chịu tác động của tải trọng gió hay kích động động đất, các công trình cảng sông, cảng biển hay những giàn khoan chịu tác động của các tải trọng sóng hay cũng có thể là dao động của các máy trong các nhà máy xí nghiệp trong quá trình hoạt động của chúng Ta có thể đặt vấn đề làm thế nào để tăng cường tuổi thọ và duy trì độ bền của các hệ cơ học nói trên Nhiều mô hình toán học được đưa ra để phục vụ thực tiễn đó. Các phương trình toán học được mô tả và giải quyết dưới nhiều phương diện khác nhau. Với hệ động lực phi tuyến, người ta bắt gặp các phương trình phi tuyến yếu và các phương trình phi tuyến mạnh. Phương trình phi tuyến yếu được quan tâm nghiên cứu và phát triển với nhiều phương pháp khác nhau trong những thập kỷ gần đây. Có thể kể đến một trong những phương pháp phổ biến nhất là phương pháp tuyến tính hóa tương đương hay phương pháp tuyến tính hóa thống kê. Đây là phương pháp được đưa ra đồng thời trong những năm 50 của thế kỷ trước bởi các tác giả Booton [1], Kazakov [2], Caughey [3, 4]. Tuy nhiên ý tưởng của phương pháp này đã được nhen nhóm từ trước đó. Ban đầu phương pháp tuyến tính hóa được trình bày cho các hệ tiền định, cơ sở toán học của nó được đề cập trong [5] bởi Krylov và Bogoluboff. Đến Caughey, ông áp dụng phương pháp tuyến tính hóa cho các hệ phi tuyến chịu kích động ngẫu nhiên. Ông gọi phương pháp này là “phương pháp tuyến tính hóa tương đương”. Còn tên gọi “phương pháp tuyến tính hóa thống kê” được trình bày bởi Booton và Kazakov. Điều thú vị là phương pháp không ngừng được cải tiến và được đóng góp bởi nhiều tác giả [6-23] sao cho nó giải quyết phù hợp với từng loại bài toán khác nhau chẳng hạn các bài toán liên quan đến các không gian trạng thái, miền các tần số, không gian các hàm đặc trưng [11]. Phương pháp dựa trên các tiêu chuẩn tuyến tính hóa để tìm ra các công thức dạng ẩn hoặc dạng hiện cho hệ số tuyến tính hóa. Hệ số này phụ thuộc vào đặc trưng đáp ứng chưa biết (như giá trị trung bình, các tương quan, các mô men bậc cao ). Tuy nhiên khi áp dụng phương pháp vào các phương trình phi tuyến mạnh thì gặp phải các sai số lớn hơn so với việc áp dụng nó vào các hệ phi tuyến yếu. Vì vậy nhu cầu cải tiến phương pháp là cần thiết cho việc giải quyết các hệ phi tuyến mạnh. Điều này là một trong những mấu chốt hình thành nhiều cải tiến mới đây [12-15]. Một trong những cải tiến đó được trình bày trong [12], trong đó các tác giả Nguyễn Đông Anh và Di Paola giải quyết bài toán cho các hệ phi tuyến bằng một phương pháp với tên gọi “phương pháp tuyến tính hóa điều chỉnh”. Phương pháp này dựa trên ý tưởng rằng thành phần phi tuyến ban đầu không được tuyến hóa trực tiếp như phương pháp tuyến tính hóa
  8. 4 kinh điển mà nó được thay thế bằng thành phần phi tuyến có bậc cao hơn, sau đó thành phần phi tuyến này được thay thế bởi thành phần phi tuyến bậc thấp hơn cùng bậc với thành phần phi tuyến ban đầu, rồi mới thay thế thành phần phi tuyến sau cùng bởi một thành phần tuyến tính. Phương pháp này sau đó được mở rộng bởi các tác giả Elishakoff, Andrimasy, Dolley [15]. Các tác giả đó đã thực hiện điều chỉnh số bước thay thế so với cách làm như ban đầu [12]. Kết quả là đối với một số hệ phi tuyến, việc thay đổi số bước thay thế như vậy dẫn đến sai số của các đáp ứng của hệ giảm đi đáng kể. Cho đến nay phương pháp mới được áp dụng cho các hệ rời rạc, còn đối với các hệ liên tục vẫn chưa có tính toán nào được thực hiện. Do đó đây là vấn đề được đặt ra trong luận văn này. Trong luận văn này, tác giả trình bày phương pháp tuyến tính hóa điểu chỉnh cho một số hệ rời rạc và một hệ liên tục điển hình là bài toán dao động của dầm Euler-Bernoulli phi tuyến chịu kích động ngoài ngẫu nhiên. Luận văn gồm 3 chương với nội dung như sau Chương 1. Tổng quan về phương pháp tuyến tính hóa tương đương. Chương này trình bày những nội dung cơ bản của phương pháp tuyến tính hóa tương đương và áp dụng phương pháp vào hai hệ một bậc tự do điển hình là hệ Duffing và hệ Van der Pol chịu kích động ngoài ngẫu nhiên ồn trắng. Chương 2. Một số mở rộng của phương pháp tuyến tính hóa tương đương. Phương pháp tuyến hóa với tên gọi “phương pháp tuyến tính hóa điều chỉnh” được trình bày trong chương này. Sau đó là một số mở rộng của phương pháp và áp dụng vào nghiên cứu một số hệ phi tuyến như Atalik-Utku, Lutes-Sarkani và Van der Pol. Chương 3. Áp dụng phương pháp tuyến tính hóa tương đương cho bài toán dao động của dầm. Chương này trình bày ứng dụng của phương pháp tuyến tính hóa điều chỉnh vào bài toán dao động của dầm Euler-Bernoulli chịu kích động ngoài ngẫu nhiên. Trong khuôn khổ luận văn, tác giả chỉ trình bày một số hệ phi tuyến điển hình với việc sử dụng phương pháp tuyến tính hóa tương đương để nghiên cứu đáp ứng của các hệ đó. Trong quá trình thực hiện luận văn này tác giả đã hết sức cố gắng nhưng chắc chắn rằng không tránh khỏi những khiếm khuyết, tác giả mong nhận được ý kiến đóng góp của quý thầy cô và các bạn.
  9. Chương 1 Tổng quan về phương pháp tuyến tính hóa tương đương 1.1. Phương pháp tuyến tính hóa tương đương cho hệ nhiều bậc tự do 1.1.1. Phương trình vi phân chuyển động Trong phần này ta xét hệ dao động nhiều bậc tự do có phương trình chuyển động được cho dưới dạng sau đây MXCX  KX  XXX,,,   Ut (1.1) T trong đó XXX,,  là các véc tơ vị trí, vận tốc và gia tốc của hệ, XXXX  1 2 n  , các ma trận vuông cấp n gồm M m , C c , K k lần lượt là ma trận ij n n ij n n ij n n khối lượng, ma trận cản và ma trận độ cứng của hệ, hàm phi tuyến T XXXXXXXXXXXX,,   ,,    ,,    ,,   là m n thành 1 2 n ột véc tơ T phần, véc tơ UUUU  1 2 n  là một quá trình ngẫu nhiên có n thành phần. Nếu hệ (1.1) không chứa thành phần phi tuyến  thì nó có dạng một hệ phương trình vi phân tuyến tính hệ số hằng số được biết trong lý thuyết dao động ngẫu nhiên tuyến tính. Vì rất khó để tìm được nghiệm chính xác của hệ phi tuyến (1.1) nên người ta sẽ tìm cách xây dựng nghiệm xấp xỉ bằng nhiều phương pháp gần đúng khác nhau [11]. Phương pháp tuyến tính hóa tương đương là một trong những phương pháp điển hình để nghiên cứu các hệ ngẫu nhiên phi tuyến. Ý tưởng của phương pháp tuyến tính hóa tương đương cho hệ nhiều bậc tự do nói chung và hệ một bậc tự do nói riêng là sự thay thế thành phần phi tuyến của hệ bằng một thành phần tuyến tính tương ứng với việc sử dụng một tiêu chuẩn tối ưu nào đó trong đó chú ý rằng kích động ngoài của hệ phi tuyến vẫn giữ cho hệ tuyến tính. Hệ tuyến tính như thế người ta gọi là hệ tuyến tính hóa tương đương. Hệ tuyến tính hóa này có một ưu điểm nổi bật là ta có thể tìm được nghiệm của nó khi biết đầy đủ các thông số của hệ. Do đó khi sử dụng hệ tuyến tính hóa tương đương, ta có thể nghiên cứu hệ phi tuyến bằng các lý thuyết đã biết của hệ tuyến tính.
  10. 6 1.1.2. Hệ tuyến tính hóa tương đương Hệ tuyến tính hóa tương đương của hệ (1.1) có dạng sau đây . MMXCCXKKXUt e e e , (1.2) trong đó MCKe,, e e là các ma trận vuông cấp n được xác định bằng cách sử dụng một tiêu chuẩn tối ưu nào đó. Có nhiều tiêu chuẩn khác nhau để tìm ra các ma trận này (xem [9]). Song một tiêu chuẩn có lẽ được sử dụng nhiều và phổ biến hơn cả là tiêu chuẩn sai số bình phương trung bình nhỏ nhất. Sai số giữa phương trình ban đầu (1.1) và phương trình tuyến tính hóa (1.2) cho bởi eMXCXKX   XXX,,   MMXe  CCX e  KKX e (1.3)  XXXMXCXKX,,.  e  e  e Các ma trận MCKe,, e e được xác định sao cho nghiệm X t của hệ tuyến tính hóa (1.2) là một xấp xỉ “tốt nhất” của nghiệm hệ phi tuyến (1.1). Vì nghiệm của hệ tuyến tính hóa phụ thuộc vào MCe, e và K e nên người ta phải thiết lập một hệ kín giữa ba ma trận này và đáp ứng X t , từ đó mới xác định các đáp ứng khác của hệ. Tiêu chuẩn sai số bình phương trung bình áp dụng cho (1.3) được biểu diễn dưới dạng: E eT e min , (1.4) e e e mij,, c ij k ij e e e e e e trong đó mij,, c ij k ij là các phần tử của các ma trận MCK,, tương ứng, ký hiệu E là toán tử kỳ vọng toán học. Điều kiện (1.4) dẫn tới các phương trình sau đây  E eT e 0 i , j 1, n , (1.5) e mij  E eT e 0 i , j 1, n , (1.6) e cij  E eT e 0 i , j 1, n . (1.7) e kij Tiêu chuẩn (1.4) có thể được viết lại là E e2 e 2 e 2 min , (1.8) 1 2 n e e e mij,, c ij k ij T trong đó e1, e 2 , , en là các phần tử của véc tơ e  e1 e 2 en  .
  11. 7 Sử dụng tính chất tuyến tính của toán tử kỳ vọng, điều kiện (1.8) dẫn tới n 2 Di min (1.9)  me,, c e k e i 1 ij ij ij với Di được định nghĩa bởi 2 2 Di E e i i 1, n , (1.10) hay 2 n 2 e e  e Di E  i-. m ijjijjijj X c X k X (1.11) j 1 2 e e e Từ (1.11), ta có thể thấy rằng đại lượng Di chỉ phụ thuộc vào mij, c ij và kij ( j 1, n ) do đó điều kiện (1.9) được thay thế bởi điều kiện đơn giản hơn sau đây [9] D2 min . (1.12) i e e e mij,, c ij k ij Tức là ta có các phương trình:  2 e Di 0 i , j 1, n , (1.13) mij  2 e Di 0 i , j 1, n , (1.14) cij  2 e Di 0 i , j 1, n . (1.15) kij Ta biến đổi vế trái của phương trình (1.13) 2  n e e  e e E i- m is X s c is X s k is X s m  ij s 1 n e e  e  2E i -  m is X s c is X s k is X s X j (1.16) s 1 n  e   e   e  2EX j  i 2 mEXX is s j cEXX is s j kEXX is s j . s 1 Do đó phương trình (1.13) tương đương với phương trình sau n e   e   e    mEXXis sj cEXX is sj kEXX is sj EX ji  . (1.17) s 1
  12. 8 Hoàn toàn tương tự, hai phương trình (1.14) và (1.15) tương đương với lần lượt hai phương trình sau đây n e   e   e    mEXXis sj cEXX is sj kEXX is sj EX ji  , (1.18) s 1 n e  e  e  mEXXis sj cEXX is sj kEXX is sj EX ji  . (1.19) s 1 Các phương trình (1.17), (1.18) và (1.19) có thể được viết dưới dạng ma trận tương ứng e  T e   T e  T  T Mi E XX C i E XX K i E XX E  i X , (1.20) e  T e   T e  T  T Mi E XX C i E XX K i E XX E  i X , (1.21) e  T e  T e T T Mi E XX C i E XX K i E XX E  i X , (1.22) e e e e e e trong đó MCKi,, i i là các véc tơ hàng thứ i của các ma trận MC, và K tương ứng. Hệ (1.20), (1.21) và (1.22) có thể viết thành dạng ma trận E XXTTT E XX E XX  Ke C e M e EXX  T EXX   T EXX   T i i i TTT   E XX E XX E XX (1.23) EXEXEX TTT    , i i i trong đó chú ý rằng ma trận đầu tiên trong vế trái của (1.23) có cỡ là 1 3n , ma trận thứ hai có cỡ là 3n 3 n , vế phải có cỡ là 1 3n . Lấy chuyển vị cả hai vế của (1.23) ta được T TT T T TTT   T E XX E XX E XX eT EX i Ki TTTTT T T eT T EXX  EXX   EXX   C EX   . i i (1.24) eT TTTT M EXX TTTT EXX  EXX  i EX   i Biến đổi các ma trận chuyển vị ta được
  13. 9 TTT  E XX E XX E XX eT EX  Ki i T   T   T eT  EXX EXX EXX Ci E  i X . eT (1.25) TTT M EXX  EXX  EXX  i EX   i Ta đưa vào ký hiệu véc tơ X   XX , (1.26) X trong đó chú ý rằng véc tơ X có dạng ma trận cỡ 3n 1. Chuyển vị của véc tơ X T  TTT XXXX   . (1.27) T Dễ thấy rằng ma trận EXX   chính là ma trận thứ nhất trong vế phải của hệ (1.25), nên hệ (1.25) được viết lại là eT Ki T EXXCEX   eT   i i (i 1, n ). (1.28) eT M i Trong hệ (1.28) với mỗi chỉ số i ta có 3n phương trình với 3n ẩn chưa biết e e e 2 2 MCKi,, i i . Do đó hệ này tất cả có 3n phương trình với 3n ẩn chưa biết là các thành T phần của các ma trận MCe, e và K e . Ma trận EXX   là ma trận vuông cấp 3n . Để T hệ (1.28) giải được với các ẩn MCKe,, e e thì điều kiện cần và đủ là ma trận EXX   i i i không suy biến, nghĩa là T det EXX  0. (1.29) T Bây giờ ta chỉ ra rằng ma trận EXX   suy biến khi và chỉ khi ít nhất một thành phần của X có thể biểu diễn bằng một tổ hợp tuyến tính qua các thành phần còn lại. Để chứng minh điều kiện đủ, ta chú rằng, sự phụ thuộc tuyến tính tương đương với sự tồn tại 3n số thực ak,, b k c k ( k 1, n ) không đồng thời bằng không sao cho n    ak X k b k X k c k X k 0. (1.30) k 1
  14. 10 Sử dụng ký hiệu véc tơ, phương trình (1.30) được viết lại là T X v 0, (1.31) trong đó T v  a1 a 2 an b 1 b 2 b n c 1 c 2 c n  0. (1.32) Nhân trái hai vế của (1.31) và áp dụng toán tử kỳ vọng vào hai vế, đồng thời chú ý tính chất tuyến tính của toán tử này ta được T E  X X v 0. (1.33) T Vì véc tơ v khác không nên từ phương trình (1.33) ta suy ra ma trận EXX   là một ma trận suy biến. T Để chứng minh điều kiện cần, ta giả sử rằng ma trận EXX   là một ma trận suy biến. Khí đó tồn tại một véc tơ v khác không sao cho phương trình (1.33) thỏa mãn. Nhân trái cả hai vế của phương trình (1.33) với vT , ta được T vT E  X X v 0. (1.34) T Sử dụng tính chất tuyến tính của toán tử kỳ vọng và chú ý đại lượng X v là một vô hướng, ta biến đổi vế trái của (1.34) TTT 2 TT vEXX  v EvXXv  E  Xv . (1.35) Từ (1.34) và (1.35) suy ra rằng T 2 E  X v 0. (1.36) Mà E0 0 , nên phương trình (1.36) dẫn tới phương trình (1.31). Điều đó có nghĩa các phần tử của véc tơ X là phụ thuộc tuyến tính. Ta có điều phải chứng minh. Trong [9] đã chỉ ra rằng các phần tử của các ma trận MCKe,, e e có thể được xác định như sau
  15. 11  me E i , (1.37) ij  X j  ce E i , (1.38) ij  X j e i kij E . (1.39) X j Các công thức (1.37), (1.38) và (1.39) được sử dụng đầu tiên cho các bài toán dao động ngẫu nhiên dừng bởi Atalik và Utku năm 1976 (xem [7]), sau đó Spanos [19] chỉ ra rằng nó vẫn có thể áp dụng được cho bài toán dao động ngẫu nhiên không dừng. 1.1.3. Ma trận mật độ phổ Để đơn giản, ta xét hệ (1.1) mà thành phần phi tuyến  không chứa tọa độ X mà chỉ phụ thuộc vào tọa độ của vị trí và vận tốc, nghĩa là MX CX  KX  X,, X  U t (1.40) trong đó U t là quá trình ngẫu nhiên trung bình không với ma trận mật độ phổ SU  SSSUUUUUU    1 1 1 2 1 n SU  (1.41) SSS    UUUUUUn1 n 2 n n v i S  là ph chéo c a các thành ph n UU, (i, j 1, n ). ớ UUi j ổ ủ ầ i j Từ hệ tuyến tính hóa tương đương (1.2), hàm truyền của hệ này là 1 e e 2 H  K K i  C C  M . (1.42) Ma trận HH  là ma trận vuông cấp n với điều kiện không suy biến, ij n n nghĩa là e e 2 det K K i C C  M 0. (1.43)
  16. 12 Các thành phần Hij  liên quan đến các xung đơn vị hij t bởi phép biến đổi Fourier sau H h t e i t dt, ij ij (1.44) 1 h t H ei t d . ij ij 2 Nghiệm Xi t của (1.2) có thể biểu diễn qua xung đơn vị hij t và kích động ngoài U t của hệ bởi tích phân Duhamel như sau X h t  U  d ,. X h t  U  d  (1.45) i ik 1 k 1 1 j jp 2 p 2 1 Từ đó ta tính được mô men bậc hai EXX ht  ht  EU  U  dd   . (1.46) i j ik 1 jp 2 k 1 p 2 1 2 Trong khi đó mô men EUU k 1 p  2 liên quan đến mật độ phổ SU  của kích động ngoài U t như sau E U  U  S  e  2  1 d . (1.47) k1 p 2 Uk U p Thay (1.47) vào (1.46) và đổi thứ tự lấy tích phân ta được  EXX htht   S  e  2  1 ddd    i j ik1 jp 2 Uk U p  1 2    h t  e i  t 1 d  h t  e i  t  2 d  S  d  (1.48) ik1 1  jp 2 1  Uk U p   H  H  S  d . ik jp Uk U p Biểu thức (1.48) có thể viết lại dưới dạng ma trận như sau E XXTT H  S  H  d . (1.49) U
  17. 13 Từ (1.49) ta có thể chỉ ra rằng ma trận phổ đầu ra SX  của đáp ứng X t của (1.2) có dạng T SHSHXU     . (1.50) Khi đó ta có thể viết mô men bậc hai (1.49) như sau E X X S  d , (1.51) i j Xi X j S  là các thành ph n c a ma tr n S  nh t (1.50). trong đó XXi j ầ ủ ậ X được xác đị ừ Như vậy từ (1.28) và (1.51) ta đã thiết lập được một hệ phương trình gồm các ẩn là các phần tử của các ma trận KCe, e . Nói chung, hệ phương trình này là một hệ phương trình đại số phi tuyến. Để đơn giản trong việc tính toán ta xét trường hợp riêng là hệ một bậc tự do chịu kích động ngoài ngẫu nhiên. 1.2. Phương pháp tuyến tính hóa tương cho hệ một bậc tự do Hệ một bậc tự do được xem là trường hợp riêng của hệ nhiều bậc tự do. Tuy nhiên trường hợp hệ một bậc tự lại rất quan trọng và rất được quan tâm nghiên cứu. Trong phần này ta đề cập tới dao động của hệ một bậc tự chịu kích động ngẫu nhiên ồn trắng có phương trình như sau 2 x 2 hx  0 x g x , x  f t , (1.52) trong đó hàm phi tuyến g x, x phụ thuộc vào hai đối số là tọa độ của vị trí x và vận tốc x , hàm f t là một kích động ngoài ngẫu nhiên ồn trắng Gaussian có trung bình không. Dựa vào (1.2), phương trình tuyến tính hóa tương đương của phương trình phi tuyến (1.52) có dạng đơn giản hơn như sau 2 x 2 h b x  0 k x f t , (1.53) trong đó hai hệ số b, k được tìm từ tiêu chuẩn tuyến tính hóa (1.12). Cụ thể trong trường hợp này ta có tiêu chuẩn sau đây
  18. 14 2 e1 E g x, x bx  kx min. (1.54) b, k Tiêu chuẩn (1.54) được gọi là tiêu chuẩn sai số bình phương trung bình đối với hệ một bậc tự do. Nó dẫn đến hệ hai phương trình với hai ẩn b, k e 1 0, (1.55) b e 1 0. (1.56) k Giải hệ (1.55)-(1.56) với chú ý rằng E xx 0 ta được E xg  b , 2 (1.57) E x E xg k . 2 (1.58) E x Từ phương trình tuyến tính hóa, các quá trình x và x là các quá trình chuẩn. Do đó các mô men bậc lẻ đều bằng không, các mô men E xg  và E xg đều biểu diễn qua được các mô men bậc hai của x và x . Hệ (1.57) và (1.58) lập thành một hệ hai 2 2 phương trình nhưng lại có bốn ẩn bao gồm b,,, k E x E x . Để đóng kín hệ này ta cần tới hai phương trình nữa. Dựa vào ma trận của phổ đầu ra được trình bày trong mục 1.1.3, áp dụng (1.50) và (1.51) cho hệ một bậc tự do, ta có hai phương trình sau đây để đóng kín hệ (1.57)-(1.58) S  d  E x2 f , (1.59) 2 2 2 2 2 2h b   0 k 2S  d  E x2 f , (1.60) 2 2 2 2 2 2h b   0 k trong đó S f  là hàm mật độ phổ của kích động ngoài f t . Với kích động ồn trắng ta có SSf  0 const . Do đó, sử dụng định lý giá trị thặng dư trong lý thuyết hàm biến phức ta có thể tính được hai tích phân suy rộng (1.59) và (1.60) dưới dạng kết quả sau 2 S E x2 0 , (1.61) 2 2 2h b 0 k 2 S E x2 0 . (1.62) 2 2h b
  19. 15 1.3. Áp dụng phương pháp tuyến tính hóa tương đương cho các hệ Duffing và Van der Pol 1.3.1. Hệ Duffing Ta xét hệ Duffing một bậc tự do có phương trình chuyển động được mô tả bởi 2 3 x 2 hx  0 x  x  t , (1.63) trong đó hàm phi tuyến là hàm bậc ba đối với tọa độ x với tham số phi tuyến  được giả sử là dương. Kích động ngoài  t là kích động ồn trắng cường độ đơn vị với hàm tương quan R t cho bởi R  E  t  t    , (1.64) với   là hàm Delta-Dirac. Áp dụng (1.57) và (1.58) với g  x3 , các hệ số b, k của phương trình tuyến tính hóa được xác định là b 0, (1.65) E x4 k  , 2 (1.66) E x trong đó chú ý rằng x, x là các quá trình Gaussian độc lập, hay các quá trình chuẩn. Ta có công thức sau đây rút ra từ quá trình chuẩn x n 2n 2 E x 2 n 1 !! E x ( n 1,2,3 ). (1.67) Sử dụng (1.67), biểu thức của k trong (1.66) có dạng 2 k 3 E x . (1.68) Từ biểu thức (1.68), phương trình tuyến tính hóa của (1.63) là 2 2 x 2 hx  0 3  E x x  t . (1.69)
  20. 16 Đáp ứng bình phương trung bình của x được tìm từ quan hệ (1.61) như sau (xem [9])  2 E x2 . (1.70) 2 2 2.2h 0 3  E x Nếu nói đến mật độ phổ hằng số S0 của kích động đầu vào thì ta chú ý quan hệ giữa 2  và S0 cho bởi 2  2 S0 . (1.71) 2 Giải phương trình (1.70) với ẩn E x ta thu được 1 3 2 E x2  2  4 . (1.72) 0 0 6 h Nghiệm (1.72) là một nghiệm xấp xỉ của đáp ứng trung bình phương của hệ phi tuyến gốc (1.63). Để đánh giá độ chính xác của nghiệm xấp xỉ này, ta sẽ so sánh giá trị của nó khi biết các tham số đầu vào ,,, 0 h  với nghiệm chính xác của (1.63) (xem [11]). Hàm mật độ chính xác của (1.63) là 4h 12 2 1 4  p x Aexp 2 0 x  x  , (1.73)  2 4  trong đó hằng số A được tìm từ điều kiện chuẩn hóa của hàm p x 4h 1 1  A 1 exp  2 x 2  x 4 dx . (1.74) 2 0   2 4  Đáp ứng bình phương trung bình chính xác là 2 4h 1 2 2 1 4  xexp 2 0 x  x   2 4 E x2  . (1.75) cx 4h 12 2 1 4  exp x  x dx 2 0   2 4  Kết quả tính toán số đáp ứng bình phương trung bình (1.72) và nghiệm chính xác (1.75) theo sự thay đổi của tham số phi tuyến  được trình bày trong Bảng 1.
  21. 17 2 Bảng 1. Đáp ứng bình phương trung bình của hệ Duffing với các tham số h 0.5,0 1,  2  E x2 E x2 cx kd Sai số (%) 0.1 0.8176 0.5054 1.4876 0.5 0.5792 0.5486 5.2861 1.0 0.4679 0.4343 7.1938 5.0 0.2543 0.2270 10.7384 10 0.1889 0.1667 11.7721 50 0.0904 0.0784 13.2721 100 0.0650 0.0561 13.6491 Quan sát bảng ta nhận thấy rằng khi tham số  ( 0.1) nhỏ thì phương pháp tuyến tính hóa tương đương (phương pháp kinh điển) cho kết quả của đáp ứng bình phương trung bình E x2 tính theo (1.72) có sai số so với nghiệm chính xác E x2 tính kd cx theo (1.75) là nhỏ (1.4876%). Khi tăng tham số  thì sai số tăng dần, nhất là khi  100 thì kết quả của phương pháp tuyến tính hóa kinh điển cho sai số lên tới 13.6491% so với nghiệm chính xác. 1.3.2. Hệ Van der Pol Xét phương trình vi phân có dạng sau đây 2 2 x  x x  0 x  t , (1.76) trong đó ,,,  0  là các hằng số dương cho trước, hàm  t là quá trình ngẫu nhiên ồn trắng trung bình không, cường độ đơn vị và có hàm tương quan dưới dạng (1.64). . Áp dụng (1.57)-(1.58) cho hàm g  x2 x ta thu được các hệ số tuyến tính hóa của phương trình tuyến tính hóa của phương trình phi tuyến (1.76) là 2 b  E x , (1.77) k 0. (1.78)
  22. 18 Do đó phương trình tuyến tính hóa của (1.76) có dạng 2 2 x  E x x  0 x  t . (1.79) 2 Từ phương trình (1.79), đáp ứng bình phương trung bình E x được tìm từ quan hệ sau (xem [9])  2 E x2 . (1.80) 2 2 2 E x 0 2 Giải phương trình (1.80) với ẩn E x ta được 1 2 2 E x2 2 . (1.81) 2  2 0 Đây là một nghiệm xấp xỉ của đáp ứng bình phương trung bình của hệ Van der Pol một bậc tự do. Khi biết các tham số đầu vào, ta hoàn toàn tính được đáp ứng bình phương trung bình này. Ta có thể quan sát sai số tương đối của nghiệm xấp xỉ này so với nghiệm mô phỏng theo Monte-Carlo khi kích động ngoài là ngẫu nhiên ồn trắng. Các nghiệm tính toán bằng số được trình bày trong Bảng 2 với sự thay đổi của  2 . Bảng 2. Đáp ứng bình phương trung bình của hệ Van der Pol với các tham số 0.2, 0 1,  2 2 E x2 E x2  MC kd error (%) 0.02 0.2081 0.1366 34.3573 0.20 0.3638 0.2791 23.2741 1.00 0.7325 0.5525 24.5742 2.00 1.0255 0.7589 25.9998 4.00 1.4525 1.0512 27.6248 Quan sát Bảng 2 ta thấy rằng sai số của đáp ứng bình phương trung bình E x2 của kd phương pháp kinh điển tính theo công thức (1.81) so với đáp ứng bình phương trung bình E x2 thu được từ mô phỏng Monte-Carlo là khá lớn, từ 23.2741% đến MC 34.3573% khi  2 thay đổi từ 0.02 đến 4.00.
  23. 19 Tổng kết chương 1 Trên đây, phương pháp tuyến tính hóa tương đương được trình bày với nội dung của tiêu chuẩn tuyến tính hóa kinh điển, giả thiết Gaussian ồn trắng. Với hai hệ điển hình được minh họa, phương pháp tuyến tính hóa tương đương cho kết quả nghiệm giải tích gần đúng tính đáp ứng bình phương trung bình của hệ. Tuy nhiên cũng dễ nhận thấy rằng khi hệ phi tuyến yếu thì nghiệm giải tích này cho một kết quả có thể nói là tốt, nhưng khi tính phi tuyến của hệ tăng lên thì sai số tăng lên. Do đó đòi hỏi chúng ta phải cải tiến phương pháp để đạt được sai số nhỏ hơn theo mong muốn. Để làm điều đó, ta có thể sử dụng các tiêu chuẩn tuyến tính hóa khác nhau, như được giới thiệu trong [11]. Trong nội dung của chương sau, tác giả trình bày một cách tiếp cận khác đối với bài toán tuyến tính hóa tương đương nhằm tăng độ chính xác của đáp ứng của hệ phi tuyến.
  24. Chương 2 Một số mở rộng của phương pháp tuyến tính hóa tương đương Như đã nói ở chương 1, trong một số hệ phi tuyến, phương pháp tuyến tính hóa tương đương Gaussian có thể cho ta sai số lớn khi tính phi tuyến của hệ tăng lên. Để giảm sai số, trong chương này, tác giả trình bày một cách tiếp cận khác của bài toán tuyến tính hóa tương đương với tên gọi “Phương pháp tuyến tính hóa điều chỉnh” [12]. Sau đó là một số mở rộng của phương pháp này. 2.1. Phương pháp tuyến tính hóa điều chỉnh Xét phương trình phương trình vi phân chuyển động của hệ một bậc tự do chịu kích động ngoài ngẫu nhiên Gaussian f t được cho dưới dạng sau đây 2 z 2 hz  0 z g z ,z  f t , (2.1) trong đó các hằng số dương h,0 cho trước, hàm phi tuyến g z, z được giả thiết có dạng đa thức của các đối số z và z n n 2p 2 q 1 2 p 2 q 1 g z,, z  pqz  z  pq zz  (2.2) p 0 q 0 ( pq,  pq ( p, q 0,1,2, ) là các hằng số). Theo phương pháp tuyến tính hóa tương đương kinh điển, thì các số hạng phi tuyến của hàm g z, z được thay thế bởi các số hạng tuyến tính 2p 2 q 1 pqz z pq z, (2.3) 2p 2 q 1 pqz z  pq z , (2.4) trong đó pq,  pq được tìm từ tiêu chuẩn sai số bình phương trung bình nhỏ nhất 2p 2 q 1 2 E pq z z pq z min, (2.5) pq
  25. 21 2p 2 q 1 2 E pq z z  pq z  min. (2.6)  pq Tuy nhiên như đã chỉ ra [12], trong một số hệ phi tuyến, phương pháp tuyến tính hóa tương đương kinh điển cho kết quả tốt đối với các hệ phi tuyến bé, còn đối với các hệ phi tuyến mạnh thì phương pháp này lại cho sai số lớn. Điều này đòi hỏi cần phải cải tiến phương pháp theo một cách nào mà sai số của các hệ phi tuyến mạnh có thể giảm đi. Theo tinh thần này, ta sẽ tuyến tính hóa các số hạng phi tuyến của g z, z bằng một con đường khác với phương pháp tuyến tính hóa kinh điển. Trước tiên số hạng phi tuyến được thay thế bằng một số hạng phi tuyến bậc cao hơn, số hạng phi tuyến bậc cao này lại được thay thế bởi số hạng phi tuyến cùng bậc với số hạng phi tuyến ban 2p 2 q 1 đầu, cuối cùng ta mới thay thế bởi một số hạng tuyến tính. Cụ thể ta có pq z z 1 2p 2 q 1 2 p 2 q 1 4 p 4 q 1 được thay thế bởi pq z z z  z pq z  z , sau đó số hạng này được thay thế 2 2p 2 q 1 2 2p 2 q 1 3 bởi pq z z , cuối cùng là số hạng phi tuyến pq z z được thay thế bởi pq z : 2p 2 q 1 1 4 p 4 q 1 2 2 p 2 q 1 3 pqz z pqz  z pqz  z pq z. (2.7) Tương tự 2p 2 q 1 1 4 p 4 q 1 2 2 p 2 q 1 3 pqz z  pq zz   pq zz   pq z , (2.8) k k trong đó các hệ số pq,  pq k 1,2,3 được tìm từ tiêu chuẩn sai số bình phương trung bình áp dụng cho mỗi bước thay thế. Tức là với mỗi bước thay thế, sai số bình phương trung bình giữa các số hạng phi tuyến được cực tiểu hóa theo các tham số tương ứng. k Với sơ đồ (2.7) ta có các tiêu chuẩn sau đây để xác định các hệ số pq : 2 2p 2 q 1 1 4 p 4 q 1 E z z z  z min, (2.9) pq pq 1 pq 2 1 4p 4 q 1 2 2 p 2 q 1 E z z z  z min, (2.10) pq pq 2 pq 2 2 2p 2 q 1 3 E z z z min. (2.11) pq pq 3 pq 1 2 3 Các tiêu chuẩn (2.9), (2.10), (2.11) dẫn đến đạo hàm riêng theo pq , pq , pq của các biểu thức ở vế phải tương ứng bằng không  2 E z2p z 2 q 1 1 z 4 p z 4 q 1 0, (2.12) 1 pq pq   pq
  26. 22  2 E 1 z4p z 4 q 1 2 z 2 p z 2 q 1 0, (2.13) 2 pq pq   pq  2 E 2 z2p z 2 q 1 3 z 0. (2.14) 3 pq pq  pq Các phương trình (2.12), (2.13), (2.14) lần lượt tương đương với E z6p z 6 q 2 1 , pq pq 8p 8 q 2 (2.15) E z z E z6p z 6 q 2 2 1 , pq pq 4p 4 q 2 (2.16) E z z E z2p z 2 q 2 3 2 . pq pq 2 (2.17) E z 3 Thay (2.15) vào (2.16), lấy kết quả thu được thay vào (2.17) ta được biểu thức của pq như sau E z2p z 2 q 2 E z  6 p z 6 q 2 E z  6 p z 6 q 2 3 . pq pq 2 4p 4 q 2 8 p 8 q 2 (2.18) E z E z z E z  z k Bây giờ ta xác định các hệ số  pq k 1,2,3 của sơ đồ (2.8). Các tiêu chuẩn sai số bình phương trung bình tương ứng với sơ đồ này là 2 2p 2 q 1 1 4 p 4 q 1 E z z  zz  min, (2.19) pq pq 1  pq 2 1 4p 4 q 1 2 2 p 2 q 1 E z z  zz  min, (2.20) pq pq 2  pq 2 2 2p 2 q 1 3 E z z  z  min. (2.21) pq pq 3  pq k Tương tự như việc tính toán đối với pq k 1,2,3 , ta thu được các biểu thức của k  pq k 1,2,3 là E z6p z 6 q 2  1  , pq pq 8p 8 q 2 (2.22) E z z
  27. 23 E z6p z 6 q 2  2  1 , pq pq 4p 4 q 2 (2.23) E z z E z2p z 2 q 2  3  2 . pq pq 2 (2.24) E z 3 Thay (2.22) và (2.23) vào (2.24) ta được biểu thức của  pq E z2p z 2 q 2 E z 6 pz  6 q 2 E z 6 pz  6 q 2  3  . pq pq 2 4p 4 q 2 8 p 8 q 2 (2.25) E z E zz  E zz  3 3 Từ các sơ đồ tuyến tính hóa (2.7), (2.8) và các biểu thức của pq ,  pq , hàm phi tuyến g z, z được thay thế bởi hàm tuyến tính của z và z n n 3 3 g z,. z  pq z  pq z  (2.26) p 0 q 0 Như vậy phương trình tuyến tính hóa của (2.1) có dạng n n 2 3 3 z 2 hz  0 z  pq z  pq z  f t , (2.27) p 0 q 0 3 3 trong đó pq ,  pq được xác định từ các biểu thức (2.18), (2.25). Từ phương trình (2.27) ta có thể xác định được một cách gần đúng đáp ứng bình phương trung bình của hệ phi tuyến (2.1) ban đầu. Các đáp ứng z, z là các quá trình Gaussian độc lập nên tất cả các mô men bậc cao đều biểu diễn được thông qua mô men bậc một và bậc hai của hệ. Để minh họa cho phương pháp tuyến tính hóa điều chỉnh ở trên, ta áp dụng phương pháp này cho các hệ Atalik-Utku, Lutes-Sarkani và Van der Pol. 2.2. Áp dụng phương pháp tuyến tính hóa điều chỉnh cho các hệ Atalik-Utku, Lutes-Sarkani và Van der Pol 2.2.1. Hệ Atalik-Utku Xét hệ phi tuyến sau z 2 hz   z3 f t , (2.28)
  28. 24 trong đó h, là các hằng số dương, kích động ngoài f t là quá trình Gaussian ồn trắng, trung bình không và có hàm tương quan Rf  2 Rf  E f t f t     , (2.29) với   là hàm Delta-Dirac. Hàm mật độ xác suất dừng chính xác của hệ (2.28) là (xem [15]) h 4  p z Aexp z  , (2.30)  2  trong đó A được tìm từ điều kiện chuẩn hóa hàm p z h  A 1 exp z 4 dz . (2.31) 2    Biểu thức chính xác của đáp ứng bình phương trung bình của hệ (2.28) là 2 h 4  zexp z  dz  2  E z2 z 2 p z dz . (2.32) cx h 4  exp z dz 2    Ta chú ý công thức tích phân s s 1 r 1 s uexp au du a r  , (2.33) 0 r r trong đó  v là ký hiệu hàm Gamma  v uv 1 exp u du . (2.34) 0 Sử dụng (2.33) để biến đổi (2.32): 3 1 h 4 3 3 2  2  2 4 4 4  E z2 . 0.337989 . 1 (2.35) cx h 1 h  1 h 4 1  2  4 4  4
  29. 25 Phương trình (2.28) có hàm phi tuyến g  z3 chỉ phụ thuộc vào z , nên áp dụng sơ đồ tuyến tính hóa (2.7) ta có sơ đồ sau đây cho  z3 3 5 3 z 1 z 2 z 3 z. (2.36) Sử dụng công thức (2.18) và chú ý (1.67) ta thu được 7  7 z3 z 5 z 3  E z 2 z. (2.37) 2 9E z 9 3 Phương trình tuyến tính hóa tương đương của (2.28) sử dụng tuyến tính hóa điều chỉnh: 7 z 2 hz   E z2 z f t . (2.38) 3 2 Đáp ứng bình phương trung bình E z được tìm từ  2 E z2 . (2.39) 7 4h E z2 3 2 Giải phương trình (2.39) với ẩn E z ta được 3 2  2 E z2 0.327327 . (2.40) 28 h h  Phương trình tuyến tính hóa tương đương của (2.28) sử dụng phương pháp tuyến tính hóa kinh điển có dạng 2 z 2 hz  3 E z z f t . (2.41) Đáp ứng bình phương trung bình ứng với phương pháp kinh điển 1 2  2 E z2 0.288675 . (2.42) kd 2 3 h h 
  30. 26 Sai số tương đối giữa đáp ứng bình phương trung bình của phương pháp tuyến tính hóa điều chỉnh E z2 (biểu thức (2.40) ) và đáp ứng bình phương trung bình chính dc xác E z2 (biểu thức (2.35)) là cx 2 2 E z E z dc cx 100% 3.15%. (2.43) 1 E z2 cx Sai số tương đối giữa đáp ứng bình phương trung bình của phương pháp kinh điển và đáp ứng bình phương trung bình chính xác là 2 2 E z E z kd cx 100% 14.59%. (2.44) 2 E z2 cx So sánh (2.43) và (2.44) ta thấy đối với phương pháp tuyến tính hóa điều chỉnh, sai số là 3.15% nhỏ hơn xấp xỉ 4.6 lần so với phương pháp tuyến tính hóa kinh điển (14.59%). Như vậy bằng việc sử dụng phương pháp tuyến tính hóa điều chỉnh kết quả của đáp ứng bình phương trung bình của hệ Atalik-Utku cho kết quả gần nghiệm chính xác hơn rất nhiều so với phương pháp tuyến tính hóa kinh điển. 2.2.2. Hệ Lutes-Sarkani Xét hệ phương trình phi tuyến sau a z  zsgn z f t , (2.45) trong đó a là một số nguyên dương,  là hằng số dương cho trước, hàm f t là quá trình Gaussian dừng, ồn trắng với hàm mật độ phổ hằng số S0 . Hàm phi tuyến a g  zsgn z được tuyến tính hóa thành keq z k eq zsgn z . Số mũ của hàm phi tuyến ban đầu giảm một lượng là a 1, do đó trong sơ đồ tuyến tính hóa điều chỉnh, số a hạng phi tuyến ban đầu  zsgn z được thay thế bởi số hạng phi tuyến có bậc lớn a a 1 2 a 1 một lượng là a 1, nghĩa là k1 zsgn z k 1 z sgn z , sau đó số hạng này được a a thay thế bởi k2 zsgn z , cuối cùng số hạng phi tuyến k2 zsgn z được thay thế bởi số hạng tuyến tính keq z . Sơ đồ tuyến tính hóa theo mô tả ở trên như sau
  31. 27 a2 a 1 a  zsgn zkz 1 sgn zkz 2 sgn zkz eq sgn zkz eq . (2.46) Với mỗi bước thay thế, tiêu chuẩn sai số bình phương trung bình được sử dụng để tìm các hệ số k1,, k 2 keq . Cụ thể là 2 E  zasgn z k z2 a 1 sgn z min, (2.47) 1 k1 2 E k z2a 1 sgn z k z a sgn z min, (2.48) 1 2 k2 2 E k za sgn z k z sgn z min. (2.49) 2 eq keq Các tiêu chuẩn (2.47), (2.48), (2.49) dẫn tới lần lượt các đạo hàm riêng theo k1,, k 2 keq :  a2 a 1 2 E  zsgn z k z sgn z 0, (2.50) 1 k1  2a 1 a 2 E k zsgn z k z sgn z 0, (2.51) 1 2 k2  a 2 E k zsgn z k z sgn z 0. (2.52) 2 eq keq Các phương trình (2.50), (2.51), (2.52) lần lượt dẫn tới các biểu thức sau đây cho k1,, k 2 keq : E z 3a 1 k1  , (2.53) E z 4a 2 E z 3a 1 k2 k 1 , (2.54) E z 2a E z a 1 keq k2 . (2.55) E z 2 Thay (2.53) và (2.54) vào (2.55) ta được E za 1 E z 3 a 1 E z 3 a 1 keq  . (2.56) E z2 E z 2a E z 4 a 2
  32. 28 Bi u th c t ng quát cho E z a là (xem [15]) ể ứ ổ a 2 a u u  E z  a exp du , (2.57) z  2 2  2 2 trong đó  z E z là ký hiệu phương sai của z (do kỳ vọng của z bằng không). Sử dụng công thức tích phân (2.33) ta biến đổi (2.57): a 2 a 2 a2 a u  2  u  E z z uexp du z ua exp du   2 0 2  2 0 2  a 1 (2.58) a a a 1 2 2z1 1 a 1  z a 1  22  . 2 2 2 2 2 2 Thay biểu thức (2.58) vào (2.56) ta được 2 a 2 3 a a 1   a 1 2 2 k  22 . (2.59) eq z 3 2a 1 4 a 1    2 2 2 Từ phương trình tuyến tính hóa của (2.45) và chú ý (2.59), đáp ứng bình phương trung bình của z được tìm từ quan hệ 2 S0  z . (2.60) keq Thay (2.59) vào (2.60) và giải phương trình với ẩn  z ta thu được 2 2 a 1 1 2 a 1 a 1  2 S0 2 2a 1 4 a 1 a 2 3 a z 2      . (2.61)  2 2 2 2  Biểu thức (2.61) là biểu thức xác định phương sai của đáp ứng của hệ Lutes-Sarkani theo phương pháp tuyến tính hóa điều chỉnh. Bây giờ ta tìm biểu thức phương sai chính xác của hệ. Hàm mật độ xác suất chính xác của đáp ứng của hệ (2.45) cho bởi [16]  za 1  p z Aexp  , (2.62) a 1 S0 
  33. 29 trong đó hằng số A được tìm từ điều kiện chuẩn hóa hàm p z a 1 1  u  A exp  du . (2.63) a 1 S0  Sử dụng công thức tích phân (2.33) ta tính được 1 1  a 1 a 1 1 A  . (2.64) a 1 S0 2 a 1 Phương sai chính xác:  ua 1   2 u 2 p u du 2 A u 2 exp z, cx  0 a 1 S0  2 (2.65) 1 a 1 2 S0 3 1 a 1 a 1   .  a 1 a 1 Ngoài ra sử dụng phương pháp tuyến tính hóa kinh điển ta thu được biểu thức của phương sai như sau 2 2 a 1 S a 1 2  2 0 . (2.66) z,kd a  2 2a a / 2 Bảng 3. Phương sai của hệ Lutes-Sarkani với các giá trị khác nhau của a 2 2 2 2 a  z, cx  z, kd 1 (%)  z, I 2 (%)  z, II 3 (%) 1 3.1416 3.1416 0.0000 3.1416 0.0000 3.1416 0 2 1.6655 1.5708 5.6877 1.6784 0.7713 1.7599 5.6693 3 1.1981 1.0233 14.5904 1.1603 3.1546 1.2614 5.2820 4 0.9761 0.7531 22.8490 0.8884 8.9861 0.9881 1.2229 5 0.8474 0.5939 29.9225 0.7201 15.0206 0.8134 4.0131 6 0.7635 0.4894 35.8981 0.6056 20.6846 0.6917 9.4038 7 0.7045 0.4159 40.9630 0.5226 25.8224 0.6019 14.5548 Quan sát Bảng 3 ta thấy rằng sai số 1 của phương pháp kinh điển so với nghiệm chính xác là khá lớn khi số mũ a của thành phần phi tuyến tăng lên, tuy nhiên khi sử 2 dụng phương pháp tuyến tính hóa điều chỉnh thì sai số 2 của phương sai  z, I tính
  34. 30 theo biểu thức (2.62) thì sai số có giảm so với tuyến tính hóa kinh điển. Đặc biệt khi a 2 thì sai số 2 khá nhỏ (0.7713%). 2.2.3. Hệ Van der Pol Ta xét hệ Van der Pol có phương trình (1.76), hàm phi tuyến g  x2 x được tuyến tính hóa theo sơ đồ (2.8) với p 1, q 0 như sau 2 4 2 x x 1 x x   2 x x   3 x . (2.67) Áp dụng kết quả (2.25) ta thu được hệ số 3 : E x2 x 2 E x 6 x  2 E x 6 x  2   . 3 2 4 2 8 2 (2.68) E x E x x  E x x  Các đáp ứng x, x là các quá trình Gaussian độc lập nên sử dụng công thức (1.67) ta có 2 2 2 2 E x x E x E x  , 2 4 2 2 2 E x x 3 E x E x  , 3 (2.69) 6 2 2 2 E x x 15 E x E x  , 4 8 2 2 2 E x x 105 E x E x  . Thay (2.69) vào (2.68), sau khi rút gọn ta được 5   E x2 . (2.70) 3 7 Phương trình tuyến tính hóa của hệ Van der Pol (1.76): 5 2 2 x  E x x  0 x  t . (2.71) 7 Đáp ứng bình phương trung bình của (2.71) được tìm từ quan hệ 2 2  E x . (2.72) 5 2 2 2 E x 0 7
  35. 31 2 Giải phương trình (2.72) với ẩn E x ta được 7 10 2 E x2 2 . (2.73) 10 7  2 0 Biểu thức (2.73) là giá trị của đáp ứng trung bình phương được tính theo phương pháp tuyến tính hóa điều chỉnh. Nó khác với biểu thức (1.81) tính theo phương pháp tuyến tính hóa kinh điển. Bảng 4. Đáp ứng bình phương trung bình của hệ Van der Pol với các phương pháp khác nhau ( 0.2, 0 1,  2 ) 2 2 2 2 E x E x (%) E x 2 (%)  MC kd 1 dc 0.02 0.2081 0.1366 34.3573 0.1791 13.9356 0.20 0.3638 0.2791 23.2741 0.3437 5.5250 1.00 0.7325 0.5525 24.5742 0.6657 9.1195 2.00 1.0255 0.7589 25.9998 0.9096 11.3018 4.00 1.4525 1.0512 27.6248 1.2553 13.5766 Trong phần này ta sẽ so sánh các giá trị số của đáp ứng bình phương trung bình E x2 tính theo (1.81) và E x2 (2.73) n u kd dc tính theo phương pháp tuyế tính hóa điề ch nh mô ph ng -Carlo E x2 . Quan sát B ng 4, ỉ vơí giá trị ỏ theo phương pháp Monte MC ả ta thấy sai số 2 của đáp ứng bình phương trung bình tính theo tuyến tính hóa điều 2 chỉnh đã giảm rõ rệt so với phương pháp kinh điển 1 . Đặc biệt khi  0.2 thì sai số 2 chỉ còn là xấp xỉ 5.5%. 2.3. Mở rộng của phương pháp tuyến tính hóa điều chỉnh: điều chỉnh hai bước Như đã biết, phương pháp tuyến tính hóa điều chỉnh cho kết quả tốt hơn so với phương pháp tuyến tính hóa kinh điển trong một số hệ phi tuyến. Tuy nhiên với nội dung của tuyến tính hóa điều chỉnh, một câu hỏi xuất hiện là điều gì xảy ra nếu ta thay đổi số bước điều chỉnh của phương pháp. Nghĩa là hiện nay ta mới thay thế số hạng phi tuyến của hệ bởi một số hạng phi tuyến khác nhưng với chỉ một lần, bây giờ ta sẽ
  36. 32 không dừng lại ở đó mà số hạng phi tuyến được thay thế bởi nhiều số hạng phi tuyến khác, sau đó mới tuyến tính hóa thành số hạng tuyến tính. Tuy nhiên để đơn giản cho việc trình bày, ở đây ta sẽ sử dụng hai lần thay thế số hạng phi tuyến, ta gọi đây là điều chỉnh hai bước. Phương pháp điều chỉnh hai bước này được trình bày với hai hệ Atalik-Utku và Lutes-Sarkani. 2.3.1. Phương pháp điều chỉnh hai bước cho hệ Atalik-Utku Số hạng phi tuyến của hệ Atalik-Utku (2.28) được tuyến tính hóa theo điều chỉnh hai bước với sơ đồ điều chỉnh như sau 3 5 7 5 3 z 1 z 2 z 3 z 4 z 5 z. (2.74) Trong sơ đồ (2.74), ta không dừng lại ở việc thay thế z3 bởi z5 mà ta còn tăng lên z7 rồi sau đó mới giảm dần bậc của z . Với mỗi bước thay thế ta cũng sử dụng tiêu chuẩn sai số bình phương trung bình. Tính toán tương tự như điều chỉnh một bước (2.36) (ta hiểu điều chỉnh số mũ của z ) ta thu được  11  77  77  z3 z 5 z 7 z 5 z 3 E z 2 z. (2.75) 22 2 9E z 2 117 E z 117 39 117 E z Phương trình tuyến tính hóa của (2.28) ứng với (2.75) là 77 z 2 hz   E z2 z f t . (2.76) 39 Đáp ứng bình phương trung bình thu được từ (2.76): 39 2  2 E z2 0.355842 . (2.77) 308 h h  Sai số tương đối của (2.77) so với nghiệm chính xác (2.35) là 3 5.28% , sai số này lớn hơn sai số của phương pháp điều chỉnh một bước (2.40). Do đó đối với hệ Atalik- Utku, phương pháp tuyến tính hóa điều chỉnh một bước cho kết quả tốt hơn cả. Bây giờ ta kiểm tra xem tuyến tính hóa điều chỉnh có cho kết quả mong muốn đối với hệ Lutes-Sarkani không.
  37. 33 2.3.2. Phương pháp tuyến tính hóa điều chỉnh hai bước cho hệ Lutes-Sarkani Ta sử dụng tuyến tính hóa điều chỉnh hai bước cho hệ được mô tả bởi phương trình (2.45) như sau  zasgn zkz 2 a 1 sgn zkz 3 a 2 sgn zkz 2 a 1 sgn z 1 2 3 (2.78) a k4 zsgn z keq . II z sgn z k eq . II z . Tương tự như tuyến tính hóa điều chỉnh một bước (2.46) ta thu được E za 1 E z 3 a 1 E z 5 a 3 E z 5 a 3 E z 3 a 1 keq. II  . (2.79) E z2 E z 2a E z 4 a 2 E z 6 a 4 E z 4 a 2 Sử dụng công thức (2.58) ta biến đổi (2.79) thành 2 2 a 2 3 a 5 a 2 a 1    a 1 2 2 2 2 keq. II  z 22 . (2.80) 3 2a 1 4 a 1 6 a 3     2 2 2 2 Biểu thức phương sai ứng với (2.80): 2 2 a 1 2 2a 1 4 a 1 6 a 3 a 1    a 1 2 S0 2 2 2 2 z, II 22 2  . (2.81)  a 2 3 a 5 a 2    2 2 2  Biểu thức (2.81) là giá trị xấp xỉ của phương sai của hệ Lutes-Sarkani sử dụng phương pháp tuyến tính hóa điều chỉnh hai bước. Số liệu tính toán cho (2.81) được trình bày trong Bảng 3 (cột thứ 7). Rõ ràng rằng với việc sử dụng tuyến tính hóa điều chỉnh hai bước thì sai số của đáp ứng bình phương trung bình giảm hơn so với tuyến tính hóa điều chỉnh một bước khi mà giá trị của a lớn dần. Điều đó cho thấy một lợi thế của điều chỉnh hai bước cho hệ Lutes-Sarkani so với điều chỉnh một bước. Còn đối với hệ Atalik-Utku như tính toán ở trên thì điều chỉnh hai bước lại cho kết quả xấu hơn so với một bước. Như vậy phương pháp tuyến tính hóa điều chỉnh cho một số hệ phi tuyến cho kết quả với sai số nhỏ hơn so với phương pháp tuyến tính hóa kinh điển. Đây là một ưu điểm của phương pháp khi áp dụng cho các hệ phi tuyến ở trên. Tuy nhiên ta cũng cần khảo sát nhiều hệ hơn để kiểm tra độ chính xác của phương pháp.
  38. 34 Tổng kết chương 2 Chương 2 trình bày phương pháp tuyến tính hóa điều chỉnh cho hệ có hàm phi tuyến là đa thức của các đối số vị trí và vận tốc. Sau đó là một số mở rộng của phương pháp với việc thay đổi số bước điều chỉnh. Phương pháp xuất phát từ ý tưởng rằng số hạng phi tuyến ban đầu không được tuyến tính hóa trực tiếp như thông thường mà được thay thế bởi số hạng phi tuyến có bậc cao hơn, sau đó lại thay thế bởi số hạng phi tuyến có bậc thấp hơn và cuối cùng mới tuyến tính hóa thành số hạng tuyến tính. Do đó phương pháp mở ra một cách thức mới để tuyến tính hóa một hệ phi tuyến chịu kích động ngẫu nhiên. Trong chương này các hệ một bậc tự do minh họa cho phương pháp tuyến tính hóa điều chỉnh. Kết quả thu được đã khẳng định tính ưu thế của phương pháp trong một số hệ phi tuyến so với việc sử dụng phương pháp tuyến tính hóa kinh điển. Trong chương tiếp theo phương pháp được áp dụng sang một trong những hệ liên tục thường gặp trong thực tế, đó là kết cấu dầm Euler-Bernoulli phi tuyến chịu kích động ngoài ngẫu nhiên.
  39. Chương 3 Áp dụng phương pháp tuyến tính hóa tương đương cho bài toán dao động của dầm Phương pháp tuyến tính hóa tương đương được quan tâm nghiên cứu đối với các hệ rời rạc một hay nhiều bậc tự do. Ngoài ra, nó cũng được áp dụng rộng rãi nghiên cứu nhiều hệ liên tục khác nhau. Phương pháp tuyến tính hóa điều chỉnh được nghiên cứu gần đây [12, 15] mới chỉ tập trung vào các hệ rời rạc, do đó trong nghiên cứu này, tác giả sẽ áp dụng phương pháp vào một trong những bài toán dao động thường gặp của hệ liên tục, đó là bài toán dao động của dầm. 3.1. Phương trình dao động của dầm Các bài toán của các hệ liên tục khá đa dạng và phức tạp. Trong chương này, tác giả chỉ giới hạn nghiên cứu trong phạm vi bài toán dầm Euler-Bernoulli. Ta xét bài toán dầm Euler-Bernoulli với nền đàn hồi bị ràng buộc ở hai đầu chịu tải trọng phân bố p x, t cho bởi phương trình dao động sau [18] 4w  2 w  2 w  w EI N A  K w p x,, t (3.1) x4  x 2  t 2  t f trong đó 2 EAL  w N dx (3.2) 2L0  x Các đại lượng EIAL,,,, lần lượt là mô đun đàn hồi, mô men quán tính của tiết diện ngang, mật độ khối lượng, diện tích mặt cắt ngang và chiều dài của dầm,  là hệ số cản nhớt, K f là hệ số cứng của nền đàn hồi, w x, t là độ võng của dầm, hàm p x, t là quá trình ergodic dừng đóng vai trò là tải trọng phân bố tác dụng lên dầm. Tải trọng này được giả thiết rằng có thể tách thành tích của hai thành phần không gian và thời gian như sau p x,, t f x q t (3.3)
  40. 36 trong đó f x là hàm tiền định phụ thuộc vào x , còn q t là quá trình ngẫu nhiên dừng với hàm tương quan theo Wiener-Khintchine: Eqtqt Rtt S ei t2 t 1 d  (3.4) 1 2 q 2 1 với S  là hàm mật độ phổ của q t . Giả sử độ võng w x, t được biểu diễn dưới dạng w x,, t  wn t  n x (3.5) n 1 trong đó n x ( n 1,2, ) là các hàm dạng trực giao thỏa mãn các quan hệ sau [18] d 4 EIn A 2  , (3.6) dx4 n n 1 0 khi n m x  d   ( ). (3.7) n m nm 0 1 khi n m L Thay (3.5) vào (3.2) ta có 2 EALL d EA d  d  N w tn dx w w n m dx n  n m 2L0 n 1 dx 2 L 0 n 1 m 1 dx dx EA L d d  EA 1 d  d  w wn m dx w w n m d (3.8) n m 2  n m 2Ln 1 m 10 dx dx 2 L n 1 m 1 0 d d  EA D w w , 2  nm n m 2L n 1 m 1 trong đó 1 d d  D n m d D . (3.9) nm mn 0 d d  Thay (3.5) và (3.8) vào vế phải của phương trình (3.1): d4 EA d 2  EI w tn D w w w t n n 4 2  ij i j  n 2 n 1dx 2 L i 1 j 1 n 1 dx A wn  n   w  n  n K f  w n  n (3.10) n 1 n 1 n 1 d4EA d 2  AwwKw   EIn Dww n w .  n  n f n n 4 2  ij i j 2 n n 1 dx2 L i 1 j 1 dx
  41. 37 Sử dụng (3.6) và (3.10), phương trình (3.1) trở thành EA d 2 Aw  wKw Aw 2  Dwwwn fxqt .  n  n f n n n n2  ij i j n 2 (3.11) n 12L n 1 i 1 j 1 dx Nhân cả hai vế của (3.11) với m rồi lấy tích phân hai vế trên đoạn 0, L và chú ý (3.7) ta được L Aw  w K w A 2 w   dx  n  n f n n n n m n 1 0 (3.12) EA LLd 2 D w w wn  dx f x  dx q t . 2 ij i j n 2 m m 2Ln 1 i 1 j 1 0 dx 0 Chú ý rằng L 1 0 khi n m  dx L   d  , (3.13) n m n m 0 0 Lkhi n m LLL2 x L dnd d  n d  n d  n d  m 2 mdx  m  m dx dx dx dx dx dx dx 0 0x 0 0 (3.14) 11 d d  1 0 n m d D , nm L0 d d  L trong đó các hàm n ( n 1,2, ) thỏa mãn d2 d 2    0,n n 0. n x 0 x L 2 2 (3.15) dxx 0 dx x L Thay (3.13) và (3.14) vào (3.12) rồi chia cả hai vế cho AL ta được  K E Q w ww  2 f w DDwww m qt , m  m m m m4  ij nm i j n (3.16) AALA 2 n 1 i 1 j 1 trong đó 1 L 1 Q f dx f  d . (3.17) m m m L 0 0 Trong chương này ta chỉ giới hạn tính toán cho dầm với hai đầu tựa đơn giản, kích động ngoài p x, t được giả thiết là chỉ có phần ngẫu nhiên q t còn phần tiền định f x 1. Hàm dạng n cho dầm tựa giản đơn [20] m  2 sin m  . (3.18)
  42. 38 Ta tính được các đại lượng Qm và Dnm 1 0 khi m ch½n 2 m Q 2 sin m  d  1 1 , (3.19) m 2 2 0 m khi m lÎ m 1d d  1 2m 2 khi m n D D n m d 2 nm 2 cos n  cos m  d  (3.20) nm mn 0d d  0 0 khi m n Sử dụng (3.20) ta biến đổi số hạng cuối cùng trong phương trình (3.16) như sau 2 2 DDwwwij nm i j n  DDww ii nm i n  Dw ii i  Dw nm n n 1 i 1 j 1 n 1 i 1 i 1 n 1 (3.21) 2 2 4 2 2 2 DwDwii i mm m  DwDw nn n mm m m  nww n m. i 1 n 1 n 1 Thay (3.21) vào (3.16) ta được phương trình vi phân chuyển động của dầm đàn hồi tựa giản đơn ở hai đầu như sau K E 4 m 2 Q w w  2 w f w n 2 w 2 w m q t . m  m m m m4  n m (3.22) AALA 2 n 1 Đặt EI 4  2 , (3.23) 0 AL4 I R , (3.24) A K f 2 . (3.25) A 0 Từ (3.6), (3.18) và (3.23) ta rút ra 2 2 4 m 0 m . (3.26) Sử dụng các ký hiệu (3.23)-(3.26), ta đưa phương trình (3.22) về dạng   2 Q w w  21 w m n 2 w 2 w m q t m  m m 4 m 2 2  n m . (3.27) A m 2 R mn 1 A ( m 1,2,3, )
  43. 39 Ta sẽ sử dụng phương trình (3.27) để nghiên cứu đáp ứng của hệ ban đầu cho bởi phương trình (3.1). Kích động ngoài ngẫu nhiên q t có hàm tương quan (3.4) và hàm mật độ phổ SS  0 const . Trong phần tiếp theo ta sẽ tuyến tính hóa phương trình phi tuyến (3.27) bằng việc sử dụng phương pháp tuyến tính hóa điều chỉnh. 3.2. Phương pháp tuyến tính hóa điều chỉnh cho bài toán dao động của dầm Phương trình phi tuyến (3.27) được tuyến tính hóa thành phương trình sau đây  Q w w  k 2 w m q t , (3.28) m AA m m, eq m m trong đó km, eq là hệ số tuyến tính hóa tương đương và là đại lượng không thứ nguyên được tìm từ một tiêu chuẩn tuyến tính hóa. Trong phương trình (3.27) ta sẽ tuyến tính 2 2 hóa thành phần phi tuyến  n wn w m . Do đó ta cần tuyến tính hóa các số hạng phi n 1 2 tuyến có dạng wn w m . Ta đề nghị sơ đồ tuyến tính hóa cho số hạng này như sau 2 4 2 wnm w a nmnm w w b nmnm w w c nmm w . (3.29) 2 Nghĩa là số hạng phi tuyến wn w m ban đầu được thay thế bởi số hạng phi tuyến bậc cao 4 2 hơn anm w n w m , tiếp theo ta thay thế số hạng này bởi bnm w n w m , cuối cùng số hạng phi 2 tuyến bnm w n w m được thay thế bởi số hạng tuyến tính cnm w m . Trong sơ đồ này ta chỉ giới hạn sử dụng phương pháp tuyến tính hóa điều chỉnh một bước để tính toán các đáp ứng của hệ. Đặt 1 2 4 enm w n w m a nm w n w m , 2 4 2 enm a nm w n w m b nm w n w m , (3.30) 3 2 enm b nm w n w m c nm w m. 1 2 3 Các đại lượng enm,, e nm e nm là các sai số giữa các số hạng phi tuyến trong mỗi lần thay thế. Đòi hỏi trung bình phương của các số hạng này phải đạt cực tiểu theo các hệ số tương ứng anm,, b nm c nm : 2 1 2 4 2 E enm E w n w m a nm w n w m min, (3.31) anm 2 2 4 2 2 E enm E a nm w n w m b nm w n w m min, (3.32) bnm
  44. 40 2 3 2 2 E enm E b nm w n w m c nm w m min. (3.33) cnm Các tiêu chuẩn sai số bình phương trung bình (3.31), (3.32), (3.33) dẫn tới lần lượt các phương trình sau đây để xác định các hệ số anm,, b nm c nm  2 E w2 w a w 4 w 0, (3.34) n m nm n m anm  2 E a w4 w b w 2 w 0, (3.35) nm n m nm n m bnm  2 E b w2 w c w 0. (3.36) nm n m nm m cnm Từ (3.34), (3.35), (3.36) ta thu được các hệ số như sau E w6 w 2 a n m , (3.37) nm 8 2 E wn w m E w6 w 2 b a n m , (3.38) nm nm 4 2 E wn w m E w2 w 2 c b n m . (3.39) nm nm 2 E wm Thế (3.37) và (3.38) vào (3.39) ta được biểu thức hệ số của số hạng tuyến tính hóa E w2 w 2 E w 6 w 2 E w 6 w 2 c n m n m n m . (3.40) nm 2 4 2 8 2 E wm E w n w m E w n w m Chú ý rằng với quá trình ngẫu nhiên Gaussian trung bình không, ta có công thức sau [20]: E z1 z 2 z 2m   E z j z k , (3.41)  j k trong đó tổng lấy theo  bao gồm tất cả các cặp độc lập của z1, z 2 , , z 2m, có tất cả số 2m ! cặp độc lập là . Sử dụng (3.41) ta biểu diễn được các đại lượng sau đây qua các 2m m ! mô men bậc hai
  45. 41 2 2 2 E wn w m  nn  mm 2  nm , 4 2 2 2 E wn w m 3 nn  mm 12  nn  nm , (3.42) 6 2 3 2 2 E wn w m 15 nn  mm 90  nn  nm , 8 2 4 2 3 3 2 E wn w m 105 nn  mm 120  nn  nm 720  nn  nm , trong đó ký hiệu nm là các mô men bậc hai 2 nm  mn E w n w m,.  nn E w n (3.43) Thế các biểu thức mô men trong (3.42) vào biểu thức (3.40): 2 2 2 5nn  nn  mm 2  nm  nn  mm 6  nm c . nm 2 2 3 2 (3.44) mm  nn  mm 4  nm 7  nn  mm 8  nm 48  nn  nm 2 2 Trong biểu thức (3.44) chú ý rằng cnm c mn . Như vậy thành phần phi tuyến  n wn w m n 1 của (3.27) được tuyến tính hóa theo phương pháp tuyến tính hóa điều chỉnh trở thành 2 2 2 2 n wn w m  n c nm w m  n c nm w m. (3.45) n 1 n 1 n 1 Và do đó phương trình tuyến tính hóa của (3.27) có dạng sau đây  1 Q w w  21 n 2 c w m q t , m  m m 4 2 2  nm m (3.46) A m2 R mn 1 A trong đó cnm được xác định bởi biểu thức (3.44). So sánh (3.28) và (3.46), biểu thức của hệ số tuyến tính hóa tương đương km, eq là 1 k 1 n2 c . m, eq4 2 2  nm (3.47) m2 R m n 1 Thay (3.44) vào (3.47) ta được 2 2 2 2 5 n nn  nn  mm 2  nm  nn  mm 6  nm k 1 . (3.48) m, eq 4 2 2  2 2 3 2 m2 R m n 1 mm  nn  mm 4  nm 7  nn  mm 8  nm 48  nn  nm
  46. 42 Như vậy hệ số tuyến tính hóa tương đương km, eq phụ thuộc vào các mô men bậc hai nm . Bây giờ ta tìm các mô men bậc hai nm khi biết hàm mật độ phổ đầu vào S  . Hàm truyền Hm  ứng với phương trình (3.28) là 1 H  . (3.49) m  k2  2 i  m, eq m A Hàm Hm  liên quan đến xung đơn vị hm t ứng với (3.28) bởi cặp biến đổi Fourier sau đây 1 h t H ei t d , (3.50) m m 2 H h t e i t dt. (3.51) m m Nghiệm wm t của (3.28) có thể biểu diễn dưới dạng Q w h t n q  d  , n n 1 1 1 A (3.52) Q w h t m q  d  . m m 2 2 2 A Sử dụng (3.52) ta có QQ EwwEht n qdht    m q  d   n m n 1 1 1 m 2 2 2 AA (3.53) QQ n m h t  h t  E q  q  d  d  . 2 n 1 m 2 1 2 1 2 A Thay biểu thức (3.4) vào (3.53) ta được QQ Eww n m htht   Seddd i t2 t 1     n m 2 n 1 m 2 1 2 A QQ n m h t  e i  t 1 d  h t  e i  t  2 d  S  d  2 n 1 1 m 2 2 A (3.54) QQ n m H  H  S  d  2 n m A QQn m Snm , A 2
  47. 43 trong đó S S H  H  S  d . (3.55) nm mn n m Trong trường hợp kích động ngoài là ồn trắng với hàm mật độ phổ SS  0 const thì (3.55) có dạng S S S H  H  d . (3.56) nm mn0 n m Nếu biết được Snm thì từ (3.54) ta sẽ tính được các mô men bậc hai E wn w m . Thay biểu thức của Hm  từ (3.49) vào (3.56) ta được d SS . (3.57) nm 0 2 2 2 2  kn,, eq n  i  k m eq  m  i  AA Sử dụng định lý thặng dư trong lý thuyết hàm biến phức để tính tích phân suy rộng (3.57) ta thu được kết quả sau  1 SSnm 4 0 . (3.58) A 2 22  2 2 kneqn,,,, k meqm  2 k neqn  k meqm  A QQn m Ta thay nm E w n w m S nm từ (3.54) vào biểu thức của km, eq từ (3.48) ta được A 2 5 QQQSSSS2 2 2 k 1 n2 n m n nn nn mm nm m, eq 4 2 2 2 2 2 m2 R mn 1 A Qm S mm S nn S mm 4 S nm (3.59) 2 2 SSSnn mm 6 nm 2 2 3 2 2 . 7QSSQQSQSSn nn mm 8 n m nm 48 n nn nm Trong biểu thức lấy tổng của (3.59) nếu ta giới hạn tổng của M số hạng đầu tiên, nghĩa là biểu thức của km, eq ( m 1,2, , M ) là một tổng của hữu hạn các số hạng:
  48. 44 5 M QQQSSSS2 2 2 k 1 n2 n m n nn nn mm nm m, eq 4 2 2 2 2 2 m2 R mn 1 A Qm S mm S nn S mm 4 S nm (3.60) 2 2 SSSnn mm 6 nm 2 2 3 2 2 . 7QSSQQSQSSn nn mm 8 n m nm 48 n nn nm Khi đó nếu thay (3.58) vào (3.60) ta thu được một hệ phi tuyến đóng kín đối với M ẩn km, eq ( m 1,2, , M ). Ta phải giải số hệ này để thu được các giá trị km, eq . Sau đó thay vào (3.58) ta thu được Snm , rồi sử dụng (3.54) để tìm các mô men bậc hai E wn w m  của phương trình tuyến tính hóa (3.28). Nếu sử dụng phương pháp tuyến tính hóa kinh điển thì biểu thức của hệ số tuyến tính hóa tương đương cho bởi [20] 2 1 M QSSS 2 2 k 1 n2 n nn mm nm . m, eq .kd 4 2 2  (3.61) m2 R mn 1 A Smm Bài toán dao động ngẫu nhiên của dầm phi tuyến (3.1) còn được nghiên cứu bằng phương pháp năng lượng, vấn đề này được trình bày trong [18]. Theo đó, các hệ số tuyến tính hóa km, eq được tìm từ hệ T 2 T k16 k Mk4 BEwUEwU 1 2 2 EwU 2 , (3.62) 1,eq 2, eq M , eq 2 1 2 M  0 trong đó ma trận B và hàm thế năng U được xác định bởi E w2 w 2 E w 2 w 2 E w 2 w 2 1 1 1 2 1 M E w2 w 2 E w 2 w 2 E w 2 w 2 B 2 1 2 2 2 M , (3.63) E w2 w 2 E w 2 w 2 E w 2 w 2 MMMM1 2 2  2 MM1 U 0 m4 w 2 m 2 w 2 .  m2  m (3.64) 2 m 1 4R m 1 Các hệ phương trình (3.60), (3.61) và (3.62) là các hệ đại số phi tuyến đối với M ẩn km, eq ( m 1,2, , M ), việc tìm nghiệm giải tích của nó là rất khó. Do đó ta phải giải chúng bằng phương pháp số. Trong phần tính toán số của nghiên cứu này tác giả chỉ xét trường hợp riêng giới hạn với giá trị M 3. Tức là biểu diễn (3.5) chỉ với 3 mode
  49. 45 là w1,, w 2 w 3 . Tuy nhiên đại lượng Q2 0 (xem (3.19)) nên từ (3.27) ta có thể thấy rằng đáp ứng w2 bằng không. Các kết quả số được thực hiện trên phần mềm MATLAB 7.2 (Xem phần Phụ lục). Tổng quát hơn thì M phương trình (3.27) sẽ rút xuống còn tất cả các phương trình ứng với giá trị m là số lẻ. Hình 1. ng bình ph ình 2 theo tham s v i . Đáp ứ ương trung b E w1 ố R ớ S0 1, 0,  0.1 Hình 2 ng bình ph ình 2 theo tham s v i . . Đáp ứ ương trung b E w1 ố R ớ S0 1, 1,  0.1
  50. 46 Hình 3 ng bình ph ình 2 theo tham s v i . . Đáp ứ ương trung b E w1 ố R ớ S0 5, 1,  0.1 Hình 4 ng bình ph ình 2 theo tham s v i . . Đáp ứ ương trung b E w1 ố ớ SR0 5, 0.1, 1 Các Hình 1, 2, 3, 4 trình bày đồ thị của đáp ứng bình phương trung bình của mode đầu tiên w1 thay đổi theo tham số R và với các phương pháp khác nhau trong trường hợp số mode bằng 3 ( A 1, 0 1). Đồ thị đường nét liền là phương pháp tuyến tính hóa điều chỉnh (3.60), đồ thị đường nét đứt là phương pháp năng lượng (3.62), đồ thị với nét chấm chấm là phương pháp tuyến tính hóa kinh điển (3.61), các dấu sao “*” chỉ các kết quả mô phỏng số bằng phương pháp Monte-Carlo. Các đồ thị cho thấy rằng
  51. 47 phương pháp tuyến tính hóa điều chỉnh cho kết quả gần với kết quả mô phỏng hơn hai phương pháp còn lại, do đó đây là phương pháp cho kết quả tốt đối với dàm Euler- Bernoulli tựa giản đơn chịu kích động ngẫu nhiên. Tổng kết chương 3 Chương 3 trình bày áp dụng phương pháp tuyến tính hóa điều chỉnh cho bài toán dao động của dầm phi tuyến chịu kích động ngẫu nhiên. Kết quả thu được là một hệ phương trình đại số phi tuyến đóng kín (3.60) đối với các ẩn là các hệ số tuyến tính hóa tương đương. Phương pháp số được sử dụng để giải hệ phi tuyến này. Kết quả tính toán số được thực hiện cho ba đáp ứng đầu tiên của hệ phi tuyến (3.27). Các đáp ứng bình phương trung bình thu được bằng phương pháp tuyến tính hóa điều chỉnh được so sánh với các phương pháp đã biết gồm phương pháp năng lượng, phương pháp tuyến tính hóa kinh điển và phương pháp mô phỏng Monte-Carlo. Các kết quả chỉ ra rằng phương pháp tuyến tính hóa điều chỉnh có thể cho kết quả tốt hơn so với phương pháp tuyến tính hóa kinh điển.
  52. 48 KẾT LUẬN Phương pháp tuyến tính hóa tương đương là một trong những phương pháp quan trọng để nghiên cứu các hệ động lực chịu kích động ngẫu nhiên. Phương pháp được phát triển và mở rộng bởi nhiều nhà nghiên cứu trong nhiều thập kỷ qua. Với phương phương pháp tuyến tính hóa điều chỉnh, tuy được đưa ra từ năm 1995 nhưng có rất ít các công trình nghiên cứu tiếp theo. Các công trình nghiên cứu tập chung vào các hệ rời rạc, còn các hệ liên tục vẫn còn là một khoảng trống. Trong luận văn này tác giả đã sử dụng phương pháp tuyến tính hóa điều chỉnh để giải quyết các bài toán cơ hệ phi tuyến rời rạc và liên tục với nội dung bao gồm: Hệ thống lại phương pháp tuyến tính hóa đương cho hệ một và nhiều bậc tự do, đồng trình bày những kết quả nổi bật của phương pháp tuyến tính hóa điều chỉnh cho một số hệ phi tuyến điển hình. Mở rộng phương pháp tuyến tính hóa điều chỉnh cho hệ liên tục là một dầm phi tuyến chịu kích động ngẫu nhiên và thu được kết quả là một hệ kín các phương trình đại số phi tuyến đối với ẩn là các hệ số của các phương trình tuyến tính hóa tương đương. Tính toán số đối với hệ đại số phi tuyến để tìm các đáp ứng bình phương trung bình đầu tiên của cơ hệ phi tuyến và so sánh chúng với các phương pháp năng lượng, phương pháp tuyến tính hóa kinh điển và phương pháp mô phỏng Monte-Carlo. Các tính toán số được thực hiện trên phần mềm Matlat. Luận văn mới chỉ dừng lại cho việc áp dụng phương pháp tuyến tính hóa điều chỉnh cho dầm phi tuyến với hai đầu tựa giản đơn. Phương pháp này còn cần được nghiên cứu để có thể áp dụng cho nhiều hệ phi tuyến khác. Đây là vấn đề cần được tiếp tục phát triển.
  53. 49 TÀI LIỆU THAM KHẢO 1. R.C. Booton (1954), “The analysis of nonlinear central systems with random inputs”, IRE Trans. Circuit Theory, 1, pp. 32-34. 2. I.E. Kazakov (1954), “An approximate method for the statistical investigation for nonlinear systems”, Trudy VVIA im Prof. N. E. Zhukovskogo, 394, pp. 1-52 (in Russian). 3. T.K. Caughey (1963), “Equivalent linearization techniques”, Journal of the Acoustical Society of America, 35, pp. 1706-1711 (Reference is made to presentations of the procedure in lectures delivered in 1953 at the California Institute of Technology). 4. T.K. Caughey (1959), “Response of Van der Pol's oscillator to random excitation”, J App Mech, 26, pp. 345-348. 5. N.M. Krylov and N.N. Bogoliubov (1937), Introduction to nonlinear mechanics, Academy of Sciences of Ukraine Publishing House (in Russian). 6. R.L. Stratonovich (1967), Topics in the Theory of Random Noise, Vol. II. Gordon and Breach. New York. 7. T.S. Atalik and S. Utku (1976), “Stochastic of linearization of multi-degree of freedom nonlinear”, J Earth Eng Struct Dynamic, 4, pp. 441-420. 8. F. Casciati and L. Faravelli (1986), “Equivalent linearization in nonlinear random vibration problems”, Proc Conf on Vibration Problems in Eng, Xian, China, pp. 986-991. 9. J.B. Roberts and P.D. Spanos (1990), Random vibration and statistical linearization, Wiley, New York. 10. L. Socha and T.T. Soong (1991), “Linearization in analysis of nonlinear stochastic systems”, Appl Mech Rev, 44, pp. 399-422. 11. L. Socha (2008), Linearization methods for stochastic dynamic systems, Springer, Lecture Notes in Physics. 12. N.D. Anh and M. Di Paola (1995), “Some extensions of Gaussian equivalent linearization”, International Conference on Nonlinear Stochastic Dynamics, Hanoi, Vietnam, pp 5-16. 13. N.D. Anh and W. Schiehlen (1997), “New criterion for Gaussian equivalent linearization”, European Journal of Mechanics, A Solids, 16(6), pp. 1025-1039. 14. N.D. Anh and W. Schiehlen (1997), “An extension to the mean square criterion of Gaussian equivalent linearization”, Vietnam J Math, 25(2), pp. 115-123. 15. I. Elishakoff, L. Andrimasy and M. Dolley (2008), “Application and extension of the stochastic linearization by Anh and Di Paola”, Acta Mech, 204, pp. 89-98. 16. L.D. Lutes and S. Sarkani, (2004), "Random vibration: Analysis of structural and mechanical systems", Elsevier, Amsterdam, pp. 423-424, 437-438. 17. I. Elishakoff (2000), "Stochastic linearization technique: A new interpretation and a selective review", Shock Vibration Dig, 32, pp. 179-188. 18. J.J. Fang, I. Elishakoff and R. Caimi (1995), “Nonlinear response of a beam under stationary random excitation by improved stochastic linearization method”, Applied Mathematical Modeling, 19, pp. 106-111. 19. P.D. Spanos (1980), “Formulation of stochastic linearization for symmetric or asymmetric MDOF nonlinear systems”, Journal of Applied Mechanics, ASME, 47, pp. 209-211. 20. P. Seide (1976), “Nonlinear stresses and deflections of beams subjected to random time dependent uniform pressure”, J. Eng. Indus, No. 75-DET-23. 21. N.D. Anh, I. Elishakoff and N.N. Hieu (2011), “Extension of the regulated stochastic linearization to beam vibrations” (to be submitted to Vietnam Journal of Mechanics). 22. N.D. Anh and N.N. Hieu (2010), “The Duffing oscillator under combined periodic and random excitations” (to be submitted to Journal of Probabilistic Engineering Mechanics, Elsevier). 23. N.D. Anh, N.N. Hieu and N.N. Linh (2011), “A dual criterion of equivalent linearization method for nonlinear systems subjected to random excitation” (to be submitted to Acta Mechanica, Springer).
  54. 50 PHỤ LỤC Các chương trình Matlab cho các hình 1, 2, 3, 4 (chương 3) và Mô phỏng số với công cụ Simulink Giải thích: S0: Mật độ phổ hằng số của kích động đầu vào  bet0: Đại lượng này được định nghĩa bởi bet0 trong đó  0.1, A 1 A ome: Đại lượng này là 0 1 rhoA: Đại lượng A 1 K f alp: Đại lượng 2 A 0 1. Chương trình Matlab cho Hình 1 % %Copyright @ 2011 Nguyen Nhu Hieu %Date: May 10, 2011 %Faculty of Engineering Mechanics and Automation, University of Engineering and %Technology, Vietnam National University % function Mean_Square_Response_W1 clc S0=1; bet0=0.1; ome=1; rhoA=1; alp=0; lam0=8*S0/(pi*rhoA*bet0*ome^2); R=0.05:0.05:2; num=max(size(R)); for i=1:1:num lam0=8*S0/(pi*rhoA*bet0*ome^2); ka=fsolve(@(ka) MyFunRL(ka, R(i), S0, bet0, ome, rhoA, alp), [0.1 9]); k=fsolve(@(k) MyFunEM(k, R(i), S0, bet0, ome, rhoA, alp), [0.1 0.1]); kc=fsolve(@(kc) MyFunCL(kc, R(i), S0, bet0, ome, rhoA, alp), [0.1 0.1 9]); lambdaRL(i)=lam0/ka(1); lambdaEM(i)=lam0/k(1); lambdaCL(i)=lam0/kc(1); end plot(R, lambdaRL,'r-'); hold on plot(R,lambdaEM,'b ');
  55. 51 hold on plot(R, lambdaCL,'k:'); hold on R_1=[0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0]; E_Smu1=[ 0.9619 1.9012 2.8651 3.7852 4.6582 5.4410 6.1916 6.9486 7.5441 8.1681] ; %Mean square responses by Monte-Carlo simulation for i=1:1:10 plot(R_1(i),E_Smu1(i),'r*'); hold on end set(gcf, 'color', 'white'); xlabel('R'); ylabel('mean square response E[W_{1}^{2}]'); legend('Present Method','Energy Method','Conventional Method','Simulation'); title('S_{0}=1, \alpha =0, \beta =0.1'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function FunRL=MyFunRL(ka, R, S0, bet0, ome, rhoA, alp); % Regulated Stochastic Linearization %ka: Linearization coefficient %R: Radius of gyration of the beam cross section %S0: Spectral density %bet0=beta/rhoA S11=pi*S0/(bet0*ome^2*ka(1)); S33=pi*S0/(81*bet0*ome^2*ka(2)); S13=4*pi*S0*bet0/((ka(1)-81*ka(2))^2*ome^4+2*bet0^2*(ka(1)+81*ka(2))*ome^2); Q1=2*sqrt(2)/pi; Q3=2*sqrt(2)/(3*pi); lam11=(Q1/rhoA)^2*S11; lam33=(Q3/rhoA)^2*S33; lam13=(Q1*Q3)/(rhoA)^2*S13; ka1eq=1+alp+1/(2*R^2)*(7/3*lam11+45*(lam11*lam33+2*lam13^2)*(lam11*lam33 ^3+6*lam13^2*lam33^2)^2/(lam11*(lam11*lam33^2+4*lam33*lam13^2)*(7*lam11* lam33^4+8*lam13^3*lam33^2+48*lam13^2*lam33^3))); ka3eq=1+alp/81+1/(18*R^2)*(21*lam33+5*(lam11*lam33+2*lam13^2)*(lam11^3*la m33+6*lam13^2*lam11^2)^2/(lam33*(lam11^2*lam33+4*lam11*lam13^2)*(7*lam3 3*lam11^4+8*lam13^3*lam11^2+48*lam13^2*lam11^3))); FunRL=[ka1eq-ka(1); ka3eq-ka(2)]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function FunEM=MyFunEM(k, R, S0, bet0, ome, rhoA, alp); %Energy Method %k: Linearization coefficient
  56. 52 S11=pi*S0/(bet0*ome^2*k(1)); S33=pi*S0/(81*bet0*ome^2*k(2)); S13=4*pi*S0*bet0/((k(1)-81*k(2))^2*ome^4+2*bet0^2*(k(1)+81*k(2))*ome^2); Q1=2*sqrt(2)/pi; Q3=2*sqrt(2)/(3*pi); lam11=(Q1/rhoA)^2*S11; lam33=(Q3/rhoA)^2*S33; lam13=(Q1*Q3)/(rhoA)^2*S13; E2W1=lam11; EW1W3=lam13; E2W3=lam33; E4W1=3*lam11^2; E4W3=3*lam33^2; E6W1=15*lam11^3; E6W3=15*lam33^3; E22W1W3=lam11*lam33+2*lam13^2; % E[w1^2*W^2] E42W1W3=3*lam11^2*lam33+12*lam11*lam13^2; E24W1W3=3*lam33^2*lam11+12*lam33*lam13^2; E2W1U=ome^2/2*((1+alp)*E4W1+(81+alp)*E22W1W3+1/(4*R^2)*(E6W1+81*E24 W1W3+18*E42W1W3)); E2W3U=ome^2/2*((1+alp)*E22W1W3+(81+alp)*E4W3+1/(4*R^2)*(E42W1W3+81 *E6W3+18*E24W1W3)); Det=E4W1*E4W3-(E22W1W3)^2; Det1=2/(ome^2)*(E2W1U*E4W3-E2W3U*E22W1W3); Det3=2/(ome^2)*(E4W1*E2W3U-E22W1W3*E2W1U); k1eq=Det1/Det; k3eq=Det3/(81*Det); FunEM=[k1eq-k(1); k3eq-k(2)]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function FunCL=MyFunCL(kc, R, S0, bet0, ome, rhoA, alp); %Conventional Linearization Technique %kc: Linearization coefficient Q1=2*sqrt(2)/pi; Q2=0; Q3=2*sqrt(2)/(3*pi); S11=pi*S0/(bet0*ome^2*kc(1)); S22=pi*S0/(16*bet0*ome^2*kc(2)); S33=pi*S0/(81*bet0*ome^2*kc(3)); S12=4*pi*S0*bet0/((kc(1)-16*kc(2))^2*ome^4+2*bet0^2*(kc(1)+16*kc(2))*ome^2); S13=4*pi*S0*bet0/((kc(1)-81*kc(3))^2*ome^4+2*bet0^2*(kc(1)+81*kc(3))*ome^2); S23=4*pi*S0*bet0/((16*kc(2)- 81*kc(3))^2*ome^4+2*bet0^2*(16*kc(2)+81*kc(3))*ome^2);
  57. 53 k1eq=1+alp+1/(2*R^2)*((Q1/rhoA)^2*(S11*S11+2*S11^2)/S11+4*(Q2/rhoA)^2*(S1 1*S22+2*S12^2)/S11+9*(Q3/rhoA)^2*(S11*S33+2*S13^2)/S11); k2eq=1+alp/16+1/(8*R^2)*((Q1/rhoA)^2*(S22*S11+2*S12^2)/S22+4*(Q2/rhoA)^2* (S22*S22+2*S22^2)/S22+9*(Q3/rhoA)^2*(S22*S33+2*S23^2)/S22); k3eq=1+alp/81+1/(18*R^2)*((Q1/rhoA)^2*(S33*S11+2*S13^2)/S33+4*(Q2/rhoA)^2 *(S33*S22+2*S23^2)/S33+9*(Q3/rhoA)^2*(S33*S33+2*S33^2)/S33); lam11=(Q1/rhoA)^2*S11; FunCL=[k1eq-kc(1); k2eq-kc(2); k3eq-kc(3)]; % END PROGRAM 2. Chương trình Matlab cho Hình 2 % %Copyright @ 2011 Nguyen Nhu Hieu %Date: May 10, 2011 %Faculty of Engineering Mechanics and Automation, University of Engineering and %Technology, Vietnam National University % function Mean_Square_Response_W1 clc S0=1; bet0=0.1; ome=1; rhoA=1; alp=1; lam0=8*S0/(pi*rhoA*bet0*ome); R=0.05:0.05:2; num=max(size(R)); for i=1:1:num lam0=8*S0/(pi*rhoA*bet0*ome^2); ka=fsolve(@(ka) MyFunRL(ka, R(i), S0, bet0, ome, rhoA, alp), [0.1 9]); k=fsolve(@(k) MyFunEM(k, R(i), S0, bet0, ome, rhoA, alp), [0.1 0.1]); kc=fsolve(@(kc) MyFunCL(kc, R(i), S0, bet0, ome, rhoA, alp), [0.1 0.1 9]); lambdaRL(i)=lam0/ka(1); lambdaEM(i)=lam0/k(1); lambdaCL(i)=lam0/kc(1); end plot(R, lambdaRL,'r-'); hold on plot(R,lambdaEM,'b '); hold on plot(R, lambdaCL,'k:'); hold on R_1=[0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0];
  58. 54 E_Smu1=[0.9391 1.8162 2.6608 3.4275 4.1112 4.7288 5.3034 5.7712 6.1482 6.5725]; %Mean square responses by Monte Carlo simulation for i=1:1:10 plot(R_1(i),E_Smu1(i),'r*'); hold on end set(gcf, 'color', 'white'); xlabel('R'); ylabel('mean square response E[W_{1}^{2}]'); legend('Present Method','Energy Method','Conventional Method','Simulation'); title('S_{0}=1, \alpha =1, \beta =0.1'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function FunRL=MyFunRL(ka, R, S0, bet0, ome, rhoA, alp); % Regulated Stochastic Linearization %ka: Linearization coefficient %R: Radius of gyration of the beam cross section %S0: Spectral density %bet0=beta/rhoA S11=pi*S0/(bet0*ome^2*ka(1)); S33=pi*S0/(81*bet0*ome^2*ka(2)); S13=4*pi*S0*bet0/((ka(1)-81*ka(2))^2*ome^4+2*bet0^2*(ka(1)+81*ka(2))*ome^2); Q1=2*sqrt(2)/pi; Q3=2*sqrt(2)/(3*pi); lam11=(Q1/rhoA)^2*S11; lam33=(Q3/rhoA)^2*S33; lam13=(Q1*Q3)/(rhoA)^2*S13; ka1eq=1+alp+1/(2*R^2)*(7/3*lam11+45*(lam11*lam33+2*lam13^2)*(lam11*lam33 ^3+6*lam13^2*lam33^2)^2/(lam11*(lam11*lam33^2+4*lam33*lam13^2)*(7*lam11* lam33^4+8*lam13^3*lam33^2+48*lam13^2*lam33^3))); ka3eq=1+alp/81+1/(18*R^2)*(21*lam33+5*(lam11*lam33+2*lam13^2)*(lam11^3*la m33+6*lam13^2*lam11^2)^2/(lam33*(lam11^2*lam33+4*lam11*lam13^2)*(7*lam3 3*lam11^4+8*lam13^3*lam11^2+48*lam13^2*lam11^3))); FunRL=[ka1eq-ka(1); ka3eq-ka(2)]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function FunEM=MyFunEM(k, R, S0, bet0, ome, rhoA, alp); %Energy Method %k: Linearization coefficient S11=pi*S0/(bet0*ome^2*k(1)); S33=pi*S0/(81*bet0*ome^2*k(2)); S13=4*pi*S0*bet0/((k(1)-81*k(2))^2*ome^4+2*bet0^2*(k(1)+81*k(2))*ome^2); Q1=2*sqrt(2)/pi; Q3=2*sqrt(2)/(3*pi); lam11=(Q1/rhoA)^2*S11; lam33=(Q3/rhoA)^2*S33;
  59. 55 lam13=(Q1*Q3)/(rhoA)^2*S13; E2W1=lam11; EW1W3=lam13; E2W3=lam33; E4W1=3*lam11^2; E4W3=3*lam33^2; E6W1=15*lam11^3; E6W3=15*lam33^3; E22W1W3=lam11*lam33+2*lam13^2; % E[w1^2*W^2] E42W1W3=3*lam11^2*lam33+12*lam11*lam13^2; E24W1W3=3*lam33^2*lam11+12*lam33*lam13^2; E2W1U=ome^2/2*((1+alp)*E4W1+(81+alp)*E22W1W3+1/(4*R^2)*(E6W1+81*E24 W1W3+18*E42W1W3)); E2W3U=ome^2/2*((1+alp)*E22W1W3+(81+alp)*E4W3+1/(4*R^2)*(E42W1W3+81 *E6W3+18*E24W1W3)); Det=E4W1*E4W3-(E22W1W3)^2; Det1=2/(ome^2)*(E2W1U*E4W3-E2W3U*E22W1W3); Det3=2/(ome^2)*(E4W1*E2W3U-E22W1W3*E2W1U); k1eq=Det1/Det; k3eq=Det3/(81*Det); FunEM=[k1eq-k(1); k3eq-k(2)]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function FunCL=MyFunCL(kc, R, S0, bet0, ome, rhoA, alp); %Conventional Linearization Technique %kc: Linearization coefficient Q1=2*sqrt(2)/pi; Q2=0; Q3=2*sqrt(2)/(3*pi); S11=pi*S0/(bet0*ome^2*kc(1)); S22=pi*S0/(16*bet0*ome^2*kc(2)); S33=pi*S0/(81*bet0*ome^2*kc(3)); S12=4*pi*S0*bet0/((kc(1)-16*kc(2))^2*ome^4+2*bet0^2*(kc(1)+16*kc(2))*ome^2); S13=4*pi*S0*bet0/((kc(1)-81*kc(3))^2*ome^4+2*bet0^2*(kc(1)+81*kc(3))*ome^2); S23=4*pi*S0*bet0/((16*kc(2)- 81*kc(3))^2*ome^4+2*bet0^2*(16*kc(2)+81*kc(3))*ome^2); k1eq=1+alp+1/(2*R^2)*((Q1/rhoA)^2*(S11*S11+2*S11^2)/S11+4*(Q2/rhoA)^2*(S1 1*S22+2*S12^2)/S11+9*(Q3/rhoA)^2*(S11*S33+2*S13^2)/S11); k2eq=1+alp/16+1/(8*R^2)*((Q1/rhoA)^2*(S22*S11+2*S12^2)/S22+4*(Q2/rhoA)^2* (S22*S22+2*S22^2)/S22+9*(Q3/rhoA)^2*(S22*S33+2*S23^2)/S22); k3eq=1+alp/81+1/(18*R^2)*((Q1/rhoA)^2*(S33*S11+2*S13^2)/S33+4*(Q2/rhoA)^2 *(S33*S22+2*S23^2)/S33+9*(Q3/rhoA)^2*(S33*S33+2*S33^2)/S33); lam11=(Q1/rhoA)^2*S11; FunCL=[k1eq-kc(1); k2eq-kc(2); k3eq-kc(3)]; % END PROGRAM
  60. 56 3. Chương trình Matlab cho Hình 3 % %Copyright @ 2011 Nguyen Nhu Hieu %Date: May 10, 2011 %Faculty of Engineering Mechanics and Automation, University of Engineering and %Technology, Vietnam National University % function Mean_Square_Response_W1 clc S0=5; bet0=0.1; ome=1; rhoA=1; alp=1; lam0=8*S0/(pi*rhoA*bet0*ome); R=0.05:0.05:2; num=max(size(R)); for i=1:1:num lam0=8*S0/(pi*rhoA*bet0*ome^2); ka=fsolve(@(ka) MyFunRL(ka, R(i), S0, bet0, ome, rhoA, alp), [0.1 9]); k=fsolve(@(k) MyFunEM(k, R(i), S0, bet0, ome, rhoA, alp), [0.1 0.1]); kc=fsolve(@(kc) MyFunCL(kc, R(i), S0, bet0, ome, rhoA, alp), [0.1 0.1 9]); lambdaRL(i)=lam0/ka(1); lambdaEM(i)=lam0/k(1); lambdaCL(i)=lam0/kc(1); end plot(R, lambdaRL,'r-'); hold on plot(R,lambdaEM,'b '); hold on plot(R, lambdaCL,'k:'); hold on %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% hold on R_1=[0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0]; E_Smu1=[ 2.1812 4.2702 6.1984 8.1410 10.0478 12.0609 13.8322 15.5734 17.2122 18.8455 ]; for i=1:1:10 plot(R_1(i),E_Smu1(i),'r*'); hold on end set(gcf, 'color', 'white'); xlabel('R'); ylabel('mean square response E[W_{1}^{2}]'); legend('Present Method','Energy Method','Conventional Method','Simulation');
  61. 57 title('S_{0}= 5, \alpha = 1, \beta = 0.1'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function FunRL=MyFunRL(ka, R, S0, bet0, ome, rhoA, alp);% Regulated Stochastic Linearization %ka: Linearization coefficient %R: Radius of gyration of the beam cross section %S0: Spectral density %bet0=beta/rhoA S11=pi*S0/(bet0*ome^2*ka(1)); S33=pi*S0/(81*bet0*ome^2*ka(2)); S13=4*pi*S0*bet0/((ka(1)-81*ka(2))^2*ome^4+2*bet0^2*(ka(1)+81*ka(2))*ome^2); Q1=2*sqrt(2)/pi; Q3=2*sqrt(2)/(3*pi); lam11=(Q1/rhoA)^2*S11; lam33=(Q3/rhoA)^2*S33; lam13=(Q1*Q3)/(rhoA)^2*S13; ka1eq=1+alp+1/(2*R^2)*(7/3*lam11+45*(lam11*lam33+2*lam13^2)*(lam11*lam33 ^3+6*lam13^2*lam33^2)^2/(lam11*(lam11*lam33^2+4*lam33*lam13^2)*(7*lam11* lam33^4+8*lam13^3*lam33^2+48*lam13^2*lam33^3))); ka3eq=1+alp/81+1/(18*R^2)*(21*lam33+5*(lam11*lam33+2*lam13^2)*(lam11^3*la m33+6*lam13^2*lam11^2)^2/(lam33*(lam11^2*lam33+4*lam11*lam13^2)*(7*lam3 3*lam11^4+8*lam13^3*lam11^2+48*lam13^2*lam11^3))); FunRL=[ka1eq-ka(1); ka3eq-ka(2)]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function FunEM=MyFunEM(k, R, S0, bet0, ome, rhoA, alp); %Energy Method %k: Linearization coefficient S11=pi*S0/(bet0*ome^2*k(1)); S33=pi*S0/(81*bet0*ome^2*k(2)); S13=4*pi*S0*bet0/((k(1)-81*k(2))^2*ome^4+2*bet0^2*(k(1)+81*k(2))*ome^2); Q1=2*sqrt(2)/pi; Q3=2*sqrt(2)/(3*pi); lam11=(Q1/rhoA)^2*S11; lam33=(Q3/rhoA)^2*S33; lam13=(Q1*Q3)/(rhoA)^2*S13; E2W1=lam11; EW1W3=lam13; E2W3=lam33; E4W1=3*lam11^2; E4W3=3*lam33^2; E6W1=15*lam11^3; E6W3=15*lam33^3; E22W1W3=lam11*lam33+2*lam13^2; % E[w1^2*W^2]
  62. 58 E42W1W3=3*lam11^2*lam33+12*lam11*lam13^2; E24W1W3=3*lam33^2*lam11+12*lam33*lam13^2; E2W1U=ome^2/2*((1+alp)*E4W1+(81+alp)*E22W1W3+1/(4*R^2)*(E6W1+81*E24 W1W3+18*E42W1W3)); E2W3U=ome^2/2*((1+alp)*E22W1W3+(81+alp)*E4W3+1/(4*R^2)*(E42W1W3+81 *E6W3+18*E24W1W3)); Det=E4W1*E4W3-(E22W1W3)^2; Det1=2/(ome^2)*(E2W1U*E4W3-E2W3U*E22W1W3); Det3=2/(ome^2)*(E4W1*E2W3U-E22W1W3*E2W1U); k1eq=Det1/Det; k3eq=Det3/(81*Det); FunEM=[k1eq-k(1); k3eq-k(2)]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function FunCL=MyFunCL(kc, R, S0, bet0, ome, rhoA, alp); %Conventional Linearization Technique %kc: Linearization coefficient Q1=2*sqrt(2)/pi; Q2=0; Q3=2*sqrt(2)/(3*pi); S11=pi*S0/(bet0*ome^2*kc(1)); S22=pi*S0/(16*bet0*ome^2*kc(2)); S33=pi*S0/(81*bet0*ome^2*kc(3)); S12=4*pi*S0*bet0/((kc(1)-16*kc(2))^2*ome^4+2*bet0^2*(kc(1)+16*kc(2))*ome^2); S13=4*pi*S0*bet0/((kc(1)-81*kc(3))^2*ome^4+2*bet0^2*(kc(1)+81*kc(3))*ome^2); S23=4*pi*S0*bet0/((16*kc(2)- 81*kc(3))^2*ome^4+2*bet0^2*(16*kc(2)+81*kc(3))*ome^2); k1eq=1+alp+1/(2*R^2)*((Q1/rhoA)^2*(S11*S11+2*S11^2)/S11+4*(Q2/rhoA)^2*(S1 1*S22+2*S12^2)/S11+9*(Q3/rhoA)^2*(S11*S33+2*S13^2)/S11); k2eq=1+alp/16+1/(8*R^2)*((Q1/rhoA)^2*(S22*S11+2*S12^2)/S22+4*(Q2/rhoA)^2* (S22*S22+2*S22^2)/S22+9*(Q3/rhoA)^2*(S22*S33+2*S23^2)/S22); k3eq=1+alp/81+1/(18*R^2)*((Q1/rhoA)^2*(S33*S11+2*S13^2)/S33+4*(Q2/rhoA)^2 *(S33*S22+2*S23^2)/S33+9*(Q3/rhoA)^2*(S33*S33+2*S33^2)/S33); lam11=(Q1/rhoA)^2*S11; FunCL=[k1eq-kc(1); k2eq-kc(2); k3eq-kc(3)]; % END PROGRAM
  63. 59 4. Chương trình Matlab cho Hình 4 % %Copyright @ 2011 Nguyen Nhu Hieu %Date: May 10, 2011 %Faculty of Engineering Mechanics and Automation, University of Engineering and %Technology, Vietnam National University % function Mean_Square_Response_W1 clc S0=5; bet0=0.1; ome=1; rhoA=1; R=1; lam0=8*S0/(pi*rhoA*bet0*ome^2); alp=0.05:0.05:10; num=max(size(alp)); for i=1:1:num ka=fsolve(@(ka) MyFunRL(ka, alp(i), S0, bet0, ome, rhoA, R), [0.1 9]); k=fsolve(@(k) MyFunEM(k, alp(i), S0, bet0, ome, rhoA, R), [0.1 0.1]); kc=fsolve(@(kc) MyFunCL(kc, alp(i), S0, bet0, ome, rhoA, R), [0.1 0.1 9]); lambdaRL(i)=lam0/ka(1); lambdaEM(i)=lam0/k(1); lambdaCL(i)=lam0/kc(1); end plot(alp, lambdaRL,'r-'); hold on plot(alp,lambdaEM,'b '); hold on plot(alp, lambdaCL,'k:'); hold on alpha=[1 2 3 4 5 6 7 8 9 10]; MMM=[10.1095 9.5256 9.0638 8.6552 8.2053 7.7941 7.4623 7.0973 6.7656 6.4760]; for i=1:1:10 plot(alpha(i), MMM(i), 'r*'); hold on end set(gcf, 'color', 'white'); xlabel('\alpha'); ylabel('mean square response E[W_{1}^{2}]');
  64. 60 legend('Present Method','Energy Method','Conventional Method','Simulation'); title('S_{0}=5, R =1, \beta =0.1'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function FunRL=MyFunRL(ka, alp, S0, bet0, ome, rhoA, R); % Regulated Stochastic Linearization %ka: Linearization coefficient %R: Radius of gyration of the beam cross section %S0: Spectral density %bet0=beta/rhoA S11=pi*S0/(bet0*ome^2*ka(1)); S33=pi*S0/(81*bet0*ome^2*ka(2)); S13=4*pi*S0*bet0/((ka(1)-81*ka(2))^2*ome^4+2*bet0^2*(ka(1)+81*ka(2))*ome^2); Q1=2*sqrt(2)/pi; Q3=2*sqrt(2)/(3*pi); lam11=(Q1/rhoA)^2*S11; lam33=(Q3/rhoA)^2*S33; lam13=(Q1*Q3)/(rhoA)^2*S13; ka1eq=1+alp+1/(2*R^2)*(7/3*lam11+45*(lam11*lam33+2*lam13^2)*(lam11*lam33 ^3+6*lam13^2*lam33^2)^2/(lam11*(lam11*lam33^2+4*lam33*lam13^2)*(7*lam11* lam33^4+8*lam13^3*lam33^2+48*lam13^2*lam33^3))); ka3eq=1+alp/81+1/(18*R^2)*(21*lam33+5*(lam11*lam33+2*lam13^2)*(lam11^3*la m33+6*lam13^2*lam11^2)^2/(lam33*(lam11^2*lam33+4*lam11*lam13^2)*(7*lam3 3*lam11^4+8*lam13^3*lam11^2+48*lam13^2*lam11^3))); FunRL=[ka1eq-ka(1); ka3eq-ka(2)]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function FunEM=MyFunEM(k, alp, S0, bet0, ome, rhoA, R); %Energy Method %k: Linearization coefficient S11=pi*S0/(bet0*ome^2*k(1)); S33=pi*S0/(81*bet0*ome^2*k(2)); S13=4*pi*S0*bet0/((k(1)-81*k(2))^2*ome^4+2*bet0^2*(k(1)+81*k(2))*ome^2); Q1=2*sqrt(2)/pi; Q3=2*sqrt(2)/(3*pi); lam11=(Q1/rhoA)^2*S11; lam33=(Q3/rhoA)^2*S33; lam13=(Q1*Q3)/(rhoA)^2*S13; E2W1=lam11; EW1W3=lam13; E2W3=lam33; E4W1=3*lam11^2; E4W3=3*lam33^2; E6W1=15*lam11^3;
  65. 61 E6W3=15*lam33^3; E22W1W3=lam11*lam33+2*lam13^2; % E[w1^2*W^2] E42W1W3=3*lam11^2*lam33+12*lam11*lam13^2; E24W1W3=3*lam33^2*lam11+12*lam33*lam13^2; E2W1U=ome^2/2*((1+alp)*E4W1+(81+alp)*E22W1W3+1/(4*R^2)*(E6W1+81*E24 W1W3+18*E42W1W3)); E2W3U=ome^2/2*((1+alp)*E22W1W3+(81+alp)*E4W3+1/(4*R^2)*(E42W1W3+81 *E6W3+18*E24W1W3)); Det=E4W1*E4W3-(E22W1W3)^2; Det1=2/(ome^2)*(E2W1U*E4W3-E2W3U*E22W1W3); Det3=2/(ome^2)*(E4W1*E2W3U-E22W1W3*E2W1U); k1eq=Det1/Det; k3eq=Det3/(81*Det); FunEM=[k1eq-k(1); k3eq-k(2)]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function FunCL=MyFunCL(kc, alp, S0, bet0, ome, rhoA, R); %Conventional Linearization Technique %kc: Linearization coefficient Q1=2*sqrt(2)/pi; Q2=0; Q3=2*sqrt(2)/(3*pi); S11=pi*S0/(bet0*ome^2*kc(1)); S22=pi*S0/(16*bet0*ome^2*kc(2)); S33=pi*S0/(81*bet0*ome^2*kc(3)); S12=4*pi*S0*bet0/((kc(1)-16*kc(2))^2*ome^4+2*bet0^2*(kc(1)+16*kc(2))*ome^2); S13=4*pi*S0*bet0/((kc(1)-81*kc(3))^2*ome^4+2*bet0^2*(kc(1)+81*kc(3))*ome^2); S23=4*pi*S0*bet0/((16*kc(2)- 81*kc(3))^2*ome^4+2*bet0^2*(16*kc(2)+81*kc(3))*ome^2); k1eq=1+alp+1/(2*R^2)*((Q1/rhoA)^2*(S11*S11+2*S11^2)/S11+4*(Q2/rhoA)^2*(S1 1*S22+2*S12^2)/S11+9*(Q3/rhoA)^2*(S11*S33+2*S13^2)/S11); k2eq=1+alp/16+1/(8*R^2)*((Q1/rhoA)^2*(S22*S11+2*S12^2)/S22+4*(Q2/rhoA)^2* (S22*S22+2*S22^2)/S22+9*(Q3/rhoA)^2*(S22*S33+2*S23^2)/S22); k3eq=1+alp/81+1/(18*R^2)*((Q1/rhoA)^2*(S33*S11+2*S13^2)/S33+4*(Q2/rhoA)^2 *(S33*S22+2*S23^2)/S33+9*(Q3/rhoA)^2*(S33*S33+2*S33^2)/S33); lam11=(Q1/rhoA)^2*S11; FunCL=[k1eq-kc(1); k2eq-kc(2); k3eq-kc(3)]; % END PROGRAM 5. Chương trình đầu vào cho Mô phỏng bằng công cụ Simulink trong Matlab % %Copyright @ 2011 Nguyen Nhu Hieu %Date: May 10, 2011
  66. 62 %Faculty of Engineering Mechanics and Automation, University of Engineering and %Technology, Vietnam National University % clear; clc; num=10000; %num=1; sum1=0; sum3=0; R=1; S0=5; alpha=4; Q1=2*sqrt(2)/pi; Q2=0; Q3=2*sqrt(2)/(3*pi); p1=(Q1*sqrt(2*pi*S0))^2; p2=0; p3=(Q3*sqrt(2*pi*S0))^2; beta0=0.1; b1=1+alpha; om1=1; b2=1+alpha/16; om2=4; b3=1+alpha/81; om3=9; A1=1/2/R^2; A2=1/8/R^2; A3=1/18/R^2; %err=0.0001; for i=1:1:num Ran0=round(100000000*rand(1)); sim('ptvp1'); w1=yout(:,1); %w2=yout(:,2); w3=yout(:,3); MM1=mean(w1.*w1); MM3=mean(w3.*w3); sum1=sum1+MM1; sum3=sum3+MM3; end sum1=sum1/num sum3=sum3/num % THE END
  67. 63 Sơ đồ mô phỏng bằng công cụ Simulink trong Matlab