Luận án Phân tích đáp ứng của Profile cánh máy bay theo cách tiếp cận đối ngẫu

884.4  27.43  6.38  -1.16  s (m)   (kg/m3)  ca (kg.m 2/s)  ka1 (N.m)  ka2 (N.m)  0.6  1.225  0.036  6.833  9.967  ka3 (N.m)  ka4 (N.m)  ka5 (N.m)  cl cm 667.685  26.569  -5087.931  3.358  -0.635  Kp  Ki (1/s)  Kd (s)      0  Thay đổi để khảo sát  0      Ví dụ 5: Thiết diện cánh phi tuyến bậc 5, điều khiển vi phân (D control).  Các số liệu được cho trên bảng 10.  Bảng 10: Số liệu của ví dụ 5  a b (m)  mT (kg)  mW (kg)  Ia (kg.m 2)  -0.6847  0.135  12.387  2.049  0.0558  xa kh (N/m)  ch (kg/s)  cla cma 0.3314  2884.4  27.43  6.38  -1.16  s (m)   (kg/m3)  ca (kg.m 2/s)  ka1 (N.m)  ka2 (N.m)  0.6  1.225  0.036  6.833  9.967  ka3 (N.m)  ka4 (N.m)  ka5 (N.m)  cl cm 667.685  26.569  -5087.931  3.358  -0.635  Kp  Ki (1/s)  Kd (s)      0  0  Thay đổi để khảo sát       87  4.4.2. Tìm vận tốc tới hạn bằng phương pháp số Lời  giải  số  thu  được  bằng  cách  giải  hệ  phương  trình  vi  phân  phi  tuyến  (111).  Hệ  này  được  giải  bằng  hàm  ode45  trong  MATLAB.  Để  sử  dụng  hàm  ode45,  phương  trình  vi  phân  được  chuyển  về  hệ  phương  trình  vi  phân  cấp  1  (114), trong đó các tham số được tính từ (115), độ cứng phi tuyến được xác định  từ (112). Quy trình tính toán số sau được sử dụng để xác định vận tốc tới hạn và  tần số dao động LCO:  - Đầu tiên, một giá trị ban đầu của góc xoáy được cố định. Các giá trị ban  đầu khác       0 , 0 , 0h h a   được gán bằng 0.  - Thay đổi dần dần vận tốc dòng khí U. Với mỗi vận tốc dòng khí, phương  trình vi phân được giải. Trong các ví dụ thì khoảng thời gian giải phương trình  vi phân từ 0 đến 120s, là đủ dài để thu được dao động ổn định (trong tài liệu Li  vcs 2011 chỉ  tính đến 30s). Nếu biên độ dao động hội  tụ  thì vận  tốc dòng khí  nhỏ hơn vận tốc tới hạn. Khi đáp ứng bắt đầu hình thành dao động LCO (biên độ  dao động  tại  các  thời điểm cuối bằng  hoặc  lớn  hơn biên độ dao động  tại  thời  điểm ban đầu) thì vận tốc đạt tới vận tốc flutter. Khi đó vận tốc dòng, biên độ  dao động và tần số dao động sẽ được ghi lại.   - Tăng giá trị ban đầu của góc xoáy a(0) và lặp lại quá trình bên trên. Quá  trình trên thực tế tính toán cho thấy khi tăng giá trị ban đầu đủ lớn thì biên độ  dao động ổn định hầu như không tăng nữa. Đây có thể xác định là vị trí mà vận  tốc tới hạn có giá trị nhỏ nhất.   Trên thực tế quy trình tính toán như trên mất nhiều thời gian ở việc “tăng  dần” vận tốc dòng. Vì vậy để tăng tốc độ tính toán, quy trình cần làm từ bước  thô đến tinh. Cụ thể, đầu tiên bước vận tốc được lấy giá trị lớn để xác định vận  tốc gần với vận  tốc tới hạn. Sau đó  ta  thay đổi xuất phát điểm của vận  tốc và  giảm bước vận tốc. Quá trình tiếp diễn đến khi đạt được kết quả đủ độ chính xác  cần thiết.   88  Để đảm bảo độ tin cậy của tính toán số, kết quả tính vận tốc flutter trong ví  dụ 1 và 2 sẽ được so sánh với Li vcs 2011. Đoạn mã thể hiện việc tìm vận tốc  flutter bằng cách giải phương trình vi phân trong ví dụ 1 được cho trong phụ lục.  Hình 39, 40 thể hiện kết quả tính vận tốc tới hạn theo các điều kiện đầu. Hình 39  là kết quả của luận án còn hình 40 là kết quả trong Li vcs 2011 [53].  -0.35 -0.25 -0.15 -0.05 0.05 0.15 0.25 0.35 7.5 8 8.5 9 9.5 10 10.5 11 11.5 12 Vận tốc tới hạn , m/s a ( 0 ), r a d h(0)=0 h(0)=0.01 h(0)=0.05 Hình 39: Biên flutter của ví dụ 1 tính theo phương pháp số được đề cập trong luận án  Hình 40: Biên flutter của ví dụ 1 (lấy từ bài báo Li vcs 2011)   89  So sánh 2 hình vẽ ta thấy có sự trùng hợp rõ ràng và điều đó thể hiện độ tin  cậy của các tính toán số. Một số so sánh định lượng cũng được chỉ ra trong bảng  11, 12.  Bảng 11: So sánh một số vận tốc flutter trong ví dụ 1  h(0)  a(0)  Vận tốc tới hạn (Li vcs 2011)  Vận tốc tới hạn (Luận án)  0.01  0.1  7.92  7.955  0  ±0.2  7.9  7.94  Bảng 12: So sánh một số vận tốc flutter trong ví dụ 2  h(0)  a(0)  Vận tốc tới hạn (Li vcs 2011)  Vận tốc tới hạn (Luận án)  0.001  0.001  11.4  11.27  0.5  0.5  10.6  10.525  0.01  0.1  10.7  10.65  Kết quả so sánh  trên các bảng cho thấy sự phù hợp và khẳng định độ  tin  cậy của các tính toán số trong luận án.  4.5. Kết quả tính toán với ví dụ 1 Hình 41 cho thấy sự so sánh giữa các đường cong biên độ-vận tốc trong ví  dụ 1. Bảng 13 cho các kết quả định lượng so sánh các vận tốc tới hạn flutter.  7.5 8 8.5 9 9.5 10 10.5 11 11.5 12 12.5 0.09 0.1 0.11 0.12 0.13 0.14 Biên độ LCO (rad) V ận  t ố c  tớ i  h ạn  ( m /s ) TTH đối ngẫu  cải tiến 3 TTH kinh điển TTH đối ngẫu  thông thường - Dấu X: kết quả số - Nét đứt: đường  cong không ổn định Hình 41: Kết quả so sánh đường cong biên độ-vận tốc trong ví dụ 1   90  Bảng 13: So sánh các vận tốc tới hạn trong ví dụ 1  Đáp  ứng  tại  biên  độ  lớn  (độ  phi  tuyến lớn)  Vận  tốc  tới  hạn  nhỏ  nhất (m/s) Biên độ LCO  (rad)  Vận  tốc  tới  hạn  (m/s)  Tính toán số  0.1484  11.66  7.94  TTH kinh điển  12.25  (sai  số  5.06%)  7.9455  TTH đối ngẫu thông  thường  10.5666  (sai  số  9.38%)  7.9455  TTH  đối  ngẫu  cải  tiến 1  11.3431  (sai  số  2.72%)  7.9456  TTH  đối  ngẫu  cải  tiến 2  11.51  (sai  số  1.29%)  7.9457  TTH  đối  ngẫu  cải  tiến 3  11.591  (sai  số  0.59%)  7.9456  Kết quả cho thấy:    - Trên đường cong biên độ  - vận  tốc,  kỹ  thuật đối ngẫu đầy đủ 3  thành  phần cho kết quả gần với tính toán số nhất.    - Trên hình 41 và  trên bảng 13,  ta  thấy  tất cả các phương pháp đều  thể  hiện sự tồn tại của một vận tốc tới hạn flutter tối thiểu. Các phương pháp tuyến  tính hóa đều dự báo  tốt vận  tốc  tối  thiểu này và các kết quả không khác nhau  nhiều giữa các phương pháp (cột cuối cùng của bảng 13).     - Các phương pháp tuyến tính hóa đối ngẫu cải tiến đều cho kết quả cải  thiện so với phương pháp tuyến tính hóa kinh điển trên miền phi tuyến (biên độ  dao động lớn) như thấy trên cột 2 của bảng 13.  4.6. Kết quả tính toán với ví dụ 2 Hình 42 cho thấy sự so sánh giữa các đường cong biên độ-vận tốc trong ví  dụ 1. Bảng 14 cho các kết quả định lượng so sánh các vận tốc tới hạn flutter.   91  10.5 10.6 10.7 10.8 10.9 11 11.1 11.2 11.3 11.4 11.5 0.09 0.11 0.13 0.15 0.17 Biên độ LCO (rad) V ận  t ố c  tớ i  h ạn  ( m /s ) TTH đối ngẫu  cải tiến 3 TTH kinh điển TTH đối ngẫu  thông thường - Dấu X: kết quả số - Nét đứt: đường  cong không ổn định Hình 42: Kết quả so sánh đường cong biên độ-vận tốc trong ví dụ 2  Bảng 14: So sánh các vận tốc tới hạn flutter trong ví dụ 2  Đáp  ứng  tại  biên  độ  lớn  (độ  phi  tuyến lớn)  Vận  tốc  tới  hạn  nhỏ  nhất (m/s) Biên độ LCO  (rad)  Vận  tốc  tới  hạn  (m/s)  Tính toán số  0.1743  11.29  10.525  TTH kinh điển  11.4404  (sai  số  1.33%)  10.5246  TTH đối ngẫu thông  thường  11.0654  (sai  số  1.99%)  10.5247  TTH  đối  ngẫu  cải  tiến 1  11.2457  (sai  số  0.39%)  10.5246  TTH  đối  ngẫu  cải  tiến 2  11.2782  (sai  số  0.1%)  10.5247  TTH  đối  ngẫu  cải  tiến 3  11.2939  (sai  số  0.03%)  10.5246   92  Kết quả tiếp tục cho thấy một số kết luận đã thu được ở ví dụ 1:   - Kỹ thuật  tuyến tính hóa đối ngẫu đầy đủ 3 thành phần cho kết quả gần  với tính toán số nhất trong miền phi tuyến lớn (biên độ LCO lớn).  - Các  kỹ  thuật  tuyến  tính  hóa đều  tìm  ra được vận  tốc  tới hạn  flutter  tối  thiểu và sai khác rất ít so với tính toán số.   4.7. Kết quả tính toán với ví dụ 3 Trong ví dụ 3 này, ta sẽ tính hai trường hợp tham số Kp khác nhau để đánh  giá ảnh hưởng của điều khiển tỷ lệ (P-control) đối với hệ trong ví dụ 1. Ta xét 2  trường hợp Kp=-0.5 và Kp=0.5. Các kết quả so sánh được cho trên hình 43, 44 và  bảng 15, 16.  9.4 9.6 9.8 10 10.2 10.4 10.6 10.8 11 0.07 0.075 0.08 0.085 0.09 0.095 Biên độ LCO (rad) V ận  t ố c  tớ i  h ạn  ( m /s ) TTH đối ngẫu  cải tiến 3 TTH kinh điển TTH đối ngẫu  thông thường - Dấu X: kết quả số - Nét đứt: đường  cong không ổn định Hình 43: Kết quả so sánh đường cong biên độ-vận tốc trong ví dụ 3, Kp=-0.5  Bảng 15: So sánh các vận tốc tới hạn flutter trong ví dụ 3, Kp=-0.5  Đáp  ứng  tại  biên  độ  lớn  (độ  phi  tuyến lớn)  Vận  tốc  tới  hạn  nhỏ  nhất (m/s) Biên độ LCO  (rad)  Vận  tốc  tới  hạn  (m/s)   93  Tính toán số  0.0972  10.465  9.415  TTH kinh điển  10.9671  (sai  số  4.8%)  9.4207  TTH đối ngẫu thông  thường  9.8347  (sai  số  6.02%)  9.421  TTH  đối  ngẫu  cải  tiến 1  10.1836  (sai  số  2.69%)  9.4209  TTH  đối  ngẫu  cải  tiến 2  10.2784  (sai  số  1.78%)  9.4207  TTH  đối  ngẫu  cải  tiến 3  10.3285  (sai  số  1.30%)  9.421  6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5 15.5 0.1 0.12 0.14 0.16 0.18 0.2 0.22 Biên độ LCO (rad) V ận  t ố c  tớ i  h ạn  ( m /s ) TTH đối ngẫu  cải tiến 3 TTH kinh điển TTH đối ngẫu  thông thường - Dấu X: kết quả số - Nét đứt: đường  cong không ổn định Hình 44: Kết quả so sánh đường cong biên độ-vận tốc trong ví dụ 3, Kp=0.5  Bảng 16: So sánh các vận tốc tới hạn flutter trong ví dụ 3, Kp=0.5  Đáp  ứng  tại  biên  độ  lớn  (độ  phi  tuyến lớn)  Vận  tốc  tới  hạn  nhỏ  nhất (m/s) Biên độ LCO  (rad)  Vận  tốc  tới  hạn  (m/s)   94  Tính toán số  0.2255  14.745  6.985  TTH kinh điển  15.1676  (sai  số  2.87%)  6.9877  TTH đối ngẫu thông  thường  13.9926  (sai  số  5.1%)  6.988  TTH  đối  ngẫu  cải  tiến 1  14.5797  (sai  số  1.12%)  6.9877  TTH  đối  ngẫu  cải  tiến 2  14.6995  (sai  số  0.31%)  6.9879  TTH  đối  ngẫu  cải  tiến 3  14.7562  (sai  số  0.08%)  6.9877  Qua các kết quả tính toán, ta có một số nhận xét sau:  - Trong các trường hợp điều khiển tỷ lệ (P-control), kỹ thuật tuyến tính hóa  đối ngẫu cải tiến 3 (sử dụng 3 thành phần) có độ chính xác tốt nhất khi tính toán  vận tốc tới hạn ở miền có biên độ LCO lớn (độ phi tuyến cao).  - Các trường hợp tuyến tính hóa đều cho kết quả vận tốc tới hạn nhỏ nhất là  gần như nhau và gần với kết quả tính toán số.  - So sánh bảng 13 và bảng 15 ta thấy hệ số tỷ lệ Kp âm sẽ làm tăng vận tốc  giới hạn tối thiểu (từ 7.94m/s lên 9.415m/s) và giảm vận tốc tới hạn ở miền biên  độ lớn (từ 11.66m/s xuống 10.465m/s). Ngược lại, so sánh bảng 13 và bảng 16  ta thấy hệ số tỷ lệ Kp dương sẽ làm giảm vận tốc giới hạn tối thiểu (từ 7.94m/s  xuống 6.985m/s) và  tăng vận  tốc  tới hạn ở miền biên độ  lớn  (từ 11.66m/s  lên  14.745m/s). Xét về mặt hiệu quả điều khiển thì hệ số Kp âm là tốt hơn.  4.8. Kết quả tính toán với ví dụ 4 Trong ví dụ 4 này, ta sẽ tính hai trường hợp tham số Ki khác nhau để đánh  giá ảnh hưởng của điều khiển tích phân (I-control) đối với hệ trong ví dụ 1. Ta  xét 2 trường hợp Ki=-2 (1/s) và Ki=2 (1/s). Trường hợp Ki âm, hệ ổn định và các   95  kết quả so sánh được cho trên hình 45 và bảng 17. Trường hợp Ki dương, hệ mất  ổn định. Phương pháp số thể hiện điều này ở việc biên độ dao động tiến ra vô  cùng. Còn các phương pháp tuyến tính hóa thể hiện điều này ở việc tồn tại giá trị  riêng của bài toán giá trị riêng (37) có phần thực dương. Các kết quả này được  thể hiện trên hình 46 và bảng 18.  7 7.5 8 8.5 9 9.5 10 0.1 0.11 0.12 0.13 0.14 0.15 Biên độ LCO (rad) V ận  t ố c  tớ i  h ạn  ( m /s ) TTH đối ngẫu  cải tiến 3 TTH kinh điển TTH đối ngẫu  thông thường - Dấu X: kết quả số - Nét đứt: đường  cong không ổn định Hình 45: Kết quả so sánh đường cong biên độ-vận tốc trong ví dụ 4, Ki=-2 (1/s)  Bảng 17: So sánh các vận tốc tới hạn flutter trong ví dụ 4, Ki=-2 (1/s)  Đáp ứng  tại biên độ  lớn  (độ  phi  tuyến lớn)  Vận  tốc  tới  hạn  nhỏ  nhất (m/s) Biên độ LCO  (rad)  Vận  tốc  tới  hạn  (m/s)  Tính toán số  0.1595  9.385  7.055  TTH kinh điển  9.5412  (sai  số  1.66%)  7.058  TTH đối ngẫu thông  thường  8.848  (sai  số  5.72%)  7.079  TTH  đối  ngẫu  cải  tiến 1  9.1912  (sai  số  2.06%)  7.0577   96  TTH  đối  ngẫu  cải  tiến 2  9.2598  (sai  số  1.33%)  7.0578  TTH  đối  ngẫu  cải  tiến 3  9.2924  (sai  số  0.99%)  7.0577  -0.01 -0.008 -0.006 -0.004 -0.002 0 0.002 0.004 0.006 0.008 0.01 0 20 40 60 80 100 120 Thời gian (s) G ó c  d ao  đ ộ n g   a  ( ra d ) Hình 46: Sự phân kỳ của dao động xoắn trong ví dụ 4 khi Ki=2 (1/s), h(0)=0(m),  a(0)=0.0125(rad), u=2.6 (m/s)  Bảng 18: Các giá trị riêng tính theo các phương pháp tuyến tính hóa, ví dụ 4, Ki=2 (1/s)    TTH kinh  điển  TTH đối  ngẫu thông  thường  TTH đối  ngẫu cải  tiến 1  TTH đối  ngẫu cải  tiến 2  TTH đối  ngẫu cải  tiến 3  Các giá trị  riêng của  (37) ứng với  trường hợp  biên độ rất  nhỏ =0.01  rad  -2.7795  +14.1352i  -2.7843  +14.1308i  -2.7819  +14.1330i  -2.7818  +14.1331i  -2.7817  +14.1331i  -2.7795  -14.1352i  -2.7843  -14.1308i  -2.7819  -14.1330i  -2.7818  -14.1331i  -2.7817  -14.1331i  14.1336i  14.1300i  14.1318i  14.1319i  14.1319i  -14.1336i  -14.1300i  -14.1318i  -14.1319i  -14.1319i  1.6001 1.6078 1.6039 1.6037 1.6037  97  Qua các kết quả tính toán, ta có một số nhận xét sau:  - Trong các trường hợp điều khiển tích phân (I-control), kỹ thuật tuyến tính  hóa đối ngẫu cải  tiến 3  (sử dụng 3  thành phần) vẫn cho  thấy độ chính xác  tốt  nhất khi tính toán vận tốc tới hạn ở miền có biên độ LCO lớn (độ phi tuyến cao).  - Các trường hợp tuyến tính hóa đều cho kết quả vận tốc tới hạn nhỏ nhất là  gần như nhau và gần với kết quả tính toán số.  - So sánh bảng 13 và bảng 17  ta  thấy hệ số  tích phân Ki âm sẽ  làm giảm  vận tốc giới hạn tối thiểu (từ 7.94m/s xuống 7.055m/s) và giảm vận tốc tới hạn ở  miền biên độ lớn (từ 11.66m/s xuống 9.385m/s).   - Hệ số tích phân Ki dương sẽ làm hệ mất ổn định được thể hiện trên cả các  tính toán số và lời giải tuyến tính hóa. Như thấy trên bảng 18, giá trị riêng là số  thực dương sẽ làm cho dao động tiến ra vô cùng ngay cả với biên độ rất nhỏ.  - Như  vậy,  lời giải  tuyến  tính  hóa  không chỉ  xác  định được  các đáp  ứng  LCO mà  còn chỉ  ra  được cả  các  trường hợp  không  tồn  tại  LCO, biên độ dao  động  tiến ra vô cùng. Ngoài ra, qua ví dụ này ta  thấy điều khiển  tích phân  (I- control) trong trường hợp này là không có hiệu quả (hoặc làm giảm vận tốc tới  hạn tối thiểu, hoặc làm hệ mất ổn định).  4.9. Kết quả tính toán với ví dụ 5 Trong ví dụ 5 này, ta sẽ tính hai trường hợp tham số Kd khác nhau để đánh  giá ảnh hưởng của điều khiển vi phân (D-control) đối với hệ trong ví dụ 1. Ta  xét 2 trường hợp Kd=-0.5 (s) và Kd=0.5 (s). Các kết quả so sánh được cho trên  hình 47, 48 và bảng 19,20.   98  2.8 2.82 2.84 2.86 2.88 2.9 0.07 0.075 0.08 0.085 0.09 Biên độ LCO (rad) V ận  t ố c  tớ i  h ạn  ( m /s ) TTH đối ngẫu  cải tiến 3 TTH kinh điển TTH đối ngẫu  thông thường - Dấu X: kết quả số - Nét đứt: đường  cong không ổn định Hình 47: Kết quả so sánh đường cong biên độ-vận tốc trong ví dụ 5, Kd=-0.5 (s)  Bảng 19: So sánh các vận tốc tới hạn flutter trong ví dụ 5, Kd=-0.5 (s)  Đáp ứng  tại biên độ  lớn  (độ  phi  tuyến lớn)  Vận  tốc  tới  hạn  nhỏ  nhất (m/s) Biên độ LCO  (rad)  Vận  tốc  tới  hạn  (m/s)  Tính toán số  0.0932  2.87  2.805  TTH kinh điển  2.8857  (sai  số  0.55%)  2.8063  TTH đối ngẫu thông  thường  2.841  (sai  số  1.01%)  2.8063  TTH  đối  ngẫu  cải  tiến 1  2.8606  (sai  số  0.33%)  2.8063  TTH  đối  ngẫu  cải  tiến 2  2.8649  (sai  số  0.18%)  2.8063  TTH  đối  ngẫu  cải  tiến 3  2.867  (sai  số  0.1%)  2.8063   99  9.45 9.5 9.55 9.6 9.65 9.7 9.75 9.8 9.85 0.05 0.055 0.06 0.065 Biên độ LCO (rad) V ận  t ố c  tớ i  h ạn  ( m /s ) TTH đối ngẫu  cải tiến 3 TTH kinh điển TTH đối ngẫu  thông thường - Dấu X: kết quả số - Nét đứt: đường  cong không ổn định Hình 48: Kết quả so sánh đường cong biên độ-vận tốc trong ví dụ 5, Kd=0.5 (s)  Bảng 20: So sánh các vận tốc flutter trong ví dụ 5, Kd=0.5 (s)  Đáp ứng  tại biên độ  lớn  (độ  phi  tuyến lớn)  Vận  tốc  tới  hạn  nhỏ  nhất (m/s) Biên độ LCO  (rad)  Vận  tốc  tới  hạn  (m/s)  Tính toán số  0.0655  9.78  9.5  TTH kinh điển  9.8426  (sai  số  0.64%)  9.4846  TTH đối ngẫu thông  thường  9.5608  (sai  số  2.24%)  9.4846  TTH  đối  ngẫu  cải  tiến 1  9.6472  (sai  số  1.36%)  9.4846  TTH  đối  ngẫu  cải  tiến 2  9.6703  (sai  số  1.12%)  9.4846  TTH  đối  ngẫu  cải  tiến 3  9.6826  (sai  số  1%)  9.4846   100  Qua các kết quả tính toán, ta có một số nhận xét sau:  - Trong các trường hợp điều khiển vi phân (D-control), kỹ thuật tuyến tính  hóa đối ngẫu cải tiến 3 (sử dụng 3 thành phần) có độ chính xác tương đối tốt khi  tính toán vận tốc tới hạn ở miền có biên độ LCO lớn (độ phi tuyến cao).  - Các trường hợp tuyến tính hóa đều cho kết quả vận tốc tới hạn nhỏ nhất là  gần như nhau và gần với kết quả tính toán số.   - Đối với dao động có biên độ nhỏ hệ phi tuyến sẽ rất gần với hệ tuyến tính  (tính phi tuyến yếu) khi đó phương pháp TTH kinh điển cho kết quả tốt hơn, tuy  nhiên sai số của các phương pháp TTH đều nhỏ do tính phi tuyến yếu.  - So sánh bảng 13 và bảng 19 ta thấy hệ số tỷ lệ Kd âm sẽ làm giảm vận tốc  giới hạn tối thiểu (từ 7.94m/s xuống 2.805m/s) và giảm vận tốc tới hạn ở miền  biên độ  lớn (từ 11.66m/s xuống 2.87m/s). Ngược  lại,  so sánh bảng 13 và bảng  20 ta thấy hệ số tỷ lệ Kd dương sẽ làm tăng vận tốc giới hạn tối thiểu (từ 7.94m/s  lên  9.5m/s)  và  giảm  vận  tốc  tới  hạn  ở  miền  biên  độ  lớn  (từ  11.66m/s  xuống  9.78m/s). Như vậy trong cả hai trường hợp Kd dương và âm thì miền biến thiên  của vận tốc giới hạn đều thu hẹp lại, tức là ít phụ thuộc vào biên độ dao động và  tính phi tuyến ít đi.  - Xét về mặt hiệu quả điều khiển thì hệ số Kd dương sẽ tốt hơn (tăng vận  tốc giới hạn tối thiểu).  Kết luận chương 4 Trong chương này đã khảo sát mô hình khí động phi tuyến của một  thiết  diện cánh dựa  trên mô hình  thí nghiệm được nhiều nhà nghiên cứu quốc tế sử  dụng. Trong mô hình được trình bày, độ cứng xoắn ka được xét là một hàm của  a dưới dạng đa thức bậc 4. Trong luận án này quan tâm tới thuật toán điều khiển  PID. Đây  là  thuật  toán điều khiển đơn giản nhất và được  sử dụng nhiều nhất.  Các kết quả của các kỹ thuật đối ngẫu được áp dụng vào bài toán ổn định flutter  của thiết diện cánh cụ thể. Bằng kỹ thuật tuyến tính thu được vận tốc U là hàm  của biên độ dao động xoắn A. Từ đó sẽ vẽ được đường đặc trưng biên độ A với   101  vận tốc tới hạn U và cũng cho ta mối quan hệ giữa tần số với vận tốc tới hạn. Để  khảo  sát  cụ  thể  đã xét  2  ví dụ  thiết  diện  cánh  phi  tuyến  đến bậc  3  và  bậc  5,  không có điều khiển, 2 ví dụ thiết diện cánh phi tuyến đến bậc 5, có điều khiển  tích phân và  điều khiển vi phân. Đã xác định vận tốc tới hạn bằng phương pháp  số giải hệ phương trình vi phân phi tuyến bằng hàm ode45 trong MATLAB. Đã  xác định quy trình tính toán số để xác định vận tốc tới hạn và tần số dao động  LCO. Quy  trình cần  làm  từ  bước  thô đến  tinh. Cụ  thể, đầu  tiên bước  vận  tốc  được lấy giá trị lớn để xác định vận tốc gần với vận tốc tới hạn. Sau đó thay đổi  xuất phát điểm của vận tốc và giảm bước vận tốc. Quá trình tiếp diễn đến khi đạt  được kết quả đủ độ chính xác cần thiết. Để đánh giá độ tin cậy của tính toán số,  kết quả tính vận tốc flutter trong ví dụ 1 và 2 được so sánh với Li vcs 2011. Kết  quả so sánh trên đồ thị và bảng số liệu cho thấy có sự phù hợp rõ ràng và điều đó  thể hiện độ tin cậy của các tính toán số trong luận án.   Kết quả phân tích cũng cho thấy về tổng thể các kỹ thuật tuyến tính hóa và  phương pháp số đều không có nhiều sự khác biệt về định tính trên đường cong  biên độ - vận tốc. Đối với dao động nhỏ hệ phi tuyến sẽ rất gần với hệ tuyến tính  (tính phi tuyến yếu) khi đó phương pháp TTH kinh điển cho kết quả tốt hơn, tuy  nhiên sai số của các phương pháp TTH đều nhỏ do tính phi tuyến yếu. Đối với  dao động có biên độ lớn (tính phi tuyến lớn) các kỹ thuật TTH đối ngẫu cho kết  quả có sai số nhỏ hơn kỹ thuật TTH kinh điển. Ngoài ra, nhìn chung kỹ thuật đối  ngẫu đầy đủ 3 thành phần cho kết quả gần với tính toán số nhất.    102  KẾT LUẬN VÀ KIẾN NGHỊ A. Trong luận án đã tiến hành nghiên cứu vấn đề sau đây.   - Đã tiến hành phân tích tổng quan một số kết quả chính ở trên thế giới và  trong nước về vấn đề ổn định flutter của thiết diện cánh chịu lực khí động.  - Trên cơ sở lý thuyết khí động học và các số liệu thực nghiệm đã có xây  dựng  mô  hình  thiết  diện  cánh  máy  bay  2  chiều  uốn-xoắn  chuyển  động  trong  dòng khí không nén được, dựa trên mô hình thí nghiệm được nhiều nhà nghiên  cứu quốc tế sử dụng. Trong mô hình được trình bày, độ cứng xoắn ka được xét  là một hàm của a dưới dạng đa thức bậc 5. Phương trình phi tuyến thu được từ  mô hình được dùng để phân tích đáp ứng cũng như các hiện tượng flutter.   - Đã áp dụng công cụ CFD và CSM để tối ưu hóa hình dạng cánh máy bay  đã được áp dụng. Phương pháp này được sử dụng cho các thiết kế của các cánh  máy bay UAV, bay ở số Reynolds thấp. Một bài toán giảm thiểu lực cản, trong  khi vẫn đạt yêu cầu về lực nâng được giải quyết bằng phương pháp SQP. Thiết  diện  cánh  Eppler  66  được  chọn  làm  cánh  máy  bay  ban  đầu.  Mô  phỏng  cũng  được thực hiện bằng cách sử dụng ANSYS Workbench phiên bản 16.0 để chứng  minh hiệu quả của cánh máy bay tối ưu. Trong trường hợp này, kết quả cho thấy  các cánh máy bay tối ưu đạt được sự giảm 20% lực cản so với các cánh máy bay  ban đầu và vẫn có thể để đảm bảo các yêu cầu lực nâng tối thiểu.  - Sau khi thiết lập phương trình dao động hai bậc tự do của thiết diện cánh  đã trình bày một số phương pháp giải tích, đặc biệt đã phát triển tiêu chuẩn đối  ngẫu có trọng số cho phương pháp tuyến tính hóa tương đương hệ dao động phi  tuyến tuần hoàn và ngẫu nhiên, trong đó tiêu chuẩn tuyến tính hóa kinh điển là  trường hợp riêng khi trọng số bằng không.  - Áp dụng cho bài toán ổn định flutter của thiết diện cánh đã thu được vận  tốc U là hàm của biên độ dao động xoắn A. Từ đó sẽ vẽ được đường đặc trưng  biên độ A với vận tốc tới hạn U và cũng cho ta mối quan hệ giữa tần số với vận  tốc tới hạn. Để khảo sát cụ thể đã xét 2 ví dụ thiết diện cánh phi tuyến đến bậc 3   103  và bậc 5, không có điều khiển, 2 ví dụ thiết diện cánh phi tuyến đến bậc 5, có  điều khiển tích phân và  điều khiển vi phân.   - Đã xác định vận tốc tới hạn bằng phương pháp số giải hệ phương trình vi  phân  phi  tuyến  bằng  hàm ode45  trong  MATLAB.  Đã xác  định quy  trình  tính  toán số để xác định vận tốc tới hạn và tần số dao động LCO. Quy trình cần làm  từ bước thô đến tinh. Cụ thể, đầu tiên bước vận tốc được lấy giá trị lớn để xác  định vận tốc gần với vận tốc tới hạn. Sau đó thay đổi xuất phát điểm của vận tốc  và giảm bước vận tốc. Quá trình tiếp diễn đến khi đạt được kết quả đủ độ chính  xác cần thiết.   - Đã đánh giá độ tin cậy của tính toán số của luận án, kết quả tính vận tốc  flutter trong ví dụ 1 và 2 được so sánh với Li vcs 2011 [53]. Kết quả so sánh trên  đồ thị và bảng số liệu cho thấy có sự phù hợp rõ ràng và điều đó thể hiện độ tin  cậy của các tính toán số trong luận án.   - Kết quả phân tích cũng cho thấy về định tính tổng thể các kỹ thuật tuyến  tính hóa và phương pháp số đều không có nhiều sự khác biệt  trên đường cong  biên độ - vận tốc. Đối với dao động có biên độ nhỏ hệ phi tuyến sẽ rất gần với  hệ  tuyến tính (tính phi  tuyến yếu) khi đó phương pháp TTH kinh điển cho kết  quả  tốt hơn,  tuy  nhiên  sai  số  của  các  phương  pháp  TTH  đều  nhỏ  do  tính  phi  tuyến  yếu. Đối  với dao động có biên độ  lớn  (tính  phi  tuyến  lớn)  các  kỹ  thuật  TTH đối ngẫu cho kết quả có sai số nhỏ hơn kỹ thuật TTH kinh điển. Ngoài ra,  kỹ thuật đối ngẫu đầy đủ 3 thành phần thường cho kết quả gần với tính toán số  nhất.            B. Những kết quả mới của luận án bao gồm:  -  Lần đầu  tiên áp  dụng và  phát  triển  sáng  tạo cách  tiếp  cận  đối  ngẫu  để  phân  tích các hiện  tượng dao động flutter xuất hiện  trong  thiết diện cánh chịu  lực khí động.  - Đối với tiêu chuẩn tuyến tính hóa có trọng số đã giải quyết vấn đề chọn  giá trị trọng số thông qua nghiên cứu đề xuất 3 cách lựa chọn tương ứng với 3  cải tiến. Đã tiến hành khảo sát tần số dao động riêng trong hệ dao động bảo toàn   104  bậc cao và hệ Duffing chịu  tải trọng ngẫu nhiên cho  thấy các cải  tiến này đều  cho kết quả tốt hơn so với phương pháp tuyến tính hóa kinh điển.  - Áp dụng cho bài toán ổn định flutter của thiết diện cánh, với các dạng  phi  tuyến đa thức hay gặp và tính phí tuyến không quá nhỏ, cho thấy các tiêu chuẩn  đối ngẫu cải tiến đều cho kết quả chính xác hơn tiêu chuẩn kinh điển. Ngoài ra  các  tiêu chuẩn đối ngẫu đầy đủ 3  thành phần  thường cho  kết quả  tốt hơn  tiêu  chuẩn  đối ngẫu  2  thành phần  trong các ví dụ khảo  sát. Điều đó cho  thấy  tiêu  chuẩn đối ngẫu đầy đủ 3  thành phần có thể áp dụng tốt cho các hệ có tính phi  tuyến rõ rệt.  Hướng nghiên cứu  tiếp  theo sau  luận án có  thể tập  trung vào một số việc  sau đây:  - Đối với cách tiếp cận đối ngẫu trong phương pháp tuyến tính hóa tương  đương cần khảo sát nhiều hệ dao động phi tuyến khác nhau để phát hiện các hiện  tượng cũng như các khả năng và giới hạn áp dụng của tiêu chuẩn tương đương  đối ngẫu.  - Nghiên  cứu  mở  rộng bài  toán ổn  định  thiết  diện  cánh 2 chiều  sang  bài  toán ổn định cánh có xét đến chiều dài cánh. Đây là bài toán rất phức tạp nhưng  rất cần thiết cho kỹ thuật hàng không. Nghiên cứu thêm các phương pháp CFD,  CSM để áp dụng cho bài toán ổn định cánh.   105  DANH SÁCH CÔNG TRÌNH ĐàĐƯỢC CÔNG BỐ CỦA LUẬN ÁN  1. Nguyễn Đông Anh, Nguyễn Minh Triết, Mở Rộng Tiêu Chuẩn Đối Ngẫu Cho Các Hệ Phi Tuyến Dao Động Tuần Hoàn, Tạp chí Khoa học Giáo dục Kỹ thuật,  Trường Đại học Sư phạm Kỹ thuật TPHCM, 2015, p.03:07.  2.  Nguyen  Minh  Triet,  Nguyen  Ngoc  Viet,  Pham  Manh  Thang,  Aerodynamic Analysis of Aircraft Wing,  VNU  Journal  of  Science,  Natural  Sciences  and  Technology, 2015, p.68:75.  3. Nguyen Minh  Triet, Extension of dual equivalent linearization technique to flutter analysis of two dimensional nonlinear airfoils,  Vietnam  Journal  of  Mechanics, vol. 37, N3, 2015, p.217:230.  4.  Nguyen  Minh  Triet,  A Full Dual Mean Square Error Criterion For The Equivalent Linearization, Journal of Science and Technology, 2015, p.557:562.  5. Nguyen Minh Triet, M.T. Pham, M. C.  Vu, D.A. Nguyen - "Design wireless control system for aircraft model "  Proceedings  of  the  3rd  International  Conference  on  Engineering  Mechanics  and  Automation,    ICEMA3,  2014,   p.283:286.  6.  Minh  Triet  Nguyen,  Ngoc  Viet  Nguyen,  Van  Manh  Hoang,  Manh  Thang  Pham - Aerodynamic shape optimization of airfoil using SQP method  -  Tuyển  tập Hội nghị Khoa học toàn quốc Cơ học Vật rắn biến dạng lần thứ XII, 2015,  p.1442:1449.  7.  Minh  Triet  Nguyen,  Van  Long  Nguyen,  Ngoc  Viet  Nguyen,  Ngoc Linh Nguyen, Manh Thang Pham “A Study on Low-Speed Wind Tunnel – Theory and Experiment”    Proceedings  of  the  4rd  International  Conference  on  Engineering Mechanics and Automation - ICEMA4, 2016.  8.  Minh  Triet  Nguyen,  Ngoc  Viet  Nguyen,  Van  Manh  Hoang,  Manh  Thang  Pham “Aerodynamic analysis and experiment of an airfoil in a low speed wind tunnel”. Proceedings  of  the  4rd  International  Conference  on  Engineering  Mechanics and Automation - ICEMA4, 2016.   106   107  TÀI LIỆU THAM KHẢO  Tiếng Việt 1. Hoàng  Thị  Bích  Ngọc,  Đinh  Văn  Phong,  Nguyễn  Hồng  Sơn  (2009),  Nghiên cứu hiện tượng đàn hồi cánh dưới tác dụng của lực khí động, Tuyển tập  công trình Hội nghị Cơ học toàn quốc, Hà Nội, Tập 1, tr. 485-494.  2. Lã Hải Dũng  (2005), Nghiên cứu quá  trình  flutter uốn  - xoắn cánh máy  bay, Kỷ yếu Hội thảo toàn quốc cơ học và khí cụ bay có điều khiển lần thứ nhất,  NXB ĐHQG Hà Nội, tr.238-244  3. Nguyễn  Ngọc  Linh  (2015),  Luận  án  tiến  sĩ:  Phân  tích  dao  động  ngẫu  nhiên phi tuyến bằng phương pháp tuyến tính hóa tương đương, Viện Cơ học.  4. Nguyễn Văn Khang (1998), Dao động kỹ thuật, Nhà xuất bản Khoa học  và Kỹ thuật, Hà Nội.  5. Nguyễn Văn Khang, Trần Ngọc An (2014), Flutter instability analysis of  bridge decks using the step-by-step method. Vietnam J. of Mechanics, vol 36, N  1, p.1-11.  6. Phạm Duy Hòa, Nguyễn Văn Mỹ, Phan Thanh Hoàng (2014), Nghiên cứu  kiểm soát ổn định flutter của kết cấu cầu hệ treo bằng khe slot, Tuyển tập công  trình  Hội nghị  cơ  học  toàn  quốc,  Tập 1, Nhà  xuất  bản Khoa học  tự  nhiên  và  công nghệ, ISBN 978-604-913-233-9.  7. Trần  Ngọc  An  (2014),  Tính  Toán Ổn Định Khí  Động  Flutter  Của Dầm  Chủ Trong Kết Cấu Cầu Hệ Dây Bằng Phương Pháp Bước Lặp, Luận án tiến sĩ  Cơ học, Hà Nội.  8.  Trần  Thế  Văn  (2012),  Luận  án  tiến  sĩ:  Nghiên  cứu  ổn  định  của  tấm  composite lớp chịu tải trọng khí động, Học viện kỹ thuật quân sự.  Tiếng Anh 9. Abdelkader M., Shaqura M., Claudel C.G. and Gueaieb W. (2013), “A UAV  based system for real time flash flood monitoring in desert environments using   108  Lagrangian microsensors”, Proceedings of the 2013 International Conference on  Unmanned Aircraft Systems, Atlanta, USA, May.  10. Abdelkefi. A, R. Vasconcellos, F.D. Marques, M.R. Hajj (2012), Modeling  and  identification  of  freeplay  nonlinearity,  Journal  of  Sound  and  Vibration,  Volume 331, Issue 8, , Pages 1898-1907.  11. Allen  CB,  Taylor  NV,  Fenwick  CL,  Gaitonde  AL,  Jones  DP  (2005).  A  comparison of full non-linear and reduced order aerodynamic models in control  law design  for  active  flutter  suppression.  Int J Numer Meth Eng;64(12):1628– 48.  12. Anderson.  J.  D.  (2010),  “Fundamentals  of  Aero-dynamics”,  McGraw-Hill  Science.  13. Anh  N.D.,  Hieu  N.N.,  Linh  N.N.  (2012a),  A  dual  criterion  of  equivalent  linearization method for nonlinear systems subjected to random excitation, Acta  Mechanica, 223:645-654.  14. Anh N.D., Nguyen N.X., Hoa L.T. (2013), Design of three-element dynamic  vibration absorber for damped linear structures, Journal of Sound and Vibration,  DOI:10.1016/j.jsv.2013.03.032.  15. Anh  N.D.,  Zakovorotny  V.L.,  Hieu  N.N.,  Diep  D.V.  (2012b),  A  dual  criterion of stochastic linearization method for multi-degree-of-freedom systems  subjected to random excitation, Acta Mechanica, 223:2667-2684.  16. Anh,  N.D.,  Hung,  L.X.,  Viet,  L.D.  (2012c)  Dual  approach  to  local  mean  square  error  criterion  for  stochastic  equivalent  linearization,  Acta  Mech.  DOI  10.1007/s00707-012-0751-8   17. Anh. N.D. (2010), Duality in the analysis of responses to nonlinear systems,  Vietnam Journal of Mechanics, VAST 32:263-266.  18. Antoniou.A and Lu. W.S.  (2003), Optimization: Methods, Algorithms, and  Applications, Kluwer Academic.   109  19. Asami T., Nishihara O. (1999), Analytical and experimental evaluation of an  air  damped dynamic  vibration  absorber:  design optimizations  of  three-element  type model, Journal of Vibration and Acoustics, 121:334-342.  20. Brat. W.V. (1971), “Flight test measurement of exterior turbulent boundary  layer pressure fluctuations on Boeing model 737 airplane”, Journal of Sound and  Vibration, vol.14, pp. 439-457.  21. Caughey  TK.  (1963):  Equivalent  linearization  techniques.  Journal  of  the  Acoustical Society of America; 35:1706–17112.   22. Chen Feixin, Liu  Jike, Chen Yanmao  (2013),  Flutter  analysis  of  an  airfoil  with nonlinear damping using equivalent  linearization,  Journal of Aeronautics,  DOI:10.1016/j.cja.2013.07.020.  23. Chen FX, Chen YM, Liu JK. (2012) Equivalent linearization method for the  flutter system of an airfoil with multiple nonlinearities. Commun Nonlinear Sci  Numer Simul;17(12):4529–35.  24. Chen  YM,  Liu  JK,  Meng  G.  (2011)  Equivalent  damping  of  aeroelastic  system of an airfoil with cubic stiffness. J Fluids Struct;27(8):1447–54.  25. Chung  KW,  Chan  CL,  Lee  BHK.  (2007)  Bifurcation  analysis  of  a  two- degree-of-freedom aeroelastic system with freeplay structural nonlinearity by a  perturbation-incremental method. J Sound Vib;299(3):520–39.  26. Collar, A.R. (1978) The first fifty years of aeroelasticity, Aerospace, 2–20.  27. Crandall  S.H.:  (2006),  "A  half-century  of  stochastic  equivalent  linearization", Struct. Control Health Monit., 13, 27–40.  28. Ding  Qian,  Wang  Dong-Li  (2006),  The  flutter  of  an  airfoil  with  cubic  structural and aerodynamic non-linearities, Aerospace Science and Technology,  10:427-434.  29. Djayapertapa  L,  Allen  CB  (2001a).  Aeroservoelastic  simulation  by  time  marching. Aeronaut J;105(1054):667–78.   110  30. Djayapertapa L, Allen CB, Fiddes SP (2001b). Two-dimensional  transonic  aeroservoelastic  computations  in  the  time  domain.  Int  J  Numer  Meth  Eng;52(12):1355–77.  31. Dowell, E.H., Clark, R., Cox, D., Curtiss, H.C., Edwards, J.W., Hall, K.C.,  Peters,  D.A.,  Scanlan,  R.,  Simiu,  E.,  Sisto,  F.,  Strganac,  T.W.  and  Tang.D  (2015),  A  modern  course  in  aeroelasticity.  5th  ed.  Springer  International  Publishing Switzerland.  32. Dowell,  E.H.,  Edwards,  J.W.  and  Strganac,  T.W.  (2003),  Nonlinear  Aeroelasticity. Journal of Aircraft, 40(5), 857–74.  33. Elishakoff I., Andrimasy L., Dolley M. (2009):Application and extension of  the stochastic linearization by Anh and Di Paola. Acta Mech. 204, 89–98  34. Flomenhoft,  H.I.  (1997)  The  Revolution  in  Structural  Dynamics,  Dynaflo  Press.  35. Friedmann, P.P. (1999) Renaissance of aeroelasticity and its future. Journal  of Aircraft, 36(1), 105–21.  36. Fung, Y. C. (1993), An Introduction to the Theory of Aeroelasticity, Dover  Edition.  37. Garrick,  I.E.  and  Reid,W.H.  (1981)  Historical  development  of  aircraft  flutter. Journal of Aircraft, 18(11), 897–912.  38. Glauert,  H.  (1959),  The  Elements  of  Aerofoil  and  Airscrew  Theory,  2nd  edition, Cambridge University Press.  39. Gupta.  S.G.,  M.M.  Ghonge,  P.M.  Jawandhiya  (2013),  “Review  of  Unmanned  Aircraft  System  (UAS)”,  International  Journal  of  Advanced  Research  in  Computer  Engineering  &  Technology,  Issue  4,  Vol.  2,  pp.  1646- 1658.  40. Hardin P.  and  Jensen  R.  (2011),  “Small-scale  unmanned  aerial  vehicles  in  environmental  remote  sensing:  Challenges  and  opportunities”,  GISci.  Remote  Sens., 48(1), 99-111.   111  41. Henshaw MJ, Badcock KJ, Vio GA, Allen AM, Chamberlain J, Kaynes I, et  al.  (2007),  Non-linear  aeroelastic  prediction  for  aircraft  applications.  Prog  Aerosp Sci;43(4–6):65–137.  42. Hodges, D.H. and Pierce, G.A. (2002)  Introduction to Structural Dynamics  and Aeroelasticity, Cambridge University Press.  43. Ira  H.  Abbott  and  Albert  E.  Von  Doenhoff  (1951),  "Theory  of  Wing  Sections", Dover Publishing, New York,.  44. James. R.M. (1997), “The theory and design of two-airfoil lifting systems”,  Computer Methods in Applied Mechanics and Engineering, vol. 10, pp. 13-43,.  45. Jin  Y.,  Yuan  X.,  Shin  B.R.  (2002),  “Numerical  Analysis  of  the  Airfoil’s  Fluid-Structure Interaction Problems at Large at Large Mean Incidence Angle”,  Proc. ICCFD, Sydney, Australia, 15–19 July.  46. Katz,  J.  and  Plotkin,  A.  (2001)  Low  Speed  Aerodynamics,  2nd  edn,  Cambridge University Press.  47. Ko J, Strganac TW, Junkins JL (2002). Structured model reference adaptive  control for a wing section with structural nonlinearity. J Vib Control; 8(5):553– 557.  48. Krylov  N.  M.,  Bogolyubov  N.  N.  (1937):  Introduction  to  non-linear  mechanics (in Russian). Kiev: Publisher AN SSSR.  49. Langley R S.  (1988): An  investigation of multiple solutions yielded by the  equivalent linearization method. J. Sound Vib.  127:271-281.  50. Lee  B.H.K.,  Gong  L.,  Wong  Y.S.  (1997),  Analysis  and  computation  of  nonlinear dynamic response of a two-degree-freedom system and its application  in aeroelasticity, Journal of Fluids and Structures, 11:225-246.  51. Lee B.H.K., Price S.J., Wong Y.S. (1999), Nonlinear aeroelastic analysis of  airfoils: bifurcation and chaos, Progress in Aerospace Sciences, 35:205-334.  52. Lee BHK, Liu L, Chung KW. (2005) Airfoil motion in subsonic flow with  strong cubic nonlinear restoring forces. J Sound Vib;281(3–5):699–717.   112  53. Li  DC,  Guo  SJ,  Xiang  JW  (2011).  Adaptive  control  of  a  nonlinear  aeroelastic system. Aerosp Sci Technol;15(5):343–352.  54. Li DC, Guo SJ, Xiang JW, Di Matteo N. (2010), Control of an aeroelastic  system  with  control  surface  nonlinearity.  In:  Proceedings  of  51st  AIAA/ASME/ASCE/AHS/ASC  structures,  structural  dynamics,  and  materials  conference.  55. Li DC, Guo SJ, Xiang JW. (2012) Study of the conditions that cause chaotic  motion  in  a  two-dimensional  airfoil  with  structural  nonlinearities  in  subsonic  flow. J Fluids Struct;33:109–26.  56. Liu J.K., Zhao L.C. (1992), Bifurcation analysis of airfoils in incompressible  flow, Journal of Sound and Vibration,154:117-124.  57. Liu  L,  Dowell  EH,  Thomas  JP.  (2007)  A  high  dimensional  harmonic  balance approach for an aeroelastic airfoil with cubic restoring forces. J Fluids  Struct;23(3):351–63.  58. Liu L, Dowell EH. (2005) Harmonic balance approach for an airfoil with a  freeplay control surface. AIAA J;43(4):802–15.  59. Liu  L,  Wong  YS,  Lee  BHK.  (2000)  Application  of  the  center  manifold  theory in nonlinear aeroelasticity. J Sound Vib;234(4):641–59.  60. Liu L, Wong YS, Lee BHK. (2002) Non-linear aeroelastic analysis using the  point  transformation method, part 1:  freeplay model.  J Sound Vib;253(2):447– 69.  61. Liu. P, Albert Y. Che, Yin-Nan Huang, Jen-Yu Han,  Jihn-Sung Lai, Shih- Chung Kang, Tzong-Hann Wu, Ming-Chang Wen and Meng-Han Tsai (2014),  “A  review  of  rotorcraft  Unmanned  Aerial  Vehicle  (UAV)  developments  and  applications in civil engineering”, Smart Structures and Systems, Vol. 13, No. 6,  pp. 1065-1094.  62. Livne, E. (2003) Future of airplane aeroelasticity. Journal of Aircraft, 40(6),  1066–92.   113  63. Ma.  K.,  H.  Wen,  T.  Hu,  T.D.  Topping,  D.  Isheim,  D.N.  Seidman,  E.J.  Laverni,  J.  M.  Schoenung  (2014),  “Mechanical  behavior  and  strengthening  mechanisms in ultrafine grain precipitation-strengthened aluminum alloy”, Acta  Materialia, Vol. 62, pp. 141–155.  64. Mohamed.  G.(2000),  “Flow  Control:  Passive,  Active  and  Reactive  Flow  Management”, Cambrigde University Press.  65. Nguyen  Dong  Anh,  Nguyen  Xuan  Nguyen,  Nguyen  Hoang  Quan  (2014),  Global-local approach  to  the design of dynamic vibration absorber for damped  structures, Journal of Vibration and Control, doi: 10.1177/1077546314561282  66. Petrila. T and Trif. D (2005), “Basics of fluid mechanics and introduction to  computational fluid dynamics”, Springer-Verlag.  67. Pines.  S  (1958),  An  elementary  explanation  of  the  flutter  mechanism.  In:  Proceedings nature  specialists  meeting  on  dynamics  and  aeroelasticit,  Institute  of theAeronautical Sciences, Ft.Worth, Texas, pp 52–58  68. Platanitis G, Strganac TW (2004). Control of a nonlinear wing section using  leading and trailing-edge surfaces. J Guidance Control Dyn;27(1):52–58.  69. Prabhakar  A.  and  Ohri  A.  (2013),  “CFD  Analysis  on  MAV  NACA  2412  Wing  in  High  Lift  Take-Off  Configuration  for  Enhanced  Lift  Generation”,  J  Aeronaut Aerospace Eng., 2: 125. doi:10.4172/2168-9792.1000125,.  70. Proppe, C., Pradlwarter, H.J., Schüller, G.I. (2003): Equivalent linearization  and  Monte-Carlo  simulation  in  stochastic  dynamics.  J.  Probab.  Eng.  Mech.  18(1), 1–15  71. R. Eppler (1990), “Airfoil Design and Data”, Springer Berlin Heidelberg.  72. Rao,  S.  S.  (2010):  Mechanical  Vibrations,  5th  ed.,  Addison-Wesley,  Reading, MA  73. Rathinam,S., Kim Z.W. and Sengupta R. (2008), “Vision-based monitoring  of  locally  linear  structures  using  an  unmanned  aerial  vehicle”,  Journal  of  Infrastructure Systems, 14(1), 52-63.   114  74. Roberts,  J.B.,  Spanos,  P.D.  (1990):  Random  Vibration  and  Statistical  Linearization. Wiley, New York.  75. Rumsey.C. L., Ying. S. X., et  al  (2002),  “A CFD Prediction of High Lift:  review  of  present  CFD  capability”,  Progress  in  Aerospace  Sciences,  vol.  38,  issue 2.  76. Shahrzad  P.,  Mahzoon  M.  (2002),  Limit  cycle  flutter  of  airfoils  in  steady  and unsteady flows, Journal of Sound and Vibration, 256:213-225.  77. Socha,  L.  (2008):  Linearization  methods  for  stochastic  dynamic  system.  Lecture Notes in Physics. Springer, Berlin.  78. Strganac. T.W., J. Ko, D.E. Thompson (2000), Identification and control of  limit cycle oscillations in aeroelastic systems, Journal of Guidance, Control, and  Dynamics 23 (6) 1127–1133.  79. T.  Petrila  and  D.  (2005),  “Basics  of  fluid  mechanics  and  introduction  to  computational fluid dynamics”, Springer-Verlag.  80. Tang DM, Dowell EH. (2010) Aeroelastic airfoil with free play at angle of  attack with gust excitation. AIAA J;48(2):427–42.  81. Tang  DM,  Henri  PG,  Dowell  EH.  (2004a)  Study  of  airfoil  gust  response  alleviation  using  an  electro-magnetic  dry  friction  damper,  part  1:  theory.  J  Sound Vib;269(3):853–74.  82. Tang  DM,  Henri  PG,  Dowell  EH.  (2004b)  Study  of  airfoil  gust  response  alleviation using an electro-magnetic dry friction damper, part 2: experiment. J  Sound Vib;269(3):875–97.  83. Thomas J. Muller and James D. DeLaurier (2003), “Aerodynamics of small  vehicles”, Annual Review in Fluid Mechanics, 35:89–111.  84. Turner D., Lucieer A. and Watson C. (2012), “An automated technique for  generating  georectified  mosaics  from  ultra-high  resolution  unmanned  aerial  vehicle  (UAV)  imagery,  based  on  structure  from  motion  (SfM) point  clouds”,  Remote Sensing, 4(5), 1392-1410.   115  85. Venter, G. (2010). “Review of optimization techniques, In: Encyclopedia of  aerospace engineering”, Wiley & Sons Ltd.  86. Wei,  X  and  Mottershead,  JE  (2014)  Aeroelastic  systems  with  softening  nonlinearity. AIAA Journal, 52 (9). pp. 1922-1927.  87. Wright,  J.  R.,  and  Cooper,  J.  E.  (2007).  Introduction  to  Aircraft  Aeroelasticity and Loads: John Wiley & Sons.  88. Yang Y.R.  (1995), KBM method of analyzing limit cycle flutter of a wing  with an external store and comparison with a wind-tunnel test, Journal of Sound  and Vibration,187:271-280.  89. Yang Z.C., Zhao L.C. (1988), Analysis of limit cycle flutter of an airfoil in  incompressible flow, Journal of Sound and Vibration, 123:1-13.  90. Zhao YH, Hu HY (2004). Aeroelastic analysis of a non-linear airfoil based  on unsteady vortex lattice model. J Sound Vib;276(3–5): 491–510.   116  PHỤ LỤC 1. Ma trận Sylvester và kết thức Trong toán học, ma trận Sylvester là ma trận gắn liền với hai đa thức một biến.  Các phần tử của ma trận Sylvester là các hệ số của đa thức. Định thức của ma  trận Sylvester là kết thức (resultant) của hai đa thức. Kết thức của 2 đa thức sẽ  bằng 0 khi 2 đa thức có nghiệm chung. Ma trận Sylvester được định nghĩa như  sau:  Cho p và q là hai đa thức khác không, có bậc tương ứng là m và n. Khi đó:      2 0 1 2 2 0 1 2 ... ... m m n n p z p p z p z p z q z q q z q z q z           Ma trận Sylvester tương ứng với p và q sẽ có kích cỡ (m+n)(m+n) và thu được  theo cách sau:  - Hàng đầu tiên là:  1 1 0... 0 ... 0m mp p p p    -  Hàng  thứ  hai  giống  hàng  đầu  tiên  nhưng  được  dịch  một  cột  sang  bên  phải,  phần tử đầu tiên của hàng bằng 0.  - n – 2 hàng tiếp theo thu được bằng cách tương tự, dịch các hệ số sang bên phải  một cột trong mỗi lần và cho các thành phần khác của hàng bằng 0.  - Hàng n + 1 bằng:  1 1 0... 0 ... 0n nq q q q   - Các hàng tiếp theo thu được bằng cách như đã nói ở trên.  Xét ví dụ 2 đa thức bậc 2 có dạng:      2 0 1 2 2 0 1 2 p z p p z p z q z q q z q z       Ma trận Sylvester tương ứng với 2 đa thức này là:  2 1 0 2 1 0 , 2 1 0 2 1 0 0 0 0 0 p q p p p p p p S q q q q q q              117  và kết thức của 2 đa thức có dạng:                            2 , 2 0 2 0 0 1 1 0 2 1 2 1p qR p q q p p q p q p q q p       2. Đoạn mã MATLAB tính ví dụ mục 3.4 syms t x real % cac bien ky hieu n=1; % gia tri cua n v1=int((cos(2*pi*t))^(2*n+2),0,1); % ham trung gian v2=int((cos(2*pi*t))^(4*n+2),0,1); % ham trung gian v3=int((cos(2*pi*t))^2,0,1); % ham trung gian r2=v1^2/v2/v3; % he so tuong quan % tan so chinh xac om_e=double(2*pi/4/sqrt(n+1)/int(1/sqrt(1-x^(2*n+2)),0,1)); % tan so theo tieu chuan kinh dien om_kd=double(sqrt(v1/v3)); % tan so theo tieu chuan doi ngau thong thuong om_dn=double(sqrt(1/(2-r2)*v1/v3)); % tan so theo tieu chuan doi ngau cai tien 1 om_dn1=double(sqrt((3-r2)/2/(2-r2)*v1/v3)); % tan so theo tieu chuan doi ngau cai tien 2 om_dn2=double(sqrt((1/r2+2*(1-r2)/r2^2*log(1-r2/2))*v1/v3)); % tan so theo tieu chuan doi ngau cai tien 3 om_dn3=double(sqrt(3/(4-r2)*v1/v3)); % hien thi ket qua [om_e om_kd om_dn om_dn1 om_dn2 om_dn3] 3. Đoạn mã MATLAB tính ví dụ mục 3.5 om=1;h=0.5;gm=0.1;s=1; syms x real % bien ky hieu % dich chuyen chinh xac temp1=double(int(exp(-4*h/s*(1/2*om^2*x^2+1/4*gm*x^4)),-inf,inf)); temp2=double(int(x^2*exp(-4*h/s*(1/2*om^2*x^2+1/4*gm*x^4)),-inf,inf)); x2_e=temp2/temp1; % dich chuyen theo tieu chuan kinh dien temp1=roots([3*gm,om^2,-s/4/h]); x2_kd=temp1(find(temp1>0)); % dich chuyen theo tieu chuan doi ngau thong thuong temp1=roots([15/7*gm,om^2,-s/4/h]); x2_dn=temp1(find(temp1>0)); % dich chuyen theo tieu chuan doi ngau cai tien 1 temp1=roots([18/7*gm,om^2,-s/4/h]); x2_dn1=temp1(find(temp1>0)); % dich chuyen theo tieu chuan doi ngau cai tien 2 temp1=roots([(log(7/10)*20/3+5)*gm,om^2,-s/4/h]); x2_dn2=temp1(find(temp1>0)); % dich chuyen theo tieu chuan doi ngau cai tien 3 temp1=roots([45/17*gm,om^2,-s/4/h]); x2_dn3=temp1(find(temp1>0)); % hien thi ket qua [x2_e x2_kd x2_dn x2_dn1 x2_dn2 x2_dn3]  118  4. Đoạn mã MATLAB lập phương trình với vận tốc tới hạn Hàm sau mô tả phương trình để xác định vận  tốc  tới hạn u và  tính tần số dao  động . Nếu biến đánh dấu flag có giá  trị 1 thì cho ra phương trình để  tìm u.  Nếu đánh dấu flag có giá trị 2 thì tính  theo các tham số và theo u.  function f=phuongtrinh(u,ka,ca,a,b,mt,mw,Ia,xa,kh,ch,clalp,cmalp,clbta,cmbta,s,ro,kp ,ki,kd,flag) m1=mt;m2=mw*xa*b;m3=Ia; k1=-kh;k2=-ro*u^2*b*s*(clalp+clbta*kp);k3=-ro*u^2*b*s*clbta*ki;k4=-ch- ro*u*b*s*clalp; k5=-ro*u*b*s*(clalp*(1/2-a)*b+clbta*u*kd);k6=0;k7=- ka+ro*u^2*b^2*s*(cmalp+clbta*kp); k8=ro*u^2*b^2*s*clbta*ki;k9=ro*u*b^2*s*cmalp;k10=- ca+ro*u*b^2*s*(cmalp*(1/2-a)*b+u*cmbta*kd); if flag==1 p2=(m1*m3-m2^2);q2=(-m1*k10+m2*k9+m2*k5-k4*m3); p1=-(m2*k6+m2*k2+k4*k10-k5*k9-k1*m3-m1*k7);q1=-(-m1*k8+k1*k10+m2*k3-k6*k5- k2*k9+k7*k4); p0=(k1*k7-k6*k2-k3*k9+k8*k4);q0=(k1*k8-k6*k3); f=(p2*q0-q2*p0)^2+(p0*q1-p1*q0)*(p2*q1-q2*p1); elseif flag==2 temp1=(-k1*k8+k6*k3)/(-m2*k9+m3*k4-m2*k5+m1*k10)-(k6*k2+k3*k9-k8*k4- k1*k7)/(m2^2-m1*m3); temp2=(k4*k10+m2*k2+m2*k6-m1*k7-k5*k9-k1*m3)/(m2^2-m1*m3)-(k7*k4+k1*k10- k6*k5-m1*k8+m2*k3-k2*k9)/(-m2*k9+m3*k4-m2*k5+m1*k10); f=sqrt(temp1/temp2); end Đoạn mã sau thể hiện một đoạn chạy thử với một tập số liệu đầu vào để tính U  và   function chaythu a=-0.6847;b=0.135;mt=12.387;mw=2.049;Ia=0.0558;xa=0.3314;kh=2884.4; ch=27.43;clalp=6.28;cmalp=-1.160;s=0.6;ro=1.225;clbta=3.358;cmbta=-0.635; kp=0;ki=0;kd=0;ka=6.833;ca=0.036; u=fzero(@phuongtrinh,[0 15],[],ka,ca,a,b,mt,mw,Ia,xa,kh,ch,clalp,cmalp,clbta,cmbta,s,ro,kp,ki,kd,1) om=phuongtrinh(u,ka,ca,a,b,mt,mw,Ia,xa,kh,ch,clalp,cmalp,clbta,cmbta,s,ro,k p,ki,kd,2) 5. Đoạn mã MATLAB giải phương trình vi phân tìm vận tốc tới hạn Hàm sau mô tả phương trình vi phân.  function dx=ptvp(t,x,u,ka1,ka2,ka3,ka4,ka5,ca,a,b,mt,mw,Ia,xa,kh,ch,clalp,cmalp,clbt a,cmbta,s,ro,kp,ki,kd) alp=x(2); ka=ka1+ka2*alp+ka3*alp^2+ka4*alp^3+ka5*alp^4; m1=mt;m2=mw*xa*b;m3=Ia; k1=-kh;k2=-ro*u^2*b*s*(clalp+clbta*kp);k3=-ro*u^2*b*s*clbta*ki;k4=-ch- ro*u*b*s*clalp; k5=-ro*u*b*s*(clalp*(1/2-a)*b+clbta*u*kd);k6=0;k7=- ka+ro*u^2*b^2*s*(cmalp+clbta*kp); k8=ro*u^2*b^2*s*clbta*ki;k9=ro*u*b^2*s*cmalp;k10=- ca+ro*u*b^2*s*(cmalp*(1/2-a)*b+u*cmbta*kd); SystemMatrix=blkdiag(eye(3),[m1 m2;m2 m3])\[0 0 0 1 0;0 0 0 0 1;0 1 0 0 0;k1 k2 k3 k4 k5;k6 k7 k8 k9 k10]; dx=SystemMatrix*x;  119  Hàm sau mô tả đoạn giải lặp phương trình vi phân để tìm vận tốc tới hạn.  % cac tham so cua canh a=-0.6847;b=0.135;mt=12.387;mw=2.049;Ia=0.0558;xa=0.3314;kh=2884.4; ch=27.43;clalp=6.28;cmalp=-1.160;s=0.6;ro=1.225;clbta=3.358;cmbta=-0.635; kp=0;ki=0;kd=0;ca=0.036; ka1=6.833;ka2=9.967;ka3=667.685;ka4=26.569;ka5=-5087.931; %cac dieu kien dau a0=0.07;h0=0; % van toc ban dau initu=7.5; % do tang van toc tu tho den tinh incu_range=[1 0.5 0.1 0.05 0.01 0.005]; % thoi gian tinh va cac vi tri de xac dinh bien do dao dong Tf=120;Tcut=[59 60 119 120]; % neu bien flutter=0 thi van toc con duoi toi han flutter=0; for lanlap=1:length(incu_range) u=initu; while flutter==0 [t,y]=ode45(@ptvp,[0 Tf],[h0 a0 0 0 0]',[],u,ka1,ka2,ka3,ka4,ka5,ca,a,b,mt,mw,Ia,xa,kh,ch,clalp,cmalp,clbta,cmb ta,s,ro,kp,ki,kd); vitri=find(t>Tcut(1) & t<Tcut(2)); A1=1/2*(max(y(vitri,2))-min(y(vitri,2))); V1=1/2*(max(y(vitri,5))-min(y(vitri,5))); vitri=find(t>Tcut(3) & t<Tcut(4)); A2=1/2*(max(y(vitri,2))-min(y(vitri,2))); V2=1/2*(max(y(vitri,5))-min(y(vitri,5))); if (abs(A2-A1)/A1A1) & (abs(A2/a0) > 0.1) flutter=1; else u=u+incu_range(lanlap); end end initu=u-incu_range(lanlap);flutter=0; end % cho ra ket qua bien do va tan so biendo=A2 tanso=sqrt(V2/A2)